IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 33651 for trunk


Ignore:
Timestamp:
Apr 1, 2012, 2:48:41 PM (14 years ago)
Author:
eugene
Message:

updates to relphot to support ubercal, distributed dvo; some optmimization (eg, image matching to area); cleanup names of static helper arrays in ImageOps, MosaicOps; cleanup rules for -reset options (relphot does not reset astrom bits); add flat-field correction; determine ubercal distance for images; better plot ranges; bcatalog limit density selecting objects by number of measurements; allow relphot to update the output catalog schema (catformat); make relphot -averages and relphot -update-objects equivalent in how low-quality objects are handled; record stage at which average is determined

Location:
trunk/Ohana/src/relphot
Files:
25 edited
5 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/Makefile

    r30616 r33651  
    1 default: relphot
     1default: relphot relphot_client
    22help:
    33        @echo "make options: relphot (default)"
    … …  
    1818
    1919relphot: $(BIN)/relphot.$(ARCH)
    20 install: $(DESTBIN)/relphot
     20relphot_client: $(BIN)/relphot_client.$(ARCH)
     21
     22install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client
    2123
    2224RELPHOT = \
    … …  
    2426$(SRC)/GridOps.$(ARCH).o         \
    2527$(SRC)/ImageOps.$(ARCH).o        \
     28$(SRC)/ImageSubset.$(ARCH).o     \
     29$(SRC)/BrightCatalog.$(ARCH).o   \
    2630$(SRC)/MosaicOps.$(ARCH).o       \
    2731$(SRC)/SetSignals.$(ARCH).o      \
    … …  
    4852$(RELPHOT): $(INC)/relphot.h
    4953$(BIN)/relphot.$(ARCH): $(RELPHOT)
     54
     55RELPHOT_CLIENT = \
     56$(SRC)/ConfigInit.$(ARCH).o      \
     57$(SRC)/GridOps.$(ARCH).o         \
     58$(SRC)/ImageOps.$(ARCH).o        \
     59$(SRC)/ImageSubset.$(ARCH).o     \
     60$(SRC)/BrightCatalog.$(ARCH).o   \
     61$(SRC)/MosaicOps.$(ARCH).o       \
     62$(SRC)/SetSignals.$(ARCH).o      \
     63$(SRC)/Shutdown.$(ARCH).o        \
     64$(SRC)/StarOps.$(ARCH).o         \
     65$(SRC)/bcatalog.$(ARCH).o        \
     66$(SRC)/args.$(ARCH).o            \
     67$(SRC)/help.$(ARCH).o            \
     68$(SRC)/plotstuff.$(ARCH).o       \
     69$(SRC)/liststats.$(ARCH).o       \
     70$(SRC)/initialize.$(ARCH).o      \
     71$(SRC)/load_catalogs.$(ARCH).o   \
     72$(SRC)/reload_catalogs.$(ARCH).o \
     73$(SRC)/relphot_objects.$(ARCH).o         \
     74$(SRC)/relphot_client.$(ARCH).o  \
     75$(SRC)/setExclusions.$(ARCH).o   \
     76$(SRC)/setMrelFinal.$(ARCH).o    \
     77$(SRC)/write_coords.$(ARCH).o
     78
     79$(RELPHOT_CLIENT): $(INC)/relphot.h
     80$(BIN)/relphot_client.$(ARCH): $(RELPHOT_CLIENT)

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

    r31668 r33651  
    77# define GRID_V2
    88# define NO_IMAGE -100
     9
     10// choose off_t or int depending on full-scale relphot analysis resources
     11// # define IDX_T off_t
     12# define IDX_T int
     13
     14typedef enum {
     15    MODE_NONE   = 0,
     16    MODE_LOAD   = 1,
     17    MODE_UPDATE = 2,
     18    MODE_UPDATE_OBJECTS = 3,
     19} ModeType;
    920
    1021typedef struct {
    … …  
    1829  short Xm;
    1930  float secz;
    20   char flags;
     31  float ubercalDist;
     32  unsigned int flags;
    2133  Coords coords;
    2234} Mosaic;
    … …  
    3042  double min;
    3143  double max;
     44  double Upper80;
     45  double Lower20;
    3246  double total;
    3347  int    Nmeas;
    3448} StatType;
    3549
     50typedef struct {
     51  AverageTiny *average;       // array of (minimal) average data
     52  MeasureTiny *measure;       // array of (minimal) measure data
     53  SecFilt     *secfilt;       // array of secfilt data (matched to average by Nsecfilt)
     54  off_t       Naverage;
     55  off_t       Nmeasure;
     56} BrightCatalog;
     57
     58typedef struct {
     59  float Mcal;
     60  float dMcal;
     61  unsigned int imageID;
     62  unsigned int photom_map_id;
     63  unsigned int flags;
     64  unsigned int tzero;
     65  unsigned char trate;
     66} ImageSubset;
     67
     68typedef struct {
     69  Catalog *catalog;           // array of catalogs generated
     70  int NCATALOG;               // number of catalogs allocated
     71  int Ncatalog;               // number of catalogs generated
     72  int Nsecfilt;               // number of catalogs generated
     73  off_t *NAVERAGE;            // allocated Averages per catalog
     74  off_t *NMEASURE;            // allocated Measures per catalog
     75  int   *index;               // lookup table catID -> catalog[i]
     76  int   *catIDs;              // lookup table catID <- catalog[i]
     77  int    maxID;               // max catID value to date
     78} CatalogSplitter;
     79
    3680/* global variables set in parameter file */
    37 char   ImageCat[256];
    38 char   ImageTemplate[256];
    39 char   CatTemplate[256];
    40 char   GSCFILE[256];
    41 char   CATDIR[256];
     81# define MAX_PATH_LENGTH 1024
     82char   ImageCat[MAX_PATH_LENGTH];
     83char   ImageTemplate[MAX_PATH_LENGTH];
     84char   CatTemplate[MAX_PATH_LENGTH];
     85char   GSCFILE[MAX_PATH_LENGTH];
     86char  *CATDIR;
    4287char   CATMODE[16];    /* raw, mef, split, mysql */
    4388char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
    44 char   CameraConfig[256];
    45 char   SKY_TABLE[256];
     89char   CameraConfig[MAX_PATH_LENGTH];
     90char   CAMERA[64];    /* eg, gpc1 */
     91char   SKY_TABLE[MAX_PATH_LENGTH];
    4692int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
     93
     94int          HOST_ID;
     95char        *HOSTDIR;
     96char        *IMAGES;
     97char        *BCATALOG;
     98ModeType     MODE;
    4799
    48100double MAG_LIM;
    … …  
    54106double MIN_ERROR;
    55107double IMFIT_SYS_SIGMA_LIM;
     108double CLOUD_TOLERANCE;
     109
     110int    PARALLEL;
     111int    PARALLEL_MANUAL;
     112int    PARALLEL_SERIAL;
    56113
    57114int    VERBOSE;
    … …  
    79136int    FREEZE_MOSAICS;
    80137int    USE_GRID;
     138int    KEEP_UBERCAL;
    81139char  *OUTROOT;
     140char  *UPDATE_CATFORMAT;
    82141int    PLOTDELAY;
    83142int    UpdateAverages;
     143
     144char  *PhotcodeList;
    84145
    85146int    RELPHOT_GRID_X;
    … …  
    101162double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
    102163
    103 int ImagSelect, ImagMin, ImagMax;
     164int ImagSelect;
     165double ImagMin, ImagMax;
    104166
    105167int DophotSelect, DophotValue;
    … …  
    112174
    113175SkyRegion UserPatch;
    114 int UserPatchSelect;
     176char *UserCatalog;
    115177
    116178int USE_BASIC_CHECK;
    … …  
    131193off_t        *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
    132194int           args                PROTO((int argc, char **argv));
     195int           args_client         PROTO((int argc, char **argv));
    133196int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
    134197void          clean_images        PROTO((void));
    135 void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
     198void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr));
    136199void          clean_mosaics       PROTO((void));
    137200void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
    … …  
    156219Coords       *getCoords           PROTO((off_t meas, int cat));
    157220off_t         getImageEntry       PROTO((off_t meas, int cat));
    158 float         getMcal             PROTO((off_t meas, int cat));
     221float         getMcal             PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
     222float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
    159223float         getMgrid            PROTO((off_t meas, int cat));
    160224float         getMmos             PROTO((off_t meas, int cat));
    161225float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
     226short         getUbercalDist      PROTO((off_t meas, int cat));
    162227Image        *getimage            PROTO((off_t N));
    163228Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
    164 void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
     229ImageSubset  *getimages_subset    PROTO((off_t *N));
     230void          global_stats        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    165231void          initGrid            PROTO((int dX, int dY));
    166232void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
    167233void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
     234void          initImagesSubset    PROTO((ImageSubset *input, off_t *line_number, off_t N));
    168235void          initImages          PROTO((Image *input, off_t *LineNumber, off_t N));
    169236void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
    … …  
    172239void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
    173240void          initialize          PROTO((int argc, char **argv));
     241void          initialize_client   PROTO((int argc, char **argv));
    174242void          initstats           PROTO((char *mode));
    175243int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
    176 Catalog       *load_catalogs      PROTO((SkyList *skylist, int *Ncatalog));
    177 SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect));
    178 Image         *select_images      PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
     244Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath));
     245Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog));
     246
     247SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
     248Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
    179249
    180250int           main                PROTO((int argc, char **argv));
    … …  
    185255void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
    186256void          plot_defaults       PROTO((Graphdata *graphdata));
    187 void          plot_grid           PROTO((Catalog *catalog));
     257void          plot_grid           PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    188258void          plot_images         PROTO((void));
    189259void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) );
    190260void          plot_mosaic_fields  PROTO((Catalog *catalog));
    191261void          plot_mosaics        PROTO((void));
    192 void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
     262void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
    193263void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
    194264void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
    195 void          reload_catalogs     PROTO((SkyList *skylist));
     265void          reload_catalogs     PROTO((SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath));
     266int           reload_catalogs_parallel PROTO((SkyList *sky));
    196267int           reload_images       PROTO((FITS_DB *db));
    197 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
    198 void          setMcal             PROTO((Catalog *catalog, int Poor));
     268int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
     269void          setMcal             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
    199270void          setMcalFinal        PROTO((void));
    200 int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
    201 void          setMgrid            PROTO((Catalog *catalog));
    202 int           setMmos             PROTO((Catalog *catalog, int Poor));
    203 int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
    204 void          setMrelFinal        PROTO((Catalog *catalog));
    205 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
     271int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
     272void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
     273int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
     274int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
     275void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
     276int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
    206277int           setMave             PROTO((Catalog *catalog, int Ncatalog));
    207278void          set_ZP              PROTO((double ZERO));
    208279int           setrefcode          PROTO((Image *image, off_t Nimage));
    209 void          skip_measurements   PROTO((Catalog *catalog, int pass));
     280void          skip_measurements   PROTO((Catalog *catalog, int pass, FlatCorrectionTable *flatcorr));
    210281void          sortA               PROTO((double *X, int N));
    211282void          sortB               PROTO((double *X, double *Y, int N));
    … …  
    213284void          sortD               PROTO((double *X, double *Y, double *Z, int N));
    214285StatType      statsImageM         PROTO((Catalog *catalog));
    215 StatType      statsImageN         PROTO((Catalog *catalog));
     286StatType      statsImageN         PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    216287StatType      statsImageX         PROTO((Catalog *catalog));
    217288StatType      statsImagedM        PROTO((Catalog *catalog));
    218289StatType      statsMosaicM        PROTO((Catalog *catalog));
    219 StatType      statsMosaicN        PROTO((Catalog *catalog));
     290StatType      statsMosaicN        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
    220291StatType      statsMosaicX        PROTO((Catalog *catalog));
    221292StatType      statsMosaicdM       PROTO((Catalog *catalog));
    222 StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode));
     293StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr));
    223294StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog, int Nsec));
    224295StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog, int Nsec));
    … …  
    226297void          wimages             PROTO((void));
    227298void          write_coords        PROTO((Header *header, Coords *coords));
    228 int relphot_objects (void);
     299int relphot_objects (int hostID, char *hostpath);
    229300
    230301void relphot_usage (void);
    231302void relphot_help (int argc, char **argv);
    232303
     304void relphot_client_usage (void);
     305void relphot_client_help (int argc, char **argv);
     306
    233307off_t getImageByID (off_t ID);
    234308
    … …  
    236310int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);
    237311
    238 int populate_tiny_values (Catalog *catalog);
    239 int free_tiny_values (Catalog *catalog);
     312BrightCatalog *BrightCatalogLoad(char *filename);
     313int BrightCatalogSave(char *filename, BrightCatalog *catalog);
     314BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog);
     315CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt);
     316int BrightCatalogSplitFree (CatalogSplitter *catalogs);
     317int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog);
     318
     319int ImageSubsetSave(char *filename, ImageSubset *image, off_t Nimage);
     320ImageSubset *ImageSubsetLoad(char *filename, off_t *nimage);

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

    r30616 r33651  
    3333  GetConfig (config, "IMAGE_GOOD_FRACTION",    "%lf", 0, &IMAGE_GOOD_FRACTION);
    3434
    35   GetConfig (config, "GSCFILE",                "%s",  0, GSCFILE);
    36   GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
     35  // force CATDIR to be absolute (so parallel mode will work)
     36  char *tmpcatdir = NULL;
     37  ALLOCATE (tmpcatdir, char, MAX_PATH_LENGTH);
     38  GetConfig (config, "CATDIR",                 "%s",  0, tmpcatdir);
     39  CATDIR = abspath (tmpcatdir, MAX_PATH_LENGTH);
     40  free (tmpcatdir);
     41
     42  GetConfig (config, "CAMERA",                  "%s",  0, CAMERA);
     43
     44  GetConfig (config, "GSCFILE",                 "%s",  0, GSCFILE);
    3745  ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
    3846  ScanConfig (config, "CATFORMAT",              "%s",  0, CATFORMAT);

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

    r31450 r33651  
    284284
    285285/* direct (non-iterative) solution for Mgrid values for all grid bins */
    286 void setMgridDirect (Catalog *catalog, int Ncatalog) {
     286void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    287287
    288288  int **gotstar, **gridmeas;
    … …  
    399399
    400400        // skip images marked as BAD
    401         Mcal = getMcal  (m, c);
     401        Mcal = getMcal  (m, c, flatcorr, catalog);
    402402        if (isnan(Mcal)) {
    403403          Ncal ++;
    … …  
    522522
    523523/* determine Mgrid values for all grid bins */
    524 void setMgrid (Catalog *catalog) {
     524void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    525525
    526526  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
    … …  
    554554        continue;
    555555      }
    556       Mcal = getMcal  (m, c);
     556      Mcal = getMcal  (m, c, flatcorr, catalog);
    557557      if (isnan(Mcal)) {
    558558        Ncal ++;
    … …  
    609609}
    610610
    611 void plot_grid (Catalog *catalog) {
     611void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    612612
    613613  int i, j, m, c, n, N, Narea;
    … …  
    641641        continue;
    642642      }
    643       Mcal  = getMcal  (m, c);
     643      Mcal  = getMcal  (m, c, flatcorr, catalog);
    644644      if (isnan(Mcal)) continue;
    645645      Mmos  = getMmos  (m, c);

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

    r31668 r33651  
    11# include "relphot.h"
    2 
    3 static off_t       **bin;     // link from catalog,measure to image         
    4 static off_t       **clist;   // catalog which supplied measurement on image
    5 static off_t       **mlist;   // measure reference for measurement on image 
    6 static off_t        *Nlist;   // number of measurements on image             
    7 static off_t        *NLIST;   // allocated number of measurements on image   
    8 
    9 static Image        *image;   // list of available images 
     2void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
     3
     4// the MeasureToImage, ImageToCatalog, and ImageToMeasure arrays are a substantial part of the memory footprint.
     5// total data volume is currently:
     6// a) Naverage*sizeof(AverageTiny)      [averages]       : 32
     7// b) Naverage*sizeof(Secfilt)*Nsecfilt [secfilt values] : 8 * 32
     8// b) Nmeasure*sizeof(Measure)          [measurements]   : 72
     9// c) Nmeasure*sizeof(IDX_T)*3          [image idx]      : 3 * 16
     10// d) Nmeasure*sizeof(IDX_T)*3          [mosaic idx]     : 3 * 16
     11// e) Nimage*sizeof(Image)              [image data]     : 360
     12
     13// for 3pi analysis in 2012, we have Nmeasure ~ 20 x Naverage
     14// so, each measurement costs ~14.5B (Ave + Sec) + 72B (Meas) + 96B (idx)!
     15// we are using off_t (64bit) to avoid the 32bit limit of an int, but
     16// if we really had 2^31 measurements in a single analysis, we would be using > 350GB of ram...
     17// until we reach that point, it is sort of silly to use IDX_T = off_t here.
     18
     19// in fact, we are safer than this, because the number of detections per table is rarely so large.
     20
     21// with IDX_T = int, each measurement costs ~14.5B (Ave + Sec) + 72B (Meas) + 48B (idx)!
     22
     23// elsewhere (not in relphot_images), we need to use off_t because a single catalog
     24
     25// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage)
    1026static off_t        Nimage;   // number of available images
     27static Image        *image;   // array of available images 
     28
     29// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row
     30// each image belongs
    1131static off_t       *LineNumber; // match of subset to full image table
    1232
    … …  
    1535static off_t        *imageIdx; // list of index for image IDs
    1636
     37// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas])
     38// each image has N_onImage[ImageIndex] measurements
     39static off_t        *N_onImage;   // actual number of measurements on image         
     40static off_t        *N_ONIMAGE;   // allocated number of measurements on image   
     41
     42// relationships between the measure,catalog set and the images:
     43static IDX_T       **MeasureToImage; // image index from measure,catalog   : MeasureToImage[cat][meas] = ImageIndex
     44static IDX_T       **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
     45static IDX_T       **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
     46
     47// MeasureToImage was 'bin'
     48// ImageToCatalog was 'clist'
     49// ImageToMeasure was 'mlist'
     50
     51// N_onImage was 'Nlist'
     52// N_ONIMAGE was 'NLIST'
     53
    1754Image *getimages (off_t *N, off_t **line_number) {
    1855
    … …  
    4279  }
    4380  llsortpair (imageIDs, imageIdx, Nimage);
     81}
     82
     83void initImagesSubset (ImageSubset *input, off_t *line_number, off_t N) {
     84
     85  off_t i;
     86
     87  // we have been given an ImageSubset array, containing a reduced set of image fields
     88  // create full a Image array and save the needed values
     89  ALLOCATE (image, Image, N);
     90  for (i = 0; i < N; i++) {
     91    image[i].imageID       = input[i].imageID      ;
     92    image[i].photom_map_id = input[i].photom_map_id;
     93    image[i].flags         = input[i].flags        ;
     94    image[i].Mcal          = input[i].Mcal         ;
     95    image[i].dMcal         = input[i].dMcal        ;
     96    image[i].tzero         = input[i].tzero        ;
     97    image[i].trate         = input[i].trate        ;
     98  }
     99  LineNumber = line_number;
     100  Nimage = N;
     101
     102  ALLOCATE (imageIDs, off_t, Nimage);
     103  ALLOCATE (imageIdx, off_t, Nimage);
     104
     105  for (i = 0; i < Nimage; i++) {
     106    imageIdx[i] = i;
     107    imageIDs[i] = image[i].imageID;
     108  }
     109  llsortpair (imageIDs, imageIdx, Nimage);
     110}
     111
     112ImageSubset *getimages_subset (off_t *N) {
     113
     114  *N = Nimage;
     115
     116  off_t i;
     117
     118  // we have been given an ImageSubset array, containing a reduced set of image fields
     119  // create full a Image array and save the needed values
     120  ImageSubset *subset = NULL;
     121  ALLOCATE (subset, ImageSubset, Nimage);
     122  for (i = 0; i < Nimage; i++) {
     123    subset[i].imageID       = image[i].imageID      ;
     124    subset[i].photom_map_id = image[i].photom_map_id;
     125    subset[i].flags         = image[i].flags        ;
     126    subset[i].Mcal          = image[i].Mcal         ;
     127    subset[i].dMcal         = image[i].dMcal        ;
     128    subset[i].tzero         = image[i].tzero        ;
     129    subset[i].trate         = image[i].trate        ;
     130  }
     131  return subset;
    44132}
    45133
    … …  
    71159void initImageBins (Catalog *catalog, int Ncatalog, int doImageList) {
    72160
    73   off_t i, j;
    74 
    75   ALLOCATE (bin, off_t *, Ncatalog);
     161  IDX_T i, j;
     162
     163  ALLOCATE (MeasureToImage, IDX_T *, Ncatalog);
    76164  for (i = 0; i < Ncatalog; i++) {
    77     ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
    78     for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
     165    ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1));
     166    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1;
    79167  }
    80168
    81169  if (doImageList) {
    82     ALLOCATE (Nlist, off_t, Nimage);
    83     ALLOCATE (NLIST, off_t, Nimage);
    84     ALLOCATE (clist, off_t *, Nimage);
    85     ALLOCATE (mlist, off_t *, Nimage);
     170    ALLOCATE (N_onImage, off_t, Nimage);
     171    ALLOCATE (N_ONIMAGE, off_t, Nimage);
     172    ALLOCATE (ImageToCatalog, IDX_T *, Nimage);
     173    ALLOCATE (ImageToMeasure, IDX_T *, Nimage);
    86174
    87175    for (i = 0; i < Nimage; i++) {
    88       Nlist[i] = 0;
    89       NLIST[i] = 100;
    90       ALLOCATE (clist[i], off_t, NLIST[i]);
    91       ALLOCATE (mlist[i], off_t, NLIST[i]);
     176      N_onImage[i] = 0;
     177      N_ONIMAGE[i] = 30;
     178      ALLOCATE (ImageToCatalog[i], IDX_T, N_ONIMAGE[i]);
     179      ALLOCATE (ImageToMeasure[i], IDX_T, N_ONIMAGE[i]);
    92180    }
    93181  }
    … …  
    99187
    100188  for (i = 0; i < Ncatalog; i++) {
    101     free (bin[i]);
    102   }
    103   free (bin);
     189    free (MeasureToImage[i]);
     190  }
     191  free (MeasureToImage);
    104192
    105193  if (doImageList) {
    106194    for (i = 0; i < Nimage; i++) {
    107       free (clist[i]);
    108       free (mlist[i]);
    109     }
    110     free (clist);
    111     free (mlist);
    112     free (Nlist);
    113     free (NLIST);
    114   }
    115 }
    116 
    117 /* select all image equivalent to the active photcode set */
     195      free (ImageToCatalog[i]);
     196      free (ImageToMeasure[i]);
     197    }
     198    free (ImageToCatalog);
     199    free (ImageToMeasure);
     200    free (N_onImage);
     201    free (N_ONIMAGE);
     202  }
     203}
     204
     205/* select all images equivalent to the active photcode set */
    118206void findImages (Catalog *catalog, int Ncatalog, int doImageList) {
    119207
    … …  
    210298
    211299  // index for (catalog, measure) -> image
    212   bin[cat][meas] = idx;
     300  MeasureToImage[cat][meas] = idx;
    213301
    214302  if (doImageList) {
    215303    // index for image, Nentry -> catalog
    216     clist[idx][Nlist[idx]] = cat;
     304    ImageToCatalog[idx][N_onImage[idx]] = cat;
    217305
    218306    // index for image, Nentry -> measure
    219     mlist[idx][Nlist[idx]] = meas;
    220     Nlist[idx] ++;
    221 
    222     if (Nlist[idx] == NLIST[idx]) {
    223       NLIST[idx] += 100;
    224       REALLOCATE (clist[idx], off_t, NLIST[idx]);
    225       REALLOCATE (mlist[idx], off_t, NLIST[idx]);
     307    ImageToMeasure[idx][N_onImage[idx]] = meas;
     308    N_onImage[idx] ++;
     309
     310    if (N_onImage[idx] == N_ONIMAGE[idx]) {
     311      N_ONIMAGE[idx] += 30;
     312      REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]);
     313      REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]);
    226314    }   
    227315  }
    … …  
    234322  off_t i;
    235323
    236   i = bin[cat][meas];
     324  i = MeasureToImage[cat][meas];
    237325  return (i);
    238326}
    239327
    240 float getMcal (off_t meas, int cat) {
     328// returns image.Mcal - ff(x,y)
     329float getMcal (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
    241330
    242331  off_t i;
    243   float value;
    244 
    245   i = bin[cat][meas];
     332  float value, offset;
     333
     334  i = MeasureToImage[cat][meas];
    246335  if (i == -1) return (NAN);
    247336
    248337  if (image[i].flags & IMAGE_BAD) return (NAN); 
    249338  value = image[i].Mcal;
     339  offset = 0.0;
     340
     341  // to do this, I need to pass in the catalog and flatcorr pointers
     342  int flat_id = image[i].photom_map_id;
     343  if (flat_id) {
     344    offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
     345  }
     346  value -= offset;
     347
    250348  return (value);
    251349}
    252350
     351// returns image.Mcal - ff(x,y)
     352short getUbercalDist (off_t meas, int cat) {
     353
     354  off_t i;
     355  short distance;
     356
     357  i = MeasureToImage[cat][meas];
     358  if (i == -1) return (1000);
     359
     360  if (image[i].flags & IMAGE_BAD) return (1000); 
     361  distance = image[i].ubercalDist; // was dummy3 in structure
     362  return (distance);
     363}
     364
     365float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
     366
     367  off_t i = MeasureToImage[cat][meas];
     368  if (i == -1) return (NAN);
     369
     370  float offset = 0.0;
     371
     372  // to do this, I need to pass in the catalog and flatcorr pointers
     373  int flat_id = image[i].photom_map_id;
     374  if (flat_id) {
     375    offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
     376  }
     377
     378  return (offset);
     379}
     380
    253381Coords *getCoords (off_t meas, int cat) {
    254382
    255383  off_t i;
    256384
    257   i = bin[cat][meas];
     385  i = MeasureToImage[cat][meas];
    258386  if (i == -1) return (NULL);
    259387  return (&image[i].coords);
    … …  
    261389
    262390/* determine Mcal values for all images */
    263 void setMcal (Catalog *catalog, int PoorImages) {
    264 
    265   off_t i, j, m, c, n, N, Nmax, mark, bad;
    266   int Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys, Nbright;
    267   float Msys, Mrel, Mmos, Mgrid, McalBright, McalBright2;
    268   double *list, *dlist;
     391void setMcal (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     392
     393  off_t i, j, m, c, n, N, Nmax;
     394  int mark, bad, Nfew, Nbad, Nmos, Nrel, Ngrid, Nsys, Nbright;
     395  float Msys, Mrel, Mmos, Mgrid, Mflat;
     396  double *list, *dlist, *Mlist, *dMlist;
    269397  StatType stats;
    270398
    271399  if (FREEZE_IMAGES) return;
    272400
     401  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
     402
    273403  int Nsecfilt = GetPhotcodeNsecfilt ();
    274404
    275405  if (PoorImages) {
    276     IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     406    // XXX use bad stars and measurements for PoorImages? or not?
     407    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
     408    IMAGE_BAD = 0;
    277409  }
    278410
    279411  Nmax = 0;
    280412  for (i = 0; i < Nimage; i++) {
    281     Nmax = MAX (Nmax, Nlist[i]);
     413    Nmax = MAX (Nmax, N_onImage[i]);
    282414  }
    283415  ALLOCATE (list, double, Nmax);
    284416  ALLOCATE (dlist, double, Nmax);
     417  ALLOCATE (Mlist, double, Nmax);
     418  ALLOCATE (dMlist, double, Nmax);
    285419
    286420  Nfew = Nbad = Nmos = Ngrid = Nrel = Nsys = 0;
    287 
    288   // counters to measure the bright-end scatter
    289   McalBright = McalBright2 = 0.0;
    290   Nbright = 0;
    291421
    292422  for (i = 0; i < Nimage; i++) {
    … …  
    298428    }     
    299429
     430    // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
     431    // we optionally do not recalibrate images with UBERCAL zero points
     432    if (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
     433      image[i].ubercalDist = 0; // was dummy3
     434      if (KEEP_UBERCAL) continue;
     435    }
     436
     437    int minUbercalDist = 1000;
     438
     439    // number of stars to measure the bright-end scatter
     440    Nbright = 0;
     441
    300442    N = 0;
    301     for (j = 0; j < Nlist[i]; j++) {
     443    for (j = 0; j < N_onImage[i]; j++) {
    302444     
    303       m = mlist[i][j];
    304       c = clist[i][j];
     445      m = ImageToMeasure[i][j];
     446      c = ImageToCatalog[i][j];
    305447     
    306448      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    … …  
    324466      }
    325467     
     468      // image.Mcal is not supposed to include the flat-field correction, so we need to
     469      // apply that offset as well here for this image (in other words, each detection is
     470      // being compared to the model, excluding the zero point, Mcal.  The model includes
     471      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
     472
     473      Mflat = getMflat (m, c, flatcorr, catalog);
     474
    326475      n = catalog[c].measureT[m].averef;
    327476      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
    … …  
    330479        continue;
    331480      }
    332       list[N] = Msys - Mrel - Mmos - Mgrid;
     481
     482      PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     483      if (!code) goto skip;
     484      if (code->equiv < 1) goto skip;
     485      int Nsec = GetPhotcodeNsec (code->equiv);
     486      if (Nsec == -1) goto skip;
     487      minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
     488
     489      // fprintf (stderr, "%1d, %3d : %3d, %3d : %10.6f %10.6f : %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n", (int) i, (int) j, (int) c, (int) m, catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mmos, Mgrid, Mflat);
     490
     491    skip:
     492      list[N] = Msys - Mrel - Mmos - Mgrid + Mflat;
    333493      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    334494      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
    335         McalBright += list[N];
    336         McalBright2 += SQ(list[N]);
     495        Mlist[Nbright] = list[N];
     496        dMlist[Nbright] = dlist[N];
    337497        Nbright ++;
    338498      }
    339499      N++;
    340500    }
    341     /* Nlist[i] is all measurements, N is good measurements */
     501    /* N_onImage[i] is all measurements, N is good measurements */
     502
     503    // if (VERBOSE2) fprintf (stderr, "meas skipped: (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ngrid, Nrel, Nsys);
    342504
    343505    /* too few good measurements or too many bad measurements */
    344506    if (!PoorImages) {
    345       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
     507      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onImage[i]);
    346508      if (mark) {
    347509        image[i].flags |= ID_IMAGE_PHOTOM_FEW;
    … …  
    351513      }     
    352514    }
    353    
     515
    354516    liststats (list, dlist, N, &stats);
    355 
    356     float CLOUD_TOLERANCE = 0.01;
    357517    image[i].Mcal  = stats.mean;
    358518    image[i].dMcal = stats.error;
    359     image[i].dMagSys = stats.sigma;
    360519    image[i].nFitPhotom = N;
    361520    image[i].Xm    = 100.0*log10(stats.chisq);
    362521
     522    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
     523
     524    // bright end scatter
     525    liststats (Mlist, dMlist, Nbright, &stats);
     526    image[i].dMagSys = stats.sigma;
     527
    363528    if (image[i].Mcal < -CLOUD_TOLERANCE) {
    364529      image[i].Mcal = 0.0;
    365530    }
     531
     532    // minUbercalDist calculated here is the min value for any star owned by this image
     533    // since this particular image is tied to that star, bump its distance by 1
     534    image[i].ubercalDist = minUbercalDist + 1;
    366535  }
    367536  free (list);
    368537  free (dlist);
     538  free (Mlist);
     539  free (dMlist);
    369540
    370541  fprintf (stderr, "%d images marked having too few measurements (Nbad: %d, Nmos: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", Nfew, Nbad, Nmos, Ngrid, Nrel, Nsys);
    … …  
    395566  ALLOCATE (dlist, double, Nimage);
    396567
     568  // measure stats for Mcal and dMcal
    397569  for (i = N = 0; i < Nimage; i++) {
    398570    if (image[i].flags & IMAGE_BAD) continue;
    399     mlist[N] = fabs (image[i].Mcal);
     571    mlist[N] = image[i].Mcal;
    400572    slist[N] = image[i].dMcal;
    401573    dlist[N] = 1;
    … …  
    412584  Nmark = 0;
    413585  for (i = 0; i < Nimage; i++) {
     586    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
     587    if (KEEP_UBERCAL && (image[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
     588
    414589    mark = FALSE;
    415590    image[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
    … …  
    430605}
    431606
     607void plot_setMcal (double *list, int Npts, StatType *stats, float clouds) {
     608
     609  off_t i;
     610  double *xlist;
     611  Graphdata graphdata;
     612
     613  if (!PLOTSTUFF) return;
     614
     615  ALLOCATE (xlist, double, Npts);
     616
     617  /**** Mcal vs seq ****/
     618  float minMcal    = +100.0;
     619  float maxMcal    = -100.0;
     620  for (i = 0; i < Npts; i++) {
     621    xlist[i] = i;
     622    minMcal = MIN (list[i], minMcal);
     623    maxMcal = MAX (list[i], maxMcal);
     624  }
     625
     626  float McalRange  = MAX(1.2*(maxMcal - minMcal), 0.21);
     627  float McalCenter = 0.5*(maxMcal + minMcal);
     628
     629  plot_defaults (&graphdata);
     630  graphdata.xmin = -1;
     631  graphdata.xmax = Npts + 1;
     632  graphdata.ymin = McalCenter - 0.5*McalRange;
     633  graphdata.ymax = McalCenter + 0.5*McalRange;
     634  plot_list (&graphdata, xlist, list, Npts, "sequence vs Mcal", "%s.seq.png", OUTROOT);
     635
     636  free (xlist);
     637}
     638
    432639void plot_images () {
    433640
    … …  
    443650
    444651  /**** dMcal vs airmass ****/
     652  float minAirmass = 1000.0;
     653  float maxAirmass =    0.0;
     654  float minMcal    = +100.0;
     655  float maxMcal    = -100.0;
     656  float mindMcal    = +100.0;
     657  float maxdMcal    = -100.0;
    445658  for (i = 0; i < Nimage; i++) {
    446659    Mlist[i] = image[i].Mcal;
    447660    dlist[i] = image[i].dMcal;
    448661    xlist[i] = image[i].secz;
    449   }
     662    minAirmass = MIN (image[i].secz, minAirmass);
     663    maxAirmass = MAX (image[i].secz, maxAirmass);
     664    minMcal = MIN (image[i].Mcal, minMcal);
     665    maxMcal = MAX (image[i].Mcal, maxMcal);
     666    mindMcal = MIN (image[i].dMcal, mindMcal);
     667    maxdMcal = MAX (image[i].dMcal, maxdMcal);
     668  }
     669
     670  float AirmassRange  = MAX(1.2*(maxAirmass - minAirmass), 0.25);
     671  float AirmassCenter = 0.5*(maxAirmass + minAirmass);
     672
     673  float McalRange  = MAX(1.2*(maxMcal - minMcal), 0.21);
     674  float McalCenter = 0.5*(maxMcal + minMcal);
     675
     676  float dMcalRange  = MAX(1.2*(maxdMcal - mindMcal), 0.21);
     677  float dMcalCenter = 0.5*(maxdMcal + mindMcal);
    450678
    451679  plot_defaults (&graphdata);
    452   graphdata.ymin = PlotdMmin;
    453   graphdata.ymax = PlotdMmax;
     680  graphdata.xmin = AirmassCenter - 0.5*AirmassRange;
     681  graphdata.xmax = AirmassCenter + 0.5*AirmassRange;
     682  graphdata.ymin = McalCenter - 0.5*McalRange;
     683  graphdata.ymax = McalCenter + 0.5*McalRange;
    454684  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
     685
    455686  plot_defaults (&graphdata);
     687  graphdata.xmin = McalCenter - 0.5*McalRange;
     688  graphdata.xmax = McalCenter + 0.5*McalRange;
     689  graphdata.ymin = dMcalCenter - 0.5*dMcalRange;
     690  graphdata.ymax = dMcalCenter + 0.5*dMcalRange;
    456691  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", "%s.Mcal.dMcal.png", OUTROOT);
    457692
    … …  
    471706  plot_defaults (&graphdata);
    472707  graphdata.style = 1;
     708  graphdata.xmin = dMcalCenter - 0.5*dMcalRange;
     709  graphdata.xmax = dMcalCenter + 0.5*dMcalRange;
    473710  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT);
    474711
    … …  
    478715}
    479716
    480 StatType statsImageN (Catalog *catalog) {
     717StatType statsImageN (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    481718
    482719  off_t i, j, m, c, n, N;
    … …  
    496733
    497734    N = 0;
    498     for (j = 0; j < Nlist[i]; j++) {
    499 
    500       m = mlist[i][j];
    501       c = clist[i][j];
    502 
    503       Mcal  = getMcal  (m, c);
     735    for (j = 0; j < N_onImage[i]; j++) {
     736
     737      m = ImageToMeasure[i][j];
     738      c = ImageToCatalog[i][j];
     739
     740      Mcal  = getMcal  (m, c, flatcorr, catalog);
    504741      if (isnan(Mcal)) continue;
    505742      Mmos  = getMmos  (m, c);

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

    r31668 r33651  
    11# include "relphot.h"
     2void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
     3
     4// see discussion in ImagesOps.c re: IDX_T
    25
    36// array of mosaic definition structures
    … …  
    58static Mosaic *mosaic;
    69
    7 // list of all images associated with a mosaic
    8 static off_t   *Nimlist;
    9 static off_t   **imlist; /* mosaic -> image[] */
    10 static off_t   **bin;    /* catalog, measure -> mosaic */
    11 
    12 // list of mosaics associated with an image
    13 static off_t   *mosimage;
    14 
    15 // list of mosaic associated with each image 
    16 static int     **clist;  /* mosaic -> catalog[] */
    17 static off_t   **mlist;  /* mosiac -> measure[] */
    18 static off_t    *Nlist;
    19 static off_t    *NLIST;
     10// relationships between the mosaics and their associated images
     11static off_t   *MosaicN_Image; // number of images associated with the given mosaic
     12static off_t  **MosaicToImage; // list of imagesa associated with the given mosaic
     13
     14// mosaic index for given image : ImageToMosaic[ImageIndex] = MosaicIndex (ImageIndex : 0 < Nimage)
     15static off_t   *ImageToMosaic;
     16
     17// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas])
     18// each mosaic has N_onMosaic[MosaicIndex] measurements
     19static off_t    *N_onMosaic;   // actual number of measurements on mosaic       
     20static off_t    *N_ONMOSAIC;   // allocated number of measurements on mosaic   
     21
     22// relationships between the measure,catalog set and the mosaics:
     23static off_t   **MeasureToMosaic; // Mosaic index from measure,catalog  : MeasureToMosaic[cat][meas] = MosaicIndex
     24static off_t   **MosaicToCatalog; // catalog for given measure on mosaic : MosaicCatalog[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex])
     25static off_t   **MosaicToMeasure; // measure for given measure on mosaic : MosaicMeasure[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex])
     26
     27// MeasureToMosaic was 'bin'
     28// MosaicToCatalog was 'clist'
     29// MosaicToMeasure was 'mlist'
     30
     31// N_onMosaic was 'Nlist'
     32// N_ONMOSAIC was 'NLIST'
     33
     34// ImageToMosaic was 'mosimage'
     35
     36// MosaicN_image was 'Nimlist'
     37// MosaicToImage was 'imlist'
    2038
    2139/* find mosaic frames (unique time periods & photcode name matches mosaic) */
    22 void initMosaics (Image *image, off_t Nimage) {
    23 
    24   off_t i, j, status, found, NMOSAIC, *NIMLIST;
     40void initMosaics_old (Image *image, off_t Nimage) {
     41
     42  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
    2543  unsigned int start, stop;
    2644  char *pname;
    … …  
    3250  ALLOCATE (mosaic, Mosaic, NMOSAIC);
    3351
    34   ALLOCATE (imlist,  off_t *, NMOSAIC);
    35   ALLOCATE (Nimlist, off_t,   NMOSAIC);
    36   ALLOCATE (NIMLIST, off_t,   NMOSAIC);
    37 
    38   ALLOCATE (mosimage, off_t, Nimage); // mosaic to which image belongs
     52  ALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
     53  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
     54  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
     55
     56  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
    3957
    4058  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
    … …  
    4563  /* generate list of unique mosaics */
    4664  for (i = 0; i < Nimage; i++) {
    47     mosimage[i] = -1;
     65    ImageToMosaic[i] = -1;
    4866
    4967    /* select valid mosaic images by photcode */
    … …  
    6482
    6583      // add reference from image to mosaic
    66       mosimage[i] = j;
     84      ImageToMosaic[i] = j;
    6785
    6886      /* add image to mosaic image list */
    69       imlist[j][Nimlist[j]] = i;
    70       Nimlist[j] ++;
    71       if (Nimlist[j] == NIMLIST[j]) {
    72         NIMLIST[j] += 10;
    73         REALLOCATE (imlist[j], off_t, NIMLIST[j]);
     87      MosaicToImage[j][MosaicN_Image[j]] = i;
     88      MosaicN_Image[j] ++;
     89      if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
     90        MosaicN_IMAGE[j] += 10;
     91        REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
    7492      }
    7593
    … …  
    88106    mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
    89107
     108    // XXX do we need to do something about flag consistency across a mosaic?
     109
    90110    /* add image to mosaic image list */
    91     NIMLIST[Nmosaic] = 10;
    92     Nimlist[Nmosaic] = 1;
    93     ALLOCATE (imlist[Nmosaic], off_t, NIMLIST[Nmosaic]);
    94     imlist[Nmosaic][0] = i;
     111    MosaicN_IMAGE[Nmosaic] = 10;
     112    MosaicN_Image[Nmosaic] = 1;
     113    ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
     114    MosaicToImage[Nmosaic][0] = i;
    95115
    96116    // add reference from image to mosaic
    97     mosimage[i] = Nmosaic;
     117    ImageToMosaic[i] = Nmosaic;
    98118   
    99119    Nmosaic ++;
    … …  
    101121      NMOSAIC += 10;
    102122      REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
    103       REALLOCATE (imlist,  off_t *, NMOSAIC);
    104       REALLOCATE (Nimlist, off_t,   NMOSAIC);
    105       REALLOCATE (NIMLIST, off_t,   NMOSAIC);
    106     }
    107   }
     123      REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
     124      REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
     125      REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
     126    }
     127  }
     128
     129  // free this or not?
     130  free (MosaicN_IMAGE);
    108131
    109132  initMosaicGrid (image, Nimage);
     133
     134  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
    110135  return;
     136}
     137
     138off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
     139
     140void sort_times (unsigned int *T, int N) {
     141
     142# define SWAPFUNC(A,B){ unsigned int tmp; \
     143  tmp = T[A]; T[A] = T[B]; T[B] = tmp; \
     144}
     145# define COMPARE(A,B)(T[A] < T[B])
     146
     147  OHANA_SORT (N, COMPARE, SWAPFUNC);
     148
     149# undef SWAPFUNC
     150# undef COMPARE
     151
     152}
     153
     154/* find mosaic frames (unique time periods & photcode name matches mosaic) */
     155void initMosaics (Image *image, off_t Nimage) {
     156
     157  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
     158  unsigned int start, stop, *startTimes, *startTimesMosaic;
     159  char *pname;
     160
     161  if (!MOSAIC_ZEROPT) return;
     162
     163  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
     164   * entry in the image table that represents this mosaic.  Instead, it is an
     165   * internal construct that defines a group of related images
     166   */
     167
     168  // generate a list of unique start times (these define the mosaics)
     169  ALLOCATE (startTimes, unsigned int, Nimage);
     170  for (i = 0; i < Nimage; i++) {
     171    startTimes[i] = image[i].tzero;
     172  }
     173  sort_times (startTimes, Nimage);
     174 
     175  Nmosaic = 0;
     176  NMOSAIC = 1000;
     177  ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
     178  startTimesMosaic[0] = startTimes[0];
     179
     180  for (i = 0; i < Nimage; i++) {
     181    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
     182      fprintf (stderr, "error?\n");
     183      abort();
     184    }
     185    if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
     186    Nmosaic ++;
     187    if (Nmosaic >= NMOSAIC) {
     188      NMOSAIC += 1000;
     189      REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
     190    }
     191    startTimesMosaic[Nmosaic] = startTimes[i];
     192  }
     193  Nmosaic ++;
     194
     195  // now I have a list of uniq start times, and they are in order
     196  // create the mosaic arrays for these times
     197  ALLOCATE (mosaic, Mosaic, Nmosaic);
     198
     199  ALLOCATE (MosaicToImage, off_t *, Nmosaic);
     200  ALLOCATE (MosaicN_Image, off_t,   Nmosaic);
     201  ALLOCATE (MosaicN_IMAGE, off_t,   Nmosaic);
     202
     203  for (i = 0; i < Nmosaic; i++) {
     204    /* a new mosaic, define ranges */
     205    mosaic[i].start = startTimesMosaic[i];
     206    mosaic[i].stop  = 0;
     207    mosaic[i].Mcal  = 0.0;
     208    mosaic[i].dMcal = 0.0;
     209    mosaic[i].dMsys = 0.0;
     210    mosaic[i].Xm    = 0.0;
     211    mosaic[i].flags = 0;
     212    mosaic[i].secz  = NAN;
     213    mosaic[i].photcode = 0;
     214
     215    MosaicN_IMAGE[i] = 10;
     216    MosaicN_Image[i] = 0;
     217    ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]);
     218    MosaicToImage[i][0] = -1;
     219  }
     220
     221  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
     222
     223  // assign each image to a mosaic
     224  for (i = 0; i < Nimage; i++) {
     225    ImageToMosaic[i] = -1;
     226
     227    /* select valid mosaic images by photcode */
     228    pname = GetPhotcodeNamebyCode (image[i].photcode);
     229    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
     230    if (status) continue;
     231
     232    /* set image time range */
     233    // start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
     234    // stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
     235
     236    start = image[i].tzero;
     237    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
     238
     239    j = findMosaic(startTimesMosaic, Nmosaic, start);
     240    if (j == -1) {
     241      fprintf (stderr, "error?\n");
     242      abort();
     243    }
     244
     245    // add reference from image to mosaic
     246    ImageToMosaic[i] = j;
     247
     248    // have we already found this mosaic?
     249    found = (MosaicN_Image[j] > 0);
     250
     251    /* add image to mosaic image list */
     252    MosaicToImage[j][MosaicN_Image[j]] = i;
     253    MosaicN_Image[j] ++;
     254    if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
     255      MosaicN_IMAGE[j] += 10;
     256      REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
     257    }
     258    if (found) continue;
     259   
     260    /* a new mosaic, define ranges */
     261    if (mosaic[j].start != start) {
     262      fprintf (stderr, "error?\n");
     263      abort();
     264    }
     265    mosaic[j].stop  = stop;
     266    mosaic[j].Mcal  = 0.0;
     267    mosaic[j].dMcal = 0.0;
     268    mosaic[j].dMsys = 0.0;
     269    mosaic[j].Xm    = 0.0;
     270    mosaic[j].flags  = image[i].flags;
     271    mosaic[j].secz  = image[i].secz;
     272    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
     273  }
     274
     275  // free this or not?
     276  free (MosaicN_IMAGE);
     277  free (startTimes);
     278  free (startTimesMosaic);
     279
     280  initMosaicGrid (image, Nimage);
     281
     282  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
     283  return;
     284}
     285
     286// use bisection to find the overlapping mosaic
     287off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) {
     288
     289  off_t Nlo, Nhi, N;
     290
     291  // find the last mosaic before start
     292  Nlo = 0; // startTimes[Nlo] guaranteed to be <= start
     293  Nhi = Nmosaic - 1; // startTimes[Nhi] guaranteed to be >= start
     294  while (Nhi - Nlo > 10) {
     295    N = 0.5*(Nlo + Nhi);
     296    if (startTimes[N] < start) {
     297      Nlo = MAX(N, 0);
     298    } else {
     299      Nhi = MIN(N, Nmosaic - 1);
     300    }
     301  }
     302
     303  // check for the matched mosaic starting from Nlo
     304  for (N = Nlo; N < Nmosaic; N++) {
     305    if (start > startTimes[N])  continue;
     306    return (N);
     307  }
     308  return (-1);
    111309}
    112310
    … …  
    129327    dS = 0.0;
    130328    Mcal = dMcal = Xm = 0;
    131     for (j = 0; j < Nimlist[i]; j++) {
    132       m = imlist[i][j];
     329    for (j = 0; j < MosaicN_Image[i]; j++) {
     330      m = MosaicToImage[i][j];
    133331      NX = image[m].NX;
    134332      NY = image[m].NY;
    … …  
    157355      dMcal += image[m].dMcal;
    158356      Xm    += image[m].Xm;
     357
     358      // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but
     359      // solve for a single offset for each exposure (Mosaic.Mcal). 
     360
     361      // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
     362
     363      // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
     364
    159365      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    160366      image[m].Mcal  = 0.0;
    … …  
    162368      image[m].Xm    = NAN_S_SHORT;
    163369    }
    164     dS /= Nimlist[i];
     370    dS /= MosaicN_Image[i];
    165371    strcpy (mosaic[i].coords.ctype, "RA---TAN");
    166372    mosaic[i].coords.crval1 = Rmin;
    … …  
    176382    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
    177383
    178     mosaic[i].Mcal  = Mcal / Nimlist[i];
    179     mosaic[i].dMcal = dMcal / Nimlist[i];
    180     mosaic[i].Xm    = Xm / Nimlist[i];
     384    mosaic[i].Mcal  = Mcal / MosaicN_Image[i];
     385    mosaic[i].dMcal = dMcal / MosaicN_Image[i];
     386    mosaic[i].Xm    = Xm / MosaicN_Image[i];
    181387  }
    182388  if (!USE_GRID) return;
    … …  
    197403  image = getimages (&Nimage, NULL);
    198404
     405  // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL
     406  // images? (no need to update them)
     407
    199408  // copy the mosaic results to the images.  set the mosaic Mcal to 0.0 since we have moved its
    200409  // impact to the images
    201410  for (i = 0; i < Nmosaic; i++) {
    202     for (j = 0; j < Nimlist[i]; j++) {
    203       im = imlist[i][j];
     411    for (j = 0; j < MosaicN_Image[i]; j++) {
     412      im = MosaicToImage[i][j];
    204413      image[im].Mcal += mosaic[i].Mcal;
    205414      image[im].dMcal = mosaic[i].dMcal;
    … …  
    221430  if (!MOSAIC_ZEROPT) return;
    222431
    223   ALLOCATE (bin, off_t *, Ncatalog);
     432  ALLOCATE (MeasureToMosaic, off_t *, Ncatalog);
    224433  for (i = 0; i < Ncatalog; i++) {
    225     ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
    226     for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
     434    ALLOCATE (MeasureToMosaic[i], off_t, MAX (catalog[i].Nmeasure, 1));
     435    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToMosaic[i][j] = -1;
    227436  }
    228437
    229438  if (doMosaicList) {
    230439    /* mosaic -> measure */
    231     ALLOCATE (Nlist, off_t,   Nmosaic);
    232     ALLOCATE (NLIST, off_t,   Nmosaic);
    233     ALLOCATE (clist, int *,  Nmosaic);
    234     ALLOCATE (mlist, off_t *, Nmosaic);
     440    ALLOCATE (N_onMosaic, off_t,   Nmosaic);
     441    ALLOCATE (N_ONMOSAIC, off_t,   Nmosaic);
     442    ALLOCATE (MosaicToCatalog, off_t *, Nmosaic);
     443    ALLOCATE (MosaicToMeasure, off_t *, Nmosaic);
    235444
    236445    for (i = 0; i < Nmosaic; i++) {
    237       Nlist[i] = 0;
    238       NLIST[i] = 100;
    239       ALLOCATE (clist[i], int,   NLIST[i]);
    240       ALLOCATE (mlist[i], off_t, NLIST[i]);
     446      N_onMosaic[i] = 0;
     447      N_ONMOSAIC[i] = 100;
     448      ALLOCATE (MosaicToCatalog[i], off_t, N_ONMOSAIC[i]);
     449      ALLOCATE (MosaicToMeasure[i], off_t, N_ONMOSAIC[i]);
    241450    }
    242451  }
    … …  
    251460
    252461  for (i = 0; i < Ncatalog; i++) {
    253     free (bin[i]);
    254   }
    255   free (bin);
     462    free (MeasureToMosaic[i]);
     463  }
     464  free (MeasureToMosaic);
    256465
    257466  if (doMosaicList) {
    258467    /* mosaic -> measure */
    259468    for (i = 0; i < Nmosaic; i++) {
    260       free (clist[i]);
    261       free (mlist[i]);
    262     }
    263     free (Nlist);
    264     free (NLIST);
    265     free (clist);
    266     free (mlist);
     469      free (MosaicToCatalog[i]);
     470      free (MosaicToMeasure[i]);
     471    }
     472    free (N_onMosaic);
     473    free (N_ONMOSAIC);
     474    free (MosaicToCatalog);
     475    free (MosaicToMeasure);
    267476  }
    268477}
    … …  
    310519  }
    311520
    312   mosID = mosimage[idx];
     521  mosID = ImageToMosaic[idx];
    313522  if (mosID < 0) {
    314523    Image *image = getimage(idx);
    … …  
    320529  {
    321530    Image *image = getimage(idx);
    322     unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
     531    // XXX we are now matching with just tzero.  be careful for cameras with drift
     532    // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
     533    unsigned int imageStart = image[0].tzero;
    323534    if (imageStart != mosaic[mosID].start) {
    324535      fprintf (stderr, "error in image to mosaic match\n");
    … …  
    327538  }
    328539
    329   bin[cat][meas] = mosID;
     540  MeasureToMosaic[cat][meas] = mosID;
    330541
    331542  if (doMosaicList) {
    332     clist[mosID][Nlist[mosID]] = cat;
    333     mlist[mosID][Nlist[mosID]] = meas;
    334     Nlist[mosID] ++;
     543    MosaicToCatalog[mosID][N_onMosaic[mosID]] = cat;
     544    MosaicToMeasure[mosID][N_onMosaic[mosID]] = meas;
     545    N_onMosaic[mosID] ++;
    335546   
    336     if (Nlist[mosID] == NLIST[mosID]) {
    337       NLIST[mosID] += 100;
    338       REALLOCATE (clist[mosID], int,   NLIST[mosID]);
    339       REALLOCATE (mlist[mosID], off_t, NLIST[mosID]);
     547    if (N_onMosaic[mosID] == N_ONMOSAIC[mosID]) {
     548      N_ONMOSAIC[mosID] += 100;
     549      REALLOCATE (MosaicToCatalog[mosID], off_t, N_ONMOSAIC[mosID]);
     550      REALLOCATE (MosaicToMeasure[mosID], off_t, N_ONMOSAIC[mosID]);
    340551    }   
    341552  }
    … …  
    350561  if (!MOSAIC_ZEROPT) return (0);
    351562
    352   i = bin[cat][meas];
     563  i = MeasureToMosaic[cat][meas];
    353564  if (i == -1) return (NAN);
    354565
    … …  
    358569}
    359570
    360 int setMmos (Catalog *catalog, int PoorImages) {
    361 
    362   off_t i, j, m, c, N, Nmax;
    363   int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt;
    364   float Msys, Mrel, Mcal, Mgrid;
     571int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {
     572
     573  off_t i, j, m, c, n, N, Nmax;
     574  int mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nbright;
     575  float Msys, Mrel, Mcal, Mgrid, Mflat;
    365576  double *list, *dlist, *Mlist, *dMlist;
    366577  StatType stats;
    … …  
    372583  image = getimages (&N, NULL);
    373584
    374   Nsecfilt = GetPhotcodeNsecfilt ();
     585  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
     586
     587  int Nsecfilt = GetPhotcodeNsecfilt ();
    375588
    376589  if (PoorImages) {
     590    // XXX use bad stars and measurements for PoorImages? or not?
     591    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    377592    IMAGE_BAD = 0;
    378593  }
    … …  
    380595  Nmax = 0;
    381596  for (i = 0; i < Nmosaic; i++) {
    382     Nmax = MAX (Nmax, Nlist[i]);
     597    Nmax = MAX (Nmax, N_onMosaic[i]);
    383598  }
    384599  ALLOCATE (list, double, Nmax);
    … …  
    397612    }     
    398613
     614    // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
     615    // we optionally do not recalibrate images with UBERCAL zero points
     616    if (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL) {
     617      mosaic[i].ubercalDist = 0;
     618      // propagate ubercalDist to the images
     619      for (j = 0; j < MosaicN_Image[i]; j++) {
     620        off_t im = MosaicToImage[i][j];
     621        image[im].ubercalDist = mosaic[i].ubercalDist;
     622        // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
     623      }
     624      if (KEEP_UBERCAL) continue;
     625    }
     626
     627    int minUbercalDist = 1000;
     628
     629    // number of stars to measure the bright-end scatter
     630    Nbright = 0;
     631
    399632    N = 0;
    400     for (j = 0; j < Nlist[i]; j++) {
     633    for (j = 0; j < N_onMosaic[i]; j++) {
    401634     
    402       m = mlist[i][j];
    403       c = clist[i][j];
     635      m = MosaicToMeasure[i][j];
     636      c = MosaicToCatalog[i][j];
    404637     
    405638      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
    … …  
    407640          continue;
    408641      }
    409       Mcal  = getMcal  (m, c);
     642      Mcal  = getMcal  (m, c, flatcorr, catalog);
    410643      if (isnan(Mcal)) {
    411644          Ncal++;
    … …  
    423656      }
    424657     
     658      // image.Mcal is not supposed to include the flat-field correction, so we need to
     659      // apply that offset as well here for this image (in other words, each detection is
     660      // being compared to the model, excluding the zero point, Mcal.  The model includes
     661      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
     662
     663      Mflat = getMflat (m, c, flatcorr, catalog);
     664
    425665      n = catalog[c].measureT[m].averef;
    426666      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
    … …  
    429669        continue;
    430670      }
    431       list[N]  = Msys - Mrel - Mcal - Mgrid;
     671
     672      PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
     673      if (!code) goto skip;
     674      if (code->equiv < 1) goto skip;
     675      int Nsec = GetPhotcodeNsec (code->equiv);
     676      if (Nsec == -1) goto skip;
     677      minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);
     678
     679    skip:
     680      list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
    432681      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    433       Mlist[N] = Msys;
    434       dMlist[N] = list[N];
     682      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
     683        Mlist[Nbright] = list[N];
     684        dMlist[Nbright] = dlist[N];
     685        Nbright ++;
     686      }
    435687      N++;
    436688    }
    437     /* Nlist[i] is all measurements, N is good measurements */
     689    /* N_onMosaic[i] is all measurements, N is good measurements */
    438690
    439691    /* too few good measurements or too many bad measurements (skip in PoorImages run) */
    440 
    441692    if (!PoorImages) {
    442       mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
     693      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[i]);
    443694      if (mark) {
    444         if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  Nlist[i]); }
     695        if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[MosaicToImage[i][0]].name,  i,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[i]); }
    445696        mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW;
    446697        Nfew ++;
    … …  
    450701    }
    451702    liststats (list, dlist, N, &stats);
    452     if (PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
    453 
    454     float CLOUD_TOLERANCE = 0.01;
     703    if (VERBOSE2 && PoorImages) fprintf (stderr, "Mmos: %f %f %d "OFF_T_FMT"\n", stats.mean, stats.sigma, stats.Nmeas,  N);
     704
    455705    mosaic[i].Mcal  = stats.mean;
    456706    mosaic[i].dMcal = stats.error;
    457     mosaic[i].dMsys = stats.sigma;
    458707    mosaic[i].nFitPhotom = N;
    459708    mosaic[i].Xm    = 100.0*log10(stats.chisq);
    460709
     710    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
     711
     712    // bright end scatter
     713    liststats (Mlist, dMlist, Nbright, &stats);
     714    mosaic[i].dMsys = stats.sigma;
     715
    461716    if (mosaic[i].Mcal < -CLOUD_TOLERANCE) {
    462717      mosaic[i].Mcal = 0.0;
     718    }
     719
     720    // minUbercalDist calculated here is the min value for any star owned by this image
     721    // since this particular image is tied to that star, bump its distance by 1
     722    mosaic[i].ubercalDist = minUbercalDist + 1;
     723
     724    // propagate ubercalDist to the images
     725    for (j = 0; j < MosaicN_Image[i]; j++) {
     726      off_t im = MosaicToImage[i][j];
     727      image[im].ubercalDist = mosaic[i].ubercalDist;
     728      // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    463729    }
    464730  }
    … …  
    481747// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
    482748// driven by the coupling of common stars.
     749// XXX this function is no longer used because we force significantly negative clouds to 0.0
    483750int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
    484751
    … …  
    515782
    516783    if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) {
    517       fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
     784      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    518785    }
    519786    if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) {
    … …  
    521788      // NOTE the negative sign: down below, we are going to add in the negative of Mcal
    522789      // to this image, and the propagated mean values for other images
    523       fprintf (stderr, "anti-clouds: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
     790      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    524791    }
    525792    if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) {
    526       fprintf (stderr, "cloudy    : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
     793      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    527794    }
    528795  }
    … …  
    579846
    580847        // find the source of this measurement (skip unassigned measurements)
    581         nMos = bin[i][m];
     848        nMos = MeasureToMosaic[i][m];
    582849        if (nMos == -1) continue;
    583850
    … …  
    641908    }
    642909    imageOffset[i] = dM / NSlist[i];
    643     // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]);
     910    // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]);
    644911  }
    645912
    646913  // for (i = 0; i < Nmosaic; i++) {
    647   //   fprintf (stderr, "correction: %s : %f\n", image[imlist[i][0]].name, imageOffset[i]);
     914  //   fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]);
    648915  // }
    649916
    … …  
    727994}
    728995
    729 StatType statsMosaicN (Catalog *catalog) {
     996StatType statsMosaicN (Catalog *catalog, FlatCorrectionTable *flatcorr) {
    730997
    731998  off_t i, j, m, c, n, N;
    … …  
    7461013
    7471014    N = 0;
    748     for (j = 0; j < Nlist[i]; j++) {
    749 
    750       m = mlist[i][j];
    751       c = clist[i][j];
    752 
    753       Mcal = getMcal  (m, c);
     1015    for (j = 0; j < N_onMosaic[i]; j++) {
     1016
     1017      m = MosaicToMeasure[i][j];
     1018      c = MosaicToCatalog[i][j];
     1019
     1020      Mcal = getMcal  (m, c, flatcorr, catalog);
    7541021      if (isnan(Mcal)) continue;
    7551022      Mgrid = getMgrid (m, c);
    … …  
    8331100  Nmark = Nscatter = Noffset = 0;
    8341101  for (i = 0; i < Nmosaic; i++) {
     1102    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
     1103    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
     1104
    8351105    mark = FALSE;
    8361106    if (mosaic[i].dMcal > MaxScatter) {
    … …  
    8721142  N = 0;
    8731143  for (i = 0; i < Nmosaic; i++)
    874     N = MAX (N, Nlist[i]);
     1144    N = MAX (N, N_onMosaic[i]);
    8751145
    8761146  ALLOCATE (xlist, double, N);
    … …  
    8811151    Xmin = Ymin = +360.0;
    8821152    Xmax = Ymax = -360.0;
    883     for (j = 0; j < Nlist[i]; j++) {
     1153    for (j = 0; j < N_onMosaic[i]; j++) {
    8841154     
    885       m = mlist[i][j];
    886       c = clist[i][j];
     1155      m = MosaicToMeasure[i][j];
     1156      c = MosaicToCatalog[i][j];
    8871157     
    8881158      if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
    … …  
    9591229
    9601230  Imax = 0;
    961   Nmax = Nimlist[0];
    962   for (i = 0; i < Nmosaic; i++) {
    963     if (Nimlist[i] > Nmax) {
     1231  Nmax = MosaicN_Image[0];
     1232  for (i = 0; i < Nmosaic; i++) {
     1233    if (MosaicN_Image[i] > Nmax) {
    9641234      Imax = i;
    965       Nmax = Nimlist[i];
     1235      Nmax = MosaicN_Image[i];
    9661236    }
    9671237  }
    … …  
    9691239  *refmosaic = &mosaic[Imax];
    9701240  *Nimage = Nmax;
    971   return (imlist[Imax]);
    972 }
     1241  return (MosaicToImage[Imax]);
     1242}

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

    r29754 r33651  
    3030  for (i = 0; i < 36; i++) {
    3131    switch (i) {
    32       /* can't redirect this signals */
     32      /* can't redirect these signals */
    3333    case SIGKILL:    /* kill -9: cannot be caught or ignored */
    3434    case SIGSTOP:    /* SIGSTOP: cannot be caught or ignored */
     35
    3536      /* ignore these signals */
    3637    case SIGCHLD:    /* child halted: ignore */

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

    r7080 r33651  
    11# include "relphot.h"
    22
    3 static FITS_DB *db;
     3static FITS_DB *db = NULL;
    44
    55void set_db (FITS_DB *in) {
    … …  
    2222
    2323  SetProtect (TRUE);
    24   gfits_db_close (db);
     24  if (db) gfits_db_close (db);
    2525  fprintf (stderr, "ERROR: addstar halted\n");
    2626  exit (1);

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

    r31668 r33651  
    4747}
    4848
    49 int setMrel (Catalog *catalog, int Ncatalog) {
     49int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    5050
    5151  off_t j, k, m;
    … …  
    5959  for (i = 0; i < Ncatalog; i++) {
    6060    for (j = 0; j < catalog[i].Naverage; j++) {
     61      // XXX accumulate all secfilt values in a single pass?
     62
     63      int minUbercalDist = 1000;
    6164
    6265      int Ns;
    … …  
    6669        int Nsec = GetPhotcodeNsec(thisCode);
    6770
    68         /* calculate the average value for a single star */
     71        /* calculate the average mag in this SEC photcode for a single star */
    6972
    7073        // skip bad stars
    7174        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
     75
     76        N = 0;
    7277        m = catalog[i].averageT[j].measureOffset;
    73 
    74         N = 0;
    7578        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
     79
    7680          // skip measurements that do not match the current photcode
    77           int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    78           if (ecode != thisCode) { continue; }
    79 
    80           if (catalog[i].measureT[m].dbFlags & MEAS_BAD) {
    81             Nbad ++;
    82             continue;
    83           }
    84           // XXX allow REF stars (no Image Entry) to be included in the calculation this
    85           // should be optionally set, and should allow for REF stars to be downweighted by
    86           // more than their reported errors.  how such info is carried is unclear...
     81          PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
     82          if (!code) continue;
     83          if (code->equiv != thisCode) { continue; }
     84
     85          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) { Nbad ++; continue; }
     86
    8787          if (getImageEntry (m, i) < 0) {
    8888            Mcal = Mmos = Mgrid = 0;
    8989          } else {
    90             Mcal  = getMcal  (m, i);
    91             if (isnan(Mcal)) {
    92               Ncal ++;
    93               continue;
    94             }
     90            Mcal  = getMcal  (m, i, flatcorr, catalog);
     91            if (isnan(Mcal)) { Ncal ++; continue; }
    9592            Mmos  = getMmos  (m, i);
    96             if (isnan(Mmos)) {
    97               Nmos ++;
    98               continue;
    99             }
     93            if (isnan(Mmos)) { Nmos ++; continue; }
    10094            Mgrid = getMgrid (m, i);
    101             if (isnan(Mgrid)) {
    102               Ngrid++;
    103               continue;
    104             }
     95            if (isnan(Mgrid)) { Ngrid++; continue; }
    10596          }
    10697
    10798          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
    108           if (isnan(Msys)) {
    109             Nsys++;
    110             continue;
    111           }
     99          if (isnan(Msys)) { Nsys++; continue; }
    112100          list[N] = Msys - Mcal - Mmos - Mgrid;
    113           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
     101
     102          int myUbercalDist = getUbercalDist(m,i);
     103          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     104
     105          // dlist gives the error, which is used as the weight in WT_MEAN.
     106          // we can modify the resulting weight in a few ways:
     107          // 1) MIN_ERROR guarantees a floor
     108          // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     109          // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     110          // 4) some reference photcode of some kind can be specified as fixed and have a high weight
     111          dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
     112
     113          // up-weight the ubercal values (or convergence can take a long time...)
     114          if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     115            dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
     116          }
    114117
    115118          // tie down reference photometry if the -refcode (code) option is selected
     119          // eg, -refcode g_SDSS
     120          // this probably makes no sense in the context of multifilter analysis
    116121          if (refPhotcode) {
    117             if (GetPhotcodeEquivCodebyCode(catalog[i].measureT[m].photcode) == refPhotcode[0].equiv) {
    118               // increase the weight by a factor of 100:
    119               dlist[N] = 0.01*catalog[i].measureT[m].dM;
     122            if (code->code == refPhotcode->code) {
     123              // tiny error -> large weight
     124              // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
     125              dlist[N] = 0.0001;
    120126            }
    121127          }
    … …  
    125131        // when performing the grid analysis, STAR_TOOFEW will be set to 1;
    126132        if (N <= STAR_TOOFEW) { /* too few measurements */
     133          // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].measureT[catalog[i].averageT[j].measureOffset].imageID, N, STAR_TOOFEW);
    127134          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
    128135          Nfew ++;
    … …  
    132139
    133140        liststats (list, dlist, N, &stats);
    134 
     141       
    135142        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    136143        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
    137         catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
     144        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
     145       
     146        catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
    138147      }
    139148    }
    … …  
    144153}
    145154
    146 int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) {
     155int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
     156
     157  off_t k;
     158
     159  int Nsecfilt = GetPhotcodeNsecfilt ();
     160
     161  off_t m = average[0].measureOffset;
     162
     163  for (k = 0; k < average[0].Nmeasure; k++, m++) {
     164    fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
     165  }
     166
     167  int Ns;
     168  for (Ns = 0; Ns < Nsecfilt; Ns++) {
     169    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
     170  }
     171  return 1;
     172}
     173
     174# define MARK_SKIP_MEAS \
     175  catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
     176  catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     177
     178// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
     179// * setMrel uses the internal Tiny structures only
     180// * setMrelOutput skips stars for which there are too few good measurements
     181// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
     182
     183// setMave is also similar to the above.  but:
     184// * setMave is called by relphot -update-objects
     185// * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2
     186// * setMave updates 2MASS average flags
     187// * setMave updates average EXT flags (PS1 and 2MASS)
     188int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
    147189
    148190  int i, N;
    149191  off_t j, k, m, Nmax;
    150   float Msys, Mcal, Mmos, Mgrid;
    151   double *list, *dlist;
    152   StatType stats;
    153   int Nsec, Nsecfilt, ecode;
    154 
    155   Nsecfilt = GetPhotcodeNsecfilt ();
     192  float Msys, Mcal, Mmos, Mgrid, Map;
     193  double *list, *dlist, *aplist, *daplist;
     194  StatType stats, apstats;
     195  int Nsec;
     196
     197  int Nsecfilt = GetPhotcodeNsecfilt ();
    156198
    157199  /* Nmeasure is now different, need to reallocate */
    … …  
    165207  ALLOCATE (dlist, double, MAX (1, Nmax));
    166208
    167   for (i = 0; i < Ncatalog; i++) {
    168     for (j = 0; j < catalog[i].Naverage; j++) {
    169       /* skip stars already calibrated */
    170       if (catalog[i].found[j]) continue; 
     209  ALLOCATE (aplist, double, MAX (1, Nmax));
     210  ALLOCATE (daplist, double, MAX (1, Nmax));
     211
     212  for (i = 0; i < Ncatalog; i++) {
     213    for (j = 0; j < catalog[i].Naverage; j++) {
     214
     215      if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) {
     216        fprintf (stderr, "test obj\n");
     217        print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure);
     218      }
     219
     220      int minUbercalDist = 1000;
    171221
    172222      int Ns;
    173223      for (Ns = 0; Ns < Nphotcodes; Ns++) {
     224
    174225        int thisCode = photcodes[Ns][0].code;
    175226        Nsec = GetPhotcodeNsec(thisCode);
     227
     228        /* star/photcodes already calibrated */
     229        if (catalog[i].found[Nsecfilt*j+Nsec]) continue; 
     230
     231        int Ncode = 0;
     232        int Next = 0;
     233        int haveSynth = FALSE;
     234        int haveUbercal = FALSE;
     235        int isRefPhot = FALSE;
    176236
    177237        N = 0;
    178238        m = catalog[i].averageT[j].measureOffset;
    179239        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
     240
    180241          // skip measurements that do not match the current photcode
    181           ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    182           if (ecode != thisCode) { continue; }
     242          PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
     243          if (!code) continue;
     244          if (code->equiv != thisCode) { continue; }
     245          Ncode ++;
    183246
    184247          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    185248
    186           // XXX allow REF stars (no Image Entry) to be included in the calculation this
    187           // should be optionally set, and should allow for REF stars to be downweighted by
    188           // more than their reported errors.  how such info is carried is unclear...
    189249          if (getImageEntry (m, i) < 0) {
     250            // these detetions have no image (eg, ref values such as 2MASS)
    190251            Mcal = Mmos = Mgrid = 0;
     252            isRefPhot = TRUE;
    191253          } else {
    192             Mcal  = getMcal  (m, i);
    193             if (isnan(Mcal)) continue;
     254            Mcal  = getMcal  (m, i, flatcorr, catalog);
     255            if (isnan(Mcal)) {
     256              MARK_SKIP_MEAS;
     257              continue;
     258            }
    194259            Mmos  = getMmos  (m, i);
    195             if (isnan(Mmos)) continue;
     260            if (isnan(Mmos)) {
     261              MARK_SKIP_MEAS;
     262              continue;
     263            }
    196264            Mgrid = getMgrid (m, i);
    197             if (isnan(Mgrid)) continue;
     265            if (isnan(Mgrid)) {
     266              MARK_SKIP_MEAS;
     267              continue;
     268            }
    198269          }
    199270
    200271          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
     272          if (isnan(Msys)) {
     273            if (!isRefPhot) {
     274              MARK_SKIP_MEAS;
     275            }
     276            continue;
     277          }
    201278          list[N] = Msys - Mcal - Mmos - Mgrid;
    202           dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
     279
     280          Map = PhotAper (&catalog[i].measure[m]);
     281          aplist[N] = Map - Mcal - Mmos - Mgrid;
     282
     283          int myUbercalDist = getUbercalDist(m,i);
     284          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
     285
     286          // count the extended detections
     287          if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
     288            if (!isnan(catalog[i].measure[m].Map)) {
     289              if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
     290                Next ++;
     291              }
     292            }
     293          }
     294          // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
     295          if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) {
     296            if (catalog[i].measure[m].photFlags & 0x00c00000) {
     297              Next ++;
     298            }
     299            if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) {
     300              // detections without one of these bits should only be used in PASS_1
     301              MARK_SKIP_MEAS;
     302              continue;
     303            }
     304          }
     305          // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
     306          if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) {
     307            if (pass < 4) {
     308              MARK_SKIP_MEAS;
     309              continue;
     310            }
     311            haveSynth = TRUE;
     312          }
     313
     314          // dlist gives the error, which is used as the weight in WT_MEAN.
     315          // we can modify the resulting weight in a few ways:
     316          // 1) MIN_ERROR guarantees a floor
     317          // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
     318          // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
     319          // 4) some reference photcode of some kind can be specified as fixed and have a high weight
     320          dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
     321
     322          // up-weight the ubercal values (or convergence can take a long time...) (XXX make this optional?)
     323          if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     324            dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
     325            haveUbercal = TRUE;
     326          }
     327
     328          // tie down reference photometry if the -refcode (code) option is selected
     329          // eg, -refcode g_SDSS
     330          // this probably makes no sense in the context of multifilter analysis
     331          if (refPhotcode) {
     332            if (code->code == refPhotcode->code) {
     333              // tiny error -> large weight
     334              // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
     335              dlist[N] = 0.0001;
     336            }
     337          }
    203338          N++;
    204339        }
    205340        if (N < 1) continue;
    206341
     342        for (m = 0; m < N; m++) {
     343          daplist[m] = dlist[m];
     344        }
     345
     346        // XXX force WT_MEAN or MEAN here?
    207347        liststats (list, dlist, N, &stats);
    208         if (mark) catalog[i].found[j] = TRUE;
     348        liststats (aplist, daplist, N, &apstats);
     349        catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
     350        switch (pass) {
     351          case 0:
     352            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
     353            break;
     354          case 1:
     355            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
     356            break;
     357          case 2:
     358            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
     359            break;
     360          case 3:
     361            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
     362            break;
     363          case 4:
     364            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
     365            break;
     366        }
     367
     368        if (haveSynth) {
     369          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
     370        }       
     371        if (haveUbercal) {
     372          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
     373        }       
    209374
    210375        /* use sigma or error in dM for output? */
    211376        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
     377        catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean;
    212378        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
    213379        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
     380        catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
     381        catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
     382
     383        catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
     384        catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
     385        catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
     386
     387        if ((Next > 0) && (Next > 0.5*N)) {
     388          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
     389        }
    214390      }
    215391    }
    … …  
    218394  free (list);
    219395  free (dlist);
     396
     397  free (aplist);
     398  free (daplist);
    220399  return (TRUE);
    221400}
     401
     402/* grab Nsec for named photcode */
     403# define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME)      \
     404  short MY_NSEC = -1;                           \
     405  {                                             \
     406  PhotCode *code = GetPhotcodebyName (NAME);    \
     407  if (code) {                                   \
     408    MY_NSEC = GetPhotcodeNsec (code->equiv);    \
     409  } }
    222410
    223411// For each average object, set the average mags based on existing equiv photometry.
    224412// NOTE: this function operates on the real Measure & Average structures, not the
    225413// MeasureTiny & AverageTiny structures
     414// NOTE: this function is called on the remote machine -- make sure any (global) options
     415// are passed to the relphot_client program
    226416int setMave (Catalog *catalog, int Ncatalog) {
    227417
    … …  
    231421  double *list, *dlist;
    232422  StatType stats;
    233   PhotCode *code;
    234423  DVOAverageFlags flagBits;
    235424
    … …  
    248437  Nsecfilt = GetPhotcodeNsecfilt ();
    249438
    250 # define PSFQUALSTATS 1
     439  // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data
     440  NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J");
     441  NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H");
     442  NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K");
     443
    251444  for (i = 0; i < Ncatalog; i++) {
    252445    for (j = 0; j < catalog[i].Naverage; j++) {
    253446
    254447      // update average photometry for each of the average filters
     448
     449      // XXX Note that this would be faster if we had an array of results and accumulated
     450      // them in a single pass
     451
     452      int Next = 0;
     453
    255454      for (Ns = 0; Ns < Nsecfilt; Ns++) {
    256455
    257         code = GetPhotcodebyNsec (Ns);
     456        PhotCode *code = GetPhotcodebyNsec (Ns);
    258457        Nc = code[0].code;
    259458       
    … …  
    267466          if (isnan(Msys)) continue;
    268467
     468          // reject POOR detections (PSF_QF < 0.85) and count extended detections
    269469          // XXX only apply this filter for psphot data from GPC1 for now...
    270           if (PSFQUALSTATS && (catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
    271               if (catalog[i].measure[m].psfQual < 0.85) continue;
     470          if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
     471            if (catalog[i].measure[m].psfQual < 0.85) continue;
     472            if (!isnan(catalog[i].measure[m].Map)) {
     473              if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
     474                Next ++;
     475              }
     476            }
    272477          }
    273478
    … …  
    287492        catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
    288493        catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
     494
     495        catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80;
     496        catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20;
     497
     498        if ((Next > 0) && (Next > 0.5*N)) {
     499          catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT;
     500        }
    289501      }
    290502
    291503      // update average flags based on the detection stats. 
    292       // XXX we need to clean this up -- this is a serious set of hacks...
    293       if (PSFQUALSTATS) {
    294         int Galaxy2MASS, GalaxySDSS, goodPS1, nEXT, nPSF, good2MASS;
    295 
    296         Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
    297         GalaxySDSS = FALSE;  // best guess for galaxy based on SDSS measurements (XXX need to fix SDSS flags)
    298         goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
    299         good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
    300         nEXT = nPSF = 0;     // number of PS1 PSF vs EXT measurements
    301 
    302         m = catalog[i].average[j].measureOffset;
    303         for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    304 
    305           // PS1 data :
    306           if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
    307             if (catalog[i].measure[m].psfQual > 0.85) {
    308               goodPS1 = TRUE;
    309               if (!isnan(catalog[i].measure[m].Map)) {
    310                 if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
    311                   nEXT ++;
    312                 } else {
    313                   nPSF ++;
    314                 }
     504
     505      int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
     506      int goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
     507      int good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
     508      int nEXT = 0;
     509      int nPSF = 0;     // number of PS1 PSF vs EXT measurements
     510      int have2MASS = FALSE;
     511     
     512      // count, flag good and extended detections
     513      m = catalog[i].average[j].measureOffset;
     514      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
     515
     516        // PS1 data :
     517        if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
     518          if (catalog[i].measure[m].psfQual > 0.85) {
     519            goodPS1 = TRUE;
     520            if (!isnan(catalog[i].measure[m].Map)) {
     521              if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
     522                nEXT ++;
     523              } else {
     524                nPSF ++;
    315525              }
    316526            }
    317527          }
     528        }
    318529         
    319           // 2MASS data:
    320           if (catalog[i].measure[m].photcode == 2011) {
    321             if (catalog[i].measure[m].photFlags & 0x00c00000) {
    322               Galaxy2MASS = TRUE;
    323             }
    324             if (catalog[i].measure[m].photFlags & 0x00000007) {
    325               good2MASS = TRUE;
    326             }
    327           } 
    328         }
    329 
    330         // we attempt to set a few flags here; reset those bits before trying:
    331         catalog[i].average[j].flags &= ~flagBits;
    332 
    333         if (nEXT && (nEXT > nPSF)) {
    334           catalog[i].average[j].flags |= ID_OBJ_EXT;
    335         }
    336         if (goodPS1) {
    337           catalog[i].average[j].flags |= ID_OBJ_GOOD;
    338         }
    339         if (Galaxy2MASS) {
    340           catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
    341         }
    342         if (good2MASS) {
    343           catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
    344         }
     530        // 2MASS data J-band flags
     531        if (catalog[i].measure[m].photcode == 2011) {
     532          // only need to do this once (always have JHK triplet; galaxy flag is same for all 3)
     533          have2MASS = TRUE;
     534          if (catalog[i].measure[m].photFlags & 0x00c00000) {
     535            Galaxy2MASS = TRUE;
     536          }
     537          if (catalog[i].measure[m].photFlags & 0x00000007) {
     538            good2MASS = TRUE;
     539            if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0;
     540          } else {
     541            if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1;
     542          }
     543        } 
     544        // 2MASS data H-band flags
     545        if (catalog[i].measure[m].photcode == 2012) {
     546          if (catalog[i].measure[m].photFlags & 0x00000007) {
     547            good2MASS = TRUE;
     548            if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0;
     549          } else {
     550            if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1;
     551          }
     552        } 
     553        // 2MASS data K-band flags
     554        if (catalog[i].measure[m].photcode == 2013) {
     555          if (catalog[i].measure[m].photFlags & 0x00000007) {
     556            good2MASS = TRUE;
     557            if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0;
     558          } else {
     559            if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1;
     560          }
     561        } 
     562      }
     563
     564      // we attempt to set a few flags here; reset those bits before trying:
     565      catalog[i].average[j].flags &= ~flagBits;
     566
     567      // XXX set the secfilt bits?
     568      if (nEXT && (nEXT > nPSF)) {
     569        catalog[i].average[j].flags |= ID_OBJ_EXT;
     570      }
     571      if (goodPS1) {
     572        catalog[i].average[j].flags |= ID_OBJ_GOOD;
     573      }
     574      if (Galaxy2MASS) {
     575        catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
     576        if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT;
     577        if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT;
     578        if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT;
     579      }
     580      if (good2MASS) {
     581        catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
    345582      }
    346583    }
    … …  
    353590
    354591/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
    355 int setMcalOutput (Catalog *catalog, int Ncatalog) {
     592int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    356593
    357594  int i;
    … …  
    368605      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    369606        if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    370         Mcal  = getMcal  (m, i);
     607        Mcal  = getMcal  (m, i, flatcorr, catalog);
    371608        if (isnan(Mcal)) continue;
    372609        Mmos  = getMmos  (m, i);
    … …  
    377614        // set the output calibration
    378615        catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
     616
     617        if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
     618          myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag");
     619        }
    379620      }
    380621    }
    … …  
    385626void clean_stars (Catalog *catalog, int Ncatalog) {
    386627
    387   int i, j, Ndel, Nave, Ntot, mark, Ns;
     628  int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
    388629  float dM, Xm;
    389630  double Chisq, MaxScatter, MaxChisq;
    … …  
    402643
    403644  int Nsecfilt = GetPhotcodeNsecfilt ();
     645
     646  // XX int oldPLOTSTUFF = PLOTSTUFF;
     647  // XX PLOTSTUFF = TRUE;
     648  // XX plot_chisq (catalog, Ncatalog);
     649  // XX PLOTSTUFF = oldPLOTSTUFF;
    404650
    405651  // eliminate bad stars using the stats for a single secfilt at a time
    … …  
    430676    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
    431677
    432     Ndel = Nave = 0;
     678    Ndel = Nave = Nscat = Nnan = Nchi = 0;
    433679    for (i = 0; i < Ncatalog; i++) {
    434680      for (j = 0; j < catalog[i].Naverage; j++) {
    … …  
    440686          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
    441687          Ndel ++;
     688          if (dM > MaxScatter) { Nscat ++; }
     689          if (Xm == NAN_S_SHORT) { Nnan ++; }
     690          if (Chisq > MaxChisq) { Nchi ++; }
    442691        } else {
    443692          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
    … …  
    446695      }
    447696    }
    448     fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
     697    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
    449698    initstats (STATMODE);
    450699  }
    … …  
    456705# define NSIGMA_CLIP 3.0
    457706# define NSIGMA_REJECT 5.0
    458 void clean_measures (Catalog *catalog, int Ncatalog, int final) {
     707void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) {
    459708
    460709  off_t j, k, m, Nmax, Ndel, Nave;
    … …  
    492741      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    493742       
    494         /* on final processing, skip stars already measured */
    495         if (final && catalog[i].found[j]) continue; 
    496 
    497743        int thisCode = photcodes[Ns][0].code;
    498744        int Nsec = GetPhotcodeNsec(thisCode);
    499745       
     746        /* on final processing, skip stars already measured */
     747        if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue; 
     748
    500749        /* skip bad stars to prevent them from becoming good (on inner sample) */
    501750        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    … …  
    509758          if (ecode != thisCode) { continue; }
    510759
    511           /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
    512           Mcal  = getMcal  (m, i);
     760          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
     761          Mcal  = getMcal  (m, i, flatcorr, catalog);
    513762          if (isnan(Mcal)) { Ncal ++; continue; }
    514763          Mmos  = getMmos  (m, i);
    … …  
    522771          N++;
    523772        }
    524         if (N <= TOOFEW) { Nfew ++; continue; }
     773        if (N <= TOOFEW) {
     774          Nfew ++;
     775          continue;
     776        }
    525777
    526778        /* 3-sigma clip based on stats of inner 50% */
    … …  
    553805          if (ecode != thisCode) { continue; }
    554806
    555           /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
    556           Mcal  = getMcal  (m, i);
     807          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
     808          Mcal  = getMcal  (m, i, flatcorr, catalog);
    557809          if (isnan(Mcal)) continue;
    558810          Mmos  = getMmos  (m, i);
    … …  
    572824        /* mark bad measures (> 3 sigma deviant) */
    573825        for (k = 0; k < N; k++) {
    574           if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
     826          // treat the scatter of the star as a systematic term.  this is a bit of an
     827          // over-estimage (a perfect Gauss distribution with perfect errors would have
     828          // mySigma = sqrt(2) too large)
     829          float mySigma = hypot (stats.sigma, dlist[k]);
     830          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
    575831            catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
    576832            if (final) {
    … …  
    591847  free (ilist);
    592848  free (tlist);
    593 }
    594 
    595 StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) {
     849
     850
     851}
     852
     853StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) {
    596854
    597855  off_t j, k, m, Ntot;
    … …  
    625883        int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    626884        if (ecode != seccode) { continue;}
    627         Mcal = getMcal  (m, i);
     885        Mcal = getMcal  (m, i, flatcorr, catalog);
    628886        if (isnan(Mcal)) { continue;}
    629887        Mmos = getMmos  (m, i);
    … …  
    712970
    713971      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
     972      if (isnan(dM)) continue;
    714973      list[n] = dM;
    715974      dlist[n] = 1;

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

    r31668 r33651  
    3838  UserPatch.Dmin = -90;
    3939  UserPatch.Dmax = +90;
    40   UserPatchSelect = FALSE;
    4140  if ((N = get_argument (argc, argv, "-region"))) {
    4241    remove_argument (N, &argc, argv);
    … …  
    4948    UserPatch.Dmax = atof (argv[N]);
    5049    remove_argument (N, &argc, argv);
    51     UserPatchSelect = TRUE;
     50  }
     51
     52  /* specify region file by name (eg n0000/0000.00) */
     53  UserCatalog = NULL;
     54  if ((N = get_argument (argc, argv, "-catalog"))) {
     55    remove_argument (N, &argc, argv);
     56    UserCatalog = strcreate (argv[N]);
     57    remove_argument (N, &argc, argv);
    5258  }
    5359
    … …  
    6874  }
    6975
     76  // XXX for the moment, make this selection manual.  it needs to be automatic
     77  // based on the state of the SkyTable
     78  PARALLEL = FALSE;
     79  if ((N = get_argument (argc, argv, "-parallel"))) {
     80    PARALLEL = TRUE;
     81    remove_argument (N, &argc, argv);
     82  }
     83  // this is a test mode : rather than launching the remote jobs and waiting for completion,
     84  // relphot will simply list the remote command and wait for the user to signal completion
     85  PARALLEL_MANUAL = FALSE;
     86  if ((N = get_argument (argc, argv, "-parallel-manual"))) {
     87    PARALLEL = TRUE; // -parallel-manual implies -parallel
     88    PARALLEL_MANUAL = TRUE;
     89    remove_argument (N, &argc, argv);
     90  }
     91  // this is a test mode : rather than launching the relphot_client jobs remotely, they are
     92  // run in serial via 'system'
     93  PARALLEL_SERIAL = FALSE;
     94  if ((N = get_argument (argc, argv, "-parallel-serial"))) {
     95    if (PARALLEL_MANUAL) {
     96      fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n");
     97      exit (1);
     98    }
     99    PARALLEL = TRUE; // -parallel-serial implies -parallel
     100    PARALLEL_SERIAL = TRUE;
     101    remove_argument (N, &argc, argv);
     102  }
     103
    70104  PLOTSTUFF = FALSE;
    71105  if ((N = get_argument (argc, argv, "-plot"))) {
    … …  
    90124  }     
    91125
    92   strcpy (STATMODE, "CHI_INNER_WTMEAN");
     126  strcpy (STATMODE, "WT_MEAN");
    93127  if ((N = get_argument (argc, argv, "-statmode"))) {
    94128    remove_argument (N, &argc, argv);
    … …  
    128162  }
    129163
     164  UPDATE_CATFORMAT = NULL;
     165  if ((N = get_argument (argc, argv, "-update-catformat"))) {
     166    remove_argument (N, &argc, argv);
     167    UPDATE_CATFORMAT = strcreate (argv[N]);
     168    remove_argument (N, &argc, argv);
     169  }
     170
    130171  SAVE_IMAGE_UPDATES = TRUE;
    131172  if ((N = get_argument (argc, argv, "-skip-image-updates"))) {
    … …  
    140181    remove_argument (N, &argc, argv);
    141182    MaxDensityUse = TRUE;
     183  }
     184
     185  CLOUD_TOLERANCE = 0.02;
     186  if ((N = get_argument (argc, argv, "-cloud-limit"))) {
     187    remove_argument (N, &argc, argv);
     188    CLOUD_TOLERANCE = atof(argv[N]);
     189    remove_argument (N, &argc, argv);
    142190  }
    143191
    … …  
    189237  }
    190238
     239  KEEP_UBERCAL = TRUE;
     240  if ((N = get_argument (argc, argv, "-reset-ubercal"))) {
     241    remove_argument (N, &argc, argv);
     242    KEEP_UBERCAL = FALSE;
     243  }
     244
    191245  MIN_ERROR = 0.001;
    192246  if ((N = get_argument (argc, argv, "-minerror"))) {
    … …  
    247301
    248302  if (UpdateAverages && (argc == 1)) return TRUE;
    249   if (UserPatchSelect && (argc == 2)) return TRUE;
    250   if (argc != 3) relphot_usage ();
     303  if (argc != 2) relphot_usage ();
    251304
    252305  return TRUE;
    253306}
     307
     308int args_client (int argc, char **argv) {
     309
     310  int N;
     311  double trange;
     312
     313  // by definition, the client is not parallel
     314  PARALLEL = FALSE;
     315  PARALLEL_MANUAL = FALSE;
     316  PARALLEL_SERIAL = FALSE;
     317
     318  HOST_ID = 0;
     319  if ((N = get_argument (argc, argv, "-hostID"))) {
     320    remove_argument (N, &argc, argv);
     321    HOST_ID = atoi (argv[N]);
     322    remove_argument (N, &argc, argv);
     323  }
     324  if (!HOST_ID) relphot_client_usage();
     325
     326  HOSTDIR = NULL;
     327  if ((N = get_argument (argc, argv, "-hostdir"))) {
     328    remove_argument (N, &argc, argv);
     329    HOSTDIR = strcreate (argv[N]);
     330    remove_argument (N, &argc, argv);
     331  }
     332  if (!HOSTDIR) relphot_client_usage();
     333
     334  // CATDIR[0] = 0;
     335  // if ((N = get_argument (argc, argv, "-catdir"))) {
     336  //   remove_argument (N, &argc, argv);
     337  //   strcpy (CATDIR, argv[N]);
     338  //   remove_argument (N, &argc, argv);
     339  // }
     340  // if (!CATDIR[0]) relphot_client_usage();
     341
     342  IMAGES = NULL; // used in -update mode
     343  BCATALOG = NULL; // used in -load mode
     344  MODE = MODE_NONE;
     345  if ((N = get_argument (argc, argv, "-load"))) {
     346    MODE = MODE_LOAD;
     347    remove_argument (N, &argc, argv);
     348    BCATALOG = strcreate (argv[N]);
     349    remove_argument (N, &argc, argv);
     350  }
     351  if ((N = get_argument (argc, argv, "-update-catalogs"))) {
     352    if (MODE) {
     353      fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n");
     354      relphot_client_usage();
     355    }
     356    MODE = MODE_UPDATE;
     357    remove_argument (N, &argc, argv);
     358    IMAGES = strcreate (argv[N]);
     359    remove_argument (N, &argc, argv);
     360  }
     361  if ((N = get_argument (argc, argv, "-update-objects"))) {
     362    if (MODE) {
     363      fprintf (stderr, "ERROR: cannot mix modes (-load, -update-catalogs, -update-objects)\n");
     364      relphot_client_usage();
     365    }
     366    MODE = MODE_UPDATE_OBJECTS;
     367    remove_argument (N, &argc, argv);
     368  }
     369  if (!MODE) relphot_client_usage();
     370
     371  strcpy (STATMODE, "WT_MEAN");
     372  if ((N = get_argument (argc, argv, "-statmode"))) {
     373    remove_argument (N, &argc, argv);
     374    strcpy (STATMODE, argv[N]);
     375    remove_argument (N, &argc, argv);
     376  }
     377
     378  /* specify portion of the sky */
     379  UserPatch.Rmin = 0;
     380  UserPatch.Rmax = 360;
     381  UserPatch.Dmin = -90;
     382  UserPatch.Dmax = +90;
     383  if ((N = get_argument (argc, argv, "-region"))) {
     384    remove_argument (N, &argc, argv);
     385    UserPatch.Rmin = atof (argv[N]);
     386    remove_argument (N, &argc, argv);
     387    UserPatch.Rmax = atof (argv[N]);
     388    remove_argument (N, &argc, argv);
     389    UserPatch.Dmin = atof (argv[N]);
     390    remove_argument (N, &argc, argv);
     391    UserPatch.Dmax = atof (argv[N]);
     392    remove_argument (N, &argc, argv);
     393  }
     394
     395  VERBOSE = VERBOSE2 = FALSE;
     396  if ((N = get_argument (argc, argv, "-v"))) {
     397    VERBOSE = TRUE;
     398    remove_argument (N, &argc, argv);
     399  }
     400  if ((N = get_argument (argc, argv, "-vv"))) {
     401    VERBOSE2 = VERBOSE = TRUE;
     402    remove_argument (N, &argc, argv);
     403  }
     404
     405  RESET = FALSE;
     406  if ((N = get_argument (argc, argv, "-reset"))) {
     407    remove_argument (N, &argc, argv);
     408    RESET = TRUE;
     409  }
     410
     411  KEEP_UBERCAL = TRUE;
     412  if ((N = get_argument (argc, argv, "-reset-ubercal"))) {
     413    remove_argument (N, &argc, argv);
     414    KEEP_UBERCAL = FALSE;
     415  }
     416
     417  /* define time */
     418  TimeSelect = FALSE;
     419  if ((N = get_argument (argc, argv, "-time"))) {
     420    remove_argument (N, &argc, argv);
     421    if (!ohana_str_to_time (argv[N], &TSTART)) {
     422      fprintf (stderr, "ERROR: syntax error\n");
     423      return (FALSE);
     424    }
     425    remove_argument (N, &argc, argv);
     426    if (!ohana_str_to_dtime (argv[N], &trange)) {
     427      if (!ohana_str_to_time (argv[N], &TSTOP)) {
     428        fprintf (stderr, "ERROR: syntax error\n");
     429        return (FALSE);
     430      }
     431    } else {
     432      if (trange < 0) {
     433        trange = fabs (trange);
     434        TSTOP = TSTART;
     435        TSTART -= trange;
     436      } else {
     437        TSTOP = TSTART + trange;
     438      }
     439    }
     440    remove_argument (N, &argc, argv);
     441    TimeSelect = TRUE;
     442  }
     443
     444
     445  MIN_ERROR = 0.001;
     446  if ((N = get_argument (argc, argv, "-minerror"))) {
     447    remove_argument (N, &argc, argv);
     448    MIN_ERROR = atof (argv[N]);
     449    remove_argument (N, &argc, argv);
     450    /* require MIN_ERROR > 0 */
     451  } 
     452
     453  UPDATE = FALSE;
     454  if ((N = get_argument (argc, argv, "-update"))) {
     455    remove_argument (N, &argc, argv);
     456    UPDATE = TRUE;
     457  }
     458
     459  UPDATE_CATFORMAT = NULL;
     460  if ((N = get_argument (argc, argv, "-update-catformat"))) {
     461    remove_argument (N, &argc, argv);
     462    UPDATE_CATFORMAT = strcreate (argv[N]);
     463    remove_argument (N, &argc, argv);
     464  }
     465
     466  AreaSelect = FALSE;
     467  if ((N = get_argument (argc, argv, "-area"))) {
     468    remove_argument (N, &argc, argv);
     469    AreaXmin = atof (argv[N]);
     470    remove_argument (N, &argc, argv);
     471    AreaXmax = atof (argv[N]);
     472    remove_argument (N, &argc, argv);
     473    AreaYmin = atof (argv[N]);
     474    remove_argument (N, &argc, argv);
     475    AreaYmax = atof (argv[N]);
     476    remove_argument (N, &argc, argv);
     477    AreaSelect = TRUE;
     478  }
     479
     480  ImagSelect = FALSE;
     481  if ((N = get_argument (argc, argv, "-instmag"))) {
     482    remove_argument (N, &argc, argv);
     483    ImagMin = atof (argv[N]);
     484    remove_argument (N, &argc, argv);
     485    ImagMax = atof (argv[N]);
     486    remove_argument (N, &argc, argv);
     487    ImagSelect = TRUE;
     488  }
     489
     490  DophotSelect = FALSE;
     491  if ((N = get_argument (argc, argv, "-dophot"))) {
     492    remove_argument (N, &argc, argv);
     493    DophotValue = atof (argv[N]);
     494    remove_argument (N, &argc, argv);
     495    DophotSelect = TRUE;
     496  }
     497
     498  MaxDensityUse = FALSE;
     499  if ((N = get_argument (argc, argv, "-max-density"))) {
     500    remove_argument (N, &argc, argv);
     501    MaxDensityValue = atof(argv[N]);
     502    remove_argument (N, &argc, argv);
     503    MaxDensityUse = TRUE;
     504  }
     505
     506  if ((MODE == MODE_UPDATE_OBJECTS)  && (argc == 1)) return TRUE;
     507  if (argc != 2) relphot_client_usage ();
     508
     509  return TRUE;
     510}
     511

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

    r31668 r33651  
    11# include "relphot.h"
    22
    3 extern double drand48();
    4 
    5 int CopyAverageTiny (AverageTiny *averageT, Average *average) {
    6 
    7   averageT[0].R             = average[0].R;
    8   averageT[0].D             = average[0].D;
    9   averageT[0].flags         = average[0].flags;
    10   averageT[0].Nmeasure      = average[0].Nmeasure;
    11   averageT[0].measureOffset = average[0].measureOffset;
    12 
    13   // make Nmeasure & measureOffset optional?
    14 
    15   return (TRUE);
    16 }
    17 
    18 int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) {
    19 
    20   measureT[0].dR       = measure[0].dR;
    21   measureT[0].dD       = measure[0].dD;
    22   measureT[0].M        = measure[0].M;
    23   measureT[0].Mcal     = measure[0].Mcal;
    24   measureT[0].dM       = measure[0].dM;
    25   measureT[0].airmass  = measure[0].airmass;
    26   measureT[0].Xccd     = measure[0].Xccd;
    27   measureT[0].Yccd     = measure[0].Yccd;
    28   measureT[0].t        = measure[0].t;
    29   measureT[0].dt       = measure[0].dt;
    30   measureT[0].averef   = measure[0].averef;
    31   measureT[0].imageID  = measure[0].imageID;
    32   measureT[0].dbFlags  = measure[0].dbFlags;
    33   measureT[0].photcode = measure[0].photcode;
    34 
    35   return (TRUE);
    36 }
     3int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog);
    374
    385int bcatalog (Catalog *subcatalog, Catalog *catalog) {
    396 
    407  off_t i, j, offset;
    41   int ecode, found, Ns;
     8  int found, Ns, *Nvalid;
    429  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
    4310  float mag;
    44   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf;
     11  int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf, Nnan, Nbad, Npoor;
    4512
    4613  int Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    5320  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
    5421  ALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE);
    55   ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     22  ALLOCATE (subcatalog[0].secfilt,  SecFilt,    NAVERAGE*Nsecfilt);
    5623  Nmeasure = Naverage = 0;
    5724
    58   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = 0;
     25  Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
     26
     27  // copy the following fields to the subcatalog:
     28  subcatalog[0].catID = catalog[0].catID;
     29  subcatalog[0].filename = catalog[0].filename;
     30
     31  // array to count measures with photcodes->equiv to each of the secfilts
     32  ALLOCATE (Nvalid, int, Nsecfilt);
     33
     34  // flags used by the photometry analysis (excluding UBERCAL)
     35  unsigned int PHOTOM_FLAGS =
     36    ID_MEAS_NOCAL |       // detection ignored for this analysis (photcode, time range)
     37    ID_MEAS_POOR_PHOTOM | // detection is photometry outlier
     38    ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement
     39    ID_MEAS_AREA;         // detetion was outside acceptable area of device
    5940
    6041  /* exclude stars not in range or with too few measurements */
    … …  
    6344
    6445    /* start with all stars good */
    65     CopyAverageTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]);
     46    CopyAverageToTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]);
    6647    subcatalog[0].averageT[Naverage].measureOffset = Nmeasure;
    6748
    6849    for (j = 0; j < Nsecfilt; j++) {
    6950      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
    70     }
    71 
     51      Nvalid[j] = 0; // reset the Nvalid array for this star
     52    }
     53
     54    // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images)
    7255    if (RESET) {
    7356      int Ns;
     57      DVOSecfiltFlags secfiltBits =
     58        ID_SECF_STAR_FEW |
     59        ID_SECF_STAR_POOR |
     60        ID_PHOTOM_PASS_0 |
     61        ID_PHOTOM_PASS_1 |
     62        ID_PHOTOM_PASS_2 |
     63        ID_PHOTOM_PASS_3 |
     64        ID_PHOTOM_PASS_4 |
     65        ID_SECF_USE_SYNTH |
     66        ID_SECF_USE_UBERCAL |
     67        ID_SECF_OBJ_EXT;
    7468      for (Ns = 0; Ns < Nphotcodes; Ns++) {
    7569
    … …  
    7771        int Nsec = GetPhotcodeNsec(thisCode);
    7872
    79         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M  = NAN;
    80         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN;
    81         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_FEW;
    82         subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_POOR;
     73        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
     74        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
     75        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
     76        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm     = NAN;
     77        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
     78        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
     79        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
     80        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
     81        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
     82        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
    8383      }
    8484    }
    … …  
    9292
    9393      /* select measurements by photcode */
    94       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode);
     94      PhotCode *code = GetPhotcodebyCode (catalog[0].measure[offset].photcode);
     95      if (!code) continue;
    9596      found = FALSE;
     97      int Nsec = -1; // Nsec equivalent for the measurement
    9698      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
    97         if (ecode == photcodes[Ns][0].code) found = TRUE;
     99        if (code->equiv != photcodes[Ns][0].code) continue;
     100        found = TRUE;
     101        Nsec = GetPhotcodeNsec(code->equiv);
    98102      }
    99103      if (!found) {
    … …  
    113117      // skip garbage measurements
    114118      if (catalog[0].measure[offset].psfQual < 0.85) { Npsfqf ++; continue; }
    115 
    116       // check for galaxies
    117       if (!isnan(catalog[0].measure[offset].Map)) {
     119      if (isnan(catalog[0].measure[offset].M)) { Nnan ++; continue; }
     120
     121      if (catalog[0].measure[offset].photFlags & code->photomBadMask) {
     122        Nbad++;
     123        continue;
     124      }
     125      if (catalog[0].measure[offset].photFlags & code->photomPoorMask) { Npoor++; continue;}
     126
     127      // check for galaxies (XXX skip for now)
     128      if (FALSE && !isnan(catalog[0].measure[offset].Map)) {
    118129        if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
    119130          nEXT ++;
    … …  
    137148      }
    138149
    139       CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
     150      // count this measurement as valid for this secfilt entry
     151      if (Nsec > -1) {
     152        assert (Nsec < Nsecfilt);
     153        Nvalid[Nsec] ++;
     154      }
     155
     156      CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
    140157      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    141158      subcatalog[0].measureT[Nmeasure].averef = Naverage;
    142159      if (RESET) {
    143         subcatalog[0].measureT[Nmeasure].Mcal = 0;
    144         subcatalog[0].measureT[Nmeasure].dbFlags &= 0xff00;
    145         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    146         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
    147         subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
     160        if (!KEEP_UBERCAL) {
     161          subcatalog[0].measureT[Nmeasure].Mcal = 0;
     162          subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     163        }
     164        subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS;
    148165      }
    149166      Nmeasure ++;
    … …  
    156173
    157174    // skip object if it is likely to be a galaxy
    158     if (nEXT >= nPSF) {
     175    if (FALSE && (nEXT >= nPSF)) {
    159176      Nmeasure -= Nm;
    160177      Ngalaxy ++;
    … …  
    163180
    164181    // XXXX test : what checks do I need to make elsewhere to avoid problems here?
    165     if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */
     182    if (Nm <= STAR_TOOFEW) { /* enough measurements total? */
    166183      Nmeasure -= Nm;
    167184      Nfew ++;
    168185      continue;
    169186    }
     187    int anySecfiltGood = FALSE;
     188    for (Ns = 0; Ns < Nsecfilt; Ns++) { /* enough measurements in at least one band? */
     189      if (Nvalid[Ns] <= STAR_TOOFEW) continue;
     190      anySecfiltGood = TRUE;
     191    }
     192    if (!anySecfiltGood) {
     193      Nmeasure -= Nm;
     194      Nfew ++;
     195      continue;
     196    }
     197
    170198    subcatalog[0].averageT[Naverage].Nmeasure = Nm;
    171199    Naverage ++;
    … …  
    186214  // limit the total number of stars in the catalog
    187215  if (MaxDensityUse) {
    188     LimitDensityCatalog (subcatalog, catalog);
     216    LimitDensityCatalog_ByNmeasure (subcatalog, catalog);
    189217  }
    190218
    … …  
    192220    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n",
    193221             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
    194     fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d galaxies\n",
    195              Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Ngalaxy);
     222    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n",
     223             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
    196224  }
    197225  return (TRUE);
    198226}
    199227
    200 int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
     228// sort by decreasing Nmeasure (X)
     229void sort_by_Nmeasure (int *X, off_t *Y, off_t N) {
     230
     231# define SWAPFUNC(A,B){ int tmpI; off_t tmpT;   \
     232  tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \
     233  tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \
     234}
     235# define COMPARE(A,B)(X[A] > X[B])
     236
     237  OHANA_SORT (N, COMPARE, SWAPFUNC);
     238
     239# undef SWAPFUNC
     240# undef COMPARE
     241
     242}
     243
     244int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) {
    201245
    202246  Catalog tmpcatalog;
    … …  
    206250  int Nsecfilt = GetPhotcodeNsecfilt ();
    207251
    208   gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
    209   gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
    210   gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
    211   gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
     252  gfits_scan (&oldcatalog[0].header, "RA0",  "%lf", 1, &Rmin);
     253  gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin);
     254  gfits_scan (&oldcatalog[0].header, "RA1",  "%lf", 1, &Rmax);
     255  gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax);
    212256
    213257  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
    … …  
    229273
    230274  // we are going to select Nmax entries by generating a random-sorted index list
    231   off_t *index, tmp, i, j, ave;
     275  int *value;
     276  off_t *index, i, j, ave;
    232277  ALLOCATE (index, off_t, Naverage);
     278  ALLOCATE (value, int, Naverage);
    233279  for (i = 0; i < Naverage; i++) {
    234280    index[i] = i;
    235   }
    236   for (i = 0; i < Naverage; i++) {
    237     j = (Naverage - i) * drand48() + i; // a number between i and Naverage
    238     tmp = index[j];
    239     index[j] = index[i];
    240     index[i] = tmp;
    241   }
     281    value[i] = subcatalog[0].averageT[i].Nmeasure;
     282  }
     283  sort_by_Nmeasure (value, index, Naverage);
    242284
    243285  // count the number of measurements this selection will yield
    … …  
    266308      Nmeasure ++;
    267309    }
     310    for (j = 0; j < Nsecfilt; j++) {
     311      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
     312    }
     313  }
     314
     315  if (VERBOSE) {
     316    char *basename = filebasename (oldcatalog[0].filename);
     317    fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n",
     318             Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  oldcatalog[0].Nmeasure, basename);
     319    free (basename);
     320  }
     321
     322  free (index);
     323  free (value);
     324  free (subcatalog[0].averageT);
     325  free (subcatalog[0].measureT);
     326  free (subcatalog[0].secfilt);
     327
     328  subcatalog[0].averageT = tmpcatalog.averageT;
     329  subcatalog[0].measureT = tmpcatalog.measureT;
     330  subcatalog[0].secfilt = tmpcatalog.secfilt;
     331  subcatalog[0].Naverage = Nmax;
     332  subcatalog[0].Nmeasure = Nmeasure;
     333  subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt;
     334  subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt;
     335
     336  return (TRUE);
     337}
     338
     339int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
     340
     341  Catalog tmpcatalog;
     342
     343  double Rmin, Rmax, Dmin, Dmax;
     344
     345  int Nsecfilt = GetPhotcodeNsecfilt ();
     346
     347  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
     348  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
     349  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
     350  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
     351
     352  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
     353
     354  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
     355  assert (AREA > 0);
     356
     357  off_t Nmax = MaxDensityValue * AREA;
     358  if (subcatalog[0].Naverage <= Nmax) {
     359    if (VERBOSE) {
     360      fprintf (stderr, "subcatalog has less than the max density\n");
     361    }
     362    return (TRUE);
     363  }
     364
     365  off_t Naverage = subcatalog[0].Naverage;
     366
     367  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
     368
     369  // we are going to select Nmax entries by generating a random-sorted index list
     370  off_t *index, tmp, i, j, ave;
     371  ALLOCATE (index, off_t, Naverage);
     372  for (i = 0; i < Naverage; i++) {
     373    index[i] = i;
     374  }
     375  for (i = 0; i < Naverage; i++) {
     376    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
     377    tmp = index[j];
     378    index[j] = index[i];
     379    index[i] = tmp;
     380  }
     381
     382  // count the number of measurements this selection will yield
     383  off_t NMEASURE = 0;
     384  for (i = 0; i < Nmax; i++) {
     385    ave = index[i];
     386    NMEASURE += subcatalog[0].averageT[ave].Nmeasure;
     387  }
     388
     389  // allocate the output data
     390  ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax);
     391  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
     392  ALLOCATE (tmpcatalog.secfilt,  SecFilt, Nmax * Nsecfilt);
     393
     394  off_t Nmeasure = 0;
     395
     396  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
     397  for (i = 0; i < Nmax; i++) {
     398    ave = index[i];
     399    tmpcatalog.averageT[i] = subcatalog[0].averageT[ave];
     400    tmpcatalog.averageT[i].measureOffset = Nmeasure;
     401    for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) {
     402      off_t offset = subcatalog[0].averageT[ave].measureOffset + j;
     403      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
     404      tmpcatalog.measureT[Nmeasure].averef = i;
     405      Nmeasure ++;
     406    }
     407    for (j = 0; j < Nsecfilt; j++) {
     408      tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j];
     409    }
    268410  }
    269411
    … …  
    287429  return (TRUE);
    288430}
    289 
    290 // for the cases where we are not using a subset of the data, we still need to have a copy of these fields
    291 int populate_tiny_values (Catalog *catalog) {
    292 
    293   off_t i;
    294 
    295   AverageTiny *averageT;
    296   Average *average;
    297 
    298   MeasureTiny *measureT;
    299   Measure *measure;
    300 
    301   ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure);
    302   ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage);
    303 
    304   average  = catalog[0].average;
    305   averageT = catalog[0].averageT;
    306 
    307   measure  = catalog[0].measure;
    308   measureT = catalog[0].measureT;
    309 
    310   for (i = 0; i < catalog[0].Naverage; i++) {
    311     // CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]);
    312     averageT[i].R             = average[i].R;
    313     averageT[i].D             = average[i].D;
    314     averageT[i].flags         = average[i].flags;
    315     averageT[i].Nmeasure      = average[i].Nmeasure;
    316     averageT[i].measureOffset = average[i].measureOffset;
    317   }
    318 
    319   for (i = 0; i < catalog[0].Nmeasure; i++) {
    320     // CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
    321     measureT[i].dR       = measure[i].dR;
    322     measureT[i].dD       = measure[i].dD;
    323     measureT[i].M        = measure[i].M;
    324     measureT[i].Mcal     = measure[i].Mcal;
    325     measureT[i].dM       = measure[i].dM;
    326     measureT[i].airmass  = measure[i].airmass;
    327     measureT[i].Xccd     = measure[i].Xccd;
    328     measureT[i].Yccd     = measure[i].Yccd;
    329     measureT[i].t        = measure[i].t;
    330     measureT[i].dt       = measure[i].dt;
    331     measureT[i].averef   = measure[i].averef;
    332     measureT[i].imageID  = measure[i].imageID;
    333     measureT[i].dbFlags  = measure[i].dbFlags;
    334     measureT[i].photcode = measure[i].photcode;
    335   }
    336 
    337   return (TRUE);
    338 }
    339 
    340 int free_tiny_values (Catalog *catalog) {
    341 
    342   if (catalog[0].averageT) free (catalog[0].averageT);
    343   if (catalog[0].measureT) free (catalog[0].measureT);
    344   return (TRUE);
    345 }
    346 

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

    r31450 r33651  
    11# include "relphot.h"
    22
    3 void global_stats (Catalog *catalog, int Ncatalog) {
     3void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    44
    55  StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD;
    … …  
    1717    int seccode = photcodes[Ns][0].code;
    1818
    19     stN = statsStarN (catalog, Ncatalog, Nsec, seccode);
     19    stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr);
    2020    stX = statsStarX (catalog, Ncatalog, Nsec);
    2121    stS = statsStarS (catalog, Ncatalog, Nsec);
    … …  
    2727  }
    2828 
    29   imN = statsImageN (catalog);
     29  imN = statsImageN (catalog, flatcorr);
    3030  imX = statsImageX (catalog);
    3131  imM = statsImageM (catalog);
    3232  imD = statsImagedM (catalog);
    3333 
    34   msN = statsMosaicN (catalog);
     34  msN = statsMosaicN (catalog, flatcorr);
    3535  msM = statsMosaicM (catalog);
    3636  msD = statsMosaicdM (catalog);

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

    r30616 r33651  
    22
    33void relphot_usage (void) {
    4   fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
    5   fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
     4  fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -region RA RA DEC DEC\n");
     5  fprintf (stderr, "       or:    relphot (photcodes) -catalog (name)\n");
    66  fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
    77  fprintf (stderr, "  use -h for more usage information\n");
    … …  
    1818
    1919show_help:
    20   fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
    21   fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
     20  fprintf (stderr, "ERROR: USAGE: relphot (photcode) -region RA RA DEC DEC\n");
     21  fprintf (stderr, "       or:    relphot (photcode) -catalog (name)\n");
    2222  fprintf (stderr, "       or:    relphot -averages -region RA RA DEC DEC\n");
    2323  fprintf (stderr, "  options: \n");
    … …  
    3737  fprintf (stderr, "  -imfreeze\n");
    3838  fprintf (stderr, "  -grid\n");
     39  fprintf (stderr, "  -reset-ubercal : also reset ubercal-ed zero points (otherwise they are sacrosanct)\n");
    3940  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
    4041  fprintf (stderr, "  -instmag min max\n");
    … …  
    4344}
    4445
     46void relphot_client_usage (void) {
     47  fprintf (stderr, "ERROR: USAGE: relphot (photcodes) -load (filename) -hostID (hostID) -hostdir (hostdir) [options]\n");
     48  fprintf (stderr, "       or:    relphot -update (filename) -hostID (hostID) -hostdir (hostdir) [options]\n");
     49  fprintf (stderr, "       or:    relphot -update-objects -hostID (hostID) -hostdir (hostdir) [options]\n");
     50  fprintf (stderr, "  use -h for more usage information\n");
     51  exit (2);
     52}
     53
     54void relphot_client_help (int argc, char **argv) {
     55
     56  /* check for help request */
     57  if (get_argument (argc, argv, "-help")) goto show_help;
     58  if (get_argument (argc, argv, "-h"))    goto show_help;
     59  if (argc == 1) relphot_client_usage();
     60  return;
     61
     62show_help:
     63  fprintf (stderr, "USAGE: relphot_client [-load / -update] (db info)\n\n");
     64  fprintf (stderr, "       relphot_client -load (bcatalog) : extract the bright catalog subset from client's tables\n");
     65  fprintf (stderr, "                            (bcatalog) : location where the bright subset is saved\n");
     66  fprintf (stderr, "       relphot_client -update (images) : apply calculated zero points to the client's tables\n\n");
     67  fprintf (stderr, "                              (images) : location of the table with the image zero points\n");
     68  fprintf (stderr, "       relphot_client -update-objects  : determine average magnitudes for objects\n\n");
     69  fprintf (stderr, "       db info : -hostID (hostID) -hostdir (hostdir) -catdir (catdir)\n");
     70  fprintf (stderr, "other options:\n");
     71  fprintf (stderr, "  -v  : verbose output\n");
     72  fprintf (stderr, "  -vv : extra verbose output\n");
     73  fprintf (stderr, "  \n");
     74  exit (2);
     75}
     76

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

    r31450 r33651  
    22
    33void initialize (int argc, char **argv) {
    4 
    5   int N;
    64
    75  relphot_help (argc, argv);
    … …  
    108
    119  if (!UpdateAverages) {
    12 
    13     N = UserPatchSelect ? 1 : 2;
    14 
    15 # if (0)
    16     // XXX DEP
    17     if ((photcode = GetPhotcodebyName (argv[N])) == NULL) {
    18       fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]);
    19       exit (1);
    20     }
    21     if (photcode[0].type != PHOT_SEC) {
    22       fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[N]);
    23       exit (1);
    24     }
    25     // PhotSec is used to select the single average photcode being processed
    26     PhotSec = GetPhotcodeNsec (photcode[0].code);
    27 # endif
    28 
    29     Nphotcodes = 0;
    30     photcodes = NULL;
    31     int NPHOTCODES = 10;
    32     ALLOCATE (photcodes, PhotCode *, NPHOTCODES);
    33 
    34     /* parse the comma-separated list of photcodesKeep */
    35     char *myList = strcreate(argv[N]);
    36     char *list = myList;
    37     char *codename = NULL;
    38     char *ptr = NULL;
    39     while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
    40       list = NULL; // pass NULL on successive strtok_r calls
    41       fprintf (stderr, "PHOTCODE LIST: %s\n", myList);
    42       fprintf (stderr, "codename: %s\n", codename);
    43       if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) {
    44         fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
    45         exit (1);
    46       }
    47       if (photcodes[Nphotcodes][0].type != PHOT_SEC) {
    48           fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename);
    49           exit (1);
    50       }
    51       Nphotcodes ++;
    52       CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10);
    53     }
     10    // load the list of photcodes into the globals (photcodes, Nphotcodes)
     11    // only a single remaining argument in this mode (the list of photcodes, eg g,r,i)
     12    PhotcodeList = strcreate (argv[1]);
     13    photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes
    5414  }
    55   // XXX DEP PhotNsec = GetPhotcodeNsecfilt ();
    5615  if (USE_GRID && (Nphotcodes > 1)) {
    5716    fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
    … …  
    9554  srand48(B);
    9655}
     56
     57void initialize_client (int argc, char **argv) {
     58
     59  // XXX need to determine which globals can affect relphot_client in either mode and pass appropriately
     60
     61  relphot_client_help (argc, argv);
     62  ConfigInit (&argc, argv);
     63  args_client (argc, argv);
     64
     65  if (MODE == MODE_UPDATE_OBJECTS) return;
     66
     67  // load the list of photcodes into the globals (photcodes, Nphotcodes)
     68  PhotcodeList = strcreate (argv[1]);
     69  photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes
     70}
     71
     72void ParsePhotcodeList_old (char *word) {
     73
     74  Nphotcodes = 0;
     75  photcodes = NULL;
     76  int NPHOTCODES = 10;
     77  ALLOCATE (photcodes, PhotCode *, NPHOTCODES);
     78
     79  /* parse the comma-separated list of photcodesKeep */
     80  char *myList = strcreate(word);
     81  char *list = myList;
     82  char *codename = NULL;
     83  char *ptr = NULL;
     84  while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
     85    list = NULL; // pass NULL on successive strtok_r calls
     86    if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) {
     87      fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
     88      exit (1);
     89    }
     90    if (photcodes[Nphotcodes][0].type != PHOT_SEC) {
     91      fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename);
     92      exit (1);
     93    }
     94    Nphotcodes ++;
     95    CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10);
     96  }
     97}
     98

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

    r16040 r33651  
    4141  stats[0].min    = value[0];
    4242  stats[0].max    = value[N-1];
     43  int N80 = MIN (N-1, 0.8*N);
     44  int N20 = MAX (0, 0.2*N);
     45  stats[0].Upper80 = value[N80];
     46  stats[0].Lower20 = value[N20];
    4347
    4448  switch (statmode) {
    … …  
    106110    dS += M;
    107111  }
    108   X2 = X2 / Nm;
     112  X2 = X2 / (Nm - 1);
    109113  dS = sqrt (dS / Nm);
    110114

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

    r31668 r33651  
    11# include "relphot.h"
    2 
    3 Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) {
     2# include <sys/types.h>
     3# include <sys/wait.h>
     4
     5// this function loops over the catalogs, loads the data, and extracts a bright subset
     6// the returned array (catalog, Ncatalog) has the same layout as the full database, but
     7// only a subset of the detetions & objects
     8
     9// if this function is called in parallel mode, it in turn calls load_catalogs_parallel,
     10// which distributes the work to the remote hosts and loads their results
     11
     12// if this function is called with a specified hostID, then only the fraction of the
     13// database hosted by that hostID is loaded
     14Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath) {
    415
    516  off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total;
    617  Catalog *catalog, tcatalog;
    718
     19  // XXX need to decide how to determine PARALLEL mode...
     20  if (PARALLEL && !hostID) {
     21    catalog = load_catalogs_parallel (skylist, Ncatalog);
     22    return catalog;
     23  }
     24
    825  if (VERBOSE2) fprintf (stderr, "loading catalog data\n");
    926
     27  // a bit of an over-alloc, since we don't load all catalogs for a region
    1028  ALLOCATE (catalog, Catalog, skylist[0].Nregions);
    1129
    … …  
    1432  // load data from each region file, only use bright stars
    1533  for (i = 0; i < skylist[0].Nregions; i++) {
     34    // XXX keep in mind that not all catalogs are loaded
    1635    dvo_catalog_init (&catalog[i], TRUE);
    1736
     37    // does this host ID match the desired location for the table?
     38    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
     39
    1840    // set up the basic catalog info
    19     tcatalog.filename = skylist[0].filename[i];
     41    char hostfile[1024];
     42    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
     43    tcatalog.filename = hostID ? hostfile : skylist[0].filename[i];
     44
    2045    tcatalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    2146    tcatalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    … …  
    3459    }
    3560
     61    if (!tcatalog.sorted) {
     62      fprintf (stderr, "this database is not sorted.  please sort using addstar -resort\n");
     63      exit (1);
     64    }
     65
    3666    Nstar_total += tcatalog.Naverage;
    3767    Nmeas_total += tcatalog.Nmeasure;
    … …  
    4878    Nstar += catalog[i].Naverage;
    4979    Nmeas += catalog[i].Nmeasure;
     80    if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) {
     81      fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n");
     82      fprintf (stderr, "  If you need to do this, and you can afford the RAM, rebuild with IDX_T set to off_t (see relphot.h, ImagesOps.c)\n");
     83      exit (3);
     84    }
    5085  }
    5186  if (Nstar < 2) {
    … …  
    5489
    5590  fprintf (stderr, "using "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measurements)\n", Nstar, Nstar_total, Nmeas, Nmeas_total);
    56   if (Nstar < 1) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n");
     91  if (!hostID && (Nstar < 1)) Shutdown ("%s", "ERROR: no stars match the minimum requirements; exiting \n");
     92  // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency)
    5793   
    5894
    … …  
    67103   need to use an XCLD lock here. 
    68104*/
     105
     106// CATDIR is supplied globally
     107# define DEBUG 1
     108Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
     109
     110  int Nsecfilt  = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
     111
     112  // launch the setphot_client jobs to the parallel hosts
     113
     114  // load the list of hosts
     115  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
     116  if (!table) {
     117    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
     118    exit (1);
     119  }   
     120
     121  int i;
     122  for (i = 0; i < table->Nhosts; i++) {
     123
     124    // ensure that the paths are absolute path names
     125    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     126    free (table->hosts[i].pathname);
     127    table->hosts[i].pathname = tmppath;
     128
     129    ALLOCATE (table->hosts[i].results, char, 1024);
     130    snprintf (table->hosts[i].results, 1024, "%s/relphot.catalog.subset.dat", table->hosts[i].pathname);
     131
     132    // options / arguments that can affect relphot_client -load:
     133    // VERBOSE, VERBOSE2
     134    // KEEP_UBERCAL
     135    // RESET (-reset)
     136    // TimeSelect -time
     137    // DophotSelect
     138    // (note that psfQual is applied rigidly at 0.85, as is the galaxy test)
     139    // MAG_LIM
     140    // SIGMA_LIM
     141    // ImagSelect, ImagMin, ImagMax
     142    // MaxDensityUse, MaxDensityValue
     143
     144    char command[1024];
     145    snprintf (command, 1024, "relphot_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D CAMERA %s -D MAG_LIM %f -D SIGMA_LIM %f",
     146              PhotcodeList, table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA, MAG_LIM, SIGMA_LIM);
     147
     148    char tmpline[1024];
     149    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
     150    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
     151    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
     152    if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal",  command);                    strcpy (command, tmpline); }
     153    if (DophotSelect)  { snprintf (tmpline, 1024, "%s -dophot %d",      command, DophotValue);       strcpy (command, tmpline); }
     154    if (ImagSelect)    { snprintf (tmpline, 1024, "%s -instmag %f %f",  command, ImagMin, ImagMax);  strcpy (command, tmpline); }
     155    if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue);   strcpy (command, tmpline); }
     156    if (TimeSelect) {
     157      char *tstart = ohana_sec_to_date (TSTART);
     158      char *tstop  = ohana_sec_to_date (TSTOP);
     159      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop);
     160      free (tstart);
     161      free (tstop);
     162      strcpy (command, tmpline);
     163    }
     164
     165    fprintf (stderr, "command: %s\n", command);
     166
     167    if (PARALLEL_MANUAL) continue;
     168
     169    if (PARALLEL_SERIAL) {
     170      int status = system (command);
     171      if (status) {
     172        fprintf (stderr, "ERROR running relphot_client\n");
     173        exit (2);
     174      }
     175    } else {
     176      // launch the job on the remote machine (no handshake)
     177      int errorInfo = 0;
     178      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
     179      if (!pid) {
     180        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
     181        continue;
     182      }
     183      table->hosts[i].pid = pid; // save for future reference
     184    }
     185  }
     186
     187  if (PARALLEL_MANUAL) {
     188    fprintf (stderr, "run the relphot_client commands above.  when these are done, hit return\n");
     189    getchar();
     190  }
     191  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
     192    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     193  }
     194
     195  // each host generates a BrightCatalog structure, with the measure, average, etc value
     196  // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt).  I need to
     197  // split out the per-catalog measurements into separate catalog entries.
     198
     199  // set up an initial array of catalogs
     200  CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt);
     201
     202  for (i = 0; i < table->Nhosts; i++) {
     203
     204    BrightCatalog *bcatalog = NULL;
     205    while ((bcatalog = BrightCatalogLoad (table->hosts[i].results)) == NULL) {
     206      // failed to get the data from this host.  This can happen for various reasons.  Give the user a chance to try again...
     207      fprintf (stderr, "failed to read data from %s, stopping operations until this can be fixed\n", table->hosts[i].hostname);
     208      fprintf (stderr, "you may run the command manually and send this process the CONT signal\n");
     209      int pid = getpid();
     210      kill (pid, SIGSTOP);
     211      fprintf (stderr, "retrying %s\n", table->hosts[i].results);
     212    }
     213    free (table->hosts[i].results);
     214    table->hosts[i].results = NULL;
     215   
     216    BrightCatalogSplit (catalogs, bcatalog);
     217
     218    free (bcatalog->average);
     219    free (bcatalog->measure);
     220    free (bcatalog->secfilt);
     221    free (bcatalog);
     222  }
     223
     224  Catalog *catalog = catalogs->catalog;
     225  *Ncatalog = catalogs->Ncatalog;
     226  BrightCatalogSplitFree (catalogs);
     227
     228  int Nmeasure = 0;
     229  int Naverage = 0;
     230  for (i = 0; i < catalogs->Ncatalog; i++) {
     231    Nmeasure += catalogs->catalog[i].Nmeasure;
     232    Naverage += catalogs->catalog[i].Naverage;
     233  }
     234
     235  fprintf (stderr, "loaded %d catalogs, using a total of %d stars (%d measures)\n", catalogs->Ncatalog, Naverage, Nmeasure);
     236
     237  return (catalog);
     238}     

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

    r32346 r33651  
    99// This function generates a subset of the images based on selections.  Input db has already
    1010// been loaded with the raw fits table data
    11 SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
     11SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region) {
    1212
    1313  Image     *image, *subset;
    … …  
    3535
    3636  // determine the populated SkyRegions overlapping the requested area
    37   if (RegionSelect) {
     37
     38  // if the user selects a specific catalog, do that one;
     39  // otherwise limit to the selection region
     40  if (regionName) {
     41    Nchar = strlen(regionName);
     42    if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
     43    skylist = SkyListByName (sky, regionName);
     44
     45    // select the images which overlap the selected sky regions
     46    // 'subset' points to a new copy of the data (different from 'image')
     47    subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
     48    MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
     49  } else {
    3850    if (region[0].Rmin > region[0].Rmax) {
    3951      SkyRegion subregion;
    … …  
    7486      MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    7587    }
    76   } else {
    77     Nchar = strlen(regionName);
    78     if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
    79     skylist = SkyListByName (sky, regionName);
    80 
    81     // select the images which overlap the selected sky regions
    82     // 'subset' points to a new copy of the data (different from 'image')
    83     subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
    84     MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
    85   }
     88  }
    8689
    8790  // save the subset of images in the static reference in ImageOps, set up indexes

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

    r31450 r33651  
    11# include "relphot.h"
    22   
    3 void plot_scatter (Catalog *catalog, int Ncatalog) {
     3void plot_scatter (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    44
    55  off_t i, j, k, m, N, Ntot;
    … …  
    1111  for (i = 0; i < Ncatalog; i++) {
    1212    for (j = 0; j < catalog[i].Naverage; j++) {
    13       Ntot += catalog[i].average[j].Nmeasure;
     13      Ntot += catalog[i].averageT[j].Nmeasure;
    1414    }
    1515  }
    … …  
    3232            /* calculate the average value for a single star */
    3333            if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue; 
    34             m = catalog[i].average[j].measureOffset;
     34            m = catalog[i].averageT[j].measureOffset;
    3535
    36             for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
     36            for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
    3737                // skip measurements that do not match the current photcode
    3838                int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
    … …  
    4040
    4141                if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
    42                 Mcal = getMcal  (m, i);
     42                Mcal = getMcal  (m, i, flatcorr, catalog);
    4343                if (isnan(Mcal)) continue;
    4444                Mmos = getMmos  (m, i);

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

    r27611 r33651  
    165165  graphdata[0].color = black;
    166166  graphdata[0].lweight = 0;
    167   graphdata[0].size = 0.5;
     167  graphdata[0].size = 1.0;
    168168
    169169  graphdata[0].xmin = dUNDEF;

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

    r31668 r33651  
    11# include "relphot.h"
     2
     3# define DEBUG 1
    24
    35# define TIMESTAMP(TIME) \
    … …  
    79    gettimeofday (&start, (void *) NULL);
    810
    9 void reload_catalogs (SkyList *skylist) {
     11void reload_catalogs (SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath) {
    1012
    1113  int i;
    … …  
    2527  double time7 = 0.0;
    2628
     29  // XXX need to decide how to determine PARALLEL mode...
     30  if (PARALLEL && !hostID) {
     31    reload_catalogs_parallel (skylist);
     32    return;
     33  }
     34
    2735  if (VERBOSE) fprintf (stderr, "re-loading catalog data\n");
    2836
    2937  /* load data from each region file */
    3038  for (i = 0; i < skylist[0].Nregions; i++) {
     39
     40    // does this host ID match the desired location for the table?
     41    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
     42
    3143    gettimeofday (&start, (void *) NULL);
    32     catalog.filename = skylist[0].filename[i];
     44
     45    // set up the basic catalog info
     46    char hostfile[1024];
     47    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
     48    catalog.filename = hostID ? hostfile : skylist[0].filename[i];
    3349
    3450    // only update existing db tables
    … …  
    5773    TIMESTAMP(time2);
    5874
    59     populate_tiny_values(&catalog);
     75    populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE);
    6076    TIMESTAMP(time3a);
    6177
     78    // XXX need to worry about the image subset data
    6279    initImageBins  (&catalog, 1, FALSE);
    6380    initMosaicBins (&catalog, 1, FALSE);
    … …  
    6986    TIMESTAMP(time4);
    7087
    71     setMrelFinal (&catalog);
     88    setMrelFinal (&catalog, flatcorr);
    7289    TIMESTAMP(time5);
    7390
     91    // modify the output format as desired (ignore current format on disk)
     92    if (UPDATE_CATFORMAT) {
     93      catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT);
     94    }
    7495    dvo_catalog_save (&catalog, VERBOSE);
    7596    dvo_catalog_unlock (&catalog);
    … …  
    85106  }
    86107
    87   fprintf (stderr, "time1  %f : find catalog\n", time1);
    88   fprintf (stderr, "time2  %f : load catalog\n", time2);
    89   fprintf (stderr, "time3a %f : init imbins\n", time3a);
    90   fprintf (stderr, "time3b %f : init imbins\n",  time3b);
    91   fprintf (stderr, "time4  %f : find images\n",  time4);
    92   fprintf (stderr, "time5  %f : set Mrel\n",     time5);
    93   fprintf (stderr, "time6  %f : save catalog\n", time6);
    94   fprintf (stderr, "time7  %f : free catalog\n", time7);
     108  fprintf (stderr, "time step 1  %10.3f sec : find catalog\n", time1);
     109  fprintf (stderr, "time step 2  %10.3f sec : load catalog\n", time2);
     110  fprintf (stderr, "time step 3  %10.3f sec : make tiny values\n", time3a);
     111  fprintf (stderr, "time step 4  %10.3f sec : init imbins\n",  time3b);
     112  fprintf (stderr, "time step 5  %10.3f sec : find images\n",  time4);
     113  fprintf (stderr, "time step 6  %10.3f sec : set Mrel\n",     time5);
     114  fprintf (stderr, "time step 7  %10.3f sec : save catalog\n", time6);
     115  fprintf (stderr, "time step 8  %10.3f sec : free catalog\n", time7);
    95116}
     117
     118// XXX Image to Image Subset
     119int reload_catalogs_parallel (SkyList *sky) {
     120
     121  off_t Nimage;
     122  ImageSubset *image = getimages_subset (&Nimage);
     123
     124  // write out the subset table of image information
     125  char imageFile[512];
     126  snprintf (imageFile, 512, "%s/Images.subset.dat", CATDIR);
     127  if (!ImageSubsetSave (imageFile, image, Nimage)) {
     128    fprintf (stderr, "failed to write image subset\n");
     129    exit (1);
     130  }
     131  free (image);
     132
     133  // now launch the relphot_client jobs to the parallel hosts
     134
     135  // load the list of hosts
     136  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
     137
     138  int i;
     139  for (i = 0; i < table->Nhosts; i++) {
     140
     141    // ensure that the paths are absolute path names
     142    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     143    free (table->hosts[i].pathname);
     144    table->hosts[i].pathname = tmppath;
     145
     146    char command[1024];
     147    snprintf (command, 1024, "relphot_client %s -update-catalogs %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f",
     148              PhotcodeList, imageFile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR);
     149
     150    // options & configs which affect relphot_client -update-catalogs
     151    // VERBOSE, VERBOSE2
     152    // RESET
     153    // TimeSelect
     154    // AreaSelect
     155    // STATMODE
     156    // STAR_TOOFEW
     157    // MIN_ERROR
     158   
     159    char tmpline[1024];
     160    if (VERBOSE)          { snprintf (tmpline, 1024, "%s -v",                   command);                                         strcpy (command, tmpline); }
     161    if (VERBOSE2)         { snprintf (tmpline, 1024, "%s -vv",                  command);                                         strcpy (command, tmpline); }
     162    if (RESET)            { snprintf (tmpline, 1024, "%s -reset",               command);                                         strcpy (command, tmpline); }
     163    if (UPDATE)           { snprintf (tmpline, 1024, "%s -update",              command);                                         strcpy (command, tmpline); }
     164    if (!KEEP_UBERCAL)    { snprintf (tmpline, 1024, "%s -reset-ubercal",       command);                                         strcpy (command, tmpline); }
     165    if (UPDATE_CATFORMAT) { snprintf (tmpline, 1024, "%s -update-catformat %s", command, UPDATE_CATFORMAT);                       strcpy (command, tmpline); }
     166    if (AreaSelect)       { snprintf (tmpline, 1024, "%s -area %f %f %f %f",    command, AreaXmin, AreaXmax, AreaYmin, AreaYmax); strcpy (command, tmpline); }
     167    if (TimeSelect) {
     168      char *tstart = ohana_sec_to_date (TSTART);
     169      char *tstop  = ohana_sec_to_date (TSTOP);
     170      snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop);
     171      free (tstart);
     172      free (tstop);
     173      strcpy (command, tmpline);
     174    }
     175
     176    fprintf (stderr, "command: %s\n", command);
     177
     178    if (PARALLEL_MANUAL) continue;
     179
     180    if (PARALLEL_SERIAL) {
     181      int status = system (command);
     182      if (status) {
     183        fprintf (stderr, "ERROR running relphot_client\n");
     184        exit (2);
     185      }
     186    } else {
     187      // launch the job on the remote machine (no handshake)
     188      int errorInfo = 0;
     189      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
     190      if (!pid) {
     191        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
     192        exit (1);
     193      }
     194      table->hosts[i].pid = pid; // save for future reference
     195    }
     196  }
     197
     198  if (PARALLEL_MANUAL) {
     199    fprintf (stderr, "run the relphot_client commands above.  when these are done, hit return\n");
     200    getchar();
     201  }
     202  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
     203    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     204  }
     205  return (TRUE);
     206}     

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

    r31668 r33651  
    11# include "relphot.h"
    2 # define USE_DIRECT 0
    3 
    4 # define TIMESTAMP(TIME) \
     2
     3# define MARKTIME(MSG,...) {                    \
     4    float dtime;                                \
    55    gettimeofday (&stop, (void *) NULL);        \
    66    dtime = DTIME (stop, start);                \
    7     TIME += dtime;                              \
    8     gettimeofday (&start, (void *) NULL);
    9 
    10 # define MARKTIME(MSG,...) { \
    11   float dtime; \
    12   gettimeofday (&stop, (void *) NULL); \
    13   dtime = DTIME (stop, start); \
    14   fprintf (stderr, MSG, __VA_ARGS__); }
     7    fprintf (stderr, MSG, __VA_ARGS__); }
    158
    169int main (int argc, char **argv) {
    1710
    1811  int i, status, Ncatalog;
    19   Catalog *catalog;
     12  Catalog *catalog = NULL;
    2013  FITS_DB db;
    2114  struct timeval start, stop;
    … …  
    2417  gettimeofday (&start, (void *) NULL);
    2518
    26   // XXX quick and stupid test:
    27   if (0) {
    28     int i, j;
    29     off_t *Nlist, *NLIST, **mlist;
    30     off_t *NlistI, *NLISTI, **mlistI;
    31     int **clist;
    32     int **clistI;
    33     float dtime;
    34     float time1 = 0.0;
    35     float time2 = 0.0;
    36 
    37     int Nmosaic = 10000;
    38     int Nimage = 600000;
    39 
    40     for (j = 0; j < 19000; j++) {
    41 
    42       // mosaic indexes
    43       ALLOCATE (Nlist, off_t,   Nmosaic);
    44       ALLOCATE (NLIST, off_t,   Nmosaic);
    45       ALLOCATE (clist, int *,   Nmosaic);
    46       ALLOCATE (mlist, off_t *, Nmosaic);
    47    
    48       for (i = 0; i < Nmosaic; i++) {
    49         Nlist[i] = 0;
    50         NLIST[i] = 100;
    51         ALLOCATE (clist[i], int,   NLIST[i]);
    52         ALLOCATE (mlist[i], off_t, NLIST[i]);
    53       }
    54 
    55       // image indexes
    56       ALLOCATE (NlistI, off_t,   Nimage);
    57       ALLOCATE (NLISTI, off_t,   Nimage);
    58       ALLOCATE (clistI, int *,   Nimage);
    59       ALLOCATE (mlistI, off_t *, Nimage);
    60    
    61       for (i = 0; i < Nimage; i++) {
    62         NlistI[i] = 0;
    63         NLISTI[i] = 100;
    64         ALLOCATE (clistI[i], int,   NLISTI[i]);
    65         ALLOCATE (mlistI[i], off_t, NLISTI[i]);
    66       }
    67 
    68       TIMESTAMP(time1);
    69 
    70       // free mosaic indexes
    71       for (i = 0; i < Nmosaic; i++) {
    72         free (clist[i]);
    73         free (mlist[i]);
    74       }
    75       free (Nlist);
    76       free (NLIST);
    77       free (clist);
    78       free (mlist);
    79 
    80       // free image indexes
    81       for (i = 0; i < Nimage; i++) {
    82         free (clistI[i]);
    83         free (mlistI[i]);
    84       }
    85       free (NlistI);
    86       free (NLISTI);
    87       free (clistI);
    88       free (mlistI);
    89 
    90       TIMESTAMP(time2);
    91 
    92       if (j % 100 == 0) {
    93         fprintf (stderr, "done with %d\n", j);
    94         fprintf (stderr, "time1  %f : init mosaic\n", time1);
    95         fprintf (stderr, "time2  %f : free mosaic\n", time2);
    96       }
    97     }
    98    
    99     fprintf (stderr, "time1  %f : init mosaic\n", time1);
    100     fprintf (stderr, "time2  %f : free mosaic\n", time2);
    101     exit (1);
    102   }
    103 
    10419  /* get configuration info, args */
    10520  initialize (argc, argv);
    … …  
    10722  /* the object analysis is a separate process iterating over catalogs */
    10823  if (UpdateAverages) {
    109     relphot_objects ();
     24    relphot_objects (0, NULL);
    11025    exit (0);
    11126  }
    … …  
    12540  // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE)
    12641  if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) {
    127     // XXX get ZERO_POINT from config
    128     dvo_image_create (&db, 25.0);
    129     // Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
     42
     43    // 25.0 is the nominal zero point for measurements in the database
     44    // for all measurements except PHOT_REF:
     45    // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0
     46    // for measurements using PHOT_REF:
     47    // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known)
     48    dvo_image_create (&db, 25.0);
     49
     50    // XXX why are we running relphot on an empty database?
    13051  } else {
    13152    if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
    … …  
    13354  MARKTIME("-- load image data: %f sec\n", dtime);
    13455
    135   /* load regions and images based on specified sky patch */
    136   // XXX need to mimic old-style load by passing patch name
    137   // XXX need to reduce number of global variables in use.
    138   // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
    139   skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
     56  /* load regions and images based on specified sky patch and/or catalog */
     57  skylist = load_images (&db, UserCatalog, &UserPatch);
    14058  MARKTIME("-- load images: %f sec\n", dtime);
    14159
    … …  
    14361  if (!UPDATE) dvo_image_unlock (&db);
    14462
    145   /* load catalog data from region files */
    146   catalog = load_catalogs (skylist, &Ncatalog);
    147   MARKTIME("-- load catalog data: %f sec\n", dtime);
     63  // load the flat correction table (if defined)
     64  char flatcorrfile[256];
     65  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
     66  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
     67
     68  if (NLOOP > 0) {
     69    /* load catalog data from region files (hostID is 0 since we are not a client */
     70    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
     71    MARKTIME("-- load catalog data: %f sec\n", dtime);
    14872 
    149   /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
    150 
    151   /* match measurements with images, mosaics */
    152   initImageBins  (catalog, Ncatalog, TRUE);
    153   MARKTIME("-- make image bins: %f sec\n", dtime);
    154 
    155   initMosaicBins (catalog, Ncatalog, TRUE);
    156   initGridBins   (catalog, Ncatalog);
    157   initMrel (catalog, Ncatalog);
    158 
    159   findImages (catalog, Ncatalog, TRUE);
    160   MARKTIME("-- set up image indexes: %f sec\n", dtime);
    161 
    162   findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
    163   MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
    164 
    165   SAVEPLOT = FALSE;
    166 
    167   setExclusions (catalog, Ncatalog);
    168 
    169   global_stats (catalog, Ncatalog);
    170 
    171   if (PLOTSTUFF) {
    172     plot_star_coords (catalog, Ncatalog);
    173     // plot_mosaic_fields (catalog);
    174   }
    175 
    176   // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
    177   if (USE_GRID) {
     73    /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
     74
     75    /* match measurements with images, mosaics */
     76    initImageBins  (catalog, Ncatalog, TRUE);
     77    MARKTIME("-- make image bins: %f sec\n", dtime);
     78
     79    initMosaicBins (catalog, Ncatalog, TRUE);
     80    initGridBins   (catalog, Ncatalog);
     81    initMrel (catalog, Ncatalog);
     82
     83    findImages (catalog, Ncatalog, TRUE);
     84    MARKTIME("-- set up image indexes: %f sec\n", dtime);
     85
     86    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
     87    MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
     88
     89    SAVEPLOT = FALSE;
     90
     91    setExclusions (catalog, Ncatalog, TRUE);
     92
     93    global_stats (catalog, Ncatalog, flatcorr);
     94
     95    if (PLOTSTUFF) {
     96      plot_star_coords (catalog, Ncatalog);
     97      // plot_mosaic_fields (catalog);
     98    }
     99
     100    // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
     101    if (USE_GRID) {
    178102      int star_toofew;
    179103
    … …  
    183107      STAR_TOOFEW = 0;
    184108      for (i = 0; i < NGRID; i++) {
    185           STAR_BAD = ID_STAR_POOR;
    186           setMrel  (catalog, Ncatalog);
    187           STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    188           setMgrid (catalog);
     109        STAR_BAD = ID_STAR_POOR;
     110        setMrel  (catalog, Ncatalog, flatcorr);
     111        STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     112        setMgrid (catalog, flatcorr);
    189113      }
    190114      STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
    191115      STAR_TOOFEW = star_toofew;
    192   }
    193 
    194   /* determine fit values */
    195   for (i = 0; i < NLOOP; i++) {
    196     setMrel  (catalog, Ncatalog);
    197     setMcal  (catalog, FALSE);
    198     setMmos  (catalog, FALSE);
    199     setMgrid (catalog);
     116    }
     117
     118    /* determine fit values */
     119    for (i = 0; i < NLOOP; i++) {
     120      if (PLOTSTUFF) {
     121        plot_scatter (catalog, Ncatalog, flatcorr);
     122      }
     123      setMrel  (catalog, Ncatalog, flatcorr);
     124      if (PLOTSTUFF) {
     125        plot_scatter (catalog, Ncatalog, flatcorr);
     126      }
     127      setMcal  (catalog, FALSE, flatcorr);
     128      setMmos  (catalog, FALSE, flatcorr);
     129      setMgrid (catalog, flatcorr);
    200130   
     131      if (PLOTSTUFF) {
     132        plot_scatter (catalog, Ncatalog, flatcorr);
     133        plot_grid (catalog, flatcorr);
     134        plot_mosaics ();
     135        plot_images ();
     136        plot_stars (catalog, Ncatalog);
     137        plot_chisq (catalog, Ncatalog);
     138      }
     139      // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
     140      // if (i % 6 == 1) rationalize_images ();
     141      if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE, flatcorr);
     142      if (i % 6 == 3) clean_stars (catalog, Ncatalog);
     143      if (i % 6 == 5) clean_mosaics ();
     144      if (i % 6 == 5) clean_images ();
     145
     146      // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
     147      // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
     148      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
     149      // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
     150      global_stats (catalog, Ncatalog, flatcorr);
     151      MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
     152    }
     153
    201154    if (PLOTSTUFF) {
    202       plot_scatter (catalog, Ncatalog);
    203       plot_grid (catalog);
     155      plot_scatter (catalog, Ncatalog, flatcorr);
     156      plot_grid (catalog, flatcorr);
    204157      plot_mosaics ();
    205158      plot_images ();
    … …  
    207160      plot_chisq (catalog, Ncatalog);
    208161    }
    209     // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog);
    210     // if (i % 6 == 1) rationalize_images ();
    211     if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE);
    212     if (i % 6 == 3) clean_stars (catalog, Ncatalog);
    213     if (i % 6 == 5) clean_mosaics ();
    214     if (i % 6 == 5) clean_images ();
    215 
    216     // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE);
    217     // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
    218     // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
    219     // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
    220     global_stats (catalog, Ncatalog);
    221     MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
    222   }
    223 
    224   if (PLOTSTUFF) {
    225     plot_scatter (catalog, Ncatalog);
    226     plot_grid (catalog);
    227     plot_mosaics ();
    228     plot_images ();
    229     plot_stars (catalog, Ncatalog);
    230     plot_chisq (catalog, Ncatalog);
    231   }
    232162 
    233   if (USE_GRID) dump_grid ();
    234 
    235   /* set Mcal & Mmos for bad images */
    236   setMcal  (catalog, TRUE);
    237   setMmos  (catalog, TRUE);
    238   MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
    239 
    240   setMcalFinal (); // copy per-mosaic calibrations to the images
    241 
    242   if (SAVE_IMAGE_UPDATES) {
    243 
    244     FITS_DB dbX;
    245     off_t *LineNumber;
    246     Image *image, *subset;
    247     off_t Nimage, Nsubset, i, Nx;
    248     char filename[1024];
    249 
    250     gfits_db_init (&dbX);
    251     dbX.lockstate = LCK_XCLD;
    252     dbX.timeout   = 60.0;
    253     dbX.mode      = db.mode;
    254     dbX.format    = db.format;
    255 
    256     snprintf (filename, 1024, "%s.bck", ImageCat);
    257     if (!gfits_db_lock (&dbX, filename)) {
    258       fprintf (stderr, "can't lock backup image image catalog\n");
    259       return (FALSE);
    260     }
     163    if (USE_GRID) dump_grid ();
     164
     165    /* set Mcal & Mmos for bad images */
     166    setMcal  (catalog, TRUE, flatcorr);
     167    setMmos  (catalog, TRUE, flatcorr);
     168    MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
     169
     170    setMcalFinal (); // copy per-mosaic calibrations to the images
     171
     172    if (SAVE_IMAGE_UPDATES) {
     173
     174      FITS_DB dbX;
     175      off_t *LineNumber;
     176      Image *image, *subset;
     177      off_t Nimage, Nsubset, i, Nx;
     178      char filename[1024];
     179
     180      gfits_db_init (&dbX);
     181      dbX.lockstate = LCK_XCLD;
     182      dbX.timeout   = 60.0;
     183      dbX.mode      = db.mode;
     184      dbX.format    = db.format;
     185
     186      snprintf (filename, 1024, "%s.bck", ImageCat);
     187      if (!gfits_db_lock (&dbX, filename)) {
     188        fprintf (stderr, "can't lock backup image image catalog\n");
     189        return (FALSE);
     190      }
    261191   
    262     // apply the changes from the image subset to the full table:
    263 
    264     // convert database table to internal structure (binary to Image)
    265     // 'image' points to the same memory as db->ftable->buffer
    266     image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
    267     if (!image) {
     192      // apply the changes from the image subset to the full table:
     193
     194      // convert database table to internal structure (binary to Image)
     195      // 'image' points to the same memory as db->ftable->buffer
     196      image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped);
     197      if (!image) {
    268198        fprintf (stderr, "ERROR: failed to read images\n");
    269199        exit (2);
    270     }
    271     gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
    272 
    273     subset = getimages (&Nsubset, &LineNumber);
    274     for (i = 0; i < Nsubset; i++) {
    275       if (LineNumber[i] == -1) continue;
    276       memcpy (&image[LineNumber[i]], &subset[i], Nx);
    277     }
    278 
    279     // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
    280     // the entire Image table, but with the updates in place I think this says: re-work the
    281     // ftable/vtable usage in this program to be more sensible...
    282     gfits_copy_header (&db.header,  &dbX.header);
    283     gfits_copy_matrix (&db.matrix,  &dbX.matrix);
    284     gfits_copy_header (&db.theader, &dbX.theader);
    285     gfits_copy_ftable (&db.ftable,  &dbX.ftable);
    286 
    287     dbX.ftable.header = &dbX.theader;
    288     dbX.virtual = FALSE;
    289 
    290     // copy Images.dat to another structure
    291     dvo_image_save (&dbX, VERBOSE);
    292     dvo_image_unlock (&dbX);
    293     MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
    294   }
    295 
    296   // only change the real database if -update is requested
     200      }
     201      gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
     202
     203      subset = getimages (&Nsubset, &LineNumber);
     204      for (i = 0; i < Nsubset; i++) {
     205        if (LineNumber[i] == -1) continue;
     206        memcpy (&image[LineNumber[i]], &subset[i], Nx);
     207      }
     208
     209      // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
     210      // the entire Image table, but with the updates in place I think this says: re-work the
     211      // ftable/vtable usage in this program to be more sensible...
     212      gfits_copy_header (&db.header,  &dbX.header);
     213      gfits_copy_matrix (&db.matrix,  &dbX.matrix);
     214      gfits_copy_header (&db.theader, &dbX.theader);
     215      gfits_copy_ftable (&db.ftable,  &dbX.ftable);
     216
     217      dbX.ftable.header = &dbX.theader;
     218      dbX.virtual = FALSE;
     219
     220      // save Images.dat using the copied structure
     221      if (UPDATE_CATFORMAT) {
     222        // ensure the db format is updated
     223        dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
     224        char photcodeFile[1024];
     225        sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR);
     226        SavePhotcodesFITS (photcodeFile);
     227      }
     228      dvo_image_save (&dbX, VERBOSE);
     229      dvo_image_unlock (&dbX);
     230      MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
     231    }
     232
     233    /* at this point, we have correct cal coeffs in the image/mosaic structures */
     234    for (i = 0; i < Ncatalog; i++) {
     235      free_tiny_values (&catalog[i]);
     236      dvo_catalog_free (&catalog[i]);
     237    }
     238    freeImageBins (Ncatalog, TRUE);
     239    freeMosaicBins (Ncatalog, TRUE);
     240    freeGridBins (Ncatalog);
     241
     242  } // (NLOOP > 0) : this loop determines the offsets per chip
     243
     244  reload_images (&db);
     245
     246  // only change the real database files if -update is requested
    297247  if (!UPDATE) exit (0);
    298248 
    299   reload_images (&db);
    300 
    301   /* at this point, we have correct cal coeffs in the image/mosaic structures */
    302   for (i = 0; i < Ncatalog; i++) {
    303     free_tiny_values (&catalog[i]);
    304     dvo_catalog_free (&catalog[i]);
    305   }
    306   freeImageBins (Ncatalog, TRUE);
    307   freeMosaicBins (Ncatalog, TRUE);
    308   freeGridBins (Ncatalog);
    309 
    310249  /* load catalog data from region files, update Mrel include all data */
    311   reload_catalogs (skylist);
     250  reload_catalogs (skylist, flatcorr, 0, NULL);
    312251  MARKTIME("-- updated all catalogs: %f sec\n", dtime);
    313252
     253  if (UPDATE_CATFORMAT) {
     254    // ensure the db format is updated
     255    db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
     256  }
    314257  dvo_image_update (&db, VERBOSE);
    315258  dvo_image_unlock (&db);

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

    r31450 r33651  
    11# include "relphot.h"
    22
    3 int relphot_objects () {
     3// apply the current calibration values to the measurements to determine the average
     4// magnitudes, applying the desired clipping analysis.  Note that this does not depend on
     5// the form of the image zero point measurement.  If there is a grid or map term, and this
     6// has been correctly propagated to the measurement, the average will respect that value.
     7
     8int relphot_objects_parallel (SkyList *sky);
     9
     10int relphot_objects (int hostID, char *hostpath) {
    411
    512  off_t i, j, k;
    … …  
    1724  skylist = SkyListByPatch (sky, -1, &UserPatch);
    1825
     26  // XXX need to decide how to determine PARALLEL mode...
     27  if (PARALLEL && !hostID) {
     28    relphot_objects_parallel (skylist);
     29    return TRUE;
     30  }
     31
    1932  // load data from each region file, only use bright stars
    2033  for (i = 0; i < skylist[0].Nregions; i++) {
    2134
     35    // does this host ID match the desired location for the table?
     36    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
     37
    2238    // set up the basic catalog info
    23     catalog.filename  = skylist[0].filename[i];
     39    // set up the basic catalog info
     40    char hostfile[1024];
     41    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
     42    catalog.filename  = hostID ? hostfile : skylist[0].filename[i];
    2443    catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    2544    catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    … …  
    3857    }
    3958
    40     // reset
     59    // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images)
    4160    if (RESET) {
    4261      Nsecfilt = catalog.Nsecfilt;
     62      DVOSecfiltFlags secfiltBits =
     63        ID_SECF_STAR_FEW |
     64        ID_SECF_STAR_POOR |
     65        ID_PHOTOM_PASS_0 |
     66        ID_PHOTOM_PASS_1 |
     67        ID_PHOTOM_PASS_2 |
     68        ID_PHOTOM_PASS_3 |
     69        ID_PHOTOM_PASS_4 |
     70        ID_SECF_USE_SYNTH |
     71        ID_SECF_USE_UBERCAL |
     72        ID_SECF_OBJ_EXT;
     73     
    4374      for (j = 0; j < catalog.Naverage; j++) {
    4475        catalog.average[j].flags = 0;
    4576        for (k = 0; k < Nsecfilt; k++) {
    46           catalog.secfilt[j*Nsecfilt + k].M     = NAN;
    47           catalog.secfilt[j*Nsecfilt + k].M_20  = NAN_S_SHORT;
    48           catalog.secfilt[j*Nsecfilt + k].M_80  = NAN_S_SHORT;
    49           catalog.secfilt[j*Nsecfilt + k].dM    = NAN;
    50           catalog.secfilt[j*Nsecfilt + k].Xm    = NAN_S_SHORT;
    51           catalog.secfilt[j*Nsecfilt + k].Ncode = 0;
    52           catalog.secfilt[j*Nsecfilt + k].Nused = 0;
    53           // XXX reset the photometry flags for secfilt entries?
     77          catalog.secfilt[j*Nsecfilt + k].M      = NAN;
     78          catalog.secfilt[j*Nsecfilt + k].Map    = NAN;
     79          catalog.secfilt[j*Nsecfilt + k].dM     = NAN;
     80          catalog.secfilt[j*Nsecfilt + k].Xm     = NAN;
     81          catalog.secfilt[j*Nsecfilt + k].M_20   = NAN_S_SHORT;
     82          catalog.secfilt[j*Nsecfilt + k].M_80   = NAN_S_SHORT;
     83          catalog.secfilt[j*Nsecfilt + k].Ncode  = 0;
     84          catalog.secfilt[j*Nsecfilt + k].Nused  = 0;
     85          catalog.secfilt[j*Nsecfilt + k].flags &= ~secfiltBits;
     86          catalog.secfilt[j*Nsecfilt + k].ubercalDist = 1000;
    5487        }
    5588      }
    … …  
    76109  return (TRUE);
    77110}
     111
     112// CATDIR is supplied globally
     113# define DEBUG 1
     114int relphot_objects_parallel (SkyList *sky) {
     115
     116  // launch the setphot_client jobs to the parallel hosts
     117
     118  // load the list of hosts
     119  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
     120  if (!table) {
     121    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
     122    exit (1);
     123  }   
     124
     125  int i;
     126  for (i = 0; i < table->Nhosts; i++) {
     127
     128    // ensure that the paths are absolute path names
     129    char *tmppath = abspath (table->hosts[i].pathname, MAX_PATH_LENGTH);
     130    free (table->hosts[i].pathname);
     131    table->hosts[i].pathname = tmppath;
     132
     133    // options / arguments that can affect relphot_client -update-objects:
     134    // VERBOSE, VERBOSE2
     135    // KEEP_UBERCAL
     136    // RESET (-reset)
     137    // TimeSelect -time
     138    // DophotSelect
     139    // (note that psfQual is applied rigidly at 0.85, as is the galaxy test)
     140    // MAG_LIM
     141    // SIGMA_LIM
     142    // ImagSelect, ImagMin, ImagMax
     143    // MaxDensityUse, MaxDensityValue
     144
     145    char command[1024];
     146    snprintf (command, 1024, "relphot_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f",
     147              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR);
     148
     149    char tmpline[1024];
     150    if (VERBOSE)       { snprintf (tmpline, 1024, "%s -v",              command);                    strcpy (command, tmpline); }
     151    if (VERBOSE2)      { snprintf (tmpline, 1024, "%s -vv",             command);                    strcpy (command, tmpline); }
     152    if (RESET)         { snprintf (tmpline, 1024, "%s -reset",          command);                    strcpy (command, tmpline); }
     153    if (UPDATE)        { snprintf (tmpline, 1024, "%s -update",         command);                    strcpy (command, tmpline); }
     154    if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal",  command);                    strcpy (command, tmpline); }
     155
     156    fprintf (stderr, "command: %s\n", command);
     157
     158    if (PARALLEL_MANUAL) continue;
     159
     160    if (PARALLEL_SERIAL) {
     161      int status = system (command);
     162      if (status) {
     163        fprintf (stderr, "ERROR running relphot_client\n");
     164        exit (2);
     165      }
     166    } else {
     167      // launch the job on the remote machine (no handshake)
     168      int errorInfo = 0;
     169      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
     170      if (!pid) {
     171        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
     172        exit (1);
     173      }
     174      table->hosts[i].pid = pid; // save for future reference
     175    }
     176  }
     177
     178  if (PARALLEL_MANUAL) {
     179    fprintf (stderr, "run the relphot_client commands above.  when these are done, hit return\n");
     180    getchar();
     181  }
     182  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
     183    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
     184  }
     185
     186  return TRUE;
     187}     

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

    r31450 r33651  
    2525 
    2626  Image *image;
    27   off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
     27  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE, D_NIMAGE;
    2828  off_t *line_number;
    2929  int InRange, ecode, found;
    … …  
    6363  DminSkyRegion = +90.0;
    6464  DmaxSkyRegion = -90.0;
     65
     66  D_NIMAGE = 6000;
    6567
    6668  // FILE *ftest = fopen ("relphot.dump.dat", "w");
    … …  
    110112
    111113  nimage = 0;
    112   NIMAGE = 100;
     114  NIMAGE = D_NIMAGE;
    113115  ALLOCATE (image, Image, NIMAGE);
    114116  ALLOCATE (line_number, off_t, NIMAGE);
    … …  
    117119  for (i = 0; i < Ntimage; i++) {
    118120     
     121    if (!(i % 300000)) fprintf (stderr, ".");
     122
    119123    /* exclude images by photcode */
    120124    ecode = GetPhotcodeEquivCodebyCode (timage[i].photcode);
    … …  
    137141      continue;
    138142    }
    139 
    140     // XXX temporary test : record center coords for each accepted and rejected chip/mosaic
    141     // double RAo, DECo;
    142143
    143144    /* define image corners - note the DIS images (mosaic phu) are special */
    … …  
    148149      Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY;
    149150      Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY;
    150       // XY_to_RD(&RAo, &DECo, 0.0, 0.0, &timage[i].coords);
    151151    } else {
    152152      Xi[0] = 0;            Yi[0] = 0;
    … …  
    155155      Xi[3] = 0;            Yi[3] = timage[i].NY;
    156156      Xi[4] = 0;            Yi[4] = 0;
    157       // XY_to_RD(&RAo, &DECo, 0.5*timage[i].NX, 0.5*timage[i].NY, &timage[i].coords);
    158157    }
    159158    found = FALSE;
    … …  
    230229
    231230  found_it:
    232     // XXX We claim this is a good image: write to a test file
    233     // fprintf (ftest, "%s : %lf %lf  : %f %f %d %x\n", timage[i].name, RAo, DECo, timage[i].Mcal, timage[i].dMcal, timage[i].nstar, timage[i].flags);
    234 
    235231    image[nimage] = timage[i];
    236     /* always allow 'few' images to succeed, if possible */
     232    /* always allow 'few' images to succeed, if possible (new images / detections may have
     233     * been added) */
    237234    if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) {
    238235      image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR);
    239236    }
    240237    if (RESET) {
    241       // XXX assignMcal (&image[nimage], (double *) NULL, -1);
    242       // XXX this needs to be thought through a bit more
    243       image[nimage].Mcal = 0.0;
    244       image[nimage].dMcal = NAN;
     238      if (!KEEP_UBERCAL) {
     239        image[nimage].Mcal = 0.0;
     240        image[nimage].dMcal = NAN;
     241        image[nimage].flags &= ~ID_IMAGE_PHOTOM_UBERCAL;
     242      }
    245243      image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR;
     244      image[nimage].ubercalDist = 1000;
    246245    }
    247246    line_number[nimage] = i;
    248247    nimage ++;
    249248    if (nimage == NIMAGE) {
    250       NIMAGE += 100;
     249      NIMAGE += D_NIMAGE;
     250      D_NIMAGE = MAX (100000, D_NIMAGE * 1.5);
    251251      REALLOCATE (image, Image, NIMAGE);
    252252      REALLOCATE (line_number, off_t, NIMAGE);
    … …  
    262262  REALLOCATE (line_number, off_t, MAX (nimage, 1));
    263263  free (skycoords);
     264  free (RmaxSky);
     265  free (index);
    264266
    265267  *Nimage  = nimage;

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

    r31450 r33651  
    66// and uses the bits read from disk as the test
    77
    8 int setExclusions (Catalog *catalog, int Ncatalog) {
     8int setExclusions (Catalog *catalog, int Ncatalog, int verbose) {
    99
    1010  off_t i, j, k, m, Narea, Nnocal, Ngood;
    … …  
    5959    }
    6060  }
    61   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements marked by area\n",    Narea);
    62   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements marked nocal\n",      Nnocal);
    63   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measurements kept for analysis\n", Ngood);
     61  if (verbose) fprintf (stderr, OFF_T_FMT" measurements marked by area\n",    Narea);
     62  if (verbose) fprintf (stderr, OFF_T_FMT" measurements marked nocal\n",      Nnocal);
     63  if (verbose) fprintf (stderr, OFF_T_FMT" measurements kept for analysis\n", Ngood);
    6464  return (TRUE);
    6565}

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

    r31450 r33651  
    55// output dbFlags values
    66
    7 void setMrelFinal (Catalog *catalog) {
    8 
    9   off_t i, j, m;
     7void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
     8
     9  off_t i;
    1010  int ecode;
     11
     12  int Nsecfilt = GetPhotcodeNsecfilt ();
    1113
    1214  /* if we reset the catalog, reset all the current measurements */
    1315  if (RESET) {
    14 
    15     int Nsecfilt = GetPhotcodeNsecfilt ();
    16 
     16    // flags used by the photometry analysis (excluding UBERCAL)
     17    unsigned int PHOTOM_FLAGS =
     18      ID_MEAS_NOCAL |       // detection ignored for this analysis (photcode, time range)
     19      ID_MEAS_POOR_PHOTOM | // detection is photometry outlier
     20      ID_MEAS_SKIP_PHOTOM | // detection was ignored for photometry measurement
     21      ID_MEAS_AREA;       // detetion was outside acceptable area of device
     22   
     23    // flags used by the photometry analysis (excluding UBERCAL)
     24    // unsigned int secfiltFlags =
     25    //   ID_PHOTOM_PASS_0 | // average measured at pass 0
     26    //   ID_PHOTOM_PASS_1 | // average measured at pass 1
     27    //   ID_PHOTOM_PASS_2 | // average measured at pass 2
     28    //   ID_PHOTOM_PASS_3 | // average measured at pass 3
     29    //   ID_PHOTOM_PASS_4 | // average measured at pass 3
     30    //   ID_SECF_USE_SYNTH | // average measured at pass 3
     31    //   ID_SECF_USE_UBERCAL | // average measured at pass 3
     32    //   ID_SECF_OBJ_EXT; // average measured at pass 3
     33   
    1734    int Ns;
    1835    for (Ns = 0; Ns < Nphotcodes; Ns++) {
    … …  
    2239
    2340      for (i = 0; i < catalog[0].Naverage; i++) {
    24         catalog[0].secfilt[Nsecfilt*i+Nsec].M  = NAN;
    25         catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN;
    26         catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT;
    27 
    28         m = catalog[0].average[i].measureOffset;
     41        off_t N = Nsecfilt*i+Nsec;
     42        catalog[0].secfilt[N].M      = NAN;
     43        catalog[0].secfilt[N].Map    = NAN;
     44        catalog[0].secfilt[N].dM     = NAN;
     45        catalog[0].secfilt[N].Xm     = NAN;
     46        catalog[0].secfilt[N].M_20   = NAN_S_SHORT;
     47        catalog[0].secfilt[N].M_80   = NAN_S_SHORT;
     48        catalog[0].secfilt[N].Ncode  = 0;
     49        catalog[0].secfilt[N].Nused  = 0;
     50        catalog[0].secfilt[N].flags  = 0; // XXX there are no *astrometry* bits in secfilt.flags
     51        catalog[0].secfilt[N].ubercalDist = 1000;
     52
     53        off_t m = catalog[0].average[i].measureOffset;
     54        off_t j;
    2955        for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    3056       
    … …  
    3965          }
    4066       
    41           catalog[0].measure[m].Mcal = 0;
    42           catalog[0].measure[m].dbFlags &= 0xff00;
    43           catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
    44           catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    45           catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA;
    46           catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
     67          if (!KEEP_UBERCAL) {
     68            catalog[0].measure[m].Mcal = 0;
     69            catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
     70          }
     71          catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS;
    4772        }
    4873      }
    … …  
    5176
    5277  // this sets flags in the measureT element, not the measure element
    53   setExclusions (catalog, 1);  /* mark by area */
     78  setExclusions (catalog, 1, VERBOSE);  /* mark by area */
    5479
    5580  /* set catalog[0].found[i] = FALSE */
    56   ALLOCATE (catalog[0].found, off_t, MAX (1, catalog[0].Naverage));
    57   for (i = 0; i < catalog[0].Naverage; i++) {
     81  ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
     82  for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) {
    5883    catalog[0].found[i] = FALSE;
    5984  }
    6085
     86  clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    6187  for (i = 0; i < 5; i++) {
    62     skip_measurements (catalog, i);       /* set ID_MEAS_SKIP for measures to be skipped */
    63     setMrelOutput  (catalog, 1, FALSE);   /* set Mrel using all measures not skipped */
    64     clean_measures (catalog, 1, TRUE);    /* mark outliers ID_MEAS_POOR_PHOTOM */
    65     setMrelOutput  (catalog, 1, TRUE);    /* set Mrel using remaining measures */
    66   }
    67   setMcalOutput (catalog, 1);
    68 
    69   int Nsecfilt = GetPhotcodeNsecfilt ();
     88    skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
     89    setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
     90  }
     91  setMcalOutput (catalog, 1, flatcorr);
    7092
    7193  /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
    7294  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
     95  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
    7396  for (i = 0; i < catalog[0].Naverage; i++) {
    7497    catalog[0].average[i].flags &= ~ID_STAR_FEW;
    7598    catalog[0].average[i].flags &= ~ID_STAR_POOR;
     99    off_t j;
    76100    for (j = 0; j < Nsecfilt; j++) {
    77101        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW;
    78102        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR;
    79103    }
    80     m = catalog[0].average[i].measureOffset;
     104    off_t m = catalog[0].average[i].measureOffset;
    81105    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
    82106      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
    … …  
    86110
    87111/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
    88 void skip_measurements (Catalog *catalog, int pass) {
     112// For each star & photcode, we are selecting the measurements to use in the average
     113// photometry value.  This function is called several times, lowering the bar to
     114// acceptance on each pass. The rejections at the different passes are:
     115
     116// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
     117// 1 : accept measurements thought to be POOR based on photflags
     118// 2 : accept the measurements marked as outliers; accept images marked as outliers
     119// 3 : accept measurements thought to be BAD based on photflags
     120// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)
     121
     122// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
     123// non-outliers with Imag range
     124
     125// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
     126//
     127//   off_t k;
     128//
     129//   int Nsecfilt = GetPhotcodeNsecfilt ();
     130//
     131//   off_t m = average[0].measureOffset;
     132//
     133//   for (k = 0; k < average[0].Nmeasure; k++, m++) {
     134//     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
     135//   }
     136//
     137//   int Ns;
     138//   for (Ns = 0; Ns < Nsecfilt; Ns++) {
     139//     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
     140//   }
     141//   return 1;
     142// }
     143
     144int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
     145
     146void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {
    89147
    90148  off_t i, k, m;
    91   int ecode, found, Ns;
    92149  off_t Ntot, Ntry, Nkeep, Nskip;
    93150  float mag;
    94151
     152  int Nsecfilt = GetPhotcodeNsecfilt ();
     153
    95154  Ntot = Ntry = Nskip = Nkeep = 0;
    96155
    97156  /* allow measures from images marked POOR and FEW */
    98   if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
     157  if (pass >= 2) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
    99158 
    100159  /* allow measures marked as outliers (POOR) and off image region (AREA) */
    101   if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
    102 
    103   /* mark measures which should be ignored on second pass */
     160  if (pass >= 2) {
     161    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
     162  } else {
     163    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
     164  }
     165
     166  /* mark measures which should be ignored on this pass */
    104167  for (i = 0; i < catalog[0].Naverage; i++) {
    105168    Ntot += catalog[0].average[i].Nmeasure;
    106     if (catalog[0].found[i]) continue;
    107 
    108     m = catalog[0].average[i].measureOffset;
    109     for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
    110       Ntry++;
    111 
    112       /* clear SKIP for all measures at first */
    113       catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
    114 
    115       /** never use these measurements (wrong photcode, bad time range) */
    116       /* skipped via NOCAL, don't mark as skipped */
    117       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
    118       found = FALSE;
    119       for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
    120         if (ecode == photcodes[Ns][0].code) found = TRUE;
     169
     170    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
     171      fprintf (stderr, "test obj\n");
     172      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
     173    }
     174
     175    // mark measurements for each secfilt separately
     176    int Ns;
     177    for (Ns = 0; Ns < Nphotcodes; Ns++) {
     178     
     179      int thisCode = photcodes[Ns][0].code;
     180      int Nsec = GetPhotcodeNsec(thisCode);
     181     
     182      /* star/photcodes already calibrated */
     183      if (catalog[0].found[Nsecfilt*i+Nsec]) continue; 
     184
     185      m = catalog[0].average[i].measureOffset;
     186      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
     187        Ntry++;
     188
     189        // skip measurements not related to this photcode
     190        PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
     191        if (!code) continue;
     192        if (code->equiv != thisCode) continue;
     193
     194        // skip measurements by time range (mark as skipped for this secfilt)
     195        // XXX note that this is a bit dangerous : some objects may never get calibrated
     196        if (TimeSelect) {
     197          if (catalog[0].measure[m].t < TSTART) goto skip;
     198          if (catalog[0].measure[m].t > TSTOP) goto skip;
     199        }
     200       
     201        // we now have a measurement relavant to this photcode (and time range)
     202
     203        // clear SKIP for all measures at first
     204        catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
     205        catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
     206
     207        // skip measurements from BAD images and mosaics (not REF mags)
     208        // do NOT skip measurements without a matching image (REF mags)
     209        off_t Nim = getImageEntry (m, 0);
     210        if (Nim > -1) {
     211          if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
     212          if (isnan(getMmos (m, 0))) goto skip;
     213        }
     214
     215        // PASS 4 : skip measurements by inst mag limit (not REF mags)
     216        if ((pass < 4) && ImagSelect) {
     217          if (Nim > -1) {
     218            mag = PhotInst (&catalog[0].measure[m]);
     219            if (mag < ImagMin) goto skip;
     220            if (mag > ImagMax) goto skip;
     221          }
     222        }
     223
     224        // PASS 3 : accept bad measurements (eg, SAT, CR)
     225        if (pass < 3) {
     226          if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
     227          if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
     228            if (catalog[0].measure[m].psfQual < 0.85) goto skip;
     229          }
     230        }
     231       
     232        // PASS 2 : internal outliers accepted
     233
     234        // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
     235        if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;
     236
     237        // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)
     238
     239        Nkeep ++;
     240       
     241        continue;
     242
     243      skip:
     244        catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     245        catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
     246        Nskip ++;
    121247      }
    122       if (!found) continue;
    123 
    124       /* skip measurements by time range */
    125       if (TimeSelect) {
    126         if (catalog[0].measure[m].t < TSTART) goto skip;
    127         if (catalog[0].measure[m].t > TSTOP) goto skip;
    128       }
    129 
    130       /* skip measurements with sat. dophot values */
    131       // XXX we need to fix this: use the photFlags to ID saturation...
    132       // if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip;
    133 
    134       /* skip measurements from BAD images and mosaics */
    135       /* do NOT skip measurements without a matching image */
    136       if (isnan(getMcal (m, 0))) goto skip;
    137       if (isnan(getMmos (m, 0))) goto skip;
    138 
    139       /* skip measurements by inst mag limit */
    140       if ((pass < 4) && ImagSelect) {
    141         mag = PhotInst (&catalog[0].measure[m]);
    142         if (mag < ImagMin) goto skip;
    143         if (mag > ImagMax) goto skip;
    144       }
    145 
    146       // d1 = (catalog[0].measure[m].dophot == 1);
    147       // d2 = (catalog[0].measure[m].dophot == 2);
    148       // fprintf (stderr, "fix the dophot / type selections!");
    149       // abort ();
    150 
    151       /* skip measurements by measurement error */
    152       // if ((pass < 2) && !(d1 || d2)) goto skip;
    153 
    154       /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */
    155       // if ((pass < 1) && !d1) goto skip;
    156       Nkeep ++;
    157 
    158       continue;
    159 
    160     skip:
    161       catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    162       catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
    163       Nskip ++;
    164     }
    165   }
    166   if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
     248    }
     249  }
     250  if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
    167251                        pass, Ntot, Ntry, Nskip, Nkeep);
    168252}
Note: See TracChangeset for help on using the changeset viewer.