IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 20, 2005, 4:27:25 AM (21 years ago)
Author:
eugene
Message:

substantial cleanup of photometry extraction functions

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/dvo/photometry.c

    r4420 r4585  
    33/* match code to measure  */
    44# define TESTCODE(C,M) \
    5   switch (C[0].type) { \
    6   case PHOT_DEP: \
    7   case PHOT_REF: \
    8     if (C[0].code != M.source) continue; \
    9     break; \
    10   case PHOT_PRI: \
    11   case PHOT_SEC: \
    12     if (C[0].code != GetPhotcodeEquivCodebyCode (M.source)) continue; \
    13     break; \
    14   default: \
    15     break; \
    16 }
     5  if ((C != NULL) && ApplySelections[SelectionParam]) { \
     6    switch (C[0].type) { \
     7    case PHOT_DEP: \
     8    case PHOT_REF: \
     9      if (C[0].code != M.source) continue; \
     10      break; \
     11    case PHOT_PRI: \
     12    case PHOT_SEC: \
     13      if (C[0].code != GetPhotcodeEquivCodebyCode (M.source)) continue; \
     14      break; \
     15    default: \
     16      break; \
     17  } }
    1718
    1819/* exclusions based on measure.params  */
    1920# define TESTMEASURE(M) \
    20   if (TimeSelect && (M.t < tzero)) continue; \
    21   if (TimeSelect && (M.t > tend)) continue; \
    22   if (ErrSelect  && (M.dM > ErrValue)) continue; \
    23   if (TypeSelect && (TypeValue != GetMeasureTypeCode (&M))) continue; \
    24   if (iMagSelect && (PhotInst (&M) < iMagMin)) continue; \
    25   if (FlagSelect && (M.flags != FlagValue)) continue;
     21  if (ApplySelections[SelectionParam]) { \
     22    if (TimeSelect && (M.t < tzero)) continue; \
     23    if (TimeSelect && (M.t > tend)) continue; \
     24    if (ErrSelect  && (M.dM > ErrValue)) continue; \
     25    if (TypeSelect && (TypeValue != GetMeasureTypeCode (&M))) continue; \
     26    if (iMagSelect && (PhotInst (&M) < iMagMin)) continue; \
     27    if (FlagSelect && (M.flags != FlagValue)) continue; \
     28  }
    2629
    2730# define SETMAG(MOUT,MEAS,MODE) \
    … …  
    3437  if (MODE == MAG_AVE)  MOUT = PhotRel  (&MEAS, average, secfilt); \
    3538  if (MODE == MAG_REF)  MOUT = PhotCal  (&MEAS, average, secfilt, measure, GetPhotcodeEquivbyCode (MEAS.source)); \
    36   if (MagSelect && (MOUT > MagMax)) continue; \
    37   if (MagSelect && (MOUT < MagMin)) continue;
     39  if (ApplySelections[SelectionParam]) { \
     40    if (MagSelect && (MOUT > MagMax)) continue; \
     41    if (MagSelect && (MOUT < MagMin)) continue; \
     42  }
    3843
    3944/* selection criteria */
    … …  
    5156static double TypefracValue;
    5257
     58/* apply selections or not */
     59static int ApplySelections[4];
     60static int SelectionParam;
     61
    5362/* applied to Average quantities */
    5463static int PhotcodeSelect;
    … …  
    6978static int TimeFormat;
    7079
    71 /* db image table */
    72 static Image *image = NULL;
    73 static int *subset = NULL;
    74 static int Nimage = 0;
    75 static int Nsubset = 0;
    76 static Coords mosaic;
     80int GetTimeSelection (time_t *tz, time_t *te) {
     81  *tz = tzero;
     82  *te = tend;
     83  return (TimeSelect);
     84}
    7785
    7886int GetPhotcodeInfo (char *string, PhotCode **Code, int *Mode) {
    … …  
    156164}
    157165 
     166int SetSelectionParam (int param) {
     167  SelectionParam = param;
     168  return (TRUE);
     169}
     170
     171int GetSelectionParam () {
     172  return (SelectionParam);
     173}
     174
    158175int GetMagMode (char *string) {
    159176
    … …  
    196213    fprintf (stderr, "value may be one of the following:\n");
    197214    fprintf (stderr, " ra dR dec dD mag dmag Mrel Mcal photcode time fwhm dophot xccd yccd xmosaic ymosaic flags\n");
     215    fprintf (stderr, "value may also be a valid photcode\n");
     216    fprintf (stderr, "photcodes or 'mag' may have optional magnitude type: mag,[Minst, Mcat, Msys, Mrel, Mcal]\n");
    198217  }
    199218  return (param);
    200219}
    201220 
    202 int TestPhotSelections (PhotCode **code) {
    203 
    204   int NeedPhotcode;
    205 
    206   /* if i've supplied a photcode (code != NULL), i'm not allowed to restrict it */
    207   if ((code[0] != NULL) && PhotcodeSelect) return (FALSE);
    208 
    209   /* if I have an average or ensemble restriction, I must have a photcode */
    210   if (ChiSelect || NphotSelect || NcodeSelect || FWHMSelect) {
    211     if (code[0] != NULL) return (TRUE);
    212     if (PhotcodeSelect) {
    213       code[0] = PhotcodeValue;
    214       return (TRUE);
    215     }
    216     return (FALSE);
    217   }
    218 
    219   /* if I have a photcode restriction, apply it to the code */
    220   if (PhotcodeSelect) {
    221     code[0] = PhotcodeValue;
    222     return (TRUE);
    223   }
    224   return (TRUE);
    225 }
    226 
    227221int GetAverageParam (char *parname) {
    228222
    … …  
    233227  if (!strcasecmp (parname, "dec"))   param = AVE_DEC;
    234228  if (!strcasecmp (parname, "dmag"))  param = AVE_dMAG;
     229  if (!strcasecmp (parname, "mag"))   param = AVE_MAG;
    235230  if (!strcasecmp (parname, "Nmeas")) param = AVE_NMEAS;
    236231  if (!strcasecmp (parname, "Nmiss")) param = AVE_NMISS;
    … …  
    251246}
    252247
     248/* I've set some selections - if these require a photcode, check if I set one */
     249int TestPhotSelections (PhotCode **code, int param) {
     250
     251  int NeedPhotcode;
     252
     253  /* if i've supplied a photcode (code != NULL), i'm not allowed to restrict it */
     254  if (code[0] != NULL) {
     255    if (PhotcodeSelect) {
     256      fprintf (stderr, "photcode selection rules violated: cannot restrict photcode with a photcode\n");
     257      return (FALSE);
     258    } else {
     259      return (TRUE);
     260    }
     261  }
     262
     263  /* if I have an average or ensemble restriction, I must have a photcode */
     264  NeedPhotcode = FALSE;
     265  NeedPhotcode |= ChiSelect;
     266  NeedPhotcode |= NphotSelect;
     267  NeedPhotcode |= NcodeSelect;
     268  NeedPhotcode |= ErrSelect;
     269  NeedPhotcode |= TypeSelect;
     270  NeedPhotcode |= TypefracSelect;
     271 
     272  /* for measure tests, supply MEAS_ZERO */
     273  NeedPhotcode |= (param == AVE_Xm);
     274  NeedPhotcode |= (param == AVE_MAG);
     275  NeedPhotcode |= (param == AVE_dMAG);
     276  NeedPhotcode |= (param == AVE_TYPE);
     277  NeedPhotcode |= (param == AVE_NPHOT);
     278  NeedPhotcode |= (param == AVE_NCODE);
     279
     280  if (NeedPhotcode || PhotcodeSelect) {
     281    if (!PhotcodeSelect) return (FALSE);
     282    code[0] = PhotcodeValue;
     283    if (code[0][0].type == PHOT_DEP) return (FALSE);
     284    if (code[0][0].type == PHOT_REF) return (FALSE);
     285  }
     286  return (TRUE);
     287}
     288
    253289void GetAverageParamHelp () {
    254290  fprintf (stderr, "value may be one of the following:\n");
    … …  
    259295}
    260296
    261 /* load images based on parameters and region, etc */
    262 int SetImageSelection (int mode, Graphdata *graphsky, int RegionSelect) {
    263 
    264   image = NULL;
    265   subset = NULL;
    266  
    267   switch (mode) {
    268     case MEAS_XCCD:
    269     case MEAS_YCCD:
    270       break;
    271     case MEAS_XMOSAIC:
    272     case MEAS_YMOSAIC:
    273       /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */
    274       mosaic.crpix1 = mosaic.crpix2 = 0.0;
    275       mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600;
    276       mosaic.pc1_1  = mosaic.pc2_2  = 1.0;
    277       mosaic.pc1_2  = mosaic.pc2_1  = 0.0;
    278       mosaic.Npolyterms = 0;
    279       strcpy (mosaic.ctype, "RA---SIN");
    280       break;
    281     default:
    282       return (TRUE);
    283   }
    284 
    285   if ((image = LoadImages (&Nimage)) == NULL) return (FALSE);
    286   BuildChipMatch (image, Nimage);
    287   image_subset (image, Nimage, &subset, &Nsubset, graphsky, RegionSelect, tzero, (double) tend - tzero, TimeSelect);
    288   sort_image_subset (image, subset, Nsubset);
    289   return (TRUE);
    290 }
    291 
    292 /* free loaded images */
    293 void FreeImageSelection () {
    294   if (image != NULL) free (image);
    295   if (subset != NULL) free (subset);
    296   image = NULL;
    297   subset = NULL;
    298   return;
    299 }
    300 
    301297/* (re)load photcodes from photcode table */
    302298int InitPhotcodes () {
    … …  
    318314/* remove standard photometry filtering options, set selections */
    319315/* not all functions respect all selections... */
    320 int SetPhotSelections (int *argc, char **argv) {
    321 
    322   int N;
     316int SetPhotSelections (int *argc, char **argv, int Nparams) {
     317
     318  int i, N;
    323319  double trange;
    324320
    … …  
    367363  }
    368364
    369   /* select on value of flag (MEASURE ONLY) */
    370   ChiSelect = FALSE;
    371   if (N = get_argument (*argc, argv, "-chisq")) {
    372     ChiSelect = TRUE;
    373     remove_argument (N, argc, argv);
    374     ChiLimit = atof (argv[N]);
     365  /* select on value of Chisq (AVERAGE ONLY) */
     366  for (i = 0; i < 4; i++) ApplySelections[i] = TRUE;
     367  if (N = get_argument (*argc, argv, "-apply")) {
     368    remove_argument (N, argc, argv);
     369    if (strlen(argv[N]) != Nparams) {
     370      fprintf (stderr, "-apply selection must define all parameter choices\n");
     371      return (FALSE);
     372    }
     373    for (i = 0; i < Nparams; i++) {
     374      if (toupper(argv[N][i]) == 'Y') {
     375        ApplySelections[i] = TRUE;
     376      } else {
     377        ApplySelections[i] = FALSE;
     378      }
     379    }
    375380    remove_argument (N, argc, argv);
    376381  }
    … …  
    384389    if (PhotcodeValue == NULL) {
    385390      fprintf (stderr, "photcode not found in photcode table\n");
    386       PhotcodeSelect = FALSE;
     391      return (FALSE);;
    387392    }
    388393    remove_argument (N, argc, argv);
    … …  
    422427  }
    423428
     429  /* select on value of Chisq (AVERAGE ONLY) */
     430  ChiSelect = FALSE;
     431  if (N = get_argument (*argc, argv, "-chisq")) {
     432    ChiSelect = TRUE;
     433    remove_argument (N, argc, argv);
     434    ChiLimit = atof (argv[N]);
     435    remove_argument (N, argc, argv);
     436  }
     437
    424438  /* select on measurement type: 1,2,3 (AVERAGE ONLY) */
    425439  TypeSelect = FALSE;
    … …  
    508522  ALLOCATE (list, double, NLIST);
    509523
    510   /* chisq, fwhm, Nphot, Ncode */
    511   /* if (!TestAverage (code, average, secfilt, measure)) return (list); */
     524  /* check selections based on averages & ensembles: chisq, Nphot, etc */
     525  if (!TestAverage (code, average, secfilt, measure)) return (list);
    512526
    513527  /* look for measures */
    514528  for (i = 0; i < average[0].Nm; i++) {
    515     if (code != NULL) {
    516       TESTCODE (code, measure[i]);
    517     }
    518     /* exclusions based on measure.params  */
    519     TESTMEASURE (measure[i]);
    520 
    521     /* need to have the mag value to test MagSelect */
    522     SETMAG (M, measure[i], mode);
     529    TESTCODE (code, measure[i]);  /* skip measurements not matching photcode */
     530    TESTMEASURE (measure[i]);     /* exclusions based on measure.params  */
     531    SETMAG (M, measure[i], mode); /* set appropriate magnitude (also does MagSelect) */
    523532
    524533    /* assign value */
    525     value = 0;
    526     switch (param) {
    527       case MEAS_MAG: /* magnitudes are already determined above */
    528         value = M;
    529         break;
    530       case MEAS_RA: /* OK */
    531         value = average[0].R - measure[i].dR / 360000.0;
    532         break;
    533       case MEAS_DEC: /* OK */
    534         value = average[0].D - measure[i].dD / 360000.0;
    535         break;
    536       case MEAS_dMAG: /* OK */
    537         value = 0.001*measure[i].dM;
    538         break;
    539       case MEAS_AIRMASS: /* OK */
    540         value = 0.001*measure[i].airmass;
    541         break;
    542       case MEAS_EXPTIME: /* OK */
    543         value = pow (10.0, measure[i].dt * 0.0004);
    544         break;
    545       case MEAS_PHOTCODE: /* OK */
    546         value = measure[i].source;
    547         break;
    548       case MEAS_TIME: /* OK */
    549         value = TimeValue (measure[i].t, TimeReference, TimeFormat);
    550         break;
    551       case MEAS_dR: /* OK */
    552         value = 0.01*measure[i].dR;
    553         break;
    554       case MEAS_dD: /* OK */
    555         value = 0.01*measure[i].dD;
    556         break;
    557       case MEAS_FWHM: /* OK */
    558         value = 0.01*measure[i].FWx;
    559         break;
    560       case MEAS_DOPHOT: /* OK */
    561         value = measure[i].dophot;
    562         break;
    563       case MEAS_FLAGS: /* ? */
    564         value = measure[i].flags;
    565         break;
    566       case MEAS_XCCD: /* OK */
    567         ra  = average[0].R - measure[i].dR / 360000.0;
    568         dec = average[0].D - measure[i].dD / 360000.0;
    569         m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
    570         if (m == -1) break;
    571         if (!FindMosaicForImage (image, Nimage, m)) break;
    572         RD_to_XY (&x, &y, ra, dec, &image[m].coords);
    573         value = x;
    574         break;
    575       case MEAS_YCCD: /* OK */
    576         ra  = average[0].R - measure[i].dR / 360000.0;
    577         dec = average[0].D - measure[i].dD / 360000.0;
    578         m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
    579         if (m == -1) break;
    580         if (!FindMosaicForImage (image, Nimage, m)) break;
    581         RD_to_XY (&x, &y, ra, dec, &image[m].coords);
    582         value = y;
    583         break;
    584       case MEAS_XMOSAIC: /* OK */
    585         ra  = average[0].R - measure[i].dR / 360000.0;
    586         dec = average[0].D - measure[i].dD / 360000.0;
    587         m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
    588         if (m == -1) break;
    589         /* this should use the mosaic assoicated with the image */
    590         /* the solution here depends on image.coords having crref1,2 set to match the boresite center */
    591         /* XXX set the mosaic.coords to have a unity WRP transformation in which x,y = ra,dec */
    592         mosaic.crval1 = image[m].coords.crval1;
    593         mosaic.crval2 = image[m].coords.crval2;
    594         RD_to_XY (&x, &y, ra, dec, &mosaic);
    595         value = x;
    596         break;
    597       case MEAS_YMOSAIC: /* OK */
    598         ra  = average[0].R - measure[i].dR / 360000.0;
    599         dec = average[0].D - measure[i].dD / 360000.0;
    600         m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
    601         if (m == -1) break;
    602         /* same comments as above */
    603         mosaic.crval1 = image[m].coords.crval1;
    604         mosaic.crval2 = image[m].coords.crval2;
    605         RD_to_XY (&x, &y, ra, dec, &mosaic);
    606         value = y;
    607         break;
    608     }
    609     list[Nlist] = value;
     534    list[Nlist] = GetMeasure (param, &average[0], &measure[i], M);
    610535    Nlist ++;
    611536  }
    … …  
    622547  value = NO_MAG;
    623548
    624   /* need valid code if param is: */
    625   switch (param) {
    626     case AVE_MAG:
    627     case AVE_dMAG:
    628     case AVE_Xm:
    629     case AVE_TYPE:
    630     case AVE_NCODE:
    631     case AVE_NPHOT:
    632       if (code == NULL) return (NO_MAG);
    633       break;
    634     default:
    635       break;
    636   }
    637 
    638   /* chisq, fwhm, Nphot, Ncode */
     549  /* this function requires code set for certain value of param. 
     550     use TestPhotSelectionsAverage to validate code/param choices */
     551
     552  /* filter by average quantities (eg, chisq, Nphot, etc) */
    639553  if (!TestAverage (code, average, secfilt, measure)) return (NO_MAG);
    640554
    … …  
    728642}
    729643
    730 /* determine the representative dophot type for this photcode */
     644/* determine the representative dophot type for this photcode (must be PRI/SEC) */
    731645int DetermineTypeCode (Average *average, Measure *measure, int code) {
    732646
    … …  
    785699}
    786700
     701/* test if this average object meets the specified selection criteria.
     702   for photcode-dependent quantities, only test for PRI/SEC photcodes */
    787703int TestAverage (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {
    788704
    789   int i, Nm, Type;
     705  int i, Nm, Type, Select;
    790706  double fwhm, typefrac, dM, Xm;
    791707
    792   /* exclusions based on average.params  */
    793   if (ChiSelect) {
    794     /* some average.params need a valid photcode for reference */
    795     if (code == NULL) return (FALSE);
    796     Xm = PhotXm (code, average, secfilt);
    797     if (Xm == -1) return (FALSE);
    798     if (Xm > ChiLimit) return (FALSE);
    799   }
    800  
    801   /* for ErrSelect, check average errors */
    802   if (ErrSelect) {
    803     if (code == NULL) return (FALSE);  /* need a valid code */
    804     dM = iPhotdM (code, average, secfilt);
    805     if (dM > ErrValue) return (NO_MAG);
    806   }
    807  
     708  if (!ApplySelections[SelectionParam]) return (TRUE);
     709
    808710  /* pass objects with more than FWHMfrac points with FWHM above / below FWHMvalue */
    809711  if (FWHMSelect) {
    … …  
    827729  }
    828730
     731  /* all selections below require a valid photcode */
     732  Select = ChiSelect || ErrSelect || NcodeSelect || NphotSelect || TypeSelect || TypefracSelect;
     733  if (!Select) return (TRUE);
     734
     735  /* must have a valid code of some kind */
     736  if (code == NULL) return (FALSE);
     737
     738  /* only PRI/SEC photcodes apply the filter */
     739  if (code[0].type == PHOT_DEP) return (TRUE);
     740  if (code[0].type == PHOT_REF) return (TRUE);
     741
     742  /* exclusions based on average.params  */
     743  if (ChiSelect) {
     744    Xm = PhotXm (code, average, secfilt);
     745    if (Xm == -1) return (FALSE);
     746    if (Xm > ChiLimit) return (FALSE);
     747  }
     748 
     749  /* for ErrSelect, check average errors */
     750  if (ErrSelect) {
     751    dM = iPhotdM (code, average, secfilt);
     752    if (dM > ErrValue) return (NO_MAG);
     753  }
     754 
    829755  /* for NcodeSelect, count Nmeas for appropriate photcode */
    830756  if (NcodeSelect) {
    … …  
    937863}
    938864
     865/* extract delta-mag pairs applying specified selections */
    939866double *ExtractDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {
    940867
    941   int i, j, A1, A2, N1, N2, Nlist, NLIST;
    942   double *M1, *M2, *list;
    943 
    944   /* check for special case of measure-measure */
     868  int i, j, A1, A2, N1, N2, Np, Nlist, NLIST;
     869  double *M1, *M2, *list, NoMag;
     870
     871  /* check for special case of measure-measure - this is needed to drop self-matches */
    945872  A1 = ((mode[0] == MAG_AVE) || (mode[0] == MAG_REF));
    946873  A2 = ((mode[1] == MAG_AVE) || (mode[1] == MAG_REF));
    … …  
    955882  ALLOCATE (list, double, NLIST);
    956883  M1 = M2 = NULL;
     884  NoMag = NO_MAG * 0.001;
    957885
    958886  /* one of the two is an average, must do independently */
    … …  
    960888  if (N1 == 0) goto skip;
    961889 
     890  Np = GetSelectionParam ();
     891  SetSelectionParam (Np + 1);
    962892  M2 = ExtractMagnitudes (code[1], mode[1], average, secfilt, measure, &N2);
    963893  if (N2 == 0) goto skip;
    964894
    965   /* average mags can return NO_MAG : skip them? */
     895  /* magnitudes may be NO_MAG : set delta to NO_MAG */
    966896  for (i = 0; i < N1; i++) {
    967897    for (j = 0; j < N2; j++) {
    968       list[Nlist] = M1[i] - M2[j];
     898      if ((M1[i] == NoMag) || (M2[j] == NoMag)) {
     899        list[Nlist] = NoMag;
     900      } else {
     901        list[Nlist] = M1[i] - M2[j];
     902      }
    969903      Nlist ++;
    970904    }
    971905  }
    972   *nlist = Nlist;
    973   return (list);
    974906
    975907 skip:
    976908  if (M1 != NULL) free (M1);
    977909  if (M2 != NULL) free (M2);
    978   M1 = M2 = NULL;
     910  *nlist = Nlist;
    979911  return (list);
    980912}
    981913 
    982 /* extract a list of delta-measure-mags from the specified average entry based on the pre-set selections */
     914/* extract a list of delta-measure-mags from the specified average entry based on the
     915   pre-set selections - does not return self-matched measurements */
    983916double *ExtractMeasuresDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {
    984917
    985   int i, j, Nlist, NLIST;
    986   double *list, M1, M2, value;
     918  int i, j, Np0, Np1, Nlist, NLIST;
     919  double *list, M1, M2, value, NoMag;
    987920 
    988921  *nlist = 0;
    … …  
    990923  NLIST = MAX (1, average[0].Nm*average[0].Nm);
    991924  ALLOCATE (list, double, NLIST);
    992 
    993   /* must have two code values */
     925  NoMag = NO_MAG * 0.001;
     926
     927  /* must have two code values - drop this test? this is programming case, not a user case */
    994928  if (code == NULL) return (list);
    995929  if (code[0] == NULL) return (list);
    996930  if (code[1] == NULL) return (list);
    997931
    998   /* exclude based on chisq for both codes  */
    999   if (ChiSelect) {
    1000     value = PhotXm (code[0], average, secfilt);
    1001     if (value == -1) return (list);
    1002     if (value > ChiLimit) return (list);
    1003     value = PhotXm (code[1], average, secfilt);
    1004     if (value == -1) return (list);
    1005     if (value > ChiLimit) return (list);
    1006   }
     932  /* exclude based on average parameters for both codes  */
     933  if (!TestAverage (code[0], average, secfilt, measure)) return (list);
     934  if (!TestAverage (code[1], average, secfilt, measure)) return (list);
     935
     936  Np0 = GetSelectionParam ();
     937  Np1 = Np0 + 1;
    1007938
    1008939  /* loop twice over all measures */
    1009940  for (i = 0; i < average[0].Nm; i++) {
     941    SetSelectionParam (Np0);
    1010942    TESTCODE (code[0], measure[i]);
    1011943    TESTMEASURE (measure[i]);
    … …  
    1013945    for (j = 0; j < average[0].Nm; j++) {
    1014946      if (i == j) continue;
     947      SetSelectionParam (Np1);
    1015948      TESTCODE (code[1], measure[j]);
    1016949      TESTMEASURE (measure[j]);
    1017950      SETMAG(M2, measure[j], mode[1]);
    1018       list[Nlist] = M1 - M2;
     951      if ((M1 == NoMag) || (M2 == NoMag)) {
     952        list[Nlist] = NoMag;
     953      } else {
     954        list[Nlist] = M1 - M2;
     955      }
    1019956      Nlist ++;
    1020957    }
    … …  
    1024961}
    1025962
     963/* extract a measurement list matching the number of dmag entries */
    1026964double *ExtractByDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {
    1027965
    … …  
    1037975  }
    1038976
     977  /* one of the two entries results in a single element. extract the other */
    1039978  if (A1) {
    1040979    list = ExtractMeasures (code[1], mode[1], average, secfilt, measure, &N1, param);
    … …  
    1049988double *ExtractMeasuresByDMag (PhotCode **code, int *mode, int use_first, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {
    1050989
    1051   int i, j, k, m, n, Nlist, NLIST;
     990  int i, j, k, n, Nlist, NLIST;
    1052991  time_t t;
    1053992  int np, found;
    1054993  double *list, M1, M2, value;
    1055   double ra, dec, x, y;
    1056994 
    1057995  *nlist = 0;
    … …  
    10821020
    10831021      /* assign value */
    1084       value = 0;
    1085       switch (param) {
    1086         case MEAS_MAG: /* magnitudes are already determined above */
    1087           value = use_first ? M1 : M2;
    1088           break;
    1089         case MEAS_RA: /* OK */
    1090           value = average[0].R - measure[n].dR / 360000.0;
    1091           break;
    1092         case MEAS_DEC: /* OK */
    1093           value = average[0].D - measure[n].dD / 360000.0;
    1094           break;
    1095         case MEAS_dMAG: /* OK */
    1096           value = 0.001*measure[n].dM;
    1097           break;
    1098         case MEAS_AIRMASS: /* OK */
    1099           value = 0.001*measure[n].airmass;
    1100           break;
    1101         case MEAS_EXPTIME: /* OK */
    1102           value = pow (10.0, measure[n].dt * 0.0004);
    1103           break;
    1104         case MEAS_PHOTCODE: /* OK */
    1105           value = measure[n].source;
    1106           break;
    1107         case MEAS_TIME: /* OK */
    1108           value = TimeValue (measure[n].t, TimeReference, TimeFormat);
    1109           break;
    1110         case MEAS_dR: /* OK */
    1111           value = 0.01*measure[n].dR;
    1112           break;
    1113         case MEAS_dD: /* OK */
    1114           value = 0.01*measure[n].dD;
    1115           break;
    1116         case MEAS_FWHM: /* OK */
    1117           value = 0.01*measure[n].FWx;
    1118           break;
    1119         case MEAS_DOPHOT: /* OK */
    1120           value = measure[n].dophot;
    1121           break;
    1122         case MEAS_FLAGS: /* ? */
    1123           value = measure[n].flags;
    1124           break;
    1125         case MEAS_XCCD: /* OK */
    1126           ra  = average[0].R - measure[n].dR / 360000.0;
    1127           dec = average[0].D - measure[n].dD / 360000.0;
    1128           m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
    1129           if (m == -1) break;
    1130           if (!FindMosaicForImage (image, Nimage, m)) break;
    1131           RD_to_XY (&x, &y, ra, dec, &image[m].coords);
    1132           value = x;
    1133           break;
    1134         case MEAS_YCCD: /* OK */
    1135           ra  = average[0].R - measure[n].dR / 360000.0;
    1136           dec = average[0].D - measure[n].dD / 360000.0;
    1137           m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
    1138           if (m == -1) break;
    1139           if (!FindMosaicForImage (image, Nimage, m)) break;
    1140           RD_to_XY (&x, &y, ra, dec, &image[m].coords);
    1141           value = y;
    1142           break;
    1143         case MEAS_XMOSAIC: /* OK */
    1144           ra  = average[0].R - measure[n].dR / 360000.0;
    1145           dec = average[0].D - measure[n].dD / 360000.0;
    1146           m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
    1147           if (m == -1) break;
    1148           mosaic.crval1 = image[m].coords.crval1;
    1149           mosaic.crval2 = image[m].coords.crval2;
    1150           RD_to_XY (&x, &y, ra, dec, &mosaic);
    1151           value = x;
    1152           break;
    1153         case MEAS_YMOSAIC: /* OK */
    1154           ra  = average[0].R - measure[n].dR / 360000.0;
    1155           dec = average[0].D - measure[n].dD / 360000.0;
    1156           m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
    1157           if (m == -1) break;
    1158           mosaic.crval1 = image[m].coords.crval1;
    1159           mosaic.crval2 = image[m].coords.crval2;
    1160           RD_to_XY (&x, &y, ra, dec, &mosaic);
    1161           value = y;
    1162           break;
    1163       }
    1164       list[Nlist] = value;
     1022      list[Nlist] = GetMeasure (param, &average[0], &measure[n], (use_first ? M1 : M2));
    11651023      Nlist ++;
    11661024    }
    … …  
    11701028}
    11711029
     1030double GetMeasure (int param, Average *average, Measure *measure, double mag) {
     1031
     1032  int m;
     1033  double ra, dec, x, y;
     1034  double value;
     1035  Image *image;
     1036  Coords *mosaic;
     1037
     1038  value = 0;
     1039  switch (param) {
     1040    case MEAS_MAG: /* magnitudes are already determined above */
     1041      value = mag;
     1042      break;
     1043    case MEAS_RA: /* OK */
     1044      value = average[0].R - measure[0].dR / 360000.0;
     1045      break;
     1046    case MEAS_DEC: /* OK */
     1047      value = average[0].D - measure[0].dD / 360000.0;
     1048      break;
     1049    case MEAS_dMAG: /* OK */
     1050      value = 0.001*measure[0].dM;
     1051      break;
     1052    case MEAS_AIRMASS: /* OK */
     1053      value = 0.001*measure[0].airmass;
     1054      break;
     1055    case MEAS_EXPTIME: /* OK */
     1056      value = pow (10.0, measure[0].dt * 0.0004);
     1057      break;
     1058    case MEAS_PHOTCODE: /* OK */
     1059      value = measure[0].source;
     1060      break;
     1061    case MEAS_TIME: /* OK */
     1062      value = TimeValue (measure[0].t, TimeReference, TimeFormat);
     1063      break;
     1064    case MEAS_dR: /* OK */
     1065      value = 0.01*measure[0].dR;
     1066      break;
     1067    case MEAS_dD: /* OK */
     1068      value = 0.01*measure[0].dD;
     1069      break;
     1070    case MEAS_FWHM: /* OK */
     1071      value = 0.01*measure[0].FWx;
     1072      break;
     1073    case MEAS_DOPHOT: /* OK */
     1074      value = measure[0].dophot;
     1075      break;
     1076    case MEAS_FLAGS: /* ? */
     1077      value = measure[0].flags;
     1078      break;
     1079    case MEAS_XCCD: /* OK */
     1080      ra  = average[0].R - measure[0].dR / 360000.0;
     1081      dec = average[0].D - measure[0].dD / 360000.0;
     1082      image = MatchImage (measure[0].t, measure[0].source);
     1083      if (image == NULL) break;
     1084      RD_to_XY (&x, &y, ra, dec, &image[0].coords);
     1085      value = x;
     1086      break;
     1087    case MEAS_YCCD: /* OK */
     1088      ra  = average[0].R - measure[0].dR / 360000.0;
     1089      dec = average[0].D - measure[0].dD / 360000.0;
     1090      image = MatchImage (measure[0].t, measure[0].source);
     1091      if (image == NULL) break;
     1092      RD_to_XY (&x, &y, ra, dec, &image[0].coords);
     1093      value = y;
     1094      break;
     1095    case MEAS_XMOSAIC: /* OK */
     1096      ra  = average[0].R - measure[0].dR / 360000.0;
     1097      dec = average[0].D - measure[0].dD / 360000.0;
     1098      mosaic = MatchMosaic (measure[0].t, measure[0].source);
     1099      if (mosaic == NULL) break;
     1100      RD_to_XY (&x, &y, ra, dec, mosaic);
     1101      value = x;
     1102      break;
     1103    case MEAS_YMOSAIC: /* OK */
     1104      ra  = average[0].R - measure[0].dR / 360000.0;
     1105      dec = average[0].D - measure[0].dD / 360000.0;
     1106      mosaic = MatchMosaic (measure[0].t, measure[0].source);
     1107      if (mosaic == NULL) break;
     1108      RD_to_XY (&x, &y, ra, dec, mosaic);
     1109      value = y;
     1110      break;
     1111  }
     1112  return (value);
     1113}
     1114
     1115/** the mosaic entries do not use the registered mosaic found
     1116    by MatchImage (via FindMosaicForImage).  Rather, they use
     1117    a coordinate frame saved by SetImageSelection
     1118**/
Note: See TracChangeset for help on using the changeset viewer.