Index: /trunk/Ohana/src/addstar/src/GetFileMode.c =================================================================== --- /trunk/Ohana/src/addstar/src/GetFileMode.c (revision 35415) +++ /trunk/Ohana/src/addstar/src/GetFileMode.c (revision 35416) @@ -20,5 +20,8 @@ havePHOT_VER = gfits_scan (header, "PHOT_VER", "%s", 1, ctmp); haveTARG_VER = gfits_scan (header, "TARG_VER", "%s", 1, ctmp); - if (havePHOT_VER && haveTARG_VER) return SDSS_OBJ; + if (havePHOT_VER && haveTARG_VER) { + if (VERBOSE) fprintf (stderr, "found SDSS objects\n"); + return SDSS_OBJ; + } if (haveNaxis && (Naxis == 2)) { @@ -28,6 +31,8 @@ if ((Nx > 0) && (Ny > 0)) { if (haveCTYPE && !strcmp (&ctype[4], "-WRP")) { + if (VERBOSE) fprintf (stderr, "found MOSAIC CMP\n"); return MOSAIC_CMP; } + if (VERBOSE) fprintf (stderr, "found SIMPLE CMP\n"); return SIMPLE_CMP; } @@ -36,6 +41,8 @@ if (haveNaxis && (TEXTMODE || !simple)) { if (haveCTYPE && !strcmp (&ctype[4], "-WRP")) { + if (VERBOSE) fprintf (stderr, "found MOSAIC CMP\n"); return MOSAIC_CMP; } + if (VERBOSE) fprintf (stderr, "found SIMPLE CMP\n"); return SIMPLE_CMP; } @@ -43,21 +50,27 @@ if (!extend && strcmp (&ctype[4], "-DIS")) { if (!strcmp (&ctype[4], "-WRP")) { + if (VERBOSE) fprintf (stderr, "found MOSAIC CMF\n"); return MOSAIC_CMF; } + if (VERBOSE) fprintf (stderr, "found SIMPLE CMF\n"); return SIMPLE_CMF; } if (!extend && !strcmp (&ctype[4], "-DIS")) { + if (VERBOSE) fprintf (stderr, "found MOSAIC PHU\n"); return MOSAIC_PHU; } if (extend && strcmp (&ctype[4], "-DIS")) { + if (VERBOSE) fprintf (stderr, "found SIMPLE MEF\n"); return SIMPLE_MEF; } if (extend && !strcmp (&ctype[4], "-DIS")) { + if (VERBOSE) fprintf (stderr, "found MOSAIC MEF\n"); return MOSAIC_MEF; } + if (VERBOSE) fprintf (stderr, "extension type is unknown\n"); return (NONE); } Index: /trunk/Ohana/src/addstar/src/MatchHeaders.c =================================================================== --- /trunk/Ohana/src/addstar/src/MatchHeaders.c (revision 35415) +++ /trunk/Ohana/src/addstar/src/MatchHeaders.c (revision 35416) @@ -61,4 +61,7 @@ if (!strcmp (exttype, "PS1_SV1")) goto keep; if (!strcmp (exttype, "PS1_SV2")) goto keep; + if (!strcmp (exttype, "PS1_DV3")) { + goto keep; + } continue; Index: /trunk/Ohana/src/addstar/src/ReadStarsFITS.c =================================================================== --- /trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 35415) +++ /trunk/Ohana/src/addstar/src/ReadStarsFITS.c (revision 35416) @@ -15,4 +15,5 @@ Stars *Convert_PS1_SV1_Alt PROTO((FTable *table, unsigned int *nstars)); Stars *Convert_PS1_SV2 PROTO((FTable *table, unsigned int *nstars)); +Stars *Convert_PS1_DV3 PROTO((FTable *table, unsigned int *nstars)); // given a file with the pointer at the start of the table block and the @@ -78,4 +79,7 @@ stars = Convert_PS1_SV2 (&table, &Nstars); } + if (!strcmp (type, "PS1_DV3")) { + stars = Convert_PS1_DV3 (&table, &Nstars); + } if (stars == NULL) { fprintf (stderr, "invalid table type %s\n", type); @@ -962,3 +966,85 @@ } - +Stars *Convert_PS1_DV3 (FTable *table, unsigned int *nstars) { + + off_t Nstars; + unsigned int i; + double ZeroPt; + Stars *stars; + CMF_PS1_DV3 *ps1data; + + ps1data = gfits_table_get_CMF_PS1_DV3 (table, &Nstars, NULL); + if (!ps1data) { + fprintf (stderr, "skipping inconsistent entry\n"); + return (NULL); + } + ZeroPt = GetZeroPoint(); + + ALLOCATE (stars, Stars, Nstars); + for (i = 0; i < Nstars; i++) { + InitStar (&stars[i]); + stars[i].measure.Xccd = ps1data[i].X; + stars[i].measure.Yccd = ps1data[i].Y; + stars[i].measure.dXccd = ToShortPixels(ps1data[i].dX); + stars[i].measure.dYccd = ToShortPixels(ps1data[i].dY); + + stars[i].measure.posangle = ToShortDegrees(ps1data[i].posangle); + stars[i].measure.pltscale = ps1data[i].pltscale; + + if ((ps1data[i].M >= 0.0) || isnan(ps1data[i].M)) { + stars[i].measure.M = NAN; + } else { + stars[i].measure.M = ps1data[i].M + ZeroPt; + } + stars[i].measure.dM = ps1data[i].dM; + stars[i].measure.dMcal = ps1data[i].dMcal; + stars[i].measure.Map = ps1data[i].Map + ZeroPt; + + stars[i].measure.Mkron = (ps1data[i].kronFlux > 0.0) ? -2.5*log10(ps1data[i].kronFlux) + ZeroPt : NAN; + stars[i].measure.dMkron = (ps1data[i].kronFlux > 0.0) ? ps1data[i].kronFluxErr / ps1data[i].kronFlux : NAN; + + // these fluxes are converted from counts to counts/sec in FilterStars.c + stars[i].measure.FluxPSF = ps1data[i].Flux; + stars[i].measure.dFluxPSF = ps1data[i].dFlux; + stars[i].measure.FluxKron = ps1data[i].kronFlux; + stars[i].measure.dFluxKron = ps1data[i].kronFluxErr; + + stars[i].measure.Sky = ps1data[i].sky; + stars[i].measure.dSky = ps1data[i].dSky; + + stars[i].measure.psfChisq = ps1data[i].psfChisq; + stars[i].measure.psfQF = ps1data[i].psfQF; + stars[i].measure.psfQFperf = ps1data[i].psfQFperf; + stars[i].measure.psfNdof = ps1data[i].psfNdof; + stars[i].measure.psfNpix = ps1data[i].psfNpix; + stars[i].measure.crNsigma = ps1data[i].crNsigma; + stars[i].measure.extNsigma = ps1data[i].extNsigma; + + stars[i].measure.FWx = ToShortPixels(ps1data[i].fx); + stars[i].measure.FWy = ToShortPixels(ps1data[i].fy); + stars[i].measure.theta = ToShortDegrees(ps1data[i].df); + + stars[i].measure.Mxx = ToShortPixels(ps1data[i].Mxx); + stars[i].measure.Mxy = ToShortPixels(ps1data[i].Mxy); + stars[i].measure.Myy = ToShortPixels(ps1data[i].Myy); + + stars[i].measure.photFlags = ps1data[i].flags; + + // this is may optionally be replaced by the internal sequence (see FilterStars.c) + stars[i].measure.detID = ps1data[i].detID; + + // the Average fields and the following Measure fields are set in FilterStars after + // the image metadata is in hand: dR, dD, Mcal, dt, airmass, az, t, imageID, extID. + + // averef is set in find_matches + + // dbFlags is zero on ingest. + + // the following fields are currently not being set anywhere: t_msec + } + *nstars = Nstars; + return (stars); +} + + + Index: /trunk/Ohana/src/addstar/src/load2mass_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/load2mass_catalog.c (revision 35415) +++ /trunk/Ohana/src/addstar/src/load2mass_catalog.c (revision 35416) @@ -20,31 +20,11 @@ // XXX for now, the output objects will have limited astrometric interpretation... // XXX every 3 stars represents 3 measurements and 1 average + dvo_average_init (&catalog[0].average[Nave]); catalog[0].average[Nave].R = stars[i].average.R; catalog[0].average[Nave].D = stars[i].average.D; - catalog[0].average[Nave].dR = 0; - catalog[0].average[Nave].dD = 0; - catalog[0].average[Nave].uR = 0; - catalog[0].average[Nave].uD = 0; - catalog[0].average[Nave].duR = 0; - catalog[0].average[Nave].duD = 0; - catalog[0].average[Nave].P = 0; - catalog[0].average[Nave].dP = 0; - - // XXX for now, set the average mag data to NULL - catalog[0].average[Nave].Nmeasure = 0; - catalog[0].average[Nave].Nmissing = 0; - catalog[0].average[Nave].ChiSqAve = NAN_S_SHORT; catalog[0].average[Nave].measureOffset = Nmeas; - catalog[0].average[Nave].missingOffset = -1; - catalog[0].average[Nave].flags = 0; for (j = 0; j < Nsec; j++) { - catalog[0].secfilt[Nave*Nsec+j].M = NAN; - catalog[0].secfilt[Nave*Nsec+j].dM = NAN; - catalog[0].secfilt[Nave*Nsec+j].Xm = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsec+j].M_20 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsec+j].M_80 = NAN_S_SHORT; - catalog[0].secfilt[Nave*Nsec+j].Ncode = 0; - catalog[0].secfilt[Nave*Nsec+j].Nused = 0; + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsec+j]); } Index: /trunk/Ohana/src/delstar/Makefile =================================================================== --- /trunk/Ohana/src/delstar/Makefile (revision 35415) +++ /trunk/Ohana/src/delstar/Makefile (revision 35416) @@ -38,4 +38,5 @@ $(SRC)/find_matches.$(ARCH).o \ $(SRC)/parse_time.$(ARCH).o \ +$(SRC)/ImageIDsOps.$(ARCH).o \ $(SRC)/check_permissions.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o @@ -49,4 +50,5 @@ $(SRC)/delete_photcodes.$(ARCH).o \ $(SRC)/delete_photcodes_single.$(ARCH).o \ +$(SRC)/ImageIDsOps.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o Index: /trunk/Ohana/src/delstar/include/delstar.h =================================================================== --- /trunk/Ohana/src/delstar/include/delstar.h (revision 35415) +++ /trunk/Ohana/src/delstar/include/delstar.h (revision 35416) @@ -117,5 +117,5 @@ int delete_duplicate_images (int hostID, char *hostpath, FITS_DB *db); int delete_duplicate_image_measures (int hostID, char *hostpath, IndexArray *imageID); -int delete_duplicate_image_measures_parallel (SkyList *sky); +int delete_duplicate_image_measures_parallel (SkyList *sky, IndexArray *imageID); int delete_duplicate_image_measures_catalog (Catalog *catalog, IndexArray *imageID); @@ -123,2 +123,4 @@ IndexArray *make_index_array (Image *image, off_t Nimage, int mode); +IndexArray *ImageIDLoad(char *filename); +int ImageIDSave(char *filename, IndexArray *imageID); Index: /trunk/Ohana/src/delstar/src/ImageIDsOps.c =================================================================== --- /trunk/Ohana/src/delstar/src/ImageIDsOps.c (revision 35416) +++ /trunk/Ohana/src/delstar/src/ImageIDsOps.c (revision 35416) @@ -0,0 +1,133 @@ +# include "delstar.h" +# define BZERO_INT16 1.0*0x8000 +# define BZERO_INT32 1.0*0x80000000 +# define BZERO_INT64 1.0*0x8000000000000000 + +IndexArray *ImageIDLoad(char *filename) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + IndexArray *imageID = NULL; + ALLOCATE (imageID, IndexArray, 1); + + gfits_scan (&header, "MIN_ID", OFF_T_FMT, 1, &imageID->minID); + gfits_scan (&header, "MAX_ID", OFF_T_FMT, 1, &imageID->maxID); + gfits_scan (&header, "RANGE", OFF_T_FMT, 1, &imageID->range); + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return NULL; + } + fclose (f); + + char type[16]; + unsigned short *value = gfits_get_bintable_column_data (&theader, &ftable, "VALUE", type, &Nrow, &Ncol); + myAssert (!strcmp(type, "short"), "wrong column type"); + myAssert (Nrow == imageID->range, "wrong number of rows in table?"); + + // XXX free the fits table data here + + ALLOCATE (imageID->value, off_t, imageID->range); + for (i = 0; i < imageID->range; i++) { + imageID->value[i] = value[i]; + } + fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); + + free (value); + + return imageID; +} + +int ImageIDSave(char *filename, IndexArray *imageID) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_modify (&header, "MIN_ID", OFF_T_FMT, 1, imageID->minID); + gfits_modify (&header, "MAX_ID", OFF_T_FMT, 1, imageID->maxID); + gfits_modify (&header, "RANGE", OFF_T_FMT, 1, imageID->range); + + gfits_create_table_header (&theader, "BINTABLE", "IMAGE_ID"); + + gfits_define_bintable_column (&theader, "I", "VALUE", "value", NULL, 1.0, BZERO_INT16); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + unsigned short *value; + + // create intermediate storage arrays + ALLOCATE (value, unsigned short, imageID->range); + + // assign the storage arrays + for (i = 0; i < imageID->range; i++) { + value[i] = imageID->value[i]; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "VALUE", value, imageID->range); + + free (value); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image ID file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + fflush (f); + + return TRUE; +} Index: /trunk/Ohana/src/delstar/src/args.c =================================================================== --- /trunk/Ohana/src/delstar/src/args.c (revision 35415) +++ /trunk/Ohana/src/delstar/src/args.c (revision 35416) @@ -328,4 +328,9 @@ remove_argument (N, &argc, argv); } + if ((N = get_argument (argc, argv, "-dup-images"))) { + if (MODE != MODE_NONE) usage(); + MODE = MODE_DUP_IMAGES; + remove_argument (N, &argc, argv); + } if (MODE == MODE_NONE) delstar_client_usage (); Index: /trunk/Ohana/src/delstar/src/delete_duplicate_images.c =================================================================== --- /trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 35415) +++ /trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 35416) @@ -84,5 +84,5 @@ // launch the remote jobs if (PARALLEL && !hostID) { - int status = delete_duplicate_image_measures_parallel (skylist); + int status = delete_duplicate_image_measures_parallel (skylist, imageID); return status; } @@ -128,5 +128,14 @@ // CATDIR is supplied globally -int delete_duplicate_image_measures_parallel (SkyList *sky) { +int delete_duplicate_image_measures_parallel (SkyList *sky, IndexArray *imageID) { + + // write out the subset table of image information + char imageFile[512]; + snprintf (imageFile, 512, "%s/ImageIDs.tmp.fits", CATDIR); + + if (!ImageIDSave (imageFile, imageID)) { + fprintf (stderr, "failed to write image ID table\n"); + exit (1); + } // launch the delstar_client jobs to the parallel hosts @@ -157,6 +166,4 @@ fprintf (stderr, "command: %s\n", command); - - fprintf (stderr, "parallel version is not ready (need to save / load the imageID array\n"); if (PARALLEL_MANUAL) continue; @@ -463,11 +470,11 @@ // XXX this is probably bad : it assumes all chips // should be deleted for a given exposure - off_t parentID = image[i].parentID; - if (parentID) { - off_t Np = parentID - imageID->minID; - myAssert (Np >= 0, "oops"); - myAssert (Np < imageID->range, "oops"); - imageID->value[Np] = TRUE; - } + // off_t parentID = image[i].parentID; + // if (parentID) { + // off_t Np = parentID - imageID->minID; + // myAssert (Np >= 0, "oops"); + // myAssert (Np < imageID->range, "oops"); + // imageID->value[Np] = TRUE; + // } } Index: /trunk/Ohana/src/delstar/src/delstar_client.c =================================================================== --- /trunk/Ohana/src/delstar/src/delstar_client.c (revision 35415) +++ /trunk/Ohana/src/delstar/src/delstar_client.c (revision 35416) @@ -6,23 +6,23 @@ int main (int argc, char **argv) { + char imageFile[512]; + SetSignals (); ConfigInit (&argc, argv); args_client (argc, argv); - // load the current sky table (layout of all SkyRegions) - SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); - if (!sky) { - fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); - exit (2); - } - SkyTableSetFilenames (sky, CATDIR, "cpt"); - - SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - if (!skylist) { - fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); - exit (2); - } - switch (MODE) { + case MODE_DUP_IMAGES: + // read the subset table of image information + snprintf (imageFile, 512, "%s/ImageIDs.tmp.fits", CATDIR); + + IndexArray *imageID = ImageIDLoad (imageFile); + if (!imageID) { + fprintf (stderr, "failed to read image ID table\n"); + exit (1); + } + if (!delete_duplicate_image_measures (HOST_ID, HOSTDIR, imageID)) exit (1); + exit (0); + break; case MODE_IMAGEFILE: break; Index: /trunk/Ohana/src/dvomerge/Makefile =================================================================== --- /trunk/Ohana/src/dvomerge/Makefile (revision 35415) +++ /trunk/Ohana/src/dvomerge/Makefile (revision 35416) @@ -23,5 +23,5 @@ dvoverify_client : $(BIN)/dvoverify_client.$(ARCH) -all: dvomerge dvomerge_client dvoconvert dvosecfilt dvoverify dvoverify_client +all: dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvoverify dvoverify_client # $(SRC)/dvomergeContinue.$(ARCH).o @@ -93,4 +93,5 @@ DVOSECFILT = \ $(SRC)/dvosecfilt.$(ARCH).o \ +$(SRC)/dvosecfilt_catalogs.$(ARCH).o \ $(SRC)/SetSignals.$(ARCH).o \ $(SRC)/ConfigInit.$(ARCH).o \ @@ -101,4 +102,16 @@ $(DVOSECFILT) : $(INC)/dvomerge.h $(BIN)/dvosecfilt.$(ARCH) : $(DVOSECFILT) + +DVOSECFILT_CLIENT = \ +$(SRC)/dvosecfilt_client.$(ARCH).o \ +$(SRC)/dvosecfilt_catalogs.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/help.$(ARCH).o \ +$(SRC)/args.$(ARCH).o + +$(DVOSECFILT_CLIENT) : $(INC)/dvomerge.h +$(BIN)/dvosecfilt_client.$(ARCH) : $(DVOSECFILT_CLIENT) DVOREPAIR = \ @@ -136,5 +149,5 @@ $(BIN)/dvoverify_client.$(ARCH) : $(DVOVERIFY_CLIENT) -INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvorepair dvoverify dvoverify_client +INSTALL = dvomerge dvomerge_client dvoconvert dvosecfilt dvosecfilt_client dvorepair dvoverify dvoverify_client # dependancy rules for binary code ######################### Index: /trunk/Ohana/src/dvomerge/include/dvomerge.h =================================================================== --- /trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 35415) +++ /trunk/Ohana/src/dvomerge/include/dvomerge.h (revision 35416) @@ -36,4 +36,5 @@ int REPLACE_BY_PHOTCODE; +char *SINGLE_CPT; SkyRegion UserPatch; // used by MODE CAT @@ -103,4 +104,11 @@ void dvosecfilt_help PROTO((int argc, char **argv)); int dvosecfilt_args PROTO((int *argc, char **argv)); + +void dvosecfilt_client_usage PROTO((void)); +void dvosecfilt_client_help PROTO((int argc, char **argv)); +int dvosecfilt_client_args PROTO((int *argc, char **argv)); + +int dvosecfilt_catalogs PROTO((int Nsecfilt)); +int dvosecfilt_parallel PROTO((SkyTable *insky, int Nsecfilt)); int SkyTablePopulatedRange PROTO((off_t *ns, off_t *ne, SkyTable *sky, off_t Nstart)); Index: /trunk/Ohana/src/dvomerge/src/args.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/args.c (revision 35415) +++ /trunk/Ohana/src/dvomerge/src/args.c (revision 35416) @@ -170,8 +170,60 @@ int N; - /* extra error messages */ - VERBOSE = FALSE; - if ((N = get_argument (*argc, argv, "-v"))) { - VERBOSE = TRUE; + HOST_ID = 0; + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + SINGLE_CPT = NULL; + if ((N = get_argument (*argc, argv, "-cpt"))) { + remove_argument (N, argc, argv); + SINGLE_CPT = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (*argc, argv, "-region"))) { + remove_argument (N, argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, argc, argv); + } + + PARALLEL = FALSE; + if ((N = get_argument (*argc, argv, "-parallel"))) { + PARALLEL = TRUE; + remove_argument (N, argc, argv); + } + // this is a test mode : rather than launching the remote jobs and waiting for completion, + // relphot will simply list the remote command and wait for the user to signal completion + PARALLEL_MANUAL = FALSE; + if ((N = get_argument (*argc, argv, "-parallel-manual"))) { + PARALLEL = TRUE; // -parallel-manual implies -parallel + PARALLEL_MANUAL = TRUE; + remove_argument (N, argc, argv); + } + // this is a test mode : rather than launching the relphot_client jobs remotely, they are + // run in serial via 'system' + PARALLEL_SERIAL = FALSE; + if ((N = get_argument (*argc, argv, "-parallel-serial"))) { + if (PARALLEL_MANUAL) { + fprintf (stderr, "ERROR: cannot mix -parallel-manual and -parallel-serial\n"); + exit (1); + } + PARALLEL = TRUE; // -parallel-serial implies -parallel + PARALLEL_SERIAL = TRUE; remove_argument (N, argc, argv); } @@ -180,2 +232,63 @@ return TRUE; } + +/*** check for command line options ***/ +int dvosecfilt_client_args (int *argc, char **argv) { + + int N; + + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + SINGLE_CPT = NULL; + if ((N = get_argument (*argc, argv, "-cpt"))) { + remove_argument (N, argc, argv); + SINGLE_CPT = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (*argc, argv, "-region"))) { + remove_argument (N, argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, argc, argv); + } + + HOST_ID = 0; + if ((N = get_argument (*argc, argv, "-hostID"))) { + remove_argument (N, argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, argc, argv); + } + if (!HOST_ID) dvosecfilt_client_usage(); + + HOSTDIR = NULL; + if ((N = get_argument (*argc, argv, "-hostdir"))) { + remove_argument (N, argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + if (!HOSTDIR) dvosecfilt_client_usage(); + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + if (*argc != 3) dvosecfilt_client_usage(); + return TRUE; +} Index: /trunk/Ohana/src/dvomerge/src/dvosecfilt.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvosecfilt.c (revision 35415) +++ /trunk/Ohana/src/dvomerge/src/dvosecfilt.c (revision 35416) @@ -3,10 +3,4 @@ // modify the number of average photcodes : change only secfilt tables in place int main (int argc, char **argv) { - - char filename[256], *input; - off_t i, j, k, Nsecfilt, NsecInput, Nstart; - SkyTable *insky; - Catalog catalog; - SecFilt *insec, *outsec; SetSignals (); @@ -15,75 +9,8 @@ dvosecfilt_args (&argc, argv); - input = argv[1]; - Nsecfilt = atoi(argv[2]); + strcpy (CATDIR, argv[1]); + int Nsecfilt = atoi(argv[2]); - // load the sky table for the existing database - insky = SkyTableLoadOptimal (input, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); - SkyTableSetFilenames (insky, input, "cpt"); - - // loop over all input catalogs - for (i = 0; i < insky[0].Nregions; i++) { - if (!insky[0].regions[i].table) continue; - if (VERBOSE) fprintf (stderr, "table: %s\n", insky[0].regions[i].name); - - // set the parameters which guide catalog open/load/create - catalog.filename = insky[0].filename[i]; - catalog.Nsecfilt = 0; - - // always load all of the data (if any exists) - catalog.catflags = LOAD_SECF; - catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data - catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - - if (!dvo_catalog_open (&catalog, &insky[0].regions[i], VERBOSE, "w")) { - fprintf (stderr, "ERROR: failure to open catalog file %s\n", filename); - exit (2); - } - - if (catalog.Naves_disk == 0) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - continue; - } - - if (Nsecfilt == catalog.Nsecfilt) { - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - continue; - } - - NsecInput = catalog.Nsecfilt; - Nstart = MIN(Nsecfilt, NsecInput); - - insec = catalog.secfilt; - ALLOCATE (outsec, SecFilt, catalog.Naves_disk * Nsecfilt); - - for (k = 0; k < catalog.Naves_disk; k++) { - for (j = 0; (j < catalog.Nsecfilt) && (j < Nsecfilt); j++) { - outsec[k*Nsecfilt + j].M = insec[k*NsecInput + j].M; - outsec[k*Nsecfilt + j].dM = insec[k*NsecInput + j].dM; - outsec[k*Nsecfilt + j].Xm = insec[k*NsecInput + j].Xm; - outsec[k*Nsecfilt + j].flags = insec[k*NsecInput + j].flags; - outsec[k*Nsecfilt + j].Ncode = insec[k*NsecInput + j].Ncode; - outsec[k*Nsecfilt + j].Nused = insec[k*NsecInput + j].Nused; - outsec[k*Nsecfilt + j].M_20 = insec[k*NsecInput + j].M_20; - outsec[k*Nsecfilt + j].M_80 = insec[k*NsecInput + j].M_80; - } - for (j = Nstart; j < Nsecfilt; j++) { - outsec[k*Nsecfilt + j].M = NAN; - outsec[k*Nsecfilt + j].dM = NAN; - outsec[k*Nsecfilt + j].Xm = NAN_S_SHORT; - } - } - free (catalog.secfilt); - catalog.secfilt = outsec; - catalog.Nsecfilt = Nsecfilt; - catalog.Nsecf_mem = Nsecfilt * catalog.Naves_disk; - catalog.Nsecf_disk = Nsecfilt * catalog.Naves_disk; - - dvo_catalog_save (&catalog, VERBOSE); - dvo_catalog_unlock (&catalog); - dvo_catalog_free (&catalog); - } + dvosecfilt_catalogs (Nsecfilt); exit (0); Index: /trunk/Ohana/src/dvomerge/src/dvosecfilt_catalogs.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvosecfilt_catalogs.c (revision 35416) +++ /trunk/Ohana/src/dvomerge/src/dvosecfilt_catalogs.c (revision 35416) @@ -0,0 +1,153 @@ +# include "dvomerge.h" + +int dvosecfilt_catalogs (int Nsecfilt) { + + off_t i, j, k, NsecInput, Nstart; + Catalog catalog; + SecFilt *insec, *outsec; + char filename[DVO_MAX_PATH]; + SkyList *skylist; + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + if (PARALLEL && !HOST_ID) { + int status = dvosecfilt_parallel (sky, Nsecfilt); + return status; + } + + // determine the populated SkyRegions overlapping the requested area (default depth) + if (SINGLE_CPT) { + skylist = SkyRegionByCPT (sky, SINGLE_CPT); + } else { + skylist = SkyListByPatch (sky, -1, &UserPatch); + } + myAssert (skylist, "ooops!"); + + // loop over all input catalogs + for (i = 0; i < skylist->Nregions; i++) { + if (!skylist->regions[i]->table) continue; + if (VERBOSE) fprintf (stderr, "table: %s\n", skylist->regions[i]->name); + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist->regions[i], HOST_ID)) continue; + + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", HOSTDIR, skylist->regions[i]->name); + catalog.filename = HOST_ID ? hostfile : skylist->filename[i]; + + // always load all of the data (if any exists) + catalog.Nsecfilt = 0; + catalog.catflags = LOAD_SECF; // XXX this will fail for MEF version + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + + if (!dvo_catalog_open (&catalog, skylist->regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", filename); + exit (2); + } + + if (catalog.Naves_disk == 0) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + if (Nsecfilt == catalog.Nsecfilt) { + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + NsecInput = catalog.Nsecfilt; + Nstart = MIN(Nsecfilt, NsecInput); + + insec = catalog.secfilt; + ALLOCATE (outsec, SecFilt, catalog.Naves_disk * Nsecfilt); + + for (k = 0; k < catalog.Naves_disk; k++) { + for (j = 0; (j < catalog.Nsecfilt) && (j < Nsecfilt); j++) { + outsec[k*Nsecfilt + j] = insec[k*NsecInput + j]; + } + for (j = Nstart; j < Nsecfilt; j++) { + dvo_secfilt_init (&outsec[k*Nsecfilt + j]); + } + } + free (catalog.secfilt); + catalog.secfilt = outsec; + catalog.Nsecfilt = Nsecfilt; + catalog.Nsecf_mem = Nsecfilt * catalog.Naves_disk; + catalog.Nsecf_disk = Nsecfilt * catalog.Naves_disk; + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + + return TRUE; +} + +// launch the dvosecfilt_client jobs to the parallel hosts +int dvosecfilt_parallel (SkyTable *sky, int Nsecfilt) { + + // ensure that the paths are absolute path names + char *abscatdir = abspath (CATDIR, DVO_MAX_PATH); + + // 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; + + // options / arguments that can affect relastro_client -update-objects: + char command[DVO_MAX_PATH]; + snprintf (command, DVO_MAX_PATH, "dvosecfilt_client %s %d -hostID %d -hostdir %s -region %f %f %f %f", + abscatdir, Nsecfilt, table->hosts[i].hostID, table->hosts[i].pathname, + UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax + ); + + char tmpline[DVO_MAX_PATH]; + if (VERBOSE) { snprintf (tmpline, DVO_MAX_PATH, "%s -v", command); strcpy (command, tmpline); } + if (SINGLE_CPT) { snprintf (tmpline, DVO_MAX_PATH, "%s -cpt %s", command, SINGLE_CPT); 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 photdbc_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) { + 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 photdbc_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/src/dvosecfilt_client.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/dvosecfilt_client.c (revision 35416) +++ /trunk/Ohana/src/dvomerge/src/dvosecfilt_client.c (revision 35416) @@ -0,0 +1,17 @@ +# include "dvomerge.h" + +// modify the number of average photcodes : change only secfilt tables in place +int main (int argc, char **argv) { + + SetSignals (); + dvosecfilt_client_help (argc, argv); + ConfigInit (&argc, argv); + dvosecfilt_client_args (&argc, argv); + + strcpy (CATDIR, argv[1]); + int Nsecfilt = atoi(argv[2]); + + dvosecfilt_catalogs (Nsecfilt); + + exit (0); +} Index: /trunk/Ohana/src/dvomerge/src/help.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/help.c (revision 35415) +++ /trunk/Ohana/src/dvomerge/src/help.c (revision 35416) @@ -25,5 +25,16 @@ void dvosecfilt_usage(void) { - fprintf (stderr, "USAGE: dvosecfilt (catdir) (Nsecfilt)\n"); + fprintf (stderr, "USAGE: dvosecfilt (catdir) (Nsecfilt)\n\n"); + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -v : verbose mode\n"); + + exit (2); +} + +void dvosecfilt_client_usage(void) { + + fprintf (stderr, "USAGE: dvosecfilt_client (catdir) (Nsecfilt)\n\n"); + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -v : verbose mode\n"); exit (2); @@ -123,4 +134,26 @@ } +void dvosecfilt_client_help (int argc, char **argv) { + + /* check for help request */ + if (get_argument (argc, argv, "-help")) goto show_help; + if (get_argument (argc, argv, "-h")) goto show_help; + return; + +show_help: + + fprintf (stderr, "USAGE\n"); + fprintf (stderr, " dvosecfilt_client (catdir) (Nsecfilt)\n\n"); + + fprintf (stderr, " change number of secfilt entries in photcode table (updates secfilt tables (cps), NSECFILT in cpt files)\n"); + fprintf (stderr, " NOTE: the user must change the photcode table to reflect the change (use photcode-table -export / -import)\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} + void dvorepair_help (int argc, char **argv) { Index: /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c =================================================================== --- /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 35415) +++ /trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c (revision 35416) @@ -299,53 +299,21 @@ // XXX should we accept the input measurements for these fields? - + dvo_average_init (&output[0].average[Nave]); output[0].average[Nave].R = input[0].average[N].R; output[0].average[Nave].D = input[0].average[N].D; - output[0].average[Nave].dR = 0; - output[0].average[Nave].dD = 0; - - output[0].average[Nave].Nmeasure = 0; // this value is update as the measurements are associated with this entry below - output[0].average[Nave].Nmissing = 0; - output[0].average[Nave].Nextend = 0; - output[0].average[Nave].measureOffset = Nmeas; - output[0].average[Nave].missingOffset = -1; - output[0].average[Nave].extendOffset = -1; - - output[0].average[Nave].uR = 0; - output[0].average[Nave].uD = 0; - output[0].average[Nave].duR = 0; - output[0].average[Nave].duD = 0; - output[0].average[Nave].P = 0; - output[0].average[Nave].dP = 0; - - output[0].average[Nave].stargal = 0; - output[0].average[Nave].ChiSqAve = 0.0; - output[0].average[Nave].ChiSqPM = 0.0; - output[0].average[Nave].ChiSqPar = 0.0; - output[0].average[Nave].Tmean = 0; - output[0].average[Nave].Trange = 0; - output[0].average[Nave].Npos = 0; - output[0].average[Nave].objID = objID; // we create objID values in the context of the output db output[0].average[Nave].catID = catID; // we create catID values in the context of the output db - output[0].average[Nave].flags = 0; // XXX why reset these? + if (PSPS_ID) { output[0].average[Nave].extID = CreatePSPSObjectID(output[0].average[Nave].R, output[0].average[Nave].D); - } else { - output[0].average[Nave].extID = 0; - } + } objID ++; + // init the new secfilt entries for (j = 0; j < NsecfiltOut; j++) { int outputIndex = (Nave * NsecfiltOut) + j; - output[0].secfilt[outputIndex].M = NAN; - output[0].secfilt[outputIndex].dM = NAN; - output[0].secfilt[outputIndex].Xm = NAN_S_SHORT; - output[0].secfilt[outputIndex].M_20 = NAN_S_SHORT; - output[0].secfilt[outputIndex].M_80 = NAN_S_SHORT; - output[0].secfilt[outputIndex].Ncode = 0; - output[0].secfilt[outputIndex].Nused = 0; + dvo_secfilt_init (&output[0].secfilt[outputIndex]); } @@ -359,11 +327,5 @@ if (isfinite(input[0].secfilt[N*NsecfiltIn+j].M)) { - output[0].secfilt[outputIndex].M = input[0].secfilt[N*NsecfiltIn+j].M; - output[0].secfilt[outputIndex].dM = input[0].secfilt[N*NsecfiltIn+j].dM; - output[0].secfilt[outputIndex].Xm = input[0].secfilt[N*NsecfiltIn+j].Xm; - output[0].secfilt[outputIndex].M_20 = input[0].secfilt[N*NsecfiltIn+j].M_20; - output[0].secfilt[outputIndex].M_80 = input[0].secfilt[N*NsecfiltIn+j].M_80; - output[0].secfilt[outputIndex].Ncode = input[0].secfilt[N*NsecfiltIn+j].Ncode; - output[0].secfilt[outputIndex].Nused = input[0].secfilt[N*NsecfiltIn+j].Nused; + output[0].secfilt[outputIndex] = input[0].secfilt[N*NsecfiltIn+j]; } } Index: /trunk/Ohana/src/dvomerge/test/dvosecfilt.parallel.sh =================================================================== --- /trunk/Ohana/src/dvomerge/test/dvosecfilt.parallel.sh (revision 35416) +++ /trunk/Ohana/src/dvomerge/test/dvosecfilt.parallel.sh (revision 35416) @@ -0,0 +1,9 @@ +#!/bin/csh -f + + rm -rf catdir.merge + rm -rf catdir.merge.p? + rsync -auv catdir.2mass/ catdir.merge/ + cp HostTable.dat catdir.merge/ + dvodist -out catdir.merge + dvomerge -parallel catdir.grizy into catdir.merge + dvosecfilt -parallel catdir.merge 10 Index: /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps.c =================================================================== --- /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps.c (revision 35415) +++ /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps.c (revision 35416) @@ -57,4 +57,5 @@ } + int status = TRUE; // select measurements for each populated catalog for (i = 0; i < skylist[0].Nregions; i++) { @@ -86,7 +87,13 @@ // NOTE: this is where the real action happens if (SAVE_REMOTE && HOST_ID) { - append_detections_dvopsps_catalog (&catalog); + if (!append_detections_dvopsps_catalog (&catalog)) { + fprintf (stderr, "ERROR: failure to append detections to output catalog\n"); + status = FALSE; + } } else { - insert_detections_dvopsps_catalog (&catalog, mysqlReal); + if (!insert_detections_dvopsps_catalog (&catalog, mysqlReal)) { + fprintf (stderr, "ERROR: failure to insert detections into mysql database\n"); + status = FALSE; + } } @@ -97,8 +104,11 @@ if (SAVE_REMOTE && HOST_ID) { - save_detections_dvopsps (); - } - - return (TRUE); + if (!save_detections_dvopsps ()) { + fprintf (stderr, "ERROR: failure to save output file with detections\n"); + status = FALSE; + } + } + + return (status); } @@ -189,13 +199,17 @@ for (i = 0; i < table->Nhosts; i++) { - while ((detections = DetectionsLoad (table->hosts[i].results, &Ndetections)) == NULL) { + if ((detections = DetectionsLoad (table->hosts[i].results, &Ndetections)) == NULL) { // failed to get the data from this host. This can happen for various reasons. Give the user a chance to try again... + // fprintf (stderr, "you may run the command manually\n"); fprintf (stderr, "failed to read data from %s\n", table->hosts[i].hostname); - fprintf (stderr, "you may run the command manually\n"); + return FALSE; } free (table->hosts[i].results); table->hosts[i].results = NULL; - insert_detections_mysql_array (mysql, detections, Ndetections); + if (!insert_detections_mysql_array (mysql, detections, Ndetections)) { + fprintf (stderr, "failed to insert data for %s\n", table->hosts[i].hostname); + return FALSE; + } free (detections); @@ -218,4 +232,6 @@ insert_detections_mysql_init (&buffer); + int status = TRUE; + gettimeofday (&start, (void *) NULL); for (i = 0; i < Ndetections; i++) { @@ -224,5 +240,8 @@ insert_detections_mysql_detvalue (&buffer, &detections[i]); if (buffer.Nbuffer > MAX_BUFFER) { - insert_detections_mysql_commit (&buffer, mysql); + if (!insert_detections_mysql_commit (&buffer, mysql)) { + fprintf (stderr, "failure to insert detections in mysql (commit)\n"); + status = FALSE; + } buffer.Nbuffer = 0; bzero (buffer.buffer, buffer.Nalloc); @@ -230,9 +249,17 @@ } } - insert_detections_mysql_commit (&buffer, mysql); + if (!insert_detections_mysql_commit (&buffer, mysql)) { + fprintf (stderr, "failure to insert detections in mysql (commit 2)\n"); + status = FALSE; + } FreeIOBuffer (&buffer); + if (!status) { + MARKTIME("-- failed to insert %d rows in %f sec\n", Ndetections, dtime); + return FALSE; + } + MARKTIME("-- inserted %d rows in %f sec\n", Ndetections, dtime); - return (TRUE); + return TRUE; } Index: /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c =================================================================== --- /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c (revision 35415) +++ /trunk/Ohana/src/dvopsps/src/insert_detections_dvopsps_catalog.c (revision 35416) @@ -74,11 +74,14 @@ int save_detections_dvopsps () { + int status = TRUE; + if (!DetectionsSave (RESULT_FILE, detections, Ndetections)) { fprintf (stderr, "failed to save detection file %s\n", RESULT_FILE); + status = FALSE; } free (detections); Ndetections = 0; - return (TRUE); + return (status); } @@ -101,4 +104,6 @@ int Ninsert = 0; + int status = TRUE; + // NOTE for testing, just do a few objects if (!mysql) { @@ -115,8 +120,14 @@ // XXX check return status - insert_detections_mysql_value (&buffer, &average[i], &measure[m+j]); + if (!insert_detections_mysql_value (&buffer, &average[i], &measure[m+j])) { + fprintf (stderr, "failure to insert detections in mysql\n"); + status = FALSE; + } if (buffer.Nbuffer > MAX_BUFFER) { - insert_detections_mysql_commit (&buffer, mysql); + if (!insert_detections_mysql_commit (&buffer, mysql)) { + fprintf (stderr, "failure to insert detections in mysql (commit)\n"); + status = FALSE; + } if (DEBUG) fprintf (stderr, "inserted %d rows\n", Ninsert); Ninsert = 0; @@ -134,7 +145,16 @@ } } - insert_detections_mysql_commit (&buffer, mysql); + if (!insert_detections_mysql_commit (&buffer, mysql)) { + fprintf (stderr, "failure to insert detections in mysql (commit)\n"); + status = FALSE; + } + if (VERBOSE) fprintf (stderr, "inserted %d rows\n", Ninsert); FreeIOBuffer (&buffer); + + if (!status) { + MARKTIME("-- failed to insert "OFF_T_FMT" rows in %f sec\n", found, dtime); + return (FALSE); + } MARKTIME("-- inserted "OFF_T_FMT" rows in %f sec, skipped %d without IDs\n", found, dtime, missingID); @@ -188,4 +208,5 @@ PrintIOBuffer (buffer, "%lu, ", average->extID); // objID PrintIOBuffer (buffer, "%u, ", measure->dbFlags); // flags + float zp = code->C * 0.001 + code->K * (measure->airmass - 1) - measure->Mcal; float exptime = pow(10.0, 0.4 * measure->dt); @@ -211,4 +232,6 @@ MYSQL_RES *result; + int status = TRUE; + // check that the last two chars are ,\n and replace with ;\n if (!strcmp(&buffer->buffer[buffer->Nbuffer-2], ",\n")) { @@ -222,9 +245,10 @@ if (mysql) { if (DEBUG) fprintf (stderr, "%s\n", buffer->buffer); - int status = mysql_query(mysql, buffer->buffer); - if (status) { + int mysqlStatus = mysql_query(mysql, buffer->buffer); + if (mysqlStatus) { fprintf (stderr, "error with insert:\n"); fprintf (stderr, "%s\n", mysql_error(mysql)); fprintf (stderr, "Nbuffer: %d\n", buffer->Nbuffer); + status = FALSE; } result = mysql_store_result (mysql); @@ -234,4 +258,4 @@ } - return TRUE; -} + return status; +} Index: /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c =================================================================== --- /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 35415) +++ /trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 35416) @@ -96,22 +96,22 @@ PrintIOBuffer (&buffer, "CREATE TABLE %s_cpt (" - "RA FLOAT, " - "DEC_ FLOAT, " - "RA_ERR REAL, " - "DEC_ERR REAL, " - "U_RA REAL, " - "U_DEC REAL, " - "V_RA_ERR REAL, " - "V_DEC_ERR REAL, " - "PAR REAL, " - "PAR_ERR REAL, " - "CHISQ_POS REAL, " - "CHISQ_PM REAL, " - "CHISQ_PAP REAL, " + "RA DOUBLE, " + "DEC_ DOUBLE, " + "RA_ERR FLOAT, " + "DEC_ERR FLOAT, " + "U_RA FLOAT, " + "U_DEC FLOAT, " + "V_RA_ERR FLOAT, " + "V_DEC_ERR FLOAT, " + "PAR FLOAT, " + "PAR_ERR FLOAT, " + "CHISQ_POS FLOAT, " + "CHISQ_PM FLOAT, " + "CHISQ_PAP FLOAT, " "MEAN_EPOCH INT, " "TIME_RANGE INT, " - "PSF_QF REAL, " - "PSF_QF_PERF REAL, " - "STARGAL_SEP REAL, " + "PSF_QF FLOAT, " + "PSF_QF_PERF FLOAT, " + "STARGAL_SEP FLOAT, " "NUMBER_POS SMALLINT, " "NMEASURE SMALLINT, " @@ -152,14 +152,14 @@ PrintIOBuffer (&buffer, "CREATE TABLE %s_cps (" - "MAG REAL, " - "MAG_AP REAL, " - "MAG_KRON REAL, " - "MAG_KRON_ERR REAL, " - "MAG_ERR REAL, " - "MAG_CHI REAL, " - "FLUX_PSF REAL, " - "FLUX_PSF_ERR REAL, " - "FLUX_KRON REAL, " - "FLUX_KRON_ERR REAL, " + "MAG FLOAT, " + "MAG_AP FLOAT, " + "MAG_KRON FLOAT, " + "MAG_KRON_ERR FLOAT, " + "MAG_ERR FLOAT, " + "MAG_CHI FLOAT, " + "FLUX_PSF FLOAT, " + "FLUX_PSF_ERR FLOAT, " + "FLUX_KRON FLOAT, " + "FLUX_KRON_ERR FLOAT, " "FLAGS INT, " "NCODE SMALLINT, " Index: /trunk/Ohana/src/dvopsps/src/mysql_dvopsps.c =================================================================== --- /trunk/Ohana/src/dvopsps/src/mysql_dvopsps.c (revision 35415) +++ /trunk/Ohana/src/dvopsps/src/mysql_dvopsps.c (revision 35416) @@ -58,11 +58,11 @@ // dump_result (connection); - sprintf (query, "show variables like 'max_allowed_packet';"); - if (mysql_query (connection, query)) { - fprintf (stderr, "failed to set max_allowed_packet\n"); - fprintf (stderr, "%s\n", mysql_error (connection)); - return NULL; - } - dump_result (connection); + // sprintf (query, "show variables like 'max_allowed_packet';"); + // if (mysql_query (connection, query)) { + // fprintf (stderr, "failed to set max_allowed_packet\n"); + // fprintf (stderr, "%s\n", mysql_error (connection)); + // return NULL; + // } + // dump_result (connection); if (0) { Index: /trunk/Ohana/src/kapa2/src/DrawObjects.c =================================================================== --- /trunk/Ohana/src/kapa2/src/DrawObjects.c (revision 35415) +++ /trunk/Ohana/src/kapa2/src/DrawObjects.c (revision 35416) @@ -633,5 +633,5 @@ float *x, *y, *dxm, *dxp; double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by; - double sx0, sy0, sx1, sy1, sz, sx10, sx11; + double sx0, sy0, sx1, sy1, sz, sx10, sx11, X0, X1, Y0, Y1; mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0); @@ -651,4 +651,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[1].dfy*0.03; + X0 = graph[0].axis[0].fx; + X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + Y0 = graph[0].axis[1].fy; + Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dxp[i]))) continue; @@ -662,9 +667,9 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLine (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (sx1, sx10, sx1, sx11); + ClipLine (sx1, sx10, sx1, sx11, X0, Y0, X1, Y1); } } @@ -679,9 +684,9 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLine (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (sx1, sx10, sx1, sx11); + ClipLine (sx1, sx10, sx1, sx11, X0, Y0, X1, Y1); } } @@ -695,5 +700,5 @@ float *x, *y, *dym, *dyp; double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by; - double sx0, sy0, sx1, sy1, sz, sx10, sx11; + double sx0, sy0, sx1, sy1, sz, sx10, sx11, X0, X1, Y0, Y1; mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0); @@ -713,4 +718,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[0].dfx*0.03; + X0 = graph[0].axis[0].fx; + X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + Y0 = graph[0].axis[1].fy; + Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dyp[i]))) continue; @@ -724,9 +734,9 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLine (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (sx10, sy1, sx11, sy1); + ClipLine (sx10, sy1, sx11, sy1, X0, Y0, X1, Y1); } } @@ -741,11 +751,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLine (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (sx10, sy1, sx11, sy1); - } - } - } -} + ClipLine (sx10, sy1, sx11, sy1, X0, Y0, X1, Y1); + } + } + } +} Index: /trunk/Ohana/src/kapa2/src/PSObjects.c =================================================================== --- /trunk/Ohana/src/kapa2/src/PSObjects.c (revision 35415) +++ /trunk/Ohana/src/kapa2/src/PSObjects.c (revision 35416) @@ -578,4 +578,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[1].dfy*0.03; + double X0 = graph[0].axis[0].fx; + double X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + double Y0 = graph[0].axis[1].fy; + double Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dxp[i]))) continue; @@ -589,9 +594,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLinePS (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx0, sy0, sx1, sy1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (sx1, sx10, sx1, sx11); + ClipLinePS (sx1, sx10, sx1, sx11, X0, Y0, X1, Y1, f); + // DrawLine (sx1, sx10, sx1, sx11); } } @@ -606,9 +613,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLinePS (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx0, sy0, sx1, sy1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (sx1, sx10, sx1, sx11); + ClipLinePS (sx1, sx10, sx1, sx11, X0, Y0, X1, Y1, f); + // DrawLine (sx1, sx10, sx1, sx11); } } @@ -641,4 +650,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[0].dfx*0.03; + double X0 = graph[0].axis[0].fx; + double X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + double Y0 = graph[0].axis[1].fy; + double Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dyp[i]))) continue; @@ -652,9 +666,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLinePS (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx0, sy0, sx1, sy1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (sx10, sy1, sx11, sy1); + ClipLinePS (sx10, sy1, sx11, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx10, sy1, sx11, sy1); } } @@ -669,9 +685,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (sx0, sy0, sx1, sy1); + ClipLinePS (sx0, sy0, sx1, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx0, sy0, sx1, sy1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (sx10, sy1, sx11, sy1); + ClipLinePS (sx10, sy1, sx11, sy1, X0, Y0, X1, Y1, f); + // DrawLine (sx10, sy1, sx11, sy1); } } Index: /trunk/Ohana/src/kapa2/src/bDrawObjects.c =================================================================== --- /trunk/Ohana/src/kapa2/src/bDrawObjects.c (revision 35415) +++ /trunk/Ohana/src/kapa2/src/bDrawObjects.c (revision 35416) @@ -578,4 +578,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[1].dfy*0.03; + double X0 = graph[0].axis[0].fx; + double X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + double Y0 = graph[0].axis[1].fy; + double Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dxp[i]))) continue; @@ -589,9 +594,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (buffer, sx0, sy0, sx1, sy1); + bDrawClipLine (buffer, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx0, sy0, sx1, sy1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (buffer, sx1, sx10, sx1, sx11); + bDrawClipLine (buffer, sx1, sx10, sx1, sx11, X0, Y0, X1, Y1); + // DrawLine (buffer, sx1, sx10, sx1, sx11); } } @@ -606,9 +613,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (buffer, sx0, sy0, sx1, sy1); + bDrawClipLine (buffer, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx0, sy0, sx1, sy1); if (bar) { sx10 = sy1 - sz; sx11 = sy1 + sz; - DrawLine (buffer, sx1, sx10, sx1, sx11); + bDrawClipLine (buffer, sx1, sx10, sx1, sx11, X0, Y0, X1, Y1); + // DrawLine (buffer, sx1, sx10, sx1, sx11); } } @@ -640,4 +649,9 @@ bar = object[0].ebar; sz = object[0].size*graph[0].axis[0].dfx*0.03; + double X0 = graph[0].axis[0].fx; + double X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx; + double Y0 = graph[0].axis[1].fy; + double Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy; + for (i = 0; i < object[0].Npts; i++) { if (!(finite(x[i]) && finite(y[i]) && finite(dyp[i]))) continue; @@ -651,9 +665,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (buffer, sx0, sy0, sx1, sy1); + bDrawClipLine (buffer, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx0, sy0, sx1, sy1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (buffer, sx10, sy1, sx11, sy1); + bDrawClipLine (buffer, sx10, sy1, sx11, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx10, sy1, sx11, sy1); } } @@ -668,9 +684,11 @@ (sy1 < graph[0].axis[1].fy) && (sy1 > graph[0].axis[1].fy + graph[0].axis[1].dfy))) { - DrawLine (buffer, sx0, sy0, sx1, sy1); + bDrawClipLine (buffer, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx0, sy0, sx1, sy1); if (bar) { sx10 = sx1 - sz; sx11 = sx1 + sz; - DrawLine (buffer, sx10, sy1, sx11, sy1); + bDrawClipLine (buffer, sx10, sy1, sx11, sy1, X0, Y0, X1, Y1); + // DrawLine (buffer, sx10, sy1, sx11, sy1); } } Index: /trunk/Ohana/src/libautocode/def/cmf-ps1-dv3.d =================================================================== --- /trunk/Ohana/src/libautocode/def/cmf-ps1-dv3.d (revision 35416) +++ /trunk/Ohana/src/libautocode/def/cmf-ps1-dv3.d (revision 35416) @@ -0,0 +1,65 @@ +# name of structure type +STRUCT CMF_PS1_V3 +EXTNAME CMF_PS1_V3 +TYPE BINTABLE +SIZE 224 + +# elements of data structure / FITS table +FIELD detID, IPP_IDET, unsigned int, detection ID +FIELD X, X_PSF, float, x coord, pixels +FIELD Y, Y_PSF, float, y coord, pixels +FIELD dX, X_PSF_SIG, float, x coord error, pixels +FIELD dY, Y_PSF_SIG, float, y coord error, pixels +FIELD posangle, POSANGLE, float, Posangle at source, degrees +FIELD pltscale, PLTSCALE, float, Plate Scale at source, arcsec/pixel +FIELD M, PSF_INST_MAG, float, inst mags, mags +FIELD dM, PSF_INST_MAG_SIG, float, inst mag error, mags +FIELD Flux, PSF_INST_FLUX, float, psf flux, counts +FIELD dFlux, PSF_INST_FLUX_SIG, float, psf flux error, counts +FIELD Map, AP_MAG, float, standard aperture mag, mags +FIELD MapRaw, AP_MAG_RAW, float, raw aperture mag, mags +FIELD apRadius, AP_MAG_RADIUS, float, radius used for fit, pixels +FIELD apFlux, AP_FLUX, float, ap flux +FIELD apFluxErr, AP_FLUX_SIG, float, ap flux err +FIELD Mpeak, PEAK_FLUX_AS_MAG, float, peak flux as a mag, mags +FIELD Mcalib, CAL_PSF_MAG, float, calibrated psf mag, mags +FIELD dMcal, CAL_PSF_MAG_SIG, float, zero point scatter, mags + +# NOTE: RA & DEC (both double) need to be on an 8-byte boundary... +FIELD RA, RA_PSF, double, PSF RA coord, degrees +FIELD DEC, DEC_PSF, double, PSF DEC coord, degrees + +FIELD sky, SKY, float, sky flux, cnts/sec +FIELD dSky, SKY_SIGMA, float, sky flux error, cnts/sec +FIELD psfChisq, PSF_CHISQ, float, psf fit chisq +FIELD crNsigma, CR_NSIGMA, float, Nsigma deviations from PSF to CF +FIELD extNsigma, EXT_NSIGMA, float, Nsigma deviations from PSF to EXT +FIELD fx, PSF_MAJOR, float, psf fit major axis, pixels +FIELD fy, PSF_MINOR, float, psf fit minor axis, pixels +FIELD df, PSF_THETA, float, ellipse angle, degrees +FIELD psfQF, PSF_QF, float, quality factor +FIELD psfQFperf, PSF_QF_PERFECT, float, quality factor perfect +FIELD psfNdof, PSF_NDOF, int, psf degrees of freedom +FIELD psfNpix, PSF_NPIX, int, psf number of pixels +FIELD Mxx, MOMENTS_XX, float, second moment X, pixels^2 +FIELD Mxy, MOMENTS_XY, float, second moment Y, pixels^2 +FIELD Myy, MOMENTS_YY, float, second moment XY, pixels^2 +FIELD Mr1, MOMENTS_R1, float, first radial moment, pixels +FIELD Mrh, MOMENTS_RH, float, half radial moment, pixels^1/2 +FIELD kronFlux, KRON_FLUX, float, kron flux, counts +FIELD kronFluxErr, KRON_FLUX_ERR, float, kron flux error, counts +FIELD kronInner, KRON_FLUX_INNER, float, kron flux 1ID) { - case MEAS_MAG: /* magnitudes are already determined above */ - equiv = GetPhotcodeEquivbyCode (measure[0].photcode); - - // we return the magnitude for this measure if: - if (field->photcode->type == PHOT_MAG) goto valid_photcode; + case MEAS_MAG: { /* magnitudes are already determined above */ + PhotCode *myEquiv = GetPhotcodeEquivbyCode (measure[0].photcode); + + // if we request mag:ave, use equiv for photcode + if (field->photcode->type == PHOT_MAG) { + equiv = myEquiv; + goto valid_photcode; + } + + // if we ask for 2MASS_K, etc (REF values), return NAN unless measure->code matches if ((field->photcode->type == PHOT_REF) && (measure[0].photcode == field->photcode->code)) goto valid_photcode; + + // if we ask for GPC1.g.XY03:rel, etc (DEP values), return NAN unless measure->code matches if ((field->photcode->type == PHOT_DEP) && (measure[0].photcode == field->photcode->code)) goto valid_photcode; - if ((equiv != NULL) && (field->photcode->type == PHOT_SEC) && (equiv[0].code == field->photcode->code)) goto valid_photcode; + // if we ask for g:ave, or other SEC-level values, return the corresponding field + if (field->photcode->type == PHOT_SEC) { + switch (field->magMode) { + // measure-like : return non-NAN if measure.equiv.photcode matches field.photcode + case MAG_INST: + case MAG_CAT: + case MAG_SYS: + case MAG_REL: + case MAG_CAL: + case MAG_APER: + case MAG_APER_INST: + case MAG_KRON: + case MAG_KRON_INST: + case MAG_KRON_ERR: + case MAG_ERR: + case MAG_PHOT_FLAGS: + equiv = myEquiv; + if (equiv && (equiv->code == field->photcode->code)) goto valid_photcode; + break; + + // mean-like : return value for the given photcode + case MAG_AVE: + case MAG_REF: + case MAG_CHISQ: + case MAG_AVE_ERR: + case MAG_NCODE: + case MAG_NPHOT: + case MAG_FLUX_PSF: + case MAG_FLUX_PSF_ERR: + case MAG_FLUX_KRON: + case MAG_FLUX_KRON_ERR: + equiv = field->photcode; + goto valid_photcode; + break; + default: + fprintf (stderr, "error"); + return value; + } + } break; @@ -182,4 +241,5 @@ } break; + } case MEAS_RA: /* OK */ value.Flt = average[0].R - measure[0].dR / 3600.0; @@ -465,19 +525,69 @@ value.Flt = FromShortPixels(measure[0].dRsys); break; - case MEAS_XMOSAIC: /* OK */ - ra = average[0].R - measure[0].dR / 3600.0; - dec = average[0].D - measure[0].dD / 3600.0; - mosaic = MatchMosaicMetadata (measure[0].imageID); - if (mosaic == NULL) break; - RD_to_XY (&x, &y, ra, dec, mosaic); - value.Flt = x; + + case MEAS_XFIELD: /* offset relative to exposure center in ra,dec space */ + if (!haveFieldMeas) { + if (REMOTE_CLIENT) { + fieldc = MatchFieldMetadata (measure[0].imageID); + } else { + fprintf (stderr, "non-parallel Xmos broken\n"); + abort(); + // fieldc = MatchField (measure[0].t, measure[0].photcode); + } + if (fieldc == NULL) break; + double Rm = average[0].R - measure[0].dR / 3600.0; + double Dm = average[0].D - measure[0].dD / 3600.0; + RD_to_XY (&XFIELD_MEAS, &YFIELD_MEAS, Rm, Dm, fieldc); + } + value.Flt = XFIELD_MEAS; + break; + case MEAS_YFIELD: /* OK */ + if (!haveFieldMeas) { + if (REMOTE_CLIENT) { + fieldc = MatchFieldMetadata (measure[0].imageID); + } else { + fprintf (stderr, "non-parallel Xmos broken\n"); + abort(); + // fieldc = MatchField (measure[0].t, measure[0].photcode); + } + if (fieldc == NULL) break; + double Rm = average[0].R - measure[0].dR / 3600.0; + double Dm = average[0].D - measure[0].dD / 3600.0; + RD_to_XY (&XFIELD_MEAS, &YFIELD_MEAS, Rm, Dm, fieldc); + } + value.Flt = YFIELD_MEAS; + break; + + case MEAS_XMOSAIC: /* offset relative to exposure center in camera coords */ + if (!haveMosaicMeas) { + if (REMOTE_CLIENT) { + mosaic = MatchMosaicMetadata (measure[0].imageID); + } else { + fprintf (stderr, "non-parallel Xmos broken\n"); + abort(); + mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + } + if (mosaic == NULL) break; + double Rm = average[0].R - measure[0].dR / 3600.0; + double Dm = average[0].D - measure[0].dD / 3600.0; + RD_to_XY (&XMOS_MEAS, &YMOS_MEAS, Rm, Dm, mosaic); + } + value.Flt = XMOS_MEAS; break; case MEAS_YMOSAIC: /* OK */ - ra = average[0].R - measure[0].dR / 3600.0; - dec = average[0].D - measure[0].dD / 3600.0; - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); - if (mosaic == NULL) break; - RD_to_XY (&x, &y, ra, dec, mosaic); - value.Flt = y; + if (!haveMosaicMeas) { + if (REMOTE_CLIENT) { + mosaic = MatchMosaicMetadata (measure[0].imageID); + } else { + fprintf (stderr, "non-parallel Xmos broken\n"); + abort(); + mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + } + if (mosaic == NULL) break; + double Rm = average[0].R - measure[0].dR / 3600.0; + double Dm = average[0].D - measure[0].dD / 3600.0; + RD_to_XY (&XMOS_MEAS, &YMOS_MEAS, Rm, Dm, mosaic); + } + value.Flt = YMOS_MEAS; break; Index: /trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dbFields.c (revision 35415) +++ /trunk/Ohana/src/libdvo/src/dbFields.c (revision 35416) @@ -87,4 +87,9 @@ strcpy (code[0].name, "MAG"); code[0].type = PHOT_MAG; + // the field call 'mag' is only valid for mextract + // it should default to REL, but other types should default to AVE + if (useDefault) { + *mode = MAG_REL; + } free (tmpstring); return (code); @@ -220,4 +225,6 @@ if (!strcasecmp (fieldName, "YFIX")) ESCAPE (MEAS_YFIX, MAG_NONE, OPIHI_FLT); if (!strcasecmp (fieldName, "POS_SYS_ERR")) ESCAPE (MEAS_POS_SYS_ERR, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "XFIELD")) ESCAPE (MEAS_XFIELD, MAG_NONE, OPIHI_FLT); + if (!strcasecmp (fieldName, "YFIELD")) ESCAPE (MEAS_YFIELD, MAG_NONE, OPIHI_FLT); if (!strcasecmp (fieldName, "XMOSAIC")) ESCAPE (MEAS_XMOSAIC, MAG_NONE, OPIHI_FLT); if (!strcasecmp (fieldName, "YMOSAIC")) ESCAPE (MEAS_YMOSAIC, MAG_NONE, OPIHI_FLT); @@ -254,5 +261,5 @@ // check for code:mode in photcode name - code = ParsePhotcodeField (fieldName, &mode, MAG_REL); + code = ParsePhotcodeField (fieldName, &mode, MAG_AVE); if (code == NULL) { gprint (GP_ERR, "unknown field '%s' for measurement table in DVO database\n", fieldName); Index: /trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- /trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 35415) +++ /trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 35416) @@ -112,4 +112,14 @@ // init all data, or just catalog data +void dvo_averageT_init (AverageTiny *average) { + average->R = 0; + average->D = 0; + average->flags = 0; + average->Nmeasure = 0; + average->measureOffset = -1; + average->catID = 0; +} + +// init all data, or just catalog data void dvo_secfilt_init (SecFilt *secfilt) { secfilt->M = NAN; @@ -142,5 +152,4 @@ measure->dR = NAN; measure->dD = NAN; - measure->M = NAN; measure->Mcal = NAN; measure->Map = NAN; @@ -203,4 +212,34 @@ measure->dbFlags = 0; measure->photFlags = 0; +} + +void dvo_measureT_init (MeasureTiny *measure) { + measure->dR = NAN; + measure->dD = NAN; + measure->M = NAN; + measure->Mcal = NAN; + measure->dM = NAN; + + measure->airmass = NAN; + measure->Xccd = NAN; + measure->Yccd = NAN; + measure->Xfix = NAN; + measure->Yfix = NAN; + + measure->t = 0; + measure->dt = NAN; + measure->averef = 0; + + measure->imageID = 0; + + measure->dbFlags = 0; + measure->photFlags = 0; + measure->photcode = 0; + + measure->catID = 0; + + measure->dXccd = 0; + measure->dYccd = 0; + measure->dRsys = 0; } Index: /trunk/Ohana/src/libkapa/src/KapaOpen.c =================================================================== --- /trunk/Ohana/src/libkapa/src/KapaOpen.c (revision 35415) +++ /trunk/Ohana/src/libkapa/src/KapaOpen.c (revision 35416) @@ -120,5 +120,16 @@ return (-1); -accepted: + accepted: + { + // we need to do some minimal handshake here. I will send out + // a 4 char message : KAPA + int Nout = write (BindSocket, "KAPA", 4); + if (Nout != 4) { + if (DEBUG) fprintf (stderr, "connection failed\n"); + close (BindSocket); + return (-1); + } + } + if (DEBUG) fprintf (stderr, "connection accepted\n"); return (BindSocket); @@ -203,4 +214,19 @@ perror ("connect: "); exit (2); + } + + // apparently, I can connect on someone else's port (eg GoogleTalkPlugin) + // do a simple handshake before we set !NONBLOCK: + char line[5]; + int Nout = read (InitSocket, line, 4); + if (Nout != 4) { + if (DEBUG) fprintf (stderr, "connection failed\n"); + close (InitSocket); + return (-1); + } + if (strncmp (line, "KAPA", 4)) { + if (DEBUG) fprintf (stderr, "connection to the wrong server\n"); + close (InitSocket); + return (-1); } @@ -290,7 +316,7 @@ if (name == NULL) { - sprintf (temp, "%s -socket %s &", kapa_exec, socket_name); + sprintf (temp, "%s -socket %s", kapa_exec, socket_name); } else { - sprintf (temp, "%s -socket %s -name %s &", kapa_exec, socket_name, name); + sprintf (temp, "%s -socket %s -name %s", kapa_exec, socket_name, name); } Index: /trunk/Ohana/src/opihi/cmd.astro/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 35415) +++ /trunk/Ohana/src/opihi/cmd.astro/Makefile (revision 35416) @@ -30,4 +30,5 @@ $(SRC)/drizzle.$(ARCH).o \ $(SRC)/flux.$(ARCH).o \ +$(SRC)/fitplx.$(ARCH).o \ $(SRC)/fixwrap.$(ARCH).o \ $(SRC)/fixcols.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/cmd.astro/fitplx.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/fitplx.c (revision 35416) +++ /trunk/Ohana/src/opihi/cmd.astro/fitplx.c (revision 35416) @@ -0,0 +1,335 @@ +# include "astro.h" +# define J2000 51544.5 /* Modified Julian date at standard epoch J2000 */ + +# define ESCAPE(MSG,...) { \ + gprint (GP_ERR, MSG, __VA_ARGS__); \ + return FALSE; } + +typedef struct { + double Ro, dRo; + double Do, dDo; + + double uR, duR; + double uD, duD; + + double p, dp; + + double chisq; + int Nfit; +} PlxFit; + +int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE); +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius); +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time); + +int fitplx (int argc, char **argv) { + + int i, N; + + Vector *rvec, *dvec, *tvec, *dRvec, *dDvec; + + Vector *mvec = NULL; // mask vector + if ((N = get_argument (argc, argv, "-mask"))) { + remove_argument (N, &argc, argv); + if ((mvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + CastVector (mvec, OPIHI_INT); + } + + int VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + remove_argument (N, &argc, argv); + VERBOSE = TRUE; + } + + if (argc != 6) { + gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n"); + // what about the errors? + return (FALSE); + } + + /* select input / output buffers */ + // if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + // Nx = buf[0].header.Naxis[0]; + // Ny = buf[0].header.Naxis[1]; + + if ((rvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[1]); + if ((dvec = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[3]); + if ((tvec = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[5]); + + if ((dRvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[2]); + if ((dDvec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", argv[4]); + + double *R = rvec->elements.Flt; + double *D = dvec->elements.Flt; + double *T = tvec->elements.Flt; + + double *dR = dRvec->elements.Flt; + double *dD = dDvec->elements.Flt; + + int *mask = NULL; + if (mvec) { + mask = mvec->elements.Int; + } + + N = tvec->Nelements; // XXX check other lengths + + // find mean values to remove + double Npts = 0; + double Tmean = 0; + double Rmean = 0; + double Dmean = 0; + double Tmin = +1000000; + double Tmax = -1000000; + for (i = 0; i < N; i++) { + if (mask && !mask[i]) continue; + Rmean += R[i]; + Dmean += D[i]; + Tmean += T[i]; + Tmin = MIN(Tmin, T[i]); + Tmax = MAX(Tmax, T[i]); + Npts += 1.0; + } + Rmean /= Npts; + Dmean /= Npts; + Tmean /= Npts; + + float Trange = Tmax - Tmin; + // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax); + + /* project coordinates to a plane centered on the object with units of arcsec */ + Coords coords; + coords.crval1 = Rmean; + coords.crval2 = Dmean; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "RA---SIN"); + + double *X, *Y, *t, *pX, *pY, *dX, *dY; + ALLOCATE (X, double, N); + ALLOCATE (Y, double, N); + ALLOCATE (dX, double, N); + ALLOCATE (dY, double, N); + ALLOCATE (t, double, N); + ALLOCATE (pX, double, N); + ALLOCATE (pY, double, N); + + float pXmin = +2.0; + float pXmax = -2.0; + float pYmin = +2.0; + float pYmax = -2.0; + + int n = 0; + for (i = 0; i < N; i++) { + if (mask && !mask[i]) continue; + RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords); + dX[n] = dR[i]; + dY[n] = dD[i]; + t[n] = (T[i] - Tmean) / 365.25; + ParFactor (&pX[n], &pY[n], R[i], D[i], T[i]); + pXmin = MIN (pXmin, pX[n]); + pXmax = MAX (pXmax, pX[n]); + pYmin = MIN (pYmin, pY[n]); + pYmax = MAX (pYmax, pY[n]); + n++; + } + float dXRange = pXmax - pXmin; + float dYRange = pYmax - pYmin; + float parRange = hypot (dXRange, dYRange); + + // fprintf (stderr, "par factor range: %f\n", parRange); + + PlxFit fit; + FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE); + + // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo); + + XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords); + if (VERBOSE) { + fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo); + fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD); + fprintf (stderr, "par: %f +/- %f\n", fit.p, fit.dp); + fprintf (stderr, "chisq: %f Nfit %d\n", fit.chisq, fit.Nfit); + } + + set_variable ("RA", Rmean); + set_variable ("DEC", Dmean); + set_variable ("dR", fit.dRo); + set_variable ("dD", fit.dDo); + set_variable ("uR", fit.uR); + set_variable ("uD", fit.uD); + set_variable ("duR", fit.duR); + set_variable ("duD", fit.duD); + set_variable ("plx", fit.p); + set_variable ("dplx", fit.dp); + + set_variable ("Tmean", Tmean); + set_variable ("Trange", Trange); + set_variable ("Prange", parRange); + + set_variable ("chisq", fit.chisq); + set_variable ("Nfit", fit.Nfit); + + return (TRUE); +} + +/* do we want an init function which does the alloc and a clear function to free? */ +int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE) { + + int i; + + static double **A = NULL; + static double **B = NULL; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2; + double chisq, Xf, Yf; + + if (A == NULL) { + ALLOCATE (A, double *, 5); + ALLOCATE (B, double *, 5); + for (i = 0; i < 5; i++) { + ALLOCATE (A[i], double, 5); + ALLOCATE (B[i], double, 1); + memset (A[i], 0, 5*sizeof(double)); + memset (B[i], 0, 1*sizeof(double)); + } + } + + PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0; + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + for (i = 0; i < Npts; i++) { + + if (VERBOSE) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], dX[i], Y[i], dY[i], T[i], pR[i], pD[i]); + + /* handle case where dX or dY = 0.0 */ + wx = 1.0 / SQ(dX[i]); + wy = 1.0 / SQ(dY[i]); + + Wx += wx; + Wy += wy; + + Tx += T[i]*wx; + Ty += T[i]*wy; + + Tx2 += SQ(T[i])*wx; + Ty2 += SQ(T[i])*wy; + + PR += pR[i]*wx; + PD += pD[i]*wy; + + PRT += pR[i]*T[i]*wx; + PDT += pD[i]*T[i]*wy; + + PRX += pR[i]*X[i]*wx; + PDY += pD[i]*Y[i]*wy; + + PR2 += SQ(pR[i])*wx; + PD2 += SQ(pD[i])*wy; + + Xs += X[i]*wx; + Ys += Y[i]*wy; + + XT += X[i]*T[i]*wx; + YT += Y[i]*T[i]*wy; + } + + A[0][0] = Wx; + A[0][1] = Tx; + A[0][4] = PR; + + A[1][0] = Tx; + A[1][1] = Tx2; + A[1][4] = PRT; + + A[2][2] = Wy; + A[2][3] = Ty; + A[2][4] = PD; + + A[3][2] = Ty; + A[3][3] = Ty2; + A[3][4] = PDT; + + A[4][0] = PR; + A[4][1] = PRT; + A[4][2] = PD; + A[4][3] = PDT; + A[4][4] = PR2 + PD2; + + B[0][0] = Xs; + B[1][0] = XT; + B[2][0] = Ys; + B[3][0] = YT; + B[4][0] = PRX + PDY; + + dgaussjordan ((double **)A, (double **)B, 5, 1); + + fit[0].Ro = B[0][0]; + fit[0].uR = B[1][0]; + fit[0].Do = B[2][0]; + fit[0].uD = B[3][0]; + fit[0].p = B[4][0]; + + fit[0].dRo = sqrt(A[0][0]); + fit[0].duR = sqrt(A[1][1]); + fit[0].dDo = sqrt(A[2][2]); + fit[0].duD = sqrt(A[3][3]); + fit[0].dp = sqrt(A[4][4]); + + // add up the chi square for the fit + chisq = 0.0; + for (i = 0; i < Npts; i++) { + Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i]; + Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i]; + chisq += SQ(X[i] - Xf) / SQ(dX[i]); + chisq += SQ(Y[i] - Yf) / SQ(dY[i]); + if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq); + } + fit[0].Nfit = Npts; + + // the reduced chisq is divided by (Ndof = 2*Npts - 5) + fit[0].chisq = chisq / (2.0*Npts - 5.0); + return (TRUE); +} + +/* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */ +int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) { + + double n = mjd - J2000; // day number relative to standard epoch + double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) + double g = (357.528 + 0.9856003 * n)*RAD_DEG; // Mean anomaly + + *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees + *beta = 0.0; // approx latitude + *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees + *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU + return TRUE; +} + +/* given RA, DEC, Time, calculate the parallax factor */ +// RA,DEC are decimal degrees +// Time is MJD +int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) { + + double lambda, beta, epsilon, Radius; + + // Time must be mjd + sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius); + + double lambda_rad = lambda*RAD_DEG; + double epsilon_rad = epsilon*RAD_DEG; + double RA_rad = RA*RAD_DEG; + double DEC_rad = DEC*RAD_DEG; + + double x = Radius*cos(lambda_rad); + double y = Radius*cos(epsilon_rad)*sin(lambda_rad); + double z = Radius*sin(epsilon_rad)*sin(lambda_rad); + + *pR = +(y*cos(RA_rad) - x*sin(RA_rad)); + *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad); + + return TRUE; +} Index: /trunk/Ohana/src/opihi/cmd.astro/init.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 35415) +++ /trunk/Ohana/src/opihi/cmd.astro/init.c (revision 35416) @@ -14,4 +14,5 @@ int drizzle PROTO((int, char **)); int flux PROTO((int, char **)); +int fitplx PROTO((int, char **)); int fixwrap PROTO((int, char **)); int fiximage PROTO((int, char **)); @@ -71,4 +72,5 @@ {1, "drizzle", drizzle, "transform image to image"}, {1, "flux", flux, "flux in a convex contour"}, + {1, "fitplx", fitplx, "fit proper motion and parallax"}, {1, "fixwrap", fixwrap, "fix megacam over-wrapped pixels"}, {1, "fiximage", fiximage, "fix pixels in an image by interpolation"}, Index: /trunk/Ohana/src/opihi/cmd.data/Makefile =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/Makefile (revision 35415) +++ /trunk/Ohana/src/opihi/cmd.data/Makefile (revision 35416) @@ -150,4 +150,5 @@ $(SRC)/vsmooth.$(ARCH).o \ $(SRC)/vstats.$(ARCH).o \ +$(SRC)/xsection.$(ARCH).o \ $(SRC)/wd.$(ARCH).o \ $(SRC)/write_vectors.$(ARCH).o \ Index: /trunk/Ohana/src/opihi/cmd.data/init.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/init.c (revision 35415) +++ /trunk/Ohana/src/opihi/cmd.data/init.c (revision 35416) @@ -141,4 +141,5 @@ int wd PROTO((int, char **)); int write_vectors PROTO((int, char **)); +int xsection PROTO((int, char **)); int zap PROTO((int, char **)); int zplot PROTO((int, char **)); @@ -299,4 +300,5 @@ {1, "wd", wd, "write an image to a file"}, {1, "write", write_vectors, "write vectors to datafile"}, + {1, "xsection", xsection, "generate cross-section histogram for an image (or region)"}, {1, "zap", zap, "assign values to pixel regions"}, {1, "zplot", zplot, "plot x y with size scaled by z"}, Index: /trunk/Ohana/src/opihi/cmd.data/xsection.c =================================================================== --- /trunk/Ohana/src/opihi/cmd.data/xsection.c (revision 35416) +++ /trunk/Ohana/src/opihi/cmd.data/xsection.c (revision 35416) @@ -0,0 +1,109 @@ +# include "data.h" + +// take an image and accumulate a histogram in some direction +int xsection (int argc, char **argv) { + + int i, j, N, Nbins; + int sx, sy, nx, ny, bin; + int *vecN; + float *V, delta; + Vector *vec1, *vec2; + Buffer *buf; + + // int Quiet = FALSE; + // if ((N = get_argument (argc, argv, "-q"))) { + // Quiet = TRUE; + // remove_argument (N, &argc, argv); + // } + // if ((N = get_argument (argc, argv, "-quiet"))) { + // Quiet = TRUE; + // remove_argument (N, &argc, argv); + // } + + delta = 1.0; + if ((N = get_argument (argc, argv, "-delta"))) { + remove_argument (N, &argc, argv); + delta = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + sx = sy = 0.0; + nx = ny = 0.0; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + sx = atof (argv[N]); + remove_argument (N, &argc, argv); + sy = atof (argv[N]); + remove_argument (N, &argc, argv); + nx = atof (argv[N]); + remove_argument (N, &argc, argv); + ny = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 5) { + gprint (GP_ERR, "USAGE: xsection (angle) [-region sx sy nx ny] [-delta binsize]\n"); + return (FALSE); + } + + if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + + int Nx = buf[0].matrix.Naxis[0]; + int Ny = buf[0].matrix.Naxis[1]; + + /* if either range is set to zero, use the rest of the chip */ + if (nx == 0) nx = Nx - sx; + if (ny == 0) ny = Ny - sy; + + if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) { + gprint (GP_ERR, "region out of range\n"); + return (FALSE); + } + + + // center coords are center of the region + float Xo = sx + 0.5*nx; + float Yo = sy + 0.5*ny; + + float angle = atof(argv[4]); + float secant = MAX(2.0*fabs(cos(angle*RAD_DEG)/Nx), 2.0*fabs(sin(angle*RAD_DEG)/Ny)); + float range = 1.1 / secant; + + Nbins = 2*range / delta; + + if ((vec1 = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((vec2 = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE); + + ResetVector (vec1, OPIHI_FLT, Nbins + 1); + ResetVector (vec2, OPIHI_FLT, Nbins + 1); + bzero (vec1[0].elements.Flt, vec1[0].Nelements*sizeof(opihi_flt)); + bzero (vec2[0].elements.Flt, vec2[0].Nelements*sizeof(opihi_flt)); + + V = (float *)buf[0].matrix.buffer; + + float dx = cos(angle*RAD_DEG); + float dy = sin(angle*RAD_DEG); + + ALLOCATE (vecN, int, Nbins + 1); + memset (vecN, 0, (Nbins + 1)*sizeof(int)); + + for (j = sy; j < sy + ny; j++) { + float dY = (j - Yo); + for (i = sx; i < sx + nx; i++) { + float value = V[i*Nx + j]; + if (!isfinite(value)) continue; + float dX = (i - Xo); + float L = (dX * dx) + (dY * dy); + bin = MAX (MIN (Nbins, (L / delta) + 0.5*Nbins), 0); + vec2[0].elements.Flt[bin] += value; + vecN[bin] ++; + } + } + for (i = 0; i < Nbins + 1; i++) { + vec1[0].elements.Flt[i] = (i - 0.5*Nbins)*delta; + vec2[0].elements.Flt[i] /= (float) vecN[i]; + } + + return (TRUE); +} + Index: /trunk/Ohana/src/opihi/dvo/avextract.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/avextract.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/avextract.c (revision 35416) @@ -137,5 +137,5 @@ // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results if (PARALLEL && !HOST_ID) { - int status = HostTableParallelOps (argc, argv, RESULT_FILE, TRUE, 0, VERBOSE); + int status = HostTableParallelOps (skylist, argc, argv, RESULT_FILE, TRUE, 0, VERBOSE); dbFreeFields (fields, Nfields); Index: /trunk/Ohana/src/opihi/dvo/avmatch.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/avmatch.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/avmatch.c (revision 35416) @@ -95,5 +95,5 @@ // I need to pass the RA & DEC vectors to the remote clients... - int status = HostTableParallelOps (argc, argv, RESULT_FILE, TRUE, RAvec->Nelements, VERBOSE); + int status = HostTableParallelOps (skylist, argc, argv, RESULT_FILE, TRUE, RAvec->Nelements, VERBOSE); if (vec) free (vec); Index: /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/dvo_host_utils.c (revision 35416) @@ -8,5 +8,5 @@ # define DIE(WHO,MSG) { perror(WHO); myAbort(MSG); } -int HostTableLaunchJobs (HostTable *table, char *basecmd, char *options, int VERBOSE) { +int HostTableLaunchJobs (SkyList *sky, HostTable *table, char *basecmd, char *options, int VERBOSE) { char uniquer[12]; @@ -35,6 +35,22 @@ int top_status = TRUE; - int i; + int i, j; for (i = 0; i < table->Nhosts; i++) { + + if (sky && (sky->Nregions < table->Nhosts)) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], table->hosts[i].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", table->hosts[i].hostname); + continue; + } + // fprintf (stderr, "not skip host %s\n", table->hosts[i].hostname); + } // ensure that the paths are absolute path names @@ -88,14 +104,14 @@ // an alternative ending step ignores the result files and instead saves the names into // the list 'result:n' for the user to access as desired -int HostTableParallelOps (int argc, char **argv, char *ResultFile, int ReadVectors, int Nelements, int VERBOSE) { +int HostTableParallelOps (SkyList *sky, int argc, char **argv, char *ResultFile, int ReadVectors, int Nelements, int VERBOSE) { int i; - // load the list of hosts - SkyTable *sky = GetSkyTable(); - if (!sky) { - gprint (GP_ERR, "failed to load sky table for database\n"); - return FALSE; - } + // XX // load the list of hosts + // XX SkyTable *sky = GetSkyTable(); + // XX if (!sky) { + // XX gprint (GP_ERR, "failed to load sky table for database\n"); + // XX return FALSE; + // XX } char *tmppath = GetCATDIR (); @@ -145,5 +161,5 @@ // launch this command remotely - HostTableLaunchJobs (table, basecmd, options, VERBOSE); + HostTableLaunchJobs (sky, table, basecmd, options, VERBOSE); free (options); free (basecmd); @@ -183,4 +199,5 @@ if (table->hosts[i].status) continue; + if (!table->hosts[i].pid) continue; if (ReadVectors) { Index: /trunk/Ohana/src/opihi/dvo/mextract.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/mextract.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/mextract.c (revision 35416) @@ -16,5 +16,5 @@ if (field->ID == MEAS_CENTER_OFFSET) return TRUE; // 0.5*NX, 0.5*NY if (field->ID == MEAS_EXPNAME_AS_INT) return TRUE; // expname (or as int) - if (field->ID == MEAS_AIRMASS) return TRUE; // airmass + if (field->ID == MEAS_MEAN_AIRMASS) return TRUE; // airmass return FALSE; } @@ -171,5 +171,5 @@ // call the remote client - int status = HostTableParallelOps (targc, targv, RESULT_FILE, TRUE, 0, VERBOSE); + int status = HostTableParallelOps (skylist, targc, targv, RESULT_FILE, TRUE, 0, VERBOSE); dbFreeFields (fields, Nfields); Index: /trunk/Ohana/src/opihi/dvo/mmatch.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/mmatch.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/mmatch.c (revision 35416) @@ -106,76 +106,4 @@ RAvec = NULL; DECvec = NULL; - if (PARALLEL && !HOST_ID) { - - // We need to copy the args to a temp array and modify them so that we send the - // correct set to the remote client. The args list looks like this: - // if (!CoordsFile) : mmatch (RA) (DEC) (RADIUS) field, ... - // if ( CoordsFile) : mmatch (RADIUS) field, ... [because we stripped off the -coords filename elements] - - // allocate the temp array and copy all but (RA) (DEC) - int targc = 0; - char **targv = NULL; - ALLOCATE (targv, char *, argc + 2); - for (i = 0; i < argc; i++) { - if (!CoordsFile && (i == 1)) continue; - if (!CoordsFile && (i == 2)) continue; - targv[targc] = strcreate (argv[i]); - targc ++; - } - - // if not specified, create the coords.fits input file - if (!CoordsFile) { - // get vectors corresponding to coordinates of interest - if ((RAvec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto help; - if ((DECvec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto help; - - ALLOCATE (vec, Vector *, 2); - vec[0] = RAvec; - vec[1] = DECvec; - - CoordsFile = abspath("coords.fits", 1024); - int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, NULL); - if (!status) goto escape; - } - - // add the coords file to the args list - targv[targc+0] = strcreate ("-coords"); - targv[targc+1] = CoordsFile; // this gets freed with targv - targc += 2; - - // if needed, add the index vector to the args list - if (IDXvec) { - REALLOCATE (targv, char *, targc + 2); - targv[targc+0] = strcreate ("-index"); - targv[targc+1] = strcreate (IDXvec[0].name); - targc += 2; - } - - if (loadImages) { - Image *image; - off_t Nimage; - if ((image = LoadImagesDVO (&Nimage)) == NULL) goto escape; - - char *filename = abspath("image.metadata.fits", DVO_MAX_PATH); - ImageMetadataSave (filename, image, Nimage); - - REALLOCATE (targv, char *, targc + 2); - targv[targc+0] = strcreate ("-image-metadata"); - targv[targc+1] = strcreate (filename); - targc += 2; - } - - // call the remote client - int status = HostTableParallelOps (targc, targv, RESULT_FILE, TRUE, 0, VERBOSE); - if (vec) free (vec); - - // free up targv - for (i = 0; i < targc; i++) { - free (targv[i]); - } - free (targv); - - return status; - } // get vectors corresponding to coordinates of interest @@ -193,4 +121,77 @@ remove_argument (1, &argc, argv); } + + /* load regions which contain all supplied RA,DEC coordinates */ + if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape; + + if (PARALLEL && !HOST_ID) { + + // We need to copy the args to a temp array and modify them so that we send the + // correct set to the remote client. The args list looks like this: + // if (!CoordsFile) : mmatch (RADIUS) field, ... [we removed RA & DEC above] + // if ( CoordsFile) : mmatch (RADIUS) field, ... [because we stripped off the -coords filename elements] + + // allocate the temp array and copy all but (RA) (DEC) + int targc = 0; + char **targv = NULL; + ALLOCATE (targv, char *, argc + 2); + for (i = 0; i < argc; i++) { + targv[targc] = strcreate (argv[i]); + targc ++; + } + + // if not specified, create the coords.fits input file + // NOTE: RAvec, DECvec were set above + if (!CoordsFile) { + ALLOCATE (vec, Vector *, 2); + vec[0] = RAvec; + vec[1] = DECvec; + + // XXX this is now set for both cases... + CoordsFile = abspath("coords.fits", 1024); + int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, NULL); + if (!status) goto escape; + } + + // add the coords file to the args list + targv[targc+0] = strcreate ("-coords"); + targv[targc+1] = CoordsFile; // this gets freed with targv + targc += 2; + + // if needed, add the index vector to the args list + if (IDXvec) { + REALLOCATE (targv, char *, targc + 2); + targv[targc+0] = strcreate ("-index"); + targv[targc+1] = strcreate (IDXvec[0].name); + targc += 2; + } + + if (loadImages) { + Image *image; + off_t Nimage; + if ((image = LoadImagesDVO (&Nimage)) == NULL) goto escape; + + char *filename = abspath("image.metadata.fits", DVO_MAX_PATH); + ImageMetadataSave (filename, image, Nimage); + + REALLOCATE (targv, char *, targc + 2); + targv[targc+0] = strcreate ("-image-metadata"); + targv[targc+1] = strcreate (filename); + targc += 2; + } + + // call the remote client + int status = HostTableParallelOps (skylist, targc, targv, RESULT_FILE, TRUE, 0, VERBOSE); + if (vec) free (vec); + + // free up targv + for (i = 0; i < targc; i++) { + free (targv[i]); + } + free (targv); + + return status; + } // END of remote call section + RADIUS = atof (argv[1]); remove_argument (1, &argc, argv); @@ -209,6 +210,4 @@ } - /* load regions which contain all supplied RA,DEC coordinates */ - if ((skylist = SelectRegionsByCoordVectors (RAvec, DECvec)) == NULL) goto escape; /* create output storage vectors */ Index: /trunk/Ohana/src/opihi/dvo/photometry.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/photometry.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/photometry.c (revision 35416) @@ -1127,5 +1127,5 @@ double value; Image *image; - Coords *mosaic; + // Coords *mosaic; value = 0; @@ -1193,8 +1193,9 @@ # endif break; +# if 0 case MEAS_XMOSAIC: /* OK */ ra = average[0].R - measure[0].dR / 3600.0; dec = average[0].D - measure[0].dD / 3600.0; - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaic (measure[0].t, measure[0].photcode); // XXX not used anymore if (mosaic == NULL) break; RD_to_XY (&x, &y, ra, dec, mosaic); @@ -1204,9 +1205,10 @@ ra = average[0].R - measure[0].dR / 3600.0; dec = average[0].D - measure[0].dD / 3600.0; - mosaic = MatchMosaic (measure[0].t, measure[0].photcode); + mosaic = MatchMosaic (measure[0].t, measure[0].photcode); // XXX not used anymore if (mosaic == NULL) break; RD_to_XY (&x, &y, ra, dec, mosaic); value = y; break; +# endif } return (value); Index: /trunk/Ohana/src/opihi/dvo/remote.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/remote.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/remote.c (revision 35416) @@ -61,6 +61,18 @@ } + // load the list of hosts + SkyTable *sky = GetSkyTable(); + if (!sky) { + gprint (GP_ERR, "failed to load sky table for database\n"); + return FALSE; + } + SkyList *skylist = NULL; + ALLOCATE (skylist, SkyList, 1); + skylist[0].Nregions = sky[0].Nregions; + strcpy (skylist[0].hosts, sky[0].hosts); + // strip of the 'remote' and send the remaining arguments to the remote machine - int status = HostTableParallelOps (argc - 1, &argv[1], NULL, ReadVectors, 0, VERBOSE); + int status = HostTableParallelOps (skylist, argc - 1, &argv[1], NULL, ReadVectors, 0, VERBOSE); + free (skylist); return status; } Index: /trunk/Ohana/src/opihi/dvo/skyregion.c =================================================================== --- /trunk/Ohana/src/opihi/dvo/skyregion.c (revision 35415) +++ /trunk/Ohana/src/opihi/dvo/skyregion.c (revision 35416) @@ -98,4 +98,5 @@ new[0].Nregions = 0; new[0].ownElements = FALSE; // this list is only holding a view to the elements + strcpy (new[0].hosts, sky[0].hosts); // output list @@ -107,4 +108,5 @@ list[0].Nregions = 0; list[0].ownElements = FALSE; // this list is only holding a view to the elements + strcpy (list[0].hosts, sky[0].hosts); for (i = 0; i < Npts; i++) { Index: /trunk/Ohana/src/opihi/include/dvoshell.h =================================================================== --- /trunk/Ohana/src/opihi/include/dvoshell.h (revision 35415) +++ /trunk/Ohana/src/opihi/include/dvoshell.h (revision 35416) @@ -99,6 +99,6 @@ dbValue dbExtractImages PROTO((Image *image, off_t Nimage, off_t N, dbField *field)); -int HostTableLaunchJobs PROTO((HostTable *table, char *basecmd, char *options, int VERBOSE)); -int HostTableParallelOps PROTO((int argc, char **argv, char *ResultFile, int ReadVectors, int Nelements, int VERBOSE)); +int HostTableLaunchJobs PROTO((SkyList *sky, HostTable *table, char *basecmd, char *options, int VERBOSE)); +int HostTableParallelOps PROTO((SkyList *sky, int argc, char **argv, char *ResultFile, int ReadVectors, int Nelements, int VERBOSE)); int HostTableReloadResults PROTO((char *uniquer, int VERBOSE)); int HostTableGetResults PROTO((char *uniquer, int VERBOSE)); Index: /trunk/Ohana/src/relastro/Makefile =================================================================== --- /trunk/Ohana/src/relastro/Makefile (revision 35415) +++ /trunk/Ohana/src/relastro/Makefile (revision 35416) @@ -64,4 +64,6 @@ $(SRC)/high_speed_objects.$(ARCH).o \ $(SRC)/high_speed_utils.$(ARCH).o \ +$(SRC)/hpm_catalogs.$(ARCH).o \ +$(SRC)/hpm_objects.$(ARCH).o \ $(SRC)/relastro_merge_source.$(ARCH).o \ $(SRC)/resort_catalog.$(ARCH).o \ @@ -106,4 +108,6 @@ $(SRC)/high_speed_objects.$(ARCH).o \ $(SRC)/high_speed_utils.$(ARCH).o \ +$(SRC)/hpm_catalogs.$(ARCH).o \ +$(SRC)/hpm_objects.$(ARCH).o \ $(SRC)/relastro_merge_source.$(ARCH).o \ $(SRC)/resort_catalog.$(ARCH).o \ Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 35415) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 35416) @@ -20,5 +20,5 @@ typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode; -typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS} FitTarget; +typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS, TARGET_HPM} FitTarget; typedef enum { @@ -442,2 +442,6 @@ PhotCode **ParsePhotcodeList (char *rawlist, int *nphotcodes, int needAve); + +int hpm_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath); +int hpm_catalogs_parallel (SkyList *skylist); +int hpm_objects (SkyRegion *region, Catalog *catalog); Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 35415) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 35416) @@ -224,5 +224,5 @@ if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { - // project all of the R,D coordinates to a plane centered on this coordinate set + // project all of the R,D coordinates to a plane centered on this coordinate. set // the times to be relative to Tmean (this is required for parallax as well) for (k = 0; k < N; k++) { @@ -270,4 +270,7 @@ catalog[i].average[j].flags |= ID_STAR_FIT_PAR; Npar ++; + } else { + // need to set mode = FIT_PM_ONLY if we do not fit for parallax + mode = FIT_PM_ONLY; } } @@ -321,5 +324,5 @@ fit.uD); - //make sure that the fit succeeded + // make sure that the fit succeeded status = TRUE; status &= finite(fit.Ro); @@ -364,6 +367,6 @@ catalog[i].average[j].ChiSqPM = fitPM.chisq; catalog[i].average[j].ChiSqPar = fitPAR.chisq; - catalog[i].average[j].Tmean = (Tmean * 86400 * 365.26) + T2000; - catalog[i].average[j].Trange = (Trange * 86400 * 365.26); + catalog[i].average[j].Tmean = (Tmean * 86400 * 365.25) + T2000; + catalog[i].average[j].Trange = (Trange * 86400 * 365.25); catalog[i].average[j].Npos = fit.Nfit; if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n", @@ -385,5 +388,5 @@ } - fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT"\n", (NaveSum + NpmSum + NparSum), NaveSum, NpmSum, NparSum, NskipSum); + if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT"\n", (NaveSum + NpmSum + NparSum), NaveSum, NpmSum, NparSum, NskipSum); return (TRUE); } Index: /trunk/Ohana/src/relastro/src/args.c =================================================================== --- /trunk/Ohana/src/relastro/src/args.c (revision 35415) +++ /trunk/Ohana/src/relastro/src/args.c (revision 35416) @@ -89,5 +89,5 @@ if ((N = get_argument (argc, argv, "-high-speed"))) { // XXX include a parallax / no-parallax option - if (N >= argc - 3) usage(); + if (N >= argc - 4) usage(); FIT_TARGET = TARGET_HIGH_SPEED; remove_argument (N, &argc, argv); @@ -95,4 +95,14 @@ remove_argument (N, &argc, argv); PHOTCODE_B_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + RADIUS = atof(argv[N]); + remove_argument (N, &argc, argv); + HIGH_SPEED_DIR = abspath(argv[N], DVO_MAX_PATH); + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-hpm"))) { + if (N >= argc - 2) usage(); + FIT_TARGET = TARGET_HPM; remove_argument (N, &argc, argv); RADIUS = atof(argv[N]); @@ -106,6 +116,6 @@ remove_argument (N, &argc, argv); PARALLEL_OUTPUT = TRUE; - if (FIT_TARGET != TARGET_HIGH_SPEED) { - fprintf (stderr, "-parallel-output only valid for -high-speed mode\n"); + if ((FIT_TARGET != TARGET_HIGH_SPEED) && (FIT_TARGET != TARGET_HPM)) { + fprintf (stderr, "-parallel-output only valid for -high-speed or -hpm modes\n"); exit (1); } @@ -456,5 +466,5 @@ if ((N = get_argument (argc, argv, "-high-speed"))) { // XXX include a parallax / no-parallax option - if (N >= argc - 3) usage_client(); + if (N >= argc - 5) usage_client(); FIT_TARGET = TARGET_HIGH_SPEED; remove_argument (N, &argc, argv); @@ -466,4 +476,14 @@ remove_argument (N, &argc, argv); HIGH_SPEED_DIR = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-hpm"))) { + if (N >= argc - 3) usage(); + FIT_TARGET = TARGET_HPM; + remove_argument (N, &argc, argv); + RADIUS = atof(argv[N]); + remove_argument (N, &argc, argv); + HIGH_SPEED_DIR = abspath(argv[N], DVO_MAX_PATH); remove_argument (N, &argc, argv); } @@ -668,4 +688,5 @@ fprintf (stderr, " OR: relastro -update-objects [options]\n"); fprintf (stderr, " OR: relastro -high-speed [options]\n"); + fprintf (stderr, " OR: relastro -hpm [options]\n"); fprintf (stderr, " OR: relastro -merge-source [options]\n\n"); @@ -679,4 +700,5 @@ fprintf (stderr, " -update-mosaics\n"); fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n"); + fprintf (stderr, " -hpm (radius) (output catdir)\n"); fprintf (stderr, " -merge-source (objID) (catID) into (objID) (catID)\n\n"); Index: /trunk/Ohana/src/relastro/src/hpm_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/hpm_catalogs.c (revision 35416) +++ /trunk/Ohana/src/relastro/src/hpm_catalogs.c (revision 35416) @@ -0,0 +1,203 @@ +# include "relastro.h" + +int hpm_catalogs_parallel (SkyList *sky); + +int hpm_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath) { + + int i; + Catalog catalog; + + // tell libdvo the CATDIR + dvo_set_catdir(CATDIR); + + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + hpm_catalogs_parallel (skylist); + goto finish; + } + + // load data from each region file, only use bright stars + 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; + + // set up the basic catalog info + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : skylist[0].filename[i]; + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "r")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (!catalog.Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + hpm_objects (skylist[0].regions[i], &catalog); + + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + } + +finish: + + // in relastro (NOT relastr_client), save the Images.dat, Photcodes.dat, SkyTable.fits files + if (!hostID) { + char photcodeFile[1024]; + snprintf (photcodeFile, 1024, "%s/Photcodes.dat", HIGH_SPEED_DIR); + if (!SavePhotcodesFITS (photcodeFile)) { + fprintf (stderr, "error saving photcode table %s\n", photcodeFile); + exit (1); + } + + // need to copy across the Images, SkyTable, and Photcode tables: + char *skyfile = SkyTableFilename (HIGH_SPEED_DIR); + SkyTableSave (sky, skyfile); + + char line[2048]; + snprintf (line, 2048, "cp %s/Images.dat %s/Images.dat", CATDIR, HIGH_SPEED_DIR); + int status = system (line); + if (status) { + fprintf (stderr, "copy of Images.dat failed\n"); + exit (3); + } + } + + return (TRUE); +} + +// CATDIR is supplied globally +# define DEBUG 1 +int hpm_catalogs_parallel (SkyList *skylist) { + + // launch the setphot_client jobs to the parallel hosts + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, skylist->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", skylist->hosts, CATDIR); + exit (1); + } + + // If we want to parallelize the output, the remote client needs to load the HostTable + // for the output catdir. Require the output HostTable to match the input (ie, same + // hostIDs must exist). Preferred, but not required, that they match in the host + // location. + HostTable *tableOut = NULL; + if (PARALLEL_OUTPUT) { + // load the list of hosts + tableOut = HostTableLoad (HIGH_SPEED_DIR, skylist->hosts); + if (!tableOut) { + fprintf (stderr, "ERROR: failure reading Host Table %s for output database %s\n", skylist->hosts, HIGH_SPEED_DIR); + exit (1); + } + if (table->Nhosts != tableOut->Nhosts) { + fprintf (stderr, "ERROR: output HostTable must have matching hosts (%d in, %d out)\n", table->Nhosts, tableOut->Nhosts); + 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; + + // use this directory to save the output catalogs (distributed if parallel) + char *outputDir = strcreate(HIGH_SPEED_DIR); + if (PARALLEL_OUTPUT) { + // find the matching host entry in the output HostTable (match by ID, not index) + int hostID = table->hosts[i].hostID; + int index = table->index[hostID]; + if (index == -1) { + fprintf (stderr, "ERROR: output HostTable must have matching hosts (host ID %d not found)\n", hostID); + exit (1); + } + + free (outputDir); + outputDir = abspath (tableOut->hosts[index].pathname, DVO_MAX_PATH); + } + + char command[1024]; + snprintf (command, 1024, "relastro_client -hpm %f %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", + RADIUS, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + free (outputDir); + + // options / arguments that can affect relastro_client -high-speed + char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command); strcpy (command, tmpline); } + + if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); strcpy (command, tmpline); } + + if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command); strcpy (command, tmpline); } + + if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } + if (PhotFlagSelect) { snprintf (tmpline, 1024, "%s +photflags", command); strcpy (command, tmpline); } + if (PhotFlagBad) { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad); strcpy (command, tmpline); } + if (PhotFlagPoor) { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor); strcpy (command, tmpline); } + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + free (tstart); + free (tstop); + 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 relastro_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 relastro_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/relastro/src/hpm_objects.c =================================================================== --- /trunk/Ohana/src/relastro/src/hpm_objects.c (revision 35416) +++ /trunk/Ohana/src/relastro/src/hpm_objects.c (revision 35416) @@ -0,0 +1,298 @@ +# include "relastro.h" +# define MAX_TRANGE 180.0 +# define MIN_PS1_DET 3 +# define MAX_CHISQ_PM 1000000.0 +# define MIN_NPOS_FRAC 1.0 + +# define NEXT_I { if (Ngroup < 2) slowMoving[ni] = TRUE; newI = TRUE; i++; continue; } +# define NEXT_J { j++; continue; } + +int hpm_objects (SkyRegion *region, Catalog *catalog) { + + off_t i, j, J, ni, nj, *N1, k, m; + int *slowMoving; + double *X1, *Y1; + double dX, dY, dR, RADIUS2; + Coords tcoords; + Catalog catalogOut, testcat; + + int Nsecfilt; + char filename[1024]; + + // we need at least 2 objects if we are going to match anything... + if (catalog[0].Naverage < 2) return (TRUE); + if (VERBOSE) fprintf (stderr, "checking "OFF_T_FMT" objects\n", catalog[0].Naverage); + + // we save the best objects in an hpm dvodb + snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name); + + dvo_catalog_init (&catalogOut, TRUE); /* init new catalog */ + catalogOut.filename = strcreate(filename); + + Nsecfilt = GetPhotcodeNsecfilt(); + catalogOut.filename = filename; + catalogOut.Nsecfilt = Nsecfilt; + catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data + + catalogOut.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalogOut.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + + if (!dvo_catalog_open (&catalogOut, region, VERBOSE2,"w")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", + filename); + exit (2); + } + + off_t NAVERAGE = 1000; + off_t NMEASURE = 10000; + REALLOCATE (catalogOut.average, Average, NAVERAGE); + REALLOCATE (catalogOut.measure, Measure, NMEASURE); + REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt); + catalogOut.Naverage = 0; + catalogOut.Nmeasure = 0; + + // testcat is used to determine the fit for a single object group + // objects which do not have a high-quality testcat fit are not kept + dvo_catalog_init (&testcat, TRUE); + testcat.Naverage = 1; // this is fixed -- only one obj in testcat + testcat.Nsecfilt = Nsecfilt; + off_t NMEASURE_TESTCAT = 1000; + ALLOCATE (testcat.average, Average, 1); + ALLOCATE (testcat.measure, Measure, NMEASURE_TESTCAT); + ALLOCATE (testcat.secfilt, SecFilt, Nsecfilt); + + // mask with which to mark objects to be ignored + ALLOCATE (slowMoving, int, catalog[0].Naverage); + memset (slowMoving, 0, catalog[0].Naverage*sizeof(int)); + off_t Nslow = 0; + + time_t T2000 = ohana_date_to_sec ("2000/01/01"); + + // mark (exclude) objects which do not meet the selection criterion + for (i = 0; i < catalog[0].Naverage; i++) { + + int XVERB = FALSE; + XVERB |= (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object\n"); + } + + // count the PS1 detections via explicit photcode ranges? + // XXX this is a total hard-wired hack... + int Nps1 = 0; + double Tmin = +100.0*365.0; // +/- 100 years + double Tmax = -100.0*365.0; + m = catalog[0].average[i].measureOffset; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + if (catalog[0].measure[m+j].photcode < 10000) continue; + if (catalog[0].measure[m+j].photcode > 10600) continue; + double To = (catalog[0].measure[m+j].t - T2000) / 86400.0 ; // time in days relative to J2000 in years + Tmin = MIN(Tmin, To); + Tmax = MAX(Tmax, To); + Nps1 ++; + } + if (Nps1 < MIN_PS1_DET) { + slowMoving[i] = TRUE; + Nslow ++; + continue; + } + // selection criteria : (Nps1 > 2) && (Trange < 180) + double Trange = Tmax - Tmin; + if (Trange > MAX_TRANGE) { + slowMoving[i] = TRUE; + Nslow ++; + continue; + } + } + + fprintf (stderr, OFF_T_FMT" slow, "OFF_T_FMT" total objects; "OFF_T_FMT" possible fast\n", Nslow, catalog[0].Naverage, catalog[0].Naverage - Nslow); + if (catalog[0].Naverage == Nslow) { + fprintf (stderr, "no possible fast objects, skipping this catalog\n"); + dvo_catalog_save (&catalogOut, VERBOSE2); + dvo_catalog_unlock (&catalogOut); + dvo_catalog_free (&catalogOut); + free (slowMoving); + return TRUE; + } + + // double loop over unmarked objects (sorted in RA / X) + // record all pairs within the desired match distance + + // create tmp local positional index + ALLOCATE (X1, double, catalog[0].Naverage); + ALLOCATE (Y1, double, catalog[0].Naverage); + ALLOCATE (N1, off_t, catalog[0].Naverage); + + // define a local projection + tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); + if (region[0].Dmax < 90) { + tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax); + } else { + tcoords.crval2 = 90.0; + } + tcoords.crpix1 = 0; + tcoords.crpix2 = 0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 1; + strcpy (tcoords.ctype, "RA---ARC"); + + /* build spatial index (RA sort) referencing input array sequence */ + for (i = 0; i < catalog[0].Naverage; i++) { + int status = RD_to_XY (&X1[i], &Y1[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N1[i] = i; + assert (status); + } + sort_coords_index (X1, Y1, N1, catalog[0].Naverage); + + RADIUS2 = SQ(RADIUS); + + // group is a list of objects that are within a clump. this set will be tested + // via a clipped fit to the measurements after the group is identified + int Ngroup = 0; + int NGROUP = 100; + int *group = NULL; + ALLOCATE (group, int, NGROUP); + + // in the loop below, we need to do a bunch of things when we go to the next main object + int newI = TRUE; + + // mark (skip) objects with both sets of target photcodes + // the loop below is attempting to find associations of multiple objects which have + // passed the cuts above. the index i is following the primary object of interest + // the index j is used to explore possible near neighbors. + // When we go to the next object 'i', Nmatch is reset + for (i = j = 0; (i < catalog[0].Naverage) && (j < catalog[0].Naverage);) { + + ni = N1[i]; + nj = N1[j]; + + if (newI) { + Ngroup = 1; + group[0] = ni; + newI = FALSE; + } + + int XVERB = FALSE; + XVERB |= (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object %d or %d\n", (int) ni, (int) nj); + } + + if (slowMoving[ni]) NEXT_I; + if (slowMoving[nj]) NEXT_J; + + if (!finite(X1[i]) || !finite(Y1[i])) NEXT_I; + if (!finite(X1[j]) || !finite(Y1[j])) NEXT_J; + + // look for pairs that are within the maximum separation + dX = X1[i] - X1[j]; + if (XVERB) { + fprintf (stderr, "%d %d : %f\n", (int) ni, (int) nj, dX); + } + + if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small + if (dX >= +1.02*RADIUS) NEXT_J; // positive dX: j is too small + + // within match range; look for valid matches & accumulate the group + for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) { + if (J == i) continue; // avoid auto-matches + nj = N1[J]; + + dX = X1[i] - X1[J]; + + if (slowMoving[nj]) continue; + + XVERB = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object pt2 %d or %d : %f\n", (int) ni, (int) nj, dX); + fprintf (stderr, ".\n"); + } + + dY = Y1[i] - Y1[J]; + dR = dX*dX + dY*dY; + if (dR > RADIUS2) continue; + + /*** a match is found ***/ + group[Ngroup] = nj; + Ngroup ++; + CHECK_REALLOCATE (group, int, NGROUP, Ngroup, 100); + } + + if (Ngroup < 2) NEXT_I; + + // we now have spatially associated group of objects. now we need to see if the set of + // measurements can be fitted reasonably with a proper motion (& parallax?) + + // the mean object will start with info from the primary object + // remember: testcat.Naverage = 1 -- does not change + ni = group[0]; + int Nmatchmeas = 0; + testcat.average[0] = catalog[0].average[ni]; + testcat.average[0].measureOffset = 0; + for (J = 0; J < Ngroup; J++) { + nj = group[J]; + m = catalog[0].average[nj].measureOffset; + for (k = 0; k < catalog[0].average[nj].Nmeasure; k++) { + testcat.measure[Nmatchmeas] = catalog[0].measure[m+k]; + /* Set offset RA and Dec wrt correct average value*/ + testcat.measure[Nmatchmeas].dR = catalog[0].measure[m+k].dR + 3600.0*(catalog[0].average[ni].R - catalog[0].average[nj].R); + testcat.measure[Nmatchmeas].dD = catalog[0].measure[m+k].dD + 3600.0*(catalog[0].average[ni].D - catalog[0].average[nj].D); + testcat.measure[Nmatchmeas].averef = 0; + Nmatchmeas++; + CHECK_REALLOCATE (testcat.measure, Measure, NMEASURE_TESTCAT, Nmatchmeas, 1000); + } + } + testcat.average[0].Nmeasure = Nmatchmeas; + testcat.Nmeasure = Nmatchmeas; + populate_tiny_values (&testcat, DVO_TV_MEASURE); + + // we have now accumulated the measurements for this group, let's try a fit + // this needs to be a (fairly robust) clipped fit or we will have a hard time + // distinguishing a bad fit from a fit with 1 or 2 bad points + FIT_MODE = FIT_PM_ONLY; + UpdateObjects (&testcat, 1); + + // logic for keeping the fit (anything else?) + int good = FALSE; + good |= (testcat.average[0].ChiSqPM < MAX_CHISQ_PM); + good |= (testcat.average[0].Npos > MIN_NPOS_FRAC * testcat.average[0].Nmeasure); + good &= ((testcat.average[0].flags & ID_STAR_FIT_PM) > 0); + + if (good) { + // save the new object on catalogOut + // (note: there is only one object in testcat and thus measureOffset = 0 + catalogOut.average[catalogOut.Naverage] = testcat.average[0]; + catalogOut.average[catalogOut.Naverage].measureOffset = catalogOut.Nmeasure; + for (k = 0; k < testcat.average[0].Nmeasure; k++) { + catalogOut.measure[catalogOut.Nmeasure] = testcat.measure[k]; + catalogOut.measure[catalogOut.Nmeasure].averef = catalogOut.Naverage; + catalogOut.Nmeasure ++; + CHECK_REALLOCATE (catalogOut.measure, Measure, NMEASURE, catalogOut.Nmeasure, 100); + } + catalogOut.Naverage ++; + CHECK_REALLOCATE (catalogOut.average, Average, NAVERAGE, catalogOut.Naverage, 100); + } + NEXT_I; + } + catalogOut.Nsecfilt = Nsecfilt; + catalogOut.Nsecf_mem = Nsecfilt * catalogOut.Naverage; + + populate_tiny_values (&catalogOut, DVO_TV_MEASURE); + UpdateObjects (&catalogOut, 1); + + fprintf (stderr, "found "OFF_T_FMT" matches with "OFF_T_FMT" detections\n", catalogOut.Naverage, catalogOut.Nmeasure); + dvo_catalog_save (&catalogOut, VERBOSE2); + dvo_catalog_unlock (&catalogOut); + dvo_catalog_free (&catalogOut); + free (slowMoving); + free (X1); + free (Y1); + free (N1); + + return TRUE; +} Index: /trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro.c (revision 35415) +++ /trunk/Ohana/src/relastro/src/relastro.c (revision 35416) @@ -19,4 +19,8 @@ /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */ high_speed_catalogs (sky, skylist, 0, NULL); + exit (0); + + case TARGET_HPM: + hpm_catalogs (sky, skylist, 0, NULL); exit (0); Index: /trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_client.c (revision 35415) +++ /trunk/Ohana/src/relastro/src/relastro_client.c (revision 35416) @@ -58,4 +58,10 @@ } + case TARGET_HPM: { + // USAGE: relastro_client -high-speed + hpm_catalogs (sky, skylist, HOST_ID, HOSTDIR); + break; + } + // XXX loading the images is fairly costly -- see if we can do an image subset case TARGET_UPDATE_OFFSETS: { Index: /trunk/Ohana/src/relphot/Makefile =================================================================== --- /trunk/Ohana/src/relphot/Makefile (revision 35415) +++ /trunk/Ohana/src/relphot/Makefile (revision 35416) @@ -40,8 +40,10 @@ $(SRC)/load_catalogs.$(ARCH).o \ $(SRC)/load_images.$(ARCH).o \ +$(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plot_scatter.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/reload_catalogs.$(ARCH).o \ $(SRC)/relphot.$(ARCH).o \ +$(SRC)/relphot_images.$(ARCH).o \ $(SRC)/relphot_objects.$(ARCH).o \ $(SRC)/select_images.$(ARCH).o \ @@ -69,4 +71,5 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/help.$(ARCH).o \ +$(SRC)/synthetic_mags.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/liststats.$(ARCH).o \ Index: /trunk/Ohana/src/relphot/include/relphot.h =================================================================== --- /trunk/Ohana/src/relphot/include/relphot.h (revision 35415) +++ /trunk/Ohana/src/relphot/include/relphot.h (revision 35416) @@ -4,4 +4,13 @@ # include # include + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); /* # define GRID_V1 */ @@ -193,4 +202,5 @@ int UpdateAverages; int ApplyOffsets; +int SyntheticPhotometry; char *PhotcodeList; @@ -353,4 +363,5 @@ void write_coords PROTO((Header *header, Coords *coords)); int relphot_objects (int hostID, char *hostpath); +int relphot_images (void); void relphot_usage (void); @@ -390,2 +401,5 @@ void SetMrelInfoInit (SetMrelInfo *results, int allocLists); + +int init_synthetic_mags (); +int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm); Index: /trunk/Ohana/src/relphot/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/MosaicOps.c (revision 35416) @@ -152,10 +152,4 @@ } -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stopTimer, (void *) NULL); \ - dtime = DTIME (stopTimer, startTimer); \ - fprintf (stderr, MSG, __VA_ARGS__); } - /* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...) this function will also identify the images NOT in the subset which belong to a selected mosaic @@ -167,9 +161,7 @@ char *pname; - struct timeval startTimer, stopTimer; - if (!MOSAIC_ZEROPT) return; - gettimeofday (&startTimer, (void *) NULL); + INITTIME; /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no @@ -264,5 +256,5 @@ } } - fprintf (stderr, "%d total images, %d overlap skyregion, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark); + fprintf (stderr, "%d total images, %d overlap skyregion & match selection criteria, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark); MARKTIME("find unselected images matching selected mosaics: %f sec\n", dtime); Index: /trunk/Ohana/src/relphot/src/args.c =================================================================== --- /trunk/Ohana/src/relphot/src/args.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/args.c (revision 35416) @@ -251,4 +251,8 @@ remove_argument (N, &argc, argv); MOSAIC_ZEROPT = TRUE; + if (!strcasecmp (MOSAICNAME, "none")) { + fprintf (stderr, "mosaic astrometry selected by MOSAICNAME not defined\n"); + exit (2); + } } @@ -315,4 +319,11 @@ remove_argument (N, &argc, argv); DophotSelect = TRUE; + } + + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = TRUE; + init_synthetic_mags(); } @@ -516,4 +527,11 @@ } + SyntheticPhotometry = FALSE; + if ((N = get_argument (argc, argv, "-synthphot"))) { + remove_argument (N, &argc, argv); + SyntheticPhotometry = FALSE; + init_synthetic_mags(); + } + ImagSelect = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { Index: /trunk/Ohana/src/relphot/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relphot/src/bcatalog.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/bcatalog.c (revision 35416) @@ -179,4 +179,5 @@ REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); } + myAssert (Nmeasure < NMEASURE, "realloc failure"); } @@ -204,4 +205,17 @@ continue; } + + // for w-band photometry (& other cases?) convert gri photometry + // to a synthetic w-band magnitude + add_synthetic_mags (&subcatalog[0].averageT[Naverage], + &subcatalog[0].secfilt[Nsecfilt*Naverage], + &subcatalog[0].measureT[Nmeasure], + &Nmeasure, + &Nm); + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + } + myAssert (Nmeasure < NMEASURE, "realloc failure"); subcatalog[0].averageT[Naverage].Nmeasure = Nm; Index: /trunk/Ohana/src/relphot/src/global_stats.c =================================================================== --- /trunk/Ohana/src/relphot/src/global_stats.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/global_stats.c (revision 35416) @@ -1,9 +1,3 @@ # include "relphot.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stopTimer, (void *) NULL); \ - dtime = DTIME (stopTimer, startTimer); \ - fprintf (stderr, MSG, __VA_ARGS__); } void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) { @@ -11,10 +5,8 @@ StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; - // struct timeval startTimer, stopTimer; + // INITTIME; fprintf (stderr, "\n"); fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); - - // gettimeofday (&startTimer, (void *) NULL); int Ns; Index: /trunk/Ohana/src/relphot/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/load_catalogs.c (revision 35416) @@ -148,12 +148,14 @@ char tmpline[1024]; - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } - if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (RESET_ZEROPTS) { snprintf (tmpline, 1024, "%s -reset-zpts", command); strcpy (command, tmpline); } - if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } - if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } - if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } - if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (RESET_ZEROPTS) { snprintf (tmpline, 1024, "%s -reset-zpts", command); strcpy (command, tmpline); } + if (!KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -reset-ubercal", command); strcpy (command, tmpline); } + if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (SyntheticPhotometry) { snprintf (tmpline, 1024, "%s -synthphot", command); strcpy (command, tmpline); } + if (TimeSelect) { char *tstart = ohana_sec_to_date (TSTART); Index: /trunk/Ohana/src/relphot/src/load_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/load_images.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/load_images.c (revision 35416) @@ -1,9 +1,3 @@ # include "relphot.h" - -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } // This function generates a subset of the images based on selections. Input db has already @@ -14,5 +8,4 @@ off_t Nimage, Nsubset, Nchar; off_t *LineNumber; - struct timeval start, stop; char *inSubset; @@ -20,5 +13,5 @@ SkyList *skylist = NULL; - gettimeofday (&start, (void *) NULL); + INITTIME; // load the current sky table (layout of all SkyRegions) @@ -93,5 +86,5 @@ } -# if (1) +# if (0) FILE *ftest = fopen ("skydump.dat", "w"); int i; Index: /trunk/Ohana/src/relphot/src/relphot.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/relphot.c (revision 35416) @@ -1,19 +1,5 @@ # include "relphot.h" -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - int main (int argc, char **argv) { - - int i, status, Ncatalog; - Catalog *catalog = NULL; - FITS_DB db; - struct timeval start, stop; - SkyList *skylist = NULL; - - gettimeofday (&start, (void *) NULL); /* get configuration info, args */ @@ -36,242 +22,5 @@ } - /* register database handle with shutdown procedure */ - set_db (&db); - db.mode = dvo_catalog_catmode (CATMODE); - db.format = dvo_catalog_catformat (CATFORMAT); - - /* lock and load the image db table */ - status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); - if (!status && UPDATE) { - fprintf (stderr, "error\n"); - Shutdown ("ERROR: failure to lock image catalog %s", db.filename); - } - // if the file is missing, db.dbstate will have a value of either: - // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) - if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { - - // 25.0 is the nominal zero point for measurements in the database - // for all measurements except PHOT_REF: - // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 - // for measurements using PHOT_REF: - // measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known) - dvo_image_create (&db, 25.0); - - // XXX why are we running relphot on an empty database? - } else { - if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); - } - MARKTIME("-- load image data: %f sec\n", dtime); - - /* load regions and images based on specified sky patch and/or catalog */ - skylist = load_images (&db, UserCatalog, &UserPatch); - MARKTIME("-- load images: %f sec\n", dtime); - - /* unlock, if we can (else, unlocked below) */ - if (!UPDATE) dvo_image_unlock (&db); - - // load the flat correction table (if defined) - char flatcorrfile[256]; - sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); - - if (NLOOP > 0) { - /* load catalog data from region files (hostID is 0 since we are not a client */ - catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); - MARKTIME("-- load catalog data: %f sec\n", dtime); - - /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ - - /* match measurements with images, mosaics */ - initImageBins (catalog, Ncatalog, TRUE); - MARKTIME("-- make image bins: %f sec\n", dtime); - - initMosaicBins (catalog, Ncatalog, TRUE); - initGridBins (catalog, Ncatalog); - initMrel (catalog, Ncatalog); - - findImages (catalog, Ncatalog, TRUE); - MARKTIME("-- set up image indexes: %f sec\n", dtime); - - findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ - MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); - - SAVEPLOT = FALSE; - - setExclusions (catalog, Ncatalog, TRUE); - - global_stats (catalog, Ncatalog, flatcorr); - - if (PLOTSTUFF) { - plot_star_coords (catalog, Ncatalog); - // plot_mosaic_fields (catalog); - } - - // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: - if (USE_GRID) { - int star_toofew; - - // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW - // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel - star_toofew = STAR_TOOFEW; - STAR_TOOFEW = 0; - for (i = 0; i < NGRID; i++) { - STAR_BAD = ID_STAR_POOR; - setMrel (catalog, Ncatalog, flatcorr); - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - setMgrid (catalog, flatcorr); - } - STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; - STAR_TOOFEW = star_toofew; - } - - /* determine fit values */ - for (i = 0; i < NLOOP; i++) { - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - setMrel (catalog, Ncatalog, flatcorr); // threaded - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - } - setMcal (catalog, FALSE, flatcorr); - setMmos (catalog, FALSE, flatcorr); - setMgrid (catalog, flatcorr); - MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); - - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - } - // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); - // if (i % 6 == 1) rationalize_images (); - - // NOTE : in the past, I was not iterating enough before cleaning. make sure we do - // at least 8 loops first -- that should get the systematic errors down to the ~1% - // level, even in cases where we have an even split between photometric data and - // data with 1 mag of extinction. - - if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); - if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); - if ((i > 8) && (i % 8 == 5)) clean_mosaics (); - if ((i > 8) && (i % 8 == 5)) clean_images (); - - // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); - // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); - // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); - // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); - if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); - MARKTIME("-- finished loop %d: %f sec\n", i, dtime); - } - - if (PLOTSTUFF) { - plot_scatter (catalog, Ncatalog, flatcorr); - plot_grid (catalog, flatcorr); - plot_mosaics (); - plot_images (); - plot_stars (catalog, Ncatalog); - plot_chisq (catalog, Ncatalog); - } - - if (USE_GRID) dump_grid (); - - /* set Mcal & Mmos for bad images */ - setMcal (catalog, TRUE, flatcorr); - setMmos (catalog, TRUE, flatcorr); - MARKTIME("-- finalize Mcal values: %f sec\n", dtime); - - setMcalFinal (); // copy per-mosaic calibrations to the images - - if (SAVE_IMAGE_UPDATES) { - - FITS_DB dbX; - off_t *LineNumber; - Image *image, *subset; - off_t Nimage, Nsubset, i, Nx; - char filename[1024]; - - gfits_db_init (&dbX); - dbX.lockstate = LCK_XCLD; - dbX.timeout = 60.0; - dbX.mode = db.mode; - dbX.format = db.format; - - snprintf (filename, 1024, "%s.bck", ImageCat); - if (!gfits_db_lock (&dbX, filename)) { - fprintf (stderr, "can't lock backup image image catalog\n"); - return (FALSE); - } - - // apply the changes from the image subset to the full table: - - // convert database table to internal structure (binary to Image) - // 'image' points to the same memory as db->ftable->buffer - image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); - if (!image) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); - } - gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); - - subset = getimages (&Nsubset, &LineNumber); - for (i = 0; i < Nsubset; i++) { - if (LineNumber[i] == -1) continue; - memcpy (&image[LineNumber[i]], &subset[i], Nx); - } - - // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump - // the entire Image table, but with the updates in place I think this says: re-work the - // ftable/vtable usage in this program to be more sensible... - gfits_copy_header (&db.header, &dbX.header); - gfits_copy_matrix (&db.matrix, &dbX.matrix); - gfits_copy_header (&db.theader, &dbX.theader); - gfits_copy_ftable (&db.ftable, &dbX.ftable); - - dbX.ftable.header = &dbX.theader; - dbX.virtual = FALSE; - - // save Images.dat using the copied structure - if (UPDATE_CATFORMAT) { - // ensure the db format is updated - dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); - char photcodeFile[1024]; - sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); - SavePhotcodesFITS (photcodeFile); - } - dvo_image_save (&dbX, VERBOSE); - dvo_image_unlock (&dbX); - MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); - } - - /* at this point, we have correct cal coeffs in the image/mosaic structures */ - for (i = 0; i < Ncatalog; i++) { - free_tiny_values (&catalog[i]); - dvo_catalog_free (&catalog[i]); - } - freeImageBins (Ncatalog, TRUE); - freeMosaicBins (Ncatalog, TRUE); - freeGridBins (Ncatalog); - - } // (NLOOP > 0) : this loop determines the offsets per chip - - reload_images (&db); - - // only change the real database files if -update is requested - if (!UPDATE) exit (0); - - /* load catalog data from region files, update Mrel include all data */ - reload_catalogs (skylist, flatcorr, 0, NULL); - MARKTIME("-- updated all catalogs: %f sec\n", dtime); - - if (UPDATE_CATFORMAT) { - // ensure the db format is updated - db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); - } - dvo_image_update (&db, VERBOSE); - dvo_image_unlock (&db); + relphot_images (); exit (0); Index: /trunk/Ohana/src/relphot/src/relphot_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/relphot_images.c (revision 35416) +++ /trunk/Ohana/src/relphot/src/relphot_images.c (revision 35416) @@ -0,0 +1,253 @@ +# include "relphot.h" + +// 25.0 is the nominal zero point for measurements in the database +// for all measurements except PHOT_REF: +// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + 25.0 +// for measurements using PHOT_REF: +// measure.M contains: -2.5*log(DN) + 2.5*log(exptime) + True Zero Point (but exptime and ZP are known) + +int relphot_images () { + + int i, status, Ncatalog; + Catalog *catalog = NULL; + FITS_DB db; + SkyList *skylist = NULL; + + INITTIME; + + /* register database handle with shutdown procedure */ + set_db (&db); + db.mode = dvo_catalog_catmode (CATMODE); + db.format = dvo_catalog_catformat (CATFORMAT); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status && UPDATE) { + fprintf (stderr, "error\n"); + Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + } + + // if the file is missing, db.dbstate will have a value of either: + // LCK_EMPTY (if UPDATE) or LCK_MISSING (if !UPDATE) + if ((db.dbstate == LCK_EMPTY) || (db.dbstate == LCK_MISSING)) { + Shutdown ("ERROR: database %s contains no image data", CATDIR); + } + + // read data from Image.dat file + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("-- load image data: %f sec\n", dtime); + + /* load regions and images based on specified sky patch and/or catalog */ + skylist = load_images (&db, UserCatalog, &UserPatch); + MARKTIME("-- load images: %f sec\n", dtime); + + /* unlock, if we can (else, unlocked below) */ + if (!UPDATE) dvo_image_unlock (&db); + + // load the flat correction table (if defined) + char flatcorrfile[256]; + sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR); + FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE); + + if (NLOOP > 0) { + /* load catalog data from region files (hostID is 0 since we are not a client */ + catalog = load_catalogs (skylist, &Ncatalog, 0, NULL); + MARKTIME("-- load catalog data: %f sec\n", dtime); + + /* add in a loop over the catalogs calling dvo_catalog_chipcoords */ + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("-- make image bins: %f sec\n", dtime); + + initMosaicBins (catalog, Ncatalog, TRUE); + initGridBins (catalog, Ncatalog); + initMrel (catalog, Ncatalog); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("-- set up image indexes: %f sec\n", dtime); + + findMosaics (catalog, Ncatalog, TRUE); /* also sets Grid values */ + MARKTIME("-- set up mosaic indexes: %f sec\n", dtime); + + SAVEPLOT = FALSE; + + setExclusions (catalog, Ncatalog, TRUE); + + global_stats (catalog, Ncatalog, flatcorr); + + if (PLOTSTUFF) { + plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // if we are measuring the flat-field correction grid, we need to perform a number of iterations first: + if (USE_GRID) { + int star_toofew; + + // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW + // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel + star_toofew = STAR_TOOFEW; + STAR_TOOFEW = 0; + for (i = 0; i < NGRID; i++) { + STAR_BAD = ID_STAR_POOR; + setMrel (catalog, Ncatalog, flatcorr); + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + setMgrid (catalog, flatcorr); + } + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + STAR_TOOFEW = star_toofew; + } + + /* determine fit values */ + for (i = 0; i < NLOOP; i++) { + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMrel (catalog, Ncatalog, flatcorr); // threaded + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + } + setMcal (catalog, FALSE, flatcorr); + setMmos (catalog, FALSE, flatcorr); + setMgrid (catalog, flatcorr); + MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime); + + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + // if (i < NLOOP - 1) rationalize_mosaics (catalog, Ncatalog); + // if (i % 6 == 1) rationalize_images (); + + // NOTE : in the past, I was not iterating enough before cleaning. make sure we do + // at least 8 loops first -- that should get the systematic errors down to the ~1% + // level, even in cases where we have an even split between photometric data and + // data with 1 mag of extinction. + + if ((i > 8) && (i % 8 == 2)) clean_measures (catalog, Ncatalog, FALSE, flatcorr); + if ((i > 8) && (i % 8 == 3)) clean_stars (catalog, Ncatalog); + if ((i > 8) && (i % 8 == 5)) clean_mosaics (); + if ((i > 8) && (i % 8 == 5)) clean_images (); + + // if ((i == 1) || (i == 5) || (i == 9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); + // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics (); + // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images (); + if (i % 3 == 2) global_stats (catalog, Ncatalog, flatcorr); + MARKTIME("-- finished loop %d: %f sec\n", i, dtime); + } + + if (PLOTSTUFF) { + plot_scatter (catalog, Ncatalog, flatcorr); + plot_grid (catalog, flatcorr); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + } + + if (USE_GRID) dump_grid (); + + /* set Mcal & Mmos for bad images */ + setMcal (catalog, TRUE, flatcorr); + setMmos (catalog, TRUE, flatcorr); + MARKTIME("-- finalize Mcal values: %f sec\n", dtime); + + setMcalFinal (); // copy per-mosaic calibrations to the images + + if (SAVE_IMAGE_UPDATES) { + + FITS_DB dbX; + off_t *LineNumber; + Image *image, *subset; + off_t Nimage, Nsubset, i, Nx; + char filename[1024]; + + gfits_db_init (&dbX); + dbX.lockstate = LCK_XCLD; + dbX.timeout = 60.0; + dbX.mode = db.mode; + dbX.format = db.format; + + snprintf (filename, 1024, "%s.bck", ImageCat); + if (!gfits_db_lock (&dbX, filename)) { + fprintf (stderr, "can't lock backup image image catalog\n"); + return (FALSE); + } + + // apply the changes from the image subset to the full table: + + // convert database table to internal structure (binary to Image) + // 'image' points to the same memory as db->ftable->buffer + image = gfits_table_get_Image (&db.ftable, &Nimage, &db.swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + gfits_scan (db.ftable.header, "NAXIS1", OFF_T_FMT, 1, &Nx); + + subset = getimages (&Nsubset, &LineNumber); + for (i = 0; i < Nsubset; i++) { + if (LineNumber[i] == -1) continue; + memcpy (&image[LineNumber[i]], &subset[i], Nx); + } + + // copy Images.dat data (all or only the subset / vtable elements?) I think I need to dump + // the entire Image table, but with the updates in place I think this says: re-work the + // ftable/vtable usage in this program to be more sensible... + gfits_copy_header (&db.header, &dbX.header); + gfits_copy_matrix (&db.matrix, &dbX.matrix); + gfits_copy_header (&db.theader, &dbX.theader); + gfits_copy_ftable (&db.ftable, &dbX.ftable); + + dbX.ftable.header = &dbX.theader; + dbX.virtual = FALSE; + + // save Images.dat using the copied structure + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + dbX.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + char photcodeFile[1024]; + sprintf (photcodeFile, "%s/Photcodes.dat", CATDIR); + SavePhotcodesFITS (photcodeFile); + } + dvo_image_save (&dbX, VERBOSE); + dvo_image_unlock (&dbX); + MARKTIME("-- save Image.dat.bck: %f sec\n", dtime); + } + + /* at this point, we have correct cal coeffs in the image/mosaic structures */ + for (i = 0; i < Ncatalog; i++) { + free_tiny_values (&catalog[i]); + dvo_catalog_free (&catalog[i]); + } + freeImageBins (Ncatalog, TRUE); + freeMosaicBins (Ncatalog, TRUE); + freeGridBins (Ncatalog); + + } // (NLOOP > 0) : this loop determines the offsets per chip + + reload_images (&db); + + // only change the real database files if -update is requested + if (!UPDATE) exit (0); + + /* load catalog data from region files, update Mrel include all data */ + reload_catalogs (skylist, flatcorr, 0, NULL); + MARKTIME("-- updated all catalogs: %f sec\n", dtime); + + if (UPDATE_CATFORMAT) { + // ensure the db format is updated + db.format = dvo_catalog_catformat (UPDATE_CATFORMAT); + } + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + + exit (0); +} + Index: /trunk/Ohana/src/relphot/src/select_images.c =================================================================== --- /trunk/Ohana/src/relphot/src/select_images.c (revision 35415) +++ /trunk/Ohana/src/relphot/src/select_images.c (revision 35416) @@ -17,10 +17,4 @@ off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); -# define MARKTIME(MSG,...) { \ - float dtime; \ - gettimeofday (&stop, (void *) NULL); \ - dtime = DTIME (stop, start); \ - fprintf (stderr, MSG, __VA_ARGS__); } - Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, char *inSubset, off_t **LineNumber, off_t *Nimage, SkyRegion *region) { @@ -31,5 +25,4 @@ double Ri[5], Di[5], Xi[5], Yi[5]; Coords tcoords; - struct timeval start, stop; if (skylist[0].Nregions < 1) { @@ -40,5 +33,5 @@ } - gettimeofday (&start, (void *) NULL); + INITTIME; // the comparison is made in the catalog local projection. below we set crval1,2 @@ -82,5 +75,5 @@ } -# if (1) +# if (0) // XXX quick test of the dec bands: FILE *fout = fopen ("dec.bands.dat", "w"); @@ -132,5 +125,5 @@ // this adds 1.3 sec for 3M images if (!FindMosaicForImage (timage, Ntimage, i)) { - fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); continue; } Index: /trunk/Ohana/src/relphot/src/synthetic_mags.c =================================================================== --- /trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 35416) +++ /trunk/Ohana/src/relphot/src/synthetic_mags.c (revision 35416) @@ -0,0 +1,71 @@ +# include "relphot.h" + +static int NSr = -1; +static int NSi = -1; +static int W_SYNTH = -1; + +int init_synthetic_mags () { + int Pr = GetPhotcodeCodebyName ("r"); + int Pi = GetPhotcodeCodebyName ("i"); + + NSr = GetPhotcodeNsec (Pr); + NSi = GetPhotcodeNsec (Pi); + + W_SYNTH = GetPhotcodeCodebyName ("SYNTH.w"); + + return TRUE; +} + +int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm) { + + if (!SyntheticPhotometry) return TRUE; + + myAssert (NSr >= 0, "failed to call init_synthetic_mags"); + myAssert (NSi >= 0, "failed to call init_synthetic_mags"); + myAssert (W_SYNTH >= 0, "failed to call init_synthetic_mags"); + + // XXX for now, hard-wire the photometric relationship + + // we want to create a synthetic w-band measurement based on the observed mean r & i photometry + + // from Tonry + // (w-r)_obs = 0.042 + 0.166 (r-i)_obs - 0.398 (r-i)_obs^2, (r-i)_obs < 0.5 + // (w-r)_obs = 0.268 - 0.435 (r-i)_obs - 0.078 (r-i)_obs^2, (r-i)_obs > 0.5 + // thus, for (r-i < 0.5): + // (r - i < 0.5) : w = r + 0.042 + 0.166*(r-i) - 0.398(r-i)^2 + // (r - i > 0.5) : w = r + 0.268 - 0.435 (r-i) - 0.078(r-i)^2 + + float Mr = secfilt[NSr].M; + float Mi = secfilt[NSi].M; + + if (!isfinite(Mr)) return FALSE; + if (!isfinite(Mi)) return FALSE; + + // apply a color correction + // XXX this is very GPC1 specific and hard-wired -- be very afraid! + // select the filter; default to fixed photcode and mag limit otherwise + + float Mri = Mr - Mi; + + // saturate at some valid range limits + Mri = MAX (MIN(Mri, 2.0), -0.2); + + float Mw = NAN; + if (Mri < 0.5) { + Mw = Mr + 0.042 + 0.166*(Mri) - 0.398*SQ(Mri); + } else { + Mw = Mr + 0.268 - 0.435*(Mri) - 0.078*SQ(Mri); + } + + dvo_measureT_init (measure); + + measure->M = Mw; + measure->dM = 0.001; + measure->photcode = W_SYNTH; + measure->dbFlags |= ID_MEAS_SYNTH_MAG; + + (*Nmeasure) ++; + (*Nm) ++; + + return TRUE; +} Index: /trunk/Ohana/src/uniphot/src/match_flatcorr_to_images.c =================================================================== --- /trunk/Ohana/src/uniphot/src/match_flatcorr_to_images.c (revision 35415) +++ /trunk/Ohana/src/uniphot/src/match_flatcorr_to_images.c (revision 35416) @@ -43,4 +43,5 @@ for (i = 0; i < Nimage; i++) { if (!image[i].photcode) continue; // skip PHU images + if (image[i].photcode > maxCode) continue; // skip PHU images int found = FALSE;