Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 37115) +++ /trunk/Ohana/src/addstar/Makefile (revision 37116) @@ -78,8 +78,8 @@ $(SRC)/parse_time.$(ARCH).o \ $(SRC)/replace_match.$(ARCH).o \ -$(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/resort_catalogs.$(ARCH).o \ $(SRC)/resort_threaded.$(ARCH).o \ $(SRC)/resort_unthreaded.$(ARCH).o \ -$(SRC)/resort_unthreaded_catalogs.$(ARCH).o \ +$(SRC)/resort_catalog.$(ARCH).o \ $(SRC)/StarOps.$(ARCH).o \ $(SRC)/ReadStarsFITS.$(ARCH).o \ @@ -105,8 +105,8 @@ $(SRC)/addstar_client.$(ARCH).o \ $(SRC)/args_parallel_client.$(ARCH).o \ -$(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/resort_catalogs.$(ARCH).o \ $(SRC)/resort_threaded.$(ARCH).o \ $(SRC)/resort_unthreaded.$(ARCH).o \ -$(SRC)/resort_unthreaded_catalogs.$(ARCH).o +$(SRC)/resort_catalog.$(ARCH).o ADDSTARD = \ Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 37115) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 37116) @@ -237,9 +237,11 @@ Stars *rd_gsc PROTO((char *filename, unsigned int *nstars)); int replace_match PROTO((Average *average, Measure *measure, Stars *star)); -int resort_threaded PROTO((AddstarClientOptions *options, SkyTable *sky)); -int resort_unthreaded PROTO((AddstarClientOptions *options, SkyTable *sky, int hostID, char *hostpath)); -int resort_unthreaded_catalogs PROTO((AddstarClientOptions *options, SkyList *skylist, int hostID, char *hostpath)); +int resort_catalogs PROTO((AddstarClientOptions *options, SkyTable *sky)); +int resort_catalogs_parallel PROTO((AddstarClientOptions *options, SkyList *sky)); +int resort_threaded PROTO((SkyList *skylist, int ForceSort)); +int resort_unthreaded PROTO((SkyList *skylist, int ForceSort)); void resort_catalog PROTO((Catalog *catalog)); void resort_catalog_old PROTO((Catalog *catalog)); + Stars *ReadStarsFITS PROTO((FILE *f, Header *header, Header *in_theader, unsigned int *nstars)); Stars *ReadStarsTEXT PROTO((FILE *f, unsigned int *nstars)); @@ -253,7 +255,7 @@ off_t *init_lensing_links PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing)); off_t *init_missing_links PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing)); -off_t add_meas_link PROTO((Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE)); -off_t add_miss_link PROTO((Average *average, off_t *next_miss, off_t Nmissing)); -off_t add_lens_link PROTO((Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING)); +int add_meas_link PROTO((Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE)); +int add_miss_link PROTO((Average *average, off_t *next_miss, off_t Nmissing)); +int add_lens_link PROTO((Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING)); off_t *build_measure_links PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure)); off_t *build_lensing_links PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing)); @@ -316,4 +318,5 @@ void *memrchr(const void *s, int c, size_t n); int addstar_create_ID (); +int strextend (char *input, char *format,...); /** Index: /trunk/Ohana/src/addstar/src/BoundaryTreeIO.c =================================================================== --- /trunk/Ohana/src/addstar/src/BoundaryTreeIO.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/BoundaryTreeIO.c (revision 37116) @@ -1,3 +1,5 @@ # include "addstar.h" + +/** this file is deprecated (functions moved to libdvo) ***/ # define GET_COLUMN_NEW(OUT,NAME,TYPE) \ Index: /trunk/Ohana/src/addstar/src/addstar.c =================================================================== --- /trunk/Ohana/src/addstar/src/addstar.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/addstar.c (revision 37116) @@ -40,9 +40,5 @@ if (options.mode == ADDSTAR_MODE_RESORT) { - if (NTHREADS == 0) { - resort_unthreaded (&options, sky, 0, NULL); - } else { - resort_threaded (&options, sky); - } + resort_catalogs (&options, sky); exit (0); } @@ -71,5 +67,4 @@ skylist = SkyListForStars (sky, -1, stars, Nstars); break; - case ADDSTAR_MODE_RESORT: case ADDSTAR_MODE_REFCAT: skylist = SkyListByPatch (sky, -1, &UserPatch); @@ -161,21 +156,4 @@ if (Nsubset) free (subset); break; - case ADDSTAR_MODE_RESORT: - if (options.nosort == 3) catalog.sorted = FALSE; - - // no need to resort empty catalogs - if (catalog.Naves_disk == 0) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - continue; - } - - if (OLD_RESORT) { - resort_catalog_old (&catalog); - } else { - resort_catalog (&catalog); - } - Nsubset = 1; - break; } /* report total updated values */ Index: /trunk/Ohana/src/addstar/src/addstar_client.c =================================================================== --- /trunk/Ohana/src/addstar/src/addstar_client.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/addstar_client.c (revision 37116) @@ -25,5 +25,5 @@ if (options.mode == ADDSTAR_MODE_RESORT) { - resort_unthreaded (&options, sky, HOST_ID, HOSTDIR); + resort_catalogs (&options, sky); exit (0); } Index: /trunk/Ohana/src/addstar/src/args.c =================================================================== --- /trunk/Ohana/src/addstar/src/args.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/args.c (revision 37116) @@ -34,4 +34,10 @@ remove_argument (N, &argc, argv); } + OLD_RESORT = FALSE; + if ((N = get_argument (argc, argv, "-old-resort"))) { + remove_argument (N, &argc, argv); + OLD_RESORT = TRUE; + } + if ((N = get_argument (argc, argv, "-create-id"))) { options.mode = ADDSTAR_MODE_CREATE_ID; @@ -350,4 +356,6 @@ // XXX for the moment, make this selection manual. it needs to be automatic // based on the state of the SkyTable + HOST_ID = 0; + HOSTDIR = NULL; PARALLEL = FALSE; if ((N = get_argument (argc, argv, "-parallel"))) { @@ -466,4 +474,17 @@ } +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + /** addstar modes: Index: /trunk/Ohana/src/addstar/src/args_parallel_client.c =================================================================== --- /trunk/Ohana/src/addstar/src/args_parallel_client.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/args_parallel_client.c (revision 37116) @@ -87,4 +87,10 @@ } + OLD_RESORT = FALSE; + if ((N = get_argument (argc, argv, "-old-resort"))) { + remove_argument (N, &argc, argv); + OLD_RESORT = TRUE; + } + /* extra error messages */ VERBOSE = FALSE; @@ -118,2 +124,16 @@ exit (2); } + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + Index: /trunk/Ohana/src/addstar/src/build_links.c =================================================================== --- /trunk/Ohana/src/addstar/src/build_links.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/build_links.c (revision 37116) @@ -102,5 +102,5 @@ /* average[].measureOffset, average[].Nmeasure are valid within an addstar run */ -off_t add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) { +int add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) { off_t k, m; @@ -133,7 +133,4 @@ return (TRUE); } - -/* Missing does not carry enough information to reconstruct the links - we must always save the missing table, if it exists */ Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) { @@ -160,4 +157,6 @@ } +/*******************************************************************************************/ + /* build the initial links assuming the table is sorted */ off_t *init_missing_links (Average *average, off_t Naverage, Missing *missing, off_t Nmissing) { @@ -186,5 +185,5 @@ } -off_t add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) { +int add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) { off_t k, m; @@ -205,4 +204,7 @@ } +/* Missing does not carry enough information to reconstruct the links + we must always save the missing table, if it exists */ + Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss) { @@ -307,5 +309,5 @@ /* average[].lensingOffset, average[].Nlensing are valid within an addstar run */ -off_t add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) { +int add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) { off_t k, m; Index: /trunk/Ohana/src/addstar/src/findskycell.c =================================================================== --- /trunk/Ohana/src/addstar/src/findskycell.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/findskycell.c (revision 37116) @@ -27,7 +27,7 @@ int apply_tree (char *treefile, char *datafile); -float SCALE = 1.0; -int NX_SUB = 1; -int NY_SUB = 1; +float SCALE = 1.0; +float NX_SUB = 1.0; +float NY_SUB = 1.0; double R_MIN = 0.0; @@ -217,6 +217,6 @@ ALLOCATE (tree.Xo, double *, tree.Nzone); ALLOCATE (tree.Yo, double *, tree.Nzone); - ALLOCATE (tree.dX, int *, tree.Nzone); - ALLOCATE (tree.dY, int *, tree.Nzone); + ALLOCATE (tree.dX, float *, tree.Nzone); + ALLOCATE (tree.dY, float *, tree.Nzone); ALLOCATE (tree.cell, int *, tree.Nzone); ALLOCATE (tree.name, char **, tree.Nzone); @@ -251,6 +251,6 @@ ALLOCATE (tree.Xo[zone], double, tree.NBAND[zone]); ALLOCATE (tree.Yo[zone], double, tree.NBAND[zone]); - ALLOCATE (tree.dX[zone], int, tree.NBAND[zone]); - ALLOCATE (tree.dY[zone], int, tree.NBAND[zone]); + ALLOCATE (tree.dX[zone], float, tree.NBAND[zone]); + ALLOCATE (tree.dY[zone], float, tree.NBAND[zone]); ALLOCATE (tree.cell[zone], int, tree.NBAND[zone]); ALLOCATE (tree.name[zone], char *, tree.NBAND[zone]); @@ -260,6 +260,6 @@ tree.Xo[zone][band] = NAN; tree.Yo[zone][band] = NAN; - tree.dX[zone][band] = -1; - tree.dY[zone][band] = -1; + tree.dX[zone][band] = NAN; + tree.dY[zone][band] = NAN; tree.cell[zone][band] = -1; ALLOCATE (tree.name[zone][band], char, BOUNDARY_TREE_NAME_LENGTH); @@ -286,6 +286,6 @@ REALLOCATE (tree.Xo[zone], double, tree.NBAND[zone]); REALLOCATE (tree.Yo[zone], double, tree.NBAND[zone]); - REALLOCATE (tree.dX[zone], int, tree.NBAND[zone]); - REALLOCATE (tree.dY[zone], int, tree.NBAND[zone]); + REALLOCATE (tree.dX[zone], float, tree.NBAND[zone]); + REALLOCATE (tree.dY[zone], float, tree.NBAND[zone]); REALLOCATE (tree.cell[zone], int, tree.NBAND[zone]); REALLOCATE (tree.name[zone], char *, tree.NBAND[zone]); @@ -295,6 +295,6 @@ tree.Xo[zone][band] = NAN; tree.Yo[zone][band] = NAN; - tree.dX[zone][band] = -1; - tree.dY[zone][band] = -1; + tree.dX[zone][band] = NAN; + tree.dY[zone][band] = NAN; tree.cell[zone][j] = -1; ALLOCATE (tree.name[zone][j], char, BOUNDARY_TREE_NAME_LENGTH); Index: /trunk/Ohana/src/addstar/src/resort_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_catalog.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/resort_catalog.c (revision 37116) @@ -7,6 +7,6 @@ void resort_catalog_old (Catalog *catalog) { - off_t *next_meas; - off_t Naves, Nmeas; + off_t *next_meas, *next_lens; + off_t Naves, Nmeas, Nlens; if (catalog[0].sorted == TRUE) return; @@ -16,11 +16,14 @@ /* internal counters */ Nmeas = catalog[0].Nmeasure; + Nlens = catalog[0].Nlensing; Naves = catalog[0].Naverage; /* set up pointers for linked list of measure, missing */ next_meas = build_measure_links (catalog[0].average, Naves, catalog[0].measure, Nmeas); + next_lens = build_lensing_links (catalog[0].average, Naves, catalog[0].lensing, Nlens); catalog[0].sorted = TRUE; catalog[0].measure = sort_measure (catalog[0].average, Naves, catalog[0].measure, Nmeas, next_meas); + catalog[0].lensing = sort_lensing (catalog[0].average, Naves, catalog[0].lensing, Nlens, next_lens); MARKTIME (" match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves); Index: /trunk/Ohana/src/addstar/src/resort_catalogs.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_catalogs.c (revision 37116) +++ /trunk/Ohana/src/addstar/src/resort_catalogs.c (revision 37116) @@ -0,0 +1,120 @@ +# include "addstar.h" + +// pass in options and skylist? +int resort_catalogs (AddstarClientOptions *options, SkyTable *sky) { + + INITTIME; + + if (options->only_images) { + fprintf (stderr, "-image (only images) makes no sense with -resort\n"); + exit (2); + } + if (options->only_match) { + fprintf (stderr, "-only-match makes no sense with -resort\n"); + exit (2); + } + if (options->nosort == 1) { + fprintf (stderr, "-nosort makes no sense with -resort\n"); + exit (2); + } + if (options->update) { + fprintf (stderr, "-update makes no sense with -resort\n"); + exit (2); + } + if (options->calibrate) { + fprintf (stderr, "-cal (calibrate) makes no sense with -resort\n"); + exit (2); + } + + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + // in these cases, limit the sky catalogs to an existing subset + if (options->existing_regions) { + SkyList *tmp; + tmp = SkyListExistingSubset (skylist, CATDIR); + SkyListFree (skylist); + skylist = tmp; + } + if (VERBOSE) fprintf (stderr, "writing to "OFF_T_FMT" regions\n", skylist[0].Nregions); + + if (PARALLEL && !HOST_ID) { + resort_catalogs_parallel (options, skylist); + return TRUE; + } + + int ForceSort = (options->nosort == 3); + + if (NTHREADS == 0) { + resort_unthreaded (skylist, ForceSort); + } else { + resort_threaded (skylist, ForceSort); + } + + MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime); + return TRUE; +} + +// CATDIR is supplied globally +# define DEBUG 1 +int resort_catalogs_parallel (AddstarClientOptions *options, SkyList *sky) { + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + char command[1024]; + snprintf (command, 1024, "addstar_client -resort"); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-hostID %d", table->hosts[i].hostID); + strextend (command, "-hostdir %s", table->hosts[i].pathname); + strextend (command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + if (VERBOSE) strextend (command, "%s -v"); + if (NTHREADS) strextend (command, "%s -threads %d", NTHREADS); + if (OLD_RESORT) strextend (command, "%s -old-resort"); + if (options->nosort == 3) strextend (command, "%s -force-sort"); + if (options->existing_regions) strextend (command, "%s -existing-regions"); + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running addstar_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the addstar_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: /trunk/Ohana/src/addstar/src/resort_threaded.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_threaded.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/resort_threaded.c (revision 37116) @@ -13,5 +13,5 @@ int iThread; int state; - int nosort; + int forcesort; char *filename; SkyRegion *region; @@ -37,6 +37,10 @@ } + // chose the catalog file (local or remote?) + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, threadData->region->name); + catalog.filename = HOST_ID ? hostfile : threadData->filename; + // set the parameters which guide catalog open/load/create - catalog.filename = threadData->filename; 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 @@ -64,5 +68,5 @@ // this is an overloaded value to mean 'force sort' - if (threadData->nosort == 3) catalog.sorted = FALSE; + if (threadData->forcesort) catalog.sorted = FALSE; if (OLD_RESORT) { @@ -88,5 +92,5 @@ } -int resort_threaded (AddstarClientOptions *options, SkyTable *sky) { +int resort_threaded (SkyList *skylist, int ForceSort) { int j, launched, done, Nwait, Nrun; @@ -97,40 +101,4 @@ ThreadData *threadData; pthread_t *threads; - - double dtime; - struct timeval start, stop; - gettimeofday (&start, NULL); - - SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - - // in these cases, limit the sky catalogs to an existing subset - if (options->existing_regions) { - SkyList *tmp; - tmp = SkyListExistingSubset (skylist, CATDIR); - SkyListFree (skylist); - skylist = tmp; - } - if (VERBOSE) fprintf (stderr, "writing to "OFF_T_FMT" regions\n", skylist[0].Nregions); - - if (options->only_images) { - fprintf (stderr, "-image (only images) makes no sense with -resort\n"); - exit (2); - } - if (options->only_match) { - fprintf (stderr, "-only-match makes no sense with -resort\n"); - exit (2); - } - if (options->nosort == 1) { - fprintf (stderr, "-nosort makes no sense with -resort\n"); - exit (2); - } - if (options->update) { - fprintf (stderr, "-update makes no sense with -resort\n"); - exit (2); - } - if (options->calibrate) { - fprintf (stderr, "-cal (calibrate) makes no sense with -resort\n"); - exit (2); - } ALLOCATE(threads, pthread_t, NTHREADS); @@ -148,4 +116,7 @@ Naverage = Nmeasure = 0; for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; launched = FALSE; @@ -166,5 +137,5 @@ if (threadData[j].state == TS_FAIL) goto failure; if (threadData[j].state != TS_WAIT) continue; - threadData[j].nosort = options->nosort; + threadData[j].forcesort = ForceSort; threadData[j].filename = skylist[0].filename[i]; threadData[j].region = skylist[0].regions[i]; @@ -203,8 +174,4 @@ } - gettimeofday (&stop, NULL); - dtime = DTIME (stop, start); - fprintf (stderr, "SUCCESS: elapsed time %9.4f sec for, "OFF_T_FMT" average, "OFF_T_FMT" measure\n", dtime, Naverage, Nmeasure); - // for (i = 0; i < NTHREADS; i++) { // pthread_cancel(threads[i]); Index: /trunk/Ohana/src/addstar/src/resort_unthreaded.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_unthreaded.c (revision 37115) +++ /trunk/Ohana/src/addstar/src/resort_unthreaded.c (revision 37116) @@ -1,50 +1,64 @@ # include "addstar.h" -int resort_unthreaded (AddstarClientOptions *options, SkyTable *sky, int hostID, char *hostpath) { +// pass in options and skylist? +int resort_unthreaded (SkyList *skylist, int ForceSort) { - double dtime; - struct timeval start, stop; + off_t i; + off_t Naverage, Nmeasure; + Catalog catalog; - gettimeofday (&start, NULL); + /* match stars to existing catalog data (or otherwise manipulate catalog data) */ + Naverage = Nmeasure = 0; + for (i = 0; i < skylist[0].Nregions; i++) { - SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - - // in these cases, limit the sky catalogs to an existing subset - if (options->existing_regions) { - SkyList *tmp; - tmp = SkyListExistingSubset (skylist, CATDIR); - SkyListFree (skylist); - skylist = tmp; + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], HOST_ID)) continue; + + // chose the catalog file (local or remote?) + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i]; + + // set the parameters which guide catalog open/load/create + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_LENSING; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + // an error exit status here is a significant error (disk I/O or file access) + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + // Naves_disk == 0 implies an empty catalog file, skip empty catalogs + if (catalog.Naves_disk == 0) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + // this is an overloaded value to mean 'force sort' + if (ForceSort) catalog.sorted = FALSE; + + if (OLD_RESORT) { + resort_catalog_old (&catalog); + } else { + resort_catalog (&catalog); + } + + // report total updated values + Naverage += catalog.Naverage; + Nmeasure += catalog.Nmeasure; + + // write out catalog, if appropriate + SetProtect (TRUE); + dvo_catalog_save (&catalog, VERBOSE); + SetProtect (FALSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); } - if (VERBOSE) fprintf (stderr, "writing to "OFF_T_FMT" regions\n", skylist[0].Nregions); - - if (options->only_images) { - fprintf (stderr, "-image (only images) makes no sense with -resort\n"); - exit (2); - } - if (options->only_match) { - fprintf (stderr, "-only-match makes no sense with -resort\n"); - exit (2); - } - if (options->nosort == 1) { - fprintf (stderr, "-nosort makes no sense with -resort\n"); - exit (2); - } - if (options->update) { - fprintf (stderr, "-update makes no sense with -resort\n"); - exit (2); - } - if (options->calibrate) { - fprintf (stderr, "-cal (calibrate) makes no sense with -resort\n"); - exit (2); - } - - resort_unthreaded_catalogs (options, skylist, hostID, hostpath); - - gettimeofday (&stop, NULL); - dtime = DTIME (stop, start); - fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); - return TRUE; } + Index: unk/Ohana/src/addstar/src/resort_unthreaded_catalogs.c =================================================================== --- /trunk/Ohana/src/addstar/src/resort_unthreaded_catalogs.c (revision 37115) +++ (revision ) @@ -1,132 +1,0 @@ -# include "addstar.h" - -int resort_unthreaded_parallel (AddstarClientOptions *options, SkyList *sky); - -// pass in options and skylist? -int resort_unthreaded_catalogs (AddstarClientOptions *options, SkyList *skylist, int hostID, char *hostpath) { - - off_t i; - off_t Naverage, Nmeasure; - Catalog catalog; - - if (PARALLEL && !hostID) { - resort_unthreaded_parallel (options, skylist); - return TRUE; - } - - /* match stars to existing catalog data (or otherwise manipulate catalog data) */ - Naverage = Nmeasure = 0; - for (i = 0; i < skylist[0].Nregions; i++) { - - // does this host ID match the desired location for the table? - if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; - - // chose the catalog file (local or remote?) - char hostfile[1024]; - snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); - catalog.filename = hostID ? hostfile : skylist[0].filename[i]; - - // set the parameters which guide catalog open/load/create - catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data - catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_LENSING; - catalog.Nsecfilt = GetPhotcodeNsecfilt (); - - // an error exit status here is a significant error (disk I/O or file access) - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { - fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); - exit (2); - } - - // Naves_disk == 0 implies an empty catalog file, skip empty catalogs - if (catalog.Naves_disk == 0) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - continue; - } - - // this is an overloaded value to mean 'force sort' - if (options->nosort == 3) catalog.sorted = FALSE; - - if (OLD_RESORT) { - resort_catalog_old (&catalog); - } else { - resort_catalog (&catalog); - } - - // report total updated values - Naverage += catalog.Naverage; - Nmeasure += catalog.Nmeasure; - - // write out catalog, if appropriate - SetProtect (TRUE); - dvo_catalog_save (&catalog, VERBOSE); - SetProtect (FALSE); - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - } - return TRUE; -} - -// CATDIR is supplied globally -# define DEBUG 1 -int resort_unthreaded_parallel (AddstarClientOptions *options, SkyList *sky) { - - // load the list of hosts - HostTable *table = HostTableLoad (CATDIR, sky->hosts); - if (!table) { - fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); - exit (1); - } - - int i; - for (i = 0; i < table->Nhosts; i++) { - - // ensure that the paths are absolute path names - char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); - free (table->hosts[i].pathname); - table->hosts[i].pathname = tmppath; - - char command[1024]; - snprintf (command, 1024, "addstar_client -resort -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", - table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, - UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); - - char tmpline[1024]; - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (options->nosort == 3) { snprintf (tmpline, 1024, "%s -force-sort", command); strcpy (command, tmpline); } - if (options->existing_regions) { snprintf (tmpline, 1024, "%s -existing-regions", command); strcpy (command, tmpline); } - - fprintf (stderr, "command: %s\n", command); - - if (PARALLEL_MANUAL) continue; - - if (PARALLEL_SERIAL) { - int status = system (command); - if (status) { - fprintf (stderr, "ERROR running addstar_client\n"); - exit (2); - } - } else { - // launch the job on the remote machine (no handshake) - int errorInfo = 0; - int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); - if (!pid) { - if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); - exit (1); - } - table->hosts[i].pid = pid; // save for future reference - } - } - - if (PARALLEL_MANUAL) { - fprintf (stderr, "run the addstar_client commands above. when these are done, hit return\n"); - getchar(); - } - if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { - HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); - } - - return TRUE; -} - Index: /trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- /trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 37115) +++ /trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 37116) @@ -129,10 +129,16 @@ int merge_catalogs_old PROTO((SkyRegion *region, Catalog *output, Catalog *input, double RADIUS, int *secflitMap)); +off_t *build_measure_links PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure)); off_t *init_measure_links PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure)); +int add_meas_link PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE)); +Measure *sort_measure PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next)); + +off_t *build_lensing_links PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing)); +off_t *init_lensing_links PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing)); +int add_lens_link PROTO((Average *average, off_t *next, off_t Nlensing, off_t NLENSING)); +Lensing *sort_lensing PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next)); + off_t *init_missing_links PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing)); -off_t *build_measure_links PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure)); -int add_meas_link PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE)); int add_miss_link PROTO((Average *average, off_t *next, off_t Nmissing)); -Measure *sort_measure PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next)); Missing *sort_missing PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss)); Index: /trunk/Ohana/src/dvomerge/src/LoadCatalog.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/LoadCatalog.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/LoadCatalog.c (revision 37116) @@ -8,5 +8,5 @@ // always load all of the data (if any exists) - catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING | LOAD_LENSOBJ; catalog[0].catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data Index: /trunk/Ohana/src/dvomerge/src/build_links.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/build_links.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/build_links.c (revision 37116) @@ -1,3 +1,26 @@ # include "dvomerge.h" + +/* + +There are two modes for the measure table: sorted and unsorted. + +In sorted mode, all measures associated with a given average are in a single block. The block +is pointed to by average->measureOffset and the range is average->Nmeasure + +In unsorted mode, it is not possible to go directly from average to measure without scanning. +In this case, it is necessary to use the value measure->averef to find the corresponding +average entry. + +Note that average->measureOffset and measure->averef are only valid for a given load of the +data: they refer to the sequence number in the data blocks. + +next_meas is a list of the equivalent sequence of the measure block as if it were sorted. + +to find the sequence of measurements for a given average: +n_0 = average->measureOffset +n_1 = next_meas[n_0] +n_i = next_meas[n_i-1] + +*/ /* build the initial links assuming the table is sorted, @@ -6,17 +29,17 @@ off_t i, j, N; - off_t *next; + off_t *next_meas; N = 0; - ALLOCATE (next, off_t, Nmeasure); + ALLOCATE (next_meas, off_t, Nmeasure); for (i = 0; i < Naverage; i++, N++) { for (j = 0; j < average[i].Nmeasure - 1; j++, N++) { - next[N] = N + 1; + next_meas[N] = N + 1; if (N >= Nmeasure) { fprintf (stderr, "WARNING: N out of bounds (1)\n"); } } - next[N] = -1; + next_meas[N] = -1; if (N >= Nmeasure) { fprintf (stderr, "WARNING: N out of bounds (2)\n"); @@ -28,59 +51,5 @@ } } - return (next); -} - -/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */ -int add_meas_link (Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE) { - - off_t k, m; - - /* if we have trouble, check validity of next[m] : m < Nmeasure */ - m = average[0].measureOffset; - - for (k = 0; k < average[0].Nmeasure - 1; k++) { - m = next[m]; - if (m >= NMEASURE) { - fprintf (stderr, "WARNING: m out of bounds (3)\n"); - } - } - - /* set up references */ - next[Nmeasure] = -1; - if (Nmeasure >= NMEASURE) { - fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n"); - } - - if (m == -1) { - average[0].measureOffset = Nmeasure; - } else { - next[m] = Nmeasure; - if (m >= NMEASURE) { - fprintf (stderr, "WARNING: m out of bounds (4)\n"); - } - } - - return (TRUE); -} - -Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next) { - - off_t i, k, n, N; - Measure *tmpmeasure; - - /* fix order of Measure (memory intensive, but fast) */ - N = 0; - ALLOCATE (tmpmeasure, Measure, Nmeasure); - for (i = 0; i < Naverage; i++) { - n = average[i].measureOffset; - average[i].measureOffset = N; - for (k = 0; k < average[i].Nmeasure; k++, N++) { - tmpmeasure[N] = measure[n]; - tmpmeasure[N].averef = i; - n = next[n]; - } - } - free (measure); - return (tmpmeasure); + return (next_meas); } @@ -95,7 +64,7 @@ off_t i, m, k, Nm, averef; - off_t *next; - - ALLOCATE (next, off_t, Nmeasure); + off_t *next_meas; + + ALLOCATE (next_meas, off_t, Nmeasure); /* reset the Nm, offset values for average */ @@ -108,5 +77,5 @@ averef = measure[Nm].averef; m = average[averef].measureOffset; - next[Nm] = -1; + next_meas[Nm] = -1; if (m == -1) { /* no links yet for source */ @@ -116,6 +85,6 @@ } - for (k = 0; next[m] != -1; k++) { - m = next[m]; + for (k = 0; next_meas[m] != -1; k++) { + m = next_meas[m]; if (m >= Nmeasure) { fprintf (stderr, "WARNING: m out of bounds (1)\n"); @@ -124,11 +93,69 @@ average[averef].Nmeasure = k + 2; - next[m] = Nm; + next_meas[m] = Nm; if (m >= Nmeasure) { fprintf (stderr, "WARNING: m out of bounds (2)\n"); } } - return (next); -} + return (next_meas); +} + +/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */ +int add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) { + + off_t k, m; + + /* if we have trouble, check validity of next_meas[m] : m < Nmeasure */ + m = average[0].measureOffset; + + for (k = 0; k < average[0].Nmeasure - 1; k++) { + m = next_meas[m]; + if (m >= NMEASURE) { + fprintf (stderr, "WARNING: m out of bounds (3)\n"); + } + } + + /* set up references */ + next_meas[Nmeasure] = -1; + if (Nmeasure >= NMEASURE) { + fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n"); + } + + if (m == -1) { + average[0].measureOffset = Nmeasure; + } else { + next_meas[m] = Nmeasure; + if (m >= NMEASURE) { + fprintf (stderr, "WARNING: m out of bounds (4)\n"); + } + } + + return (TRUE); +} + +Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) { + + off_t i, k, n, N; + Measure *tmpmeasure; + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeasure); + for (i = 0; i < Naverage; i++) { + n = average[i].measureOffset; + average[i].measureOffset = N; + for (k = 0; k < average[i].Nmeasure; k++, N++) { + if (n == -1) abort(); + tmpmeasure[N] = measure[n]; + if (measure[n].averef != i) abort(); + tmpmeasure[N].averef = i; + n = next_meas[n]; + } + } + free (measure); + return (tmpmeasure); +} + +/*******************************************************************************************/ /* build the initial links assuming the table is sorted */ @@ -136,15 +163,15 @@ off_t i, j, N; - off_t *next; + off_t *next_miss; N = 0; - ALLOCATE (next, off_t, Nmissing); + ALLOCATE (next_miss, off_t, Nmissing); for (i = 0; i < Naverage; i++) { for (j = 0; j < average[i].Nmissing - 1; j++, N++) { - next[N] = N + 1; + next_miss[N] = N + 1; } if (average[i].Nmissing > 0) { - next[N] = -1; + next_miss[N] = -1; if (N >= Nmissing) { fprintf (stderr, "overflow in init_missing_links"); @@ -155,8 +182,8 @@ } - return (next); -} - -int add_miss_link (Average *average, off_t *next, off_t Nmissing) { + return (next_miss); +} + +int add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) { off_t k, m; @@ -165,13 +192,13 @@ if (average[0].Nmissing < 1) { average[0].missingOffset = Nmissing; - next[Nmissing] = -1; + next_miss[Nmissing] = -1; return (TRUE); } m = average[0].missingOffset; - for (k = 0; k < average[0].Nmissing - 1; k++) m = next[m]; + for (k = 0; k < average[0].Nmissing - 1; k++) m = next_miss[m]; /* set up references */ - next[Nmissing] = -1; - next[m] = Nmissing; + next_miss[Nmissing] = -1; + next_miss[m] = Nmissing; return (TRUE); } @@ -180,5 +207,5 @@ we must always save the missing table, if it exists */ -Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next) { +Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss) { off_t i, k, n, N; @@ -193,5 +220,5 @@ for (k = 0; k < average[i].Nmissing; k++, N++) { tmpmissing[N] = missing[n]; - n = next[n]; + n = next_miss[n]; } } @@ -199,2 +226,154 @@ return (tmpmissing); } + +/*******************************************************************************************/ + +/* build the initial links assuming the table is sorted, + not partial, and has a correct set of average[].lensingOffset,Nlensing values */ +off_t *init_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) { + + off_t i, j, N; + off_t *next_lens; + + N = 0; + + ALLOCATE (next_lens, off_t, Nlensing); + for (i = 0; i < Naverage; i++) { + if (!average[i].Nlensing) continue; + for (j = 0; j < average[i].Nlensing - 1; j++, N++) { + next_lens[N] = N + 1; + if (N >= Nlensing) { + fprintf (stderr, "WARNING: N out of bounds (1)\n"); + } + } + next_lens[N] = -1; + if (N >= Nlensing) { + fprintf (stderr, "WARNING: N out of bounds (2)\n"); + } + + if (N >= Nlensing) { + fprintf (stderr, "overflow in init_lensing_links\n"); + abort (); + } + N++; + } + return (next_lens); +} + +/* construct lensing links which are valid FOR THIS LOAD + * - if we have a full load, we will get links which can + * be used by other programs (eg, relphot, etc) + * - if we have a partial load, the links are only valid + * for that partial load + */ + +off_t *build_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) { + + off_t i, m, k, Nm, averef; + off_t *next_lens; + + ALLOCATE (next_lens, off_t, Nlensing); + + /* reset the Nm, offset values for average */ + for (i = 0; i < Naverage; i++) { + average[i].lensingOffset = -1; + average[i].Nlensing = 0; + } + + for (Nm = 0; Nm < Nlensing; Nm++) { + averef = lensing[Nm].averef; + m = average[averef].lensingOffset; + next_lens[Nm] = -1; + + if (m == -1) { /* no links yet for source */ + average[averef].lensingOffset = Nm; + average[averef].Nlensing = 1; + continue; + } + + for (k = 0; next_lens[m] != -1; k++) { + m = next_lens[m]; + if (m >= Nlensing) { + fprintf (stderr, "WARNING: m out of bounds (1)\n"); + } + } + + average[averef].Nlensing = k + 2; + next_lens[m] = Nm; + if (m >= Nlensing) { + fprintf (stderr, "WARNING: m out of bounds (2)\n"); + } + } + return (next_lens); +} + +/* average[].lensingOffset, average[].Nlensing are valid within an addstar run */ +int add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) { + + off_t k, m; + + /* if we have trouble, check validity of next_lens[m] : m < Nlensing */ + m = average[0].lensingOffset; + + for (k = 0; k < average[0].Nlensing - 1; k++) { + m = next_lens[m]; + if (m >= NLENSING) { + fprintf (stderr, "WARNING: m out of bounds (3)\n"); + } + } + + /* set up references */ + next_lens[Nlensing] = -1; + if (Nlensing >= NLENSING) { + fprintf (stderr, "WARNING: Nlensing out of bounds (1)\n"); + } + + if (m == -1) { + average[0].lensingOffset = Nlensing; + } else { + next_lens[m] = Nlensing; + if (m >= NLENSING) { + fprintf (stderr, "WARNING: m out of bounds (4)\n"); + } + } + + return (TRUE); +} + +Lensing *sort_lensing (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next_lens) { + + off_t i, k, n, np, N; + Lensing *tmplensing; + + /* + for (i = 0; i < Naverage; i++) { + if (average[i].Nlensing != 4) { + fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nlensing, (int) average[i].lensingOffset); + } + } + */ + + /* fix order of Lensing (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmplensing, Lensing, Nlensing); + for (i = 0; i < Naverage; i++) { + if (!average[i].Nlensing) continue; + n = average[i].lensingOffset; + average[i].lensingOffset = N; + for (k = 0; k < average[i].Nlensing; k++, N++) { + if (n == -1) { + fprintf (stderr, "entry after %d has a problem\n", (int) np); + abort(); + } + tmplensing[N] = lensing[n]; + if (lensing[n].averef != i) abort(); + tmplensing[N].averef = i; + np = n; + n = next_lens[n]; + } + } + free (lensing); + return (tmplensing); +} + + Index: /trunk/Ohana/src/dvomerge/src/dvomergeContinue.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeContinue.c (revision 37116) @@ -131,6 +131,4 @@ merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap); - 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); Index: /trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (revision 37116) @@ -74,6 +74,4 @@ merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap); - 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); Index: /trunk/Ohana/src/dvomerge/src/dvomergeCreate.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeCreate.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeCreate.c (revision 37116) @@ -178,5 +178,4 @@ SkyListFree (inlist); - outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; dvo_catalog_save (&outcatalog, VERBOSE); dvo_catalog_unlock (&outcatalog); Index: /trunk/Ohana/src/dvomerge/src/dvomergeFromList.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeFromList.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeFromList.c (revision 37116) @@ -143,6 +143,4 @@ merge_catalogs_old (&skyregion, &outcatalog, &incatalog, RADIUS, secfiltMap); - outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; - dmhObjectAdd (outstat[0].history, &outcatalog.header, inStats); Index: /trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeUpdate.c (revision 37116) @@ -13,7 +13,5 @@ int NsecfiltInput, NsecfiltOutput; - double dtime; - struct timeval start, stop; - gettimeofday (&start, NULL); + INITTIME; CONTINUE = FALSE; @@ -121,13 +119,10 @@ } - gettimeofday (&stop, NULL); - dtime = DTIME (stop, start); - if (!status) { - fprintf (stderr, "ERROR: elapsed time %9.4f sec\n", dtime); + MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime); exit (3); } - fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); + MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime); exit (0); } Index: /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (revision 37116) @@ -163,7 +163,10 @@ LoadCatalog (&outcatalog, outlist[0].regions[j], outcatalog.filename, "w", NsecfiltOutput); + // if no catalog already exists, use the input catalog to define the format + if (outcatalog.Naves_disk == 0) { + outcatalog.catformat = incatalog.catformat; + } + merge_catalogs_old (outlist[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap); - - outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; if (outstat[j].missed) { Index: /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c (revision 37116) @@ -72,6 +72,4 @@ merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap); - 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); @@ -115,7 +113,5 @@ int NsecfiltInput, NsecfiltOutput; - double dtime; - struct timeval start, stop; - gettimeofday (&start, NULL); + INITTIME; CONTINUE = FALSE; @@ -134,4 +130,5 @@ fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does not handle merge tracking yet\n"); + fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does work in a parallel context\n"); exit (2); @@ -293,8 +290,5 @@ } - gettimeofday (&stop, NULL); - dtime = DTIME (stop, start); - fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime); - + MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime); exit (0); @@ -318,4 +312,5 @@ } } + MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime); exit (1); } Index: /trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c (revision 37116) @@ -8,10 +8,8 @@ // input entries always define new objects -#define notyet 1 - int merge_catalogs_new (SkyRegion *region, Catalog *output, Catalog *input, int *secfiltMap) { off_t i, j, offset; - off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, NsecfiltIn, NsecfiltOut, Nm; + off_t NAVERAGE, NMEASURE, NLENSING, Naverage, Nmeasure, Nlensing, NsecfiltIn, NsecfiltOut, Nm; Naverage = output[0].Naverage; @@ -63,4 +61,21 @@ } output[0].average[Naverage].Nmeasure = Nm; + + Nm = 0; + for (j = 0; j < input[0].average[i].Nlensing; j++) { + offset = input[0].average[i].lensingOffset + j; + + output[0].lensing[Nlensing] = input[0].lensing[offset]; + output[0].lensing[Nlensing].averef = Naverage; + + Nlensing ++; + Nm ++; + if (Nlensing == NLENSING) { + NLENSING += 1000; + REALLOCATE (output[0].lensing, Lensing, NLENSING); + } + } + output[0].average[Naverage].Nlensing = Nm; + Naverage ++; if (Naverage == NAVERAGE) { @@ -72,7 +87,8 @@ REALLOCATE (output[0].average, Average, MAX (Naverage, 1)); REALLOCATE (output[0].measure, Measure, MAX (Nmeasure, 1)); + REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1)); REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1)); output[0].Naverage = Naverage; - output[0].Nmeasure = Nmeasure; + output[0].Nlensing = Nlensing; output[0].Nsecf_mem = Naverage * NsecfiltOut; @@ -84,8 +100,9 @@ if (VERBOSE) { - fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n", + fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing) for catalog\n", i, output[0].Naverage, - output[0].Nmeasure); + output[0].Nmeasure, + output[0].Nlensing); } return (TRUE); Index: /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 37115) +++ /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 37116) @@ -14,6 +14,6 @@ double *X1, *Y1, *X2, *Y2; double dX, dY, dR; - off_t *N1, *N2, *next_meas; - off_t Nave, NAVE, Nmeas, NMEAS, Nmatch; + off_t *N1, *N2, *next_meas, *next_lens; + off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS; int NsecfiltIn; int NsecfiltOut; @@ -52,4 +52,5 @@ Nmatch = 0; NMEAS = Nmeas = output[0].Nmeasure; + NLENS = Nlens = output[0].Nlensing; // current max obj ID for this catalog @@ -113,5 +114,4 @@ } - /* set up pointers for linked list of measure */ if (output[0].sorted && (output[0].Nmeasure >= output[0].Nmeas_disk)) { @@ -119,6 +119,8 @@ // is sorted while processed next_meas = init_measure_links (output[0].average, Nave, output[0].measure, Nmeas); + next_lens = init_lensing_links (output[0].average, Nave, output[0].lensing, Nlens); } else { next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas); + next_lens = build_lensing_links (output[0].average, Nave, output[0].lensing, Nlens); } @@ -182,4 +184,9 @@ REALLOCATE (output[0].measure, Measure, NMEAS); } + if (Nlens + input[0].average[N].Nmeasure >= NLENS) { + NLENS = Nlens + input[0].average[N].Nmeasure + 1000; + REALLOCATE (next_lens, off_t, NLENS); + REALLOCATE (output[0].lensing, Lensing, NLENS); + } // 4) average properties from the input and the output db need to be properly merged. @@ -191,4 +198,5 @@ if (REPLACE_BY_PHOTCODE) { // index to first measure for this object + // XXX this does not support lensing measurements int Mout = output[0].average[n].measureOffset; if (replace_match(&output[0].average[n], &output[0].measure[Mout], &input[0].average[N], &input[0].measure[offset])) { @@ -248,4 +256,22 @@ } + // if lensing measurements exist, add them too + if (output[0].lensing) { + for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin++) { + /* add to end of lensing list */ + add_lens_link (&output[0].average[n], next_lens, Nlens, NLENS); + + // set the new lensing + off_t lensoff = input[0].average[N].lensingOffset + Nin; + output[0].lensing[Nlens] = input[0].lensing[lensoff]; + + output[0].lensing[Nlens].averef = n; + output[0].lensing[Nlens].objID = output[0].average[n].objID; + output[0].lensing[Nlens].catID = output[0].catID; + output[0].average[n].Nlensing ++; + Nlens ++; + } + } + // update the average properties to reflect the incoming entries: // if the original value is NAN but the input value is not, accept the input: @@ -288,4 +314,9 @@ REALLOCATE (output[0].measure, Measure, NMEAS); } + if (Nlens + input[0].average[N].Nlensing >= NLENS) { + NLENS = Nlens + input[0].average[N].Nlensing + 1000; + REALLOCATE (next_lens, off_t, NLENS); + REALLOCATE (output[0].lensing, Lensing, NLENS); + } if (Nave >= NAVE) { NAVE = Nave + 1000; @@ -355,4 +386,31 @@ next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1; } + + /** add lensing for this input average object **/ + if (output[0].lensing) { + output[0].average[Nave].lensingOffset = Nlens; + for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin ++) { + // supply the lensing values from this detection + off_t lensoff = input[0].average[N].lensingOffset + Nin; + output[0].lensing[Nlens] = input[0].lensing[lensoff]; + + // the following lensing elements cannot be set until here: + output[0].lensing[Nlens].averef = Nave; + output[0].lensing[Nlens].objID = output[0].average[Nave].objID; + output[0].lensing[Nlens].catID = output[0].catID; + + // as we add lensing, update Nlensing to match + output[0].average[Nave].Nlensing ++; + + /* we set next[Nlens] to -1 here, and update correctly below */ + next_lens[Nlens] = -1; + Nlens ++; + } + int Ngroup = input[0].average[N].Nlensing; + for (j = 0; j < Ngroup - 1; j++) { + next_lens[Nlens - Ngroup + j] = Nlens - Ngroup + j + 1; + } + } + Nave ++; } @@ -362,4 +420,5 @@ REALLOCATE (output[0].average, Average, Nave); REALLOCATE (output[0].measure, Measure, Nmeas); + REALLOCATE (output[0].lensing, Lensing, Nlens); # define NOSORT 0 @@ -369,4 +428,5 @@ output[0].sorted = TRUE; output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas); + output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens); } @@ -384,8 +444,10 @@ output[0].Naverage = Nave; output[0].Nmeasure = Nmeas; + output[0].Nlensing = Nlens; output[0].Nsecf_mem = Nave*NsecfiltOut; - if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nmatch); + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nlens, Nmatch); free (next_meas); + free (next_lens); free (X2); Index: /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c =================================================================== --- /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 37115) +++ /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 37116) @@ -352,5 +352,5 @@ PrintIOBuffer (sec_buffer, "%hd, ", secfilt->ubercalDist); PrintIOBuffer (sec_buffer, "%.6f, ", secfilt->Mstdev); - PrintIOBuffer (sec_buffer, "%lu ", secfilt->stackDetectID); + PrintIOBuffer (sec_buffer, "%lu ", 0L); // secfilt->stackDetectID PrintIOBuffer (sec_buffer, "),\n"); secfilt ++; Index: /trunk/Ohana/src/libautocode/def/secfilt-ps1-v4.d =================================================================== --- /trunk/Ohana/src/libautocode/def/secfilt-ps1-v4.d (revision 37115) +++ /trunk/Ohana/src/libautocode/def/secfilt-ps1-v4.d (revision 37116) @@ -23,3 +23,4 @@ FIELD ubercalDist, UBERCAL_DIST, short, number of images from an ubercal-image FIELD Mstdev, MAG_STDEV, short, standard deviation of measurements, millimags -FIELD stackDetectID, STACK_DETECT_ID, uint64_t, detection ID on stack used for fluxes (if any) +FIELD stackPrmryOff, STACK_PRIMARY_OFF, int, measure entry which is primary stack detection +FIELD stackBestOff, STACK_BEST_OFF, int, measure entry which is best stack detection Index: /trunk/Ohana/src/libautocode/def/secfilt-ps1-v5.d =================================================================== --- /trunk/Ohana/src/libautocode/def/secfilt-ps1-v5.d (revision 37115) +++ /trunk/Ohana/src/libautocode/def/secfilt-ps1-v5.d (revision 37116) @@ -42,6 +42,6 @@ FIELD dFapStk, FLUX_AP_STK_ERR, float, mean flux ap error -## need an 8-byte boundary here ## -FIELD stackDetectID, STACK_DETECT_ID, uint64_t, detection ID on stack used for fluxes (if any) +FIELD stackPrmryOff, STACK_PRIMARY_OFF, int, measure entry which is primary stack detection +FIELD stackBestOff, STACK_BEST_OFF, int, measure entry which is best stack detection ## *** this section is for forced-warp mean values *** Index: /trunk/Ohana/src/libautocode/def/secfilt.d =================================================================== --- /trunk/Ohana/src/libautocode/def/secfilt.d (revision 37115) +++ /trunk/Ohana/src/libautocode/def/secfilt.d (revision 37116) @@ -42,6 +42,6 @@ FIELD dFapStk, FLUX_AP_STK_ERR, float, mean flux ap error -## need an 8-byte boundary here ## -FIELD stackDetectID, STACK_DETECT_ID, uint64_t, detection ID on stack used for fluxes (if any) +FIELD stackPrmryOff, STACK_PRIMARY_OFF, int, measure entry which is primary stack detection +FIELD stackBestOff, STACK_BEST_OFF, int, measure entry which is best stack detection ## *** this section is for forced-warp mean values *** Index: /trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- /trunk/Ohana/src/libdvo/include/dvo.h (revision 37115) +++ /trunk/Ohana/src/libdvo/include/dvo.h (revision 37116) @@ -128,5 +128,5 @@ 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_WARP_USED = 0x00000400, // measurement used to find mean warp photometry ID_MEAS_UNDEF_4 = 0x00000800, // unused ID_MEAS_BLEND_MEAS_X = 0x00001000, // detection is within radius of multiple objects across catalogs @@ -383,12 +383,12 @@ char ***name; // projection cell name - int NX_SUB; - int NY_SUB; + float NX_SUB; + float NY_SUB; double dPix; double **Xo; double **Yo; - int **dX; - int **dY; + float **dX; + float **dY; } BoundaryTree; @@ -411,9 +411,9 @@ double Do; double dPix; - int dX; - int dY; - - int NX_SUB; - int NY_SUB; + float dX; + float dY; + + float NX_SUB; + float NY_SUB; char *basename; @@ -710,9 +710,11 @@ float PhotFluxInstErr (Measure *measure, dvoMagClassType class); float PhotFluxCatErr (Measure *measure, dvoMagClassType class); +float PhotFluxAveErr (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source); float PhotXm (PhotCode *code, Average *average, SecFilt *secfilt); float PhotZeroPoint (Measure *measure, Average *average, SecFilt *secfilt); -float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt); +int PhotNphot (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source); +float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source); float PhotMmin (PhotCode *code, Average *average, SecFilt *secfilt); float PhotMmax (PhotCode *code, Average *average, SecFilt *secfilt); Index: /trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- /trunk/Ohana/src/libdvo/include/dvodb.h (revision 37115) +++ /trunk/Ohana/src/libdvo/include/dvodb.h (revision 37116) @@ -54,7 +54,9 @@ MAG_OPTION_NPHOT, // Nused MAG_OPTION_UC_DIST, - MAG_OPTION_STACK_DET_ID, MAG_OPTION_FLAGS, } dvoMagOptionType; + +// MAG_OPTION_STACK_PRIMARY_OFF, +// MAG_OPTION_STACK_BEST_OFF, typedef enum { Index: /trunk/Ohana/src/libdvo/src/BoundaryTree.c =================================================================== --- /trunk/Ohana/src/libdvo/src/BoundaryTree.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/BoundaryTree.c (revision 37116) @@ -84,6 +84,6 @@ gfits_scan (headerPHU, "DEC_ORI", "%lf", 1, &tree->DEC_origin); gfits_scan (headerPHU, "DEC_OFF", "%lf", 1, &tree->DEC_offset); - gfits_scan (headerPHU, "NX_SUB", "%d", 1, &tree->NX_SUB); - gfits_scan (headerPHU, "NY_SUB", "%d", 1, &tree->NY_SUB); + gfits_scan (headerPHU, "NX_SUB", "%f", 1, &tree->NX_SUB); + gfits_scan (headerPHU, "NY_SUB", "%f", 1, &tree->NY_SUB); gfits_scan (headerPHU, "PIXSCALE", "%lf", 1, &tree->dPix); @@ -111,6 +111,6 @@ ALLOCATE (tree->Xo, double *, tree->Nzone); ALLOCATE (tree->Yo, double *, tree->Nzone); - ALLOCATE (tree->dX, int *, tree->Nzone); - ALLOCATE (tree->dY, int *, tree->Nzone); + ALLOCATE (tree->dX, float *, tree->Nzone); + ALLOCATE (tree->dY, float *, tree->Nzone); ALLOCATE (tree->cell, int *, tree->Nzone); ALLOCATE (tree->projID, int *, tree->Nzone); @@ -121,6 +121,6 @@ ALLOCATE (tree->Xo[i], double, tree->Nband[i]); ALLOCATE (tree->Yo[i], double, tree->Nband[i]); - ALLOCATE (tree->dX[i], int, tree->Nband[i]); - ALLOCATE (tree->dY[i], int, tree->Nband[i]); + ALLOCATE (tree->dX[i], float, tree->Nband[i]); + ALLOCATE (tree->dY[i], float, tree->Nband[i]); ALLOCATE (tree->cell[i], int, tree->Nband[i]); ALLOCATE (tree->projID[i], int, tree->Nband[i]); @@ -147,6 +147,6 @@ GET_COLUMN_NEW(&headerCell, &ftableCell, Xo, "X_CENT", double); GET_COLUMN_NEW(&headerCell, &ftableCell, Yo, "Y_CENT", double); - GET_COLUMN_NEW(&headerCell, &ftableCell, dX, "X_GRID", int); - GET_COLUMN_NEW(&headerCell, &ftableCell, dY, "Y_GRID", int); + GET_COLUMN_NEW(&headerCell, &ftableCell, dX, "X_GRID", float); + GET_COLUMN_NEW(&headerCell, &ftableCell, dY, "Y_GRID", float); GET_COLUMN_NEW(&headerCell, &ftableCell, name, "NAME", char); // XXX how is this done? gfits_free_header (&headerCell); @@ -214,6 +214,6 @@ gfits_modify (&header, "DEC_OFF", "%lf", 1, tree->DEC_offset); - gfits_modify (&header, "NX_SUB", "%d", 1, tree->NX_SUB); - gfits_modify (&header, "NY_SUB", "%d", 1, tree->NY_SUB); + gfits_modify (&header, "NX_SUB", "%f", 1, tree->NX_SUB); + gfits_modify (&header, "NY_SUB", "%f", 1, tree->NY_SUB); gfits_modify (&header, "PIXSCALE", "%lf", 1, tree->dPix); @@ -289,6 +289,6 @@ gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, fmt, "NAME", "cell name", "none", 1.0, 0.0); @@ -310,6 +310,6 @@ double *Xo ; ALLOCATE (Xo, double, Ncell); double *Yo ; ALLOCATE (Yo, double, Ncell); - int *dX ; ALLOCATE (dX, int, Ncell); - int *dY ; ALLOCATE (dY, int, Ncell); + float *dX ; ALLOCATE (dX, float, Ncell); + float *dY ; ALLOCATE (dY, float, Ncell); char *name ; ALLOCATE (name, char, Ncell*BOUNDARY_TREE_NAME_LENGTH); Index: /trunk/Ohana/src/libdvo/src/TessellationTable.c =================================================================== --- /trunk/Ohana/src/libdvo/src/TessellationTable.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/TessellationTable.c (revision 37116) @@ -28,9 +28,9 @@ tess[i].dPix = NAN; - tess[i].dX = 0; - tess[i].dY = 0; - - tess[i].NX_SUB = 0; - tess[i].NY_SUB = 0; + tess[i].dX = NAN; + tess[i].dY = NAN; + + tess[i].NX_SUB = NAN; + tess[i].NY_SUB = NAN; tess[i].basename = NULL; @@ -123,6 +123,6 @@ GET_COLUMN_NEW(Xo, "X_CENT", double); GET_COLUMN_NEW(Yo, "Y_CENT", double); - GET_COLUMN_NEW(dX, "X_GRID", int); - GET_COLUMN_NEW(dY, "Y_GRID", int); + GET_COLUMN_NEW(dX, "X_GRID", float); + GET_COLUMN_NEW(dY, "Y_GRID", float); GET_COLUMN_NEW(Rmin, "R_MIN", double); GET_COLUMN_NEW(Rmax, "R_MAX", double); @@ -130,6 +130,6 @@ GET_COLUMN_NEW(Dmax, "D_MAX", double); GET_COLUMN_NEW(dPix, "SCALE", double); - GET_COLUMN_NEW(NX_SUB, "NX_SUB", int); - GET_COLUMN_NEW(NY_SUB, "NY_SUB", int); + GET_COLUMN_NEW(NX_SUB, "NX_SUB", float); + GET_COLUMN_NEW(NY_SUB, "NY_SUB", float); GET_COLUMN_NEW(TYPE, "TYPE", int); @@ -280,6 +280,6 @@ gfits_modify (&header, "DEC_OFF", "%lf", 1, rings->DEC_offset); - gfits_modify (&header, "NX_SUB", "%d", 1, rings->NX_SUB); - gfits_modify (&header, "NY_SUB", "%d", 1, rings->NY_SUB); + gfits_modify (&header, "NX_SUB", "%f", 1, rings->NX_SUB); + gfits_modify (&header, "NY_SUB", "%f", 1, rings->NY_SUB); gfits_modify (&header, "PIXSCALE", "%lf", 1, rings->dPix); } @@ -300,6 +300,6 @@ gfits_define_bintable_column (&theader, "D", "X_CENT", "projection cell center pixel", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "Y_CENT", "projection cell center pixel", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_GRID", "skycell grid spacing", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_GRID", "skycell grid spacing", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "R_MIN", "RA limit (lower)", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "R_MAX", "RA limit (upper)", "none", 1.0, 0.0); @@ -307,6 +307,6 @@ gfits_define_bintable_column (&theader, "D", "D_MAX", "DEC limit (upper)", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "D", "SCALE", "pixel scale for projection cell", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "NX_SUB", "skycell subdivision in x", "none", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "NY_SUB", "skycell subdivision in y", "none", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "TYPE", "type of tessellation", "none", 1.0, 0.0); @@ -330,6 +330,6 @@ double *Xo ; ALLOCATE (Xo , double, Nout); double *Yo ; ALLOCATE (Yo , double, Nout); - int *dX ; ALLOCATE (dX , int, Nout); - int *dY ; ALLOCATE (dY , int, Nout); + float *dX ; ALLOCATE (dX , float, Nout); + float *dY ; ALLOCATE (dY , float, Nout); double *Rmin ; ALLOCATE (Rmin , double, Nout); double *Rmax ; ALLOCATE (Rmax , double, Nout); @@ -337,6 +337,6 @@ double *Dmax ; ALLOCATE (Dmax , double, Nout); double *dPix ; ALLOCATE (dPix , double, Nout); - int *NX_SUB ; ALLOCATE (NX_SUB, int, Nout); - int *NY_SUB ; ALLOCATE (NY_SUB, int, Nout); + float *NX_SUB ; ALLOCATE (NX_SUB, float, Nout); + float *NY_SUB ; ALLOCATE (NY_SUB, float, Nout); int *TYPE ; ALLOCATE (TYPE, int, Nout); @@ -559,4 +559,8 @@ int xi = x / tess[myTess].tree->dX[zone][band]; int yi = y / tess[myTess].tree->dY[zone][band]; + + xi = MAX(MIN(xi, tess[myTess].tree->NX_SUB - 1.0), 0); + yi = MAX(MIN(yi, tess[myTess].tree->NY_SUB - 1.0), 0); + int N = xi + tess[myTess].tree->NX_SUB * yi; Index: /trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 37116) @@ -200,4 +200,14 @@ } + /* in some cases below we need a measurement which matches the requested photcode. + here is some sample code to find this. + + for (i = 0; i < average[0].Nmeasure; i++) { + if (field->photcode->code != measure[i].photcode) continue; + value.Flt = measure[i].M; + break; + } + */ + switch (field->magOption) { case MAG_OPTION_MAG: @@ -207,5 +217,5 @@ break; case MAG_LEVEL_REF: - // why is measure needed here? + // which measure is needed here? // PhotRef (field->photcode, average, secfilt, measure, field->magClass, field->magSource); break; @@ -242,5 +252,6 @@ break; case MAG_LEVEL_REF: - // value.Flt = PhotFluxRef (equiv, average, secfilt, measure, field->magClass, field->magSource); + // XXX which measure is needed here? + // value.Flt = PhotFluxRef (field->photcode, average, secfilt, measure, field->magClass, field->magSource); break; case MAG_LEVEL_INST: @@ -258,5 +269,5 @@ case MAG_LEVEL_AVE: case MAG_LEVEL_REF: - // value.Flt = PhotFluxAveErr (equiv, average, secfilt, field->magClass, field->magSource); + value.Flt = PhotFluxAveErr (field->photcode, average, secfilt, field->magClass, field->magSource); break; case MAG_LEVEL_INST: @@ -271,7 +282,15 @@ case MAG_OPTION_STDEV: + value.Flt = PhotMstdev (field->photcode, average, secfilt, field->magClass, field->magSource); + break; case MAG_OPTION_CHISQ: + value.Flt = PhotXm (field->photcode, average, secfilt); + break; case MAG_OPTION_MIN: + value.Flt = PhotMmin (field->photcode, average, secfilt); + break; case MAG_OPTION_MAX: + value.Flt = PhotMmax (field->photcode, average, secfilt); + break; case MAG_OPTION_NCODE: { int Nsec = GetPhotcodeNsec (field->photcode->code); @@ -281,13 +300,13 @@ } case MAG_OPTION_NPHOT: { - int Nsec = GetPhotcodeNsec (field->photcode->code); - if (Nsec == -1) break; - value.Int = secfilt[Nsec].Nused; + value.Int = PhotNphot (field->photcode, average, secfilt, field->magClass, field->magSource); break; } case MAG_OPTION_UC_DIST: - case MAG_OPTION_STACK_DET_ID: - value.Int = PhotStackID (field->photcode, average, secfilt); - break; + value.Flt = PhotUCdist (field->photcode, average, secfilt); + break; + // XX case MAG_OPTION_STACK_DET_ID: + // XX value.Int = PhotStackID (field->photcode, average, secfilt); + // XX break; case MAG_OPTION_FLAGS: { int Nsec = GetPhotcodeNsec (field->photcode->code); @@ -301,124 +320,4 @@ break; -# if (0) - switch (field->magMode) { - case MAG_AVE: - value.Flt = PhotAve (field->photcode, average, secfilt); - break; - case MAG_REF: - value.Flt = PhotRef (field->photcode, average, secfilt, measure); - break; - case MAG_INST: - case MAG_CAT: - case MAG_SYS: - case MAG_REL: - case MAG_CAL: - // XXX need to code this correctly: this returns just the first matching value - value.Flt = NAN; - for (i = 0; i < average[0].Nmeasure; i++) { - if (field->photcode->code != measure[i].photcode) continue; - value.Flt = measure[i].M; - break; - } - break; - case MAG_CHISQ: - value.Flt = PhotXm (field->photcode, average, secfilt); - break; - case MAG_ERR: - if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) { - for (i = 0; i < average[0].Nmeasure; i++) { - if (field->photcode->code != measure[i].photcode) continue; - value.Flt = measure[i].dM; - break; - } - } else { - value.Flt = PhotAveErr (field->photcode, average, secfilt); - } - break; - case MAG_STDEV: - value.Flt = PhotMstdev (field->photcode, average, secfilt); - break; - case MAG_PHOT_FLAGS: - if ((field->photcode->type == PHOT_REF) || (field->photcode->type == PHOT_DEP)) { - for (i = 0; i < average[0].Nmeasure; i++) { - if (field->photcode->code != measure[i].photcode) continue; - value.Int = measure[i].photFlags; - break; - } - } else { - int Nsec = GetPhotcodeNsec (field->photcode->code); - if (Nsec == -1) break; - value.Int = secfilt[Nsec].flags; - } - break; - case MAG_NCODE: - // XXX push these into dvo_photcode_ops APIs - // XXX do I need to allow for conversion to equiv? - int Nsec = GetPhotcodeNsec (field->photcode->code); - if (Nsec == -1) break; - value.Int = secfilt[Nsec].Ncode; - break; - case MAG_NPHOT: - int Nsec = GetPhotcodeNsec (field->photcode->code); - if (Nsec == -1) break; - value.Int = secfilt[Nsec].Nused; - break; - - case MAG_APER: - case MAG_APER_AVE: - value.Flt = PhotAperAve (field->photcode, average, secfilt); - break; - - case MAG_APER_REF: - value.Flt = PhotAperRef (field->photcode, average, secfilt, measure); - break; - - case MAG_KRON: - case MAG_KRON_AVE: - value.Flt = PhotKronAve (field->photcode, average, secfilt); - break; - - case MAG_KRON_REF: - value.Flt = PhotKronRef (field->photcode, average, secfilt, measure); - break; - - case MAG_KRON_ERR: - value.Flt = PhotKronAveErr (field->photcode, average, secfilt); - break; - - case MAG_MIN: - value.Flt = PhotMmin (field->photcode, average, secfilt); - break; - case MAG_MAX: - value.Flt = PhotMmax (field->photcode, average, secfilt); - break; - case MAG_UC_DIST: - value.Flt = PhotUCdist (field->photcode, average, secfilt); - break; - case MAG_STACK_DET_ID: - value.Int = PhotStackID (field->photcode, average, secfilt); - break; - - case MAG_FLUX_PSF: - value.Flt = PhotAveFluxPSF (field->photcode, average, secfilt); - break; - case MAG_FLUX_PSF_ERR: - value.Flt = PhotAvedFluxPSF (field->photcode, average, secfilt); - break; - case MAG_FLUX_KRON: - value.Flt = PhotAveFluxKron (field->photcode, average, secfilt); - break; - case MAG_FLUX_KRON_ERR: - value.Flt = PhotAvedFluxKron (field->photcode, average, secfilt); - break; - } - break; - case AVE_dMAG: - value.Flt = PhotAveErr (field->photcode, average, secfilt); - break; - case AVE_Xm: - value.Flt = PhotXm (field->photcode, average, secfilt); - break; -# endif case AVE_TYPE: break; Index: /trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37116) @@ -248,5 +248,4 @@ case MAG_OPTION_NPHOT: case MAG_OPTION_UC_DIST: - case MAG_OPTION_STACK_DET_ID: case MAG_OPTION_FLAGS: case MAG_OPTION_NONE: Index: /trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dbFields.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dbFields.c (revision 37116) @@ -78,5 +78,4 @@ if (!strcasecmp (string, "nphot")) return (MAG_OPTION_NPHOT); if (!strcasecmp (string, "uc_dist")) return (MAG_OPTION_UC_DIST); - if (!strcasecmp (string, "stack_det_id")) return (MAG_OPTION_STACK_DET_ID); if (!strcasecmp (string, "flags")) return (MAG_OPTION_FLAGS); @@ -145,5 +144,4 @@ if (!strcasecmp (string, "max")) return (MAG_MAX); if (!strcasecmp (string, "ucdist")) return (MAG_UC_DIST); - if (!strcasecmp (string, "stackDetectID")) return (MAG_STACK_DET_ID); if (!strcasecmp (string, "fluxpsf")) return (MAG_FLUX_PSF); if (!strcasecmp (string, "fluxpsferr")) return (MAG_FLUX_PSF_ERR); @@ -296,5 +294,4 @@ switch (field->magOption) { - case MAG_OPTION_STACK_DET_ID: case MAG_OPTION_NCODE: case MAG_OPTION_NPHOT: Index: /trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 37116) @@ -197,5 +197,6 @@ secfilt->dFapStk = NAN; - secfilt->stackDetectID = 0; + secfilt->stackPrmryOff = -1; + secfilt->stackBestOff = -1; secfilt->MpsfWrp = NAN; @@ -860,11 +861,7 @@ catalog[0].lensobj = NULL; } - // if (catalog[0].X_t) free (catalog[0].X_t); - // if (catalog[0].Y_t) free (catalog[0].Y_t); - if (catalog[0].nOwn_t) free (catalog[0].nOwn_t); - if (catalog[0].found_t) free (catalog[0].found_t); - // if (catalog[0].image_t) free (catalog[0].image_t); - // if (catalog[0].mosaic_t) free (catalog[0].mosaic_t); - if (catalog[0].foundWarp_t) free (catalog[0].foundWarp_t); + if (catalog[0].nOwn_t) { free (catalog[0].nOwn_t); catalog[0].nOwn_t = NULL; } + if (catalog[0].found_t) { free (catalog[0].found_t); catalog[0].found_t = NULL; } + if (catalog[0].foundWarp_t) { free (catalog[0].foundWarp_t); catalog[0].foundWarp_t = NULL; } } Index: /trunk/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 37116) @@ -392,4 +392,11 @@ return (lensing); } +# define SKIPPING_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, NAME)) { \ + fprintf (stderr, "ERROR: format %s not defined for lensing, skipping\n", NAME); \ + *Nlensing = 0; \ + return NULL; \ + } + if (!strcmp (extname, "DVO_LENSING")) { lensing = gfits_table_get_Lensing (ftable, Nlensing, NULL); @@ -402,16 +409,16 @@ } -// CONVERT_FORMAT ("DVO_LENSING_PS1_REF", PS1_REF, PS1_REF); -// CONVERT_FORMAT ("DVO_LENSING_ELIXIR", ELIXIR, Elixir); -// CONVERT_FORMAT ("DVO_LENSING_LONEOS", LONEOS, Loneos); -// CONVERT_FORMAT ("DVO_LENSING_PANSTARRS_DEV_0", PANSTARRS_DEV_0, Panstarrs_DEV_0); -// CONVERT_FORMAT ("DVO_LENSING_PANSTARRS_DEV_1", PANSTARRS_DEV_1, Panstarrs_DEV_1); -// CONVERT_FORMAT ("DVO_LENSING_PS1_DEV_1", PS1_DEV_1, PS1_DEV_1); -// CONVERT_FORMAT ("DVO_LENSING_PS1_DEV_2", PS1_DEV_2, PS1_DEV_2); -// CONVERT_FORMAT ("DVO_LENSING_PS1_V1", PS1_V1, PS1_V1); -// CONVERT_FORMAT ("DVO_LENSING_PS1_V2", PS1_V2, PS1_V2); -// CONVERT_FORMAT ("DVO_LENSING_PS1_V3", PS1_V3, PS1_V3); -// CONVERT_FORMAT ("DVO_LENSING_PS1_V4", PS1_V4, PS1_V4); - CONVERT_FORMAT ("DVO_LENSING_PS1_V5", PS1_V5, PS1_V5); + SKIPPING_FORMAT ("DVO_LENSING_PS1_REF", PS1_REF, PS1_REF); + SKIPPING_FORMAT ("DVO_LENSING_ELIXIR", ELIXIR, Elixir); + SKIPPING_FORMAT ("DVO_LENSING_LONEOS", LONEOS, Loneos); + SKIPPING_FORMAT ("DVO_LENSING_PANSTARRS_DEV_0", PANSTARRS_DEV_0, Panstarrs_DEV_0); + SKIPPING_FORMAT ("DVO_LENSING_PANSTARRS_DEV_1", PANSTARRS_DEV_1, Panstarrs_DEV_1); + SKIPPING_FORMAT ("DVO_LENSING_PS1_DEV_1", PS1_DEV_1, PS1_DEV_1); + SKIPPING_FORMAT ("DVO_LENSING_PS1_DEV_2", PS1_DEV_2, PS1_DEV_2); + SKIPPING_FORMAT ("DVO_LENSING_PS1_V1", PS1_V1, PS1_V1); + SKIPPING_FORMAT ("DVO_LENSING_PS1_V2", PS1_V2, PS1_V2); + SKIPPING_FORMAT ("DVO_LENSING_PS1_V3", PS1_V3, PS1_V3); + SKIPPING_FORMAT ("DVO_LENSING_PS1_V4", PS1_V4, PS1_V4); + CONVERT_FORMAT ("DVO_LENSING_PS1_V5", PS1_V5, PS1_V5); # undef CONVERT_FORMAT Index: /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 37116) @@ -251,5 +251,6 @@ out[i].Mstdev = in[i].Mstdev; out[i].ubercalDist = in[i].ubercalDist; - out[i].stackDetectID = in[i].stackDetectID; + out[i].stackPrmryOff = in[i].stackPrmryOff; + out[i].stackBestOff = in[i].stackBestOff; } return (out); @@ -281,5 +282,6 @@ out[i].Mstdev = in[i].Mstdev; out[i].ubercalDist = in[i].ubercalDist; - out[i].stackDetectID = in[i].stackDetectID; + out[i].stackPrmryOff = in[i].stackPrmryOff; + out[i].stackBestOff = in[i].stackBestOff; } return (out); Index: /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 37116) @@ -319,5 +319,6 @@ out[i].dFapStk = in[i].dFapStk; - out[i].stackDetectID = in[i].stackDetectID; + out[i].stackPrmryOff = in[i].stackPrmryOff; + out[i].stackBestOff = in[i].stackBestOff; out[i].MpsfWrp = in[i].MpsfWrp; @@ -387,5 +388,6 @@ out[i].dFapStk = in[i].dFapStk; - out[i].stackDetectID = in[i].stackDetectID; + out[i].stackPrmryOff = in[i].stackPrmryOff; + out[i].stackBestOff = in[i].stackBestOff; out[i].MpsfWrp = in[i].MpsfWrp; Index: /trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 37115) +++ /trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 37116) @@ -687,14 +687,94 @@ } -float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - float Mstdev; - - if (code == NULL) return NAN; - - Ns = photcodes[0].hashNsec[code[0].code]; - Mstdev = (Ns == -1) ? NAN : 0.001*secfilt[Ns].Mstdev; +float PhotMstdev (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) { + + if (code == NULL) return NAN; + + int Ns = photcodes[0].hashNsec[code[0].code]; + if (Ns == -1) return (NAN); + + float Mstdev = NAN; + switch (source) { + case MAG_SRC_CHP: + switch (class) { + case MAG_CLASS_PSF: + Mstdev = secfilt[Ns].Mstdev; + break; + case MAG_CLASS_KRON: + Mstdev = secfilt[Ns].sMkron; + break; + case MAG_CLASS_APER: + Mstdev = secfilt[Ns].sMap; + break; + default: + break; + } + break; + case MAG_SRC_WRP: + switch (class) { + case MAG_CLASS_PSF: + Mstdev = secfilt[Ns].sFpsfWrp; + break; + case MAG_CLASS_KRON: + Mstdev = secfilt[Ns].sFkronWrp; + break; + case MAG_CLASS_APER: + Mstdev = secfilt[Ns].sFapWrp; + break; + default: + break; + } + break; + case MAG_SRC_STK: + default: + break; + } return (Mstdev); +} + +int PhotNphot (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) { + + if (code == NULL) return 0; + + int Ns = photcodes[0].hashNsec[code[0].code]; + if (Ns == -1) return 0; + + int Nphot = 0; + switch (source) { + case MAG_SRC_CHP: + switch (class) { + case MAG_CLASS_PSF: + Nphot = secfilt[Ns].Nused; + break; + case MAG_CLASS_KRON: + Nphot = secfilt[Ns].NusedKron; + break; + case MAG_CLASS_APER: + Nphot = secfilt[Ns].NusedAp; + break; + default: + break; + } + break; + case MAG_SRC_WRP: + switch (class) { + case MAG_CLASS_PSF: + Nphot = secfilt[Ns].NusedWrp; + break; + case MAG_CLASS_KRON: + Nphot = secfilt[Ns].NusedKronWrp; + break; + case MAG_CLASS_APER: + Nphot = secfilt[Ns].NusedApWrp; + break; + default: + break; + } + break; + case MAG_SRC_STK: + default: + break; + } + return (Nphot); } @@ -783,15 +863,15 @@ } -unsigned int PhotStackID (PhotCode *code, Average *average, SecFilt *secfilt) { - - int Ns; - unsigned int ID; - - if (code == NULL) return 0; - - Ns = photcodes[0].hashNsec[code[0].code]; - ID = (Ns == -1) ? 0 : secfilt[Ns].stackDetectID; - return (ID); -} +// XX unsigned int PhotStackID (PhotCode *code, Average *average, SecFilt *secfilt) { +// XX +// XX int Ns; +// XX unsigned int ID; +// XX +// XX if (code == NULL) return 0; +// XX +// XX Ns = photcodes[0].hashNsec[code[0].code]; +// XX ID = (Ns == -1) ? 0 : secfilt[Ns].stackDetectID; +// XX return (ID); +// XX } // Xm is now (2014.07.03) stored as the chisq except in dvo formats which use as short @@ -1112,4 +1192,64 @@ } return (Fave); +} + +float PhotFluxAveErr (PhotCode *code, Average *average, SecFilt *secfilt, dvoMagClassType class, dvoMagSourceType source) { + + if (code == NULL) return NAN; + + int Ns = photcodes[0].hashNsec[code[0].code]; + if (Ns == -1) return (NAN); + + float dFave = NAN; + switch (source) { + case MAG_SRC_CHP: + switch (class) { + case MAG_CLASS_PSF: + dFave = secfilt[Ns].dM * MagToFlux(secfilt[Ns].M); + break; + case MAG_CLASS_KRON: + dFave = secfilt[Ns].dMkron * MagToFlux(secfilt[Ns].Mkron); + break; + case MAG_CLASS_APER: + dFave = secfilt[Ns].dMap * MagToFlux(secfilt[Ns].Map); + break; + default: + break; + } + break; + case MAG_SRC_WRP: + switch (class) { + case MAG_CLASS_PSF: + dFave = secfilt[Ns].dFpsfWrp; + break; + case MAG_CLASS_KRON: + dFave = secfilt[Ns].dFkronWrp; + break; + case MAG_CLASS_APER: + dFave = secfilt[Ns].dFapWrp; + break; + default: + break; + } + break; + case MAG_SRC_STK: + switch (class) { + case MAG_CLASS_PSF: + dFave = secfilt[Ns].dFpsfStk; + break; + case MAG_CLASS_KRON: + dFave = secfilt[Ns].dFkronStk; + break; + case MAG_CLASS_APER: + dFave = secfilt[Ns].dFapStk; + break; + default: + break; + } + break; + default: + break; + } + return (dFave); } Index: /trunk/Ohana/src/opihi/dvo/gstar.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/gstar.c (revision 37115) +++ /trunk/Ohana/src/opihi/dvo/gstar.c (revision 37116) @@ -23,4 +23,24 @@ GSTAR_STACK_FLUX_KRON, GSTAR_STACK_FLUX_KRON_ERR, + GSTAR_STACK_FLUX_APER, + GSTAR_STACK_FLUX_APER_ERR, + GSTAR_WARP_FLUX_PSF, + GSTAR_WARP_FLUX_PSF_ERR, + GSTAR_WARP_FLUX_KRON, + GSTAR_WARP_FLUX_KRON_ERR, + GSTAR_WARP_FLUX_APER, + GSTAR_WARP_FLUX_APER_ERR, + GSTAR_STACK_MAG_PSF, + GSTAR_STACK_MAG_PSF_ERR, + GSTAR_STACK_MAG_KRON, + GSTAR_STACK_MAG_KRON_ERR, + GSTAR_STACK_MAG_APER, + GSTAR_STACK_MAG_APER_ERR, + GSTAR_WARP_MAG_PSF, + GSTAR_WARP_MAG_PSF_ERR, + GSTAR_WARP_MAG_KRON, + GSTAR_WARP_MAG_KRON_ERR, + GSTAR_WARP_MAG_APER, + GSTAR_WARP_MAG_APER_ERR, } GSTAR_SECF_CODES; @@ -31,8 +51,8 @@ void sort_index_float (float *X, int *IDX, int N) { -# define SWAPFUNC(A,B){ float tmp; int itmp; \ - tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ - itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ -} +# define SWAPFUNC(A,B){ float tmp; int itmp; \ + tmp = X[A]; X[A] = X[B]; X[B] = tmp; \ + itmp = IDX[A]; IDX[A] = IDX[B]; IDX[B] = itmp; \ + } # define COMPARE(A,B)(X[A] < X[B]) @@ -52,5 +72,5 @@ off_t i, Nstars, *N1; off_t j, k, m, N, Nlo, Nhi; - int Nsecfilt, NPTS, QUIET, FULL_OUTPUT, STACK_OUTPUT, INST, SHOW_MASKS; + int Nsecfilt, NPTS, QUIET, INST, SHOW_MASKS; int found, GetMeasures; int SaveVectors; @@ -109,5 +129,5 @@ } - FULL_OUTPUT = FALSE; + int FULL_OUTPUT = FALSE; if ((N = get_argument (argc, argv, "-full"))) { FULL_OUTPUT = TRUE; @@ -115,7 +135,43 @@ } - STACK_OUTPUT = FALSE; - if ((N = get_argument (argc, argv, "-stack"))) { + int STACK_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+stack"))) { STACK_OUTPUT = TRUE; + remove_argument (N, &argc, argv); + } + + int EXTRA_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+extras"))) { + EXTRA_OUTPUT = TRUE; + remove_argument (N, &argc, argv); + } + + int WARP_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+warp"))) { + WARP_OUTPUT = TRUE; + remove_argument (N, &argc, argv); + } + + int FLUX_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+flux"))) { + FLUX_OUTPUT = TRUE; + remove_argument (N, &argc, argv); + } + + int MAG_OUTPUT = TRUE; + if ((N = get_argument (argc, argv, "-mag"))) { + MAG_OUTPUT = FALSE; + remove_argument (N, &argc, argv); + } + + int KRON_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+kron"))) { + KRON_OUTPUT = TRUE; + remove_argument (N, &argc, argv); + } + + int APER_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "+aper"))) { + APER_OUTPUT = TRUE; remove_argument (N, &argc, argv); } @@ -185,6 +241,6 @@ // an error exit status here is a significant error if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) { - fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); - exit (2); + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); } dvo_catalog_unlock (&catalog); @@ -248,5 +304,5 @@ k = N1[i]; if (!QUIET) { - gprint (GP_LOG, "star: "OFF_T_FMT"\n", i); + gprint (GP_LOG, "star: "OFF_T_FMT"\n", k); gprint (GP_LOG, "%11.7f ", catalog.average[k].R); gprint (GP_LOG, "%11.7f ", catalog.average[k].D); @@ -259,88 +315,228 @@ if (FULL_OUTPUT) { - gprint (GP_LOG, "%f ", catalog.average[k].dR); - gprint (GP_LOG, "%f ", catalog.average[k].dD); - gprint (GP_LOG, "%f ", catalog.average[k].uR); - gprint (GP_LOG, "%f ", catalog.average[k].uD); - gprint (GP_LOG, "%f ", catalog.average[k].duR); - gprint (GP_LOG, "%f ", catalog.average[k].duD); - gprint (GP_LOG, "%f ", catalog.average[k].P); - gprint (GP_LOG, "%f ", catalog.average[k].dP); - - gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPM); - gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPar); - - date = ohana_sec_to_date (catalog.average[k].Tmean); - gprint (GP_LOG, "%20s ", date); - gprint (GP_LOG, "%f ", catalog.average[k].Trange / 86400.0); + gprint (GP_LOG, "%f ", catalog.average[k].dR); + gprint (GP_LOG, "%f ", catalog.average[k].dD); + gprint (GP_LOG, "%f ", catalog.average[k].uR); + gprint (GP_LOG, "%f ", catalog.average[k].uD); + gprint (GP_LOG, "%f ", catalog.average[k].duR); + gprint (GP_LOG, "%f ", catalog.average[k].duD); + gprint (GP_LOG, "%f ", catalog.average[k].P); + gprint (GP_LOG, "%f ", catalog.average[k].dP); + + gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPM); + gprint (GP_LOG, "%f ", catalog.average[k].ChiSqPar); + + date = ohana_sec_to_date (catalog.average[k].Tmean); + gprint (GP_LOG, "%20s ", date); + gprint (GP_LOG, "%f ", catalog.average[k].Trange / 86400.0); } - gprint (GP_LOG, "\n"); + gprint (GP_LOG, "\n\n"); /* filter names */ + gprint (GP_LOG, "filter : "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_FILTER_NAMES); gprint (GP_LOG, "\n"); /* average mags */ + gprint (GP_LOG, "chp_psf_ave: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG); gprint (GP_LOG, "\n"); /* average mag errors */ + gprint (GP_LOG, "chp_psf_err: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_ERR); gprint (GP_LOG, "\n"); /* average mag chisq */ + gprint (GP_LOG, "chp_psf_chi: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_CHISQ); gprint (GP_LOG, "\n"); - if (FULL_OUTPUT) { + + if (APER_OUTPUT) { /* Map */ + gprint (GP_LOG, "chp_ap__ave: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_AP_MAG); gprint (GP_LOG, "\n"); - + } + + if (KRON_OUTPUT) { /* Mkron */ + gprint (GP_LOG, "chp_krn_ave: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG); gprint (GP_LOG, "\n"); /* dMkron */ + gprint (GP_LOG, "chp_krn_err: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_KRON_MAG_ERR); gprint (GP_LOG, "\n"); - + } + + if (EXTRA_OUTPUT) { /* Mmin */ + gprint (GP_LOG, "chp_psf_min: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_MIN); gprint (GP_LOG, "\n"); /* Mmax */ + gprint (GP_LOG, "chp_psf_max: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_AVE_MAG_MAX);; gprint (GP_LOG, "\n"); /* UCDIST */ + gprint (GP_LOG, "ubercal_dst: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_UCDIST);; gprint (GP_LOG, "\n"); /* secfilt flags */ + gprint (GP_LOG, "filtflags: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_SECF_FLAGS); gprint (GP_LOG, "\n"); - } - if (STACK_OUTPUT) { + if (STACK_OUTPUT && MAG_OUTPUT) { /* FluxPSF */ + gprint (GP_LOG, "stk_psf_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_PSF); + gprint (GP_LOG, "\n"); + + /* dFluxPSF */ + gprint (GP_LOG, "stk_psf_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_PSF_ERR); + gprint (GP_LOG, "\n"); + + if (KRON_OUTPUT) { + /* MagKron */ + gprint (GP_LOG, "stk_krn_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_KRON); + gprint (GP_LOG, "\n"); + + /* dMagKron */ + gprint (GP_LOG, "stk_krn_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_KRON_ERR); + gprint (GP_LOG, "\n"); + } + + if (APER_OUTPUT) { + /* MagAper */ + gprint (GP_LOG, "stk_ap__ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_APER); + gprint (GP_LOG, "\n"); + + /* dMagAper */ + gprint (GP_LOG, "stk_ap__err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_MAG_APER_ERR); + gprint (GP_LOG, "\n"); + } + } + if (WARP_OUTPUT && MAG_OUTPUT) { + /* FluxPSF */ + gprint (GP_LOG, "wrp_psf_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_PSF); + gprint (GP_LOG, "\n"); + + /* dFluxPSF */ + gprint (GP_LOG, "wrp_psf_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_PSF_ERR); + gprint (GP_LOG, "\n"); + + if (KRON_OUTPUT) { + /* MagKron */ + gprint (GP_LOG, "wrp_krn_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_KRON); + gprint (GP_LOG, "\n"); + + /* dMagKron */ + gprint (GP_LOG, "wrp_krn_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_KRON_ERR); + gprint (GP_LOG, "\n"); + + } + + if (APER_OUTPUT) { + /* MagAper */ + gprint (GP_LOG, "wrp_ap__ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_APER); + gprint (GP_LOG, "\n"); + + /* dMagAper */ + gprint (GP_LOG, "wrp_ap__err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_MAG_APER_ERR); + gprint (GP_LOG, "\n"); + } + } + if (STACK_OUTPUT && FLUX_OUTPUT) { + /* FluxPSF */ + gprint (GP_LOG, "stk_psf_ave: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF); gprint (GP_LOG, "\n"); - + /* dFluxPSF */ + gprint (GP_LOG, "stk_psf_err: "); for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_PSF_ERR); gprint (GP_LOG, "\n"); - - /* FluxKron */ - for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON); - gprint (GP_LOG, "\n"); - - /* dFluxKron */ - for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON_ERR); - gprint (GP_LOG, "\n"); - + + if (KRON_OUTPUT) { + /* FluxKron */ + gprint (GP_LOG, "stk_krn_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON); + gprint (GP_LOG, "\n"); + + /* dFluxKron */ + gprint (GP_LOG, "stk_krn_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_KRON_ERR); + gprint (GP_LOG, "\n"); + } + + if (APER_OUTPUT) { + /* FluxAper */ + gprint (GP_LOG, "stk_ap_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_APER); + gprint (GP_LOG, "\n"); + + /* dFluxAper */ + gprint (GP_LOG, "stk_ap_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_STACK_FLUX_APER_ERR); + gprint (GP_LOG, "\n"); + } } - } + if (WARP_OUTPUT && FLUX_OUTPUT) { + /* FluxPSF */ + gprint (GP_LOG, "wrp_psf_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_PSF); + gprint (GP_LOG, "\n"); + + /* dFluxPSF */ + gprint (GP_LOG, "wrp_psf_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_PSF_ERR); + gprint (GP_LOG, "\n"); + + if (KRON_OUTPUT) { + /* FluxKron */ + gprint (GP_LOG, "wrp_krn_ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_KRON); + gprint (GP_LOG, "\n"); + + /* dFluxKron */ + gprint (GP_LOG, "wrp_krn_err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_KRON_ERR); + gprint (GP_LOG, "\n"); + } + + if (APER_OUTPUT) { + /* FluxAper */ + gprint (GP_LOG, "wrp_ap__ave: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_APER); + gprint (GP_LOG, "\n"); + + /* dFluxAper */ + gprint (GP_LOG, "wrp_ap__err: "); + for (j = 0; j < Nsecfilt; j++) printPhotcodeSequence (&catalog.average[k], &catalog.secfilt[Nsecfilt*k], j, GSTAR_WARP_FLUX_APER_ERR); + gprint (GP_LOG, "\n"); + } + } + } + gprint (GP_LOG, "\n"); + if (GetMeasures || SaveVectors) { @@ -413,31 +609,31 @@ if (FULL_OUTPUT) { - gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mcal); - gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Map); - gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mkron); - gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].dMkron); - gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[Nv].dt)); - gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].airmass); - gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].az); - gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Xccd); - gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Yccd); - gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dXccd)); - gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dYccd)); - gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Sky); - gprint (GP_LOG, "%5.1f ", catalog.measure[Nv].dSky); - gprint (GP_LOG, "%8d ", catalog.measure[Nv].averef); - gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].detID); - gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].imageID); - gprint (GP_LOG, "%.3f ", catalog.measure[Nv].psfQF); - gprint (GP_LOG, "%.1f ", catalog.measure[Nv].psfChisq); - gprint (GP_LOG, "%.1f ", catalog.measure[Nv].crNsigma); - gprint (GP_LOG, "%.1f ", catalog.measure[Nv].extNsigma); - gprint (GP_LOG, "%.1f ", FromShortDegrees(catalog.measure[Nv].theta)); + gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mcal); + gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Map); + gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].Mkron); + gprint (GP_LOG, "%6.3f ", catalog.measure[Nv].dMkron); + gprint (GP_LOG, "%5.1f ", pow(10.0, 0.4*catalog.measure[Nv].dt)); + gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].airmass); + gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].az); + gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Xccd); + gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Yccd); + gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dXccd)); + gprint (GP_LOG, "%3.1f ", FromShortPixels(catalog.measure[Nv].dYccd)); + gprint (GP_LOG, "%6.1f ", catalog.measure[Nv].Sky); + gprint (GP_LOG, "%5.1f ", catalog.measure[Nv].dSky); + gprint (GP_LOG, "%8d ", catalog.measure[Nv].averef); + gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].detID); + gprint (GP_LOG, "0x%08x ", catalog.measure[Nv].imageID); + gprint (GP_LOG, "%5.3f ", catalog.measure[Nv].psfQF); + gprint (GP_LOG, "%7.1f ", catalog.measure[Nv].psfChisq); + gprint (GP_LOG, "%3.1f ", catalog.measure[Nv].crNsigma); + gprint (GP_LOG, "%4.1f ", catalog.measure[Nv].extNsigma); + gprint (GP_LOG, "%5.1f ", FromShortDegrees(catalog.measure[Nv].theta)); } - if (STACK_OUTPUT) { - gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxPSF); - gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxPSF); - gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxKron); - gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxKron); + if (FLUX_OUTPUT) { + gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxPSF); + gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxPSF); + gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].FluxKron); + gprint (GP_LOG, "%10.3e ", catalog.measure[Nv].dFluxKron); } gprint (GP_LOG, "\n"); @@ -547,17 +743,17 @@ switch (type) { - case GSTAR_AVE_MAG: /* average mags */ - if (seq == -1) { - print_double (NAN); - } else { - print_double (secfilt[seq].M); - } - break; - - case GSTAR_AVE_MAG_ERR: /* average mags errors */ - if (seq == -1) { - print_double (NAN); - } else { - print_double (secfilt[seq].dM); + case GSTAR_AVE_MAG_MIN: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].Mmin); + } + break; + + case GSTAR_AVE_MAG_MAX: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].Mmax); } break; @@ -592,4 +788,22 @@ break; + /*** CHIP : MAG ***/ + + case GSTAR_AVE_MAG: /* average mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].M); + } + break; + + case GSTAR_AVE_MAG_ERR: /* average mags errors */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dM); + } + break; + case GSTAR_AVE_AP_MAG: /* average ap mags */ if (seq == -1) { @@ -600,20 +814,4 @@ break; - case GSTAR_AVE_MAG_MIN: /* average ap mags */ - if (seq == -1) { - print_double (NAN); - } else { - print_double (secfilt[seq].Mmin); - } - break; - - case GSTAR_AVE_MAG_MAX: /* average ap mags */ - if (seq == -1) { - print_double (NAN); - } else { - print_double (secfilt[seq].Mmax); - } - break; - case GSTAR_AVE_KRON_MAG: /* average ap mags */ if (seq == -1) { @@ -632,4 +830,6 @@ break; + /******************* STACK : FLUX **********************/ + case GSTAR_STACK_FLUX_PSF: /* average ap mags */ if (seq == -1) { @@ -664,4 +864,170 @@ break; + case GSTAR_STACK_FLUX_APER: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].FapStk); + } + break; + + case GSTAR_STACK_FLUX_APER_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].dFapStk); + } + break; + + /******************* STACK : MAG **********************/ + + case GSTAR_STACK_MAG_PSF: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MpsfStk); + } + break; + + case GSTAR_STACK_MAG_PSF_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFpsfStk / secfilt[seq].FpsfStk); + } + break; + + case GSTAR_STACK_MAG_KRON: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MkronStk); + } + break; + + case GSTAR_STACK_MAG_KRON_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFkronStk / secfilt[seq].FkronStk); + } + break; + + case GSTAR_STACK_MAG_APER: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MapStk); + } + break; + + case GSTAR_STACK_MAG_APER_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFapStk / secfilt[seq].FapStk); + } + break; + + /******************* WARP : FLUX **********************/ + + case GSTAR_WARP_FLUX_PSF: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].FpsfWrp); + } + break; + + case GSTAR_WARP_FLUX_PSF_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].dFpsfWrp); + } + break; + + case GSTAR_WARP_FLUX_KRON: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].FkronWrp); + } + break; + + case GSTAR_WARP_FLUX_KRON_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].dFkronWrp); + } + break; + + case GSTAR_WARP_FLUX_APER: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].FapWrp); + } + break; + + case GSTAR_WARP_FLUX_APER_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double_exp (secfilt[seq].dFapWrp); + } + break; + + /******************* WARP : MAG **********************/ + + case GSTAR_WARP_MAG_PSF: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MpsfWrp); + } + break; + + case GSTAR_WARP_MAG_PSF_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFpsfWrp / secfilt[seq].FpsfWrp); + } + break; + + case GSTAR_WARP_MAG_KRON: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MkronWrp); + } + break; + + case GSTAR_WARP_MAG_KRON_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFkronWrp / secfilt[seq].FkronWrp); + } + break; + + case GSTAR_WARP_MAG_APER: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].MapWrp); + } + break; + + case GSTAR_WARP_MAG_APER_ERR: /* average ap mags */ + if (seq == -1) { + print_double (NAN); + } else { + print_double (secfilt[seq].dFapWrp / secfilt[seq].FapWrp); + } + break; + } } Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 37115) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 37116) @@ -440,4 +440,5 @@ int MatchImageName (off_t meas, int cat, char *name); int MatchImageSkycellID (off_t meas, int cat, int myTessID, int myProjectionID, int mySkycellID); +int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID); int load_tess (char *treefile); @@ -451,5 +452,5 @@ int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt); int setMrelAverageExposure (off_t meas, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found); -int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell); +int setMrelAverageStack (off_t meas, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure); int setGlobalObjStats (Average *average, Measure *measure); Index: /trunk/Ohana/src/relphot/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/ImageOps.c (revision 37115) +++ /trunk/Ohana/src/relphot/src/ImageOps.c (revision 37116) @@ -491,4 +491,27 @@ if (projectID[i] != myProjectionID) return FALSE; if (skycellID[i] != mySkycellID) return FALSE; + + return TRUE; +} + +// returns image.Mcal - ff(x,y) +int FindImageSkycellID (off_t meas, int cat, int *myTessID, int *myProjectionID, int *mySkycellID) { + + off_t i; + + if (!MeasureToImage) return FALSE; + + if (meas < 0) return FALSE; + + i = MeasureToImage[cat][meas]; + if (i == -1) return FALSE; + + if (tessID[i] == -1) return FALSE; + if (projectID[i] == -1) return FALSE; + if (skycellID[i] == -1) return FALSE; + + *myTessID = tessID[i]; + *myProjectionID = projectID[i]; + *mySkycellID = skycellID[i]; return TRUE; Index: /trunk/Ohana/src/relphot/src/StarOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/StarOps.c (revision 37115) +++ /trunk/Ohana/src/relphot/src/StarOps.c (revision 37116) @@ -1086,4 +1086,10 @@ if (ecode != thisCode) { continue; } + // SKIP gpc1 stack data + if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); @@ -1131,4 +1137,10 @@ if (ecode != thisCode) { continue; } + // SKIP gpc1 stack data + if (isGPC1stack(catalog[i].measureT[m].photcode)) continue; + + // SKIP gpc1 forced-warp data + if (isGPC1warp(catalog[i].measureT[m].photcode)) continue; + // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers Mcal = getMcal (m, i, flatcorr, catalog); Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 37115) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 37116) @@ -78,5 +78,5 @@ int Nsec = GetPhotcodeNsec(thisCode); - fprintf (stderr, "really use dvo_secfilt_init?"); + // fprintf (stderr, "really use dvo_secfilt_init?"); dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]); Index: /trunk/Ohana/src/relphot/src/setMrelCatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 37115) +++ /trunk/Ohana/src/relphot/src/setMrelCatalog.c (revision 37116) @@ -50,8 +50,8 @@ int isSetMrelFinal = (pass >= 0); - char *primaryCell = NULL; - if (isSetMrelFinal) { - ALLOCATE (primaryCell, char, DVO_MAX_PATH); - } + // XX char *primaryCell = NULL; + // XX if (isSetMrelFinal) { + // XX ALLOCATE (primaryCell, char, DVO_MAX_PATH); + // XX } for (j = 0; j < catalog[Nc].Naverage; j++) { @@ -64,5 +64,5 @@ // only apply Stack operation on setMrelFinal in first pass if (isSetMrelFinal && (pass == 0)) { - setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell); + setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m]); setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]); } @@ -74,5 +74,5 @@ } } - if (primaryCell) free (primaryCell); + // if (primaryCell) free (primaryCell); return (TRUE); } @@ -441,9 +441,13 @@ // only apply Stack operation on setMrelFinal in first pass -int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) { - - off_t k, ID; - - float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0; +// this function has 3 goals, not to be confused: +// 1) find and mark the PRIMARY detections (regardless of the quality of the detection) +// 2) select the BEST detections per filter (regardless of PRIMARY) +// 3) apply the zero point and AB->Jy transformations +int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure) { + + off_t k; + + float Mcal= 0, Mmos = 0, Mgrid = 0, Finst = 0; // set the primary projection cell and skycell for this coordinate @@ -467,12 +471,8 @@ int haveStack = FALSE; - // need to find the measurement closest to the center of its skycell, as well as the - // closest for the subset of primary projection cells - - float stackCenterOffsetMin = 1e9; - off_t stackCenterMeasureMin = -1; - - float stackPrimaryOffsetMin = 1e9; - off_t stackPrimaryMeasureMin = -1; + float psfQFbest = 0.0; + + off_t stackBestMeasure = -1; + off_t stackPrimaryMeasure = -1; int isBad = FALSE; @@ -481,4 +481,7 @@ off_t meas = measureOffset; for (k = 0; k < average[0].Nmeasure; k++, meas++) { + + // only examine gpc1 stack data + if (!isGPC1stack(measure[k].photcode)) continue; // skip measurements that do not match the current photcode @@ -489,14 +492,30 @@ // clear this bit for all measurements measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY; - - // only examine gpc1 stack data - if (!isGPC1stack(measure[k].photcode)) continue; + measure[k].dbFlags &= ~ID_MEAS_STACK_PHOT_SRC; haveStack = TRUE; haveStackObject = TRUE; + // match measurement to its image (this is just a check, right?) + if (getImageEntry (meas, cat) < 0) { + fprintf (stderr, "*** warning : this should be impossible\n"); + abort(); + continue; + } + + // find the PRIMARY measurement + + // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the + // primary skycell for this position + if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) { + stackPrimaryMeasure = k; + measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + secfilt[Nsec].stackPrmryOff = meas; + } + + // now choose the BEST measurements (may also be PRIMARY) + if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); - if (getImageEntry (meas, cat) < 0) continue; // measurements without an image are either external reference photometry or // data for which the associated image has not been loaded (probably because of @@ -510,26 +529,11 @@ if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid); - Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF); - if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys); - - unsigned int stackImageID; - - // which stack image should we use for the mean value? - // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell - if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) { - float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); - if (stackPrimaryOffset < stackPrimaryOffsetMin) { - stackPrimaryOffsetMin = stackPrimaryOffset; - // stackPrimaryIDmin = stackImageID; - stackPrimaryMeasureMin = k; - } - } - - // get the center distance for the generic case: - float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID); - if (stackCenterOffset < stackCenterOffsetMin) { - stackCenterOffsetMin = stackCenterOffset; - // stackCenterIDmin = stackImageID; - stackCenterMeasureMin = k; + // NOTE: negative and insignificant vlues values are allowed, but not NAN flux values + Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF); + if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst); + + if (measure[k].psfQFperf > psfQFbest) { + psfQFbest = measure[k].psfQFperf; + stackBestMeasure = k; } @@ -554,11 +558,10 @@ // measurements which are bad will not have a valid stack entry and are skipped - k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin; + k = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure; if (k < 0) continue; - // we are now populating stackDetectID not stack Image ID in secfilt // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin : stackCenterIDmin; - ID = measure[k].extID; // for the stack, this is the stackDetectID + // ID = measure[k].extID; // for the stack, this is the stackDetectID // get the zero point for the selected image @@ -605,13 +608,14 @@ secfilt[Nsec].dFapStk = zpFactor * measure[k].dFluxAp; + // Jy to AB mags secfilt[Nsec].MpsfStk = (measure[k].FluxPSF > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN; secfilt[Nsec].MkronStk = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN; secfilt[Nsec].MapStk = (measure[k].FluxAp > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN; - secfilt[Nsec].stackDetectID = ID; - + secfilt[Nsec].stackBestOff = k + measureOffset; + + // this is the measurement used by secfilt[] measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC; - if (stackPrimaryMeasureMin >= 0) { - measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY; + if (k == stackPrimaryMeasure) { secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY; } @@ -642,5 +646,5 @@ off_t k; - float Fsys = 0, Mcal= 0; + float Mcal= 0; // we are measuring means for 3 types of FLUXes: psf, ap, kron. I am using the psf mag @@ -664,5 +668,5 @@ // option for a test print - if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) { + if ((average[0].objID == 0xd6e) && (average[0].catID == 0x4984)) { fprintf (stderr, "test obj\n"); print_measure_set_alt (average, secfilt, measure); @@ -680,4 +684,9 @@ if (found[Nsec]) continue; + // I need to restrict the measurements coming only from the same skycell as the secfilt values + // if we cannot make that association, skip this set of warps + int stkTessID, stkProjID, stkSkycellID; + if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &stkTessID, &stkProjID, &stkSkycellID)) continue; + off_t meas = measureOffset; @@ -694,4 +703,8 @@ if (!code) continue; if (code->equiv != thisCode) { continue; } + + measure[k].dbFlags &= ~ID_MEAS_WARP_USED; + + if (!MatchImageSkycellID (meas, cat, stkTessID, stkProjID, stkSkycellID)) continue; if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); @@ -720,9 +733,18 @@ // are small (and not measured) float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF); - if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys); + if (isnan(Fpsf)) SKIP_THIS_MEAS(Nsys); // if (Msys < 0.0) SKIP_THIS_MEAS(Nsys); // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys); float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF); + if (isnan(dFpsf)) SKIP_THIS_MEAS (Nsys); + + int isBad = (measure[k].psfQF < 0.85); + isBad |= isnan(measure[k].psfQF); + if (isBad) SKIP_THIS_MEAS (Nsys); + + int isPoor = (measure[k].psfQFperf < 0.85); + if ((pass == 0) && isPoor) SKIP_THIS_MEAS (Nsys); + dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error Fpsflist[Npsf] = Fpsf * Fcal; @@ -730,4 +752,6 @@ wpsflist[Npsf] = 1.0; Npsf ++; + + measure[k].dbFlags |= ID_MEAS_WARP_USED; float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);