IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 23, 2007, 4:39:53 PM (19 years ago)
Author:
eugene
Message:

remove the concept of primary / secondary photcodes: all average mags go into the secfilt table

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/dvo-mods-2007-02/Ohana/src/libdvo/src/LoadPhotcodes.c

    r12011 r12031  
    11# include <dvo.h>
    22
    3 # define NCTERMS 4
    4 # define F 0.001
    5 /* # define NO_MAG 100.0 */
     3int LoadPhotcodes (char *catdir_file, char *master_file) {
    64
    7 static PhotCodeData *photcodes = NULL;
    8 static double ZERO_POINT;
    9 static short int Nseclist[0x10000];
     5  /* first try to load the photcodes from the specified CATDIR location */
     6  if (LoadPhotcodesFITS (catdir_file)) return TRUE;
     7 
     8  /* next try to load the photcodes from the master text photcode file */
     9  if (LoadPhotcodesFITS (master_file)) return TRUE;
    1010
    11 /* static short int iZERO_POINT; */
    12 
    13 void SetZeroPoint (double ZP) {
    14   ZERO_POINT = ZP;
    15   /* iZERO_POINT = 1000 * ZP; */
     11  return FALSE;
    1612}
    17 
    18 int LoadPhotcodes (char *filename) {
    19  
    20   FILE *f;
    21   int i, Ns, Np, NPHOTCODE, Npri, Nfield;
    22   PhotCode *photcode;
    23   int code;
    24   char *c;
    25   char line[256];
    26   char name[32], type[32], Zero[32], Airmass[32], Offset[32],
    27     C1[32], C2[32], Slope[32], Color[32], Primary[32];
    28   int c1, c2;
    29 
    30   /* allocate space to photcode table, free existing data */
    31   if (photcodes == NULL) {
    32     ALLOCATE (photcodes, PhotCodeData, 1);
    33     photcodes[0].code = NULL;
    34   }
    35   if (photcodes[0].code != NULL) free (photcodes[0].code);
    36 
    37   f = fopen (filename, "r");
    38   if (f == (FILE *) NULL) {
    39     photcodes[0].Ncode    = 0;
    40     photcodes[0].Nsecfilt = 0;
    41     photcodes[0].code     = (PhotCode *) NULL;
    42     return (FALSE);
    43   }
    44 
    45   Np = 0;
    46   NPHOTCODE = 10;
    47   ALLOCATE (photcode, PhotCode, NPHOTCODE);
    48   photcodes[0].Nsecfilt = 0;
    49   for (i = 0; i < 0x10000; i++) {
    50     photcodes[0].hashcode[i] = -1;
    51     photcodes[0].hashNsec[i] = -1;
    52   }
    53 
    54   while (scan_line (f, line) != EOF) {
    55     for (c = line; isspace (*c); c++);
    56     if (*c == '#') continue;
    57     Nfield = sscanf (c, "%d %s %s %s %s %s %s %s %s %s %s",
    58                      &code, name, type, Zero, Airmass, Offset, C1, C2, Slope, Color, Primary);
    59     if (Nfield != 11) { continue; }
    60    
    61     c1 = atof (C1);
    62     c2 = atof (C2);
    63     if (!strcmp (C1, "-")) { c1 = 0; }
    64     if (!strcmp (C2, "-")) { c2 = 0; }
    65 
    66     photcode[Np].type = 0;
    67     photcode[Np].code = code;
    68     strcpy (photcode[Np].name, name);
    69     if (!strncasecmp (type, "pri", 3)) {
    70       photcode[Np].type  = PHOT_PRI;
    71       photcode[Np].C     = 1000*atof (Zero);
    72       photcode[Np].K     = atof (Airmass);
    73       photcode[Np].dC    = 1000*atof (Offset);
    74       photcode[Np].dX    = 1000*atof (Color);
    75       photcode[Np].c1    = c1;
    76       photcode[Np].c2    = c2;
    77       photcode[Np].equiv = atoi (Primary);
    78       ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    79       Nseclist[0] = Np;
    80     }     
    81     if (!strncasecmp (type, "sec", 3)) {
    82       photcode[Np].type  = PHOT_SEC;
    83       photcode[Np].C     = 1000*atof (Zero);
    84       photcode[Np].K     = atof (Airmass);
    85       photcode[Np].dC    = 1000*atof (Offset);
    86       photcode[Np].dX    = 1000*atof (Color);
    87       photcode[Np].c1    = c1;
    88       photcode[Np].c2    = c2;
    89       photcode[Np].equiv = atoi (Primary);
    90       photcodes[0].Nsecfilt ++;
    91       ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    92       Nseclist[photcodes[0].Nsecfilt] = Np;
    93     }     
    94     if (!strncasecmp (type, "dep", 3)) {
    95       photcode[Np].type  = PHOT_DEP;
    96       photcode[Np].C     = 1000*atof (Zero);    /* zero point in millimags */
    97       photcode[Np].K     = atof (Airmass);      /* airmass coeff (millimag / millimag) */
    98       photcode[Np].dC    = 1000*atof (Offset);  /* color ref z.p. (millimag) */
    99       photcode[Np].dX    = 1000*atof (Color);   /* average color (millimag) */
    100       photcode[Np].c1    = c1;
    101       photcode[Np].c2    = c2;
    102       photcode[Np].equiv = atoi (Primary);
    103       ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    104     }     
    105     if (!strncasecmp (type, "ref", 3)) {
    106       photcode[Np].type  = PHOT_REF;
    107       photcode[Np].C     = 0;
    108       photcode[Np].K     = 0;
    109       photcode[Np].dC    = 0;
    110       photcode[Np].dX    = 0;
    111       photcode[Np].c1    = 0;
    112       photcode[Np].c2    = 0;
    113       photcode[Np].equiv = atoi (Primary);
    114       photcode[Np].X[0]  = 0;
    115       photcode[Np].Nc    = 0;
    116     }
    117 
    118     /* alt photcodes are a little different: they have the SAME photcode as an existing
    119        pri/sec photcode, but define an alternate transformation for that code */
    120     if (!strncasecmp (type, "alt", 3)) {
    121       photcode[Np].type  = PHOT_ALT;
    122       photcode[Np].C     = 1000*atof (Zero);    /* zero point in millimags */
    123       photcode[Np].K     = atof (Airmass);      /* airmass coeff (millimag / millimag) */
    124       photcode[Np].dC    = 1000*atof (Offset);  /* color ref z.p. (millimag) */
    125       photcode[Np].dX    = 1000*atof (Color);   /* average color (millimag) */
    126       photcode[Np].c1    = c1;
    127       photcode[Np].c2    = c2;
    128       photcode[Np].equiv = atoi (Primary);
    129       ParseColorTerms (Slope, photcode[Np].X, &photcode[Np].Nc);
    130     }
    131     if (!photcode[Np].type) {
    132       fprintf (stderr, "error in Photfile: unknown type %s\n", type);
    133     }
    134 
    135     Np++;
    136     if (Np == NPHOTCODE) {
    137       NPHOTCODE += 10;
    138       REALLOCATE (photcode, PhotCode, NPHOTCODE);
    139     }
    140   }
    141   fclose (f);
    142  
    143   /* set up hashcode for photcode refs:
    144    * the hashcode gives the structure sequence for a given photcode:
    145    * photcode[i].hashcode[photcode[i].code] == i
    146    */
    147 
    148   Ns = 0;
    149   for (i = 0; i < Np; i++) {
    150     if (photcode[i].type == PHOT_ALT) continue; /* no hashcode for ALT codes */
    151     if (photcodes[0].hashcode[photcode[i].code] != -1) {
    152       fprintf (stderr, "duplicate photcodes in file\n");
    153       code = photcodes[0].hashcode[photcode[i].code];
    154       fprintf (stderr, "conflict between %s (%d) and %s (%d)\n",
    155                photcode[i].name, photcode[i].code, photcode[code].name, photcode[code].code);
    156       free (photcode);
    157       return (FALSE);
    158     }
    159     photcodes[0].hashcode[photcode[i].code] = i;
    160     if (photcode[i].type == PHOT_SEC) {
    161       photcodes[0].hashNsec[photcode[i].code] = Ns;
    162       Ns ++;
    163     }
    164   }
    165   /* validity check for references */
    166   for (i = 0; i < Np; i++) {
    167     if (photcode[i].type == PHOT_DEP) {
    168       Npri = photcodes[0].hashcode[photcode[i].equiv];
    169       if ((Npri >= Np) || (Npri < 0)) {
    170         fprintf (stderr, "reference for dependent photcode is not in photcodes\n");
    171         free (photcode);
    172         return (FALSE);
    173       }
    174       if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) {
    175         fprintf (stderr, "reference for dependent photcode is not a primary or secondary code\n");
    176         free (photcode);
    177         return (FALSE);
    178       }
    179     }
    180     if (photcode[i].type == PHOT_ALT) {
    181       Npri = photcodes[0].hashcode[photcode[i].code];
    182       if ((Npri >= Np) || (Npri < 0)) {
    183         fprintf (stderr, "reference for alternate photcode is not in photcodes\n");
    184         free (photcode);
    185         return (FALSE);
    186       }
    187       if ((photcode[Npri].type != PHOT_PRI) && (photcode[Npri].type != PHOT_SEC)) {
    188         fprintf (stderr, "reference for alternate photcode is not a primary or secondary code\n");
    189         free (photcode);
    190         return (FALSE);
    191       }
    192     }
    193   }
    194   photcodes[0].code = photcode;
    195   photcodes[0].Ncode = Np;
    196 
    197   return (TRUE);
    198 
    199 }
    200 
    201 void ParseColorTerms (char *terms, float *X, int *N) {
    202 
    203   int i;
    204   char *p;
    205 
    206   p = terms;
    207 
    208   for (i = 0; (p != NULL) && (i < NCTERMS); i++) {
    209     X[i] = atof (p);
    210     p = strchr (p, ',');
    211     if (p == (char *) NULL) continue;
    212     p ++;
    213   }
    214   *N = i;
    215 }
    216 
    217 /********** photcode lookups **********/
    218 
    219 /* return photcode for given name */
    220 PhotCode *GetPhotcodebyName (char *name) {
    221  
    222   int i;
    223 
    224   if (name == (char *) NULL ) return (NULL);
    225  
    226   for (i = 0; i < photcodes[0].Ncode; i++) {
    227     if (!strcmp (photcodes[0].code[i].name, name)) {
    228       return (&photcodes[0].code[i]);
    229     }
    230   }
    231   return (NULL);
    232 }
    233 /* return photcode.code for given name */
    234 int GetPhotcodeCodebyName (char *name) {
    235  
    236   int i;
    237  
    238   if (name == (char *) NULL ) return (0);
    239 
    240   for (i = 0; i < photcodes[0].Ncode; i++) {
    241     if (!strcmp (photcodes[0].code[i].name, name)) {
    242       return (photcodes[0].code[i].code);
    243     }
    244   }
    245   return (0);
    246 }
    247 /* return equivalent photcode for given name */
    248 PhotCode *GetPhotcodeEquivbyName (char *name) {
    249  
    250   int i, equiv;
    251  
    252   if (name == (char *) NULL ) return (NULL);
    253 
    254   for (i = 0; i < photcodes[0].Ncode; i++) {
    255     if (!strcmp (photcodes[0].code[i].name, name)) {
    256       if (photcodes[0].code[i].equiv == 0) return (NULL);
    257       equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv];
    258       if (equiv == -1) return (NULL);
    259       return (&photcodes[0].code[equiv]);
    260     }
    261   }
    262   return (NULL);
    263 }
    264 /* return equivalent photcode.code for given name */
    265 int GetPhotcodeEquivCodebyName (char *name) {
    266  
    267   int i, equiv;
    268  
    269   if (name == (char *) NULL ) return (0);
    270 
    271   for (i = 0; i < photcodes[0].Ncode; i++) {
    272     if (!strcmp (photcodes[0].code[i].name, name)) {
    273       if (photcodes[0].code[i].equiv == 0) return (0);
    274       equiv = photcodes[0].hashcode[photcodes[0].code[i].equiv];
    275       if (equiv == -1) return (0);
    276       return (photcodes[0].code[equiv].code);
    277     }
    278   }
    279   return (0);
    280 }
    281 
    282 /* return photcode for given code */
    283 PhotCode *GetPhotcodebyCode (int code) {
    284  
    285   int entry;
    286  
    287   if (code < 0) return (NULL);
    288   if (code > 0x10000) return (NULL);
    289 
    290   entry = photcodes[0].hashcode[code];
    291   if (entry == -1) return (NULL);
    292 
    293   return (&photcodes[0].code[entry]);
    294 }
    295 /* return photcode.code for given code */
    296 char *GetPhotcodeNamebyCode (int code) {
    297  
    298   int entry;
    299  
    300   if (code < 0) return (NULL);
    301   if (code > 0x10000) return (NULL);
    302 
    303   entry = photcodes[0].hashcode[code];
    304   if (entry == -1) return (NULL);
    305 
    306   return (photcodes[0].code[entry].name);
    307 }
    308 /* return equivalent photcode for given code */
    309 PhotCode *GetPhotcodeEquivbyCode (int code) {
    310  
    311   int entry, equiv;
    312  
    313   if (code < 0) return (NULL);
    314   if (code > 0x10000) return (NULL);
    315 
    316   entry = photcodes[0].hashcode[code];
    317   if (entry == -1) return (NULL);
    318 
    319   if (photcodes[0].code[entry].equiv == 0) return (NULL);
    320   equiv = photcodes[0].hashcode[photcodes[0].code[entry].equiv];
    321 
    322   if (equiv == -1) return (NULL);
    323   return (&photcodes[0].code[equiv]);
    324 }
    325 /* return equivalent photcode.code for given code */
    326 int GetPhotcodeEquivCodebyCode (int code) {
    327  
    328   int entry;
    329  
    330   if (code < 0) return (0);
    331   if (code > 0x10000) return (0);
    332 
    333   entry = photcodes[0].hashcode[code];
    334   if (entry == -1) return (0);
    335   return (photcodes[0].code[entry].equiv);
    336 }
    337 
    338 int GetPhotcodeNsec (int code) {
    339  
    340   int Nsec;
    341  
    342   if (code < 0) return (-1);
    343   if (code > 0x10000) return (-1);
    344 
    345   Nsec = photcodes[0].hashNsec[code];
    346   return (Nsec);
    347 }
    348 
    349 /* Nsec of 0 is PRI */
    350 PhotCode *GetPhotcodebyNsec (int Nsec) {
    351  
    352   if (Nsec > photcodes[0].Nsecfilt) return (NULL);
    353   if (Nsec < 0) return (NULL);
    354  
    355   return (&photcodes[0].code[Nseclist[Nsec]]);
    356 }
    357 
    358 int GetPhotcodeNsecfilt () {
    359   return (photcodes[0].Nsecfilt);
    360 }
    361 
    362 /* ALLOCATE and return list of all photcodes
    363    with photcode.equiv == code */
    364 int *GetPhotcodeEquivList (int code, int *nlist) {
    365 
    366   int i, Nlist;
    367   int *list;
    368 
    369   ALLOCATE (list, int, MAX (1, photcodes[0].Ncode));
    370   Nlist = 0;
    371   for (i = 0; i < photcodes[0].Ncode; i++) {
    372     if (photcodes[0].code[i].equiv != code) continue;
    373     list[Nlist] = photcodes[0].code[i].code;
    374     Nlist ++;
    375   }
    376   REALLOCATE (list, int, MAX (1, Nlist));
    377 
    378   *nlist = Nlist;
    379   return (list);
    380 }
    381 
    382 /******** internal photometry conversions (keeps values in short int millimags) *********/
    383 float PhotInst (Measure *measure) {
    384 
    385   int Np;
    386   float M;
    387 
    388   Np = photcodes[0].hashcode[measure[0].photcode];
    389   if (Np == -1) return (NO_MAG);
    390 
    391   if (photcodes[0].code[Np].type == PHOT_REF) {
    392     M = measure[0].M;
    393     return (M);
    394   }
    395 
    396   M = measure[0].M - measure[0].dt - ZERO_POINT;
    397          
    398   return (M);
    399 
    400 }
    401 
    402 float PhotCat (Measure *measure) {
    403 
    404   int Np;
    405   float Mcat;
    406   PhotCode *code;
    407 
    408   Np = photcodes[0].hashcode[measure[0].photcode];
    409   if (Np == -1) return (NO_MAG);
    410 
    411   if (photcodes[0].code[Np].type == PHOT_REF) {
    412     Mcat = measure[0].M;
    413     return (Mcat);
    414   }
    415   code = &photcodes[0].code[Np];
    416   Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + F*code[0].C;
    417  
    418   return (Mcat);
    419 }
    420 
    421 float PhotSys (Measure *measure, Average *average, SecFilt *secfilt) {
    422 
    423   int i, Np;
    424   float Mcat, Mcol, Msys, mc, Mc;
    425   PhotCode *code;
    426 
    427   Np = photcodes[0].hashcode[measure[0].photcode];
    428   if (Np == -1) return (NO_MAG);
    429 
    430   if (photcodes[0].code[Np].type == PHOT_REF) {
    431     Msys = measure[0].M;
    432     return (Msys);
    433   }
    434   code = &photcodes[0].code[Np];
    435   Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + F*code[0].C;
    436 
    437   /* for DEP, color must be made of PRI/SEC */
    438   mc = PhotColorForCode (average, secfilt, NULL, code);
    439   if (mc == NO_MAG) return (Mcat);
    440   mc = mc - F*code[0].dX;
    441 
    442   Mc = mc;
    443   Mcol = 0;
    444   for (i = 0; i < code[0].Nc; i++) {
    445     Mcol += code[0].X[i]*Mc;
    446     Mc *= mc;
    447   }
    448   Msys = Mcat + Mcol;
    449   return (Msys);
    450 }
    451 
    452 float PhotRel (Measure *measure, Average *average, SecFilt *secfilt) {
    453 
    454   int i, Np;
    455   float Mcat, Mcol, Mrel, mc, Mc;
    456   PhotCode *code;
    457 
    458   Np = photcodes[0].hashcode[measure[0].photcode];
    459   if (Np == -1) return (NO_MAG);
    460 
    461   if (photcodes[0].code[Np].type == PHOT_REF) {
    462     Mcat = measure[0].M;
    463     return (Mcat);
    464   }
    465   code = &photcodes[0].code[Np];
    466   Mrel = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + F*code[0].C - measure[0].Mcal;
    467 
    468   /* for DEP, color must be made of PRI/SEC */
    469   mc = PhotColorForCode (average, secfilt, NULL, code);
    470   if (mc == NO_MAG) return (Mrel);
    471   mc = mc - F*code[0].dX;
    472 
    473   Mc = mc;
    474   Mcol = 0;
    475   for (i = 0; i < code[0].Nc; i++) {
    476     Mcol += code[0].X[i]*Mc;
    477     Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
    478   }
    479   Mrel += Mcol;
    480   return (Mrel);
    481 }
    482 
    483 /* return calibrated magnitude from measure for given photcode */
    484 float PhotCal (Measure *thisone, Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {
    485 
    486   int i, Np;
    487   float Mcal, Mrel, Mcol, mc, Mc;
    488 
    489   /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */
    490   Np = photcodes[0].hashcode[thisone[0].photcode];
    491   if (Np == -1) {
    492     return (NO_MAG);
    493   }
    494 
    495   if (photcodes[0].code[Np].type == PHOT_REF) {
    496     Mrel = thisone[0].M;
    497     return (Mrel);
    498   }
    499   if (code[0].code != photcodes[0].code[Np].equiv) {
    500     return (NO_MAG);
    501   }
    502 
    503   Mcal = PhotRel (thisone, average, secfilt) + F*code[0].C;
    504 
    505   mc = PhotColorForCode (average, secfilt, measure, code);
    506   if (mc == NO_MAG) return (Mcal);
    507   mc = mc - F*code[0].dX;
    508 
    509   Mc = mc;
    510   Mcol = 0;
    511   for (i = 0; i < code[0].Nc; i++) {
    512     Mcol += code[0].X[i]*Mc;
    513     Mc *= mc;
    514   }
    515   Mcal += Mcol;
    516   return (Mcal);
    517 }
    518 
    519 /* color term may not use DEP magnitude */
    520 float PhotColorForCode (Average *average, SecFilt *secfilt, Measure *measure, PhotCode *code) {
    521 
    522   int i, Ns1, Ns2, Ns;
    523   float m1, m2, mc;
    524   PhotCode *color;
    525 
    526   m1 = m2 = NO_MAG;
    527 
    528   if (measure == NULL) {
    529     Ns1 = photcodes[0].hashNsec[code[0].c1];
    530     Ns2 = photcodes[0].hashNsec[code[0].c2];
    531  
    532     m1 = (Ns1 == -1) ? NO_MAG : secfilt[Ns1].M;
    533     m2 = (Ns2 == -1) ? NO_MAG : secfilt[Ns2].M;
    534     mc = ((m1 == NO_MAG) || (m2 == NO_MAG)) ? NO_MAG : (m1 - m2);
    535     return (mc);
    536   }
    537 
    538   /* find magnitude matching first color term */
    539   color = GetPhotcodebyCode (code[0].c1);
    540   if (color == NULL) return (NO_MAG);
    541   if (color[0].type == PHOT_REF) {
    542     for (i = 0; (i < average[0].Nm) && (m1 == NO_MAG); i++) {
    543       if (measure[i].photcode == color[0].code) {
    544         m1 = measure[i].M;
    545       }
    546     }
    547   } else {
    548     Ns = photcodes[0].hashNsec[color[0].code];
    549     m1 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    550   }     
    551 
    552   /* find magnitude matching second color term */
    553   color = GetPhotcodebyCode (code[0].c2);
    554   if (color == NULL) return (NO_MAG);
    555   if (color[0].type == PHOT_REF) {
    556     for (i = 0; (i < average[0].Nm) && (m2 == NO_MAG); i++) {
    557       if (measure[i].photcode == color[0].code) {
    558         m2 = measure[i].M;
    559       }
    560     }
    561   } else {
    562     Ns = photcodes[0].hashNsec[color[0].code];
    563     m2 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    564   }     
    565   mc = ((m1 == NO_MAG) || (m2 == NO_MAG)) ? NO_MAG : (m1 - m2);
    566   return (mc);
    567 }
    568 
    569 /* return calibrated magnitude from average/secfilt for given photcode */
    570 float PhotRef (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {
    571 
    572   int i, Ns;
    573   float Mave, Mref, Mcol, mc;
    574   double Mc;
    575 
    576   Ns = photcodes[0].hashNsec[code[0].code];
    577   Mave = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    578   Mref = Mave + F*code[0].C;
    579 
    580   mc = PhotColorForCode (average, secfilt, measure, code);
    581   if (mc == NO_MAG) return (Mref);
    582   mc = mc - F*code[0].dX;
    583 
    584   Mc = mc;
    585   Mcol = 0;
    586   for (i = 0; i < code[0].Nc; i++) {
    587     Mcol += code[0].X[i]*Mc;
    588     Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
    589   }
    590   Mref += Mcol;
    591   return (Mref);
    592 }
    593 
    594 /***/
    595 float PhotAve (PhotCode *code, Average *average, SecFilt *secfilt) {
    596 
    597   int Ns;
    598   float Mave;
    599 
    600   Ns = photcodes[0].hashNsec[code[0].code];
    601   Mave = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    602   return (Mave);
    603 }
    604 
    605 float PhotdM (PhotCode *code, Average *average, SecFilt *secfilt) {
    606 
    607   int Ns;
    608   float dM;
    609 
    610   Ns = photcodes[0].hashNsec[code[0].code];
    611   dM  = (Ns == -1) ? NO_MAG : secfilt[Ns].dM;
    612   return (dM);
    613 }
    614 
    615 float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt) {
    616 
    617   int Ns;
    618   short Mi;
    619   float Xm;
    620 
    621   Ns = photcodes[0].hashNsec[code[0].code];
    622   Mi = (Ns == -1) ? NO_MAG : secfilt[Ns].Xm;
    623   Xm = (Mi == NO_MAG) ? -1.0 : pow (10.0, 0.01*Mi);
    624   return (Xm);
    625 }
    626 
    627 /* given a photcode pair c1 & c2, return the color of this star (NaN if not found) */
    628 int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color) {
    629 
    630   int i, Ns;
    631   double M1, M2, dM;
    632   PhotCode *code;
    633  
    634   code = GetPhotcodebyCode (c1);
    635   if (code == NULL) return (FALSE);
    636   if (code[0].type == PHOT_REF) {
    637     for (i = 0; i < average[0].Nm; i++) {
    638       if (measure[i].photcode == c1) {
    639         M1 = measure[i].M;
    640         goto filter1;
    641       }
    642     }   
    643     return (FALSE);
    644   } else {
    645     Ns = photcodes[0].hashNsec[code[0].code];
    646     M1 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    647   }     
    648 
    649 filter1:
    650   code = GetPhotcodebyCode (c2);
    651   if (code == NULL) return (FALSE);
    652   if (code[0].type == PHOT_REF) {
    653     for (i = 0; i < average[0].Nm; i++) {
    654       if (measure[i].photcode == c2) {
    655         M2 = measure[i].M;
    656         goto filter2;
    657       }
    658     }   
    659     return (FALSE);
    660   } else {
    661     Ns = photcodes[0].hashNsec[code[0].code];
    662     M2 = (Ns == -1) ? NO_MAG : secfilt[Ns].M;
    663   }     
    664  
    665 filter2:
    666 
    667   dM = M1 - M2;
    668   *color = dM;
    669  
    670   return (TRUE);
    671 }
    672 
    673 /* photcode table should have the following format:
    674 
    675 # code name     type  zero  airmass  offset  c1 c2  slope  <color>  primary
    676 1    B        pri   24.0  0.15     -       -  -   -      -        -
    677 2    B        pri   24.0  0.15     -       -  -   -      -        -
    678 3    B1       sec   22.5  0.18     0.15    1  2   0.10   0.50     1
    679 1000 USNO_B   ref   -     -        -       -  -   -      -        -
    680 
    681 */
    682 
    683 
    684 /*
    685   Nc1 = photcodes[0].code[Np].c1;
    686   Ns1 = photcodes[0].hashNsec[Nc1];
    687 
    688   Nc2 = photcodes[0].code[Np].c2;
    689   Ns2 = photcodes[0].hashNsec[Nc2];
    690 
    691   Xlam = photcodes[0].code[Np].X[0];
    692   Klam = photcodes[0].code[Np].K;
    693  
    694   m1 = (Ns1 == -1) ? average[0].M : secfilt[Ns1].M;
    695   m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M;
    696 */
Note: See TracChangeset for help on using the changeset viewer.