IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33654


Ignore:
Timestamp:
Apr 1, 2012, 3:05:06 PM (14 years ago)
Author:
eugene
Message:

add support for ubercal in setphot; add support for parallel dvo in setphot

Location:
trunk/Ohana/src/uniphot
Files:
2 deleted
13 edited
11 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/uniphot/Makefile

    r31160 r33654  
    1 default: uniphot setphot setfwhm
     1default: uniphot setphot setphot_client setfwhm
    22help:
    33@echo "make options: uniphot setphot setfwhm default help install default (uniphot setphot setfwhm)"
    … …  
    1818
    1919uniphot: $(BIN)/uniphot.$(ARCH)
     20setfwhm: $(BIN)/setfwhm.$(ARCH)
    2021setphot: $(BIN)/setphot.$(ARCH)
    21 setfwhm: $(BIN)/setfwhm.$(ARCH)
    22 install: $(DESTBIN)/uniphot $(DESTBIN)/setphot $(DESTBIN)/setfwhm
     22setphot_client: $(BIN)/setphot_client.$(ARCH)
     23
     24install: $(DESTBIN)/uniphot $(DESTBIN)/setfwhm $(DESTBIN)/setphot $(DESTBIN)/setphot_client
    2325
    2426UNIPHOT = \
    … …  
    3335$(SRC)/find_image_sgroups.$(ARCH).o \
    3436$(SRC)/fit_groups.$(ARCH).o         \
    35 $(SRC)/update.$(ARCH).o            \
    36 $(SRC)/update_catalog.$(ARCH).o    \
     37$(SRC)/update_dvo_uniphot.$(ARCH).o \
     38$(SRC)/update_catalog_uniphot.$(ARCH).o \
    3739$(SRC)/convert.$(ARCH).o            \
    3840$(SRC)/SetSignals.$(ARCH).o         \
    … …  
    4547SETPHOT = \
    4648$(SRC)/setphot.$(ARCH).o            \
    47 $(SRC)/initialize.$(ARCH).o         \
    48 $(SRC)/ConfigInit.$(ARCH).o         \
    49 $(SRC)/args.$(ARCH).o               \
    50 $(SRC)/liststats.$(ARCH).o          \
     49$(SRC)/initialize_setphot.$(ARCH).o \
    5150$(SRC)/load_zpt_table.$(ARCH).o     \
    52 $(SRC)/load_images.$(ARCH).o        \
     51$(SRC)/load_images_setphot.$(ARCH).o \
     52$(SRC)/ImageSubset.$(ARCH).o        \
    5353$(SRC)/match_zpts_to_images.$(ARCH).o       \
     54$(SRC)/match_flatcorr_to_images.$(ARCH).o           \
     55$(SRC)/update_dvo_setphot.$(ARCH).o \
    5456$(SRC)/update_catalog_setphot.$(ARCH).o \
    5557$(SRC)/SetSignals.$(ARCH).o         \
    5658$(SRC)/Shutdown.$(ARCH).o           
    5759
    58 $(SETPHOT): $(INC)/uniphot.h
     60$(SETPHOT): $(INC)/setphot.h
    5961$(BIN)/setphot.$(ARCH): $(SETPHOT)
     62
     63SETPHOT_CLIENT = \
     64$(SRC)/setphot_client.$(ARCH).o     \
     65$(SRC)/update_dvo_setphot.$(ARCH).o \
     66$(SRC)/update_catalog_setphot.$(ARCH).o \
     67$(SRC)/initialize_setphot_client.$(ARCH).o \
     68$(SRC)/ImageSubset.$(ARCH).o
     69
     70$(SETPHOT_CLIENT): $(INC)/setphot.h
     71$(BIN)/setphot_client.$(ARCH): $(SETPHOT_CLIENT)
    6072
    6173SETFWHM =                           \

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

    r31635 r33654  
    6666char         STATMODE[64];
    6767int          VERBOSE;
     68int          UBERCAL; // load the supplied ubercal zero point fits table (with flat-field corrections)
     69int          NO_METADATA; // the supplied ubercal data has no descriptive metadata
    6870int          NLOOP;
    6971int          TimeSelect;
    … …  
    9597void          XDead               PROTO((void));
    9698int           args_uniphot        PROTO((int argc, char **argv));
    97 int           args_setphot        PROTO((int argc, char **argv));
    9899void          dumpresult          PROTO((void));
    99100Group        *find_image_sgroups  PROTO((FITS_DB *db, ImageLink **imlink, int *Nsgroup));
    … …  
    103104int           gcatalog            PROTO((Catalog *catalog));
    104105void          initialize_uniphot  PROTO((int argc, char **argv));
    105 void          initialize_setphot  PROTO((int argc, char **argv));
    106106void          initstats           PROTO((char *mode));
    107107int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
    108108int           load_images_uniphot PROTO((FITS_DB *db));
     109Image        *load_images_setfwhm PROTO((FITS_DB *db, off_t *Nimage));
    109110int           main                PROTO((int argc, char **argv));
    110111int           open_graph          PROTO((int N));
    … …  
    112113void          sortB               PROTO((double *X, double *Y, int N));
    113114void          sortD               PROTO((double *X, double *Y, double *Z, int N));
    114 void          update_dvo_catalog  PROTO((Catalog *catalog, Group *sgroup, int warn));
    115115void          wcatalog            PROTO((Catalog *catalog));
    116116void          wimages             PROTO((Image *image, int Nimage));
    … …  
    125125int           SetSignals          PROTO((void));
    126126int           subset_images       PROTO((FITS_DB *db));
    127 void          update              PROTO((FITS_DB *db, Group *sgroup, int Nsgroup));
     127void          update_dvo_uniphot  PROTO((FITS_DB *db, Group *sgroup, int Nsgroup));
     128void          update_catalog_uniphot PROTO((Catalog *catalog, Group *sgroup, int warn));
    128129void          sort_time           PROTO((unsigned int *value, int N));
    129130
    130 // setphot-specific prototypes
    131 ZptTable     *load_zpt_table         PROTO((char *filename, int *nzpts));
    132 Image        *load_images_setphot    PROTO((FITS_DB *db, off_t *Nimage));
    133 int           match_zpts_to_images   PROTO((Image *image, off_t Nimage, ZptTable *zpts, int Nzpts));
    134 int           update_setphot         PROTO((Image *image, off_t Nimage));
    135 void          update_catalog_setphot PROTO((Catalog *catalog, Image *image, off_t *index, off_t Nimage));
    136 
    137 /*** time/coord conversion functions not supplied by libohana ***/
    138 time_t        TimeRef               PROTO((double time, time_t TimeReference, int TimeFormat));
    139 double        TimeValue             PROTO((time_t time, time_t TimeReference, int TimeFormat));
    140 
    141 int           hh_hms                PROTO((double hh, int *hr, int *mn, double *sc));
    142 int           dd_dms                PROTO((double dd, int *dg, int *mn, double *sc));
    143 int           hms_format            PROTO((char *line, double value));
    144 int           dms_format            PROTO((char *line, double value));
    145 int           hh_hm                 PROTO((double hh, int *hr, double *mn));
    146 int           day_to_sec            PROTO((char *string, time_t *second));
    147 int           hms_to_sec            PROTO((char *string, time_t *second));
    148 char         *ohana_sec_to_hms      PROTO((time_t second));
    149 char         *ohana_sec_to_day      PROTO((time_t second));
    150 
    151 char         *meade_deg_to_str      PROTO((double deg));
    152 char         *meade_ra_to_str       PROTO((double deg));
    153 char         *meade_dec_to_str      PROTO((double deg));
    154 char         *strptime              PROTO((const char *s, const char *format, struct tm *tm));
    155131time_t        GetTimeReference      PROTO((char *reference));
    156132int           GetTimeUnits          PROTO((char *name));

  • trunk/Ohana/src/uniphot/src/ConfigInit.c

    r27790 r33654  
    77  double ZERO_POINT;
    88  char  *config, *file;
    9   char CatdirPhotcodeFile[256];
    10   char MasterPhotcodeFile[256];
    119
    1210  /*** load configuration info ***/
    … …  
    2321  ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
    2422  ScanConfig (config, "CATFORMAT",              "%s",  0, CATFORMAT);
    25   ScanConfig (config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
    2623
    2724  sprintf (ImageCat, "%s/Images.dat", CATDIR);
    … …  
    4340  if (*CATMODE == 0) strcpy (CATMODE, "RAW");
    4441  if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
    45 
    46   /* XXX this does not yet write out the master photcode table */
    47   sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR);
    48   if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile, TRUE)) {
    49     fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile);
    50     exit (1);
    51   }
    52 
    5342}
    5443

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

    r31160 r33654  
    7272}
    7373
    74 int args_setphot (int argc, char **argv) {
    75 
    76   int N;
    77 
    78   VERBOSE = FALSE;
    79   if ((N = get_argument (argc, argv, "-v"))) {
    80     VERBOSE = TRUE;
    81     remove_argument (N, &argc, argv);
    82   }
    83 
    84   UPDATE = FALSE;
    85   if ((N = get_argument (argc, argv, "-update"))) {
    86     remove_argument (N, &argc, argv);
    87     UPDATE = TRUE;
    88   }
    89 
    90   if (argc != 2) {
    91     fprintf (stderr, "ERROR: USAGE: setphot (zptfile) [options]\n");
    92     exit (2);
    93   }
    94 
    95   return (TRUE);
    96 }
    97 
    9874int args_setfwhm (int argc, char **argv) {
    9975

  • trunk/Ohana/src/uniphot/src/convert.c

    r31160 r33654  
    11# include "uniphot.h"
    2 # define _XOPEN_SOURCE /* glibc2 (strptime) needs this */
    3 # include <time.h>
    4 
    5 /** additional time / coordinate conversions not supplied by libohana **/
    62
    73time_t GetTimeReference (char *reference) {
    … …  
    3632  return (Units);
    3733}
    38 
    39 int hh_hms (double hh, int *hr, int *mn, double *sc) {
    40 
    41   int N, flag;
    42 
    43   flag = SIGN(hh);
    44   hh = fabs(hh);
    45 
    46   // rationalize hh to range -24.0 < hh < 24.0
    47   if (hh >= 24.0) {
    48     N = (int)(hh/24.0);
    49     hh -= 24.0*N;
    50   }
    51 
    52   *hr = (int) hh;
    53   *mn = (int) 60*(hh - *hr);
    54   *sc = 3600.0*(hh - *hr - *mn / 60.0);
    55   if (*sc > 59.99) {
    56     *sc = 0.0;
    57     *mn += 1.0;
    58   }
    59   *hr *= flag;
    60   return (TRUE);
    61 }
    62  
    63 int dd_dms (double dd, int *dg, int *mn, double *sc) {
    64 
    65   int flag;
    66 
    67   flag = SIGN(dd);
    68   dd = fabs (dd);
    69   *dg = (int) dd;
    70   *mn = (int) 60*(dd - *dg);
    71   *sc = 3600.0*(dd - *dg - *mn/60.0);
    72   if (*sc > 59.99) {
    73     *sc = 0;
    74     *mn += 1.0;
    75   }
    76   *dg *= flag;
    77   return (TRUE);
    78 }
    79  
    80 int hms_format (char *line, double value) {
    81 
    82   int hr, mn;
    83   double sc;
    84 
    85   hh_hms (value, &hr, &mn, &sc);
    86   hr = (int) value;
    87   if (isnan (value))
    88     sprintf (line, "xx:xx:xx.xx");
    89   else {
    90     if (value < 0) {
    91       sprintf (line, "-%02d:%02d:%05.2f", abs(hr), mn, sc);
    92     } else {
    93       sprintf (line, "+%02d:%02d:%05.2f", hr, mn, sc);
    94     }
    95   }     
    96   return (TRUE);
    97 }
    98 
    99 int dms_format (char *line, double value) {
    100 
    101   int dg, mn;
    102   double sc;
    103 
    104   dd_dms (value, &dg, &mn, &sc);
    105   if (value < 0) {
    106     sprintf (line, "-%02d:%02d:%05.2f", abs(dg), mn, sc);
    107   } else {
    108     sprintf (line, "+%02d:%02d:%05.2f", dg, mn, sc);
    109   }
    110   return (TRUE);
    111 }
    112 
    113 /***** convert 00:00:00 or 00:00 to 0 - 86400 ****/
    114 int hms_to_sec (char *string, time_t *second) {
    115  
    116   char *p;
    117   struct tm time;
    118 
    119   p = strptime (string, "%H:%M:%S", &time);
    120   if (p != NULL) goto valid;
    121 
    122   p = strptime (string, "%H:%M", &time);
    123   if (p != NULL) goto valid;
    124 
    125   return (FALSE);
    126    
    127 valid:
    128   if (*p) return (FALSE);
    129   *second = time.tm_hour*3600 + time.tm_min*60 + time.tm_sec;
    130   return (TRUE);
    131 }
    132 
    133 /***** convert Mon[@00:00:00] or 00:00 to 0 - 86400*7 ****/
    134 int day_to_sec (char *string, time_t *second) {
    135  
    136   char *p;
    137   struct tm time;
    138 
    139   bzero (&time, sizeof(time));
    140   p = strptime (string, "%A@%H:%M:%S", &time);
    141   if (p != NULL) goto valid;
    142 
    143   p = strptime (string, "%A@%H:%M", &time);
    144   if (p != NULL) goto valid;
    145 
    146   p = strptime (string, "%A@%H", &time);
    147   if (p != NULL) goto valid;
    148 
    149   p = strptime (string, "%A", &time);
    150   if (p != NULL) goto valid;
    151 
    152   return (FALSE);
    153 
    154 valid:
    155   if (*p) return (FALSE);
    156   *second = time.tm_wday*86400 + time.tm_hour*3600 + time.tm_min*60 + time.tm_sec;
    157   return (TRUE);
    158 }
    159 
    160 /***** convert seconds to HH:MM:SS ****/
    161 char *ohana_sec_to_hms (time_t second) {
    162  
    163   struct tm *gmt;
    164   char *line;
    165 
    166   ALLOCATE (line, char, 64);
    167   gmt   = gmtime (&second);
    168   sprintf (line, "%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    169   return (line);
    170 }
    171 
    172 /***** convert seconds to Day@HH:MM:SS ****/
    173 char *ohana_sec_to_day (time_t second) {
    174  
    175   struct tm *gmt;
    176   char *line;
    177 
    178   ALLOCATE (line, char, 64);
    179   gmt   = gmtime (&second);
    180   switch (gmt[0].tm_wday) {
    181     case 0:
    182       sprintf (line, "Sun@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    183       break;
    184     case 1:
    185       sprintf (line, "Mon@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    186       break;
    187     case 2:
    188       sprintf (line, "Tue@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    189       break;
    190     case 3:
    191       sprintf (line, "Wed@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    192       break;
    193     case 4:
    194       sprintf (line, "Thu@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    195       break;
    196     case 5:
    197       sprintf (line, "Fri@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    198       break;
    199     case 6:
    200       sprintf (line, "Sat@%02d:%02d:%02d", gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec);
    201       break;
    202   }
    203   return (line);
    204 }
    205 
    206 int hh_hm (double hh, int *hr, double *mn) {
    207 
    208   int flag;
    209 
    210   flag = SIGN(hh);
    211   hh = fabs (hh);
    212 
    213   *mn = 60.0*(hh - (int)hh);
    214   *hr = (int) hh;
    215   *hr *= flag;
    216   return (TRUE);
    217 }
    218 
    219 char *meade_deg_to_str (double deg) {
    220 
    221   int hr;
    222   double mn;
    223   char *line;
    224 
    225   ALLOCATE (line, char, 16);
    226 
    227   hh_hm (deg, &hr, &mn);
    228 
    229   sprintf (line, "%03d:%04.1f", abs(hr), mn);
    230   return (line);
    231 }
    232 
    233 char *meade_ra_to_str (double deg) {
    234 
    235   int hr;
    236   double mn;
    237   char *line;
    238 
    239   ALLOCATE (line, char, 16);
    240 
    241   hh_hm (deg/15.0, &hr, &mn);
    242 
    243   sprintf (line, "%02d:%04.1f", abs(hr), mn);
    244   return (line);
    245 }
    246 
    247 char *meade_dec_to_str (double deg) {
    248 
    249   int hr;
    250   double mn;
    251   char *line;
    252 
    253   ALLOCATE (line, char, 16);
    254 
    255   hh_hm (deg, &hr, &mn);
    256 
    257   if (deg < 0) {
    258     sprintf (line, "-%02d:%04.1f", abs(hr), mn);
    259   } else {
    260     sprintf (line, "+%02d:%04.1f", hr, mn);
    261   }     
    262   return (line);
    263 }
    264 
    265 /* convert UNIX time to a value referenced to the TimeReference in the given unit */
    266 double TimeValue (time_t time, time_t TimeReference, int TimeFormat) {
    267 
    268   double value, dt;
    269 
    270   dt = (time > TimeReference) ? (time - TimeReference) : -1 * (double)(TimeReference - time);
    271   switch (TimeFormat) {
    272   case TIME_JD:
    273     value = time / 86400.0 + 2440587.5;
    274     break;
    275   case TIME_MJD:
    276     value = time / 86400.0 + 40587.0;
    277     break;
    278   case TIME_DAYS:
    279     value = dt / 86400.0;
    280     break;
    281   case TIME_HOURS:
    282     value = dt / 3600.0;
    283     break;
    284   case TIME_MINUTES:
    285     value = dt / 60.0;
    286     break;
    287   case TIME_SECONDS:
    288   default:
    289     value = dt;
    290     break;
    291   }
    292   return (value);
    293 }
    294  
    295 /* convert time value referenced to the TimeReference in the given unit to UNIX time */
    296 time_t TimeRef (double value, time_t TimeReference, int TimeFormat) {
    297 
    298   int dt;
    299   time_t time;
    300 
    301   switch (TimeFormat) {
    302   case TIME_JD:
    303     time = (value -  2440587.5) * 86400.0;
    304     return (time);
    305     break;
    306   case TIME_MJD:
    307     time = (value -  40587.0) * 86400.0;
    308     return (time);
    309     break;
    310   case TIME_DAYS:
    311     dt = value * 86400.0;
    312     break;
    313   case TIME_HOURS:
    314     dt = value * 3600.0;
    315     break;
    316   case TIME_MINUTES:
    317     dt = value * 60.0;
    318     break;
    319   case TIME_SECONDS:
    320   default:
    321     dt = value;
    322     break;
    323   }
    324 
    325   time = TimeReference + dt;
    326   return (time);
    327 }
    328 
    329 /* times may be in forms as:
    330  * 20040200450s (N seconds since 1970.0)
    331  * 2440900.232j (julian date)
    332  * 99/02/23,03:22:18 (date string)
    333  * (separators may be anything except space, +, -)
    334  * 99:02:15:12:23:30
    335  * 99:02:15:12h23m30s
    336  */
    337 

  • trunk/Ohana/src/uniphot/src/initialize.c

    r31160 r33654  
    2121}
    2222
    23 void initialize_setphot (int argc, char **argv) {
    24 
    25   /* are these set correctly? */
    26   ConfigInit (&argc, argv);
    27   args_setphot (argc, argv);
    28 }
    29 
    3023void initialize_setfwhm (int argc, char **argv) {
    3124

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

    r29001 r33654  
    1111}
    1212
    13 Image *load_images_setphot (FITS_DB *db, off_t *Nimage) {
     13Image *load_images_setfwhm (FITS_DB *db, off_t *Nimage) {
    1414
    1515  Image *image;

  • trunk/Ohana/src/uniphot/src/load_zpt_table.c

    r28717 r33654  
    1 # include "uniphot.h"
     1# include "setphot.h"
    22
    33ZptTable *load_zpt_table (char *filename, int *nzpts) {
    … …  
    4141  return zpts;
    4242}
     43
     44/* Slightly more generic loader than the ubercal version.  Still assumes Nfilters x Nseasons
     45   the input file must contain the following:
     46   PHU Header : NFILTER, NSEASON, NCHIP_X, NCHIP_Y, NCELL_X, NCELL_Y, TS0_nnnn (season nnnn start mjd), TS1_nnnn (season nnnn end mjd)
     47   
     48   NSEASON * 2 extensions with
     49   TABLE : mjd, zpt; Header: FILTER
     50   IMAGE : 1D array of flat offsets
     51*/
     52
     53ZptTable *load_zpt_ubercal(char *filename, int *nzpts, FlatCorrectionTable *flatcorrTable) {
     54
     55  int i, nfilter, nseason, ix, iy, ixc, iyc, Ncol;
     56  off_t Nrow;
     57  char type[16], filter[80];
     58  int Nzpts, NZPTS;
     59  ZptTable *zpts;
     60  Header header;
     61  Header theader;
     62  Matrix matrix;
     63  FTable ftable;
     64
     65  // parameters describing the flat-field correction
     66  int NSEASON;
     67  int NFILTER;
     68  int NCHIP_X;
     69  int NCHIP_Y;
     70  int NCELL_X;
     71  int NCELL_Y;
     72  int CHIP_DX;
     73  int CHIP_DY;
     74
     75  *nzpts = 0;
     76
     77  FILE *f = fopen (filename, "r");
     78  if (!f) {
     79    fprintf (stderr, "ERROR: cannot open zpt table file %s\n", filename);
     80    exit (1);
     81  }
     82
     83  /* load in PHU segment (ignore) */
     84  if (!gfits_fread_header (f, &header)) {
     85    if (VERBOSE) fprintf (stderr, "can't read Flat Correction header\n");
     86    fclose (f);
     87    return (NULL);
     88  }
     89  if (!gfits_fread_matrix (f, &matrix, &header)) {
     90    if (VERBOSE) fprintf (stderr, "can't read Flat Correction matrix\n");
     91    gfits_free_header (&header);
     92    fclose (f);
     93    return (NULL);
     94  }
     95
     96  Nzpts = 0;
     97  NZPTS = 0;
     98  ALLOCATE (zpts, ZptTable, NZPTS);
     99
     100  // this function would be better if we read the list of filters, seasons, and the dimensions from the header
     101  // for current testing, make fake smfs that correspond to specific chips, filter, and mjd ranges?
     102
     103  // the simple files from Eddie have no internal metadata describing the corrections,
     104  // so they must be manually encoded
     105
     106  // hard-wired values which describe the ubercal analysis
     107  if (NO_METADATA) {
     108    NFILTER = 5;
     109    NSEASON = 4;
     110    NCHIP_X = 8;
     111    NCHIP_Y = 8;
     112    NCELL_X = 2;
     113    NCELL_Y = 2;
     114    CHIP_DX = 2*2424;  // != 4880
     115    CHIP_DY = 2*2430;  // != 4864
     116    // Note that Eddie has identified the center of the boundary between cells xy3n and
     117    // xy4n and used that to split the chip.  I am setting CHIP_DX,DY to force the same split
     118  } else {
     119    if (!gfits_scan (&header, "NSEASON", "%d", 1, &NSEASON)) {
     120      fprintf (stderr, "cannot find NSEASON in header of %s\n", filename);
     121      fclose (f);
     122      return NULL;
     123    }
     124    if (!gfits_scan (&header, "NFILTER", "%d", 1, &NFILTER)) {
     125      fprintf (stderr, "cannot find NFILTER in header of %s\n", filename);
     126      fclose (f);
     127      return NULL;
     128    }
     129    if (!gfits_scan (&header, "NCHIP_X", "%d", 1, &NCHIP_X)) {
     130      fprintf (stderr, "cannot find NCHIP_X in header of %s\n", filename);
     131      fclose (f);
     132      return NULL;
     133    }
     134    if (!gfits_scan (&header, "NCHIP_Y", "%d", 1, &NCHIP_Y)) {
     135      fprintf (stderr, "cannot find NCHIP_Y in header of %s\n", filename);
     136      fclose (f);
     137      return NULL;
     138    }
     139    if (!gfits_scan (&header, "NCELL_X", "%d", 1, &NCELL_X)) {
     140      fprintf (stderr, "cannot find NCELL_X in header of %s\n", filename);
     141      fclose (f);
     142      return NULL;
     143    }
     144    if (!gfits_scan (&header, "NCELL_Y", "%d", 1, &NCELL_Y)) {
     145      fprintf (stderr, "cannot find NCELL_Y in header of %s\n", filename);
     146      fclose (f);
     147      return NULL;
     148    }
     149    if (!gfits_scan (&header, "CHIP_DX", "%d", 1, &CHIP_DX)) {
     150      fprintf (stderr, "cannot find CHIP_DX in header of %s\n", filename);
     151      fclose (f);
     152      return NULL;
     153    }
     154    if (!gfits_scan (&header, "CHIP_DY", "%d", 1, &CHIP_DY)) {
     155      fprintf (stderr, "cannot find CHIP_DY in header of %s\n", filename);
     156      fclose (f);
     157      return NULL;
     158    }
     159  }
     160
     161  flatcorrTable->Nseason = NSEASON;
     162  ALLOCATE (flatcorrTable->tstart, e_time, NSEASON);
     163  ALLOCATE (flatcorrTable->tstop,  e_time, NSEASON);
     164
     165  // ubercal hard-coded values
     166  char filters_uc[5][3] = {"g", "r", "i", "z", "y"};
     167  double tstart_uc[] = {55000.0, 55296.0, 55327.0, 55662.0};
     168  double tstop_uc[]  = {55296.0, 55327.0, 55662.0, 60000.0};
     169
     170  double mjdstart, mjdstop;
     171  for (i = 0; i < NSEASON; i++) {
     172    if (NO_METADATA) {
     173      mjdstart = tstart_uc[i];
     174      mjdstop  = tstop_uc[i];
     175  } else {
     176      char name[9];
     177      snprintf (name, 9, "TS0_%04d", i);
     178      if (!gfits_scan (&header, name, "%lf", 1, &mjdstart)) {
     179        fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     180        fclose (f);
     181        return NULL;
     182      }
     183
     184      snprintf (name, 9, "TS1_%04d", i);
     185      if (!gfits_scan (&header, name, "%lf", 1, &mjdstop)) {
     186        fprintf (stderr, "cannot find %s in header of %s\n", name, filename);
     187        fclose (f);
     188        return NULL;
     189      }
     190    }
     191    flatcorrTable->tstart[i] = ohana_mjd_to_sec(mjdstart);
     192    flatcorrTable->tstop[i]  = ohana_mjd_to_sec(mjdstop);
     193  }
     194
     195  // we have 5 filters, and 4 flat-field correction sets for each
     196  flatcorrTable->Ncorr = NFILTER*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y;
     197  flatcorrTable->Nimage = NFILTER*NSEASON*NCHIP_X*NCHIP_Y;
     198
     199  ALLOCATE (flatcorrTable->corr, FlatCorrection, flatcorrTable->Ncorr);
     200  ALLOCATE (flatcorrTable->image, FlatCorrectionImage, flatcorrTable->Nimage);
     201  memset (flatcorrTable->corr, 0, flatcorrTable->Ncorr*sizeof(FlatCorrection));
     202
     203  int corrID = 1;
     204  int Nimage = 0;
     205  ftable.header = &theader;
     206  for (nfilter = 0; nfilter < NFILTER; nfilter++) {
     207    // *** load the ZERO POINT table ***
     208
     209    // load data for this header
     210    if (!gfits_load_header (f, &theader)) return (NULL);
     211
     212    // read the fits table bytes
     213    if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL);
     214 
     215    // skip over remaining bytes in data segment
     216    fseeko (f, ftable.datasize - ftable.validsize, SEEK_CUR);
     217
     218    // need to create and assign to flat-field correction
     219    double *mjd = gfits_get_bintable_column_data (&theader, &ftable, "mjd_obs", type, &Nrow, &Ncol);
     220    assert (!strcmp(type, "double"));
     221
     222    double *zp = gfits_get_bintable_column_data (&theader, &ftable, "zp", type, &Nrow, &Ncol);
     223    assert (!strcmp(type, "double"));
     224     
     225    // float *zperr = gfits_get_bintable_column_data (&theader, &ftable, "resid", type, &Nrow, &Ncol);
     226    // assert (!strcmp(type, "float"));
     227     
     228    NZPTS += Nrow;
     229    REALLOCATE (zpts, ZptTable, NZPTS);
     230    for (i = 0; i < Nrow; i++) {
     231      zpts[i+Nzpts].time = ohana_mjd_to_sec(mjd[i]);
     232      zpts[i+Nzpts].zpt = zp[i];
     233      zpts[i+Nzpts].zpt_err = 0.0;
     234      zpts[i+Nzpts].found = FALSE;
     235    }
     236    Nzpts += Nrow;
     237
     238    // *** load the flat-field correction image ***
     239
     240    // the image contains the flat-field corrections for a specific filter
     241
     242    // load data for this header
     243    if (!gfits_load_header (f, &header)) return (NULL);
     244
     245    if (NO_METADATA) {
     246      strcpy (filter, filters_uc[nfilter]);
     247    } else {
     248      if (!gfits_scan (&header, "FILTER", "%s", 1, filter)) {
     249        fprintf (stderr, "cannot find FILTER in header of %s\n", filename);
     250        fclose (f);
     251        return NULL;
     252      }
     253    }
     254   
     255    // read the fits table bytes
     256    double *offset64 = NULL;
     257    float *offset32 = NULL;
     258    int use32 = TRUE;
     259   
     260    if (!gfits_fread_matrix (f, &matrix, &header)) return (NULL);
     261    switch (header.bitpix) {
     262      case -32: // float
     263        offset32 = (float *) matrix.buffer;
     264        use32 = TRUE;
     265        break;
     266      case -64: // double
     267        offset64 = (double *) matrix.buffer;
     268        use32 = FALSE;
     269        break;
     270      default:
     271        fprintf (stderr, "invalid bitpix for flat-field correction image: %d\n", header.bitpix);
     272        exit (2);
     273    }
     274
     275    // XXX the initial hacked-together table from Eddie is missing the last 2 elements.  they should be zero
     276    if (NO_METADATA) {
     277      matrix.Naxis[0] += 2;
     278      REALLOCATE (offset64, double, matrix.Naxis[0]);
     279      matrix.buffer = (char *) offset64;
     280    }
     281
     282    for (nseason = 0; nseason < NSEASON; nseason++) { // seasons
     283      for (iy = 0; iy < NCHIP_Y; iy++) { // y-chip
     284        for (ix = 0; ix < NCHIP_X; ix++) { // x-chip
     285
     286          // photcode name
     287          char photname[64];
     288          snprintf (photname, 64, "GPC1.%s.XY%d%d", filter, ix, iy);
     289          // note that the XY00, XY07, etc, chips will have photcode values of 0
     290
     291          flatcorrTable->image[Nimage].photcode = GetPhotcodeCodebyName(photname);
     292          flatcorrTable->image[Nimage].Nx = NCELL_X;
     293          flatcorrTable->image[Nimage].Ny = NCELL_Y;
     294          flatcorrTable->image[Nimage].ID = corrID;
     295          flatcorrTable->image[Nimage].DX = CHIP_DX;
     296          flatcorrTable->image[Nimage].DY = CHIP_DY;
     297          flatcorrTable->image[Nimage].tstart = flatcorrTable->tstart[nseason];
     298          flatcorrTable->image[Nimage].tstop  = flatcorrTable->tstop[nseason];
     299         
     300          int seq_full = -1; // sequence number within the full table (all filters concatenated together)
     301          int seq_filt = -1; // sequence number for just this filter
     302
     303          // XXX we should have a NCHIP_X * NCHIP_Y array of x and y parity values
     304
     305          // This enforces a 180 chip rotation for XY3n - XY7n & is only known to be valid for GPC1 (the XYnn names as well)
     306          for (iyc = 0; iyc < NCELL_Y; iyc++) {
     307            for (ixc = 0; ixc < NCELL_X; ixc++) {
     308              if (ix < 4) {
     309                // chips ix < 4 should be flipped
     310                seq_full = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     311                seq_filt =                                                   nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + NCELL_Y - 1 - iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + NCELL_X - 1 - ixc);
     312                // = .... + y_parity*iyc - (y_parity - 1)*(NCELL_Y - 1) / 2
     313              } else {
     314                seq_full = nfilter*NSEASON*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
     315                seq_filt =                                                   nseason*NCHIP_X*NCHIP_Y*NCELL_X*NCELL_Y + (iy*NCELL_Y + iyc)*NCHIP_X*NCELL_X + (ix*NCELL_X + ixc);
     316              }
     317              assert (seq_full > -1);
     318              assert (seq_filt > -1);
     319              assert (seq_full < flatcorrTable->Ncorr);
     320              if (matrix.Naxes == 1) {
     321                assert (seq_filt < matrix.Naxis[0]);
     322              }
     323              if (matrix.Naxes == 2) {
     324                assert (seq_filt < matrix.Naxis[0]*matrix.Naxis[1]);
     325              }
     326              assert (!flatcorrTable->corr[seq_full].ID);
     327              flatcorrTable->corr[seq_full].x = ixc;
     328              flatcorrTable->corr[seq_full].y = iyc;
     329              if (use32) {
     330                flatcorrTable->corr[seq_full].offset = offset32[seq_filt];
     331              } else {
     332                flatcorrTable->corr[seq_full].offset = offset64[seq_filt];
     333              }
     334              flatcorrTable->corr[seq_full].ID = corrID;
     335            }
     336          }
     337          corrID ++;
     338          Nimage ++;
     339        }
     340      }
     341    }         
     342  }
     343
     344  /*** convert from corr,image format to offsets ***/
     345  FlatCorrectionInternal (flatcorrTable);
     346
     347  fprintf (stderr, "loaded %d zero points\n", Nzpts);
     348
     349  *nzpts = Nzpts;
     350  return zpts;
     351}
     352
     353static short maxCode = 0;
     354static short *zpt_codeval = NULL;
     355static float *zpt_offsets = NULL;
     356static short *zpt_index = NULL;
     357
     358int parse_zpt_offsets (char *ZPT_OFFSET_FILTERS, char *ZPT_OFFSET_VALUES) {
     359
     360  if (!ZPT_OFFSET_FILTERS) return TRUE;
     361  assert (ZPT_OFFSET_FILTERS);
     362  assert (ZPT_OFFSET_VALUES);
     363
     364  char *filters = strcreate(ZPT_OFFSET_FILTERS);
     365  char *values = strcreate(ZPT_OFFSET_VALUES);
     366
     367  char *p1 = filters;
     368  char *p2 = values;
     369
     370  char *filter = NULL;
     371  char *value = NULL;
     372
     373  char *s1 = NULL;
     374  char *s2 = NULL;
     375
     376  // save an array of the zero point offsets and their photcodes
     377  int Ncode = 0;
     378  int NCODE = 8;
     379  ALLOCATE (zpt_codeval, short, NCODE);
     380  ALLOCATE (zpt_offsets, float, NCODE);
     381
     382  while (1) {
     383    filter = strtok_r (p1, ",", &s1);
     384    if (!filter) break;
     385
     386    value  = strtok_r (p2, ",", &s2);
     387    if (!value) {
     388      fprintf (stderr, "ERROR: mismatch between list of photcodes and list of offsets\n");
     389      exit (1);
     390    }
     391
     392    // do something here
     393    PhotCode *code = GetPhotcodebyName (filter);
     394    if (!code) {
     395      fprintf (stderr, "ERROR: unknown photcode %s\n", filter);
     396      exit (1);
     397    }
     398    if (code->type != PHOT_SEC) {
     399      fprintf (stderr, "ERROR: photcode %s is not an average (SEC) type\n", filter);
     400      exit (1);
     401    }
     402
     403    zpt_codeval[Ncode] = code->code;
     404    zpt_offsets[Ncode] = atof(value);
     405    Ncode ++;
     406
     407    if (Ncode >= NCODE) {
     408      NCODE += 8;
     409      REALLOCATE (zpt_codeval, short, NCODE);
     410      REALLOCATE (zpt_offsets, float, NCODE);
     411    }
     412   
     413    maxCode = MAX(code->code, maxCode);
     414
     415    p1 = NULL;
     416    p2 = NULL;
     417  }
     418   
     419  // generate an index to the zpt offsets by the photcode
     420  ALLOCATE (zpt_index, short, maxCode + 1);
     421
     422  int i;
     423  for (i = 0; i < maxCode + 1; i++) {
     424    zpt_index[i] = -1;
     425  }
     426
     427  for (i = 0; i < Ncode; i++) {
     428    short seq = zpt_codeval[i];
     429    if (zpt_index[seq] != -1) {
     430      fprintf (stderr, "duplicate photcode in -zpt-offset list\n");
     431      exit (2);
     432    }
     433    zpt_index[seq] = i;
     434  }
     435  return TRUE;
     436}
     437
     438float apply_zpt_offset (short code) {
     439
     440  if (!zpt_offsets) return 0.0;
     441
     442  // code is a primary photcode (but may be out of range, in which case no correction is supplied)
     443  if (code <= 0) return 0.0;
     444  if (code > maxCode) return 0.0;
     445
     446  short index = zpt_index[code];
     447  if (index == -1) return 0.0;
     448
     449  float offset = zpt_offsets[index];
     450  return offset;
     451}

  • trunk/Ohana/src/uniphot/src/match_zpts_to_images.c

    r29001 r33654  
    1 # include "uniphot.h"
     1# include "setphot.h"
    22
    33/* sort a coordinate pair (X,Y) and the associated index (S) */
    … …  
    3939  sort_zpts_by_time (zpts, Nzpts);
    4040
     41  if (RESET) {
     42    for (i = 0; i < Nimage; i++) {
     43      image[i].Mcal = 0.0;
     44      image[i].dMcal = NAN;
     45      image[i].flags &= ~ID_IMAGE_PHOTOM_NOCAL; // clear the NOCAL flag
     46      if (UBERCAL) {
     47        image[i].flags &= ~ID_IMAGE_PHOTOM_UBERCAL; // clear the UBERCAL flag
     48      }
     49    }
     50  }
     51
    4152  NImatch = 0; // matched images
    4253  for (i = Nz = 0; (i < Nimage) && (Nz < Nzpts); ) {
    4354
    44     if (i % 1000 == 0) fprintf (stderr, ".");
    45     if (Nz % 100 == 0) fprintf (stderr, "!");
     55    // if (i % 1000 == 0) fprintf (stderr, ".");
     56    // if (Nz % 100 == 0) fprintf (stderr, "!");
     57
     58    if (!isfinite(zpts[Nz].zpt)) {
     59      Nz++;
     60      continue;
     61    }
    4662
    4763    Ni = index[i];
    … …  
    6985
    7086    // we have a match: set zpt and record the match
    71     image[Ni].Mcal = SCALE*code[0].C - zpts[Nz].zpt;
     87    // is the zero point supplied nominally corrected or not?
     88    // UBERCAL includes 2.5log(exptime) + K*airmass in the zero point
     89    if (UBERCAL) {
     90      image[Ni].Mcal = SCALE*code[0].C - zpts[Nz].zpt + 2.5*log10(image[Ni].exptime) + code[0].K*(image[Ni].secz - 1.000);
     91      myAssert (isfinite(image[Ni].Mcal), "oops, ubercal made a nan image");
     92    } else {
     93      image[Ni].Mcal = SCALE*code[0].C - zpts[Nz].zpt;
     94    }
     95
     96    // if we have defined zero point offsets, then apply them here
     97    float offset = apply_zpt_offset (code[0].equiv);
     98    assert (isfinite(offset));
     99    image[Ni].Mcal += offset;
     100
    72101    image[Ni].dMcal = zpts[Nz].zpt_err;
    73102    image[Ni].flags &= ~ID_IMAGE_PHOTOM_NOCAL; // clear the NOCAL flag
     103    // image[Ni].flags |=  ID_IMAGE_PHOTOM_EXTERN; XXX do we want some flag like this?
     104    if (UBERCAL) {
     105      image[Ni].flags |=  ID_IMAGE_PHOTOM_UBERCAL;
     106    }
    74107    zpts[Nz].found = TRUE;
    75108    NImatch ++;

  • trunk/Ohana/src/uniphot/src/setfwhm.c

    r31160 r33654  
    2020
    2121  // load images
    22   image  = load_images_setphot (&db, &Nimage);
     22  image  = load_images_setfwhm (&db, &Nimage);
    2323  if (!UPDATE) dvo_image_unlock (&db);
    2424 

  • trunk/Ohana/src/uniphot/src/setphot.c

    r27790 r33654  
    1 # include "uniphot.h"
     1# include "setphot.h"
    22
    33int main (int argc, char **argv) {
    … …  
    88  ZptTable *zpts;
    99  Image *image;
     10  FlatCorrectionTable flatcorrTable;
    1011
    1112  /* get configuration info, args, lockfile */
    … …  
    1718  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
    1819
    19   zpts = load_zpt_table (argv[1], &Nzpts);
     20  if (UBERCAL) {
     21    zpts = load_zpt_ubercal (argv[1], &Nzpts, &flatcorrTable);
     22
     23    char flatcorrfile[MAX_PATH_LENGTH];
     24    int size = snprintf (flatcorrfile, MAX_PATH_LENGTH, "%s/flatcorr.fits", CATDIR);
     25    assert (size < MAX_PATH_LENGTH);
     26    FlatCorrectionSave(&flatcorrTable, flatcorrfile);
     27    // XXX should this program save any old copy of this file?
     28  } else {
     29    zpts = load_zpt_table (argv[1], &Nzpts);
     30  }
     31
     32  if (!zpts) Shutdown ("failed to load zero points, or empty table");
    2033
    2134  // load images
    … …  
    2538  match_zpts_to_images (image, Nimage, zpts, Nzpts);
    2639
    27   update_setphot(image, Nimage);
     40  if (UBERCAL) {
     41    match_flatcorr_to_images (image, Nimage, &flatcorrTable);
     42  }
    2843
    29   // write image table
     44  if (!IMAGES_ONLY) {
     45    status = update_dvo_setphot (image, Nimage, &flatcorrTable);
     46  }
     47
     48  // write image table (even if some remote clients failed)
    3049  if (UPDATE) {
    3150    dvo_image_save (&db, VERBOSE);
    3251  }
    3352
     53  if (!status) exit (1);
    3454  exit (0);
    3555}
    … …  
    3959   setphot (zpt_table)
    4060
    41  * load text table of zpts, time (photcode?)
    42  * load images
    43  * match images to zpts
    44  * set zpts
    45  * load catalogs
    46  * update detection (& averages?)
     61   setphot has 2 modes : basic & ubercal
     62
     63   * outline of basic mode:
     64
     65   ** load text table of zpts, time (photcode?)
     66   ** load images
     67   ** match images to zpts
     68   ** set zpts
     69   ** load catalogs
     70   ** update detection (& averages?)
     71
     72   * outline of ubercal mode:
     73
     74   ** load fits table of zpts & time + table of flat-field corrections
     75   ** load images
     76   ** match images to zpts
     77   ** set zpts
     78   ** load catalogs
     79   ** update detection (& averages?)
    4780
    4881 */

  • trunk/Ohana/src/uniphot/src/uniphot.c

    r31160 r33654  
    7474  if (!UPDATE) exit (0);
    7575
    76   update (&db, sgroup, Nsgroup);
     76  update_dvo_uniphot (&db, sgroup, Nsgroup);
    7777  dvo_image_unlock (&db);
    7878

  • trunk/Ohana/src/uniphot/src/update_catalog_setphot.c

    r28717 r33654  
    1 # include "uniphot.h"
    2 # include <glob.h>
     1# include "setphot.h"
    32
    4 int update_setphot (Image *image, off_t Nimage) {
     3void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage, FlatCorrectionTable *flatcorr) {
    54
    6   SkyRegion UserPatch;
    7   SkyTable *sky = NULL;
    8   SkyList *skylist = NULL;
    9   Catalog catalog;
    10   off_t i, maxID, *index;
     5  off_t i, j, found;
    116
    12   // create an index for the image IDs
    13   maxID = 0;
    14   for (i = 0; i < Nimage; i++) {
    15     maxID = MAX(maxID, image[i].imageID);
    16   }
    17   ALLOCATE (index, off_t, maxID + 1);
    18   for (i = 0; i < Nimage; i++) {
    19     if (!image[i].imageID) continue; // images with ID == 0 are virtual
    20     index[image[i].imageID] = i;
    21   }
    22  
    23 
    24   // load the current sky table (layout of all SkyRegions)
    25   sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
    26   SkyTableSetFilenames (sky, CATDIR, "cpt");
    27  
    28   // determine the populated SkyRegions overlapping the requested area (default depth)
    29   UserPatch.Rmin = 0;
    30   UserPatch.Rmax = 360;
    31   UserPatch.Dmin = -90;
    32   UserPatch.Dmax = +90;
    33   skylist = SkyListByPatch (sky, -1, &UserPatch);
    34 
    35   // update measurements for each populated catalog
    36   for (i = 0; i < skylist[0].Nregions; i++) {
    37 
    38     // set up the basic catalog info
    39     catalog.filename  = skylist[0].filename[i];
    40     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    41     catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    42     catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    43     catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    44 
    45     if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
    46       fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
    47       exit (1);
    48     }
    49     if (!catalog.Naves_disk) {
    50       if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
    51       dvo_catalog_unlock (&catalog);
    52       dvo_catalog_free (&catalog);
    53       continue;
    54     }
    55 
    56     update_catalog_setphot (&catalog, image, index, Nimage);
    57 
    58     if (!UPDATE) {
    59       dvo_catalog_unlock (&catalog);
    60       dvo_catalog_free (&catalog);
    61       continue;
    62     }
    63    
    64     if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename);
    65    
    66     dvo_catalog_save (&catalog, VERBOSE);
    67     dvo_catalog_unlock (&catalog);
    68     dvo_catalog_free (&catalog);
    69   }
    70   return (TRUE);
    71 }     
    72 
    73 void update_catalog_setphot (Catalog *catalog, Image *image, off_t *index, off_t Nimage) {
    74 
    75   off_t i, j, m, id, idx, found;
    76   float Mcal, dMcal;
     7  // if we are resetting, reset all flags
     8  DVOMeasureFlags photomFlags = ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_PHOTOM_UBERCAL;
    779
    7810  found = 0;   
    7911  for (i = 0; i < catalog[0].Naverage; i++) {
    8012
    81     m = catalog[0].average[i].measureOffset;
     13    off_t m = catalog[0].average[i].measureOffset;
    8214    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    83       idx = catalog[0].measure[m].imageID;
     15      off_t idx = catalog[0].measure[m].imageID;
    8416      if (idx <= 0) continue; // detections with imageID == 0 do not have a valid image (eg, ref photcode)
    8517
    86       id = index[idx];
    87       Mcal = image[id].Mcal;
    88       dMcal = image[id].dMcal;
     18      off_t id = index[idx];
     19      float Mcal = image[id].Mcal;
     20      float dMcal = image[id].dMcal;
     21      float Mcal_offset = 0.0;
    8922
    90       catalog[0].measure[m].Mcal = Mcal;
     23      // if we know about a flat-field correction, then we need to apply the sub-chip correction
     24      int flat_id = image[id].photom_map_id;
     25      if (flat_id > 0) {
     26          Mcal_offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[0].measure[m].Xccd, catalog[0].measure[m].Yccd);
     27      }
     28
     29      catalog[0].measure[m].Mcal = Mcal - Mcal_offset;
    9130      catalog[0].measure[m].dMcal = dMcal;
     31      myAssert(isfinite(catalog[0].measure[m].Mcal), "oops: ubercal made a nan");
     32
     33      if (RESET) {
     34        catalog[0].measure[m].dbFlags &= ~photomFlags;
     35      }
     36
     37      // if we are setting the zero points from an UBERCAL database, and this detection is from one of those images,
     38      // then tag the measurement as well
     39      if (UBERCAL && (image[id].flags & ID_IMAGE_PHOTOM_UBERCAL)) {
     40        catalog[0].measure[m].dbFlags |=  ID_MEAS_PHOTOM_UBERCAL;
     41      }
     42
    9243      found ++;
    9344    }
Note: See TracChangeset for help on using the changeset viewer.