Index: /trunk/Ohana/src/relastro/Makefile =================================================================== --- /trunk/Ohana/src/relastro/Makefile (revision 33651) +++ /trunk/Ohana/src/relastro/Makefile (revision 33652) @@ -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 +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,6 +66,72 @@ $(SRC)/relastro_merge_source.$(ARCH).o \ $(SRC)/resort_catalog.$(ARCH).o \ +$(SRC)/BrightCatalog.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o -$(RELASTRO): $(INC)/relastro.h $(KAPA_INCS) -$(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS) +#$(SRC)/write_coords.$(ARCH).o \ + +$(RELASTRO): $(INC)/relastro.h +$(BIN)/relastro.$(ARCH): $(RELASTRO) + +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 +$(BIN)/relastro_client.$(ARCH): $(RELASTRO_CLIENT) + +testparallax: $(BIN)/testparallax.$(ARCH) + +TESTPAR = \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/testparallax.$(ARCH).o + +$(TESTPAR): $(INC)/relastro.h +$(BIN)/testparallax.$(ARCH): $(TESTPAR) Index: /trunk/Ohana/src/relastro/doc/parallel.txt =================================================================== --- /trunk/Ohana/src/relastro/doc/parallel.txt (revision 33652) +++ /trunk/Ohana/src/relastro/doc/parallel.txt (revision 33652) @@ -0,0 +1,91 @@ + +2012.03.05 + + Notes on parallel relastro modifications + + relastro has 4 major modes: + + * relastro_objects + * high_speed_catalogs + * relastro_merge_source + * standard relastro (update chips, etc) + + Of these: + * relastro_merge_source can be skipped for not (not interesting for parallel analysis) + * relastro_objects, high_speed_catalogs, and standard relastro all call UpdateObjects, one way or another. + + The big question is if / how we can use MeasureTiny (or if we need a separate astrometry version) + + * average.d fields we are NOT using in UpdateObjects (28 / 120 bytes) + FIELD Nmissing, NMISSING, unsigned short, number of missings + FIELD Nextend, NEXTEND, unsigned short, number of extended measurements + FIELD missingOffset, OFF_MISSING, int, offset to first missing obs + FIELD extendOffset, OFF_EXTEND, int, offset to first extended measurement + FIELD photFlagsUpper, PHOTFLAGS_U, uint32_t, upper bit of 2 bit summary of per-measure photflags + FIELD photFlagsLower, PHOTFLAGS_L, uint32_t, lower bit of 2 bit summary of per-measure photflags + FIELD extID, EXT_ID, uint64_t, external ID for object (eg PSPS objID) + + * measure.d fields in use in UpdateObjects +* measure[m].t +* measure[m].dbFlags +* measure[m].dt +* measure[0].dR +* measure[0].dD +* measure[0].M +* measure[0].dM +* measure[0].photcode +* measure[0].photFlags +* measure[m].Xccd; +* measure[m].Yccd; +* measure[m].averef; +* measure[0].imageID; +* measure[0].airmass +* measure[0].Mcal (only used in PhotRelTiny) ++ measure[0].dXccd ++ measure[0].dYccd ++ measure[m].dRsys + + float dR; + float dD; + float M; + float Mcal; + float dM; + float airmass; + float Xccd; + float Yccd; + float dt; + int t; + unsigned int averef; + unsigned int imageID; + unsigned int dbFlags; + unsigned int photFlags; +- int catID; // unsigned int? + unsigned short photcode; + +maybe +FIELD Map, M_APER, float, aperture mag, mag +FIELD psfQual, PSF_QF, float, psf coverage/quality factor + +not used +FIELD dMcal, MAG_CAL_ERR, float, systematic calibration error, mag +FIELD az, AZ, float, telescope azimuth +FIELD Sky, SKY_FLUX, float, local estimate of sky flux, counts/sec +FIELD dSky, SKY_FLUX_ERR, float, local estimate of sky flux, counts/sec +FIELD t_msec, TIME_MSEC, unsigned short, time fraction of second, milliseconds +FIELD detID, DET_ID, unsigned int, detection ID +FIELD objID, OBJ_ID, unsigned int, unique ID for object in table +FIELD catID, CAT_ID, unsigned int, unique ID for table in which object was first realized +FIELD psfChisq, PSF_CHISQ, float, psf fit chisq +FIELD psfNdof, PSF_NDOF, int, psf degrees of freedom +FIELD psfNpix, PSF_NPIX, int, psf number of pixels +FIELD crNsigma, CR_NSIGMA, float, Nsigma deviation towards CR +FIELD extNsigma, EXT_NSIGMA, float, Nsigma deviation towards EXT +FIELD FWx, FWHM_MAJOR, short, object fwhm major axis, 1/100 of pixels +FIELD FWy, FWHM_MINOR, short, object fwhm minor axis, 1/100 of pixels +FIELD theta, PSF_THETA, short, angle wrt ccd X dir, (0xffff/360) deg +FIELD Mxx, MXX, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD Mxy, MXY, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD Myy, MYY, short, second moments in pixel coords, 1/100 of pixels^2 +FIELD posangle, POSANGLE, short, position angle sky to chip, (0xffff/360) deg +FIELD pltscale, PLTSCALE, float, plate scale, arcsec/pixel + Index: /trunk/Ohana/src/relastro/doc/plx_factor.pro =================================================================== --- /trunk/Ohana/src/relastro/doc/plx_factor.pro (revision 33652) +++ /trunk/Ohana/src/relastro/doc/plx_factor.pro (revision 33652) @@ -0,0 +1,116 @@ +;+ +; +; PLX_FACTOR +; +; Given a Julian date (scalar or vector), this returns the parallax +; factor for a set of equatorial coordinates (RA,Dec). User must have +; an appropriate Earth ephemeris savefile in the path specified at the +; top of the procedure. The savefile should return the variables: +; +; jd, x, y, z +; +; that specify the position of the Earth in the coordinate system of +; solar system's barycenter. For my savefile, I computed this from +; the JPL ephemeris DE405. +; +; INPUT: +; +; epoch (JD) -- julian date(s) for which to compute the parallax +; factor +; +; ra (deg) -- "catalog" RA of object (unperturbed by parallax, +; aberration, etc.) +; +; de (deg) -- "catalog" Dec of object (unperturbed by parallax, +; aberration, etc.) +; +; NOTE: Epoch input can be an array, with RA and Dec input as +; scalars. However, if RA and Dec are input as arrays (i.e., the +; parallax factor is to be computed for multiple objects), then they +; must have the same dimensions as the epoch array. +; +; OUTPUT: +; +; The output produced by this program are two arrays of parallax +; factors (f_ra, f_de) such that: +; +; RA = RA_0 + f_ra * parallax +; Dec = Dec_0 + f_de * parallax +; +; Note that f_ra and f_de have dimensionless units. Also note that +; the RA correction factor includes a factor of cos(Dec) such that: +; f_ra * parallax = delta(RA) * cos(Dec) +; +; Written by Trent J. Dupuy -- 2010 Jul 14 +; +; REVISIONS +; --------- +; 2011 Jan 25: Added /SPITZER keyword that uses the Spitzer ephemeris +; instead of the Earth ephemeris. +;- +pro plx_factor, epoch, ra_in, de_in, f_ra, f_de, spitzer=spitzer + +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;;;;;;;;;;; SPECIFY PATH FOR EPHEMERIS SAVEFILE HERE ;;;;;;;;;;;;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +savefile = '/Users/tdupuy/data/jpl/'+$ + 'JPL_DE405_EARTH-SOLARBARYCENTER.sav' +if keyword_set(spitzer) then begin + savefile = '/Users/tdupuy/data/jpl/'+$ + 'JPL_SPITZER-SOLARBARYCENTER.sav' + message,'Using spitzer ephemeris.',/continue +endif +restore, savefile +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + +; +; check inputs +; +if n_params() lt 5 then begin + print, 'PLX_FACTOR, epoch (JD), ra (deg), dec (deg), f_ra, f_de' + return +endif +if total(size(ra_in,/dim) ne size(de_in,/dim)) ne 0 then $ + message, 'input RA and Dec dimensions must match each other' +if (size(ra_in,/dim))[0] ne 0 and $ + total(size(ra_in,/dim) ne size(epoch,/dim)) ne 0 then $ + message, '(RA, Dec) not specified as scalars, so their '+$ + 'dimensions must match the input epoch array' + +; +; convert to radians and compute trig functions +; +radeg = 180d/!dpi +ra = ra_in/radeg +de = de_in/radeg +sina = sin(ra) +cosa = cos(ra) +sind = sin(de) +cosd = cos(de) + +; +; interpolate the Earth ephemeris at the input epochs +; +x_in = interpol( x, jd, double(epoch) ) +y_in = interpol( y, jd, double(epoch) ) +z_in = interpol( z, jd, double(epoch) ) +if min(epoch) lt min(jd) or max(epoch) gt max(jd) then $ + message, /continue, 'inputs epochs extend beyond the range '+$ + 'of the ephemeris being used.' + +; +; could also apply an offset here due to the position of the observer +; on the earth being slightly different than the geocenter. as +; described in Sec 6.2.3 of Kovalesky & Seidelmann this effect is +; about 4 uas for a parallax of 0.1" -- i.e., negligible +; + +; +; compute parallax factor (from pg. B28 of Astronomical Almanac 2010, +; also see Sec 6.2.1 of Kovalesky & Seidelmann's Fundamentals +; of Astrometry for a derivation) +; +f_ra = x_in*sina - y_in*cosa +f_de = x_in*cosa*sind + y_in*sina*sind - z_in*cosd + +end Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 33651) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 33652) @@ -4,4 +4,9 @@ # include # include +# include + +// choose off_t or int depending on full-scale relphot analysis resources +// # define IDX_T off_t +# define IDX_T int typedef enum { @@ -15,5 +20,34 @@ 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 enum { + MARK_MEAS_DEFAULT = 0x0000, + MARK_TOO_FEW_MEAS = 0x0001, + MARK_NAN_POS_ERROR = 0x0002, + MARK_NAN_MAG_ERROR = 0x0004, + MARK_BIG_MAG_ERROR = 0x0008, + MARK_BIG_OFFSET = 0x0010, +} MeasurementMask; + +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 +120,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; @@ -110,4 +149,10 @@ double RADIUS; // match radius for high-speed objects +int PARALLEL; +int PARALLEL_MANUAL; +int PARALLEL_SERIAL; + +int PARALLEL_OUTPUT; + int VERBOSE; int VERBOSE2; @@ -115,4 +160,5 @@ int RESET; int NLOOP; +int NTHREADS; int UPDATE; int PLOTSTUFF; @@ -123,4 +169,5 @@ int PM_TOOFEW; double PM_DT_MIN; +double PAR_FACTOR_MIN; int PLOTDELAY; int CHIPORDER; @@ -128,7 +175,4 @@ int MaxDensityUse; double MaxDensityValue; - -int UserCatalog; -double UserCatalogRA, UserCatalogDEC; char *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST; @@ -142,12 +186,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 +222,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 +268,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 +357,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 +377,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 +409,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,7 +417,9 @@ 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); -int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog); +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *catalog); +int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog); int initializeConstraints(); @@ -378,4 +428,13 @@ 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); + +PhotCode **ParsePhotcodeList (char *rawlist, int *nphotcodes, int needAve); Index: /trunk/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- /trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 33652) +++ /trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 33652) @@ -0,0 +1,890 @@ +# include "relastro.h" + +// 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, Average, and Secfilt tables, but with cat_id appended to each row so we can +// reassign correctly on read + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// XXX double check the header extname values for each table +// XXX make sure we free things as we can +// XXX make sure we close files as we can +// XXX handle free and close on error return as well +BrightCatalog *BrightCatalogLoad(char *filename) { + + int i, Ncol; + off_t Nrow; + char type[16]; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + header.buffer = NULL; + matrix.buffer = NULL; + ftable.buffer = NULL; + theader.buffer = NULL; + BrightCatalog *catalog = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + goto escape; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + goto escape; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + goto escape; + } + + ALLOCATE (catalog, BrightCatalog, 1); + catalog->Naverage = 0; + catalog->Nmeasure = 0; + + ftable.header = &theader; + + // *** Measure Data *** + { + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(dR, "RA_OFF", float); + GET_COLUMN(dD, "DEC_OFF", float); + GET_COLUMN(M, "MAG_SYS", float); + GET_COLUMN(Mcal, "MAG_CAL", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(airmass, "AIRMASS", float); + GET_COLUMN(Xccd, "X_CCD", float); + GET_COLUMN(Yccd, "Y_CCD", float); + GET_COLUMN(dt, "EXPTIME", float); + GET_COLUMN(t, "TIME", int); + GET_COLUMN(averef, "AVE_REF", int); // XXX signed vs unsigned? + GET_COLUMN(imageID, "IMAGE_ID", int); + GET_COLUMN(dbFlags, "DB_FLAGS", int); + GET_COLUMN(photFlags, "PHOT_FLAGS", int); + GET_COLUMN(catID, "CAT_ID", int); + GET_COLUMN(photcode, "PHOTCODE", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + MeasureTiny *measure = NULL; + ALLOCATE (measure, MeasureTiny, Nrow); + for (i = 0; i < Nrow; i++) { + measure[i].dR = dR[i]; + measure[i].dD = dD[i]; + measure[i].M = M[i]; + measure[i].Mcal = Mcal[i]; + measure[i].dM = dM[i]; + measure[i].airmass = airmass[i]; + measure[i].Xccd = Xccd[i]; + measure[i].Yccd = Yccd[i]; + measure[i].dt = dt[i]; + measure[i].t = t[i]; + measure[i].averef = averef[i]; + measure[i].imageID = imageID[i]; + measure[i].dbFlags = dbFlags[i]; + measure[i].photFlags = photFlags[i]; + measure[i].catID = catID[i]; + measure[i].photcode = photcode[i]; + } + // fprintf (stderr, "loaded data for %lld measures\n", (long long) Nrow); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + catalog->measure = measure; + catalog->Nmeasure = Nrow; + } + + /*** load Average values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + GET_COLUMN(R, "RA", double); + GET_COLUMN(D, "DEC", double); + GET_COLUMN(dR, "RA_ERR", float); + GET_COLUMN(dD, "DEC_ERR", float); + GET_COLUMN(uR, "U_RA", float); + GET_COLUMN(uD, "U_DEC", float); + GET_COLUMN(duR, "V_RA_ERR", float); + GET_COLUMN(duD, "V_DEC_ERR", float); + GET_COLUMN(P, "PAR", float); + GET_COLUMN(dP, "PAR_ERR", float); + GET_COLUMN(ChiSqAve, "CHISQ_POS", float); + GET_COLUMN(ChiSqPM, "CHISQ_PM", float); + GET_COLUMN(ChiSqPar, "CHISQ_PAP", float); + GET_COLUMN(Tmean, "MEAN_EPOCH", int); + GET_COLUMN(Trange, "TIME_RANGE", int); + GET_COLUMN(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); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + Average *average = NULL; + ALLOCATE (average, Average, Nrow); + for (i = 0; i < Nrow; i++) { + average[i].R = R[i] ; + average[i].D = D[i] ; + average[i].dR = dR[i] ; + average[i].dD = dD[i] ; + average[i].uR = uR[i] ; + average[i].uD = uD[i] ; + average[i].duR = duR[i] ; + average[i].duD = duD[i] ; + average[i].P = P[i] ; + average[i].dP = dP[i] ; + average[i].ChiSqAve = ChiSqAve[i] ; + average[i].ChiSqPM = ChiSqPM[i] ; + average[i].ChiSqPar = ChiSqPar[i] ; + average[i].Tmean = Tmean[i] ; + average[i].Trange = Trange[i] ; + average[i].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 (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; + catalog->Naverage = Nrow; + } + + /*** load Secfilt values ***/ + { + + // load data for this header + if (!gfits_load_header (f, &theader)) goto escape; + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) goto escape; + + // need to create and assign to flat-field correction + GET_COLUMN(M, "MAG", float); + GET_COLUMN(dM, "MAG_ERR", float); + GET_COLUMN(Xm, "MAG_CHI", float); + GET_COLUMN(flags, "FLAGS", int); + GET_COLUMN(Ncode, "NCODE", short); + GET_COLUMN(Nused, "NUSED", short); + GET_COLUMN(M_20, "MAG_20", short); + GET_COLUMN(M_80, "MAG_80", short); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + SecFilt *secfilt = NULL; + ALLOCATE (secfilt, SecFilt, Nrow); + for (i = 0; i < Nrow; i++) { + secfilt[i].M = M[i]; + secfilt[i].dM = dM[i]; + secfilt[i].Xm = Xm[i]; + secfilt[i].flags = flags[i]; + secfilt[i].Ncode = Ncode[i]; + secfilt[i].Nused = Nused[i]; + secfilt[i].M_20 = M_20[i]; + secfilt[i].M_80 = M_80[i]; + } + fprintf (stderr, "loaded data for %lld averages\n", (long long) Nrow); + + free (M ); + free (dM ); + free (Xm ); + free (flags); + free (Ncode); + free (Nused); + free (M_20 ); + free (M_80 ); + catalog->secfilt = secfilt; + // assert Nsecfilt * Naverage = Nrow? + } + + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + + return catalog; + +escape: + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + if (catalog) free (catalog); + + fclose (f); + return NULL; +} + +// we are passed a BrightCatalog structure, write it to a FITS table (3 ext) +int BrightCatalogSave(char *filename, BrightCatalog *catalog) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file for output %s\n", filename); + return FALSE; + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + /*** MeasureTiny ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); + + gfits_define_bintable_column (&theader, "E", "RA_OFF", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_OFF", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "AIRMASS", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "X_CCD", "ccd x coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "Y_CCD", "ccd y coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "EXPTIME", "-2.5 * log (exposure time)", "sec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "TIME", "time of exp", "sec", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "AVE_REF", "pointer to average table", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog", NULL, 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode", NULL, 1.0, 1.0*0x80); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + float *dR ; ALLOCATE (dR , float, catalog->Nmeasure); + float *dD ; ALLOCATE (dD , float, catalog->Nmeasure); + float *M ; ALLOCATE (M , float, catalog->Nmeasure); + float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); + float *dM ; ALLOCATE (dM , float, catalog->Nmeasure); + float *airmass ; ALLOCATE (airmass , float, catalog->Nmeasure); + float *Xccd ; ALLOCATE (Xccd , float, catalog->Nmeasure); + float *Yccd ; ALLOCATE (Yccd , float, catalog->Nmeasure); + float *dt ; ALLOCATE (dt , float, catalog->Nmeasure); + int *t ; ALLOCATE (t , int , catalog->Nmeasure); + int *averef ; ALLOCATE (averef , int , catalog->Nmeasure); + int *imageID ; ALLOCATE (imageID , int , catalog->Nmeasure); + int *dbFlags ; ALLOCATE (dbFlags , int , catalog->Nmeasure); + int *photFlags ; ALLOCATE (photFlags, int , catalog->Nmeasure); + int *catID ; ALLOCATE (catID , int , catalog->Nmeasure); + short *photcode ; ALLOCATE (photcode , short, catalog->Nmeasure); + + // assign the storage arrays + MeasureTiny *measure = catalog->measure; + for (i = 0; i < catalog->Nmeasure; i++) { + dR[i] = measure[i].dR ; + dD[i] = measure[i].dD ; + M[i] = measure[i].M ; + Mcal[i] = measure[i].Mcal ; + dM[i] = measure[i].dM ; + airmass[i] = measure[i].airmass ; + Xccd[i] = measure[i].Xccd ; + Yccd[i] = measure[i].Yccd ; + dt[i] = measure[i].dt ; + t[i] = measure[i].t ; + averef[i] = measure[i].averef ; + catID[i] = measure[i].catID ; + imageID[i] = measure[i].imageID ; + dbFlags[i] = measure[i].dbFlags ; + photFlags[i]= measure[i].photFlags; + photcode[i] = measure[i].photcode ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA_OFF", dR, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DEC_OFF", dD, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AIRMASS", airmass, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "X_CCD", Xccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "Y_CCD", Yccd, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "EXPTIME", dt, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "TIME", t, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "AVE_REF", averef, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DB_FLAGS", dbFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOT_FLAGS", photFlags, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "PHOTCODE", photcode, catalog->Nmeasure); + + free (dR ); + free (dD ); + free (M ); + free (Mcal ); + free (dM ); + free (airmass ); + free (Xccd ); + free (Yccd ); + free (dt ); + free (t ); + free (averef ); + free (imageID ); + free (dbFlags ); + free (photFlags ); + free (catID ); + free (photcode); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** Average ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "AVERAGE"); + + gfits_define_bintable_column (&theader, "D", "RA", "RA", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "DEC", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "RA_ERR", "RA error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DEC_ERR", "DEC error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_RA", "RA*cos(D) proper-motion", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "U_DEC", "DEC proper-motion", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "V_RA_ERR", "RA*cos(D) p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "V_DEC_ERR", "DEC p-m error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PAR", "parallax", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PAR_ERR", "parallax error", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_POS", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_PM", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CHISQ_PAP", "astrometry analysis chisq", "", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "MEAN_EPOCH", "mean epoch (PM-PAR ref)", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "J", "TIME_RANGE", "mean epoch (PM-PAR ref)", "", 1.0, 1.0*0x8000); + gfits_define_bintable_column (&theader, "E", "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 + gfits_create_table (&theader, &ftable); + + // create intermediate storage arrays + double *R ; ALLOCATE (R , double , catalog->Naverage); + double *D ; ALLOCATE (D , double , catalog->Naverage); + float *dR ; ALLOCATE (dR , float , catalog->Naverage); + float *dD ; ALLOCATE (dD , float , catalog->Naverage); + float *uR ; ALLOCATE (uR , float , catalog->Naverage); + float *uD ; ALLOCATE (uD , float , catalog->Naverage); + float *duR ; ALLOCATE (duR , float , catalog->Naverage); + float *duD ; ALLOCATE (duD , float , catalog->Naverage); + float *P ; ALLOCATE (P , float , catalog->Naverage); + float *dP ; ALLOCATE (dP , float , catalog->Naverage); + float *ChiSqAve ; ALLOCATE (ChiSqAve , float , catalog->Naverage); + float *ChiSqPM ; ALLOCATE (ChiSqPM , float , catalog->Naverage); + float *ChiSqPar ; ALLOCATE (ChiSqPar , float , catalog->Naverage); + int *Tmean ; ALLOCATE (Tmean , int , catalog->Naverage); + int *Trange ; ALLOCATE (Trange , int , catalog->Naverage); + float *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 + Average *average = catalog->average; + for (i = 0; i < catalog->Naverage; i++) { + R[i] = average[i].R ; + D[i] = average[i].D ; + dR[i] = average[i].dR ; + dD[i] = average[i].dD ; + uR[i] = average[i].uR ; + uD[i] = average[i].uD ; + duR[i] = average[i].duR ; + duD[i] = average[i].duD ; + P[i] = average[i].P ; + dP[i] = average[i].dP ; + ChiSqAve[i] = average[i].ChiSqAve ; + ChiSqPM[i] = average[i].ChiSqPM ; + ChiSqPar[i] = average[i].ChiSqPar ; + Tmean[i] = average[i].Tmean ; + Trange[i] = average[i].Trange ; + 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, "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); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + + /*** SecFilt ***/ + { + gfits_create_table_header (&theader, "BINTABLE", "SECFILT"); + + gfits_define_bintable_column (&theader, "E", "MAG", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_ERR", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "MAG_CHI", "magnitude (sys)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "magnitude (cal)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NCODE", "magnitude (err)", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "NUSED", "airmass", NULL, 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_20", "ccd x coord", "pix", 1.0, 0.0); + gfits_define_bintable_column (&theader, "I", "MAG_80", "ccd y coord", "pix", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // Nsecfilt is number of average filters, Nsec is total number of values + int Nsecfilt = GetPhotcodeNsecfilt(); + int Nsec = Nsecfilt * catalog->Naverage; + + // create intermediate storage arrays + float *M ; ALLOCATE (M , float, Nsec); + float *dM ; ALLOCATE (dM , float, Nsec); + float *Xm ; ALLOCATE (Xm , float, Nsec); + int *flags ; ALLOCATE (flags , int, Nsec); + short *Ncode ; ALLOCATE (Ncode , short, Nsec); + short *Nused ; ALLOCATE (Nused , short, Nsec); + short *M_20 ; ALLOCATE (M_20 , short, Nsec); + short *M_80 ; ALLOCATE (M_80 , short, Nsec); + + // assign the storage arrays + SecFilt *secfilt = catalog->secfilt; + for (i = 0; i < Nsec; i++) { + M [i] = secfilt[i]. M ; + dM [i] = secfilt[i]. dM ; + Xm [i] = secfilt[i]. Xm ; + flags [i] = secfilt[i]. flags ; + Ncode [i] = secfilt[i]. Ncode ; + Nused [i] = secfilt[i]. Nused ; + M_20 [i] = secfilt[i]. M_20 ; + M_80 [i] = secfilt[i]. M_80 ; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "MAG", M , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_ERR", dM , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_CHI", Xm , Nsec); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NCODE", Ncode, Nsec); + gfits_set_bintable_column (&theader, &ftable, "NUSED", Nused, Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_20", M_20 , Nsec); + gfits_set_bintable_column (&theader, &ftable, "MAG_80", M_80 , Nsec); + + free (M ); + free (dM ); + free (Xm ); + free (flags ); + free (Ncode ); + free (Nused ); + free (M_20 ); + free (M_80 ); + + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + gfits_free_header (&theader); + gfits_free_table (&ftable); + } + return TRUE; +} + +// merge a list of catalogs into a single BrightCatalog array set +BrightCatalog *BrightCatalogMerge (Catalog *catalog, int Ncatalog) { + + off_t i, j, k; + + BrightCatalog *bcatalog = NULL; + ALLOCATE (bcatalog, BrightCatalog, 1); + + int Nmeas = 0; + int Naves = 0; + for (i = 0; i < Ncatalog; i++) { + Nmeas += catalog[i].Nmeasure; + Naves += catalog[i].Naverage; + } + + // XXX this prevents different catalogs from having different Nsecfilt values + int Nsecfilt = GetPhotcodeNsecfilt(); + + ALLOCATE (bcatalog[0].measure, MeasureTiny, Nmeas); + ALLOCATE (bcatalog[0].average, Average, Naves); + ALLOCATE (bcatalog[0].secfilt, SecFilt, Naves*Nsecfilt); + + int Nm = 0; + int Na = 0; + for (i = 0; i < Ncatalog; i++) { + if (!catalog[i].Naverage) continue; + for (j = 0; j < catalog[i].Naverage; j++) { + bcatalog[0].average[Na] = catalog[i].average[j]; + for (k = 0; k < Nsecfilt; k++) { + bcatalog[0].secfilt[Nsecfilt*Na + k] = catalog[i].secfilt[Nsecfilt*j + k]; + } + Na++; + assert (Na <= Naves); + } + for (j = 0; j < catalog[i].Nmeasure; j++) { + bcatalog[0].measure[Nm] = catalog[i].measureT[j]; + Nm++; + assert (Nm <= Nmeas); + } + } + bcatalog->Naverage = Na; + bcatalog->Nmeasure = Nm; + return bcatalog; +} + +// distribute a bright catalog across separate catalogs +CatalogSplitter *BrightCatalogSplitInit (int Nsecfilt) { + + int i; + + CatalogSplitter *catalogs = NULL; + + ALLOCATE (catalogs, CatalogSplitter, 1); + + // as we see new BrightCatalogs, we will update maxID and extend this array + catalogs->maxID = 0; + ALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = 0; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + catalogs->Nsecfilt = Nsecfilt; + + catalogs->Ncatalog = 0; + catalogs->NCATALOG = 16; + catalogs->catalog = NULL; + ALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + + ALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + ALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + ALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + + for (i = 0; i < catalogs->NCATALOG; i++) { + dvo_catalog_init (&catalogs->catalog[i], TRUE); + catalogs->catIDs[i] = 0; + catalogs->NAVERAGE[i] = 100; + catalogs->NMEASURE[i] = 100; + catalogs->catalog[i].Naverage = 0; + catalogs->catalog[i].Nmeasure = 0; + ALLOCATE (catalogs->catalog[i].average, Average, catalogs->NAVERAGE[i]); + ALLOCATE (catalogs->catalog[i].measureT, MeasureTiny, catalogs->NMEASURE[i]); + ALLOCATE (catalogs->catalog[i].secfilt, SecFilt, catalogs->NAVERAGE[i]*Nsecfilt); + } + return catalogs; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplitFree (CatalogSplitter *catalogs) { + + free (catalogs->catIDs); + free (catalogs->NAVERAGE); + free (catalogs->NMEASURE); + free (catalogs->index); + free (catalogs); + return TRUE; +} + +// distribute a bright catalog across separate catalogs +int BrightCatalogSplit (CatalogSplitter *catalogs, BrightCatalog *bcatalog) { + + int i; + + int Nsecfilt = catalogs->Nsecfilt; + + // find the max value of catID in this BrightCatalog + int catIDmax = 0; + for (i = 0; i < bcatalog->Naverage; i++) { + catIDmax = MAX(catIDmax, bcatalog->average[i].catID); + } + // XXX validate the measure ID range here + + int maxIDold = catalogs->maxID; + catalogs->maxID = MAX (maxIDold, catIDmax); + + // extend the index array and init + REALLOCATE (catalogs->index, int, catalogs->maxID + 1); + for (i = maxIDold + 1; i <= catalogs->maxID; i++) catalogs->index[i] = -1; + + // identify the new catID values + for (i = 0; i < bcatalog->Naverage; i++) { + int catID = bcatalog->average[i].catID; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + int idx = catalogs->index[catID]; + if (idx != -1) continue; // already have seen this one + + // a new catID value + int Ncat = catalogs->Ncatalog; // the next available slot + catalogs->catIDs[Ncat] = catID; + assert (Ncat >= 0); + assert (Ncat < catalogs->NCATALOG); + + catalogs->catalog[Ncat].Nsecfilt = Nsecfilt; + catalogs->catalog[Ncat].catID = catID; + + catalogs->index[catID] = Ncat; + assert (catID > 0); + assert (catID < catalogs->maxID + 1); + + catalogs->Ncatalog ++; + + if (catalogs->Ncatalog >= catalogs->NCATALOG) { + catalogs->NCATALOG += 16; + // fprintf (stderr, "realloc catalogs->catalog: old: %llx ", (long long) catalogs->catalog); + REALLOCATE (catalogs->catalog, Catalog, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->catalog, (long long) (catalogs->catalog + sizeof(Catalog)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NAVERAGE: old: %llx ", (long long) catalogs->NAVERAGE); + REALLOCATE (catalogs->NAVERAGE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NAVERAGE, (long long) (catalogs->NAVERAGE + sizeof(off_t)*catalogs->NCATALOG)); + + // fprintf (stderr, "realloc catalogs->NMEASURE: old: %llx ", (long long) catalogs->NMEASURE); + REALLOCATE (catalogs->NMEASURE, off_t, catalogs->NCATALOG); + // fprintf (stderr, "new: %llx - %llx\n", (long long) catalogs->NMEASURE, (long long) (catalogs->NMEASURE + sizeof(off_t)*catalogs->NCATALOG)); + + REALLOCATE (catalogs->catIDs, int, catalogs->NCATALOG); + + int j; + for (j = catalogs->NCATALOG - 16; j < catalogs->NCATALOG; j++) { + dvo_catalog_init (&catalogs->catalog[j], TRUE); + catalogs->catIDs[j] = 0; + catalogs->NAVERAGE[j] = 100; + catalogs->NMEASURE[j] = 100; + catalogs->catalog[j].Naverage = 0; + catalogs->catalog[j].Nmeasure = 0; + ALLOCATE (catalogs->catalog[j].average, Average, catalogs->NAVERAGE[j]); + ALLOCATE (catalogs->catalog[j].measureT, MeasureTiny, catalogs->NMEASURE[j]); + ALLOCATE (catalogs->catalog[j].secfilt, SecFilt, catalogs->NAVERAGE[j]*Nsecfilt); + } + } + } + + // each bcatalog (from each host) has a different set of catID ranges + // they also (probably) have contiguous ranges. But, it is a bit tricky to be sure I + // can use that info, so perhaps ignore it + + // assign the averages to the corresponding catalog + for (i = 0; i < bcatalog->Naverage; i++) { + int ID = bcatalog->average[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + + int Na = catalogs->catalog[Nc].Naverage; + catalogs->catalog[Nc].average[Na] = bcatalog->average[i]; + + // secfilt entries are grouped in blocks for each average + int k; + for (k = 0; k < Nsecfilt; k++) { + catalogs->catalog[Nc].secfilt[Na*Nsecfilt + k] = bcatalog->secfilt[i*Nsecfilt + k]; + } + + catalogs->catalog[Nc].Naverage ++; + if (catalogs->catalog[Nc].Naverage >= catalogs->NAVERAGE[Nc]) { + catalogs->NAVERAGE[Nc] += 100; + REALLOCATE (catalogs->catalog[Nc].average, Average, catalogs->NAVERAGE[Nc]); + REALLOCATE (catalogs->catalog[Nc].secfilt, SecFilt, catalogs->NAVERAGE[Nc]*Nsecfilt); + } + } + + // assign the measures to the corresponding catalog + // XXX what about averef and related links? Do I need them? + for (i = 0; i < bcatalog->Nmeasure; i++) { + int ID = bcatalog->measure[i].catID; + int Nc = catalogs->index[ID]; + assert (Nc > -1); + assert (Nc < catalogs->NCATALOG); + int Na = catalogs->catalog[Nc].Nmeasure; + catalogs->catalog[Nc].measureT[Na] = bcatalog->measure[i]; + catalogs->catalog[Nc].Nmeasure ++; + if (catalogs->catalog[Nc].Nmeasure >= catalogs->NMEASURE[Nc]) { + catalogs->NMEASURE[Nc] += 100; + REALLOCATE (catalogs->catalog[Nc].measureT, MeasureTiny, catalogs->NMEASURE[Nc]); + } + } + + return TRUE; +} Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 33652) @@ -26,14 +26,9 @@ 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); GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); + + GetConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN); GetConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP); @@ -42,6 +37,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: /trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitChip.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/FitChip.c (revision 33652) @@ -4,5 +4,5 @@ int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) { - int i, NstatFull, Nstat, Niter, skip; + int i, NstatFull, Nstat, Niter; float dLsig, dMsig, dRsig; float dLsigFull, dMsigFull, dRsigFull; @@ -27,9 +27,9 @@ if (raw[i].mask) continue; if (isnan(raw[i].dMag)) { - raw[i].mask |= 0x0004; + raw[i].mask |= MARK_NAN_MAG_ERROR; continue; } if (raw[i].dMag > SIGMA_LIM) { - raw[i].mask |= 0x0008; + raw[i].mask |= MARK_BIG_MAG_ERROR; continue; } // is this redundant with ImageOps.c:915? @@ -39,5 +39,5 @@ dR = hypot (dL, dM); if (dR > dRmax) { - raw[i].mask |= 0x0010; + raw[i].mask |= MARK_BIG_OFFSET; continue; } @@ -52,6 +52,15 @@ // fit the requested order polynomial if (CHIPORDER > 0) { - image[0].coords.Npolyterms = CHIPORDER; - } + int Ndof_min = 3; + int order_max = 0.5*(sqrt(4*Nmatch - 4*Ndof_min + 1) - 3); + int order_use = MIN (CHIPORDER, order_max); + if (order_use < 1) { + if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nmatch); + image[0].flags |= ID_IMAGE_ASTROM_FEW; + return FALSE; + } + image[0].coords.Npolyterms = order_use; + } + if (fit) fit_free (fit); fit = fit_init (image[0].coords.Npolyterms); @@ -63,4 +72,5 @@ } +# if (0) // check if the fit has enough data points for the polynomial order skip = FALSE; @@ -86,4 +96,5 @@ return FALSE; } +# endif // measure the fit, update the coords & object coordinates @@ -91,4 +102,5 @@ fprintf (stderr, "failed to fit new model\n"); image[0].flags |= ID_IMAGE_ASTROM_FAIL; + if (fit) fit_free (fit); return FALSE; } @@ -97,4 +109,5 @@ fprintf (stderr, "failed to fit new model\n"); image[0].flags |= ID_IMAGE_ASTROM_FAIL; + if (fit) fit_free (fit); return FALSE; } @@ -114,13 +127,13 @@ for (i = 0; i < Nmatch; i++) { if (!raw[i].mask) continue; - if (raw[i].mask & 0x01) nMask1 ++; - if (raw[i].mask & 0x02) nMask2 ++; - if (raw[i].mask & 0x04) nMask3 ++; - if (raw[i].mask & 0x08) nMask4 ++; - if (raw[i].mask & 0x10) nMask5 ++; + if (raw[i].mask & MARK_TOO_FEW_MEAS ) nMask1 ++; + if (raw[i].mask & MARK_NAN_POS_ERROR) nMask2 ++; + if (raw[i].mask & MARK_NAN_MAG_ERROR) nMask3 ++; + if (raw[i].mask & MARK_BIG_MAG_ERROR) nMask4 ++; + if (raw[i].mask & MARK_BIG_OFFSET ) nMask5 ++; } GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM); - GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM); + GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM); int Nm = 0; @@ -133,5 +146,5 @@ image[0].nLinkAstrom = (Nm / Ns); - if (VERBOSE) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5); + if (VERBOSE2) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5); image[0].dXpixSys = dLsig; Index: /trunk/Ohana/src/relastro/src/FixProblemImages.c =================================================================== --- /trunk/Ohana/src/relastro/src/FixProblemImages.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/FixProblemImages.c (revision 33652) @@ -5,4 +5,7 @@ int FixProblemImages (SkyList *skylist) { + + fprintf (stderr, "this function is currently not used : check on change from measure -> measureT in load_catalogs\n"); + abort(); int Nbad; @@ -64,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: /trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- /trunk/Ohana/src/relastro/src/GetAstromError.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/GetAstromError.c (revision 33652) @@ -1,4 +1,46 @@ # 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); + + // do not raise an exception, just send back the result + if (isnan(code[0].astromErrSys)) return NAN; + + 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) { @@ -30,4 +72,8 @@ code = GetPhotcodebyCode (measure[0].photcode); + + // do not raise an exception, just send back the result + if (isnan(code[0].astromErrSys)) return NAN; + AS = code[0].astromErrScale; MS = code[0].astromErrMagScale; Index: /trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relastro/src/ImageOps.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/ImageOps.c (revision 33652) @@ -3,17 +3,22 @@ # define USE_IMAGE_ID 1 -static off_t **bin; // link from catalog,measure to image -static int **clist; // catalog which supplied measurement on image -static off_t **mlist; // measure reference for measurement on image -static off_t *Nlist; // number of measurements on image -static off_t *NLIST; // allocated number of measurements on image - +// we have an array of images: image[ImageIndex] (ImageIndex : 0 < Nimage) static Image *image; // list of available images static off_t Nimage; // number of available images + +// if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row +// each image belongs static off_t *LineNumber; // match of subset to full image table -static int *Ncatlist; // catalogs associated with each image -static int *NCATLIST; // catalogs associated with each image +static off_t *N_onImage; // number of measurements on image +static off_t *N_ONIMAGE; // allocated number of measurements on image + +static int *Ncatlist; // catalogs associated with each image +static int *NCATLIST; // catalogs associated with each image static int **catlist; // catalogs associated with each image + +static IDX_T **MeasureToImage; // link from catalog,measure to image +static IDX_T **ImageToCatalog; // catalog which supplied measurement on image +static IDX_T **ImageToMeasure; // measure reference for measurement on image // if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted @@ -28,4 +33,11 @@ # endif +// MeasureToImage was 'bin' +// ImageToCatalog was 'clist' +// ImageToMeasure was 'mlist' + +// N_onImage was 'Nlist' +// N_ONIMAGE was 'NLIST' + Image *getimages (off_t *N, off_t **line_number) { @@ -105,20 +117,20 @@ off_t i, j; - ALLOCATE (bin, off_t *, Ncatalog); + ALLOCATE (MeasureToImage, IDX_T *, Ncatalog); for (i = 0; i < Ncatalog; i++) { - ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1)); - for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; - } - - ALLOCATE (Nlist, off_t, Nimage); - ALLOCATE (NLIST, off_t, Nimage); - ALLOCATE (clist, int *, Nimage); - ALLOCATE (mlist, off_t *, Nimage); + ALLOCATE (MeasureToImage[i], IDX_T, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToImage[i][j] = -1; + } + + ALLOCATE (N_onImage, off_t, Nimage); + ALLOCATE (N_ONIMAGE, off_t, Nimage); + ALLOCATE (ImageToCatalog, IDX_T *, Nimage); + ALLOCATE (ImageToMeasure, IDX_T *, Nimage); for (i = 0; i < Nimage; i++) { - Nlist[i] = 0; - NLIST[i] = 100; - clist[i] = NULL; // we allocate these iff they are needed in matchImage - mlist[i] = NULL; // we allocate these iff they are needed in matchImage + N_onImage[i] = 0; + N_ONIMAGE[i] = 30; + ImageToCatalog[i] = NULL; // we allocate these iff they are needed in matchImage + ImageToMeasure[i] = NULL; // we allocate these iff they are needed in matchImage } @@ -141,15 +153,15 @@ for (i = 0; i < Ncatalog; i++) { - free (bin[i]); - } - free (bin); + free (MeasureToImage[i]); + } + free (MeasureToImage); for (i = 0; i < Nimage; i++) { - if (clist[i]) { free (clist[i]); } - if (mlist[i]) { free (mlist[i]); } - } - free (clist); - free (mlist); - free (Nlist); - free (NLIST); + if (ImageToCatalog[i]) { free (ImageToCatalog[i]); } + if (ImageToMeasure[i]) { free (ImageToMeasure[i]); } + } + free (ImageToCatalog); + free (ImageToMeasure); + free (N_onImage); + free (N_ONIMAGE); } @@ -169,9 +181,26 @@ } - for (i = 0; VERBOSE2 && (i < Nimage); i++) { - name = GetPhotcodeNamebyCode (image[i].photcode); - fprintf (stderr, "image "OFF_T_FMT" has "OFF_T_FMT" measures (%s, %s)\n", i, Nlist[i], - ohana_sec_to_date(image[i].tzero), name); - } + // watch for chips with few stars + int Nfew = 0; + int Nbad = 0; + for (i = 0; i < Nimage; i++) { + // ignore the PHU entries + if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) { + name = GetPhotcodeNamebyCode (image[i].photcode); + fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)", i, image[i].imageID, image[i].name, N_onImage[i], image[i].nstar, + ohana_sec_to_date(image[i].tzero), name); + if (N_onImage[i] < 20) { + fprintf (stderr, "*"); + Nfew ++; + } + if (N_onImage[i] < 15) { + fprintf (stderr, "*** warning ***"); + Nbad ++; + } + fprintf (stderr, "\n"); + } + } + fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad); } @@ -181,8 +210,8 @@ off_t idx, ID; - Measure *measure; + MeasureTiny *measure; int i, found; - measure = &catalog[cat].measure[meas]; + measure = &catalog[cat].measureT[meas]; ID = measure[0].imageID; @@ -194,25 +223,25 @@ // index for (catalog, measure) -> image - bin[cat][meas] = idx; + MeasureToImage[cat][meas] = idx; // if we need to allocate an image index table, do so here - if (!clist[idx]) { - ALLOCATE (clist[idx], int, NLIST[idx]); - } - if (!mlist[idx]) { - ALLOCATE (mlist[idx], off_t, NLIST[idx]); + if (!ImageToCatalog[idx]) { + ALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]); + } + if (!ImageToMeasure[idx]) { + ALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]); } // index for image, Nentry -> catalog - clist[idx][Nlist[idx]] = cat; + ImageToCatalog[idx][N_onImage[idx]] = cat; // index for image, Nentry -> measure - mlist[idx][Nlist[idx]] = meas; - Nlist[idx] ++; - - if (Nlist[idx] == NLIST[idx]) { - NLIST[idx] += 100; - REALLOCATE (clist[idx], int, NLIST[idx]); - REALLOCATE (mlist[idx], off_t, NLIST[idx]); + ImageToMeasure[idx][N_onImage[idx]] = meas; + N_onImage[idx] ++; + + if (N_onImage[idx] == N_ONIMAGE[idx]) { + N_ONIMAGE[idx] += 30; + REALLOCATE (ImageToCatalog[idx], IDX_T, N_ONIMAGE[idx]); + REALLOCATE (ImageToMeasure[idx], IDX_T, N_ONIMAGE[idx]); } @@ -241,7 +270,7 @@ off_t i; - Measure *measure; - - measure = &catalog[cat].measure[meas]; + MeasureTiny *measure; + + measure = &catalog[cat].measureT[meas]; /* find the image that supplied this measurement */ @@ -255,17 +284,17 @@ // index for (catalog, measure) -> image - bin[cat][meas] = i; + MeasureToImage[cat][meas] = i; // index for image, Nentry -> catalog - clist[i][Nlist[i]] = cat; + ImageToCatalog[i][N_onImage[i]] = cat; // index for image, Nentry -> measure - mlist[i][Nlist[i]] = meas; - Nlist[i] ++; - - if (Nlist[i] == NLIST[i]) { - NLIST[i] += 100; - REALLOCATE (clist[i], int, NLIST[i]); - REALLOCATE (mlist[i], off_t, NLIST[i]); + ImageToMeasure[i][N_onImage[i]] = meas; + N_onImage[i] ++; + + if (N_onImage[i] == N_ONIMAGE[i]) { + N_ONIMAGE[i] += 30; + REALLOCATE (ImageToCatalog[i], IDX_T, N_ONIMAGE[i]); + REALLOCATE (ImageToMeasure[i], IDX_T, N_ONIMAGE[i]); } return; @@ -281,5 +310,5 @@ off_t i; - i = bin[cat][meas]; + i = MeasureToImage[cat][meas]; if (i == -1) return (NULL); return (&image[i].coords); @@ -333,5 +362,5 @@ } -void dump_measures(Average *average, Measure *measure) { +void dump_measures (Average *average, Measure *measure) { off_t j, off; @@ -402,11 +431,16 @@ } - for (i = 0; i < Nlist[im]; i++) { - m = mlist[im][i]; - c = clist[im][i]; - - X = catalog[c].measure[m].Xccd; - Y = catalog[c].measure[m].Yccd; - n = catalog[c].measure[m].averef; + int NoffRAave = 0; int NoffRAori = 0; + int NoffDECave = 0; int NoffDECori = 0; + + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + Measure *measure = &catalog[c].measure[m]; + + X = measure[0].Xccd; + Y = measure[0].Yccd; + n = measure[0].averef; if (moscoords == NULL) { @@ -428,47 +462,60 @@ // NOTE: This should never happen, or our StarMap tests are not working 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].measure); + NoffRAave ++; + if (VERBOSE2) { + 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].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].measure); + NoffDECave ++; + if (VERBOSE2) { + 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].measure); + } // abort (); } // complain if the new location is far from the old location - if (fabs(catalog[c].measure[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].measure); + if (fabs(catalog[c].measureT[m].dR - dR) > DPOS_MAX_ASEC) { + NoffRAori ++; + if (VERBOSE2) { + 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].measure); + } // abort(); } - if (fabs(catalog[c].measure[m].dD - dD) > DPOS_MAX_ASEC) { - fprintf (stderr, "measurement is far from original 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].measure); - // abort(); - } - - dPos += SQ(catalog[c].measure[m].dR - dR) + SQ(catalog[c].measure[m].dD - dD); + if (fabs(measure[0].dD - dD) > DPOS_MAX_ASEC) { + NoffDECori ++; + if (VERBOSE2) { + fprintf (stderr, "measurement is far from original 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].measure); + } // abort(); + } + + dPos += SQ(measure[0].dR - dR) + SQ(measure[0].dD - dD); nPos ++; - catalog[c].measure[m].dR = dR; - catalog[c].measure[m].dD = dD; + measure[0].dR = dR; + measure[0].dD = dD; - if (catalog[c].measure[m].dR > +180.0*3600.0) { + if (measure[0].dR > +180.0*3600.0) { // average on high end of boundary, move star up R += 360.0; - catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); - } - if (catalog[c].measure[m].dR < -180.0*3600.0) { + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (measure[0].dR < -180.0*3600.0) { // average on low end of boundary, move star down R -= 360.0; - catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); } // set the systematic error for this image: - catalog[c].measure[m].dRsys = ToShortPixels(dPosSys); - } - + measure[0].dRsys = ToShortPixels(dPosSys); + } + + int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori; + if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); saveOffsets (dPos, nPos, im); @@ -487,4 +534,7 @@ Mosaic *mosaic; Coords *moscoords, *imcoords, *oldcoords; + + fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n"); + abort(); // check if this image is bad and should be skipped @@ -507,11 +557,14 @@ imcoords = &image[im].coords; - for (i = 0; i < Nlist[im]; i++) { - m = mlist[im][i]; - c = clist[im][i]; - - X = catalog[c].measure[m].Xccd; - Y = catalog[c].measure[m].Yccd; - n = catalog[c].measure[m].averef; + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated + Measure *measure = &catalog[c].measure[m]; + + X = measure[0].Xccd; + Y = measure[0].Yccd; + n = measure[0].averef; dR = dD = 0.0; @@ -527,16 +580,16 @@ } - catalog[c].measure[m].dR = dR; - catalog[c].measure[m].dD = dD; - - if (catalog[c].measure[m].dR > +180.0*3600.0) { + measure[0].dR = dR; + measure[0].dD = dD; + + if (measure[0].dR > +180.0*3600.0) { // average on high end of boundary, move star up R += 360.0; - catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); - } - if (catalog[c].measure[m].dR < -180.0*3600.0) { + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); + } + if (measure[0].dR < -180.0*3600.0) { // average on low end of boundary, move star down R -= 360.0; - catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R); + measure[0].dR = 3600.0*(catalog[c].average[n].R - R); } } @@ -567,28 +620,30 @@ } - ALLOCATE (raw, StarData, Nlist[im]); - - for (i = 0; i < Nlist[im]; i++) { - m = mlist[im][i]; - c = clist[im][i]; + ALLOCATE (raw, StarData, N_onImage[im]); + + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + MeasureTiny *measure = &catalog[c].measureT[m]; /* apply the current image transformation or use the current value of R+dR, D+dD? */ - raw[i].X = catalog[c].measure[m].Xccd; - raw[i].Y = catalog[c].measure[m].Yccd; - - raw[i].Mag = catalog[c].measure[m].M; - raw[i].dMag = catalog[c].measure[m].dM; - raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS); - - n = catalog[c].measure[m].averef; + raw[i].X = measure[0].Xccd; + raw[i].Y = measure[0].Yccd; + + raw[i].Mag = measure[0].M; + raw[i].dMag = measure[0].dM; + raw[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); + + n = measure[0].averef; // an object with only one detection provides no information about the image calibration // XXX this is already taken care of in bcatalog - raw[i].mask = 0x0000; + raw[i].mask = MARK_MEAS_DEFAULT; if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) { - raw[i].mask |= 0x0001; - } - if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) { - raw[i].mask |= 0x0002; + raw[i].mask |= MARK_TOO_FEW_MEAS; + } + if (!finite(measure[0].dR) || !finite(measure[0].dD)) { + raw[i].mask |= MARK_NAN_POS_ERROR; } raw[i].Nmeas = catalog[c].average[n].Nmeasure; // record so we can check how well connected an image is @@ -613,5 +668,5 @@ } - *Nstars = Nlist[im]; + *Nstars = N_onImage[im]; return (raw); } @@ -639,10 +694,12 @@ } - ALLOCATE (ref, StarData, Nlist[im]); - - for (i = 0; i < Nlist[im]; i++) { - m = mlist[im][i]; - c = clist[im][i]; - n = catalog[c].measure[m].averef; + ALLOCATE (ref, StarData, N_onImage[im]); + + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + 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? */ @@ -650,7 +707,7 @@ ref[i].D = catalog[c].average[n].D; - ref[i].Mag = catalog[c].measure[m].M; - ref[i].dMag = catalog[c].measure[m].dM; - ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS); + ref[i].Mag = measure[0].M; + ref[i].dMag = measure[0].dM; + ref[i].dPos = GetAstromErrorTiny (&measure[0], ERROR_MODE_POS); ref[i].mask = FALSE; @@ -675,5 +732,5 @@ } - *Nstars = Nlist[im]; + *Nstars = N_onImage[im]; return (ref); } @@ -682,6 +739,8 @@ void FlagOutliers2D(Catalog *catalog); +// operates on Full values (not tiny) void FlagOutliers (Catalog *catalog) { - // XXX FlagOutliers is just using FlagOutliers2D + + // XXX FlagOutliers is now just using FlagOutliers2D FlagOutliers2D(catalog); return; @@ -696,5 +755,5 @@ assert(catalog[0].Nsecfilt == Nsecfilt); - if (VERBOSE) fprintf (stderr, "marking poor measures\n"); + if (VERBOSE2) fprintf (stderr, "marking poor measures\n"); Nmax = 0; for (i = 0; i < catalog[0].Naverage; i++) { @@ -716,15 +775,17 @@ for (j = 0; j < catalog[0].Naverage; j++) { + // pointer to this set of measurements + m = catalog[0].average[j].measureOffset; + Measure *measure = &catalog[0].measure[m]; + /* accumulate list of valid measurements */ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - R[N] = catalog[0].measure[m].dR; - D[N] = catalog[0].measure[m].dD; - dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); - dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (!MeasFilterTest(&measure[k], FALSE)) continue; + R[N] = measure[k].dR; + D[N] = measure[k].dD; + dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA); + dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC); if (isnan(R[N]) || isnan(D[N])) continue; N++; @@ -741,19 +802,18 @@ /* compare per-object distance to this standard deviation, and flag outliers*/ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { //reset flag on each invocation - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; + if (!MeasFilterTest(&measure[k], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta)))) { - catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; Ndel++; } @@ -787,5 +847,5 @@ assert(catalog[0].Nsecfilt == Nsecfilt); - if (VERBOSE) fprintf (stderr, "marking poor measures\n"); + if (VERBOSE2) fprintf (stderr, "marking poor measures\n"); Nmax = 0; for (i = 0; i < catalog[0].Naverage; i++) { @@ -809,19 +869,21 @@ for (j = 0; j < catalog[0].Naverage; j++) { + // pointer to this set of measurements + m = catalog[0].average[j].measureOffset; + Measure *measure = &catalog[0].measure[m]; + /* accumulate list of valid measurements */ - m = catalog[0].average[j].measureOffset; N = 0; - - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - - //reset flag on each invocation - catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { + + // reset flag on each invocation + measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; // skip measurements based on user selected criteria - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - R[N] = catalog[0].measure[m].dR; - D[N] = catalog[0].measure[m].dD; - dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA); - dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC); + if (!MeasFilterTest(&measure[k], FALSE)) continue; + R[N] = measure[k].dR; + D[N] = measure[k].dD; + dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA); + dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC); if (isnan(R[N]) || isnan(D[N])) continue; N++; @@ -837,15 +899,14 @@ /* calculate deviations of all points*/ - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { - //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { + // skip bad measurements + if (!MeasFilterTest(&measure[k], FALSE)) continue; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta))); - index[N] = m; + index[N] = k; N++; } @@ -856,10 +917,10 @@ // recalculate image center, sigma based on closest 50% of points - for(k = 0; k < N; k++) { + for (k = 0; k < N; k++) { off_t ind = (off_t) index[k]; - R[k] = catalog[0].measure[ind].dR; - D[k] = catalog[0].measure[ind].dD; - dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA); - dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC); + R[k] = measure[ind].dR; + D[k] = measure[ind].dD; + dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA); + dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC); } liststats (R, dR, N, &statsR); @@ -869,26 +930,25 @@ // use these new statistics to flag outliers - m = catalog[0].average[j].measureOffset; N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) { + for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { //skip bad measurements - if (!MeasFilterTest(&catalog[0].measure[m], FALSE)) continue; - x = catalog[0].measure[m].dR - statsR.median; - y = catalog[0].measure[m].dD - statsD.median; + if (!MeasFilterTest(&measure[k], FALSE)) continue; + x = measure[k].dR - statsR.median; + y = measure[k].dD - statsD.median; theta = atan2(y,x); d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + SQR(statsD.sigma * Ns * sin(theta))); if ((d2[N]) > 1) { - catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM; + measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; Ndel ++; } N++; Nave++; - } //done rejecting outliers + } // done rejecting outliers // examine results // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); - } //done looping over objects + } // done looping over objects if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); @@ -901,5 +961,65 @@ } - +/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ +// we only optionally apply the sigma limit: for object averages, this should not be used (should it?) +int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) { + int found, k; + long mask; + PhotCode *code; + float mag; + + if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE; + if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all 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 = PhotInstTiny (measure); + if (mag < ImagMin || mag > ImagMax) return FALSE; + } + + return TRUE; +} /** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ Index: /trunk/Ohana/src/relastro/src/ParFactor.c =================================================================== --- /trunk/Ohana/src/relastro/src/ParFactor.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/ParFactor.c (revision 33652) @@ -1,3 +1,4 @@ # include "relastro.h" +# define J2000 2451545. /* Julian date at standard epoch */ # if (0) @@ -31,5 +32,4 @@ double n, L, g; -# define J2000 2451545. /* Julian date at standard epoch */ n = jd - J2000; @@ -46,6 +46,4 @@ int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius) { -# define J2000 2451545. /* Julian date at standard epoch */ - double n = jd - J2000; // day number double L = 280.460 + 0.9856474 * n; // mean solar longitute (corr. for aberration) @@ -60,4 +58,5 @@ /* given RA, DEC, Time, calculate the parallax factor */ +// Time is relative to Tmean, Tmean is years relative to J2000 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time, double Tmean) { @@ -66,5 +65,8 @@ /* given a time T in UNIX seconds, determine the solar longitude S */ - jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); + // jd = ohana_sec_to_jd (365.25*86400.0*(Time + Tmean)); + jd = 365.25*(Time + Tmean) + J2000; + // fprintf (stderr, "Time: %f, jd: %f\n", Time, jd); + sun_ecliptic (jd, &lambda, &beta, &epsilon, &Radius); Index: /trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- /trunk/Ohana/src/relastro/src/StarMaps.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/StarMaps.c (revision 33652) @@ -68,9 +68,11 @@ for (i = 0; i < catalog[0].Nmeasure; i++) { - N = getImageByID(catalog[0].measure[i].imageID); + MeasureTiny *measure = &catalog[0].measureT[i]; + + N = getImageByID(measure[0].imageID); if (N < 0) continue; - xbin = catalog[0].measure[i].Xccd / starmap[N].Nx; - ybin = catalog[0].measure[i].Yccd / starmap[N].Ny; + xbin = measure[0].Xccd / starmap[N].Nx; + ybin = measure[0].Yccd / starmap[N].Ny; xbin = MAX(0, MIN(NX_MAP-1, xbin)); @@ -82,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: /trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateChips.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/UpdateChips.c (revision 33652) @@ -1,13 +1,38 @@ # include "relastro.h" + int plotChipFits (double *Ro, double *Do, char *mode, int Nimage); int saveCenter (Image *image, double *Ro, double *Do, int im); - +off_t getNextImageForThread (); +int UpdateChips_threaded (Catalog *catalog, int Ncatalog); +void *UpdateChips_worker (void *data); + +enum {THREAD_RUN, THREAD_DONE}; + +typedef struct { + int entry; + int state; + double *Ro; + double *Do; + char *mode; + off_t Nskip; + off_t Nmosaic; + off_t NnewFit; + off_t NoldFit; + Catalog *catalog; + int Ncatalog; +} ThreadInfo; + +static Image *image = NULL; +static off_t Nimage = 0; +static off_t nextImage = 0; + +// update astrometry of all chips relative to the average positions +// if NTHREADS is non-zero, call the threaded version of this function int UpdateChips (Catalog *catalog, int Ncatalog) { - int Nskip, Nmosaic, NnewFit, NoldFit; + off_t Nskip, Nmosaic, NnewFit, NoldFit; /* we can measure new image parameters for each non-mosaic chip independently */ - off_t i, Nimage, Nraw, Nref, nFitAstr; - Image *image; + off_t i, Nraw, Nref, nFitAstr; StarData *raw, *ref; Coords *oldCoords; @@ -16,4 +41,9 @@ char *mode; + if (NTHREADS) { + UpdateChips_threaded (catalog, Ncatalog); + return TRUE; + } + Nskip = Nmosaic = NnewFit = NoldFit = 0; @@ -25,5 +55,10 @@ ALLOCATE (mode, char, Nimage); + // XXX for faster processing in the future, this can be easily run in parallel + // each chip is fitted independently, so we could do N at once in parallel for (i = 0; i < Nimage; i++) { + + // XXX looks like everything below is thread safe : we can unroll this into a set of + // helper functions that grab the next available chip.... /* skip all except WRP images */ @@ -47,4 +82,5 @@ Nskip ++; mode[i] = 0; + free (raw); continue; } @@ -62,5 +98,5 @@ // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n", i, Nraw); if (!FitChip (raw, ref, Nraw, &image[i])) { - if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); // restore status quo ante @@ -105,9 +141,216 @@ plotChipFits (Ro, Do, mode, Nimage); - - fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic); + free (Ro); + free (Do); + free (mode); + + fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", + (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic); return (TRUE); } +int UpdateChips_threaded (Catalog *catalog, int Ncatalog) { + + int i; + off_t Nskip, Nmosaic, NnewFit, NoldFit; + double *Ro, *Do; + char *mode; + + image = getimages (&Nimage, NULL); + nextImage = 0; + + // save fit results for summary plot + ALLOCATE (Ro, double, Nimage); + ALLOCATE (Do, double, Nimage); + ALLOCATE (mode, char, Nimage); + + pthread_attr_t attr; + pthread_attr_init (&attr); + pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED); + + pthread_t *threads; + ALLOCATE (threads, pthread_t, NTHREADS); + + ThreadInfo *threadinfo; + ALLOCATE (threadinfo, ThreadInfo, NTHREADS); + + // launch N worker threads + for (i = 0; i < NTHREADS; i++) { + threadinfo[i].entry = i; + threadinfo[i].state = THREAD_RUN; + threadinfo[i].Ro = Ro; + threadinfo[i].Do = Do; + threadinfo[i].mode = mode; + threadinfo[i].Nskip = 0; + threadinfo[i].Nmosaic = 0; + threadinfo[i].NnewFit = 0; + threadinfo[i].NoldFit = 0; + threadinfo[i].catalog = catalog; + threadinfo[i].Ncatalog = Ncatalog; + pthread_create (&threads[i], NULL, UpdateChips_worker, &threadinfo[i]); + } + pthread_attr_destroy (&attr); + + // wait until all threads have finished + while (1) { + int allDone = TRUE; + for (i = 0; i < NTHREADS; i++) { + if (threadinfo[i].state == THREAD_RUN) allDone = FALSE; + } + if (allDone) { + break; + } + usleep (500000); + } + + // all threads are done, free the threads array and grab the info + free (threads); + + plotChipFits (Ro, Do, mode, Nimage); + free (Ro); + free (Do); + free (mode); + + // report stats & summary from the threads + Nskip = Nmosaic = NnewFit = NoldFit = 0; + for (i = 0; i < NTHREADS; i++) { + fprintf (stderr, "UpdateChips thread %d : %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", + i, (int) threadinfo[i].NnewFit, (int) threadinfo[i].NoldFit, (int) threadinfo[i].Nskip, (int) threadinfo[i].Nmosaic); + + NnewFit += threadinfo[i].NnewFit; + NoldFit += threadinfo[i].NoldFit; + Nskip += threadinfo[i].Nskip; + Nmosaic += threadinfo[i].Nmosaic; + } + free (threadinfo); + + fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", + (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic); + return (TRUE); +} + +void *UpdateChips_worker (void *data) { + + /* we can measure new image parameters for each non-mosaic chip independently */ + off_t Nraw, Nref, nFitAstr; + StarData *raw, *ref; + Coords *oldCoords; + float dXpixSys, dYpixSys; + + ThreadInfo *threadinfo = data; + + while (1) { + + off_t i = getNextImageForThread(); + if (i == -1) { + threadinfo->state = THREAD_DONE; + return NULL; + } + + /* skip all except WRP images */ + if (strcmp(&image[i].coords.ctype[4], "-WRP")) { + threadinfo->Nmosaic ++; + threadinfo->mode[i] = 0; + continue; + } + + /* convert measure coordinates to raw entries */ + raw = getImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, &Nraw, MODE_MOSAIC); + if (!raw) { + threadinfo->Nskip ++; + threadinfo->mode[i] = 0; + continue; + } + + /* convert average coordinates to ref entries */ + ref = getImageRef (threadinfo->catalog, threadinfo->Ncatalog, i, &Nref, MODE_MOSAIC); + if (!ref) { + threadinfo->Nskip ++; + threadinfo->mode[i] = 0; + free (raw); + continue; + } + + // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions. + assert (Nraw == Nref); + + // save these in case of failure + saveCoords (&image[i].coords, i); + dXpixSys = image[i].dXpixSys; + dYpixSys = image[i].dYpixSys; + nFitAstr = image[i].nFitAstrom; + + // FitChip does iterative, clipped fitting + // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n", i, Nraw); + if (!FitChip (raw, ref, Nraw, &image[i])) { + if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); + + // restore status quo ante + oldCoords = getCoords (i); + memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i); + threadinfo->mode[i] = 1; + threadinfo->NoldFit ++; + free (raw); + free (ref); + continue; + } + + if (!checkStarMap (i)) { + if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); + // restore status quo ante + oldCoords = getCoords (i); + memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i); + threadinfo->mode[i] = 2; + image[i].flags |= ID_IMAGE_ASTROM_POOR; + threadinfo->NoldFit ++; + free (raw); + free (ref); + continue; + } + + saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i); + threadinfo->mode[i] = 3; + threadinfo->NnewFit ++; + free (raw); + free (ref); + } + + // we should never reach here... + return NULL; +} + +// mutex to lock UpdateChips_worker operations +static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER; + +// we have an array of chips (image, Nimage). we need to hand out images one at a time to +// the worker threads as they need +off_t getNextImageForThread () { + + pthread_mutex_lock (&UpdateChips_mutex); + if (nextImage >= Nimage) { + pthread_mutex_unlock (&UpdateChips_mutex); + return (-1); + } + off_t thisImage = nextImage; + nextImage ++; + + pthread_mutex_unlock (&UpdateChips_mutex); + return (thisImage); +} + +// This function uses the mosaic array in MosaicOps.c to associate PHU + CHIP this is +// thread safe : the lookup (getMosaicForImage) is selecting an element of a previously +// allocated and assigned static array. This function does NOT use the thread-unsafe (and +// somewhat slower) functions in libdvo/src/mosaic_astrom.c to associate CHIP and PHU int saveCenter (Image *image, double *Ro, double *Do, int im) { @@ -136,17 +379,12 @@ // note that for a Simple image, L,M = P,Q XY_to_LM (&L, &M, X, Y, imcoords); - LM_to_RD (&R, &D, L, M, imcoords); + LM_to_RD (&R, &D, L, M, imcoords); // because of the block above, ctype is not -WRP } else { XY_to_LM (&L, &M, X, Y, imcoords); XY_to_LM (&P, &Q, L, M, moscoords); - LM_to_RD (&R, &D, P, Q, moscoords); - } - - double Rmid; - if (UserCatalog) { - Rmid = UserCatalogRA; - } else { - Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax); - } + LM_to_RD (&R, &D, P, Q, moscoords); // moscoords.ctype is -DIS + } + + double Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax); R = ohana_normalize_angle_to_midpoint (R, Rmid); Index: /trunk/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 33652) @@ -1,4 +1,5 @@ # include "relastro.h" +// this function operates on Measure, not MeasureTiny int UpdateMeasures (Catalog *catalog, int Ncatalog) { Index: /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 33652) @@ -7,14 +7,29 @@ // 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++) { + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], 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 @@ -22,14 +37,17 @@ catalog.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); exit (1); } if (!catalog.Naves_disk) { - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); dvo_catalog_unlock (&catalog); dvo_catalog_free (&catalog); continue; } + + // set the values in MeasureTiny needed by UpdateObjects + populate_tiny_values(&catalog, DVO_TV_MEASURE); // match measurements with images @@ -40,4 +58,5 @@ UpdateMeasures (&catalog, 1); + // uses MeasureTiny values, but will also update Measure values when populated UpdateObjects (&catalog, 1); @@ -50,2 +69,97 @@ 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 (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 (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) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 33652) @@ -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) { @@ -72,11 +74,14 @@ T2000 = ohana_date_to_sec ("2000/01/01"); // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? + Nsecfilt = GetPhotcodeNsecfilt (); - assert (catalog[0].Nsecfilt == Nsecfilt); + if (Ncatalog) { + assert (catalog[0].Nsecfilt == Nsecfilt); + } NaveSum = NparSum = NpmSum = NskipSum = 0; for (i = 0; i < Ncatalog; i++) { - if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); + if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); Nave = Npar = Npm = Nskip = 0; @@ -86,13 +91,4 @@ XVERB = FALSE; - // skip objects which are known to be problematic - // XXX include this code or not? -# if (0) - if (catalog[i].average[j].code & STAR_BAD) { - Nskip ++; - continue; - } -# endif - if (catalog[i].average[j].Nmeasure == 0) { continue; @@ -101,50 +97,56 @@ N = 0; m = catalog[i].average[j].measureOffset; - Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25); + 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 + mode = FIT_MODE; // find the basic properties of the detections for this object (Tmin, Tmax, Tmean) - Tmean = 0; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - - //does the measurement pass the supplied filtering constraints? - if (!MeasFilterTest(&catalog[i].measure[m], FALSE)) { - catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + for (k = 0; k < catalog[i].average[j].Nmeasure; k++) { + + // does the measurement pass the supplied filtering constraints? + if (!MeasFilterTestTiny(&measure[k], FALSE)) { + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; } continue; } //outlier rejection - if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) { - catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ; + 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; } - // exclude measurements by previous outlier detection - // XXX include this code or not? -# if (0) - if (catalog[i].measure[m].dbFlags & MEAS_BAD) { - catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM; + 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]); + D[N] = getMeanD (&measure[k], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + + // XXX I think this is a problem: T[] is time in years relative to J2000, but ParFactor expects + // to get Time in years relative to UNIX Tzero (1970/01/01) + T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years + + // dX, dY : error in arcsec -- + dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA); + dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC); + + // allow a given photcode or measurement to be + // ignored if the error is NAN (for photcode, set astromErrSys to NaN) + if (isnan(dX[N])) { + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; continue; } -# endif - - catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ; - - R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - T[N] = (catalog[i].measure[m].t - T2000) / (86400*365.25) ; // time relative to J2000 in years - - Tmin = MIN(Tmin, T[N]); - Tmax = MAX(Tmax, T[N]); - Tmean += T[N]; - - // dX, dY : error in arcsec -- - dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA); - dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC); + if (isnan(dY[N])) { + measure[k].dbFlags &= ~ID_MEAS_USED_OBJ; + continue; + } // add systematic error in quadrature, if desired // only do this after the fit has converged (or you will never improve the poor images) // if (INCLUDE_SYS_ERR) { - // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys); + // float dRsys = FromShortPixels(measure[k].dRsys); // dX[N] = hypot(dX[N], dRsys); // dY[N] = hypot(dY[N], dRsys); @@ -154,5 +156,5 @@ // dY[N] = 0.1; - dT[N] = catalog[i].measure[m].dt; + dT[N] = measure[k].dt; // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in @@ -172,8 +174,22 @@ catalog[i].average[j].flags &= ~ID_STAR_FEW; + // find Tmin & Tmax from the list of accepted measurements + Tmean = 0; + Tmin = Tmax = T[0]; + for (k = 0; k < N; k++) { + Tmin = MIN(Tmin, T[k]); + Tmax = MAX(Tmax, T[k]); + Tmean += T[k]; + } // XXX add the parallax factor range as a criterion as well Trange = Tmax - Tmin; if (Trange < PM_DT_MIN) mode = FIT_AVERAGE; if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE; + + if (FIT_TARGET == TARGET_HIGH_SPEED) { + Tmean = 0.5*(Tmax - Tmin); + } else { + Tmean /= (float) N; + } // too few measurements for average position (require 2 values) @@ -192,14 +208,7 @@ coords.crval2 = D[0]; - if (FIT_TARGET == TARGET_HIGH_SPEED) { - Tmean = 0.5*(Tmax - Tmin); - } else { - Tmean /= (float) N; - } - // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e); XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC); XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST); - // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY); // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq @@ -232,11 +241,25 @@ // fprintf (stderr, "parallax fitting is still untested (%s, %d)\n", __FILE__, __LINE__); + float pXmin = +2.0; + float pXmax = -2.0; + float pYmin = +2.0; + float pYmax = -2.0; for (k = 0; k < N; k++) { ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean); - } - FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB); - XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); - catalog[i].average[j].flags |= ID_STAR_FIT_PAR; - Npar ++; + pXmin = MIN (pXmin, pX[k]); + pXmax = MAX (pXmax, pX[k]); + pYmin = MIN (pYmin, pY[k]); + pYmax = MAX (pYmax, pY[k]); + } + float dXRange = pXmax - pXmin; + float dYRange = pYmax - pYmin; + float parRange = hypot (dXRange, dYRange); + + if (parRange >= PAR_FACTOR_MIN) { + FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB); + XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords); + catalog[i].average[j].flags |= ID_STAR_FIT_PAR; + Npar ++; + } } @@ -307,8 +330,11 @@ // the measure fields must be updated before the average fields - m = catalog[i].average[j].measureOffset; - for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { - setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); - setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]); + 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: /trunk/Ohana/src/relastro/src/args.c =================================================================== --- /trunk/Ohana/src/relastro/src/args.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/args.c (revision 33652) @@ -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,10 @@ if ((N = get_argument (argc, argv, "-update-objects"))) { remove_argument (N, &argc, argv); - FIT_TARGET = TARGET_OBJECTS; + FIT_TARGET = TARGET_UPDATE_OBJECTS; + } + + if ((N = get_argument (argc, argv, "-update-offsets"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_UPDATE_OFFSETS; } @@ -92,6 +98,16 @@ RADIUS = atof(argv[N]); remove_argument (N, &argc, argv); - HIGH_SPEED_DIR = strcreate(argv[N]); - remove_argument (N, &argc, argv); + HIGH_SPEED_DIR = abspath(argv[N], MAX_PATH_LENGTH); + remove_argument (N, &argc, argv); + } + + PARALLEL_OUTPUT = FALSE; + if ((N = get_argument (argc, argv, "-parallel-output"))) { + remove_argument (N, &argc, argv); + PARALLEL_OUTPUT = TRUE; + if (FIT_TARGET != TARGET_HIGH_SPEED) { + fprintf (stderr, "-parallel-output only valid for -high-speed mode\n"); + exit (1); + } } @@ -118,14 +134,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 +150,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 +329,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,17 +361,15 @@ } - /* 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"))) { remove_argument (N, &argc, argv); NLOOP = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + NTHREADS = 0; + if ((N = get_argument (argc, argv, "-threads"))) { + remove_argument (N, &argc, argv); + NTHREADS = atof (argv[N]); remove_argument (N, &argc, argv); } @@ -356,9 +379,284 @@ } +int args_client (int argc, char **argv) { + + int N; + double trange; + char *endptr; + + /* 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); + } + + if ((N = get_argument (argc, argv, "-testobj1"))) { + if (N > argc - 3) usage (); + remove_argument (N, &argc, argv); + OBJ_ID_SRC = strtol(argv[N], &endptr, 0); + if (*endptr) usage (); + remove_argument (N, &argc, argv); + CAT_ID_SRC = strtol(argv[N], &endptr, 0); + if (*endptr) usage (); + remove_argument (N, &argc, argv); + } + + if ((N = get_argument (argc, argv, "-testobj2"))) { + if (N > argc - 3) usage (); + remove_argument (N, &argc, argv); + OBJ_ID_DST = strtol(argv[N], &endptr, 0); + if (*endptr) usage (); + remove_argument (N, &argc, argv); + CAT_ID_DST = strtol(argv[N], &endptr, 0); + if (*endptr) usage (); + 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; + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + // do not hold the lock on the Image table for relastro_client -update-offsets + if (FIT_TARGET == TARGET_UPDATE_OFFSETS) { + 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 +667,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 +685,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 +694,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: /trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /trunk/Ohana/src/relastro/src/bcatalog.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/bcatalog.c (revision 33652) @@ -1,3 +1,8 @@ # include "relastro.h" + +static int Nskip1 = 0; +static int Nskip2 = 0; +static unsigned int Tref = 1323353985; +static short Cref = 10355; int bcatalog (Catalog *subcatalog, Catalog *catalog) { @@ -12,10 +17,21 @@ /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ + memset(subcatalog, 0, sizeof(Catalog)); NAVERAGE = 50; NMEASURE = 1000; - ALLOCATE (subcatalog[0].average, Average, NAVERAGE); - ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); - ALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + ALLOCATE (subcatalog[0].average, Average, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); Nmeasure = Naverage = 0; + + DVOAverageFlags averageBits = + ID_STAR_FEW | + ID_STAR_POOR | + ID_STAR_FIT_AVE | + ID_STAR_FIT_PM | + ID_STAR_FIT_PAR | + ID_STAR_USE_AVE | + ID_STAR_USE_PM | + ID_STAR_USE_PAR; /* exclude stars not in range or with too few measurements */ @@ -23,4 +39,5 @@ if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) { continue; + // XXX can't (easily) count this for a given chip } @@ -33,9 +50,5 @@ if (RESET) { - // XXX reset the ra,dec coords? - // XXX define astrometry vs photometry average.dbFlags - // XXX put photometry flags in secfilt (ie, per average mag?) - subcatalog[0].average[Naverage].flags &= ~ID_STAR_FEW; - subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR; + subcatalog[0].average[Naverage].flags &= ~averageBits; } @@ -48,4 +61,7 @@ if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; + if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nskip1 ++; + } continue; } @@ -54,4 +70,7 @@ if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) { catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP; + if (FALSE && (abs(catalog[0].measure[offset].t - Tref) < 10) && (catalog[0].measure[offset].photcode == Cref)) { + Nskip2 ++; + } continue; } @@ -70,13 +89,12 @@ // allowed. - // CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); - - subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; - subcatalog[0].measure[Nmeasure].averef = Naverage; + 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; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; - subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } @@ -85,5 +103,5 @@ if (Nmeasure == NMEASURE) { NMEASURE += 1000; - REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); } } @@ -96,26 +114,50 @@ } } - REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); - REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1)); - REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt); + REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, MAX (Naverage, 1)*Nsecfilt); subcatalog[0].Naverage = Naverage; subcatalog[0].Nmeasure = Nmeasure; subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; - subcatalog[0].Nsecf_mem = Naverage * Nsecfilt; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + subcatalog[0].catID = catalog[0].catID; assert (Nsecfilt == catalog[0].Nsecfilt); // limit the total number of stars in the catalog if (MaxDensityUse) { - LimitDensityCatalog (subcatalog, catalog); - } - - if (VERBOSE) { - fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n", i, subcatalog[0].Naverage, subcatalog[0].Nmeasure); - } + LimitDensityCatalog_ByNmeasure (subcatalog, catalog); + } else { + if (VERBOSE2) { + char *basename = filebasename (catalog[0].filename); + fprintf (stderr, "subset of "OFF_T_FMT" ("OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" total) measures for catalog %s\n", + subcatalog[0].Naverage, catalog[0].Naverage, subcatalog[0].Nmeasure, catalog[0].Nmeasure, basename); + free (basename); + } + } return (TRUE); } +void bcatalog_show_skips () { + fprintf (stderr, "Nskip: %d, %d\n", Nskip1, Nskip2); +} + +// sort by decreasing Nmeasure (X) +void sort_by_Nmeasure (int *X, off_t *Y, off_t N) { + +# define SWAPFUNC(A,B){ int tmpI; off_t tmpT; \ + tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \ + tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \ +} +# define COMPARE(A,B)(X[A] > X[B]) + + OHANA_SORT (N, COMPARE, SWAPFUNC); + +# undef SWAPFUNC +# undef COMPARE + +} + /* this version does NOT use AverageTiny, MeasureTiny */ -int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) { +int LimitDensityCatalog_ByNmeasure (Catalog *subcatalog, Catalog *oldcatalog) { Catalog tmpcatalog; @@ -125,8 +167,8 @@ int Nsecfilt = GetPhotcodeNsecfilt (); - gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); - gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); - gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); - gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + gfits_scan (&oldcatalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&oldcatalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&oldcatalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&oldcatalog[0].header, "DEC1", "%lf", 1, &Dmax); if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); @@ -147,16 +189,14 @@ // select a random subset of Nmax stars from subcatalog using Fisher-Yates - // we are going to select Nmax entries by generating a random-sorted index list - off_t *index, tmp, i, j, ave; + // we are going to select Nmax entries by choosing the brightest objects + int *value; + off_t *index, i, j, ave; ALLOCATE (index, off_t, Naverage); + ALLOCATE (value, int, Naverage); for (i = 0; i < Naverage; i++) { index[i] = i; - } - for (i = 0; i < Naverage; i++) { - j = (Naverage - i) * drand48() + i; // a number between i and Naverage - tmp = index[j]; - index[j] = index[i]; - index[i] = tmp; - } + value[i] = subcatalog[0].average[i].Nmeasure; + } + sort_by_Nmeasure (value, index, Naverage); // count the number of measurements this selection will yield @@ -168,7 +208,7 @@ // allocate the output data - ALLOCATE (tmpcatalog.average, Average, Nmax); - ALLOCATE (tmpcatalog.measure, Measure, NMEASURE); - ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); off_t Nmeasure = 0; @@ -181,11 +221,112 @@ for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { off_t offset = subcatalog[0].average[ave].measureOffset + j; - tmpcatalog.measure[Nmeasure] = subcatalog[0].measure[offset]; - tmpcatalog.measure[Nmeasure].averef = i; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; Nmeasure ++; } - } - - if (VERBOSE) { + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { + char *basename = filebasename (oldcatalog[0].filename); + fprintf (stderr, "limited to "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) stars, "OFF_T_FMT" ("OFF_T_FMT" subset, "OFF_T_FMT" total) measures for catalog %s\n", + Nmax, subcatalog[0].Naverage, oldcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, oldcatalog[0].Nmeasure, basename); + free (basename); + } + + free (index); + free (value); + free (subcatalog[0].average); + free (subcatalog[0].measureT); + free (subcatalog[0].secfilt); + + subcatalog[0].average = tmpcatalog.average; + subcatalog[0].measureT = tmpcatalog.measureT; + subcatalog[0].secfilt = tmpcatalog.secfilt; + subcatalog[0].Naverage = Nmax; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = oldcatalog[0].Nsecfilt; + subcatalog[0].Nsecf_mem = Naverage * oldcatalog[0].Nsecfilt; + + return (TRUE); +} + +/* this version does NOT use AverageTiny, MeasureTiny */ +int LimitDensityCatalog_RandomSample (Catalog *subcatalog, Catalog *catalog) { + + Catalog tmpcatalog; + + double Rmin, Rmax, Dmin, Dmax; + + int Nsecfilt = GetPhotcodeNsecfilt (); + + gfits_scan (&catalog[0].header, "RA0", "%lf", 1, &Rmin); + gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin); + gfits_scan (&catalog[0].header, "RA1", "%lf", 1, &Rmax); + gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax); + + if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax); + + float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin)); + assert (AREA > 0); + + off_t Nmax = MaxDensityValue * AREA; + if (subcatalog[0].Naverage <= Nmax) { + if (VERBOSE2) { + fprintf (stderr, "subcatalog has less than the max density\n"); + } + return (TRUE); + } + + off_t Naverage = subcatalog[0].Naverage; + + // select a random subset of Nmax stars from subcatalog using Fisher-Yates + + // we are going to select Nmax entries by generating a random-sorted index list + off_t *index, tmp, i, j, ave; + ALLOCATE (index, off_t, Naverage); + for (i = 0; i < Naverage; i++) { + index[i] = i; + } + for (i = 0; i < Naverage; i++) { + j = (Naverage - i) * drand48() + i; // a number between i and Naverage + tmp = index[j]; + index[j] = index[i]; + index[i] = tmp; + } + + // count the number of measurements this selection will yield + off_t NMEASURE = 0; + for (i = 0; i < Nmax; i++) { + ave = index[i]; + NMEASURE += subcatalog[0].average[ave].Nmeasure; + } + + // allocate the output data + ALLOCATE (tmpcatalog.average, Average, Nmax); + ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE); + ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt); + + off_t Nmeasure = 0; + + // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links) + for (i = 0; i < Nmax; i++) { + ave = index[i]; + tmpcatalog.average[i] = subcatalog[0].average[ave]; + tmpcatalog.average[i].measureOffset = Nmeasure; + for (j = 0; j < tmpcatalog.average[i].Nmeasure; j++) { + off_t offset = subcatalog[0].average[ave].measureOffset + j; + tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset]; + tmpcatalog.measureT[Nmeasure].averef = i; + Nmeasure ++; + } + for (j = 0; j < Nsecfilt; j++) { + tmpcatalog.secfilt[i*Nsecfilt + j] = subcatalog[0].secfilt[ave*Nsecfilt + j]; + } + } + + if (VERBOSE2) { fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure, catalog[0].filename); @@ -193,9 +334,9 @@ free (subcatalog[0].average); - free (subcatalog[0].measure); + free (subcatalog[0].measureT); free (subcatalog[0].secfilt); subcatalog[0].average = tmpcatalog.average; - subcatalog[0].measure = tmpcatalog.measure; + subcatalog[0].measureT = tmpcatalog.measureT; subcatalog[0].secfilt = tmpcatalog.secfilt; subcatalog[0].Naverage = Nmax; Index: /trunk/Ohana/src/relastro/src/dvo_astrom_ops.c =================================================================== --- /trunk/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 33652) @@ -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: /trunk/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitpoly.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/fitpoly.c (revision 33652) @@ -149,4 +149,6 @@ if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) { + array_free (matrix, fit[0].Nelems); + array_free (vector, fit[0].Nelems); return (FALSE); } Index: /trunk/Ohana/src/relastro/src/high_speed_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/high_speed_catalogs.c (revision 33652) @@ -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,24 +11,22 @@ 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_parallel (skylist); + goto finish; + } // load data from each region file, only use bright stars for (i = 0; i < skylist[0].Nregions; i++) { + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + // set up the basic catalog info - catalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : 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 @@ -37,10 +34,10 @@ catalog.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) { + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "r")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); exit (1); } if (!catalog.Naves_disk) { - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); dvo_catalog_unlock (&catalog); dvo_catalog_free (&catalog); @@ -54,4 +51,152 @@ } +finish: + + // in relastro (NOT relastr_client), save the Images.dat, Photcodes.dat, SkyTable.fits files + if (!hostID) { + 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); } + +// CATDIR is supplied globally +# define DEBUG 1 +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, skylist->hosts); + if (!table) { + fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", skylist->hosts, CATDIR); + exit (1); + } + + // If we want to parallelize the output, the remote client needs to load the HostTable + // for the output catdir. Require the output HostTable to match the input (ie, same + // hostIDs must exist). Preferred, but not required, that they match in the host + // location. + HostTable *tableOut = NULL; + if (PARALLEL_OUTPUT) { + // load the list of hosts + tableOut = HostTableLoad (HIGH_SPEED_DIR, skylist->hosts); + if (!tableOut) { + fprintf (stderr, "ERROR: failure reading Host Table %s for output database %s\n", skylist->hosts, HIGH_SPEED_DIR); + exit (1); + } + if (table->Nhosts != tableOut->Nhosts) { + fprintf (stderr, "ERROR: output HostTable must have matching hosts (%d in, %d out)\n", table->Nhosts, tableOut->Nhosts); + 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; + + // use this directory to save the output catalogs (distributed if parallel) + char *outputDir = strcreate(HIGH_SPEED_DIR); + if (PARALLEL_OUTPUT) { + // find the matching host entry in the output HostTable (match by ID, not index) + int hostID = table->hosts[i].hostID; + int index = table->index[hostID]; + if (index == -1) { + fprintf (stderr, "ERROR: output HostTable must have matching hosts (host ID %d not found)\n", hostID); + exit (1); + } + + free (outputDir); + outputDir = abspath (tableOut->hosts[index].pathname, MAX_PATH_LENGTH); + } + + char command[1024]; + 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, outputDir, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); + + free (outputDir); + + // options / arguments that can affect relastro_client -high-speed + 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 (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); + + 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) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: /trunk/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- /trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 33652) @@ -6,36 +6,36 @@ 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; - int foundA, foundB; + 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, XVERB; 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, VERBOSE2,"w")) { fprintf (stderr, "ERROR: failure to open catalog file %s\n", filename); @@ -45,7 +45,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): @@ -77,4 +77,11 @@ // mark (exclude) objects with both sets of target photcodes for (i = 0; i < catalog[0].Naverage; i++) { + + XVERB = (catalog[0].average[i].objID == OBJ_ID_SRC) && (catalog[0].average[i].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[i].objID == OBJ_ID_DST) && (catalog[0].average[i].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object\n"); + } + // do any of the measures for this object match group A? m = catalog[0].average[i].measureOffset; @@ -114,7 +121,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 +131,5 @@ NgroupBbad++; } + continue; } // this object does not have a detection matching the contraints from either gorupA or groupB -- skip as if slow @@ -173,4 +181,10 @@ nj = N1[j]; + XVERB = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object %d or %d\n", (int) ni, (int) nj); + } + if (slowMoving[ni]) NEXT_I; if (slowMoving[nj]) NEXT_J; @@ -185,4 +199,7 @@ // look for pairs that are within the maximum separation dX = X1[i] - X1[j]; + if (XVERB) { + fprintf (stderr, "%d %d : %f\n", (int) ni, (int) nj, dX); + } if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small @@ -197,4 +214,11 @@ nj = N1[J]; if (!groupB[nj]) continue; + + XVERB = (catalog[0].average[ni].objID == OBJ_ID_SRC) && (catalog[0].average[ni].catID == CAT_ID_SRC); + XVERB |= (catalog[0].average[nj].objID == OBJ_ID_DST) && (catalog[0].average[nj].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "test object pt2 %d or %d : %f\n", (int) ni, (int) nj, dX); + fprintf (stderr, ".\n"); + } dY = Y1[i] - Y1[J]; @@ -209,16 +233,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;lphotcodes match the photcodeSet +// return TRUE if measure->photcode match any in the photcodeSet +// (but, return FALSE if the measurement is bad) int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) { int found, k; - assert (Nset > 0); - - if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE; - if (!finite(measure[0].M)) return FALSE; + if (!Nset) return TRUE; + // XXX require a set or not? assert (Nset > 0); + + if (!finite(measure[0].dR)) return FALSE; + if (!finite(measure[0].dD)) return FALSE; + if (!finite(measure[0].M)) return FALSE; + float dX = GetAstromError (measure, ERROR_MODE_RA); + if (isnan(dX)) return FALSE; + + float dY = GetAstromError (measure, ERROR_MODE_DEC); + if (isnan(dY)) return FALSE; + /* select measurements by photcode, or equiv photcode, if specified */ - if (Nset > 0) { - found = FALSE; - for (k = 0; (k < Nset) && !found; k++) { - if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE; - if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; - } - if (!found) return FALSE; - } - - return TRUE; + found = FALSE; + for (k = 0; (k < Nset) && !found; k++) { + if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE; + if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; + } + + return found; } Index: /trunk/Ohana/src/relastro/src/high_speed_utils.c =================================================================== --- /trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 33652) @@ -17,8 +17,10 @@ static dbStack *stackA = NULL; static int NstackA = 0; + static dbField *fieldsB = NULL; static int NfieldsB = 0; static dbStack *stackB = NULL; static int NstackB = 0; + static dbValue *valuesA; // static int NvaluesA; @@ -41,5 +43,5 @@ ALLOCATE(words, char *, NwordsAllocated); - char *copy_of_str = strdup(string); + char *copy_of_str = strcreate(string); char *saveptr = NULL; @@ -60,4 +62,7 @@ // Create the data strucutres required to apply avextract like cuts int setupConstraints(char *whereString, dbField **pFields, int *pNfields, dbStack **pStack, int *pNstack) { + + // an empty where is allowed (no further filtering) + if (!whereString[0]) return TRUE; // split the string into words @@ -89,11 +94,18 @@ int initializeConstraints() { + // if constaints haven't been set by the user use 2mass - if (!WHERE_A[0]) { - strcpy(WHERE_A, DEFAULT_WHERE_A); - } - if (!WHERE_B[0]) { - strcpy(WHERE_B, DEFAULT_WHERE_B); - } + // XXX if (!WHERE_A[0]) { + // XXX strcpy(WHERE_A, DEFAULT_WHERE_A); + // XXX } + // XXX if (!WHERE_B[0]) { + // XXX strcpy(WHERE_B, DEFAULT_WHERE_B); + // XXX } + + fprintf (stderr, "where A: %s\n", WHERE_A); + fprintf (stderr, "where B: %s\n", WHERE_B); + + // pull TIMEFORMAT & TIMEREF from the Config system + dbExtractAveragesInit(); if (! setupConstraints(WHERE_A, &fieldsA, &NfieldsA, &stackA, &NstackA)) { @@ -102,5 +114,8 @@ exit(1); } - ALLOCATE(valuesA, dbValue, NfieldsA); + if (NfieldsA) { + ALLOCATE(valuesA, dbValue, NfieldsA); + } + if (! setupConstraints(WHERE_B, &fieldsB, &NfieldsB, &stackB, &NstackB)) { print_error(); @@ -108,5 +123,8 @@ exit(1); } - ALLOCATE(valuesB, dbValue, NfieldsB); + if (NfieldsB) { + ALLOCATE(valuesB, dbValue, NfieldsB); + } + return TRUE; } @@ -116,8 +134,10 @@ int applyConstraintsA(Catalog *catalog, off_t i) { off_t m = catalog[0].average[i].measureOffset; - off_t n; int Nsecfilt = GetPhotcodeNsecfilt(); - dbExtractAveragesInit(); + // an empty WHERE matches all objects + if (!NfieldsA) return TRUE; + + off_t n; for (n = 0; n < NfieldsA; n++) { valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]); @@ -130,5 +150,6 @@ int Nsecfilt = GetPhotcodeNsecfilt(); - dbExtractAveragesInit(); + // an empty WHERE matches all objects + if (!NfieldsB) return TRUE; off_t n; Index: /trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- /trunk/Ohana/src/relastro/src/initialize.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/initialize.c (revision 33652) @@ -2,9 +2,4 @@ void initialize (int argc, char **argv) { - - int NPHOTCODES; - char *codename, *ptr, *list; - - ptr = NULL; ConfigInit (&argc, argv); @@ -13,90 +8,13 @@ if (FIT_TARGET == TARGET_MERGE_SOURCE) return; - /* build a list of accepted photcodes. these will be used by bcatalog to accept or - reject loaded data */ + fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); + fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); - NphotcodesKeep = 0; - photcodesKeep = NULL; - if (PHOTCODE_KEEP_LIST != NULL) { - NPHOTCODES = 10; - ALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES); - - /* parse the comma-separated list of photcodesKeep */ - list = PHOTCODE_KEEP_LIST; - while ((codename = strtok_r (list, ",", &ptr)) != NULL) { - list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_KEEP_LIST); - fprintf (stderr, "codename: %s\n", codename); - if ((photcodesKeep[NphotcodesKeep] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); - } - NphotcodesKeep ++; - CHECK_REALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES, NphotcodesKeep, 10); - } - } - - NphotcodesSkip = 0; - photcodesSkip = NULL; - if (PHOTCODE_SKIP_LIST != NULL) { - NPHOTCODES = 10; - ALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES); - - /* parse the comma-separated list of photcodesSkip */ - list = PHOTCODE_SKIP_LIST; - while ((codename = strtok_r (list, ",", &ptr)) != NULL) { - list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_SKIP_LIST); - fprintf (stderr, "codename: %s\n", codename); - if ((photcodesSkip[NphotcodesSkip] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); - } - NphotcodesSkip ++; - CHECK_REALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES, NphotcodesSkip, 10); - } - } - - NphotcodesGroupA = 0; - photcodesGroupA = NULL; - if (PHOTCODE_A_LIST != NULL) { - NPHOTCODES = 10; - ALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES); - - /* parse the comma-separated list of photcodesGroupA */ - fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); - list = PHOTCODE_A_LIST; - while ((codename = strtok_r (list, ",", &ptr)) != NULL) { - list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "codename: %s\n", codename); - if ((photcodesGroupA[NphotcodesGroupA] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); - } - NphotcodesGroupA ++; - CHECK_REALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES, NphotcodesGroupA, 10); - } - } - - NphotcodesGroupB = 0; - photcodesGroupB = NULL; - if (PHOTCODE_B_LIST != NULL) { - NPHOTCODES = 10; - ALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES); - - fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); - /* parse the comma-separated list of photcodesGroupB */ - list = PHOTCODE_B_LIST; - while ((codename = strtok_r (list, ",", &ptr)) != NULL) { - list = NULL; // pass NULL on successive strtok_r calls - fprintf (stderr, "codename: %s\n", codename); - if ((photcodesGroupB[NphotcodesGroupB] = GetPhotcodebyName (codename)) == NULL) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); - exit (1); - } - NphotcodesGroupB ++; - CHECK_REALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES, NphotcodesGroupB, 10); - } - } + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); + photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); initstats (STATMODE); @@ -123,2 +41,24 @@ } +void initialize_client (int argc, char **argv) { + + ConfigInit (&argc, argv); + args_client (argc, argv); + + fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); + fprintf (stderr, "PHOTCODE_SKIP_LIST: %s\n", PHOTCODE_SKIP_LIST); + fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST); + fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST); + + photcodesKeep = ParsePhotcodeList (PHOTCODE_KEEP_LIST, &NphotcodesKeep, FALSE); + photcodesSkip = ParsePhotcodeList (PHOTCODE_SKIP_LIST, &NphotcodesSkip, FALSE); + photcodesGroupA = ParsePhotcodeList (PHOTCODE_A_LIST, &NphotcodesGroupA, TRUE); + photcodesGroupB = ParsePhotcodeList (PHOTCODE_B_LIST, &NphotcodesGroupB, TRUE); + + initstats (STATMODE); + + // IMAGE_BAD = ID_IMAGE_ASTROM_POOR | ID_IMAGE_ASTROM_FEW | ID_IMAGE_ASTROM_SKIP; + // STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + // MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR_ASTROM | ID_MEAS_SKIP_ASTROM | ID_MEAS_AREA; +} + Index: /trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/load_catalogs.c (revision 33652) @@ -1,17 +1,28 @@ # include "relastro.h" -Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) { +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog); +void bcatalog_show_skips (); + +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath) { int i, j, k, m, Nstar; Catalog *catalog, *pcatalog, tcatalog; + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + catalog = load_catalogs_parallel (skylist, Ncatalog); + return catalog; + } + if (VERBOSE) fprintf (stderr, "loading catalog data\n"); ALLOCATE (catalog, Catalog, skylist[0].Nregions); - initStarMaps(); - // load data from each region file, only use bright stars + int Ncat = 0; for (i = 0; i < skylist[0].Nregions; i++) { + + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; // we only allow output if we do not use a subset. in this case, @@ -19,46 +30,68 @@ pcatalog = subselect ? &tcatalog : &catalog[i]; + // define the catalog file name + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + pcatalog->filename = hostID ? hostfile : skylist[0].filename[i]; + // set up the basic catalog info - pcatalog[0].filename = skylist[0].filename[i]; - pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data - pcatalog[0].catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data - pcatalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point - pcatalog[0].Nsecfilt = GetPhotcodeNsecfilt (); - + pcatalog->catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + pcatalog->catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + pcatalog->catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + pcatalog->Nsecfilt = GetPhotcodeNsecfilt (); + + // loads Average, Measure, SecFilt if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE2, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename); exit (1); } - if (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; } } } } - } - - // XXX keep this test? + Ncat ++; + } + + // XXX TEST : bcatalog_show_skips(); + Nstar = 0; - for (i = 0; i < skylist[0].Nregions; i++) { + for (i = 0; i < Ncat; i++) { Nstar += catalog[i].Naverage; + if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) { + fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n"); + fprintf (stderr, " If you need to do this, and you can afford the RAM, rebuild with IDX_T set to off_t (see relastro.h, ImagesOps.c)\n"); + exit (3); + } } if (Nstar < 2) { @@ -66,6 +99,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); } @@ -76,2 +110,132 @@ need to use an XCLD lock here. */ + +// CATDIR is supplied globally +# define DEBUG 1 +Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) { + + // 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; + + ALLOCATE (table->hosts[i].results, char, 1024); + snprintf (table->hosts[i].results, 1024, "%s/relastro.catalog.subset.dat", table->hosts[i].pathname); + + // 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 -load-objects %s -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f -D RELASTRO_SIGMA_LIM %f", + table->hosts[i].results, table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR, SIGMA_LIM); + + char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (ImagSelect) { snprintf (tmpline, 1024, "%s -instmag %f %f", command, ImagMin, ImagMax); strcpy (command, tmpline); } + if (MaxDensityUse) { snprintf (tmpline, 1024, "%s -max-density %f", command, MaxDensityValue); strcpy (command, tmpline); } + + if (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) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + // each host generates a BrightCatalog structure, with the measure, average, etc value + // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt). I need to + // split out the per-catalog measurements into separate catalog entries. + + // set up an initial array of catalogs + int Nsecfilt = GetPhotcodeNsecfilt (); + CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt); + + for (i = 0; i < table->Nhosts; i++) { + + BrightCatalog *bcatalog = BrightCatalogLoad (table->hosts[i].results); + assert (bcatalog); + + BrightCatalogSplit (catalogs, bcatalog); + + free (bcatalog->average); + free (bcatalog->measure); + free (bcatalog->secfilt); + free (bcatalog); + } + + Catalog *catalog = catalogs->catalog; + *Ncatalog = catalogs->Ncatalog; + BrightCatalogSplitFree (catalogs); + + int Nmeasure = 0; + int Naverage = 0; + for (i = 0; i < catalogs->Ncatalog; i++) { + Nmeasure += catalogs->catalog[i].Nmeasure; + Naverage += catalogs->catalog[i].Naverage; + } + + fprintf (stderr, "loaded %d catalogs, using a total of %d stars (%d measures)\n", catalogs->Ncatalog, Naverage, Nmeasure); + + return (catalog); +} Index: /trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- /trunk/Ohana/src/relastro/src/load_images.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/load_images.c (revision 33652) @@ -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: /trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/relastro.c (revision 33652) @@ -3,118 +3,41 @@ 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); - /* the object analysis is a separate process iterating over catalogs */ - 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); - /* the object analysis is a separate process iterating over catalogs */ - 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); + case TARGET_SIMPLE: + case TARGET_CHIPS: + case TARGET_MOSAICS: + relastro_images (skylist); + exit (0); - /* 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; + case TARGET_UPDATE_OFFSETS: + // iterate over catalogs to make detection coordinates consistant + UpdateObjectOffsets (skylist, 0, NULL); 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: /trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_client.c (revision 33652) +++ /trunk/Ohana/src/relastro/src/relastro_client.c (revision 33652) @@ -0,0 +1,102 @@ +# include "relastro.h" + +// relastro_client is run on a remote host and is responsible for updating the catalogs +// owned by that host. + +// there are four modes: + +// 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_client -update-offsets : load image table containing the updated astrometry +// information and determine per-object astrometry + +// 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) { + + // get configuration info, args, lockfile (set CATDIR, HOST_ID, HOSTDIR, etc) + initialize_client (argc, argv); + + SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + SkyList *skylist = SkyListByPatch (sky, -1, &UserPatch); + + switch (FIT_TARGET) { + + case TARGET_LOAD_OBJECTS: { + // USAGE: relastro_client -load-objects + int Ncatalog; + Catalog *catalog = load_catalogs (skylist, &Ncatalog, TRUE, HOST_ID, HOSTDIR); + if (!catalog) { + fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR); + exit (2); + } + BrightCatalog *bcatalog = BrightCatalogMerge (catalog, Ncatalog); + if (!BrightCatalogSave (BCATALOG, bcatalog)) { + fprintf (stderr, "ERROR saving bright catalog from %s\n", CATDIR); + exit (2); + } + break; + } + + case TARGET_UPDATE_OBJECTS: { + // USAGE: relastro_client -update-objects + relastro_objects (skylist, HOST_ID, HOSTDIR); + break; + } + + 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; + } + + default: + fprintf (stderr, "impossible!"); + abort(); + } + + exit (0); +} Index: /trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_images.c (revision 33652) +++ /trunk/Ohana/src/relastro/src/relastro_images.c (revision 33652) @@ -0,0 +1,112 @@ +# 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); + + if (Ncatalog == 0) { + fprintf (stderr, "ERROR: no valid data for relastro, exiting\n"); + exit (2); + } + + /* 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: /trunk/Ohana/src/relastro/src/relastro_merge_source.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_merge_source.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/relastro_merge_source.c (revision 33652) @@ -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 @@ -52,5 +47,5 @@ catalog_src.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catalog_src, region_src, VERBOSE, "w")) { + if (!dvo_catalog_open (&catalog_src, region_src, VERBOSE2, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog_src.filename); exit (1); @@ -127,5 +122,5 @@ resort_catalog (&catalog_src); - dvo_catalog_save (&catalog_src, VERBOSE); + dvo_catalog_save (&catalog_src, VERBOSE2); dvo_catalog_unlock (&catalog_src); dvo_catalog_free (&catalog_src); Index: /trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro_objects.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/relastro_objects.c (revision 33652) @@ -1,22 +1,32 @@ # include "relastro.h" -int relastro_objects () { +int relastro_objects_parallel (SkyList *sky); + +int relastro_objects (SkyList *skylist, int hostID, char *hostpath) { int i, j, k, m; - SkyTable *sky = NULL; - SkyList *skylist = NULL; Catalog catalog; + DVOMeasureFlags measureBits = + ID_MEAS_POOR_ASTROM | + ID_MEAS_SKIP_ASTROM | + ID_MEAS_USED_OBJ | + ID_MEAS_USED_CHIP; - // 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); + DVOAverageFlags averageBits = + ID_STAR_FEW | + ID_STAR_POOR | + ID_STAR_FIT_AVE | + ID_STAR_FIT_PM | + ID_STAR_FIT_PAR | + ID_STAR_USE_AVE | + ID_STAR_USE_PM | + ID_STAR_USE_PAR; + + // XXX need to decide how to determine PARALLEL mode... + if (PARALLEL && !hostID) { + relastro_objects_parallel (skylist); + return TRUE; } @@ -24,6 +34,11 @@ for (i = 0; i < skylist[0].Nregions; i++) { + // does this host ID match the desired location for the table? + if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue; + // set up the basic catalog info - catalog.filename = skylist[0].filename[i]; + char hostfile[1024]; + snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name); + catalog.filename = hostID ? hostfile : 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 @@ -31,10 +46,11 @@ catalog.Nsecfilt = GetPhotcodeNsecfilt (); - if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + // loads Average, Measure, SecFilt + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE2, "w")) { fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); exit (1); } if (!catalog.Naves_disk) { - if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + if (VERBOSE2) fprintf (stderr, "no data in %s, skipping\n", catalog.filename); dvo_catalog_unlock (&catalog); dvo_catalog_free (&catalog); @@ -42,20 +58,25 @@ } + // operates on full tables (Average, Measure, Secfilt) if (FlagOutlier) { FlagOutliers(&catalog); } - // XXX consider what gets reset (only ASTROM flags) + // reset only the astrometry-related average and measure bits if (RESET) { for (j = 0; j < catalog.Naverage; j++) { - catalog.average[j].flags = 0; + catalog.average[j].flags &= ~averageBits; m = catalog.average[j].measureOffset; for (k = 0; k < catalog.average[j].Nmeasure; k++) { - catalog.measure[m+k].dbFlags = 0; + catalog.measure[m+k].dbFlags &= ~measureBits; } } } + populate_tiny_values(&catalog, DVO_TV_MEASURE); + UpdateObjects (&catalog, 1); + + free_tiny_values(&catalog); if (!UPDATE) { @@ -70,2 +91,104 @@ return (TRUE); } + +// CATDIR is supplied globally +# define DEBUG 1 +int relastro_objects_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 -update-objects: + // 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 + + // FIT_MODE + // PM_TOOFEW + // SRC_MEAS_TOOFEW + // FIT_TARGET + + char command[1024]; + snprintf (command, 1024, "relastro_client -update-objects -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); + + char tmpline[1024]; + if (FIT_MODE == FIT_PM_ONLY) { snprintf (tmpline, 1024, "%s -pm", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PAR_ONLY) { snprintf (tmpline, 1024, "%s -par", command); strcpy (command, tmpline); } + if (FIT_MODE == FIT_PM_AND_PAR) { snprintf (tmpline, 1024, "%s -pmpar", command); strcpy (command, tmpline); } + + if (VERBOSE) { snprintf (tmpline, 1024, "%s -v", command); strcpy (command, tmpline); } + if (VERBOSE2) { snprintf (tmpline, 1024, "%s -vv", command); strcpy (command, tmpline); } + if (RESET) { snprintf (tmpline, 1024, "%s -reset", command); strcpy (command, tmpline); } + if (UPDATE) { snprintf (tmpline, 1024, "%s -update", command); strcpy (command, tmpline); } + + if (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) { + HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} Index: /trunk/Ohana/src/relastro/src/resort_catalog.c =================================================================== --- /trunk/Ohana/src/relastro/src/resort_catalog.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/resort_catalog.c (revision 33652) @@ -1,6 +1,3 @@ # include "relastro.h" - -# define myAbort(MSG) { fprintf (stderr, "%s\n", MSG); abort(); } -# define myAssert(LOGIC,MSG) { if (!(LOGIC)) { fprintf (stderr, "%s\n", MSG); abort(); } } // sort the measure Sequence based on the average Sequence entries Index: /trunk/Ohana/src/relastro/src/save_catalogs.c =================================================================== --- /trunk/Ohana/src/relastro/src/save_catalogs.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/save_catalogs.c (revision 33652) @@ -8,6 +8,6 @@ for (i = 0; i < Ncatalog; i++) { - if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename); - dvo_catalog_save (&catalog[i], VERBOSE); + if (VERBOSE2) fprintf (stderr, "saving catalog %s\n", catalog[i].filename); + dvo_catalog_save (&catalog[i], VERBOSE2); dvo_catalog_unlock (&catalog[i]); dvo_catalog_free (&catalog[i]); Index: /trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- /trunk/Ohana/src/relastro/src/select_images.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/select_images.c (revision 33652) @@ -25,5 +25,5 @@ Image *image; - off_t i, j, k, m, nStart, iSky, nimage, NIMAGE; + off_t i, j, k, m, nStart, iSky, nimage, NIMAGE, D_NIMAGE; off_t *line_number; int InRange, found; @@ -68,4 +68,6 @@ DminSkyRegion = +90.0; DmaxSkyRegion = -90.0; + + D_NIMAGE = 6000; /* compare with each region file */ @@ -113,5 +115,5 @@ nimage = 0; - NIMAGE = 100; + NIMAGE = D_NIMAGE; ALLOCATE (image, Image, NIMAGE); ALLOCATE (line_number, off_t, NIMAGE); @@ -167,6 +169,6 @@ /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */ - double RminImage = 360.0; - double RmaxImage = 0.0; + double RminImage = RmidSkyRegion + 180.0; + double RmaxImage = RmidSkyRegion - 180.0; double DminImage = +90.0; double DmaxImage = -90.0; @@ -194,4 +196,19 @@ if (RmaxImage < RminSkyRegion) continue; + if (!strncmp(timage[i].name, "o5903g0638o", 10) && (timage[i].photcode == 10355)) { + fprintf (stderr, "test image\n"); + } + + // require that the full image be inside the region of interest (only for non PHU + // images) XXX : if we calibrate the mosaic, require that the full mosaic be inside + // the region and all of its chips.. +# define FULL_OVERLAP 1 + if (FULL_OVERLAP && strcmp(&timage[i].coords.ctype[4], "-DIS")) { + if (RmaxImage > UserPatch.Rmax) continue; + if (RminImage < UserPatch.Rmin) continue; + if (DmaxImage > UserPatch.Dmax) continue; + if (DminImage < UserPatch.Dmin) continue; + } + // image overlaps region, keep it if (USE_BASIC_CHECK) goto found_it; @@ -243,8 +260,5 @@ } if (RESET) { - // XXX do we need / want to do this in relastro? - // assignMcal (&image[nimage], (double *) NULL, -1); - // image[nimage].Mcal = NAN; - // image[nimage].dMcal = NAN; + // this only resets the astrometry image flags, not the photometry ones image[nimage].flags &= ~badImage; } @@ -252,5 +266,6 @@ nimage ++; if (nimage == NIMAGE) { - NIMAGE += 100; + NIMAGE += D_NIMAGE; + D_NIMAGE = MAX (100000, D_NIMAGE * 1.5); REALLOCATE (image, Image, NIMAGE); REALLOCATE (line_number, off_t, NIMAGE); @@ -264,4 +279,6 @@ REALLOCATE (line_number, off_t, MAX (nimage, 1)); free (skycoords); + free (RmaxSky); + free (index); *Nimage = nimage; Index: /trunk/Ohana/src/relastro/src/testparallax.c =================================================================== --- /trunk/Ohana/src/relastro/src/testparallax.c (revision 33651) +++ /trunk/Ohana/src/relastro/src/testparallax.c (revision 33652) @@ -1,31 +1,114 @@ # include "relastro.h" +# define NPTS 64 +# define J2000 2451545. /* Julian date at standard epoch */ +# define T1970 2440587.500000 /* JD at UNIX ref time */ int main (int argc, char **argv) { + + int i, Npts; + char line[1024]; + double R[NPTS], D[NPTS], Time[NPTS], JD[NPTS]; + double X[NPTS], Y[NPTS], dX[NPTS], dY[NPTS], pX[NPTS], pY[NPTS]; + double Tref[NPTS], Tjyrs[NPTS], TrefS, TjyrsS, TrefMean, TjyrsMean, Ro, Do; + Coords coords; + PMFit fitPAR; - if (argc != N) { - fprintf (stderr, "USAGE: testparallax (file.dat)\n"); - exit (1); + if (argc != 2) { + fprintf (stderr, "USAGE: testparallax (file.dat)\n"); + exit (1); + } + + // test parallax program + FILE *f = fopen (argv[1], "r"); + if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); } + + // XXX uncomment to skip first line scan_line(f, line); + + Npts = 0; + for (i = 0; scan_line(f, line) != EOF; i++) { + if (Npts == NPTS) { + fprintf (stderr, "too many point: use dynamic alloc\n"); + exit (2); } + if (line[0] == '#') continue; + dparse (&R[Npts], 1, line); + dparse (&D[Npts], 2, line); + dparse (&JD[Npts], 3, line); + dparse (&Time[Npts], 3, line); - // test parallax program - FILE *f = fopen (argv[1], "r"); - if (!f) { fprintf (stderr, "failed to open %s\n", argv[1]); exit (2); } + dparse (&dX[Npts], 4, line); + dparse (&dY[Npts], 4, line); - char buffer[1024]; - scan_line(f, line); - - while (scan_line(f, line) != EOF) { - dparse (&jd, 1, line); - dparse (&time, 2, line); - dparse (&ra, 3, line); - dparse (&dec, 4, line); + // trent's file + // dparse (&JD[Npts], 1, line); + // dparse (&Time[Npts], 2, line); + // dparse (&R[Npts], 3, line); + // dparse (&D[Npts], 4, line); + Npts ++; + } - } + /* project coordinates to a plane centered on the object with units of arcsec */ + coords.crval1 = 0; + coords.crval2 = 0; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "RA---SIN"); - // project to local coords + // use one point as a local reference + coords.crval1 = R[0]; + coords.crval2 = D[0]; - // run fitter - FitPMandPar + TrefS = TjyrsS = 0.0; + // project to local coords + for (i = 0; i < Npts; i++) { + RD_to_XY (&X[i], &Y[i], R[i], D[i], &coords); + Tjyrs[i] = Time[i] / 365.25; + Tref[i] = Time[i] / 365.25; + // Tunix[i] = (JD[i] - T1970) / 365.25; // time relative to T1970 in years + // Tjyrs[i] = (JD[i] - J2000) / 365.25; // time relative to J2000 in years + + TrefS += Tref[i]; + TjyrsS += Tjyrs[i]; + + // nominal for PS1 + // dX[i] = 0.020; + // dY[i] = 0.020; + } + TrefMean = TrefS / Npts; + TjyrsMean = TjyrsS / Npts; + + for (i = 0; i < Npts; i++) { + Tref[i] -= TrefMean; + ParFactor (&pX[i], &pY[i], R[i], D[i], Tjyrs[i], 0.0); + fprintf (stderr, "%f %f : %f %f : %f %f : %f %f\n", R[i], D[i], X[i], Y[i], Tref[i], Tjyrs[i], pX[i], pY[i]); + } + + // run fitter + // FitPMandPar (&fitPAR, X, dX, Y, dY, Tref, pX, pY, Npts, TRUE); + FitPMandPar (&fitPAR, X, dX, Y, dY, Tjyrs, pX, pY, Npts, TRUE); + XY_to_RD (&Ro, &Do, fitPAR.Ro, fitPAR.Do, &coords); + + fprintf (stderr, "Rx : %f\n", fitPAR.Ro); + fprintf (stderr, "Dx : %f\n", fitPAR.Do); + fprintf (stderr, "Ro : %f\n", Ro); + fprintf (stderr, "Do : %f\n", Do); + fprintf (stderr, "dRo : %f\n", fitPAR.dRo); + fprintf (stderr, "dDo : %f\n", fitPAR.dDo); + fprintf (stderr, "uR : %f\n", fitPAR.uR); + fprintf (stderr, "uD : %f\n", fitPAR.uD); + fprintf (stderr, "duR : %f\n", fitPAR.duR); + fprintf (stderr, "duD : %f\n", fitPAR.duD); + fprintf (stderr, "p : %f\n", fitPAR.p); + fprintf (stderr, "dp : %f\n", fitPAR.dp); + fprintf (stderr, "uTot : %f\n", hypot(fitPAR.uR,fitPAR.uD)); + fprintf (stderr, "PA : %f\n", DEG_RAD*atan2(fitPAR.uR,fitPAR.uD) + 360.0); + + exit (0); } +