Index: trunk/Ohana/src/addstar/Makefile =================================================================== --- trunk/Ohana/src/addstar/Makefile (revision 41145) +++ trunk/Ohana/src/addstar/Makefile (revision 41151) @@ -27,4 +27,5 @@ loadgaia : $(BIN)/loadgaia.$(ARCH) loadgaia_dr2 : $(BIN)/loadgaia_dr2.$(ARCH) +loadukirt_uhs : $(BIN)/loadukirt_uhs.$(ARCH) loadstarpar : $(BIN)/loadstarpar.$(ARCH) loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH) @@ -44,5 +45,5 @@ # programs in 'SERVER' use the client-server concept and are out of date -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 +INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadgaia_dr2 loadukirt_uhs loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos SERVER = addstarc addstard addstart @@ -305,4 +306,26 @@ $(SRC)/psps_ids.$(ARCH).o +LOAD-UKIRT_UHS = \ +$(SRC)/loadukirt_uhs.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args_loadukirt_uhs.$(ARCH).o \ +$(SRC)/find_matches_ukirt_uhs.$(ARCH).o \ +$(SRC)/loadukirt_uhs_catalog.$(ARCH).o \ +$(SRC)/loadukirt_uhs_make_subset.$(ARCH).o \ +$(SRC)/loadukirt_uhs_readstars.$(ARCH).o \ +$(SRC)/loadukirt_uhs_readstars_uhs.$(ARCH).o \ +$(SRC)/loadukirt_uhs_readstars_ugcs.$(ARCH).o \ +$(SRC)/loadukirt_uhs_readstars_ugps.$(ARCH).o \ +$(SRC)/loadukirt_uhs_readstars_ulas.$(ARCH).o \ +$(SRC)/loadukirt_uhs_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 \ @@ -528,4 +551,5 @@ $(LOAD-GAIA) : $(INC)/addstar.h $(INC)/gaia.h $(LOAD-GAIA_DR2) : $(INC)/addstar.h $(INC)/gaia_dr2.h +$(LOAD-UKIRT_UHS) : $(INC)/addstar.h $(INC)/ukirt_uhs.h $(LOAD-STARPAR) : $(INC)/addstar.h $(INC)/loadstarpar.h $(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h @@ -551,4 +575,5 @@ $(BIN)/loadgaia.$(ARCH) : $(LOAD-GAIA) $(BIN)/loadgaia_dr2.$(ARCH) : $(LOAD-GAIA_DR2) +$(BIN)/loadukirt_uhs.$(ARCH) : $(LOAD-UKIRT_UHS) $(BIN)/loadstarpar.$(ARCH) : $(LOAD-STARPAR) $(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT) Index: trunk/Ohana/src/addstar/include/ukirt_uhs.h =================================================================== --- trunk/Ohana/src/addstar/include/ukirt_uhs.h (revision 41151) +++ trunk/Ohana/src/addstar/include/ukirt_uhs.h (revision 41151) @@ -0,0 +1,65 @@ + +typedef enum {UKIRT_MODE_NONE, UKIRT_MODE_UHS, UKIRT_MODE_UGCS, UKIRT_MODE_UGPS, UKIRT_MODE_ULAS} UkirtMode; + +# define NSTARS_MAX 10000000 +// # define NSTARS_MAX 10 + +# define BUFFER_SIZE 0x100000 +// # define BUFFER_SIZE 900 + +typedef struct { + Average average; + Measure *measure; + int flag; // in a subset? + int found; // assigned to an object? +} UKIRT_Stars; + +int HOST_ID; +char *HOSTDIR; +char *CPT_FILE; +char *INPUT; +int UKIRT_NFILTER; + +UkirtMode UKIRT_MODE; + +AddstarClientOptions args_loadukirt_uhs (int *argc, char **argv, AddstarClientOptions options); +// AddstarClientOptions args_loadukirt_uhs_client (int *argc, char **argv, AddstarClientOptions options); + +int loadukirt_uhs_table (SkyList *skylistInput, char *filename, AddstarClientOptions *options); + +UKIRT_Stars *loadukirt_uhs_make_subset (UKIRT_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset); + +// int loadukirt_uhs_save_remote (UKIRT_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 loadukirt_uhs_catalog (UKIRT_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_ukirt_uhs (SkyRegion *region, UKIRT_Stars *stars, int Nstars, Catalog *catalog, AddstarClientOptions *options); + +int loadukirt_uhs_save_stars (char *filename, UKIRT_Stars *stars, int Nstars); +// UKIRT_Stars *loadukirt_uhs_load_stars (char *filename, int *nstars); + +UKIRT_Stars *loadukirt_uhs_readstars (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars); + +UKIRT_Stars *loadukirt_uhs_readstars_ugps (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars); +UKIRT_Stars *loadukirt_uhs_readstars_ugcs (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars); +UKIRT_Stars *loadukirt_uhs_readstars_ulas (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars); +UKIRT_Stars *loadukirt_uhs_readstars_uhs (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars); + + +int loadukirt_uhs_sortStars (UKIRT_Stars *stars, int Nstars); + +// int loadukirt_uhs_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/ReadImageHeader.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 41145) +++ trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 41151) @@ -118,5 +118,5 @@ if (photcode == 0) { if (!gfits_scan (header, "PHOTCODE", "%s", 1, photname)) { - fprintf (stderr, "photcode not supplied in header\n"); + fprintf (stderr, "failure: photcode not supplied in header\n"); return (FALSE); } Index: trunk/Ohana/src/addstar/src/args_loadukirt_uhs.c =================================================================== --- trunk/Ohana/src/addstar/src/args_loadukirt_uhs.c (revision 41151) +++ trunk/Ohana/src/addstar/src/args_loadukirt_uhs.c (revision 41151) @@ -0,0 +1,165 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +static void help (void); + +AddstarClientOptions args_loadukirt_uhs (int *argc, char **argv, AddstarClientOptions options) { + + int N; + + /* check for help request */ + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + help (); + } + + // UKIRT_MODE is global + UKIRT_MODE = UKIRT_MODE_NONE; + if ((N = get_argument (*argc, argv, "-uhs"))) { + UKIRT_MODE = UKIRT_MODE_UHS; + UKIRT_NFILTER = 1; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-ugcs"))) { + UKIRT_MODE = UKIRT_MODE_UGCS; + UKIRT_NFILTER = 6; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-ugps"))) { + UKIRT_MODE = UKIRT_MODE_UGPS; + UKIRT_NFILTER = 5; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-ulas"))) { + UKIRT_MODE = UKIRT_MODE_ULAS; + UKIRT_NFILTER = 5; + remove_argument (N, argc, argv); + } + if (UKIRT_MODE == UKIRT_MODE_NONE) { + fprintf (stderr, "missing or invalid ukirt mode\n"); + 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 loadukirt_uhs */ + // 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: loadukirt_uhs [options] (fitsfile)\n"); + exit (2); + } + return (options); +} + +static void help () { + + fprintf (stderr, "USAGE: loadukirt_uhs [options] (file) [-uhs | -ugcs | -ulas | -ugps] [..more files]\n"); + fprintf (stderr, " add data from UKIRT CSV file to DVO\n"); + fprintf (stderr, " mode is required : -uhs | -ugcs | -ulas | -ugps\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_ukirt_uhs.c =================================================================== --- trunk/Ohana/src/addstar/src/find_matches_ukirt_uhs.c (revision 41151) +++ trunk/Ohana/src/addstar/src/find_matches_ukirt_uhs.c (revision 41151) @@ -0,0 +1,346 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +int find_matches_ukirt_uhs (SkyRegion *region, UKIRT_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/loadukirt_uhs.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs.c (revision 41151) @@ -0,0 +1,40 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +/* This is the DVO program to upload UKIRT_UHS detections from cvs files supplied by Mike + Read into a DVO database. It is modeled on the loadgaia_dr2 program with column modifications + determined by Zhoujian Zhang + + USAGE: loadukirt_uhs -D CATDIR (catdir) (ukirt_uhsfile) [...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_loadukirt_uhs (&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); + + loadukirt_uhs_table (skylist, argv[1], &options); + + SkyTableFree (sky); + SkyListFree(skylist); + FreePhotcodeTable(); + free (CATDIR); + + ohana_memcheck (VERBOSE); + ohana_memdump (VERBOSE); + exit (0); +} + Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_catalog.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_catalog.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_catalog.c (revision 41151) @@ -0,0 +1,32 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +int loadukirt_uhs_catalog (UKIRT_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_ukirt_uhs (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/loadukirt_uhs_make_subset.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_make_subset.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_make_subset.c (revision 41151) @@ -0,0 +1,58 @@ +# include "addstar.h" +# include "ukirt_uhs.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) + +UKIRT_Stars *loadukirt_uhs_make_subset (UKIRT_Stars *stars, int Nstars, int start, SkyRegion *region, int *nsubset) { + + int i; + + // collect array of (Stars *) stars in a new output catalog + int Nsubset = 0; + int NSUBSET = 3000; + + ALLOCATE_PTR (subset, UKIRT_Stars, NSUBSET); + // for (int i = 0; i < NSUBSET; i++) { + // ALLOCATE (subset[i].measure, Measure, UKIRT_NFILTER); + // } + + // 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; + + // the *measure value is a copy of the version allocated to stars[i] + // do NOT free subset[i].measure + subset[Nsubset] = stars[i]; + Nsubset ++; + + stars[i].flag = TRUE; + + CHECK_REALLOCATE (subset, UKIRT_Stars, NSUBSET, Nsubset, 10000); + + /* + if (NSTARS >= Nstars) { + NSTARS += 10000; + REALLOCATE (stars, UKIRT_Stars, NSTARS); + for (int j = Nstars; j < NSTARS; j++) { + ALLOCATE (stars[j].measure, Measure, NFILTER); + } + } + */ + } + + *nsubset = Nsubset; + return subset; +} Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars.c (revision 41151) @@ -0,0 +1,46 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +UKIRT_Stars *loadukirt_uhs_readstars (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars) { + + switch (UKIRT_MODE) { + case UKIRT_MODE_UHS: { + UKIRT_Stars *stars = loadukirt_uhs_readstars_uhs (f, buffer, nstart, options, nstars); + return stars; + } + + case UKIRT_MODE_UGCS: { + UKIRT_Stars *stars = loadukirt_uhs_readstars_ugcs (f, buffer, nstart, options, nstars); + return stars; + } + + case UKIRT_MODE_UGPS: { + UKIRT_Stars *stars = loadukirt_uhs_readstars_ugps (f, buffer, nstart, options, nstars); + return stars; + } + + case UKIRT_MODE_ULAS: { + UKIRT_Stars *stars = loadukirt_uhs_readstars_ulas (f, buffer, nstart, options, nstars); + return stars; + } + + default: + fprintf (stderr, "programming error: invalid mode\n"); + exit (1); + } + return NULL; +} + +int loadukirt_uhs_sortStars (UKIRT_Stars *stars, int Nstars) { + +# define SWAPFUNC(A,B){ UKIRT_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/loadukirt_uhs_readstars_ugcs.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ugcs.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ugcs.c (revision 41151) @@ -0,0 +1,802 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +/* this function reads the values of interest from the UKIRT UGCS CSV files: + + data model description: doc/mapping_UKIDSS_unWISE_DVO.pdf + + N -- field (column) number in CSV file (1 counting) + | UKIRT field name : DVO measure field + 1 sourceID : extID + 2 cuEventID : + 3 frameSetID : imageID (assumes we do not try to go backwards for these photcodes) + 4 ra : R + 5 dec : D + + 6 sigRa : + 7 sigDec : + 8 epoch : t ???? -- how is there a single epoch? + 9 muRa : + 10 muDec : + 11 sigMuRa : + 12 sigMuDec : + 13 chi2 : + 14 nFrames : + + 15 cx : + 16 cy : + 17 cz : + 18 htmID : + 19 l : + 20 b : + 21 lambda : + 22 eta : + 23 priOrSec : primary bit in dbFlags + 24 zmyPnt : + 25 zmyPntErr : + 26 ymjPnt : + 27 ymjPntErr : + 28 jmhPnt : + 29 jmhPntErr : + 30 hmk_1Pnt : + 31 hmk_1PntErr : + 32 zmyExt : + 33 zmyExtErr : + 34 ymjExt : + 35 ymjExtErr : + 36 jmhExt : + 37 jmhExtErr : + 38 hmk_1Ext : + 39 hmk_1ExtErr : + + 40 mergedClassStat : + 41 mergedClass : psfQFperf or photFlags2 + 42 pStar : psfChisq + 43 pGalaxy : extNsigma + 44 pNoise : + 45 pSaturated : + + 46 zHallMag : M + 47 zHallMagErr : dM + 48 zPetroMag : + 49 zPetroMagErr : + 50 zAperMag3 : Map + 51 zAperMag3Err : dMap + 52 zAperMag4 : Mkron + 53 zAperMag4Err : dMron + 54 zAperMag6 : + 55 zAperMag6Err : + 56 zGausig : + 57 zEll : + 58 zPA : posangle + 59 zErrBits : + 60 zDeblend : + 61 zClass : psfQF + 62 zClassStat : + 63 zppErrBits : photFlags + 64 zSeqNum : detID + 65 zObjID : + 66 zXi : + 67 zEta : + + 68 yHallMag : M + 69 yHallMagErr : dM + 70 yPetroMag : + 71 yPetroMagErr : + 72 yAperMag3 : Map + 73 yAperMag3Err : dMap + 74 yAperMag4 : Mkron + 75 yAperMag4Err : dMron + 76 yAperMag6 : + 77 yAperMag6Err : + 78 yGausig : + 79 yEll : + 80 yPA : posangle + 81 yErrBits : + 82 yDeblend : + 83 yClass : psfQF + 84 yClassStat : + 85 yppErrBits : photFlags + 86 ySeqNum : detID + 87 yObjID : + 88 yXi : + 89 yEta : + + 90 jHallMag : M + 91 jHallMagErr : dM + 92 jPetroMag : + 93 jPetroMagErr : + 94 jAperMag3 : Map + 95 jAperMag3Err : dMap + 96 jAperMag4 : Mkron + 97 jAperMag4Err : dMron + 98 jAperMag6 : + 99 jAperMag6Err : +100 jGausig : +101 jEll : +102 jPA : posangle +103 jErrBits : +104 jDeblend : +105 jClass : psfQF +106 jClassStat : +107 jppErrBits : photFlags +108 jSeqNum : detID +109 jObjID : +110 jXi : +111 jEta : + +112 hHallMag : M +113 hHallMagErr : dM +114 hPetroMag : +115 hPetroMagErr : +116 hAperMag3 : Map +117 hAperMag3Err : dMap +118 hAperMag4 : Mkron +119 hAperMag4Err : dMron +120 hAperMag6 : +121 hAperMag6Err : +122 hGausig : +123 hEll : +124 hPA : posangle +125 hErrBits : +126 hDeblend : +127 hClass : psfQF +128 hClassStat : +129 hppErrBits : photFlags +130 hSeqNum : detID +131 hObjID : +132 hXi : +133 hEta : + +134 k_1HallMag : M +135 k_1HallMagErr : dM +136 k_1PetroMag : +137 k_1PetroMagErr : +138 k_1AperMag3 : Map +139 k_1AperMag3Err : dMap +140 k_1AperMag4 : Mkron +141 k_1AperMag4Err : dMron +142 k_1AperMag6 : +143 k_1AperMag6Err : +144 k_1Gausig : +145 k_1Ell : +146 k_1PA : posangle +147 k_1ErrBits : +148 k_1Deblend : +149 k_1Class : psfQF +150 k_1ClassStat : +151 k_1ppErrBits : photFlags +152 k_1SeqNum : detID +153 k_1ObjID : +154 k_1Xi : +155 k_1Eta : + +156 k_2HallMag : M +157 k_2HallMagErr : dM +158 k_2PetroMag : +159 k_2PetroMagErr : +160 k_2AperMag3 : Map +161 k_2AperMag3Err : dMap +162 k_2AperMag4 : Mkron +163 k_2AperMag4Err : dMron +164 k_2AperMag6 : +165 k_2AperMag6Err : +166 k_2Gausig : +167 k_2Ell : +168 k_2PA : posangle +169 k_2ErrBits : +170 k_2Deblend : +171 k_2Class : psfQF +172 k_2ClassStat : +173 k_2ppErrBits : photFlags +174 k_2SeqNum : detID +175 k_2ObjID : +176 k_2Xi : +177 k_2Eta : + +*/ + +# define iPARSE(NS,NE,FIELD,NAME) { \ + cA = iparse_csv_rpt (&ivalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = ivalue; } + +# define jPARSE(NS,NE,FIELD,NAME) { \ + cA = jparse_csv_rpt (&jvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = jvalue; } + +# define dPARSE(NS,NE,FIELD,NAME) { \ + cA = dparse_csv_rpt (&dvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = dvalue; } + +static int Nline_read = 0; // track number of lines read so far (use to skip lines as well) +static int Nskip = 0; // UKIRT UHS CSV files have a single header row (and no special character to mark) + +UKIRT_Stars *loadukirt_uhs_readstars_ugcs (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars) { + + int codeZ = GetPhotcodeCodebyName ("UKIRT_Z"); if (!codeZ) Shutdown ("missing photcode UKIRT_Z"); + int codeY = GetPhotcodeCodebyName ("UKIRT_Y"); if (!codeY) Shutdown ("missing photcode UKIRT_Y"); + int codeJ = GetPhotcodeCodebyName ("UKIRT_J"); if (!codeJ) Shutdown ("missing photcode UKIRT_J"); + int codeH = GetPhotcodeCodebyName ("UKIRT_H"); if (!codeH) Shutdown ("missing photcode UKIRT_H"); + int codeK = GetPhotcodeCodebyName ("UKIRT_K"); if (!codeK) Shutdown ("missing photcode UKIRT_K"); + + // XXX I need a UKIRT UHS Mean Epoch + // time_t UKIRT_UHS_EPOCH = ohana_date_to_sec ("2016/01/01,00:00:00"); + // fprintf (stderr, "WARNING: using an invalid UKIRT_UHS_EPOCH (see loadukirt_uhs_readstars.c:60)\n"); + + int Nelem = 0; // number of valid rows read (vector elements) + int NELEM = 10000; // currently-allocated number of output rows + + // vectors to hold the data loaded from the csv file (names are UKIRT native names) + ALLOCATE_PTR (sourceID, uint64_t, NELEM); + ALLOCATE_PTR (frameSetID, uint64_t, NELEM); // XXX doc says 8-bytes for this value! + ALLOCATE_PTR (ra, double, NELEM); + ALLOCATE_PTR (dec, double, NELEM); + ALLOCATE_PTR (epoch, double, NELEM); + ALLOCATE_PTR (priOrSec, uint64_t, NELEM); + ALLOCATE_PTR (mergedClass, int, NELEM); // XXX doc says 2-bytes (could use a short) + ALLOCATE_PTR (pStar, double, NELEM); + ALLOCATE_PTR (pGalaxy, double, NELEM); + + ALLOCATE_PTR (zHallMag, double, NELEM); + ALLOCATE_PTR (zHallMagErr, double, NELEM); + ALLOCATE_PTR (zAperMag3, double, NELEM); + ALLOCATE_PTR (zAperMag3Err, double, NELEM); + ALLOCATE_PTR (zAperMag4, double, NELEM); + ALLOCATE_PTR (zAperMag4Err, double, NELEM); + ALLOCATE_PTR (zPA, double, NELEM); + ALLOCATE_PTR (zClass, int, NELEM); + ALLOCATE_PTR (zppErrBits, int, NELEM); + ALLOCATE_PTR (zSeqNum, int, NELEM); + + ALLOCATE_PTR (yHallMag, double, NELEM); + ALLOCATE_PTR (yHallMagErr, double, NELEM); + ALLOCATE_PTR (yAperMag3, double, NELEM); + ALLOCATE_PTR (yAperMag3Err, double, NELEM); + ALLOCATE_PTR (yAperMag4, double, NELEM); + ALLOCATE_PTR (yAperMag4Err, double, NELEM); + ALLOCATE_PTR (yPA, double, NELEM); + ALLOCATE_PTR (yClass, int, NELEM); + ALLOCATE_PTR (yppErrBits, int, NELEM); + ALLOCATE_PTR (ySeqNum, int, NELEM); + + ALLOCATE_PTR (jHallMag, double, NELEM); + ALLOCATE_PTR (jHallMagErr, double, NELEM); + ALLOCATE_PTR (jAperMag3, double, NELEM); + ALLOCATE_PTR (jAperMag3Err, double, NELEM); + ALLOCATE_PTR (jAperMag4, double, NELEM); + ALLOCATE_PTR (jAperMag4Err, double, NELEM); + ALLOCATE_PTR (jPA, double, NELEM); + ALLOCATE_PTR (jClass, int, NELEM); + ALLOCATE_PTR (jppErrBits, int, NELEM); + ALLOCATE_PTR (jSeqNum, int, NELEM); + + ALLOCATE_PTR (hHallMag, double, NELEM); + ALLOCATE_PTR (hHallMagErr, double, NELEM); + ALLOCATE_PTR (hAperMag3, double, NELEM); + ALLOCATE_PTR (hAperMag3Err, double, NELEM); + ALLOCATE_PTR (hAperMag4, double, NELEM); + ALLOCATE_PTR (hAperMag4Err, double, NELEM); + ALLOCATE_PTR (hPA, double, NELEM); + ALLOCATE_PTR (hClass, int, NELEM); + ALLOCATE_PTR (hppErrBits, int, NELEM); + ALLOCATE_PTR (hSeqNum, int, NELEM); + + ALLOCATE_PTR (k1HallMag, double, NELEM); + ALLOCATE_PTR (k1HallMagErr, double, NELEM); + ALLOCATE_PTR (k1AperMag3, double, NELEM); + ALLOCATE_PTR (k1AperMag3Err, double, NELEM); + ALLOCATE_PTR (k1AperMag4, double, NELEM); + ALLOCATE_PTR (k1AperMag4Err, double, NELEM); + ALLOCATE_PTR (k1PA, double, NELEM); + ALLOCATE_PTR (k1Class, int, NELEM); + ALLOCATE_PTR (k1ppErrBits, int, NELEM); + ALLOCATE_PTR (k1SeqNum, int, NELEM); + + ALLOCATE_PTR (k2HallMag, double, NELEM); + ALLOCATE_PTR (k2HallMagErr, double, NELEM); + ALLOCATE_PTR (k2AperMag3, double, NELEM); + ALLOCATE_PTR (k2AperMag3Err, double, NELEM); + ALLOCATE_PTR (k2AperMag4, double, NELEM); + ALLOCATE_PTR (k2AperMag4Err, double, NELEM); + ALLOCATE_PTR (k2PA, double, NELEM); + ALLOCATE_PTR (k2Class, int, NELEM); + ALLOCATE_PTR (k2ppErrBits, int, NELEM); + ALLOCATE_PTR (k2SeqNum, int, NELEM); + + // 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 = *nstart; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + while (!EndOfFile && (Nelem < NSTARS_MAX)) { + int Nbytes = BUFFER_SIZE - 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; } + + // for UKIRT UHS, 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 : TBD for UHS + + // the start of the line is the 1st element (fields are 1-counting) + jPARSE ( 1, 1, sourceID, "sourceID"); + iPARSE ( 3, 1, frameSetID, "frameSetID"); + dPARSE ( 4, 3, ra, "ra"); + dPARSE ( 5, 4, dec, "dec"); + dPARSE ( 8, 5, epoch, "epoch"); + jPARSE ( 23, 8, priOrSec, "priOrSec"); // watch for epoch + iPARSE ( 41, 23, mergedClass, "mergedClass"); + dPARSE ( 42, 41, pStar, "pStar"); + dPARSE ( 43, 42, pGalaxy, "pGalaxy"); + dPARSE ( 46, 43, zHallMag, "zHallMag"); + dPARSE ( 47, 46, zHallMagErr, "zHallMagErr"); + dPARSE ( 50, 47, zAperMag3, "zAperMag3"); + dPARSE ( 51, 50, zAperMag3Err, "zAperMag3Err"); + dPARSE ( 52, 51, zAperMag4, "zAperMag4"); + dPARSE ( 53, 52, zAperMag4Err, "zAperMag4Err"); + dPARSE ( 58, 53, zPA, "zPA"); + iPARSE ( 61, 58, zClass, "zClass"); + iPARSE ( 63, 61, zppErrBits, "zppErrBits"); + iPARSE ( 64, 63, zSeqNum, "zSeqNum"); + dPARSE ( 68, 64, yHallMag, "yHallMag"); + dPARSE ( 69, 68, yHallMagErr, "yHallMagErr"); + dPARSE ( 72, 69, yAperMag3, "yAperMag3"); + dPARSE ( 73, 72, yAperMag3Err, "yAperMag3Err"); + dPARSE ( 74, 73, yAperMag4, "yAperMag4"); + dPARSE ( 75, 74, yAperMag4Err, "yAperMag4Err"); + dPARSE ( 80, 75, yPA, "yPA"); + iPARSE ( 83, 80, yClass, "yClass"); + iPARSE ( 85, 83, yppErrBits, "yppErrBits"); + iPARSE ( 86, 85, ySeqNum, "ySeqNum"); + dPARSE ( 90, 86, jHallMag, "jHallMag"); + dPARSE ( 91, 90, jHallMagErr, "jHallMagErr"); + dPARSE ( 94, 91, jAperMag3, "jAperMag3"); + dPARSE ( 95, 94, jAperMag3Err, "jAperMag3Err"); + dPARSE ( 96, 95, jAperMag4, "jAperMag4"); + dPARSE ( 97, 96, jAperMag4Err, "jAperMag4Err"); + dPARSE (102, 97, jPA, "jPA"); + iPARSE (105, 102, jClass, "jClass"); + iPARSE (107, 105, jppErrBits, "jppErrBits"); + iPARSE (108, 107, jSeqNum, "jSeqNum"); + dPARSE (112, 108, hHallMag, "hHallMag"); + dPARSE (113, 112, hHallMagErr, "hHallMagErr"); + dPARSE (116, 113, hAperMag3, "hAperMag3"); + dPARSE (117, 116, hAperMag3Err, "hAperMag3Err"); + dPARSE (118, 117, hAperMag4, "hAperMag4"); + dPARSE (119, 118, hAperMag4Err, "hAperMag4Err"); + dPARSE (124, 119, hPA, "hPA"); + iPARSE (127, 124, hClass, "hClass"); + iPARSE (129, 127, hppErrBits, "hppErrBits"); + iPARSE (130, 129, hSeqNum, "hSeqNum"); + dPARSE (134, 130, k1HallMag, "k1HallMag"); + dPARSE (135, 134, k1HallMagErr, "k1HallMagErr"); + dPARSE (138, 135, k1AperMag3, "k1AperMag3"); + dPARSE (139, 138, k1AperMag3Err, "k1AperMag3Err"); + dPARSE (140, 139, k1AperMag4, "k1AperMag4"); + dPARSE (141, 140, k1AperMag4Err, "k1AperMag4Err"); + dPARSE (146, 141, k1PA, "k1PA"); + iPARSE (149, 146, k1Class, "k1Class"); + iPARSE (151, 149, k1ppErrBits, "k1ppErrBits"); + iPARSE (152, 151, k1SeqNum, "k1SeqNum"); + dPARSE (156, 152, k2HallMag, "k2HallMag"); + dPARSE (157, 156, k2HallMagErr, "k2HallMagErr"); + dPARSE (160, 157, k2AperMag3, "k2AperMag3"); + dPARSE (161, 160, k2AperMag3Err, "k2AperMag3Err"); + dPARSE (162, 161, k2AperMag4, "k2AperMag4"); + dPARSE (163, 162, k2AperMag4Err, "k2AperMag4Err"); + dPARSE (168, 163, k2PA, "k2PA"); + iPARSE (171, 168, k2Class, "k2Class"); + iPARSE (173, 171, k2ppErrBits, "k2ppErrBits"); + iPARSE (174, 173, k2SeqNum, "k2SeqNum"); + + 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 (sourceID, uint64_t, NELEM); + REALLOCATE (frameSetID, uint64_t, NELEM); + REALLOCATE (ra, double, NELEM); + REALLOCATE (dec, double, NELEM); + REALLOCATE (priOrSec, uint64_t, NELEM); + REALLOCATE (mergedClass, int, NELEM); + REALLOCATE (pStar, double, NELEM); + REALLOCATE (pGalaxy, double, NELEM); + REALLOCATE (zHallMag, double, NELEM); + REALLOCATE (zHallMagErr, double, NELEM); + REALLOCATE (zAperMag3, double, NELEM); + REALLOCATE (zAperMag3Err, double, NELEM); + REALLOCATE (zAperMag4, double, NELEM); + REALLOCATE (zAperMag4Err, double, NELEM); + REALLOCATE (zPA, double, NELEM); + REALLOCATE (zClass, int, NELEM); + REALLOCATE (zppErrBits, int, NELEM); + REALLOCATE (zSeqNum, int, NELEM); + REALLOCATE (yHallMag, double, NELEM); + REALLOCATE (yHallMagErr, double, NELEM); + REALLOCATE (yAperMag3, double, NELEM); + REALLOCATE (yAperMag3Err, double, NELEM); + REALLOCATE (yAperMag4, double, NELEM); + REALLOCATE (yAperMag4Err, double, NELEM); + REALLOCATE (yPA, double, NELEM); + REALLOCATE (yClass, int, NELEM); + REALLOCATE (yppErrBits, int, NELEM); + REALLOCATE (ySeqNum, int, NELEM); + REALLOCATE (jHallMag, double, NELEM); + REALLOCATE (jHallMagErr, double, NELEM); + REALLOCATE (jAperMag3, double, NELEM); + REALLOCATE (jAperMag3Err, double, NELEM); + REALLOCATE (jAperMag4, double, NELEM); + REALLOCATE (jAperMag4Err, double, NELEM); + REALLOCATE (jPA, double, NELEM); + REALLOCATE (jClass, int, NELEM); + REALLOCATE (jppErrBits, int, NELEM); + REALLOCATE (jSeqNum, int, NELEM); + REALLOCATE (hHallMag, double, NELEM); + REALLOCATE (hHallMagErr, double, NELEM); + REALLOCATE (hAperMag3, double, NELEM); + REALLOCATE (hAperMag3Err, double, NELEM); + REALLOCATE (hAperMag4, double, NELEM); + REALLOCATE (hAperMag4Err, double, NELEM); + REALLOCATE (hPA, double, NELEM); + REALLOCATE (hClass, int, NELEM); + REALLOCATE (hppErrBits, int, NELEM); + REALLOCATE (hSeqNum, int, NELEM); + REALLOCATE (k1HallMag, double, NELEM); + REALLOCATE (k1HallMagErr, double, NELEM); + REALLOCATE (k1AperMag3, double, NELEM); + REALLOCATE (k1AperMag3Err, double, NELEM); + REALLOCATE (k1AperMag4, double, NELEM); + REALLOCATE (k1AperMag4Err, double, NELEM); + REALLOCATE (k1PA, double, NELEM); + REALLOCATE (k1Class, int, NELEM); + REALLOCATE (k1ppErrBits, int, NELEM); + REALLOCATE (k1SeqNum, int, NELEM); + REALLOCATE (k2HallMag, double, NELEM); + REALLOCATE (k2HallMagErr, double, NELEM); + REALLOCATE (k2AperMag3, double, NELEM); + REALLOCATE (k2AperMag3Err, double, NELEM); + REALLOCATE (k2AperMag4, double, NELEM); + REALLOCATE (k2AperMag4Err, double, NELEM); + REALLOCATE (k2PA, double, NELEM); + REALLOCATE (k2Class, int, NELEM); + REALLOCATE (k2ppErrBits, int, NELEM); + REALLOCATE (k2SeqNum, 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; + + ALLOCATE_PTR (stars, UKIRT_Stars, NSTARS); + + for (int i = 0; i < NstarsIn; i++) { + + Rmin = MIN (Rmin, ra[i]); + Rmax = MAX (Rmax, ra[i]); + Dmin = MIN (Dmin, dec[i]); + Dmax = MAX (Dmax, dec[i]); + + // only allocate the measures for stars as we define them + ALLOCATE (stars[Nstars].measure, Measure, UKIRT_NFILTER); + + // we have NFILTER (== 6) UKIRT measurements per object + dvo_average_init (&stars[Nstars].average); + for (int j = 0; j < UKIRT_NFILTER; j++) { + dvo_measure_init (&stars[Nstars].measure[j]); + } + + stars[Nstars].average.R = ra[i]; + stars[Nstars].average.D = dec[i]; + stars[Nstars].average.dR = NAN; + stars[Nstars].average.dD = NAN; + + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + int isPrimary = (!priOrSec[i] || (priOrSec[i] == frameSetID[i])) ? ID_MEAS_STACK_PRIMARY : 0x00; + int photFlags2 = isPrimary; + switch (mergedClass[i]) { + case 1: photFlags2 |= 0x01; break; // Galaxy + case 0: photFlags2 |= 0x02; break; // Noise + case -1: photFlags2 |= 0x04; break; // Star + case -2: photFlags2 |= 0x08; break; // probably star + case -3: photFlags2 |= 0x10; break; // probably galaxy + case -9: photFlags2 |= 0x20; break; // saturated + default: break; + } + + // y measure + stars[Nstars].measure[0].extID = sourceID[i]; + stars[Nstars].measure[0].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[0].R = ra[i]; + stars[Nstars].measure[0].D = dec[i]; + stars[Nstars].measure[0].psfChisq = pStar[i]; + stars[Nstars].measure[0].extNsigma = pGalaxy[i]; + stars[Nstars].measure[0].M = yHallMag[i]; + stars[Nstars].measure[0].dM = yHallMagErr[i]; + stars[Nstars].measure[0].Map = yAperMag3[i]; + stars[Nstars].measure[0].dMap = yAperMag3Err[i]; + stars[Nstars].measure[0].Mkron = yAperMag4[i]; + stars[Nstars].measure[0].dMkron = yAperMag4Err[i]; + stars[Nstars].measure[0].posangle = yPA[i]; // XXX units + stars[Nstars].measure[0].psfQF = yClass[i]; + stars[Nstars].measure[0].photFlags = yppErrBits[i]; + stars[Nstars].measure[0].photFlags2 = photFlags2; + stars[Nstars].measure[0].detID = ySeqNum[i]; + stars[Nstars].measure[0].photcode = codeY; + + // z measure + stars[Nstars].measure[1].extID = sourceID[i]; + stars[Nstars].measure[1].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[1].R = ra[i]; + stars[Nstars].measure[1].D = dec[i]; + stars[Nstars].measure[1].psfChisq = pStar[i]; + stars[Nstars].measure[1].extNsigma = pGalaxy[i]; + stars[Nstars].measure[1].M = zHallMag[i]; + stars[Nstars].measure[1].dM = zHallMagErr[i]; + stars[Nstars].measure[1].Map = zAperMag3[i]; + stars[Nstars].measure[1].dMap = zAperMag3Err[i]; + stars[Nstars].measure[1].Mkron = zAperMag4[i]; + stars[Nstars].measure[1].dMkron = zAperMag4Err[i]; + stars[Nstars].measure[1].posangle = zPA[i]; // XXX units + stars[Nstars].measure[1].psfQF = zClass[i]; + stars[Nstars].measure[1].photFlags = zppErrBits[i]; + stars[Nstars].measure[1].photFlags2 = photFlags2; + stars[Nstars].measure[1].detID = zSeqNum[i]; + stars[Nstars].measure[1].photcode = codeZ; + + // J measure + stars[Nstars].measure[2].extID = sourceID[i]; + stars[Nstars].measure[2].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[2].R = ra[i]; + stars[Nstars].measure[2].D = dec[i]; + stars[Nstars].measure[2].psfChisq = pStar[i]; + stars[Nstars].measure[2].extNsigma = pGalaxy[i]; + stars[Nstars].measure[2].M = jHallMag[i]; + stars[Nstars].measure[2].dM = jHallMagErr[i]; + stars[Nstars].measure[2].Map = jAperMag3[i]; + stars[Nstars].measure[2].dMap = jAperMag3Err[i]; + stars[Nstars].measure[2].Mkron = jAperMag4[i]; + stars[Nstars].measure[2].dMkron = jAperMag4Err[i]; + stars[Nstars].measure[2].posangle = jPA[i]; // XXX units + stars[Nstars].measure[2].psfQF = jClass[i]; + stars[Nstars].measure[2].photFlags = jppErrBits[i]; + stars[Nstars].measure[2].photFlags2 = photFlags2; + stars[Nstars].measure[2].detID = jSeqNum[i]; + stars[Nstars].measure[2].photcode = codeJ; + + // H measure + stars[Nstars].measure[3].extID = sourceID[i]; + stars[Nstars].measure[3].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[3].R = ra[i]; + stars[Nstars].measure[3].D = dec[i]; + stars[Nstars].measure[3].psfChisq = pStar[i]; + stars[Nstars].measure[3].extNsigma = pGalaxy[i]; + stars[Nstars].measure[3].M = hHallMag[i]; + stars[Nstars].measure[3].dM = hHallMagErr[i]; + stars[Nstars].measure[3].Map = hAperMag3[i]; + stars[Nstars].measure[3].dMap = hAperMag3Err[i]; + stars[Nstars].measure[3].Mkron = hAperMag4[i]; + stars[Nstars].measure[3].dMkron = hAperMag4Err[i]; + stars[Nstars].measure[3].posangle = hPA[i]; // XXX units + stars[Nstars].measure[3].psfQF = hClass[i]; + stars[Nstars].measure[3].photFlags = hppErrBits[i]; + stars[Nstars].measure[3].photFlags2 = photFlags2; + stars[Nstars].measure[3].detID = hSeqNum[i]; + stars[Nstars].measure[3].photcode = codeH; + + // K1 measure + stars[Nstars].measure[4].extID = sourceID[i]; + stars[Nstars].measure[4].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[4].R = ra[i]; + stars[Nstars].measure[4].D = dec[i]; + stars[Nstars].measure[4].psfChisq = pStar[i]; + stars[Nstars].measure[4].extNsigma = pGalaxy[i]; + stars[Nstars].measure[4].M = k1HallMag[i]; + stars[Nstars].measure[4].dM = k1HallMagErr[i]; + stars[Nstars].measure[4].Map = k1AperMag3[i]; + stars[Nstars].measure[4].dMap = k1AperMag3Err[i]; + stars[Nstars].measure[4].Mkron = k1AperMag4[i]; + stars[Nstars].measure[4].dMkron = k1AperMag4Err[i]; + stars[Nstars].measure[4].posangle = k1PA[i]; // XXX units + stars[Nstars].measure[4].psfQF = k1Class[i]; + stars[Nstars].measure[4].photFlags = k1ppErrBits[i]; + stars[Nstars].measure[4].photFlags2 = photFlags2; + stars[Nstars].measure[4].detID = k1SeqNum[i]; + stars[Nstars].measure[4].photcode = codeK; + + // K1 measure + stars[Nstars].measure[5].extID = sourceID[i]; + stars[Nstars].measure[5].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[5].R = ra[i]; + stars[Nstars].measure[5].D = dec[i]; + stars[Nstars].measure[5].psfChisq = pStar[i]; + stars[Nstars].measure[5].extNsigma = pGalaxy[i]; + stars[Nstars].measure[5].M = k2HallMag[i]; + stars[Nstars].measure[5].dM = k2HallMagErr[i]; + stars[Nstars].measure[5].Map = k2AperMag3[i]; + stars[Nstars].measure[5].dMap = k2AperMag3Err[i]; + stars[Nstars].measure[5].Mkron = k2AperMag4[i]; + stars[Nstars].measure[5].dMkron = k2AperMag4Err[i]; + stars[Nstars].measure[5].posangle = k2PA[i]; // XXX units + stars[Nstars].measure[5].psfQF = k2Class[i]; + stars[Nstars].measure[5].photFlags = k2ppErrBits[i]; + stars[Nstars].measure[5].photFlags2 = photFlags2; + stars[Nstars].measure[5].detID = k2SeqNum[i]; + stars[Nstars].measure[5].photcode = codeK; + + stars[Nstars].average.Nmeasure = UKIRT_NFILTER; + Nstars ++; + + if (NSTARS >= Nstars) { + NSTARS += 10000; + REALLOCATE (stars, UKIRT_Stars, NSTARS); + } + } + + FREE (sourceID); + FREE (frameSetID); + FREE (ra); + FREE (dec); + FREE (epoch); + FREE (priOrSec); + FREE (mergedClass); + FREE (pStar); + FREE (pGalaxy); + + FREE (zHallMag); + FREE (zHallMagErr); + FREE (zAperMag3); + FREE (zAperMag3Err); + FREE (zAperMag4); + FREE (zAperMag4Err); + FREE (zPA); + FREE (zClass); + FREE (zppErrBits); + FREE (zSeqNum); + + FREE (yHallMag); + FREE (yHallMagErr); + FREE (yAperMag3); + FREE (yAperMag3Err); + FREE (yAperMag4); + FREE (yAperMag4Err); + FREE (yPA); + FREE (yClass); + FREE (yppErrBits); + FREE (ySeqNum); + + FREE (jHallMag); + FREE (jHallMagErr); + FREE (jAperMag3); + FREE (jAperMag3Err); + FREE (jAperMag4); + FREE (jAperMag4Err); + FREE (jPA); + FREE (jClass); + FREE (jppErrBits); + FREE (jSeqNum); + + FREE (hHallMag); + FREE (hHallMagErr); + FREE (hAperMag3); + FREE (hAperMag3Err); + FREE (hAperMag4); + FREE (hAperMag4Err); + FREE (hPA); + FREE (hClass); + FREE (hppErrBits); + FREE (hSeqNum); + + FREE (k1HallMag); + FREE (k1HallMagErr); + FREE (k1AperMag3); + FREE (k1AperMag3Err); + FREE (k1AperMag4); + FREE (k1AperMag4Err); + FREE (k1PA); + FREE (k1Class); + FREE (k1ppErrBits); + FREE (k1SeqNum); + + FREE (k2HallMag); + FREE (k2HallMagErr); + FREE (k2AperMag3); + FREE (k2AperMag3Err); + FREE (k2AperMag4); + FREE (k2AperMag4Err); + FREE (k2PA); + FREE (k2Class); + FREE (k2ppErrBits); + FREE (k2SeqNum); + + *nstars = Nstars; + *nstart = Nstart; + return (stars); +} Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ugps.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ugps.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ugps.c (revision 41151) @@ -0,0 +1,682 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +/* this function reads the values of interest from the UKIRT UGCS CSV files: + + data model description: doc/mapping_UKIDSS_unWISE_DVO.pdf + + N -- field (column) number in CSV file (1 counting) + | UKIRT field name : DVO measure field + 1 sourceID : extID + 2 cuEventID : + 3 frameSetID : imageID (assumes we do not try to go backwards for these photcodes) + 4 ra : R + 5 dec : D + + 6 sigRa : + 7 sigDec : + 8 epoch : t ???? -- how is there a single epoch? + 9 muRa : + 10 muDec : + 11 sigMuRa : + 12 sigMuDec : + 13 chi2 : + 14 nFrames : + + 15 cx : + 16 cy : + 17 cz : + 18 htmID : + 19 l : + 20 b : + 21 lambda : + 22 eta : + 23 priOrSec : primary bit in dbFlags + + 24 jmhPnt : + 25 jmhPntErr : + 26 hmk1Pnt : + 27 hmk1PntErr : + 28 h2mk1Pnt : + 29 h2mk1PntErr : + + 30 mergedClassStat : + 31 mergedClass : psfQFperf or photFlags2 + 32 pStar : psfChisq + 33 pGalaxy : extNsigma + 34 pNoise : + 35 pSaturated : + + 36 jAperMag1 : M + 37 jAperMag1Err : dM + 38 jAperMag3 : Map + 39 jAperMag3Err : dMap + 40 jAperMag4 : Mkron + 41 jAperMag4Err : dMron + 42 jGausig : + 43 jEll : + 44 jPA : posangle + 45 jErrBits : + 46 jDeblend : + 47 jClass : psfQF + 48 jClassStat : + 49 jppErrBits : photFlags + 50 jSeqNum : detID + 51 jObjID : + 52 jXi : + 53 jEta : + 54 hAperMag1 : M + 55 hAperMag1Err : dM + 56 hAperMag3 : Map + 57 hAperMag3Err : dMap + 58 hAperMag4 : Mkron + 59 hAperMag4Err : dMron + 60 hGausig : + 61 hEll : + 62 hPA : posangle + 63 hErrBits : + 64 hDeblend : + 65 hClass : psfQF + 66 hClassStat : + 67 hppErrBits : photFlags + 68 hSeqNum : detID + 69 hObjID : + 70 hXi : + 71 hEta : + 72 k_1AperMag1 : M + 73 k_1AperMag1Err : dM + 74 k_1AperMag3 : Map + 75 k_1AperMag3Err : dMap + 76 k_1AperMag4 : Mkron + 77 k_1AperMag4Err : dMron + 78 k_1Gausig : + 79 k_1Ell : + 80 k_1PA : posangle + 81 k_1ErrBits : + 82 k_1Deblend : + 83 k_1Class : psfQF + 84 k_1ClassStat : + 85 k_1ppErrBits : photFlags + 86 k_1SeqNum : detID + 87 k_1ObjID : + 88 k_1Xi : + 89 k_1Eta : + 90 k_2AperMag1 : M + 91 k_2AperMag1Err : dM + 92 k_2AperMag3 : Map + 93 k_2AperMag3Err : dMap + 94 k_2AperMag4 : Mkron + 95 k_2AperMag4Err : dMron + 96 k_2Gausig : + 97 k_2Ell : + 98 k_2PA : posangle + 99 k_2ErrBits : +100 k_2Deblend : +101 k_2Class : psfQF +102 k_2ClassStat : +103 k_2ppErrBits : photFlags +104 k_2SeqNum : detID +105 k_2ObjID : +106 k_2Xi : +107 k_2Eta : +108 h_2AperMag1 : M +109 h_2AperMag1Err : dM +110 h_2AperMag3 : Map +111 h_2AperMag3Err : dMap +112 h_2AperMag4 : Mkron +113 h_2AperMag4Err : dMron +114 h_2Gausig : +115 h_2Ell : +116 h_2PA : posangle +117 h_2ErrBits : +118 h_2Deblend : +119 h_2Class : psfQF +120 h_2ClassStat : +121 h_2ppErrBits : photFlags +122 h_2SeqNum : detID +123 h_2ObjID : +124 h_2Xi : +125 h_2Eta : + +*/ + +# define iPARSE(NS,NE,FIELD,NAME) { \ + cA = iparse_csv_rpt (&ivalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = ivalue; } + +# define jPARSE(NS,NE,FIELD,NAME) { \ + cA = jparse_csv_rpt (&jvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = jvalue; } + +# define dPARSE(NS,NE,FIELD,NAME) { \ + cA = dparse_csv_rpt (&dvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = dvalue; } + +static int Nline_read = 0; // track number of lines read so far (use to skip lines as well) +static int Nskip = 0; // UKIRT UHS CSV files have a single header row (and no special character to mark) + +UKIRT_Stars *loadukirt_uhs_readstars_ugps (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars) { + + int codeJ = GetPhotcodeCodebyName ("UKIRT_J"); if (!codeJ) Shutdown ("missing photcode UKIRT_J"); + int codeH = GetPhotcodeCodebyName ("UKIRT_H"); if (!codeH) Shutdown ("missing photcode UKIRT_H"); + int codeK = GetPhotcodeCodebyName ("UKIRT_K"); if (!codeK) Shutdown ("missing photcode UKIRT_K"); + + // XXX I need a UKIRT UHS Mean Epoch + // time_t UKIRT_UHS_EPOCH = ohana_date_to_sec ("2016/01/01,00:00:00"); + // fprintf (stderr, "WARNING: using an invalid UKIRT_UHS_EPOCH (see loadukirt_uhs_readstars.c:60)\n"); + + int Nelem = 0; // number of valid rows read (vector elements) + int NELEM = 10000; // currently-allocated number of output rows + + // vectors to hold the data loaded from the csv file (names are UKIRT native names) + ALLOCATE_PTR (sourceID, uint64_t, NELEM); + ALLOCATE_PTR (frameSetID, uint64_t, NELEM); // XXX doc says 8-bytes for this value! + ALLOCATE_PTR (ra, double, NELEM); + ALLOCATE_PTR (dec, double, NELEM); + ALLOCATE_PTR (epoch, double, NELEM); + ALLOCATE_PTR (priOrSec, uint64_t, NELEM); + ALLOCATE_PTR (mergedClass, int, NELEM); // XXX doc says 2-bytes (could use a short) + ALLOCATE_PTR (pStar, double, NELEM); + ALLOCATE_PTR (pGalaxy, double, NELEM); + + ALLOCATE_PTR (jAperMag1, double, NELEM); + ALLOCATE_PTR (jAperMag1Err, double, NELEM); + ALLOCATE_PTR (jAperMag3, double, NELEM); + ALLOCATE_PTR (jAperMag3Err, double, NELEM); + ALLOCATE_PTR (jAperMag4, double, NELEM); + ALLOCATE_PTR (jAperMag4Err, double, NELEM); + ALLOCATE_PTR (jPA, double, NELEM); + ALLOCATE_PTR (jClass, int, NELEM); + ALLOCATE_PTR (jppErrBits, int, NELEM); + ALLOCATE_PTR (jSeqNum, int, NELEM); + + ALLOCATE_PTR (hAperMag1, double, NELEM); + ALLOCATE_PTR (hAperMag1Err, double, NELEM); + ALLOCATE_PTR (hAperMag3, double, NELEM); + ALLOCATE_PTR (hAperMag3Err, double, NELEM); + ALLOCATE_PTR (hAperMag4, double, NELEM); + ALLOCATE_PTR (hAperMag4Err, double, NELEM); + ALLOCATE_PTR (hPA, double, NELEM); + ALLOCATE_PTR (hClass, int, NELEM); + ALLOCATE_PTR (hppErrBits, int, NELEM); + ALLOCATE_PTR (hSeqNum, int, NELEM); + + ALLOCATE_PTR (k1AperMag1, double, NELEM); + ALLOCATE_PTR (k1AperMag1Err, double, NELEM); + ALLOCATE_PTR (k1AperMag3, double, NELEM); + ALLOCATE_PTR (k1AperMag3Err, double, NELEM); + ALLOCATE_PTR (k1AperMag4, double, NELEM); + ALLOCATE_PTR (k1AperMag4Err, double, NELEM); + ALLOCATE_PTR (k1PA, double, NELEM); + ALLOCATE_PTR (k1Class, int, NELEM); + ALLOCATE_PTR (k1ppErrBits, int, NELEM); + ALLOCATE_PTR (k1SeqNum, int, NELEM); + + ALLOCATE_PTR (k2AperMag1, double, NELEM); + ALLOCATE_PTR (k2AperMag1Err, double, NELEM); + ALLOCATE_PTR (k2AperMag3, double, NELEM); + ALLOCATE_PTR (k2AperMag3Err, double, NELEM); + ALLOCATE_PTR (k2AperMag4, double, NELEM); + ALLOCATE_PTR (k2AperMag4Err, double, NELEM); + ALLOCATE_PTR (k2PA, double, NELEM); + ALLOCATE_PTR (k2Class, int, NELEM); + ALLOCATE_PTR (k2ppErrBits, int, NELEM); + ALLOCATE_PTR (k2SeqNum, int, NELEM); + + ALLOCATE_PTR (h2AperMag1, double, NELEM); + ALLOCATE_PTR (h2AperMag1Err, double, NELEM); + ALLOCATE_PTR (h2AperMag3, double, NELEM); + ALLOCATE_PTR (h2AperMag3Err, double, NELEM); + ALLOCATE_PTR (h2AperMag4, double, NELEM); + ALLOCATE_PTR (h2AperMag4Err, double, NELEM); + ALLOCATE_PTR (h2PA, double, NELEM); + ALLOCATE_PTR (h2Class, int, NELEM); + ALLOCATE_PTR (h2ppErrBits, int, NELEM); + ALLOCATE_PTR (h2SeqNum, int, NELEM); + + // 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 = *nstart; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + while (!EndOfFile && (Nelem < NSTARS_MAX)) { + int Nbytes = BUFFER_SIZE - 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; } + + // for UKIRT UHS, 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 : TBD for UHS + + // the start of the line is the 1st element (fields are 1-counting) + jPARSE ( 1, 1, sourceID, "sourceID"); + iPARSE ( 3, 1, frameSetID, "frameSetID"); + dPARSE ( 4, 3, ra, "ra"); + dPARSE ( 5, 4, dec, "dec"); + dPARSE ( 8, 5, epoch, "epoch"); + jPARSE ( 23, 8, priOrSec, "priOrSec"); // watch for epoch + iPARSE ( 31, 23, mergedClass, "mergedClass"); + dPARSE ( 32, 31, pStar, "pStar"); + dPARSE ( 33, 32, pGalaxy, "pGalaxy"); + dPARSE ( 36, 33, jAperMag1, "jAperMag1"); + dPARSE ( 37, 36, jAperMag1Err, "jAperMag1Err"); + dPARSE ( 38, 37, jAperMag3, "jAperMag3"); + dPARSE ( 39, 38, jAperMag3Err, "jAperMag3Err"); + dPARSE ( 40, 39, jAperMag4, "jAperMag4"); + dPARSE ( 41, 40, jAperMag4Err, "jAperMag4Err"); + dPARSE ( 44, 41, jPA, "jPA"); + iPARSE ( 47, 44, jClass, "jClass"); + iPARSE ( 49, 47, jppErrBits, "jppErrBits"); + iPARSE ( 50, 49, jSeqNum, "jSeqNum"); + dPARSE ( 54, 50, hAperMag1, "hAperMag1"); + dPARSE ( 55, 54, hAperMag1Err, "hAperMag1Err"); + dPARSE ( 56, 55, hAperMag3, "hAperMag3"); + dPARSE ( 57, 56, hAperMag3Err, "hAperMag3Err"); + dPARSE ( 58, 57, hAperMag4, "hAperMag4"); + dPARSE ( 59, 58, hAperMag4Err, "hAperMag4Err"); + dPARSE ( 62, 59, hPA, "hPA"); + iPARSE ( 65, 62, hClass, "hClass"); + iPARSE ( 67, 65, hppErrBits, "hppErrBits"); + iPARSE ( 68, 67, hSeqNum, "hSeqNum"); + dPARSE ( 72, 68, k1AperMag1, "k1AperMag1"); + dPARSE ( 73, 72, k1AperMag1Err, "k1AperMag1Err"); + dPARSE ( 74, 73, k1AperMag3, "k1AperMag3"); + dPARSE ( 75, 74, k1AperMag3Err, "k1AperMag3Err"); + dPARSE ( 76, 75, k1AperMag4, "k1AperMag4"); + dPARSE ( 77, 76, k1AperMag4Err, "k1AperMag4Err"); + dPARSE ( 80, 77, k1PA, "k1PA"); + iPARSE ( 83, 80, k1Class, "k1Class"); + iPARSE ( 85, 83, k1ppErrBits, "k1ppErrBits"); + iPARSE ( 86, 85, k1SeqNum, "k1SeqNum"); + dPARSE ( 90, 86, k2AperMag1, "k2AperMag1"); + dPARSE ( 91, 90, k2AperMag1Err, "k2AperMag1Err"); + dPARSE ( 92, 91, k2AperMag3, "k2AperMag3"); + dPARSE ( 93, 92, k2AperMag3Err, "k2AperMag3Err"); + dPARSE ( 94, 93, k2AperMag4, "k2AperMag4"); + dPARSE ( 95, 94, k2AperMag4Err, "k2AperMag4Err"); + dPARSE ( 98, 95, k2PA, "k2PA"); + iPARSE (101, 98, k2Class, "k2Class"); + iPARSE (103, 101, k2ppErrBits, "k2ppErrBits"); + iPARSE (104, 103, k2SeqNum, "k2SeqNum"); + dPARSE (108, 104, h2AperMag1, "hAperMag1"); + dPARSE (109, 108, h2AperMag1Err, "hAperMag1Err"); + dPARSE (110, 109, h2AperMag3, "hAperMag3"); + dPARSE (111, 110, h2AperMag3Err, "hAperMag3Err"); + dPARSE (112, 111, h2AperMag4, "hAperMag4"); + dPARSE (113, 112, h2AperMag4Err, "hAperMag4Err"); + dPARSE (116, 113, h2PA, "hPA"); + iPARSE (119, 116, h2Class, "hClass"); + iPARSE (121, 119, h2ppErrBits, "hppErrBits"); + iPARSE (122, 121, h2SeqNum, "hSeqNum"); + + 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 (sourceID, uint64_t, NELEM); + REALLOCATE (frameSetID, uint64_t, NELEM); + REALLOCATE (ra, double, NELEM); + REALLOCATE (dec, double, NELEM); + REALLOCATE (priOrSec, uint64_t, NELEM); + REALLOCATE (mergedClass, int, NELEM); + REALLOCATE (pStar, double, NELEM); + REALLOCATE (pGalaxy, double, NELEM); + REALLOCATE (jAperMag1, double, NELEM); + REALLOCATE (jAperMag1Err, double, NELEM); + REALLOCATE (jAperMag3, double, NELEM); + REALLOCATE (jAperMag3Err, double, NELEM); + REALLOCATE (jAperMag4, double, NELEM); + REALLOCATE (jAperMag4Err, double, NELEM); + REALLOCATE (jPA, double, NELEM); + REALLOCATE (jClass, int, NELEM); + REALLOCATE (jppErrBits, int, NELEM); + REALLOCATE (jSeqNum, int, NELEM); + REALLOCATE (hAperMag1, double, NELEM); + REALLOCATE (hAperMag1Err, double, NELEM); + REALLOCATE (hAperMag3, double, NELEM); + REALLOCATE (hAperMag3Err, double, NELEM); + REALLOCATE (hAperMag4, double, NELEM); + REALLOCATE (hAperMag4Err, double, NELEM); + REALLOCATE (hPA, double, NELEM); + REALLOCATE (hClass, int, NELEM); + REALLOCATE (hppErrBits, int, NELEM); + REALLOCATE (hSeqNum, int, NELEM); + REALLOCATE (k1AperMag1, double, NELEM); + REALLOCATE (k1AperMag1Err, double, NELEM); + REALLOCATE (k1AperMag3, double, NELEM); + REALLOCATE (k1AperMag3Err, double, NELEM); + REALLOCATE (k1AperMag4, double, NELEM); + REALLOCATE (k1AperMag4Err, double, NELEM); + REALLOCATE (k1PA, double, NELEM); + REALLOCATE (k1Class, int, NELEM); + REALLOCATE (k1ppErrBits, int, NELEM); + REALLOCATE (k1SeqNum, int, NELEM); + REALLOCATE (k2AperMag1, double, NELEM); + REALLOCATE (k2AperMag1Err, double, NELEM); + REALLOCATE (k2AperMag3, double, NELEM); + REALLOCATE (k2AperMag3Err, double, NELEM); + REALLOCATE (k2AperMag4, double, NELEM); + REALLOCATE (k2AperMag4Err, double, NELEM); + REALLOCATE (k2PA, double, NELEM); + REALLOCATE (k2Class, int, NELEM); + REALLOCATE (k2ppErrBits, int, NELEM); + REALLOCATE (k2SeqNum, int, NELEM); + REALLOCATE (h2AperMag1, double, NELEM); + REALLOCATE (h2AperMag1Err, double, NELEM); + REALLOCATE (h2AperMag3, double, NELEM); + REALLOCATE (h2AperMag3Err, double, NELEM); + REALLOCATE (h2AperMag4, double, NELEM); + REALLOCATE (h2AperMag4Err, double, NELEM); + REALLOCATE (h2PA, double, NELEM); + REALLOCATE (h2Class, int, NELEM); + REALLOCATE (h2ppErrBits, int, NELEM); + REALLOCATE (h2SeqNum, 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; + + ALLOCATE_PTR (stars, UKIRT_Stars, NSTARS); + + for (int i = 0; i < NstarsIn; i++) { + + Rmin = MIN (Rmin, ra[i]); + Rmax = MAX (Rmax, ra[i]); + Dmin = MIN (Dmin, dec[i]); + Dmax = MAX (Dmax, dec[i]); + + // only allocate the measures for stars as we define them + ALLOCATE (stars[Nstars].measure, Measure, UKIRT_NFILTER); + + // we have NFILTER (== 6) UKIRT measurements per object + dvo_average_init (&stars[Nstars].average); + for (int j = 0; j < UKIRT_NFILTER; j++) { + dvo_measure_init (&stars[Nstars].measure[j]); + } + + stars[Nstars].average.R = ra[i]; + stars[Nstars].average.D = dec[i]; + stars[Nstars].average.dR = NAN; + stars[Nstars].average.dD = NAN; + + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + int isPrimary = (!priOrSec[i] || (priOrSec[i] == frameSetID[i])) ? ID_MEAS_STACK_PRIMARY : 0x00; + int photFlags2 = isPrimary; + switch (mergedClass[i]) { + case 1: photFlags2 |= 0x01; break; // Galaxy + case 0: photFlags2 |= 0x02; break; // Noise + case -1: photFlags2 |= 0x04; break; // Star + case -2: photFlags2 |= 0x08; break; // probably star + case -3: photFlags2 |= 0x10; break; // probably galaxy + case -9: photFlags2 |= 0x20; break; // saturated + default: break; + } + + // J measure + stars[Nstars].measure[0].extID = sourceID[i]; + stars[Nstars].measure[0].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[0].R = ra[i]; + stars[Nstars].measure[0].D = dec[i]; + stars[Nstars].measure[0].psfChisq = pStar[i]; + stars[Nstars].measure[0].extNsigma = pGalaxy[i]; + stars[Nstars].measure[0].M = jAperMag1[i]; + stars[Nstars].measure[0].dM = jAperMag1Err[i]; + stars[Nstars].measure[0].Map = jAperMag3[i]; + stars[Nstars].measure[0].dMap = jAperMag3Err[i]; + stars[Nstars].measure[0].Mkron = jAperMag4[i]; + stars[Nstars].measure[0].dMkron = jAperMag4Err[i]; + stars[Nstars].measure[0].posangle = jPA[i]; // XXX units + stars[Nstars].measure[0].psfQF = jClass[i]; + stars[Nstars].measure[0].photFlags = jppErrBits[i]; + stars[Nstars].measure[0].photFlags2 = photFlags2; + stars[Nstars].measure[0].detID = jSeqNum[i]; + stars[Nstars].measure[0].photcode = codeJ; + + // H measure + stars[Nstars].measure[1].extID = sourceID[i]; + stars[Nstars].measure[1].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[1].R = ra[i]; + stars[Nstars].measure[1].D = dec[i]; + stars[Nstars].measure[1].psfChisq = pStar[i]; + stars[Nstars].measure[1].extNsigma = pGalaxy[i]; + stars[Nstars].measure[1].M = hAperMag1[i]; + stars[Nstars].measure[1].dM = hAperMag1Err[i]; + stars[Nstars].measure[1].Map = hAperMag3[i]; + stars[Nstars].measure[1].dMap = hAperMag3Err[i]; + stars[Nstars].measure[1].Mkron = hAperMag4[i]; + stars[Nstars].measure[1].dMkron = hAperMag4Err[i]; + stars[Nstars].measure[1].posangle = hPA[i]; // XXX units + stars[Nstars].measure[1].psfQF = hClass[i]; + stars[Nstars].measure[1].photFlags = hppErrBits[i]; + stars[Nstars].measure[1].photFlags2 = photFlags2; + stars[Nstars].measure[1].detID = hSeqNum[i]; + stars[Nstars].measure[1].photcode = codeH; + + // K1 measure + stars[Nstars].measure[2].extID = sourceID[i]; + stars[Nstars].measure[2].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[2].R = ra[i]; + stars[Nstars].measure[2].D = dec[i]; + stars[Nstars].measure[2].psfChisq = pStar[i]; + stars[Nstars].measure[2].extNsigma = pGalaxy[i]; + stars[Nstars].measure[2].M = k1AperMag1[i]; + stars[Nstars].measure[2].dM = k1AperMag1Err[i]; + stars[Nstars].measure[2].Map = k1AperMag3[i]; + stars[Nstars].measure[2].dMap = k1AperMag3Err[i]; + stars[Nstars].measure[2].Mkron = k1AperMag4[i]; + stars[Nstars].measure[2].dMkron = k1AperMag4Err[i]; + stars[Nstars].measure[2].posangle = k1PA[i]; // XXX units + stars[Nstars].measure[2].psfQF = k1Class[i]; + stars[Nstars].measure[2].photFlags = k1ppErrBits[i]; + stars[Nstars].measure[2].photFlags2 = photFlags2; + stars[Nstars].measure[2].detID = k1SeqNum[i]; + stars[Nstars].measure[2].photcode = codeK; + + // K1 measure + stars[Nstars].measure[3].extID = sourceID[i]; + stars[Nstars].measure[3].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[3].R = ra[i]; + stars[Nstars].measure[3].D = dec[i]; + stars[Nstars].measure[3].psfChisq = pStar[i]; + stars[Nstars].measure[3].extNsigma = pGalaxy[i]; + stars[Nstars].measure[3].M = k2AperMag1[i]; + stars[Nstars].measure[3].dM = k2AperMag1Err[i]; + stars[Nstars].measure[3].Map = k2AperMag3[i]; + stars[Nstars].measure[3].dMap = k2AperMag3Err[i]; + stars[Nstars].measure[3].Mkron = k2AperMag4[i]; + stars[Nstars].measure[3].dMkron = k2AperMag4Err[i]; + stars[Nstars].measure[3].posangle = k2PA[i]; // XXX units + stars[Nstars].measure[3].psfQF = k2Class[i]; + stars[Nstars].measure[3].photFlags = k2ppErrBits[i]; + stars[Nstars].measure[3].photFlags2 = photFlags2; + stars[Nstars].measure[3].detID = k2SeqNum[i]; + stars[Nstars].measure[3].photcode = codeK; + + // H2 measure + stars[Nstars].measure[4].extID = sourceID[i]; + stars[Nstars].measure[4].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[4].R = ra[i]; + stars[Nstars].measure[4].D = dec[i]; + stars[Nstars].measure[4].psfChisq = pStar[i]; + stars[Nstars].measure[4].extNsigma = pGalaxy[i]; + stars[Nstars].measure[4].M = h2AperMag1[i]; + stars[Nstars].measure[4].dM = h2AperMag1Err[i]; + stars[Nstars].measure[4].Map = h2AperMag3[i]; + stars[Nstars].measure[4].dMap = h2AperMag3Err[i]; + stars[Nstars].measure[4].Mkron = h2AperMag4[i]; + stars[Nstars].measure[4].dMkron = h2AperMag4Err[i]; + stars[Nstars].measure[4].posangle = h2PA[i]; // XXX units + stars[Nstars].measure[4].psfQF = h2Class[i]; + stars[Nstars].measure[4].photFlags = h2ppErrBits[i]; + stars[Nstars].measure[4].photFlags2 = photFlags2; + stars[Nstars].measure[4].detID = h2SeqNum[i]; + stars[Nstars].measure[4].photcode = codeH; + + stars[Nstars].average.Nmeasure = UKIRT_NFILTER; + Nstars ++; + + if (NSTARS >= Nstars) { + NSTARS += 10000; + REALLOCATE (stars, UKIRT_Stars, NSTARS); + } + } + + FREE (sourceID); + FREE (frameSetID); + FREE (ra); + FREE (dec); + FREE (epoch); + FREE (priOrSec); + FREE (mergedClass); + FREE (pStar); + FREE (pGalaxy); + + FREE (jAperMag1); + FREE (jAperMag1Err); + FREE (jAperMag3); + FREE (jAperMag3Err); + FREE (jAperMag4); + FREE (jAperMag4Err); + FREE (jPA); + FREE (jClass); + FREE (jppErrBits); + FREE (jSeqNum); + + FREE (hAperMag1); + FREE (hAperMag1Err); + FREE (hAperMag3); + FREE (hAperMag3Err); + FREE (hAperMag4); + FREE (hAperMag4Err); + FREE (hPA); + FREE (hClass); + FREE (hppErrBits); + FREE (hSeqNum); + + FREE (k1AperMag1); + FREE (k1AperMag1Err); + FREE (k1AperMag3); + FREE (k1AperMag3Err); + FREE (k1AperMag4); + FREE (k1AperMag4Err); + FREE (k1PA); + FREE (k1Class); + FREE (k1ppErrBits); + FREE (k1SeqNum); + + FREE (k2AperMag1); + FREE (k2AperMag1Err); + FREE (k2AperMag3); + FREE (k2AperMag3Err); + FREE (k2AperMag4); + FREE (k2AperMag4Err); + FREE (k2PA); + FREE (k2Class); + FREE (k2ppErrBits); + FREE (k2SeqNum); + + FREE (h2AperMag1); + FREE (h2AperMag1Err); + FREE (h2AperMag3); + FREE (h2AperMag3Err); + FREE (h2AperMag4); + FREE (h2AperMag4Err); + FREE (h2PA); + FREE (h2Class); + FREE (h2ppErrBits); + FREE (h2SeqNum); + + *nstars = Nstars; + *nstart = Nstart; + return (stars); +} Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_uhs.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_uhs.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_uhs.c (revision 41151) @@ -0,0 +1,338 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +/* this function reads the values of interest from the UKIRT UHS CSV files: + + data model description: doc/mapping_UKIDSS_unWISE_DVO.pdf + + N -- field (column) number in CSV file (1 counting) + | UKIRT field name : DVO measure field + 1 sourceID : extID + 2 cuEventID : + 3 frameSetID : imageID (assumes we do not try to go backwards for these photcodes) + 4 ra : R + 5 dec : D + 6 cx : + 7 cy : + 8 cz : + 9 htmID : + 10 l : + 11 b : + 12 lambda : + 13 eta : + 14 priOrSec : primary bit in dbFlags + 15 mergedClassStat : + 16 mergedClass : psfQFperf or photFlags2 + 17 pStar : psfChisq + 18 pGalaxy : extNsigma + 19 pNoise : + 20 pSaturated : + 21 eBV : + 22 aJ : + 23 jHallMag : M + 24 jHallMagErr : dM + 25 jPetroMag : + 26 jPetroMagErr : + 27 jAperMag3 : Map + 28 jAperMag3Err : dMap + 29 jAperMag4 : Mkron + 30 jAperMag4Err : dMron + 31 jAperMag6 : + 32 jAperMag6Err : + 33 jGausig : + 34 jEll : + 35 jPA : posangle + 36 jErrBits : + 37 jDeblend : + 38 jClass : psfQF + 39 jClassStat : + 40 jppErrBits : photFlags + 41 jSeqNum : detID + +*/ + +# define iPARSE(NS,NE,FIELD,NAME) { \ + cA = iparse_csv_rpt (&ivalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = ivalue; } + +# define jPARSE(NS,NE,FIELD,NAME) { \ + cA = jparse_csv_rpt (&jvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = jvalue; } + +# define dPARSE(NS,NE,FIELD,NAME) { \ + cA = dparse_csv_rpt (&dvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = dvalue; } + +static int Nline_read = 0; // track number of lines read so far (use to skip lines as well) +static int Nskip = 1; // UKIRT UHS CSV files have a single header row (and no special character to mark) + +UKIRT_Stars *loadukirt_uhs_readstars_uhs (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars) { + + int codeJ = GetPhotcodeCodebyName ("UKIRT_J"); if (!codeJ) Shutdown ("missing photcode UKIRT_UHS_J"); + + // XXX I need a UKIRT UHS Mean Epoch + // time_t UKIRT_UHS_EPOCH = ohana_date_to_sec ("2016/01/01,00:00:00"); + // fprintf (stderr, "WARNING: using an invalid UKIRT_UHS_EPOCH (see loadukirt_uhs_readstars.c:60)\n"); + + int Nelem = 0; // number of valid rows read (vector elements) + int NELEM = 10000; // currently-allocated number of output rows + + // vectors to hold the data loaded from the csv file (names are UKIRT native names) + ALLOCATE_PTR (sourceID, uint64_t, NELEM); + ALLOCATE_PTR (frameSetID, uint64_t, NELEM); // XXX doc says 8-bytes for this value! + ALLOCATE_PTR (ra, double, NELEM); + ALLOCATE_PTR (dec, double, NELEM); + ALLOCATE_PTR (priOrSec, uint64_t, NELEM); + ALLOCATE_PTR (mergedClass, int, NELEM); // XXX doc says 2-bytes (could use a short) + ALLOCATE_PTR (pStar, double, NELEM); + ALLOCATE_PTR (pGalaxy, double, NELEM); + ALLOCATE_PTR (jHallMag, double, NELEM); + ALLOCATE_PTR (jHallMagErr, double, NELEM); + ALLOCATE_PTR (jAperMag3, double, NELEM); + ALLOCATE_PTR (jAperMag3Err, double, NELEM); + ALLOCATE_PTR (jAperMag4, double, NELEM); + ALLOCATE_PTR (jAperMag4Err, double, NELEM); + ALLOCATE_PTR (jPA, double, NELEM); + ALLOCATE_PTR (jClass, int, NELEM); + ALLOCATE_PTR (jppErrBits, int, NELEM); + ALLOCATE_PTR (jSeqNum, int, NELEM); + + // 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 = *nstart; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + while (!EndOfFile && (Nelem < NSTARS_MAX)) { + int Nbytes = BUFFER_SIZE - 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; } + + // for UKIRT UHS, 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 : TBD for UHS + + // the start of the line is the 1st element (fields are 1-counting) + jPARSE ( 1, 1, sourceID, "sourceID"); + iPARSE ( 3, 1, frameSetID, "frameSetID"); + dPARSE ( 4, 3, ra, "ra"); + dPARSE ( 5, 4, dec, "dec"); + jPARSE (14, 5, priOrSec, "priOrSec"); + iPARSE (16, 14, mergedClass, "mergedClass"); + dPARSE (17, 16, pStar, "pStar"); + dPARSE (18, 17, pGalaxy, "pGalaxy"); + dPARSE (23, 18, jHallMag, "jHallMag"); + dPARSE (24, 23, jHallMagErr, "jHallMagErr"); + dPARSE (27, 24, jAperMag3, "jAperMag3"); + dPARSE (28, 27, jAperMag3Err, "jAperMag3Err"); + dPARSE (29, 28, jAperMag4, "jAperMag4"); + dPARSE (30, 29, jAperMag4Err, "jAperMag4Err"); + dPARSE (35, 30, jPA, "jPA"); + iPARSE (38, 35, jClass, "jClass"); + iPARSE (40, 38, jppErrBits, "jppErrBits"); + iPARSE (41, 40, jSeqNum, "jSeqNum"); + + 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 (sourceID, uint64_t, NELEM); + REALLOCATE (frameSetID, uint64_t, NELEM); + REALLOCATE (ra, double, NELEM); + REALLOCATE (dec, double, NELEM); + REALLOCATE (priOrSec, uint64_t, NELEM); + REALLOCATE (mergedClass, int, NELEM); + REALLOCATE (pStar, double, NELEM); + REALLOCATE (pGalaxy, double, NELEM); + REALLOCATE (jHallMag, double, NELEM); + REALLOCATE (jHallMagErr, double, NELEM); + REALLOCATE (jAperMag3, double, NELEM); + REALLOCATE (jAperMag3Err, double, NELEM); + REALLOCATE (jAperMag4, double, NELEM); + REALLOCATE (jAperMag4Err, double, NELEM); + REALLOCATE (jPA, double, NELEM); + REALLOCATE (jClass, int, NELEM); + REALLOCATE (jppErrBits, int, NELEM); + REALLOCATE (jSeqNum, 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; + + ALLOCATE_PTR (stars, UKIRT_Stars, NSTARS); + + for (int i = 0; i < NstarsIn; i++) { + + Rmin = MIN (Rmin, ra[i]); + Rmax = MAX (Rmax, ra[i]); + Dmin = MIN (Dmin, dec[i]); + Dmax = MAX (Dmax, dec[i]); + + // only allocate the measures for stars as we define them + ALLOCATE (stars[Nstars].measure, Measure, UKIRT_NFILTER); + + // we have one UKIRT UHS (J-band) measurement per object + dvo_average_init (&stars[Nstars].average); + for (int j = 0; j < UKIRT_NFILTER; j++) { + dvo_measure_init (&stars[Nstars].measure[j]); + } + + stars[Nstars].average.R = ra[i]; + stars[Nstars].average.D = dec[i]; + stars[Nstars].average.dR = NAN; + stars[Nstars].average.dD = NAN; + + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + int isPrimary = (!priOrSec[i] || (priOrSec[i] == frameSetID[i])) ? ID_MEAS_STACK_PRIMARY : 0x00; + int photFlags2 = isPrimary; + switch (mergedClass[i]) { + case 1: photFlags2 |= 0x01; break; // Galaxy + case 0: photFlags2 |= 0x02; break; // Noise + case -1: photFlags2 |= 0x04; break; // Star + case -2: photFlags2 |= 0x08; break; // probably star + case -3: photFlags2 |= 0x10; break; // probably galaxy + case -9: photFlags2 |= 0x20; break; // saturated + default: break; + } + + stars[Nstars].measure[0].extID = sourceID[i]; + stars[Nstars].measure[0].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[0].R = ra[i]; + stars[Nstars].measure[0].D = dec[i]; + stars[Nstars].measure[0].psfChisq = pStar[i]; + stars[Nstars].measure[0].extNsigma = pGalaxy[i]; + stars[Nstars].measure[0].M = jHallMag[i]; + stars[Nstars].measure[0].dM = jHallMagErr[i]; + stars[Nstars].measure[0].Map = jAperMag3[i]; + stars[Nstars].measure[0].dMap = jAperMag3Err[i]; + stars[Nstars].measure[0].Mkron = jAperMag4[i]; + stars[Nstars].measure[0].dMkron = jAperMag4Err[i]; + stars[Nstars].measure[0].posangle = jPA[i]; // XXX units + stars[Nstars].measure[0].psfQF = jClass[i]; + stars[Nstars].measure[0].photFlags = jppErrBits[i]; + stars[Nstars].measure[0].photFlags2 = photFlags2; + stars[Nstars].measure[0].detID = jSeqNum[i]; + stars[Nstars].measure[0].photcode = codeJ; + + stars[Nstars].average.Nmeasure = UKIRT_NFILTER; + Nstars ++; + + if (NSTARS >= Nstars) { + NSTARS += 10000; + REALLOCATE (stars, UKIRT_Stars, NSTARS); + } + } + + FREE (sourceID); + FREE (frameSetID); + FREE (ra); + FREE (dec); + FREE (priOrSec); + FREE (mergedClass); + FREE (pStar); + FREE (pGalaxy); + FREE (jHallMag); + FREE (jHallMagErr); + FREE (jAperMag3); + FREE (jAperMag3Err); + FREE (jAperMag4); + FREE (jAperMag4Err); + FREE (jPA); + FREE (jClass); + FREE (jppErrBits); + FREE (jSeqNum); + + *nstars = Nstars; + *nstart = Nstart; + return (stars); +} + Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ulas.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ulas.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_readstars_ulas.c (revision 41151) @@ -0,0 +1,822 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +/* this function reads the values of interest from the UKIRT UGCS CSV files: + + data model description: doc/mapping_UKIDSS_unWISE_DVO.pdf + + N -- field (column) number in CSV file (1 counting) + | UKIRT field name : DVO measure field + + 1 sourceID : extID + 2 frameSetID : imageID (assumes we do not try to go backwards for these photcodes) + 3 ra : R + 4 dec : D + 5 sigRa : + 6 sigDec : + 7 epoch : t ???? -- how is there a single epoch? + 8 muRa : + 9 muDec : + 10 sigMuRa : + 11 sigMuDec : + 12 chi2 : + 13 nFrames : + 14 cx : + 15 cy : + 16 cz : + 17 htmID : + 18 l : + 19 b : + 20 lambda : + 21 eta : + 22 priOrSec : primary bit in dbFlags + 23 ymj_1Pnt : + 24 ymj_1PntErr : + 25 j_1mhPnt : + 26 j_1mhPntErr : + 27 hmkPnt : + 28 hmkPntErr : + 29 ymj_1Ext : + 30 ymj_1ExtErr : + 31 j_1mhExt : + 32 j_1mhExtErr : + 33 hmkExt : + 34 hmkExtErr : + 35 mergedClassStat : + 36 mergedClass : psfQFperf or photFlags2 + 37 pStar : psfChisq + 38 pGalaxy : extNsigma + 39 pNoise : + 40 pSaturated : + 41 eBV : + 42 aY : + 43 aJ + 44 aH + 45 aK + 46 yHallMag : M + 47 yHallMagErr : dM + 48 yPetroMag : + 49 yPetroMagErr : + 50 yAperMag3 : Map + 51 yAperMag3Err : dMap + 52 yAperMag4 : Mkron + 53 yAperMag4Err : dMron + 54 yAperMag6 : + 55 yAperMag6Err : + 56 yGausig : + 57 yEll : + 58 yPA : posangle + 59 yErrBits : + 60 yDeblend : + 61 yClass : psfQF + 62 yClassStat : + 63 yppErrBits : photFlags + 64 ySeqNum : detID + 65 yObjID : + 66 yXi : + 67 yEta : + 68 j_1HallMag : M + 69 j_1HallMagErr : dM + 70 j_1PetroMag : + 71 j_1PetroMagErr : + 72 j_1AperMag3 : Map + 73 j_1AperMag3Err : dMap + 74 j_1AperMag4 : Mkron + 75 j_1AperMag4Err : dMron + 76 j_1AperMag6 : + 77 j_1AperMag6Err : + 78 j_1Gausig : + 79 j_1Ell : + 80 j_1PA : posangle + 81 j_1ErrBits : + 82 j_1Deblend : + 83 j_1Class : psfQF + 84 j_1ClassStat : + 85 j_1ppErrBits : photFlags + 86 j_1SeqNum : detID + 87 j_1ObjID : + 88 j_1Xi : + 89 j_1Eta : + 90 j_2HallMag : M + 91 j_2HallMagErr : dM + 92 j_2PetroMag : + 93 j_2PetroMagErr : + 94 j_2AperMag3 : Map + 95 j_2AperMag3Err : dMap + 96 j_2AperMag4 : Mkron + 97 j_2AperMag4Err : dMron + 98 j_2AperMag6 : + 99 j_2AperMag6Err : +100 j_2Gausig : +101 j_2Ell : +102 j_2PA : posangle +103 j_2ErrBits : +104 j_2Deblend : +105 j_2Class : psfQF +106 j_2ClassStat : +107 j_2ppErrBits : photFlags +108 j_2SeqNum : detID +109 j_2ObjID : +110 j_2Xi : +111 j_2Eta : +112 hHallMag : M +113 hHallMagErr : dM +114 hPetroMag : +115 hPetroMagErr : +116 hAperMag3 : Map +117 hAperMag3Err : dMap +118 hAperMag4 : Mkron +119 hAperMag4Err : dMron +120 hAperMag6 : +121 hAperMag6Err : +122 hGausig : +123 hEll : +124 hPA : posangle +125 hErrBits : +126 hDeblend : +127 hClass : psfQF +128 hClassStat : +129 hppErrBits : photFlags +130 hSeqNum : detID +131 hObjID : +132 hXi : +133 hEta : +134 kHallMag : M +135 kHallMagErr : dM +136 kPetroMag : +137 kPetroMagErr : +138 kAperMag3 : Map +139 kAperMag3Err : dMap +140 kAperMag4 : Mkron +141 kAperMag4Err : dMron +142 kAperMag6 : +143 kAperMag6Err : +144 kGausig : +145 kEll : +146 kPA : posangle +147 kErrBits : +148 kDeblend : +149 kClass : psfQF +150 kClassStat : +151 kppErrBits : photFlags +152 kSeqNum : detID +153 kObjID : +154 kXi : +155 kEta : + +*/ + +# define iPARSE(NS,NE,FIELD,NAME) { \ + cA = iparse_csv_rpt (&ivalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = ivalue; } + +# define jPARSE(NS,NE,FIELD,NAME) { \ + cA = jparse_csv_rpt (&jvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = jvalue; } + +# define dPARSE(NS,NE,FIELD,NAME) { \ + cA = dparse_csv_rpt (&dvalue, (NS), (NE), cA, &readStatus); \ + if (!readStatus && VERBOSE) { \ + gprint (GP_ERR, "suspect field: %d (%s) in %s\n", (NS), NAME, c0); \ + } \ + lineStatus &= readStatus; FIELD[Nelem] = dvalue; } + +static int Nline_read = 0; // track number of lines read so far (use to skip lines as well) +static int Nskip = 0; // UKIRT UHS CSV files have a single header row (and no special character to mark) + +UKIRT_Stars *loadukirt_uhs_readstars_ulas (FILE *f, char *buffer, int *nstart, AddstarClientOptions *options, int *nstars) { + + int codeY = GetPhotcodeCodebyName ("UKIRT_Y"); if (!codeY) Shutdown ("missing photcode UKIRT_Y"); + int codeJ = GetPhotcodeCodebyName ("UKIRT_J"); if (!codeJ) Shutdown ("missing photcode UKIRT_J"); + int codeH = GetPhotcodeCodebyName ("UKIRT_H"); if (!codeH) Shutdown ("missing photcode UKIRT_H"); + int codeK = GetPhotcodeCodebyName ("UKIRT_K"); if (!codeK) Shutdown ("missing photcode UKIRT_K"); + + // XXX I need a UKIRT UHS Mean Epoch + // time_t UKIRT_UHS_EPOCH = ohana_date_to_sec ("2016/01/01,00:00:00"); + // fprintf (stderr, "WARNING: using an invalid UKIRT_UHS_EPOCH (see loadukirt_uhs_readstars.c:60)\n"); + + int Nelem = 0; // number of valid rows read (vector elements) + int NELEM = 10000; // currently-allocated number of output rows + + // vectors to hold the data loaded from the csv file (names are UKIRT native names) + ALLOCATE_PTR (sourceID, uint64_t, NELEM); + ALLOCATE_PTR (frameSetID, uint64_t, NELEM); // XXX doc says 8-bytes for this value! + ALLOCATE_PTR (ra, double, NELEM); + ALLOCATE_PTR (dec, double, NELEM); + ALLOCATE_PTR (epoch, double, NELEM); + ALLOCATE_PTR (priOrSec, uint64_t, NELEM); + ALLOCATE_PTR (mergedClass, int, NELEM); // XXX doc says 2-bytes (could use a short) + ALLOCATE_PTR (pStar, double, NELEM); + ALLOCATE_PTR (pGalaxy, double, NELEM); + + ALLOCATE_PTR (yHallMag, double, NELEM); + ALLOCATE_PTR (yHallMagErr, double, NELEM); + ALLOCATE_PTR (yAperMag3, double, NELEM); + ALLOCATE_PTR (yAperMag3Err, double, NELEM); + ALLOCATE_PTR (yAperMag4, double, NELEM); + ALLOCATE_PTR (yAperMag4Err, double, NELEM); + ALLOCATE_PTR (yPA, double, NELEM); + ALLOCATE_PTR (yClass, int, NELEM); + ALLOCATE_PTR (yppErrBits, int, NELEM); + ALLOCATE_PTR (ySeqNum, int, NELEM); + + ALLOCATE_PTR (j1HallMag, double, NELEM); + ALLOCATE_PTR (j1HallMagErr, double, NELEM); + ALLOCATE_PTR (j1AperMag3, double, NELEM); + ALLOCATE_PTR (j1AperMag3Err, double, NELEM); + ALLOCATE_PTR (j1AperMag4, double, NELEM); + ALLOCATE_PTR (j1AperMag4Err, double, NELEM); + ALLOCATE_PTR (j1PA, double, NELEM); + ALLOCATE_PTR (j1Class, int, NELEM); + ALLOCATE_PTR (j1ppErrBits, int, NELEM); + ALLOCATE_PTR (j1SeqNum, int, NELEM); + + ALLOCATE_PTR (j2HallMag, double, NELEM); + ALLOCATE_PTR (j2HallMagErr, double, NELEM); + ALLOCATE_PTR (j2AperMag3, double, NELEM); + ALLOCATE_PTR (j2AperMag3Err, double, NELEM); + ALLOCATE_PTR (j2AperMag4, double, NELEM); + ALLOCATE_PTR (j2AperMag4Err, double, NELEM); + ALLOCATE_PTR (j2PA, double, NELEM); + ALLOCATE_PTR (j2Class, int, NELEM); + ALLOCATE_PTR (j2ppErrBits, int, NELEM); + ALLOCATE_PTR (j2SeqNum, int, NELEM); + + ALLOCATE_PTR (hHallMag, double, NELEM); + ALLOCATE_PTR (hHallMagErr, double, NELEM); + ALLOCATE_PTR (hAperMag3, double, NELEM); + ALLOCATE_PTR (hAperMag3Err, double, NELEM); + ALLOCATE_PTR (hAperMag4, double, NELEM); + ALLOCATE_PTR (hAperMag4Err, double, NELEM); + ALLOCATE_PTR (hPA, double, NELEM); + ALLOCATE_PTR (hClass, int, NELEM); + ALLOCATE_PTR (hppErrBits, int, NELEM); + ALLOCATE_PTR (hSeqNum, int, NELEM); + + ALLOCATE_PTR (kHallMag, double, NELEM); + ALLOCATE_PTR (kHallMagErr, double, NELEM); + ALLOCATE_PTR (kAperMag3, double, NELEM); + ALLOCATE_PTR (kAperMag3Err, double, NELEM); + ALLOCATE_PTR (kAperMag4, double, NELEM); + ALLOCATE_PTR (kAperMag4Err, double, NELEM); + ALLOCATE_PTR (kPA, double, NELEM); + ALLOCATE_PTR (kClass, int, NELEM); + ALLOCATE_PTR (kppErrBits, int, NELEM); + ALLOCATE_PTR (kSeqNum, int, NELEM); + + // 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 = *nstart; // location of the last valid byte in the buffer (start filling here) + int EndOfFile = FALSE; + while (!EndOfFile && (Nelem < NSTARS_MAX)) { + int Nbytes = BUFFER_SIZE - 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; } + + // for UKIRT UHS, 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 : TBD for UHS + + // the start of the line is the 1st element (fields are 1-counting) + jPARSE ( 1, 1, sourceID, "sourceID"); + iPARSE ( 2, 1, frameSetID, "frameSetID"); + dPARSE ( 3, 2, ra, "ra"); + dPARSE ( 4, 3, dec, "dec"); + dPARSE ( 7, 4, epoch, "epoch"); + jPARSE ( 22, 7, priOrSec, "priOrSec"); // watch for epoch + iPARSE ( 36, 22, mergedClass, "mergedClass"); + dPARSE ( 37, 36, pStar, "pStar"); + dPARSE ( 38, 37, pGalaxy, "pGalaxy"); + dPARSE ( 46, 38, yHallMag, "yHallMag"); + dPARSE ( 47, 46, yHallMagErr, "yHallMagErr"); + dPARSE ( 50, 47, yAperMag3, "yAperMag3"); + dPARSE ( 51, 50, yAperMag3Err, "yAperMag3Err"); + dPARSE ( 52, 51, yAperMag4, "yAperMag4"); + dPARSE ( 53, 52, yAperMag4Err, "yAperMag4Err"); + dPARSE ( 58, 53, yPA, "yPA"); + iPARSE ( 61, 58, yClass, "yClass"); + iPARSE ( 63, 61, yppErrBits, "yppErrBits"); + iPARSE ( 64, 63, ySeqNum, "ySeqNum"); + dPARSE ( 68, 64, j1HallMag, "j1HallMag"); + dPARSE ( 69, 68, j1HallMagErr, "j1HallMagErr"); + dPARSE ( 72, 69, j1AperMag3, "j1AperMag3"); + dPARSE ( 73, 72, j1AperMag3Err, "j1AperMag3Err"); + dPARSE ( 74, 73, j1AperMag4, "j1AperMag4"); + dPARSE ( 75, 74, j1AperMag4Err, "j1AperMag4Err"); + dPARSE ( 80, 75, j1PA, "j1PA"); + iPARSE ( 83, 80, j1Class, "j1Class"); + iPARSE ( 85, 83, j1ppErrBits, "j1ppErrBits"); + iPARSE ( 86, 85, j1SeqNum, "j1SeqNum"); + dPARSE ( 90, 86, j2HallMag, "j2HallMag"); + dPARSE ( 91, 90, j2HallMagErr, "j2HallMagErr"); + dPARSE ( 94, 91, j2AperMag3, "j2AperMag3"); + dPARSE ( 95, 94, j2AperMag3Err, "j2AperMag3Err"); + dPARSE ( 96, 95, j2AperMag4, "j2AperMag4"); + dPARSE ( 97, 96, j2AperMag4Err, "j2AperMag4Err"); + dPARSE (102, 97, j2PA, "j2PA"); + iPARSE (105, 102, j2Class, "j2Class"); + iPARSE (107, 105, j2ppErrBits, "j2ppErrBits"); + iPARSE (108, 107, j2SeqNum, "j2SeqNum"); + dPARSE (112, 108, hHallMag, "hHallMag"); + dPARSE (113, 112, hHallMagErr, "hHallMagErr"); + dPARSE (116, 113, hAperMag3, "hAperMag3"); + dPARSE (117, 116, hAperMag3Err, "hAperMag3Err"); + dPARSE (118, 117, hAperMag4, "hAperMag4"); + dPARSE (119, 118, hAperMag4Err, "hAperMag4Err"); + dPARSE (124, 119, hPA, "hPA"); + iPARSE (127, 124, hClass, "hClass"); + iPARSE (129, 127, hppErrBits, "hppErrBits"); + iPARSE (130, 129, hSeqNum, "hSeqNum"); + dPARSE (134, 130, kHallMag, "kHallMag"); + dPARSE (135, 134, kHallMagErr, "kHallMagErr"); + dPARSE (138, 135, kAperMag3, "kAperMag3"); + dPARSE (139, 138, kAperMag3Err, "kAperMag3Err"); + dPARSE (140, 139, kAperMag4, "kAperMag4"); + dPARSE (141, 140, kAperMag4Err, "kAperMag4Err"); + dPARSE (146, 141, kPA, "kPA"); + iPARSE (149, 146, kClass, "kClass"); + iPARSE (151, 149, kppErrBits, "kppErrBits"); + iPARSE (152, 151, kSeqNum, "kSeqNum"); + + 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 (sourceID, uint64_t, NELEM); + REALLOCATE (frameSetID, uint64_t, NELEM); + REALLOCATE (ra, double, NELEM); + REALLOCATE (dec, double, NELEM); + REALLOCATE (priOrSec, uint64_t, NELEM); + REALLOCATE (mergedClass, int, NELEM); + REALLOCATE (pStar, double, NELEM); + REALLOCATE (pGalaxy, double, NELEM); + REALLOCATE (yHallMag, double, NELEM); + REALLOCATE (yHallMagErr, double, NELEM); + REALLOCATE (yAperMag3, double, NELEM); + REALLOCATE (yAperMag3Err, double, NELEM); + REALLOCATE (yAperMag4, double, NELEM); + REALLOCATE (yAperMag4Err, double, NELEM); + REALLOCATE (yPA, double, NELEM); + REALLOCATE (yClass, int, NELEM); + REALLOCATE (yppErrBits, int, NELEM); + REALLOCATE (ySeqNum, int, NELEM); + REALLOCATE (j1HallMag, double, NELEM); + REALLOCATE (j1HallMagErr, double, NELEM); + REALLOCATE (j1AperMag3, double, NELEM); + REALLOCATE (j1AperMag3Err, double, NELEM); + REALLOCATE (j1AperMag4, double, NELEM); + REALLOCATE (j1AperMag4Err, double, NELEM); + REALLOCATE (j1PA, double, NELEM); + REALLOCATE (j1Class, int, NELEM); + REALLOCATE (j1ppErrBits, int, NELEM); + REALLOCATE (j1SeqNum, int, NELEM); + REALLOCATE (j2HallMag, double, NELEM); + REALLOCATE (j2HallMagErr, double, NELEM); + REALLOCATE (j2AperMag3, double, NELEM); + REALLOCATE (j2AperMag3Err, double, NELEM); + REALLOCATE (j2AperMag4, double, NELEM); + REALLOCATE (j2AperMag4Err, double, NELEM); + REALLOCATE (j2PA, double, NELEM); + REALLOCATE (j2Class, int, NELEM); + REALLOCATE (j2ppErrBits, int, NELEM); + REALLOCATE (j2SeqNum, int, NELEM); + REALLOCATE (hHallMag, double, NELEM); + REALLOCATE (hHallMagErr, double, NELEM); + REALLOCATE (hAperMag3, double, NELEM); + REALLOCATE (hAperMag3Err, double, NELEM); + REALLOCATE (hAperMag4, double, NELEM); + REALLOCATE (hAperMag4Err, double, NELEM); + REALLOCATE (hPA, double, NELEM); + REALLOCATE (hClass, int, NELEM); + REALLOCATE (hppErrBits, int, NELEM); + REALLOCATE (hSeqNum, int, NELEM); + REALLOCATE (kHallMag, double, NELEM); + REALLOCATE (kHallMagErr, double, NELEM); + REALLOCATE (kAperMag3, double, NELEM); + REALLOCATE (kAperMag3Err, double, NELEM); + REALLOCATE (kAperMag4, double, NELEM); + REALLOCATE (kAperMag4Err, double, NELEM); + REALLOCATE (kPA, double, NELEM); + REALLOCATE (kClass, int, NELEM); + REALLOCATE (kppErrBits, int, NELEM); + REALLOCATE (kSeqNum, 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; + + ALLOCATE_PTR (stars, UKIRT_Stars, NSTARS); + + for (int i = 0; i < NstarsIn; i++) { + + Rmin = MIN (Rmin, ra[i]); + Rmax = MAX (Rmax, ra[i]); + Dmin = MIN (Dmin, dec[i]); + Dmax = MAX (Dmax, dec[i]); + + // only allocate the measures for stars as we define them + ALLOCATE (stars[Nstars].measure, Measure, UKIRT_NFILTER); + + // we have NFILTER (== 6) UKIRT measurements per object + dvo_average_init (&stars[Nstars].average); + for (int j = 0; j < UKIRT_NFILTER; j++) { + dvo_measure_init (&stars[Nstars].measure[j]); + } + + stars[Nstars].average.R = ra[i]; + stars[Nstars].average.D = dec[i]; + stars[Nstars].average.dR = NAN; + stars[Nstars].average.dD = NAN; + + stars[Nstars].flag = FALSE; + stars[Nstars].found = FALSE; + + int isPrimary = (!priOrSec[i] || (priOrSec[i] == frameSetID[i])) ? ID_MEAS_STACK_PRIMARY : 0x00; + int photFlags2 = isPrimary; + switch (mergedClass[i]) { + case 1: photFlags2 |= 0x01; break; // Galaxy + case 0: photFlags2 |= 0x02; break; // Noise + case -1: photFlags2 |= 0x04; break; // Star + case -2: photFlags2 |= 0x08; break; // probably star + case -3: photFlags2 |= 0x10; break; // probably galaxy + case -9: photFlags2 |= 0x20; break; // saturated + default: break; + } + + // y measure + stars[Nstars].measure[0].extID = sourceID[i]; + stars[Nstars].measure[0].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[0].R = ra[i]; + stars[Nstars].measure[0].D = dec[i]; + stars[Nstars].measure[0].psfChisq = pStar[i]; + stars[Nstars].measure[0].extNsigma = pGalaxy[i]; + stars[Nstars].measure[0].M = yHallMag[i]; + stars[Nstars].measure[0].dM = yHallMagErr[i]; + stars[Nstars].measure[0].Map = yAperMag3[i]; + stars[Nstars].measure[0].dMap = yAperMag3Err[i]; + stars[Nstars].measure[0].Mkron = yAperMag4[i]; + stars[Nstars].measure[0].dMkron = yAperMag4Err[i]; + stars[Nstars].measure[0].posangle = yPA[i]; // XXX units + stars[Nstars].measure[0].psfQF = yClass[i]; + stars[Nstars].measure[0].photFlags = yppErrBits[i]; + stars[Nstars].measure[0].photFlags2 = photFlags2; + stars[Nstars].measure[0].detID = ySeqNum[i]; + stars[Nstars].measure[0].photcode = codeY; + + // J_1 measure + stars[Nstars].measure[1].extID = sourceID[i]; + stars[Nstars].measure[1].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[1].R = ra[i]; + stars[Nstars].measure[1].D = dec[i]; + stars[Nstars].measure[1].psfChisq = pStar[i]; + stars[Nstars].measure[1].extNsigma = pGalaxy[i]; + stars[Nstars].measure[1].M = j1HallMag[i]; + stars[Nstars].measure[1].dM = j1HallMagErr[i]; + stars[Nstars].measure[1].Map = j1AperMag3[i]; + stars[Nstars].measure[1].dMap = j1AperMag3Err[i]; + stars[Nstars].measure[1].Mkron = j1AperMag4[i]; + stars[Nstars].measure[1].dMkron = j1AperMag4Err[i]; + stars[Nstars].measure[1].posangle = j1PA[i]; // XXX units + stars[Nstars].measure[1].psfQF = j1Class[i]; + stars[Nstars].measure[1].photFlags = j1ppErrBits[i]; + stars[Nstars].measure[1].photFlags2 = photFlags2; + stars[Nstars].measure[1].detID = j1SeqNum[i]; + stars[Nstars].measure[1].photcode = codeJ; + + // J measure + stars[Nstars].measure[2].extID = sourceID[i]; + stars[Nstars].measure[2].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[2].R = ra[i]; + stars[Nstars].measure[2].D = dec[i]; + stars[Nstars].measure[2].psfChisq = pStar[i]; + stars[Nstars].measure[2].extNsigma = pGalaxy[i]; + stars[Nstars].measure[2].M = j2HallMag[i]; + stars[Nstars].measure[2].dM = j2HallMagErr[i]; + stars[Nstars].measure[2].Map = j2AperMag3[i]; + stars[Nstars].measure[2].dMap = j2AperMag3Err[i]; + stars[Nstars].measure[2].Mkron = j2AperMag4[i]; + stars[Nstars].measure[2].dMkron = j2AperMag4Err[i]; + stars[Nstars].measure[2].posangle = j2PA[i]; // XXX units + stars[Nstars].measure[2].psfQF = j2Class[i]; + stars[Nstars].measure[2].photFlags = j2ppErrBits[i]; + stars[Nstars].measure[2].photFlags2 = photFlags2; + stars[Nstars].measure[2].detID = j2SeqNum[i]; + stars[Nstars].measure[2].photcode = codeJ; + + // H measure + stars[Nstars].measure[3].extID = sourceID[i]; + stars[Nstars].measure[3].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[3].R = ra[i]; + stars[Nstars].measure[3].D = dec[i]; + stars[Nstars].measure[3].psfChisq = pStar[i]; + stars[Nstars].measure[3].extNsigma = pGalaxy[i]; + stars[Nstars].measure[3].M = hHallMag[i]; + stars[Nstars].measure[3].dM = hHallMagErr[i]; + stars[Nstars].measure[3].Map = hAperMag3[i]; + stars[Nstars].measure[3].dMap = hAperMag3Err[i]; + stars[Nstars].measure[3].Mkron = hAperMag4[i]; + stars[Nstars].measure[3].dMkron = hAperMag4Err[i]; + stars[Nstars].measure[3].posangle = hPA[i]; // XXX units + stars[Nstars].measure[3].psfQF = hClass[i]; + stars[Nstars].measure[3].photFlags = hppErrBits[i]; + stars[Nstars].measure[3].photFlags2 = photFlags2; + stars[Nstars].measure[3].detID = hSeqNum[i]; + stars[Nstars].measure[3].photcode = codeH; + + // K measure + stars[Nstars].measure[4].extID = sourceID[i]; + stars[Nstars].measure[4].imageID = frameSetID[i]; // XXX check for frameSetID > 31-bit int + stars[Nstars].measure[4].R = ra[i]; + stars[Nstars].measure[4].D = dec[i]; + stars[Nstars].measure[4].psfChisq = pStar[i]; + stars[Nstars].measure[4].extNsigma = pGalaxy[i]; + stars[Nstars].measure[4].M = kHallMag[i]; + stars[Nstars].measure[4].dM = kHallMagErr[i]; + stars[Nstars].measure[4].Map = kAperMag3[i]; + stars[Nstars].measure[4].dMap = kAperMag3Err[i]; + stars[Nstars].measure[4].Mkron = kAperMag4[i]; + stars[Nstars].measure[4].dMkron = kAperMag4Err[i]; + stars[Nstars].measure[4].posangle = kPA[i]; // XXX units + stars[Nstars].measure[4].psfQF = kClass[i]; + stars[Nstars].measure[4].photFlags = kppErrBits[i]; + stars[Nstars].measure[4].photFlags2 = photFlags2; + stars[Nstars].measure[4].detID = kSeqNum[i]; + stars[Nstars].measure[4].photcode = codeK; + + stars[Nstars].average.Nmeasure = UKIRT_NFILTER; + Nstars ++; + + if (NSTARS >= Nstars) { + NSTARS += 10000; + REALLOCATE (stars, UKIRT_Stars, NSTARS); + } + } + + FREE (sourceID); + FREE (frameSetID); + FREE (ra); + FREE (dec); + FREE (epoch); + FREE (priOrSec); + FREE (mergedClass); + FREE (pStar); + FREE (pGalaxy); + + FREE (yHallMag); + FREE (yHallMagErr); + FREE (yAperMag3); + FREE (yAperMag3Err); + FREE (yAperMag4); + FREE (yAperMag4Err); + FREE (yPA); + FREE (yClass); + FREE (yppErrBits); + FREE (ySeqNum); + + FREE (j1HallMag); + FREE (j1HallMagErr); + FREE (j1AperMag3); + FREE (j1AperMag3Err); + FREE (j1AperMag4); + FREE (j1AperMag4Err); + FREE (j1PA); + FREE (j1Class); + FREE (j1ppErrBits); + FREE (j1SeqNum); + + FREE (j2HallMag); + FREE (j2HallMagErr); + FREE (j2AperMag3); + FREE (j2AperMag3Err); + FREE (j2AperMag4); + FREE (j2AperMag4Err); + FREE (j2PA); + FREE (j2Class); + FREE (j2ppErrBits); + FREE (j2SeqNum); + + FREE (hHallMag); + FREE (hHallMagErr); + FREE (hAperMag3); + FREE (hAperMag3Err); + FREE (hAperMag4); + FREE (hAperMag4Err); + FREE (hPA); + FREE (hClass); + FREE (hppErrBits); + FREE (hSeqNum); + + FREE (kHallMag); + FREE (kHallMagErr); + FREE (kAperMag3); + FREE (kAperMag3Err); + FREE (kAperMag4); + FREE (kAperMag4Err); + FREE (kPA); + FREE (kClass); + FREE (kppErrBits); + FREE (kSeqNum); + + *nstars = Nstars; + *nstart = Nstart; + return (stars); +} + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Index: trunk/Ohana/src/addstar/src/loadukirt_uhs_table.c =================================================================== --- trunk/Ohana/src/addstar/src/loadukirt_uhs_table.c (revision 41151) +++ trunk/Ohana/src/addstar/src/loadukirt_uhs_table.c (revision 41151) @@ -0,0 +1,90 @@ +# include "addstar.h" +# include "ukirt_uhs.h" + +// one pass reads a single file at a time +int loadukirt_uhs_table (SkyList *skylistInput, char *filename, AddstarClientOptions *options) { + + // open the input file here and pass the pointer below + FILE *f = fopen (filename, "r"); + if (f == NULL) Shutdown ("can't read ukirt_uhs file: %s", filename); + + // we allocate one extra byte into which we never read so there will always be a NULL terminating the string + ALLOCATE_PTR (buffer, char, BUFFER_SIZE + 1); + bzero (buffer, BUFFER_SIZE + 1); + + // starting point of valid data in buffer; this is updated on each pass to readstars + int Nstart = 0; + + while (1) { + + // read the chunk from the file + fprintf (stderr, "loading %s\n", filename); + + // on each pass, we read a new chunk of data and receive Nstars + + int Nstars = 0; + UKIRT_Stars *stars = loadukirt_uhs_readstars (f, buffer, &Nstart, options, &Nstars); + if (!Nstars) { + fclose (f); + free (buffer); + free (stars); + return TRUE; // end of file reached + } + + fprintf (stderr, "writing %d stars to dvo\n", Nstars); + + // sort the stars by RA (Nmeasure agnostic) + loadukirt_uhs_sortStars (stars, Nstars); + + // scan through the stars, loading the containing catalogs + // skip through table for unsaved stars + for (int 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 (Nmeasure agnostic) + int Nsubset; + UKIRT_Stars *subset = loadukirt_uhs_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: + // loadukirt_uhs_save_remote (subset, Nsubset, hosts, region, skylist[0].filename[0], options); + // loadukirt_uhs_catalog (Nmeasure agnostic) + loadukirt_uhs_catalog (subset, Nsubset, region, skylist[0].filename[0], options); + free (subset); + SkyListFree (skylist); + } + for (int i = 0; i < Nstars; i++) { + free (stars[i].measure); + } + free (stars); + } + + // wait for last remote clients to finish + // NOTE: disable parallel mode for now: + // harvest_all (); + + // we should not actually reach this point + return FALSE; +} + +/* modifying to read in smaller blocks of the input file at a time + + each pass on readstars needs to ... + +*/ +