IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41390 for trunk


Ignore:
Timestamp:
Jul 27, 2020, 11:06:35 AM (6 years ago)
Author:
eugene
Message:

add tgroups concept: periods of time which are photometric

Location:
trunk/Ohana/src/relphot
Files:
14 edited
8 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot

  • trunk/Ohana/src/relphot/Makefile

    r39636 r41390  
    3232$(SRC)/BrightCatalog.$(ARCH).o   \
    3333$(SRC)/MosaicOps.$(ARCH).o       \
     34$(SRC)/TGroupOps.$(ARCH).o       \
     35$(SRC)/ReadTGroupFile.$(ARCH).o  \
    3436$(SRC)/SetSignals.$(ARCH).o      \
    3537$(SRC)/Shutdown.$(ARCH).o        \
    … …  
    8991$(SRC)/BrightCatalog.$(ARCH).o   \
    9092$(SRC)/MosaicOps.$(ARCH).o       \
     93$(SRC)/TGroupOps.$(ARCH).o       \
     94$(SRC)/ReadTGroupFile.$(ARCH).o  \
    9195$(SRC)/SetSignals.$(ARCH).o      \
    9296$(SRC)/Shutdown.$(ARCH).o        \

  • trunk/Ohana/src/relphot/include/relphot.h

    r40386 r41390  
    9090} Mosaic;
    9191
     92typedef struct {
     93  unsigned int start;
     94  unsigned int stop;
     95  short photcode;
     96  float McalPSF;
     97  float McalAPER;
     98  float dKlam;
     99  float dMcal;
     100  float dMsys;
     101  float McalChiSq;
     102  unsigned short nFitPhotom;
     103  unsigned int flags;
     104
     105  off_t Nimage;
     106  off_t NIMAGE;
     107  off_t *image;
     108
     109  off_t Nmeasure;
     110  off_t NMEASURE;
     111  off_t *measure;
     112
     113  off_t Ncatalog; // XX not needed
     114  off_t NCATALOG; // XX not needed
     115  off_t *catalog;
     116} TGroup;
     117
    92118typedef enum {
    93119  STATS_NONE,
    … …  
    145171  int Ncal;
    146172  int Nmos;
     173  int Ngrp;
    147174  int Ngrid;
    148175
    … …  
    305332int    VERBOSE_IMAGE;
    306333int    MOSAIC_ZEROPT;
     334int    TGROUP_ZEROPT;
     335int    TGROUP_FIT_AIRMASS;
    307336
    308337int    TEST_IMAGE1;
    … …  
    403432void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr));
    404433void          clean_mosaics       PROTO((void));
     434void          clean_tgroups       PROTO((void));
    405435void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
    406436int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
    … …  
    410440void          findImages          PROTO((Catalog *catalog, int Ncatalog, int doImageList));
    411441int           findMosaics         PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
     442int           findTGroups         PROTO((Catalog *catalog, int Ncatalog));
    412443
    413444void clearImages (void);
    … …  
    428459void          freeImageBins       PROTO((int Ncatalog, int doImageList));
    429460void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
     461void          freeTGroupBins      PROTO((int Ncatalog));
     462
    430463void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
    431464int           gcatalog            PROTO((Catalog *catalog, int FINAL));
    … …  
    437470float         getMgrid            PROTO((off_t meas, int cat));
    438471float         getMmos             PROTO((off_t meas, int cat));
     472float         getMgrp             PROTO((off_t meas, int cat, float airmass));
    439473float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat, dvoMagClassType class, dvoMagSourceType source));
    440474short         getUbercalDist      PROTO((off_t meas, int cat));
    … …  
    452486void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
    453487void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
     488void          initTGroups         PROTO((Image *subset, off_t *LineNumber, off_t Nsubset));
     489void          freeTGroups         PROTO((void));
     490void          initTGroupBins      PROTO((Catalog *catalog, int Ncatalog));
    454491void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
    455492RelphotMode   initialize          PROTO((int argc, char **argv));
    … …  
    459496int           liststats_init      PROTO((StatType *stats));
    460497int           liststats_irls      PROTO((StatDataSet *dataset, int Npoints, StatType *stats));
     498int           liststats_fit1d     PROTO((double *value, double *err, double *x, int Npts, StatType *stats, double *dk));
     499
     500unsigned int *ReadTGroupFile      PROTO((FILE *f, int *nelem));
     501void          loadTGroups         PROTO((char *filename));
     502void          initTGroupsMcal     PROTO((void));
     503TGroup       *getTGroupForImage   PROTO((off_t im));
     504off_t         findTGroup          PROTO((unsigned int start, int photcode));
    461505
    462506Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
    … …  
    470514void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat, int doImageList));
    471515void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat, int doMosaicList));
     516void          matchTGroups        PROTO((Catalog *catalog, off_t meas, int cat));
     517
    472518double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
    473519void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
    … …  
    491537int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
    492538void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
    493 void          setMcalFinal        PROTO((void));
     539void          setMcalFromMosaics  PROTO((void));
     540void          setMcalFromTGroups  PROTO((void));
    494541int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    495542void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    496543int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
     544int           setMgrp             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
    497545int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    498546void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
    … …  
    677725void RepairWarpMeasuresOpenLogfile ();
    678726void RepairWarpMeasuresCloseLogfile ();
     727
     728int save_images_backup (FITS_DB *db);

  • trunk/Ohana/src/relphot/src/GridOps.c

    r40291 r41390  
    289289
    290290  int **gotstar, **gridmeas;
    291   int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys;
     291  int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp;
    292292  double **A, **B, *Mjx, *Wjx;
    293   float Msys, Mcal, Mmos, Merr, Wsys;
     293  float Msys, Mcal, Mmos, Mgrp, Merr, Wsys;
    294294  double Mj, Wj;
    295295 
    … …  
    305305  }
    306306
    307   Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = 0;
     307  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0;
    308308
    309309  ALLOCATE (gotstar, int *, Ncatalog);
    … …  
    411411        if (isnan(Mmos)) {
    412412          Nmos ++;
     413          continue;
     414        }
     415
     416        // skip mosaics marked as BAD
     417        Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     418        if (isnan(Mgrp)) {
     419          Ngrp ++;
    413420          continue;
    414421        }
    … …  
    526533void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    527534
    528   int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
     535  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit;
    529536  double *list, *dlist;
    530   float Msys, Mrel, Mcal, Mmos;
     537  float Msys, Mrel, Mcal, Mmos, Mgrp;
    531538
    532539  StatType stats;
    … …  
    544551  ALLOCATE (dlist, double, Nmax);
    545552
    546   Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0;
     553  Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0;
    547554
    548555  for (i = 0; i < Ngrid; i++) {
    … …  
    568575        continue;
    569576      }
     577      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     578      if (isnan(Mgrp)) {
     579        Ngrp ++;
     580        continue;
     581      }
    570582      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    571583      if (isnan(Mrel)) {
    … …  
    580592        continue;
    581593      }
    582       list[N] = Msys - Mrel - Mcal - Mmos;
     594      list[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
    583595      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    584596      N++;
    … …  
    616628
    617629  int i, j, m, c, n, N, Narea;
    618   float Msys, Mrel, Mcal, Mmos;
     630  float Msys, Mrel, Mcal, Mmos, Mgrp;
    619631  double *xlist, *Mlist, *dlist, *ylist;
    620632  Graphdata graphdata;
    … …  
    649661      Mmos  = getMmos  (m, c);
    650662      if (isnan(Mmos)) continue;
     663      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     664      if (isnan(Mgrp)) continue;
    651665      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    652666      if (isnan(Mrel)) continue;
    … …  
    657671      xlist[N] = Xmeas[c][m];
    658672      ylist[N] = Ymeas[c][m];
    659       Mlist[N] = Msys - Mrel - Mcal - Mmos;
    660       dlist[N] = Msys - Mrel - Mcal - Mmos - gridM[i];
     673      Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
     674      dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i];
    661675      N++;
    662676    }

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r40291 r41390  
    586586
    587587  off_t i, j, m, c, n;
    588   int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys;
     588  int mark, bad, Nfew, Nbad, Nmos, Ngrp, Nrel, Ngrid, Nsys;
    589589
    590590  StatType stats;
    … …  
    612612  StatDataSet *brightStars = StatDataSetAlloc (1, Nmax);
    613613
    614   Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
     614  Nfew = Nbad = Nmos = Ngrp = Ngrid = Nrel = Nsys = 0;
    615615
    616616  int Ncalibrated = 0;
    … …  
    631631
    632632    // only freeze the gpc1 chips
    633     if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
     633    // XXX : when do we FREEZE_IMAGES in a mix of GPC1, GPC2, etc?
     634    // if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
     635    if (FREEZE_IMAGES) continue;
    634636    if (image[i].photcode == 0) continue; // skip the PHU images
    635637
    … …  
    666668          continue;
    667669      }
     670      float Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     671      if (isnan(Mgrp)) {
     672          Ngrp ++;
     673          continue;
     674      }
    668675      float Mgrid = getMgrid (m, c);
    669676      if (isnan(Mgrid)) {
    … …  
    706713
    707714    skip:
    708       psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrid + Mflat;
     715      psfStars->flxlist[Nref] = MsysPSF - MrelPSF - Mmos - Mgrp - Mgrid + Mflat;
    709716      psfStars->errlist[Nref] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    710717      psfStars->wgtlist[Nref] = 1;
    … …  
    718725
    719726      if (isfinite(MrelKron) && isfinite(MsysKron)) {
    720         kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrid + Mflat;
     727        kronStars->flxlist[Nkron] = MsysKron - MrelKron - Mmos - Mgrp - Mgrid + Mflat;
    721728        kronStars->errlist[Nkron] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    722729        kronStars->wgtlist[Nkron] = 1;
    … …  
    10231030  off_t i, j, m, c, n, N;
    10241031  double *list, *dlist;
    1025   float Mcal, Mmos, Mgrid;
     1032  float Mcal, Mmos, Mgrp, Mgrid;
    10261033
    10271034  StatType stats;
    … …  
    10501057      Mmos  = getMmos  (m, c);
    10511058      if (isnan(Mmos)) continue;
     1059      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     1060      if (isnan(Mgrp)) continue;
    10521061      Mgrid = getMgrid (m, c);
    10531062      if (isnan(Mgrid)) continue;

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r40291 r41390  
    720720}
    721721
    722 void setMcalFinal () {
     722void setMcalFromMosaics () {
    723723
    724724  off_t i, j, im, Nimage;
    … …  
    919919  int Nbad;
    920920  int Ncal;
     921  int Ngrp;
    921922  int Ngrid;
    922923  int Nrel;
    … …  
    11541155    if (fout) {
    11551156      float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
     1157      // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass);
    11561158      float Mgrid    = getMgrid (m, c);
    11571159      float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    … …  
    11751177    if (isnan(Mcal)) {
    11761178      info->Ncal++;
     1179      continue;
     1180    }
     1181    float Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     1182    if (isnan(Mgrp)) {
     1183      info->Ngrp++;
    11771184      continue;
    11781185    }
    … …  
    16921699  off_t i, j, m, c, n, N;
    16931700  double *list, *dlist;
    1694   float Mcal, Mgrid, Mrel;
     1701  float Mcal, Mgrp, Mgrid, Mrel;
    16951702  StatType stats;
    16961703
    … …  
    17161723      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
    17171724      if (isnan(Mcal)) continue;
     1725      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass);
     1726      if (isnan(Mgrp)) continue;
    17181727      Mgrid = getMgrid (m, c);
    17191728      if (isnan(Mgrid)) continue;

  • trunk/Ohana/src/relphot/src/StarOps.c

    r40291 r41390  
    105105  results->Ncal  = 0;
    106106  results->Nmos  = 0;
     107  results->Ngrp  = 0;
    107108  results->Ngrid = 0;
    108109}
    … …  
    264265  summary->Ncal  += results->Ncal ;
    265266  summary->Nmos  += results->Nmos ;
     267  summary->Ngrp  += results->Ngrp ;
    266268  summary->Ngrid += results->Ngrid;
    267269}
    … …  
    313315    SetMrelInfoAccum (&summary, &results);
    314316  }
    315   if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     317  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys);
    316318
    317319  SetMrelInfoFree (&results);
    … …  
    336338    SetMrelInfoAccum (&summary, &results);
    337339  }
    338   if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     340  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys);
    339341
    340342  SetMrelInfoFree (&results);
    … …  
    393395  // report stats & summary from the threads
    394396  for (i = 0; i < NTHREADS; i++) {
    395     if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n",
     397    if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n",
    396398             i,
    397399             threadinfo[i].summary.Ncode,
    … …  
    400402             threadinfo[i].summary.Ncal,
    401403             threadinfo[i].summary.Nmos,
     404             threadinfo[i].summary.Ngrp,
    402405             threadinfo[i].summary.Ngrid,
    403406             threadinfo[i].summary.Nsys);
    … …  
    405408    SetMrelInfoFree (&threadinfo[i].summary);
    406409  }
    407   if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
     410  if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrp: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrp, summary.Ngrid, summary.Nsys);
    408411  free (threadinfo);
    409412
    … …  
    488491        float Mmos  = getMmos  (m, i);
    489492        if (isnan(Mmos)) continue;
     493        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
     494        if (isnan(Mgrp)) continue;
    490495        float Mgrid = getMgrid (m, i);
    491496        if (isnan(Mgrid)) continue;
    … …  
    502507
    503508        // set the output calibration
    504         catalog[i].measure[m].McalPSF  = McalPSF  + Mmos + Mgrid;
    505         catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Mmos + Mgrid : McalPSF  + Mmos + Mgrid;
     509        // float Moff = Mmos + Mgrp + Mgrid;
     510        float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
     511        catalog[i].measure[m].McalPSF  = McalPSF + Moff;
     512        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
    506513
    507514        if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    … …  
    662669  off_t *ilist;
    663670  double *tlist, *list, *dlist;
    664   float Msys, Mcal, Mmos, Mgrid;
    665   int Ncal, Nmos, Ngrid, Nfew;
     671  float Msys, Mcal, Mmos, Mgrp, Mgrid;
     672  int Ncal, Nmos, Ngrp, Ngrid, Nfew;
    666673
    667674  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    689696
    690697  Ndel = Nave = 0;
    691   Ncal = Nmos = Ngrid = Nfew = 0;
     698  Ncal = Nmos = Ngrp = Ngrid = Nfew = 0;
    692699  for (i = 0; i < Ncatalog; i++) {
    693700    for (j = 0; j < catalog[i].Naverage; j++) {
    … …  
    724731          Mmos  = getMmos  (m, i);
    725732          if (isnan(Mmos)) { Nmos ++; continue; }
     733          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass);
     734          if (isnan(Mgrp)) { Ngrp ++; continue; }
    726735          Mgrid = getMgrid (m, i);
    727736          if (isnan(Mgrid)) { Ngrid ++; continue; }
    728737
    729738          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
    730           list[N] = Msys - Mcal - Mmos - Mgrid;
     739          list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
    731740          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
    732741          N++;

  • trunk/Ohana/src/relphot/src/args.c

    r40386 r41390  
    336336    remove_argument (N, &argc, argv);
    337337    USE_MCAL_PSF_FOR_STACK_APER = TRUE;
     338  }
     339
     340  // calibrate the photometry by night / time period
     341  TGROUP_ZEROPT = FALSE;
     342  char *TGROUP_FILENAME = NULL;
     343  if ((N = get_argument (argc, argv, "-tgroup"))) {
     344    TGROUP_ZEROPT = TRUE;
     345    remove_argument (N, &argc, argv);
     346    TGROUP_FILENAME = strcreate(argv[N]);
     347    remove_argument (N, &argc, argv);
     348    FREEZE_IMAGES = TRUE;
     349  }
     350  if ((N = get_argument (argc, argv, "-tgroups"))) {
     351    if (TGROUP_ZEROPT) {
     352      fprintf (stderr, "use only one of -tgroups and -tgroup (same meaning)\n");
     353      exit (2);
     354    }
     355    TGROUP_ZEROPT = TRUE;
     356    remove_argument (N, &argc, argv);
     357    TGROUP_FILENAME = strcreate(argv[N]);
     358    remove_argument (N, &argc, argv);
     359    FREEZE_IMAGES = TRUE;
     360  }
     361  TGROUP_FIT_AIRMASS = FALSE;
     362  if ((N = get_argument (argc, argv, "-tgroup-fit-airmass"))) {
     363    if (!TGROUP_ZEROPT) {
     364      fprintf (stderr, "-tgroup-fit-airmass requires -tgroup (filename) option\n");
     365      exit (2);
     366    }
     367    TGROUP_FIT_AIRMASS = TRUE;
     368    remove_argument (N, &argc, argv);
    338369  }
    339370
    … …  
    519550  if (argc != 1) relphot_usage ();
    520551
     552  if (TGROUP_ZEROPT) {
     553    loadTGroups (TGROUP_FILENAME);
     554    free (TGROUP_FILENAME);
     555  }
     556
    521557  return mode;
    522558}
    … …  
    542578  FreePhotcodeTable();
    543579 
     580  freeTGroups();
     581
    544582  ohana_memcheck (VERBOSE);
    545583  ohana_memdump (VERBOSE);

  • trunk/Ohana/src/relphot/src/liststats.c

    r40291 r41390  
    172172  for (i = ks; i < ke; i++) {
    173173    M  = SQ (value[i] - Mo);
    174     dM = SQ (dvalue[i]);
    175     X2 += M / dM;
     174    if (dvalue) {
     175      dM = SQ (dvalue[i]);
     176      X2 += M / dM;
     177    }
    176178    dS += M;
    177179  }
    178   X2 = X2 / (Nm - 1);
     180  X2 = dvalue ? X2 / (Nm - 1) : NAN;
    179181  dS = sqrt (dS / (Nm - 1));
    180182
    … …  
    187189  return (TRUE);
    188190}
     191
     192int liststats_fit1d (double *value, double *err, double *x, int Npts, StatType *stats, double *dk) {
     193
     194  double M = 0.0, MZ = 0.0, Z1 = 0.0, Z2 = 0.0, R = 0.0;
     195
     196  for (int i = 0; i < Npts; i++) {
     197    double wt = 1.0 / (err[i]*err[i]);
     198    double V = value[i];
     199    double Z = x[i];
     200    M  += V   * wt;
     201    MZ += V*Z * wt;
     202    Z1 += Z   * wt;
     203    Z2 += Z*Z * wt;
     204    R  += 1.0 * wt;
     205  }
     206
     207  double DetInv = (R*Z2 - Z1*Z1);
     208  double Det = (fabs(DetInv) < -1e6) ? 0.0 : 1.0 / DetInv; // do not allow unstable solutions
     209
     210  double zp = (M*Z2 - Z1*MZ)*Det;
     211  *dk = (MZ*R - M*Z1)*Det;
     212
     213  stats->mean = zp;
     214  stats->chisq = NAN;
     215  return TRUE;
     216}
     217
    189218
    190219// These should probably be tunable:

  • trunk/Ohana/src/relphot/src/load_images.c

    r40291 r41390  
    110110  MARKTIME("init images: %f sec\n", dtime);
    111111
     112  initTGroups (subset, LineNumber, Nsubset);
     113
    112114  return TRUE;
    113115}

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r41341 r41390  
    8686    initImageBins  (&catalog, 1, FALSE);
    8787    initMosaicBins (&catalog, 1, FALSE);
     88    initTGroupBins (&catalog, 1);
     89
    8890    initGridBins   (&catalog, 1);
    8991    TIMESTAMP(time3b);
    … …  
    9193    findImages (&catalog, 1, FALSE);
    9294    findMosaics (&catalog, 1, FALSE);
     95    findTGroups (&catalog, 1);
     96
    9397    TIMESTAMP(time4);
    9498
    … …  
    117121    TIMESTAMP(time6);
    118122
     123    // XXX freeTgroupBins here:
     124
    119125    freeImageBins (1, FALSE);
    120126    freeMosaicBins (1, FALSE);
     127    freeTGroupBins (1);
    121128    freeGridBins (1);
    122129    TIMESTAMP(time7);

  • trunk/Ohana/src/relphot/src/relphot_images.c

    r40291 r41390  
    7070
    7171    initMosaicBins (catalog, Ncatalog, TRUE);
     72    initTGroupBins (catalog, Ncatalog);
     73
    7274    initGridBins   (catalog, Ncatalog);
    7375    initMrel (catalog, Ncatalog);
    … …  
    7779
    7880    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
     81    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
     82
     83    findTGroups (catalog, Ncatalog);
    7984    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    8085
    … …  
    129134
    130135      setMmos  (catalog, FALSE, flatcorr);
     136      setMgrp  (catalog, FALSE, flatcorr);
     137
    131138      setMgrid (catalog, flatcorr);
    132139      MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    … …  
    151158      if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog);
    152159      if ((i > 8) && (i % 8 == 5)) clean_mosaics ();
     160      if ((i > 8) && (i % 8 == 5)) clean_tgroups ();
    153161      if ((i > 8) && (i % 8 == 5)) clean_images ();
    154162
    … …  
    172180    setMcal  (catalog, onlyPoorImages, flatcorr);
    173181    setMmos  (catalog, onlyPoorImages, flatcorr);
     182    setMgrp  (catalog, onlyPoorImages, flatcorr);
    174183    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    175184
    176     setMcalFinal (); // copy per-mosaic calibrations to the images
     185    setMcalFromMosaics (); // copy per-mosaic calibrations to the images
     186    setMcalFromTGroups (); // copy per-tgroup calibrations to the images
    177187
    178188    if (SAVE_IMAGE_UPDATES) {
    … …  
    249259    freeImageBins (Ncatalog, TRUE);
    250260    freeMosaicBins (Ncatalog, TRUE);
     261    freeTGroupBins (Ncatalog);
    251262    freeGridBins (Ncatalog);
    252263
    … …  
    255266    // If nloop == 0, the above pass is not performed, in which case
    256267    // the Mcal values passed to the mosaics are not returned to the images...
    257     setMcalFinal (); // copy per-mosaic calibrations to the images
     268    setMcalFromMosaics (); // copy per-mosaic calibrations to the images
     269    setMcalFromTGroups (); // copy per-tgroup calibrations to the images
    258270  }
    259271

  • trunk/Ohana/src/relphot/src/relphot_parallel_images.c

    r40291 r41390  
    116116    setMcal  (catalog, FALSE, flatcorr);
    117117    setMmos  (catalog, FALSE, flatcorr);
     118    // setMgrp  (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups
    118119    MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
    119120   
    … …  
    152153  setMcal  (catalog, onlyPoorImages, flatcorr);
    153154  setMmos  (catalog, onlyPoorImages, flatcorr);
     155  // setMgrp  (catalog, onlyPoorImages, flatcorr); XXX see above
    154156  MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    155157
    156   setMcalFinal (); // copy per-mosaic calibrations to the images
     158  setMcalFromMosaics (); // copy per-mosaic calibrations to the images
     159  // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images
    157160
    158161  share_image_mags (regionHosts, -1);

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r40373 r41390  
    221221
    222222    // ** Choose the calibration (depends on the mode : do I have an image reference or not?)
    223     float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0;
     223    float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
    224224    off_t meas = measureOffset + k;
    225225    if (getImageEntry (meas, cat) < 0) {
    … …  
    235235      Mmos  = getMmos  (meas, cat);
    236236      if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
     237      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass);
     238      if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
    237239      Mgrid = getMgrid (meas, cat);
    238240      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
    … …  
    261263      if (CHECK_VALID_MAG(Map, dMap)) {
    262264        int Nap = results->aperData[Nsec].Nlist;
    263         results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
     265        results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    264266        results->aperData[Nsec].errlist[Nap] = dMap;
    265267        results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
    … …  
    281283    if (CHECK_VALID_MAG(Mkron, dMkron)) {
    282284      int Nkron = results->kronData[Nsec].Nlist;
    283       results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
     285      results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    284286      results->kronData[Nsec].errlist[Nkron] = dMkron;
    285287      results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    … …  
    294296    if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
    295297      int Npsf = results->psfData[Nsec].Nlist;
    296       results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
     298      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrp - Mgrid; // this is consistent with PhotRel
    297299      results->psfData[Nsec].errlist[Npsf] = dMpsf;
    298300      results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
    … …  
    711713      Mmos      = getMmos  (measSeq, cat);
    712714      Mgrid     = getMgrid (measSeq, cat);
     715      // XXX can Mmos and Mgrid exist for stacks?
    713716    }
    714717   
Note: See TracChangeset for help on using the changeset viewer.