Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 33652) +++ /trunk/Ohana/src/addstar/Makefile (revision 33653) @@ -24,8 +24,12 @@ sedstar : $(BIN)/sedstar.$(ARCH) load2mass : $(BIN)/load2mass.$(ARCH) +loadwise : $(BIN)/loadwise.$(ARCH) +loadsupercos: $(BIN)/loadsupercos.$(ARCH) gztest : $(BIN)/gztest.$(ARCH) mkcmf : $(BIN)/mkcmf.$(ARCH) -all: addstar sedstar load2mass skycells mkcmf +all: addstar sedstar load2mass skycells mkcmf loadwise loadsupercos + +INSTALL = addstar sedstar load2mass skycells mkcmf loadwise loadsupercos # I need to fix the client/server version of addstar now that I have dropped Stars @@ -214,4 +218,40 @@ $(SRC)/psps_ids.$(ARCH).o +LOAD-WISE = \ +$(SRC)/loadwise.$(ARCH).o \ +$(SRC)/loadwise_prelim_rawdata.$(ARCH).o \ +$(SRC)/loadwise_ops.$(ARCH).o \ +$(SRC)/loadwise_star_full.$(ARCH).o \ +$(SRC)/args_loadwise.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/find_matches_refstars.$(ARCH).o \ +$(SRC)/build_links.$(ARCH).o \ +$(SRC)/replace_match.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/update_coords.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + +LOAD-SUPERCOS = \ +$(SRC)/loadsupercos.$(ARCH).o \ +$(SRC)/loadsupercos_ops.$(ARCH).o \ +$(SRC)/loadsupercos_survey.$(ARCH).o \ +$(SRC)/loadsupercos_plates.$(ARCH).o \ +$(SRC)/loadsupercos_rawdata.$(ARCH).o \ +$(SRC)/args_loadsupercos.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o \ +$(SRC)/find_matches_refstars.$(ARCH).o \ +$(SRC)/build_links.$(ARCH).o \ +$(SRC)/replace_match.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/update_coords.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/parse_time.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + SEDSTAR = \ $(SRC)/sedstar.$(ARCH).o \ @@ -241,4 +281,5 @@ $(SKYCELLS) : $(INC)/addstar.h $(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h +$(LOAD-WISE) : $(INC)/addstar.h $(INC)/WISE.h $(MKCMF) : $(INC)/mkcmf.h @@ -250,8 +291,8 @@ $(BIN)/sedstar.$(ARCH) : $(SEDSTAR) $(BIN)/load2mass.$(ARCH) : $(LOAD-2MASS) +$(BIN)/loadwise.$(ARCH) : $(LOAD-WISE) +$(BIN)/loadsupercos.$(ARCH): $(LOAD-SUPERCOS) $(BIN)/skycells.$(ARCH) : $(SKYCELLS) $(BIN)/mkcmf.$(ARCH) : $(MKCMF) - -INSTALL = addstar sedstar load2mass skycells mkcmf # I need to fix the client/server version of addstar now that I have dropped Stars Index: /trunk/Ohana/src/addstar/include/WISE.h =================================================================== --- /trunk/Ohana/src/addstar/include/WISE.h (revision 33652) +++ /trunk/Ohana/src/addstar/include/WISE.h (revision 33653) @@ -21,35 +21,25 @@ short WISE_W1, WISE_W2, WISE_W3, WISE_W4; -SkyTable *get2mass_acc (SkyRegion *patch, char *path, char *accel); -Stars *get2mass_2DR_data (SkyRegion *region, char *filename, SkyRegion *patch, int photcode, int *nstars); -Stars *get2mass_AS_data (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars); -Stars *get2mass_AS_rawdata (SkyRegion *region, char *filename, SkyRegion *patch, int phocode, int *nstars); +AddstarClientOptions args_loadwise (int argc, char **argv, AddstarClientOptions options); -SkyTable *scan2mass_acc (char *path, char *accel); -int scan2mass_as_data (char *filename); +int loadwise_prelim_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options); -char *skipNbounds (char *line, char bound, int Nbound, int Nbyte); -e_time get2mass_time (char *ptr, int Nbound, int Nbyte); -e_time get2mass_date (char *ptr, int Nbound, int Nmax); +int getWISE_setup (); +int getWISE_coords (char *line, double *R, double *D, int Nmax); +e_time getWISE_date (char *ptr, int Nbound, int Nmax); +e_time getWISE_time (char *ptr, int Nbound, int Nmax); -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); -int load2mass_catalog (Catalog *catalog, Stars *stars, int Nstars); +int getWISE_sortStars (WISE_Stars *tstars, int Ntstars); +int loadwise_star_full (Stars **star, char *line, int Nmax); -int get2mass_setup (int photcode); -int get2mass_coords (char *line, double *R, double *D, int Nmax); -int get2mass_star (Stars *star, char *line, int Nmax); -int get2mass_3star (Stars *star, char *line, int Nmax); +char *nextWISEfield (char *line); +int setWISE_cc_flag (Stars *star, char qual); +int setWISE_ph_qual (Stars *star, char qual); +int setWISE_rd_flag (Stars *star, char qual); +int setWISE_bl_flag (Stars *star, char qual); +int setWISE_gal_flag (Stars *star, char qual); +int setWISE_mp_flag (Stars *star, char qual); +int setWISE_dup_flag (Stars *star, char qual); +int setWISE_use_flag (Stars *star, char qual); -int get2mass_3star_full (Stars **star, char *line, int Nmax); -char *next2MASSfield (char *line); -int set2MASS_ph_qual (Stars *star, char qual); -int set2MASS_rd_flag (Stars *star, char qual); -int set2MASS_cc_flag (Stars *star, char qual); -int set2MASS_bl_flag (Stars *star, char qual); -int set2MASS_gal_flag (Stars *star, char qual); -int set2MASS_mp_flag (Stars *star, char qual); -int set2MASS_dup_flag (Stars *star, char qual); -int set2MASS_use_flag (Stars *star, char qual); -int get2mass_sortStars (WISE_Stars *tstars, int Ntstars); +char *skipNbounds (char *line, char bound, int Nbound, int Nbyte); Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 33652) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 33653) @@ -28,4 +28,7 @@ fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", NAME); \ exit (0); } + +# define dCOS(A) ((double) cos ((double)RAD_DEG*A)) +# define dSIN(A) ((double) sin ((double)RAD_DEG*A)) typedef struct { Index: /trunk/Ohana/src/addstar/include/supercos.h =================================================================== --- /trunk/Ohana/src/addstar/include/supercos.h (revision 33653) +++ /trunk/Ohana/src/addstar/include/supercos.h (revision 33653) @@ -0,0 +1,85 @@ +/* header for items specific to loadsupercos.c */ + +typedef struct { + double longitude; + double latitude; + double plateScale; + int survey_id; +} Survey; + +typedef struct{ + int64_t objID; // 0 + int8_t surveyID; // 8 + int32_t plateID; // 9 + double ra; // 33 + double dec; // 41 + double xCen; // 129 + double yCen; // 137 + float aU; // 145 + float bU; // 149 + int16_t thetaU; // 154 + int8_t class; // 165 + int32_t quality; // 104 + float prfStat; // 208 + float prfMag; // 212 + float gMag; // 216 + float sMag; // 220 +} Detection; + +typedef struct{ + int64_t objID; // 0 + int8_t surveyID; // 8 + int32_t plateID; // 9 + int64_t parentID; // 13 + int64_t sourceID; // 21 + int32_t recNum; // 29 + double ra; // 33 + double dec; // 41 + int64_t htmId; // 49 + double cx; // 57 + double cy; // 65 + double cz; // 77 + double xmin; // 81 + double xmax; // 89 + double ymin; // 97 + double ymax; // 105 + int32_t area; // 113 + float ipeak; // 117 + float cosmag; // 121 + float isky; // 125 + double xCen; // 129 + double yCen; // 137 + float aU; // 145 + float bU; // 149 + int16_t thetaU; // 154 + float aI; // 155 + float bI; // 159 + int16_t thetaI; // 163 + int8_t class; // 165 + int16_t pa; // 166 + int32_t ap1; // 168 + int32_t ap2; // 172 + int32_t ap3; // 176 + int32_t ap4; // 180 + int32_t ap5; // 184 + int32_t ap6; // 188 + int32_t ap7; // 192 + int32_t ap8; // 196 + int32_t blend; // 200 + int32_t quality; // 204 + float prfStat; // 208 + float prfMag; // 212 + float gMag; // 216 + float sMag; // 220 + int16_t SSAfield; // 224 + int16_t seam; // 226 (228 total) +} FullDetection; + +AddstarClientOptions args_loadsupercos (int *argc, char **argv, AddstarClientOptions options); +Survey *loadsupercos_survey (char *filename, int *nsurvey); +Image *loadsupercos_plates (Survey *survey, int Nsurvey, char *filename, int *nimage); +int loadsupercos_getFilterInfo (char *line, char *emulsion, char *filterID); +int loadsupercos_rawdata (Image *image, int *imlist, int Nimage, SkyList *skytable, char *filename, AddstarClientOptions options); +int loadsupercos_sortStars (Stars *tstars, int Ntstars); +int *loadsupercos_image_index (Image *image, int Nimage); +int loadsupercos_getST (char *line, double *st); Index: /trunk/Ohana/src/addstar/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/addstar/src/ConfigInit.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/ConfigInit.c (revision 33653) @@ -213,7 +213,6 @@ /* get detection filtering mask */ - if(!ScanConfig (config, "DETECTIONFILTER", "%d", 0, &options.detectionFilter)) { + if (!ScanConfig (config, "DETECTIONFILTER", "%d", 0, &options.detectionFilter)) { options.detectionFilter = 0; - fprintf (stderr, "Could not find 'DETECTIONFILTER' in config, using '%x'\n", options.detectionFilter); } Index: /trunk/Ohana/src/addstar/src/ReadStarsSDSS.c =================================================================== --- /trunk/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 33653) @@ -4,6 +4,4 @@ # define NFILTER 5 -# define dCOS(A) ((double) cos ((double)RAD_DEG*A)) -# define dSIN(A) ((double) sin ((double)RAD_DEG*A)) # define GET_COLUMN_5(NAME,TYPE) \ @@ -279,24 +277,4 @@ } -// ha/dec -> alt/az -int altaz (double *alt, double *az, double ha, double dec, double latitude) { - - double sind, sinh, cosh; - - sind = dSIN (dec) * dSIN (latitude) + dCOS (dec) * dCOS (ha) * dCOS (latitude); - *alt = DEG_RAD * asin (sind); - - sinh = - dCOS (dec) * dSIN (ha); - cosh = dSIN (dec) * dCOS (latitude) - dCOS (dec) * dCOS (ha) * dSIN (latitude); - - *az = DEG_RAD * atan2 (sinh, cosh); - - // I may need to use the parallactic angle to get the right plate rotation... - // sinh = -dCOS(az) * dSIN(alt) * dSIN(ha) * dSIN(latitude) + dSIN(az) * dSIN(alt) * dCOS(ha) - dSIN(ha) * dCOS(alt) * dCOS(latitude); - // cosh = -dSIN(az) * dSIN(ha) * dSIN(latitude) - dCOS(az) * dCOS(ha); - // rot = -DEG_RAD * atan2 (sinh, cosh); - return TRUE; -} - int SetSDSSFlags (Stars *star, unsigned int flags1, unsigned int flags2) { Index: /trunk/Ohana/src/addstar/src/args_loadsupercos.c =================================================================== --- /trunk/Ohana/src/addstar/src/args_loadsupercos.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/args_loadsupercos.c (revision 33653) @@ -0,0 +1,90 @@ +# include "addstar.h" +static void help (void); + +AddstarClientOptions args_loadsupercos (int *argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + help (); + } + + // a global used by find_matches_refstars.c (value is 1 except for load2mass & loadwise) + NREFSTAR_GROUP = 1; + + /*** check for command line options ***/ + + /* basic mode: image, list, refcat */ + options.mode = M_REFCAT; + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + /* 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; + if ((N = get_argument (*argc, argv, "-only-match"))) { + options.only_match = TRUE; + remove_argument (N, argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (*argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + /* other addstar options which cannot be used in load2mass */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (*argc < 4) { + fprintf (stderr, "USAGE: loadsupercos [options] (surveys.csv) (plates.csv) (detections.bin) [..more files]\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadsupercos [options] (surveys.csv) (plates.csv) (detections.bin) [..more files]\n"); + fprintf (stderr, " add data from Supercosmos to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -only-match : only add measurements to existing objects\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} Index: /trunk/Ohana/src/addstar/src/args_loadwise.c =================================================================== --- /trunk/Ohana/src/addstar/src/args_loadwise.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/args_loadwise.c (revision 33653) @@ -0,0 +1,90 @@ +# include "addstar.h" +static void help (void); + +AddstarClientOptions args_loadwise (int argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + // a global used by find_matches_refstars.c (value is 1 except for load2mass & loadwise) + NREFSTAR_GROUP = 4; + + /*** check for command line options ***/ + + /* basic mode: image, list, refcat */ + options.mode = M_REFCAT; + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + /* 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; + if ((N = get_argument (argc, argv, "-only-match"))) { + options.only_match = TRUE; + remove_argument (N, &argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, &argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + /* other addstar options which cannot be used in load2mass */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (argc < 2) { + fprintf (stderr, "USAGE: load2mass [options] (wisefile) [..more files]\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadwise [options] (wisefile) [..more files]\n"); + fprintf (stderr, " add data from WISE catalog to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -only-match : only add measurements to existing objects\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} Index: /trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/find_matches.c (revision 33653) @@ -291,11 +291,16 @@ for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].Map = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000; + catalog[0].secfilt[Nave*Nsecfilt+j].flags = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].dummy[0] = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].dummy[1] = 0; } Index: /trunk/Ohana/src/addstar/src/find_matches_closest.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 33653) @@ -292,11 +292,16 @@ for (j = 0; j < Nsecfilt; j++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].M = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].Map = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN; + catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].M_20 = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].M_80 = NAN_S_SHORT; + catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].ubercalDist = 1000; + catalog[0].secfilt[Nave*Nsecfilt+j].flags = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].dummy[0] = 0; + catalog[0].secfilt[Nave*Nsecfilt+j].dummy[1] = 0; } Index: /trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 33653) @@ -13,5 +13,10 @@ int Nsecfilt; - if ((NREFSTAR_GROUP != 1) && (NREFSTAR_GROUP != 3)) { + switch (NREFSTAR_GROUP) { + case 1: + case 3: + case 4: + break; + default: fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n"); exit (1); Index: /trunk/Ohana/src/addstar/src/loadsupercos.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadsupercos.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/loadsupercos.c (revision 33653) @@ -0,0 +1,73 @@ +# include "addstar.h" +# include "supercos.h" + +/* This is the DVO program to upload Supercosmos detections into a DVO database. It is modeled + on the loadwise program. Like that case, it is expected to be run only rarely (once?). The + supercosmos data are delivered as *.bin files. The format of the file is described by the + SuperCosDetection structure. It does not allow a subset of the sky to be uploaded; entire + Supercosmos files are loaded if supplied on the command line. + + USAGE: loadsupercos -D CATDIR (catdir) (surveys.csv) (plates.csv) (scfile) [...more files] + */ + +int main (int argc, char **argv) { + + int i, status; + SkyTable *sky; + SkyList *skylist = NULL; + AddstarClientOptions options; + FITS_DB db; + + // need to construct these options with args_loadWISE... + options = ConfigInit (&argc, argv); + options = args_loadsupercos (&argc, argv, options); + + // load the full sky description table: + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, 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; + } + + // load the images table (needed regardless) + int Nsurvey = 0; + Survey *survey = loadsupercos_survey(argv[1], &Nsurvey); + + // load the images table (needed regardless) + int Nimage = 0; + Image *image = loadsupercos_plates(survey, Nsurvey, argv[2], &Nimage); + + // attempt to open the existing images table. if it does not exist, we will save the loaded images + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + status = dvo_image_lock (&db, ImageCat, 3600.0, LCK_XCLD); // 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) { + if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); + dvo_image_create (&db, GetZeroPoint()); + dvo_image_addrows (&db, image, Nimage); + SetProtect (TRUE); + dvo_image_update (&db, VERBOSE); + SetProtect (FALSE); + if (VERBOSE) fprintf (stderr, "created image table for Supercosmos\n"); + } + + int *imlist = loadsupercos_image_index (image, Nimage); + + for (i = 3; i < argc; i++) { + fprintf (stderr, "loading %s\n", argv[i]); + loadsupercos_rawdata (image, imlist, Nimage, skylist, argv[i], options); + } + exit (0); +} + Index: /trunk/Ohana/src/addstar/src/loadsupercos_ops.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadsupercos_ops.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/loadsupercos_ops.c (revision 33653) @@ -0,0 +1,122 @@ +# include "addstar.h" +# include "supercos.h" + +// this is defined in libohana/src/string.c but not generally exposed +char *_parse_nextword_csv (char *string); + +int loadsupercos_getFilterInfo (char *line, char *emulsion, char *filterID) { + + int i; + int Nchar = 0; + char *p1 = line; + + // emulsion is field 11, filterID is field 12 + for (i = 1; i < 11; i++) { + p1 = _parse_nextword_csv (p1); + if (!p1) { + fprintf (stderr, "error parsing filter info for line %s\n", line); + abort(); + } + } + if (*p1 == ',') { + fprintf (stderr, "missing emulsion info for line %s\n", line); + abort(); + } + + char *p2 = _parse_nextword_csv (p1); + if (!p2) { + fprintf (stderr, "error parsing filter info for line %s\n", line); + abort(); + } + if (*p2 == ',') { + fprintf (stderr, "missing filter ID for line %s\n", line); + abort(); + } + + char *p3 = _parse_nextword_csv (p2); + if (!p3) { + fprintf (stderr, "error parsing filter info for line %s\n", line); + abort(); + } + + Nchar = p2 - p1 - 1; + strncpy (emulsion, p1, Nchar); + emulsion[Nchar] = 0; + stripwhite (emulsion); + + Nchar = p3 - p2 - 1; + strncpy (filterID, p2, Nchar); + filterID[Nchar] = 0; + stripwhite (filterID); + + return TRUE; +} + +int loadsupercos_getST (char *line, double *st) { + + int i; + char *p1 = line; + + // lstObs is field 20 + for (i = 1; i < 20; i++) { + p1 = _parse_nextword_csv (p1); + if (!p1) { + fprintf (stderr, "error parsing filter info for line %s\n", line); + abort(); + } + } + if (*p1 == ',') { + fprintf (stderr, "missing sidereal time info for line %s\n", line); + abort(); + } + + char hour[3], minute[3]; + strncpy(hour, p1, 2); hour[2] = 0; + strncpy(minute, p1+2, 2); minute[2] = 0; + + double minuteF = atof (minute); + double hourF = atof(hour); + *st = hourF + minuteF / 60.0; + + return TRUE; +} + +int loadsupercos_sortStars (Stars *tstars, int Ntstars) { + +# define SWAPFUNC(A,B){ Stars temp = tstars[A]; tstars[A] = tstars[B]; tstars[B] = temp; } +# define COMPARE(A,B)(tstars[A].average.R < tstars[B].average.R) + + OHANA_SORT (Ntstars, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + + return TRUE; +} + +int *loadsupercos_image_index (Image *image, int Nimage) { + + int i, *index; + + // note that externID is unsigned, so we cannot set max to -1 + int maxIndex = 0; + for (i = 0; i < Nimage; i++) { + maxIndex = MAX (maxIndex, image[i].externID); + } + + ALLOCATE (index, int, maxIndex + 1); + for (i = 0; i < maxIndex + 1; i++) { + index[i] = -1; + } + + for (i = 0; i < Nimage; i++) { + int plateID = image[i].externID; + if (plateID < 0) abort(); + if (plateID > maxIndex) abort(); + if (index[plateID] != -1) abort(); + index[plateID] = i; + } + + return index; + +} Index: /trunk/Ohana/src/addstar/src/loadsupercos_plates.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadsupercos_plates.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/loadsupercos_plates.c (revision 33653) @@ -0,0 +1,347 @@ +# include "addstar.h" +# include "supercos.h" + +Image *loadsupercos_plates (Survey *survey, int Nsurvey, char *filename, int *nimage) { + + // these are fields read directly from the file. some are supplied unchanged to images, + // some require adjustments. + + double mjd, exptime, sidtime, RAo, DECo, aXmin, aXmax, aYmin, aYmax, alt, az, stepsize; + int nstar, survey_id, plate_id; + char line[3000], emulsion[128], filterID[128]; + + FILE *f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "failed to open SuperCOSMOS plate table file %s\n", filename); + exit (1); + } + if (scan_line (f, line) == EOF) { + fprintf (stderr, "error reading header line of %s\n", filename); + exit (1); + } + if (scan_line (f, line) == EOF) { + fprintf (stderr, "error reading empty plate line of %s\n", filename); + exit (1); + } + + int NIMAGE = 1000; + int Nimage = 0; + Image *image = NULL; + + ALLOCATE (image, Image, NIMAGE); + + int status = TRUE; + while (scan_line (f, line) != EOF) { + + iparse_csv (&plate_id, 1, line); + iparse_csv (&survey_id, 2, line); + if (survey_id >= Nsurvey) { + fprintf (stderr, "unknown survey %d\n", survey_id); + abort(); + } + + dparse_csv (&mjd, 28, line); // NOTE dparse is 1 counting (but 0 == 1?) + iparse_csv (&nstar, 106, line); + dparse_csv (&aXmin, 39, line); + dparse_csv (&aXmax, 40, line); + dparse_csv (&aYmin, 41, line); + dparse_csv (&aYmax, 42, line); + + dparse_csv (&stepsize, 37, line); // microns / pixel + dparse_csv (&exptime, 21, line); + dparse_csv (&sidtime, 20, line); + loadsupercos_getST (line, &sidtime); + + dparse_csv (&RAo, 14, line); + dparse_csv (&DECo, 15, line); + + loadsupercos_getFilterInfo (line, emulsion, filterID); + + // need to pass sidereal time in degrees + altaz (&alt, &az, 15.0*sidtime - RAo, DECo, survey[survey_id].latitude); + + float scale = stepsize * survey[survey_id].plateScale / 1000.0 / 3600.0; // scale in degrees / pixel + + image[Nimage].tzero = ohana_mjd_to_sec (mjd); + image[Nimage].nstar = nstar; + image[Nimage].secz = 1.0 / cos(RAD_DEG*alt); + + image[Nimage].NX = 5.00 / scale; + image[Nimage].NY = 4.75 / scale; + + image[Nimage].apmifit = 0.0; + image[Nimage].dapmifit = 0.0; + image[Nimage].Mcal = 0.0; + image[Nimage].dMcal = 0.0; + image[Nimage].Xm = 0.0; + + char photname[128]; + snprintf (photname, 128, "SCOS.%s.%s", emulsion, filterID); + image[Nimage].photcode = GetPhotcodeCodebyName (photname); + if (!image[Nimage].photcode) { + fprintf (stderr, "unknown photcode %s\n", photname); + fprintf (stderr, "line %d: %s\n", Nimage, line); + abort(); + } + image[Nimage].exptime = exptime * 60.0; + image[Nimage].sidtime = sidtime; // sidereal time in hours + image[Nimage].latitude = survey[survey_id].latitude; + + image[Nimage].RAo = RAo; + image[Nimage].DECo = DECo; + image[Nimage].Radius = NAN; + + snprintf (image[Nimage].name, 121, "SCOS.%02d.%06d", survey_id, plate_id); + image[Nimage].imageID = Nimage + 1; + image[Nimage].externID = plate_id; + + image[Nimage].detection_limit = 0; + image[Nimage].saturation_limit = 0; + image[Nimage].cerror = 0; + image[Nimage].fwhm_x = 0; + image[Nimage].fwhm_y = 0; + image[Nimage].trate = 0.0; + image[Nimage].ccdnum = 0; + image[Nimage].flags = 0; + image[Nimage].parentID = 0; + image[Nimage].sourceID = 0; + image[Nimage].nLinkAstrom = 0; + image[Nimage].nLinkPhotom = 0; + image[Nimage].dXpixSys = NAN; + image[Nimage].dYpixSys = NAN; + image[Nimage].dMagSys = NAN; + image[Nimage].nFitAstrom = 0; + image[Nimage].nFitPhotom = 0; + image[Nimage].photom_map_id = 0; + image[Nimage].astrom_map_id = 0; + + + // for now, we define a fake coordinate system based on the boresite center + strcpy (image[Nimage].coords.ctype, "RA---TAN"); + + image[Nimage].coords.crval1 = RAo; + image[Nimage].coords.crval2 = DECo; + + image[Nimage].coords.crpix1 = 0.5*image[Nimage].NX; + image[Nimage].coords.crpix2 = 0.5*image[Nimage].NY; + image[Nimage].coords.cdelt1 = image[Nimage].coords.cdelt2 = scale; + + image[Nimage].coords.pc1_1 = 0.0; + image[Nimage].coords.pc1_2 = 1.0; + image[Nimage].coords.pc2_1 = 1.0; + image[Nimage].coords.pc2_2 = 0.0; + + image[Nimage].coords.Npolyterms = 0; + memset (image[Nimage].coords.polyterms, 0, 2*7*sizeof(float)); + + Nimage ++; + CHECK_REALLOCATE (image, Image, NIMAGE, Nimage, 1000); + } + fclose (f); + + if (!status) abort(); + + *nimage = Nimage; + return image; +} + +/*** Image fields and mappings from supercosmos ***/ + +// tzero : mjd [27] +// nstar : objNum [105] +// secz : need to get telescope lattitude & longitude numbers +// NX : from aXmin, aXmax? [38,39] +// NY : from aXmin, aXmax? [40,41] +// apmifit : 0.0 +// dapmifit : 0.0 +// Mcal : 0.0 +// dMcal : 0.0 +// Xm : NAN +// photcode : from filterID [11] +// exptime : expLength (min) [20] +// sidtime : lstObs [19] +// latitude : (from survey_id) [01] +// RAo : raPnt [13] +// DECo : decPnt [14] +// Radius : (calculate) +// DUMMY : +// name : generate from plateID? [00] +// detection_limit : NAN +// saturation_limit : NAN +// cerror : NAN +// fwhm_x : NAN +// fwhm_y : +// trate : +// ccdnum : +// flags : +// imageID : +// parentID : +// externID : +// sourceID : +// nLinkAstrom : +// nLinkPhotom : +// dummy3 : +// dXpixSys : +// dYpixSys : +// dMagSys : +// nFitAstrom : +// nFitPhotom : +// photom_map_id : +// astrom_map_id : + +/*** SuperCOSMOS plate table (downloaded from the SSA). CSV file contains the following + * fields: + ***/ + +// 1 plateID +// 2 surveyID +// 3 fieldID +// 4 plateNum +// 5 dateMeas +// 6 timeMeas +// 7 instrument +// 8 softVersion +// 9 operator +// 10 scanMaterial +// 11 emulsion +// 12 filterID +// 13 filterComm +// 14 raPnt +// 15 decPnt +// 16 radecSys +// 17 equinox +// 18 eqTsys +// 19 utDateObs +// 20 lstObs +// 21 expLength +// 22 epoch +// 23 domeTemp +// 24 domePressure +// 25 domeHumid +// 26 waveEffect +// 27 tropl +// 28 mjd +// 29 focusNX +// 30 focusNY +// 31 calType +// 32 stepWdg +// 33 nSteps +// 34 orientation +// 35 emulpos +// 36 sosp +// 37 stepsize +// 38 scanlen +// 39 aXmin +// 40 aXmax +// 41 aYmin +// 42 aYmax +// 43 xPnt +// 44 yPnt +// 45 areaCut +// 46 apParam +// 47 dbParam +// 48 dbAmin +// 49 dbAmax +// 50 dbAcut +// 51 dbLevel +// 52 skySquare +// 53 skyDefn +// 54 skyFilter +// 55 fThresh +// 56 fSclen +// 57 pCut +// 58 starCat +// 59 brightLim +// 60 faintLim +// 61 equinRef +// 62 tSysRef +// 63 maxIter +// 64 rCritIni +// 65 rCritAbs +// 66 rCritRel +// 67 rCritFin +// 68 distType +// 69 scdGrid +// 70 raCol +// 71 decCol +// 72 raPMCol +// 73 decPMCol +// 74 plxCol +// 75 rvCol +// 76 magCol +// 77 cEquinox +// 78 cEqTSys +// 79 cEpoch +// 80 cEpTSys +// 81 starsInit +// 82 starsSelec +// 83 starsUsed +// 84 coeffs1 +// 85 coeffs2 +// 86 coeffs3 +// 87 coeffs4 +// 88 coeffs5 +// 89 coeffs6 +// 90 cubicDist +// 91 obsRaPnt +// 92 obsDecPnt +// 93 shiftPnt +// 94 areaMin +// 95 areaMax +// 96 magMin +// 97 magMax +// 98 ellipMin +// 99 ellipMax +//100 ellipMode +//101 orientMode +//102 ellipMed +//103 orientMed +//104 ellipMean +//105 orientMean +//106 objNum +//107 nParents +//108 lane1Count +//109 lane2Count +//110 lane3Count +//111 lane4Count +//112 lane5Count +//113 lane6Count +//114 lane7Count +//115 lane8Count +//116 lane9Count +//117 lane10Count +//118 lane11Count +//119 lane12Count +//120 lane13Count +//121 lane14Count +//122 lane15Count +//123 lane16Count +//124 lane17Count +//125 lane18Count +//126 lane19Count +//127 lane20Count +//128 lane21Count +//129 lane22Count +//130 lane23Count +//131 lane24Count +//132 lane25Count +//133 lane26Count +//134 lane27Count +//135 lane28Count +//136 lane29Count +//137 lane30Count +//138 astResidX +//139 astResidY +//140 nCalCoeffs +//141 CalXK +//142 gradient +//143 getCal1 +//144 getCal2 +//145 getCal3 +//146 getCal4 +//147 getCal5 +//148 getCal6 +//149 ffZp +//150 galZp +//151 galGrad +//152 directory Index: /trunk/Ohana/src/addstar/src/loadsupercos_rawdata.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadsupercos_rawdata.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/loadsupercos_rawdata.c (revision 33653) @@ -0,0 +1,222 @@ +# include "addstar.h" +# include "supercos.h" + +# define NBYTE 228 +# define NRECORDS 1000000 +# define DEBUG 0 + +int loadsupercos_rawdata (Image *image, int *imlist, int Nimage, SkyList *skytable, char *filename, AddstarClientOptions options) { + + int i, j, verbose; + int Nstars, NSTARS, Nrecords; + + FILE *f; + char *buffer; + + Detection *sstars; + Stars **stars, *tstars; + SkyList *skylist; + SkyRegion *region; + Catalog catalog; + + // the data files are written as an integer number of 228 byte records + ALLOCATE (buffer, char, NRECORDS*NBYTE); + + f = fopen (filename, "r"); + if (f == NULL) Shutdown ("can't read the Supercosmos data file: %s", filename); + + ALLOCATE (tstars, Stars, NRECORDS); + ALLOCATE (sstars, Detection, NRECORDS); + + double Rmin = 360.0; + double Rmax = 0.0; + double Dmin = +90.0; + double Dmax = -90.0; + + while ((Nrecords = fread (buffer, NBYTE, NRECORDS, f)) != 0) { + if (Nrecords == -1) Shutdown ("error reading from raw file %s", filename); + if (DEBUG) fprintf (stderr, "read %d bytes (%d bytes x %d rows)", NBYTE * Nrecords, NBYTE, Nrecords); + + if (VERBOSE) fprintf (stderr, "read .. "); + + // extract the basic values from the buffer + for (i = 0; i < Nrecords; i++) { + sstars[i].objID = *(int64_t *) &buffer[i*NBYTE + 0]; + sstars[i].surveyID = *(int8_t *) &buffer[i*NBYTE + 8]; + sstars[i].plateID = *(int32_t *) &buffer[i*NBYTE + 9]; + sstars[i].ra = *(double *) &buffer[i*NBYTE + 33]; + sstars[i].dec = *(double *) &buffer[i*NBYTE + 41]; + sstars[i].xCen = *(double *) &buffer[i*NBYTE + 129]; + sstars[i].yCen = *(double *) &buffer[i*NBYTE + 137]; + sstars[i].aU = *(float *) &buffer[i*NBYTE + 145]; + sstars[i].bU = *(float *) &buffer[i*NBYTE + 149]; + sstars[i].thetaU = *(int16_t *) &buffer[i*NBYTE + 154]; + sstars[i].class = *(int8_t *) &buffer[i*NBYTE + 165]; + sstars[i].quality = *(int32_t *) &buffer[i*NBYTE + 204]; + sstars[i].prfStat = *(float *) &buffer[i*NBYTE + 208]; + sstars[i].prfMag = *(float *) &buffer[i*NBYTE + 212]; + sstars[i].gMag = *(float *) &buffer[i*NBYTE + 216]; + sstars[i].sMag = *(float *) &buffer[i*NBYTE + 220]; + } + + // convert to the Stars format + for (i = 0; i < Nrecords; i++) { + InitStar (&tstars[i]); + tstars[i].measure.Xccd = sstars[i].xCen; + tstars[i].measure.Yccd = sstars[i].yCen; + + tstars[i].measure.M = sstars[i].sMag; + tstars[i].measure.Map = sstars[i].gMag; + tstars[i].measure.photFlags = sstars[i].class; + + tstars[i].average.R = sstars[i].ra; + tstars[i].average.D = sstars[i].dec; + + int Ni = imlist[sstars[i].plateID]; + if (Ni == -1) abort(); + + // XXX fix these + tstars[i].measure.Mcal = 0.0; + tstars[i].measure.dt = image[Ni].exptime; + + double sidtime = image[Ni].sidtime; + double latitude = image[Ni].latitude; + + double alt, az; + altaz (&alt, &az, 15.0*sidtime - tstars[i].average.R, tstars[i].average.D, latitude); + + tstars[i].measure.airmass = 1.0 / cos(RAD_DEG*alt); + tstars[i].measure.az = az; + tstars[i].measure.t = image[Ni].tzero; + tstars[i].measure.imageID = image[Ni].imageID; + tstars[i].measure.photcode = image[Ni].photcode; + + tstars[i].measure.FWx = ToShortPixels(sstars[i].aU); + tstars[i].measure.FWy = ToShortPixels(sstars[i].bU); + tstars[i].measure.theta = ToShortDegrees(sstars[i].thetaU); + + Rmin = MIN (Rmin, tstars[i].average.R); + Rmax = MAX (Rmax, tstars[i].average.R); + Dmin = MIN (Dmin, tstars[i].average.D); + Dmax = MAX (Dmax, tstars[i].average.D); + } + + if (VERBOSE) fprintf (stderr, "scan %d stars (%10.6f - %10.6f, %10.6f - %10.6f) .. ", Nrecords, Rmin, Rmax, Dmin, Dmax); + + // sort the tstars by RA + loadsupercos_sortStars (tstars, Nrecords); + + // scan through the stars, loading the containing catalogs + // skip through table for unsaved stars + for (i = 0; i < Nrecords; i++) { + if (tstars[i].found != -1) continue; + + // scan forward until we reach the UserPatch (stars sorted by RA) + if (tstars[i].average.R < UserPatch.Rmin) continue; + if (tstars[i].average.R > UserPatch.Rmax) break; + if (tstars[i].average.D < UserPatch.Dmin) continue; + if (tstars[i].average.D > UserPatch.Dmax) continue; + + // identify the relevant catalog + skylist = SkyRegionByPoint_List (skytable, -1, tstars[i].average.R, tstars[i].average.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]); + + int this_image_id = tstars[i].measure.imageID; + + // loop over stars in this WISE region that are also in this output region + Nstars = 0; + NSTARS = 1000; + ALLOCATE (stars, Stars *, NSTARS); + + for (j = i; j < Nrecords; j++) { + if (tstars[j].found != -1) continue; + + if (tstars[j].measure.imageID != this_image_id) continue; + + // check if in skyregion + if (tstars[j].average.R < region[0].Rmin) continue; + if (tstars[j].average.R > region[0].Rmax) break; + if (tstars[j].average.D < region[0].Dmin) continue; + if (tstars[j].average.D > region[0].Dmax) continue; + + // check if in UserPatch + if (tstars[j].average.R < UserPatch.Rmin) continue; + if (tstars[j].average.R > UserPatch.Rmax) break; + if (tstars[j].average.D < UserPatch.Dmin) continue; + if (tstars[j].average.D > UserPatch.Dmax) continue; + + stars[Nstars] = &tstars[j]; + Nstars ++; + + CHECK_REALLOCATE (stars, Stars *, NSTARS, Nstars, 1000); + } + + 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); + + if (1) { + verbose = VERBOSE; + VERBOSE = FALSE; + + // now we have all of the loaded stars in this catalog + catalog.filename = skylist[0].filename[0]; + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[0], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + // assume no input star matches an existing star + // simply add to the existing table + // loadWISE_catalog (&catalog, stars, Nstars); + find_matches_refstars (skylist[0].regions[0], stars, Nstars, &catalog, options); + // loadWISE_catalog (&catalog, stars, Nstars); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + // free (catalog.filename); + // XXX don't free this! it points to an element of the skytable + } + + SkyListFree (skylist); + + free (stars); + VERBOSE = verbose; + } + if (VERBOSE) fprintf (stderr, "done\n"); + } + + if (VERBOSE) fprintf (stderr, "\n"); + + free (tstars); + free (sstars); + + fclose (f); + free (buffer); + return (TRUE); +} + +/* + for each WISE file: + for each data block + generate a table of: R, D, byte, flag + for each unsaved star + find containing catalog + load catalog + find all contained stars + add to catalog + save catalog + mark all contained stars +*/ + Index: /trunk/Ohana/src/addstar/src/loadsupercos_survey.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadsupercos_survey.c (revision 33653) +++ /trunk/Ohana/src/addstar/src/loadsupercos_survey.c (revision 33653) @@ -0,0 +1,69 @@ +# include "addstar.h" +# include "supercos.h" + +Survey *loadsupercos_survey (char *filename, int *nsurvey) { + + char line[1024]; + + FILE *f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "failed to open SuperCOSMOS survey table file %s\n", filename); + exit (1); + } + if (scan_line (f, line) == EOF) { + fprintf (stderr, "error reading header line of %s\n", filename); + exit (1); + } + + int NSURVEY = 20; + int Nsurvey = 0; + Survey *survey = NULL; + ALLOCATE (survey, Survey, NSURVEY); + while (scan_line (f, line) != EOF) { + + iparse_csv (&survey[Nsurvey].survey_id, 1, line); + dparse_csv (&survey[Nsurvey].longitude, 10, line); + dparse_csv (&survey[Nsurvey].latitude, 11, line); + dparse_csv (&survey[Nsurvey].plateScale, 13, line); // arcsec / mm + + if (Nsurvey != survey[Nsurvey].survey_id) { + fprintf (stderr, "NEED to SORT surveys\n"); + abort(); + } + + Nsurvey ++; + + CHECK_REALLOCATE (survey, Survey, NSURVEY, Nsurvey, 10); + } + fclose (f); + + *nsurvey = Nsurvey; + return survey; +} + +// fields in the survey.csv file (this could be read from the header line, but it is too much work) +// 1 surveyID +// 2 surveyName +// 3 systemID +// 4 fieldOfView +// 5 decMin +// 6 decMax +// 7 numFields +// 8 telescope +// 9 telAperture +// 10 telLong +// 11 telLat +// 12 telHeight +// 13 plateScale +// 14 colour +// 15 waveMin +// 16 waveMax +// 17 waveEff +// 18 magLimit +// 19 epochMin +// 20 epochMax +// 21 epTsys +// 22 equinox +// 23 eqTsys +// 24 surveyRef + Index: /trunk/Ohana/src/addstar/src/loadwise.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadwise.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/loadwise.c (revision 33653) @@ -2,16 +2,24 @@ # include "WISE.h" +/* This is the DVO program to upload WISE detections into a DVO database. + It is modeled on the load2mass program, but unlike that case, it is expected to be run only + rarely (once?). The WISE data are delivered as *.bz2 files. This program assumes the files + have been decompressed. It does not allow a subset of the sky to be uploaded; entire WISE + files are loaded if supplied on the command line. + + USAGE: loadwise -D CATDIR (catdir) (wisefile) [...more files] + + */ + int main (int argc, char **argv) { - char *path; int i; - SkyTable *sky, *skyWISE; + SkyTable *sky; SkyList *skylist = NULL; - SkyList *overlap = NULL; AddstarClientOptions options; // need to construct these options with args_loadWISE... options = ConfigInit (&argc, argv); - options = args_loadWISE (argc, argv, options); + options = args_loadwise (argc, argv, options); // load the full sky description table: @@ -30,23 +38,9 @@ } - path = WISE_DIR_PRELIM; - - // the accel.dat file has the raw filenames - // test if the file exists, or else try the .gz version - // XXX can we build an accelerator? - skyWISE = loadWISE_acc (path, "accel.dat"); - - for (i = 0; i < skyWISE[0].Nregions; i++) { - // check if any of the skylist entries overlap this WISE catalog: - overlap = SkyListByBounds_List (skylist, -1, skyWISE[0].regions[i].Rmin, skyWISE[0].regions[i].Rmax, skyWISE[0].regions[i].Dmin, skyWISE[0].regions[i].Dmax); - if (overlap[0].Nregions == 0) { - SkyListFree (overlap); - continue; - } - - fprintf (stderr, "loading %s\n", skyWISE[0].filename[i]); - loadWISE_as_rawdata (overlap, skyWISE[0].filename[i], options); - SkyListFree (overlap); + for (i = 1; i < argc; i++) { + fprintf (stderr, "loading %s\n", argv[i]); + loadwise_prelim_rawdata (skylist, argv[i], options); } exit (0); } + Index: /trunk/Ohana/src/addstar/src/loadwise_ops.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadwise_ops.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/loadwise_ops.c (revision 33653) @@ -1,9 +1,4 @@ # include "addstar.h" # include "WISE.h" - -static int FilterSkip; -static int TimeSkip; -static int Qentry; -static int Photcode; int getWISE_setup () { Index: /trunk/Ohana/src/addstar/src/loadwise_prelim_rawdata.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadwise_prelim_rawdata.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/loadwise_prelim_rawdata.c (revision 33653) @@ -10,5 +10,5 @@ # define DEBUG 0 -int loadWISE_as_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) { +int loadwise_prelim_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) { int i, j, verbose; @@ -159,5 +159,5 @@ stars[Nstars+3][0].average.R = tstars[j].R; stars[Nstars+3][0].average.D = tstars[j].D; - getWISE_star_full (&stars[Nstars], &buffer[offset], Nbyte - offset); + loadwise_star_full (&stars[Nstars], &buffer[offset], Nbyte - offset); tstars[j].flag = TRUE; Index: /trunk/Ohana/src/addstar/src/loadwise_star_full.c =================================================================== --- /trunk/Ohana/src/addstar/src/loadwise_star_full.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/loadwise_star_full.c (revision 33653) @@ -2,88 +2,96 @@ # include "WISE.h" -// XXX check to see if desired output format is PS1_V1 or later? (use 16bit version if not?) +int setWISE_ph_qual (Stars *star, char qual); +int setWISE_var_flag (Stars *star, char qual); +int setWISE_ext_flag (Stars **star, char *ptr); +int setWISE_cc_flag (Stars *star, char qual); +int setWISE_sat_flag (Stars *star, char *ptr); +int setWISE_blend_flag (Stars **star, char *ptr); // fill in the data for a WISE quad star. takes a pointer to the start of the line the // RA and DEC have already been set -int getWISE_star_full (Stars **star, char *line, int Nmax) { - +int loadwise_star_full (Stars **star, char *line, int Nmax) { + + int i; char *ptr; - double dMfull; - double jd; if (line == NULL) Shutdown ("format error in WISE"); + ptr = line; + + // I can assign dRA and dDEC to dX and dY if I can consistently set posangle and pltscale + star[0][0].measure.posangle = 0.0; + star[0][0].measure.pltscale = 1.0; + ptr = skipNbounds (ptr, '|', 3, Nmax); // skip: desig, ra, dec, - star[0][0].measure.FWx = ToShortPixels(strtod (ptr, NULL)); - ptr = nextWISEfield (ptr); // sig_dec - star[0][0].measure.FWy = ToShortPixels(strtod (ptr, NULL)); - ptr = nextWISEfield (ptr); // sig_radec - star[0][0].measure.theta = ToShortPixels(strtod (ptr, NULL)); - - star[1][0].measure.FWx = star[0][0].measure.FWx; - star[1][0].measure.FWy = star[0][0].measure.FWy; - star[1][0].measure.theta = star[0][0].measure.theta; - star[2][0].measure.FWx = star[0][0].measure.FWx; - star[2][0].measure.FWy = star[0][0].measure.FWy; - star[2][0].measure.theta = star[0][0].measure.theta; - star[3][0].measure.FWx = star[0][0].measure.FWx; - star[3][0].measure.FWy = star[0][0].measure.FWy; - star[3][0].measure.theta = star[0][0].measure.theta; - - ptr = skipNbounds (ptr, '|', 10, Nmax); // skip: glon, glat, elon, elat, wx, wy, cntr, source_id, coadd_id, src + star[0][0].measure.dXccd = ToShortPixels(strtod (ptr, NULL)); // sig_ra + ptr = nextWISEfield (ptr); + star[0][0].measure.dYccd = ToShortPixels(strtod (ptr, NULL)); // sig_dec + ptr = nextWISEfield (ptr); + + // we only know a single set of values for all 4 bands + for (i = 1; i < 4; i++) { + star[i][0].measure.dXccd = star[0][0].measure.dXccd; + star[i][0].measure.dYccd = star[0][0].measure.dYccd; + star[i][0].measure.posangle = star[0][0].measure.posangle; + star[i][0].measure.pltscale = star[0][0].measure.pltscale; + } + ptr = skipNbounds (ptr, '|', 5, Nmax); // skip: sig_radec, glon, glat, elon, elat + + star[0][0].measure.Xccd = strtod (ptr, NULL); // wx + ptr = nextWISEfield (ptr); // skip wx + star[0][0].measure.Yccd = strtod (ptr, NULL); // wy + ptr = nextWISEfield (ptr); // skip wy + + for (i = 1; i < 4; i++) { + star[i][0].measure.Xccd = star[0][0].measure.Xccd; + star[i][0].measure.Yccd = star[0][0].measure.Yccd; + } + ptr = skipNbounds (ptr, '|', 4, Nmax); // skip: cntr, source_id, coadd_id, src // W1 - star[0][0].measure.M = strtod (ptr, NULL); - ptr = nextWISEfield (ptr); // j_cmsig - star[0][0].measure.dM = strtod (ptr, NULL); - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w1snr, w1rchi2 - - // W2 - star[1][0].measure.M = strtod (ptr, NULL); - ptr = nextWISEfield (ptr); // j_cmsig - star[1][0].measure.dM = strtod (ptr, NULL); - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w2snr, w2rchi2 - - // W3 - star[2][0].measure.M = strtod (ptr, NULL); - ptr = nextWISEfield (ptr); // j_cmsig - star[2][0].measure.dM = strtod (ptr, NULL); - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w3snr, w3rchi2 - - // W4 - star[3][0].measure.M = strtod (ptr, NULL); - ptr = nextWISEfield (ptr); // j_cmsig - star[3][0].measure.dM = strtod (ptr, NULL); - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w4snr, w4rchi2 - - star[0][0].measure.photFlags = 0; - star[1][0].measure.photFlags = 0; - star[2][0].measure.photFlags = 0; - star[3][0].measure.photFlags = 0; - ptr = skipNbounds (ptr, '|', 8, Nmax); // skip fields - - setWISE_cc_flag (star[0], ptr[0]); // cc_flg - setWISE_cc_flag (star[1], ptr[1]); // cc_flg - setWISE_cc_flag (star[2], ptr[2]); // cc_flg - setWISE_cc_flag (star[3], ptr[2]); // cc_flg - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: ext_flg, var_flg - - setWISE_ph_qual (star[0], ptr[0]); // ph_qual - setWISE_ph_qual (star[1], ptr[1]); // ph_qual - setWISE_ph_qual (star[2], ptr[2]); // ph_qual - setWISE_ph_qual (star[3], ptr[3]); // ph_qual - ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: ext_flg, var_flg - - // setWISE_rd_flag (star[0], ptr[0]); - // setWISE_bl_flag (star[0], ptr[0]); - // setWISE_gal_flag (star[0], ptr[0]); - // setWISE_mp_flag (star[0], ptr[0]); - - // ptr = nextWISEfield (ptr); // x_scan - // star[0][0].measure.Xccd = strtod (ptr, NULL); - // star[2][0].measure.Xccd = star[1][0].measure.Xccd = star[0][0].measure.Xccd; - // ptr = nextWISEfield (ptr); // jdate (julian date) - - jd = 2455263.0; // NOTE : WISE prelim release does not contain per-detection time info. + for (i = 0; i < 4; i++) { + star[0][0].measure.M = strtod (ptr, NULL); // w?mpro + ptr = nextWISEfield (ptr); + star[0][0].measure.dM = strtod (ptr, NULL); // w?sigmpro + ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: w?sigmpro, w?snr + star[0][0].measure.psfChisq = strtod (ptr, NULL); // w?rchi2 + ptr = nextWISEfield (ptr); // skip : w1rchi2 + + // init the photFlags field + star[i][0].measure.photFlags = 0; + } + + ptr = nextWISEfield (ptr); // skip: rchi2 + + // set blend flags for all 4 measures + setWISE_blend_flag (star, ptr); // nb & na both used here + ptr = skipNbounds (ptr, '|', 2, Nmax); // skip: nb, na + + for (i = 0; i < 4; i++) { + setWISE_sat_flag (star[i], ptr); // w1sat + ptr = nextWISEfield (ptr); + } + ptr = nextWISEfield (ptr); // skip satnum + + for (i = 0; i < 4; i++) { + setWISE_cc_flag (star[i], ptr[i]); // cc_flg + } + ptr = nextWISEfield (ptr); // skip cc_flags + + // set ext flags for all 4 measures + setWISE_ext_flag (star, ptr); // ext_flg + ptr = nextWISEfield (ptr); // skip ext_flags + + for (i = 0; i < 4; i++) { + setWISE_var_flag (star[i], ptr[i]); // var_flg + } + ptr = nextWISEfield (ptr); // skip var_flags + + for (i = 0; i < 4; i++) { + setWISE_ph_qual (star[i], ptr[i]); // ph_qual + } + + double jd = 2455263.0; // NOTE : WISE prelim release does not contain per-detection time info. // the release is based on data taken in the period 14 January 2010 to 29 April 2010 star[0][0].measure.t = ohana_jd_to_sec (jd); @@ -92,48 +100,19 @@ star[3][0].measure.t = star[0][0].measure.t; - // ptr = nextWISEfield (ptr); // j_psfchi - // star[0][0].measure.psfChisq = strtod (ptr, NULL); - // ptr = nextWISEfield (ptr); // h_psfchi - // star[1][0].measure.psfChisq = strtod (ptr, NULL); - // ptr = nextWISEfield (ptr); // k_psfchi - // star[2][0].measure.psfChisq = strtod (ptr, NULL); - // - // ptr = nextWISEfield (ptr); // j_m_stdap - // star[0][0].measure.Map = strtod (ptr, NULL); - // ptr = nextWISEfield (ptr); // j_msig_stdap (skip?) - // - // ptr = nextWISEfield (ptr); // h_m_stdap - // star[1][0].measure.Map = strtod (ptr, NULL); - // ptr = nextWISEfield (ptr); // h_msig_stdap (skip?) - // - // ptr = nextWISEfield (ptr); // k_m_stdap - // star[2][0].measure.Map = strtod (ptr, NULL); - // ptr = nextWISEfield (ptr); // k_msig_stdap (skip?) - // - // ptr = nextWISEfield (ptr); // dist_edge_ns (skip) - // ptr = nextWISEfield (ptr); // dist_edge_ew (skip) - // ptr = nextWISEfield (ptr); // dist_edge_flg (skip) - // - // ptr = nextWISEfield (ptr); // dup_src (one flag for all filters) - // setWISE_dup_flag (star[0], ptr[0]); - // setWISE_dup_flag (star[1], ptr[0]); - // setWISE_dup_flag (star[2], ptr[0]); - // - // ptr = nextWISEfield (ptr); // use_src (one flag for all filters) - // setWISE_use_flag (star[0], ptr[0]); - // setWISE_use_flag (star[1], ptr[0]); - // setWISE_use_flag (star[2], ptr[0]); - - star[0][0].measure.photcode = TM_J; + star[0][0].measure.photcode = WISE_W1; star[0][0].measure.detID = 0; star[0][0].measure.imageID = 0; - star[1][0].measure.photcode = TM_H; + star[1][0].measure.photcode = WISE_W2; star[1][0].measure.detID = 0; star[1][0].measure.imageID = 0; - star[2][0].measure.photcode = TM_K; + star[2][0].measure.photcode = WISE_W3; star[2][0].measure.detID = 0; star[2][0].measure.imageID = 0; + + star[3][0].measure.photcode = WISE_W4; + star[3][0].measure.detID = 0; + star[3][0].measure.imageID = 0; return TRUE; @@ -153,15 +132,71 @@ } +# define FLAG_PH_A 0x00000001 // quality flag 'A' +# define FLAG_PH_B 0x00000002 // quality flag 'B' +# define FLAG_PH_C 0x00000004 // quality flag 'C' +# define FLAG_PH_U 0x00000008 // quality flag 'U' +# define FLAG_PH_X 0x00000010 // quality flag 'X' + +# define FLAG_SATURATED_PIX 0x00000100 // sat > 0.0 + +# define FLAG_CC_PERSIST 0x00010000 // 'p' or 'P' +# define FLAG_CC_HALO 0x00020000 // 'h' or 'H' +# define FLAG_CC_GHOST 0x00020000 // 'o' or 'O' (note overloading) +# define FLAG_CC_SPIKE 0x00040000 // 'd' or 'D' +# define FLAG_CC_SPURIOUS 0x00080000 // Upper Case letters + +# define FLAG_BLEND_ACTIVE 0x00100000 // nb > 1, na == 0 +# define FLAG_BLEND_PASSIVE 0x00200000 // nb > 1, na == 1 + +# define FLAG_EXTENDED 0x01000000 // ext == 1 + +# define FLAG_VARIABLE_LEVEL1 0x10000000 // var_flg == 1 to 4 +# define FLAG_VARIABLE_LEVEL2 0x20000000 // var_flg == 5 to 7 +# define FLAG_VARIABLE_LEVEL3 0x40000000 // var_flg == 8 or 9 + +int setWISE_blend_flag (Stars **star, char *ptr) { + + int nb = atoi (ptr); + if (nb == 1) return TRUE; + + ptr = nextWISEfield (ptr); // skip to na + int na = atoi (ptr); + + if (na == 0) { + star[0][0].measure.photFlags |= FLAG_BLEND_ACTIVE; + star[1][0].measure.photFlags |= FLAG_BLEND_ACTIVE; + star[2][0].measure.photFlags |= FLAG_BLEND_ACTIVE; + star[3][0].measure.photFlags |= FLAG_BLEND_ACTIVE; + } else { + star[0][0].measure.photFlags |= FLAG_BLEND_PASSIVE; + star[1][0].measure.photFlags |= FLAG_BLEND_PASSIVE; + star[2][0].measure.photFlags |= FLAG_BLEND_PASSIVE; + star[3][0].measure.photFlags |= FLAG_BLEND_PASSIVE; + } + return TRUE; +} + +int setWISE_sat_flag (Stars *star, char *ptr) { + + float sat = strtod (ptr, NULL); + if (sat > 0.0) { + star[0].measure.photFlags |= FLAG_SATURATED_PIX; + } + return TRUE; +} + int setWISE_cc_flag (Stars *star, char qual) { switch (qual) { - case 'p': star[0].measure.photFlags |= 0x00010000; break; - case 'P': star[0].measure.photFlags |= 0x00010000; break; - case 'h': star[0].measure.photFlags |= 0x00020000; break; - case 'H': star[0].measure.photFlags |= 0x00020000; break; - case 'd': star[0].measure.photFlags |= 0x00040000; break; - case 'D': star[0].measure.photFlags |= 0x00040000; break; - case 'o': star[0].measure.photFlags |= 0x00020000; break; - case 'O': star[0].measure.photFlags |= 0x00020000; break; + case 'p': star[0].measure.photFlags |= FLAG_CC_PERSIST; break; + case 'h': star[0].measure.photFlags |= FLAG_CC_HALO; break; + case 'd': star[0].measure.photFlags |= FLAG_CC_SPIKE; break; + case 'o': star[0].measure.photFlags |= FLAG_CC_GHOST; break; + + case 'P': star[0].measure.photFlags |= FLAG_CC_PERSIST | FLAG_CC_SPURIOUS; break; + case 'H': star[0].measure.photFlags |= FLAG_CC_HALO | FLAG_CC_SPURIOUS; break; + case 'D': star[0].measure.photFlags |= FLAG_CC_SPIKE | FLAG_CC_SPURIOUS; break; + case 'O': star[0].measure.photFlags |= FLAG_CC_GHOST | FLAG_CC_SPURIOUS; break; + case '0': break; default: @@ -172,12 +207,42 @@ } +int setWISE_ext_flag (Stars **star, char *ptr) { + + int ext = atoi (ptr); + if (ext == 1) { + star[0][0].measure.photFlags |= FLAG_EXTENDED; + star[1][0].measure.photFlags |= FLAG_EXTENDED; + star[2][0].measure.photFlags |= FLAG_EXTENDED; + star[3][0].measure.photFlags |= FLAG_EXTENDED; + } + return TRUE; +} + +int setWISE_var_flag (Stars *star, char qual) { + + int value = qual - '0'; + assert (value >= 0); + assert (value <= 9); + if (value == 0) return TRUE; + if (value < 5) { + star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL1; + return TRUE; + } + if (value < 8) { + star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL2; + return TRUE; + } + star[0].measure.photFlags |= FLAG_VARIABLE_LEVEL3; + return TRUE; +} + int setWISE_ph_qual (Stars *star, char qual) { switch (qual) { - case 'A': star[0].measure.photFlags |= 0x00000001; break; // was: 0x0004 - case 'B': star[0].measure.photFlags |= 0x00000002; break; // was: 0x0005 - case 'C': star[0].measure.photFlags |= 0x00000004; break; // was: 0x0006 - case 'U': star[0].measure.photFlags |= 0x00000040; break; // was: 0x0001 - case 'X': star[0].measure.photFlags |= 0x00000080; break; // was: 0x0000 + case 'A': star[0].measure.photFlags |= FLAG_PH_A; break; + case 'B': star[0].measure.photFlags |= FLAG_PH_B; break; + case 'C': star[0].measure.photFlags |= FLAG_PH_C; break; + case 'U': star[0].measure.photFlags |= FLAG_PH_U; break; + case 'X': star[0].measure.photFlags |= FLAG_PH_X; break; default: fprintf (stderr, "error!\n"); @@ -187,78 +252,27 @@ } -int setWISE_rd_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags |= 0x00000100; break; // was: 0x0000 - case '1': star[0].measure.photFlags |= 0x00000200; break; // was: 0x0010 - case '2': star[0].measure.photFlags |= 0x00000400; break; // was: 0x0020 - case '3': star[0].measure.photFlags |= 0x00000800; break; // was: 0x0030 - case '4': star[0].measure.photFlags |= 0x00001000; break; // was: 0x0040 - case '6': star[0].measure.photFlags |= 0x00002000; break; // was: 0x0050 - case '9': star[0].measure.photFlags |= 0x00004000; break; // was: 0x0060 - default: - fprintf (stderr, "error!\n"); - exit (2); - } - return (TRUE); -} - -int setWISE_bl_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags &= ~0x00300000; break; // was: ~0x0008 - case '1': star[0].measure.photFlags |= 0x00100000; break; // was: ~0x0008 - default: star[0].measure.photFlags |= 0x00200000; break; // was: 0x0008 - } - return (TRUE); -} - -int setWISE_gal_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags &= ~0x00c00000; break; // was: ~0x0080 - case '1': star[0].measure.photFlags |= 0x00400000; break; // was: ~0x0080 - default: star[0].measure.photFlags |= 0x00800000; // was: 0x0080 - star[0].measure.extNsigma = 100.0; - break; - } - return (TRUE); -} - -int setWISE_mp_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags &= ~0x03000000; break; // was: ~0x0800 - case '1': star[0].measure.photFlags |= 0x01000000; break; // was: ~0x0800 - default: star[0].measure.photFlags |= 0x02000000; break; // was: 0x0800 - } - return (TRUE); -} - -int setWISE_dup_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags &= ~0x0c000000; break; // was: ~0x1000 - case '1': star[0].measure.photFlags |= 0x04000000; break; // was: ~0x1000 - default: star[0].measure.photFlags |= 0x08000000; break; // was: 0x1000 - } - return (TRUE); -} - -int setWISE_use_flag (Stars *star, char qual) { - - switch (qual) { - case '0': star[0].measure.photFlags &= ~0x10000000; break; // was: ~0x2000 - case '1': star[0].measure.photFlags |= 0x10000000; break; // was: 0x2000 - default: abort(); - } - return (TRUE); -} - -// unused photFlags: -// 0x0000.8000 -// 0x0004.0000 -// 0x0008.0000 -// 0x2000.0000 -// 0x4000.0000 -// 0x8000.0000 +/* unused photFlag bits: + +0x0000.0020 +0x0000.0040 +0x0000.0080 + +0x0000.0200 +0x0000.0400 +0x0000.0800 + +0x0000.1000 +0x0000.2000 +0x0000.4000 +0x0000.8000 + +0x0040.0000 +0x0080.0000 + +0x0200.0000 +0x0400.0000 +0x0800.0000 + +0x8000.0000 + +*/ Index: /trunk/Ohana/src/addstar/src/mkcmf.c =================================================================== --- /trunk/Ohana/src/addstar/src/mkcmf.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/mkcmf.c (revision 33653) @@ -14,9 +14,12 @@ void gauss_init (int Nbin); double rnd_gauss (double mean, double sigma); -void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, int Nstars); +void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars); void writeStars_PS1_V1 (FTable *ftable, double *X, double *Y, double *M, int Nstars); void writeStars_PS1_DEV_1 (FTable *ftable, double *X, double *Y, double *M, int Nstars); void writeStars_PS1_DEV_0 (FTable *ftable, double *X, double *Y, double *M, int Nstars); +int ADDNOISE = TRUE; +float BAD_PSFQF_FRAC = 0.0; + int main (int argc, char **argv) { @@ -25,4 +28,5 @@ int N, Nstars, NSTARS, found; + unsigned int *Flag; double *X, *Y, *M, Xmax, Ymax; @@ -34,4 +38,19 @@ Coords coords; + static char *photcode = "SIMTEST.r.Chip"; + if ((N = get_argument (argc, argv, "-photcode"))) { + remove_argument (N, &argc, argv); + photcode = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + static double mjd; + mjd = NAN; + if ((N = get_argument (argc, argv, "-mjd"))) { + remove_argument (N, &argc, argv); + mjd = strtod (argv[N], NULL); + remove_argument (N, &argc, argv); + } + static char *date = "2001-01-01"; if ((N = get_argument (argc, argv, "-date"))) { @@ -58,4 +77,39 @@ } + double CRPIX1 = 0.0; + double CRPIX2 = 0.0; + if ((N = get_argument (argc, argv, "-crpix"))) { + remove_argument (N, &argc, argv); + CRPIX1 = atof (argv[N]); + remove_argument (N, &argc, argv); + CRPIX2 = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + float exptime = 1.0; + if ((N = get_argument (argc, argv, "-exptime"))) { + remove_argument (N, &argc, argv); + exptime = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + // XXX note that the airmass and ra,dec,mjd can be inconsistent (for a given observatory location) + float airmass = 1.0; + if ((N = get_argument (argc, argv, "-airmass"))) { + remove_argument (N, &argc, argv); + airmass = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + int NX = 0; + int NY = 0; + if ((N = get_argument (argc, argv, "-size"))) { + remove_argument (N, &argc, argv); + NX = atof (argv[N]); + remove_argument (N, &argc, argv); + NY = atof (argv[N]); + remove_argument (N, &argc, argv); + } + // add support for all cmf types static char *type = "PS1_V2"; @@ -66,4 +120,32 @@ } + // add support for all cmf types + int FROM_COORDS = FALSE; + if ((N = get_argument (argc, argv, "-coords"))) { + remove_argument (N, &argc, argv); + FROM_COORDS = TRUE; + } + + // expect a field with the photom flags for each detection + int READ_FLAGS = FALSE; + if ((N = get_argument (argc, argv, "-flags"))) { + remove_argument (N, &argc, argv); + READ_FLAGS = TRUE; + } + + // add support for all cmf types + ADDNOISE = TRUE; + if ((N = get_argument (argc, argv, "-no-noise"))) { + remove_argument (N, &argc, argv); + ADDNOISE = FALSE; + } + + // random bad PSF_QF values + if ((N = get_argument (argc, argv, "-bad-psfqf-frac"))) { + remove_argument (N, &argc, argv); + BAD_PSFQF_FRAC = atof(argv[N]); + remove_argument (N, &argc, argv); + } + if (argc != 3) { fprintf (stderr, "USAGE mkcmf (input) (output) [-date date] [-time time] [-radec ra dec] [-cmftype type]\n"); @@ -71,4 +153,21 @@ } + /* bore site center guess */ + strcpy (coords.ctype, "DEC--TAN"); + coords.crval1 = RA; + coords.crval2 = DEC; + coords.crpix1 = CRPIX1; + coords.crpix2 = CRPIX2; + + coords.cdelt1 = 0.25/3600.0; + coords.cdelt2 = 0.25/3600.0; + + coords.pc1_1 = 1; + coords.pc2_2 = 1; + coords.pc1_2 = 0; + coords.pc2_1 = 0; + coords.Npolyterms = 1; + + // load stars and generate complete output fields f = fopen (argv[1], "r"); if (f == NULL) { @@ -83,7 +182,33 @@ ALLOCATE (Y, double, NSTARS); ALLOCATE (M, double, NSTARS); + ALLOCATE (Flag, unsigned int, NSTARS); Xmax = Ymax = 0; - while (fscanf (f, "%lf %lf %lf", &X[Nstars], &Y[Nstars], &M[Nstars]) != EOF) { + while (1) { + int status; + double ra, dec, mag, xobs, yobs, xraw, yraw, mraw; + unsigned int flags; + if (FROM_COORDS) { + status = fscanf (f, "%lf %lf %lf %lf %lf %lf", &ra, &dec, &mag, &xraw, &yraw, &mraw); + RD_to_XY (&xobs, &yobs, ra, dec, &coords); + } else { + status = fscanf (f, "%lf %lf %lf", &xobs, &yobs, &mraw); + } + if (status == EOF) break; + if (READ_FLAGS) { + status = fscanf (f, "%x", &flags); + if (status == EOF) { + fprintf (stderr, "error: missing flag for last star?\n"); + exit (1); + } + Flag[Nstars] = flags; + } else { + Flag[Nstars] = 0; + } + + X[Nstars] = xobs; + Y[Nstars] = yobs; + M[Nstars] = mraw; + Xmax = MAX(Xmax, X[Nstars]); Ymax = MAX(Ymax, Y[Nstars]); @@ -94,6 +219,11 @@ REALLOCATE (Y, double, NSTARS); REALLOCATE (M, double, NSTARS); + REALLOCATE (Flag, unsigned int, NSTARS); } Nstars ++; + } + if (NX && NY) { + Xmax = NX; + Ymax = NY; } @@ -106,34 +236,22 @@ // XXX add minimum needed header fields - gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax + 50)); - gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax + 50)); + gfits_print (&header, "IMNAXIS1", "%d", 1, (int)(Xmax)); + gfits_print (&header, "IMNAXIS2", "%d", 1, (int)(Ymax)); gfits_modify (&header, "NSTARS", "%d", 1, Nstars); - gfits_modify (&header, "PHOTCODE", "%s", 1, "SIMTEST.r.Chip"); - gfits_modify (&header, "DATE-OBS", "%s", 1, date); - gfits_modify (&header, "UTC-OBS", "%s", 1, time); + gfits_modify (&header, "PHOTCODE", "%s", 1, photcode); + if (isfinite(mjd)) { + gfits_modify (&header, "MJD-OBS", "%lf", 1, mjd); + } else { + gfits_modify (&header, "DATE-OBS", "%s", 1, date); + gfits_modify (&header, "UTC-OBS", "%s", 1, time); + } gfits_modify (&header, "ZERO_PT", "%lf", 1, 25.0); - gfits_modify (&header, "EXPTIME", "%lf", 1, 1.0); - gfits_modify (&header, "AIRMASS", "%lf", 1, 1.0); + gfits_modify (&header, "EXPTIME", "%lf", 1, exptime); + gfits_modify (&header, "AIRMASS", "%lf", 1, airmass); gfits_modify (&header, "NASTRO", "%d", 1, 10); gfits_modify (&header, "IMAGEID", "%d", 1, (int)(1000*drand48())); gfits_modify (&header, "SOURCEID", "%d", 1, (int)(100*drand48())); - - /* bore site center guess */ - strcpy (coords.ctype, "DEC--TAN"); - coords.crval1 = RA; - coords.crval2 = DEC; - coords.crpix1 = 0.0; - coords.crpix2 = 0.0; - - coords.cdelt1 = 0.25/3600.0; - coords.cdelt2 = 0.25/3600.0; - - coords.pc1_1 = 1; - coords.pc2_2 = 1; - coords.pc1_2 = 0; - coords.pc2_1 = 0; - coords.Npolyterms = 1; PutCoords (&coords, &header); @@ -157,5 +275,5 @@ } if (!strcmp(type, "PS1_V2")) { - writeStars_PS1_V2 (&ftable, X, Y, M, Nstars); + writeStars_PS1_V2 (&ftable, X, Y, M, Flag, Nstars); found = TRUE; } @@ -211,5 +329,7 @@ A = time(NULL); - for (B = 0; A == time(NULL); B++); + // XXX this is expensive if called a lot (1 sec min) + // for (B = 0; A == time(NULL); B++); + B = A + 10000; srand48(B); @@ -273,7 +393,13 @@ fSN = 1.0 / sqrt(flux); - stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0); + stars[i].X = X[i]; + stars[i].Y = Y[i]; + stars[i].M = M[i]; + + if (ADDNOISE) { + stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M += fSN*rnd_gauss(0.0, 1.0); + } stars[i].dX = FX * fSN; @@ -312,7 +438,13 @@ fSN = 1.0 / sqrt(flux); - stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0); + stars[i].X = X[i]; + stars[i].Y = Y[i]; + stars[i].M = M[i]; + + if (ADDNOISE) { + stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M += fSN*rnd_gauss(0.0, 1.0); + } stars[i].dX = FX * fSN; @@ -354,7 +486,13 @@ fSN = 1.0 / sqrt(flux); - stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0); + stars[i].X = X[i]; + stars[i].Y = Y[i]; + stars[i].M = M[i]; + + if (ADDNOISE) { + stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M += fSN*rnd_gauss(0.0, 1.0); + } stars[i].dX = FX * fSN; @@ -380,5 +518,5 @@ } -void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, int Nstars) { +void writeStars_PS1_V2 (FTable *ftable, double *X, double *Y, double *M, unsigned int *Flag, int Nstars) { int i; @@ -396,8 +534,21 @@ fSN = 1.0 / sqrt(flux); - stars[i].X = X[i] + FX * fSN * rnd_gauss(0.0, 1.0); - stars[i].Y = Y[i] + FY * fSN * rnd_gauss(0.0, 1.0); - stars[i].M = M[i] + fSN*rnd_gauss(0.0, 1.0); - + stars[i].X = X[i]; + stars[i].Y = Y[i]; + stars[i].M = M[i]; + + if (ADDNOISE) { + stars[i].X += FX * fSN * rnd_gauss(0.0, 1.0); + stars[i].Y += FY * fSN * rnd_gauss(0.0, 1.0); + stars[i].M += fSN*rnd_gauss(0.0, 1.0); + } + + // randomly give poor PSFQF values + if ((BAD_PSFQF_FRAC > 0.0) && (drand48() < BAD_PSFQF_FRAC)) { + stars[i].psfQual = 0.25; + } else { + stars[i].psfQual = PSFQUAL; + } + stars[i].dX = FX * fSN; stars[i].dY = FY * fSN; @@ -413,7 +564,6 @@ stars[i].fy = FY; stars[i].df = DF; - stars[i].psfQual = PSFQUAL; stars[i].nFrames = 1; - stars[i].flags = FLAGS; + stars[i].flags = Flag[i]; } Index: /trunk/Ohana/src/addstar/src/parse_time.c =================================================================== --- /trunk/Ohana/src/addstar/src/parse_time.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/parse_time.c (revision 33653) @@ -145,2 +145,23 @@ } + +// ha/dec -> alt/az +int altaz (double *alt, double *az, double ha, double dec, double latitude) { + + double sind, sinh, cosh; + + sind = dSIN (dec) * dSIN (latitude) + dCOS (dec) * dCOS (ha) * dCOS (latitude); + *alt = DEG_RAD * asin (sind); + + sinh = - dCOS (dec) * dSIN (ha); + cosh = dSIN (dec) * dCOS (latitude) - dCOS (dec) * dCOS (ha) * dSIN (latitude); + + *az = DEG_RAD * atan2 (sinh, cosh); + + // I may need to use the parallactic angle to get the right plate rotation... + // sinh = -dCOS(az) * dSIN(alt) * dSIN(ha) * dSIN(latitude) + dSIN(az) * dSIN(alt) * dCOS(ha) - dSIN(ha) * dCOS(alt) * dCOS(latitude); + // cosh = -dSIN(az) * dSIN(ha) * dSIN(latitude) - dCOS(az) * dCOS(ha); + // rot = -DEG_RAD * atan2 (sinh, cosh); + return TRUE; +} + Index: /trunk/Ohana/src/addstar/src/resort_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_catalog.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/resort_catalog.c (revision 33653) @@ -28,7 +28,4 @@ return; } - -# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); } -# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } } // sort the measure Sequence based on the average Sequence entries Index: /trunk/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 33652) +++ /trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 33653) @@ -546,5 +546,9 @@ } -// define the parameters of a single sky projection center +// define the parameters of a projection centers for this ring +// dec : ~ center of ring in Dec +// dDEC : approximate height +// nring : number of cells generated for this ring +// format : guide to generate the filenames (c-type string format) SkyRectangle *sky_rectangle_ring (float dec, float dDEC, int *nring, char *format) { Index: /trunk/Ohana/src/addstar/test/relphot.flatcorr.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/relphot.flatcorr.dvo (revision 33653) +++ /trunk/Ohana/src/addstar/test/relphot.flatcorr.dvo (revision 33653) @@ -0,0 +1,563 @@ +# -*-sh-*- + +# relphot can apply a flat-field correction to each image (image->photom_map_id specifies the correction to use) +# we would like to test that the flat-field is correctly included in the analysis + +# define a flat-field correction +# write a fake flat-field correction file in the right format +# create a set of cmf files with the flat-field correction applied +# load cmfs into DVO +# apply the flat-field correction with setphot? + +# problems: +# * setphot only applies the very-specific ubercal flat-corrections (including specific MJD ranges for different corrections) +# * no other mechanism to tell DVO about a flat-field correction + +input tap.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +# set various global variables +macro init + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + + # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters + $zpt_nominal:g = 24.58 + $zpt_nominal:r = 24.80 + $zpt_nominal:i = 24.74 + $zpt_nominal:z = 24.26 + $zpt_nominal:y = 23.41 + + # klam is negative, so klam*(airmass - 1) increase the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + $klam_nominal:i = -0.04 + $klam_nominal:z = -0.03 + $klam_nominal:y = -0.03 + + # we have two sets of images: ubercaled and not-ubercaled + + # sequence to count images (only used in this function) + create tmpseq 0 9 + + $filt_uc:n = 9 + $filt_uc:0 = g + $filt_uc:1 = g + $filt_uc:2 = g + $filt_uc:3 = r + $filt_uc:4 = r + $filt_uc:5 = r + $filt_uc:6 = i + $filt_uc:7 = i + $filt_uc:8 = i + + # airmass slopes for these sequnece + set klam_uc = zero(tmpseq) + for i 0 tmpseq[] + klam_uc[$i] = $klam_nominal:$filt_uc:$i + end + + # mjd and zpt values for ubercal'ed imaged + set exptime_uc = 10.0 + zero(tmpseq) + set airmass_uc = 1.3 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + # note that K is defined as a negative value (is this sensible?) + set zpt_uc = 25.0 + 2.5*log(exptime_uc) + klam_uc*(airmass_uc - 1.0) + tmpseq*0.0050 - 0.0025 + + set mjd_uc = zero(zpt_uc) + mjd_uc[0] = 55000.01 + mjd_uc[1] = 55000.02 + mjd_uc[2] = 55000.03 + mjd_uc[3] = 55015.01 + mjd_uc[4] = 55015.02 + mjd_uc[5] = 55015.03 + mjd_uc[6] = 55025.01 + mjd_uc[7] = 55025.02 + mjd_uc[8] = 55025.03 + + # mjd and zpt values for not-ubercal'ed images + $filt_nc:n = 9 + $filt_nc:0 = g + $filt_nc:1 = g + $filt_nc:2 = g + $filt_nc:3 = r + $filt_nc:4 = r + $filt_nc:5 = r + $filt_nc:6 = i + $filt_nc:7 = i + $filt_nc:8 = i + + # airmass slopes for these sequnece + set klam_nc = zero(tmpseq) + for i 0 tmpseq[] + klam_nc[$i] = $klam_nominal:$filt_nc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 15.0 + zero(tmpseq) + set airmass_nc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - tmpseq*0.0050 + 0.0025 + + set mjd_nc = zero(zpt_nc) + mjd_nc[0] = 55000.11 + mjd_nc[1] = 55000.12 + mjd_nc[2] = 55000.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 +end + +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + exec rm -rf $catdir + + init + mkzptfile + mkstars 500 + + tapPLAN {((mjd_uc[]*4 + mjd_nc[]*4) + (mjd_uc[] + mjd_nc[]) + (mjd_uc[] + mjd_nc[])) / 3} + + # generate the basic images and check they were correctly ingested by dvo + mksequence $fileroot $catdir + for i 0 mjd_uc[] 3 + ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw + end + for i 0 mjd_nc[] 3 + ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # run setphot on the db and check that the images now match the expected values + exec setphot -update -ubercal testzpt.fits -D CATDIR $catdir + for i 0 mjd_uc[] 3 + ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] 3 + ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the db and check that the images now match the expected values + exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR catdir.test -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& tmp.log + for i 0 mjd_uc[] 3 + ckexposure catdir.test mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] 3 + ckexposure catdir.test mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + tapDONE +end + +# go testdata/t1 catdir.test +# ckexposure catdir.test mjd_uc[0] zpt_uc[0] raw : result is mean of 0.0 (since test corrects for ZPT_NOMINAL vs REAL) +# ckexposure catdir.test mjd_nc[0] zpt_nc[0] raw : result is mean of 0.0 +# setphot -update -ubercal testzpt.fits -D CATDIR catdir.test +# ckexposure catdir.test mjd_uc[0] zpt_uc[0] corr : result is mean of 0.0 (since image is corrected to ZPT_REAL) +# ckexposure catdir.test mjd_nc[0] zpt_nc[0] corr : result is mean of -0.42 (since image is NOT corrected to ZPT_REAL) +# relphot g -region 9.5 10.5 19.5 20.5 -update -D CATDIR catdir.test -D STAR_TOOFEW 1 -D SIGMA_LIM 0.075 -statmode MEAN +# dvo: ckexposure catdir.test mjd_uc[0] zpt_uc[0] corr +# mean is now -0.42 +# dvo: ckexposure catdir.test mjd_nc[0] zpt_nc[0] corr +# mean is now -0.42 +# ** setphot assigned the zero point of the uc exposure to an effective -0.4 (ie, negative clouds), relphot was not told to keep this so it dragged the negative clouds to 0.0 -- this seems to be the right thing to do + +# re-run setphot, then relphot with -keep-ubercal + +# the test case has only 2 exposures, so the STAR_TOOFEW of 1 is necessary. Also, the mag and sigma distributions are somewhat artificial +# relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR catdir.test -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update +# dvo: ckexposure catdir.test mjd_uc[0] zpt_uc[0] corr +# mean is now 0.0 +# dvo: ckexposure catdir.test mjd_nc[0] zpt_nc[0] corr +# mean is now 0.0 +# ** setphot correctly assigned the zero point of UC & relphot adjusted NC to match + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + if ($PLOT) + dev -n 1 + lim mag dm; clear; box; plot mag dm + end + + set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X)) + set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y)) + subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0) + subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1) + subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0) + subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1) + + # uncorrected values behave like this: + # mag_DVO = m_inst + zpt_nominal + 2.5*log(exptime) + K*(airmass - 1.0) + # mag_real = m_inst + zpt_real + cell_offset + # dm = mag_DVO - mag_real = zpt_nominal - zpt_real - cell_offset + # - zpt_nominal + zpt_real + cell_offset ~ 0.0 + # zpt_real (in this case) = 25.0 + 2.5*log(exptime) + K*(airmass - 1.0) + # (actually, it is the value in the vector 'zpt' for this entry + # zpt_nominal = 24.58 + + # setphot-corrected values behave like this: + # mag_DVO = mag_real (because m_inst has cell_offset applied) + # dm = mag_DVO - mag_real ~ 0.0 + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + + if ("$MODE" == "raw") + for ix 0 $NCELL_X + for iy 0 $NCELL_Y + # vstat -q dm$ix\$iy + # echo cell_off[$ix][$iy] {$MEDIAN - $ZPT_NOMINAL + $ZPT_REAL + cell_off[$ix][$iy]} $MEDIAN $MEAN $SIGMA + + set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + cell_off[$ix][$iy] + vstat -q dm_adjust + tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER" + end + end + return + end + + if ("$MODE" == "setphot_uc") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "setphot_uc ZP $MJD_IMAGE $FILTER" + return + end + + if ("$MODE" == "setphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc ZP $MJD_IMAGE $FILTER" + return + end + + if ("$MODE" == "relphot") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "relphot ZP $MJD_IMAGE $FILTER" + return + end + +end + +macro mksequence + if ($0 != 3) + echo "mksequence (fileroot) (catdir)" + break + end + + local i + + # we have defined a set of mjd and zero point values for ubercal'ed images + # for i 0 mjd_uc[] + for i 0 mjd_uc[] 3 + mkexposure $1.uc.$i $RA_CENTER $DEC_CENTER zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + + # we have a second set of mjd and zero point values for non-ubercal'ed images + # for i 0 mjd_nc[] + for i 0 mjd_nc[] 3 + mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (fileroot) (ra) (dec) (zpt) (exptime) (airmass) (mjd) (filter) (catdir)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR + + $ROOT = $1 + $RAo = $2 + $DECo = $3 + $ZPT = $4 + $EXPTIME = $5 + $AIRMASS = $6 + $MJD = $7 + $FILTER = $8 + $CATDIR = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + # XX ctimes -ref $MJD -var datetime + # XX substr $datetime 0 10 date + # XX substr $datetime 11 8 time + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + mkinput test.in.txt $ra $dec $ZPT + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + $options = $options -type PS1_V2 + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + + # the fake images have inconsistent ra,dec and airmass,sidtime values + echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& tmp.log + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X)) + set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X)) + + # eddie says: M_real = M_inst + zpt + ff_offset + # thus, stars_mi = stars_mag - zpt - ff_offset = stars_m - offset + set stars_mi = zero(stars_m) + for i 0 stars_mi[] + $nx = cell_x[$i] + $ny = cell_y[$i] + $offset = cell_off[$nx][$ny] + stars_mi[$i] = stars_m[$i] - $offset + end + + write -f "%10.6f %10.6f %7.3f %6.1f %6.1f %7.3f" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 15.0 + 4.0 * rnd(tmp) +end + +# create a fake ubercal-style zpt/flatcorr table with metadata +macro mkzptfile + + mcreate phu 0 0 + keyword phu NSEASON -wd $NSEASON + keyword phu NFILTER -wd $NFILTER + keyword phu NCHIP_X -wd $NCHIP_X + keyword phu NCHIP_Y -wd $NCHIP_Y + keyword phu NCELL_X -wd $NCELL_X + keyword phu NCELL_Y -wd $NCELL_Y + keyword phu CHIP_DX -wd $CHIP_DX + keyword phu CHIP_DY -wd $CHIP_DY + + # define the season boundaries + keyword phu TS0_0000 -wf 55000.0 + keyword phu TS1_0000 -wf 55010.0 + keyword phu TS0_0001 -wf 55010.0 + keyword phu TS1_0001 -wf 55020.0 + keyword phu TS0_0002 -wf 55020.0 + keyword phu TS1_0002 -wf 55030.0 + + # create a single flat-field correction map for a single chip + mcreate cell_off 2 2 + cell_off[0][0] = 0.01 + cell_off[1][0] = -0.01 + cell_off[0][1] = 0.02 + cell_off[1][1] = -0.02 + + # in the correction table, we have an image of the full exposure + # (NCHIP_X x NCHIP_Y) for each season, unwrapped as a single linear + # vector (with the last two missing -- how is that hard wired?) + + local ix iy ixc iyc ns + + # unroll the chip & cell level corrections into a single vector for a single image + delete -q image_offset + for iy 0 $NCHIP_Y + for iyc 0 $NCELL_Y + for ix 0 $NCHIP_X + for ixc 0 $NCELL_X + concat cell_off[$ixc][$iyc] image_offset + end + end + end + end + + # generate the full set of corrections for a single filter + delete -q offset + for ns 0 $NSEASON + concat image_offset offset + end + + # the output fits table uses the vector names for the field names, + # so set them here to the desired names + set mjd_obs = mjd_uc + set zp = zpt_uc + set resid = 0.02*rnd(zpt_uc) - 0.01 + + delete -q flatcorr + dimenup offset flatcorr {$NSEASON*$NCHIP_X*$NCHIP_Y*$NCELL_X*$NCELL_Y} 1 + + wd phu testzpt.fits + + # filter 1 + write testzpt.fits -fits ZPTS_1 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w g + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 2 (?) + write testzpt.fits -fits ZPTS_2 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w r + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 3 + write testzpt.fits -fits ZPTS_3 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w i + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits +end + +# create a populated catdir with a couple of cmf files +macro test.relphot + if ($0 != 3) + echo "test.relphot (cmftype) (dvotype)" + break + end + + tapPLAN 4 + + exec rm -rf catdir.test + + $RA = 10.0 + $DEC = 20.0 + + for i 0 $offset:n + mkinput test.in.txt $offset:$i + exec mkcmf test.in.txt test.cmf -date 2008/1/1 -time $i\:00:00 -radec $RA $DEC -type $1 + exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 + end + + exec relphot -D CATDIR catdir.test r -region {$RA-1.0} {$RA+1.0} {$DEC-1.0} {$DEC+1.0} -update -nloop 10.0 >& /dev/null + + catdir catdir.test + skyregion {$RA-1.0} {$RA+1.0} {$DEC-1.0} {$DEC+1.0} + + imextract time Mcal + + for i 0 $offset:n + tapOK {abs(Mcal[$i] - Mcal[0] - $offset:$i) < 0.001} "Mcal $i" + end + + exec rm test.cmf test.in.txt + exec rm -r catdir.test + + tapDONE +end + +# running setphot: setphot testzpt.fits -update -ubercal -D CATDIR catdir.test Index: /trunk/Ohana/src/addstar/test/relphot.parallel.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/relphot.parallel.dvo (revision 33653) +++ /trunk/Ohana/src/addstar/test/relphot.parallel.dvo (revision 33653) @@ -0,0 +1,446 @@ +# -*-sh-*- + +# test suite for the parallel analysis mode for relphot +# this does not include / test flat-field correction options + +input tap.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +# set various global variables +macro init + echo "starting init..." + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + if (not($?VERBOSE)) set VERBOSE = 0 + + # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters + $zpt_nominal:g = 24.58 + $zpt_nominal:r = 24.80 + $zpt_nominal:i = 24.74 + $zpt_nominal:z = 24.26 + $zpt_nominal:y = 23.41 + + # klam is negative, so klam*(airmass - 1) increase the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + $klam_nominal:i = -0.04 + $klam_nominal:z = -0.03 + $klam_nominal:y = -0.03 + + # we have two sets of images: ubercaled and not-ubercaled + + # sequence to count images (only used in this function) + create tmpseq 0 9 + + # mjd and zpt values for not-ubercal'ed images + $filt_nc:n = 9 + create filtN_nc 0 $filt_nc:n + $filt_nc:0 = g; filtN_nc[0] = 1 + $filt_nc:1 = g; filtN_nc[1] = 1 + $filt_nc:2 = g; filtN_nc[2] = 1 + $filt_nc:3 = r; filtN_nc[3] = 2 + $filt_nc:4 = r; filtN_nc[4] = 2 + $filt_nc:5 = r; filtN_nc[5] = 2 + $filt_nc:6 = i; filtN_nc[6] = 3 + $filt_nc:7 = i; filtN_nc[7] = 3 + $filt_nc:8 = i; filtN_nc[8] = 3 + + # airmass slopes for these sequnece + set klam_nc = zero(tmpseq) + for i 0 tmpseq[] + klam_nc[$i] = $klam_nominal:$filt_nc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 15.0 + zero(tmpseq) + set airmass_nc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - tmpseq*0.0050 + 0.0025 + + # without an ubercal or other tie, relphot should set zpts to the mean of the zpts for a filter + subset tmp = zpt_nc where filtN_nc == 1; vstat -q tmp; set zpt_nc_g = $MEAN + subset tmp = zpt_nc where filtN_nc == 2; vstat -q tmp; set zpt_nc_r = $MEAN + subset tmp = zpt_nc where filtN_nc == 3; vstat -q tmp; set zpt_nc_i = $MEAN + + set mjd_nc = zero(zpt_nc) + mjd_nc[0] = 55000.11 + mjd_nc[1] = 55000.12 + mjd_nc[2] = 55000.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 +end + +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + exec rm -rf $catdir + + init + mkstars 500 + + tapPLAN {((mjd_nc[]*4 + mjd_nc[]) + mjd_nc[])} + + # generate the basic images and check they were correctly ingested by dvo + mksequence $fileroot $catdir + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # XXX tmp + exec rsync -aub $catdir/ $catdir.raw/ + + # parallel-ize the database + mkparallel $catdir + + # run relphot on the serial db and check that the images now match the expected values + exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + + for i 0 mjd_nc[] + # this is a little ill-defined : if we run relphot with no tied + # down data, we should convert on a solution with zero points + # matching the mean zp of photometric data, which we calculate in + # 'init' + ckexposure $catdir mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + end + + # run relphot on the parallel db and check that the images now match the expected values + exec relphot -parallel g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir.p0 -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.p0 + + exec dvodist -in $catdir.p0 >& log.dvodist.in + + for i 0 mjd_nc[] + # this is a little ill-defined : if we run relphot with no tied + # down data, we should convert on a solution with zero points + # matching the mean zp of photometric data, which we calculate in + # 'init' + ckexposure $catdir.p0 mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + end + + tapDONE +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + if (not(ra[])) + echo "no matching data for MJD = $MJD_IMAGE" + if ("$MODE" == "raw") + tapSKIP 4 + else + tapSKIP 1 + end + return + end + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + if (ra[] != stars_ra_m[]) + echo "failed to match stars in catdir to stars in memory" + break + end + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + if ($PLOT) + dev -n 0 + cplot -pt 7 -c red ra dec + dev -n 1 + lim mag dm; clear; box; plot mag dm + end + + set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X)) + set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y)) + subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0) + subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1) + subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0) + subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1) + + # uncorrected values behave like this: + # mag_DVO = m_inst + zpt_nominal + 2.5*log(exptime) + K*(airmass - 1.0) + # mag_real = m_inst + zpt_real + cell_offset + # dm = mag_DVO - mag_real = zpt_nominal - zpt_real - cell_offset + # - zpt_nominal + zpt_real + cell_offset ~ 0.0 + # zpt_real (in this case) = 25.0 + 2.5*log(exptime) + K*(airmass - 1.0) + # (actually, it is the value in the vector 'zpt' for this entry + # zpt_nominal = 24.58 + + # setphot-corrected values behave like this: + # mag_DVO = mag_real (because m_inst has cell_offset applied) + # dm = mag_DVO - mag_real ~ 0.0 + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + if ($VERBOSE) echo "$ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)" + + if ("$MODE" == "raw") + for ix 0 $NCELL_X + for iy 0 $NCELL_Y + # vstat -q dm$ix\$iy + # echo cell_off[$ix][$iy] {$MEDIAN - $ZPT_NOMINAL + $ZPT_REAL + cell_off[$ix][$iy]} $MEDIAN $MEAN $SIGMA + + set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + vstat -q dm_adjust + tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + end + end + return + end + + if ("$MODE" == "setphot_uc") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "setphot_uc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "setphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "relphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "relphot ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + +end + +macro mksequence + if ($0 != 3) + echo "mksequence (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (fileroot) (ra) (dec) (zpt) (exptime) (airmass) (mjd) (filter) (catdir)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR + + $ROOT = $1 + $RAo = $2 + $DECo = $3 + $ZPT = $4 + $EXPTIME = $5 + $AIRMASS = $6 + $MJD = $7 + $FILTER = $8 + $CATDIR = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + if ($VERBOSE) echo mkinput test.in.txt $ra $dec $ZPT + mkinput test.in.txt $ra $dec $ZPT + # costs 7/36 sec / chip + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + $options = $options -type PS1_V2 + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + if ($VERBOSE) echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + # costs 8/36 sec / chip + + # the fake images have inconsistent ra,dec and airmass,sidtime values + if ($VERBOSE) echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& /dev/null + # costs 10/36 sec / chip + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X)) + set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X)) + + if (1) + # FLAT-FIELD CORRECTION + # eddie says: M_real = M_inst + zpt + ff_offset + # thus, stars_mi = stars_mag - zpt - ff_offset = stars_m - offset + set stars_mi = zero(stars_m) + for i 0 stars_mi[] + $nx = cell_x[$i] + $ny = cell_y[$i] + # $offset = cell_off[$nx][$ny] + $offset = 0.0 + stars_mi[$i] = stars_m[$i] - $offset + end + else + set stars_mi = stars_m + end + + # delete the output file before writing + exec rm -f $1 + write -f "%10.6f %10.6f %7.3f %6.1f %6.1f %7.3f" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + echo "starting mkstars..." + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 15.0 + 4.0 * rnd(tmp) +end + +# copy the catdir to catdir.p0 and parallelize to local directories catdir.d1 - catdir.d3 +macro mkparallel + if ($0 != 2) + echo "USAGE: mkparallel (catdir)" + break + end + + local hostname catdir + + $hostname = `hostname -s` + + $tmp1 = `dirname $1` + $tmp2 = `basename $1` + $catdir = $tmp1/$tmp2 + + exec rsync -auv $catdir/ $catdir.p0/ >& log.rsync + mkdir $catdir.d1 + mkdir $catdir.d2 + mkdir $catdir.d3 + + exec rm -f $catdir.p0/HostTable.dat + exec echo "1 $hostname $catdir.d1" > $catdir.p0/HostTable.dat + exec echo "2 $hostname $catdir.d2" >> $catdir.p0/HostTable.dat + exec echo "3 $hostname $catdir.d3" >> $catdir.p0/HostTable.dat + + exec dvodist -out $catdir.p0 >& log.dvodist.out +end + Index: /trunk/Ohana/src/addstar/test/relphot.reject.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/relphot.reject.dvo (revision 33653) +++ /trunk/Ohana/src/addstar/test/relphot.reject.dvo (revision 33653) @@ -0,0 +1,645 @@ +# -*-sh-*- + +# test suite for the parallel analysis mode for relphot & setphot +# in this test suite, we add in sources with outlier measurements. these are supposed to trigger the various rejection flags in relphot. does it work? + +# measurement rejections that we perform: +# measurement is an outlier for the star +# measurement is flagged as bad by the photometry +# measurement is flagged as poor by the photometry +# Area & Time (ignore for now -- probably not used by GPC1) + +# also: +# bad images: chisq to high, zpt scatter too high +# bad objects : chisq to high, too few measurements? + +input tap.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +# set various global variables +macro init + echo "starting init..." + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + $DEFECT_FRAC = 0.03 + $OFFSET_FRAC_UC = 0.0 + $OFFSET_FRAC_NC = 0.5 + if (not($?VERBOSE)) set VERBOSE = 0 + + # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters + $zpt_nominal:g = 24.58 + $zpt_nominal:r = 24.80 + $zpt_nominal:i = 24.74 + $zpt_nominal:z = 24.26 + $zpt_nominal:y = 23.41 + + # klam is negative, so klam*(airmass - 1) increase the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + $klam_nominal:i = -0.04 + $klam_nominal:z = -0.03 + $klam_nominal:y = -0.03 + + # we have two sets of images: ubercaled and not-ubercaled + + # sequence to count images (only used in this function) + # XXX TEST create tmpseq 0 9 + create tmpseq 0 9 + + ##### mjd and zpt values for ubercal'ed images + + $filt_uc:n = tmpseq[] + create filtN_uc 0 $filt_uc:n + $filt_uc:0 = g; filtN_uc[0] = 1 + $filt_uc:1 = g; filtN_uc[1] = 1 + $filt_uc:2 = g; filtN_uc[2] = 1 + $filt_uc:3 = g; filtN_uc[3] = 1 + $filt_uc:4 = g; filtN_uc[4] = 1 + $filt_uc:5 = g; filtN_uc[5] = 1 + $filt_uc:6 = g; filtN_uc[6] = 1 + $filt_uc:7 = g; filtN_uc[7] = 1 + $filt_uc:8 = g; filtN_uc[8] = 1 + + # airmass slopes for this sequnece + set klam_uc = zero(tmpseq) + for i 0 tmpseq[] + klam_uc[$i] = $klam_nominal:$filt_uc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_uc = 15.0 + zero(tmpseq) + set airmass_uc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_uc = 25.0 + 2.5*log(exptime_uc) + klam_uc*(airmass_uc - 1.0) + tmpseq*0.0050 - 0.0025 + + set mjd_uc = zero(zpt_uc) + mjd_uc[0] = 55000.01 + mjd_uc[1] = 55000.02 + mjd_uc[2] = 55000.03 + mjd_uc[3] = 55015.01 + mjd_uc[4] = 55015.02 + mjd_uc[5] = 55015.03 + mjd_uc[6] = 55025.01 + mjd_uc[7] = 55025.02 + mjd_uc[8] = 55025.03 + + ##### mjd and zpt values for not-ubercal'ed images + + $filt_nc:n = tmpseq[] + create filtN_nc 0 $filt_nc:n + $filt_nc:0 = g; filtN_nc[0] = 1 + $filt_nc:1 = g; filtN_nc[1] = 1 + $filt_nc:2 = g; filtN_nc[2] = 1 + $filt_nc:3 = g; filtN_nc[3] = 1 + $filt_nc:4 = g; filtN_nc[4] = 1 + $filt_nc:5 = g; filtN_nc[5] = 1 + $filt_nc:6 = g; filtN_nc[6] = 1 + $filt_nc:7 = g; filtN_nc[7] = 1 + $filt_nc:8 = g; filtN_nc[8] = 1 + + # airmass slopes for these sequnece + set klam_nc = zero(tmpseq) + for i 0 tmpseq[] + klam_nc[$i] = $klam_nominal:$filt_nc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 15.0 + zero(tmpseq) + set airmass_nc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - tmpseq*0.0050 + 0.0025 + + # with an ubercal tie, relphot should set zpts to the inserted values above + set mjd_nc = zero(zpt_nc) + mjd_nc[0] = 55000.11 + mjd_nc[1] = 55000.12 + mjd_nc[2] = 55000.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 +end + +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + exec rm -rf $catdir + + init + mkzptfile + mkstars 10000 + + tapPLAN {(((mjd_uc[] + mjd_nc[])*4 + (mjd_uc[] + mjd_nc[])) + (mjd_uc[] + mjd_nc[]))} + + # generate the basic images and check they were correctly ingested by dvo + mksequence.drift $fileroot $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # run setphot on the serial db and check that the images now match the expected values + exec setphot -reset -update -ubercal testzpt.fits -D CATDIR $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the serial db and check that the images now match the expected values + # if ($VERBOSE) echo relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + if ($VERBOSE) echo relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0 + exec relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0 + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + tapDONE +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + if (not(ra[])) + echo "no matching data for MJD = $MJD_IMAGE" + if ("$MODE" == "raw") + tapSKIP 4 + else + tapSKIP 1 + end + return + end + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + if (ra[] != stars_ra_m[]) + echo "failed to match stars in catdir to stars in memory" + break + end + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + if ($PLOT) + dev -n 0 + cplot -pt 7 -c red ra dec + dev -n 1 + lim mag dm; clear; box; plot mag dm + end + + set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X)) + set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y)) + subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0) + subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1) + subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0) + subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1) + + # uncorrected values behave like this: + # mag_DVO = m_inst + zpt_nominal + 2.5*log(exptime) + K*(airmass - 1.0) + # mag_real = m_inst + zpt_real + cell_offset + # dm = mag_DVO - mag_real = zpt_nominal - zpt_real - cell_offset + # - zpt_nominal + zpt_real + cell_offset ~ 0.0 + # zpt_real (in this case) = 25.0 + 2.5*log(exptime) + K*(airmass - 1.0) + # (actually, it is the value in the vector 'zpt' for this entry + # zpt_nominal = 24.58 + + # setphot-corrected values behave like this: + # mag_DVO = mag_real (because m_inst has cell_offset applied) + # dm = mag_DVO - mag_real ~ 0.0 + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + if ($VERBOSE) echo "$ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)" + + if ("$MODE" == "raw") + for ix 0 $NCELL_X + for iy 0 $NCELL_Y + # vstat -q dm$ix\$iy + # echo cell_off[$ix][$iy] {$MEDIAN - $ZPT_NOMINAL + $ZPT_REAL + cell_off[$ix][$iy]} $MEDIAN $MEAN $SIGMA + + set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + vstat -q dm_adjust + tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + end + end + return + end + + if ("$MODE" == "setphot_uc") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "setphot_uc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "setphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "relphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "relphot ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + +end + +macro mksequence.drift + if ($0 != 3) + echo "mksequence.drift (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + mkexposure $1.uc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_UC/3600.0} zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_NC/3600.0} zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mksequence.1pt + if ($0 != 3) + echo "mksequence (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + mkexposure $1.uc.$i $RA_CENTER $DEC_CENTER zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (fileroot) (ra) (dec) (zpt) (exptime) (airmass) (mjd) (filter) (catdir)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR + + $ROOT = $1 + $RAo = $2 + $DECo = $3 + $ZPT = $4 + $EXPTIME = $5 + $AIRMASS = $6 + $MJD = $7 + $FILTER = $8 + $CATDIR = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + if ($VERBOSE) echo mkinput test.in.txt $ra $dec $ZPT + mkinput test.in.txt $ra $dec $ZPT + # costs 7/36 sec / chip + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + $options = $options -type PS1_V2 + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -flags + $options = $options -bad-psfqf-frac $DEFECT_FRAC + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + if ($VERBOSE) echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + # costs 8/36 sec / chip + + # the fake images have inconsistent ra,dec and airmass,sidtime values + if ($VERBOSE) echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& /dev/null + # costs 10/36 sec / chip + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X)) + set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X)) + + if (1) + # FLAT-FIELD CORRECTION + # eddie says: M_real = M_inst + zpt + ff_offset + # thus, stars_mi = stars_mag - zpt - ff_offset = stars_m - offset + set stars_mi = zero(stars_m) + for i 0 stars_mi[] + $nx = cell_x[$i] + $ny = cell_y[$i] + # $offset = cell_off[$nx][$ny] + $offset = 0.0 + stars_mi[$i] = stars_m[$i] - $offset + end + else + set stars_mi = stars_m + end + + # create bad detections / detections to ignore or avoid, save to a file? + + # create 3 types of things to ignore / avoid + # * badPhot (deviate the photometry) + # * badFlag (set a bad flag bit) + # * poorFlag (set a bad flag bit) + set badPhot = rnd(stars_r) < $DEFECT_FRAC + set badFlag = rnd(stars_r) < $DEFECT_FRAC + set poorFlag = rnd(stars_r) < $DEFECT_FRAC + + # create deviant photometry (1%) + # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated. This is a simple way to pass the info + set stars_mi = stars_mi + 0.2*badPhot + set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot + + # bad flags that we sould respect: + # PM_SOURCE_MODE_FAIL = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero) + # PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation + # PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF + # PM_SOURCE_MODE_DEFECT = 0x00000800, ///< Source is thought to be a defect + # PM_SOURCE_MODE_SATURATED = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail) + # PM_SOURCE_MODE_CR_LIMIT = 0x00002000, ///< Source has crNsigma above limit + # PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000, ///< could not measure the moments + # PM_SOURCE_MODE_SKY_FAILURE = 0x00010000, ///< could not measure the local sky + # PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000, ///< could not measure the local sky variance + # PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + # TOTAL_BAD = 0x1003bc88 + + # poor flags that we should respect: + # PM_SOURCE_MODE_POOR = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) + # PM_SOURCE_MODE_PAIR = 0x00000020, ///< Source fitted with a double psf + # PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources + # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N + # PM_SOURCE_MODE_BLEND_FIT = 0x00400000, ///< source was fitted as a blend + # PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + # PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + # PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip + # TOTAL_POOR = 0xe0440130 + + # delete the output file before writing + exec rm -f $1 + write -f "%10.6f %10.6f %7.3f %6.1f %6.1f %7.3f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + echo "starting mkstars..." + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE * (1 + tmpseq[]*$OFFSET_FRAC_NC) + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 15.0 + 4.0 * rnd(tmp) +end + +# copy the catdir to catdir.p0 and parallelize to local directories catdir.d1 - catdir.d3 +macro mkparallel + if ($0 != 2) + echo "USAGE: mkparallel (catdir)" + break + end + + local hostname catdir + + $hostname = `hostname -s` + + $tmp1 = `dirname $1` + $tmp2 = `basename $1` + $catdir = $tmp1/$tmp2 + + exec rsync -auv $catdir/ $catdir.p0/ >& log.rsync + mkdir $catdir.d1 + mkdir $catdir.d2 + mkdir $catdir.d3 + + exec rm -f $catdir.p0/HostTable.dat + exec echo "1 $hostname $catdir.d1" > $catdir.p0/HostTable.dat + exec echo "2 $hostname $catdir.d2" >> $catdir.p0/HostTable.dat + exec echo "3 $hostname $catdir.d3" >> $catdir.p0/HostTable.dat + + exec dvodist -out $catdir.p0 >& log.dvodist.out +end + +# create a fake ubercal-style zpt/flatcorr table with metadata +macro mkzptfile + + mcreate phu 0 0 + keyword phu NSEASON -wd $NSEASON + keyword phu NFILTER -wd $NFILTER + keyword phu NCHIP_X -wd $NCHIP_X + keyword phu NCHIP_Y -wd $NCHIP_Y + keyword phu NCELL_X -wd $NCELL_X + keyword phu NCELL_Y -wd $NCELL_Y + keyword phu CHIP_DX -wd $CHIP_DX + keyword phu CHIP_DY -wd $CHIP_DY + + # define the season boundaries + keyword phu TS0_0000 -wf 55000.0 + keyword phu TS1_0000 -wf 55010.0 + keyword phu TS0_0001 -wf 55010.0 + keyword phu TS1_0001 -wf 55020.0 + keyword phu TS0_0002 -wf 55020.0 + keyword phu TS1_0002 -wf 55030.0 + + # create a single flat-field correction map for a single chip + mcreate cell_off 2 2 + cell_off[0][0] = 0.01 + cell_off[1][0] = -0.01 + cell_off[0][1] = 0.02 + cell_off[1][1] = -0.02 + + # in the correction table, we have an image of the full exposure + # (NCHIP_X x NCHIP_Y) for each season, unwrapped as a single linear + # vector (with the last two missing -- how is that hard wired?) + + local ix iy ixc iyc ns + + # unroll the chip & cell level corrections into a single vector for a single image + delete -q image_offset + for iy 0 $NCHIP_Y + for iyc 0 $NCELL_Y + for ix 0 $NCHIP_X + for ixc 0 $NCELL_X + concat cell_off[$ixc][$iyc] image_offset + end + end + end + end + + # generate the full set of corrections for a single filter + delete -q offset + for ns 0 $NSEASON + concat image_offset offset + end + + # the output fits table uses the vector names for the field names, + # so set them here to the desired names + set mjd_obs = mjd_uc + set zp = zpt_uc + set resid = 0.02*rnd(zpt_uc) - 0.01 + + delete -q flatcorr + dimenup offset flatcorr {$NSEASON*$NCHIP_X*$NCHIP_Y*$NCELL_X*$NCELL_Y} 1 + + wd phu testzpt.fits + + # filter 1 + write testzpt.fits -fits ZPTS_1 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w g + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 2 (?) + write testzpt.fits -fits ZPTS_2 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w r + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 3 + write testzpt.fits -fits ZPTS_3 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w i + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits +end + Index: /trunk/Ohana/src/addstar/test/relphot.setphot.parallel.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/relphot.setphot.parallel.dvo (revision 33653) +++ /trunk/Ohana/src/addstar/test/relphot.setphot.parallel.dvo (revision 33653) @@ -0,0 +1,639 @@ +# -*-sh-*- + +# test suite for the parallel analysis mode for relphot & setphot + +input tap.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +# set various global variables +macro init + echo "starting init..." + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + if (not($?VERBOSE)) set VERBOSE = 0 + + # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters + $zpt_nominal:g = 24.58 + $zpt_nominal:r = 24.80 + $zpt_nominal:i = 24.74 + $zpt_nominal:z = 24.26 + $zpt_nominal:y = 23.41 + + # klam is negative, so klam*(airmass - 1) increase the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + $klam_nominal:i = -0.04 + $klam_nominal:z = -0.03 + $klam_nominal:y = -0.03 + + # we have two sets of images: ubercaled and not-ubercaled + + # sequence to count images (only used in this function) + create tmpseq 0 9 + + ##### mjd and zpt values for ubercal'ed images + + $filt_uc:n = 9 + create filtN_uc 0 $filt_uc:n + $filt_uc:0 = g; filtN_uc[0] = 1 + $filt_uc:1 = g; filtN_uc[1] = 1 + $filt_uc:2 = g; filtN_uc[2] = 1 + $filt_uc:3 = r; filtN_uc[3] = 2 + $filt_uc:4 = r; filtN_uc[4] = 2 + $filt_uc:5 = r; filtN_uc[5] = 2 + $filt_uc:6 = i; filtN_uc[6] = 3 + $filt_uc:7 = i; filtN_uc[7] = 3 + $filt_uc:8 = i; filtN_uc[8] = 3 + + # airmass slopes for this sequnece + set klam_uc = zero(tmpseq) + for i 0 tmpseq[] + klam_uc[$i] = $klam_nominal:$filt_uc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_uc = 15.0 + zero(tmpseq) + set airmass_uc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_uc = 25.0 + 2.5*log(exptime_uc) + klam_uc*(airmass_uc - 1.0) + tmpseq*0.0050 - 0.0025 + + set mjd_uc = zero(zpt_uc) + mjd_uc[0] = 55000.01 + mjd_uc[1] = 55000.02 + mjd_uc[2] = 55000.03 + mjd_uc[3] = 55015.01 + mjd_uc[4] = 55015.02 + mjd_uc[5] = 55015.03 + mjd_uc[6] = 55025.01 + mjd_uc[7] = 55025.02 + mjd_uc[8] = 55025.03 + + ##### mjd and zpt values for not-ubercal'ed images + + $filt_nc:n = 9 + create filtN_nc 0 $filt_nc:n + $filt_nc:0 = g; filtN_nc[0] = 1 + $filt_nc:1 = g; filtN_nc[1] = 1 + $filt_nc:2 = g; filtN_nc[2] = 1 + $filt_nc:3 = r; filtN_nc[3] = 2 + $filt_nc:4 = r; filtN_nc[4] = 2 + $filt_nc:5 = r; filtN_nc[5] = 2 + $filt_nc:6 = i; filtN_nc[6] = 3 + $filt_nc:7 = i; filtN_nc[7] = 3 + $filt_nc:8 = i; filtN_nc[8] = 3 + + # airmass slopes for these sequnece + set klam_nc = zero(tmpseq) + for i 0 tmpseq[] + klam_nc[$i] = $klam_nominal:$filt_nc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 15.0 + zero(tmpseq) + set airmass_nc = 1.6 + zero(tmpseq) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - tmpseq*0.0050 + 0.0025 + + # with an ubercal tie, relphot should set zpts to the inserted values above + set mjd_nc = zero(zpt_nc) + mjd_nc[0] = 55000.11 + mjd_nc[1] = 55000.12 + mjd_nc[2] = 55000.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 +end + +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + exec rm -rf $catdir + + init + mkzptfile + mkstars 500 + + tapPLAN {(((mjd_uc[] + mjd_nc[])*4 + (mjd_uc[] + mjd_nc[])) + (mjd_uc[] + mjd_nc[]))} + + # generate the basic images and check they were correctly ingested by dvo + mksequence $fileroot $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # parallel-ize the database (goes to $catdir.p0) + mkparallel $catdir + + ### SERIAL database tests + + # run setphot on the serial db and check that the images now match the expected values + exec setphot -reset -update -ubercal testzpt.fits -D CATDIR $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the serial db and check that the images now match the expected values + exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + ### PARALLEL database tests + + # run setphot on the parallel db and check that the images now match the expected values + exec setphot -reset -parallel -update -ubercal testzpt.fits -D CATDIR $catdir.p0 + + # bring the detections back to the local db + # exec dvodist -in $catdir.p0 >& log.dvodist.in + # XXX this should now not be needed + + # check the results + for i 0 mjd_uc[] + ckexposure $catdir.p0 mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir.p0 mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the parallel db and check that the images now match the expected values + exec relphot -parallel g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir.p0 -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.p0 + + for i 0 mjd_uc[] + ckexposure $catdir.p0 mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir.p0 mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + tapDONE +end + +macro qt + + ## SER + + # run setphot on the serial db and check that the images now match the expected values + exec setphot -reset -update -ubercal testzpt.fits -D CATDIR $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the serial db and check that the images now match the expected values + exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + ## PAR + + # run setphot on the parallel db and check that the images now match the expected values + exec setphot -reset -parallel -update -ubercal testzpt.fits -D CATDIR $catdir.p0 + + # check the results + for i 0 mjd_uc[] + ckexposure $catdir.p0 mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir.p0 mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the parallel db and check that the images now match the expected values + exec relphot -parallel g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir.p0 -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.p0 + + for i 0 mjd_uc[] + ckexposure $catdir.p0 mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir.p0 mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + if (not(ra[])) + echo "no matching data for MJD = $MJD_IMAGE" + if ("$MODE" == "raw") + tapSKIP 4 + else + tapSKIP 1 + end + return + end + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + if (ra[] != stars_ra_m[]) + echo "failed to match stars in catdir to stars in memory" + break + end + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + if ($PLOT) + dev -n 0 + cplot -pt 7 -c red ra dec + dev -n 1 + lim mag dm; clear; box; plot mag dm + end + + set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X)) + set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y)) + subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0) + subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1) + subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0) + subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1) + + # uncorrected values behave like this: + # mag_DVO = m_inst + zpt_nominal + 2.5*log(exptime) + K*(airmass - 1.0) + # mag_real = m_inst + zpt_real + cell_offset + # dm = mag_DVO - mag_real = zpt_nominal - zpt_real - cell_offset + # - zpt_nominal + zpt_real + cell_offset ~ 0.0 + # zpt_real (in this case) = 25.0 + 2.5*log(exptime) + K*(airmass - 1.0) + # (actually, it is the value in the vector 'zpt' for this entry + # zpt_nominal = 24.58 + + # setphot-corrected values behave like this: + # mag_DVO = mag_real (because m_inst has cell_offset applied) + # dm = mag_DVO - mag_real ~ 0.0 + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + if ($VERBOSE) echo "$ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)" + + if ("$MODE" == "raw") + for ix 0 $NCELL_X + for iy 0 $NCELL_Y + # vstat -q dm$ix\$iy + # echo cell_off[$ix][$iy] {$MEDIAN - $ZPT_NOMINAL + $ZPT_REAL + cell_off[$ix][$iy]} $MEDIAN $MEAN $SIGMA + + set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + vstat -q dm_adjust + tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + end + end + return + end + + if ("$MODE" == "setphot_uc") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "setphot_uc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "setphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "relphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "relphot ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + +end + +macro mksequence + if ($0 != 3) + echo "mksequence (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + mkexposure $1.uc.$i $RA_CENTER $DEC_CENTER zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (fileroot) (ra) (dec) (zpt) (exptime) (airmass) (mjd) (filter) (catdir)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR + + $ROOT = $1 + $RAo = $2 + $DECo = $3 + $ZPT = $4 + $EXPTIME = $5 + $AIRMASS = $6 + $MJD = $7 + $FILTER = $8 + $CATDIR = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + if ($VERBOSE) echo mkinput test.in.txt $ra $dec $ZPT + mkinput test.in.txt $ra $dec $ZPT + # costs 7/36 sec / chip + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + $options = $options -type PS1_V2 + $options = $options -coords + $options = $options -photcode GPC1.$FILTER.XY$ix\$iy + $options = $options -no-noise + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + if ($VERBOSE) echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + # costs 8/36 sec / chip + + # the fake images have inconsistent ra,dec and airmass,sidtime values + if ($VERBOSE) echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& /dev/null + # costs 10/36 sec / chip + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X)) + set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X)) + + if (1) + # FLAT-FIELD CORRECTION + # eddie says: M_real = M_inst + zpt + ff_offset + # thus, stars_mi = stars_mag - zpt - ff_offset = stars_m - offset + set stars_mi = zero(stars_m) + for i 0 stars_mi[] + $nx = cell_x[$i] + $ny = cell_y[$i] + # $offset = cell_off[$nx][$ny] + $offset = 0.0 + stars_mi[$i] = stars_m[$i] - $offset + end + else + set stars_mi = stars_m + end + + # delete the output file before writing + exec rm -f $1 + write -f "%10.6f %10.6f %7.3f %6.1f %6.1f %7.3f" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + echo "starting mkstars..." + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 15.0 + 4.0 * rnd(tmp) +end + +# copy the catdir to catdir.p0 and parallelize to local directories catdir.d1 - catdir.d3 +macro mkparallel + if ($0 != 2) + echo "USAGE: mkparallel (catdir)" + break + end + + local hostname catdir + + $hostname = `hostname -s` + + $tmp1 = `dirname $1` + $tmp2 = `basename $1` + $catdir = $tmp1/$tmp2 + + exec rsync -auv $catdir/ $catdir.p0/ >& log.rsync + mkdir $catdir.d1 + mkdir $catdir.d2 + mkdir $catdir.d3 + + exec rm -f $catdir.p0/HostTable.dat + exec echo "1 $hostname $catdir.d1" > $catdir.p0/HostTable.dat + exec echo "2 $hostname $catdir.d2" >> $catdir.p0/HostTable.dat + exec echo "3 $hostname $catdir.d3" >> $catdir.p0/HostTable.dat + + exec dvodist -out $catdir.p0 >& log.dvodist.out +end + +# create a fake ubercal-style zpt/flatcorr table with metadata +macro mkzptfile + + mcreate phu 0 0 + keyword phu NSEASON -wd $NSEASON + keyword phu NFILTER -wd $NFILTER + keyword phu NCHIP_X -wd $NCHIP_X + keyword phu NCHIP_Y -wd $NCHIP_Y + keyword phu NCELL_X -wd $NCELL_X + keyword phu NCELL_Y -wd $NCELL_Y + keyword phu CHIP_DX -wd $CHIP_DX + keyword phu CHIP_DY -wd $CHIP_DY + + # define the season boundaries + keyword phu TS0_0000 -wf 55000.0 + keyword phu TS1_0000 -wf 55010.0 + keyword phu TS0_0001 -wf 55010.0 + keyword phu TS1_0001 -wf 55020.0 + keyword phu TS0_0002 -wf 55020.0 + keyword phu TS1_0002 -wf 55030.0 + + # create a single flat-field correction map for a single chip + mcreate cell_off 2 2 + cell_off[0][0] = 0.01 + cell_off[1][0] = -0.01 + cell_off[0][1] = 0.02 + cell_off[1][1] = -0.02 + + # in the correction table, we have an image of the full exposure + # (NCHIP_X x NCHIP_Y) for each season, unwrapped as a single linear + # vector (with the last two missing -- how is that hard wired?) + + local ix iy ixc iyc ns + + # unroll the chip & cell level corrections into a single vector for a single image + delete -q image_offset + for iy 0 $NCHIP_Y + for iyc 0 $NCELL_Y + for ix 0 $NCHIP_X + for ixc 0 $NCELL_X + concat cell_off[$ixc][$iyc] image_offset + end + end + end + end + + # generate the full set of corrections for a single filter + delete -q offset + for ns 0 $NSEASON + concat image_offset offset + end + + # the output fits table uses the vector names for the field names, + # so set them here to the desired names + set mjd_obs = mjd_uc + set zp = zpt_uc + set resid = 0.02*rnd(zpt_uc) - 0.01 + + delete -q flatcorr + dimenup offset flatcorr {$NSEASON*$NCHIP_X*$NCHIP_Y*$NCELL_X*$NCELL_Y} 1 + + wd phu testzpt.fits + + # filter 1 + write testzpt.fits -fits ZPTS_1 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w g + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 2 (?) + write testzpt.fits -fits ZPTS_2 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w r + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 3 + write testzpt.fits -fits ZPTS_3 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w i + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits +end + Index: /trunk/Ohana/src/addstar/test/relphot.variable.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/relphot.variable.dvo (revision 33653) +++ /trunk/Ohana/src/addstar/test/relphot.variable.dvo (revision 33653) @@ -0,0 +1,633 @@ +# -*-sh-*- + +# test suite for the parallel analysis mode for relphot & setphot +# in this test suite, we add in sources with outlier measurements. these are supposed to trigger the various rejection flags in relphot. does it work? + +# generate a realistic set of data: cover + +# measurement rejections that we perform: +# measurement is an outlier for the star +# measurement is flagged as bad by the photometry +# measurement is flagged as poor by the photometry +# Area & Time (ignore for now -- probably not used by GPC1) + +# also: +# bad images: chisq to high, zpt scatter too high +# bad objects : chisq to high, too few measurements? + +input tap.dvo + +# set globals +if (not($?PLOT)) set PLOT = 0 + +# set various global variables +macro init + echo "starting init..." + $RA_CENTER = 10.0 + $DEC_CENTER = 20.0 + $PLATE_SCALE = 0.25 + $NSEASON = 3 + $NFILTER = 3 + $NCHIP_X = 2 + $NCHIP_Y = 2 + $NCELL_X = 2 + $NCELL_Y = 2 + $CHIP_DX = 1000 + $CHIP_DY = 1000 + $DEFECT_FRAC = 0.03 + $OFFSET_FRAC_UC = 0.0 + $OFFSET_FRAC_NC = 0.5 + $LONGITUDE = 10.4170608521 + $LATITUDE = 20.7070999146 + if (not($?VERBOSE)) set VERBOSE = 0 + + # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters + $zpt_nominal:g = 24.58 + $zpt_nominal:r = 24.80 + $zpt_nominal:i = 24.74 + $zpt_nominal:z = 24.26 + $zpt_nominal:y = 23.41 + + # klam is negative, so klam*(airmass - 1) increase the zero point + $klam_nominal:g = -0.15 + $klam_nominal:r = -0.10 + $klam_nominal:i = -0.04 + $klam_nominal:z = -0.03 + $klam_nominal:y = -0.03 + + # we have two sets of images: ubercaled and not-ubercaled + + ##### parameters for UBERCAL'ed images + + $filt_uc:n = 3 + create filtN_uc 0 $filt_uc:n + $filt_uc:0 = g; filtN_uc[0] = 1 + $filt_uc:1 = g; filtN_uc[1] = 1 + $filt_uc:2 = g; filtN_uc[2] = 1 + + # airmass slopes for this sequnece + create klam_uc 0 $filt_uc:n + for i 0 klam_uc[] + klam_uc[$i] = $klam_nominal:$filt_uc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_uc = 15.0 + zero(klam_uc) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_uc = 25.0 + 2.5*log(exptime_uc) + klam_uc*(airmass_uc - 1.0) + ramp(klam_uc)*0.0050 - 0.0025 + + set mjd_uc = zero(zpt_uc) + mjd_uc[0] = 55000.01 + mjd_uc[1] = 55000.02 + mjd_uc[2] = 55000.03 + + create airmass_uc 0 klam_uc[] + set airmass_uc = 1.6 + zero(klam_uc) + + + ##### parameters for non-UBERCAL'ed images + + $filt_nc:n = 27 + create filtN_nc 0 $filt_nc:n + for i 0 $filt_nc + $filt_nc:$i = g + filtN_nc[$i] = 1 + end + + # airmass slopes for these sequnece + create klam_nc 0 $filt_nc:n + for i 0 klam_nc[] + klam_nc[$i] = $klam_nominal:$filt_nc:$i + end + + # place this within a valid season (55000.0 - 55010.0 - 55020.0 - 55030.0) + set exptime_nc = 15.0 + zero(klam_nc) + + # ubercal zero points are defined as ZP_nominal + 2.5log(exptime) + K*(airmass - 1.0) + set zpt_nc = 25.0 + 2.5*log(exptime_nc) + klam_nc*(airmass_nc - 1.0) - ramp(klam_nc)*0.0050 + 0.0025 + + # with an ubercal tie, relphot should set zpts to the inserted values above + set mjd_nc = zero(zpt_nc) + mjd_nc[0] = 55000.11 + mjd_nc[1] = 55000.12 + mjd_nc[2] = 55000.13 + mjd_nc[3] = 55015.11 + mjd_nc[4] = 55015.12 + mjd_nc[5] = 55015.13 + mjd_nc[6] = 55025.11 + mjd_nc[7] = 55025.12 + mjd_nc[8] = 55025.13 + + + for i + set airmass_nc = 1.6 + zero(klam_nc) +end + +macro go + if ($0 != 3) + echo "go (fileroot) (catdir)" + break + end + + local i fileroot rootdir catdir + + $fileroot = $1 + $catdir = $2 + + dirname $1 -var rootdir + mkdir $rootdir + + exec rm -rf $catdir + + init + mkzptfile + mkstars 10000 + + tapPLAN {(((mjd_uc[] + mjd_nc[])*4 + (mjd_uc[] + mjd_nc[])) + (mjd_uc[] + mjd_nc[]))} + + # generate the basic images and check they were correctly ingested by dvo + mksequence.drift $fileroot $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i raw + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + end + + # run setphot on the serial db and check that the images now match the expected values + exec setphot -reset -update -ubercal testzpt.fits -D CATDIR $catdir + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i setphot_uc + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i setphot_nc + end + + # run relphot on the serial db and check that the images now match the expected values + # if ($VERBOSE) echo relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + if ($VERBOSE) echo relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0 + exec relphot g -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update -reset -D STAR_SCATTER 0.05>& log.relphot.s0 + for i 0 mjd_uc[] + ckexposure $catdir mjd_uc[$i] zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] $filt_uc:$i relphot + end + for i 0 mjd_nc[] + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot + end + + tapDONE +end + +macro ckexposure + if ($0 != 8) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + echo " mode == raw or corr" + break + end + + local CATDIR MJD_IMAGE ZPT_REAL ZPT_NOMINAL EXPTIME FILTER AIRMASS + + $CATDIR = $1 + $MJD_IMAGE = $2 + $ZPT_REAL = $3 + $EXPTIME = $4 + $AIRMASS = $5 + $FILTER = $6 + $MODE = $7 + + # XXX need a function to extract the nominal zpt for a given filter / photcode from the db + $ZPT_NOMINAL = $zpt_nominal:$FILTER + $KLAM_NOMINAL = $klam_nominal:$FILTER + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + catdir $CATDIR + + # assume we still have stars_ra, stars_dec, stars_mag in hand + if ($PLOT) + dev -n 0 + region $RA_CENTER $DEC_CENTER 0.2 + images + pmeasure -all -m 15 20 + end + + skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} + mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + if (not(ra[])) + echo "no matching data for MJD = $MJD_IMAGE" + if ("$MODE" == "raw") + tapSKIP 4 + else + tapSKIP 1 + end + return + end + match2d -closest ra dec stars_ra stars_dec 0.001 -index1 index1 -index2 index2 + + reindex stars_ra_m = stars_ra using index1 + reindex stars_dec_m = stars_dec using index1 + reindex stars_mag_m = stars_mag using index1 + if (ra[] != stars_ra_m[]) + echo "failed to match stars in catdir to stars in memory" + break + end + set dr = 3600*(ra - stars_ra_m) + set dd = 3600*(dec - stars_dec_m) + set dm = mag - stars_mag_m + + if ($PLOT) + dev -n 0 + cplot -pt 7 -c red ra dec + dev -n 1 + lim mag dm; clear; box; plot mag dm + end + + set cell_xbin = int(xccd / ($CHIP_DX / $NCELL_X)) + set cell_ybin = int(yccd / ($CHIP_DY / $NCELL_Y)) + subset dm00 = dm if (cell_xbin == 0) && (cell_ybin == 0) + subset dm01 = dm if (cell_xbin == 0) && (cell_ybin == 1) + subset dm10 = dm if (cell_xbin == 1) && (cell_ybin == 0) + subset dm11 = dm if (cell_xbin == 1) && (cell_ybin == 1) + + # uncorrected values behave like this: + # mag_DVO = m_inst + zpt_nominal + 2.5*log(exptime) + K*(airmass - 1.0) + # mag_real = m_inst + zpt_real + cell_offset + # dm = mag_DVO - mag_real = zpt_nominal - zpt_real - cell_offset + # - zpt_nominal + zpt_real + cell_offset ~ 0.0 + # zpt_real (in this case) = 25.0 + 2.5*log(exptime) + K*(airmass - 1.0) + # (actually, it is the value in the vector 'zpt' for this entry + # zpt_nominal = 24.58 + + # setphot-corrected values behave like this: + # mag_DVO = mag_real (because m_inst has cell_offset applied) + # dm = mag_DVO - mag_real ~ 0.0 + + $ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0) + if ($VERBOSE) echo "$ZPT_REAL_NORM = $ZPT_REAL - 2.5*log($EXPTIME) - $KLAM_NOMINAL*($AIRMASS - 1.0)" + + if ("$MODE" == "raw") + for ix 0 $NCELL_X + for iy 0 $NCELL_Y + # vstat -q dm$ix\$iy + # echo cell_off[$ix][$iy] {$MEDIAN - $ZPT_NOMINAL + $ZPT_REAL + cell_off[$ix][$iy]} $MEDIAN $MEAN $SIGMA + + set dm_adjust = dm$ix\$iy - $ZPT_NOMINAL + $ZPT_REAL_NORM + vstat -q dm_adjust + tapOK {abs($MEAN) < 0.005} "addstar raw ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + end + end + return + end + + if ("$MODE" == "setphot_uc") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "setphot_uc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "setphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "setphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot_nc") + vstat -q dm + tapOK {abs($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM) < 0.005} "relphot_nc ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + + if ("$MODE" == "relphot") + vstat -q dm + tapOK {abs($MEAN) < 0.005} "relphot ZP $MJD_IMAGE $FILTER : $MEAN $ZPT_NOMINAL $ZPT_REAL_NORM {($MEAN - $ZPT_NOMINAL + $ZPT_REAL_NORM)}" + return + end + +end + +macro mksequence.drift + if ($0 != 3) + echo "mksequence.drift (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + mkexposure $1.uc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_UC/3600.0} zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER {$DEC_CENTER + $i*$CHIP_DY*$PLATE_SCALE*$NCHIP_Y*$OFFSET_FRAC_NC/3600.0} zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mksequence.1pt + if ($0 != 3) + echo "mksequence (fileroot) (catdir)" + break + end + + echo "starting mksequence..." + local i T1 T2 + + ctimes -abs now -var T1 + # we have a second set of mjd and zero point values for non-ubercal'ed images + echo -no-return "making mjd_uc[] UC images: " + for i 0 mjd_uc[] + echo -no-return "$i.. " + mkexposure $1.uc.$i $RA_CENTER $DEC_CENTER zpt_uc[$i] exptime_uc[$i] airmass_uc[$i] mjd_uc[$i] $filt_uc:$i $2 + end + echo "" + echo -no-return "making mjd_nc[] NC images: " + for i 0 mjd_nc[] + echo -no-return "$i.. " + mkexposure $1.nc.$i $RA_CENTER $DEC_CENTER zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] mjd_nc[$i] $filt_nc:$i $2 + end + ctimes -abs now -var T2 + echo "done ({86400*($T2 - $T1)} sec)" +end + +macro mkexposure + if ($0 != 10) + echo "mkexposure (fileroot) (ra) (dec) (zpt) (exptime) (airmass) (mjd) (filter) (catdir)" + break + end + + local ix iy date time datetime ra dec ROOT RAo DECo ZPT MJD FILTER EXPTIME AIRMASS CATDIR + + $ROOT = $1 + $RAo = $2 + $DECo = $3 + $ZPT = $4 + $EXPTIME = $5 + $AIRMASS = $6 + $MJD = $7 + $FILTER = $8 + $CATDIR = $9 + + $TIMEFORMAT = mjd + $TIMEREF = 0.0 + + for ix 0 $NCHIP_X + for iy 0 $NCHIP_Y + $dx = $CHIP_DX * $PLATE_SCALE * ($ix - 0.5*$NCHIP_X + 0.5) + $dy = $CHIP_DY * $PLATE_SCALE * ($iy - 0.5*$NCHIP_Y + 0.5) + $ra = $RAo - $dx / 3600.0 / dcos($DECo) + $dec = $DECo - $dy / 3600.0 + # echo $ra $dec $dx $dy + if ($VERBOSE) echo mkinput test.in.txt $ra $dec $ZPT + mkinput test.in.txt $ra $dec $ZPT + # costs 7/36 sec / chip + + # ra,dec is the center of this chip + local options + $options = -mjd $MJD + $options = $options -radec $ra $dec + $options = $options -type PS1_V2 + $options = $options -coords + # use photcodes XY01, XY11, XY02, XY12 to avoid the corners + $options = $options -photcode GPC1.$FILTER.XY$ix\{$iy+1} + $options = $options -no-noise + $options = $options -flags + $options = $options -bad-psfqf-frac $DEFECT_FRAC + $options = $options -size $CHIP_DX $CHIP_DY + $options = $options -crpix {0.5*$CHIP_DX} {0.5*$CHIP_DY} + $options = $options -airmass $AIRMASS + $options = $options -exptime $EXPTIME + if ($VERBOSE) echo mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + exec mkcmf test.in.txt $ROOT.$ix.$iy.cmf $options + # costs 8/36 sec / chip + + # the fake images have inconsistent ra,dec and airmass,sidtime values + if ($VERBOSE) echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& /dev/null + # costs 10/36 sec / chip + end + end +end + +# make a simple input file for mkcmf +macro mkinput + if ($0 != 5) + echo "mkinput (filename) (ra) (dec) (zpt)" + break + end + + # chip coordinate of the stars + # ra,dec is at the chip center, which corresponds to pixel (CHIP_DX,CHIP_DY)/2 + set stars_X = (stars_ra - $2) * 3600.0 * dcos($3) / $PLATE_SCALE + $CHIP_DX*0.5 + set stars_Y = (stars_dec - $3) * 3600.0 / $PLATE_SCALE + $CHIP_DY*0.5 + set stars_M = (stars_mag - $4) + + set valid = (stars_X > 0) && (stars_X < $CHIP_DX) && (stars_Y > 0) && (stars_Y < $CHIP_DX) + subset stars_x = stars_X if valid + subset stars_y = stars_Y if valid + subset stars_m = stars_M if valid + subset stars_r = stars_ra if valid + subset stars_d = stars_dec if valid + subset stars_mr = stars_mag if valid + + set cell_x = int(stars_x / ($CHIP_DX / $NCELL_X)) + set cell_y = int(stars_y / ($CHIP_DX / $NCELL_X)) + + if (1) + # FLAT-FIELD CORRECTION + # eddie says: M_real = M_inst + zpt + ff_offset + # thus, stars_mi = stars_mag - zpt - ff_offset = stars_m - offset + set stars_mi = zero(stars_m) + for i 0 stars_mi[] + $nx = cell_x[$i] + $ny = cell_y[$i] + # $offset = cell_off[$nx][$ny] + $offset = 0.0 + stars_mi[$i] = stars_m[$i] - $offset + end + else + set stars_mi = stars_m + end + + # create bad detections / detections to ignore or avoid, save to a file? + + # create 3 types of things to ignore / avoid + # * badPhot (deviate the photometry) + # * badFlag (set a bad flag bit) + # * poorFlag (set a bad flag bit) + set badPhot = rnd(stars_r) < $DEFECT_FRAC + set badFlag = rnd(stars_r) < $DEFECT_FRAC + set poorFlag = rnd(stars_r) < $DEFECT_FRAC + + # create deviant photometry (1%) + # I'm setting the flag value 0x02 (EXTMODEL) if the object is deviated. This is a simple way to pass the info + set stars_mi = stars_mi + 0.2*badPhot + set stars_fl = zero(stars_mi) + 0x00000080*badFlag + 0x00000100*poorFlag + 0x00000002*badPhot + + # bad flags that we sould respect: + # PM_SOURCE_MODE_FAIL = 0x00000008, ///< Fit (non-linear) failed (non-converge, off-edge, run to zero) + # PM_SOURCE_MODE_SATSTAR = 0x00000080, ///< Source model peak is above saturation + # PM_SOURCE_MODE_BADPSF = 0x00000400, ///< Failed to get good estimate of object's PSF + # PM_SOURCE_MODE_DEFECT = 0x00000800, ///< Source is thought to be a defect + # PM_SOURCE_MODE_SATURATED = 0x00001000, ///< Source is thought to be saturated pixels (bleed trail) + # PM_SOURCE_MODE_CR_LIMIT = 0x00002000, ///< Source has crNsigma above limit + # PM_SOURCE_MODE_MOMENTS_FAILURE = 0x00008000, ///< could not measure the moments + # PM_SOURCE_MODE_SKY_FAILURE = 0x00010000, ///< could not measure the local sky + # PM_SOURCE_MODE_SKYVAR_FAILURE = 0x00020000, ///< could not measure the local sky variance + # PM_SOURCE_MODE_SIZE_SKIPPED = 0x10000000, ///< size could not be determined + + # poor flags that we should respect: + # PM_SOURCE_MODE_POOR = 0x00000010, ///< Fit succeeds, but low-SN, high-Chisq, or large (for PSF -- drop?) + # PM_SOURCE_MODE_PAIR = 0x00000020, ///< Source fitted with a double psf + # PM_SOURCE_MODE_BLEND = 0x00000100, ///< Source is a blend with other sources + # PM_SOURCE_MODE_BELOW_MOMENTS_SN = 0x00040000, ///< moments not measured due to low S/N + # PM_SOURCE_MODE_BLEND_FIT = 0x00400000, ///< source was fitted as a blend + # PM_SOURCE_MODE_ON_SPIKE = 0x20000000, ///< peak lands on diffraction spike + # PM_SOURCE_MODE_ON_GHOST = 0x40000000, ///< peak lands on ghost or glint + # PM_SOURCE_MODE_OFF_CHIP = 0x80000000, ///< peak lands off edge of chip + + # delete the output file before writing + exec rm -f $1 + write -f "%10.6f %10.6f %7.3f %6.1f %6.1f %7.3f 0x%x" $1 stars_r stars_d stars_mr stars_x stars_y stars_mi stars_fl +end + +# generate a set of stars with raw RA, DEC, MAG values to use for all of the output cmf files +macro mkstars + if ($0 != 2) + echo "USAGE: mkstars (Nstars)" + break + end + + echo "starting mkstars..." + local RA_RANGE DEC_RANGE + + # the images are oriented along N-S, E-W lines + + # size of region of interest in linear arcseconds + $RA_RANGE = 1.25 * $CHIP_DX * $NCHIP_X * $PLATE_SCALE + $DEC_RANGE = 1.25 * $CHIP_DY * $NCHIP_Y * $PLATE_SCALE * (1 + tmpseq[]*$OFFSET_FRAC_NC) + + create tmp 0 $1 + set stars_ra = $RA_CENTER + $RA_RANGE * (rnd(tmp) - 0.5) / 3600 / dcos ($DEC_CENTER) + set stars_dec = $DEC_CENTER + $DEC_RANGE * (rnd(tmp) - 0.5) / 3600 + set stars_mag = 15.0 + 4.0 * rnd(tmp) +end + +# copy the catdir to catdir.p0 and parallelize to local directories catdir.d1 - catdir.d3 +macro mkparallel + if ($0 != 2) + echo "USAGE: mkparallel (catdir)" + break + end + + local hostname catdir + + $hostname = `hostname -s` + + $tmp1 = `dirname $1` + $tmp2 = `basename $1` + $catdir = $tmp1/$tmp2 + + exec rsync -auv $catdir/ $catdir.p0/ >& log.rsync + mkdir $catdir.d1 + mkdir $catdir.d2 + mkdir $catdir.d3 + + exec rm -f $catdir.p0/HostTable.dat + exec echo "1 $hostname $catdir.d1" > $catdir.p0/HostTable.dat + exec echo "2 $hostname $catdir.d2" >> $catdir.p0/HostTable.dat + exec echo "3 $hostname $catdir.d3" >> $catdir.p0/HostTable.dat + + exec dvodist -out $catdir.p0 >& log.dvodist.out +end + +# create a fake ubercal-style zpt/flatcorr table with metadata +macro mkzptfile + + mcreate phu 0 0 + keyword phu NSEASON -wd $NSEASON + keyword phu NFILTER -wd $NFILTER + keyword phu NCHIP_X -wd $NCHIP_X + keyword phu NCHIP_Y -wd $NCHIP_Y + keyword phu NCELL_X -wd $NCELL_X + keyword phu NCELL_Y -wd $NCELL_Y + keyword phu CHIP_DX -wd $CHIP_DX + keyword phu CHIP_DY -wd $CHIP_DY + + # define the season boundaries + keyword phu TS0_0000 -wf 55000.0 + keyword phu TS1_0000 -wf 55010.0 + keyword phu TS0_0001 -wf 55010.0 + keyword phu TS1_0001 -wf 55020.0 + keyword phu TS0_0002 -wf 55020.0 + keyword phu TS1_0002 -wf 55030.0 + + # create a single flat-field correction map for a single chip + mcreate cell_off 2 2 + cell_off[0][0] = 0.01 + cell_off[1][0] = -0.01 + cell_off[0][1] = 0.02 + cell_off[1][1] = -0.02 + + # in the correction table, we have an image of the full exposure + # (NCHIP_X x NCHIP_Y) for each season, unwrapped as a single linear + # vector (with the last two missing -- how is that hard wired?) + + local ix iy ixc iyc ns + + # unroll the chip & cell level corrections into a single vector for a single image + delete -q image_offset + for iy 0 $NCHIP_Y + for iyc 0 $NCELL_Y + for ix 0 $NCHIP_X + for ixc 0 $NCELL_X + concat cell_off[$ixc][$iyc] image_offset + end + end + end + end + + # generate the full set of corrections for a single filter + delete -q offset + for ns 0 $NSEASON + concat image_offset offset + end + + # the output fits table uses the vector names for the field names, + # so set them here to the desired names + set mjd_obs = mjd_uc + set zp = zpt_uc + set resid = 0.02*rnd(zpt_uc) - 0.01 + + delete -q flatcorr + dimenup offset flatcorr {$NSEASON*$NCHIP_X*$NCHIP_Y*$NCELL_X*$NCELL_Y} 1 + + wd phu testzpt.fits + + # filter 1 + write testzpt.fits -fits ZPTS_1 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w g + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 2 (?) + write testzpt.fits -fits ZPTS_2 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w r + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits + + # filter 3 + write testzpt.fits -fits ZPTS_3 -f DDE -append mjd_obs zp resid + keyword flatcorr FILTER -w i + keyword flatcorr EXTNAME -w FLATCORR + wd -extend flatcorr testzpt.fits +end + Index: /trunk/Ohana/src/addstar/test/run.relastro.sh =================================================================== --- /trunk/Ohana/src/addstar/test/run.relastro.sh (revision 33653) +++ /trunk/Ohana/src/addstar/test/run.relastro.sh (revision 33653) @@ -0,0 +1,55 @@ +#!/bin/csh -f + +set region = "-region 248 251 -2 1.2" +set cat2M = catdir.2mass +set catPS = catdir.PS1 +set raw2M = /data/ipp022.0/alala.0/ipp/ippRefs/catdir.2mass + +# create the database +foreach file (smallset/*.smf) + addstar -D PHOTCODE_FILE dvo.photcodes.3pi -D CAMERA gpc1 -D CATDIR $catPS $file -update +end +addstar -D CAMERA gpc1 -D CATDIR $catPS -resort + +# NOTE: if the photcode table does not match the desire, do this: +# dvosecfilt 8 catdir.small.merge.v3 +# photcode-table -import dvo.photcodes.3pi -D CATDIR catdir.small.merge.v3 + +# make a local 2MASS database chunk +dvosplit $raw2M 4 -outdir $cat2M -set-mode split -set-format PS1_REF $region +addstar -resort -D CATDIR $cat2M $region +# cp $raw2M/Images.dat catdir.2mass.d4/ + +# dvomerge A and B to C is broken, so copy A (serial and parallel) +rsync -auv $catPS/ $catPS.ser/ +rsync -auv $catPS/ $catPS.par/ + +# merge 2MASS and PS1 (serial and parallel) +dvomerge $cat2M into $catPS.ser $region +dvomerge $cat2M into $catPS.par $region + +# ensure g,r,i,z,y,J,H,K ave values are set +relphot -averages -update -D CATDIR $catPS.ser +relphot -averages -update -D CATDIR $catPS.par + +# now parallelize the db +sed s/CATDIR/$catPS.par/ < HostTable.dat > $catPS.par/HostTable.dat +dvodist -out $catPS.par + +relastro -parallel -D CATDIR $catPS.par -high-speed r J 20.0 catdir.highspeed.par $region -D WHERE_A "(r:nphot > 3) && (r:err > 0.05)" -D WHERE_B "(J:err > 0.05)" +relastro -D CATDIR $catPS.ser -high-speed r J 20.0 catdir.highspeed.ser $region -D WHERE_A "(r:nphot > 3) && (r:err > 0.05)" -D WHERE_B "(J:err > 0.05)" + +relphot -reset -averages -update -D CATDIR catdir.highspeed.ser +relphot -reset -averages -update -D CATDIR catdir.highspeed.par + +# test relastro runs +if (0) then + relastro -D CATDIR $catPS.ser $region -update-chips -update + relastro -D CATDIR $catPS.ser $region -update-objects -update + relastro -D CATDIR $catPS.ser $region -update-objects -update -pm + + relastro -parallel -D CATDIR $catPS.par $region -update-chips -update + relastro -parallel -D CATDIR $catPS.par $region -update-objects -update + relastro -parallel -D CATDIR $catPS.par $region -update-objects -update -pm +endif + Index: /trunk/Ohana/src/addstar/test/simple.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/simple.dvo (revision 33652) +++ /trunk/Ohana/src/addstar/test/simple.dvo (revision 33653) @@ -8,4 +8,5 @@ test.fields PS1_V1 PS1_V1 test.fields PS1_V2 PS1_V1 + test.fields PS1_V3 PS1_V1 test.fields PS1_DEV_0 PS1_V2 @@ -13,4 +14,11 @@ test.fields PS1_V1 PS1_V2 test.fields PS1_V2 PS1_V2 + test.fields PS1_V3 PS1_V2 + + test.fields PS1_DEV_0 PS1_V3 + test.fields PS1_DEV_1 PS1_V3 + test.fields PS1_V1 PS1_V3 + test.fields PS1_V2 PS1_V3 + test.fields PS1_V3 PS1_V3 end @@ -31,5 +39,5 @@ mkinput exec mkcmf test.in.txt test.cmf -date 2008/1/1 -time 01:00:00 -radec $RA $DEC -type $1 - exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 + exec addstar -D CATDIR catdir.test -D CAMERA simtest test.cmf -D CATFORMAT $2 -quick-airmass list testfields -copy testfields_$1 @@ -220,2 +228,38 @@ N_FRAMES : SKIP # no accessor end + +# this list is good for PS1_V3 +list testfields_PS1_V3 + IPP_IDET : detid + X_PSF : xccd + Y_PSF : yccd + X_PSF_SIG : xccd:err # FAIL + Y_PSF_SIG : yccd:err # FAIL + RA_PSF : SKIP # astrometry is not calibrated in the cmf + DEC_PSF : SKIP # astrometry is not calibrated in the cmf + POSANGLE : SKIP # astrometry is not calibrated in the cmf + PLTSCALE : SKIP # astrometry is not calibrated in the cmf + PSF_INST_MAG : mag:inst + PSF_INST_MAG_SIG : mag:err + AP_MAG_STANDARD : mag:ap # FAIL + AP_MAG_RADIUS : SKIP # no accessor + PEAK_FLUX_AS_MAG : SKIP # no accessor + CAL_PSF_MAG : SKIP # photometry is not calibrated in the cmf + CAL_PSF_MAG_SIG : SKIP # photometry is not calibrated in the cmf + SKY : sky + SKY_SIG : sky_err + PSF_CHISQ : psf_chisq + CR_NSIGMA : cr_nsigma + EXT_NSIGMA : ext_nsigma + PSF_MAJOR : FWHM_MAJ + PSF_MINOR : FWHM_MIN + PSF_THETA : THETA # FAIL + PSF_QF : PSF_QF + PSF_NDOF : SKIP # no accessor + PSF_NPIX : SKIP # no accessor + MOMENTS_XX : SKIP # no accessor + MOMENTS_XY : SKIP # no accessor + MOMENTS_YY : SKIP # no accessor + FLAGS : phot_flags + N_FRAMES : SKIP # no accessor +end Index: /trunk/Ohana/src/addstar/test/tap.dvo =================================================================== --- /trunk/Ohana/src/addstar/test/tap.dvo (revision 33652) +++ /trunk/Ohana/src/addstar/test/tap.dvo (revision 33653) @@ -1,3 +1,5 @@ # -*-sh-*- + +if (not($?TAP_BREAK)) set TAP_BREAK = 0 macro tapOK @@ -29,6 +31,17 @@ $TAP_NTEST = $1 $TAP_NFAIL = 0 + $TAP_NSKIP = 0 $TAP_NDONE = 0 if (not($?TAP_BREAK)) set TAP_BREAK = 0 +end + +macro tapSKIP + if ($0 != 2) + echo "USAGE: tapSKIP (nskip)" + break + end + + $TAP_NSKIP ++ + $TAP_NDONE ++ end @@ -47,5 +60,9 @@ end - if ($TAP_NFAIL == 0) + if ($TAP_NSKIP) + echo "skipped $TAP_NSKIP of $TAP_NDONE" + end + + if ($TAP_NFAIL + $TAP_NSKIP == 0) echo "passed $TAP_NDONE tests" end