Index: trunk/Ohana/src/addstar/src/ReadImageHeader.c =================================================================== --- trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 36681) +++ trunk/Ohana/src/addstar/src/ReadImageHeader.c (revision 36833) @@ -44,4 +44,6 @@ image[0].NX = Nx; image[0].NY = Ny; + + image[0].refColor = NAN; if (!gfits_scan (header, "TZERO", "%d", 1, &image[0].tzero) && !ACCEPT_TIME) { Index: trunk/Ohana/src/checkastro/Makefile =================================================================== --- trunk/Ohana/src/checkastro/Makefile (revision 36833) +++ trunk/Ohana/src/checkastro/Makefile (revision 36833) @@ -0,0 +1,57 @@ +default: checkastro checkastro_client +help: + @echo "make options: checkastro (default)" + +include ../../Makefile.System +HOME = $(ROOT)/src/checkastro +BIN = $(HOME)/bin +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +include ../../Makefile.Common + +# programs may add their own internal requirements here +FULL_CFLAGS = $(BASE_CFLAGS) +FULL_CPPFLAGS = $(BASE_CPPFLAGS) +FULL_LDFLAGS = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS) + +checkastro: $(BIN)/checkastro.$(ARCH) +checkastro_client: $(BIN)/checkastro_client.$(ARCH) + +install: $(DESTBIN)/checkastro $(DESTBIN)/checkastro_client $(DESTBIN)/testparallax + +CHECKASTRO = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ +$(SRC)/checkastro.$(ARCH).o \ +$(SRC)/checkastro_images.$(ARCH).o \ +$(SRC)/MeasFilterTest.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o + +$(CHECKASTRO): $(INC)/checkastro.h +$(BIN)/checkastro.$(ARCH): $(CHECKASTRO) + +CHECKASTRO_CLIENT = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/checkastro_client.$(ARCH).o \ +$(SRC)/MeasFilterTest.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o + +$(CHECKASTRO_CLIENT): $(INC)/checkastro.h +$(BIN)/checkastro_client.$(ARCH): $(CHECKASTRO_CLIENT) + Index: trunk/Ohana/src/checkastro/include/checkastro.h =================================================================== --- trunk/Ohana/src/checkastro/include/checkastro.h (revision 36833) +++ trunk/Ohana/src/checkastro/include/checkastro.h (revision 36833) @@ -0,0 +1,154 @@ +# include +# include +# include +# include +# include +# include + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); + +// choose off_t or int depending on full-scale relphot analysis resources +// # define IDX_T off_t +# define IDX_T int + +typedef struct { + Average *average; // array of (minimal) average data + MeasureTiny *measure; // array of (minimal) measure data + SecFilt *secfilt; // array of secfilt data (matched to average by Nsecfilt) + off_t Naverage; + off_t Nmeasure; +} BrightCatalog; + +typedef struct { + Catalog *catalog; // array of catalogs generated + int NCATALOG; // number of catalogs allocated + int Ncatalog; // number of catalogs generated + int Nsecfilt; // number of catalogs generated + off_t *NAVERAGE; // allocated Averages per catalog + off_t *NMEASURE; // allocated Measures per catalog + int *index; // lookup table catID -> catalog[i] + int *catIDs; // lookup table catID <- catalog[i] + int maxID; // max catID value to date +} CatalogSplitter; + +/* global variables set in parameter file */ +# define DVO_MAX_PATH 1024 +char ImageCat[DVO_MAX_PATH]; +char GSCFILE[DVO_MAX_PATH]; +char CATDIR[DVO_MAX_PATH]; + +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char SKY_TABLE[DVO_MAX_PATH]; +int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +int HOST_ID; +char *HOSTDIR; +char *BCATALOG; + +unsigned int OBJ_ID_SRC; +unsigned int CAT_ID_SRC; +unsigned int OBJ_ID_DST; +unsigned int CAT_ID_DST; + +double SIGMA_LIM; +int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored +double MIN_ERROR; + +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +int VERBOSE; +int VERBOSE2; + +int MaxDensityUse; +double MaxDensityValue; + +char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST, *PHOTCODE_RESET_LIST; +int NphotcodesKeep, NphotcodesSkip, NphotcodesReset; +PhotCode **photcodesKeep, **photcodesSkip, **photcodesReset; + +int ImagSelect; +double ImagMin, ImagMax; + +int PhotFlagSelect, PhotFlagPoor, PhotFlagBad; + +float MinBadQF; + +int TimeSelect; +time_t TSTART, TSTOP; + +SkyRegion UserPatch; + +/*** relphot prototypes ***/ +void ConfigInit PROTO((int *argc, char **argv)); +void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); + +int args PROTO((int argc, char **argv)); +int args_client PROTO((int argc, char **argv)); +int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); + +void findImages PROTO((Catalog *catalog, int Ncatalog)); +int findMosaics PROTO((Catalog *catalog, int Ncatalog)); + +void set_db PROTO((FITS_DB *in)); +int Shutdown PROTO((char *format, ...)) OHANA_FORMAT(printf, 1, 2); +void TrapSignal PROTO((int sig)); +void SetProtect PROTO((int mode)); +int SetSignals PROTO((void)); + +void freeImageBins PROTO((int Ncatalog)); +void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); + +void initImageBins PROTO((Catalog *catalog, int Ncatalog)); +void initImages PROTO((Image *input, off_t *line_number, off_t N)); + +void initialize PROTO((int argc, char **argv)); +void initialize_client PROTO((int argc, char **argv)); + +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath)); +int load_images PROTO((FITS_DB *db, SkyList *skylist)); +Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage)); + +void check_permissions (char *basefile); +void lock_image_db (FITS_DB *db, char *filename); +void unlock_image_db (FITS_DB *db); +void create_image_db (FITS_DB *db); +void save_catalogs (Catalog *catalog, int Ncatalog); + +int reload_images (FITS_DB *db); + +int main PROTO((int argc, char **argv)); + +void matchImage PROTO((Catalog *catalog, off_t meas, int cat)); + +void set_ZP PROTO((double ZERO)); + +off_t getImageByID (off_t ID); + +int checkastro_objects (SkyList *skylist, int hostID, char *hostpath); +int checkastro_images (SkyList *skylist); + +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog); +int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog); + +BrightCatalog *BrightCatalogLoad(char *filename); +int BrightCatalogSave(char *filename, BrightCatalog *catalog); +BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog); +CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt); +int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog); +int BrightCatalogSplitFree (CatalogSplitter *catalogs); + +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog); +void bcatalog_show_skips (); + +int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim); +int MeasFilterTest(Measure *measure, int applySigmaLim); Index: trunk/Ohana/src/checkastro/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/checkastro/src/BrightCatalog.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/BrightCatalog.c (revision 36833) @@ -0,0 +1,906 @@ +# include "checkastro.h" + +// NOTE: this file is nearly identical to relphot/src/BrightCatalog.c, but here we use Average, not AverageTiny + +// checkastro_client -load reads from the local catalogs and generates a set of bright, subset catalogs +// we are going to save these as a single big FITS file, with only 3 extensions: + +// measure, average, secfilt. these only need to contain the fields relevant to the +// MeasureTiny, Average, and Secfilt tables, but with cat_id appended to each row so we can +// reassign correctly on read + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// XXX double check the header extname values for each table +// XXX make sure we free things as we can +// XXX make sure we close files as we can +// XXX handle free and close on error return as well +BrightCatalog *BrightCatalogLoad(char *filename) { + + int i, Ncol; + off_t Nrow; + char type[16]; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + BrightCatalog *catalog = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + goto escape; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + goto escape; + } + + ALLOCATE (catalog, BrightCatalog, 1); + catalog->Naverage = 0; + catalog->Nmeasure = 0; + + ftable.header = &theader; + + // *** Measure Data *** + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(dR, "RA_OFF", float); + GET_COLUMN(dD, "DEC_OFF", float); + GET_COLUMN(M, "MAG_SYS", float); + GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(airmass, "AIRMASS", float); + GET_COLUMN(Xccd, "X_CCD", float); + GET_COLUMN(Yccd, "Y_CCD", float); + GET_COLUMN(Xfix, "X_FIX", float); + GET_COLUMN(Yfix, "Y_FIX", float); + GET_COLUMN(dt, "EXPTIME", float); + GET_COLUMN(t, "TIME", int); + GET_COLUMN(averef, "AVE_REF", int); // XXX signed vs unsigned? + GET_COLUMN(imageID, "IMAGE_ID", int); + GET_COLUMN(dbFlags, "DB_FLAGS", int); + GET_COLUMN(photFlags, "PHOT_FLAGS", int); + GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(photcode, "PHOTCODE", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + MeasureTiny *measure = NULL; + ALLOCATE (measure, MeasureTiny, Nrow); + for (i = 0; i < Nrow; i++) { + measure[i].dR = dR[i]; + measure[i].dD = dD[i]; + measure[i].M = M[i]; + measure[i].Mcal = Mcal[i]; + measure[i].dM = dM[i]; + measure[i].airmass = airmass[i]; + measure[i].Xccd = Xccd[i]; + measure[i].Yccd = Yccd[i]; + measure[i].Xfix = Xfix[i]; + measure[i].Yfix = Yfix[i]; + measure[i].dt = dt[i]; + measure[i].t = t[i]; + measure[i].averef = averef[i]; + measure[i].imageID = imageID[i]; + measure[i].dbFlags = dbFlags[i]; + measure[i].photFlags = photFlags[i]; + measure[i].catID = catID[i]; + measure[i].photcode = photcode[i]; + } + // fprintf (stderr, "loaded data for %lld measures\n", (long long) Nrow); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (Xfix ); + free (Yfix ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + catalog->measure = measure; + catalog->Nmeasure = Nrow; + } + + /*** load Average values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + GET_COLUMN(R, "RA", double); + GET_COLUMN(D, "DEC", double); + GET_COLUMN(dR, "RA_ERR", float); + GET_COLUMN(dD, "DEC_ERR", float); + GET_COLUMN(uR, "U_RA", float); + GET_COLUMN(uD, "U_DEC", float); + GET_COLUMN(duR, "V_RA_ERR", float); + GET_COLUMN(duD, "V_DEC_ERR", float); + GET_COLUMN(P, "PAR", float); + GET_COLUMN(dP, "PAR_ERR", float); + GET_COLUMN(ChiSqAve, "CHISQ_POS", float); + GET_COLUMN(ChiSqPM, "CHISQ_PM", float); + GET_COLUMN(ChiSqPar, "CHISQ_PAP", float); + GET_COLUMN(Tmean, "MEAN_EPOCH", int); + GET_COLUMN(Trange, "TIME_RANGE", int); + GET_COLUMN(stargal, "STARGAL_SEP", float); + GET_COLUMN(Npos, "NUMBER_POS", short); + GET_COLUMN(Nmeasure, "NMEASURE", short); + GET_COLUMN(Nmissing, "NMISSING", short); + GET_COLUMN(Nextend, "NEXTEND", short); + GET_COLUMN(measureOffset, "OFF_MEASURE", int); + GET_COLUMN(missingOffset, "OFF_MISSING", int); + GET_COLUMN(refColor, "REF_COLOR", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(photFlagsUpper, "PHOTFLAGS_U", int); + GET_COLUMN(photFlagsLower, "PHOTFLAGS_L", int); + GET_COLUMN(objID, "OBJ_ID", int); + GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(extID, "EXT_ID", int64_t); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + Average *average = NULL; + ALLOCATE (average, Average, Nrow); + for (i = 0; i < Nrow; i++) { + average[i].R = R[i] ; + average[i].D = D[i] ; + average[i].dR = dR[i] ; + average[i].dD = dD[i] ; + average[i].uR = uR[i] ; + average[i].uD = uD[i] ; + average[i].duR = duR[i] ; + average[i].duD = duD[i] ; + average[i].P = P[i] ; + average[i].dP = dP[i] ; + average[i].ChiSqAve = ChiSqAve[i] ; + average[i].ChiSqPM = ChiSqPM[i] ; + average[i].ChiSqPar = ChiSqPar[i] ; + average[i].Tmean = Tmean[i] ; + average[i].Trange = Trange[i] ; + average[i].stargal = stargal[i] ; + average[i].Npos = Npos[i] ; + average[i].Nmeasure = Nmeasure[i] ; + average[i].Nmissing = Nmissing[i] ; + average[i].Nextend = Nextend[i] ; + average[i].measureOffset = measureOffset[i] ; + average[i].missingOffset = missingOffset[i] ; + average[i].refColor = refColor[i] ; + average[i].flags = flags[i] ; + average[i].photFlagsUpper = photFlagsUpper[i] ; + average[i].photFlagsLower = photFlagsLower[i] ; + average[i].objID = objID[i] ; + average[i].catID = catID[i] ; + average[i].extID = extID[i] ; + } + fprintf (stderr, "loaded data for %lld averages\n", (long long) Nrow); + + free (R); + free (D); + free (dR); + free (dD); + free (uR); + free (uD); + free (duR); + free (duD); + free (P); + free (dP); + free (ChiSqAve); + free (ChiSqPM); + free (ChiSqPar); + free (Tmean); + free (Trange); + free (stargal); + free (Npos); + free (Nmeasure); + free (Nmissing); + free (Nextend); + free (measureOffset); + free (missingOffset); + free (refColor); + free (flags); + free (photFlagsUpper); + free (photFlagsLower); + free (objID); + free (catID); + free (extID); + + catalog->average = average; + catalog->Naverage = Nrow; + } + + /*** load Secfilt values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(M, "MAG", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(Xm, "MAG_CHI", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(Ncode, "NCODE", short); + GET_COLUMN(Nused, "NUSED", short); + GET_COLUMN(M_20, "MAG_20", short); + GET_COLUMN(M_80, "MAG_80", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + SecFilt *secfilt = NULL; + ALLOCATE (secfilt, SecFilt, Nrow); + for (i = 0; i < Nrow; i++) { + secfilt[i].M = M[i]; + secfilt[i].dM = dM[i]; + secfilt[i].Xm = Xm[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].M_20 = M_20[i]; + secfilt[i].M_80 = M_80[i]; + } + fprintf (stderr, "loaded data for %lld averages\n", (long long) Nrow); + + free (M ); + free (dM ); + free (Xm ); + free (flags); + free (Ncode); + free (Nused); + free (M_20 ); + free (M_80 ); + catalog->secfilt = secfilt; + // assert Nsecfilt * Naverage = Nrow? + } + + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return catalog; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (catalog) free (catalog); + + fclose (f); + return NULL; +} + +// we are passed a BrightCatalog structure, write it to a FITS table (3 ext) +int BrightCatalogSave(char *filename, BrightCatalog *catalog) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + /*** MeasureTiny ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); + + gfits_define_bintable_column (&theader, "E", "RA_OFF", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_OFF", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "AIRMASS", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_CCD", "ccd x raw coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_CCD", "ccd y raw coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_FIX", "ccd x fixed coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_FIX", "ccd y fiex coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "EXPTIME", "-2.5 * log (exposure time)", "sec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "TIME", "time of exp", "sec", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "AVE_REF", "pointer to average table", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode", NULL, 1.0, 1.0*0x80); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); + float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); + float *M ; ALLOCATE (M , float, catalog->Nmeasure); + float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); + float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); + float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); + float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); + float *Xfix ; ALLOCATE (Xfix , float, catalog->Nmeasure); + float *Yfix ; ALLOCATE (Yfix , float, catalog->Nmeasure); + float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); + int *t ; ALLOCATE (t , int , catalog->Nmeasure); + int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); + int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); + int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); + int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); + int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); + short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); + + // assign the storage arrays + MeasureTiny *measure = catalog->measure; + for (i = 0; i < catalog->Nmeasure; i++) { + dR[i] = measure[i].dR ; + dD[i] = measure[i].dD ; + M[i] = measure[i].M ; + Mcal[i] = measure[i].Mcal ; + dM[i] = measure[i].dM ; + airmass[i] = measure[i].airmass ; + Xccd[i] = measure[i].Xccd ; + Yccd[i] = measure[i].Yccd ; + Xfix[i] = measure[i].Xfix ; + Yfix[i] = measure[i].Yfix ; + dt[i] = measure[i].dt ; + t[i] = measure[i].t ; + averef[i] = measure[i].averef ; + catID[i] = measure[i].catID ; + imageID[i] = measure[i].imageID ; + dbFlags[i] = measure[i].dbFlags ; + photFlags[i]= measure[i].photFlags; + photcode[i] = measure[i].photcode ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "X_CCD", Xccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "Y_CCD", Yccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "X_FIX", Xfix, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "Y_FIX", Yfix, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", dt, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "TIME", t, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AVE_REF", averef, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (Xfix ); + free (Yfix ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** Average ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "AVERAGE"); + + gfits_define_bintable_column (&theader, "D", "RA", "RA", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "DEC", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "RA_ERR", "RA error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_ERR", "DEC error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_RA", "RA*cos(D) proper-motion", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_DEC", "DEC proper-motion", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "V_RA_ERR", "RA*cos(D) p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "V_DEC_ERR", "DEC p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PAR", "parallax", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PAR_ERR", "parallax error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_POS", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_PM", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_PAP", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "MEAN_EPOCH", "mean epoch (PM-PAR ref)", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "TIME_RANGE", "mean epoch (PM-PAR ref)", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "E", "STARGAL_SEP", "star/galaxy separator", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NUMBER_POS", "number of detections used for astrometry", "", 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "I", "NMEASURE", "number of psf measurements", "", 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "I", "NMISSING", "number of missings", "", 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "I", "NEXTEND", "number of extended measurements", "", 1.0, 1.0*0x80); + gfits_define_bintable_column (&theader, "J", "OFF_MEASURE", "offset to first psf measurement", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "OFF_MISSING", "offset to first missing obs", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "OFF_EXTEND", "offset to first extended measurement", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "FLAGS", "average object flags (star; ghost; etc)", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOTFLAGS_U", "upper bit of 2 bit summary of per-measure photflags", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOTFLAGS_L", "lower bit of 2 bit summary of per-measure photflags", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "unique ID for object in table", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "unique ID for table in which object was first realized", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "K", "EXT_ID", "external ID for object (eg PSPS objID)", "", 1.0, 1.0*0x80000000); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + double *R ; ALLOCATE (R , double , catalog->Naverage); + double *D ; ALLOCATE (D , double , catalog->Naverage); + float *dR ; ALLOCATE (dR , float , catalog->Naverage); + float *dD ; ALLOCATE (dD , float , catalog->Naverage); + float *uR ; ALLOCATE (uR , float , catalog->Naverage); + float *uD ; ALLOCATE (uD , float , catalog->Naverage); + float *duR ; ALLOCATE (duR , float , catalog->Naverage); + float *duD ; ALLOCATE (duD , float , catalog->Naverage); + float *P ; ALLOCATE (P , float , catalog->Naverage); + float *dP ; ALLOCATE (dP , float , catalog->Naverage); + float *ChiSqAve ; ALLOCATE (ChiSqAve , float , catalog->Naverage); + float *ChiSqPM ; ALLOCATE (ChiSqPM , float , catalog->Naverage); + float *ChiSqPar ; ALLOCATE (ChiSqPar , float , catalog->Naverage); + int *Tmean ; ALLOCATE (Tmean , int , catalog->Naverage); + int *Trange ; ALLOCATE (Trange , int , catalog->Naverage); + float *stargal ; ALLOCATE (stargal , float , catalog->Naverage); + short *Npos ; ALLOCATE (Npos , short , catalog->Naverage); + short *Nmeasure ; ALLOCATE (Nmeasure , short , catalog->Naverage); + short *Nmissing ; ALLOCATE (Nmissing , short , catalog->Naverage); + short *Nextend ; ALLOCATE (Nextend , short , catalog->Naverage); + int *measureOffset ; ALLOCATE (measureOffset , int , catalog->Naverage); + int *missingOffset ; ALLOCATE (missingOffset , int , catalog->Naverage); + int *refColor ; ALLOCATE (refColor , float , catalog->Naverage); + int *flags ; ALLOCATE (flags , int , catalog->Naverage); + int *photFlagsUpper; ALLOCATE (photFlagsUpper, int , catalog->Naverage); + int *photFlagsLower; ALLOCATE (photFlagsLower, int , catalog->Naverage); + int *objID ; ALLOCATE (objID , int , catalog->Naverage); + int *catID ; ALLOCATE (catID , int , catalog->Naverage); + uint64_t *extID ; ALLOCATE (extID , uint64_t, catalog->Naverage); + + // assign the storage arrays + Average *average = catalog->average; + for (i = 0; i < catalog->Naverage; i++) { + R[i] = average[i].R ; + D[i] = average[i].D ; + dR[i] = average[i].dR ; + dD[i] = average[i].dD ; + uR[i] = average[i].uR ; + uD[i] = average[i].uD ; + duR[i] = average[i].duR ; + duD[i] = average[i].duD ; + P[i] = average[i].P ; + dP[i] = average[i].dP ; + ChiSqAve[i] = average[i].ChiSqAve ; + ChiSqPM[i] = average[i].ChiSqPM ; + ChiSqPar[i] = average[i].ChiSqPar ; + Tmean[i] = average[i].Tmean ; + Trange[i] = average[i].Trange ; + stargal[i] = average[i].stargal ; + Npos[i] = average[i].Npos ; + Nmeasure[i] = average[i].Nmeasure ; + Nmissing[i] = average[i].Nmissing ; + Nextend[i] = average[i].Nextend ; + measureOffset[i] = average[i].measureOffset ; + missingOffset[i] = average[i].missingOffset ; + refColor[i] = average[i].refColor ; + flags[i] = average[i].flags ; + photFlagsUpper[i] = average[i].photFlagsUpper ; + photFlagsLower[i] = average[i].photFlagsLower ; + objID[i] = average[i].objID ; + catID[i] = average[i].catID ; + extID[i] = average[i].extID ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "RA_ERR", dR, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "DEC_ERR", dD, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_RA", uR, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "U_DEC", uD, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "V_RA_ERR", duR, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "V_DEC_ERR", duD, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "PAR", P, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "PAR_ERR", dP, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "CHISQ_POS", ChiSqAve, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "CHISQ_PM", ChiSqPM, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "CHISQ_PAP", ChiSqPar, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "MEAN_EPOCH", Tmean, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "TIME_RANGE", Trange, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "STARGAL_SEP", stargal, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NUMBER_POS", Npos, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NMEASURE", Nmeasure, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NMISSING", Nmissing, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "NEXTEND", Nextend, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "REF_COLOR", refColor, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_U", photFlagsUpper, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_L", photFlagsLower, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "EXT_ID", extID, catalog->Naverage); + + free (R); + free (D); + free (dR); + free (dD); + free (uR); + free (uD); + free (duR); + free (duD); + free (P); + free (dP); + free (ChiSqAve); + free (ChiSqPM); + free (ChiSqPar); + free (Tmean); + free (Trange); + free (stargal); + free (Npos); + free (Nmeasure); + free (Nmissing); + free (Nextend); + free (measureOffset); + free (missingOffset); + free (refColor); + free (flags); + free (photFlagsUpper); + free (photFlagsLower); + free (objID); + free (catID); + free (extID); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** SecFilt ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "SECFILT"); + + gfits_define_bintable_column (&theader, "E", "MAG", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHI", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NCODE", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NUSED", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_20", "ccd x coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_80", "ccd y coord", "pix", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // Nsecfilt is number of average filters, Nsec is total number of values + int Nsecfilt = GetPhotcodeNsecfilt(); + int Nsec = Nsecfilt * catalog->Naverage; + + // create intermediate storage arrays + float *M ; ALLOCATE (M , float, Nsec); + float *dM ; ALLOCATE (dM , float, Nsec); + float *Xm ; ALLOCATE (Xm , float, Nsec); + int *flags ; ALLOCATE (flags , int, Nsec); + short *Ncode ; ALLOCATE (Ncode , short, Nsec); + short *Nused ; ALLOCATE (Nused , short, Nsec); + short *M_20 ; ALLOCATE (M_20 , short, Nsec); + short *M_80 ; ALLOCATE (M_80 , short, Nsec); + + // assign the storage arrays + SecFilt *secfilt = catalog->secfilt; + for (i = 0; i < Nsec; i++) { + M [i] = secfilt[i]. M ; + dM [i] = secfilt[i]. dM ; + Xm [i] = secfilt[i]. Xm ; + flags [i] = secfilt[i]. flags ; + Ncode [i] = secfilt[i]. Ncode ; + Nused [i] = secfilt[i]. Nused ; + M_20 [i] = secfilt[i]. M_20 ; + M_80 [i] = secfilt[i]. M_80 ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Xm , Nsec); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused, Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_20", M_20 , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_80", M_80 , Nsec); + + free (M ); + free (dM ); + free (Xm ); + free (flags ); + free (Ncode ); + free (Nused ); + free (M_20 ); + free (M_80 ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +// merge a list of catalogs into a single BrightCatalog array set +BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog) { + + off_t i, j, k; + + BrightCatalog *bcatalog = NULL; + ALLOCATE (bcatalog, BrightCatalog, 1); + + int Nmeas = 0; + int Naves = 0; + for (i = 0; i < Ncatalog; i++) { + Nmeas += catalog[i].Nmeasure; + Naves += catalog[i].Naverage; + } + + // XXX this prevents different catalogs from having different Nsecfilt values + int Nsecfilt = GetPhotcodeNsecfilt(); + + ALLOCATE (bcatalog[0].measure, MeasureTiny, Nmeas); + ALLOCATE (bcatalog[0].average, Average, Naves); + ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt); + + int Nm = 0; + int Na = 0; + for (i = 0; i < Ncatalog; i++) { + if (!catalog[i].Naverage) continue; + for (j = 0; j < catalog[i].Naverage; j++) { + bcatalog[0].average[Na] = catalog[i].average[j]; + for (k = 0; k < Nsecfilt; k++) { + bcatalog[0].secfilt[Nsecfilt*Na + k] = catalog[i].secfilt[Nsecfilt*j + k]; + } + Na++; + assert (Na <= Naves); + } + for (j = 0; j < catalog[i].Nmeasure; j++) { + bcatalog[0].measure[Nm] = catalog[i].measureT[j]; + Nm++; + assert (Nm <= Nmeas); + } + } + bcatalog->Naverage = Na; + bcatalog->Nmeasure = Nm; + return bcatalog; +} + +// distribute a bright catalog across separate catalogs +CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt) { + + int i; + + CatalogSplitter *catalogs = NULL; + + ALLOCATE (catalogs, CatalogSplitter, 1); + + // as we see new BrightCatalogs, we will update maxID and extend this array + catalogs->maxID = 0; + ALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = 0; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + catalogs->Nsecfilt = Nsecfilt; + + catalogs->Ncatalog = 0; + catalogs->NCATALOG = 16; + catalogs->catalog = NULL; + ALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + + ALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + ALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + ALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + + for (i = 0; i < catalogs->NCATALOG; i++) { + dvo_catalog_init (&catalogs->catalog[i], TRUE); + catalogs->catIDs[i] = 0; + catalogs->NAVERAGE[i] = 100; + catalogs->NMEASURE[i] = 100; + catalogs->catalog[i].Naverage = 0; + catalogs->catalog[i].Nmeasure = 0; + ALLOCATE (catalogs->catalog[i].average, Average, catalogs->NAVERAGE[i]); + ALLOCATE (catalogs->catalog[i].measureT, MeasureTiny, catalogs->NMEASURE[i]); + ALLOCATE (catalogs->catalog[i].secfilt, SecFilt, catalogs->NAVERAGE[i]*Nsecfilt); + } + return catalogs; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplitFree (CatalogSplitter *catalogs) { + + free (catalogs->catIDs); + free (catalogs->NAVERAGE); + free (catalogs->NMEASURE); + free (catalogs->index); + free (catalogs); + return TRUE; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog) { + + int i; + + int Nsecfilt = catalogs->Nsecfilt; + + // find the max value of catID in this BrightCatalog + int catIDmax = 0; + for (i = 0; i < bcatalog->Naverage; i++) { + catIDmax = MAX(catIDmax, bcatalog->average[i].catID); + } + // XXX validate the measure ID range here + + int maxIDold = catalogs->maxID; + catalogs->maxID = MAX (maxIDold, catIDmax); + + // extend the index array and init + REALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = maxIDold + 1; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + // identify the new catID values + for (i = 0; i < bcatalog->Naverage; i++) { + int catID = bcatalog->average[i].catID; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + int idx = catalogs->index[catID]; + if (idx != -1) continue; // already have seen this one + + // a new catID value + int Ncat = catalogs->Ncatalog; // the next available slot + catalogs->catIDs[Ncat] = catID; + assert (Ncat >= 0); + assert (Ncat < catalogs->NCATALOG); + + catalogs->catalog[Ncat].Nsecfilt = Nsecfilt; + catalogs->catalog[Ncat].catID = catID; + + catalogs->index[catID] = Ncat; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + + catalogs->Ncatalog ++; + + if (catalogs->Ncatalog >= catalogs->NCATALOG) { + catalogs->NCATALOG += 16; + // fprintf (stderr, "realloc catalogs->catalog: old: %llx ", (long long) catalogs->catalog); + REALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->catalog, (long long) (catalogs->catalog + sizeof(Catalog)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NAVERAGE: old: %llx ", (long long) catalogs->NAVERAGE); + REALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NAVERAGE, (long long) (catalogs->NAVERAGE + sizeof(off_t)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NMEASURE: old: %llx ", (long long) catalogs->NMEASURE); + REALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NMEASURE, (long long) (catalogs->NMEASURE + sizeof(off_t)*catalogs->NCATALOG)); + + REALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + + int j; + for (j = catalogs->NCATALOG - 16; j < catalogs->NCATALOG; j++) { + dvo_catalog_init (&catalogs->catalog[j], TRUE); + catalogs->catIDs[j] = 0; + catalogs->NAVERAGE[j] = 100; + catalogs->NMEASURE[j] = 100; + catalogs->catalog[j].Naverage = 0; + catalogs->catalog[j].Nmeasure = 0; + ALLOCATE (catalogs->catalog[j].average, Average, catalogs->NAVERAGE[j]); + ALLOCATE (catalogs->catalog[j].measureT, MeasureTiny, catalogs->NMEASURE[j]); + ALLOCATE (catalogs->catalog[j].secfilt, SecFilt, catalogs->NAVERAGE[j]*Nsecfilt); + } + } + } + + // each bcatalog (from each host) has a different set of catID ranges + // they also (probably) have contiguous ranges. But, it is a bit tricky to be sure I + // can use that info, so perhaps ignore it + + // assign the averages to the corresponding catalog + for (i = 0; i < bcatalog->Naverage; i++) { + int ID = bcatalog->average[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + + int Na = catalogs->catalog[Nc].Naverage; + catalogs->catalog[Nc].average[Na] = bcatalog->average[i]; + + // secfilt entries are grouped in blocks for each average + int k; + for (k = 0; k < Nsecfilt; k++) { + catalogs->catalog[Nc].secfilt[Na*Nsecfilt + k] = bcatalog->secfilt[i*Nsecfilt + k]; + } + + catalogs->catalog[Nc].Naverage ++; + if (catalogs->catalog[Nc].Naverage >= catalogs->NAVERAGE[Nc]) { + catalogs->NAVERAGE[Nc] += 100; + REALLOCATE (catalogs->catalog[Nc].average, Average, catalogs->NAVERAGE[Nc]); + REALLOCATE (catalogs->catalog[Nc].secfilt, SecFilt, catalogs->NAVERAGE[Nc]*Nsecfilt); + } + } + + // assign the measures to the corresponding catalog + // XXX what about averef and related links? Do I need them? + for (i = 0; i < bcatalog->Nmeasure; i++) { + int ID = bcatalog->measure[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + int Na = catalogs->catalog[Nc].Nmeasure; + catalogs->catalog[Nc].measureT[Na] = bcatalog->measure[i]; + catalogs->catalog[Nc].Nmeasure ++; + if (catalogs->catalog[Nc].Nmeasure >= catalogs->NMEASURE[Nc]) { + catalogs->NMEASURE[Nc] += 100; + REALLOCATE (catalogs->catalog[Nc].measureT, MeasureTiny, catalogs->NMEASURE[Nc]); + } + } + + return TRUE; +} Index: trunk/Ohana/src/checkastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/checkastro/src/ConfigInit.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/ConfigInit.c (revision 36833) @@ -0,0 +1,72 @@ +# include "checkastro.h" + +void ConfigInit (int *argc, char **argv) { + + char *config, *file; + char CatdirPhotcodeFile[256]; + struct stat filestat; + int status; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + // set defaults for all of these if they are not used by parallel / remote clients + if (!ScanConfig (config, "RELASTRO_SIGMA_LIM", "%lf", 0, &SIGMA_LIM)) SIGMA_LIM = 0.01; + if (!ScanConfig (config, "RELASTRO_SRC_MEAS_TOOFEW", "%d", 0, &SRC_MEAS_TOOFEW)) SRC_MEAS_TOOFEW = 3; + + // force CATDIR to be absolute (so parallel mode will work) + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + char *tmpcatdir = abspath (CATDIR, DVO_MAX_PATH); + strcpy (CATDIR, tmpcatdir); + free (tmpcatdir); + + // GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig(config, "CATMODE", "%s", 0, CATMODE); + ScanConfig(config, "CATFORMAT", "%s", 0, CATFORMAT); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) SKY_DEPTH = 2; + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) SKY_TABLE[0] = 0; + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + // check for existence of CATDIR + status = stat (CATDIR, &filestat); + if (status == -1) { + fprintf (stderr, "directory %s does not exist, giving up\n", CATDIR); + exit (1); + } + + /* update master photcode table if not defined */ + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, NULL, FALSE)) { + fprintf (stderr, "error loading photcode table %s\n", CatdirPhotcodeFile); + exit (1); + } + SetZeroPoint (25.0); + + free (config); + free (file); + +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: trunk/Ohana/src/checkastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/checkastro/src/ImageOps.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/ImageOps.c (revision 36833) @@ -0,0 +1,253 @@ +# include "checkastro.h" + +// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage) +static Image *image; // list of available images +static off_t Nimage; // number of available images + +// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row +// each image belongs +static off_t *LineNumber; // match of subset to full image table + +static off_t *N_onImage; // number of measurements on image +static off_t *N_ONIMAGE; // allocated number of measurements on image +static off_t *Nref_onImage; // number of reference measurements on image + +static IDX_T **MeasureToImage; // link from catalog,measure to image +static IDX_T **ImageToCatalog; // catalog which supplied measurement on image +static IDX_T **ImageToMeasure; // measure reference for measurement on image + +// to search by image ID, we sort (imageIDs, imageIdx) by imageIDs for an index +static off_t *imageIDs; // list of all image IDs +static off_t *imageIdx; // list of index for image IDs + +// generate the image ID index (for imageID -> image[i] lookups) +void initImages (Image *input, off_t *line_number, off_t N) { + + off_t i; + + image = input; + LineNumber = line_number; + Nimage = N; + + ALLOCATE (imageIDs, off_t, Nimage); + ALLOCATE (imageIdx, off_t, Nimage); + + for (i = 0; i < Nimage; i++) { + imageIdx[i] = i; + imageIDs[i] = image[i].imageID; + } + llsortpair (imageIDs, imageIdx, Nimage); +} + +off_t getImageByID (off_t ID) { + + // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs + // use bisection to find the specified image ID + + off_t Nlo, Nhi, N; + + Nlo = 0; Nhi = Nimage; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (imageIDs[N] < ID) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N + 1, Nimage); + } + } + + for (N = Nlo; N < Nhi; N++) { + if (imageIDs[N] == ID) + return (imageIdx[N]); + } + + return (-1); +} + +// these are really image & catalog indexes +void initImageBins (Catalog *catalog, int Ncatalog) { + + off_t i, j; + + ALLOCATE (MeasureToImage, IDX_T *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1; + } + + ALLOCATE (N_onImage, off_t, Nimage); + ALLOCATE (N_ONIMAGE, off_t, Nimage); + ALLOCATE (Nref_onImage, off_t, Nimage); + ALLOCATE (ImageToCatalog, IDX_T *, Nimage); + ALLOCATE (ImageToMeasure, IDX_T *, Nimage); + + for (i = 0; i < Nimage; i++) { + N_onImage[i] = 0; + N_ONIMAGE[i] = 30; + Nref_onImage[i] = 0; + ImageToCatalog[i] = NULL; // we allocate these iff they are needed in matchImage + ImageToMeasure[i] = NULL; // we allocate these iff they are needed in matchImage + } +} + +void freeImageBins (int Ncatalog) { + + off_t i; + + for (i = 0; i < Ncatalog; i++) { + free (MeasureToImage[i]); + } + free (MeasureToImage); + for (i = 0; i < Nimage; i++) { + if (ImageToCatalog[i]) { free (ImageToCatalog[i]); } + if (ImageToMeasure[i]) { free (ImageToMeasure[i]); } + } + free (ImageToCatalog); + free (ImageToMeasure); + free (N_onImage); + free (N_ONIMAGE); +} + +/* match measurements to images */ +void findImages (Catalog *catalog, int Ncatalog) { + + off_t i, j; + char *name; + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Nmeasure; j++) { + if (NphotcodesSkip > 0) { + int k; + int found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == catalog[i].measureT[j].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[i].measureT[j].photcode)) found = TRUE; + } + if (found) continue; + } + matchImage (catalog, j, i); + } + } + + BuildChipMatch (image, Nimage); + + char output[128]; + snprintf (output, 128, "checkastrom.%+3.0f.%+3.0f.%03.0f.%03.0f.dat", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + FILE *f = fopen (output, "w"); + if (!f) { + fprintf (stderr, "problem opening %s for output\n", output); + } + + // watch for chips with few stars + int Nfew = 0; + int Nbad = 0; + for (i = 0; i < Nimage; i++) { + // ignore the PHU entries + if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + + name = GetPhotcodeNamebyCode (image[i].photcode); + + off_t myMosaic = FindMosaicForImage (image, Nimage, i); + + /* only check exposure center */ + double Rexp = NAN, Dexp = NAN, Rccd = NAN, Dccd = NAN; + + if (myMosaic) { + myMosaic --; + XY_to_RD (&Rccd, &Dccd, 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); + XY_to_RD (&Rexp, &Dexp, 0.0, 0.0, &image[myMosaic].coords); + } + + char *date = ohana_sec_to_date(image[i].tzero); + FILE *outfile = f ? f : stderr; + fprintf (outfile, "%d %s : %10.6f %10.6f : %10.6f %10.6f "OFF_T_FMT" "OFF_T_FMT" %d %s %s\n", + image[i].imageID, image[i].name, Rccd, Dccd, Rexp, Dexp, N_onImage[i], Nref_onImage[i], image[i].nstar, + date, name); + free (date); + + if (N_onImage[i] < 20) { + Nfew ++; + } + if (N_onImage[i] < 15) { + Nbad ++; + } + } + if (f) { fclose (f); } + fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad); +} + +// this is the imageID-based match +void matchImage (Catalog *catalog, off_t meas, int cat) { + + off_t idx, ID; + MeasureTiny *measure; + + measure = &catalog[cat].measureT[meas]; + unsigned int averef = measure->averef; + Average *average = &catalog[cat].average[averef]; + assert (meas >= average->measureOffset); + assert (meas < average->measureOffset + average->Nmeasure); + + ID = measure[0].imageID; + + if (catalog[cat].measure) { + Measure *measureBig = &catalog[cat].measure[meas]; + int TESTPT = FALSE; + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measureBig->imageID == CAT_ID_SRC) && (measureBig->detID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measureBig->imageID == CAT_ID_DST) && (measureBig->detID == OBJ_ID_DST); + if (TESTPT) { + fprintf (stderr, "got test det\n"); + } + } + + idx = getImageByID (ID); + if (idx == -1) { + if (VERBOSE2) fprintf (stderr, "can't match detection to image?\n"); + return; + } + + // index for (catalog, measure) -> image + MeasureToImage[cat][meas] = idx; + + // if we need to allocate an image index table, do so here + if (!ImageToCatalog[idx]) { + ALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]); + } + if (!ImageToMeasure[idx]) { + ALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]); + } + + // index for image, Nentry -> catalog + ImageToCatalog[idx][N_onImage[idx]] = cat; + + // index for image, Nentry -> measure + ImageToMeasure[idx][N_onImage[idx]] = meas; + + // in bcatalog, we set this flag if the object includes a reference photcode + if (average->flags & 0x10) { + Nref_onImage[idx] ++; + } + + // count the number of measurements on this image: + N_onImage[idx] ++; + + if (N_onImage[idx] == N_ONIMAGE[idx]) { + N_ONIMAGE[idx] += 30; + REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]); + REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]); + } + + return; +} + +void dump_measures (Average *average, Measure *measure) { + + off_t j, off; + + for (j = 0; j < average[0].Nmeasure; j++) { + off = average[0].measureOffset + j; + fprintf (stderr, "dR, dD, mag, dMag: %f, %f, %f, %f\n", measure[off].dR, measure[off].dD, measure[off].M, measure[off].dM); + } + return; +} + Index: trunk/Ohana/src/checkastro/src/MeasFilterTest.c =================================================================== --- trunk/Ohana/src/checkastro/src/MeasFilterTest.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/MeasFilterTest.c (revision 36833) @@ -0,0 +1,133 @@ +# include "checkastro.h" + +/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ +// we only optionally apply the sigma limit: for object averages, this should not be used (should it?) +int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) { + int found, k; + long mask; + PhotCode *code; + float mag; + + if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE; + if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all checkastro tasks? + if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all checkastro tasks? + + /* select measurements by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (!found) return FALSE; + } + + if (NphotcodesSkip > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (found) return FALSE; + } + + // if (MinBadQF > 0.0) { + // if (measure[0].psfQF < MinBadQF) return FALSE; + // } + + /* select measurements by time */ + if (TimeSelect) { + if (measure[0].t < TSTART) return FALSE; + if (measure[0].t > TSTOP) return FALSE; + } + + /* select measurements by quality */ + if (PhotFlagSelect) { + if (PhotFlagBad) { + mask = PhotFlagBad; + } else { + code = GetPhotcodebyCode (measure[0].photcode); + mask = code[0].astromBadMask; + } + if (mask & measure[0].photFlags) return FALSE; + } + + /* select measurements by measurement error */ + if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) { + return FALSE; + } + + /* select measurements by mag limit */ + if (ImagSelect) { + mag = PhotInstTiny (measure); + if (mag < ImagMin || mag > ImagMax) return FALSE; + } + + return TRUE; +} + +/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ +// we only optionally apply the sigma limit: for object averages, this should not be used (should it?) +int MeasFilterTest(Measure *measure, int applySigmaLim) { + int found, k; + long mask; + PhotCode *code; + float mag; + + if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE; + if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all checkastro tasks? + if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all checkastro tasks? + + /* select measurements by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (!found) return FALSE; + } + + if (NphotcodesSkip > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + if (found) return FALSE; + } + + if (MinBadQF > 0.0) { + if (measure[0].psfQF < MinBadQF) return FALSE; + } + + /* select measurements by time */ + if (TimeSelect) { + if (measure[0].t < TSTART) return FALSE; + if (measure[0].t > TSTOP) return FALSE; + } + + /* select measurements by quality */ + if (PhotFlagSelect) { + if (PhotFlagBad) { + mask = PhotFlagBad; + } else { + code = GetPhotcodebyCode (measure[0].photcode); + mask = code[0].astromBadMask; + } + if (mask & measure[0].photFlags) return FALSE; + } + + /* select measurements by measurement error */ + if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) { + return FALSE; + } + + /* select measurements by mag limit */ + if (ImagSelect) { + mag = PhotInst (measure); + if (mag < ImagMin || mag > ImagMax) return FALSE; + } + + return TRUE; +} Index: trunk/Ohana/src/checkastro/src/SetSignals.c =================================================================== --- trunk/Ohana/src/checkastro/src/SetSignals.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/SetSignals.c (revision 36833) @@ -0,0 +1,122 @@ +# include "checkastro.h" + +static int Protect = FALSE; +static int Trapped = FALSE; + +void TrapSignal (int sig) { + fprintf (stderr, "trapped signal %d\n", sig); + if (sig == 11) { + fprintf (stderr, "seg fault\n"); + exit (1); + } + if (Protect) { + Trapped = TRUE; + fprintf (stderr, "blocking until protected sections are clear\n"); + return; + } + Shutdown ("halted by signal (trapped)"); +} + +void SetProtect (int mode) { + Protect = mode; + if (Trapped && !Protect) Shutdown ("halted by signal (protect)"); +} + +int SetSignals () { + + int i; + + /* disable almost all signal interrupts */ + for (i = 0; i < 36; i++) { + switch (i) { + /* can't redirect this signals */ + case SIGKILL: /* kill -9: cannot be caught or ignored */ + case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ + case SIGCHLD: /* child halted: ignore */ + case SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ +# ifdef SIGPWR + case SIGPWR: /* power failure - why ignore this? (Sys V) */ +# endif +# ifdef SIGWINCH + case SIGWINCH: /* window resized (4.3BSD) */ +# endif + break; + + default: + signal (i, TrapSignal); + } + } + return (TRUE); +} +/* + + Signal Value Action Comment + ------------------------------------------------------------------------- + SIGHUP 1 A Hangup detected on controlling terminal + or death of controlling process + SIGINT 2 A Interrupt from keyboard + SIGQUIT 3 A Quit from keyboard + SIGILL 4 A Illegal Instruction + SIGABRT 6 C Abort signal from abort(3) + SIGFPE 8 C Floating point exception + SIGKILL 9 AEF Kill signal + SIGSEGV 11 C Invalid memory reference + SIGPIPE 13 A Broken pipe: write to pipe with no readers + SIGALRM 14 A Timer signal from alarm(2) + SIGTERM 15 A Termination signal + SIGUSR1 30,10,16 A User-defined signal 1 + SIGUSR2 31,12,17 A User-defined signal 2 + SIGCHLD 20,17,18 B Child stopped or terminated + SIGCONT 19,18,25 Continue if stopped + SIGSTOP 17,19,23 DEF Stop process + SIGTSTP 18,20,24 D Stop typed at tty + SIGTTIN 21,21,26 D tty input for background process + SIGTTOU 22,22,27 D tty output for background process + + Next various other signals. + + Signal Value Action Comment + --------------------------------------------------------------------- + SIGTRAP 5 CG Trace/breakpoint trap + SIGIOT 6 CG IOT trap. A synonym for SIGABRT + SIGEMT 7,-,7 G + SIGBUS 10,7,10 AG Bus error + SIGSYS 12,-,12 G Bad argument to routine (SVID) + SIGSTKFLT -,16,- AG Stack fault on coprocessor + SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) + SIGIO 23,29,22 AG I/O now possible (4.2 BSD) + SIGPOLL AG A synonym for SIGIO (System V) + SIGCLD -,-,18 G A synonym for SIGCHLD + SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) + SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) + SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) + SIGPROF 27,27,29 AG Profile alarm clock + SIGPWR 29,30,19 AG Power failure (System V) + SIGINFO 29,-,- G A synonym for SIGPWR + SIGLOST -,-,- AG File lock lost + SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) + SIGUNUSED -,31,- AG Unused signal + (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and + sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a + sparc.) + + The letters in the "Action" column have the following meanings: + + A Default action is to terminate the process. + + B Default action is to ignore the signal. + + C Default action is to dump core. + + D Default action is to stop the process. + + E Signal cannot be caught. + + F Signal cannot be ignored. + + G Not a POSIX.1 conformant signal. + +*/ Index: trunk/Ohana/src/checkastro/src/Shutdown.c =================================================================== --- trunk/Ohana/src/checkastro/src/Shutdown.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/Shutdown.c (revision 36833) @@ -0,0 +1,30 @@ +# include "checkastro.h" + +static FITS_DB *db; + +void set_db (FITS_DB *in) { + db = in; +} + +/* clean up open / locked ImageCat before shutting down */ +int Shutdown (char *format, ...) { + va_list argp; + char *formatplus; + + ALLOCATE (formatplus, char, strlen(format)); + strcpy (formatplus, format); + strcat (formatplus, "\n"); + + va_start (argp, format); + vfprintf (stderr, formatplus, argp); + free (formatplus); + va_end (argp); + + SetProtect (TRUE); + gfits_db_close (db); + fprintf (stderr, "ERROR: addstar halted\n"); + exit (1); +} + + +/* XXX this is probably not needed anymore : just protect the write statements */ Index: trunk/Ohana/src/checkastro/src/args.c =================================================================== --- trunk/Ohana/src/checkastro/src/args.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/args.c (revision 36833) @@ -0,0 +1,433 @@ +# include "checkastro.h" +void usage (void); +void usage_client (void); +void usage_merge_source (void); +void usage_merge_source_id (char *name); + +int args (int argc, char **argv) { + + int N; + double trange; + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + UserPatch.Rmax = UserPatch.Rmin + 0.001; + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + UserPatch.Dmax = UserPatch.Dmin + 0.001; + remove_argument (N, &argc, argv); + } + + // XXX for the moment, make this selection manual. it needs to be automatic + // based on the state of the SkyTable + 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, + // checkastro 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 checkastro_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); + } + + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + if (!ohana_str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + PHOTCODE_RESET_LIST = NULL; + if ((N = get_argument (argc, argv, "-reset-to-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_RESET_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "+photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_SKIP_LIST = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9"); + if ((N = get_argument (argc, argv, "-photcode"))) { + remove_argument (N, &argc, argv); + char *tmp1 = strcreate(argv[N]); + + int Ntotal = strlen(tmp1) + strlen(PHOTCODE_SKIP_LIST) + 5; + + char *tmp2 = NULL; + ALLOCATE (tmp2, char, Ntotal); + snprintf (tmp2, Ntotal, "%s,%s", PHOTCODE_SKIP_LIST, tmp1); + + free (tmp1); + free (PHOTCODE_SKIP_LIST); + + PHOTCODE_SKIP_LIST = tmp2; + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE = VERBOSE2 = TRUE; + remove_argument (N, &argc, argv); + } + + /* XXX update this */ + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + // for now, make the default to ignore the photflags + // XXX make it true by default instead? + PhotFlagSelect = FALSE; + if ((N = get_argument (argc, argv, "+photflags"))) { + remove_argument (N, &argc, argv); + PhotFlagSelect = TRUE; + } + + PhotFlagBad = 0; + if ((N = get_argument (argc, argv, "-photflagbad"))) { + remove_argument (N, &argc, argv); + PhotFlagBad = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + PhotFlagPoor = 0; + if ((N = get_argument (argc, argv, "-photflagpoor"))) { + remove_argument (N, &argc, argv); + PhotFlagPoor = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + MinBadQF = 0.0; + if ((N = get_argument (argc, argv, "-min-bad-psfqf"))) { + remove_argument (N, &argc, argv); + MinBadQF = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 1) usage (); + return TRUE; +} + +int args_client (int argc, char **argv) { + + int N; + double trange; + + /* possible operations */ + // by definition, the client is not parallel + PARALLEL = FALSE; + PARALLEL_MANUAL = FALSE; + PARALLEL_SERIAL = FALSE; + + BCATALOG = NULL; + + HOST_ID = 0; + if ((N = get_argument (argc, argv, "-hostID"))) { + remove_argument (N, &argc, argv); + HOST_ID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOST_ID) usage_client(); + + HOSTDIR = NULL; + if ((N = get_argument (argc, argv, "-hostdir"))) { + remove_argument (N, &argc, argv); + HOSTDIR = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + if (!HOSTDIR) usage_client(); + + if ((N = get_argument (argc, argv, "-load-objects"))) { + remove_argument (N, &argc, argv); + BCATALOG = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-catalog"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + UserPatch.Rmax = UserPatch.Rmin + 0.001; + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + UserPatch.Dmax = UserPatch.Dmin + 0.001; + remove_argument (N, &argc, argv); + } + + MaxDensityUse = FALSE; + if ((N = get_argument (argc, argv, "-max-density"))) { + remove_argument (N, &argc, argv); + MaxDensityValue = atof(argv[N]); + remove_argument (N, &argc, argv); + MaxDensityUse = TRUE; + } + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + if (!ohana_str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + PHOTCODE_RESET_LIST = NULL; + if ((N = get_argument (argc, argv, "-reset-to-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_RESET_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "+photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_SKIP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_SKIP_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = VERBOSE2 = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-vv"))) { + VERBOSE = VERBOSE2 = TRUE; + remove_argument (N, &argc, argv); + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + // for now, make the default to ignore the photflags + // XXX make it true by default instead? + PhotFlagSelect = FALSE; + if ((N = get_argument (argc, argv, "+photflags"))) { + remove_argument (N, &argc, argv); + PhotFlagSelect = TRUE; + } + + PhotFlagBad = 0; + if ((N = get_argument (argc, argv, "-photflagbad"))) { + remove_argument (N, &argc, argv); + PhotFlagBad = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + PhotFlagPoor = 0; + if ((N = get_argument (argc, argv, "-photflagpoor"))) { + remove_argument (N, &argc, argv); + PhotFlagPoor = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + MinBadQF = 0.0; + if ((N = get_argument (argc, argv, "-min-bad-psfqf"))) { + remove_argument (N, &argc, argv); + MinBadQF = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 1) usage_client (); + return TRUE; +} + +void usage () { + fprintf (stderr, "ERROR: USAGE: checkastro [options]\n"); + fprintf (stderr, " options: \n"); + fprintf (stderr, " -region RA RA DEC DEC\n"); + fprintf (stderr, " -catalog (ra) (dec)\n\n"); + fprintf (stderr, " -time (start)(stop)\n"); + fprintf (stderr, " +photcode (code)[,code,code...]\n"); + fprintf (stderr, " -photcode (code)[,code,code...]\n"); + fprintf (stderr, " -reset-to-photcode (code)[,code,code...]\n"); + fprintf (stderr, " -max-density (rho)\n"); + fprintf (stderr, " -minerror\n"); + fprintf (stderr, " -minerror\n"); + fprintf (stderr, " -instmag min max\n\n"); + fprintf (stderr, " +photflags\n"); + fprintf (stderr, " -photflagbad\n"); + fprintf (stderr, " -photflagpoor\n"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " -vv\n"); + fprintf (stderr, " -parallel\n"); + fprintf (stderr, " -parallel-manual\n"); + fprintf (stderr, " -parallel-serial\n"); + fprintf (stderr, " \n"); + exit (2); +} + +void usage_client () { + fprintf (stderr, "ERROR: USAGE: checkastro_client -load\n"); + fprintf (stderr, " OR: checkastro_client -update-offsets\n"); + fprintf (stderr, " OR: checkastro_client -update-objects\n"); + fprintf (stderr, " OR: checkastro_client -high-speed\n\n"); + + fprintf (stderr, " specify one of the following modes: \n"); + fprintf (stderr, " -load : load the bright source detections for analysis\n"); + fprintf (stderr, " -update-offsets : apply the updated image parameters\n"); + fprintf (stderr, " -update-objects : calculate average astrometric properties of objects\n"); + fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n"); + fprintf (stderr, " support for the special high-speed mode\n"); + + fprintf (stderr, " additional options: \n"); + fprintf (stderr, " -pm : calculate proper motions\n"); + fprintf (stderr, " -par : calculate parallaxes\n"); + fprintf (stderr, " -pmpar : calculate motions and parallaxes\n"); + fprintf (stderr, " -region Rmin Rmax Dmin Dmax"); + fprintf (stderr, " -catalog RA DEC"); + fprintf (stderr, " -time (start)(stop)\n"); + fprintf (stderr, " +photcode (code)[,code,code...]\n"); + fprintf (stderr, " -photcode (code)[,code,code...]\n"); + fprintf (stderr, " -statmode (mode)\n"); + fprintf (stderr, " -reset\n"); + fprintf (stderr, " -minerror\n"); + fprintf (stderr, " -instmag min max\n\n"); + fprintf (stderr, " +photflags\n"); + fprintf (stderr, " -photflagbad\n"); + fprintf (stderr, " -photflagpoor\n"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " \n"); + exit (2); +} + +void usage_merge_source_id (char *name) { + + fprintf (stderr, "ERROR: invalid ID %s (remember to prefix 0x to hex IDs)\n", name); + exit (2); +} + +void usage_merge_source () { + fprintf (stderr, "ERROR: USAGE: checkastro -merge-source (objID) (catID) into (objID) (catID)\n"); + exit (2); +} Index: trunk/Ohana/src/checkastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/checkastro/src/bcatalog.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/bcatalog.c (revision 36833) @@ -0,0 +1,620 @@ +# include "checkastro.h" + +static int Nkeep1 = 0; +static int Nkeep2 = 0; +static int Nskip1 = 0; +static int Nskip2 = 0; + +// test image: 2013/06/15,13:25:51, GPC1.r.XY50 +static int CHECK_TEST_IMAGE = TRUE; +// static unsigned int Tref = 1378812312; +// static short Cref = 10001; + +static unsigned int Tref = 1379570672; +static short Cref = 10341; + +int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog); + +int bcatalog (Catalog *subcatalog, Catalog *catalog) { + + off_t i, j, offset; + off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; + int Nsecfilt; + + // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? + Nsecfilt = GetPhotcodeNsecfilt (); + assert (catalog[0].Nsecfilt == Nsecfilt); + + /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ + memset(subcatalog, 0, sizeof(Catalog)); + NAVERAGE = 50; + NMEASURE = 1000; + ALLOCATE (subcatalog[0].average, Average, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + Nmeasure = Naverage = 0; + + /* exclude stars not in range or with too few measurements */ + for (i = 0; i < catalog[0].Naverage; i++) { + if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) continue; + + /* start with all stars good */ + subcatalog[0].average[Naverage] = catalog[0].average[i]; + subcatalog[0].average[Naverage].measureOffset = Nmeasure; + for (j = 0; j < Nsecfilt; j++) { + subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j]; + } + + int foundRef = FALSE; + + Nm = 0; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + + offset = catalog[0].average[i].measureOffset + j; + + // filter objects based on user supplied criteria, including SIGMA_LIM + if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) { + catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; + if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nskip1 ++; + } + continue; + } + if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nkeep1 ++; + } + + // does this measurement match any of the reference photcodes? + if (!foundRef && NphotcodesReset) { + int m, found = FALSE; + for (m = 0; (m < NphotcodesReset) && !found; m++) { + if (catalog[0].measure[offset].photcode == photcodesReset[m][0].code) { + found = TRUE; + } + } + if (found) { + foundRef = TRUE; + } + } + + // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; + CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + subcatalog[0].measureT[Nmeasure].averef = Naverage; + + Nmeasure ++; + Nm ++; + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + } + } + if (foundRef) { + subcatalog[0].average[Naverage].flags |= 0x10; + } + subcatalog[0].average[Naverage].Nmeasure = Nm; + Naverage ++; + if (Naverage == NAVERAGE) { + NAVERAGE += 50; + REALLOCATE (subcatalog[0].average, Average, NAVERAGE); + REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + } + } + REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt); + subcatalog[0].Naverage = Naverage; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].catID = catalog[0].catID; + assert (Nsecfilt == catalog[0].Nsecfilt); + + // limit the total number of stars in the catalog + if (MaxDensityUse) { + LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog); + } else { + if (VERBOSE2) { + char *basename = filebasename (catalog[0].filename); + fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", + subcatalog[0].Naverage, catalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Nmeasure, basename); + free (basename); + } + } + if (CHECK_TEST_IMAGE && (Nkeep1 + Nkeep2 + Nskip1 + Nskip2 > 0)) { + fprintf (stderr, "kept %d %d, skipped %d %d\n", Nkeep1, Nkeep2, Nskip1, Nskip2); + } + + return (TRUE); +} + +void bcatalog_show_skips () { + fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2); + fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2); +} + +// sort by decreasing Nmeasure (X) +void sort_by_Nmeasure (int *X, off_t *Y, off_t N) { + +# define SWAPFUNC(A,B){ int tmpI; off_t tmpT; \ + tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \ + tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \ +} +# define COMPARE(A,B)(X[A] > X[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +/* generate a grid in a locally-projected coordinate system, try to select average entries + from each grid cell in decending Nmeasure order. + */ +int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) { + + off_t i, j; + int ix, iy; + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&oldcatalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&oldcatalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // generate a grid in locally projected space + double Rc = 0.5*(Rmin + Rmax); + double Dc = 0.5*(Dmin + Dmax); + + /* project coordinates to a plane centered on the object with units of arcsec */ + Coords coords; + coords.crval1 = Rc; + coords.crval2 = Dc; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "DEC--SIN"); + + // convert all average R,D values to X,Y: + double *X, *Y; + ALLOCATE (X, double, Naverage); + ALLOCATE (Y, double, Naverage); + float Xmin = +10000.0, Ymin = +10000.0; + float Xmax = -10000.0, Ymax = -10000.0; + for (i = 0; i < Naverage; i++) { + X[i] = NAN; + Y[i] = NAN; + // skip any stars which are outside of nominal catalog range + if (subcatalog[0].average[i].R < Rmin) continue; + if (subcatalog[0].average[i].R > Rmax) continue; + if (subcatalog[0].average[i].D < Dmin) continue; + if (subcatalog[0].average[i].D > Dmax) continue; + RD_to_XY (&X[i], &Y[i], subcatalog[0].average[i].R, subcatalog[0].average[i].D, &coords); + Xmin = MIN (Xmin, X[i]); + Xmax = MAX (Xmax, X[i]); + Ymin = MIN (Ymin, Y[i]); + Ymax = MAX (Ymax, Y[i]); + } + + // how many grid cells? what is the grid spacing? + float dX = Xmax - Xmin; + float dY = Ymax - Ymin; + + // *** XXX for the moment, I'm using a hard-wired cell size (200 arcsec ~ 3.3 arcmin) + int NX = (int)(dX / 200) + 1; + int NY = (int)(dY / 200) + 1; + // fprintf (stderr, "Density Grid: %d x %d\n", NX, NY); + // XXX check that NX,NY are sensible (5 degrees / 200 arcsec seems like the absolute max) + if (NX > 1000) { + fprintf (stderr, "serious problem with %s: NX = %d\n", subcatalog[0].filename, NX); + exit (3); + } + if (NY > 1000) { + fprintf (stderr, "serious problem with %s: NY = %d\n", subcatalog[0].filename, NY); + exit (3); + } + + // kind of ugly : generate a grid of index, Nmeasure arrays + // to be filled below (I also need NN and Nn to track the number of + // entries in each). + int **NN_grid; + int **Nn_grid; + int ***Nm_grid; + off_t ***idxgrid; + ALLOCATE (NN_grid, int *, NX); + ALLOCATE (Nn_grid, int *, NX); + ALLOCATE (Nm_grid, int **, NX); + ALLOCATE (idxgrid, off_t **, NX); + + for (ix = 0; ix < NX; ix++) { + ALLOCATE (NN_grid[ix], int, NY); + ALLOCATE (Nn_grid[ix], int, NY); + ALLOCATE (Nm_grid[ix], int *, NY); + ALLOCATE (idxgrid[ix], off_t *, NY); + for (iy = 0; iy < NY; iy++) { + Nn_grid[ix][iy] = 0; + NN_grid[ix][iy] = 100; + ALLOCATE (Nm_grid[ix][iy], int, NN_grid[ix][iy]); + ALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]); + } + } + + // assign all of the average entries to a grid cell + for (i = 0; i < Naverage; i++) { + if (isnan(X[i])) continue; + if (isnan(Y[i])) continue; + ix = (X[i] - Xmin) / 200.0; + iy = (Y[i] - Ymin) / 200.0; + int Nn = Nn_grid[ix][iy]; + Nm_grid[ix][iy][Nn] = subcatalog[0].average[i].Nmeasure; + + // if we are resetting to a given photcode, we need to have that photcode... + // XXX this is not relevant for the check... + if (FALSE && NphotcodesReset) { + int k; + int foundReset = FALSE; + int m = subcatalog[0].average[i].measureOffset; + MeasureTiny *measure = &subcatalog[0].measureT[m]; + for (j = 0; (j < subcatalog[0].average[i].Nmeasure) && !foundReset; j++) { + if (CHECK_TEST_IMAGE && (abs(measure[j].t - Tref) < 10) && (measure[j].photcode == Cref)) { + fprintf (stderr, "."); + } + for (k = 0; (k < NphotcodesReset) && !foundReset; k++) { + if (photcodesReset[k][0].code == measure[j].photcode) foundReset = TRUE; + } + } + if (!foundReset) { + Nm_grid[ix][iy][Nn] = 0; + } + } + + idxgrid[ix][iy][Nn] = i; + Nn_grid[ix][iy] ++; + if (Nn_grid[ix][iy] >= NN_grid[ix][iy]) { + NN_grid[ix][iy] += 100; + REALLOCATE (Nm_grid[ix][iy], int, NN_grid[ix][iy]); + REALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]); + } + } + + // sort all of the grid cells + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + sort_by_Nmeasure (Nm_grid[ix][iy], idxgrid[ix][iy], Nn_grid[ix][iy]); + NN_grid[ix][iy] = 0; // I'm going to use this array to track which element I've already selected + } + } + + // cycle over the grid until we ready Nmax + off_t *keepidx = NULL; + ALLOCATE (keepidx, off_t, Naverage); + memset (keepidx, 0, Naverage*sizeof(off_t)); + int Nkeep = 0; + + for (i = 0; (i < 20) && (Nkeep < Nmax); i++) { + for (ix = 0; (ix < NX) && (Nkeep < Nmax); ix++) { + for (iy = 0; (iy < NY) && (Nkeep < Nmax); iy++) { + if (NN_grid[ix][iy] >= Nn_grid[ix][iy]) continue; // all used up! + int Nn = NN_grid[ix][iy]; + keepidx[Nkeep] = idxgrid[ix][iy][Nn]; + Nkeep ++; + NN_grid[ix][iy] ++; + } + } + } + + // count the number of measurements this selection will yield + off_t ave, NMEASURE = 0; + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + + // test catID : 37262 37261 37257 37258 + int dumpit = FALSE; + dumpit |= (oldcatalog[0].catID == 37007); + // dumpit |= (oldcatalog[0].catID == 37261); + // dumpit |= (oldcatalog[0].catID == 37257); + // dumpit |= (oldcatalog[0].catID == 37258); + if (dumpit) { + char name[64]; + snprintf (name, 64, "cat.%05d.dump.dat", oldcatalog[0].catID); + FILE *fdump = fopen (name, "w"); + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure); + } + fclose (fdump); + } + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nkeep); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nkeep * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nkeep selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { + char *basename = filebasename (oldcatalog[0].filename); + fprintf (stderr, "limited to %d ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", + Nkeep, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, oldcatalog[0].Nmeasure, basename); + free (basename); + } + + free (X); + free (Y); + + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + free (Nm_grid[ix][iy]); + free (idxgrid[ix][iy]); + } + free (NN_grid[ix]); + free (Nn_grid[ix]); + free (Nm_grid[ix]); + free (idxgrid[ix]); + } + free (NN_grid); + free (Nn_grid); + free (Nm_grid); + free (idxgrid); + + free (keepidx); + + free (subcatalog[0].average); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nkeep; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt; + + return (TRUE); +} + +/* this version does NOT use AverageTiny, MeasureTiny */ +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&oldcatalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&oldcatalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // select a random subset of Nmax stars from subcatalog using Fisher-Yates + + // we are going to select Nmax entries by choosing the brightest objects + int *value; + off_t *index, i, j, ave; + ALLOCATE (index, off_t, Naverage); + ALLOCATE (value, int, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + value[i] = subcatalog[0].average[i].Nmeasure; + } + sort_by_Nmeasure (value, index, Naverage); + + // count the number of measurements this selection will yield + off_t NMEASURE = 0; + for (i = 0; i < Nmax; i++) { + ave = index[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + +# if (0) + if (oldcatalog[0].catID == 59962) { + FILE *fdump = fopen ("cat.dump.dat", "w"); + for (i = 0; i < Nmax; i++) { + ave = index[i]; + fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure); + } + fclose (fdump); + } +# endif + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nmax; i++) { + ave = index[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { + char *basename = filebasename (oldcatalog[0].filename); + fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", + Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, oldcatalog[0].Nmeasure, basename); + free (basename); + } + + free (index); + free (value); + free (subcatalog[0].average); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt; + + return (TRUE); +} + +/* this version does NOT use AverageTiny, MeasureTiny */ +int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE2) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // select a random subset of Nmax stars from subcatalog using Fisher-Yates + + // we are going to select Nmax entries by generating a random-sorted index list + off_t *index, tmp, i, j, ave; + ALLOCATE (index, off_t, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + } + for (i = 0; i < Naverage; i++) { + j = (Naverage - i) * drand48() + i; // a number between i and Naverage + tmp = index[j]; + index[j] = index[i]; + index[i] = tmp; + } + + // count the number of measurements this selection will yield + off_t NMEASURE = 0; + for (i = 0; i < Nmax; i++) { + ave = index[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nmax; i++) { + ave = index[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { + fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", + Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, catalog[0].filename); + } + + free (subcatalog[0].average); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt; + + return (TRUE); +} + Index: trunk/Ohana/src/checkastro/src/checkastro.c =================================================================== --- trunk/Ohana/src/checkastro/src/checkastro.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/checkastro.c (revision 36833) @@ -0,0 +1,14 @@ +# include "checkastro.h" + +int main (int argc, char **argv) { + + /* get configuration info, args */ + initialize (argc, argv); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + checkastro_images (skylist); + exit (0); +} Index: trunk/Ohana/src/checkastro/src/checkastro_client.c =================================================================== --- trunk/Ohana/src/checkastro/src/checkastro_client.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/checkastro_client.c (revision 36833) @@ -0,0 +1,36 @@ +# include "checkastro.h" + +// checkastro_client is run on a remote host and is responsible for updating the catalogs +// owned by that host. + +// there are four modes: + +// checkastro_client -load-objects : extract the bright catalog subset from the client's +// tables and save to a local FITS table for checkastro to +// load and used + +int main (int argc, char **argv) { + + // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) + initialize_client (argc, argv); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + // USAGE: checkastro_client -load-objects + int Ncatalog; + Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); + if (!catalog) { + fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); + exit (2); + } + + BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); + if (!BrightCatalogSave (BCATALOG, bcatalog)) { + fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); + exit (2); + } + + exit (0); +} Index: trunk/Ohana/src/checkastro/src/checkastro_images.c =================================================================== --- trunk/Ohana/src/checkastro/src/checkastro_images.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/checkastro_images.c (revision 36833) @@ -0,0 +1,48 @@ +# include "checkastro.h" + +int checkastro_images (SkyList *skylist) { + + int status, Ncatalog; + Catalog *catalog; + FITS_DB db; + + INITTIME; + + /* register database handle with shutdown procedure */ + set_db (&db); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, LCK_SOFT); + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + MARKTIME("load image data: %f sec\n", dtime); + + /* load regions and images based on specified sky patch (default depth) (require full overlap) */ + load_images (&db, skylist); + MARKTIME("load images: %f sec\n", dtime); + + // photcodesKeep is used here to allow measurements from the images being calibrated + // note if -reset-to-photcode is selected, photocodesKeep is replaced with below with photcodesReset + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + MARKTIME("load catalog data: %f sec\n", dtime); + + if (photcodesReset) { + photcodesKeep = photcodesReset; + NphotcodesKeep = NphotcodesReset; + } + + if (Ncatalog == 0) { + fprintf (stderr, "ERROR: no valid data for checkastro, exiting\n"); + exit (2); + } + + /* match measurements with images */ + initImageBins (catalog, Ncatalog); + MARKTIME("make image bins: %f sec\n", dtime); + + findImages (catalog, Ncatalog); + MARKTIME("set up image indexes: %f sec\n", dtime); + + exit (0); +} Index: trunk/Ohana/src/checkastro/src/initialize.c =================================================================== --- trunk/Ohana/src/checkastro/src/initialize.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/initialize.c (revision 36833) @@ -0,0 +1,29 @@ +# include "checkastro.h" + +void initialize (int argc, char **argv) { + + ConfigInit (&argc, argv); + args (argc, argv); + + if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST); + + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesReset = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset, FALSE); +} + +void initialize_client (int argc, char **argv) { + + ConfigInit (&argc, argv); + args_client (argc, argv); + + if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST); + + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesReset = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset, FALSE); +} Index: trunk/Ohana/src/checkastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/checkastro/src/load_catalogs.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/load_catalogs.c (revision 36833) @@ -0,0 +1,215 @@ +# include "checkastro.h" + +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) { + + int i, Nstar; + Catalog *catalog, *pcatalog, tcatalog; + + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + catalog = load_catalogs_parallel (skylist, Ncatalog); + return catalog; + } + + if (VERBOSE) fprintf (stderr, "loading catalog data\n"); + + ALLOCATE (catalog, Catalog, skylist[0].Nregions); + + // load data from each region file, only use bright stars + int Ncat = 0; + for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + + // we only allow output if we do not use a subset. in this case, + // the output parameters are correctly set for catalog[i] via pcatalog + pcatalog = subselect ? &tcatalog : &catalog[i]; + + // define the catalog file name + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + pcatalog->filename = hostID ? hostfile : skylist[0].filename[i]; + + // set up the basic catalog info + pcatalog->catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + pcatalog->catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + pcatalog->catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + pcatalog->Nsecfilt = GetPhotcodeNsecfilt (); + + // loads Average, Measure, SecFilt + if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE2, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename); + exit (1); + } + if (!pcatalog[0].Naves_disk) { + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); + dvo_catalog_unlock (pcatalog); + dvo_catalog_free (pcatalog); + continue; + } + + if (!pcatalog->sorted) { + fprintf (stderr, "this database is not sorted. please sort using addstar -resort\n"); + exit (1); + } + + // select only the brighter stars + if (subselect) { + // results are in Average, Secfilt, and MeasureTiny + // Ncat tracks the actually used catalogs + bcatalog (&catalog[Ncat], &tcatalog); + dvo_catalog_unlock (&tcatalog); + dvo_catalog_free (&tcatalog); + } + Ncat ++; + } + + // XXX TEST : bcatalog_show_skips(); + + Nstar = 0; + for (i = 0; i < Ncat; i++) { + Nstar += catalog[i].Naverage; + if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) { + fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n"); + fprintf (stderr, " If you need to do this, and you can afford the RAM, rebuild with IDX_T set to off_t (see checkastro.h, ImagesOps.c)\n"); + exit (3); + } + } + if (Nstar < 2) { + fprintf (stderr, "warning: insufficient stars %d\n", Nstar); + } + + // only return the populated catalogs + REALLOCATE (catalog, Catalog, Ncat); + *Ncatalog = Ncat; + return (catalog); +} + +/* this function loads all relevant catalog files for the first pass. it currently loads the data + read only (SOFT lock) since it assumes the image table has been locked. if we go to the new + addstar locking paradigm, in which the images and catalogs are updated independently, then we may + need to use an XCLD lock here. +*/ + +// CATDIR is supplied globally +# define DEBUG 1 +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + + char uniquer[12]; + int TIME = time(NULL); + int PID = getpid(); + snprintf (uniquer, 12, "%05d.%05d", PID, TIME % 100000); + + // load the list of hosts + HostTable *table = HostTableLoad (CATDIR, sky->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + exit (1); + } + + int i; + for (i = 0; i < table->Nhosts; i++) { + + // ensure that the paths are absolute path names + char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/checkastro.catalog.%s.dat", table->hosts[i].pathname, uniquer); + + char command[1024]; + snprintf (command, 1024, "checkastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D RELASTRO_SIGMA_LIM %f", + table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, SIGMA_LIM); + + char tmpline[1024]; + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + + if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_RESET_LIST) { snprintf (tmpline, 1024, "%s -reset-to-photcode %s", command, PHOTCODE_RESET_LIST); strcpy (command, tmpline); } + if (PhotFlagSelect) { snprintf (tmpline, 1024, "%s +photflags", command); strcpy (command, tmpline); } + if (PhotFlagBad) { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad); strcpy (command, tmpline); } + if (PhotFlagPoor) { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor); strcpy (command, tmpline); } + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + free (tstart); + free (tstop); + strcpy (command, tmpline); + } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running checkastro_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo); + exit (1); + } + table->hosts[i].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the checkastro_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + // each host generates a BrightCatalog structure, with the measure, average, etc value + // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt). I need to + // split out the per-catalog measurements into separate catalog entries. + + // set up an initial array of catalogs + int Nsecfilt = GetPhotcodeNsecfilt (); + CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt); + + for (i = 0; i < table->Nhosts; i++) { + + BrightCatalog *bcatalog = BrightCatalogLoad (table->hosts[i].results); + assert (bcatalog); + + BrightCatalogSplit (catalogs, bcatalog); + + free (bcatalog->average); + free (bcatalog->measure); + free (bcatalog->secfilt); + free (bcatalog); + } + + Catalog *catalog = catalogs->catalog; + *Ncatalog = catalogs->Ncatalog; + BrightCatalogSplitFree (catalogs); + + int Nmeasure = 0; + int Naverage = 0; + for (i = 0; i < catalogs->Ncatalog; i++) { + Nmeasure += catalogs->catalog[i].Nmeasure; + Naverage += catalogs->catalog[i].Naverage; + } + + fprintf (stderr, "loaded %d catalogs, using a total of %d stars (%d measures)\n", catalogs->Ncatalog, Naverage, Nmeasure); + + return (catalog); +} Index: trunk/Ohana/src/checkastro/src/load_images.c =================================================================== --- trunk/Ohana/src/checkastro/src/load_images.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/load_images.c (revision 36833) @@ -0,0 +1,34 @@ +# include "checkastro.h" + +int load_images (FITS_DB *db, SkyList *skylist) { + + Image *image, *subset; + off_t Nimage, Nsubset; + off_t *LineNumber; + + INITTIME; + + // convert database table to internal structure + image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + if (!image) { + fprintf (stderr, "ERROR: failed to read images\n"); + exit (2); + } + MARKTIME(" convert image table: %f sec\n", dtime); + + // select the images which overlap the selected sky regions + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + MARKTIME(" select images: %f sec\n", dtime); + + fprintf (stderr, "%llx\n", (long long) subset); + + initImages (subset, LineNumber, Nsubset); + MARKTIME(" init images: %f sec\n", dtime); + + // unlock, if we can (else, unlocked below) + int unlockImages = TRUE; + if (unlockImages) dvo_image_unlock (db); + + return TRUE; +} + Index: trunk/Ohana/src/checkastro/src/select_images.c =================================================================== --- trunk/Ohana/src/checkastro/src/select_images.c (revision 36833) +++ trunk/Ohana/src/checkastro/src/select_images.c (revision 36833) @@ -0,0 +1,132 @@ +# include "checkastro.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). All images in the image catalog are tested once, so there is no check that an + image already has been included. LineNum stores the locations in the Image database of the + list of images */ + +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) { + + Image *image; + off_t i, k, nimage, NIMAGE, D_NIMAGE; + off_t *line_number; + int found; + double Rexp, Dexp; + + double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + INITTIME; + + RminSkyRegion = +360.0; + RmaxSkyRegion = -360.0; + DminSkyRegion = +90.0; + DmaxSkyRegion = -90.0; + + D_NIMAGE = 6000; + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); + RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); + DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); + DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); + } + RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + MARKTIME("define sky region: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + BuildChipMatch (timage, Ntimage); + MARKTIME("build chip match: %f sec\n", dtime); + + nimage = 0; + NIMAGE = D_NIMAGE; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, off_t, NIMAGE); + + // we now have a region of interest: RminSkyRegion - RmaxSkyRegion, DminSkyRegion - DmaxSkyRegion, + // we now want to find all gpc1 chips (WRP) for which the exposure center (DIS) is in this region, + + // go through the complete list of images, selecting ones which overlap the region + for (i = 0; i < Ntimage; i++) { + + /* select images by photcode, or equiv photcode, if specified */ + if (NphotcodesKeep > 0) { + found = FALSE; + // we have to keep DIS mosaics explicitly (photcode = 0) + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) found = TRUE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == timage[i].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE; + } + if (!found) continue; + } + if (NphotcodesSkip > 0) { + found = FALSE; + for (k = 0; (k < NphotcodesSkip) && !found; k++) { + if (photcodesSkip[k][0].code == timage[i].photcode) found = TRUE; + if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE; + } + if (found) continue; + } + + /* exclude images by time */ + if (TimeSelect) { + if (timage[i].tzero < TSTART) continue; + if (timage[i].tzero > TSTOP) continue; + } + + // this adds 1.3 sec for 3M images + off_t myMosaic = FindMosaicForImage (timage, Ntimage, i); + if (!myMosaic) { + fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + myMosaic --; + + /* only check exposure center */ + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { + XY_to_RD (&Rexp, &Dexp, 0.0, 0.0, &timage[i].coords); + } else { + XY_to_RD (&Rexp, &Dexp, 0.0, 0.0, &timage[myMosaic].coords); + } + Rexp = ohana_normalize_angle_to_midpoint (Rexp, RmidSkyRegion); + + // check that this image is even in range of the searched region + if (Dexp + 1.5 > UserPatch.Dmax) continue; + if (Dexp - 1.5 < UserPatch.Dmin) continue; + + // the sky region RA is defined to be 0 - 360.0 + if ((Dexp < 88) && (Dexp > -88)) { + if (Rexp + 1.5/cos(Dexp*RAD_DEG) > UserPatch.Rmax) continue; + if (Rexp - 1.5/cos(Dexp*RAD_DEG) < UserPatch.Rmin) continue; + } + + image[nimage] = timage[i]; + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += D_NIMAGE; + D_NIMAGE = MAX (100000, D_NIMAGE * 1.5); + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); + } + } + MARKTIME("finish image selection: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, off_t, MAX (nimage, 1)); + + *Nimage = nimage; + *LineNumber = line_number; + return (image); +} + Index: trunk/Ohana/src/delstar/include/delstar.h =================================================================== --- trunk/Ohana/src/delstar/include/delstar.h (revision 36681) +++ trunk/Ohana/src/delstar/include/delstar.h (revision 36833) @@ -2,4 +2,13 @@ # include # include + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); // options for generating the IndexArray used to select images for deletion (delete_duplicate_images.c) @@ -67,4 +76,5 @@ int UPDATE; int IMAGE_DETAILS; +int IMAGE_DUPLICATES_BY_OBSTIME; int IMAGE_ONLY; int ORPHAN; Index: trunk/Ohana/src/delstar/src/args.c =================================================================== --- trunk/Ohana/src/delstar/src/args.c (revision 36681) +++ trunk/Ohana/src/delstar/src/args.c (revision 36833) @@ -18,4 +18,6 @@ fprintf (stderr, " -image-details : list info about the deleted images (-dup-images only)\n"); fprintf (stderr, " -image-only : only examine the image table (changes are NOT saved; -dup-images only)\n"); + fprintf (stderr, " -image-only-force : modify only the image table (-dup-images only)\n"); + fprintf (stderr, " -image-by-obstime : use date/time and photcode (not externID) to find duplicates\n"); fprintf (stderr, " -region Rmin Rmax Dmin Dmax : apply changes to this part of the sky\n"); @@ -68,4 +70,10 @@ if ((N = get_argument (argc, argv, "-image-details"))) { IMAGE_DETAILS = TRUE; + remove_argument (N, &argc, argv); + } + + IMAGE_DUPLICATES_BY_OBSTIME = FALSE; + if ((N = get_argument (argc, argv, "-image-by-obstime"))) { + IMAGE_DUPLICATES_BY_OBSTIME = TRUE; remove_argument (N, &argc, argv); } Index: trunk/Ohana/src/delstar/src/delete_duplicate_images.c =================================================================== --- trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 36681) +++ trunk/Ohana/src/delstar/src/delete_duplicate_images.c (revision 36833) @@ -1,3 +1,8 @@ # include "delstar.h" + +IndexArray *find_duplicates_obstime (Image *image, off_t Nimage, off_t *Nduplicates); +off_t find_obstime_range (Image *image, off_t Nimage, off_t firstEntry); +void sort_by_obstime (e_time *T, short *P, off_t *I, off_t N); +void sort_by_photcode (short *P, off_t *I, off_t N); // this function identifies the images to be deleted based on duplication of the @@ -18,5 +23,11 @@ off_t Nduplicates = 0; - IndexArray *imageID = find_duplicates (image, Nimage, &Nduplicates); + IndexArray *imageID = NULL; + + if (IMAGE_DUPLICATES_BY_OBSTIME) { + imageID = find_duplicates_obstime (image, Nimage, &Nduplicates); + } else { + imageID = find_duplicates (image, Nimage, &Nduplicates); + } if (Nduplicates == 0) { @@ -494,4 +505,209 @@ } +// sort by increasing obstime +void sort_by_obstime (e_time *T, short *P, off_t *I, off_t N) { + +# define SWAPFUNC(A,B){ e_time tmpT; short tmpP; off_t tmpI; \ + tmpT = T[A]; T[A] = T[B]; T[B] = tmpT; \ + tmpP = P[A]; P[A] = P[B]; P[B] = tmpP; \ + tmpI = I[A]; I[A] = I[B]; I[B] = tmpI; \ +} +# define COMPARE(A,B)(T[A] < T[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +// sort by increasing photcode +void sort_by_photcode (short *P, off_t *I, off_t N) { + +# define SWAPFUNC(A,B){ short tmpP; off_t tmpI; \ + tmpP = P[A]; P[A] = P[B]; P[B] = tmpP; \ + tmpI = I[A]; I[A] = I[B]; I[B] = tmpI; \ +} +# define COMPARE(A,B)(P[A] < P[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +static e_time *obstime = NULL; +static short *photcode = NULL; +static off_t *primary = NULL; +static off_t *idx = NULL; +static char *keep = NULL; + +static short *photcode_subset = NULL; +static off_t *idx_subset = NULL; + +static off_t Nsubset = 0; +static off_t NSUBSET = 300; + +off_t find_obstime_range (Image *image, off_t Nimage, off_t firstEntry) { + + Nsubset = 0; + + // find all entries with the same obstime: + e_time firstTime = obstime[firstEntry]; + + idx_subset[Nsubset] = idx[firstEntry]; + photcode_subset[Nsubset] = photcode[firstEntry]; + Nsubset++; + + off_t i = firstEntry + Nsubset; + while ((i < Nimage) && (obstime[i] == firstTime)) { + idx_subset[Nsubset] = idx[i]; + photcode_subset[Nsubset] = photcode[i]; + Nsubset++; + i++; + if (Nsubset >= NSUBSET) { + NSUBSET += 1000; + REALLOCATE (photcode_subset, short, NSUBSET); + REALLOCATE (idx_subset, off_t, NSUBSET); + } + } + + sort_by_photcode (photcode_subset, idx_subset, Nsubset); + + return i; +} + +// alternative version to find duplicates based on obstime and photcode +IndexArray *find_duplicates_obstime (Image *image, off_t Nimage, off_t *Nduplicates) { + + if (Nimage < 1) return NULL; + + // how to find duplicates: + // generate a set of arrays (idx, obstime, photcode, keep) + // sort the arrays by obstime + // loop over obstime. + // for a given new obstime, scan through the entries with the same value + // create a sub-array of photcodes, idx + // sort by photcode + // loop over entries + // find matching photcodes + // mark duplicate enties + + ALLOCATE (obstime, e_time, Nimage); + ALLOCATE (photcode, short, Nimage); + ALLOCATE (primary, off_t, Nimage); + ALLOCATE (idx, off_t, Nimage); + ALLOCATE (keep, char, Nimage); + + off_t i; + + INITTIME; + + // skip entries with photcode == 0? + for (i = 0; i < Nimage; i++) { + idx[i] = i; + primary[i] = -1; // only duplicates get a value for primary + keep[i] = TRUE; + obstime[i] = image[i].tzero; + photcode[i] = image[i].photcode; + } + MARKTIME(" generate index arrays: %f sec\n", dtime); + + // sort the 4 arrays + sort_by_obstime (obstime, photcode, idx, Nimage); + MARKTIME(" sort index arrays: %f sec\n", dtime); + + // image[idx[i]].tzero = obstime[i] + // keep[i] -> keep image[i] ('keep' is NOT resorted) + + // allocate arrays to store the subsets (these are global static) + // These get reallocated if necessary in find_obstime_range() + ALLOCATE (photcode_subset, short, NSUBSET); + ALLOCATE (idx_subset, off_t, NSUBSET); + + // entries of idx_subset correspond to the original image sequence: + // image[idx_subset[i]].tzero = obstime[i] + + off_t firstEntry = 0; + off_t nextEntry = 0; + + while (nextEntry < Nimage) { + if (firstEntry >= Nimage) { + fprintf (stderr, "error, too far?\n"); + } + + // generate photcode_subset, idx_subset in order of photcode for this unique obstime[firstEntry] + // returned value is the first value of the next entry + nextEntry = find_obstime_range (image, Nimage, firstEntry); + + // step through the photcodes and find duplicates + int j; + int firstCodeEntry = 0; + short firstCode = photcode_subset[firstCodeEntry]; + for (j = 1; j < Nsubset; j++) { + if (photcode_subset[j] == firstCode) { + // mark as duplicate + off_t dupIndex = idx_subset[j]; + keep[dupIndex] = FALSE; + primary[dupIndex] = idx_subset[firstCodeEntry]; + } else { + // new value of photcode, call if the first one + firstCodeEntry = j; + firstCode = photcode_subset[firstCodeEntry]; + } + } + firstEntry = nextEntry; + } + MARKTIME(" find duplicates: %f sec\n", dtime); + + IndexArray *imageID = make_index_array (image, Nimage, IMAGE_ID); + MARKTIME(" make index array: %f sec\n", dtime); + + // set imageID->value to TRUE for images we want to delete + off_t Ndup = 0; + for (i = 0; i < Nimage; i++) { + if (keep[i]) continue; + off_t Ni = image[i].imageID - imageID->minID; + myAssert (Ni >= 0, "oops"); + myAssert (Ni < imageID->range, "oops"); + imageID->value[Ni] = TRUE; + Ndup ++; + } + MARKTIME(" mark duplicates: %f sec\n", dtime); + + for (i = 0; IMAGE_DETAILS && (i < Nimage); i++) { + off_t Ni = image[i].imageID - imageID->minID; + if (!imageID->value[Ni]) continue; + + char *date = NULL; + date = ohana_sec_to_date (image[i].tzero); + fprintf (stderr, "delete image : (" OFF_T_FMT "), extID = %d : %30s : %20s %5d == ", i, image[i].externID, image[i].name, date, image[i].photcode); + free (date); + + off_t myPrimary = primary[i]; + if (myPrimary < 0) { + fprintf (stderr, "ERROR: this should never happen\n"); + abort(); + } + date = ohana_sec_to_date (image[myPrimary].tzero); + fprintf (stderr, "parent image : (" OFF_T_FMT "), extID = %d : %30s : %20s %5d\n", myPrimary, image[myPrimary].externID, image[myPrimary].name, date, image[myPrimary].photcode); + free (date); + } + + *Nduplicates = Ndup; + + free (photcode_subset); + free (idx_subset); + + free (obstime); + free (photcode); + free (primary); + free (idx); + free (keep); + + return imageID;; +} + // find the min & max values of the given ID (externID or imageID) // construct an empty array with length needed to fit IDs Index: trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 36681) +++ trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (revision 36833) @@ -168,4 +168,9 @@ if (!IDmap->old) { + fprintf (stderr, "input database has image IDs, but no Image table\n"); + return FALSE; + } + + if (!IDmap->Nmap) { fprintf (stderr, "input database has image IDs, but no Image table\n"); return FALSE; Index: trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c (revision 36681) +++ trunk/Ohana/src/dvomerge/src/dvomergeImageIDs.c (revision 36833) @@ -18,4 +18,6 @@ if (inDB.dbstate == LCK_EMPTY) { dvo_image_unlock (&inDB); // unlock input + IDmap->old = NULL; + IDmap->new = NULL; IDmap->Nmap = 0; return TRUE; @@ -99,4 +101,6 @@ if (inDB.dbstate == LCK_EMPTY) { dvo_image_unlock (&inDB); // unlock input + IDmap->old = NULL; + IDmap->new = NULL; IDmap->Nmap = 0; return TRUE; Index: trunk/Ohana/src/dvomerge/src/dvorepairCPT.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 36681) +++ trunk/Ohana/src/dvomerge/src/dvorepairCPT.c (revision 36833) @@ -126,5 +126,5 @@ average[Nave].measureOffset = -1; average[Nave].missingOffset = -1; - average[Nave].extendOffset = -1; + average[Nave].refColor = NAN; average[Nave].objID = measure[i].objID; Index: trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 36681) +++ trunk/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (revision 36833) @@ -411,5 +411,5 @@ average[Nave].measureOffset = -1; average[Nave].missingOffset = -1; - average[Nave].extendOffset = -1; + average[Nave].refColor = NAN; average[Nave].objID = measure[i].objID; Index: trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c =================================================================== --- trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 36681) +++ trunk/Ohana/src/dvomerge/src/dvorepairFixTables.c (revision 36833) @@ -266,5 +266,5 @@ average[Nave].measureOffset = -1; average[Nave].missingOffset = -1; - average[Nave].extendOffset = -1; + average[Nave].refColor = NAN; average[Nave].objID = measure[i].objID; Index: trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c =================================================================== --- trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 36681) +++ trunk/Ohana/src/dvopsps/src/insert_objects_dvopsps_catalog.c (revision 36833) @@ -292,5 +292,5 @@ PrintIOBuffer (ave_buffer, "%d, ", average->measureOffset); PrintIOBuffer (ave_buffer, "%d, ", average->missingOffset); - PrintIOBuffer (ave_buffer, "%d, ", average->extendOffset); + PrintIOBuffer (ave_buffer, "%d, ", 0); // was : average->extendOffset); PrintIOBuffer (ave_buffer, "%u, ", average->flags); PrintIOBuffer (ave_buffer, "%u, ", average->photFlagsUpper); Index: trunk/Ohana/src/kapa2/src/CheckPipe.c =================================================================== --- trunk/Ohana/src/kapa2/src/CheckPipe.c (revision 36681) +++ trunk/Ohana/src/kapa2/src/CheckPipe.c (revision 36833) @@ -1,3 +1,4 @@ # include "Ximage.h" +# include # define STRCONST(A) ((int)(0x1000000*A[0] + 0x10000*A[1] + 0x100*A[2] + 0x1*A[3])) @@ -41,6 +42,11 @@ word[4] = 0; switch (status) { - case -1: /* no input from pipe: continue */ - return (TRUE); + case -1: + if (errno == EAGAIN) { + /* no input from pipe: continue */ + return (TRUE); + } + perror ("exiting due to problem with socket connection in CheckPipe"); + return (FALSE); break; Index: trunk/Ohana/src/kapa2/src/EventLoop.c =================================================================== --- trunk/Ohana/src/kapa2/src/EventLoop.c (revision 36681) +++ trunk/Ohana/src/kapa2/src/EventLoop.c (revision 36833) @@ -1,3 +1,4 @@ # include "Ximage.h" +# define DEBUG 0 /* list events being selected below, all other masks are ignored */ @@ -6,9 +7,27 @@ int LastEvent (Display *display, int type, XEvent *event) { - int found; + int found, Nfound; found = FALSE; + Nfound = 0; while (XCheckTypedEvent (display, type, event)) { + // If the link is slow, then I should wait a little while for some config events to + // build up (typically, the window is being dragged around, so we get a whole series of + // config events. the flush / sync time is slow on a slow link, so we should try to wait + // for the config events to stop, but on a timescale deteremined by the flush/sync time + // if (type == ConfigureNotify) fprintf (stderr, "config (%d)\n", Nfound); + Nfound ++; found = TRUE; + } + if (DEBUG && found) { + if (type == ConfigureNotify) fprintf (stderr, "config (%d)\n", Nfound); + if (type == CirculateNotify) fprintf (stderr, "circul (%d)\n", Nfound); + if (type == Expose) fprintf (stderr, "expose (%d)\n", Nfound); + } + + // if I have a ConfigureNotify event, I should purge all Expose events as well: + if (found && (type == ConfigureNotify)) { + XEvent discard; + while (XCheckTypedEvent (display, Expose, &discard)); } return (found); @@ -38,13 +57,14 @@ if (XEventsQueued (display, QueuedAfterFlush) < 1) { - /* fprintf (stderr, "."); */ usleep (50000); continue; } + // If I have a config event, I want to also purge all expose events + /* grab the last entry for these events */ - if (LastEvent (display, ConfigureNotify, &event)) Reconfig (&event); - if (LastEvent (display, CirculateNotify, &event)) Reconfig (&event); - if (LastEvent (display, Expose, &event)) Refresh (); + if (LastEvent (display, ConfigureNotify, &event)) { Reconfig (&event); continue; } + if (LastEvent (display, CirculateNotify, &event)) { Reconfig (&event); continue; } + if (LastEvent (display, Expose, &event)) { Refresh (); continue; } if (LastEvent (display, MappingNotify, &event)) XRefreshKeyboardMapping ((XMappingEvent *) &event); if (LastEvent (display, MotionNotify, &event)) UpdatePointer (graphic, (XMotionEvent *) &event); Index: trunk/Ohana/src/kapa2/src/FlushDisplay.c =================================================================== --- trunk/Ohana/src/kapa2/src/FlushDisplay.c (revision 36681) +++ trunk/Ohana/src/kapa2/src/FlushDisplay.c (revision 36833) @@ -30,5 +30,7 @@ if (flush) { - XFlush (graphic->display); + // I changed XFlush to XSync to avoid lag + // of the display (build up of config events) + XSync (graphic->display, FALSE); reftime = now; } Index: trunk/Ohana/src/kapa2/src/Refresh.c =================================================================== --- trunk/Ohana/src/kapa2/src/Refresh.c (revision 36681) +++ trunk/Ohana/src/kapa2/src/Refresh.c (revision 36833) @@ -47,5 +47,4 @@ Y0 = graphic[0].dy - graphic[0].dy * section[0].y - dY; X0 = graphic[0].dx * section[0].x; - // fprintf (stderr, "section: (%s) %d %d - %d %d\n", section[0].name, X0, Y0, dX, dY); XFillRectangle (graphic[0].display, graphic[0].window, graphic[0].gc, X0, Y0, dX, dY); Index: trunk/Ohana/src/libautocode/def/average-ps1-dev-2.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-dev-2.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average-ps1-dev-2.d (revision 36833) @@ -25,5 +25,5 @@ FIELD measureOffset, OFF_MEASURE, int, offset to first psf measurement FIELD missingOffset, OFF_MISSING, int, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, int, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of astrometry ref stars FIELD code, code, unsigned short, ID code (star; ghost; etc) Index: trunk/Ohana/src/libautocode/def/average-ps1-v1.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-v1.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average-ps1-v1.d (revision 36833) @@ -30,5 +30,5 @@ FIELD measureOffset, OFF_MEASURE, uint32_t, offset to first psf measurement FIELD missingOffset, OFF_MISSING, uint32_t, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, uint32_t, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of astrometry ref stars # 'flags' was called 'code' prior to 2009.02.07 Index: trunk/Ohana/src/libautocode/def/average-ps1-v2.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-v2.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average-ps1-v2.d (revision 36833) @@ -35,5 +35,5 @@ FIELD measureOffset, OFF_MEASURE, uint32_t, offset to first psf measurement FIELD missingOffset, OFF_MISSING, uint32_t, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, uint32_t, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of astrometry ref stars # 'flags' was called 'code' prior to 2009.02.07 Index: trunk/Ohana/src/libautocode/def/average-ps1-v3.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-v3.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average-ps1-v3.d (revision 36833) @@ -35,5 +35,5 @@ FIELD measureOffset, OFF_MEASURE, uint32_t, offset to first psf measurement FIELD missingOffset, OFF_MISSING, uint32_t, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, uint32_t, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of astrometry ref stars # 'flags' was called 'code' prior to 2009.02.07 Index: trunk/Ohana/src/libautocode/def/average-ps1-v4.d =================================================================== --- trunk/Ohana/src/libautocode/def/average-ps1-v4.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average-ps1-v4.d (revision 36833) @@ -38,5 +38,5 @@ FIELD measureOffset, OFF_MEASURE, uint32_t, offset to first psf measurement FIELD missingOffset, OFF_MISSING, uint32_t, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, uint32_t, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of reference stars # 'flags' was called 'code' prior to 2009.02.07 Index: trunk/Ohana/src/libautocode/def/average.d =================================================================== --- trunk/Ohana/src/libautocode/def/average.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/average.d (revision 36833) @@ -38,5 +38,5 @@ FIELD measureOffset, OFF_MEASURE, int, offset to first psf measurement FIELD missingOffset, OFF_MISSING, int, offset to first missing obs -FIELD extendOffset, OFF_EXTEND, int, offset to first extended measurement +FIELD refColor, REF_COLOR, float, color of astrometry ref stars # 'flags' was called 'code' prior to 2009.02.07 Index: trunk/Ohana/src/libautocode/def/image-ps1-v2.d =================================================================== --- trunk/Ohana/src/libautocode/def/image-ps1-v2.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/image-ps1-v2.d (revision 36833) @@ -42,5 +42,5 @@ FIELD DECo, DEC_CENTER, float, image center, degrees FIELD Radius, RADIUS, float, image radius, degrees -FIELD DUMMY, DUMMY, float, dummy +FIELD refColor, REF_COLOR, float, dummy # should we define the max length of name as a macro? Index: trunk/Ohana/src/libautocode/def/image-ps1-v3.d =================================================================== --- trunk/Ohana/src/libautocode/def/image-ps1-v3.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/image-ps1-v3.d (revision 36833) @@ -42,5 +42,5 @@ FIELD DECo, DEC_CENTER, float, image center, degrees FIELD Radius, RADIUS, float, image radius, degrees -FIELD DUMMY, DUMMY, float, dummy +FIELD refColor, REF_COLOR, float, dummy # should we define the max length of name as a macro? Index: trunk/Ohana/src/libautocode/def/image-ps1-v4.d =================================================================== --- trunk/Ohana/src/libautocode/def/image-ps1-v4.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/image-ps1-v4.d (revision 36833) @@ -42,5 +42,5 @@ FIELD DECo, DEC_CENTER, float, image center, degrees FIELD Radius, RADIUS, float, image radius, degrees -FIELD DUMMY, DUMMY, float, dummy +FIELD refColor, REF_COLOR, float, dummy # should we define the max length of name as a macro? Index: trunk/Ohana/src/libautocode/def/image.d =================================================================== --- trunk/Ohana/src/libautocode/def/image.d (revision 36681) +++ trunk/Ohana/src/libautocode/def/image.d (revision 36833) @@ -44,5 +44,5 @@ FIELD DECo, DEC_CENTER, float, image center, degrees FIELD Radius, RADIUS, float, image radius, degrees -FIELD DUMMY, DUMMY, float, dummy +FIELD refColor, REF_COLOR, float, dummy FIELD name, NAME, char[121], name of original image Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 36681) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 36833) @@ -171,4 +171,5 @@ ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured + ID_STAR_BAD_PM = 0x00800000, // mean astrometry could not be measured ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) @@ -481,5 +482,5 @@ uint32_t measureOffset; // offset to first psf measurement uint32_t missingOffset; // offset to first missing obs - uint32_t extendOffset; // offset to first extended measurement + float refColor; // offset to first extended measurement uint32_t flags; // average object flags (star; ghost; etc) uint32_t photFlagsUpper; // upper bit of 2 bit summary of per-measure photflags Index: trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvodb.h (revision 36681) +++ trunk/Ohana/src/libdvo/include/dvodb.h (revision 36833) @@ -222,4 +222,5 @@ AVE_EXTID_HI, AVE_EXTID_LO, + AVE_REF_COLOR, }; @@ -289,5 +290,6 @@ IMAGE_NFIT_ASTROM, IMAGE_NLINK_PHOTOM, - IMAGE_NLINK_ASTROM + IMAGE_NLINK_ASTROM, + IMAGE_REF_COLOR }; Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 36833) @@ -161,4 +161,7 @@ case AVE_STARGAL: value.Flt = average[0].stargal; + break; + case AVE_REF_COLOR: + value.Flt = average[0].refColor; break; Index: trunk/Ohana/src/libdvo/src/dbExtractImages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 36833) @@ -69,6 +69,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -79,6 +84,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -90,6 +100,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -104,6 +119,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -118,6 +138,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -132,6 +157,11 @@ if (!haveCelestial) { if (!FindMosaicForImage (image, Nimage, N)) return value; - x = 0.5*image[N].NX; - y = 0.5*image[N].NY; + if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { + x = 0.0; + y = 0.0; + } else { + x = 0.5*image[N].NX; + y = 0.5*image[N].NY; + } XY_to_RD (&RAo, &DECo, x, y, &image[N].coords); haveCelestial = TRUE; @@ -344,4 +374,8 @@ value.Int = image[N].nLinkAstrom; break; + + case IMAGE_REF_COLOR: + value.Flt = image[N].refColor; + break; } return (value); Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 36833) @@ -381,4 +381,5 @@ if (!strcasecmp (fieldName, "EXTID_HI")) ESCAPE (AVE_EXTID_HI, MAG_NONE, OPIHI_INT); if (!strcasecmp (fieldName, "EXTID_LO")) ESCAPE (AVE_EXTID_LO, MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "REF_COLOR")) ESCAPE (AVE_REF_COLOR, MAG_NONE, OPIHI_FLT); // check for code:mode in photcode name @@ -522,8 +523,9 @@ if (!strcasecmp (fieldName, "UCDIST")) ESCAPE (IMAGE_UBERCAL_DIST,MAG_NONE, OPIHI_INT); - if (!strcasecmp (fieldName, "NFIT_PHOTOM")) ESCAPE (IMAGE_NFIT_PHOTOM,MAG_NONE, OPIHI_INT); - if (!strcasecmp (fieldName, "NFIT_ASTROM")) ESCAPE (IMAGE_NFIT_ASTROM,MAG_NONE, OPIHI_INT); - if (!strcasecmp (fieldName, "NLINK_PHOTOM")) ESCAPE (IMAGE_NLINK_PHOTOM,MAG_NONE, OPIHI_INT); - if (!strcasecmp (fieldName, "NLINK_ASTROM")) ESCAPE (IMAGE_NLINK_ASTROM,MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "NFIT_PHOTOM")) ESCAPE (IMAGE_NFIT_PHOTOM, MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "NFIT_ASTROM")) ESCAPE (IMAGE_NFIT_ASTROM, MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "NLINK_PHOTOM")) ESCAPE (IMAGE_NLINK_PHOTOM, MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "NLINK_ASTROM")) ESCAPE (IMAGE_NLINK_ASTROM, MAG_NONE, OPIHI_INT); + if (!strcasecmp (fieldName, "REF_COLOR")) ESCAPE (IMAGE_REF_COLOR, MAG_NONE, OPIHI_FLT); // for words that don't parse, try a photcode Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 36833) @@ -100,5 +100,5 @@ average->measureOffset = -1; average->missingOffset = -1; - average->extendOffset = -1; + average->refColor = NAN; average->flags = 0; Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 36833) @@ -119,5 +119,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].objID = in[i].objID; out[i].catID = in[i].catID; @@ -153,5 +153,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].objID = in[i].objID; out[i].catID = in[i].catID; Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 36833) @@ -140,5 +140,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].objID = in[i].objID; @@ -175,5 +175,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].objID = in[i].objID; Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 36833) @@ -145,5 +145,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].objID = in[i].objID; @@ -184,5 +184,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].objID = in[i].objID; @@ -269,4 +269,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; @@ -331,4 +332,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 36833) @@ -145,5 +145,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].photFlagsUpper = in[i].photFlagsUpper; @@ -186,5 +186,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].photFlagsUpper = in[i].photFlagsUpper; @@ -279,4 +279,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; @@ -341,4 +342,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 36833) @@ -166,5 +166,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].photFlagsUpper = in[i].photFlagsUpper; @@ -210,5 +210,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].photFlagsUpper = in[i].photFlagsUpper; @@ -317,4 +317,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; @@ -379,4 +380,5 @@ out[i].DECo = in[i].DECo; out[i].Radius = in[i].Radius; + out[i].refColor = in[i].refColor; out[i].detection_limit = in[i].detection_limit; @@ -680,5 +682,5 @@ out[i].measureOffset = in[i].measureOffset; out[i].missingOffset = in[i].missingOffset; - out[i].extendOffset = in[i].extendOffset; + out[i].refColor = in[i].refColor; out[i].flags = in[i].flags; out[i].photFlagsUpper = in[i].photFlagsUpper; Index: trunk/Ohana/src/libdvo/src/skyregion_io.c =================================================================== --- trunk/Ohana/src/libdvo/src/skyregion_io.c (revision 36681) +++ trunk/Ohana/src/libdvo/src/skyregion_io.c (revision 36833) @@ -78,4 +78,5 @@ gfits_free_matrix (&matrix); gfits_free_header (&theader); + fclose (f); return (skytable); Index: trunk/Ohana/src/libkapa/src/KapaOpen.c =================================================================== --- trunk/Ohana/src/libkapa/src/KapaOpen.c (revision 36681) +++ trunk/Ohana/src/libkapa/src/KapaOpen.c (revision 36833) @@ -218,15 +218,33 @@ // apparently, I can connect on someone else's port (eg GoogleTalkPlugin) // do a simple handshake before we set !NONBLOCK: - char line[5]; - int Nout = read (InitSocket, line, 4); - if (Nout != 4) { - if (DEBUG) fprintf (stderr, "connection failed\n"); - close (InitSocket); - return (-1); - } - if (strncmp (line, "KAPA", 4)) { - if (DEBUG) fprintf (stderr, "connection to the wrong server\n"); - close (InitSocket); - return (-1); + + // ensure the socket is NONBLOCK first... + fcntl (InitSocket, F_SETFL, O_NONBLOCK); + + int Ntry = 0; +retry_message: + { + char line[5]; + int Nout = read (InitSocket, line, 4); + if ((Nout == -1) && (errno == EAGAIN)) { + Ntry ++; + if (Ntry > 500) { + if (DEBUG) fprintf (stderr, "handshake failure\n"); + close (InitSocket); + return (-1); + } + usleep (10000); + goto retry_message; + } + if (Nout != 4) { + if (DEBUG) fprintf (stderr, "connection failed\n"); + close (InitSocket); + return (-1); + } + if (strncmp (line, "KAPA", 4)) { + if (DEBUG) fprintf (stderr, "connection to the wrong server\n"); + close (InitSocket); + return (-1); + } } @@ -265,11 +283,13 @@ } -# define NTRY 500 +# define NTRY 1000 Ntry = 0; while (Ntry < NTRY) { sock = KapaClientSocket ("localhost"); if (sock != -1) break; - // if (errno != EAGAIN) break; - if (errno != ECONNREFUSED) break; + if (errno != ECONNREFUSED) { + perror ("KapaOpen"); + break; + } // no connection yet. try again, but first check // if the kapa job has exited @@ -330,5 +350,5 @@ fcntl (InitSocket, F_SETFL, O_NONBLOCK); -# define NTRY 500 +# define NTRY 1000 Ntry = 0; while (Ntry < NTRY) { Index: trunk/Ohana/src/libkapa/src/KiiPicture.c =================================================================== --- trunk/Ohana/src/libkapa/src/KiiPicture.c (revision 36681) +++ trunk/Ohana/src/libkapa/src/KiiPicture.c (revision 36833) @@ -21,5 +21,6 @@ int KiiNewPicture1D (int fd, KiiImage *image, KapaImageData *data, Coords *coords) { - int Nwrite, Npix, Ncolors, size; + int Ncolors; + off_t Nwrite, Npix, size; float min, max; @@ -36,5 +37,5 @@ KiiSendMessage (fd, "%8d %8d", image[0].Nx, image[0].Ny); KiiSendMessage (fd, "%f %f %s %s", data[0].zero, data[0].range, data[0].name, data[0].file); - KiiSendMessage (fd, "%f %f %d ", min, max, size); + KiiSendMessage (fd, "%f %f "OFF_T_FMT" ", min, max, size); KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval1, coords[0].crpix1, coords[0].cdelt1, coords[0].pc1_1, coords[0].pc1_2); KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval2, coords[0].crpix2, coords[0].cdelt2, coords[0].pc2_1, coords[0].pc2_2); @@ -42,8 +43,18 @@ /* send the image data */ - Nwrite = write (fd, image[0].data1d, size); - if (Nwrite != size) { - fprintf (stderr, "error sending picture to kapa\n"); - return (FALSE); + off_t bytes_left = size; + while (bytes_left > 0) { + Nwrite = write (fd, image[0].data1d, bytes_left); + if (Nwrite == 0) { + fprintf (stderr, "unable to send more data to kapa?\n"); + return FALSE; + } + if (Nwrite < 0) { + if (errno == EAGAIN) continue; + if (errno == EWOULDBLOCK) continue; + perror ("KiiNewPicture1D:"); + return FALSE; + } + bytes_left -= Nwrite; } @@ -114,6 +125,6 @@ int KiiNewPicture2D (int fd, KiiImage *image, KapaImageData *data, Coords *coords) { - int j; - int Nwrite, Nbytes, Npix, Ncolors, size; + int Ncolors; + off_t j, Nwrite, Nbytes, Npix, size; float min, max; @@ -130,5 +141,5 @@ KiiSendMessage (fd, "%8d %8d", image[0].Nx, image[0].Ny); KiiSendMessage (fd, "%f %f %s %s", data[0].zero, data[0].range, data[0].name, data[0].file); - KiiSendMessage (fd, "%f %f %d ", min, max, size); + KiiSendMessage (fd, "%f %f "OFF_T_FMT" ", min, max, size); KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval1, coords[0].crpix1, coords[0].cdelt1, coords[0].pc1_1, coords[0].pc1_2); KiiSendMessage (fd, "%f %f %g %g %g ", coords[0].crval2, coords[0].crpix2, coords[0].cdelt2, coords[0].pc2_1, coords[0].pc2_2); Index: trunk/Ohana/src/opihi/dvo/Makefile =================================================================== --- trunk/Ohana/src/opihi/dvo/Makefile (revision 36681) +++ trunk/Ohana/src/opihi/dvo/Makefile (revision 36833) @@ -41,4 +41,5 @@ $(SRC)/cmpread.$(ARCH).o \ $(SRC)/coordimage.$(ARCH).o \ +$(SRC)/coordmosaic.$(ARCH).o \ $(SRC)/ddmags.$(ARCH).o \ $(SRC)/detrend.$(ARCH).o \ Index: trunk/Ohana/src/opihi/dvo/coordmosaic.c =================================================================== --- trunk/Ohana/src/opihi/dvo/coordmosaic.c (revision 36833) +++ trunk/Ohana/src/opihi/dvo/coordmosaic.c (revision 36833) @@ -0,0 +1,524 @@ +# include "dvoshell.h" + +int coordmosaic (int argc, char **argv) { + + int N, ix, iy; + off_t i, j, nt, Nimage, *subset, Nsubset; + time_t tzero; + double trange; + Buffer *bufX, *bufY; + SkyRegionSelection *selection; + Image *image; + + if (!InitPhotcodes ()) return (FALSE); + + // parse skyregion options + if ((selection = SetRegionSelection (&argc, argv)) == NULL) { + gprint (GP_ERR, "invalid sky region selection\n"); + return FALSE; + } + + int TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!ohana_str_to_time (argv[N], &tzero)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!ohana_str_to_dtime (argv[N], &trange)) { + gprint (GP_ERR, "syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + gprint (GP_ERR, "plotting in range %ds - %ds (%f seconds)\n", (int)tzero, (int)(tzero + trange), trange); + } + + /* + PhotCode *Photcode = NULL; + if ((N = get_argument (argc, argv, "-photcode"))) { + remove_argument (N, &argc, argv); + Photcode = GetPhotcodebyName (argv[N]); + if (Photcode == NULL) { + gprint (GP_ERR, "photcode not found in photcode table\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-Nphotcode"))) { + remove_argument (N, &argc, argv); + Photcode = GetPhotcodebyCode (atoi(argv[N])); + if (Photcode == NULL) { + gprint (GP_ERR, "photcode not found in photcode table\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + } + */ + + char *Name = NULL; + if ((N = get_argument (argc, argv, "-name"))) { + remove_argument (N, &argc, argv); + Name = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + /* + int chipID = -1; + if ((N = get_argument (argc, argv, "-chip"))) { + remove_argument (N, &argc, argv); + chipID = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + */ + + // save the image plate scales + Vector *cd1v = NULL, *cd2v = NULL; + if ((cd1v = SelectVector ("cd1", ANYVECTOR, TRUE)) == NULL) return (FALSE); + if ((cd2v = SelectVector ("cd2", ANYVECTOR, TRUE)) == NULL) return (FALSE); + + if (argc != 5) goto syntax; + + if ((bufX = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) goto escape; + if ((bufY = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) goto escape; + int Nx = atoi (argv[3]); + int Ny = atoi (argv[4]); + + gfits_free_matrix (&bufX[0].matrix); + gfits_free_header (&bufX[0].header); + CreateBuffer (bufX, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + strcpy (bufX[0].file, "(empty)"); + + gfits_free_matrix (&bufY[0].matrix); + gfits_free_header (&bufY[0].header); + CreateBuffer (bufY, Nx, Ny, -32, 0.0, 1.0); // initialized to 0.0 here + strcpy (bufY[0].file, "(empty)"); + + if ((image = LoadImagesDVO (&Nimage)) == NULL) return (FALSE); + image_subset (image, Nimage, &subset, &Nsubset, selection, tzero, trange, TimeSelect); + + float *vX = (float *) bufX[0].matrix.buffer; + float *vY = (float *) bufY[0].matrix.buffer; + + int Nsum = 0; + int NSUM = 1000; + opihi_flt *cd1 = NULL; + opihi_flt *cd2 = NULL; + ALLOCATE (cd1, opihi_flt, NSUM); + ALLOCATE (cd2, opihi_flt, NSUM); + + int NX = 0; + int NY = 0; + float dX = NAN; + float dY = NAN; + + // only select the DIS (mosaic) images + for (j = 0; j < Nsubset; j++) { + i = subset[j]; + if (Name && (strstr (image[i].name, Name) == (char *) NULL)) continue; + + if (strcmp (&image[i].coords.ctype[4], "-DIS")) continue; + + if (!isfinite(dX)) { + NX = image[i].NX; + NY = image[i].NY; + dX = NX / (float) Nx; + dY = NY / (float) Ny; + } else { + if ((NX != image[i].NX) || (NY != image[i].NY)) { + // fprintf (stderr, "mosaic size mismatch? %d,%d vs %d,%d\n", NX, NY, image[i].NX, image[i].NY); + } + } + + // DIS image pixels range fom -0.5*NX to +0.5*NX + // output pixels range from 0 to Nx + // dX = NX / Nx + // fx = ix * dX - 0.5*NX + + Coords coords = image[i].coords; + // coords.Npolyterms = 0; + + // the coords have an arbitrary rotation. I want derotate: + double T = 1.0 / (image[i].coords.pc1_1*image[i].coords.pc2_2 - image[i].coords.pc2_1*image[i].coords.pc1_2); + double A10 = +image[i].coords.pc2_2 * T; + double A01 = -image[i].coords.pc1_2 * T; + double B10 = -image[i].coords.pc2_1 * T; + double B01 = +image[i].coords.pc1_1 * T; + + coords.pc1_1 = image[i].coords.pc1_1 * A10 + image[i].coords.pc2_1 * A01; + coords.pc1_2 = image[i].coords.pc1_2 * A10 + image[i].coords.pc2_2 * A01; + coords.pc2_1 = image[i].coords.pc1_1 * B10 + image[i].coords.pc2_1 * B01; + coords.pc2_2 = image[i].coords.pc1_2 * B10 + image[i].coords.pc2_2 * B01; + + if (coords.Npolyterms > 1) { + for (nt = 0; nt < 3; nt++) { + coords.polyterms[nt][0] = image[i].coords.polyterms[nt][0] * A10 + image[i].coords.polyterms[nt][1] * A01; + coords.polyterms[nt][1] = image[i].coords.polyterms[nt][0] * B10 + image[i].coords.polyterms[nt][1] * B01; + } + } + + if (coords.Npolyterms > 2) { + for (nt = 3; nt < 7; nt++) { + coords.polyterms[nt][0] = image[i].coords.polyterms[nt][0] * A10 + image[i].coords.polyterms[nt][1] * A01; + coords.polyterms[nt][1] = image[i].coords.polyterms[nt][0] * B10 + image[i].coords.polyterms[nt][1] * B01; + } + } + + Coords linear = coords; + linear.Npolyterms = 0; + + for (iy = 0; iy < Ny; iy++) { + double fy = iy * dY - 0.5*NY; // coordinate in full image + for (ix = 0; ix < Nx; ix++) { + double fx = ix * dX - 0.5*NX; // coordinate in full image + + double Lo, Lm; //, dL; + double Mo, Mm; //, dM; + + XY_to_LM (&Lo, &Mo, fx, fy, &coords); + XY_to_LM (&Lm, &Mm, fx, fy, &linear); + // dL = Lo - Lm; + // dM = Mo - Mm; + + // vX[ix + iy*Nx] += dL; + // vY[ix + iy*Nx] += dM; + vX[ix + iy*Nx] += Lo; + vY[ix + iy*Nx] += Mo; + } + } + + cd1[Nsum] = image[i].coords.cdelt1; + cd2[Nsum] = image[i].coords.cdelt2; + Nsum ++; + if (Nsum == NSUM) { + NSUM += 1000; + REALLOCATE (cd1, opihi_flt, NSUM); + REALLOCATE (cd2, opihi_flt, NSUM); + } + } + + free (cd1v->elements.Flt); cd1v->elements.Flt = cd1; cd1v->Nelements = Nsum; + free (cd2v->elements.Flt); cd2v->elements.Flt = cd2; cd2v->Nelements = Nsum; + + for (iy = 0; iy < Ny; iy++) { + for (ix = 0; ix < Nx; ix++) { + vX[ix + iy*Nx] /= (float) Nsum; + vY[ix + iy*Nx] /= (float) Nsum; + } + } + + gfits_modify (&bufX->header, "MOS_NX", "%d", 1, NX); + gfits_modify (&bufX->header, "MOS_NY", "%d", 1, NY); + gfits_modify (&bufX->header, "MOS_DX", "%f", 1, dX); + gfits_modify (&bufX->header, "MOS_DY", "%f", 1, dY); + gfits_modify (&bufX->header, "MOS_NSUM", "%d", 1, Nsum); + + fprintf (stderr, "Nsum: %d\n", Nsum); + return TRUE; + + syntax: + gprint (GP_ERR, "USAGE: coordmosaic buffX buffY Nx Ny\n"); + escape: + return (FALSE); +} + +# define ESCAPE(MSG) { \ + if (src) fclose (src); \ + if (tgt) fclose (tgt); \ + fprintf (stderr, "%s\n", MSG); return FALSE; } + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&header, &table, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +int psastro_model (int argc, char **argv) { + + // I need to write out the full psastro input model based on my new fits + + // astrometry model consists of an empty PHU + 4 extensions: + // CHIPS : + // FP + // TP + // SKY + + int Ncol; + off_t Nrow, Nout; + char type[16]; + + FTable table; + Header header; + table.header = &header; + + FILE *src = NULL, *tgt = NULL; + + if (argc != 3) { + fprintf (stderr, "USAGE: psastro_model (src) (tgt)\n"); + return FALSE; + } + + // I am going to read in the tables from one file and save them (with updates) in another + src = fopen (argv[1], "r"); + if (!src) ESCAPE ("failed to open source file"); + + tgt = fopen (argv[2], "w"); + if (!tgt) ESCAPE ("failed to open output file"); + + // read the PHU + if (!gfits_load_header (src, &header)) ESCAPE ("error loading PHU"); + int Nbytes = gfits_data_size (&header); + if (Nbytes) ESCAPE ("non-zero PHU matrix??"); + + // save the PHU + if (!gfits_save_header (tgt, &header)) ESCAPE ("error saving PHU"); + gfits_free_header (&header); + + { + // read the CHIPS + if (!gfits_load_header (src, &header)) ESCAPE ("error loading CHIPS header"); + if (!gfits_fread_ftable_data (src, &table, FALSE)) ESCAPE ("error reading table for CHIPS"); + + // these set the value of Nrow + GET_COLUMN (SEGMENT , "SEGMENT", char ); int Nsegment = Ncol; + GET_COLUMN (PARENT , "PARENT" , char ); + GET_COLUMN (MINX , "MINX" , float); + GET_COLUMN (MAXX , "MAXX" , float); + GET_COLUMN (MINY , "MINY" , float); + GET_COLUMN (MAXY , "MAXY" , float); + GET_COLUMN (XORDER , "XORDER" , int ); + GET_COLUMN (YORDER , "YORDER" , int ); + GET_COLUMN (NXORDER , "NXORDER", int ); + GET_COLUMN (NYORDER , "NYORDER", int ); + GET_COLUMN (POLY_X , "POLY_X" , float); + GET_COLUMN (POLY_Y , "POLY_Y" , float); + GET_COLUMN (ERROR_X , "ERROR_X", float); + GET_COLUMN (ERROR_Y , "ERROR_Y", float); + GET_COLUMN (MASK_X , "MASK_X" , char ); + GET_COLUMN (MASK_Y , "MASK_Y" , char ); + + // I am keeping the same number of entries, but changing the values of + // POLY_X, POLY_Y + int i; + for (i = 0; i < Nrow; i++) { + // validate that NX, NY orders are 1 + assert (NXORDER[i] == 1); + assert (NYORDER[i] == 1); + + char name[16], segment[16]; + + int found, xmask, ymask; + + memset (segment, 0, 16); memcpy (segment, &SEGMENT[i*Nsegment], 4); + snprintf (name, 16, "%s_L_CX%dY%d", segment, XORDER[i], YORDER[i]); + double xvalue = get_double_variable (name, &found); + xmask = !found; + + memset (segment, 0, 16); memcpy (segment, &SEGMENT[i*Nsegment], 4); + snprintf (name, 16, "%s_M_CX%dY%d", segment, XORDER[i], YORDER[i]); + double yvalue = get_double_variable (name, &found); + ymask = !found; + + POLY_X[i] = xvalue; + POLY_Y[i] = yvalue; + ERROR_X[i] = 0.0; + ERROR_Y[i] = 0.0; + MASK_X[i] = xmask; + MASK_Y[i] = ymask; + } + + // reset the table size to 0 rows -- this forces the data array to be reallocate by gfits_set_bintable_column + gfits_modify (&header, "NAXIS2", "%d", 1, 0); + + gfits_set_bintable_column (&header, &table, "SEGMENT", SEGMENT , Nrow); + gfits_set_bintable_column (&header, &table, "PARENT" , PARENT , Nrow); + gfits_set_bintable_column (&header, &table, "MINX" , MINX , Nrow); + gfits_set_bintable_column (&header, &table, "MAXX" , MAXX , Nrow); + gfits_set_bintable_column (&header, &table, "MINY" , MINY , Nrow); + gfits_set_bintable_column (&header, &table, "MAXY" , MAXY , Nrow); + gfits_set_bintable_column (&header, &table, "XORDER" , XORDER , Nrow); + gfits_set_bintable_column (&header, &table, "YORDER" , YORDER , Nrow); + gfits_set_bintable_column (&header, &table, "NXORDER", NXORDER , Nrow); + gfits_set_bintable_column (&header, &table, "NYORDER", NYORDER , Nrow); + gfits_set_bintable_column (&header, &table, "POLY_X" , POLY_X , Nrow); + gfits_set_bintable_column (&header, &table, "POLY_Y" , POLY_Y , Nrow); + gfits_set_bintable_column (&header, &table, "ERROR_X", ERROR_X , Nrow); + gfits_set_bintable_column (&header, &table, "ERROR_Y", ERROR_Y , Nrow); + gfits_set_bintable_column (&header, &table, "MASK_X" , MASK_X , Nrow); + gfits_set_bintable_column (&header, &table, "MASK_Y" , MASK_Y , Nrow); + + free (SEGMENT ); + free (PARENT ); + free (MINX ); + free (MAXX ); + free (MINY ); + free (MAXY ); + free (XORDER ); + free (YORDER ); + free (NXORDER ); + free (NYORDER ); + free (POLY_X ); + free (POLY_Y ); + free (ERROR_X ); + free (ERROR_Y ); + free (MASK_X ); + free (MASK_Y ); + + // save the CHIPS + gfits_fwrite_Theader (tgt, &header); + gfits_fwrite_table (tgt, &table); + gfits_free_header (&header); + gfits_free_table (&table); + } + + { + // read the FP + if (!gfits_load_header (src, &header)) ESCAPE ("error loading FP header"); + if (!gfits_fread_ftable_data (src, &table, FALSE)) ESCAPE ("error reading table for FP"); + + // these set the value of Nrow + GET_COLUMN (SEGMENT , "SEGMENT", char ); int Nsegment = Ncol; + GET_COLUMN (PARENT , "PARENT" , char ); int Nparent = Ncol; + GET_COLUMN (MINX , "MINX" , float); + GET_COLUMN (MAXX , "MAXX" , float); + GET_COLUMN (MINY , "MINY" , float); + GET_COLUMN (MAXY , "MAXY" , float); + GET_COLUMN (XORDER , "XORDER" , int ); + GET_COLUMN (YORDER , "YORDER" , int ); + GET_COLUMN (NXORDER , "NXORDER", int ); + GET_COLUMN (NYORDER , "NYORDER", int ); + GET_COLUMN (POLY_X , "POLY_X" , float); + GET_COLUMN (POLY_Y , "POLY_Y" , float); + GET_COLUMN (ERROR_X , "ERROR_X", float); + GET_COLUMN (ERROR_Y , "ERROR_Y", float); + GET_COLUMN (MASK_X , "MASK_X" , char ); + GET_COLUMN (MASK_Y , "MASK_Y" , char ); + + Nout = 16; + + REALLOCATE (SEGMENT , char , Nout*Nsegment); + REALLOCATE (PARENT , char , Nout*Nparent); + REALLOCATE (MINX , float, Nout); + REALLOCATE (MAXX , float, Nout); + REALLOCATE (MINY , float, Nout); + REALLOCATE (MAXY , float, Nout); + REALLOCATE (XORDER , int , Nout); + REALLOCATE (YORDER , int , Nout); + REALLOCATE (NXORDER , int , Nout); + REALLOCATE (NYORDER , int , Nout); + REALLOCATE (POLY_X , float, Nout); + REALLOCATE (POLY_Y , float, Nout); + REALLOCATE (ERROR_X , float, Nout); + REALLOCATE (ERROR_Y , float, Nout); + REALLOCATE (MASK_X , char , Nout); + REALLOCATE (MASK_Y , char , Nout); + + float minX = MINX[0]; + float maxX = MAXX[0]; + float minY = MINY[0]; + float maxY = MAXY[0]; + + int ix, iy, nrow; + nrow = 0; + for (ix = 0; ix < 4; ix++) { + for (iy = 0; iy < 4; iy++) { + memcpy (&SEGMENT[nrow*Nsegment], "FOCAL_PLANE", strlen("FOCAL_PLANE")); + memcpy (&PARENT[nrow*Nparent], "TANGENT_PLANE", strlen("TANGENT_PLANE")); + MINX[nrow] = minX; + MAXX[nrow] = maxX; + MINY[nrow] = minY; + MAXY[nrow] = maxY; + XORDER[nrow] = ix; + YORDER[nrow] = iy; + NXORDER[nrow] = 3; + NYORDER[nrow] = 3; + + char name[16]; + + int found, xmask, ymask; + snprintf (name, 16, "L_X%dY%d", ix, iy); + double xvalue = get_double_variable (name, &found); + xmask = !found; + + snprintf (name, 16, "M_X%dY%d", ix, iy); + double yvalue = get_double_variable (name, &found); + ymask = !found; + + POLY_X[nrow] = xvalue; + POLY_Y[nrow] = yvalue; + ERROR_X[nrow] = 0.0; + ERROR_Y[nrow] = 0.0; + MASK_X[nrow] = xmask; + MASK_Y[nrow] = ymask; + nrow ++; + } + } + + // reset the table size to 0 rows -- this forces the data array to be reallocate by gfits_set_bintable_column + gfits_modify (&header, "NAXIS2", "%d", 1, 0); + + gfits_set_bintable_column (&header, &table, "SEGMENT", SEGMENT , Nout); + gfits_set_bintable_column (&header, &table, "PARENT" , PARENT , Nout); + gfits_set_bintable_column (&header, &table, "MINX" , MINX , Nout); + gfits_set_bintable_column (&header, &table, "MAXX" , MAXX , Nout); + gfits_set_bintable_column (&header, &table, "MINY" , MINY , Nout); + gfits_set_bintable_column (&header, &table, "MAXY" , MAXY , Nout); + gfits_set_bintable_column (&header, &table, "XORDER" , XORDER , Nout); + gfits_set_bintable_column (&header, &table, "YORDER" , YORDER , Nout); + gfits_set_bintable_column (&header, &table, "NXORDER", NXORDER , Nout); + gfits_set_bintable_column (&header, &table, "NYORDER", NYORDER , Nout); + gfits_set_bintable_column (&header, &table, "POLY_X" , POLY_X , Nout); + gfits_set_bintable_column (&header, &table, "POLY_Y" , POLY_Y , Nout); + gfits_set_bintable_column (&header, &table, "ERROR_X", ERROR_X , Nout); + gfits_set_bintable_column (&header, &table, "ERROR_Y", ERROR_Y , Nout); + gfits_set_bintable_column (&header, &table, "MASK_X" , MASK_X , Nout); + gfits_set_bintable_column (&header, &table, "MASK_Y" , MASK_Y , Nout); + + free (SEGMENT ); + free (PARENT ); + free (MINX ); + free (MAXX ); + free (MINY ); + free (MAXY ); + free (XORDER ); + free (YORDER ); + free (NXORDER ); + free (NYORDER ); + free (POLY_X ); + free (POLY_Y ); + free (ERROR_X ); + free (ERROR_Y ); + free (MASK_X ); + free (MASK_Y ); + + // save the FP + gfits_fwrite_Theader (tgt, &header); + gfits_fwrite_table (tgt, &table); + gfits_free_header (&header); + gfits_free_table (&table); + } + + // read the TP + if (!gfits_load_header (src, &header)) ESCAPE ("error loading TP header"); + if (!gfits_fread_ftable_data (src, &table, FALSE)) ESCAPE ("error reading table for TP"); + + // save the TP + gfits_fwrite_Theader (tgt, &header); + gfits_fwrite_table (tgt, &table); + gfits_free_header (&header); + gfits_free_table (&table); + + // read the SKY + if (!gfits_load_header (src, &header)) ESCAPE ("error loading SKY header"); + if (!gfits_fread_ftable_data (src, &table, FALSE)) ESCAPE ("error reading table for SKY"); + + // save the SKY + gfits_fwrite_Theader (tgt, &header); + gfits_fwrite_table (tgt, &table); + gfits_free_header (&header); + gfits_free_table (&table); + + fclose (tgt); + fclose (src); + + return TRUE; +} Index: trunk/Ohana/src/opihi/dvo/init.c =================================================================== --- trunk/Ohana/src/opihi/dvo/init.c (revision 36681) +++ trunk/Ohana/src/opihi/dvo/init.c (revision 36833) @@ -14,4 +14,6 @@ int cmpread PROTO((int, char **)); int coordimage PROTO((int, char **)); +int coordmosaic PROTO((int, char **)); +int psastro_model PROTO((int, char **)); int ddmags PROTO((int, char **)); int detrend PROTO((int, char **)); @@ -73,4 +75,6 @@ {1, "cmpread", cmpread, "read data from cmp format files"}, {1, "coordimage", coordimage, "generate a map of the transformation residuals"}, + {1, "coordmosaic", coordmosaic, "generate a map of the distortion"}, + {1, "psastro_model", psastro_model, "save psastro-format astrometry model"}, {1, "ddmags", ddmags, "plot magnitude differences"}, {1, "detrend", detrend, "extract from detrend database?"}, Index: trunk/Ohana/src/photdbc/Makefile =================================================================== --- trunk/Ohana/src/photdbc/Makefile (revision 36681) +++ trunk/Ohana/src/photdbc/Makefile (revision 36833) @@ -31,4 +31,5 @@ $(SRC)/args.$(ARCH).o \ $(SRC)/copy_images.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/join_stars.$(ARCH).o \ @@ -41,5 +42,4 @@ $(SRC)/ConfigInit.$(ARCH).o \ $(SRC)/args.$(ARCH).o \ -$(SRC)/copy_images.$(ARCH).o \ $(SRC)/Shutdown.$(ARCH).o \ $(SRC)/join_stars.$(ARCH).o \ Index: trunk/Ohana/src/photdbc/include/photdbc.h =================================================================== --- trunk/Ohana/src/photdbc/include/photdbc.h (revision 36681) +++ trunk/Ohana/src/photdbc/include/photdbc.h (revision 36833) @@ -20,4 +20,20 @@ int Nmeas; } StatType; + +typedef struct { + double Xc[5]; + double Yc[5]; + double Rc; + double Dc; +} SkyRegionCoords; + +# define MARKTIME(MSG,...) { \ + gettimeofday (&stopTimer, (void *) NULL); \ + float dtime = DTIME (stopTimer, startTimer); \ + fprintf (stderr, MSG, __VA_ARGS__); } + +# define INITTIME \ + struct timeval startTimer, stopTimer; \ + gettimeofday (&startTimer, (void *) NULL); int PARALLEL; @@ -80,7 +96,7 @@ PhotCodeData photcodes; -char *PHOTCODE_DROP_LIST, *PHOTCODE_SKIP_LIST; -int NphotcodesDrop, NphotcodesSkip; -PhotCode **photcodesDrop, **photcodesSkip; +char *PHOTCODE_DROP_LIST, *PHOTCODE_KEEP_LIST; +int NphotcodesDrop, NphotcodesKeep; +PhotCode **photcodesDrop, **photcodesKeep; # define FLAG_AREA 0X0001 @@ -131,6 +147,11 @@ void SetProtect (int mode); int SetSignals (void); -int copy_images (char *outdir); +int copy_images (char *outdir, SkyList *skylist); void usage(); + +void dsortindex (double *X, off_t *Y, int N); +off_t getRegionStartByRA (double R, double *Rref, off_t Nregions); + +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap); int photdbc_catalogs (char *outroot, SkyList *skylist, int hostID); Index: trunk/Ohana/src/photdbc/src/args.c =================================================================== --- trunk/Ohana/src/photdbc/src/args.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/args.c (revision 36833) @@ -85,4 +85,5 @@ } + // measurements with these photcodes are not copied to the output PHOTCODE_DROP_LIST = NULL; if ((N = get_argument (argc, argv, "-photcode-drop"))) { @@ -92,8 +93,11 @@ } - PHOTCODE_SKIP_LIST = NULL; - if ((N = get_argument (argc, argv, "-photcode-skip"))) { - remove_argument (N, &argc, argv); - PHOTCODE_SKIP_LIST = strcreate(argv[N]); + // measurements with these photcodes are kept ***regardless of quality*** + // -photcode-keep J will keep all J-band measurements of all kinds + // -photcode-keep GPC1.02.g will keep all g-band measurements from chip XY02 + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-keep"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); remove_argument (N, &argc, argv); } @@ -254,8 +258,8 @@ } - PHOTCODE_SKIP_LIST = NULL; - if ((N = get_argument (argc, argv, "-photcode-skip"))) { - remove_argument (N, &argc, argv); - PHOTCODE_SKIP_LIST = strcreate(argv[N]); + PHOTCODE_KEEP_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode-keep"))) { + remove_argument (N, &argc, argv); + PHOTCODE_KEEP_LIST = strcreate(argv[N]); remove_argument (N, &argc, argv); } @@ -290,5 +294,5 @@ fprintf (stderr, " -photcode-drop : remove these photcodes from the output (REF or DEP only)\n"); - fprintf (stderr, " -photcode-skip : ignore these photcodes when assessing the validity (keep unless object is dropped)\n"); + fprintf (stderr, " -photcode-keep : ignore these photcodes when assessing the validity (keep unless object is dropped)\n"); fprintf (stderr, " -instmag (min) (max) : range of valid instrumental magnitudes (or measurements are dropped)\n"); Index: trunk/Ohana/src/photdbc/src/copy_images.c =================================================================== --- trunk/Ohana/src/photdbc/src/copy_images.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/copy_images.c (revision 36833) @@ -1,13 +1,14 @@ # include "photdbc.h" -int copy_images (char *outdir) { +int copy_images (char *outdir, SkyList *skylist) { int status; - off_t Nimage; + off_t Nimage, Nsubset; + off_t *LineNumber; char *ImageOut; unsigned int imageID; FITS_DB in; FITS_DB out; - Image *image; + Image *image, *subset; struct stat filestat; char *path; @@ -50,5 +51,8 @@ exit (2); } - dvo_image_addrows (&out, image, Nimage); + + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, FALSE); + + dvo_image_addrows (&out, subset, Nsubset); // note that imageID is unsigned int Index: trunk/Ohana/src/photdbc/src/initialize.c =================================================================== --- trunk/Ohana/src/photdbc/src/initialize.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/initialize.c (revision 36833) @@ -13,5 +13,5 @@ photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); if (SHOW_PARAMS) { @@ -54,5 +54,5 @@ photcodesDrop = ParsePhotcodeList (PHOTCODE_DROP_LIST, &NphotcodesDrop, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); if (SHOW_PARAMS) { Index: trunk/Ohana/src/photdbc/src/make_subcatalog.c =================================================================== --- trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/make_subcatalog.c (revision 36833) @@ -9,5 +9,5 @@ off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt; double mag, minMag, minSigma; - int keep, *secSkip; + int keep, *secKeep; PhotCode *photcode; @@ -16,11 +16,11 @@ // set up a list of SEC entries to ignore when evaluating a source - ALLOCATE (secSkip, int, Nsecfilt); + ALLOCATE (secKeep, int, Nsecfilt); for (i = 0; i < Nsecfilt; i++) { - secSkip[i] = FALSE; + secKeep[i] = FALSE; photcode = GetPhotcodebyNsec(i); - for (k = 0; k < NphotcodesSkip; k++) { - if (photcodesSkip[k][0].code != photcode[0].code) continue; - secSkip[i] = TRUE; + for (k = 0; k < NphotcodesKeep; k++) { + if (photcodesKeep[k][0].code != photcode[0].code) continue; + secKeep[i] = TRUE; } } @@ -67,5 +67,5 @@ keep = FALSE; for (j = 0; !keep && (j < Nsecfilt); j++) { - if (secSkip[j]) continue; + if (secKeep[j]) continue; if (catalog[0].secfilt[Nsecfilt*i+j].Ncode >= NCODE_MIN) { keep = TRUE; @@ -101,7 +101,7 @@ if (NphotcodesDrop > 0) { found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + for (k = 0; (k < NphotcodesDrop) && !found; k++) { + if (photcodesDrop[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesDrop[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; } if (found) continue; @@ -109,9 +109,9 @@ // ignore certain photcodes to assess the measurements - if (NphotcodesSkip > 0) { + if (NphotcodesKeep > 0) { found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; + for (k = 0; (k < NphotcodesKeep) && !found; k++) { + if (photcodesKeep[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(catalog[0].measure[offset].photcode)) found = TRUE; } if (found) goto keep; Index: trunk/Ohana/src/photdbc/src/photdbc.c =================================================================== --- trunk/Ohana/src/photdbc/src/photdbc.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/photdbc.c (revision 36833) @@ -9,11 +9,11 @@ initialize (argc, argv); - // load and copy the image table - copy_images (argv[1]); - // the output catalog needs to inherit the SKY_DEPTH of the input catalog sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, 0, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); skylist = SkyListByPatch (sky, -1, ®ION); + + // load and copy the image table + copy_images (argv[1], skylist); // hostID is 0 for master program Index: trunk/Ohana/src/photdbc/src/photdbc_catalogs.c =================================================================== --- trunk/Ohana/src/photdbc/src/photdbc_catalogs.c (revision 36681) +++ trunk/Ohana/src/photdbc/src/photdbc_catalogs.c (revision 36833) @@ -120,5 +120,5 @@ if (ExcludeByMaxMinMag) { snprintf (tmpline, DVO_MAX_PATH, "%s -maxminmag %f", command, MAX_MIN_MAG); strcpy (command, tmpline); } if (PHOTCODE_DROP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-drop %s", command, PHOTCODE_DROP_LIST); strcpy (command, tmpline); } - if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-skip %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, DVO_MAX_PATH, "%s -photcode-keep %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } if (CATFORMAT) { snprintf (tmpline, DVO_MAX_PATH, "%s -set-format %s", command, CATFORMAT); strcpy (command, tmpline); } if (CATMODE) { snprintf (tmpline, DVO_MAX_PATH, "%s -set-mode %s", command, CATMODE); strcpy (command, tmpline); } Index: trunk/Ohana/src/photdbc/src/select_images.c =================================================================== --- trunk/Ohana/src/photdbc/src/select_images.c (revision 36833) +++ trunk/Ohana/src/photdbc/src/select_images.c (revision 36833) @@ -0,0 +1,234 @@ +# include "photdbc.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). All images in the image catalog are tested once, so there is no check that an + image already has been included. LineNum stores the locations in the Image database of the + list of images */ + +Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap) { + + Image *image; + off_t i, j, k, nimage, NIMAGE, D_NIMAGE; + off_t *line_number; + int found; + double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + Coords tcoords; + SkyRegionCoords *skycoords; + + double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion; + + double *RmaxSky; + off_t *index; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + // if no region is selected, we are choosing the whole sky; I should sort-circuit + // this function in that case + + INITTIME; + + // the comparison is made in the catalog local projection. below we set crval1,2 + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "DEC--TAN"); + + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + ALLOCATE (RmaxSky, double, skylist[0].Nregions); + ALLOCATE (index, off_t, skylist[0].Nregions); + + RminSkyRegion = +360.0; + RmaxSkyRegion = -360.0; + DminSkyRegion = +90.0; + DmaxSkyRegion = -90.0; + + D_NIMAGE = 6000; + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + skycoords[i].Yc[4] = skycoords[i].Yc[0]; + + RmaxSky[i] = skylist[0].regions[i][0].Rmax; + index[i] = i; + + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); + dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + + RminSkyRegion = MIN(RminSkyRegion, skylist[0].regions[i][0].Rmin); + RmaxSkyRegion = MAX(RmaxSkyRegion, skylist[0].regions[i][0].Rmax); + DminSkyRegion = MIN(DminSkyRegion, skylist[0].regions[i][0].Dmin); + DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax); + } + RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion); + MARKTIME("create sky region coords: %f sec\n", dtime); + + dsortindex (RmaxSky, index, skylist[0].Nregions); + MARKTIME("sort sky coords: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "finding images\n"); + BuildChipMatch (timage, Ntimage); + MARKTIME("build chip match: %f sec\n", dtime); + + nimage = 0; + NIMAGE = D_NIMAGE; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, off_t, NIMAGE); + + // go through the complete list of images, selecting ones which overlap any region + for (i = 0; i < Ntimage; i++) { + + /* exclude images by photcode, or equiv photcode, if specified */ + if (NphotcodesDrop > 0) { + if (timage[i].photcode == 0) goto keep_image; // PHU image can be immediately passed in here + found = FALSE; // here 'found' refers to a match to a drop photcode + for (k = 0; (k < NphotcodesDrop) && !found; k++) { + if (photcodesDrop[k][0].code == timage[i].photcode) found = TRUE; + if (photcodesDrop[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE; + } + if (found) continue; + } +keep_image: + + // this adds 1.3 sec for 3M images + if (!FindMosaicForImage (timage, Ntimage, i)) { + fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); + continue; + } + + /* define image corners - note the DIS images (mosaic phu) are special */ + if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { + Xi[0] = -0.5*timage[i].NX; Yi[0] = -0.5*timage[i].NY; + Xi[1] = +0.5*timage[i].NX; Yi[1] = -0.5*timage[i].NY; + Xi[2] = +0.5*timage[i].NX; Yi[2] = +0.5*timage[i].NY; + Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY; + Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY; + } else { + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; + } + found = FALSE; + + /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */ + double RminImage = RmidSkyRegion + 180.0; + double RmaxImage = RmidSkyRegion - 180.0; + double DminImage = +90.0; + double DmaxImage = -90.0; + for (j = 0; j < 5; j++) { + XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + Ri[j] = ohana_normalize_angle_to_midpoint (Ri[j], RmidSkyRegion); + + RminImage = MIN(RminImage, Ri[j]); + RmaxImage = MAX(RmaxImage, Ri[j]); + DminImage = MIN(DminImage, Di[j]); + DmaxImage = MAX(DmaxImage, Di[j]); + } + if (RmaxImage - RminImage > 180.0) { + double tmp = RminImage; + RmaxImage = RminImage; + RminImage = tmp - 360.0; + } + + // check that this image is even in range of the searched region + if (DminImage > DmaxSkyRegion) continue; + if (DmaxImage < DminSkyRegion) continue; + + // the sky region RA is defined to be 0 - 360.0 + if (RminImage > RmaxSkyRegion) continue; + if (RmaxImage < RminSkyRegion) continue; + + // require that the full image be inside the region of interest (only for non PHU + // images) XXX : if we calibrate the mosaic, require that the full mosaic be inside + // the region and all of its chips.. + if (UseFullOverlap && strcmp(&timage[i].coords.ctype[4], "-DIS")) { + if (RmaxImage > REGION.Rmax) continue; + if (RminImage < REGION.Rmin) continue; + if (DmaxImage > REGION.Dmax) continue; + if (DminImage < REGION.Dmin) continue; + } + + image[nimage] = timage[i]; + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += D_NIMAGE; + D_NIMAGE = MAX (100000, D_NIMAGE * 1.5); + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, off_t, NIMAGE); + } + } + MARKTIME("finish image selection: %f sec\n", dtime); + + if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, off_t, MAX (nimage, 1)); + free (skycoords); + free (RmaxSky); + free (index); + + *Nimage = nimage; + *LineNumber = line_number; + return (image); +} + +void dsortindex (double *X, off_t *Y, int N) { + +# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \ + tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \ + tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \ +} +# define COMPARE(A,B)(X[A] < X[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + +off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) { + + // use bisection to find the overlapping mosaic + + off_t Nlo, Nhi, N; + + // find the last mosaic before start + Nlo = 0; Nhi = Nregions; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (Rref[N] < R) { + Nlo = MAX(N, 0); + } else { + Nhi = MIN(N, Nregions); + } + } + return (Nlo); +} Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 36681) +++ trunk/Ohana/src/relastro/Makefile (revision 36833) @@ -127,4 +127,5 @@ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o \ +$(SRC)/syncfile.$(ARCH).o \ $(SRC)/BrightCatalog.$(ARCH).o Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 36681) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 36833) @@ -83,7 +83,8 @@ double L, M; /* Focal Plane - pixels */ double X, Y; /* Chip Coords - pixels */ - double Mag; - double dMag; - double dPos; + float Mag; + float Color; + float dMag; + float dPos; int mask; int Nmeas; @@ -207,7 +208,7 @@ double MaxDensityValue; -char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST; -int NphotcodesKeep, NphotcodesSkip; -PhotCode **photcodesKeep, **photcodesSkip; +char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST, *PHOTCODE_RESET_LIST; +int NphotcodesKeep, NphotcodesSkip, NphotcodesReset; +PhotCode **photcodesKeep, **photcodesSkip, **photcodesReset; char *PHOTCODE_A_LIST, *PHOTCODE_B_LIST; @@ -217,4 +218,8 @@ SkyRegionSelection SELECTION; +char *DCR_COLOR_POS, *DCR_COLOR_NEG; +PhotCode *DCR_PHOTCODE_POS, *DCR_PHOTCODE_NEG; +int DCR_NSEC_POS, DCR_NSEC_NEG; + int ImagSelect; double ImagMin, ImagMax; @@ -233,4 +238,7 @@ int CLIP_THRESH; int USE_BASIC_CHECK; + +int ExcludeBogus; +double ExcludeBogusRadius; FitMode FIT_MODE; @@ -307,5 +315,5 @@ int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); int liststats_pos PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB)); -Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile)); int load_images PROTO((FITS_DB *db, SkyList *skylist, int UseFullOverlap)); Image *select_images PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage, int UseFullOverlap)); @@ -417,4 +425,6 @@ int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt); int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt); +double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt); +double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt); int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt); int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt); @@ -521,2 +531,9 @@ int ImageTableSave (char *filename, Image *images, off_t Nimages); int select_mosaics_hostregion (RegionHostTable *regionHosts, Image *image, off_t Nimage); + +float getColor (off_t meas, int cat); + +int strextend (char *input, char *format,...); + +int areImagesLoaded (); +int areImagesMatched (); Index: trunk/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 36681) +++ trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 36833) @@ -167,5 +167,5 @@ GET_COLUMN(measureOffset, "OFF_MEASURE", int); GET_COLUMN(missingOffset, "OFF_MISSING", int); - GET_COLUMN(extendOffset, "OFF_EXTEND", int); + GET_COLUMN(refColor, "REF_COLOR", float); GET_COLUMN(flags, "FLAGS", int); GET_COLUMN(photFlagsUpper, "PHOTFLAGS_U", int); @@ -202,5 +202,5 @@ average[i].measureOffset = measureOffset[i] ; average[i].missingOffset = missingOffset[i] ; - average[i].extendOffset = extendOffset[i] ; + average[i].refColor = refColor[i] ; average[i].flags = flags[i] ; average[i].photFlagsUpper = photFlagsUpper[i] ; @@ -234,5 +234,5 @@ free (measureOffset); free (missingOffset); - free (extendOffset); + free (refColor); free (flags); free (photFlagsUpper); @@ -476,5 +476,5 @@ gfits_define_bintable_column (&theader, "J", "OFF_MEASURE", "offset to first psf measurement", "", 1.0, 1.0*0x8000); gfits_define_bintable_column (&theader, "J", "OFF_MISSING", "offset to first missing obs", "", 1.0, 1.0*0x8000); - gfits_define_bintable_column (&theader, "J", "OFF_EXTEND", "offset to first extended measurement", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "E", "REF_COLOR", "reference color", "", 1.0, 0.0); gfits_define_bintable_column (&theader, "J", "FLAGS", "average object flags (star; ghost; etc)", "", 1.0, 1.0*0x8000); gfits_define_bintable_column (&theader, "J", "PHOTFLAGS_U", "upper bit of 2 bit summary of per-measure photflags", "", 1.0, 1.0*0x8000); @@ -510,5 +510,5 @@ int *measureOffset ; ALLOCATE (measureOffset , int , catalog->Naverage); int *missingOffset ; ALLOCATE (missingOffset , int , catalog->Naverage); - int *extendOffset ; ALLOCATE (extendOffset , int , catalog->Naverage); + float *refColor ; ALLOCATE (refColor , float , catalog->Naverage); int *flags ; ALLOCATE (flags , int , catalog->Naverage); int *photFlagsUpper; ALLOCATE (photFlagsUpper, int , catalog->Naverage); @@ -543,5 +543,5 @@ measureOffset[i] = average[i].measureOffset ; missingOffset[i] = average[i].missingOffset ; - extendOffset[i] = average[i].extendOffset ; + refColor[i] = average[i].refColor ; flags[i] = average[i].flags ; photFlagsUpper[i] = average[i].photFlagsUpper ; @@ -575,5 +575,5 @@ gfits_set_bintable_column (&theader, &ftable, "OFF_MEASURE", measureOffset, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "OFF_MISSING", missingOffset, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "OFF_EXTEND", extendOffset, catalog->Naverage); + gfits_set_bintable_column (&theader, &ftable, "REF_COLOR", refColor, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); gfits_set_bintable_column (&theader, &ftable, "PHOTFLAGS_U", photFlagsUpper, catalog->Naverage); @@ -605,5 +605,5 @@ free (measureOffset); free (missingOffset); - free (extendOffset); + free (refColor); free (flags); free (photFlagsUpper); Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 36681) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 36833) @@ -4,5 +4,5 @@ int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) { - int i, NstatFull, Nstat, Niter; + int i, NstatFull, Nstat, Niter, Nkeep; float dLsig, dMsig, dRsig; float dLsigFull, dMsigFull, dRsigFull; @@ -24,4 +24,5 @@ // measure the scatter distribution (use only the bright end detections) + Nkeep = 0; for (i = 0; i < Nmatch; i++) { if (raw[i].mask) continue; @@ -42,5 +43,9 @@ continue; } - } + Nkeep ++; + } + + // I'm rejecting some points from the fit above; I count the remainders and actually + // use that count to set order_use below // figures to assess the fitting process: @@ -50,16 +55,15 @@ relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch); - // fit the requested order polynomial - if (CHIPORDER > 0) { - int Ndof_min = 3; - int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3); - int order_use = MIN (CHIPORDER, order_max); - if (order_use < 1) { - if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch); - image[0].flags |= ID_IMAGE_ASTROM_FEW; - return FALSE; - } - image[0].coords.Npolyterms = order_use; - } + // set the maximum order for the polynomial (based on number of stars kept above) + int order_use = 0; + if (Nkeep > 5) order_use = 1; // 4 stars per polynomial term (per dimension) + if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension) + if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension) + if (order_use < 1) { + if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep); + image[0].flags |= ID_IMAGE_ASTROM_FEW; + return FALSE; + } + image[0].coords.Npolyterms = order_use; if (fit) fit_free (fit); @@ -71,30 +75,4 @@ fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos); } - -# if (0) - // check if the fit has enough data points for the polynomial order - skip = FALSE; - switch (image[0].coords.Npolyterms) { - case 0: - case 1: - skip = (fit[0].Npts < 8); - break; - case 2: - skip = (fit[0].Npts < 11); - break; - case 3: - skip = (fit[0].Npts < 15); - break; - default: - fprintf (stderr, "invalid chip order %d\n", image[0].coords.Npolyterms); - skip = TRUE; - } - if (skip) { - if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms); - fit_free (fit); - image[0].flags |= ID_IMAGE_ASTROM_FEW; - return FALSE; - } -# endif // measure the fit, update the coords & object coordinates @@ -134,4 +112,20 @@ } + int Ncolor = 0; + float *colorList = NULL; + ALLOCATE (colorList, float, Nmatch); + + for (i = 0; i < Nmatch; i++) { + if (!raw[i].mask) continue; + if (!isnan(ref[i].Color)) { + colorList[Ncolor] = ref[i].Color; + Ncolor ++; + } + } + fsort (colorList, Ncolor); + float colorMedian = (Ncolor > 0) ? colorList[(int)(0.5*Ncolor)] : NAN; + image[0].refColor = colorMedian; + free (colorList); + GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM); GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM); @@ -151,4 +145,6 @@ image[0].dYpixSys = dMsig; image[0].nFitAstrom = fit[0].Npts; + + // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom); if (fit) fit_free (fit); Index: trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- trunk/Ohana/src/relastro/src/GetAstromError.c (revision 36681) +++ trunk/Ohana/src/relastro/src/GetAstromError.c (revision 36833) @@ -30,4 +30,5 @@ code = GetPhotcodebyCode (measure[0].photcode); + if (!code) return NAN; // do not raise an exception, just send back the result @@ -72,4 +73,5 @@ code = GetPhotcodebyCode (measure[0].photcode); + if (!code) return NAN; // do not raise an exception, just send back the result Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 36681) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 36833) @@ -39,4 +39,16 @@ // N_onImage was 'Nlist' // N_ONIMAGE was 'NLIST' + +int areImagesLoaded () { + + if (image) return TRUE; + return FALSE; +} + +int areImagesMatched () { + + if (MeasureToImage) return TRUE; + return FALSE; +} Image *getimages (off_t *N, off_t **line_number) { @@ -333,4 +345,13 @@ */ +float getColor (off_t meas, int cat) { + + off_t i; + + i = MeasureToImage[cat][meas]; + if (i == -1) return (NAN); + return (image[i].refColor); +} + void plot_images () { @@ -468,10 +489,15 @@ } + X = measure[0].Xccd; + Y = measure[0].Yccd; if (USE_FIXED_PIXCOORDS) { - X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd; - Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd; - } else { - X = measure[0].Xccd; - Y = measure[0].Yccd; + if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) { + float dX = measure[0].Xfix - measure[0].Xccd; + float dY = measure[0].Yfix - measure[0].Yccd; + if (hypot(dX,dY) < 2.0) { + X = measure[0].Xfix; + Y = measure[0].Yfix; + } + } } n = measure[0].averef; @@ -560,5 +586,5 @@ int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori; - if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); + if (VERBOSE2 && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); saveOffsets (dPos, nPos, im); @@ -609,10 +635,15 @@ Measure *measure = &catalog[c].measure[m]; + X = measure[0].Xccd; + Y = measure[0].Yccd; if (USE_FIXED_PIXCOORDS) { - X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd; - Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd; - } else { - X = measure[0].Xccd; - Y = measure[0].Yccd; + if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) { + float dX = measure[0].Xfix - measure[0].Xccd; + float dY = measure[0].Yfix - measure[0].Yccd; + if (hypot(dX,dY) < 2.0) { + X = measure[0].Xfix; + Y = measure[0].Yfix; + } + } } n = measure[0].averef; @@ -679,12 +710,16 @@ /* apply the current image transformation or use the current value of R+dR, D+dD? */ + raw[i].X = measure[0].Xccd; + raw[i].Y = measure[0].Yccd; if (USE_FIXED_PIXCOORDS) { - raw[i].X = isfinite(measure[0].Xfix) ? measure[0].Xfix : measure[0].Xccd; - raw[i].Y = isfinite(measure[0].Yfix) ? measure[0].Yfix : measure[0].Yccd; - } else { - raw[i].X = measure[0].Xccd; - raw[i].Y = measure[0].Yccd; - } - + if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) { + float dX = measure[0].Xfix - measure[0].Xccd; + float dY = measure[0].Yfix - measure[0].Yccd; + if (hypot(dX,dY) < 2.0) { + raw[i].X = measure[0].Xfix; + raw[i].Y = measure[0].Yfix; + } + } + } raw[i].Mag = measure[0].M; raw[i].dMag = measure[0].dM; @@ -738,4 +773,6 @@ StarData *ref; + int Nsecfilt = GetPhotcodeNsecfilt(); + mosaic = NULL; moscoords = NULL; @@ -765,4 +802,8 @@ ref[i].dMag = measure[0].dM; ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); + + if ((DCR_NSEC_POS >= 0) && (DCR_NSEC_NEG >= -1)) { + ref[i].Color = catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_POS].M - catalog[c].secfilt[n*Nsecfilt + DCR_NSEC_NEG].M; + } ref[i].mask = FALSE; @@ -859,5 +900,5 @@ N = 0; for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - //reset flag on each invocation + // reset flag on each invocation measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; @@ -898,4 +939,7 @@ double *R, *D, *dR, *dD, *d2; StatType statsR, statsD; + + // XXX we are not going to use this for now + return; Nsecfilt = GetPhotcodeNsecfilt(); @@ -1063,4 +1107,5 @@ } else { code = GetPhotcodebyCode (measure[0].photcode); + if (!code) return FALSE; mask = code[0].astromBadMask; } @@ -1129,4 +1174,5 @@ } else { code = GetPhotcodebyCode (measure[0].photcode); + if (!code) return FALSE; mask = code[0].astromBadMask; } Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 36681) +++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 36833) @@ -126,41 +126,44 @@ table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); - char tmpline[1024]; - if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } - if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } - if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } - - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } - if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } - - if (RESET_BAD_IMAGES) { snprintf (tmpline, 1024, "%s -reset-bad-images", command); strcpy (command, tmpline); } - - if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } - if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_ONLY) strextend (command, "-pm"); + if (FIT_MODE == FIT_PAR_ONLY) strextend (command, "-par"); + if (FIT_MODE == FIT_PM_AND_PAR) strextend (command, "-pmpar"); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (UPDATE) strextend (command, "-update"); + + if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images"); + + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); - if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command); strcpy (command, tmpline); } - if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); strcpy (command, tmpline); } + if (USE_BASIC_CHECK) strextend (command, "-basic-image-search"); + if (FlagOutlier) strextend (command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); - if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command); strcpy (command, tmpline); } - - if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } - if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } - if (PhotFlagSelect) { snprintf (tmpline, 1024, "%s +photflags", command); strcpy (command, tmpline); } - if (PhotFlagBad) { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad); strcpy (command, tmpline); } - if (PhotFlagPoor) { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor); strcpy (command, tmpline); } + if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); + + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (MinBadQF > 0.0) { snprintf (tmpline, 1024, "%s -min-bad-psfqf %f", command, MinBadQF); strcpy (command, tmpline); } - if (MaxMeanOffset != 10.0) { snprintf (tmpline, 1024, "%s -max-mean-offset %f", command, MaxMeanOffset); strcpy (command, tmpline); } + if (DCR_COLOR_POS && DCR_COLOR_NEG) { + strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); + } + + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); if (TimeSelect) { char *tstart = ohana_sec_to_date (TSTART); char *tstop = ohana_sec_to_date (TSTOP); - snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + strextend (command, "-time %s %s", tstart, tstop); free (tstart); free (tstop); - strcpy (command, tmpline); } fprintf (stderr, "command: %s\n", command); Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36681) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 36833) @@ -10,4 +10,5 @@ static double *T; static double *dT; +static double *C; void initObjectData (Catalog *catalog, int Ncatalog) { @@ -36,4 +37,6 @@ ALLOCATE (pX, double, MAX (1, Nmax)); ALLOCATE (pY, double, MAX (1, Nmax)); + + ALLOCATE (C, double, MAX (1, Nmax)); } @@ -54,4 +57,6 @@ free (pX); free (pY); + + free (C); } @@ -61,5 +66,5 @@ off_t j, k, m; - int i, N, Nsecfilt, mode, result, status, XVERB; + int i, N, Nc, Nsecfilt, mode, result, status, XVERB; StatType statsR, statsD; Coords coords; @@ -75,4 +80,6 @@ memset (&fitPAR, 0, sizeof(fitPAR)); initObjectData (catalog, Ncatalog); + + int setRefColor = areImagesMatched(); /* project coordinates to a plane centered on the object with units of arcsec */ @@ -112,8 +119,16 @@ fitPAR.chisq = NAN; + // if we fail to fit the astrometry for some reason, we need to set/reset these + catalog[i].average[j].flags |= ID_STAR_NO_ASTROM; + catalog[i].average[j].ChiSqAve = NAN; + catalog[i].average[j].ChiSqPM = NAN; + catalog[i].average[j].ChiSqPar = NAN; + catalog[i].average[j].Npos = 0; + if (catalog[i].average[j].Nmeasure == 0) { - continue; - } - + continue; + } + + Nc = 0; N = 0; m = catalog[i].average[j].measureOffset; @@ -138,22 +153,38 @@ } + // reset the bit to note that a detection was used (or not) + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } + // does the measurement pass the supplied filtering constraints? // MeasFilterTestTiny does not test psfQF + // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE); if (!keepMeasure) { - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } continue; } + // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius + if (ExcludeBogus) { + double Ri = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + double Di = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + coords.crval1 = catalog[i].average[j].R; + coords.crval2 = catalog[i].average[j].D; + double Xi, Yi; + RD_to_XY (&Xi, &Yi, Ri, Di, &coords); + double radius = hypot(Xi, Yi); + if (radius > ExcludeBogusRadius) { + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; } + continue; + } + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; } + } + // outlier rejection - if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } + if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { continue; } - - measure[k].dbFlags |= ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); @@ -170,14 +201,6 @@ // allow a given photcode or measurement to be // ignored if the error is NAN (for photcode, set astromErrSys to NaN) - if (isnan(dX[N])) { - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } - continue; - } - if (isnan(dY[N])) { - measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; - if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } - continue; - } + if (isnan(dX[N])) continue; + if (isnan(dY[N])) continue; // add systematic error in quadrature, if desired @@ -201,12 +224,21 @@ dD[N] = dY[N] / 3600.0; + if (setRefColor) { + float color = getColor (m+k, i); + if (!isnan(color)) { + C[Nc] = color; + Nc++; + } + } + + measure[k].dbFlags |= ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } + N++; - } + } // loop over measurements : catalog[i].average[j].Nmeasure // if we have too few good detections for the desired fit, or too limited a // baseline, use a fit with fewer parameters. XXX if we have too few measurements // for even the average position, consider including the lower-quality detections? - - catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM; // find Tmin & Tmax from the list of accepted measurements @@ -231,10 +263,4 @@ // too few measurements for average position (require 2 values) if (N < SRC_MEAS_TOOFEW) { - // XXX need to define PHOTOM and ASTROM object flags - // XXX reset the average value fields? - catalog[i].average[j].flags |= ID_STAR_NO_ASTROM; - catalog[i].average[j].ChiSqAve = NAN; - catalog[i].average[j].ChiSqPM = NAN; - catalog[i].average[j].ChiSqPar = NAN; if (N < 2) continue; } @@ -268,9 +294,13 @@ catalog[i].average[j].flags |= ID_STAR_FIT_PM; Npm ++; + + // XXX a hard-wired hack... + if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) { + mode = FIT_AVERAGE; + catalog[i].average[j].flags |= ID_STAR_BAD_PM; + } } if (mode == FIT_PM_AND_PAR) { - // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__); - float pXmin = +2.0; float pXmax = -2.0; @@ -296,4 +326,10 @@ catalog[i].average[j].flags |= ID_STAR_FIT_PAR; Npar ++; + + // XXX a hard-wired hack... + if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { + mode = FIT_AVERAGE; + catalog[i].average[j].flags |= ID_STAR_BAD_PM; + } } else { // need to set mode = FIT_PM_ONLY if we do not fit for parallax @@ -321,4 +357,10 @@ catalog[i].average[j].flags |= ID_STAR_FIT_AVE; Nave ++; + } + + if (setRefColor) { + dsort (C, Nc); + float colorMedian = (Nc > 0) ? C[(int)(0.5*Nc)] : NAN; + catalog[i].average[j].refColor = colorMedian; } @@ -411,4 +453,8 @@ catalog[i].average[j].Trange = (Trange * 86400 * 365.25); catalog[i].average[j].Npos = fit.Nfit; + + // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) + catalog[i].average[j].flags &= ~ID_STAR_NO_ASTROM; + if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n", catalog[i].average[j].R, Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 36681) +++ trunk/Ohana/src/relastro/src/args.c (revision 36833) @@ -179,8 +179,18 @@ FlagOutlier = FALSE; if ((N = get_argument (argc, argv, "-clip"))) { + fprintf (stderr, "-clip is currently disabled\n"); remove_argument (N, &argc, argv); CLIP_THRESH = atof (argv[N]); remove_argument (N, &argc, argv); FlagOutlier = TRUE; + } + + ExcludeBogus = FALSE; + ExcludeBogusRadius = 0.0; + if ((N = get_argument (argc, argv, "-exclude-bogus"))) { + remove_argument (N, &argc, argv); + ExcludeBogusRadius = atof (argv[N]); + remove_argument (N, &argc, argv); + ExcludeBogus = TRUE; } @@ -284,4 +294,21 @@ } + DCR_COLOR_POS = NULL; + DCR_COLOR_NEG = NULL; + if ((N = get_argument (argc, argv, "-dcr-color"))) { + remove_argument (N, &argc, argv); + DCR_COLOR_POS = strcreate(argv[N]); + remove_argument (N, &argc, argv); + DCR_COLOR_NEG = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_RESET_LIST = NULL; + if ((N = get_argument (argc, argv, "-reset-to-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_RESET_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + PHOTCODE_KEEP_LIST = NULL; if ((N = get_argument (argc, argv, "+photcode"))) { @@ -291,19 +318,21 @@ } - PHOTCODE_SKIP_LIST = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9"); + char *SuperCOSMOS_SKIP = strcreate("SCOS.103a.E,SCOS.4414.OG590,SCOS.4415.OG590,SCOS.IIIaF.OG590,SCOS.IIIaF.RG610,SCOS.IIIaF.RG630,SCOS.IIIaJ.GG385,SCOS.IIIaJ.GG395,SCOS.IVN.RG715,SCOS.IVN.RG9"); + + PHOTCODE_SKIP_LIST = NULL; if ((N = get_argument (argc, argv, "-photcode"))) { remove_argument (N, &argc, argv); - char *tmp1 = strcreate(argv[N]); - - int Ntotal = strlen(tmp1) + strlen(PHOTCODE_SKIP_LIST) + 5; - - char *tmp2 = NULL; - ALLOCATE (tmp2, char, Ntotal); - snprintf (tmp2, Ntotal, "%s,%s", PHOTCODE_SKIP_LIST, tmp1); - - free (tmp1); - free (PHOTCODE_SKIP_LIST); - - PHOTCODE_SKIP_LIST = tmp2; + char *RawSkip = strcreate(argv[N]); + + char *GotSkip = strstr (RawSkip, SuperCOSMOS_SKIP); + if (!GotSkip) { + int Ntotal = strlen(RawSkip) + strlen(SuperCOSMOS_SKIP) + 5; + ALLOCATE (PHOTCODE_SKIP_LIST, char, Ntotal); + snprintf (PHOTCODE_SKIP_LIST, Ntotal, "%s,%s", SuperCOSMOS_SKIP, RawSkip); + free (RawSkip); + free (SuperCOSMOS_SKIP); + } else { + PHOTCODE_SKIP_LIST = RawSkip; + } remove_argument (N, &argc, argv); } @@ -582,8 +611,18 @@ FlagOutlier = FALSE; if ((N = get_argument (argc, argv, "-clip"))) { + fprintf (stderr, "-clip is currently disabled\n"); remove_argument (N, &argc, argv); CLIP_THRESH = atof (argv[N]); remove_argument (N, &argc, argv); FlagOutlier = TRUE; + } + + ExcludeBogus = FALSE; + ExcludeBogusRadius = 0.0; + if ((N = get_argument (argc, argv, "-exclude-bogus"))) { + remove_argument (N, &argc, argv); + ExcludeBogusRadius = atof (argv[N]); + remove_argument (N, &argc, argv); + ExcludeBogus = TRUE; } @@ -655,4 +694,21 @@ remove_argument (N, &argc, argv); TimeSelect = TRUE; + } + + DCR_COLOR_POS = NULL; + DCR_COLOR_NEG = NULL; + if ((N = get_argument (argc, argv, "-dcr-color"))) { + remove_argument (N, &argc, argv); + DCR_COLOR_POS = strcreate(argv[N]); + remove_argument (N, &argc, argv); + DCR_COLOR_NEG = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + PHOTCODE_RESET_LIST = NULL; + if ((N = get_argument (argc, argv, "-reset-to-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_RESET_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); } @@ -861,2 +917,15 @@ exit (2); } + +int strextend (char *input, char *format,...) { + + char tmpextra[1024], tmpline[1024]; + va_list argp; + + va_start (argp, format); + vsnprintf (tmpextra, 1024, format, argp); + snprintf (tmpline, 1024, "%s %s", input, tmpextra); + strcpy (input, tmpline); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 36681) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 36833) @@ -1,8 +1,21 @@ # include "relastro.h" +static int Nkeep1 = 0; +static int Nkeep2 = 0; static int Nskip1 = 0; static int Nskip2 = 0; -static unsigned int Tref = 1323353985; -static short Cref = 10355; + +FILE *fbogus = NULL; +static int NskipBogus = 0; + +// test image: 2013/06/15,13:25:51, GPC1.r.XY50 +static int CHECK_TEST_IMAGE = FALSE; +// static unsigned int Tref = 1378812312; +// static short Cref = 10001; + +static unsigned int Tref = 1379570672; +static short Cref = 10341; + +int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog); int bcatalog (Catalog *subcatalog, Catalog *catalog) { @@ -11,4 +24,16 @@ off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; int Nsecfilt; + Coords coords; + + /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */ + coords.crval1 = 0; + coords.crval2 = 0; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "DEC--SIN"); // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? @@ -34,4 +59,13 @@ ID_STAR_NO_ASTROM ; + if (VERBOSE2 && ExcludeBogus && (fbogus == NULL)) { + char name[1024]; + snprintf (name, 1024, "%s/bogus.%02d.dat", CATDIR, HOST_ID); + fbogus = fopen (name, "w"); + if (!fbogus) { + fprintf (stderr, "trouble opening bogus detection dump : %s\n", name); + } + } + /* exclude stars not in range or with too few measurements */ for (i = 0; i < catalog[0].Naverage; i++) { @@ -60,12 +94,16 @@ if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; - if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { Nskip1 ++; } continue; } - - // filter out outliers - if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) { + if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nkeep1 ++; + } + + // filter out outliers - these are detections inconsistent with the offset distribution + // XXX disable this for now + if (FALSE && FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { @@ -76,4 +114,25 @@ catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP; + // exclude bogus + if (ExcludeBogus) { + double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); + double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); + coords.crval1 = catalog[0].average[i].R; + coords.crval2 = catalog[0].average[i].D; + double Xi, Yi; + RD_to_XY (&Xi, &Yi, Ri, Di, &coords); + double radius = hypot(Xi, Yi); + if (radius > ExcludeBogusRadius) { + NskipBogus ++; + if (VERBOSE2) { + FILE *foutput = fbogus ? fbogus : stderr; + char *date = ohana_sec_to_date(catalog[0].measure[offset].t); + fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date); + free (date); + } + continue; + } + } + // re-assess on each run of relastro if a measurement should be used @@ -87,4 +146,8 @@ // but before the final average properties are calculated, these measurements may be // allowed. + + if (CHECK_TEST_IMAGE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nkeep2 ++; + } CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); @@ -125,5 +188,5 @@ // limit the total number of stars in the catalog if (MaxDensityUse) { - LimitDensityCatalog_ByNmeasure (subcatalog, catalog); + LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog); } else { if (VERBOSE2) { @@ -134,9 +197,18 @@ } } + if (CHECK_TEST_IMAGE && (Nkeep1 + Nkeep2 + Nskip1 + Nskip2 > 0)) { + fprintf (stderr, "kept %d %d, skipped %d %d\n", Nkeep1, Nkeep2, Nskip1, Nskip2); + } + return (TRUE); } void bcatalog_show_skips () { - fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2); + if (ExcludeBogus) { + if (fbogus) fclose (fbogus); + fprintf (stderr, "NskipBogus: %d\n", NskipBogus); + // fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2); + // fprintf (stderr, "Nkeep: %d, %d\n", Nkeep1, Nkeep2); + } } @@ -155,4 +227,270 @@ # undef COMPARE +} + +/* generate a grid in a locally-projected coordinate system, try to select average entries + from each grid cell in decending Nmeasure order. + */ +int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) { + + off_t i, j; + int ix, iy; + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&oldcatalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&oldcatalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // generate a grid in locally projected space + double Rc = 0.5*(Rmin + Rmax); + double Dc = 0.5*(Dmin + Dmax); + + /* project coordinates to a plane centered on the object with units of arcsec */ + Coords coords; + coords.crval1 = Rc; + coords.crval2 = Dc; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "DEC--SIN"); + + // convert all average R,D values to X,Y: + double *X, *Y; + ALLOCATE (X, double, Naverage); + ALLOCATE (Y, double, Naverage); + float Xmin = +10000.0, Ymin = +10000.0; + float Xmax = -10000.0, Ymax = -10000.0; + for (i = 0; i < Naverage; i++) { + X[i] = NAN; + Y[i] = NAN; + // skip any stars which are outside of nominal catalog range + if (subcatalog[0].average[i].R < Rmin) continue; + if (subcatalog[0].average[i].R > Rmax) continue; + if (subcatalog[0].average[i].D < Dmin) continue; + if (subcatalog[0].average[i].D > Dmax) continue; + RD_to_XY (&X[i], &Y[i], subcatalog[0].average[i].R, subcatalog[0].average[i].D, &coords); + Xmin = MIN (Xmin, X[i]); + Xmax = MAX (Xmax, X[i]); + Ymin = MIN (Ymin, Y[i]); + Ymax = MAX (Ymax, Y[i]); + } + + // how many grid cells? what is the grid spacing? + float dX = Xmax - Xmin; + float dY = Ymax - Ymin; + + // *** XXX for the moment, I'm using a hard-wired cell size (200 arcsec ~ 3.3 arcmin) + int NX = (int)(dX / 200) + 1; + int NY = (int)(dY / 200) + 1; + // fprintf (stderr, "Density Grid: %d x %d\n", NX, NY); + // XXX check that NX,NY are sensible (5 degrees / 200 arcsec seems like the absolute max) + if (NX > 1000) { + fprintf (stderr, "serious problem with %s: NX = %d\n", subcatalog[0].filename, NX); + exit (3); + } + if (NY > 1000) { + fprintf (stderr, "serious problem with %s: NY = %d\n", subcatalog[0].filename, NY); + exit (3); + } + + // kind of ugly : generate a grid of index, Nmeasure arrays + // to be filled below (I also need NN and Nn to track the number of + // entries in each). + int **NN_grid; + int **Nn_grid; + int ***Nm_grid; + off_t ***idxgrid; + ALLOCATE (NN_grid, int *, NX); + ALLOCATE (Nn_grid, int *, NX); + ALLOCATE (Nm_grid, int **, NX); + ALLOCATE (idxgrid, off_t **, NX); + + for (ix = 0; ix < NX; ix++) { + ALLOCATE (NN_grid[ix], int, NY); + ALLOCATE (Nn_grid[ix], int, NY); + ALLOCATE (Nm_grid[ix], int *, NY); + ALLOCATE (idxgrid[ix], off_t *, NY); + for (iy = 0; iy < NY; iy++) { + Nn_grid[ix][iy] = 0; + NN_grid[ix][iy] = 100; + ALLOCATE (Nm_grid[ix][iy], int, NN_grid[ix][iy]); + ALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]); + } + } + + // assign all of the average entries to a grid cell + for (i = 0; i < Naverage; i++) { + if (isnan(X[i])) continue; + if (isnan(Y[i])) continue; + ix = (X[i] - Xmin) / 200.0; + iy = (Y[i] - Ymin) / 200.0; + int Nn = Nn_grid[ix][iy]; + Nm_grid[ix][iy][Nn] = subcatalog[0].average[i].Nmeasure; + + // if we are resetting to a given photcode, we need to have that photcode... + if (NphotcodesReset) { + int k; + int foundReset = FALSE; + int m = subcatalog[0].average[i].measureOffset; + MeasureTiny *measure = &subcatalog[0].measureT[m]; + for (j = 0; (j < subcatalog[0].average[i].Nmeasure) && !foundReset; j++) { + if (CHECK_TEST_IMAGE && (abs(measure[j].t - Tref) < 10) && (measure[j].photcode == Cref)) { + fprintf (stderr, "."); + } + for (k = 0; (k < NphotcodesReset) && !foundReset; k++) { + if (photcodesReset[k][0].code == measure[j].photcode) foundReset = TRUE; + } + } + if (!foundReset) { + Nm_grid[ix][iy][Nn] = 0; + } + } + + idxgrid[ix][iy][Nn] = i; + Nn_grid[ix][iy] ++; + if (Nn_grid[ix][iy] >= NN_grid[ix][iy]) { + NN_grid[ix][iy] += 100; + REALLOCATE (Nm_grid[ix][iy], int, NN_grid[ix][iy]); + REALLOCATE (idxgrid[ix][iy], off_t, NN_grid[ix][iy]); + } + } + + // sort all of the grid cells + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + sort_by_Nmeasure (Nm_grid[ix][iy], idxgrid[ix][iy], Nn_grid[ix][iy]); + NN_grid[ix][iy] = 0; // I'm going to use this array to track which element I've already selected + } + } + + // cycle over the grid until we ready Nmax + off_t *keepidx = NULL; + ALLOCATE (keepidx, off_t, Naverage); + memset (keepidx, 0, Naverage*sizeof(off_t)); + int Nkeep = 0; + + for (i = 0; (i < 20) && (Nkeep < Nmax); i++) { + for (ix = 0; (ix < NX) && (Nkeep < Nmax); ix++) { + for (iy = 0; (iy < NY) && (Nkeep < Nmax); iy++) { + if (NN_grid[ix][iy] >= Nn_grid[ix][iy]) continue; // all used up! + int Nn = NN_grid[ix][iy]; + keepidx[Nkeep] = idxgrid[ix][iy][Nn]; + Nkeep ++; + NN_grid[ix][iy] ++; + } + } + } + + // count the number of measurements this selection will yield + off_t ave, NMEASURE = 0; + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + + // test catID : 37262 37261 37257 37258 + int dumpit = FALSE; + dumpit |= (oldcatalog[0].catID == 37007); + // dumpit |= (oldcatalog[0].catID == 37261); + // dumpit |= (oldcatalog[0].catID == 37257); + // dumpit |= (oldcatalog[0].catID == 37258); + if (dumpit) { + char name[64]; + snprintf (name, 64, "cat.%05d.dump.dat", oldcatalog[0].catID); + FILE *fdump = fopen (name, "w"); + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure); + } + fclose (fdump); + } + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nkeep); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nkeep * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nkeep selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nkeep; i++) { + ave = keepidx[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { + char *basename = filebasename (oldcatalog[0].filename); + fprintf (stderr, "limited to %d ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", + Nkeep, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, oldcatalog[0].Nmeasure, basename); + free (basename); + } + + free (X); + free (Y); + + for (ix = 0; ix < NX; ix++) { + for (iy = 0; iy < NY; iy++) { + free (Nm_grid[ix][iy]); + free (idxgrid[ix][iy]); + } + free (NN_grid[ix]); + free (Nn_grid[ix]); + free (Nm_grid[ix]); + free (idxgrid[ix]); + } + free (NN_grid); + free (Nn_grid); + free (Nm_grid); + free (idxgrid); + + free (keepidx); + + free (subcatalog[0].average); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nkeep; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt; + + return (TRUE); } @@ -205,4 +543,15 @@ NMEASURE += subcatalog[0].average[ave].Nmeasure; } + +# if (0) + if (oldcatalog[0].catID == 59962) { + FILE *fdump = fopen ("cat.dump.dat", "w"); + for (i = 0; i < Nmax; i++) { + ave = index[i]; + fprintf (fdump, "%10.6f %10.6f %d\n", subcatalog[0].average[ave].R, subcatalog[0].average[ave].D, subcatalog[0].average[ave].Nmeasure); + } + fclose (fdump); + } +# endif // allocate the output data Index: trunk/Ohana/src/relastro/src/dvo_astrom_ops.c =================================================================== --- trunk/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 36681) +++ trunk/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 36833) @@ -91,4 +91,44 @@ } +double getMeanR_Big (Measure *measure, Average *average, SecFilt *secfilt) { + + double ra; + + /* the measure carries the instantaneous mean position at the epoch t */ + ra = average[0].R - measure[0].dR / 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (ra); +} + +double getMeanD_Big (Measure *measure, Average *average, SecFilt *secfilt) { + + double dec; + + /* the measure carries the instantaneous mean position at the epoch t */ + dec = average[0].D - measure[0].dD / 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (dec); +} + int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) { @@ -144,2 +184,3 @@ return (TRUE); } + Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 36681) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 36833) @@ -8,13 +8,36 @@ if (RELASTRO_OP == OP_MERGE_SOURCE) return; - fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); - fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); - fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); - fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); + if (DCR_COLOR_POS) fprintf (stderr, "DCR_COLOR_POS: %s\n", DCR_COLOR_POS); + if (DCR_COLOR_NEG) fprintf (stderr, "DCR_COLOR_NEG: %s\n", DCR_COLOR_NEG); - photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); - photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); - photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); + if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST); + if (PHOTCODE_A_LIST) fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); + if (PHOTCODE_B_LIST) fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); + + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesReset = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset, FALSE); + photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); + photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); + + DCR_NSEC_POS = DCR_NSEC_NEG = -1; + if (DCR_COLOR_POS) { + DCR_PHOTCODE_POS = GetPhotcodebyName (DCR_COLOR_POS); + if (!DCR_PHOTCODE_POS) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_POS); + exit (1); + } + DCR_NSEC_POS = GetPhotcodeNsec (DCR_PHOTCODE_POS[0].code); + } + if (DCR_COLOR_NEG) { + DCR_PHOTCODE_NEG = GetPhotcodebyName (DCR_COLOR_NEG); + if (!DCR_PHOTCODE_NEG) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", DCR_COLOR_NEG); + exit (1); + } + DCR_NSEC_NEG = GetPhotcodeNsec (DCR_PHOTCODE_NEG[0].code); + } initstats (STATMODE); @@ -42,13 +65,15 @@ args_client (argc, argv); - fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); - fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); - fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); - fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); + if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + if (PHOTCODE_RESET_LIST) fprintf (stderr, "PHOTCODE_RESET_LIST: %s\n", PHOTCODE_RESET_LIST); + if (PHOTCODE_A_LIST) fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); + if (PHOTCODE_B_LIST) fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); - photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); - photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); - photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); - photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesReset = ParsePhotcodeList (PHOTCODE_RESET_LIST, &NphotcodesReset, FALSE); + photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); + photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); initstats (STATMODE); Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36681) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 36833) @@ -1,17 +1,4 @@ # include "relastro.h" # define DEBUG 0 - -int strextend (char *input, char *format,...) { - - char tmpextra[1024], tmpline[1024]; - va_list argp; - - va_start (argp, format); - vsnprintf (tmpextra, 1024, format, argp); - snprintf (tmpline, 1024, "%s %s", input, tmpextra); - strcpy (input, tmpline); - - return TRUE; -} int launch_region_hosts (RegionHostTable *regionHosts) { @@ -88,7 +75,14 @@ if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PHOTCODE_RESET_LIST) strextend (command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST); if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); + + if (DCR_COLOR_POS && DCR_COLOR_NEG) { + strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); + } + if (PhotFlagSelect) strextend (command, "+photflags"); if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); @@ -97,4 +91,11 @@ if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); + + strextend (command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); + strextend (command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX); + strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS); + + if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); if (TimeSelect) { Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36681) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 36833) @@ -1,15 +1,16 @@ # include "relastro.h" -Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog); +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile); void bcatalog_show_skips (); -Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) { - - int i, j, k, m, Nstar; +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile) { + + int i, j, Nstar; + // int k, m; Catalog *catalog, *pcatalog, tcatalog; // XXX need to decide how to determine PARALLEL mode... if (PARALLEL && !hostID) { - catalog = load_catalogs_parallel (skylist, Ncatalog); + catalog = load_catalogs_parallel (skylist, Ncatalog, syncfile); return catalog; } @@ -73,4 +74,5 @@ if (RESET) { for (j = 0; j < catalog[Ncat].Naverage; j++) { +# if (0) catalog[Ncat].average[j].flags = 0; m = catalog[Ncat].average[j].measureOffset; @@ -78,4 +80,5 @@ catalog[Ncat].measure[m+k].dbFlags = 0; } +# endif } } @@ -84,5 +87,5 @@ } - // XXX TEST : bcatalog_show_skips(); + bcatalog_show_skips(); Nstar = 0; @@ -113,5 +116,5 @@ // CATDIR is supplied globally # define DEBUG 1 -Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile) { char uniquer[12]; @@ -150,24 +153,26 @@ table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM); - char tmpline[1024]; - if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } - if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } - if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } - - if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } - if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } - if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } - if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_ONLY) strextend (command, "-pm"); + if (FIT_MODE == FIT_PAR_ONLY) strextend (command, "-par"); + if (FIT_MODE == FIT_PM_AND_PAR) strextend (command, "-pmpar"); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + if (FlagOutlier) strextend (command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); - if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); strcpy (command, tmpline); } - - if (USE_FIXED_PIXCOORDS) { snprintf (tmpline, 1024, "%s -D USE_FIXED_PIXCOORDS 1", command); strcpy (command, tmpline); } - - if (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } - if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } - if (PhotFlagSelect) { snprintf (tmpline, 1024, "%s +photflags", command); strcpy (command, tmpline); } - if (PhotFlagBad) { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad); strcpy (command, tmpline); } - if (PhotFlagPoor) { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor); strcpy (command, tmpline); } + if (DCR_COLOR_POS && DCR_COLOR_NEG) { + strextend (command, "-dcr-color %s %s", DCR_COLOR_POS, DCR_COLOR_NEG); + } + + if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values @@ -175,8 +180,7 @@ char *tstart = ohana_sec_to_date (TSTART); char *tstop = ohana_sec_to_date (TSTOP); - snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + strextend (command, "-time %s %s", tstart, tstop); free (tstart); free (tstop); - strcpy (command, tmpline); } @@ -211,4 +215,13 @@ } + // update syncfile here (save lots of I/O time) + + // NOTE: if I let all hosts load blindly, I saturate the data clients with too many + // relastro_client requests. I need to have the master mediate this. the master + // will not launch the next remote job until this one says it is done + if (syncfile) { + update_sync_file (syncfile, 1); + } + // each host generates a BrightCatalog structure, with the measure, average, etc value // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt). I need to Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 36681) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 36833) @@ -33,5 +33,5 @@ // USAGE: relastro_client -load-objects int Ncatalog; - Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); + Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR, NULL); if (!catalog) { fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 36681) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 36833) @@ -28,6 +28,14 @@ // XXX pass in the image table // XXX who carries the image grid? - catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + + // photcodesKeep is used here to allow measurements from the images being calibrated + // note if -reset-to-photcode is selected, photocodesKeep is replaced with below with photcodesReset + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, NULL); MARKTIME("load catalog data: %f sec\n", dtime); + + if (photcodesReset) { + photcodesKeep = photcodesReset; + NphotcodesKeep = NphotcodesReset; + } if (Ncatalog == 0) { @@ -50,4 +58,6 @@ // set test points based on the starmap createStarMap (catalog, Ncatalog); + + // XXX NOTE : for 2mass reset, photcodesKeep should now limit to 2MASS measurements /* major modes */ Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36681) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 36833) @@ -40,4 +40,6 @@ snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); catalog.filename = hostID ? hostfile : skylist[0].filename[i]; + + // set up the basic catalog info 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 Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36681) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 36833) @@ -45,13 +45,14 @@ /* load catalog data from region files (hostID is 0 since we are not a client */ - catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, syncfile); MARKTIME("-- load catalog data: %f sec\n", dtime); + free (syncfile); - // NOTE: if I let all hosts load blindly, I saturate the data clients with too many - // relastro_client requests. I need to have the master mediate this. the master - // will not launch the next remote job until this one says it is done - char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); - update_sync_file (syncfile, 1); - + if (photcodesReset) { + photcodesKeep = photcodesReset; + NphotcodesKeep = NphotcodesReset; + } + // generate tables go from catID,objID -> catSeq,objSeq indexCatalogs (catalog, Ncatalog); Index: trunk/Ohana/src/uniphot/src/update_catalog_fiximids.c =================================================================== --- trunk/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 36681) +++ trunk/Ohana/src/uniphot/src/update_catalog_fiximids.c (revision 36833) @@ -1,3 +1,8 @@ # include "fiximids.h" + +// test image: 2013/06/15,13:25:51, GPC1.r.XY50 +static int CHECK_TEST_IMAGE = FALSE; +static unsigned int Tref = 1378812312; +static short Cref = 10001; void update_catalog_fiximids (Catalog *catalog) { @@ -18,4 +23,8 @@ short photcode = catalog[0].measure[i].photcode; e_time time = catalog[0].measure[i].t; + + if (CHECK_TEST_IMAGE && (abs(time - Tref) < 10) && (photcode == Cref)) { + fprintf (stderr, "."); + } // skip detections with no valid imageID (eg, ref photcode)