Index: /branches/eam_branches/20090715/Ohana/src/addstar/Makefile =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/Makefile (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/Makefile (revision 25020) @@ -58,4 +58,5 @@ $(SRC)/greference.$(ARCH).o \ $(SRC)/grefstars.$(ARCH).o \ +$(SRC)/GetZeroPointExposure.$(ARCH).o \ $(SRC)/LoadStars.$(ARCH).o \ $(SRC)/LoadHeaders.$(ARCH).o \ @@ -199,4 +200,5 @@ $(SRC)/replace_match.$(ARCH).o \ $(SRC)/update_coords.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ $(SRC)/StarOps.$(ARCH).o \ $(SRC)/ConfigInit.$(ARCH).o \ Index: /branches/eam_branches/20090715/Ohana/src/addstar/doc/notes.txt =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/doc/notes.txt (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/doc/notes.txt (revision 25020) @@ -1,2 +1,49 @@ + +2009.08.02 + + ZERO POINTS: addstar currently has a couple of inconsistent ways that it + measures or uses image zero points: + + * the standard zero point: data in the database is saved as M_inst + + 25.0. This value is written to the image header and checked when + the image is loaded. This is completely archaic, but removing + this concept from the DVO system will require some work and some + care. This value is NOT reflected in the image->Mcal values + + * on-the-fly calibration: addstar was used by CFHT skyprobe to + measure the sky transparency. This calculation was done by + measuring the zero point between the Tycho stars in the database + and the stars in the image. Star data is saved in + AddToCalibration as they are found (in find_matches, or + equivalent). At the end of the match processing, the zero point + for the image was calculated with FindCalibration and saved to the + image table but NOT the measure table. For skyprobe, the stars + are normally not saved to the database -- the Tycho stars are kept + in the measure & average tables, and used as the reference for the + zero point calculation above. + + For Pan-STARRS, the zero point analysis is performed by psastro + before we run addstar. However, there are a few ways we could + choose to apply the zero point calculation: + + * use the per-chip zero point reported in the individual chip + headers (easy, but fairly wrong) + + * calculate the per-exposure zero point from the chip headers and + supply. + + * calculate the per-exposure zero point from the MATCHED_REFS table + + We also need to consider the zero point error. We currently have + only dM (poisson error) + dMcal (calibration error). We need to + divide dMsys into dMcal (photometry calibration error) and dMap + (aperture correction error). + + * zero point options: + * NOMINAL : set Mcal = 0.0 + * CHIP_HEADER : each chip has zero point in header + * PHU_HEADER : a single exposure-wide value in the header + * CHIP_AVERAGE : determine zero point for exposure from chip zero points + * MATCHED_REFS : recalculate from table 2009.02.11 Index: /branches/eam_branches/20090715/Ohana/src/addstar/include/2mass.h =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/include/2mass.h (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/include/2mass.h (revision 25020) @@ -33,5 +33,5 @@ e_time get2mass_date (char *ptr, int Nbound, int Nmax); -int load2mass_as_rawdata (SkyTable *skytable, char *filename, AddstarClientOptions options); +int load2mass_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options); SkyTable *load2mass_acc (char *path, char *accel); int get2mass_3star (Stars *star, char *line, int Nmax); Index: /branches/eam_branches/20090715/Ohana/src/addstar/include/addstar.h =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/include/addstar.h (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/include/addstar.h (revision 25020) @@ -73,4 +73,5 @@ char TYCHO_DIR[256]; char SubpixDatafile[256]; +char *USE_NAME; char PASSWORD[80]; char HOSTNAME[80]; @@ -118,4 +119,7 @@ double FAKE_THETA; // boresite angle for fake images +char ZERO_POINT_OPTION[64]; +float ZERO_POINT_OFFSET; + // carries the mosaic into gstars @@ -187,4 +191,5 @@ HeaderSet *MatchHeaders PROTO((int **extsize, int *nimage, int mode, Header **headers, int Nheaders)); int LoadData PROTO((FILE *f, char *file, Image **images, int *nvalid, Stars **stars, int *Nstars, Header **headers, int *extsize, HeaderSet *headerSets, int NheaderSets)); +int GetZeroPointExposure PROTO((Header **headers, HeaderSet *headerSets, int Nimages)); int in_image PROTO((double r, double d, Image *image)); Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/ConfigInit.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/ConfigInit.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/ConfigInit.c (revision 25020) @@ -77,4 +77,9 @@ if (!ScanConfig (config, "EXTNAME-KEYWORD", "%s", 0, ExtnameKeyword)) { strcpy (ExtnameKeyword, "EXTNAME"); + } + + // if this config variable is not set, do not set any zero point offset + if (!ScanConfig (config, "ZERO_POINT_OPTION", "%s", 0, ZERO_POINT_OPTION)) { + strcpy (ZERO_POINT_OPTION, "NOMINAL"); } Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/GetZeroPointExposure.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/GetZeroPointExposure.c (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/GetZeroPointExposure.c (revision 25020) @@ -0,0 +1,96 @@ +# include "addstar.h" + +// XXX this needs to handle the nominal zero point for the photcode(s) and the airmass term +int GetZeroPointExposure (Header **headers, HeaderSet *headerSets, int Nimages) { + + if (!strcasecmp(ZERO_POINT_OPTION, "NOMINAL")) { + ZERO_POINT_OFFSET = 0.0; + return (TRUE); + } + + if (!strcasecmp(ZERO_POINT_OPTION, "CHIP_HEADER")) { + ZERO_POINT_OFFSET = 0.0; + return (TRUE); + } + + if (!strcasecmp(ZERO_POINT_OPTION, "CHIP_AVERAGE")) { + int i, Nzpt, Nmid, Nhead; + float *zpt, ZPT_OBS; + PhotCode *photcode; + char photname[80]; + + Nzpt = 0; + ALLOCATE (zpt, float, Nimages); + + for (i = 0; i < Nimages; i++) { + if (!strcmp(headerSets[i].exthead, "PHU")) continue; + Nhead = headerSets[i].extnum_head; + + if (!gfits_scan (headers[Nhead], "ZPT_OBS", "%f", 1, &ZPT_OBS)) { + fprintf (stderr, "zero point not supplied in header\n"); + continue; + } + + /* get photcode from header */ + if (!gfits_scan (headers[Nhead], "PHOTCODE", "%s", 1, photname)) { + fprintf (stderr, "photcode not supplied in header\n"); + continue; + } + photcode = GetPhotcodebyName (photname); + if (photcode == NULL) { + fprintf (stderr, "photcode %s not found in photcode table\n", photname); + continue; + } + // photcode table stores zero point in milli-mags : + // Mrel = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal; + // ZERO_POINT above is the standard 25.0 + // measure[0].M = Mobs + 2.5log(exptime) + // Mrel should be equivalent to Mref. thus Mref - Mobs - 2.5log(exptime) = -Mcal + + zpt[Nzpt] = 0.001*photcode[0].C - ZPT_OBS; + Nzpt ++; + } + + if (Nzpt == 0) { + fprintf (stderr, "WARNING: zero point is not measured, no valid entries in headers\n"); + ZERO_POINT_OFFSET = 0.0; + free (zpt); + return (FALSE); + } + + fsort (zpt, Nzpt); + Nmid = Nzpt / 2; // Nzpt = 1,2,3,4 -> Nmid = 0,1,1,2 2: (0 1), 4: 0 (1 2) 3 + ZERO_POINT_OFFSET = (Nzpt % 2) ? zpt[Nmid] : 0.5*(zpt[Nmid] + zpt[Nmid-1]); + free (zpt); + return (TRUE); + } + + // this option is not consistent with the options for photcodes + if (!strcasecmp(ZERO_POINT_OPTION, "PHU_HEADER")) { + int i, Nhead; + float ZPT_OBS; + + ZERO_POINT_OFFSET = 0.0; + for (i = 0; i < Nimages; i++) { + if (strcmp(headerSets[i].exthead, "PHU")) continue; + Nhead = headerSets[i].extnum_head; + + if (!gfits_scan (headers[Nhead], "ZPT_OBS", "%f", 1, &ZPT_OBS)) { + fprintf (stderr, "WARNING: zero point is not measured, no valid entries in headers\n"); + return (FALSE); + } + ZERO_POINT_OFFSET = ZPT_OBS; + return (TRUE); + } + fprintf (stderr, "WARNING: zero point is not measured, no valid entries in headers\n"); + return (FALSE); + } + + if (!strcasecmp(ZERO_POINT_OPTION, "MATCHED_REFS")) { + fprintf (stderr, "ERROR: MATCHED_REFS zero point analysis is not implemented\n"); + exit (1); + } + + fprintf (stderr, "ERROR: invalid value for ZERO_POINT_OPTION: %s\n", ZERO_POINT_OPTION); + exit (1); +} Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/LoadData.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/LoadData.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/LoadData.c (revision 25020) @@ -23,5 +23,13 @@ // find image rootname - name = filebasename (file); + if (USE_NAME) { + name = filebasename (USE_NAME); + } else { + name = filebasename (file); + } + + // if zero points are calculated for the full exposure using more than just the matched chip header, + // we need to perform that analysis here + GetZeroPointExposure (headers, headerSets, Nimages); // now run through the images, interpret the headers and read the stars Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/ReadImageHeader.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/ReadImageHeader.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/ReadImageHeader.c (revision 25020) @@ -8,4 +8,5 @@ double tmp, sec, Cerror, ZeroPt, FWHM_X, FWHM_Y; char *c, photname[64], line[80]; + PhotCode *photcodeData = NULL; // zero out the entire image structure @@ -104,9 +105,10 @@ return (FALSE); } - photcode = GetPhotcodeCodebyName (photname); - if (photcode == 0) { + photcodeData = GetPhotcodebyName (photname); + if (photcodeData == NULL) { fprintf (stderr, "photcode %s not found in photcode table\n", photname); return (FALSE); } + photcode = photcodeData[0].code; } if (photcode == 0) { @@ -191,4 +193,8 @@ } + // XXX this is archaic: we used to set a fixed zero point of 25 to shift data into the range + // 0 - 32 so it would fit in an unsigned int. This is also needed because some programs like + // sextractor will put in an arbitrary zero point that we need to understand to get back to + // instrumental mags. gfits_scan (header, "ZERO_PT", "%lf", 1, &ZeroPt); if (ZeroPt != GetZeroPoint()) { @@ -197,6 +203,17 @@ } + // if it exists, find the + if (!strcasecmp(ZERO_POINT_OPTION, "CHIP_HEADER")) { + float ZPT_OBS; + if (!photcodeData || !gfits_scan (header, "ZPT_OBS", "%f", 1, &ZPT_OBS)) { + fprintf (stderr, "zero point not supplied in header\n"); + ZERO_POINT_OFFSET = 0.0; + } else { + ZERO_POINT_OFFSET = 0.001*photcodeData[0].C - ZPT_OBS; + } + } + /* secz is in units milli-airmass */ - image[0].Mcal = 0.0; + image[0].Mcal = ZERO_POINT_OFFSET; image[0].Xm = NAN_S_SHORT; image[0].flags = 0; Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/SkyRegionUtils.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/SkyRegionUtils.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/SkyRegionUtils.c (revision 25020) @@ -11,4 +11,5 @@ ALLOCATE (subset, SkyList, 1); ALLOCATE (subset[0].regions, SkyRegion *, NSUBSET); + ALLOCATE (subset[0].filename, char *, NSUBSET); subset[0].ownElements = FALSE; // free these elements when freeing the list @@ -20,6 +21,11 @@ subset[0].regions[Nsubset] = input[0].regions[i]; + subset[0].filename[Nsubset] = input[0].filename[i]; Nsubset ++; - CHECK_REALLOCATE (subset[0].regions, SkyRegion *, NSUBSET, Nsubset, 100); + if (Nsubset >= NSUBSET) { + NSUBSET += 100; + REALLOCATE (subset[0].regions, SkyRegion *, NSUBSET); + REALLOCATE (subset[0].filename, char *, NSUBSET); + } subset[0].Nregions = Nsubset; } Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/args.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/args.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/args.c (revision 25020) @@ -46,4 +46,11 @@ if ((N = get_argument (argc, argv, "-list"))) { options.filelist = TRUE; + remove_argument (N, &argc, argv); + } + + USE_NAME = NULL; + if ((N = get_argument (argc, argv, "-use-name"))) { + remove_argument (N, &argc, argv); + USE_NAME = strcreate (argv[N]); remove_argument (N, &argc, argv); } Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/args_load2mass.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/args_load2mass.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/args_load2mass.c (revision 25020) @@ -38,4 +38,10 @@ } + /* only add to existing regions */ + options.existing_regions = FALSE; + if ((N = get_argument (argc, argv, "-existing-regions"))) { + options.existing_regions = TRUE; + remove_argument (N, &argc, argv); + } /* only add to existing objects */ options.only_match = FALSE; @@ -62,5 +68,4 @@ options.timeref = 0; options.mosaic = FALSE; - options.existing_regions = FALSE; options.skip_missed = FALSE; options.closest = FALSE; Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass.c (revision 25020) @@ -7,4 +7,6 @@ int i; SkyTable *sky, *sky2mass; + SkyList *skylist = NULL; + SkyList *overlap = NULL; AddstarClientOptions options; @@ -13,7 +15,19 @@ options = args_load2mass (argc, argv, options); + // load the full sky description table: sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); + // generate the subset matching the user-selected region + skylist = SkyListByPatch (sky, -1, &UserPatch); + + // if we only match to existing (already populated) regions, limit the select to those regions: + if (options.existing_regions) { + SkyList *tmp; + tmp = SkyListExistingSubset (skylist, CATDIR); + SkyListFree (skylist); + skylist = tmp; + } + path = TWO_MASS_DIR_AS; @@ -23,6 +37,14 @@ for (i = 0; i < sky2mass[0].Nregions; i++) { + // check if any of the skylist entries overlap this 2mass catalog: + overlap = SkyListByBounds_List (skylist, -1, sky2mass[0].regions[i].Rmin, sky2mass[0].regions[i].Rmax, sky2mass[0].regions[i].Dmin, sky2mass[0].regions[i].Dmax); + if (overlap[0].Nregions == 0) { + SkyListFree (overlap); + continue; + } + fprintf (stderr, "loading %s\n", sky2mass[0].filename[i]); - load2mass_as_rawdata (sky, sky2mass[0].filename[i], options); + load2mass_as_rawdata (overlap, sky2mass[0].filename[i], options); + SkyListFree (overlap); } exit (0); Index: /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass_as_rawdata.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass_as_rawdata.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/addstar/src/load2mass_as_rawdata.c (revision 25020) @@ -10,5 +10,5 @@ # define DEBUG 0 -int load2mass_as_rawdata (SkyTable *skytable, char *filename, AddstarClientOptions options) { +int load2mass_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) { int i, j, verbose; @@ -109,4 +109,13 @@ if (tstars[i].D > UserPatch.Dmax) continue; + // identify the relevant catalog + skylist = SkyRegionByPoint_List (skytable, -1, tstars[i].R, tstars[i].D); + if (skylist[0].Nregions == 0) { + SkyListFree (skylist); + continue; + } + region = skylist[0].regions[0]; + if (DEBUG) fprintf (stderr, "writing to %s\n", skylist[0].filename[0]); + // collect array of (Stars *) stars in a new output catalog Nstars = 0; @@ -114,9 +123,5 @@ ALLOCATE (stars, Stars *, NSTARS); - // identify the relevant catalog - skylist = SkyRegionByPoint (skytable, -1, tstars[i].R, tstars[i].D); - region = skylist[0].regions[0]; - if (DEBUG) fprintf (stderr, "writing to %s\n", skylist[0].filename[0]); - + // loop over stars in this 2mass region that are also in this output region for (j = i; j < Ntstars; j++) { if (tstars[j].flag) continue; @@ -166,5 +171,5 @@ } - if (DEBUG) fprintf (stderr, "selected %d stars (%10.6f - %10.6f, %10.6f - %10.6f\n", Nstars, + if (DEBUG) fprintf (stderr, "selected %d stars (%10.6f - %10.6f, %10.6f - %10.6f)\n", Nstars, region[0].Rmin, region[0].Rmax, region[0].Dmin, region[0].Dmax); Index: /branches/eam_branches/20090715/Ohana/src/getstar/Makefile =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/Makefile (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/getstar/Makefile (revision 25020) @@ -23,9 +23,12 @@ dvoImageOverlaps.install: $(DESTBIN)/dvoImageOverlaps +dvoImagesAtCoords: $(BIN)/dvoImagesAtCoords.$(ARCH) +dvoImagesAtCoords.install: $(DESTBIN)/dvoImagesAtCoords + dvoImageExtract: $(BIN)/dvoImageExtract.$(ARCH) dvoImageExtract.install: $(DESTBIN)/dvoImageExtract -all: getstar dvoImageOverlaps dvoImageExtract -install: getstar.install dvoImageOverlaps.install dvoImageExtract.install +all: getstar dvoImageOverlaps dvoImageExtract dvoImagesAtCoords +install: getstar.install dvoImageOverlaps.install dvoImageExtract.install dvoImagesAtCoords.install GETSTAR = \ @@ -61,4 +64,18 @@ $(BIN)/dvoImageOverlaps.$(ARCH): $(OVERLAPS) +COVERLAPS = \ +$(SRC)/dvoImagesAtCoords.$(ARCH).o \ +$(SRC)/args_coords.$(ARCH).o \ +$(SRC)/ConfigInit_coords.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/ReadImageFiles.$(ARCH).o \ +$(SRC)/ReadImageHeader.$(ARCH).o \ +$(SRC)/GetFileMode.$(ARCH).o \ +$(SRC)/MatchCoords.$(ARCH).o + +$(COVERLAPS): $(INC)/dvoImagesAtCoords.h $(DVO_INCS) +$(BIN)/dvoImagesAtCoords.$(ARCH): $(COVERLAPS) + EXTRACT = \ $(SRC)/dvoImageExtract.$(ARCH).o \ Index: /branches/eam_branches/20090715/Ohana/src/getstar/include/dvoImagesAtCoords.h =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/include/dvoImagesAtCoords.h (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/getstar/include/dvoImagesAtCoords.h (revision 25020) @@ -0,0 +1,39 @@ +# include +# include +# include +# include + +enum {NONE, SIMPLE_CMP, SIMPLE_CMF, SIMPLE_MEF, MOSAIC_CMP, MOSAIC_CMF, MOSAIC_MEF, MOSAIC_PHU}; + +int VERBOSE; + +char OUTPUT[256]; +char GSCFILE[256]; +char CATDIR[256]; +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char SKY_TABLE[256]; +int SKY_DEPTH; +char ImageCat[256]; +Coords *MOSAIC; // carries the mosaic into ReadImageHeader + +int WITH_PHU; +int SOLO_PHU; +int ACCEPT_ASTROM; +char *astromFile; + +// These functions belong to dvoImageOverap +int args_coords PROTO((int argc, char **argv)); +int ConfigInit_coords PROTO((int *argc, char **argv)); +int Shutdown PROTO((char *format, ...)); +void TrapSignal PROTO((int sig)); +void SetProtect PROTO((int mode)); +int SetSignals PROTO((void)); + +Image *ReadImageFiles (char *filename, int *Nimages); +int ReadImageHeader (Header *header, Image *image); +int *MatchCoords(Image *, int, double, double, int *); + +int GetFileMode (Header *header); +// int edge_check (double *x1, double *y1, double *x2, double *y2); +// double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); Index: /branches/eam_branches/20090715/Ohana/src/getstar/src/ConfigInit_coords.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/src/ConfigInit_coords.c (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/getstar/src/ConfigInit_coords.c (revision 25020) @@ -0,0 +1,46 @@ +# include "dvoImagesAtCoords.h" + +int ConfigInit_coords (int *argc, char **argv) { + + char *config, *file; + char CatdirPhotcodeFile[256]; + char MasterPhotcodeFile[256]; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (1); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "CATMODE", "%s", 0, CATMODE); + ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { + SKY_DEPTH = 2; + } + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { + SKY_TABLE[0] = 0; + } + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + /* XXX this does not yet write out the master photcode table */ + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) { + fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile); + exit (1); + } + + free (config); + free (file); + + return (TRUE); +} Index: /branches/eam_branches/20090715/Ohana/src/getstar/src/MatchCoords.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/src/MatchCoords.c (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/getstar/src/MatchCoords.c (revision 25020) @@ -0,0 +1,86 @@ +# include "dvoImagesAtCoords.h" +# ifndef FLT_MAX +# define FLT_MAX 1e32 +# endif + +/* given coordinate, find images in list that contain the point */ +int *MatchCoords (Image *dbImages, int NdbImages, double ra, double dec, int *Nmatch) { + + int i, j, N, addtolist, status; + int NMATCH, nmatch, *match; + Coords tcoords; + double r, d; + double Xi[4], Yi[4], Xo[4], Yo[4]; /* image and original corners */ + double Xmin, Xmax, Ymin, Ymax; + double xmin, xmax, ymin, ymax; + + *Nmatch = 0; + + /* match represents the subset of overlapping images */ + nmatch = 0; + NMATCH = 20; + ALLOCATE (match, int, NMATCH); + + /* setup links for mosaic WRP and DIS entries */ + BuildChipMatch (dbImages, NdbImages); + + /* run through image table and search for overlaps + also define the vtable entries for the images we keep */ + + for (i = 0; i < NdbImages; i++) { + + if (!WITH_PHU && !strcmp (&dbImages[i].coords.ctype[4], "-DIS")) continue; + if ( SOLO_PHU && strcmp (&dbImages[i].coords.ctype[4], "-DIS")) continue; + + /* if any of these images are WRP images, need to find matching DIS */ + if (!FindMosaicForImage (dbImages, NdbImages, i)) continue; + + /* define image corners */ + SetImageCorners (Xi, Yi, &dbImages[i]); + // Xi[4] = Xi[0]; Yi[4] = Yi[0]; + + + ymin = xmin = +FLT_MAX; + ymax = xmax = -FLT_MAX; + for (j = 0; j < 4; j++) { + xmin = MIN(xmin, Xi[j]); + xmax = MAX(xmax, Xi[j]); + ymin = MIN(ymin, Yi[j]); + ymax = MAX(ymax, Yi[j]); + } + + // transform input point to image coords + double x, y; + RD_to_XY(&x, &y, ra, dec, &dbImages[i].coords); + + if ((x >= xmin && x <= xmax) && (y >= ymin && y <= ymax)) { + + match[nmatch] = i; + nmatch ++; + if (nmatch == NMATCH) { + NMATCH += 20; + REALLOCATE (match, int, NMATCH); + } + } + + } + if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", nmatch); + + *Nmatch = nmatch; + return (match); +} + +void SetImageCorners (double *X, double *Y, Image *image) { + + if (!strcmp(&image[0].coords.ctype[4], "-DIS")) { + X[0] = -0.5*image[0].NX; Y[0] = -0.5*image[0].NY; + X[1] = +0.5*image[0].NX; Y[1] = -0.5*image[0].NY; + X[2] = +0.5*image[0].NX; Y[2] = +0.5*image[0].NY; + X[3] = -0.5*image[0].NX; Y[3] = +0.5*image[0].NY; + } else { + X[0] = 0; Y[0] = 0; + X[1] = image[0].NX; Y[1] = 0; + X[2] = image[0].NX; Y[2] = image[0].NY; + X[3] = 0; Y[3] = image[0].NY; + } +} Index: /branches/eam_branches/20090715/Ohana/src/getstar/src/args_coords.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/src/args_coords.c (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/getstar/src/args_coords.c (revision 25020) @@ -0,0 +1,86 @@ +# include "dvoImagesAtCoords.h" + +void help () { + fprintf (stderr, "USAGE: \n" + "dvoImagesAtCoords [-astrom astrom_file] (coords_file)\n" + ); + exit (2); +} + +int args_coords (int argc, char **argv) { + + int N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + astromFile = NULL; + if ((N = get_argument(argc, argv, "-astrom"))) { + remove_argument (N, &argc, argv); + if (argv[N] == NULL) { + fprintf (stderr, "ERROR: no file provided with -astrom\n"); + exit (1); + } + astromFile = strdup(argv[N]); + remove_argument (N, &argc, argv); + } + + /* check for command line options */ + WITH_PHU = FALSE; + if ((N = get_argument (argc, argv, "+phu"))) { + WITH_PHU = TRUE; + remove_argument (N, &argc, argv); + } + SOLO_PHU = FALSE; + if ((N = get_argument (argc, argv, "-phu"))) { + WITH_PHU = TRUE; + SOLO_PHU = TRUE; + remove_argument (N, &argc, argv); + } + + /* check for command line options */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + /* accept bad header astrometry */ + ACCEPT_ASTROM = FALSE; + if ((N = get_argument (argc, argv, "-accept"))) { + ACCEPT_ASTROM = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-accept-astrom"))) { + ACCEPT_ASTROM = TRUE; + remove_argument (N, &argc, argv); + } + /* provide a mosaic for distortion */ + MOSAIC = NULL; + if ((N = get_argument (argc, argv, "-mosaic"))) { + Header header; + ALLOCATE (MOSAIC, Coords, 1); + + remove_argument (N, &argc, argv); + if (!gfits_read_header (argv[N], &header)) { + fprintf (stderr, "ERROR: can't read header for mosaic %s\n", argv[N]); + exit (1); + } + if (!GetCoords (MOSAIC, &header)) { + fprintf (stderr, "ERROR: no astrometric solution in header\n"); + exit (1); + } + if (strcmp(&MOSAIC[0].ctype[4], "-DIS")) { + fprintf (stderr, "ERROR: not a mosaic distortion header\n"); + exit (1); + } + remove_argument (N, &argc, argv); + gfits_free_header (&header); + } + + if (argc != 2) help(); + + return (TRUE); +} Index: /branches/eam_branches/20090715/Ohana/src/getstar/src/dvoImagesAtCoords.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/getstar/src/dvoImagesAtCoords.c (revision 25020) +++ /branches/eam_branches/20090715/Ohana/src/getstar/src/dvoImagesAtCoords.c (revision 25020) @@ -0,0 +1,125 @@ +# include "dvoImagesAtCoords.h" + +typedef struct { + int id; + double ra; + double dec; +} Point; + +static int readPoints(char *filename, Point **pointsOut); +static int ListImagesAtCoords (Image *dbImages, int id, double ra, double dec, int *matches, int Nmatches); + +int main (int argc, char **argv) { + + int i, Nimages, NdbImages, status; + int Nmatches, *matches; + Image *dbImages; + int Npoints; + FITS_DB db; + + SetSignals (); + ConfigInit_coords (&argc, argv); + args_coords (argc, argv); + + Point *points; + Npoints = readPoints(argv[1], &points); + if (!Npoints) { + exit(0); + } + + if (astromFile) { + dbImages = ReadImageFiles(astromFile, &NdbImages); + } else { + /*** update the image table ***/ + /* setup image table format and lock */ + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, ImageCat, 3600.0, LCK_SOFT); // shorter timeout? + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + + /* load or create the image table */ + if (db.dbstate == LCK_EMPTY) { + fprintf (stderr, "no images in database (%s)\n", ImageCat); + exit (1); + } else { + if (!dvo_image_load (&db, VERBOSE, FALSE)) { + Shutdown ("can't read image catalog %s", db.filename); + } + } + dvo_image_unlock (&db); + + // convert database table to internal structure + dbImages = gfits_table_get_Image (&db.ftable, &NdbImages, &db.swapped); + } + + Point *pt; + for (i = 0, pt = points; i < Npoints; i++, pt++) { + matches = MatchCoords (dbImages, NdbImages, pt->ra, pt->dec, &Nmatches); + ListImagesAtCoords(dbImages, pt->id, pt->ra, pt->dec, matches, Nmatches); + } + + exit (0); +} + +static int readPoints(char *filename, Point **pointsOut) +{ + FILE *f = fopen(filename, "r"); + if (!f) { + fprintf (stderr, "failed to open coordinate file %s", filename); + exit(1); + } + + int LEN = 20; + Point *pts; + ALLOCATE(pts, Point, LEN); + int Npoints = 0; + + int Nread; + // File Format is simple text file. + // Each Line has 3 values separated by spaces that represent the coordinates to match. + // + // 'unique integer id' 'ra in degrees' 'declination in degrees' + // we don't check the id for uniqueness (garbage in garbage out). + + while ((Nread = fscanf(f, "%d %lf %lf\n", &pts[Npoints].id, &pts[Npoints].ra, &pts[Npoints].dec)) == 3) { + Npoints++; + if (Npoints >= LEN) { + LEN += 20; + REALLOCATE(pts, Point, LEN); + } + } + if (Nread != -1) { + fprintf (stderr, "unexpected data on line %d of points file %s", Npoints, filename); + exit(1); + } + + *pointsOut = pts; + + return Npoints; +} + +static int ListImagesAtCoords (Image *dbImages, int id, double ra, double dec, int *matches, int Nmatches) +{ + int i; + for (i = 0; i < Nmatches; i++) { + int N = matches[i]; + + char *name; + // output of lookup is filename[class_id.hdr] for astrometry files + char *left_bracket = rindex(dbImages[N].name, '['); + if (left_bracket) { + name = strdup(left_bracket + 1); + // zap the .hdr] + char *dot = index(name, '.'); + if (dot) { + *dot = 0; + } + } else { + name = dbImages[N].name; + } + fprintf (stdout, "%d %lf %lf %s\n", id, ra, dec, name); + } + + return (TRUE); +} + Index: /branches/eam_branches/20090715/Ohana/src/libdvo/include/dvo.h =================================================================== --- /branches/eam_branches/20090715/Ohana/src/libdvo/include/dvo.h (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/libdvo/include/dvo.h (revision 25020) @@ -311,4 +311,5 @@ float PhotInst (Measure *measure); float PhotCat (Measure *measure); +float PhotAper (Measure *measure); float PhotSys (Measure *measure, Average *average, SecFilt *secfilt); float PhotRel (Measure *measure, Average *average, SecFilt *secfilt); @@ -441,4 +442,5 @@ SkyList *SkyListByBounds PROTO((SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); SkyList *SkyListChildrenByBounds PROTO((SkyTable *table, int No, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); + int SkyListMerge PROTO((SkyList **outlist, SkyList *newlist)); int SkyListFree PROTO((SkyList *list)); @@ -447,4 +449,8 @@ int SkyTableSetFilenames PROTO((SkyTable *sky, char *path, char *ext)); +SkyList *SkyRegionByPoint_List PROTO((SkyList *inList, int depth, double ra, double dec)); +SkyList *SkyListByBounds_List PROTO((SkyList *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); +SkyList *SkyListChildrenByBounds_List PROTO((SkyList *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax)); + /* dvo-specific sorting functions */ void sortave (Average *ave, int N); Index: /branches/eam_branches/20090715/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 25020) @@ -249,4 +249,23 @@ code = &photcodes[0].code[Np]; Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; + + return (Mcat); +} + +float PhotAper (Measure *measure) { + + int Np; + float Mcat; + PhotCode *code; + + Np = photcodes[0].hashcode[measure[0].photcode]; + if (Np == -1) return (NAN); + + if (photcodes[0].code[Np].type == PHOT_REF) { + Mcat = measure[0].Map; + return (Mcat); + } + code = &photcodes[0].code[Np]; + Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C; return (Mcat); Index: /branches/eam_branches/20090715/Ohana/src/libdvo/src/skyregion_ops.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/libdvo/src/skyregion_ops.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/libdvo/src/skyregion_ops.c (revision 25020) @@ -300,4 +300,128 @@ REALLOCATE (list[0].filename, char *, NNEW); } + } + } + + list[0].Nregions = Nnew; + return (list); +} + +/* find region which overlaps c at given depth (-1 : populated ) */ +SkyList *SkyRegionByPoint_List (SkyList *inList, int depth, double ra, double dec) { + + int i; + SkyRegion **region; + SkyList *list; + + ALLOCATE (list, SkyList, 1); + ALLOCATE (list[0].regions, SkyRegion *, 1); + ALLOCATE (list[0].filename, char *, 1); + list[0].Nregions = 0; + list[0].ownElements = FALSE; // this list is only holding a view to the elements + + region = inList[0].regions; + + for (i = 0; i < inList[0].Nregions; i++) { + + // once we pass our desired depth, quit + if ((depth > -1) && (depth < region[i][0].depth)) break; + + // skip tables that do not overlap + if (ra < region[i][0].Rmin) continue; + if (ra > region[i][0].Rmax) continue; + if (dec < region[i][0].Dmin) continue; + if (dec > region[i][0].Dmax) continue; + + // skip tables that are not at our depth + if ((depth > -1) && (depth > region[i][0].depth)) continue; + if ((depth == -1) && (region[i][0].table == FALSE)) continue; + + list[0].regions[0] = region[i]; + list[0].filename[0] = inList[0].filename[i]; + list[0].Nregions = 1; + break; + } + return (list); +} + +SkyList *SkyListByBounds_List (SkyList *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax) { + + int i, j, Ns; + SkyList *list, *extra; + + Rmin = ohana_normalize_angle (Rmin); + Rmax = ohana_normalize_angle (Rmax); + + /* handle 0,360 boundary requests */ + /* this is probably wrong: I will get duplicates for all Dec bands... */ + if (Rmin > Rmax) { + list = SkyListChildrenByBounds_List (table, depth, 0.0, Rmax, Dmin, Dmax); + + extra = SkyListChildrenByBounds_List (table, depth, Rmin, 360.0, Dmin, Dmax); + REALLOCATE (list[0].regions, SkyRegion *, list[0].Nregions + extra[0].Nregions); + REALLOCATE (list[0].filename, char *, list[0].Nregions + extra[0].Nregions); + Ns = list[0].Nregions; + for (i = 0; i < extra[0].Nregions; i++) { + // search for pre-existing match + for (j = 0; j < list[0].Nregions; j++) { + if (list[0].regions[j] == extra[0].regions[i]) { + goto skip; + } + } + list[0].regions[Ns] = extra[0].regions[i]; + list[0].filename[Ns] = extra[0].filename[i]; + Ns ++; + skip: + continue; + } + list[0].Nregions = Ns; + SkyListFree (extra); + } else { + list = SkyListChildrenByBounds_List (table, depth, Rmin, Rmax, Dmin, Dmax); + } + + return (list); +} + +SkyList *SkyListChildrenByBounds_List (SkyList *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax) { + + int i, j, Ns, Ne, Nnew, NNEW; + int append; + SkyList *children; + SkyList *list; + SkyRegion **region; + + Nnew = 0; + NNEW = 50; + ALLOCATE (list, SkyList, 1); + ALLOCATE (list[0].regions, SkyRegion *, NNEW); + ALLOCATE (list[0].filename, char *, NNEW); + list[0].ownElements = FALSE; // this list is only holding a view to the elements + + region = table[0].regions; + + // can we assume the skylist is globally sorted by depth? i think so... + for (i = 0; i < table[0].Nregions; i++) { + + // once we pass our desired depth, quit + if ((depth > -1) && (depth < region[i][0].depth)) break; + + // skip tables that do not overlap + if (Rmax <= region[i][0].Rmin) continue; + if (Rmin >= region[i][0].Rmax) continue; + if (Dmax <= region[i][0].Dmin) continue; + if (Dmin >= region[i][0].Dmax) continue; + + // skip tables that are not at our depth + if ((depth > -1) && (depth > region[i][0].depth)) continue; + if ((depth == -1) && (region[i][0].table == FALSE)) continue; + + list[0].regions[Nnew] = region[i]; + list[0].filename[Nnew] = table[0].filename[i]; + Nnew ++; + if (Nnew >= NNEW) { + NNEW += 50; + REALLOCATE (list[0].regions, SkyRegion *, NNEW); + REALLOCATE (list[0].filename, char *, NNEW); } } Index: /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbExtractMeasures.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbExtractMeasures.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbExtractMeasures.c (revision 25020) @@ -117,4 +117,7 @@ case MAG_REF: value.Flt = PhotRef (equiv, average, secfilt, measure); + break; + case MAG_APER: + value.Flt = PhotAper (measure); break; case MAG_ERR: Index: /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbFields.c =================================================================== --- /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbFields.c (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/opihi/dvo/dbFields.c (revision 25020) @@ -37,4 +37,6 @@ if (!strcasecmp (string, "ave")) return (MAG_AVE); if (!strcasecmp (string, "ref")) return (MAG_REF); + if (!strcasecmp (string, "ap")) return (MAG_APER); + if (!strcasecmp (string, "aper")) return (MAG_APER); if (!strcasecmp (string, "err")) return (MAG_ERR); if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); Index: /branches/eam_branches/20090715/Ohana/src/opihi/include/dvoshell.h =================================================================== --- /branches/eam_branches/20090715/Ohana/src/opihi/include/dvoshell.h (revision 25019) +++ /branches/eam_branches/20090715/Ohana/src/opihi/include/dvoshell.h (revision 25020) @@ -27,4 +27,5 @@ MAG_INST, MAG_CAT, + MAG_APER, MAG_SYS, MAG_REL,