Index: trunk/Ohana/src/addstar/Makefile =================================================================== --- trunk/Ohana/src/addstar/Makefile (revision 37729) +++ trunk/Ohana/src/addstar/Makefile (revision 37807) @@ -24,4 +24,8 @@ sedstar : $(BIN)/sedstar.$(ARCH) load2mass : $(BIN)/load2mass.$(ARCH) +loadstarpar : $(BIN)/loadstarpar.$(ARCH) +loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH) +loadICRF : $(BIN)/loadICRF.$(ARCH) +loadICRF_client : $(BIN)/loadICRF_client.$(ARCH) loadwise : $(BIN)/loadwise.$(ARCH) dumpskycells : $(BIN)/dumpskycells.$(ARCH) @@ -31,7 +35,7 @@ mkcmf : $(BIN)/mkcmf.$(ARCH) -all: addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells findskycell - -INSTALL = addstar addstar_client sedstar load2mass skycells mkcmf loadwise loadsupercos dumpskycells findskycell +all: addstar addstar_client sedstar load2mass loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadsupercos dumpskycells findskycell + +INSTALL = addstar addstar_client sedstar load2mass loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf loadwise loadsupercos dumpskycells findskycell # I need to fix the client/server version of addstar now that I have dropped Stars @@ -241,4 +245,72 @@ $(SRC)/psps_ids.$(ARCH).o +LOAD-STARPAR = \ +$(SRC)/loadstarpar.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/args_loadstarpar.$(ARCH).o \ +$(SRC)/find_matches_starpar.$(ARCH).o \ +$(SRC)/loadstarpar_catalog.$(ARCH).o \ +$(SRC)/loadstarpar_io.$(ARCH).o \ +$(SRC)/loadstarpar_make_subset.$(ARCH).o \ +$(SRC)/loadstarpar_readstars.$(ARCH).o \ +$(SRC)/loadstarpar_remote_hosts.$(ARCH).o \ +$(SRC)/loadstarpar_save_remote.$(ARCH).o \ +$(SRC)/loadstarpar_table.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + +LOAD-STARPAR-CLIENT = \ +$(SRC)/loadstarpar_client.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/args_loadstarpar.$(ARCH).o \ +$(SRC)/find_matches_starpar.$(ARCH).o \ +$(SRC)/loadstarpar_catalog.$(ARCH).o \ +$(SRC)/loadstarpar_io.$(ARCH).o \ +$(SRC)/loadstarpar_make_subset.$(ARCH).o \ +$(SRC)/loadstarpar_readstars.$(ARCH).o \ +$(SRC)/loadstarpar_remote_hosts.$(ARCH).o \ +$(SRC)/loadstarpar_save_remote.$(ARCH).o \ +$(SRC)/loadstarpar_table.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + +LOAD-ICRF = \ +$(SRC)/loadICRF.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/args_loadICRF.$(ARCH).o \ +$(SRC)/find_matches_ICRF.$(ARCH).o \ +$(SRC)/loadICRF_catalog.$(ARCH).o \ +$(SRC)/loadICRF_io.$(ARCH).o \ +$(SRC)/loadICRF_make_subset.$(ARCH).o \ +$(SRC)/loadICRF_readstars.$(ARCH).o \ +$(SRC)/loadICRF_remote_hosts.$(ARCH).o \ +$(SRC)/loadICRF_save_remote.$(ARCH).o \ +$(SRC)/loadICRF_table.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + +LOAD-ICRF-CLIENT = \ +$(SRC)/loadICRF_client.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/SkyRegionUtils.$(ARCH).o \ +$(SRC)/args_loadICRF.$(ARCH).o \ +$(SRC)/find_matches_ICRF.$(ARCH).o \ +$(SRC)/loadICRF_catalog.$(ARCH).o \ +$(SRC)/loadICRF_io.$(ARCH).o \ +$(SRC)/loadICRF_make_subset.$(ARCH).o \ +$(SRC)/loadICRF_readstars.$(ARCH).o \ +$(SRC)/loadICRF_remote_hosts.$(ARCH).o \ +$(SRC)/loadICRF_save_remote.$(ARCH).o \ +$(SRC)/loadICRF_table.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + LOAD-WISE = \ $(SRC)/loadwise.$(ARCH).o \ @@ -309,15 +381,19 @@ $(SRC)/Shutdown.$(ARCH).o -$(ADDSTARC) : $(INC)/addstar.h -$(ADDSTARD) : $(INC)/addstar.h -$(ADDSTART) : $(INC)/addstar.h -$(ADDSTAR) : $(INC)/addstar.h -$(ADDSTAR_CLIENT) : $(INC)/addstar.h -$(SKYCELLS) : $(INC)/addstar.h -$(DUMPSKYCELLS) : $(INC)/addstar.h -$(FINDSKYCELL) : $(INC)/addstar.h -$(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h -$(LOAD-WISE) : $(INC)/addstar.h $(INC)/WISE.h -$(MKCMF) : $(INC)/mkcmf.h +$(ADDSTARC) : $(INC)/addstar.h +$(ADDSTARD) : $(INC)/addstar.h +$(ADDSTART) : $(INC)/addstar.h +$(ADDSTAR) : $(INC)/addstar.h +$(ADDSTAR_CLIENT) : $(INC)/addstar.h +$(SKYCELLS) : $(INC)/addstar.h +$(DUMPSKYCELLS) : $(INC)/addstar.h +$(FINDSKYCELL) : $(INC)/addstar.h +$(LOAD-2MASS) : $(INC)/addstar.h $(INC)/2mass.h +$(LOAD-STARPAR) : $(INC)/addstar.h $(INC)/loadstarpar.h +$(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h +$(LOAD-ICRF) : $(INC)/addstar.h $(INC)/loadICRF.h +$(LOAD-ICRF-CLIENT) : $(INC)/addstar.h $(INC)/loadICRF.h +$(LOAD-WISE) : $(INC)/addstar.h $(INC)/WISE.h +$(MKCMF) : $(INC)/mkcmf.h $(BIN)/addstar.$(ARCH) : $(ADDSTAR) @@ -330,4 +406,8 @@ $(BIN)/load2mass.$(ARCH) : $(LOAD-2MASS) $(BIN)/loadwise.$(ARCH) : $(LOAD-WISE) +$(BIN)/loadstarpar.$(ARCH) : $(LOAD-STARPAR) +$(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT) +$(BIN)/loadICRF.$(ARCH) : $(LOAD-ICRF) +$(BIN)/loadICRF_client.$(ARCH): $(LOAD-ICRF-CLIENT) $(BIN)/loadsupercos.$(ARCH) : $(LOAD-SUPERCOS) $(BIN)/skycells.$(ARCH) : $(SKYCELLS) Index: trunk/Ohana/src/addstar/include/WISE.h =================================================================== --- trunk/Ohana/src/addstar/include/WISE.h (revision 37729) +++ trunk/Ohana/src/addstar/include/WISE.h (revision 37807) @@ -21,9 +21,10 @@ short WISE_W1, WISE_W2, WISE_W3, WISE_W4; -int ALLSKY; +enum {MODE_NONE, MODE_ALLWISE, MODE_ALLSKY, MODE_PRELIM}; +int MODE; AddstarClientOptions args_loadwise (int argc, char **argv, AddstarClientOptions options); -int loadwise_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options, int ALLSKY); +int loadwise_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options); int getWISE_setup (); Index: trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- trunk/Ohana/src/addstar/include/addstar.h (revision 37729) +++ trunk/Ohana/src/addstar/include/addstar.h (revision 37807) @@ -20,5 +20,5 @@ # define IN_REGION(R,D) ( \ ((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \ -((R) >= region[0].Rmin) && ((R) < region[0].Rmax)) +((R) >= region[0].Rmin) && ((R) < region[0].Rmax)) /* grab named photcode */ @@ -31,4 +31,6 @@ # define dCOS(A) ((double) cos ((double)RAD_DEG*A)) # define dSIN(A) ((double) sin ((double)RAD_DEG*A)) + +# define myAbortF(FORMAT,...) { fprintf (stderr, FORMAT, __VA_ARGS__); abort(); } typedef struct { @@ -326,4 +328,8 @@ int strextend (char *input, char *format,...); +void initMosaicCoords (); +void saveMosaicCoords (Coords *input); + + /** there is an inconsistency to be resolved: fixed structures (like Image) Index: trunk/Ohana/src/addstar/include/loadICRF.h =================================================================== --- trunk/Ohana/src/addstar/include/loadICRF.h (revision 37807) +++ trunk/Ohana/src/addstar/include/loadICRF.h (revision 37807) @@ -0,0 +1,42 @@ + +typedef struct { + double R, D; + Measure measure; + int found; // TRUE == matched to an input star + int flag; // TRUE == included in a subset +} ICRF_Stars; + +int HOST_ID; +char *HOSTDIR; +char *CPT_FILE; +char *INPUT; + +AddstarClientOptions args_loadICRF (int *argc, char **argv, AddstarClientOptions options); +AddstarClientOptions args_loadICRF_client (int *argc, char **argv, AddstarClientOptions options); + +int loadICRF_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options); + +ICRF_Stars *loadICRF_make_subset (ICRF_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset); + +int loadICRF_save_remote (ICRF_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options); + +int save_remote_host (HostInfo *host); + +int init_remote_hosts (); +int find_empty_slot (); +int harvest_all (); +int harvest_host (); + +int loadICRF_catalog (ICRF_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options); + +int galactic_to_celestial (double *R, double *D, double l, double b); + +int find_matches_ICRF (SkyRegion *region, ICRF_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options); + +int loadICRF_save_stars (char *filename, ICRF_Stars *stars, int Nstars); +ICRF_Stars *loadICRF_load_stars (char *filename, int *nstars); + +ICRF_Stars *loadICRF_readstars (char *filename, int *nstars); + +int loadICRF_sortStars (ICRF_Stars *stars, int Nstars); +int InitICRF_Star (ICRF_Stars *star); Index: trunk/Ohana/src/addstar/include/loadstarpar.h =================================================================== --- trunk/Ohana/src/addstar/include/loadstarpar.h (revision 37807) +++ trunk/Ohana/src/addstar/include/loadstarpar.h (revision 37807) @@ -0,0 +1,41 @@ + +typedef struct { + double R, D; + StarPar starpar; + int flag; // in a subset? + int found; // assigned to an object? +} StarPar_Stars; + +int HOST_ID; +char *HOSTDIR; +char *CPT_FILE; +char *INPUT; + +AddstarClientOptions args_loadstarpar (int argc, char **argv, AddstarClientOptions options); +AddstarClientOptions args_loadstarpar_client (int argc, char **argv, AddstarClientOptions options); + +int loadstarpar_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options); + +StarPar_Stars *loadstarpar_make_subset (StarPar_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset); + +int loadstarpar_save_remote (StarPar_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options); + +int save_remote_host (HostInfo *host); + +int init_remote_hosts (); +int find_empty_slot (); +int harvest_all (); +int harvest_host (); + +int loadstarpar_catalog (StarPar_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options); + +int galactic_to_celestial (double *R, double *D, double l, double b); + +int find_matches_starpar (SkyRegion *region, StarPar_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options); + +int loadstarpar_save_stars (char *filename, StarPar_Stars *stars, int Nstars); +StarPar_Stars *loadstarpar_load_stars (char *filename, int *nstars); + +StarPar_Stars *loadstarpar_readstars (char *filename, int *nstars); + +int loadstarpar_sortStars (StarPar_Stars *stars, int Nstars); Index: trunk/Ohana/src/addstar/src/LoadData.c =================================================================== --- trunk/Ohana/src/addstar/src/LoadData.c (revision 37729) +++ trunk/Ohana/src/addstar/src/LoadData.c (revision 37807) @@ -33,4 +33,10 @@ // we need to perform that analysis here GetZeroPointExposure (headers, headerSets, Nimages); + + if (!options[0].mosaic) { + // we are requiring a single set of N x WRPs + 1 DIS per file set. + // NOTE: if the -mosaic was supplied, do not reset the mosaic coords here + initMosaicCoords (); + } // now run through the images, interpret the headers and read the stars Index: trunk/Ohana/src/addstar/src/LoadDataPMM.c =================================================================== --- trunk/Ohana/src/addstar/src/LoadDataPMM.c (revision 37729) +++ trunk/Ohana/src/addstar/src/LoadDataPMM.c (revision 37807) @@ -193,5 +193,5 @@ // XXX for now, we define a totally fake coordinate system centered on the plate center - strcpy (image[0].coords.ctype, "DEC--TAN"); + InitCoords (&image[0].coords, "DEC--TAN"); image[0].coords.crval1 = pmm_get_ra (RA); @@ -203,12 +203,4 @@ image[0].coords.crpix2 = 0.5*PLATE_NY; image[0].coords.cdelt1 = image[0].coords.cdelt2 = 0.9 / 3600.0; - - image[0].coords.pc1_1 = 0.0; - image[0].coords.pc1_2 = 1.0; - image[0].coords.pc2_1 = 1.0; - image[0].coords.pc2_2 = 0.0; - - image[0].coords.Npolyterms = 0; - memset (image[0].coords.polyterms, 0, 2*7*sizeof(float)); image[0].NX = PLATE_NX; Index: trunk/Ohana/src/addstar/src/LoadStars.c =================================================================== --- trunk/Ohana/src/addstar/src/LoadStars.c (revision 37729) +++ trunk/Ohana/src/addstar/src/LoadStars.c (revision 37807) @@ -88,4 +88,9 @@ exit (1); } + /* supplied photcode is incompatible with multi-chip images */ + if ((NheaderSets > 1) && options[0].mosaic) { + fprintf (stderr, "ERROR: -mosaic cannot be supplied to multiple images\n"); + exit (1); + } // if these are SDSS data, load with SDSS-specific wrapper Index: trunk/Ohana/src/addstar/src/NewImage.c =================================================================== --- trunk/Ohana/src/addstar/src/NewImage.c (revision 37729) +++ trunk/Ohana/src/addstar/src/NewImage.c (revision 37807) @@ -33,4 +33,5 @@ goto reject; } + saveMosaicCoords (mosaic); } Index: trunk/Ohana/src/addstar/src/ReadImageHeader.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 37729) +++ trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 37807) @@ -1,3 +1,15 @@ # include "addstar.h" + +// for DIS/WRP sets, we need to save the DIS set and supply it to the WRP entries +static Coords *mosaic = NULL; + +void initMosaicCoords () { + mosaic = NULL; +} + +void saveMosaicCoords (Coords *input) { + ALLOCATE (mosaic, Coords, 1); + memcpy (mosaic, input, sizeof(Coords)); +} /* read an image header corresponding to a CMF / CMP data block */ @@ -54,5 +66,5 @@ /* only load astrometry, NAXIS1,2, and time if this is a MOSAIC_PHU (ctype is ....-DIS) */ if (!strcmp (&image[0].coords.ctype[4], "-DIS")) { - RegisterMosaic (&image[0].coords); + saveMosaicCoords (&image[0].coords); return (TRUE); } @@ -79,8 +91,9 @@ } if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { - if (!isRegisteredMosaic()) { + if (!mosaic) { fprintf (stderr, "no mosaic for WRP image (use -mosaic)\n"); return (FALSE); } + image[0].coords.mosaic = mosaic; } else { /* force image to lie in 0-360 range */ Index: trunk/Ohana/src/addstar/src/ReadStarsSDSS.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 37729) +++ trunk/Ohana/src/addstar/src/ReadStarsSDSS.c (revision 37807) @@ -201,20 +201,10 @@ // XXX for now, we define a totally fake coordinate system centered on the first listed star - strcpy (images[N].coords.ctype, "DEC--TAN"); - + InitCoords (&images[N].coords, "DEC--TAN"); images[N].coords.crval1 = stars[0].average.R; images[N].coords.crval2 = stars[0].average.D; - images[N].coords.crpix1 = stars[0].measure.Xccd; images[N].coords.crpix2 = stars[0].measure.Yccd; images[N].coords.cdelt1 = images[N].coords.cdelt2 = 0.4 / 3600.0; - - images[N].coords.pc1_1 = 0.0; - images[N].coords.pc1_2 = 1.0; - images[N].coords.pc2_1 = 1.0; - images[N].coords.pc2_2 = 0.0; - - images[N].coords.Npolyterms = 0; - memset (images[N].coords.polyterms, 0, 2*7*sizeof(float)); images[N].NX = 2048; Index: trunk/Ohana/src/addstar/src/UpdateDatabase_Image.c =================================================================== --- trunk/Ohana/src/addstar/src/UpdateDatabase_Image.c (revision 37729) +++ trunk/Ohana/src/addstar/src/UpdateDatabase_Image.c (revision 37807) @@ -21,5 +21,4 @@ /* find correpsonding regions for image */ - RegisterMosaic (mosaic); skylist = NULL; newlist = NULL; Index: trunk/Ohana/src/addstar/src/addstar.c =================================================================== --- trunk/Ohana/src/addstar/src/addstar.c (revision 37729) +++ trunk/Ohana/src/addstar/src/addstar.c (revision 37807) @@ -132,5 +132,5 @@ } } - MARKTIME ("load cpt: %f sec\n", dtime); RESETTIME; + if (VERBOSE) MARKTIME ("load cpt: %f sec\n", dtime); RESETTIME; // Naves_disk == 0 implies an empty catalog file @@ -163,5 +163,5 @@ break; } - MARKTIME ("match stars: %f sec\n", dtime); RESETTIME; + if (VERBOSE) MARKTIME ("match stars: %f sec\n", dtime); RESETTIME; /* report total updated values */ @@ -183,5 +183,5 @@ dvo_catalog_unlock (&catalog); dvo_catalog_free (&catalog); - MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; + if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME; if (options.mode == ADDSTAR_MODE_REFCAT) free (stars); Index: trunk/Ohana/src/addstar/src/args.c =================================================================== --- trunk/Ohana/src/addstar/src/args.c (revision 37729) +++ trunk/Ohana/src/addstar/src/args.c (revision 37807) @@ -175,5 +175,5 @@ exit (1); } - RegisterMosaic (&MOSAIC); + saveMosaicCoords (&MOSAIC); remove_argument (N, &argc, argv); gfits_free_header (&header); Index: trunk/Ohana/src/addstar/src/args_client.c =================================================================== --- trunk/Ohana/src/addstar/src/args_client.c (revision 37729) +++ trunk/Ohana/src/addstar/src/args_client.c (revision 37807) @@ -83,5 +83,5 @@ exit (1); } - RegisterMosaic (&MOSAIC); + saveMosaicCoords (&MOSAIC); remove_argument (N, &argc, argv); gfits_free_header (&header); Index: trunk/Ohana/src/addstar/src/args_loadICRF.c =================================================================== --- trunk/Ohana/src/addstar/src/args_loadICRF.c (revision 37807) +++ trunk/Ohana/src/addstar/src/args_loadICRF.c (revision 37807) @@ -0,0 +1,219 @@ +# include "addstar.h" +# include "loadICRF.h" + +static void help (void); +static void help_client (void); + +AddstarClientOptions args_loadICRF (int *argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + help (); + } + + /* basic mode: image, list, refcat */ + options.mode = ADDSTAR_MODE_REFCAT; + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + HOST_ID = 0; + 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); + } + + /* only add to existing regions */ + options.existing_regions = FALSE; + if ((N = get_argument (*argc, argv, "-existing-regions"))) { + options.existing_regions = TRUE; + remove_argument (N, argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (*argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + /* other addstar options which cannot be used in loadICRF */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (*argc < 2) { + fprintf (stderr, "USAGE: loadICRF [options] (file) [..more files]\n"); + exit (2); + } + return (options); +} + +AddstarClientOptions args_loadICRF_client (int *argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + help (); + } + + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + 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) help_client(); + + 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) help_client(); + + CPT_FILE = NULL; + if ((N = get_argument (*argc, argv, "-cpt"))) { + remove_argument (N, argc, argv); + CPT_FILE = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + if (!CPT_FILE) help_client(); + + INPUT = NULL; + if ((N = get_argument (*argc, argv, "-input"))) { + remove_argument (N, argc, argv); + INPUT = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + if (!INPUT) help_client(); + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + /* only add to existing objects */ + options.only_match = FALSE; + if ((N = get_argument (*argc, argv, "-only-match"))) { + options.only_match = TRUE; + remove_argument (N, argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (*argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + /* other addstar options which cannot be used in loadICRF */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (*argc != 1) { + fprintf (stderr, "USAGE: loadICRF_client -cpt (file) -input (file)\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadICRF [options] (file) [..more files]\n"); + fprintf (stderr, " add data from ICRF list to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} + +static void help_client () { + + fprintf (stderr, "USAGE: loadICRF_client -D CATDIR catdir -hostID ID -hostdir (dir) -cpt (filename) -input (input) [options]\n"); + fprintf (stderr, " add data from stellar parameter file to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} Index: trunk/Ohana/src/addstar/src/args_loadstarpar.c =================================================================== --- trunk/Ohana/src/addstar/src/args_loadstarpar.c (revision 37807) +++ trunk/Ohana/src/addstar/src/args_loadstarpar.c (revision 37807) @@ -0,0 +1,227 @@ +# include "addstar.h" +# include "loadstarpar.h" + +static void help (void); +static void help_client (void); + +AddstarClientOptions args_loadstarpar (int argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + /* basic mode: image, list, refcat */ + options.mode = ADDSTAR_MODE_REFCAT; + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + HOST_ID = 0; + 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); + } + + /* only add to existing regions */ + options.existing_regions = FALSE; + if ((N = get_argument (argc, argv, "-existing-regions"))) { + options.existing_regions = TRUE; + remove_argument (N, &argc, argv); + } + /* only add to existing objects */ + options.only_match = FALSE; + if ((N = get_argument (argc, argv, "-only-match"))) { + options.only_match = TRUE; + remove_argument (N, &argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, &argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + /* other addstar options which cannot be used in loadstarpar */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (argc < 2) { + fprintf (stderr, "USAGE: loadstarpar [options] (fitsfile) [..more files]\n"); + exit (2); + } + return (options); +} + +AddstarClientOptions args_loadstarpar_client (int argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + 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) help_client(); + + 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) help_client(); + + CPT_FILE = NULL; + if ((N = get_argument (argc, argv, "-cpt"))) { + remove_argument (N, &argc, argv); + CPT_FILE = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!CPT_FILE) help_client(); + + INPUT = NULL; + if ((N = get_argument (argc, argv, "-input"))) { + remove_argument (N, &argc, argv); + INPUT = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!INPUT) help_client(); + + /* we do not allow a subset to be extracted -- all or nothing, babe */ + UserPatch.Rmin = 0; + UserPatch.Rmax= 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + + /* only add to existing objects */ + options.only_match = FALSE; + if ((N = get_argument (argc, argv, "-only-match"))) { + options.only_match = TRUE; + remove_argument (N, &argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + options.replace = FALSE; + if ((N = get_argument (argc, argv, "-replace"))) { + options.replace = TRUE; + remove_argument (N, &argc, argv); + } + + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + /* other addstar options which cannot be used in loadstarpar */ + options.photcode = 0; + options.timeref = 0; + options.mosaic = FALSE; + options.skip_missed = FALSE; + options.closest = FALSE; + options.nosort = FALSE; + options.update = FALSE; + options.only_images = FALSE; + options.calibrate = FALSE; + options.quality_airmass = FALSE; + ACCEPT_ASTROM = FALSE; + FORCE_READ = FALSE; + TEXTMODE = FALSE; + SUBPIX = FALSE; + DUMP = NULL; + + if (argc != 1) { + fprintf (stderr, "USAGE: loadstarpar_client -cpt (file) -input (file)\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadstarpar [options] (file) [..more files]\n"); + fprintf (stderr, " add data from stellar parameter file to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -only-match : only add measurements to existing objects\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} + +static void help_client () { + + fprintf (stderr, "USAGE: loadstarpar_client -D CATDIR catdir -hostID ID -hostdir (dir) -cpt (filename) -input (input) [options]\n"); + fprintf (stderr, " add data from stellar parameter file to DVO\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -only-match : only add measurements to existing objects\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} Index: trunk/Ohana/src/addstar/src/args_loadwise.c =================================================================== --- trunk/Ohana/src/addstar/src/args_loadwise.c (revision 37729) +++ trunk/Ohana/src/addstar/src/args_loadwise.c (revision 37807) @@ -53,8 +53,19 @@ } +enum {MODE_NONE, MODE_ALLWISE, MODE_ALLSKY, MODE_PRELIM}; +int MODE; + /* load the prelim data dump */ - ALLSKY = TRUE; - if ((N = get_argument (argc, argv, "-prelim"))) { - ALLSKY = FALSE; + MODE = MODE_ALLWISE; + if ((N = get_argument (argc, argv, "-mode"))) { + MODE = MODE_NONE; + remove_argument (N, &argc, argv); + if (!strcasecmp(argv[N], "allwise")) MODE = MODE_ALLWISE; + if (!strcasecmp(argv[N], "allsky")) MODE = MODE_ALLSKY; + if (!strcasecmp(argv[N], "prelim")) MODE = MODE_PRELIM; + if (MODE == MODE_NONE) { + fprintf (stderr, "invalid mode: %s\n", argv[N]); + exit (2); + } remove_argument (N, &argc, argv); } Index: trunk/Ohana/src/addstar/src/build_links.c =================================================================== --- trunk/Ohana/src/addstar/src/build_links.c (revision 37729) +++ trunk/Ohana/src/addstar/src/build_links.c (revision 37807) @@ -1,5 +1,3 @@ # include "addstar.h" - -# define myAbortF(FORMAT,...) { fprintf (stderr, FORMAT, __VA_ARGS__); abort(); } /* Index: trunk/Ohana/src/addstar/src/coord_transforms.c =================================================================== --- trunk/Ohana/src/addstar/src/coord_transforms.c (revision 37807) +++ trunk/Ohana/src/addstar/src/coord_transforms.c (revision 37807) @@ -0,0 +1,76 @@ +# include "addstar.h" +# include "loadstarpar.h" + +// values used by libdvo/coord_systems (from Liu et al 2011, A&A 526, A16): +// Liu etal and Reid etal disagree at the 0.3 arcsec level +// transform->phi = 62.8717488056*RAD_DEG; +// transform->Xo = 32.9319185700*RAD_DEG; +// transform->xo = 282.8594812080; + +int galactic_to_celestial (double *R, double *D, double l, double b) { + + fprintf (stderr, "define me..."); + abort(); +} + +int celestial_to_galactic (double *l, double *b, double R, double D) { + +_angp = np.radians(192.859508333) # 12h 51m 26.282s (J2000) +_dngp = np.radians(27.128336111) # +27d 07' 42.01" (J2000) +_l0 = np.radians(32.932) +_ce = np.cos(_dngp) +_se = np.sin(_dngp) + + + +} + + + + +/*** python version from LSD (Mario Juric) + +# Appendix of Reid et al. (http://adsabs.harvard.edu/cgi-bin/bib_query?2004ApJ...616..872R) +# This convention is also used by LAMBDA/WMAP (http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm) +_angp = np.radians(192.859508333) # 12h 51m 26.282s (J2000) +_dngp = np.radians(27.128336111) # +27d 07' 42.01" (J2000) +_l0 = np.radians(32.932) +_ce = np.cos(_dngp) +_se = np.sin(_dngp) + +def equgal(ra, dec): + ra = np.radians(ra) + dec = np.radians(dec) + + cd, sd = np.cos(dec), np.sin(dec) + ca, sa = np.cos(ra - _angp), np.sin(ra - _angp) + + sb = cd*_ce*ca + sd*_se + l = np.arctan2(sd - sb*_se, cd*sa*_ce) + _l0 + b = np.arcsin(sb) + + l = np.where(l < 0, l + 2.*np.pi, l) + + l = np.degrees(l) + b = np.degrees(b) + + return NamedList(('l', l), ('b', b)) + +def galequ(l, b): + l = np.radians(l) + b = np.radians(b) + + cb, sb = np.cos(b), np.sin(b) + cl, sl = np.cos(l-_l0), np.sin(l-_l0) + + ra = np.arctan2(cb*cl, sb*_ce-cb*_se*sl) + _angp + dec = np.arcsin(cb*_ce*sl + sb*_se) + + ra = np.where(ra < 0, ra + 2.*np.pi, ra) + + ra = np.degrees(ra) + dec = np.degrees(dec) + + return NamedList(('ra', ra), ('dec', dec)) + + ***/ Index: trunk/Ohana/src/addstar/src/fakeimage.c =================================================================== --- trunk/Ohana/src/addstar/src/fakeimage.c (revision 37729) +++ trunk/Ohana/src/addstar/src/fakeimage.c (revision 37807) @@ -3,5 +3,5 @@ Image *fakeimage (char *rootname, off_t *Nimage, int photcode) { - int i, j, Nx, Ny, Nchips; + int i, Nx, Ny, Nchips; double pltscale, pixscale; double Rmin, Rmax, Dmin, Dmax; @@ -33,9 +33,7 @@ /* create a mosaic distortion structure */ - strcpy (MOSAIC.ctype, "DEC--DIS"); + InitCoords (&MOSAIC, "DEC--DIS"); MOSAIC.crval1 = FAKE_RA; MOSAIC.crval2 = FAKE_DEC; - - MOSAIC.crpix1 = MOSAIC.crpix2 = 0.0; /* mosaic 'pixels' are millimeters */ @@ -49,8 +47,4 @@ MOSAIC.Npolyterms = 3; - for (i = 0; i < 7; i++) { - MOSAIC.polyterms[i][0] = 0; - MOSAIC.polyterms[i][1] = 0; - } ScanConfig (config, "DPLATE_X", "%lf", 0, &pltscale); MOSAIC.polyterms[3][0] = pltscale; // L : X^3 Y^0 @@ -59,6 +53,4 @@ MOSAIC.polyterms[4][1] = pltscale; // M : X^0 Y^0 MOSAIC.polyterms[6][1] = pltscale; // M : X^2 Y^0 - - RegisterMosaic (&MOSAIC); /* some basic data about the chisp */ @@ -85,22 +77,12 @@ strcpy (image[i+1].name, name); - strcpy (image[i+1].coords.ctype, "DEC--WRP"); + InitCoords (&image[i+1].coords, "DEC--WRP"); image[i+1].coords.crval1 = dX*pixscale; image[i+1].coords.crval2 = dY*pixscale; - - image[i+1].coords.crpix1 = image[i+1].coords.crpix2 = 0.0; + image[i+1].coords.cdelt1 = image[i+1].coords.cdelt2 = pixscale; - image[i+1].coords.pc1_1 = 1.0; - image[i+1].coords.pc1_2 = 0.0; - image[i+1].coords.pc2_1 = 0.0; - image[i+1].coords.pc2_2 = 1.0; - - image[i+1].coords.Npolyterms = 0; - for (j = 0; j < 7; j++) { - image[i+1].coords.polyterms[j][0] = 0; - image[i+1].coords.polyterms[j][1] = 0; - } + image[i+1].coords.mosaic = &MOSAIC; image[i+1].sidtime = 0.0; Index: trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches.c (revision 37729) +++ trunk/Ohana/src/addstar/src/find_matches.c (revision 37807) @@ -53,4 +53,5 @@ * coordops.update.c). We use the center of the region (catalog) for crval1,2. */ + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -59,11 +60,5 @@ 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, "DEC--ARC"); /* build spatial index (RA sort) referencing input array sequence */ Index: trunk/Ohana/src/addstar/src/find_matches_ICRF.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_ICRF.c (revision 37807) +++ trunk/Ohana/src/addstar/src/find_matches_ICRF.c (revision 37807) @@ -0,0 +1,177 @@ +# include "addstar.h" +# include "loadICRF.h" + +// we find matching stars by (R,D). we add the ICRF to the array of ICRFs with appropriate links. +int find_matches_ICRF (SkyRegion *region, ICRF_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options) { + + off_t i, j, n, N, J, Jmin; + double RADIUS, RADIUS2, Rmin; + double *X1, *Y1, *X2, *Y2; + double dX, dY, dR; + off_t *N1, *N2; + off_t Nave, NAVE, Nmeas, NMEAS, Nmatch; + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, double, Nstars); + ALLOCATE (Y1, double, Nstars); + ALLOCATE (N1, off_t, Nstars); + + /** allocate local arrays (catalog) **/ + NAVE = Nave = catalog[0].Naverage; + ALLOCATE (X2, double, NAVE); + ALLOCATE (Y2, double, NAVE); + ALLOCATE (N2, off_t, NAVE); + + /* internal counters */ + Nmatch = 0; + + // current max obj ID for this catalog + unsigned int catID = catalog[0].catID; + + int Nsecfilt = catalog[0].Nsecfilt; + + NMEAS = Nmeas = catalog[0].Nmeasure; + + /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has + * the advantage that every point in R,D has a mapping to a unique X,Y. However, note + * that not all possible X,Y points map back to R,D and the local plate scale changes + * far from the projection pole. We use the center of the region (catalog) for crval1,2. + */ + Coords tcoords; + InitCoords (&tcoords, "DEC--ARC"); + tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); + tcoords.crval2 = (region[0].Dmax < 90) ? 0.5*(region[0].Dmin + region[0].Dmax) : 90.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + + /* build spatial index (RA sort) referencing input array sequence */ + for (i = 0; i < Nstars; i++) { + RD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); + N1[i] = i; + } + sort_coords_index (X1, Y1, N1, Nstars); + + /* build spatial index (RA sort) */ + for (i = 0; i < Nave; i++) { + RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + } + sort_coords_index (X2, Y2, N2, Nave); + + /* choose a radius for matches */ + RADIUS = options->radius; /* provided by config */ + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + if (!finite(X1[i]) || !finite(Y1[i])) { + i++; + continue; + } + if (!finite(X2[j]) || !finite(Y2[j])) { + j++; + continue; + } + + /* negative dX: j is too large */ + dX = X1[i] - X2[j]; + if (dX <= -1.02*RADIUS) { + i++; + continue; + } + /* positive dX, i is too large */ + if (dX >= 1.02*RADIUS) { + j++; + continue; + } + + if (stars[N1[i]].found) { + /* this star has already been assigned to an object in this or another catalog */ + i++; + continue; + } + + /* this block will match a given detection to the closest object within range of that + detection. Note that this matches ALL detections within range of the single object + to that same object. In the case of ICRF, this is what I want because I am + attaching each radio band as a separate measurement, with a single concatenated + list. This is not super efficient, but ICRF is a small list + */ + + /* within match range; look for matches */ + Jmin = -1; + Rmin = RADIUS2; + for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) { + /* find closest match for this detection */ + dX = X1[i] - X2[J]; + dY = Y1[i] - Y2[J]; + dR = dX*dX + dY*dY; + if (dR > RADIUS2) continue; + if (dR < Rmin) { + Rmin = dR; + Jmin = J; + } + } + + /* no match, try next detection */ + if (Jmin == -1) { + i++; + continue; + } + + /*** a match is found, add to measure with links to average ***/ + Nmatch ++; + n = N2[Jmin]; + N = N1[i]; + + /* make sure there is space for next entry */ + if (Nmeas >= NMEAS) { + NMEAS = Nmeas + 1000; + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + + /** add new measurement for this star **/ + + // set the new measurements + catalog[0].measure[Nmeas] = stars[N].measure; + + // measure now carries R,D (not dR,dD) + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].objID = catalog[0].average[n].objID; + catalog[0].measure[Nmeas].catID = catID; + + catalog[0].average[n].flags |= ID_ICRF_QSO; + + stars[N].found = TRUE; + catalog[0].average[n].Nmeasure ++; + Nmeas ++; + i++; + } + + // list missed QSOs, do not insert + for (i = 0; i < Nstars; i++) { + if (stars[i].found) continue; + if (!IN_REGION (stars[i].R, stars[i].D)) continue; + fprintf (stderr, "missed ICRF QSO @ %10.6f %10.6f\n", stars[i].R, stars[i].D); + } + + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + catalog[0].sorted = FALSE; + + /* check if the catalog has changed? if no change, no need to write */ + // catalog[0].objID = objID; -- I am not adding new entries, so objID does not change + catalog[0].Naverage = Nave; + catalog[0].Nmeasure = Nmeas; + catalog[0].Nsecf_mem = Nave*Nsecfilt; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nmatch); + + free (X1); + free (Y1); + free (N1); + free (N2); + free (X2); + free (Y2); + + return (Nmatch); +} Index: trunk/Ohana/src/addstar/src/find_matches_closest.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 37729) +++ trunk/Ohana/src/addstar/src/find_matches_closest.c (revision 37807) @@ -53,4 +53,5 @@ * far from the projection pole. We use the center of the region (catalog) for crval1,2. */ + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -59,11 +60,5 @@ 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, "DEC--ARC"); /* build spatial index (RA sort) referencing input array sequence */ Index: trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c (revision 37729) +++ trunk/Ohana/src/addstar/src/find_matches_closest_refstars.c (revision 37807) @@ -47,4 +47,5 @@ * far from the projection pole. We use the center of the region (catalog) for crval1,2. */ + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -53,11 +54,5 @@ tcoords.crval2 = 90.0; } - tcoords.crpix1 = 0.0; - tcoords.crpix2 = 0.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 = 0; - strcpy (tcoords.ctype, "DEC--ARC"); /* build spatial index (RA sort) referencing input array sequence */ Index: trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 37729) +++ trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 37807) @@ -51,12 +51,8 @@ /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ /* reference for coords is catalog center */ + InitCoords (&tcoords, "DEC--ZEA"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax); - tcoords.crpix1 = tcoords.crpix2 = 0.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; - strcpy (tcoords.ctype, "DEC--ZEA"); - tcoords.Npolyterms = 0; for (i = 0; i < Nstars; i++) { Index: trunk/Ohana/src/addstar/src/find_matches_starpar.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_starpar.c (revision 37807) +++ trunk/Ohana/src/addstar/src/find_matches_starpar.c (revision 37807) @@ -0,0 +1,218 @@ +# include "addstar.h" +# include "loadstarpar.h" + +// we find matching stars by (R,D). we add the starpar to the array of starpars with appropriate links. +int find_matches_starpar (SkyRegion *region, StarPar_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options) { + + off_t i, j, n, N, J, Jmin; + double RADIUS, RADIUS2, Rmin; + double *X1, *Y1, *X2, *Y2; + double dX, dY, dR; + off_t *N1, *N2; + off_t Nave, NAVE, Nstarpar, NSTARPAR, Nmatch; + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, double, Nstars); + ALLOCATE (Y1, double, Nstars); + ALLOCATE (N1, off_t, Nstars); + + /** allocate local arrays (catalog) **/ + NAVE = Nave = catalog[0].Naverage; + ALLOCATE (X2, double, NAVE); + ALLOCATE (Y2, double, NAVE); + ALLOCATE (N2, off_t, NAVE); + + /* internal counters */ + Nmatch = 0; + NSTARPAR = Nstarpar = catalog[0].Nstarpar; + + // current max obj ID for this catalog + unsigned int objID = catalog[0].objID; + unsigned int catID = catalog[0].catID; + + int Nsecfilt = catalog[0].Nsecfilt; + + /* project onto rectilinear grid with 1 arcsec pixels. the choice of ARC projection has + * the advantage that every point in R,D has a mapping to a unique X,Y. However, note + * that not all possible X,Y points map back to R,D and the local plate scale changes + * far from the projection pole. We use the center of the region (catalog) for crval1,2. + */ + Coords tcoords; + InitCoords (&tcoords, "DEC--ARC"); + tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); + tcoords.crval2 = (region[0].Dmax < 90) ? 0.5*(region[0].Dmin + region[0].Dmax) : 90.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + + /* build spatial index (RA sort) referencing input array sequence */ + for (i = 0; i < Nstars; i++) { + RD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); + N1[i] = i; + } + sort_coords_index (X1, Y1, N1, Nstars); + + /* build spatial index (RA sort) */ + for (i = 0; i < Nave; i++) { + RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + } + sort_coords_index (X2, Y2, N2, Nave); + + /* choose a radius for matches */ + RADIUS = options->radius; /* provided by config */ + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + if (!finite(X1[i]) || !finite(Y1[i])) { + i++; + continue; + } + if (!finite(X2[j]) || !finite(Y2[j])) { + j++; + continue; + } + + /* negative dX: j is too large */ + dX = X1[i] - X2[j]; + if (dX <= -1.02*RADIUS) { + i++; + continue; + } + /* positive dX, i is too large */ + if (dX >= 1.02*RADIUS) { + j++; + continue; + } + + if (stars[N1[i]].found) { + /* this star has already been assigned to an object in this or another catalog */ + i++; + continue; + } + + /* this block will match a given detection to the closest object within range of that detection. + XXX note that this matches ALL detections within range of the single object to that same object + this is bad, but I cannot just go in linear order (ie, mark off each object as they are + used). I should make a list of all Nobj * Ndet pairs in range and choose the matches + based on their separations. UGH + */ + + /* within match range; look for matches */ + Jmin = -1; + Rmin = RADIUS2; + for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) { + /* find closest match for this detection */ + dX = X1[i] - X2[J]; + dY = Y1[i] - Y2[J]; + dR = dX*dX + dY*dY; + if (dR > RADIUS2) continue; + if (dR < Rmin) { + Rmin = dR; + Jmin = J; + } + } + + /* no match, try next detection */ + if (Jmin == -1) { + i++; + continue; + } + + /*** a match is found, add to starpar with links to average ***/ + Nmatch ++; + n = N2[Jmin]; + N = N1[i]; + + /* make sure there is space for next entry */ + if (Nstarpar >= NSTARPAR) { + NSTARPAR = Nstarpar + 1000; + REALLOCATE (catalog[0].starpar, StarPar, NSTARPAR); + } + + /** add starpar for this star **/ + + // set the new starparments + catalog[0].starpar[Nstarpar] = stars[N].starpar; + catalog[0].starpar[Nstarpar].averef = n; + catalog[0].starpar[Nstarpar].objID = catalog[0].average[n].objID; + catalog[0].starpar[Nstarpar].catID = catID; + + // NOTE: include R,D in starpar? + // XXX if I add R,D to starpar, I should rationalize R,D to the same boundary? + + stars[N].found = TRUE; + catalog[0].average[n].Nstarpar ++; + Nstarpar ++; + i++; + } + +# if (1) + /** incorporate unmatched image stars? **/ + for (i = 0; (i < Nstars) && !options->only_match; i++) { + + // skip already matched stars + if (stars[i].found) continue; + if (!IN_REGION (stars[i].R, stars[i].D)) continue; + + /* make sure there is space for next entry */ + if (Nstarpar >= NSTARPAR) { + NSTARPAR = Nstarpar + 1000; + REALLOCATE (catalog[0].starpar, StarPar, NSTARPAR); + } + if (Nave >= NAVE) { + NAVE = Nave + 1000; + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*Nsecfilt); + } + + dvo_average_init (&catalog[0].average[Nave]); + catalog[0].average[Nave].R = stars[i].R; + catalog[0].average[Nave].D = stars[i].D; + + catalog[0].average[Nave].Nstarpar = 1; + catalog[0].average[Nave].starparOffset = Nstarpar; + catalog[0].average[Nave].objID = objID; + catalog[0].average[Nave].catID = catID; + + if (PSPS_ID) { + catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R, catalog[0].average[Nave].D); + } + + objID ++; + + for (j = 0; j < Nsecfilt; j++) { + dvo_secfilt_init (&catalog[0].secfilt[Nave*Nsecfilt+j]); + } + + catalog[0].starpar[Nstarpar] = stars[i].starpar; + catalog[0].starpar[Nstarpar].averef = Nave; + catalog[0].starpar[Nstarpar].objID = catalog[0].average[Nave].objID; + catalog[0].starpar[Nstarpar].catID = catID; + + stars[i].found = TRUE; + Nstarpar ++; + Nave ++; + } +# endif + + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].starpar, StarPar, Nstarpar); + + catalog[0].sorted = FALSE; + + /* check if the catalog has changed? if no change, no need to write */ + catalog[0].objID = objID; // new max value, save on catalog close + catalog[0].Naverage = Nave; + catalog[0].Nstarpar = Nstarpar; + catalog[0].Nsecf_mem = Nave*Nsecfilt; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nstarpar: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n", Nstars, Nave, Nstarpar, Nmatch); + + free (X1); + free (Y1); + free (N1); + free (N2); + free (X2); + free (Y2); + + return (Nmatch); +} Index: trunk/Ohana/src/addstar/src/find_missing.c =================================================================== --- trunk/Ohana/src/addstar/src/find_missing.c (revision 37729) +++ trunk/Ohana/src/addstar/src/find_missing.c (revision 37807) @@ -14,5 +14,4 @@ REALLOCATE (catalog[0].missing, Missing, NMISS); } - if (!FindMosaicForImage (overlap, Noverlap, j)) continue; if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue; add_miss_link (&catalog[0].average[Nave], next_miss, Nmiss); @@ -27,5 +26,5 @@ /* add reference for undetected catalog stars */ /* XXX allow this option only for single images? */ -if (!strcmp (&image[0].coords.ctype[4], "-WRP")) RegisterMosaic (mosaic); +if (!strcmp (&image[0].coords.ctype[4], "-WRP")) image[0].coords = mosaic; for (j = 0; (j < Nave) && !options.skip_missed; j++) { n = N2[j]; Index: trunk/Ohana/src/addstar/src/loadICRF.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF.c (revision 37807) @@ -0,0 +1,54 @@ +# include "addstar.h" +# include "loadICRF.h" + +/* This is the DVO program to upload ICRF QSO data (eg, from Leonid Petrov) into a DVO + database. It is modeled on the loadstarpar program and is expected to be run only + rarely (once?). The ICRF QSO parameter data are delivered as text file. It does not + allow a subset of the sky to be uploaded; entire QSO parameter files are loaded if + supplied on the command line. + + USAGE: loadICRF -D CATDIR (catdir) (file.fits) [...more files] +*/ + +int main (int argc, char **argv) { + + int i; + AddstarClientOptions options; + + options = ConfigInit (&argc, argv); + options = args_loadICRF (&argc, argv, options); // XXX fix args after rest is done + + // load the full sky description table (dvodb must exist) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // load the list of hosts + HostTable *hosts = NULL; + if (PARALLEL) { + hosts = HostTableLoad (CATDIR, sky->hosts); + if (!hosts) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + // ensure that the paths are absolute path names + for (i = 0; i < hosts->Nhosts; i++) { + char *tmppath = abspath (hosts->hosts[i].pathname, DVO_MAX_PATH); + free (hosts->hosts[i].pathname); + hosts->hosts[i].pathname = tmppath; + } + + // set up the array of active hosts + init_remote_hosts (); + } + + // generate the subset matching the user-selected region + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + // argv[i] are a list of input tables; load one at a time + for (i = 1; i < argc; i++) { + fprintf (stderr, "loading %s\n", argv[i]); + loadICRF_table (skylist, hosts, argv[i], &options); + } + exit (0); +} Index: trunk/Ohana/src/addstar/src/loadICRF_catalog.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_catalog.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_catalog.c (revision 37807) @@ -0,0 +1,28 @@ +# include "addstar.h" +# include "loadICRF.h" + +int loadICRF_catalog (ICRF_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options) { + + Catalog catalog; + + // now we have all of the loaded stars in this catalog + catalog.filename = filename; + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_MEAS; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, region, VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + find_matches_ICRF (region, stars, Nstars, &catalog, options); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + return (TRUE); +} Index: trunk/Ohana/src/addstar/src/loadICRF_client.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_client.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_client.c (revision 37807) @@ -0,0 +1,30 @@ +# include "addstar.h" +# include "loadICRF.h" + +int main (int argc, char **argv) { + + AddstarClientOptions options; + + // need to construct these options with args_loadICRF... + options = ConfigInit (&argc, argv); + options = args_loadICRF_client (&argc, argv, options); + + // client is called with a pointer to the file to be loaded + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + SkyList *skylist = SkyRegionByCPT (sky, CPT_FILE); + + char cptfile[1024]; + snprintf (cptfile, 1024, "%s/%s.cpt", HOSTDIR, skylist[0].regions[0]->name); + + int Nstars; + ICRF_Stars *stars = loadICRF_load_stars (INPUT, &Nstars); + + loadICRF_catalog (stars, Nstars, skylist->regions[0], cptfile, &options); + + free (stars); + + exit (0); +} Index: trunk/Ohana/src/addstar/src/loadICRF_io.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_io.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_io.c (revision 37807) @@ -0,0 +1,189 @@ +# include "addstar.h" +# include "loadICRF.h" + +int loadICRF_save_stars (char *filename, ICRF_Stars *stars, int Nstars) { + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + myAbortF ("ERROR: cannot open file for output %s\n", filename); + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + gfits_create_table_header (&theader, "BINTABLE", "ICRF"); + + gfits_define_bintable_column (&theader, "D", "RA", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "RA_ERR", "", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_ERR", "", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "F_PSF", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dF_PSF", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "F_AP", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dF_AP", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "M_PSF", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dM_PSF", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "M_AP", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dM_AP", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "PHOTCODE", "", "", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "FLAGS", "", "", 1.0, FT_BZERO_INT32); + + int i; + double *ra = NULL; ALLOCATE (ra , double , Nstars); for (i = 0; i < Nstars; i++) ra [i] = stars[i].R ; + double *dec = NULL; ALLOCATE (dec , double , Nstars); for (i = 0; i < Nstars; i++) dec [i] = stars[i].D ; + float *dR = NULL; ALLOCATE (dR , float , Nstars); for (i = 0; i < Nstars; i++) dR [i] = stars[i].measure.dXccd ; + float *dD = NULL; ALLOCATE (dD , float , Nstars); for (i = 0; i < Nstars; i++) dD [i] = stars[i].measure.dYccd ; + float *Fpsf = NULL; ALLOCATE (Fpsf , float , Nstars); for (i = 0; i < Nstars; i++) Fpsf [i] = stars[i].measure.FluxPSF ; + float *dFpsf = NULL; ALLOCATE (dFpsf , float , Nstars); for (i = 0; i < Nstars; i++) dFpsf [i] = stars[i].measure.dFluxPSF; + float *Fap = NULL; ALLOCATE (Fap , float , Nstars); for (i = 0; i < Nstars; i++) Fap [i] = stars[i].measure.FluxAp ; + float *dFap = NULL; ALLOCATE (dFap , float , Nstars); for (i = 0; i < Nstars; i++) dFap [i] = stars[i].measure.dFluxAp ; + float *Mpsf = NULL; ALLOCATE (Mpsf , float , Nstars); for (i = 0; i < Nstars; i++) Mpsf [i] = stars[i].measure.M ; + float *dMpsf = NULL; ALLOCATE (dMpsf , float , Nstars); for (i = 0; i < Nstars; i++) dMpsf [i] = stars[i].measure.dM ; + float *Map = NULL; ALLOCATE (Map , float , Nstars); for (i = 0; i < Nstars; i++) Map [i] = stars[i].measure.Map ; + float *dMap = NULL; ALLOCATE (dMap , float , Nstars); for (i = 0; i < Nstars; i++) dMap [i] = stars[i].measure.dMap ; + unsigned int *photcode = NULL; ALLOCATE (photcode, unsigned int, Nstars); for (i = 0; i < Nstars; i++) photcode[i] = stars[i].measure.photcode; + unsigned int *flags = NULL; ALLOCATE (flags , unsigned int, Nstars); for (i = 0; i < Nstars; i++) flags [i] = stars[i].measure.dbFlags ; + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", ra , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DEC", dec , Nstars); + gfits_set_bintable_column (&theader, &ftable, "RA_ERR", dR , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DEC_ERR", dD , Nstars); + gfits_set_bintable_column (&theader, &ftable, "F_PSF", Fpsf , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dF_PSF", dFpsf , Nstars); + gfits_set_bintable_column (&theader, &ftable, "F_AP", Fap , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dF_AP", dFap , Nstars); + gfits_set_bintable_column (&theader, &ftable, "M_PSF", Mpsf , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dM_PSF", dMpsf , Nstars); + gfits_set_bintable_column (&theader, &ftable, "M_AP", Map , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dM_AP", dMap , Nstars); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, Nstars); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nstars); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + + return TRUE; +} + + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +ICRF_Stars *loadICRF_load_stars (char *filename, int *nstars) { + + int i, Ncol; + off_t Nrow; + char type[16]; + + ICRF_Stars *stars = NULL; + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + + 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 header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read matrix\n"); + goto escape; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + GET_COLUMN(ra , "RA", double); + GET_COLUMN(dec , "DEC", double); + GET_COLUMN(dR , "RA_ERR", float); + GET_COLUMN(dD , "DEC_ERR", float); + GET_COLUMN(Fpsf , "F_PSF", float); + GET_COLUMN(dFpsf , "dF_PSF", float); + GET_COLUMN(Fap , "F_AP", float); + GET_COLUMN(dFap , "dF_AP", float); + GET_COLUMN(Mpsf , "M_PSF", float); + GET_COLUMN(dMpsf , "dM_PSF", float); + GET_COLUMN(Map , "M_AP", float); + GET_COLUMN(dMap , "dM_AP", float); + GET_COLUMN(photcode, "PHOTCODE", int); + GET_COLUMN(flags , "FLAGS", int); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + ALLOCATE (stars, ICRF_Stars, Nrow); + + for (i = 0; i < Nrow; i++) { + InitICRF_Star (&stars[i]); + stars[i].R = ra [i]; + stars[i].D = dec [i]; + stars[i].measure.R = ra [i]; + stars[i].measure.D = dec [i]; + stars[i].measure.dXccd = dR [i]; + stars[i].measure.dYccd = dD [i]; + stars[i].measure.FluxPSF = Fpsf [i]; + stars[i].measure.dFluxPSF= dFpsf [i]; + stars[i].measure.FluxAp = Fap [i]; + stars[i].measure.dFluxAp = dFap [i]; + stars[i].measure.M = Mpsf [i]; + stars[i].measure.dM = dMpsf [i]; + stars[i].measure.Map = Map [i]; + stars[i].measure.dMap = dMap [i]; + stars[i].measure.photcode= photcode[i]; + stars[i].measure.dbFlags = flags [i]; + } + + fclose (f); + *nstars = Nrow; + return stars; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + return NULL; +} Index: trunk/Ohana/src/addstar/src/loadICRF_make_subset.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_make_subset.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_make_subset.c (revision 37807) @@ -0,0 +1,44 @@ +# include "addstar.h" +# include "loadICRF.h" + +// assign stars in the given region to the subset (NOTE: stars are sorted by RA, start is +// first entry in stars array in this region) + +ICRF_Stars *loadICRF_make_subset (ICRF_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) { + + int i; + + ICRF_Stars *subset = NULL; + + // collect array of (Stars *) stars in a new output catalog + int Nsubset = 0; + int NSUBSET = 3000; + ALLOCATE (subset, ICRF_Stars, NSUBSET); + + // find the rest of the stars in this output region + for (i = start; i < Nstars; i++) { + if (stars[i].flag) continue; + + // check if in skyregion + if (stars[i].R < region[0].Rmin) continue; + if (stars[i].R > region[0].Rmax) break; + if (stars[i].D < region[0].Dmin) continue; + if (stars[i].D > region[0].Dmax) continue; + + // check if in UserPatch (a GLOBAL) + if (stars[i].R < UserPatch.Rmin) continue; + if (stars[i].R > UserPatch.Rmax) break; + if (stars[i].D < UserPatch.Dmin) continue; + if (stars[i].D > UserPatch.Dmax) continue; + + subset[Nsubset] = stars[i]; + Nsubset ++; + + stars[i].flag = TRUE; // set flag for stars we have already assigned to a subset + + CHECK_REALLOCATE (subset, ICRF_Stars, NSUBSET, Nsubset, 10000); + } + + *nsubset = Nsubset; + return subset; +} Index: trunk/Ohana/src/addstar/src/loadICRF_readstars.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_readstars.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_readstars.c (revision 37807) @@ -0,0 +1,180 @@ +# include "addstar.h" +# include "loadICRF.h" + +// ICRF table is ASCII text, with the following interesting columns (fixed bytes) +// 0 : type (C - calibrator, N - non-calibrator, U - unreliable) +// 3 - 21 : names (make into an extID?) +// 24 - 39 : RA (ICRF) +// 40 - 55 : DEC (ICRF) +// 57 - 62 : dR (mas) [use dX,dY in milliarcsec 'pixels'] +// 64 - 69 : dD (mas) +// 79 - 84 : Number of observations used (put in dt?) +// 87 - 92 : S-band total flux density integrated over entire map, Jy (note 87 may be <) : +// 94 - 99 : S-band unresolved flux density at long VLBA baselines, Jy +// 102 - 107 : C-band total flux density integrated over entire map, Jy +// 109 - 114 : C-band unresolved flux density at long VLBA baselines, Jy +// 117 - 122 : X-band total flux density integrated over entire map, Jy +// 124 - 129 : X-band unresolved flux density at long VLBA baselines, Jy +// 132 - 137 : U-band total flux density integrated over entire map, Jy +// 139 - 144 : U-band unresolved flux density at long VLBA baselines, Jy +// 147 - 152 : K-band total flux density integrated over entire map, Jy +// 154 - 159 : K-band unresolved flux density at long VLBA baselines, Jy + +// I'm going to save these are FluxPSF (total), FluxAp (unresolved), and mags for +// the signficant ones (AB_m from Jy). Use new photcodes to represent C,X,U,K + +/* +# IVS name J2000 name Right ascension Declination D_alp D_Del Corr #Obs S-band flux C-band Flux X-band Flux U-band flux K-band Flux Type Cat +# Total Unres Total Unres Total Unres Total Unres Total Unres +# hr mn seconds deg mn seconds mas mas Jy Jy Jy Jy Jy Jy Jy Jy Jy Jy +C 2357+080 J0000+0816 00 00 07.031141 +08 16 45.05175 0.46 0.85 0.758 41 -1.00 -1.00 -1.00 -1.00 0.020 <0.014 -1.00 -1.00 -1.00 -1.00 X rfc_2014c +01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789 +0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 +*/ + + +int loadICRF_addstar(ICRF_Stars *star, float fluxPSF, float fluxAp, double R, double D, float dX, float dY, int photcode); + +ICRF_Stars *loadICRF_readstars (char *filename, int *nstars) { + + // read in the full FITS files ('cause I don't have a partial read option) + FILE *f = fopen (filename, "r"); + if (f == NULL) Shutdown ("can't read stellar parameter file: %s", filename); + + int Nstars = 0; + int NSTARS = 1000; + ICRF_Stars *stars = NULL; + ALLOCATE (stars, ICRF_Stars, NSTARS); + + double R, D; + + // ICRF QSO photcodes (we store flux in Jy and AB mags) + PhotCode *codeS = GetPhotcodebyName ("ICRF_S"); if (!codeS) { fprintf (stderr, "missing ICRF_S, exiting\n"); exit (1); } + PhotCode *codeC = GetPhotcodebyName ("ICRF_C"); if (!codeC) { fprintf (stderr, "missing ICRF_C, exiting\n"); exit (1); } + PhotCode *codeX = GetPhotcodebyName ("ICRF_X"); if (!codeX) { fprintf (stderr, "missing ICRF_X, exiting\n"); exit (1); } + PhotCode *codeU = GetPhotcodebyName ("ICRF_U"); if (!codeU) { fprintf (stderr, "missing ICRF_U, exiting\n"); exit (1); } + PhotCode *codeK = GetPhotcodebyName ("ICRF_K"); if (!codeK) { fprintf (stderr, "missing ICRF_K, exiting\n"); exit (1); } + + /* read in stars line-by-line */ + char line[1024]; + while (scan_line (f, line) != EOF) { + stripwhite (line); + if (line[0] == 0) continue; + if (line[0] == '#') continue; + + // objects with values other than "C" in the first column are not good for + // calibration, so just skip. + if (line[0] != 'C') continue; + + // terminate the partial strings containing RA and DEC: + char cr = line[39]; line[39] = 0; + char cd = line[55]; line[55] = 0; + + if (!ohana_str_to_radec (&R, &D, &line[24], &line[40])) { + line[39] = cr; + line[55] = cd; + fprintf (stderr, "problem with coords: %s\n", line); + exit (1); + } + R = ohana_normalize_angle (R); + + line[63] = 0; + float dR = atof(&line[57]); + line[70] = 0; + float dD = atof(&line[64]); + + int dX = 1000.0*100*dR; + int dY = 1000.0*100*dD; + + // XXX set this somewhere? + // line[85] = 0; + // int Nmeas = atoi(&line[79]); + + line[ 93] = 0; float fluxSt = (line[ 87] == '<') ? NAN : atof(&line[ 87]); if (fluxSt == -1.0) fluxSt = NAN; + line[100] = 0; float fluxSu = (line[ 94] == '<') ? NAN : atof(&line[ 94]); if (fluxSu == -1.0) fluxSu = NAN; + line[108] = 0; float fluxCt = (line[102] == '<') ? NAN : atof(&line[102]); if (fluxCt == -1.0) fluxCt = NAN; + line[115] = 0; float fluxCu = (line[109] == '<') ? NAN : atof(&line[109]); if (fluxCu == -1.0) fluxCu = NAN; + line[123] = 0; float fluxXt = (line[117] == '<') ? NAN : atof(&line[117]); if (fluxXt == -1.0) fluxXt = NAN; + line[130] = 0; float fluxXu = (line[124] == '<') ? NAN : atof(&line[124]); if (fluxXu == -1.0) fluxXu = NAN; + line[138] = 0; float fluxUt = (line[132] == '<') ? NAN : atof(&line[132]); if (fluxUt == -1.0) fluxUt = NAN; + line[145] = 0; float fluxUu = (line[139] == '<') ? NAN : atof(&line[139]); if (fluxUu == -1.0) fluxUu = NAN; + line[153] = 0; float fluxKt = (line[147] == '<') ? NAN : atof(&line[147]); if (fluxKt == -1.0) fluxKt = NAN; + line[160] = 0; float fluxKu = (line[154] == '<') ? NAN : atof(&line[154]); if (fluxKu == -1.0) fluxKu = NAN; + + // we treat each valid flux as a new star + if (isfinite(fluxSt)) { + loadICRF_addstar (&stars[Nstars], fluxSt, fluxSu, R, D, dX, dY, codeS->code); + Nstars++; + CHECK_REALLOCATE (stars, ICRF_Stars, NSTARS, Nstars, 1000); + } + if (isfinite(fluxCt)) { + loadICRF_addstar (&stars[Nstars], fluxCt, fluxCu, R, D, dX, dY, codeC->code); + Nstars++; + CHECK_REALLOCATE (stars, ICRF_Stars, NSTARS, Nstars, 1000); + } + if (isfinite(fluxXt)) { + loadICRF_addstar (&stars[Nstars], fluxXt, fluxXu, R, D, dX, dY, codeX->code); + Nstars++; + CHECK_REALLOCATE (stars, ICRF_Stars, NSTARS, Nstars, 1000); + } + if (isfinite(fluxUt)) { + loadICRF_addstar (&stars[Nstars], fluxUt, fluxUu, R, D, dX, dY, codeU->code); + Nstars++; + CHECK_REALLOCATE (stars, ICRF_Stars, NSTARS, Nstars, 1000); + } + if (isfinite(fluxKt)) { + loadICRF_addstar (&stars[Nstars], fluxKt, fluxKu, R, D, dX, dY, codeK->code); + Nstars++; + CHECK_REALLOCATE (stars, ICRF_Stars, NSTARS, Nstars, 1000); + } + + } + *nstars = Nstars; + return (stars); +} + +int loadICRF_addstar(ICRF_Stars *star, float fluxPSF, float fluxAp, double R, double D, float dX, float dY, int photcode) { + + InitICRF_Star (star); + + star->R = R; + star->D = D; + star->measure.R = R; + star->measure.D = D; + star->measure.dXccd = dX; // in centi-pixels, where 1 pixel = 1 mas + star->measure.dYccd = dY; // in centi-pixels, where 1 pixel = 1 mas + + star->measure.FluxPSF = fluxPSF; // flux in Jy + star->measure.dFluxPSF = 0.01; // this is just a guess + star->measure.FluxAp = fluxAp; // flux in Jy + star->measure.dFluxAp = 0.01; // this is just a guess + + star->measure.M = (fluxPSF > 0.0) ? 8.9 - 2.5*log10(fluxPSF) : NAN; + star->measure.dM = (fluxPSF > 0.0) ? 0.01 / fluxPSF : NAN; + star->measure.Map = (fluxAp > 0.0) ? 8.9 - 2.5*log10(fluxAp) : NAN; + star->measure.dMap = (fluxAp > 0.0) ? 0.01 / fluxAp: NAN; + star->measure.photcode = photcode; + star->measure.dbFlags |= ID_MEAS_ICRF_QSO; + return TRUE; +} + +int loadICRF_sortStars (ICRF_Stars *stars, int Nstars) { + +# define SWAPFUNC(A,B){ ICRF_Stars temp = stars[A]; stars[A] = stars[B]; stars[B] = temp; } +# define COMPARE(A,B)(stars[A].R < stars[B].R) + + OHANA_SORT (Nstars, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + + return TRUE; +} + +int InitICRF_Star (ICRF_Stars *star) { + + dvo_measure_init (&star[0].measure); + star[0].found = FALSE; // FALSE = not yet assigned to an object + star[0].flag = FALSE; // FALSE = not yet assigned to a subset + return TRUE; +} Index: trunk/Ohana/src/addstar/src/loadICRF_remote_hosts.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_remote_hosts.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_remote_hosts.c (revision 37807) @@ -0,0 +1,144 @@ +# include "addstar.h" +# include "loadICRF.h" +# include +# include + +# define DEBUG 1 + +// we are running N parallel remote jobs on the remote hosts. In the usual remote +// processing (eg, relphot or relastro), we launch one job per remote host. In this case, +// we are launching one job foreach table (catalog.cpt) being touched. individual hosts +// may have more than one job active at a time. + +int Nhosts = 0; +HostInfo **hosts = NULL; + +int init_remote_hosts () { + + Nhosts = 10; + ALLOCATE (hosts, HostInfo *, Nhosts); + + int i; + for (i = 0; i < Nhosts; i++) { + hosts[i] = NULL; + } + return TRUE; +} + +int find_empty_slot () { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) return i; + } + return -1; +} + +int harvest_all () { + + int slot = -1; + + while (slot != -2) { + slot = harvest_host (); + usleep (50000); + } + return TRUE; +} + +// put the host in the list in a free slot. NOTE: this should never be called if we do not +// already have a free slot. +int save_remote_host (HostInfo *host) { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) { + hosts[i] = host; + return TRUE; + } + } + myAbort ("failed to find an empty slot: this is a programming error"); + return FALSE; +} + +// wait for all children to complete, report output to stdout +// possible states when function is finished: +// * no child was harvested, but children still exist +// * no outstanding children +// * a child was harvested +int harvest_host () { + + // check if any children have finished... + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + + if (!pid) return -1; // there are outstanding children, but none have exited, no slot was opened + if ((pid == -1) && (errno == ECHILD)) return -2; // there are no outstanding children + if ((pid == -1) && (errno != ECHILD)) myAbort ("programming error?"); + + if (DEBUG) fprintf (stdout, "PID %d done\n", pid); + + // find the host which has finished + int i; + int found = FALSE; + for (i = 0; (i < Nhosts) && !found; i++) { + if (!hosts[i]) continue; // unassigned slot + if (hosts[i][0].pid != pid) continue; + found = TRUE; + break; + } + myAssert (found, "Programming error: failed to matched finished job to known host!"); + int slot = i; + + HostInfo *host = hosts[slot]; + + // check on the status of this job and report any output + if (DEBUG) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + int Nout, printHead; + + // read stdout + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (DEBUG) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + } + + // close opened connections + close (host->stdio[HOST_STDIN]); + close (host->stdio[HOST_STDOUT]); + close (host->stdio[HOST_STDERR]); + + // free data associated with the host + free (host->hostname); + free (host->pathname); + FreeIOBuffer (&host->stdout); + FreeIOBuffer (&host->stderr); + free (host); + + // free the slot + hosts[slot] = NULL; + return slot; +} + Index: trunk/Ohana/src/addstar/src/loadICRF_save_remote.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_save_remote.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_save_remote.c (revision 37807) @@ -0,0 +1,90 @@ +# include "addstar.h" +# include "loadICRF.h" +# define DEBUG 1 + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + +int loadICRF_save_remote (ICRF_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // if this region is a parallel thing, save and launch remote + if (!PARALLEL) { + loadICRF_catalog (stars, Nstars, region, fullname, options); + } else { + int N = hosts->index[region->hostID]; + HostInfo *hostMach = &hosts->hosts[N]; + + // save to a unique filename + char filename[1024]; // CATDIR/tmpdir/ICRF.PID.index.fits + snprintf (filename, 1024, "%s/tmpdir/ICRF.%s.%05d.fits", CATDIR, uniquer, region->index); + if (!check_file_access (filename, FALSE, TRUE, TRUE)) exit (2); + + // write the data to the given FITS file + loadICRF_save_stars (filename, stars, Nstars); + + int slot = -1; + while (slot == -1) { + slot = find_empty_slot (); + if (slot == -1) { + usleep (50000); + slot = harvest_host(); + myAssert (slot != -2, "we should not call harvest_host here if we have open slots"); + } + } + + // allocate a host for this job + HostInfo *host = NULL; + ALLOCATE (host, HostInfo, 1); + InitHost (host); + + // we want to run this job on the host described by hostMach. copy + // immutable data from hostMach to a locally allocated host: + host->hostID = hostMach->hostID; + host->hostname = strcreate (hostMach->hostname); + host->pathname = strcreate (hostMach->pathname); + InitIOBuffer (&host->stdout, 1000); + InitIOBuffer (&host->stderr, 1000); + + // got a valid slot, so launch a new host + + // need to generate the remote command + char command[1024]; + snprintf (command, 1024, "loadICRF_client"); + strextend (command, "-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-hostdir %s", host->pathname); + strextend (command, "-cpt %s", region->name); + strextend (command, "-input %s", filename); + strextend (command, "-D ADDSTAR_RADIUS %f", options->radius); + + fprintf (stderr, "command: %s\n", command); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + fprintf (stderr, "launch cpt %s on %s @ %d\n", region->name, host->pathname, pid); + + save_remote_host (host); + } + return TRUE; +} Index: trunk/Ohana/src/addstar/src/loadICRF_table.c =================================================================== --- trunk/Ohana/src/addstar/src/loadICRF_table.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadICRF_table.c (revision 37807) @@ -0,0 +1,50 @@ +# include "addstar.h" +# include "loadICRF.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +int loadICRF_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options) { + + int i, Nstars; + + ICRF_Stars *stars = loadICRF_readstars (filename, &Nstars); + + // sort the stars by RA + loadICRF_sortStars (stars, Nstars); + + // scan through the stars, loading the containing catalogs + // skip through table for unsaved stars + for (i = 0; i < Nstars; i++) { + if (stars[i].flag) continue; + + // scan forward until we read the UserPatch + if (stars[i].R < UserPatch.Rmin) continue; + if (stars[i].R > UserPatch.Rmax) break; + if (stars[i].D < UserPatch.Dmin) continue; + if (stars[i].D > UserPatch.Dmax) continue; + + // identify the relevant catalog + SkyList *skylist = SkyRegionByPoint_List (skylistInput, -1, stars[i].R, stars[i].D); + if (skylist[0].Nregions == 0) { + SkyListFree (skylist); + continue; + } + SkyRegion *region = skylist[0].regions[0]; + + // select stars matching this region + int Nsubset; + ICRF_Stars *subset = loadICRF_make_subset (stars, Nstars, i, region, &Nsubset); + + // In parallel mode, write out the subset to a disk file. Block until a remote host + // is available. In serial mode, just match against the appropriate region and save + loadICRF_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options); + } + + // wait for last remote clients to finish + harvest_all (); + + return TRUE; +} + Index: trunk/Ohana/src/addstar/src/loadstarpar.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar.c (revision 37807) @@ -0,0 +1,54 @@ +# include "addstar.h" +# include "loadstarpar.h" + +/* This is the DVO program to upload Greg Green's stellar parametersq into a DVO database. + It is modeled on the loadwise program and is expected to be run only rarely (once?). + The stellar parameter data are delivered as *.fits files. It does not allow a subset + of the sky to be uploaded; entire stellar parameter files are loaded if supplied on the + command line. + + USAGE: loadstarpar -D CATDIR (catdir) (file.fits) [...more files] +*/ + +int main (int argc, char **argv) { + + int i; + AddstarClientOptions options; + + // need to construct these options with args_loadstarpar... + options = ConfigInit (&argc, argv); + options = args_loadstarpar (argc, argv, options); + + // load the full sky description table (dvodb must exist) + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // load the list of hosts + HostTable *hosts = NULL; + if (PARALLEL) { + hosts = HostTableLoad (CATDIR, sky->hosts); + if (!hosts) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + // ensure that the paths are absolute path names + for (i = 0; i < hosts->Nhosts; i++) { + char *tmppath = abspath (hosts->hosts[i].pathname, DVO_MAX_PATH); + free (hosts->hosts[i].pathname); + hosts->hosts[i].pathname = tmppath; + } + + // set up the array of active hosts + init_remote_hosts (); + } + + // generate the subset matching the user-selected region + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + for (i = 1; i < argc; i++) { + fprintf (stderr, "loading %s\n", argv[i]); + loadstarpar_table (skylist, hosts, argv[i], &options); + } + exit (0); +} Index: trunk/Ohana/src/addstar/src/loadstarpar_catalog.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_catalog.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_catalog.c (revision 37807) @@ -0,0 +1,28 @@ +# include "addstar.h" +# include "loadstarpar.h" + +int loadstarpar_catalog (StarPar_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options) { + + Catalog catalog; + + // now we have all of the loaded stars in this catalog + catalog.filename = filename; + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_SECF | LOAD_STARPAR; + catalog.Nsecfilt = GetPhotcodeNsecfilt (); + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, region, VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open/create catalog file %s\n", catalog.filename); + exit (2); + } + + find_matches_starpar (region, stars, Nstars, &catalog, options); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + return (TRUE); +} Index: trunk/Ohana/src/addstar/src/loadstarpar_client.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_client.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_client.c (revision 37807) @@ -0,0 +1,27 @@ +# include "addstar.h" +# include "loadstarpar.h" + +int main (int argc, char **argv) { + + AddstarClientOptions options; + + // need to construct these options with args_loadstarpar... + options = ConfigInit (&argc, argv); + options = args_loadstarpar_client (argc, argv, options); + + // client is called with a pointer to the file to be loaded + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + SkyList *skylist = SkyRegionByCPT (sky, CPT_FILE); + + int Nstars; + StarPar_Stars *stars = loadstarpar_load_stars (INPUT, &Nstars); + + loadstarpar_catalog (stars, Nstars, skylist->regions[0], CPT_FILE, &options); + + free (stars); + + exit (0); +} Index: trunk/Ohana/src/addstar/src/loadstarpar_io.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_io.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_io.c (revision 37807) @@ -0,0 +1,191 @@ +# include "addstar.h" +# include "loadstarpar.h" + +int loadstarpar_save_stars (char *filename, StarPar_Stars *stars, int Nstars) { + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + myAbortF ("ERROR: cannot open file for output %s\n", filename); + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + gfits_create_table_header (&theader, "BINTABLE", "STARPAR"); + + gfits_define_bintable_column (&theader, "D", "RA", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "GLON", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "GLAT", "", "degree", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Ebv", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dEbv", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DistMag", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dDistMag", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "M_r", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dM_r", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "FeH", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "dFeH", "", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "averef", "", "", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "objID", "", "", 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "catID", "", "", 1.0, FT_BZERO_INT32); + + int i; + double *ra = NULL; ALLOCATE (ra , double , Nstars); for (i = 0; i < Nstars; i++) ra [i] = stars[i].starpar.R ; + double *dec = NULL; ALLOCATE (dec , double , Nstars); for (i = 0; i < Nstars; i++) dec [i] = stars[i].starpar.D ; + float *glon = NULL; ALLOCATE (glon , float , Nstars); for (i = 0; i < Nstars; i++) glon [i] = stars[i].starpar.galLon ; + float *glat = NULL; ALLOCATE (glat , float , Nstars); for (i = 0; i < Nstars; i++) glat [i] = stars[i].starpar.galLat ; + float *Ebv = NULL; ALLOCATE (Ebv , float , Nstars); for (i = 0; i < Nstars; i++) Ebv [i] = stars[i].starpar.Ebv ; + float *dEbv = NULL; ALLOCATE (dEbv , float , Nstars); for (i = 0; i < Nstars; i++) dEbv [i] = stars[i].starpar.dEbv ; + float *DistMag = NULL; ALLOCATE (DistMag , float , Nstars); for (i = 0; i < Nstars; i++) DistMag [i] = stars[i].starpar.DistMag ; + float *dDistMag = NULL; ALLOCATE (dDistMag, float , Nstars); for (i = 0; i < Nstars; i++) dDistMag[i] = stars[i].starpar.dDistMag; + float *M_r = NULL; ALLOCATE (M_r , float , Nstars); for (i = 0; i < Nstars; i++) M_r [i] = stars[i].starpar.M_r ; + float *dM_r = NULL; ALLOCATE (dM_r , float , Nstars); for (i = 0; i < Nstars; i++) dM_r [i] = stars[i].starpar.dM_r ; + float *FeH = NULL; ALLOCATE (FeH , float , Nstars); for (i = 0; i < Nstars; i++) FeH [i] = stars[i].starpar.FeH ; + float *dFeH = NULL; ALLOCATE (dFeH , float , Nstars); for (i = 0; i < Nstars; i++) dFeH [i] = stars[i].starpar.dFeH ; + unsigned int *averef = NULL; ALLOCATE (averef , unsigned int, Nstars); for (i = 0; i < Nstars; i++) averef [i] = stars[i].starpar.averef ; + unsigned int *objID = NULL; ALLOCATE (objID , unsigned int, Nstars); for (i = 0; i < Nstars; i++) objID [i] = stars[i].starpar.objID ; + unsigned int *catID = NULL; ALLOCATE (catID , unsigned int, Nstars); for (i = 0; i < Nstars; i++) catID [i] = stars[i].starpar.catID ; + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", ra , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DEC", dec , Nstars); + gfits_set_bintable_column (&theader, &ftable, "GLON", glon , Nstars); + gfits_set_bintable_column (&theader, &ftable, "GLAT", glat , Nstars); + gfits_set_bintable_column (&theader, &ftable, "Ebv", Ebv , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dEbv", dEbv , Nstars); + gfits_set_bintable_column (&theader, &ftable, "DistMag", DistMag , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dDistMag", dDistMag, Nstars); + gfits_set_bintable_column (&theader, &ftable, "M_r", M_r , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dM_r", dM_r , Nstars); + gfits_set_bintable_column (&theader, &ftable, "FeH", FeH , Nstars); + gfits_set_bintable_column (&theader, &ftable, "dFeH", dFeH , Nstars); + gfits_set_bintable_column (&theader, &ftable, "averef", averef , Nstars); + gfits_set_bintable_column (&theader, &ftable, "objID", objID , Nstars); + gfits_set_bintable_column (&theader, &ftable, "catID", catID , Nstars); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + return TRUE; +} + + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +StarPar_Stars *loadstarpar_load_stars (char *filename, int *nstars) { + + int i, Ncol; + off_t Nrow; + char type[16]; + + StarPar_Stars *stars = NULL; + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + + 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 header\n"); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read matrix\n"); + goto escape; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + GET_COLUMN(ra , "RA" , double); + GET_COLUMN(dec , "DEC" , double); + GET_COLUMN(glon , "GLON" , float); + GET_COLUMN(glat , "GLAT" , float); + GET_COLUMN(Ebv , "Ebv" , float); + GET_COLUMN(dEbv , "dEbv" , float); + GET_COLUMN(DistMag , "DistMag" , float); + GET_COLUMN(dDistMag, "dDistMag" , float); + GET_COLUMN(M_r , "M_r" , float); + GET_COLUMN(dM_r , "dM_r" , float); + GET_COLUMN(FeH , "FeH" , float); + GET_COLUMN(dFeH , "dFeH" , float); + GET_COLUMN(averef , "averef" , unsigned int); + GET_COLUMN(objID , "objID" , unsigned int); + GET_COLUMN(catID , "catID" , unsigned int); + + gfits_free_header (&theader); + gfits_free_table (&ftable); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + ALLOCATE (stars, StarPar_Stars, Nrow); + + for (i = 0; i < Nrow; i++) { + stars[i].R = ra [i]; + stars[i].D = dec [i]; + stars[i].starpar.R = ra [i]; + stars[i].starpar.D = dec [i]; + stars[i].starpar.galLon = glon [i]; + stars[i].starpar.galLat = glat [i]; + stars[i].starpar.Ebv = Ebv [i]; + stars[i].starpar.dEbv = dEbv [i]; + stars[i].starpar.DistMag = DistMag [i]; + stars[i].starpar.dDistMag = dDistMag[i]; + stars[i].starpar.M_r = M_r [i]; + stars[i].starpar.dM_r = dM_r [i]; + stars[i].starpar.FeH = FeH [i]; + stars[i].starpar.dFeH = dFeH [i]; + stars[i].starpar.averef = averef [i]; + stars[i].starpar.objID = objID [i]; + stars[i].starpar.catID = catID [i]; + } + + fclose (f); + *nstars = Nrow; + return stars; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + fclose (f); + return NULL; +} Index: trunk/Ohana/src/addstar/src/loadstarpar_make_subset.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_make_subset.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_make_subset.c (revision 37807) @@ -0,0 +1,44 @@ +# include "addstar.h" +# include "loadstarpar.h" + +// assign stars in the given region to the subset (NOTE: stars are sorted by RA, start is +// first entry in stars array in this region) + +StarPar_Stars *loadstarpar_make_subset (StarPar_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) { + + int i; + + StarPar_Stars *subset = NULL; + + // collect array of (Stars *) stars in a new output catalog + int Nsubset = 0; + int NSUBSET = 3000; + ALLOCATE (subset, StarPar_Stars, NSUBSET); + + // find the rest of the stars in this output region + for (i = start; i < Nstars; i++) { + if (stars[i].flag) continue; + + // check if in skyregion + if (stars[i].R < region[0].Rmin) continue; + if (stars[i].R > region[0].Rmax) break; + if (stars[i].D < region[0].Dmin) continue; + if (stars[i].D > region[0].Dmax) continue; + + // check if in UserPatch (a GLOBAL) + if (stars[i].R < UserPatch.Rmin) continue; + if (stars[i].R > UserPatch.Rmax) break; + if (stars[i].D < UserPatch.Dmin) continue; + if (stars[i].D > UserPatch.Dmax) continue; + + subset[Nsubset] = stars[i]; + Nsubset ++; + + stars[i].flag = TRUE; + + CHECK_REALLOCATE (subset, StarPar_Stars, NSUBSET, Nsubset, 10000); + } + + *nsubset = Nsubset; + return subset; +} Index: trunk/Ohana/src/addstar/src/loadstarpar_readstars.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_readstars.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_readstars.c (revision 37807) @@ -0,0 +1,164 @@ +# include "addstar.h" +# include "loadstarpar.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +StarPar_Stars *loadstarpar_readstars (char *filename, int *nstars) { + + // read in the full FITS files ('cause I don't have a partial read option) + FILE *f = fopen (filename, "r"); + if (f == NULL) Shutdown ("can't read stellar parameter file: %s", filename); + + int i, Ncol; + off_t Nrow; + + Header header; + Matrix matrix; + Header theader; + FTable ftable; + + // 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; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + + char type[16]; + + GET_COLUMN (glat, "l", float); + GET_COLUMN (glon, "b", float); + GET_COLUMN (conv, "conv", byte); + GET_COLUMN (lnZ, "lnZ", float); + GET_COLUMN (DistMag, "DM", float); + GET_COLUMN (Ebv, "EBV", float); + GET_COLUMN (M_r, "Mr", float); + GET_COLUMN (FeH, "FeH", float); + + // free the memory associated with the FITS files + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + // the next FITS extension contains the error information + 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; + } + fclose (f); + + double *errImage = (double *) matrix.buffer; + + myAssert (Nrow == matrix.Naxis[2], "size mismatch?"); + + int NstarsIn = matrix.Naxis[2]; + + double Rmin = +360.0; + double Rmax = -360.0; + double Dmin = +360.0; + double Dmax = -360.0; + + int Nstars = 0; + int NSTARS = 0.1*NstarsIn; + + StarPar_Stars *stars = NULL; + ALLOCATE (stars, StarPar_Stars, NSTARS); + + // libdvo uses Liu et al 2011 (A&A 526, A16) for default galactic coords, + // but Greg Green / LSD use the older Reid et al 2004 (ApJ 616, 872) definition + CoordTransform *transform = InitTransform (COORD_GALACTIC_REID_2004, COORD_CELESTIAL); + + for (i = 0; i < NstarsIn; i++) { + + // skip stars based on conv and lnZ + if (!conv[i]) continue; + if (lnZ[i] < -15.0) continue; + + double R, D; + ApplyTransform (&R, &D, glon[i], glat[i], transform); + + Rmin = MIN (Rmin, R); + Rmax = MAX (Rmax, R); + Dmin = MIN (Dmin, D); + Dmax = MAX (Dmax, D); + + float dEbv = 0.5*(errImage[i*20 + 5*0 + 3] - errImage[i*20 + 5*0 + 1]); + float dDistMag = 0.5*(errImage[i*20 + 5*1 + 3] - errImage[i*20 + 5*1 + 1]); + float dM_r = 0.5*(errImage[i*20 + 5*2 + 3] - errImage[i*20 + 5*2 + 1]); + float dFeH = 0.5*(errImage[i*20 + 5*3 + 3] - errImage[i*20 + 5*3 + 1]); + + stars[Nstars].R = R; + stars[Nstars].D = D; + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + // NOTE that we have both stars.R,D and stars.starpar.R,D. stars.R,D are used + // here (in parallel with addstar) to locate the objects. BUT, in the database, + // starpar.R,D is the authoratative position for the object (either that coming + // from Greg / Eddie's table or that coming from the simulation) + stars[Nstars].starpar.R = R; + stars[Nstars].starpar.D = D; + stars[Nstars].starpar.galLon = glon[i]; + stars[Nstars].starpar.galLat = glat[i]; + + stars[Nstars].starpar.Ebv = Ebv[i] ; + stars[Nstars].starpar.dEbv = dEbv ; + stars[Nstars].starpar.DistMag = DistMag[i]; + stars[Nstars].starpar.dDistMag = dDistMag ; + stars[Nstars].starpar.M_r = M_r[i] ; + stars[Nstars].starpar.dM_r = dM_r ; + stars[Nstars].starpar.FeH = FeH[i] ; + stars[Nstars].starpar.dFeH = dFeH ; + + Nstars ++; + + CHECK_REALLOCATE (stars, StarPar_Stars, NSTARS, Nstars, 10000); + } + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + *nstars = Nstars; + return (stars); +} + +int loadstarpar_sortStars (StarPar_Stars *stars, int Nstars) { + +# define SWAPFUNC(A,B){ StarPar_Stars temp = stars[A]; stars[A] = stars[B]; stars[B] = temp; } +# define COMPARE(A,B)(stars[A].R < stars[B].R) + + OHANA_SORT (Nstars, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + + return TRUE; +} + Index: trunk/Ohana/src/addstar/src/loadstarpar_remote_hosts.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_remote_hosts.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_remote_hosts.c (revision 37807) @@ -0,0 +1,142 @@ +# include "addstar.h" +# include "loadstarpar.h" +# include +# include + +# define DEBUG 1 + +// we are running N parallel remote jobs on the remote hosts. In the usual remote +// processing (eg, relphot or relastro), we launch one job per remote host. In this case, +// we are launching one job foreach table (catalog.cpt) being touched. individual hosts +// may have more than one job active at a time. + +int Nhosts = 0; +HostInfo **hosts = NULL; + +int init_remote_hosts () { + + Nhosts = 10; + ALLOCATE (hosts, HostInfo *, Nhosts); + + int i; + for (i = 0; i < Nhosts; i++) { + hosts[i] = NULL; + } + return TRUE; +} + +int find_empty_slot () { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) return i; + } + return -1; +} + +int harvest_all () { + + int slot = -1; + + while (slot != -2) { + slot = harvest_host (); + usleep (50000); + } + return TRUE; +} + +// put the host in the list in a free slot. NOTE: this should never be called if we do not +// already have a free slot. +int save_remote_host (HostInfo *host) { + + int i; + for (i = 0; i < Nhosts; i++) { + if (hosts[i] == NULL) { + hosts[i] = host; + return TRUE; + } + } + myAbort ("failed to find an empty slot: this is a programming error"); + return FALSE; +} + +// wait for all children to complete, report output to stdout +// possible states when function is finished: +// * no child was harvested, but children still exist +// * no outstanding children +// * a child was harvested +int harvest_host () { + + // check if any children have finished... + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + + if (!pid) return -1; // there are outstanding children, but none have exited, no slot was opened + if ((pid == -1) && (errno == ECHILD)) return -2; // there are no outstanding children + if ((pid == -1) && (errno != ECHILD)) myAbort ("programming error?"); + + // find the host which has finished + int i; + int found = FALSE; + for (i = 0; (i < Nhosts) && !found; i++) { + if (!hosts[i]) continue; // unassigned slot + if (hosts[i][0].pid != pid) continue; + found = TRUE; + break; + } + myAssert (found, "Programming error: failed to matched finished job to known host!"); + int slot = i; + + HostInfo *host = hosts[slot]; + + // check on the status of this job and report any output + if (DEBUG) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + int Nout, printHead; + + // read stdout + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (DEBUG) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + } + + // close opened connections + close (host->stdio[HOST_STDIN]); + close (host->stdio[HOST_STDOUT]); + close (host->stdio[HOST_STDERR]); + + // free data associated with the host + free (host->hostname); + free (host->pathname); + FreeIOBuffer (&host->stdout); + FreeIOBuffer (&host->stderr); + free (host); + + // free the slot + hosts[slot] = NULL; + return slot; +} + Index: trunk/Ohana/src/addstar/src/loadstarpar_save_remote.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_save_remote.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_save_remote.c (revision 37807) @@ -0,0 +1,84 @@ +# include "addstar.h" +# include "loadstarpar.h" +# define DEBUG 1 + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} + +int loadstarpar_save_remote (StarPar_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // if this region is a parallel thing, save and launch remote + if (!PARALLEL) { + loadstarpar_catalog (stars, Nstars, region, fullname, options); + } else { + int N = hosts->index[region->hostID]; + HostInfo *hostMach = &hosts->hosts[N]; + + // save to a unique filename + char filename[1024]; // CATDIR/tmpdir/starpar.PID.index.fits + snprintf (filename, 1024, "%s/tmpdir/starpar.%s.%05d.fits", CATDIR, uniquer, region->index); + + // write the data to the given FITS file + loadstarpar_save_stars (filename, stars, Nstars); + + int slot = -1; + while (slot == -1) { + slot = find_empty_slot (); + if (slot == -1) { + usleep (50000); + slot = harvest_host(); + myAssert (slot != -2, "we should not call harvest_host here if we have open slots"); + } + } + + // allocate a host for this job + HostInfo *host = NULL; + ALLOCATE (host, HostInfo, 1); + + // we want to run this job on the host described by hostMach. copy + // immutable data from hostMach to a locally allocated host: + host->hostID = hostMach->hostID; + host->hostname = strcreate (hostMach->hostname); + host->pathname = strcreate (hostMach->pathname); + InitIOBuffer (&host->stdout, 1000); + InitIOBuffer (&host->stderr, 1000); + + // got a valid slot, so launch a new host + + // need to generate the remote command + char command[1024]; + snprintf (command, 1024, "loadstarpar_client"); + strextend (command, "-hostID %d", host->hostID); + strextend (command, "-D CATDIR %s", CATDIR); + strextend (command, "-hostdir %s", host->pathname); + strextend (command, "-cpt %s", region->name); + strextend (command, "-input %s", filename); + + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", host->hostname, command, host->stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", host->hostname, errorInfo); + exit (1); + } + host->pid = pid; // save for future reference + + save_remote_host (host); + } + return TRUE; +} Index: trunk/Ohana/src/addstar/src/loadstarpar_table.c =================================================================== --- trunk/Ohana/src/addstar/src/loadstarpar_table.c (revision 37807) +++ trunk/Ohana/src/addstar/src/loadstarpar_table.c (revision 37807) @@ -0,0 +1,50 @@ +# include "addstar.h" +# include "loadstarpar.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +int loadstarpar_table (SkyList *skylistInput, HostTable *hosts, char *filename, AddstarClientOptions *options) { + + int i, Nstars; + + StarPar_Stars *stars = loadstarpar_readstars (filename, &Nstars); + + // sort the stars by RA + loadstarpar_sortStars (stars, Nstars); + + // scan through the stars, loading the containing catalogs + // skip through table for unsaved stars + for (i = 0; i < Nstars; i++) { + if (stars[i].flag) continue; + + // scan forward until we read the UserPatch + if (stars[i].R < UserPatch.Rmin) continue; + if (stars[i].R > UserPatch.Rmax) break; + if (stars[i].D < UserPatch.Dmin) continue; + if (stars[i].D > UserPatch.Dmax) continue; + + // identify the relevant catalog + SkyList *skylist = SkyRegionByPoint_List (skylistInput, -1, stars[i].R, stars[i].D); + if (skylist[0].Nregions == 0) { + SkyListFree (skylist); + continue; + } + SkyRegion *region = skylist[0].regions[0]; + + // select stars matching this region + int Nsubset; + StarPar_Stars *subset = loadstarpar_make_subset (stars, Nstars, i, region, &Nsubset); + + // In parallel mode, write out the subset to a disk file. Block until a remote host + // is available. In serial mode, just match against the appropriate region and save + loadstarpar_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options); + } + + // wait for last remote clients to finish + harvest_all (); + + return TRUE; +} + Index: trunk/Ohana/src/addstar/src/loadsupercos_plates.c =================================================================== --- trunk/Ohana/src/addstar/src/loadsupercos_plates.c (revision 37729) +++ trunk/Ohana/src/addstar/src/loadsupercos_plates.c (revision 37807) @@ -118,5 +118,5 @@ // for now, we define a fake coordinate system based on the boresite center - strcpy (image[Nimage].coords.ctype, "DEC--TAN"); + InitCoords (&image[Nimage].coords, "DEC--TAN"); image[Nimage].coords.crval1 = RAo; @@ -126,12 +126,4 @@ image[Nimage].coords.crpix2 = 0.5*image[Nimage].NY; image[Nimage].coords.cdelt1 = image[Nimage].coords.cdelt2 = scale; - - image[Nimage].coords.pc1_1 = 0.0; - image[Nimage].coords.pc1_2 = 1.0; - image[Nimage].coords.pc2_1 = 1.0; - image[Nimage].coords.pc2_2 = 0.0; - - image[Nimage].coords.Npolyterms = 0; - memset (image[Nimage].coords.polyterms, 0, 2*7*sizeof(float)); Nimage ++; Index: trunk/Ohana/src/addstar/src/loadwise.c =================================================================== --- trunk/Ohana/src/addstar/src/loadwise.c (revision 37729) +++ trunk/Ohana/src/addstar/src/loadwise.c (revision 37807) @@ -40,5 +40,5 @@ for (i = 1; i < argc; i++) { fprintf (stderr, "loading %s\n", argv[i]); - loadwise_rawdata (skylist, argv[i], options, ALLSKY); + loadwise_rawdata (skylist, argv[i], options); } exit (0); Index: trunk/Ohana/src/addstar/src/loadwise_rawdata.c =================================================================== --- trunk/Ohana/src/addstar/src/loadwise_rawdata.c (revision 37729) +++ trunk/Ohana/src/addstar/src/loadwise_rawdata.c (revision 37807) @@ -10,5 +10,5 @@ # define DEBUG 0 -int loadwise_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options, int ALLSKY) { +int loadwise_rawdata (SkyList *skytable, char *filename, AddstarClientOptions options) { int i, j, verbose; @@ -159,8 +159,15 @@ stars[Nstars+3][0].average.D = tstars[j].D; - if (ALLSKY) { - loadwise_star_allsky (&stars[Nstars], &buffer[offset], Nbyte - offset); - } else { - loadwise_star_prelim (&stars[Nstars], &buffer[offset], Nbyte - offset); + switch (MODE) { + case MODE_PRELIM: + loadwise_star_prelim (&stars[Nstars], &buffer[offset], Nbyte - offset); + break; + case MODE_ALLSKY: + loadwise_star_allsky (&stars[Nstars], &buffer[offset], Nbyte - offset); + break; + case MODE_ALLWISE: + default: + // loadwise_star_allwise (&stars[Nstars], &buffer[offset], Nbyte - offset); + break; } Index: trunk/Ohana/src/addstar/src/mkcmf.c =================================================================== --- trunk/Ohana/src/addstar/src/mkcmf.c (revision 37729) +++ trunk/Ohana/src/addstar/src/mkcmf.c (revision 37807) @@ -200,5 +200,5 @@ /* bore site center guess */ - strcpy (coords.ctype, "DEC--TAN"); + InitCoords (&coords, "DEC--TAN"); coords.crval1 = RA; coords.crval2 = DEC; @@ -208,10 +208,4 @@ coords.cdelt1 = 0.25/3600.0; coords.cdelt2 = 0.25/3600.0; - - coords.pc1_1 = 1; - coords.pc2_2 = 1; - coords.pc1_2 = 0; - coords.pc2_1 = 0; - coords.Npolyterms = 1; // load stars and generate complete output fields @@ -300,6 +294,6 @@ gfits_modify (&header, "NASTRO", "%d", 1, 100); - if (imageID == -1) imageID = 1000.0*drand48(); - if (sourceID == -1) sourceID = 100.0*drand48(); + if (imageID == -1) imageID = 100000.0*drand48(); + if (sourceID == -1) sourceID = 1000.0*drand48(); gfits_modify (&header, "IMAGEID", "%d", 1, imageID); gfits_modify (&header, "SOURCEID", "%d", 1, sourceID); Index: trunk/Ohana/src/addstar/src/resort_catalog.c =================================================================== --- trunk/Ohana/src/addstar/src/resort_catalog.c (revision 37729) +++ trunk/Ohana/src/addstar/src/resort_catalog.c (revision 37807) @@ -4,4 +4,5 @@ void resort_catalog_measure (Catalog *catalog); void resort_catalog_lensing (Catalog *catalog); +void resort_catalog_starpar (Catalog *catalog); void resort_catalog_old (Catalog *catalog) { @@ -38,4 +39,5 @@ resort_catalog_measure (catalog); resort_catalog_lensing (catalog); + resort_catalog_starpar (catalog); } @@ -335,4 +337,141 @@ } +void resort_catalog_starpar (Catalog *catalog) { + + off_t Naverage, Nstarpar; + StarPar *starpar; + Average *average; + off_t i, j, N, currentAve; + + off_t *starparSeq = NULL; + off_t *averageSeq = NULL; + StarPar *starparTMP = NULL; + + // struct timeval start, stop; + // gettimeofday (&start, NULL); + + /* internal counters */ + Nstarpar = catalog[0].Nstarpar; + Naverage = catalog[0].Naverage; + + if (!Nstarpar) return; + + starpar = catalog[0].starpar; + average = catalog[0].average; + + // we have a table of average objects and an unsorted table of measurements. each measurement + // has a reference to the average object sequence (as well as an ID) + // starpar[i].averef -> average[averef] + // starpar[i].objID = average[averef].objID + // starpar[i].catID = average[averef].catID + + // we want a sorted starpar array with all averef entries in sequence + + ALLOCATE (starparSeq, off_t, Nstarpar); + ALLOCATE (averageSeq, off_t, Nstarpar); + + for (i = 0; i < Nstarpar; i++) { + starparSeq[i] = i; + averageSeq[i] = starpar[i].averef; + + if (catalog[0].catformat >= DVO_FORMAT_PS1_V1) { + // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't) + myAssert(average[averageSeq[i]].catID == starpar[starparSeq[i]].catID, "object / detection mismatch"); +# if (1) + myAssert(average[averageSeq[i]].objID == starpar[starparSeq[i]].objID, "object / detection mismatch"); +# else + // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID. + // this happened for 28 detections in the dbs on /data/stsci1?.0/eugene/dvo3pi.20130616, but not at all (as far as I know) in the rest of LAP DVO + if (average[averageSeq[i]].objID != starpar[starparSeq[i]].objID) { + fprintf (stderr, "R"); + starpar[starparSeq[i]].objID = average[averageSeq[i]].objID; // XXX I don't really like this... + } +# endif + } + } + + // check that averageSeq is now in order + // for (i = 1; i < Nstarpar; i++) { + // if (averageSeq[i] < averageSeq[i-1]) { + // fprintf (stderr, "%d ", (int) i); + // } + // } + // fprintf (stderr, "\n"); + + SortAveMatch(starparSeq, averageSeq, Nstarpar); + // MARKTIME("sort : %f sec\n", dtime); + + // check that averageSeq is now in order + // for (i = 1; i < Nstarpar; i++) { + // if (averageSeq[i] < averageSeq[i-1]) { + // fprintf (stderr, "%d ", (int) i); + // } + // } + // fprintf (stderr, "\n"); + + // copy the starpar entries in the sorted order + ALLOCATE (starparTMP, StarPar, Nstarpar); + for (i = 0; i < Nstarpar; i++) { + j = starparSeq[i]; + starparTMP[i] = starpar[j]; + } + // MARKTIME("assign starpar : %f sec\n", dtime); + + // update the values of average.starparOffset and average.Nstarpar + FREE(starpar); + catalog[0].starpar = starparTMP; + + N = 0; + currentAve = averageSeq[0]; + average[currentAve].starparOffset = 0; + for (i = 0; i < Nstarpar; i++) { + if (averageSeq[i] != currentAve) { + // we have hit the next entry in the list + average[currentAve].Nstarpar = N; + N = 0; + currentAve = averageSeq[i]; + average[currentAve].starparOffset = i; + } + N++; + } + // N++; + average[currentAve].Nstarpar = N; + // MARKTIME("update Nstarpar : %f sec\n", dtime); + + int NstarparTotal = 0; + int starparOffsetOK = TRUE; + for (i = 0; i < Naverage; i++) { + NstarparTotal += catalog[0].average[i].Nstarpar; + if (VERBOSE && !(NstarparTotal <= catalog[0].Nstarpar)) { + fprintf (stderr, "too few starpar: %d %d %d\n", (int) i, NstarparTotal, (int) catalog[0].Nstarpar); + } + starparOffsetOK &= (catalog[0].average[i].starparOffset < catalog[0].Nstarpar); + if (VERBOSE && !(catalog[0].average[i].starparOffset < catalog[0].Nstarpar)) { + fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Nstarpar, (int) catalog[0].Nstarpar); + } + starparOffsetOK &= (catalog[0].average[i].starparOffset + catalog[0].average[i].Nstarpar <= catalog[0].Nstarpar); + if (VERBOSE && !(catalog[0].average[i].starparOffset + catalog[0].average[i].Nstarpar <= catalog[0].Nstarpar)) { + fprintf (stderr, "orrset + Nstarpar too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].starparOffset, catalog[0].average[i].Nstarpar, (int) catalog[0].Nstarpar); + } + } + + if (!starparOffsetOK) { + fprintf (stderr, "ERROR: catalog %s has an invalid starparOffset\n", catalog[0].filename); + } + + if (NstarparTotal != catalog[0].Nstarpar) { + fprintf (stderr, "ERROR: catalog %s has an invalid Nstarpar\n", catalog[0].filename); + } + + // MARKTIME(" match time %9.4f sec for %7lld starpars, %6lld average\n", dtime, (long long) Nstarpar, (long long) Naverage); + + catalog[0].sorted = TRUE; + + FREE (starparSeq); + FREE (averageSeq); + + return; +} + // sort the measure or lensing Sequence based on the average Sequence entries void SortAveMatch (off_t *MEAS, off_t *AVE, off_t N) { Index: trunk/Ohana/src/addstar/src/sky_tessalation.c =================================================================== --- trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 37729) +++ trunk/Ohana/src/addstar/src/sky_tessalation.c (revision 37807) @@ -581,4 +581,7 @@ memset (rectangle, 0, sizeof(SkyRectangle)); + + InitCoords (&rectangle[0].coords, "DEC--TAN"); + rectangle[0].coords.crval1 = CENTER_RA; rectangle[0].coords.crval2 = CENTER_DEC; @@ -609,6 +612,4 @@ rectangle[0].coords.cdelt2 = SCALE / 3600.0; - strcpy (rectangle[0].coords.ctype, "DEC--TAN"); - rectangle[0].NX = NX; rectangle[0].NY = NY; @@ -699,12 +700,10 @@ for (i = 0; i < nRA; i++) { memset (&ring[i], 0, sizeof(SkyRectangle)); - memset (&ring[i].coords, 0, sizeof(Coords)); + + InitCoords (&ring[i].coords, "DEC--TAN"); ring[i].coords.crval1 = i*dRA; ring[i].coords.crval2 = dec; ring[i].coords.pc1_1 = +1.0 * X_PARITY; - ring[i].coords.pc1_2 = +0.0; - ring[i].coords.pc2_1 = -0.0; - ring[i].coords.pc2_2 = +1.0; // range values are in projected degrees @@ -718,6 +717,4 @@ ring[i].coords.cdelt1 = SCALE / 3600.0; ring[i].coords.cdelt2 = SCALE / 3600.0; - - strcpy (ring[i].coords.ctype, "DEC--TAN"); ring[i].NX = NX*(1.0 + PADDING); @@ -771,5 +768,5 @@ memset (&ring[N], 0, sizeof(SkyRectangle)); - memset (&ring[N].coords, 0, sizeof(Coords)); + InitCoords (&ring[N].coords, "DEC--TAN"); ring[N].coords.crval1 = ra / d2r; @@ -779,7 +776,4 @@ ring[N].coords.pc1_1 = +1.0 * X_PARITY; - ring[N].coords.pc1_2 = +0.0; - ring[N].coords.pc2_1 = -0.0; - ring[N].coords.pc2_2 = +1.0; // range values are in projected degrees @@ -793,6 +787,4 @@ ring[N].coords.cdelt1 = SCALE / 3600.0; ring[N].coords.cdelt2 = SCALE / 3600.0; - - strcpy (ring[N].coords.ctype, "DEC--TAN"); ring[N].NX = NX*(1.0 + PADDING); @@ -1099,12 +1091,6 @@ ALLOCATE (refcoords, Coords, 1); - refcoords[0].crval1 = refcoords[0].crval2 = 0.0; - refcoords[0].crpix1 = refcoords[0].crpix2 = 0.0; + InitCoords (refcoords, "DEC--TAN"); refcoords[0].cdelt1 = refcoords[0].cdelt2 = scale / 3600; - refcoords[0].pc1_1 = refcoords[0].pc2_2 = 1.0; - refcoords[0].pc1_2 = refcoords[0].pc2_1 = 0.0; - refcoords[0].Npolyterms = 0; - memset (refcoords[0].polyterms, 0, 14*sizeof(float)); - strcpy (refcoords[0].ctype, "DEC--TAN"); return (TRUE); } Index: trunk/Ohana/src/addstar/test/relphot.parallel.dvo =================================================================== --- trunk/Ohana/src/addstar/test/relphot.parallel.dvo (revision 37729) +++ trunk/Ohana/src/addstar/test/relphot.parallel.dvo (revision 37807) @@ -25,17 +25,23 @@ if (not($?VERBOSE)) set VERBOSE = 0 - # images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters - $zpt_nominal:g = 24.58 - $zpt_nominal:r = 24.80 - $zpt_nominal:i = 24.74 - $zpt_nominal:z = 24.26 - $zpt_nominal:y = 23.41 - - # klam is negative, so klam*(airmass - 1) increase the zero point - $klam_nominal:g = -0.15 - $klam_nominal:r = -0.10 - $klam_nominal:i = -0.04 - $klam_nominal:z = -0.03 - $klam_nominal:y = -0.03 +# images are loaded into dvo with GPC1 photcodes, so we need to get the nominal zps for those filters +# $zpt_nominal:g = 24.58 +# $zpt_nominal:r = 24.80 +# $zpt_nominal:i = 24.74 +# $zpt_nominal:z = 24.26 +# $zpt_nominal:y = 23.41 + +# klam is negative, so klam*(airmass - 1) increase the zero point +# $klam_nominal:g = -0.15 +# $klam_nominal:r = -0.10 +# $klam_nominal:i = -0.04 +# $klam_nominal:z = -0.03 +# $klam_nominal:y = -0.03 + + $zpt_nominal:g = 24.563; $klam_nominal:g = -0.147 + $zpt_nominal:r = 24.750; $klam_nominal:r = -0.085 + $zpt_nominal:i = 24.611; $klam_nominal:i = -0.044 + $zpt_nominal:z = 24.250; $klam_nominal:z = -0.033 + $zpt_nominal:y = 23.320; $klam_nominal:y = -0.073 # we have two sets of images: ubercaled and not-ubercaled @@ -111,5 +117,5 @@ mksequence $fileroot $catdir for i 0 mjd_nc[] - ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw + ckexposure $catdir mjd_nc[$i] zpt_nc[$i] exptime_nc[$i] airmass_nc[$i] $filt_nc:$i raw 0 end @@ -121,5 +127,5 @@ # run relphot on the serial db and check that the images now match the expected values - exec relphot g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 + exec relphot -images g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -D IMAGE_OFFSET 0.5 -update >& log.relphot.s0 for i 0 mjd_nc[] @@ -128,11 +134,9 @@ # matching the mean zp of photometric data, which we calculate in # 'init' - ckexposure $catdir mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + ckexposure $catdir mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc 0 end # run relphot on the parallel db and check that the images now match the expected values - exec relphot -parallel g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir.p0 -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -keep-ubercal -D IMAGE_OFFSET 0.5 -update >& log.relphot.p0 - - exec dvodist -in $catdir.p0 >& log.dvodist.in + exec relphot -images -parallel g,r,i -v -region 9.5 10.5 19.5 20.5 -D CATDIR $catdir.p0 -D STAR_TOOFEW 1 -D SIGMA_LIM 0.07 -statmode WT_MEAN -cloud-limit 0.5 -D IMAGE_OFFSET 0.5 -update >& log.relphot.p0 for i 0 mjd_nc[] @@ -141,5 +145,5 @@ # matching the mean zp of photometric data, which we calculate in # 'init' - ckexposure $catdir.p0 mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc + ckexposure $catdir.p0 mjd_nc[$i] $zpt_nc_$filt_nc:$i exptime_nc[$i] airmass_nc[$i] $filt_nc:$i relphot_nc 1 end @@ -148,6 +152,6 @@ macro ckexposure - if ($0 != 8) - echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode)" + if ($0 != 9) + echo "ckexposure (catdir) (mjd) (zpt) (exptime) (airmass) (filter) (mode) (isParallel)" echo " mode == raw or corr" break @@ -163,4 +167,8 @@ $FILTER = $6 $MODE = $7 + $PARALLEL = + if ($8) + $PARALLEL = "-parallel" + end # XXX need a function to extract the nominal zpt for a given filter / photcode from the db @@ -182,5 +190,5 @@ skyregion {$RA_CENTER - 0.2/dcos($DEC_CENTER)} {$RA_CENTER + 0.2/dcos($DEC_CENTER)} {$DEC_CENTER - 0.2} {$DEC_CENTER + 0.2} - mextract ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) + mextract $PARALLEL ra dec mag xccd yccd where (abs(time - $MJD_IMAGE) < 0.0001) if (not(ra[])) echo "no matching data for MJD = $MJD_IMAGE" @@ -344,5 +352,5 @@ # the fake images have inconsistent ra,dec and airmass,sidtime values if ($VERBOSE) echo addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass - exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& /dev/null + exec addstar $ROOT.$ix.$iy.cmf -D CATDIR $CATDIR -D CAMERA gpc1 -quick-airmass >& log.addstar # costs 10/36 sec / chip end @@ -432,4 +440,5 @@ $catdir = $tmp1/$tmp2 + if ($VERBOSE) echo rsync -auv $catdir/ $catdir.p0/ exec rsync -auv $catdir/ $catdir.p0/ >& log.rsync mkdir $catdir.d1 @@ -442,4 +451,5 @@ exec echo "3 $hostname $catdir.d3" >> $catdir.p0/HostTable.dat + if ($VERBOSE) echo dvodist -out $catdir.p0 exec dvodist -out $catdir.p0 >& log.dvodist.out end