Index: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/Makefile (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/Makefile (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/include/checkastro.h =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/include/checkastro.h (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/include/checkastro.h (revision 36793) @@ -0,0 +1,158 @@ +# 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; +float MaxMeanOffset; + +int TimeSelect; +time_t TSTART, TSTOP; + +int FlagOutlier; +int CLIP_THRESH; + +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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/.mana =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/.mana (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/.mana (revision 36793) @@ -0,0 +1,14 @@ +echo {213983744.99843 / 261.0 / 86400.0} +echo {7 + 2*8 + 2*4 + 3*36} +echo {261*12*3600 + 139*12*3600} +$cpusec_day = 261*12*3600 + 139*12*3600 +echo {213983744.99843 / $cpusec_day} +fprintf %e {261*12*3600 + 139*12*3600} +echo {3x3} +echo {3*3} +echo {3*3 / 7} +echo {2*2 / 7} +echo {sqrt(7)} +echo {4800*0.257/60} +echo {4800*0.257/60 / sqrt(7)} +pwd Index: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/BrightCatalog.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/BrightCatalog.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/BrightCatalog.c (revision 36793) @@ -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(extendOffset, "OFF_EXTEND", int); + 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].extendOffset = extendOffset[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 (extendOffset); + 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 *extendOffset ; ALLOCATE (extendOffset , int , 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 ; + extendOffset[i] = average[i].extendOffset ; + 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, "OFF_EXTEND", extendOffset, 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 (extendOffset); + 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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ConfigInit.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ConfigInit.c (revision 36793) @@ -0,0 +1,71 @@ +# 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; + + // 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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ImageOps.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/ImageOps.c (revision 36793) @@ -0,0 +1,233 @@ +# 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); + } + } + + // 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); + + if (VERBOSE2) { + char *date = ohana_sec_to_date(image[i].tzero); + fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" total and "OFF_T_FMT" ref of %d measures (%s, %s)\n", + i, image[i].imageID, image[i].name, 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 ++; + } + } + 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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/MeasFilterTest.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/MeasFilterTest.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/MeasFilterTest.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/SetSignals.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/SetSignals.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/SetSignals.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/Shutdown.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/Shutdown.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/Shutdown.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/args.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/args.c (revision 36793) @@ -0,0 +1,466 @@ +# 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); + } + + MaxMeanOffset = 10.0; + if ((N = get_argument (argc, argv, "-max-mean-offset"))) { + remove_argument (N, &argc, argv); + MaxMeanOffset = 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 -images -update-simple [options]\n"); + fprintf (stderr, " OR: checkastro -images -update-chips [options]\n"); + fprintf (stderr, " OR: checkastro -images -update-mosaic [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-regions -update-simple [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-regions -update-chips [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-regions -update-mosaic [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-images -update-simple [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-images -update-chips [options]\n"); + fprintf (stderr, " OR: checkastro -parallel-images -update-mosaic [options]\n"); + fprintf (stderr, " OR: checkastro -update-objects [options]\n"); + fprintf (stderr, " OR: checkastro -high-speed [options]\n"); + fprintf (stderr, " OR: checkastro -hpm [options]\n"); + fprintf (stderr, " OR: checkastro -merge-source [options]\n\n"); + + fprintf (stderr, " specify one of the following modes: \n"); + fprintf (stderr, " -update-objects\n"); + fprintf (stderr, " -pm\n"); + fprintf (stderr, " -par\n"); + fprintf (stderr, " -pmpar\n"); + fprintf (stderr, " -update-simple\n"); + fprintf (stderr, " -update-chips\n"); + fprintf (stderr, " -update-mosaics\n"); + fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n"); + fprintf (stderr, " -hpm (radius) (output catdir)\n"); + fprintf (stderr, " -merge-source (objID) (catID) into (objID) (catID)\n\n"); + + fprintf (stderr, " additional 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, " -plot\n"); + fprintf (stderr, " -plotdelay (seconds)\n"); + fprintf (stderr, " -statmode (mode)\n"); + fprintf (stderr, " -reset\n"); + fprintf (stderr, " -nloop (N) : number of image-fit iterations\n"); + fprintf (stderr, " -update : apply new fit to database\n"); + fprintf (stderr, " -params\n"); + fprintf (stderr, " -plrange\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_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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/bcatalog.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/bcatalog.c (revision 36793) @@ -0,0 +1,622 @@ +# 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; + // XXX can't (easily) count this for a given chip + } + + /* 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... + 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); +} + +/* 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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_client.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_client.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_client.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_images.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_images.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/checkastro_images.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/initialize.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/initialize.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/initialize.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_catalogs.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_catalogs.c (revision 36793) @@ -0,0 +1,212 @@ +# 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 (PHOTCODE_KEEP_LIST) { snprintf (tmpline, 1024, "%s +photcode %s", command, PHOTCODE_KEEP_LIST); strcpy (command, tmpline); } + if (PHOTCODE_SKIP_LIST) { snprintf (tmpline, 1024, "%s -photcode %s", command, PHOTCODE_SKIP_LIST); strcpy (command, tmpline); } + if (PhotFlagSelect) { snprintf (tmpline, 1024, "%s +photflags", command); strcpy (command, tmpline); } + if (PhotFlagBad) { snprintf (tmpline, 1024, "%s +photflagbad %d", command, PhotFlagBad); strcpy (command, tmpline); } + if (PhotFlagPoor) { snprintf (tmpline, 1024, "%s +photflagpoor %d", command, PhotFlagPoor); strcpy (command, tmpline); } + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + snprintf (tmpline, 1024, "%s -time %s %s", command, tstart, tstop); + free (tstart); + free (tstop); + strcpy (command, tmpline); + } + + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running 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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_images.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_images.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/load_images.c (revision 36793) @@ -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: /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/select_images.c =================================================================== --- /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/select_images.c (revision 36793) +++ /branches/eam_branches/ipp-20140423/Ohana/src/checkastro/src/select_images.c (revision 36793) @@ -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 > DmaxSkyRegion) continue; + if (Dexp + 1.5 < DminSkyRegion) 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) > RmaxSkyRegion) continue; + if (Rexp + 1.5/cos(Dexp*RAD_DEG) < RminSkyRegion) 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); +} +