Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/Makefile (revision 33448) @@ -1,3 +1,3 @@ -default: relastro +default: relastro relastro_client help: @echo "make options: relastro (default)" @@ -18,5 +18,7 @@ relastro: $(BIN)/relastro.$(ARCH) -install: $(DESTBIN)/relastro $(DESTBIN)/testparallax +relastro_client: $(BIN)/relastro_client.$(ARCH) + +install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client $(DESTBIN)/testparallax RELASTRO = \ @@ -53,7 +55,7 @@ $(SRC)/select_images.$(ARCH).o \ $(SRC)/relastro.$(ARCH).o \ +$(SRC)/relastro_images.$(ARCH).o \ $(SRC)/relastro_objects.$(ARCH).o \ $(SRC)/save_catalogs.$(ARCH).o \ -$(SRC)/write_coords.$(ARCH).o \ $(SRC)/CoordOps.$(ARCH).o \ $(SRC)/FixProblemImages.$(ARCH).o \ @@ -64,8 +66,62 @@ $(SRC)/relastro_merge_source.$(ARCH).o \ $(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o + +#$(SRC)/write_coords.$(ARCH).o \ $(RELASTRO): $(INC)/relastro.h $(KAPA_INCS) $(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS) + +RELASTRO_CLIENT = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPar.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/UpdateObjects.$(ARCH).o \ +$(SRC)/UpdateObjectOffsets.$(ARCH).o \ +$(SRC)/UpdateMeasures.$(ARCH).o \ +$(SRC)/GetAstromError.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/dvo_astrom_ops.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ +$(SRC)/relastro_client.$(ARCH).o \ +$(SRC)/relastro_objects.$(ARCH).o \ +$(SRC)/save_catalogs.$(ARCH).o \ +$(SRC)/CoordOps.$(ARCH).o \ +$(SRC)/high_speed_catalogs.$(ARCH).o \ +$(SRC)/high_speed_objects.$(ARCH).o \ +$(SRC)/high_speed_utils.$(ARCH).o \ +$(SRC)/relastro_merge_source.$(ARCH).o \ +$(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/plot_scatter.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/relastroVisual.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o + +#$(SRC)/FitChip.$(ARCH).o \ +#$(SRC)/FitMosaic.$(ARCH).o \ +#$(SRC)/UpdateChips.$(ARCH).o \ +#$(SRC)/UpdateMosaic.$(ARCH).o \ +#$(SRC)/UpdateSimple.$(ARCH).o \ +#$(SRC)/write_coords.$(ARCH).o \ +#$(SRC)/FixProblemImages.$(ARCH).o \ +#$(SRC)/StarMaps.$(ARCH).o \ + +$(RELASTRO_CLIENT): $(INC)/relastro.h $(KAPA_INCS) +$(BIN)/relastro_client.$(ARCH): $(RELASTRO_CLIENT) $(KAPA_LIBS) testparallax: $(BIN)/testparallax.$(ARCH) Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h (revision 33448) @@ -15,5 +15,25 @@ typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode; -typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE} FitTarget; +typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED, TARGET_MERGE_SOURCE, TARGET_UPDATE_OBJECTS, TARGET_UPDATE_OFFSETS, TARGET_LOAD_OBJECTS} FitTarget; + +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; typedef struct { @@ -86,12 +106,17 @@ /* global variables set in parameter file */ -char ImageCat[256]; -char GSCFILE[256]; -char CATDIR[256]; +# define MAX_PATH_LENGTH 1024 +char ImageCat[MAX_PATH_LENGTH]; +char GSCFILE[MAX_PATH_LENGTH]; +char CATDIR[MAX_PATH_LENGTH]; char *HIGH_SPEED_DIR; char CATMODE[16]; /* raw, mef, split, mysql */ char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ -char SKY_TABLE[256]; +char SKY_TABLE[MAX_PATH_LENGTH]; int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +int HOST_ID; +char *HOSTDIR; +char *BCATALOG; unsigned int OBJ_ID_SRC; @@ -109,4 +134,8 @@ double RADIUS; // match radius for high-speed objects + +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; int VERBOSE; @@ -129,7 +158,4 @@ double MaxDensityValue; -int UserCatalog; -double UserCatalogRA, UserCatalogDEC; - char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST; int NphotcodesKeep, NphotcodesSkip; @@ -142,12 +168,9 @@ SkyRegionSelection SELECTION; -int AreaSelect; -double AreaXmin, AreaXmax, AreaYmin, AreaYmax; - -int ImagSelect, ImagMin, ImagMax; +int ImagSelect; +double ImagMin, ImagMax; double PlotMmin, PlotMmax, PlotdMmin, PlotdMmax; -int DophotSelect, DophotValue; int PhotFlagSelect, PhotFlagPoor, PhotFlagBad; @@ -181,6 +204,7 @@ char *GetPhotnamebyCode PROTO((PhotCodeData *photcodes, int code)); void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); -off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); +off_t *SelectRefMosaic PROTO((Mosaic **refmosaic, off_t *Nimage)); int args PROTO((int argc, char **argv)); +int args_client PROTO((int argc, char **argv)); int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); void clean_images PROTO((void)); @@ -226,8 +250,9 @@ void initMrel PROTO((Catalog *catalog, int Ncatalog)); void initialize PROTO((int argc, char **argv)); +void initialize_client PROTO((int argc, char **argv)); void initstats PROTO((char *mode)); int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); -Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect)); -SkyList *load_images PROTO((FITS_DB *db, SkyRegion *region)); +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)); @@ -314,4 +339,5 @@ void FlagOutliers(Catalog *catalog); int MeasFilterTest(Measure *measure, int applySigmaLim); +int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim); int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius); @@ -333,12 +359,16 @@ Mosaic *getMosaicForImage (off_t im); -double getMeanR (Measure *measure, Average *average, SecFilt *secfilt); -double getMeanD (Measure *measure, Average *average, SecFilt *secfilt); -int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt); -int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt); +double getMeanR (MeasureTiny *measure, Average *average, SecFilt *secfilt); +double getMeanD (MeasureTiny *measure, Average *average, SecFilt *secfilt); +int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt); +int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt); +int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt); +int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt); float GetAstromError (Measure *measure, int mode); -int relastro_objects (void); -int UpdateObjectOffsets (SkyList *skylist); +float GetAstromErrorTiny (MeasureTiny *measure, int mode); +int relastro_objects (SkyList *skylist, int hostID, char *hostpath); +int relastro_images (SkyList *skylist); +int UpdateObjectOffsets (SkyList *skylist, int hostID, char *hostpath); int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj); @@ -361,5 +391,5 @@ void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im); -int high_speed_catalogs (); +int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath); int high_speed_objects (SkyRegion *region, Catalog *catalog); int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset); @@ -369,4 +399,5 @@ int createStarMapPoints(); int checkStarMap(int N); +int createStarMap (Catalog *catalog, int Ncatalog); int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit); @@ -378,4 +409,11 @@ void setupAreaSelection(SkyRegion *region); -int relastro_merge_source (); +int relastro_merge_source (SkyTable *sky); void resort_catalog (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); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/BrightCatalog.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/BrightCatalog.c (revision 33448) @@ -1,9 +1,11 @@ # include "relastro.h" -// relphot_client reads from the local catalogs and generates a set of bright, subset catalogs +// NOTE: this file is nearly identical to relphot/src/BrightCatalog.c, but here we use Average, not AverageTiny + +// relastro_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, AverageTiny, and Secfilt tables, but with cat_id appended to each row so we can +// MeasureTiny, Average, and Secfilt tables, but with cat_id appended to each row so we can // reassign correctly on read @@ -131,32 +133,100 @@ // read the fits table bytes if (!gfits_fread_ftable_data (f, &ftable, FALSE)) return (NULL); - - // need to create and assign to flat-field correction - GET_COLUMN(R, "RA", double); - GET_COLUMN(D, "DEC", double); - GET_COLUMN(Nmeasure, "NMEAS", int); - GET_COLUMN(measureOffset, "MEASURE_OFF", int); - GET_COLUMN(flags, "FLAGS", int); - GET_COLUMN(catID, "CAT_ID", int); + + 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(Xp, "SIGMA_POS", 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); // XXX free the fits table data here - AverageTiny *average = NULL; - ALLOCATE (average, AverageTiny, Nrow); + 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].Nmeasure = Nmeasure[i]; - average[i].measureOffset = measureOffset[i]; - average[i].flags = flags[i]; - average[i].catID = catID[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].Xp = Xp[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 (Nmeasure ); - free (measureOffset ); - free (flags ); - free (catID ); + 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 (Xp); + 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; @@ -257,4 +327,5 @@ 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); @@ -343,14 +414,37 @@ } - /*** AverageTiny ***/ + /*** Average ***/ { - gfits_create_table_header (&theader, "BINTABLE", "AVERAGE_TINY"); - - gfits_define_bintable_column (&theader, "D", "RA", "ra (J2000)", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "DEC", "dec (J2000)", "degree", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NMEAS", "number of measures", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "MEASURE_OFF", "index to measurements", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ref", NULL, 1.0, 0.0); + 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", "SIGMA_POS", "position scatter", "", 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 @@ -358,36 +452,128 @@ // create intermediate storage arrays - double *R ; ALLOCATE (R, double, catalog->Naverage); - double *D ; ALLOCATE (D, double, catalog->Naverage); - int *Nmeasure ; ALLOCATE (Nmeasure, int, catalog->Naverage); - int *measureOffset ; ALLOCATE (measureOffset, int, catalog->Naverage); - int *flags ; ALLOCATE (flags, int, catalog->Naverage); - int *catID ; ALLOCATE (catID, int, catalog->Naverage); + 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 *Xp ; ALLOCATE (Xp , 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 - AverageTiny *average = catalog->average; + Average *average = catalog->average; for (i = 0; i < catalog->Naverage; i++) { - R[i] = average[i].R ; - D[i] = average[i].D ; - Nmeasure[i] = average[i].Nmeasure; - measureOffset[i] = average[i].measureOffset; - flags[i] = average[i].flags; - catID[i] = average[i].catID; + 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 ; + Xp[i] = average[i].Xp ; + 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, "NMEAS", Nmeasure, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "MEASURE_OFF", measureOffset, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, catalog->Naverage); - gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Naverage); - - free (R ); - free (D ); - free (Nmeasure ); - free (measureOffset ); - free (flags ); - free (catID ); + 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, "SIGMA_POS", Xp, 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 (Xp); + 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); @@ -486,7 +672,7 @@ ALLOCATE (bcatalog[0].measure, MeasureTiny, Nmeas); - ALLOCATE (bcatalog[0].average, AverageTiny, Naves); - ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt); - + ALLOCATE (bcatalog[0].average, Average, Naves); + ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt); + int Nm = 0; int Na = 0; @@ -494,7 +680,5 @@ if (!catalog[i].Naverage) continue; for (j = 0; j < catalog[i].Naverage; j++) { - // CopyAverageTiny (&bcatalog[0].average[Na], &catalog[i].average[j]); - bcatalog[0].average[Na] = catalog[i].averageT[j]; - bcatalog[0].average[Na].catID = catalog[i].catID; + 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]; @@ -504,7 +688,5 @@ } for (j = 0; j < catalog[i].Nmeasure; j++) { - // CopyMeasureTiny (&bcatalog[0].measure[Nm], &catalog[i].measure[j]); bcatalog[0].measure[Nm] = catalog[i].measureT[j]; - bcatalog[0].measure[Nm].catID = catalog[i].catID; Nm++; assert (Nm <= Nmeas); @@ -548,5 +730,5 @@ catalogs->catalog[i].Naverage = 0; catalogs->catalog[i].Nmeasure = 0; - ALLOCATE (catalogs->catalog[i].averageT, AverageTiny, catalogs->NAVERAGE[i]); + 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); @@ -557,12 +739,4 @@ // distribute a bright catalog across separate catalogs int BrightCatalogSplitFree (CatalogSplitter *catalogs) { - - // XXX don't free the catalogs : - // for (i = 0; i < NCATALOG; i++) { - // free (catalogs->catalog[i].averageT); - // free (catalogs->catalog[i].measureT); - // free (catalogs->catalog[i].secfilt); - // } - // free (catalogs->catalog); free (catalogs->catIDs); @@ -609,4 +783,7 @@ assert (Ncat < catalogs->NCATALOG); + catalogs->catalog[Ncat].Nsecfilt = Nsecfilt; + catalogs->catalog[Ncat].catID = catID; + catalogs->index[catID] = Ncat; assert (catID > 0); @@ -639,5 +816,5 @@ catalogs->catalog[j].Naverage = 0; catalogs->catalog[j].Nmeasure = 0; - ALLOCATE (catalogs->catalog[j].averageT, AverageTiny, catalogs->NAVERAGE[j]); + 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); @@ -658,5 +835,5 @@ int Na = catalogs->catalog[Nc].Naverage; - catalogs->catalog[Nc].averageT[Na] = bcatalog->average[i]; + catalogs->catalog[Nc].average[Na] = bcatalog->average[i]; // secfilt entries are grouped in blocks for each average @@ -669,6 +846,6 @@ if (catalogs->catalog[Nc].Naverage >= catalogs->NAVERAGE[Nc]) { catalogs->NAVERAGE[Nc] += 100; - REALLOCATE (catalogs->catalog[Nc].averageT, AverageTiny, catalogs->NAVERAGE[Nc]); - REALLOCATE (catalogs->catalog[Nc].secfilt, SecFilt, catalogs->NAVERAGE[Nc]*Nsecfilt); + REALLOCATE (catalogs->catalog[Nc].average, Average, catalogs->NAVERAGE[Nc]); + REALLOCATE (catalogs->catalog[Nc].secfilt, SecFilt, catalogs->NAVERAGE[Nc]*Nsecfilt); } } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ConfigInit.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ConfigInit.c (revision 33448) @@ -26,11 +26,4 @@ if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; - // XXX these are used in relphot to identify poor stars / images -- define - // an equivalent concept here? - // GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER); - // GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER); - // GetConfig (config, "IMAGE_OFFSET", "%lf", 0, &IMAGE_OFFSET); - // GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ); - GetConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN); GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); @@ -42,6 +35,11 @@ GetConfig (config, "ADDSTAR_RADIUS", "%lf", 0, &ADDSTAR_RADIUS); + // force CATDIR to be absolute (so parallel mode will work) + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + char *tmpcatdir = abspath (CATDIR, MAX_PATH_LENGTH); + strcpy (CATDIR, tmpcatdir); + free (tmpcatdir); + GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); - GetConfig (config, "CATDIR", "%s", 0, CATDIR); ScanConfig(config, "CATMODE", "%s", 0, CATMODE); ScanConfig(config, "CATFORMAT", "%s", 0, CATFORMAT); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/FixProblemImages.c (revision 33448) @@ -67,5 +67,6 @@ } - catalog = load_catalogs (&sublist, &Ncat, FALSE); + // XXX use a different function here + // catalog = load_catalogs (&sublist, &Ncat, FALSE); assert (Ncat == Ncatlist); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/GetAstromError.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/GetAstromError.c (revision 33448) @@ -1,4 +1,42 @@ # include "relastro.h" # define WEIGHTED_ERRORS 1 + +float GetAstromErrorTiny (MeasureTiny *measure, int mode) { + + PhotCode *code; + float dPobs, dPsys, dPtotal, dM, AS, MS, dX, dY; + + if (!WEIGHTED_ERRORS) { + // if we don't understand the errors at all, this at least lets us get things roughly + // right: + return 0.1; + } + + switch (mode) { + case ERROR_MODE_RA: + dPobs = FromShortPixels(measure[0].dXccd); // dXccd is a value in pixels + break; + case ERROR_MODE_DEC: + dPobs = FromShortPixels(measure[0].dYccd); // dYccd is a value in pixels + break; + case ERROR_MODE_POS: + dX = FromShortPixels(measure[0].dXccd); // dXccd is a value in pixels + dY = FromShortPixels(measure[0].dYccd); // dYccd is a value in pixels + dPobs = hypot (dX, dY); + break; + default: + abort(); + } + + code = GetPhotcodebyCode (measure[0].photcode); + AS = code[0].astromErrScale; + MS = code[0].astromErrMagScale; + dPsys = code[0].astromErrSys; + dM = measure[0].dM; + dPtotal = sqrt(SQ(dPsys) + SQ(AS*dPobs) + SQ(MS*dM)); + + dPtotal = MAX (dPtotal, MIN_ERROR); + return (dPtotal); +} float GetAstromError (Measure *measure, int mode) { Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/ImageOps.c (revision 33448) @@ -333,5 +333,5 @@ } -void dump_measures(Average *average, MeasureTiny *measure) { +void dump_measures (Average *average, Measure *measure) { off_t j, off; @@ -406,5 +406,5 @@ c = clist[im][i]; - MeasureTiny *measure = &catalog[c].measureT[m]; + Measure *measure = &catalog[c].measure[m]; X = measure[0].Xccd; @@ -431,10 +431,10 @@ if (fabs(dR) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measureT); + dump_measures (&catalog[c].average[n], catalog[c].measure); // abort (); } if (fabs(dD) > 3.0*ADDSTAR_RADIUS) { fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measureT); + dump_measures (&catalog[c].average[n], catalog[c].measure); // abort (); } @@ -443,5 +443,5 @@ if (fabs(catalog[c].measureT[m].dR - dR) > DPOS_MAX_ASEC) { fprintf (stderr, "measurement is far from original location (R): %f %f (%f %f)\n", catalog[c].average[n].R, catalog[c].average[n].D, dR, dD); - dump_measures (&catalog[c].average[n], catalog[c].measureT); + dump_measures (&catalog[c].average[n], catalog[c].measure); // abort(); } @@ -581,5 +581,5 @@ c = clist[im][i]; - MeasureTiny *measure = &catalog[c].measure[m]; + MeasureTiny *measure = &catalog[c].measureT[m]; /* apply the current image transformation or use the current value of R+dR, D+dD? */ @@ -589,5 +589,5 @@ raw[i].Mag = measure[0].M; raw[i].dMag = measure[0].dM; - raw[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS); + raw[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); n = measure[0].averef; @@ -654,7 +654,7 @@ m = mlist[im][i]; c = clist[im][i]; - n = catalog[c].measure[m].averef; - - MeasureTiny *measure = &catalog[c].measure[m]; + + MeasureTiny *measure = &catalog[c].measureT[m]; + n = measure[0].averef; /* apply the current image transformation or use the current value of R+dR, D+dD? */ @@ -664,5 +664,5 @@ ref[i].Mag = measure[0].M; ref[i].dMag = measure[0].dM; - ref[i].dPos = GetAstromError (&measure[0], ERROR_MODE_POS); + ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); ref[i].mask = FALSE; @@ -694,4 +694,5 @@ void FlagOutliers2D(Catalog *catalog); +// operates on Full values (not tiny) void FlagOutliers (Catalog *catalog) { @@ -731,5 +732,5 @@ // pointer to this set of measurements m = catalog[0].average[j].measureOffset; - MeasureTiny *measure = &catalog[0].measureT[m]; + Measure *measure = &catalog[0].measure[m]; /* accumulate list of valid measurements */ @@ -825,5 +826,5 @@ // pointer to this set of measurements m = catalog[0].average[j].measureOffset; - MeasureTiny *measure = &catalog[0].measureT[m]; + Measure *measure = &catalog[0].measure[m]; /* accumulate list of valid measurements */ @@ -915,9 +916,7 @@ } - - /** 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(MeasureTiny *measure, int applySigmaLim) { +int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) { int found, k; long mask; @@ -978,2 +977,64 @@ 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 relastro tasks? + if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro 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; + } + + /* 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-20111122/Ohana/src/relastro/src/StarMaps.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/StarMaps.c (revision 33448) @@ -68,5 +68,5 @@ for (i = 0; i < catalog[0].Nmeasure; i++) { - MeasureTiny *measure = &catalog[0].measureT[i] + MeasureTiny *measure = &catalog[0].measureT[i]; N = getImageByID(measure[0].imageID); @@ -84,4 +84,19 @@ return (TRUE); +} + +int createStarMap (Catalog *catalog, int Ncatalog) { + + int i; + + initStarMaps(); + + for (i = 0; i < Ncatalog; i++) { + // check coverage per chip (operates on (Average, MeasureTiny, Secfilt) + updateStarMaps (&catalog[i]); + } + + createStarMapPoints(); + return TRUE; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c (revision 33448) @@ -143,10 +143,5 @@ } - double Rmid; - if (UserCatalog) { - Rmid = UserCatalogRA; - } else { - Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax); - } + double Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax); R = ohana_normalize_angle_to_midpoint (R, Rmid); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateMeasures.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateMeasures.c (revision 33448) @@ -1,4 +1,5 @@ # include "relastro.h" +// this function operates on Measure, not MeasureTiny int UpdateMeasures (Catalog *catalog, int Ncatalog) { Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33448) @@ -7,14 +7,28 @@ // original coordinates -int UpdateObjectOffsets (SkyList *skylist) { +int UpdateObjectOffsets_parallel (SkyList *sky); + +int UpdateObjectOffsets (SkyList *skylist, int hostID, char *hostpath) { int i; Catalog catalog; + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL & !hostID) { + UpdateObjectOffsets_parallel (skylist); + return TRUE; + } + // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { + if (hostID && (skylist[0].regions[i]->hostID != hostID)) continue; + + // define the catalog file name + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : skylist[0].filename[i]; + // set up the basic catalog info - catalog.filename = skylist[0].filename[i]; catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data @@ -33,4 +47,7 @@ } + // set the values in MeasureTiny needed by UpdateObjects + populate_tiny_values(&catalog, DVO_TV_MEASURE); + // match measurements with images initImageBins (&catalog, 1, FALSE); @@ -40,5 +57,5 @@ UpdateMeasures (&catalog, 1); - XXX: note that this version of UpdateObjects is / should be expecting catalog.measure, not catalog.measureT + // uses MeasureTiny values, but will also update Measure values when populated UpdateObjects (&catalog, 1); @@ -51,2 +68,96 @@ return (TRUE); } + +// CATDIR is supplied globally +# define DEBUG 1 +int UpdateObjectOffsets_parallel (SkyList *sky) { + + // launch the setphot_client jobs to the parallel hosts + + // 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, MAX_PATH_LENGTH); + free (table->hosts[i].pathname); + table->hosts[i].pathname = tmppath; + + // options / arguments that can affect relastro_client -load: + // VERBOSE, VERBOSE2 + // RESET (-reset) + // TimeSelect -time + // (note that psfQual is applied rigidly at 0.85, as is the galaxy test) + // ImagSelect, ImagMin, ImagMax + // MaxDensityUse, MaxDensityValue + + char command[1024]; + snprintf (command, 1024, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", + table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); + + char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); 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 relastro_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 relastro_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableWaitJobs (table, __FILE__, __LINE__); + } + + return TRUE; +} Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateObjects.c (revision 33448) @@ -38,4 +38,6 @@ } +// This function operates on both Measure and MeasureTiny. In the big stages, this should +// be called with just MeasureTiny set and Measure == NULL int UpdateObjects (Catalog *catalog, int Ncatalog) { @@ -93,4 +95,6 @@ m = catalog[i].average[j].measureOffset; MeasureTiny *measure = &catalog[i].measureT[m]; + Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL; + // when we update the output measure values, we need to do it here Tmin = Tmax = (measure[0].t - T2000) / (86400*365.25); @@ -102,6 +106,7 @@ // does the measurement pass the supplied filtering constraints? - if (!MeasFilterTest(&measure[k], FALSE)) { + if (!MeasFilterTestTiny(&measure[k], FALSE)) { measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } continue; } @@ -110,8 +115,10 @@ if (FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) { measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } continue; } measure[k].dbFlags |= ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; } R[N] = getMeanR (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); @@ -127,6 +134,6 @@ // dX, dY : error in arcsec -- - dX[N] = GetAstromError (&measure[k], ERROR_MODE_RA); - dY[N] = GetAstromError (&measure[k], ERROR_MODE_DEC); + dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); + dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); // add systematic error in quadrature, if desired @@ -293,9 +300,11 @@ // the measure fields must be updated before the average fields - m = catalog[i].average[j].measureOffset; - MeasureTiny *measure = &catalog[i].measureT[m]; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { setMeanR (fit.Ro, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); setMeanD (fit.Do, &measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + if (measureBig) { + setMeanR_Big (fit.Ro, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + setMeanD_Big (fit.Do, &measureBig[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + } } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c (revision 33448) @@ -1,4 +1,5 @@ # include "relastro.h" void usage (void); +void usage_client (void); void usage_merge_source (void); void usage_merge_source_id (char *name); @@ -42,5 +43,5 @@ if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_OBJECTS; + FIT_TARGET = TARGET_UPDATE_OBJECTS; } @@ -92,5 +93,5 @@ RADIUS = atof(argv[N]); remove_argument (N, &argc, argv); - HIGH_SPEED_DIR = strcreate(argv[N]); + HIGH_SPEED_DIR = abspath(argv[N], MAX_PATH_LENGTH); remove_argument (N, &argc, argv); } @@ -118,14 +119,4 @@ if (FIT_TARGET == TARGET_NONE) usage(); - - UserCatalog = FALSE; - if ((N = get_argument (argc, argv, "-catalog"))) { - remove_argument (N, &argc, argv); - UserCatalogRA = atof(argv[N]); - remove_argument (N, &argc, argv); - UserCatalogDEC = atof(argv[N]); - remove_argument (N, &argc, argv); - UserCatalog = TRUE; - } /* specify portion of the sky : allow default of all sky? */ @@ -144,8 +135,41 @@ UserPatch.Dmax = atof (argv[N]); remove_argument (N, &argc, argv); - } else { - if (!UserCatalog) { - usage (); + } + 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, + // relastro 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 relastro_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); } @@ -290,18 +314,4 @@ } - AreaSelect = FALSE; - if ((N = get_argument (argc, argv, "-area"))) { - remove_argument (N, &argc, argv); - AreaXmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaXmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmin = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaYmax = atof (argv[N]); - remove_argument (N, &argc, argv); - AreaSelect = TRUE; - } - ImagSelect = FALSE; if ((N = get_argument (argc, argv, "-instmag"))) { @@ -336,13 +346,4 @@ } - /* XXX drop this? */ - DophotSelect = FALSE; - if ((N = get_argument (argc, argv, "-dophot"))) { - remove_argument (N, &argc, argv); - DophotValue = atof (argv[N]); - remove_argument (N, &argc, argv); - DophotSelect = TRUE; - } - NLOOP = 4; if ((N = get_argument (argc, argv, "-nloop"))) { @@ -356,9 +357,253 @@ } +int args_client (int argc, char **argv) { + + int N; + double trange; + + /* possible operations */ + FIT_TARGET = TARGET_NONE; + FIT_MODE = FIT_AVERAGE; + + // 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); + FIT_TARGET = TARGET_LOAD_OBJECTS; + } + + if ((N = get_argument (argc, argv, "-update-objects"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_UPDATE_OBJECTS; + } + + if ((N = get_argument (argc, argv, "-update-offsets"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_UPDATE_OFFSETS; + } + + // check for object fitting modes + if ((N = get_argument (argc, argv, "-pm"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_ONLY; + } + if ((N = get_argument (argc, argv, "-par"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PAR_ONLY; + } + if ((N = get_argument (argc, argv, "-pmpar"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_AND_PAR; + } + + if ((N = get_argument (argc, argv, "-high-speed"))) { + // XXX include a parallax / no-parallax option + if (N >= argc - 3) usage_client(); + FIT_TARGET = TARGET_HIGH_SPEED; + remove_argument (N, &argc, argv); + PHOTCODE_A_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + PHOTCODE_B_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + RADIUS = atof(argv[N]); + remove_argument (N, &argc, argv); + HIGH_SPEED_DIR = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + FlagOutlier = FALSE; + if ((N = get_argument (argc, argv, "-clip"))) { + remove_argument (N, &argc, argv); + CLIP_THRESH = atof (argv[N]); + remove_argument (N, &argc, argv); + FlagOutlier = TRUE; + } + + if (FIT_TARGET == TARGET_NONE) usage_client(); + + /* 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); + } + + USE_BASIC_CHECK = FALSE; + if ((N = get_argument (argc, argv, "-basic-image-search"))) { + remove_argument (N, &argc, argv); + USE_BASIC_CHECK = TRUE; + } + + 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_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); + } + + // by default, require > 10pts to clip + strcpy (STATMODE, "CHI_INNER_80_WTMEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + // do not hold the lock on the Image table + UPDATE = FALSE; + + /* 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); + } + + if (argc != 1) usage_client (); + return TRUE; +} + void usage () { - fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n"); - fprintf (stderr, " OR: relastro -catalog (ra) (dec)\n\n"); - fprintf (stderr, " OR: relastro -high-speed [options]\n\n"); + fprintf (stderr, "ERROR: USAGE: relastro -update-simple [options]\n"); + fprintf (stderr, " OR: relastro -update-chips [options]\n"); + fprintf (stderr, " OR: relastro -update-mosaic [options]\n"); + fprintf (stderr, " OR: relastro -update-objects [options]\n"); + fprintf (stderr, " OR: relastro -high-speed [options]\n"); fprintf (stderr, " OR: relastro -merge-source [options]\n\n"); + fprintf (stderr, " specify one of the following modes: \n"); fprintf (stderr, " -update-objects\n"); @@ -369,6 +614,10 @@ fprintf (stderr, " -update-chips\n"); fprintf (stderr, " -update-mosaics\n"); - fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n\n"); - fprintf (stderr, " additional options: \n"); + fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (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"); @@ -383,6 +632,4 @@ fprintf (stderr, " -plrange\n"); fprintf (stderr, " -minerror\n"); - fprintf (stderr, " -area\n"); - fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); fprintf (stderr, " -instmag min max\n\n"); fprintf (stderr, " +photflags\n"); @@ -394,4 +641,38 @@ } +void usage_client () { + fprintf (stderr, "ERROR: USAGE: relastro_client -load\n"); + fprintf (stderr, " OR: relastro_client -update-offsets\n"); + fprintf (stderr, " OR: relastro_client -update-objects\n"); + fprintf (stderr, " OR: relastro_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) { Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/bcatalog.c (revision 33448) @@ -71,5 +71,5 @@ // allowed. - CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; @@ -102,5 +102,6 @@ subcatalog[0].Nmeasure = Nmeasure; subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; - subcatalog[0].Nsecf_mem = Naverage * Nsecfilt; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].catID = catalog[0].catID; assert (Nsecfilt == catalog[0].Nsecfilt); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/dvo_astrom_ops.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 33448) @@ -1,5 +1,5 @@ # include "relastro.h" -double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) { +double getMeanR (MeasureTiny *measure, Average *average, SecFilt *secfilt) { double ra; @@ -21,5 +21,5 @@ } -double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) { +double getMeanD (MeasureTiny *measure, Average *average, SecFilt *secfilt) { double dec; @@ -41,5 +41,5 @@ } -int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) { +int setMeanR (double ra_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt) { /* math to get from new fitted position to new measure offset @@ -66,5 +66,5 @@ } -int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) { +int setMeanD (double dec_fit, MeasureTiny *measure, Average *average, SecFilt *secfilt) { /* math to get from new fitted position to new measure offset @@ -90,2 +90,56 @@ return (TRUE); } + +int setMeanR_Big (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) { + + if (!measure) return TRUE; + + /* math to get from new fitted position to new measure offset + ra_obs = average[0].R - measure[0].dR / 3600.0; + measure[0].dR = (ra_fit - ra_obs) * 3600.0; + measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0 + measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR; + */ + + /* the measure carries the instantaneous mean position at the epoch t */ + measure[0].dR += (ra_fit - average[0].R) * 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (TRUE); +} + +int setMeanD_Big (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) { + + if (!measure) return TRUE; + + /* math to get from new fitted position to new measure offset + dec_obs = average[0].D - measure[0].dD / 3600.0; + measure[0].dD = (dec_fit - dec_obs) * 3600.0; + measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0 + measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD; + */ + + /* the measure carries the instantaneous mean position at the epoch t */ + measure[0].dD += (dec_fit - average[0].D) * 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (TRUE); +} Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_catalogs.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_catalogs.c (revision 33448) @@ -1,10 +1,9 @@ # include "relastro.h" -int high_speed_catalogs () { +int high_speed_catalogs_parallel (SkyList *sky); + +int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath) { int i; - - SkyTable *sky = NULL; - SkyList *skylist = NULL; Catalog catalog; @@ -12,23 +11,10 @@ dvo_set_catdir(CATDIR); - // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, FALSE, SKY_DEPTH, VERBOSE); - SkyTableSetFilenames (sky, CATDIR, "cpt"); - - // determine the populated SkyRegions overlapping the requested area (default depth) - if (UserCatalog) { - skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); - setupAreaSelection(skylist->regions[0]); - } else { - skylist = SkyListByPatch (sky, -1, &UserPatch); - setupAreaSelection(&UserPatch); - } - initializeConstraints(); // XXX need to decide how to determine PARALLEL mode... if (PARALLEL & !hostID) { - high_speed_catalogs (skylist); - return TRUE; + high_speed_catalogs_parallel (skylist); + goto finish; } @@ -64,4 +50,21 @@ } +finish: { + char photcodeFile[1024]; + snprintf (photcodeFile, 1024, "%s/Photcodes.dat", HIGH_SPEED_DIR); + if (!SavePhotcodesFITS (photcodeFile)) { + fprintf (stderr, "error saving photcode table %s\n", photcodeFile); + exit (1); + } + + // need to copy across the Images, SkyTable, and Photcode tables: + char *skyfile = SkyTableFilename (HIGH_SPEED_DIR); + SkyTableSave (sky, skyfile); + + char line[2048]; + snprintf (line, 2048, "cp %s/Images.dat %s/Images.dat", CATDIR, HIGH_SPEED_DIR); + system (line); + } + return (TRUE); } @@ -69,12 +72,12 @@ // CATDIR is supplied globally # define DEBUG 1 -int high_speed_catalogs_parallel (SkyList *sky) { +int high_speed_catalogs_parallel (SkyList *skylist) { // launch the setphot_client jobs to the parallel hosts // load the list of hosts - HostTable *table = HostTableLoad (CATDIR, sky->hosts); + HostTable *table = HostTableLoad (CATDIR, skylist->hosts); if (!table) { - fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR); + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", skylist->hosts, CATDIR); exit (1); } @@ -88,18 +91,43 @@ table->hosts[i].pathname = tmppath; - char catalogFile[512]; - snprintf (catalogFile, 512, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname); - - // options / arguments that can affect relastro_client -update-objects: + // options / arguments that can affect relastro_client -high-speed char command[1024]; - snprintf (command, 1024, "relastro_client -high-speed -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", - table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + snprintf (command, 1024, "relastro_client -high-speed %s %s %f %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", + PHOTCODE_A_LIST, PHOTCODE_B_LIST, RADIUS, HIGH_SPEED_DIR, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); 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 (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -keep-ubercal", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (USE_BASIC_CHECK) { snprintf (tmpline, 1024, "%s -basic-image-search", command); strcpy (command, tmpline); } + + if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); strcpy (command, tmpline); } + + if (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 (WHERE_A[0]) { snprintf (tmpline, 1024, "%s -D WHERE_A \"%s\"", command, WHERE_A); strcpy (command, tmpline); } + if (WHERE_B[0]) { snprintf (tmpline, 1024, "%s -D WHERE_B \"%s\"", command, WHERE_B); strcpy (command, tmpline); } + + 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); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_objects.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/high_speed_objects.c (revision 33448) @@ -6,36 +6,35 @@ int high_speed_objects (SkyRegion *region, Catalog *catalog) { - // XXX seems silly to require region here just to find the catalog bounds / center - - off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE; - int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj; + off_t i, j, m, J, ni, nj, *N1; + off_t Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad; + off_t l, i1, NAVERAGE, NMEASURE; + int *slowMoving, *groupA, *groupB, status, Nmatch, Nmatchmeas, Nepoch, nv[2], Nmatchmeasobj; int foundA, foundB; double *X1, *Y1; double dX, dY, dR, RADIUS2; Coords tcoords; - Catalog catalog1; + Catalog catalogOut; int Nsecfilt; char filename[1024]; - Noff = strlen(CATDIR); - sprintf (filename, "%s%s", HIGH_SPEED_DIR, &catalog[0].filename[Noff]); - fprintf(stderr, "%s\n",filename); - dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/ - catalog1.filename = strcreate(filename); + snprintf (filename, 1024, "%s/%s.cpt", HIGH_SPEED_DIR, region[0].name); + fprintf (stderr, "%s\n",filename); + + dvo_catalog_init (&catalogOut, TRUE); /*initialise new catalogue*/ + catalogOut.filename = strcreate(filename); SIGMA_LIM = 0.0; Nsecfilt = GetPhotcodeNsecfilt(); - Naverage=catalog[0].Naverage; - Nmeasure=catalog[0].Nmeasure; - catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion - catalog1.Nsecfilt = Nsecfilt; - // always load all of the data (if any exists) - catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + // off_t Naverage = catalog[0].Naverage; + // off_t Nmeasure = catalog[0].Nmeasure; + catalogOut.filename = filename; // based on the input name, need to keep everything below the catdir portion + catalogOut.Nsecfilt = Nsecfilt; + catalogOut.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; // load all data - catalog1.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data - catalog1.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - - if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) { + catalogOut.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalogOut.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + + if (!dvo_catalog_open (&catalogOut, region, VERBOSE,"w")) { fprintf (stderr, "ERROR: failure to open catalog file %s\n", filename); @@ -45,7 +44,7 @@ NAVERAGE = 1000; NMEASURE = 10000; - REALLOCATE (catalog1.average, Average, NAVERAGE); - REALLOCATE (catalog1.measure, Measure, NMEASURE); - REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt); + REALLOCATE (catalogOut.average, Average, NAVERAGE); + REALLOCATE (catalogOut.measure, Measure, NMEASURE); + REALLOCATE (catalogOut.secfilt, SecFilt, NAVERAGE*Nsecfilt); // high-speed between different surveys (easier case): @@ -114,7 +113,7 @@ // fprintf (stderr, "skip "OFF_T_FMT" (group A)\n", i); } - - } else if (foundB && !foundA) { - + continue; + } + if (foundB && !foundA) { if (applyConstraintsB(catalog, i)) { groupB[i] = TRUE; @@ -124,4 +123,5 @@ NgroupBbad++; } + continue; } // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow @@ -209,16 +209,16 @@ nv[1]=nj; - catalog1.average[Nmatch]=catalog[0].average[nv[0]]; + catalogOut.average[Nmatch]=catalog[0].average[nv[0]]; /*Loop over index values and set measurements in new catalogue*/ Nmatchmeasobj=0; - catalog1.average[Nmatch].measureOffset=Nmatchmeas; + catalogOut.average[Nmatch].measureOffset=Nmatchmeas; for(l=0;lNsecfilt = 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 (VERBOSE2 && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); - - //outlier rejection + 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); + } + + // outlier rejection (operates on full tables (Average, Measure, Secfilt)) if (FlagOutlier) { FlagOutliers(&tcatalog); } - updateStarMaps (&tcatalog); - // select only the brighter stars if (subselect) { - bcatalog (&catalog[i], &tcatalog); + // 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); } else { if (RESET) { - for (j = 0; j < catalog[i].Naverage; j++) { - catalog[i].average[j].flags = 0; - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { - catalog[i].measure[m+k].dbFlags = 0; + for (j = 0; j < catalog[Ncat].Naverage; j++) { + catalog[Ncat].average[j].flags = 0; + m = catalog[Ncat].average[j].measureOffset; + for (k = 0; k < catalog[Ncat].average[j].Nmeasure; k++) { + catalog[Ncat].measure[m+k].dbFlags = 0; } } } } + Ncat ++; } // XXX keep this test? Nstar = 0; - for (i = 0; i < skylist[0].Nregions; i++) { + for (i = 0; i < Ncat; i++) { Nstar += catalog[i].Naverage; } @@ -78,6 +91,7 @@ } - // XXX consider only returning the populated catalogs - *Ncatalog = skylist[0].Nregions; + // only return the populated catalogs + REALLOCATE (catalog, Catalog, Ncat); + *Ncatalog = Ncat; return (catalog); } @@ -92,8 +106,4 @@ # define DEBUG 1 Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { - - int Nsecfilt = GetPhotcodeNsecfilt (); // set the desired number in case we need to create the catalog - - // launch the setphot_client jobs to the parallel hosts // load the list of hosts @@ -117,26 +127,34 @@ // options / arguments that can affect relastro_client -load: // VERBOSE, VERBOSE2 - // KEEP_UBERCAL // RESET (-reset) // TimeSelect -time - // DophotSelect // (note that psfQual is applied rigidly at 0.85, as is the galaxy test) - // MAG_LIM - // SIGMA_LIM // ImagSelect, ImagMin, ImagMax // MaxDensityUse, MaxDensityValue char command[1024]; - snprintf (command, 1024, "relastro_client %s -load %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -D CAMERA %s -D MAG_LIM %f -D SIGMA_LIM %f", - PhotcodeList, catalogFile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, CAMERA, MAG_LIM, SIGMA_LIM); + snprintf (command, 1024, "relastro_client -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", + catalogFile, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -keep-ubercal", command); strcpy (command, tmpline); } - if (DophotSelect) { snprintf (tmpline, 1024, "%s -dophot %d", command, DophotValue); strcpy (command, tmpline); } if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); 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); @@ -183,4 +201,5 @@ // set up an initial array of catalogs + int Nsecfilt = GetPhotcodeNsecfilt (); CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt); Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/load_images.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/load_images.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/load_images.c (revision 33448) @@ -7,5 +7,5 @@ fprintf (stderr, MSG, __VA_ARGS__); } -SkyList *load_images (FITS_DB *db, SkyRegion *region) { +int load_images (FITS_DB *db, SkyList *skylist) { Image *image, *subset; @@ -14,20 +14,5 @@ struct timeval start, stop; - SkyTable *sky = NULL; - SkyList *skylist = NULL; - gettimeofday (&start, (void *) NULL); - - // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); - SkyTableSetFilenames (sky, CATDIR, "cpt"); - - // determine the populated SkyRegions overlapping the requested area - if (UserCatalog) { - skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); - } else { - skylist = SkyListByPatch (sky, -1, region); - } - MARKTIME(" setup sky: %f sec\n", dtime); // convert database table to internal structure @@ -37,5 +22,4 @@ exit (2); } - MARKTIME(" convert image table: %f sec\n", dtime); @@ -53,5 +37,5 @@ if (!UPDATE) dvo_image_unlock (db); - return (skylist); + return TRUE; } @@ -71,4 +55,4 @@ memcpy (vtable[0].buffer[i], &image[i], Nx); } - return (TRUE); + return TRUE; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro.c (revision 33448) @@ -3,118 +3,37 @@ int main (int argc, char **argv) { - int i, status, Ncatalog; - Catalog *catalog; - FITS_DB db; - struct timeval start, stop; - SkyList *skylist = NULL; - - gettimeofday (&start, (void *) NULL); - /* get configuration info, args */ initialize (argc, argv); - /* the object analysis is a separate process iterating over catalogs */ - if (FIT_TARGET == TARGET_OBJECTS) { - relastro_objects (); - exit (0); - } + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */ - if (FIT_TARGET == TARGET_HIGH_SPEED) { - high_speed_catalogs (); - exit (0); - } + switch (FIT_TARGET) { + case TARGET_UPDATE_OBJECTS: + /* the object analysis is a separate process iterating over catalogs */ + relastro_objects (skylist, 0, NULL); + exit (0); - /* a special method to manually merge unlinked detections of sources togther (not parallel) */ - if (FIT_TARGET == TARGET_MERGE_SOURCE) { - relastro_merge_source (); - exit (0); - } + case TARGET_HIGH_SPEED: + /* high-speed is a 2pt cross-correlation process for linking moving objects (high PM) */ + high_speed_catalogs (sky, skylist, 0, NULL); + exit (0); - /* register database handle with shutdown procedure */ - set_db (&db); + case TARGET_MERGE_SOURCE: + /* a special method to manually merge unlinked detections of sources togther (not parallel) */ + relastro_merge_source (sky); + exit (0); - /* lock and load the image db table */ - status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : 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) */ - skylist = load_images (&db, &UserPatch); - MARKTIME("load images: %f sec\n", dtime); - - initCoords(); - - /* load catalog data from region files : subselect high-quality measurements */ - // XXX pass in the image table - // XXX who carries the image grid? - catalog = load_catalogs (skylist, &Ncatalog, TRUE); - MARKTIME("load catalog data: %f sec\n", dtime); - - /* match measurements with images */ - initImageBins (catalog, Ncatalog, TRUE); - MARKTIME("make image bins: %f sec\n", dtime); - - findImages (catalog, Ncatalog, TRUE); - MARKTIME("set up image indexes: %f sec\n", dtime); - - if (PLOTSTUFF) { - // plot_star_coords (catalog, Ncatalog); - // plot_mosaic_fields (catalog); - } - - // set test points based on the starmap - createStarMapPoints(); - - /* major modes */ - switch (FIT_TARGET) { case TARGET_SIMPLE: - for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); - UpdateSimple (catalog, Ncatalog); - } - break; - case TARGET_CHIPS: - for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); - UpdateChips (catalog, Ncatalog); - MARKTIME("update chips: %f sec\n", dtime); - } - // create summary plots of the process - // relastroVisualSummaryChips(); - break; - case TARGET_MOSAICS: - for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); - UpdateMosaic (catalog, Ncatalog); - } - break; + relastro_images (skylist); + exit (0); default: - fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__); - exit (2); + fprintf (stderr, "impossible!\n"); + abort(); } - - if (!UPDATE) exit (0); - - // free the image / measurement pointers - freeImageBins (1); - - reload_images (&db); - - // iterate over catalogs to make detection coordinates consistant - UpdateObjectOffsets (skylist); - - // iterate over catalogs to make detection coordinates consistant - // FixProblemImages (skylist); - - // save the updated image parameters - dvo_image_update (&db, VERBOSE); - dvo_image_unlock (&db); - - exit (0); + exit (1); } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_client.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_client.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_client.c (revision 33448) @@ -4,14 +4,18 @@ // owned by that host. -// there are two modes: +// there are four modes: -// relastro_client -load : extract the bright catalog subset from the client's tables and -// save to a local FITS table for relastro to load and used +// relastro_client -load-objects : extract the bright catalog subset from the client's +// tables and save to a local FITS table for relastro to +// load and used -// (relastro loads the bright catalog subsets from all clients, then uses this to determine -// the per-image zero points (and potentially the flat-field corrections) +// relastro_client -update-offsets : load image table containing the updated astrometry +// information and determine per-object astrometry -// relastro_client -update : load image table containing the update zero point information -// and apply this to the client's tables +// relastro_client -update-objects : load image table containing the updated astrometry +// information and determine per-object astrometry + +// relastro loads the bright catalog subsets from all clients, then determines the +// astrometry for chips/mosaics, while iteratively improving the per-object mean positions int main (int argc, char **argv) { @@ -20,22 +24,14 @@ initialize_client (argc, argv); - // load the current sky table (layout of all SkyRegions) - SkyTable *sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, TRUE, -1, VERBOSE); - if (!sky) { - fprintf (stderr, "ERROR running loading sky table from %s\n", CATDIR); - exit (2); - } + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); - if (!skylist) { - fprintf (stderr, "ERROR setting up skylist for %s\n", CATDIR); - exit (2); - } - - switch (MODE) { - case MODE_LOAD: { + switch (FIT_TARGET) { + + case TARGET_LOAD_OBJECTS: { + // USAGE: relastro_client -load-objects int Ncatalog; - Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR); + Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); if (!catalog) { fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); @@ -50,27 +46,49 @@ } - case MODE_UPDATE: { - // load the image subset table from the specified location - off_t Nimage; - ImageSubset *image = ImageSubsetLoad (IMAGES, &Nimage); - if (!image) { - fprintf (stderr, "ERROR loading image subset %s\n", CATDIR); - exit (2); - } - - // save the available image information in the static array in ImageOps.c - initImagesSubset (image, NULL, Nimage); - - // load the flat-field correction table from CATDIR - char flatcorrFile[1024]; - snprintf (flatcorrFile, 1024, "%s/flatcorr.fits", CATDIR); - FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrFile, VERBOSE); - - reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR); + case TARGET_UPDATE_OBJECTS: { + // USAGE: relastro_client -update-objects + relastro_objects (skylist, HOST_ID, HOSTDIR); break; } - case MODE_UPDATE_OBJECTS: { - relastro_objects (HOST_ID, HOSTDIR); + case TARGET_HIGH_SPEED: { + // USAGE: relastro_client -high-speed + high_speed_catalogs (sky, skylist, HOST_ID, HOSTDIR); + break; + } + + // XXX loading the images is fairly costly -- see if we can do an image subset + case TARGET_UPDATE_OFFSETS: { + FITS_DB db; + + set_db (&db); + + /* lock and load the image db table */ + int 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); + + /* load regions and images based on specified sky patch (default depth) */ + load_images (&db, skylist); + + // // load the image subset table from the specified location + // off_t Nimage; + // ImageSubset *image = ImageSubsetLoad (IMAGES, &Nimage); + // if (!image) { + // fprintf (stderr, "ERROR loading image subset %s\n", CATDIR); + // exit (2); + // } + // + // // save the available image information in the static array in ImageOps.c + // initImagesSubset (image, NULL, Nimage); + // + // // load the flat-field correction table from CATDIR + // char flatcorrFile[1024]; + // snprintf (flatcorrFile, 1024, "%s/flatcorr.fits", CATDIR); + // FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrFile, VERBOSE); + + UpdateObjectOffsets (skylist, HOST_ID, HOSTDIR); + break; } Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_images.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_images.c (revision 33448) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_images.c (revision 33448) @@ -0,0 +1,107 @@ +# include "relastro.h" + +int relastro_images (SkyList *skylist) { + + int i, status, Ncatalog; + Catalog *catalog; + FITS_DB db; + struct timeval start, stop; + + gettimeofday (&start, (void *) NULL); + + /* register database handle with shutdown procedure */ + set_db (&db); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : 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) */ + load_images (&db, skylist); + MARKTIME("load images: %f sec\n", dtime); + + initCoords(); + + /* load catalog data from region files : subselect high-quality measurements */ + // XXX pass in the image table + // XXX who carries the image grid? + catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL); + MARKTIME("load catalog data: %f sec\n", dtime); + + /* match measurements with images */ + initImageBins (catalog, Ncatalog, TRUE); + MARKTIME("make image bins: %f sec\n", dtime); + + findImages (catalog, Ncatalog, TRUE); + MARKTIME("set up image indexes: %f sec\n", dtime); + + if (PLOTSTUFF) { + // plot_star_coords (catalog, Ncatalog); + // plot_mosaic_fields (catalog); + } + + // set test points based on the starmap + createStarMap (catalog, Ncatalog); + + /* major modes */ + switch (FIT_TARGET) { + case TARGET_SIMPLE: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateSimple (catalog, Ncatalog); + } + break; + + case TARGET_CHIPS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateChips (catalog, Ncatalog); + MARKTIME("update chips: %f sec\n", dtime); + } + // create summary plots of the process + // relastroVisualSummaryChips(); + break; + + case TARGET_MOSAICS: + for (i = 0; i < NLOOP; i++) { + UpdateObjects (catalog, Ncatalog); + UpdateMosaic (catalog, Ncatalog); + } + break; + + default: + fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__); + exit (2); + } + + if (!UPDATE) exit (0); + + // free the image / measurement pointers + freeImageBins (1); + + // XXX why do I do this exactly? + reload_images (&db); + + if (PARALLEL) { + // save the updated image parameters + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + } + + // iterate over catalogs to make detection coordinates consistant + UpdateObjectOffsets (skylist, 0, NULL); + + // iterate over catalogs to make detection coordinates consistant + // FixProblemImages (skylist); + + if (!PARALLEL) { + // save the updated image parameters + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + } + + exit (0); +} Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_merge_source.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_merge_source.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_merge_source.c (revision 33448) @@ -1,8 +1,7 @@ # include "relastro.h" -int relastro_merge_source () { +int relastro_merge_source (SkyTable *sky) { int i, m; - SkyTable *sky = NULL; Catalog catalog_src, catalog_dst; @@ -13,8 +12,4 @@ dvo_set_catdir(CATDIR); - // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, NULL, NULL, FALSE, SKY_DEPTH, VERBOSE); - SkyTableSetFilenames (sky, CATDIR, "cpt"); - // we are merging the detections of the src into the dst : OBJ_ID_SRC,CAT_ID_SRC -> OBJ_ID_DST,CAT_ID_DST Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_objects.c (revision 33447) +++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/relastro_objects.c (revision 33448) @@ -3,26 +3,13 @@ int relastro_objects_parallel (SkyList *sky); -int relastro_objects (int hostID, char *hostpath) { +int relastro_objects (SkyList *skylist, int hostID, char *hostpath) { int i, j, k, m; - SkyTable *sky = NULL; - SkyList *skylist = NULL; Catalog catalog; - - // load the current sky table (layout of all SkyRegions) - sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); - SkyTableSetFilenames (sky, CATDIR, "cpt"); - - // determine the populated SkyRegions overlapping the requested area (default depth) - if (UserCatalog) { - skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC); - } else { - skylist = SkyListByPatch (sky, -1, &UserPatch); - } // XXX need to decide how to determine PARALLEL mode... if (PARALLEL & !hostID) { - relphot_objects_parallel (skylist); + relastro_objects_parallel (skylist); return TRUE; } @@ -42,4 +29,5 @@ catalog.Nsecfilt = GetPhotcodeNsecfilt (); + // loads Average, Measure, SecFilt if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); @@ -53,4 +41,5 @@ } + // operates on full tables (Average, Measure, Secfilt) if (FlagOutlier) { FlagOutliers(&catalog); @@ -68,5 +57,9 @@ } + populate_tiny_values(&catalog, DVO_TV_MEASURE); + UpdateObjects (&catalog, 1); + + free_tiny_values(&catalog); if (!UPDATE) { @@ -103,15 +96,9 @@ table->hosts[i].pathname = tmppath; - char catalogFile[512]; - snprintf (catalogFile, 512, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname); - // options / arguments that can affect relastro_client -update-objects: // VERBOSE, VERBOSE2 // RESET (-reset) // TimeSelect -time - // DophotSelect // (note that psfQual is applied rigidly at 0.85, as is the galaxy test) - // MAG_LIM - // SIGMA_LIM // ImagSelect, ImagMin, ImagMax // MaxDensityUse, MaxDensityValue @@ -127,8 +114,32 @@ char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } - if (KEEP_UBERCAL) { snprintf (tmpline, 1024, "%s -keep-ubercal", command); strcpy (command, tmpline); } + + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (FlagOutlier) { snprintf (tmpline, 1024, "%s -clip %d", command, CLIP_THRESH); 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);