Index: trunk/Ohana/src/addstar/Makefile =================================================================== --- trunk/Ohana/src/addstar/Makefile (revision 40542) +++ trunk/Ohana/src/addstar/Makefile (revision 40543) @@ -26,4 +26,5 @@ loadgalphot : $(BIN)/loadgalphot.$(ARCH) loadgaia : $(BIN)/loadgaia.$(ARCH) +loadgaia_dr2 : $(BIN)/loadgaia_dr2.$(ARCH) loadstarpar : $(BIN)/loadstarpar.$(ARCH) loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH) @@ -43,5 +44,5 @@ # programs in 'SERVER' use the client-server concept and are out of date -INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos +INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadgaia_dr2 loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos SERVER = addstarc addstard addstart @@ -286,4 +287,22 @@ $(SRC)/psps_ids.$(ARCH).o +LOAD-GAIA_DR2 = \ +$(SRC)/loadgaia_dr2.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args_loadgaia_dr2.$(ARCH).o \ +$(SRC)/find_matches_gaia_dr2.$(ARCH).o \ +$(SRC)/loadgaia_dr2_catalog.$(ARCH).o \ +$(SRC)/loadgaia_dr2_make_subset.$(ARCH).o \ +$(SRC)/loadgaia_dr2_readstars.$(ARCH).o \ +$(SRC)/loadgaia_dr2_table.$(ARCH).o \ +$(SRC)/parse_csv.$(ARCH).o \ +$(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/build_links.$(ARCH).o \ +$(SRC)/strhash.$(ARCH).o \ +$(SRC)/sortIDs.$(ARCH).o \ +$(SRC)/psps_ids.$(ARCH).o + LOAD-STARPAR = \ $(SRC)/loadstarpar.$(ARCH).o \ @@ -508,4 +527,5 @@ $(LOAD-GALPHOT) : $(INC)/addstar.h $(INC)/loadgalphot.h $(LOAD-GAIA) : $(INC)/addstar.h $(INC)/gaia.h +$(LOAD-GAIA_DR2) : $(INC)/addstar.h $(INC)/gaia_dr2.h $(LOAD-STARPAR) : $(INC)/addstar.h $(INC)/loadstarpar.h $(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h @@ -530,4 +550,5 @@ $(BIN)/loadgalphot.$(ARCH) : $(LOAD-GALPHOT) $(BIN)/loadgaia.$(ARCH) : $(LOAD-GAIA) +$(BIN)/loadgaia_dr2.$(ARCH) : $(LOAD-GAIA_DR2) $(BIN)/loadstarpar.$(ARCH) : $(LOAD-STARPAR) $(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT) Index: trunk/Ohana/src/addstar/include/gaia_dr2.h =================================================================== --- trunk/Ohana/src/addstar/include/gaia_dr2.h (revision 40543) +++ trunk/Ohana/src/addstar/include/gaia_dr2.h (revision 40543) @@ -0,0 +1,50 @@ + +// measure[0] = g, measure[1] = b, measure[2] = r, + +typedef struct { + Average average; + Measure measure[3]; + int flag; // in a subset? + int found; // assigned to an object? +} Gaia_DR2_Stars; + +int HOST_ID; +char *HOSTDIR; +char *CPT_FILE; +char *INPUT; + +AddstarClientOptions args_loadgaia_dr2 (int *argc, char **argv, AddstarClientOptions options); +// AddstarClientOptions args_loadgaia_dr2_client (int *argc, char **argv, AddstarClientOptions options); + +int loadgaia_dr2_table (int Nstart, int Nend, SkyList *skylistInput, HostTable *hosts, char **filename, AddstarClientOptions *options); + +Gaia_DR2_Stars *loadgaia_dr2_make_subset (Gaia_DR2_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset); + +// int loadgaia_dr2_save_remote (Gaia_DR2_Stars *stars, int Nstars, HostTable *hosts, SkyRegion *region, char *fullname, AddstarClientOptions *options); +// int save_remote_host (HostInfo *host); + +// int init_remote_hosts (void); +// void free_remote_hosts (void); +// int find_empty_slot (void); +// int harvest_all (void); +// int harvest_host (void); + +int loadgaia_dr2_catalog (Gaia_DR2_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_gaia_dr2 (SkyRegion *region, Gaia_DR2_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options); + +int loadgaia_dr2_save_stars (char *filename, Gaia_DR2_Stars *stars, int Nstars); +// Gaia_DR2_Stars *loadgaia_dr2_load_stars (char *filename, int *nstars); + +Gaia_DR2_Stars *loadgaia_dr2_readstars (char *filename, Gaia_DR2_Stars *stars, int *nstars, AddstarClientOptions *options); + +int loadgaia_dr2_sortStars (Gaia_DR2_Stars *stars, int Nstars); + +// int loadgaia_dr2_tmpdir (void); + +// roll these into a function with a structure to carry the Nlast, etc? +char *dparse_csv_rpt (double *X, int Nwant, int Nlast, char *line, int *status); +char *iparse_csv_rpt (int *X, int Nwant, int Nlast, char *line, int *status); +char *jparse_csv_rpt (uint64_t *X, int Nwant, int Nlast, char *line, int *status); Index: trunk/Ohana/src/addstar/src/args_loadgaia_dr2.c =================================================================== --- trunk/Ohana/src/addstar/src/args_loadgaia_dr2.c (revision 40543) +++ trunk/Ohana/src/addstar/src/args_loadgaia_dr2.c (revision 40543) @@ -0,0 +1,137 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +static void help (void); + +AddstarClientOptions args_loadgaia_dr2 (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); + } + /* override any header PHOTCODE values */ + options.photcode = 0; + + /* accept proper-motion (& parallax) data from reference */ + ACCEPT_MOTION = FALSE; + if ((N = get_argument (*argc, argv, "-accept-motion"))) { + ACCEPT_MOTION = TRUE; + remove_argument (N, argc, argv); + } + + /* provide a time for dataset */ + options.timeref = 0; + if ((N = get_argument (*argc, argv, "-time"))) { + time_t tmp; + remove_argument (N, argc, argv); + if (!ohana_str_to_time (argv[N], &tmp)) { + fprintf (stderr, "syntax error in time\n"); + exit (1); + } + options.timeref = tmp; + 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 loadgaia_dr2 */ + // 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: loadgaia_dr2 [options] (fitsfile) [..more files]\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadgaia_dr2 [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); +} Index: trunk/Ohana/src/addstar/src/find_matches_gaia_dr2.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_gaia_dr2.c (revision 40543) +++ trunk/Ohana/src/addstar/src/find_matches_gaia_dr2.c (revision 40543) @@ -0,0 +1,346 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +int find_matches_gaia_dr2 (SkyRegion *region, Gaia_DR2_Stars *stars, int NstarsIn, Catalog *tgtcat, AddstarClientOptions *options) { + + int Nsecfilt = GetPhotcodeNsecfilt (); + int Nsec = GetPhotcodeNsec (options->photcode); + + /** allocate local arrays (stars) **/ + ALLOCATE_PTR (X1, double, NstarsIn); + ALLOCATE_PTR (Y1, double, NstarsIn); + ALLOCATE_PTR (N1, off_t, NstarsIn); + + /** allocate local arrays (tgtcat) **/ + off_t NAVE = tgtcat[0].Naverage; + off_t Nave = tgtcat[0].Naverage; + ALLOCATE_PTR (X2, double, NAVE); + ALLOCATE_PTR (Y2, double, NAVE); + ALLOCATE_PTR (N2, off_t, NAVE); + + /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */ + + /* internal counters */ + off_t Nmatch = 0; + off_t NMEAS = tgtcat[0].Nmeasure; + off_t Nmeas = tgtcat[0].Nmeasure; + + off_t *next_meas = NULL; + + // current max obj ID for this tgtcat + unsigned int objID = tgtcat[0].objID; + unsigned int catID = tgtcat[0].catID; + + /* 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 + * substantially far from the projection pole. We use the center of the region (tgtcat) + * for crval1,2. + */ + + Coords tcoords; + InitCoords (&tcoords, "DEC--ARC"); + tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); + if ((region[0].Dmax < 90) && (region[0].Dmin > -90)) { + tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax); + } else { + tcoords.crval2 = (region[0].Dmax >= 90) ? 90.0 : -90.0; + } + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + + /* build spatial index (RA sort) referencing input array sequence */ + off_t Nstars = 0; + for (off_t i = 0; i < NstarsIn; i++) { + int status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].average.R, stars[i].average.D, &tcoords); + if (!status) continue; + N1[Nstars] = i; + Nstars ++; + } + if (Nstars < 1) { + if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", tgtcat[0].filename); + free (X1); + free (Y1); + free (N1); + free (X2); + free (Y2); + free (N2); + return 0; + } + if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars); + + /* build spatial index (RA sort) */ + for (off_t i = 0; i < Nave; i++) { + RD_to_XY (&X2[i], &Y2[i], tgtcat[0].average[i].R, tgtcat[0].average[i].D, &tcoords); + N2[i] = i; + } + if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave); + + /* set up pointers for linked list of measure */ + if (tgtcat[0].sorted && (tgtcat[0].Nmeasure == tgtcat[0].Nmeasure_disk)) { + // this version is only valid if we have done a full tgtcat load, and if the tgtcat + // is sorted while processed + next_meas = init_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas); + } else { + next_meas = build_measure_links (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas); + } + + /* choose a radius for matches */ + double RADIUS = (options->radius == 0) ? 2.0 : options->radius; /* provided by config */ + double RADIUS2 = RADIUS*RADIUS; + + /****************** find matched stars ********************/ + + for (off_t i = 0, 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 */ + double 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; + } + + // skip this star if already assigned to an object + if (stars[N1[i]].found) { + 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 */ + off_t Jmin = -1; + double Rmin = RADIUS2; + for (off_t J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) { + /* find closest match for this detection */ + dX = X1[i] - X2[J]; + double dY = Y1[i] - Y2[J]; + double 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 average, measure ***/ + Nmatch ++; + off_t n = N2[Jmin]; + off_t N = N1[i]; + + /* make sure there is space for next entry */ + if (Nmeas >= NMEAS - stars[N].average.Nmeasure) { + NMEAS = Nmeas + 1000 + stars[N].average.Nmeasure; + REALLOCATE (next_meas, off_t, NMEAS); + REALLOCATE (tgtcat[0].measure, Measure, NMEAS); + } + + /** add measurements for this star **/ + for (int mSrc = 0; mSrc < stars[N].average.Nmeasure; mSrc++) { + + /* add to end of measurement list */ + add_meas_link (&tgtcat[0].average[n], next_meas, Nmeas, NMEAS); + + /** add measurements for this star **/ + + // set the new measurements + tgtcat[0].measure[Nmeas] = stars[N].measure[mSrc]; + + // measure now carries R,D (not dR,dD) + tgtcat[0].measure[Nmeas].dbFlags = 0; + tgtcat[0].measure[Nmeas].averef = n; + tgtcat[0].measure[Nmeas].objID = tgtcat[0].average[n].objID; + tgtcat[0].measure[Nmeas].catID = tgtcat[0].catID; + + float dRoff = dvoOffsetR(&tgtcat[0].measure[Nmeas], &tgtcat[0].average[n]); + + // rationalize dR + if (dRoff > +180.0*3600.0) { + // average on high end of boundary, move star up + tgtcat[0].measure[Nmeas].R += 360.0; + dRoff -= 360.0*3600.0; + } + if (dRoff < -180.0*3600.0) { + // average on low end of boundary, move star down + tgtcat[0].measure[Nmeas].R -= 360.0; + dRoff += 360.0*3600.0; + } + if (fabs(dRoff) > 10*RADIUS) { + // take declination into account and check again. + double cosD = cos(RAD_DEG*tgtcat[0].average[n].D); + if (fabs(dRoff*cosD) > 10*RADIUS) { + fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", + tgtcat[0].average[n].R, tgtcat[0].average[n].D, + stars[N].average.R, stars[N].average.D, + X1[i], X2[Jmin], + Y1[i], Y2[Jmin]); + // XXX abort on this? -- this is a bad failure... + } + } + + /* Nm is updated, but not written out in -update mode (for existing entries) + Nm is recalculated in build_meas_links if loaded table is not sorted */ + tgtcat[0].average[n].Nmeasure ++; + Nmeas ++; + + /* if we choose to flag close encounters, see find_matches.c */ + /* if we choose to calculate RA,DEC averages, see update_coords.c */ + } + + if (ACCEPT_MOTION) { + tgtcat[0].average[n].uR = stars[N].average.uR; + tgtcat[0].average[n].uD = stars[N].average.uD; + tgtcat[0].average[n].duR = stars[N].average.duR; + tgtcat[0].average[n].duD = stars[N].average.duD; + tgtcat[0].average[n].P = stars[N].average.P; + tgtcat[0].average[n].dP = stars[N].average.dP; + tgtcat[0].average[n].Tmean = stars[N].average.Tmean; + } + + stars[N].found = TRUE; + i++; + } + + /*************** add unmatched stars *************************/ + + /** incorporate unmatched image stars, if this star is in field of this tgtcat **/ + /* these new entries are all written out in UPDATE mode */ + for (off_t i = 0; (i < Nstars) && !options->only_match; i++) { + /* make sure there is space for next entry */ + if (Nmeas >= NMEAS - stars[i].average.Nmeasure) { + NMEAS = Nmeas + 1000 + stars[i].average.Nmeasure; + REALLOCATE (next_meas, off_t, NMEAS); + REALLOCATE (tgtcat[0].measure, Measure, NMEAS); + } + if (Nave >= NAVE) { + NAVE = Nave + 1000; + REALLOCATE (tgtcat[0].average, Average, NAVE); + if (tgtcat[0].secfilt) { + // we only update the secfilt table if it has been allocated for output + REALLOCATE (tgtcat[0].secfilt, SecFilt, NAVE*tgtcat[0].Nsecfilt); + } + } + + if (stars[i].found) continue; + if (!IN_REGION (stars[i].average.R, stars[i].average.D)) continue; + + dvo_average_init (&tgtcat[0].average[Nave]); + tgtcat[0].average[Nave].R = stars[i].average.R; + tgtcat[0].average[Nave].D = stars[i].average.D; + tgtcat[0].average[Nave].dR = stars[i].average.dR; + tgtcat[0].average[Nave].dD = stars[i].average.dD; + + tgtcat[0].average[Nave].measureOffset = Nmeas; + tgtcat[0].average[Nave].objID = objID; + tgtcat[0].average[Nave].catID = catID; + + if (PSPS_ID) { + tgtcat[0].average[Nave].extID = CreatePSPSObjectID(tgtcat[0].average[Nave].R, tgtcat[0].average[Nave].D); + } + + if (ACCEPT_MOTION) { + tgtcat[0].average[Nave].uR = stars[i].average.uR; + tgtcat[0].average[Nave].uD = stars[i].average.uD; + tgtcat[0].average[Nave].duR = stars[i].average.duR; + tgtcat[0].average[Nave].duD = stars[i].average.duD; + tgtcat[0].average[Nave].P = stars[i].average.P; + tgtcat[0].average[Nave].dP = stars[i].average.dP; + tgtcat[0].average[Nave].Tmean = stars[i].average.Tmean; + } + + objID ++; + + // we only update the secfilt table if it has been allocated for output + for (int j = 0; tgtcat[0].secfilt && (j < Nsecfilt); j++) { + dvo_secfilt_init (&tgtcat[0].secfilt[Nave*Nsecfilt+j], SECFILT_RESET_ALL); + } + + // supply the measurements from this detection + for (int mSrc = 0; mSrc < stars[i].average.Nmeasure; mSrc++) { + dvo_measure_init (&tgtcat[0].measure[Nmeas]); + tgtcat[0].measure[Nmeas] = stars[i].measure[mSrc]; + + // the following measure elements cannot be set until here: + tgtcat[0].measure[Nmeas].dbFlags = 0; + tgtcat[0].measure[Nmeas].averef = Nave; + tgtcat[0].measure[Nmeas].objID = tgtcat[0].average[Nave].objID; + tgtcat[0].measure[Nmeas].catID = tgtcat[0].catID; + + /* set the average magnitude if not already set and the photcode.equiv is not 0 */ + /* in UPDATE mode, this value is not saved; use relphot to recalculate */ + if (Nsec > -1) { + tgtcat[0].secfilt[Nave*Nsecfilt+Nsec].MpsfChp = PhotCat (&tgtcat[0].measure[Nmeas], MAG_CLASS_PSF); + } + + tgtcat[0].average[Nave].Nmeasure ++; + Nmeas ++; + + // update the next_meas pointer for this entry (last one for this star is -1) + next_meas[Nmeas-1] = (mSrc < stars[i].average.Nmeasure - 1) ? Nmeas : -1; + } + stars[i].found = TRUE; + Nave ++; + } + + REALLOCATE (tgtcat[0].average, Average, Nave); + REALLOCATE (tgtcat[0].measure, Measure, Nmeas); + + if (options->nosort) { + tgtcat[0].sorted = FALSE; + } else { + tgtcat[0].sorted = TRUE; + tgtcat[0].measure = sort_measure (tgtcat[0].average, Nave, tgtcat[0].measure, Nmeas, next_meas); + } + + /* check if the tgtcat has changed? if no change, no need to write */ + tgtcat[0].objID = objID; // new max value, save on tgtcat close + tgtcat[0].Naverage = Nave; + tgtcat[0].Nmeasure = Nmeas; + tgtcat[0].Nsecfilt_mem = tgtcat[0].secfilt ? Nave*Nsecfilt : 0; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n", Nstars, Nave, Nmeas, Nmatch); + + free (X1); + free (Y1); + free (N1); + free (X2); + free (Y2); + free (N2); + free (next_meas); + + return (Nmatch); +} + +/* + notes: + + for finding if a tgtcat star is in an image or an image star is in the tgtcat: + + tgtcats have boundaries defined by RA and DEC, but they may curve in projection + images have boundaries which are lines in pixels coords, but curve in RA and DEC + + tgtcat[0].found[Ncat] but stars[Nstar].found + +*/ Index: trunk/Ohana/src/addstar/src/loadgaia_dr2.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgaia_dr2.c (revision 40543) +++ trunk/Ohana/src/addstar/src/loadgaia_dr2.c (revision 40543) @@ -0,0 +1,43 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +/* This is the DVO program to upload gaia_dr2 detections from gaia cvs files into a DVO database. + It is modeled on the loadgaia program but the source files are CSV tables. + + USAGE: loadgaia_dr2 -D CATDIR (catdir) (gaia_dr2file) [...more files] +*/ + +int main (int argc, char **argv) { + + SkyTable *sky; + SkyList *skylist = NULL; + AddstarClientOptions options; + + // need to construct these options with args_loadtycho... + SetSignals (); + options = ConfigInit (&argc, argv); + options = args_loadgaia_dr2 (&argc, argv, options); + + // load the full sky description table: + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // generate the subset matching the user-selected region + skylist = SkyListByPatch (sky, -1, &UserPatch); + + int Nstart = 1; + int Nend = argc; + while (Nstart < Nend) { + Nstart = loadgaia_dr2_table (Nstart, Nend, skylist, NULL, argv, &options); + } + + SkyTableFree (sky); + SkyListFree(skylist); + FreePhotcodeTable(); + free (CATDIR); + + ohana_memcheck (VERBOSE); + ohana_memdump (VERBOSE); + exit (0); +} + Index: trunk/Ohana/src/addstar/src/loadgaia_dr2_catalog.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgaia_dr2_catalog.c (revision 40543) +++ trunk/Ohana/src/addstar/src/loadgaia_dr2_catalog.c (revision 40543) @@ -0,0 +1,32 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +int loadgaia_dr2_catalog (Gaia_DR2_Stars *stars, int Nstars, SkyRegion *region, char *filename, AddstarClientOptions *options) { + + Catalog catalog; + + // now we have all of the loaded stars in this catalog + dvo_catalog_init (&catalog, TRUE); + 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 = DVO_LOAD_AVERAGE | DVO_LOAD_SECFILT | DVO_LOAD_MEASURE; + 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_gaia_dr2 (region, stars, Nstars, &catalog, options); + + SetProtect (TRUE); + if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); } + if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); } + SetProtect (FALSE); + + dvo_catalog_free (&catalog); + + return (TRUE); +} Index: trunk/Ohana/src/addstar/src/loadgaia_dr2_make_subset.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgaia_dr2_make_subset.c (revision 40543) +++ trunk/Ohana/src/addstar/src/loadgaia_dr2_make_subset.c (revision 40543) @@ -0,0 +1,44 @@ +# include "addstar.h" +# include "gaia_dr2.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) + +Gaia_DR2_Stars *loadgaia_dr2_make_subset (Gaia_DR2_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) { + + int i; + + Gaia_DR2_Stars *subset = NULL; + + // collect array of (Stars *) stars in a new output catalog + int Nsubset = 0; + int NSUBSET = 3000; + ALLOCATE (subset, Gaia_DR2_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].average.R < region[0].Rmin) continue; + if (stars[i].average.R >= region[0].Rmax) break; + if (stars[i].average.D < region[0].Dmin) continue; + if (stars[i].average.D >= region[0].Dmax) continue; + + // check if in UserPatch (a GLOBAL) + if (stars[i].average.R < UserPatch.Rmin) continue; + if (stars[i].average.R > UserPatch.Rmax) break; + if (stars[i].average.D < UserPatch.Dmin) continue; + if (stars[i].average.D > UserPatch.Dmax) continue; + + subset[Nsubset] = stars[i]; + Nsubset ++; + + stars[i].flag = TRUE; + + CHECK_REALLOCATE (subset, Gaia_DR2_Stars, NSUBSET, Nsubset, 10000); + } + + *nsubset = Nsubset; + return subset; +} Index: trunk/Ohana/src/addstar/src/loadgaia_dr2_readstars.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgaia_dr2_readstars.c (revision 40543) +++ trunk/Ohana/src/addstar/src/loadgaia_dr2_readstars.c (revision 40543) @@ -0,0 +1,436 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +/* this function reads the values of interest from the Gaia CSV files: + + data model description: + http://gea.esac.esa.int/archive/documentation/GDR2/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html + + 3 source_id (map to extID) +- 5 ref_epoch + 6 ra + 7 ra_error + 8 dec + 9 dec_error + 10 parallax + 11 parallax_error +-12 parallax_over_error + 13 pmra + 14 pmra_error + 15 pmdec + 16 pmdec_error + + 27 astrometric_n_obs_al + 29 astrometric_n_good_obs_al + 31 astrometric_gof_al + 32 astrometric_chi2_al + 33 astrometric_excess_noise + 34 astrometric_excess_noise_sig + + 42 visibility_periods_used + 46 duplicated_source + + 50 phot_g_mean_flux_over_error + 51 phot_g_mean_mag + 55 phot_bp_mean_flux_over_error + 56 phot_bp_mean_mag + 60 phot_rp_mean_flux_over_error + 61 phot_rp_mean_mag + + 63 phot_proc_mode : use for a flag + +// - sqrt(astrometric_chi2_al / (astrometric_n_good_obs_al - 5)) < 1.2 * max(1.0, exp(-0.2 * (phot_g_mean_mag - 19.5))) + + we have 32 x 2 bits of photFlags potentially available for use. Gaia does have have a single flag + value, so I will generate bits to represent some data quality cuts proposed in the Gaia docs. + + I'm also going to overload a couple of values as needed: + + measure.psfNpix : astrometric_n_obs_al (nObs) + measure.psfNdof : astrometric_n_good_obs_al (nGood) + measure.extNsigma : astrometric_gof_al (goFit) + measure.psfChisq : astrometric_chi2_al (chiSq) + measure.psfQF : astrometric_excess_noise (exNoise) + measure.psfQFPerf : astrometric_excess_noise_sig (exNoiseSig) + measure.FWx : visibility_periods_used (nPeriods) + + photFlags : + 0x0000.0001 : visibility_periods_used > 6 (nPeriods > 6) + 0x0000.0002 : visibility_periods_used > 8 (nPeriods > 8) + 0x0000.0004 : astrometric_excess_noise < 1.0 (exNoise < 1) + 0x0000.0008 : ChiSq < limit ( + + 0x0000.0010 : duplicated_source + 0x0000.0020 : phot_proc_mode == 0 + 0x0000.0040 : phot_proc_mode == 1 + 0x0000.0080 : phot_proc_mode == 2 + +*/ + +Gaia_DR2_Stars *loadgaia_dr2_readstars (char *filename, Gaia_DR2_Stars *stars, int *nstars, AddstarClientOptions *options) { + + int codeG = GetPhotcodeCodebyName ("GAIA_G_DR2"); if (!codeG) Shutdown ("missing photocde GAIA_G_DR2"); + int codeB = GetPhotcodeCodebyName ("GAIA_B_DR2"); if (!codeB) Shutdown ("missing photocde GAIA_B_DR2"); + int codeR = GetPhotcodeCodebyName ("GAIA_R_DR2"); if (!codeR) Shutdown ("missing photocde GAIA_R_DR2"); + + // GAIA DR2 Epoch is 2015.5 == 2015/07/02,15:00:00 + time_t GAIA_DR2_EPOCH = ohana_date_to_sec ("2015/07/02,15:00:00"); + + // read in the full FITS files + FILE *f = fopen (filename, "r"); + if (f == NULL) Shutdown ("can't read gaia_dr2 file: %s", filename); + + int Nskip = 1; // Gaia CSV files have a single header row (and now special character to mark) + + int Nelem = 0; // number of valid rows read (vector elements) + int NELEM = 10000; // currently-allocated number of output rows + + // for initial testing, a minimal subset: + ALLOCATE_PTR (srcID, uint64_t, NELEM); + + ALLOCATE_PTR (Rg, double, NELEM); + ALLOCATE_PTR (dRg, double, NELEM); + ALLOCATE_PTR (Dg, double, NELEM); + ALLOCATE_PTR (dDg, double, NELEM); + + ALLOCATE_PTR (plx, double, NELEM); + ALLOCATE_PTR (dplx, double, NELEM); + + ALLOCATE_PTR (uR, double, NELEM); + ALLOCATE_PTR (duR, double, NELEM); + + ALLOCATE_PTR (uD, double, NELEM); + ALLOCATE_PTR (duD, double, NELEM); + + ALLOCATE_PTR (nObs, int, NELEM); + ALLOCATE_PTR (nGood, int, NELEM); + + ALLOCATE_PTR (goFit, double, NELEM); + ALLOCATE_PTR (chiSq, double, NELEM); + + ALLOCATE_PTR (exNoise, double, NELEM); + ALLOCATE_PTR (exNoiseSig, double, NELEM); + + ALLOCATE_PTR (nPeriods, int, NELEM); + ALLOCATE_PTR (dupSource, int, NELEM); + + ALLOCATE_PTR (gg, double, NELEM); + ALLOCATE_PTR (dgg, double, NELEM); + + ALLOCATE_PTR (bg, double, NELEM); + ALLOCATE_PTR (dbg, double, NELEM); + + ALLOCATE_PTR (rg, double, NELEM); + ALLOCATE_PTR (drg, double, NELEM); + + ALLOCATE_PTR (procMode, int, NELEM); + + // we allocate one extra byte into which we never read so there will always be a NULL terminating the string + ALLOCATE_PTR (buffer, char, 0x10001); + bzero (buffer, 0x10001); + + int Nline_read = 0; // track number of lines read so far (use to skip lines as well) + + // we have a working buffer read from the file. we parse the lines in the working buffer + // until we reach the last chunk without an EOL char. at that point, we shift the start + // of the last (partial) line to the start of the buffer and re-fill. + + // we treat \n\r pair as a single EOL char to handle mac files: + + int Nstart = 0; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + while (!EndOfFile) { + int Nbytes = 0x10000 - Nstart; + // we have allocated one extra byte into which we never read so there will always be a NULL terminating the string + bzero (&buffer[Nstart], Nbytes + 1); + int Nread = fread (&buffer[Nstart], 1, Nbytes, f); + if (ferror (f)) { + perror ("error reading data file"); + break; + } + // we still need to parse the rest of the buffer, but there might not be an EOL on the last line + if (Nread == 0) { + EndOfFile = TRUE; + } + + int bufferStatus = TRUE; + char *c0 = buffer; // c0 always marks the start of a line + char *cA = NULL; // cA will carry the curr point within the line + while (bufferStatus) { + + // find the end of this current line (\n or \r). if we hit the end of the buffer (NULL), + // attempt to read more data. finish up when we hit the end of the file + char *c1 = strchr (c0, '\n'); // find the end of this current line (also valid for a Mac: \r\n) + if (!c1) { + c1 = strchr (c0, '\r'); // try \r for Windows files + } + if (!c1) { + Nstart = strlen (c0); + if (EndOfFile) { + // if we have reached EOF, we need to do one last pass in case there is a line without a return + c1 = c0 + Nstart; + bufferStatus = FALSE; + if (Nstart == 0) continue; // if we have reached EOF and c0 points at the last valid character, we are done + } else { + // if we have not reached EOF, we need to shift the buffer to the start of this line and read more data + memmove (buffer, c0, Nstart); + bufferStatus = FALSE; + continue; + } + } + *c1 = 0; // mark the end of the line + Nline_read ++; + + // skip to the next line (but if EOF, do not overrun buffer) + if (Nline_read <= Nskip) { if (!EndOfFile) { c0 = c1 + 1; } continue; } + + // these are not needed: Gaia CSV files do not have any commented-out lines + if (*c0 == '#') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + if (*c0 == '!') { if (!EndOfFile) { c0 = c1 + 1; } continue; } + + // XXX for Gaia, we know which columns we want in advance + + int lineStatus = TRUE; + int readStatus; + double dvalue; + int ivalue; + uint64_t jvalue; + + // cA will follow the currently extracted field, c0 points to the start of the line + cA = c0; + + // Tref == 2015.5 for DR2 + + // the start of the line is the 1st element (fields are 1-counting) + cA = jparse_csv_rpt (&jvalue, 3, 1, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 3, "source_id", c0); } lineStatus &= readStatus; srcID[Nelem] = jvalue; + cA = dparse_csv_rpt (&dvalue, 6, 3, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 6, "Ra", c0); } lineStatus &= readStatus; Rg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 7, 6, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 7, "dRa", c0); } lineStatus &= readStatus; dRg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 8, 7, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 8, "Dec", c0); } lineStatus &= readStatus; Dg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 9, 8, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 9, "dDec", c0); } lineStatus &= readStatus; dDg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 10, 9, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 10, "plx", c0); } lineStatus &= readStatus; plx[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 11, 10, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 11, "dplx", c0); } lineStatus &= readStatus; dplx[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 13, 11, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 13, "uR", c0); } lineStatus &= readStatus; uR[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 14, 13, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 14, "duR", c0); } lineStatus &= readStatus; duR[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 15, 14, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 15, "uD", c0); } lineStatus &= readStatus; uD[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 16, 15, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 16, "duD", c0); } lineStatus &= readStatus; duD[Nelem] = dvalue; + cA = iparse_csv_rpt (&ivalue, 27, 16, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 27, "n_obs", c0); } lineStatus &= readStatus; nObs[Nelem] = ivalue; + cA = iparse_csv_rpt (&ivalue, 29, 27, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 29, "n_good", c0); } lineStatus &= readStatus; nGood[Nelem] = ivalue; + cA = dparse_csv_rpt (&dvalue, 31, 29, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 31, "GOF", c0); } lineStatus &= readStatus; goFit[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 32, 31, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 32, "chiSq", c0); } lineStatus &= readStatus; chiSq[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 33, 32, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 33, "exNoise", c0); } lineStatus &= readStatus; exNoise[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 34, 33, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 34, "exNoiseSig", c0); } lineStatus &= readStatus; exNoiseSig[Nelem]= dvalue; + cA = iparse_csv_rpt (&ivalue, 42, 34, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 42, "visibility_periods", c0); } lineStatus &= readStatus; nPeriods[Nelem] = ivalue; + + // the duplicated_source field has the value false (or true) + for (int i = 42; i < 46; i++) cA = parse_nextword_csv (cA); + dupSource[Nelem] = strncmp (cA, "false", 5) ? 1 : 0; +// cA = dparse_csv_rpt (&dvalue, 46, 42, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 46, "phot_r", c0); } lineStatus &= readStatus; dupSource[Nelem] = dvalue; + + cA = dparse_csv_rpt (&dvalue, 50, 46, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 50, "dphot_g", c0); } lineStatus &= readStatus; dgg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 51, 50, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 51, "phot_g", c0); } lineStatus &= readStatus; gg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 55, 51, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 55, "dphot_b", c0); } lineStatus &= readStatus; dbg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 56, 55, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 56, "phot_b", c0); } lineStatus &= readStatus; bg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 60, 56, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 60, "dphot_r", c0); } lineStatus &= readStatus; drg[Nelem] = dvalue; + cA = dparse_csv_rpt (&dvalue, 61, 60, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 61, "phot_r", c0); } lineStatus &= readStatus; rg[Nelem] = dvalue; + cA = iparse_csv_rpt (&ivalue, 63, 61, cA, &readStatus); if (!readStatus && VERBOSE) { gprint (GP_ERR, "suspect field: %d (%s) in %s\n", 63, "phot_mode", c0); } lineStatus &= readStatus; procMode[Nelem] = ivalue; + + if (!lineStatus && VERBOSE) { + // why do I need to copy temp here, does gprint modify the value of temp? + char temp[32]; + strncpy (temp, c0, 32); + temp[31] = 0; + gprint (GP_ERR, "skip line %s\n\n", temp); + } + + Nelem ++; + if (Nelem == NELEM) { + NELEM += 1000; + REALLOCATE (srcID, uint64_t, NELEM); + REALLOCATE (Rg, double, NELEM); + REALLOCATE (dRg, double, NELEM); + REALLOCATE (Dg, double, NELEM); + REALLOCATE (dDg, double, NELEM); + REALLOCATE (plx, double, NELEM); + REALLOCATE (dplx, double, NELEM); + REALLOCATE (uR, double, NELEM); + REALLOCATE (duR, double, NELEM); + REALLOCATE (uD, double, NELEM); + REALLOCATE (duD, double, NELEM); + REALLOCATE (nObs, int, NELEM); + REALLOCATE (nGood, int, NELEM); + REALLOCATE (goFit, double, NELEM); + REALLOCATE (chiSq, double, NELEM); + REALLOCATE (exNoise, double, NELEM); + REALLOCATE (exNoiseSig, double, NELEM); + REALLOCATE (nPeriods, int, NELEM); + REALLOCATE (dupSource, int, NELEM); + REALLOCATE (gg, double, NELEM); + REALLOCATE (dgg, double, NELEM); + REALLOCATE (bg, double, NELEM); + REALLOCATE (dbg, double, NELEM); + REALLOCATE (rg, double, NELEM); + REALLOCATE (drg, double, NELEM); + REALLOCATE (procMode, int, NELEM); + } + if (!EndOfFile) { + c0 = c1 + 1; + } + } + } + + // Nelem is now the number of items (objects,stars) read from the Gaia CSV file + int NstarsIn = Nelem; + + double Rmin = +360.0; + double Rmax = -360.0; + double Dmin = +360.0; + double Dmax = -360.0; + + // start off where we finished on a previous read + int Nstars = *nstars; + int NSTARS = Nstars + 0.1*NstarsIn; + + if (!stars) { + ALLOCATE (stars, Gaia_DR2_Stars, NSTARS); + } else { + REALLOCATE (stars, Gaia_DR2_Stars, NSTARS); + } + + for (int i = 0; i < NstarsIn; i++) { + + Rmin = MIN (Rmin, Rg[i]); + Rmax = MAX (Rmax, Rg[i]); + Dmin = MIN (Dmin, Dg[i]); + Dmax = MAX (Dmax, Dg[i]); + + // each Gaia DR2 source corresponds to 3 measurements: Gg, Gr, Gb: + + dvo_average_init (&stars[Nstars].average); + for (int j = 0; j < 3; j++) { + dvo_measure_init (&stars[Nstars].measure[j]); + } + + stars[Nstars].average.R = Rg[i]; + stars[Nstars].average.D = Dg[i]; + stars[Nstars].average.dR = 0.001*dRg[i]; + stars[Nstars].average.dD = 0.001*dDg[i]; + + if (ACCEPT_MOTION) { + // GAIA values are reported in milliarcseconds, but DVO wants arcseconds + stars[Nstars].average.uR = uR[i] * 0.001; + stars[Nstars].average.uD = uD[i] * 0.001; + stars[Nstars].average.duR = duR[i] * 0.001; + stars[Nstars].average.duD = duD[i] * 0.001; + stars[Nstars].average.P = plx[i] * 0.001; + stars[Nstars].average.dP = dplx[i] * 0.001; + stars[Nstars].average.Tmean = GAIA_DR2_EPOCH; + } + + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + // for the following cut, see https://arxiv.org/pdf/1804.09378.pdf (Gaia Collaboration et al 2018) + double ChiSqLimit = (nGood[i] > 5) ? sqrt(chiSq[i] / (nGood[i] - 5)) : NAN; + double ChiCompare = 1.2 * MAX (1.0, exp(-0.2*(gg[i] - 19.5))); + + int NVis6 = (nPeriods[i] > 6) ? 0x00000001 : 0x0; + int NVis8 = (nPeriods[i] > 8) ? 0x00000002 : 0x0; + int ExNoise = (exNoise[i] < 1.0) ? 0x00000004 : 0x0; + int ChiGood = (ChiSqLimit < ChiCompare) ? 0x00000008 : 0x0; + + int isDup = dupSource[i] ? 0x00000010 : 0x0; + int Mode0 = (procMode[i] == 0) ? 0x00000020 : 0x0; + int Mode1 = (procMode[i] == 1) ? 0x00000040 : 0x0; + int Mode2 = (procMode[i] == 2) ? 0x00000080 : 0x0; + + int photFlags = NVis6 | NVis8 | ExNoise | ChiGood | isDup | Mode0 | Mode1 | Mode2; + + for (int j = 0; j < 3; j++) { + stars[Nstars].measure[j].R = Rg[i]; + stars[Nstars].measure[j].D = Dg[i]; + + stars[Nstars].measure[j].psfNpix = nObs[i]; + stars[Nstars].measure[j].psfNdof = nGood[i]; + stars[Nstars].measure[j].extNsigma = goFit[i]; + stars[Nstars].measure[j].psfChisq = chiSq[i]; + stars[Nstars].measure[j].psfQF = exNoise[i]; + stars[Nstars].measure[j].psfQFperf = exNoiseSig[i]; + stars[Nstars].measure[j].FWx = nPeriods[i]; // this one is not great + + stars[Nstars].measure[j].photFlags = photFlags; + stars[Nstars].measure[j].t = GAIA_DR2_EPOCH; + } + + // the flux reported by Gaia is in electrons / second; we just want to use the mags: + + // dgg, dbg, drg store (flux / error), i.e., S/N. mag error is 2.5/ln(10) ~ 1.0857 / (S/N) + stars[Nstars].measure[0].M = gg[i]; + stars[Nstars].measure[0].dM = 1.08573620476 / dgg[i]; + + stars[Nstars].measure[1].M = bg[i]; + stars[Nstars].measure[1].dM = 1.08573620476 / dbg[i]; + + stars[Nstars].measure[2].M = rg[i]; + stars[Nstars].measure[2].dM = 1.08573620476 / drg[i]; + + stars[Nstars].measure[0].photcode = codeG; + stars[Nstars].measure[1].photcode = codeB; + stars[Nstars].measure[2].photcode = codeR; + + // stars[Nstars].measure.dM = dmgFlux[i] / flux; + // stars[Nstars].measure.FluxPSF = flux; + // stars[Nstars].measure.dFluxPSF = dgFlux[i]; + // stars[Nstars].measure.t = options->timeref; + + stars[Nstars].average.Nmeasure = 3; + Nstars ++; + + CHECK_REALLOCATE (stars, Gaia_DR2_Stars, NSTARS, Nstars, 10000); + } + + free (srcID); + free (Rg); + free (dRg); + free (Dg); + free (dDg); + free (plx); + free (dplx); + free (uR); + free (duR); + free (uD); + free (duD); + free (nObs); + free (nGood); + free (goFit); + free (chiSq); + free (exNoise); + free (exNoiseSig); + free (nPeriods); + free (dupSource); + free (gg); + free (dgg); + free (bg); + free (dbg); + free (rg); + free (drg); + free (procMode); + + free (buffer); + + *nstars = Nstars; + return (stars); +} + +int loadgaia_dr2_sortStars (Gaia_DR2_Stars *stars, int Nstars) { + +# define SWAPFUNC(A,B){ Gaia_DR2_Stars temp = stars[A]; stars[A] = stars[B]; stars[B] = temp; } +# define COMPARE(A,B)(stars[A].average.R < stars[B].average.R) + + OHANA_SORT (Nstars, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + + return TRUE; +} + Index: trunk/Ohana/src/addstar/src/loadgaia_dr2_table.c =================================================================== --- trunk/Ohana/src/addstar/src/loadgaia_dr2_table.c (revision 40543) +++ trunk/Ohana/src/addstar/src/loadgaia_dr2_table.c (revision 40543) @@ -0,0 +1,64 @@ +# include "addstar.h" +# include "gaia_dr2.h" +# define NSTARS_MAX 1000000 + +int loadgaia_dr2_table (int Nstart, int Nend, SkyList *skylistInput, HostTable *hosts, char **filename, AddstarClientOptions *options) { + OHANA_UNUSED_PARAM(hosts); + + int i; + + int Nstars = 0; + Gaia_DR2_Stars *stars = NULL; + for (i = Nstart; (Nstars < NSTARS_MAX) && (i < Nend); i++) { + // read the next file and append to the current set of stars + fprintf (stderr, "loading %s\n", filename[i]); + stars = loadgaia_dr2_readstars (filename[i], stars, &Nstars, options); + } + Nstart = i; // we pass back the entry for the next file to be read + if (!stars) return Nstart; + + fprintf (stderr, "writing %d stars to dvo\n", Nstars); + + // sort the stars by RA + loadgaia_dr2_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].average.R < UserPatch.Rmin) continue; + if (stars[i].average.R > UserPatch.Rmax) break; + if (stars[i].average.D < UserPatch.Dmin) continue; + if (stars[i].average.D > UserPatch.Dmax) continue; + + // identify the relevant catalog + SkyList *skylist = SkyRegionByPoint_List (skylistInput, -1, stars[i].average.R, stars[i].average.D); + if (skylist[0].Nregions == 0) { + SkyListFree (skylist); + continue; + } + SkyRegion *region = skylist[0].regions[0]; + + // select stars matching this region + int Nsubset; + Gaia_DR2_Stars *subset = loadgaia_dr2_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 + // NOTE: disable parallel mode for now: + // loadgaia_dr2_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options); + loadgaia_dr2_catalog (subset, Nsubset, region, skylist[0].filename[0], options); + free (subset); + SkyListFree (skylist); + } + + // wait for last remote clients to finish + // NOTE: disable parallel mode for now: + // harvest_all (); + + free (stars); + return Nstart; +} + Index: trunk/Ohana/src/addstar/src/parse_csv.c =================================================================== --- trunk/Ohana/src/addstar/src/parse_csv.c (revision 40543) +++ trunk/Ohana/src/addstar/src/parse_csv.c (revision 40543) @@ -0,0 +1,116 @@ +# include "addstar.h" +# include "gaia_dr2.h" + +// I want a function which takes a pointer to the start of a line +// and an entry number N and returns the Nth element, along with the pointer +// to the following (or to th + +char *dparse_csv_rpt (double *X, int Nwant, int Nlast, char *line, int *status) { + + int i; + char *word; + char *ptr; + + // only scan forward for NX from this point + int NX = Nwant - Nlast + 1; + + // scan for the NX word + word = line; + for (i = 0; i < NX - 1; i++) { + word = parse_nextword_csv (word); + } + + if (word[0] == '"') word[0] = ' '; + if (word[0] == ',') { + *X = NAN; + *status = 1; + return word; + } + + // return FALSE if this field is not consistent with a double + *X = strtod (word, &ptr); + + // default result is success + *status = 1; + + // if the leading value is '-', return -1 so we can handle -00 + if (word[0] == '-') { *status = -1; } + + // if ptr = word, then there is not actually a valid number here + if (ptr == word) { *status = FALSE; *X = NAN; } + + return word; +} + +char *iparse_csv_rpt (int *X, int Nwant, int Nlast, char *line, int *status) { + + int i; + char *word; + char *ptr; + + // only scan forward for NX from this point + int NX = Nwant - Nlast + 1; + + // scan for the NX word + word = line; + for (i = 0; i < NX - 1; i++) { + word = parse_nextword_csv (word); + } + + if (word[0] == '"') word[0] = ' '; + if (word[0] == ',') { + *X = 0; + *status = 1; + return word; + } + + *X = strtol (word, &ptr, 0); + + // default result is success + *status = 1; + + // if the leading value is '-', return -1 so we can handle -00 + if (word[0] == '-') { *status = -1; } + + // if ptr = word, then there is not actually a valid number here + if (ptr == word) { *status = FALSE; *X = 0; } + + return word; +} + +char *jparse_csv_rpt (uint64_t *X, int Nwant, int Nlast, char *line, int *status) { + + int i; + char *word; + char *ptr; + + // only scan forward for NX from this point + int NX = Nwant - Nlast + 1; + + // scan for the NX word + word = line; + for (i = 0; i < NX - 1; i++) { + word = parse_nextword_csv (word); + } + + if (word[0] == '"') word[0] = ' '; + if (word[0] == ',') { + *X = 0; + *status = 1; + return word; + } + + *X = strtoll (word, &ptr, 0); + + // default result is success + *status = 1; + + // if the leading value is '-', return -1 so we can handle -00 + if (word[0] == '-') { *status = -1; } + + // if ptr = word, then there is not actually a valid number here + if (ptr == word) { *status = FALSE; *X = 0; } + + return word; +} +