Index: trunk/Ohana/src/addstar/src/FilterStars.c =================================================================== --- trunk/Ohana/src/addstar/src/FilterStars.c (revision 28968) +++ trunk/Ohana/src/addstar/src/FilterStars.c (revision 29001) @@ -1,5 +1,3 @@ # include "addstar.h" -#include "pslib/psMetadata.h" -#include "psmodules/pmSourceMasks.h" # define EXTERNAL_ID TRUE Index: trunk/Ohana/src/addstar/src/LoadData.c =================================================================== --- trunk/Ohana/src/addstar/src/LoadData.c (revision 28968) +++ trunk/Ohana/src/addstar/src/LoadData.c (revision 29001) @@ -80,4 +80,8 @@ inStars = ReadStarsFITS (f, headers[Nhead], headers[Ndata], &images[0][Nvalid].nstar); + if (!inStars) { + // XXX need to free the data here or in ReadStarsFITS on error + continue; + } inStars = FilterStars (inStars, &images[0][Nvalid], Nvalid, options); *stars = MergeStars (*stars, Nstars, inStars, images[0][Nvalid].nstar); Index: trunk/Ohana/src/addstar/src/ReadStarsFITS.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 28968) +++ trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 29001) @@ -49,5 +49,5 @@ if (stars == NULL) { fprintf (stderr, "ERROR: invalid table type %s\n", type); - exit (1); + return (NULL); } // Nstars is not necessarily == *nstars (The former is the number of detections, the @@ -69,4 +69,8 @@ swapped = FALSE; smpdata = gfits_table_get_SMPData (table, &Nstars, &swapped); + if (!smpdata) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); @@ -108,4 +112,8 @@ ps1data = gfits_table_get_PS1_DEV_0 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); @@ -153,4 +161,8 @@ ps1data = gfits_table_get_PS1_DEV_1 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); @@ -212,4 +224,8 @@ ps1data = gfits_table_get_CMF_PS1_V1 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); @@ -284,4 +300,8 @@ ps1data = gfits_table_get_CMF_PS1_V1_Alt (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); @@ -348,4 +368,8 @@ ps1data = gfits_table_get_CMF_PS1_V2 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } ZeroPt = GetZeroPoint(); Index: trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 28968) +++ trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 29001) @@ -218,5 +218,5 @@ } - N = N1[i]; + N = i; if (stars[N][0].found > -1) continue; if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue; @@ -286,5 +286,5 @@ for (j = 0; j < NREFSTAR_GROUP; j++) { - N = N1[i + j]; + N = i + j; catalog[0].measure[Nmeas] = stars[N][0].measure; Index: trunk/Ohana/src/addstar/test/dvomerge.dvo =================================================================== --- trunk/Ohana/src/addstar/test/dvomerge.dvo (revision 28968) +++ trunk/Ohana/src/addstar/test/dvomerge.dvo (revision 29001) @@ -95,5 +95,5 @@ subset T = t if (id == imageID[$i]) set dT = T - time[$i] - vstat -q dT + vstat dT tapOK {abs($MEAN) < 0.00001} "time for measure ID $i (MEAN)" tapOK {abs($SIGMA) < 0.00001} "time for measure ID $i (SIGMA)" Index: trunk/Ohana/src/delstar/src/delete_imagename.c =================================================================== --- trunk/Ohana/src/delstar/src/delete_imagename.c (revision 28968) +++ trunk/Ohana/src/delstar/src/delete_imagename.c (revision 29001) @@ -19,4 +19,8 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } /* find image in db by name */ Index: trunk/Ohana/src/delstar/src/delete_times.c =================================================================== --- trunk/Ohana/src/delstar/src/delete_times.c (revision 28968) +++ trunk/Ohana/src/delstar/src/delete_times.c (revision 29001) @@ -26,4 +26,8 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } /* find images for time range, delete each image */ Index: trunk/Ohana/src/delstar/src/find_image_db.c =================================================================== --- trunk/Ohana/src/delstar/src/find_image_db.c (revision 28968) +++ trunk/Ohana/src/delstar/src/find_image_db.c (revision 29001) @@ -17,4 +17,8 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } Nlist = 0; Index: trunk/Ohana/src/dvomerge/Makefile =================================================================== --- trunk/Ohana/src/dvomerge/Makefile (revision 28968) +++ trunk/Ohana/src/dvomerge/Makefile (revision 29001) @@ -20,5 +20,5 @@ dvomerge : $(BIN)/dvoconvert.$(ARCH) -all: dvomerge dvoconvert +all: dvomerge dvoconvert dvosecfilt DVOMERGE = \ @@ -58,5 +58,17 @@ $(BIN)/dvoconvert.$(ARCH) : $(DVOCONVERT) -INSTALL = dvomerge dvoconvert +DVOSECFILT = \ +$(SRC)/dvosecfilt.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/help.$(ARCH).o \ +$(SRC)/args.$(ARCH).o + +$(DVOSECFILT) : $(INC)/dvomerge.h + +$(BIN)/dvosecfilt.$(ARCH) : $(DVOSECFILT) + +INSTALL = dvomerge dvoconvert dvosecfilt # dependancy rules for binary code ######################### Index: trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 28968) +++ trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 29001) @@ -55,5 +55,10 @@ int dvoconvert_args PROTO((int *argc, char **argv)); +void dvosecfilt_usage PROTO((void)); +void dvosecfilt_help PROTO((int argc, char **argv)); +int dvosecfilt_args PROTO((int *argc, char **argv)); + int SkyTablePopulatedRange PROTO((off_t *ns, off_t *ne, SkyTable *sky, off_t Nstart)); +int SkyListPopulatedRange PROTO((off_t *ns, off_t *ne, SkyList *sky, off_t Nstart)); SkyList *SkyTablePopulatedList PROTO((SkyTable *sky)); Index: trunk/Ohana/src/dvomerge/src/SkyTablePopulatedRange.c =================================================================== --- trunk/Ohana/src/dvomerge/src/SkyTablePopulatedRange.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/SkyTablePopulatedRange.c (revision 29001) @@ -12,4 +12,27 @@ Ns = sky[0].regions[Nstart].childS; Ne = sky[0].regions[Nstart].childE; + if (Ns == 0) { + fprintf (stderr, "no populated tables at an appropriate depth\n"); + exit (1); + } + Nstart = Ns; + } + + *ns = Ns; + *ne = Ne; + return (TRUE); +} + +int SkyListPopulatedRange (off_t *ns, off_t *ne, SkyList *sky, off_t Nstart) { + + off_t Ns, Ne; + + // given the starting sky region, find the populated range at or below this region + + Ns = Nstart; + Ne = Nstart + 1; + while (!sky[0].regions[Nstart][0].table) { + Ns = sky[0].regions[Nstart][0].childS; + Ne = sky[0].regions[Nstart][0].childE; if (Ns == 0) { fprintf (stderr, "no populated tables at an appropriate depth\n"); Index: trunk/Ohana/src/dvomerge/src/args.c =================================================================== --- trunk/Ohana/src/dvomerge/src/args.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/args.c (revision 29001) @@ -15,4 +15,22 @@ remove_argument (N, argc, argv); ALTERNATE_PHOTCODE_FILE = strdup(argv[N]); + remove_argument (N, argc, argv); + } + + /*** provide additional data ***/ + /* restrict to a portion of the sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (*argc, argv, "-region"))) { + remove_argument (N, argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmax = atof (argv[N]); remove_argument (N, argc, argv); } @@ -37,2 +55,18 @@ return TRUE; } + +/*** check for command line options ***/ +int dvosecfilt_args (int *argc, char **argv) { + + int N; + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + if (*argc != 3) dvosecfilt_usage(); + return TRUE; +} Index: trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 29001) @@ -13,4 +13,9 @@ images = gfits_table_get_Image (&in[0].ftable, &Nimages, &in[0].swapped); + if (!images) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + ALLOCATE (IDmap->old, off_t, Nimages); ALLOCATE (IDmap->new, off_t, Nimages); @@ -59,4 +64,6 @@ off_t Nlo, Nhi, N; + + if (!IDmap->Nmap) return 0; // find the a close entry below desired ID Index: trunk/Ohana/src/dvomerge/src/dvoconvert.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoconvert.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/dvoconvert.c (revision 29001) @@ -148,4 +148,8 @@ // convert the raw image table to Image type (byteswap if needed) images = gfits_table_get_Image (&inDB.ftable, &Nimages, &inDB.swapped); + if (!images) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } // update additional metadata Index: trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 29001) @@ -6,5 +6,5 @@ off_t i, j, Ns, Ne; SkyTable *outsky, *insky; - SkyList *inlist; + SkyList *outlist, *inlist; Catalog incatalog, outcatalog; char filename[256], *input, *output; @@ -70,4 +70,7 @@ SkyTableSetFilenames (insky, input, "cpt"); + // XXX apply this...generate the subset matching the user-selected region + inlist = SkyListByPatch (insky, -1, &UserPatch); + // generate an output table populated at the desired depth outsky = SkyTableLoadOptimal (output, NULL, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); @@ -79,15 +82,21 @@ // loop over the populatable output tables; check for data in input in the corresponding regions - SkyTablePopulatedRange (&Ns, &Ne, insky, 0); - depth = insky[0].regions[Ns].depth; + SkyListPopulatedRange (&Ns, &Ne, inlist, 0); + depth = inlist[0].regions[Ns][0].depth; - // loop over the populatable output regions - for (i = 0; i < outsky[0].Nregions; i++) { - if (!outsky[0].regions[i].table) continue; - if (VERBOSE) fprintf (stderr, "output: %s\n", outsky[0].regions[i].name); + // loop over the populated input regions + for (i = 0; i < inlist[0].Nregions; i++) { + if (!inlist[0].regions[i][0].table) continue; + if (VERBOSE) fprintf (stderr, "input: %s\n", inlist[0].regions[i][0].name); - SetPhotcodeTable(outputPhotcodes); + SetPhotcodeTable(inputPhotcodes); // load / create output catalog (if catalog does not exist, it will be created) - LoadCatalog (&outcatalog, &outsky[0].regions[i], outsky[0].filename[i], "w"); + LoadCatalog (&incatalog, &inlist[0].regions[i][0], inlist[0].filename[i], "r"); + // skip empty input catalogs + if (!incatalog.Naves_disk) { + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); + continue; + } // combine only tables at equal or larger depth @@ -97,31 +106,31 @@ // compare to a slightly reduced footprint float dPos = 2.0/3600.0; - inlist = SkyListByBounds (insky, depth, outsky[0].regions[i].Rmin + dPos, outsky[0].regions[i].Rmax - dPos, outsky[0].regions[i].Dmin + dPos, outsky[0].regions[i].Dmax - dPos); - for (j = 0; j < inlist[0].Nregions; j++) { - if (VERBOSE) fprintf (stderr, "input : %s\n", inlist[0].regions[j][0].name); + outlist = SkyListByBounds (outsky, depth, inlist[0].regions[i][0].Rmin + dPos, inlist[0].regions[i][0].Rmax - dPos, inlist[0].regions[i][0].Dmin + dPos, inlist[0].regions[i][0].Dmax - dPos); + for (j = 0; j < outlist[0].Nregions; j++) { + if (VERBOSE) fprintf (stderr, "output : %s\n", outlist[0].regions[j][0].name); // load input catalog - SetPhotcodeTable(inputPhotcodes); - LoadCatalog (&incatalog, inlist[0].regions[j], inlist[0].filename[j], "r"); + SetPhotcodeTable(outputPhotcodes); + LoadCatalog (&outcatalog, outlist[0].regions[j], outlist[0].filename[j], "w"); - // skip empty input catalogs - if (!incatalog.Naves_disk) { - dvo_catalog_unlock (&incatalog); - dvo_catalog_free (&incatalog); - continue; - } dvo_update_image_IDs (&IDmap, &incatalog); - merge_catalogs_old (&outsky[0].regions[i], &outcatalog, &incatalog, RADIUS, secfiltMap); - dvo_catalog_unlock (&incatalog); - dvo_catalog_free (&incatalog); + merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap); - fprintf (stderr, "merged %s into %s\n", outsky[0].regions[i].name, inlist[0].regions[j][0].name); + outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + // if we receive a signal which would cause us to exit, wait until the full catalog is written + SetProtect (TRUE); + dvo_catalog_save (&outcatalog, VERBOSE); + SetProtect (FALSE); + + dvo_catalog_unlock (&outcatalog); + dvo_catalog_free (&outcatalog); + + fprintf (stderr, "merged %s into %s\n", inlist[0].regions[i][0].name, outlist[0].regions[j][0].name); } - SkyListFree (inlist); + SkyListFree (outlist); - outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; - dvo_catalog_save (&outcatalog, VERBOSE); - dvo_catalog_unlock (&outcatalog); - dvo_catalog_free (&outcatalog); + dvo_catalog_unlock (&incatalog); + dvo_catalog_free (&incatalog); } @@ -153,5 +162,7 @@ // load the image table if (inDB.dbstate == LCK_EMPTY) { - Shutdown ("can't find input image catalog %s", inDB.filename); + // Shutdown ("can't find input image catalog %s", inDB.filename); + IDmap->Nmap = 0; + return TRUE; } if (!dvo_image_load (&inDB, VERBOSE, TRUE)) { Index: trunk/Ohana/src/dvomerge/src/dvosecfilt.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvosecfilt.c (revision 29001) +++ trunk/Ohana/src/dvomerge/src/dvosecfilt.c (revision 29001) @@ -0,0 +1,90 @@ +# include "dvomerge.h" + +// modify the number of average photcodes : change only secfilt tables in place +int main (int argc, char **argv) { + + char filename[256], *input; + off_t i, j, k, Nsecfilt, NsecInput, Nstart; + SkyTable *insky; + Catalog catalog; + SecFilt *insec, *outsec; + + SetSignals (); + dvosecfilt_help (argc, argv); + ConfigInit (&argc, argv); + dvosecfilt_args (&argc, argv); + + input = argv[1]; + Nsecfilt = atoi(argv[2]); + + // load the sky table for the existing database + insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (insky, input, "cpt"); + + // loop over all input catalogs + for (i = 0; i < insky[0].Nregions; i++) { + if (!insky[0].regions[i].table) continue; + if (VERBOSE) fprintf (stderr, "table: %s\n", insky[0].regions[i].name); + + // set the parameters which guide catalog open/load/create + catalog.filename = insky[0].filename[i]; + catalog.Nsecfilt = 0; + + // always load all of the data (if any exists) + catalog.catflags = LOAD_SECF; + 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 + + if (!dvo_catalog_open (&catalog, &insky[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", filename); + exit (2); + } + + if (catalog.Naves_disk == 0) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + if (Nsecfilt == catalog.Nsecfilt) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + NsecInput = catalog.Nsecfilt; + Nstart = MIN(Nsecfilt, NsecInput); + + insec = catalog.secfilt; + ALLOCATE (outsec, SecFilt, catalog.Naves_disk * Nsecfilt); + + for (k = 0; k < catalog.Naves_disk; k++) { + for (j = 0; (j < catalog.Nsecfilt) && (j < Nsecfilt); j++) { + outsec[k*Nsecfilt + j].M = insec[k*NsecInput + j].M; + outsec[k*Nsecfilt + j].dM = insec[k*NsecInput + j].dM; + outsec[k*Nsecfilt + j].Xm = insec[k*NsecInput + j].Xm; + outsec[k*Nsecfilt + j].flags = insec[k*NsecInput + j].flags; + outsec[k*Nsecfilt + j].Ncode = insec[k*NsecInput + j].Ncode; + outsec[k*Nsecfilt + j].Nused = insec[k*NsecInput + j].Nused; + outsec[k*Nsecfilt + j].M_20 = insec[k*NsecInput + j].M_20; + outsec[k*Nsecfilt + j].M_80 = insec[k*NsecInput + j].M_80; + } + for (j = Nstart; j < Nsecfilt; j++) { + outsec[k*Nsecfilt + j].M = NAN; + outsec[k*Nsecfilt + j].dM = NAN; + outsec[k*Nsecfilt + j].Xm = NAN_S_SHORT; + } + } + free (catalog.secfilt); + catalog.secfilt = outsec; + catalog.Nsecfilt = Nsecfilt; + catalog.Nsecf_mem = Nsecfilt * catalog.Naves_disk; + catalog.Nsecf_disk = Nsecfilt * catalog.Naves_disk; + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + + exit (0); +} Index: trunk/Ohana/src/dvomerge/src/dvoverify.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/dvoverify.c (revision 29001) @@ -154,4 +154,8 @@ // convert the raw image table to Image type (byteswap if needed) images = gfits_table_get_Image (&inDB.ftable, &Nimages, &inDB.swapped); + if (!images) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } // update additional metadata Index: trunk/Ohana/src/dvomerge/src/help.c =================================================================== --- trunk/Ohana/src/dvomerge/src/help.c (revision 28968) +++ trunk/Ohana/src/dvomerge/src/help.c (revision 29001) @@ -4,4 +4,5 @@ fprintf (stderr, "USAGE: dvomerge (input1) and (input2) to (output)\n"); fprintf (stderr, " OR: dvomerge (input) into (output)\n"); + fprintf (stderr, " [-region Rmin Rmax Dmin Dmax]\n"); exit (2); } @@ -10,4 +11,11 @@ fprintf (stderr, "USAGE: dvoconvert (input) to (output)\n"); + + exit (2); +} + +void dvosecfilt_usage(void) { + + fprintf (stderr, "USAGE: dvosecfilt (input) -photcodes (photcodes.txt)\n"); exit (2); @@ -57,2 +65,23 @@ } +void dvosecfilt_help (int argc, char **argv) { + + /* check for help request */ + if (get_argument (argc, argv, "-help")) goto show_help; + if (get_argument (argc, argv, "-h")) goto show_help; + return; + +show_help: + + fprintf (stderr, "USAGE\n"); + fprintf (stderr, " dvosecfilt (input) (Nsecfilt)\n\n"); + + fprintf (stderr, " change number of secfilt entries in photcode table (updates secfilt tables only)\n"); + + fprintf (stderr, " optional flags:\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/gastro2/src/rfits.c =================================================================== --- trunk/Ohana/src/gastro2/src/rfits.c (revision 28968) +++ trunk/Ohana/src/gastro2/src/rfits.c (revision 29001) @@ -15,4 +15,8 @@ stars = gfits_table_get_SMPData (&table, &Nstars, NULL); + if (!stars) { + fprintf (stderr, "ERROR: failed to read stars\n"); + exit (2); + } ALLOCATE (stardata, StarData, Nstars); Index: trunk/Ohana/src/getstar/src/dvoImageExtract.c =================================================================== --- trunk/Ohana/src/getstar/src/dvoImageExtract.c (revision 28968) +++ trunk/Ohana/src/getstar/src/dvoImageExtract.c (revision 29001) @@ -32,4 +32,8 @@ // convert database table to internal structure images = gfits_table_get_Image (&db.ftable, &Nimages, &db.swapped); + if (!images) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } matches = SelectImages (argv[1], images, Nimages, &Nmatches); Index: trunk/Ohana/src/getstar/src/dvoImageOverlaps.c =================================================================== --- trunk/Ohana/src/getstar/src/dvoImageOverlaps.c (revision 28968) +++ trunk/Ohana/src/getstar/src/dvoImageOverlaps.c (revision 29001) @@ -35,4 +35,8 @@ // convert database table to internal structure dbImages = gfits_table_get_Image (&db.ftable, &NdbImages, &db.swapped); + if (!dbImages) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } // for (i = 1; i < 2; i++) { Index: trunk/Ohana/src/getstar/src/dvoImagesAtCoords.c =================================================================== --- trunk/Ohana/src/getstar/src/dvoImagesAtCoords.c (revision 28968) +++ trunk/Ohana/src/getstar/src/dvoImagesAtCoords.c (revision 29001) @@ -55,4 +55,8 @@ // convert database table to internal structure dbImages = gfits_table_get_Image (&db.ftable, &NdbImages, &db.swapped); + if (!dbImages) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } } Index: trunk/Ohana/src/imregister/src/detsearch.c =================================================================== --- trunk/Ohana/src/imregister/src/detsearch.c (revision 28968) +++ trunk/Ohana/src/imregister/src/detsearch.c (revision 29001) @@ -36,4 +36,8 @@ if (!output.Modify && !output.Delete && !output.Altpath) gfits_db_close (&db); detrend = gfits_table_get_DetReg (&db.ftable, &Ndetrend, &db.swapped); + if (!detrend) { + fprintf (stderr, "ERROR: failed to read detrend info\n"); + exit (2); + } match = MatchCriteria (detrend, Ndetrend, &Nmatch); /* match basic criteria */ Index: trunk/Ohana/src/imregister/src/imphotsearch.c =================================================================== --- trunk/Ohana/src/imregister/src/imphotsearch.c (revision 28968) +++ trunk/Ohana/src/imregister/src/imphotsearch.c (revision 29001) @@ -33,4 +33,8 @@ image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } match = subset (image, Nimage, &Nmatch); Index: trunk/Ohana/src/imregister/src/imregtable.c =================================================================== --- trunk/Ohana/src/imregister/src/imregtable.c (revision 28968) +++ trunk/Ohana/src/imregister/src/imregtable.c (revision 29001) @@ -34,4 +34,8 @@ image = gfits_table_get_RegImage (&ftable, &Nimage, NULL); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } /* load database table */ Index: trunk/Ohana/src/imregister/src/imsearch.c =================================================================== --- trunk/Ohana/src/imregister/src/imsearch.c (revision 28968) +++ trunk/Ohana/src/imregister/src/imsearch.c (revision 29001) @@ -30,4 +30,8 @@ image = gfits_table_get_RegImage (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } match = match_criteria (image, Nimage, &Nmatch); Index: trunk/Ohana/src/imregister/src/imstatreg.c =================================================================== --- trunk/Ohana/src/imregister/src/imstatreg.c (revision 28968) +++ trunk/Ohana/src/imregister/src/imstatreg.c (revision 29001) @@ -67,4 +67,8 @@ } subset = gfits_table_get_RegImage (&temp_db.ftable, &Nsubset, &temp_db.swapped); + if (!subset) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } fprintf (stderr, "temporary database read\n"); @@ -93,4 +97,8 @@ } image = gfits_table_get_RegImage (&image_db.ftable, &Nimage, &image_db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } fprintf (stderr, "main database read\n"); Index: trunk/Ohana/src/imregister/src/photsearch.c =================================================================== --- trunk/Ohana/src/imregister/src/photsearch.c (revision 28968) +++ trunk/Ohana/src/imregister/src/photsearch.c (revision 29001) @@ -42,7 +42,15 @@ if (output.convert) { photpars_old = gfits_table_get_PhotParsOld (&db.ftable, &Nphotpars, &db.swapped); + if (!photpars_old) { + fprintf (stderr, "ERROR: failed to read photometry parameters\n"); + exit (2); + } photpars = PhotParsOld_to_PhotPars (photpars_old, Nphotpars); } else { photpars = gfits_table_get_PhotPars (&db.ftable, &Nphotpars, &db.swapped); + if (!photpars) { + fprintf (stderr, "ERROR: failed to read photometry parameters\n"); + exit (2); + } } Index: trunk/Ohana/src/imregister/src/spsearch.c =================================================================== --- trunk/Ohana/src/imregister/src/spsearch.c (revision 28968) +++ trunk/Ohana/src/imregister/src/spsearch.c (revision 29001) @@ -30,4 +30,8 @@ spectrum = gfits_table_get_Spectrum (&db.ftable, &Nspectrum, &db.swapped); + if (!spectrum) { + fprintf (stderr, "ERROR: failed to read spectrum\n"); + exit (2); + } match = match_criteria (spectrum, Nspectrum, &Nmatch); Index: trunk/Ohana/src/libautocode/def/autocode.c =================================================================== --- trunk/Ohana/src/libautocode/def/autocode.c (revision 28968) +++ trunk/Ohana/src/libautocode/def/autocode.c (revision 29001) @@ -66,5 +66,5 @@ if (Ncols != $SIZE) { fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, $SIZE); - exit (1); + return NULL; } Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 28968) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 29001) @@ -85,29 +85,36 @@ # define PHOT_MAG 0x06 /* generic magnitude; never stored */ -/* Image.code values. these are codes to note bad images */ -# define ID_IMAGE_NEW 0x0000 /* no nrphot attempted */ -# define ID_IMAGE_NOCAL 0x0001 /* used within nrphot to mean "don't apply fit" */ -# define ID_IMAGE_POOR 0x0002 /* relphot says image is bad */ -# define ID_IMAGE_SKIP 0x0004 /* external information image is bad */ -# define ID_IMAGE_FEW 0x0008 /* currently too few measurements for good value */ +/* Image.code values -- these values are 32 bit (as of PS1_V1) */ +typedef enum { + ID_IMAGE_NEW = 0x00000000, /* no calibrations yet attempted */ + ID_IMAGE_PHOTOM_NOCAL = 0x00000001, /* user-set value used within relphot: ignore */ + ID_IMAGE_PHOTOM_POOR = 0x00000002, /* relphot says image is bad (dMcal > limit) */ + ID_IMAGE_PHOTOM_SKIP = 0x00000004, /* user-set value: assert that this image has bad photometry */ + ID_IMAGE_PHOTOM_FEW = 0x00000008, /* currently too few measurements for photometry */ + ID_IMAGE_ASTROM_NOCAL = 0x00000010, /* user-set value used within relastro: ignore */ + ID_IMAGE_ASTROM_POOR = 0x00000020, /* relastro says image is bad (dR,dD > limit) */ + ID_IMAGE_ASTROM_FAIL = 0x00000040, /* relastro fit diverged, fit not applied */ + ID_IMAGE_ASTROM_SKIP = 0x00000080, /* user-set value: assert that this image has bad astrometry */ + ID_IMAGE_ASTROM_FEW = 0x00000100, /* currently too few measurements for astrometry */ +} DVOImageFlags; /* Measure.flags values -- these values are 32 bit (as of PS1_V1) */ typedef enum { - ID_MEAS_NOCAL = 0x0001, // detection ignored for this analysis (photcode, time range) -- internal only - ID_MEAS_POOR_PHOTOM = 0x0002, // detection is photometry outlier - ID_MEAS_SKIP_PHOTOM = 0x0004, // detection was ignored for photometry measurement - ID_MEAS_AREA = 0x0008, // detection near image edge - ID_MEAS_POOR_ASTROM = 0x0010, // detection is astrometry outlier - ID_MEAS_SKIP_ASTROM = 0x0020, // detection was ignored for astrometry measurement - ID_MEAS_USED_OBJ = 0x0040, // detection was used during opdate objects - ID_MEAS_USED_CHIP = 0x0080, // detection was used during update chips - ID_MEAS_BLEND_MEAS = 0x0100, // detection is within radius of multiple objects - ID_MEAS_BLEND_OBJ = 0x0200, // multiple detections within radius of object - ID_MEAS_UNDEF_3 = 0x0400, // unused - ID_MEAS_UNDEF_4 = 0x0800, // unused - ID_MEAS_BLEND_MEAS_X = 0x1000, // detection is within radius of multiple objects across catalogs - ID_MEAS_ARTIFACT = 0x2000, // detection is thought to be non-astronomical - ID_MEAS_UNDEF_5 = 0x4000, // unused - ID_MEAS_UNDEF_6 = 0x8000, // unused + ID_MEAS_NOCAL = 0x00000001, // detection ignored for this analysis (photcode, time range) -- internal only + ID_MEAS_POOR_PHOTOM = 0x00000002, // detection is photometry outlier + ID_MEAS_SKIP_PHOTOM = 0x00000004, // detection was ignored for photometry measurement + ID_MEAS_AREA = 0x00000008, // detection near image edge + ID_MEAS_POOR_ASTROM = 0x00000010, // detection is astrometry outlier + ID_MEAS_SKIP_ASTROM = 0x00000020, // detection was ignored for astrometry measurement + ID_MEAS_USED_OBJ = 0x00000040, // detection was used during opdate objects + ID_MEAS_USED_CHIP = 0x00000080, // detection was used during update chips + ID_MEAS_BLEND_MEAS = 0x00000100, // detection is within radius of multiple objects + ID_MEAS_BLEND_OBJ = 0x00000200, // multiple detections within radius of object + ID_MEAS_UNDEF_3 = 0x00000400, // unused + ID_MEAS_UNDEF_4 = 0x00000800, // unused + ID_MEAS_BLEND_MEAS_X = 0x00001000, // detection is within radius of multiple objects across catalogs + ID_MEAS_ARTIFACT = 0x00002000, // detection is thought to be non-astronomical + ID_MEAS_UNDEF_5 = 0x00004000, // unused + ID_MEAS_UNDEF_6 = 0x00008000, // unused } DVOMeasureFlags; @@ -147,4 +154,8 @@ ID_STAR_USE_PM = 0x00100000, // proper motion used (not AVE or PAR) ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) + ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) + ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) + ID_OBJ_GOOD = 0x04000000, // good-quality measurement in our data (eg, PS) + ID_OBJ_GOOD_ALT = 0x08000000, // good-quality measurement in external data (eg, 2MASS) } DVOAverageFlags; Index: trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c =================================================================== --- trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/LoadPhotcodesFITS.c (revision 29001) @@ -43,4 +43,8 @@ if (!strcmp (extname, "DVO_PHOTCODE")) { PhotCode_Elixir *photcode_elixir = gfits_table_get_PhotCode_Elixir (&db.ftable, &Ncode, &db.swapped); + if (!photcode_elixir) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_Elixir_To_Internal (photcode_elixir, Ncode); free (photcode_elixir); @@ -49,4 +53,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_ELIXIR")) { PhotCode_Elixir *photcode_elixir = gfits_table_get_PhotCode_Elixir (&db.ftable, &Ncode, &db.swapped); + if (!photcode_elixir) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_Elixir_To_Internal (photcode_elixir, Ncode); free (photcode_elixir); @@ -55,4 +63,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_DEV_1")) { PhotCode_PS1_DEV_1 *photcode_ps1_dev_1 = gfits_table_get_PhotCode_PS1_DEV_1 (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_dev_1) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_DEV_1_To_Internal (photcode_ps1_dev_1, Ncode); free (photcode_ps1_dev_1); @@ -61,4 +73,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_DEV_2")) { PhotCode_PS1_DEV_2 *photcode_ps1_dev_2 = gfits_table_get_PhotCode_PS1_DEV_2 (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_dev_2) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_DEV_2_To_Internal (photcode_ps1_dev_2, Ncode); free (photcode_ps1_dev_2); @@ -67,4 +83,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_DEV_3")) { PhotCode_PS1_DEV_3 *photcode_ps1_dev_3 = gfits_table_get_PhotCode_PS1_DEV_3 (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_dev_3) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_DEV_3_To_Internal (photcode_ps1_dev_3, Ncode); free (photcode_ps1_dev_3); @@ -73,4 +93,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_V1")) { PhotCode_PS1_V1 *photcode_ps1_v1 = gfits_table_get_PhotCode_PS1_V1 (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_v1) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_V1_To_Internal (photcode_ps1_v1, Ncode); free (photcode_ps1_v1); @@ -79,4 +103,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_V2")) { PhotCode_PS1_V2 *photcode_ps1_v2 = gfits_table_get_PhotCode_PS1_V2 (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_v2) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_V2_To_Internal (photcode_ps1_v2, Ncode); free (photcode_ps1_v2); @@ -85,4 +113,8 @@ if (!strcmp (extname, "DVO_PHOTCODE_PS1_REF")) { PhotCode_PS1_REF *photcode_ps1_ref = gfits_table_get_PhotCode_PS1_REF (&db.ftable, &Ncode, &db.swapped); + if (!photcode_ps1_ref) { + fprintf (stderr, "ERROR: failed to read photcodes\n"); + exit (2); + } photcode = PhotCode_PS1_REF_To_Internal (photcode_ps1_ref, Ncode); free (photcode_ps1_ref); Index: trunk/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c =================================================================== --- trunk/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/cmf-ps1-v1-alt.c (revision 29001) @@ -14,5 +14,5 @@ if (ftable[0].header[0].Naxis[0] != 136) { fprintf (stderr, "ERROR: wrong format for CMF_PS1_V1_Alt: "OFF_T_FMT" vs %d\n", ftable[0].header[0].Naxis[0], 136); - exit (2); + return (NULL); } Index: trunk/Ohana/src/libdvo/src/coordops.c =================================================================== --- trunk/Ohana/src/libdvo/src/coordops.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/coordops.c (revision 29001) @@ -20,5 +20,10 @@ int XY_to_LM (double *L, double *M, double x, double y, Coords *coords) { + OhanaProjection proj; + OhanaProjectionMode mode; double X, Y, X2, XY, Y2, X3, Y3; + + proj = GetProjection (coords[0].ctype); + mode = GetProjectionMode (proj); /** convert pixel coordinates to cartesian system **/ @@ -45,4 +50,9 @@ } + if (mode == PROJ_MODE_CARTESIAN) { + *L += coords[0].crval1; + *M += coords[0].crval2; + } + return (TRUE); } @@ -64,11 +74,8 @@ /** Locally Cartesian Projections **/ if (mode == PROJ_MODE_CARTESIAN) { - *ra = L + coords[0].crval1; - *dec = M + coords[0].crval2; - /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */ if (proj == PROJ_WRP) { if (!gotMosaic) return (FALSE); - XY_to_RD (ra, dec, *ra, *dec, &mosaic); + XY_to_RD (ra, dec, L, M, &mosaic); } return (TRUE); @@ -183,5 +190,4 @@ double phi, theta; - double Lo, Mo; double sphi, cphi, stht, ctht; double salp, calp, sdel, cdel, sdp, cdp; @@ -201,11 +207,9 @@ if (proj == PROJ_WRP) { if (!gotMosaic) return (FALSE); - RD_to_XY (&Lo, &Mo, ra, dec, &mosaic); - *L = (Lo - coords[0].crval1); - *M = (Mo - coords[0].crval2); + RD_to_XY (L, M, ra, dec, &mosaic); return (TRUE); } - *L = (ra - coords[0].crval1); - *M = (dec - coords[0].crval2); + *L = ra; + *M = dec; return (TRUE); } @@ -310,7 +314,17 @@ double dX, dY, Lo, Mo, dL, dM; double dLdX, dLdY, dMdX, dMdY, Do; + OhanaProjection proj; + OhanaProjectionMode mode; + + proj = GetProjection (coords[0].ctype); + mode = GetProjectionMode (proj); *x = 0; *y = 0; + + if (mode == PROJ_MODE_CARTESIAN) { + L -= coords[0].crval1; + M -= coords[0].crval2; + } /* start with linear solution for X,Y */ @@ -422,4 +436,6 @@ } +enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN}; + int GetCoords (Coords *coords, Header *header) { @@ -428,4 +444,5 @@ double equinox; char *ctype; + int mode; rotate = 0.0; @@ -437,25 +454,45 @@ strcpy (coords[0].ctype, "NONE"); - status = FALSE; - if (gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) { - status = gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); - status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); - status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); - status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); - - if (gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1)) { - status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2); - if (gfits_scan (header, "CROTA2", "%lf", 1, &rotate)) { - Lambda = coords[0].cdelt2 / coords[0].cdelt1; - coords[0].pc1_1 = cos(rotate*RAD_DEG); - coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda; - coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda; - coords[0].pc2_2 = cos(rotate*RAD_DEG); - } - if (gfits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1)) { - status &= gfits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2); - status &= gfits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1); - status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2); - } + mode = COORD_TYPE_NONE; + { + int haveCTYPE, haveCDELT, haveCROTA, haveCDij, havePCij, haveRAo; + float tmp; + char stmp[80]; + + // there are a few different representations for scale and rotation. choose an appropriate + // set: (CDELTi + CROTAi), (CDELTi + PCij), (CDij), + + haveCTYPE = gfits_scan (header, "CTYPE2", "%s", 1, stmp); + haveCDELT = gfits_scan (header, "CDELT1", "%f", 1, &tmp); + haveCROTA = gfits_scan (header, "CROTA1", "%f", 1, &tmp); + haveCDij = gfits_scan (header, "CD1_1", "%f", 1, &tmp); + havePCij = gfits_scan (header, "PC001001", "%f", 1, &tmp); + haveRAo = gfits_scan (header, "RA_O", "%f", 1, &tmp); + + if (haveCTYPE && havePCij && haveCDELT) { mode = COORD_TYPE_PC; goto gotit; } + if (haveCTYPE && haveCROTA && haveCDELT) { mode = COORD_TYPE_ROT; goto gotit; } + if (haveCTYPE && haveCDij) { mode = COORD_TYPE_CD; goto gotit; } + if (haveRAo) { mode = COORD_TYPE_LIN; goto gotit; } + // fprintf (stderr, "no valid WCS keywords\n"); + return (FALSE); + } + +gotit: + + status = TRUE; + switch (mode) { + case COORD_TYPE_PC: + status &= gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype); + status &= gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); + status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); + + status &= gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1); + status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2); + status &= gfits_scan (header, "PC001001", "%f", 1, &coords[0].pc1_1); + status &= gfits_scan (header, "PC001002", "%f", 1, &coords[0].pc1_2); + status &= gfits_scan (header, "PC002001", "%f", 1, &coords[0].pc2_1); + status &= gfits_scan (header, "PC002002", "%f", 1, &coords[0].pc2_2); /* set NPLYTERM based on header. if NPLYTERM is missing, it should have a @@ -491,24 +528,48 @@ break; } - } else { - if (gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1)) { - status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2); - status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1); - status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2); - /* renormalize */ - scale = hypot (coords[0].pc1_1, coords[0].pc1_2); - coords[0].cdelt1 = coords[0].cdelt2 = scale; - coords[0].pc1_1 /= scale; - coords[0].pc1_2 /= scale; - coords[0].pc2_1 /= scale; - coords[0].pc2_2 /= scale; - } else { - status = FALSE; - } - } - } else { - /* some of my thesis data uses this simple linear model - convert on read? */ - if (gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) { - status = gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1); + break; + + case COORD_TYPE_ROT: + status &= gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype); + status &= gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); + status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); + + status &= gfits_scan (header, "CDELT1", "%f", 1, &coords[0].cdelt1); + status &= gfits_scan (header, "CDELT2", "%f", 1, &coords[0].cdelt2); + + status &= gfits_scan (header, "CROTA2", "%lf", 1, &rotate); + Lambda = coords[0].cdelt2 / coords[0].cdelt1; + coords[0].pc1_1 = cos(rotate*RAD_DEG); + coords[0].pc1_2 = -sin(rotate*RAD_DEG) * Lambda; + coords[0].pc2_1 = sin(rotate*RAD_DEG) / Lambda; + coords[0].pc2_2 = cos(rotate*RAD_DEG); + break; + + case COORD_TYPE_CD: + status &= gfits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype); + status &= gfits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= gfits_scan (header, "CRPIX1", "%f", 1, &coords[0].crpix1); + status &= gfits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= gfits_scan (header, "CRPIX2", "%f", 1, &coords[0].crpix2); + + status &= gfits_scan (header, "CD1_1", "%f", 1, &coords[0].pc1_1); + status &= gfits_scan (header, "CD1_2", "%f", 1, &coords[0].pc1_2); + status &= gfits_scan (header, "CD2_1", "%f", 1, &coords[0].pc2_1); + status &= gfits_scan (header, "CD2_2", "%f", 1, &coords[0].pc2_2); + /* renormalize */ + scale = hypot (coords[0].pc1_1, coords[0].pc1_2); + coords[0].cdelt1 = coords[0].cdelt2 = scale; + coords[0].pc1_1 /= scale; + coords[0].pc1_2 /= scale; + coords[0].pc2_1 /= scale; + coords[0].pc2_2 /= scale; + break; + + case COORD_TYPE_LIN: + /* some of my thesis data uses this simple linear model - convert on read? */ + status &= gfits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1); + status &= gfits_scan (header, "RA_X", "%f", 1, &coords[0].pc1_1); status &= gfits_scan (header, "RA_Y", "%f", 1, &coords[0].pc1_2); status &= gfits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); @@ -518,6 +579,7 @@ coords[0].cdelt1 = coords[0].cdelt2 = 1.0; strcpy (coords[0].ctype, "GENE"); - } - } + break; + } + if (status) { if (!gfits_scan (header, "EQUINOX", "%lf", 1, &equinox)) { @@ -530,6 +592,7 @@ } } + if (!status) { - fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype); + // fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype); coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; Index: trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 29001) @@ -125,4 +125,8 @@ /* no conversions currently defined */ catalog[0].missing = gfits_table_get_Missing (&ftable, &catalog[0].Nmissing, NULL); + if (!catalog[0].missing) { + fprintf (stderr, "ERROR: failed to read missing\n"); + exit (2); + } if (Nmissing != catalog[0].Nmiss_disk) { fprintf (stderr, "Warning: mismatch between Nmissing in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nmissing, catalog[0].Nmiss_disk); Index: trunk/Ohana/src/libdvo/src/dvo_catalog_split.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 29001) @@ -285,4 +285,8 @@ /* no conversions currently defined : this just does the byte swap */ catalog[0].missing = gfits_table_get_Missing (&ftable, &Nmissing, NULL); + if (!catalog[0].missing) { + fprintf (stderr, "ERROR: failed to read missing\n"); + exit (2); + } if (Nmissing != catalog[0].Nmiss_disk) { fprintf (stderr, "Warning: mismatch between Nmissing in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nmissing, catalog[0].Nmiss_disk); @@ -473,4 +477,8 @@ /* no conversions currently defined : this just does the byte swap */ catalog[0].missing = gfits_table_get_Missing (&ftable, &Nmissing, NULL); + if (!catalog[0].missing) { + fprintf (stderr, "ERROR: failed to read missing\n"); + exit (2); + } if (Nmissing != Nrows) { fprintf (stderr, "Warning: mismatch between Nmissing in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nmissing, Nrows); Index: trunk/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 29001) @@ -32,4 +32,8 @@ Average_##TYPE *tmpAverage; \ tmpAverage = gfits_table_get_Average_##TYPE (ftable, Naverage, NULL); \ + if (!tmpAverage) { \ + fprintf (stderr, "ERROR: failed to read averages\n"); \ + exit (2); \ + } \ average = Average_##TYPE##_ToInternal (tmpAverage, *Naverage, primary); \ free (tmpAverage); \ @@ -39,4 +43,8 @@ if (!strcmp (extname, "DVO_AVERAGE")) { average = gfits_table_get_Average (ftable, Naverage, NULL); + if (!average) { + fprintf (stderr, "ERROR: failed to read averages\n"); + exit (2); + } *format = DVO_FORMAT_INTERNAL; return (average); @@ -111,4 +119,8 @@ Measure_##TYPE *tmpMeasure; \ tmpMeasure = gfits_table_get_Measure_##TYPE (ftable, Nmeasure, NULL); \ + if (!tmpMeasure) { \ + fprintf (stderr, "ERROR: failed to read measures\n"); \ + exit (2); \ + } \ measure = Measure_##TYPE##_ToInternal (tmpMeasure, *Nmeasure); \ free (tmpMeasure); \ @@ -118,4 +130,8 @@ if (!strcmp (extname, "DVO_MEASURE")) { measure = gfits_table_get_Measure (ftable, Nmeasure, NULL); + if (!measure) { + fprintf (stderr, "ERROR: failed to read measures\n"); + exit (2); + } *format = DVO_FORMAT_INTERNAL; return (measure); @@ -190,4 +206,8 @@ SecFilt_##TYPE *tmpSecFilt; \ tmpSecFilt = gfits_table_get_SecFilt_##TYPE (ftable, Nsecfilt, NULL); \ + if (!tmpSecFilt) { \ + fprintf (stderr, "ERROR: failed to read secfilts\n"); \ + exit (2); \ + } \ secfilt = SecFilt_##TYPE##_ToInternal (tmpSecFilt, *Nsecfilt); \ free (tmpSecFilt); \ @@ -197,4 +217,8 @@ if (!strcmp (extname, "DVO_SECFILT")) { secfilt = gfits_table_get_SecFilt (ftable, Nsecfilt, NULL); + if (!secfilt) { + fprintf (stderr, "ERROR: failed to read secfilts\n"); + exit (2); + } *format = DVO_FORMAT_INTERNAL; return (secfilt); @@ -263,4 +287,8 @@ Image_Elixir *tmpimage; tmpimage = gfits_table_get_Image_Elixir (ftable, &Nimage, NULL); + if (!tmpimage) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } ftable[0].buffer = (char *) Image_Elixir_ToInternal (tmpimage, Nimage); free (tmpimage); @@ -284,4 +312,8 @@ *format = DVO_FORMAT_##FORMAT; \ tmpimage = gfits_table_get_Image_##TYPE (ftable, &Nimage, NULL); \ + if (!tmpimage) { \ + fprintf (stderr, "ERROR: failed to read images\n"); \ + exit (2); \ + } \ ftable[0].buffer = (char *) Image_##TYPE##_ToInternal (tmpimage, Nimage); \ free (tmpimage); \ @@ -296,4 +328,8 @@ Image *image; image = gfits_table_get_Image (ftable, &Nimage, NULL); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } *format = DVO_FORMAT_INTERNAL; return (TRUE); Index: trunk/Ohana/src/libdvo/src/skyregion_io.c =================================================================== --- trunk/Ohana/src/libdvo/src/skyregion_io.c (revision 28968) +++ trunk/Ohana/src/libdvo/src/skyregion_io.c (revision 29001) @@ -40,4 +40,9 @@ ALLOCATE (skytable, SkyTable, 1); skytable[0].regions = gfits_table_get_SkyRegion (&ftable, &skytable[0].Nregions, NULL); + if (!skytable[0].regions) { + fprintf (stderr, "ERROR: failed to read sky regions\n"); + exit (2); + } + ALLOCATE (skytable[0].filename, char *, skytable[0].Nregions); for (i = 0; i < skytable[0].Nregions; i++) { Index: trunk/Ohana/src/markstar/src/find_images.c =================================================================== --- trunk/Ohana/src/markstar/src/find_images.c (revision 28968) +++ trunk/Ohana/src/markstar/src/find_images.c (revision 29001) @@ -26,4 +26,8 @@ timage = gfits_table_get_Image (&db[0].ftable, &Ntimage, &db[0].swapped); + if (!timage) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } /* set up buffers for images, temporary storage */ Index: trunk/Ohana/src/mosastro/src/rfits.c =================================================================== --- trunk/Ohana/src/mosastro/src/rfits.c (revision 28968) +++ trunk/Ohana/src/mosastro/src/rfits.c (revision 29001) @@ -32,4 +32,9 @@ off_t Nstars; stars = gfits_table_get_SMPData (&table, &Nstars, NULL); + if (!stars) { + fprintf (stderr, "ERROR: failed to read stars\n"); + exit (2); + } + mychip[0].Nstars = Nstars; gfits_scan (table.header, "NAXIS1", "%d", 1, &Nx); Index: trunk/Ohana/src/opihi/cmd.astro/cgrid.c =================================================================== --- trunk/Ohana/src/opihi/cmd.astro/cgrid.c (revision 28968) +++ trunk/Ohana/src/opihi/cmd.astro/cgrid.c (revision 29001) @@ -153,4 +153,12 @@ SetVector (&Yvec, OPIHI_FLT, NELEMENTS); + { // sanity check + float Nelem; + Nelem = 180.0 / minorRA; + if (!isfinite(Nelem) || (fabs(Nelem) > 10000)) { fprintf (stderr, "absurd cgrid spacing %f\n", minorRA); return FALSE; } + Nelem = 180.0 / minorDEC; + if (!isfinite(Nelem) || (fabs(Nelem) > 10000)) { fprintf (stderr, "absurd cgrid spacing %f\n", minorDEC); return FALSE; } + } + /*** do consecutive RA lines, first increasing **/ OnPic = TRUE; Index: trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- trunk/Ohana/src/opihi/cmd.data/Makefile (revision 28968) +++ trunk/Ohana/src/opihi/cmd.data/Makefile (revision 29001) @@ -19,6 +19,7 @@ $(SRC)/init.$(ARCH).o \ $(SRC)/accum.$(ARCH).o \ +$(SRC)/applyfit1d.$(ARCH).o \ $(SRC)/applyfit2d.$(ARCH).o \ -$(SRC)/applyfit1d.$(ARCH).o \ +$(SRC)/applyfit3d.$(ARCH).o \ $(SRC)/box.$(ARCH).o \ $(SRC)/book.$(ARCH).o \ @@ -46,4 +47,5 @@ $(SRC)/fit1d.$(ARCH).o \ $(SRC)/fit2d.$(ARCH).o \ +$(SRC)/fit3d.$(ARCH).o \ $(SRC)/gaussj.$(ARCH).o \ $(SRC)/gaussdeviate.$(ARCH).o \ Index: trunk/Ohana/src/opihi/cmd.data/applyfit2d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 28968) +++ trunk/Ohana/src/opihi/cmd.data/applyfit2d.c (revision 29001) @@ -10,5 +10,5 @@ if (argc != 4) { - gprint (GP_ERR, "USAGE: applyfit x y z\n"); + gprint (GP_ERR, "USAGE: applyfit2d x y z\n"); return (FALSE); } Index: trunk/Ohana/src/opihi/cmd.data/applyfit3d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/applyfit3d.c (revision 29001) +++ trunk/Ohana/src/opihi/cmd.data/applyfit3d.c (revision 29001) @@ -0,0 +1,94 @@ +# include "data.h" + +int applyfit3d (int argc, char **argv) { + + int i, j, k, io, jo, n, order, nterm; + char *c, name[64]; + double ***C, X, Y, Z; + opihi_flt *F; + Vector *xvec, *yvec, *zvec, *Fvec; + + if (argc != 5) { + gprint (GP_ERR, "USAGE: applyfit x y z F\n"); + return (FALSE); + } + + c = get_variable ("Cnnn"); + if (c == NULL) { + gprint (GP_ERR, "no fit available\n"); + return (FALSE); + } + order = atof (c); + free (c); + + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((Fvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + // this is a bit crude, but perhaps needed + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + CastVector (zvec, OPIHI_FLT); + CastVector (Fvec, OPIHI_FLT); + + nterm = order + 1; + ALLOCATE (C, double **, nterm); + for (i = 0; i < order + 1; i++) { + ALLOCATE (C[i], double *, nterm); + for (j = 0; j < nterm - i; j++) { + ALLOCATE (C[i][j], double, nterm); + for (k = 0; k < nterm - i - j; k++) { + sprintf (name, "CX%dY%dZ%d", i, j, k); + c = get_variable (name); + if (c == NULL) { + gprint (GP_ERR, "missing fit term %d,%d,%d\n", i, j, k); + for (io = 0; io < i; io++) { + for (jo = 0; jo < j; jo++) { + free (C[io][jo]); + } + free (C[io]); + } + free (C); + return (FALSE); + } + C[i][j][k] = atof (c); + free (c); + } + } + } + + ResetVector (Fvec, OPIHI_FLT, xvec[0].Nelements); + bzero (Fvec[0].elements.Flt, sizeof(opihi_flt)*zvec[0].Nelements); + F = Fvec[0].elements.Flt; + + opihi_flt *x = xvec[0].elements.Flt; + opihi_flt *y = yvec[0].elements.Flt; + opihi_flt *z = zvec[0].elements.Flt; + + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++, F++) { + Z = Y = X = 1; + for (k = 0; k < nterm; k++) { + Y = Z; + for (j = 0; j < nterm - k; j++) { + X = Y; + for (i = 0; i < nterm - k - j; i++) { + *F += C[i][j][k]*X; + X = X * (*x); + } + Y = Y * (*y); + } + Z = Z * (*z); + } + } + + for (i = 0; i < nterm; i++) { + for (j = 0; j < nterm - i; j++) { + free (C[i][j]); + } + free (C[i]); + } + free (C); + + return (TRUE); +} Index: trunk/Ohana/src/opihi/cmd.data/fit2d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/fit2d.c (revision 28968) +++ trunk/Ohana/src/opihi/cmd.data/fit2d.c (revision 29001) @@ -168,4 +168,15 @@ } + /** test print **/ + if (VERBOSE) { + for (i = 0; i < wterm; i++) { + for (j = 0; j < wterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); + } + gprint (GP_ERR, "-----\n"); + } + dgaussjordan (c, b, wterm, 1); @@ -176,4 +187,14 @@ gprint (GP_ERR, "x^%d y^%d: %g\n", nx, ny, b[i][0]); } + } + } + + /** test print **/ + if (VERBOSE) { + for (i = 0; i < wterm; i++) { + for (j = 0; j < wterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); } } Index: trunk/Ohana/src/opihi/cmd.data/fit3d.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/fit3d.c (revision 29001) +++ trunk/Ohana/src/opihi/cmd.data/fit3d.c (revision 29001) @@ -0,0 +1,352 @@ +# include "data.h" + +int fit3d (int argc, char **argv) { + + double **c, **b, ***s, X, Y, Z, dF1, dF2; + double ClipNSigma, mean, sigma, maxsigma; + int ix, iy, i, j, k, n, Npt, Nmask, nx, ny, nz, nterm, mterm, wterm, order; + int N, Weight, Quiet, ClipNiter, VERBOSE; + opihi_flt *x, *y, *z, *F, *dFv, *Ffit, *Ff; + char name[64], *mask; + Vector *xvec, *yvec, *zvec, *Fvec, *dFvec; + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + Quiet = FALSE; + if ((N = get_argument (argc, argv, "-q"))) { + Quiet = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-quiet"))) { + Quiet = TRUE; + remove_argument (N, &argc, argv); + } + + ClipNSigma = 0; + ClipNiter = 1; + if ((N = get_argument (argc, argv, "-clip"))) { + remove_argument (N, &argc, argv); + ClipNSigma = atof(argv[N]); + remove_argument (N, &argc, argv); + ClipNiter = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + dFv = NULL; + dFvec = NULL; + Weight = FALSE; + if ((N = get_argument (argc, argv, "-dF"))) { + remove_argument (N, &argc, argv); + if ((dFvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + Weight = TRUE; + } + + if (argc != 6) { + gprint (GP_ERR, "USAGE: fit x y z F order [-dF wt]\n"); + return (FALSE); + } + + if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((zvec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE); + if ((Fvec = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE); + + if (xvec[0].Nelements != yvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + if (xvec[0].Nelements != zvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + if (xvec[0].Nelements != Fvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + CastVector (zvec, OPIHI_FLT); + CastVector (Fvec, OPIHI_FLT); + + if (Weight) { + if (xvec[0].Nelements != dFvec[0].Nelements) { + gprint (GP_ERR, "vectors must have same length\n"); + return (FALSE); + } + CastVector (dFvec, OPIHI_FLT); + } + + order = atof (argv[5]); + nterm = order + 1; + wterm = (order + 3)*(order + 2)*(order + 1)/6; + mterm = (order + 2)*(order + 1)/2; + + ALLOCATE (Ffit, opihi_flt, xvec[0].Nelements); + ALLOCATE (mask, char, xvec[0].Nelements); + memset (mask, 0, xvec[0].Nelements); + + /* allocate the summation matrices */ + ALLOCATE (b, double *, wterm); + ALLOCATE (c, double *, wterm); + for (i = 0; i < wterm; i++) { + ALLOCATE (c[i], double, wterm); + ALLOCATE (b[i], double, 1); + } + + /* we initially sum the values in a cube */ + ALLOCATE (s, double **, mterm); + for (i = 0; i < mterm; i++) { + ALLOCATE (s[i], double *, mterm); + for (j = 0; j < mterm; j++) { + ALLOCATE (s[i][j], double, mterm); + } + } + + for (N = 0; N < ClipNiter; N++) { + + /* init registers for current pass */ + for (i = 0; i < wterm; i++) { + memset (c[i], 0, wterm*sizeof(double)); + memset (b[i], 0, sizeof(double)); + } + for (i = 0; i < mterm; i++) { + for (j = 0; j < mterm; j++) { + memset (s[i][j], 0, mterm*sizeof(double)); + } + } + + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + F = Fvec[0].elements.Flt; + if (Weight) dFv = dFvec[0].elements.Flt; + + /* add up the x,y,z values */ + for (i = 0; i < xvec[0].Nelements; i++, x++, y++, z++, F++) { + if (mask[i]) continue; + if (!finite(*x)) continue; + if (!finite(*y)) continue; + if (!finite(*z)) continue; + if (!finite(*F)) continue; + dF1 = 1.0; + if (Weight) { + dF1 = 1.0 / SQ(*dFv); + dFv ++; + } + X = Y = Z = dF1; + for (nz = 0; nz < mterm; nz++) { + Y = Z; + for (ny = 0; ny < mterm - nz; ny++) { + X = Y; + for (nx = 0; nx < mterm - nz - ny; nx++) { + s[nx][ny][nz] += X; + X = X * (*x); + } + Y = Y * (*y); + } + Z = Z * (*z); + } + } + + /* add up the F values */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + F = Fvec[0].elements.Flt; + for (i = 0; i < xvec[0].Nelements; i++, x++, y++, z++, F++) { + if (mask[i]) continue; + if (!finite(*x)) continue; + if (!finite(*y)) continue; + if (!finite(*z)) continue; + if (!finite(*F)) continue; + dF1 = 1.0; + if (Weight) { + dF1 = 1.0 / SQ(*dFv); + dFv ++; + } + X = Y = Z = *F*dF1; + j = 0; + for (nz = 0; nz < nterm; nz++) { + Y = Z; + for (ny = 0; ny < nterm - nz; ny++) { + X = Y; + for (nx = 0; nx < nterm - ny - nz; nx++, j++) { + b[j][0] += X; + X = X * (*x); + } + Y = Y * (*y); + } + Z = Z * (*z); + } + } + + /* re-sort mterm x mterm matrix to wterm matrix */ + ix = 0; + for (k = 0; k < nterm; k++) { + for (j = 0; j < nterm - k; j++) { + for (i = 0; i < nterm - k - j; i++, ix++) { + iy = 0; + for (nz = 0; nz < nterm; nz++) { + for (ny = 0; ny < nterm - nz; ny++) { + for (nx = 0; nx < nterm - ny - nz; nx++, iy++) { + c[ix][iy] = s[nx+i][ny+j][nz+k]; + } + } + } + } + } + } + + /** test print **/ + if (VERBOSE) { + for (i = 0; i < wterm; i++) { + for (j = 0; j < wterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); + } + gprint (GP_ERR, "-----\n"); + } + + dgaussjordan (c, b, wterm, 1); + + /** test print **/ + if (VERBOSE) { + i = 0; + for (nz = 0; nz < nterm; nz++) { + for (ny = 0; ny < nterm - nz; ny++) { + for (nx = 0; nx < nterm - nz - ny; nx++, i++) { + gprint (GP_ERR, "x^%d y^%d z^%d: %g\n", nx, ny, nz, b[i][0]); + } + } + } + } + + /** test print **/ + if (VERBOSE) { + for (i = 0; i < wterm; i++) { + for (j = 0; j < wterm; j++) { + gprint (GP_ERR, "%g ", c[i][j]); + } + gprint (GP_ERR, "\n"); + } + } + + /* the b[][0] terms are in the following order: + y^0 x^0, y^0 x^1, ... y^0 x^N + y^1 x^0, y^1 x^1, ... y^1 x^N + ... + y^N x^0, y^N x^1, ... y^N x^N + */ + /* generate fitted values */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + Ff = Ffit; + for (n = 0; n < xvec[0].Nelements; n++, x++, y++, z++, Ff++) { + if (!finite(*x)) continue; + if (!finite(*y)) continue; + if (!finite(*z)) continue; + if (!finite(*F)) continue; + *Ff = 0; + Z = Y = X = 1; + i = 0; + for (nz = 0; nz < nterm; nz++) { + Y = Z; + for (ny = 0; ny < nterm - nz; ny++) { + X = Y; + for (nx = 0; nx < nterm - nz - ny; nx++, i++) { + *Ff += b[i][0]*X; + X = X * (*x); + } + Y = Y * (*y); + } + Z = Z * (*z); + } + } + + /* measure fit residual scatter */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + F = Fvec[0].elements.Flt; + Ff = Ffit; + dF1 = dF2 = 0; + for (i = Npt = 0; i < xvec[0].Nelements; i++, x++, y++, z++, F++, Ff++) { + if (mask[i]) continue; + if (!finite(*x)) continue; + if (!finite(*y)) continue; + if (!finite(*z)) continue; + if (!finite(*F)) continue; + dF1 += (*F - *Ff); + dF2 += SQ(*F - *Ff); + Npt ++; + } + mean = dF1 / Npt; + sigma = sqrt (fabs(dF2/Npt - SQ(mean))); + maxsigma = ClipNSigma * sigma; + + if (VERBOSE) gprint (GP_ERR, "mean: %g, sigma: %g, maxsigma: %g\n", mean, sigma, maxsigma); + + /* mask outlier points */ + x = xvec[0].elements.Flt; + y = yvec[0].elements.Flt; + z = zvec[0].elements.Flt; + F = Fvec[0].elements.Flt; + Ff = Ffit; + Nmask = 0; + for (i = 0; ClipNSigma && (i < xvec[0].Nelements); i++, x++, y++, z++, Ff++) { + dF1 = (*F - *Ff); + if (fabs(dF1) > maxsigma) { + mask[i] = TRUE; + Nmask ++; + } else { + mask[i] = FALSE; + } + } + if (VERBOSE) gprint (GP_ERR, "pass: %d, Nmask: %d\n", N, Nmask); + } + + if (!Quiet) gprint (GP_ERR, "F = "); + + N = 0; + for (nz = 0; nz < nterm; nz++) { + for (ny = 0; ny < nterm - nz; ny++) { + for (nx = 0; nx < nterm - nz - ny; nx++, N++) { + sprintf (name, "CX%dY%dZ%d", nx, ny, nz); + set_variable (name, b[N][0]); + if (!Quiet) gprint (GP_ERR, "%f x^%d y^%d z^%d ", b[N][0], nx, ny, nz); + } + } + } + sprintf (name, "Cnnn"); + set_variable (name, (double) order); + + if (!Quiet) gprint (GP_ERR, "\n"); + + /* free internal data */ + free (Ffit); + free (mask); + + for (i = 0; i < wterm; i++) { + free (c[i]); + free (b[i]); + } + free (b); + free (c); + + for (i = 0; i < mterm; i++) { + for (j = 0; j < mterm; j++) { + free (s[i][j]); + } + free (s[i]); + } + free (s); + + return (TRUE); +} Index: trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/init.c (revision 28968) +++ trunk/Ohana/src/opihi/cmd.data/init.c (revision 29001) @@ -5,4 +5,5 @@ int applyfit1d PROTO((int, char **)); int applyfit2d PROTO((int, char **)); +int applyfit3d PROTO((int, char **)); int box PROTO((int, char **)); int book_command PROTO((int, char **)); @@ -34,4 +35,5 @@ int fit1d PROTO((int, char **)); int fit2d PROTO((int, char **)); +int fit3d PROTO((int, char **)); int gaussjordan PROTO((int, char **)); int gaussdeviate PROTO((int, char **)); @@ -138,4 +140,5 @@ {1, "applyfit1d", applyfit1d, "apply 1-d fit to new vector"}, {1, "applyfit2d", applyfit2d, "apply 2-d fit to new vector"}, + {1, "applyfit3d", applyfit3d, "apply 3-d fit to new vector"}, {1, "book", book_command, "commands to manipulate book/page/word data"}, {1, "box", box, "draw a box on the plot"}, @@ -168,4 +171,5 @@ {1, "fit1d", fit1d, "fit polynomial to vector pair"}, {1, "fit2d", fit2d, "fit 2-d polynomial to vector triplet"}, + {1, "fit3d", fit3d, "fit 3-d polynomial to vector quad"}, {1, "gaussdev", gaussdeviate, "generate a gaussian deviate vector"}, {1, "gaussint", gaussintegral, "return the integrated gaussian vector"}, Index: trunk/Ohana/src/opihi/cmd.data/vgauss.c =================================================================== --- trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 28968) +++ trunk/Ohana/src/opihi/cmd.data/vgauss.c (revision 29001) @@ -42,7 +42,7 @@ if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE); - REQUIRE_VECTOR_FLT (xvec, FALSE); - REQUIRE_VECTOR_FLT (yvec, FALSE); - REQUIRE_VECTOR_FLT (svec, FALSE); + CastVector (xvec, OPIHI_FLT); + CastVector (yvec, OPIHI_FLT); + CastVector (svec, OPIHI_FLT); Npts = xvec[0].Nelements; Index: trunk/Ohana/src/opihi/dvo/LoadImages.c =================================================================== --- trunk/Ohana/src/opihi/dvo/LoadImages.c (revision 28968) +++ trunk/Ohana/src/opihi/dvo/LoadImages.c (revision 29001) @@ -65,4 +65,8 @@ image = gfits_table_get_Image (&db.ftable, Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + return (NULL); + } if (dvoUseImageCache && image) { cacheNimage = *Nimage; @@ -71,5 +75,4 @@ lastModified = getLastModified(filename); } - return (image); } Index: trunk/Ohana/src/opihi/dvo/cmpReadFile.c =================================================================== --- trunk/Ohana/src/opihi/dvo/cmpReadFile.c (revision 28968) +++ trunk/Ohana/src/opihi/dvo/cmpReadFile.c (revision 29001) @@ -20,4 +20,8 @@ smpdata = gfits_table_get_SMPData (&table, &Nstars, NULL); + if (!smpdata) { + fprintf (stderr, "ERROR: failed to read stars\n"); + exit (2); + } ALLOCATE (stars, CMPstars, Nstars); Index: trunk/Ohana/src/photdbc/src/copy_images.c =================================================================== --- trunk/Ohana/src/photdbc/src/copy_images.c (revision 28968) +++ trunk/Ohana/src/photdbc/src/copy_images.c (revision 29001) @@ -46,4 +46,8 @@ image = gfits_table_get_Image (&in.ftable, &Nimage, &in.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } dvo_image_addrows (&out, image, Nimage); Index: trunk/Ohana/src/photdbc/src/find_images.c =================================================================== --- trunk/Ohana/src/photdbc/src/find_images.c (revision 28968) +++ trunk/Ohana/src/photdbc/src/find_images.c (revision 29001) @@ -28,4 +28,8 @@ timage = gfits_table_get_Image (&db[0].ftable, &Ntimage, &db[0].swapped); + if (!timage) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } nimage = 0; @@ -97,6 +101,6 @@ if (found) { image[nimage] = timage[i]; - if (image[nimage].code == ID_IMAGE_NOCAL) { - image[nimage].code &= ~ID_IMAGE_NOCAL; + if (image[nimage].code == ID_IMAGE_PHOTOM_NOCAL) { + image[nimage].code &= ~ID_IMAGE_PHOTOM_NOCAL; } line_number[nimage] = i; Index: trunk/Ohana/src/photdbc/src/flag_measures.c =================================================================== --- trunk/Ohana/src/photdbc/src/flag_measures.c (revision 28968) +++ trunk/Ohana/src/photdbc/src/flag_measures.c (revision 29001) @@ -16,4 +16,8 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } PhotNsec = GetPhotcodeNsecfilt (); Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 28968) +++ trunk/Ohana/src/relastro/Makefile (revision 29001) @@ -58,4 +58,5 @@ $(SRC)/CoordOps.$(ARCH).o \ $(SRC)/FixProblemImages.$(ARCH).o \ +$(SRC)/StarMaps.$(ARCH).o \ $(SRC)/high_speed_catalogs.$(ARCH).o \ $(SRC)/high_speed_objects.$(ARCH).o \ Index: trunk/Ohana/src/relastro/doc/notes.txt =================================================================== --- trunk/Ohana/src/relastro/doc/notes.txt (revision 28968) +++ trunk/Ohana/src/relastro/doc/notes.txt (revision 29001) @@ -1,2 +1,28 @@ + +2010.07.15 + + Some issues: + + * when I apply the fitted image parameters (UpdateObjectOffsets), I + currently check for good/bad fits and warn if the image is bad. + But the correction requires multiple re-loads of the catalog + tables. + + * I would like to check for bad fits by defining a grid of test + points and seeing how much they move. However, for some chips, + the region of valid data is much smaller than the regino with any + data. + + * I would like to define a grid (say, Nx * Ny with Nx >> polynomial + order) and only use those grid regions in testing the solution + + * This means knowing which detections belong on a specific image. + I'm not sure if there is a good way to do this. at the moment, it + is easy to do detection -> image, but not so easy to do + image->detection. + + * That might argue for generating the grid for each image, then + running through all detections and, as detections are examined, + populate the grid elements for their own image. 2010.03.24 Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 28968) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 29001) @@ -3,4 +3,5 @@ # include # include +# include typedef enum { @@ -76,4 +77,10 @@ int Nmeas; } StatType; + +# define MARKTIME(MSG,...) { \ + float dtime; \ + gettimeofday (&stop, (void *) NULL); \ + dtime = DTIME (stop, start); \ + fprintf (stderr, MSG, __VA_ARGS__); } /* global variables set in parameter file */ @@ -145,4 +152,10 @@ int DoUpdateChips; int DoUpdateMosaics; + +// StarMap parameters: +int NX_MAP; +int NY_MAP; +double DPOS_MAX; +double ADDSTAR_RADIUS; /*** relphot prototypes ***/ @@ -207,4 +220,6 @@ void create_image_db (FITS_DB *db); void save_catalogs (Catalog *catalog, int Ncatalog); + +int reload_images (FITS_DB *db); int main PROTO((int argc, char **argv)); @@ -267,8 +282,8 @@ double **poly2d_copy (double **poly, int Nx, int Ny); double poly2d_eval (double **poly, int Nx, int Ny, double x, double y); -int fit_apply_coords (CoordFit *fit, Coords *coords); +int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef); int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo); CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo); -int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image); void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords); void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords); @@ -293,4 +308,5 @@ StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); +off_t getImageByID (off_t ID); Mosaic *getmosaics (off_t *N); @@ -331,2 +347,7 @@ int high_speed_objects (SkyRegion *region, Catalog *catalog); int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset); + +int initStarMaps (); +int updateStarMaps(Catalog *catalog); +int createStarMapPoints(); +int checkStarMap(int N); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 28968) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 29001) @@ -32,4 +32,9 @@ GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); + + GetConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP); + GetConfig (config, "RELASTRO_MAP_NY", "%d", 0, &NY_MAP); + GetConfig (config, "RELASTRO_DPOS_MAX", "%lf", 0, &DPOS_MAX); + GetConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS); GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 28968) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 29001) @@ -3,13 +3,6 @@ // XXX make these user parameters -# define FIT_CHIP_MAX_ERROR 0.05 # define FIT_CHIP_NITER 3 # define FIT_CHIP_NSIGMA 3.0 - - -// XXX we should test if the fit is sufficiently constrained across the chip, or if the -// new positions deviate too much from the old positions - -// XXX add visualization tools: per chip residual plots // XXX save the fit[0].Npts value in the image table? @@ -17,5 +10,5 @@ // XXX save measurements of the fit quality (scatter, chisq) in the image table -int FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { +int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) { int i, Nscatter, Niter, skip; @@ -29,5 +22,5 @@ for (i = Nscatter = 0; i < Nmatch; i++) { if (raw[i].mask) continue; - if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue; + if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) continue; dL = raw[i].L - ref[i].L; @@ -39,15 +32,20 @@ } - // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma - dsort (values, Nscatter); - dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)]; - relastroVisualPlotScatter(values, dRmax, Nscatter); - relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch); + if (Nscatter > 5) { + // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma + // XXX this test is not sensible for Nscatter < XXX (5?) + dsort (values, Nscatter); + dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)]; + relastroVisualPlotScatter(values, dRmax, Nscatter); + relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch); + } else { + dRmax = 0; + } // fit the requested order polynomial if (CHIPORDER > 0) { - coords[0].Npolyterms = CHIPORDER; + image[0].coords.Npolyterms = CHIPORDER; } - fit = fit_init (coords[0].Npolyterms); + fit = fit_init (image[0].coords.Npolyterms); // generate the fit matches @@ -56,5 +54,5 @@ continue; } - if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) { + if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) { continue; } @@ -64,5 +62,8 @@ dM = raw[i].M - ref[i].M; dR = hypot (dL, dM); - if (dR > dRmax) continue; + + // fprintf (stderr, "fit %f %f -> %f %f : %f %f (%f vs %f) wt: %f\n", raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].L, raw[i].M, dR, dRmax, raw[i].dPos); + + if ((dRmax > 0.0) && (dR > dRmax)) continue; fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); @@ -71,5 +72,5 @@ // check if the fit has enough data points for the polynomial order skip = FALSE; - switch (coords[0].Npolyterms) { + switch (image[0].coords.Npolyterms) { case 0: case 1: @@ -83,11 +84,12 @@ break; default: - fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms); + fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms); skip = TRUE; } if (skip) { - fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, coords[0].Npolyterms); + if (VERBOSE) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms); fit_free (fit); free (values); + image[0].flags |= ID_IMAGE_ASTROM_FEW; return FALSE; } @@ -98,9 +100,11 @@ if (!fit_eval (fit)) { fprintf (stderr, "failed to fit new model\n"); + image[0].flags |= ID_IMAGE_ASTROM_FAIL; return FALSE; } - if (!fit_apply_coords (fit, coords)) { + if (!fit_apply_coords (fit, &image[0].coords, FALSE)) { fprintf (stderr, "failed to fit new model\n"); + image[0].flags |= ID_IMAGE_ASTROM_FAIL; return FALSE; } @@ -109,5 +113,5 @@ for (i = 0; i < Nmatch; i++) { - XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords); + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords); } Index: trunk/Ohana/src/relastro/src/FitMosaic.c =================================================================== --- trunk/Ohana/src/relastro/src/FitMosaic.c (revision 28968) +++ trunk/Ohana/src/relastro/src/FitMosaic.c (revision 29001) @@ -26,5 +26,5 @@ } fit_eval (fit); - fit_apply_coords (fit, coords); + fit_apply_coords (fit, coords, TRUE); fit_free (fit); Index: trunk/Ohana/src/relastro/src/FitSimple.c =================================================================== --- trunk/Ohana/src/relastro/src/FitSimple.c (revision 28968) +++ trunk/Ohana/src/relastro/src/FitSimple.c (revision 29001) @@ -26,5 +26,5 @@ } fit_eval (fit); - fit_apply_coords (fit, coords); + fit_apply_coords (fit, coords, TRUE); fit_free (fit); Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 28968) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 29001) @@ -319,4 +319,15 @@ } +void dump_measures(Average *average, Measure *measure) { + + off_t j, off; + + for (j = 0; j < average[0].Nmeasure; j++) { + off = average[0].measureOffset + j; + fprintf (stderr, "%f, %f\n", measure[off].dR, measure[off].dD); + } + return; +} + // return StarData values for detections in the specified image, converting coordinates from the // chip positions: X,Y -> L,M -> P,Q -> R,D @@ -325,5 +336,5 @@ off_t i, m, c, n, nPos; double X, Y, L, M, P, Q, R, D, dR, dD; - double dPos; + double dPos, DPOS_MAX_ASEC; Mosaic *mosaic; @@ -341,4 +352,10 @@ } imcoords = &image[im].coords; + + if (moscoords) { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords[0].cdelt1, moscoords[0].cdelt2); + } else { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords[0].cdelt1, imcoords[0].cdelt2); + } // accumulate the rms position offsets. if this value, or any specific entry, is too @@ -372,25 +389,26 @@ // complain if the new location is far from the average location - if (fabs(dR) > 2.0) { - fprintf (stderr, "!"); - setBadCoords (im); // report a failure for this image - return; - } - if (fabs(dD) > 2.0) { - fprintf (stderr, "*"); - setBadCoords (im); // report a failure for this image - return; + // NOTE: This should never happen, or our StarMap tests are not working + if (fabs(dR) > 1.5*ADDSTAR_RADIUS) { + fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); + dump_measures (&catalog[c].average[n], catalog[c].measure); + // abort (); + } + if (fabs(dD) > 1.5*ADDSTAR_RADIUS) { + fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); + dump_measures (&catalog[c].average[n], catalog[c].measure); + // abort (); } // complain if the new location is far from the old location - if (fabs(catalog[c].measure[m].dR - dR) > 2.0) { - fprintf (stderr, "@"); - setBadCoords (im); // report a failure for this image - return; - } - if (fabs(catalog[c].measure[m].dD - dD) > 2.0) { - fprintf (stderr, "#"); - setBadCoords (im); // report a failure for this image - return; + if (fabs(catalog[c].measure[m].dR - dR) > DPOS_MAX_ASEC) { + fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); + dump_measures (&catalog[c].average[n], catalog[c].measure); + // abort(); + } + if (fabs(catalog[c].measure[m].dD - dD) > DPOS_MAX_ASEC) { + fprintf (stderr, "measurement is far from original location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); + dump_measures (&catalog[c].average[n], catalog[c].measure); + // abort(); } @@ -414,4 +432,6 @@ saveOffsets (dPos, nPos, im); + + // XXX this function should probably update the average positions as well. return; Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 29001) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 29001) @@ -0,0 +1,153 @@ +# include "relastro.h" + +// XXX mark the image grid based on the loaded detections +// XXX modify tcatalog to exclude any photcode or other choices that +// only affect image selection? + +// XXX this stuff is only going to work for SIMPLE and CHIP; need to +// use parentID to get from imageID to the mosaic pseudo-image + +typedef struct { + double L; + double M; + double X; + double Y; + double dL; + double dM; +} StarMapPoint; + +typedef struct { + int Nx; + int Ny; + int *stars; // arrays to count the number of stars in each map bin + StarMapPoint *points; // test points generated based on map + int Npoints; +} StarMap; + +static StarMap *starmap; +// static int Nstarmap; -- is it always == Nimages? + +/** we have two possible options for how the user defines the map bins: + a) the user specifies the number of cells in the map, leading to varying binning + b) the user specifies the binning factor, leading to variable grid sizes. + + We want to use (b) since we know the order of the model +**/ + +int initStarMaps () { + + Image *images; + off_t i, Nimages; + + images = getimages(&Nimages); + + ALLOCATE (starmap, StarMap, Nimages); + + for (i = 0; i < Nimages; i++) { + starmap[i].Npoints = 0; + starmap[i].points = NULL; + starmap[i].Nx = images[i].NX / NX_MAP; + starmap[i].Ny = images[i].NY / NY_MAP; + ALLOCATE (starmap[i].stars, int, NX_MAP*NY_MAP); + memset (starmap[i].stars, 0, sizeof(int)*NX_MAP*NY_MAP); + } + return (TRUE); +} + +int updateStarMaps(Catalog *catalog) { + + Image *images; + off_t i, N, Nimages; + int xbin, ybin; + struct timeval start, stop; + + gettimeofday (&start, (void *) NULL); + + images = getimages(&Nimages); + + for (i = 0; i < catalog[0].Nmeasure; i++) { + + N = getImageByID(catalog[0].measure[i].imageID); + if (N < 0) continue; + + xbin = catalog[0].measure[i].Xccd / starmap[N].Nx; + ybin = catalog[0].measure[i].Yccd / starmap[N].Ny; + + xbin = MAX(0, MIN(NX_MAP-1, xbin)); + ybin = MAX(0, MIN(NY_MAP-1, ybin)); + + starmap[N].stars[ybin*NX_MAP + xbin] ++; + } + MARKTIME("assign stars to starmap bins: %f sec\n", dtime); + + return (TRUE); +} + +int createStarMapPoints() { + + Image *images; + off_t i, Nimages; + int ix, iy; + + images = getimages(&Nimages); + + for (i = 0; i < Nimages; i++) { + + assert (!starmap[i].points); + + ALLOCATE (starmap[i].points, StarMapPoint, NX_MAP*NY_MAP); + starmap[i].Npoints = 0; + + for (ix = 0; ix < NX_MAP; ix++) { + for (iy = 0; iy < NY_MAP; iy++) { + if (starmap[i].stars[iy*NX_MAP + ix] < 1) continue; + StarMapPoint *point = &starmap[i].points[starmap[i].Npoints]; + starmap[i].Npoints++; + + // set the pixel coordinates + point[0].X = ix * starmap[i].Nx; // XXX fix 0.5 pixel offset + point[0].Y = iy * starmap[i].Nx; // XXX fix 0.5 pixel offset + + // set the transformed coordinates + XY_to_LM (&point[0].L, &point[0].M, point[0].X, point[0].Y, &images[i].coords); + + point[0].dL = 0.0; + point[0].dM = 0.0; + } + } + + if (VERBOSE) fprintf (stderr, "starmap: %d points for image %s\n", starmap[i].Npoints, images[i].name); + } + + return (TRUE); +} + +int checkStarMap(int N) { + + Image *images; + off_t i, Nimages; + double L, M, dLmax, dMmax; + + images = getimages(&Nimages); + + dLmax = dMmax = 0.0; + + for (i = 0; i < starmap[N].Npoints; i++) { + + // set the transformed coordinates + XY_to_LM (&L, &M, starmap[N].points[i].X, starmap[N].points[i].Y, &images[N].coords); + + starmap[N].points[i].dL = starmap[N].points[i].L - L; + starmap[N].points[i].dM = starmap[N].points[i].M - M; + + dLmax = MAX(fabs(starmap[N].points[i].dL), dLmax); + dMmax = MAX(fabs(starmap[N].points[i].dM), dMmax); + } + + if (VERBOSE) fprintf (stderr, "max deviations for %s using %d pts : %f, %f\n", images[N].name, starmap[N].Npoints, dLmax, dMmax); + + if (dLmax > DPOS_MAX) return (FALSE); + if (dMmax > DPOS_MAX) return (FALSE); + return (TRUE); +} + Index: trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateChips.c (revision 28968) +++ trunk/Ohana/src/relastro/src/UpdateChips.c (revision 29001) @@ -31,8 +31,18 @@ // FitChip does iterative, clipped fitting // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n", i, Nraw); - if (!FitChip (raw, ref, Nraw, &image[i].coords)) { - fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + if (!FitChip (raw, ref, Nraw, &image[i])) { + if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); oldCoords = getCoords (i); memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + free (raw); + free (ref); + continue; + } + + if (!checkStarMap (i)) { + if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + oldCoords = getCoords (i); + memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + image[i].flags |= ID_IMAGE_ASTROM_POOR; } Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 28968) +++ trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 29001) @@ -5,4 +5,11 @@ off_t i, Nimage; Image *image; + + int badImage = + ID_IMAGE_ASTROM_NOCAL | + ID_IMAGE_ASTROM_POOR | + ID_IMAGE_ASTROM_FAIL | + ID_IMAGE_ASTROM_SKIP | + ID_IMAGE_ASTROM_FEW; image = getimages (&Nimage); @@ -13,4 +20,7 @@ if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + // skip images that have failed solutions (divergent or otherwise) + if (image[i].flags & badImage) continue; + /* convert measure coordinates to raw entries */ fixImageRaw (catalog, Ncatalog, i); Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 28968) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 29001) @@ -45,5 +45,5 @@ offset = catalog[0].average[i].measureOffset + j; - //filter objects based on user supplied criteria + // filter objects based on user supplied criteria if (!MeasFilterTest(&catalog[0].measure[offset])) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; @@ -51,5 +51,5 @@ } - //filter out outliers + // filter out outliers if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; Index: trunk/Ohana/src/relastro/src/fitpoly.c =================================================================== --- trunk/Ohana/src/relastro/src/fitpoly.c (revision 28968) +++ trunk/Ohana/src/relastro/src/fitpoly.c (revision 29001) @@ -110,4 +110,5 @@ if (fit[0].Npts == 0) { fprintf (stderr, "warning: no valid pts\n"); + return (FALSE); } @@ -140,9 +141,9 @@ } - for (i = 0; i < fit[0].Nelems; i++) { + for (i = 0; FALSE && (i < fit[0].Nelems); i++) { ix = i % fit[0].Nterms; iy = i / fit[0].Nterms; - // fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", - // ix, iy, vector[i][0], ix, iy, vector[i][1]); + fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + ix, iy, vector[i][0], ix, iy, vector[i][1]); } @@ -151,9 +152,9 @@ } - for (i = 0; i < fit[0].Nelems; i++) { + for (i = 0; FALSE && i < fit[0].Nelems; i++) { ix = i % fit[0].Nterms; iy = i / fit[0].Nterms; - // fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", - // ix, iy, vector[i][0], ix, iy, vector[i][1]); + fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + ix, iy, vector[i][0], ix, iy, vector[i][1]); } @@ -274,25 +275,35 @@ /* this should only apply to the polynomial, not the projection terms */ /* compare with psastro supporting code */ -int fit_apply_coords (CoordFit *fit, Coords *coords) { +int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef) { double Xo, Yo, R1, R2; CoordFit *modfit; - /* I have L,M = fit(X,Y). set corresponding terms for coords */ - - // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... - // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... - - if (!CoordsGetCenter (fit, 0.001, &Xo, &Yo)) { - fprintf (stderr, "failed to modify model\n"); - return (FALSE); - } - coords[0].crpix1 = Xo; - coords[0].crpix2 = Yo; - - // resulting fit should have zero constant terms - modfit = CoordsSetCenter (fit, Xo, Yo); - - /* we do not modify crval1,2: these are kept at the default values */ + if (keepRef) { + // adjust crpix1,2 as needed: + // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... + // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... + + if (!CoordsGetCenter (fit, 0.001, &Xo, &Yo)) { + fprintf (stderr, "failed to modify model\n"); + return (FALSE); + } + coords[0].crpix1 = Xo; + coords[0].crpix2 = Yo; + /* we do not modify crval1,2: these are kept at the default values */ + + // resulting fit should have zero constant terms + modfit = CoordsSetCenter (fit, Xo, Yo); + } else { + // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... + // P = cd1*x, Q = cd2*y + // L = pc1_1*P + pc1_2*Q + ... + + /* modify crval1,2: these are kept at the default values */ + + coords[0].crpix1 = 0.0; + coords[0].crpix2 = 0.0; + modfit = fit; + } // set cdelt1, cdelt2 @@ -331,6 +342,10 @@ coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2; } - - fit_free (modfit); + + if (keepRef) { + fit_free (modfit); + } else { + fit_apply (fit, &coords[0].crval1, &coords[0].crval2, coords[0].crpix1, coords[0].crpix2); + } /* keep the order and type from initial values */ Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 28968) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 29001) @@ -100,5 +100,5 @@ initstats (STATMODE); - // IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP; + // IMAGE_BAD = ID_IMAGE_ASTROM_POOR | ID_IMAGE_ASTROM_FEW | ID_IMAGE_ASTROM_SKIP; // STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; // MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_ASTROM | ID_MEAS_SKIP_ASTROM | ID_MEAS_AREA; Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 28968) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 29001) @@ -9,4 +9,6 @@ ALLOCATE (catalog, Catalog, skylist[0].Nregions); + + initStarMaps(); // load data from each region file, only use bright stars @@ -34,4 +36,6 @@ FlagOutliers(&tcatalog); } + + updateStarMaps (&tcatalog); // select only the brighter stars Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 28968) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 29001) @@ -33,4 +33,9 @@ // convert database table to internal structure image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + MARKTIME(" convert image table: %f sec\n", dtime); @@ -52,2 +57,19 @@ return (skylist); } + +int reload_images (FITS_DB *db) { + + Image *image; + off_t Nimage, Nx, i; + VTable *vtable; + + image = getimages (&Nimage); + + vtable = &db[0].vtable; + + gfits_scan (vtable[0].header, "NAXIS1", OFF_T_FMT, 1, &Nx); + for (i = 0; i < Nimage; i++) { + memcpy (vtable[0].buffer[i], &image[i], Nx); + } + return (TRUE); +} Index: trunk/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- trunk/Ohana/src/relastro/src/mkpolyterm.c (revision 28968) +++ trunk/Ohana/src/relastro/src/mkpolyterm.c (revision 29001) @@ -66,4 +66,5 @@ Xo -= beta[0][0]; Yo -= beta[1][0]; + dPos = hypot(beta[0][0], beta[1][0]); if (i == 0) { Index: trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro.c (revision 28968) +++ trunk/Ohana/src/relastro/src/relastro.c (revision 29001) @@ -1,9 +1,3 @@ # include "relastro.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } int main (int argc, char **argv) { @@ -49,4 +43,6 @@ /* load catalog data from region files : subselect high-quality measurements */ + // XXX pass in the image table + // XXX who carries the image grid? catalog = load_catalogs (skylist, &Ncatalog, TRUE); MARKTIME("load catalog data: %f sec\n", dtime); @@ -63,4 +59,7 @@ // plot_mosaic_fields (catalog); } + + // set test points based on the starmap + createStarMapPoints(); /* major modes */ @@ -98,9 +97,11 @@ freeImageBins (1); + reload_images (&db); + // iterate over catalogs to make detection coordinates consistant UpdateObjectOffsets (skylist); // iterate over catalogs to make detection coordinates consistant - FixProblemImages (skylist); + // FixProblemImages (skylist); // save the updated image parameters Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 28968) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 29001) @@ -28,4 +28,9 @@ double *RmaxSky; off_t *index; + + int badImage = + ID_IMAGE_ASTROM_POOR | + ID_IMAGE_ASTROM_FAIL | + ID_IMAGE_ASTROM_FEW; if (skylist[0].Nregions < 1) { @@ -171,6 +176,6 @@ image[nimage] = timage[i]; /* always allow 'few' images to succeed, if possible */ - if (image[nimage].flags & ID_IMAGE_FEW) { - image[nimage].flags &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR); + if (image[nimage].flags & ID_IMAGE_ASTROM_FEW) { + image[nimage].flags &= ~ID_IMAGE_ASTROM_FEW; } if (RESET) { @@ -178,5 +183,5 @@ assignMcal (&image[nimage], (double *) NULL, -1); image[nimage].dMcal = NAN; - image[nimage].flags &= ~ID_IMAGE_POOR; + image[nimage].flags &= ~badImage; } line_number[nimage] = i; Index: trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relphot/src/ImageOps.c (revision 28968) +++ trunk/Ohana/src/relphot/src/ImageOps.c (revision 29001) @@ -323,5 +323,5 @@ /* on PoorImages run, skip good images */ if (PoorImages) { - bad = image[i].flags & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP); + bad = image[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); if (!bad) continue; } @@ -369,8 +369,8 @@ mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]); if (mark) { - image[i].flags |= ID_IMAGE_FEW; + image[i].flags |= ID_IMAGE_PHOTOM_FEW; Nfew ++; } else { - image[i].flags &= ~ID_IMAGE_FEW; + image[i].flags &= ~ID_IMAGE_PHOTOM_FEW; } } @@ -387,5 +387,5 @@ if (PoorImages) { - IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP; + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; @@ -429,11 +429,11 @@ for (i = 0; i < Nimage; i++) { mark = FALSE; - image[i].flags &= ~ID_IMAGE_POOR; + image[i].flags &= ~ID_IMAGE_PHOTOM_POOR; mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].Mcal - MedOffset) > MaxOffset); if (mark) { Nmark ++; - image[i].flags |= ID_IMAGE_POOR; + image[i].flags |= ID_IMAGE_PHOTOM_POOR; } else { - image[i].flags &= ~ID_IMAGE_POOR; + image[i].flags &= ~ID_IMAGE_PHOTOM_POOR; } } Index: trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- trunk/Ohana/src/relphot/src/MosaicOps.c (revision 28968) +++ trunk/Ohana/src/relphot/src/MosaicOps.c (revision 29001) @@ -189,6 +189,6 @@ image[im].dMcal = mosaic[i].dMcal; image[im].Xm = mosaic[i].Xm; - image[im].flags |= (mosaic[i].flags & ID_IMAGE_FEW); - image[im].flags |= (mosaic[i].flags & ID_IMAGE_POOR); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW); + image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR); } } @@ -351,5 +351,5 @@ /* on PoorImages run, skip good images */ if (PoorImages) { - bad = mosaic[i].flags & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP); + bad = mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP); if (!bad) continue; } @@ -401,8 +401,8 @@ if (mark) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), ("OFF_T_FMT" < %d) || ("OFF_T_FMT" < %f*"OFF_T_FMT")\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]); - mosaic[i].flags |= ID_IMAGE_FEW; + mosaic[i].flags |= ID_IMAGE_PHOTOM_FEW; Nfew ++; } else { - mosaic[i].flags &= ~ID_IMAGE_FEW; + mosaic[i].flags &= ~ID_IMAGE_PHOTOM_FEW; } } @@ -421,5 +421,5 @@ if (PoorImages) { - IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP; + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; @@ -593,7 +593,7 @@ if (mark) { Nmark ++; - mosaic[i].flags |= ID_IMAGE_POOR; + mosaic[i].flags |= ID_IMAGE_PHOTOM_POOR; } else { - mosaic[i].flags &= ~ID_IMAGE_POOR; + mosaic[i].flags &= ~ID_IMAGE_PHOTOM_POOR; } } Index: trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- trunk/Ohana/src/relphot/src/StarOps.c (revision 28968) +++ trunk/Ohana/src/relphot/src/StarOps.c (revision 29001) @@ -188,4 +188,7 @@ StatType stats; PhotCode *code; + DVOAverageFlags flagBits; + + flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT; // pre-allocate a list for stats purposes @@ -281,15 +284,18 @@ } + // we attempt to set a few flags here; reset those bits before trying: + catalog[i].average[j].flags &= ~flagBits; + if (nEXT && (nEXT > nPSF)) { - catalog[i].average[j].flags |= 0x00010000; + catalog[i].average[j].flags |= ID_OBJ_EXT; } if (goodPS1) { - catalog[i].average[j].flags |= 0x00020000; + catalog[i].average[j].flags |= ID_OBJ_GOOD; } if (Galaxy2MASS) { - catalog[i].average[j].flags |= 0x00040000; + catalog[i].average[j].flags |= ID_OBJ_EXT_ALT; } if (good2MASS) { - catalog[i].average[j].flags |= 0x00080000; + catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT; } } @@ -302,5 +308,5 @@ } -/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_NOCAL */ +/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */ int setMcalOutput (Catalog *catalog, int Ncatalog) { @@ -310,5 +316,5 @@ MEAS_BAD = ID_MEAS_NOCAL; - IMAGE_BAD = ID_IMAGE_NOCAL; + IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; for (i = 0; i < Ncatalog; i++) { @@ -470,5 +476,5 @@ /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */ image_bad = IMAGE_BAD; - IMAGE_BAD = ID_IMAGE_NOCAL; + IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; m = catalog[i].average[j].measureOffset; N = 0; Index: trunk/Ohana/src/relphot/src/initialize.c =================================================================== --- trunk/Ohana/src/relphot/src/initialize.c (revision 28968) +++ trunk/Ohana/src/relphot/src/initialize.c (revision 29001) @@ -27,5 +27,5 @@ initstats (STATMODE); - IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP; + IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP; STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA; Index: trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- trunk/Ohana/src/relphot/src/load_images.c (revision 28968) +++ trunk/Ohana/src/relphot/src/load_images.c (revision 29001) @@ -25,4 +25,8 @@ // convert database table to internal structure image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } // select the images which overlap the selected sky regions Index: trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- trunk/Ohana/src/relphot/src/select_images.c (revision 28968) +++ trunk/Ohana/src/relphot/src/select_images.c (revision 29001) @@ -157,11 +157,11 @@ image[nimage] = timage[i]; /* always allow 'few' images to succeed, if possible */ - if (image[nimage].flags & ID_IMAGE_FEW) { - image[nimage].flags &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR); + if (image[nimage].flags & ID_IMAGE_PHOTOM_FEW) { + image[nimage].flags &= ~(ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR); } if (RESET) { assignMcal (&image[nimage], (double *) NULL, -1); image[nimage].dMcal = NAN; - image[nimage].flags &= ~ID_IMAGE_POOR; + image[nimage].flags &= ~ID_IMAGE_PHOTOM_POOR; } line_number[nimage] = i; Index: trunk/Ohana/src/relphot/src/setMrelFinal.c =================================================================== --- trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 28968) +++ trunk/Ohana/src/relphot/src/setMrelFinal.c (revision 29001) @@ -75,5 +75,5 @@ /* allow measures from images marked POOR and FEW */ - if (pass >= 3) IMAGE_BAD = ID_IMAGE_NOCAL; + if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL; /* allow measures marked as outliers (POOR) and off image region (AREA) */ Index: trunk/Ohana/src/uniphot/src/find_image_sgroups.c =================================================================== --- trunk/Ohana/src/uniphot/src/find_image_sgroups.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/find_image_sgroups.c (revision 29001) @@ -13,4 +13,8 @@ imlink = *Imlink; image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } coords.crpix1 = coords.crpix2 = 0.0; @@ -28,5 +32,5 @@ for (i = 0; i < Nimage; i++) { if (imlink[i].sgroup != NULL) continue; - if (image[i].flags & ID_IMAGE_NOCAL) continue; + if (image[i].flags & ID_IMAGE_PHOTOM_NOCAL) continue; /* new sgroup, set ref coords */ @@ -53,5 +57,5 @@ for (j = 0; j < Nimage; j++) { - if (image[j].flags & ID_IMAGE_NOCAL) continue; + if (image[j].flags & ID_IMAGE_PHOTOM_NOCAL) continue; if (imlink[j].sgroup != NULL) continue; Index: trunk/Ohana/src/uniphot/src/find_image_tgroups.c =================================================================== --- trunk/Ohana/src/uniphot/src/find_image_tgroups.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/find_image_tgroups.c (revision 29001) @@ -12,4 +12,9 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + ALLOCATE (imlink, ImageLink, Nimage); @@ -64,5 +69,5 @@ if (image[j].tzero < tmin[i]) continue; if (image[j].tzero > tmax[i]) continue; - if (image[j].flags & ID_IMAGE_NOCAL) continue; + if (image[j].flags & ID_IMAGE_PHOTOM_NOCAL) continue; group[i].image[Nentry] = &image[j]; Index: trunk/Ohana/src/uniphot/src/initialize.c =================================================================== --- trunk/Ohana/src/uniphot/src/initialize.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/initialize.c (revision 29001) @@ -16,5 +16,5 @@ } - IMAGE_BAD = ID_IMAGE_NOCAL | ID_IMAGE_POOR | ID_IMAGE_SKIP | ID_IMAGE_FEW; + IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW; initstats (STATMODE); Index: trunk/Ohana/src/uniphot/src/load_images.c =================================================================== --- trunk/Ohana/src/uniphot/src/load_images.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/load_images.c (revision 29001) @@ -22,4 +22,8 @@ /* use a vtable to keep the images to be calibrated */ image = gfits_table_get_Image (&db[0].ftable, Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } fprintf (stderr, "loaded "OFF_T_FMT" images\n", *Nimage); Index: trunk/Ohana/src/uniphot/src/match_zpts_to_images.c =================================================================== --- trunk/Ohana/src/uniphot/src/match_zpts_to_images.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/match_zpts_to_images.c (revision 29001) @@ -71,5 +71,5 @@ image[Ni].Mcal = SCALE*code[0].C - zpts[Nz].zpt; image[Ni].dMcal = zpts[Nz].zpt_err; - image[Ni].flags &= ~ID_IMAGE_NOCAL; // clear the NOCAL flag + image[Ni].flags &= ~ID_IMAGE_PHOTOM_NOCAL; // clear the NOCAL flag zpts[Nz].found = TRUE; NImatch ++; Index: trunk/Ohana/src/uniphot/src/subset_images.c =================================================================== --- trunk/Ohana/src/uniphot/src/subset_images.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/subset_images.c (revision 29001) @@ -9,4 +9,8 @@ /* use a vtable to keep the images to be calibrated */ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } Nkeep = 0; @@ -16,5 +20,5 @@ for (i = 0; i < Nimage; i++) { - image[i].flags |= ID_IMAGE_NOCAL; + image[i].flags |= ID_IMAGE_PHOTOM_NOCAL; /* select images by photcode */ @@ -27,5 +31,5 @@ if (image[i].tzero > TSTOP) continue; } - image[i].flags &= ~ID_IMAGE_NOCAL; + image[i].flags &= ~ID_IMAGE_PHOTOM_NOCAL; keep[Nkeep] = i; Nkeep ++; Index: trunk/Ohana/src/uniphot/src/update.c =================================================================== --- trunk/Ohana/src/uniphot/src/update.c (revision 28968) +++ trunk/Ohana/src/uniphot/src/update.c (revision 29001) @@ -14,8 +14,12 @@ image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } /* clear the NOCAL flags */ for (i = 0; i < Nimage; i++) { - image[i].flags &= ~ID_IMAGE_NOCAL; + image[i].flags &= ~ID_IMAGE_PHOTOM_NOCAL; }