Index: trunk/Ohana/src/libdvo/Makefile =================================================================== --- trunk/Ohana/src/libdvo/Makefile (revision 37729) +++ trunk/Ohana/src/libdvo/Makefile (revision 37807) @@ -38,4 +38,5 @@ $(DESTINC)/ps1_v5_defs.h \ $(DESTINC)/ps1_ref_defs.h \ +$(DESTINC)/ps1_sim_defs.h \ $(DESTINC)/cmf-ps1-dv3.h \ $(DESTINC)/cmf-ps1-v5.h \ @@ -50,4 +51,5 @@ $(SRC)/version.$(ARCH).o \ $(SRC)/coordops.$(ARCH).o \ +$(SRC)/galaxy_model.$(ARCH).o \ $(SRC)/dvosorts.$(ARCH).o \ $(SRC)/dvo_photcode_ops.$(ARCH).o \ @@ -84,8 +86,12 @@ $(SRC)/dvo_convert_PS1_V5.$(ARCH).o \ $(SRC)/dvo_convert_PS1_REF.$(ARCH).o \ +$(SRC)/dvo_convert_PS1_SIM.$(ARCH).o \ $(SRC)/flatcorr_io.$(ARCH).o \ $(SRC)/skyregion_io.$(ARCH).o \ $(SRC)/skyregion_gsc.$(ARCH).o \ $(SRC)/skyregion_ops.$(ARCH).o \ +$(SRC)/AstromOffsetMapIO.$(ARCH).o \ +$(SRC)/AstromOffsetMapOps.$(ARCH).o \ +$(SRC)/AstromOffsetMapUtils.$(ARCH).o \ $(SRC)/dvo_set_skyregion.$(ARCH).o \ $(SRC)/dvo_set_catdir.$(ARCH).o \ Index: trunk/Ohana/src/libdvo/include/dvo.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvo.h (revision 37729) +++ trunk/Ohana/src/libdvo/include/dvo.h (revision 37807) @@ -24,4 +24,5 @@ DVO_FORMAT_PS1_DEV_3, DVO_FORMAT_PS1_REF, + DVO_FORMAT_PS1_SIM, DVO_FORMAT_PS1_V1, DVO_FORMAT_PS1_V2, @@ -36,29 +37,4 @@ typedef enum {DVO_CAT_OPEN_FAIL, DVO_CAT_OPEN_OK, DVO_CAT_OPEN_EMPTY} DVOCatalogOpenModes; - -typedef enum { - PROJ_NONE, // undefined - PROJ_ZEA, // zenithal - PROJ_ZPL, // zenithal - PROJ_ARC, // zenithal - PROJ_STG, // zenithal - PROJ_SIN, // zenithal - PROJ_TAN, // zenithal - PROJ_TNX, // zenithal - PROJ_DIS, // zenithal (TAN + polyterms) - PROJ_LIN, // cartesian - PROJ_PLY, // cartesian - PROJ_WRP, // cartesian - PROJ_AIT, // pseudocyl - PROJ_GLS, // pseudocyl - PROJ_PAR, // pseudocyl -} OhanaProjection; - -typedef enum { - PROJ_MODE_NONE, - PROJ_MODE_CARTESIAN, - PROJ_MODE_ZENITHAL, - PROJ_MODE_PSEUDOCYL, -} OhanaProjectionMode; // these are used as NAN for types of int values @@ -79,17 +55,19 @@ /* catalog values to be loaded */ -# define LOAD_NONE 0x000 -# define LOAD_AVES 0x001 -# define LOAD_MEAS 0x002 -# define LOAD_MISS 0x004 -# define LOAD_SECF 0x008 -# define SKIP_AVES 0x010 -# define SKIP_MEAS 0x020 -# define SKIP_MISS 0x040 -# define SKIP_SECF 0x080 -# define LOAD_LENSING 0x100 -# define LOAD_LENSOBJ 0x200 -# define SKIP_LENSING 0x400 -# define SKIP_LENSOBJ 0x800 +# define LOAD_NONE 0x0000 +# define LOAD_AVES 0x0001 +# define LOAD_MEAS 0x0002 +# define LOAD_MISS 0x0004 +# define LOAD_SECF 0x0008 +# define SKIP_AVES 0x0010 +# define SKIP_MEAS 0x0020 +# define SKIP_MISS 0x0040 +# define SKIP_SECF 0x0080 +# define LOAD_LENSING 0x0100 +# define LOAD_LENSOBJ 0x0200 +# define SKIP_LENSING 0x0400 +# define SKIP_LENSOBJ 0x0800 +# define LOAD_STARPAR 0x1000 +# define SKIP_STARPAR 0x2000 /* photometry code types */ @@ -114,4 +92,5 @@ ID_IMAGE_ASTROM_FEW = 0x00000100, /* currently too few measurements for astrometry */ ID_IMAGE_PHOTOM_UBERCAL = 0x00000200, /* externally-supplied photometry zero point from ubercal analysis */ + ID_IMAGE_ASTROM_GMM = 0x00000400, /* image was fitted to positions corrected by the galaxy motion model */ } DVOImageFlags; @@ -136,4 +115,6 @@ ID_MEAS_STACK_PRIMARY = 0x00010000, // this stack measurement is in the primary skycell ID_MEAS_STACK_PHOT_SRC = 0x00020000, // this measurement supplied the stack photometry + ID_MEAS_ICRF_QSO = 0x00040000, // this measurement is an ICRF reference position + ID_MEAS_IMAGE_EPOCH = 0x00080000, // this measurement is registered to the image epoch (not tied to ref catalog epoch) } DVOMeasureFlags; @@ -152,34 +133,40 @@ */ +/** these names were previously used for Average flags in old (LONEOS-era) dvo versions + ID_BAD_OBJECT = 0x00004000, // if all measurements are bad, set this bit + ID_MOVING = 0x00008000, // is a moving object + ID_ROCK = 0x0000a000, // 0x8000 + 0x2000 + ID_GHOST = 0x0000c001, // 0x8000 + 0x4000 + 0x0001 + ID_TRAIL = 0x0000c002, // 0x8000 + 0x4000 + 0x0002 + ID_BLEED = 0x0000c003, // 0x8000 + 0x4000 + 0x0003 + ID_COSMIC = 0x0000c004, // 0x8000 + 0x4000 + 0x0004 +**/ + /* Average.flags values -- these values are 32 bit (as of PS1_V1) */ typedef enum { - ID_STAR_FEW = 0x00000001, // used within relphot: skip star - ID_STAR_POOR = 0x00000002, // used within relphot: skip star - ID_PROPER = 0x00000400, // star with large proper motion - ID_TRANSIENT = 0x00001000, // is this mutually exclusive with USNO? - ID_VARIABLE = 0x00002000, // not currently set? - ID_ASTEROID = 0x00002000, // identified with an asteroid - ID_BAD_OBJECT = 0x00004000, // if all measurements are bad, set this bit - ID_MOVING = 0x00008000, // is a moving object - ID_ROCK = 0x0000a000, // 0x8000 + 0x2000 - ID_GHOST = 0x0000c001, // 0x8000 + 0x4000 + 0x0001 - ID_TRAIL = 0x0000c002, // 0x8000 + 0x4000 + 0x0002 - ID_BLEED = 0x0000c003, // 0x8000 + 0x4000 + 0x0003 - ID_COSMIC = 0x0000c004, // 0x8000 + 0x4000 + 0x0004 - ID_STAR_FIT_AVE = 0x00010000, // average position fitted - ID_STAR_FIT_PM = 0x00020000, // proper motion fitted - ID_STAR_FIT_PAR = 0x00040000, // parallax fitted - ID_STAR_USE_AVE = 0x00080000, // average position used (not PM or PAR) - ID_STAR_USE_PM = 0x00100000, // proper motion used (not AVE or PAR) - ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) - ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured - ID_STAR_BAD_PM = 0x00800000, // mean astrometry could not be measured - ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) - ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) - ID_OBJ_GOOD = 0x04000000, // good-quality measurement in our data (eg,PS) - ID_OBJ_GOOD_ALT = 0x08000000, // good-quality measurement in external data (eg, 2MASS) - ID_OBJ_GOOD_STACK = 0x10000000, // good-quality object in the stack (> 1 good stack) + ID_STAR_FEW = 0x00000001, // used within relphot: skip star + ID_STAR_POOR = 0x00000002, // used within relphot: skip star + ID_ICRF_QSO = 0x00000004, // object IDed with known ICRF quasar (may have ICRF position measurement) + ID_OTHEF_QSO = 0x00000008, // object IDed with known or suspected quasar, without ICRF reference data + ID_PROPER = 0x00000010, // star with large proper motion + ID_TRANSIENT = 0x00000020, // identified as a non-periodic (stationary) transient + ID_VARIABLE = 0x00000040, // identified as a period variable + ID_ASTEROID = 0x00000080, // identified with a known solar-system object (asteroid or other) + // bits 0x00000100 - 0x00008000 are currently unused + ID_STAR_FIT_AVE = 0x00010000, // average position fitted + ID_STAR_FIT_PM = 0x00020000, // proper motion fitted + ID_STAR_FIT_PAR = 0x00040000, // parallax fitted + ID_STAR_USE_AVE = 0x00080000, // average position used (not PM or PAR) + ID_STAR_USE_PM = 0x00100000, // proper motion used (not AVE or PAR) + ID_STAR_USE_PAR = 0x00200000, // parallax used (not AVE or PM) + ID_STAR_NO_ASTROM = 0x00400000, // mean astrometry could not be measured + ID_STAR_BAD_PM = 0x00800000, // mean astrometry could not be measured + ID_OBJ_EXT = 0x01000000, // extended in our data (eg, PS) + ID_OBJ_EXT_ALT = 0x02000000, // extended in external data (eg, 2MASS) + ID_OBJ_GOOD = 0x04000000, // good-quality measurement in our data (eg,PS) + ID_OBJ_GOOD_ALT = 0x08000000, // good-quality measurement in external data (eg, 2MASS) + ID_OBJ_GOOD_STACK = 0x10000000, // good-quality object in the stack (> 1 good stack) ID_OBJ_SUSPECT_STACK = 0x20000000, // suspect object in the stack (> 1 good or suspect stack, < 2 good) - ID_OBJ_BAD_STACK = 0x40000000, // good-quality object in the stack (> 1 good stack) + ID_OBJ_BAD_STACK = 0x40000000, // good-quality object in the stack (> 1 good stack) } DVOAverageFlags; @@ -202,4 +189,22 @@ } DVOSecfiltFlags; +/* definitions for parallel dvo host information + XXX : need better names (safer namespace) +*/ + +typedef enum { + DATA_ON_TGT = 0x01, + DATA_ON_BCK = 0x02, + DATA_USE_BCK = 0x04, + + DATA_COPY_FAILURE = 0x80, +} SkyTableDataFlags; + +typedef enum { + HOST_STDIN = 0, + HOST_STDOUT = 1, + HOST_STDERR = 2, +} HostInfoIOfd; + /*** general dvo structures (internal use only / not IO) ***/ @@ -253,22 +258,4 @@ } FlatCorrectionTable; -/* definitions for parallel dvo host information - XXX : need better names (safer namespace) -*/ - -typedef enum { - DATA_ON_TGT = 0x01, - DATA_ON_BCK = 0x02, - DATA_USE_BCK = 0x04, - - DATA_COPY_FAILURE = 0x80, -} SkyTableDataFlags; - -typedef enum { - HOST_STDIN = 0, - HOST_STDOUT = 1, - HOST_STDERR = 2, -} HostInfoIOfd; - typedef struct { char *hostname; // name of remote machine @@ -289,4 +276,9 @@ } HostTable; +typedef struct { + int Nhosts; + HostInfo **hosts; +} HostTableGroup; + // A RegionHost processes data for some region in parallel with other regions typedef struct RegionHostInfo { @@ -314,4 +306,6 @@ Image *image; off_t *imseq; + + AstromOffsetTable *astromTable; int *neighbors; // list of neighbor index values @@ -444,4 +438,6 @@ double R; double D; + float RoffGAL; + float DoffGAL; float M; float Mcal; @@ -645,5 +641,5 @@ int lensobjOffset; // offset to mean lensing data int extendOffset; // offset to extended object entry - int paramsOffset; // offset to stellar parameter data + int starparOffset; // offset to stellar parameter data float refColorBlue; // color of astrometry ref stars float refColorRed; // color of astrometry ref stars @@ -679,9 +675,10 @@ Lensing *lensing; Lensobj *lensobj; + StarPar *starpar; int Nsecfilt; /* number of secfilt entries for each average entry */ - off_t Naverage, Nmeasure, Nmissing, Nlensing, Nlensobj, Nsecf_mem; /* current number of each component in memory */ - off_t Naves_disk, Nmeas_disk, Nmiss_disk, Nlensing_disk, Nlensobj_disk, Nsecf_disk; /* current number of each component on disk */ - off_t Naves_off, Nmeas_off, Nmiss_off, Nlensing_off, Nlensobj_off, Nsecf_off; /* index of first loaded data value */ + off_t Naverage, Nmeasure, Nmissing, Nlensing, Nlensobj, Nstarpar, Nsecf_mem; /* current number of each component in memory */ + off_t Naves_disk, Nmeas_disk, Nmiss_disk, Nlensing_disk, Nlensobj_disk, Nstarpar_disk, Nsecf_disk; /* current number of each component on disk */ + off_t Naves_off, Nmeas_off, Nmiss_off, Nlensing_off, Nlensobj_off, Nstarpar_off, Nsecf_off; /* index of first loaded data value */ // note that we use these for the full-sky relphot analysis @@ -704,4 +701,5 @@ struct Catalog *lensing_catalog; /* lensing catalog data (split) */ struct Catalog *lensobj_catalog; /* lensobj catalog data (split) */ + struct Catalog *starpar_catalog; /* starpar catalog data (split) */ unsigned int objID; @@ -739,21 +737,4 @@ int gfits_db_close PROTO((FITS_DB *db)); int gfits_db_free PROTO((FITS_DB *db)); - -/* in coords.c, using libautocode/def/coords.d */ -int XY_to_LM (double *L, double *M, double x, double y, Coords *coords); -int LM_to_XY (double *x, double *y, double L, double M, Coords *coords); -int RD_to_LM (double *L, double *M, double ra, double dec, Coords *coords); -int LM_to_RD (double *ra, double *dec, double L, double M, Coords *coords); -int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords); -int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords); -int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords); -int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords); -int GetCoords (Coords *coords, Header *header); -int PutCoords (Coords *coords, Header *header); -void RegisterMosaic (Coords *coords); -void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); -OhanaProjection GetProjection (char *ctype); -int SetProjection (char *ctype, OhanaProjection proj); -OhanaProjectionMode GetProjectionMode (OhanaProjection proj); char *libdvo_version (void); @@ -913,4 +894,5 @@ Lensing *FtableToLensing (FTable *ftable, off_t *Nlensing, char *format); Lensobj *FtableToLensobj (FTable *ftable, off_t *Nlensobj, char *format); +StarPar *FtableToStarPar (FTable *ftable, off_t *Nstarpar, char *format); int FtableToImage (FTable *ftable, Header *theader, char *format); @@ -920,4 +902,5 @@ int LensingToFtable (FTable *ftable, Lensing *lensing, off_t Nlensing, char format); int LensobjToFtable (FTable *ftable, Lensobj *lensobj, off_t Nlensobj, char format); +int StarParToFtable (FTable *ftable, StarPar *starpar, off_t Nstarpar, char format); int ImageToFtable (FTable *ftable, Header *theader, char format); int ImageToVtable (VTable *vtable, Header *theader, char format); @@ -936,4 +919,5 @@ # include "ps1_v5_defs.h" # include "ps1_ref_defs.h" +# include "ps1_sim_defs.h" /*** DVO image db I/O Functions ***/ @@ -950,4 +934,8 @@ void dvo_image_create (FITS_DB *db, double ZeroPoint); +int gfits_table_set_Image (FTable *ftable); +int gfits_table_mkheader_Image (Header *header); +Image *gfits_table_get_Image (FTable *ftable, off_t *Ndata, char *swapped); + /* flatcorr APIs */ FlatCorrectionTable *FlatCorrectionLoad (char *filename, int VERBOSE); @@ -963,4 +951,5 @@ SkyTable *SkyTableLoadOptimal PROTO((char *catdir, char *SKYFILE, char *GSCFILE, int readwrite, int depth, int VERBOSE)); int SkyTableSetDepth PROTO((SkyTable *sky, int depth)); +SkyList *SkyRegionByIndex PROTO((SkyTable *table, int index)); SkyList *SkyRegionByCPT PROTO((SkyTable *table, char *filename)); SkyList *SkyRegionByPoint PROTO((SkyTable *table, int depth, double ra, double dec)); @@ -1009,4 +998,9 @@ int HostTableTestHost (SkyRegion *region, int hostID); +void InitHost (HostInfo *host); + +HostTableGroup *HostTableGroups (HostTable *table, int *ngroups); +int HostTableGroupWaitJobsGetIO (HostTableGroup *table, char *file, int lineno, int VERBOSE); + // functions to support tiny versions of Average and Measure void CopyAverageToTiny (AverageTiny *averageT, Average *average); @@ -1042,4 +1036,5 @@ void dvo_lensing_init (Lensing *lensing); void dvo_lensobj_init (Lensobj *lensobj, int toZero); +void dvo_starpar_init (StarPar *starpar); void InitRegionHosts (RegionHostInfo *hosts, int Nhosts, int NHOSTS); @@ -1051,3 +1046,40 @@ int RegionHostFindNeighbors (RegionHostTable *table, int Nhost); +// galaxy_model: +int TransformProperMotion_radians (double *uR, double *uD, double uL, double uB, double Rrad, double Drad, CoordTransform *transform); +int TransformProperMotion (double *uR, double *uD, double uL, double uB, double R, double D, CoordTransform *transform); +int SolarMotionModel_radians (double *uL_sol, double *uB_sol, double Lrad, double Brad, double distance); +int SolarMotionModel (double *uL_sol, double *uB_sol, double L, double B, double distance); +int GalaxyMotionModel_radians (double *uL_gal, double *uB_gal, double Lrad, double Brad); +int GalaxyMotionModel (double *uL_gal, double *uB_gal, double L, double B); +int InitGalaxyModel (char *version); + +# define LENSFIELD(NAME) float LensValue_##NAME (PhotCode *code, Lensobj *lensobj); + +LENSFIELD(X11_sm_obj) +LENSFIELD(X12_sm_obj) +LENSFIELD(X22_sm_obj) +LENSFIELD(E1_sm_obj) +LENSFIELD(E2_sm_obj) + +LENSFIELD(X11_sh_obj) +LENSFIELD(X12_sh_obj) +LENSFIELD(X22_sh_obj) +LENSFIELD(E1_sh_obj) +LENSFIELD(E2_sh_obj) + +LENSFIELD(X11_sm_psf) +LENSFIELD(X12_sm_psf) +LENSFIELD(X22_sm_psf) +LENSFIELD(E1_sm_psf) +LENSFIELD(E2_sm_psf) + +LENSFIELD(X11_sh_psf) +LENSFIELD(X12_sh_psf) +LENSFIELD(X22_sh_psf) +LENSFIELD(E1_sh_psf) +LENSFIELD(E2_sh_psf) + +# undef LENSFIELD + # endif // DVO_H Index: trunk/Ohana/src/libdvo/include/dvodb.h =================================================================== --- trunk/Ohana/src/libdvo/include/dvodb.h (revision 37729) +++ trunk/Ohana/src/libdvo/include/dvodb.h (revision 37807) @@ -55,4 +55,38 @@ MAG_OPTION_UC_DIST, MAG_OPTION_FLAGS, + + MAG_OPTION_X11_SM_OBJ, + MAG_OPTION_X12_SM_OBJ, + MAG_OPTION_X22_SM_OBJ, + MAG_OPTION_E1_SM_OBJ, + MAG_OPTION_E2_SM_OBJ, + + MAG_OPTION_X11_SH_OBJ, + MAG_OPTION_X12_SH_OBJ, + MAG_OPTION_X22_SH_OBJ, + MAG_OPTION_E1_SH_OBJ, + MAG_OPTION_E2_SH_OBJ, + + MAG_OPTION_X11_SM_PSF, + MAG_OPTION_X12_SM_PSF, + MAG_OPTION_X22_SM_PSF, + MAG_OPTION_E1_SM_PSF, + MAG_OPTION_E2_SM_PSF, + + MAG_OPTION_X11_SH_PSF, + MAG_OPTION_X12_SH_PSF, + MAG_OPTION_X22_SH_PSF, + MAG_OPTION_E1_SH_PSF, + MAG_OPTION_E2_SH_PSF, + + MAG_OPTION_F_AP_R5, + MAG_OPTION_F_ERR_AP_R5, + MAG_OPTION_F_STDEV_AP_R5, + MAG_OPTION_F_FILL_AP_R5, + + MAG_OPTION_F_AP_R6, + MAG_OPTION_F_ERR_AP_R6, + MAG_OPTION_F_STDEV_AP_R6, + MAG_OPTION_F_FILL_AP_R6, } dvoMagOptionType; @@ -139,4 +173,10 @@ MEAS_XFIX, MEAS_YFIX, + MEAS_XOFF_KH, + MEAS_YOFF_KH, + MEAS_XOFF_DCR, + MEAS_YOFF_DCR, + MEAS_ROFF_GAL, + MEAS_DOFF_GAL, MEAS_POS_SYS_ERR, MEAS_XFIELD, @@ -162,5 +202,38 @@ MEAS_FLAT, MEAS_CENTER_OFFSET, - MEAS_REF_COLOR, + MEAS_REF_COLOR_BLUE, + MEAS_REF_COLOR_RED, + MEAS_X11_SM_OBJ, + MEAS_X12_SM_OBJ, + MEAS_X22_SM_OBJ, + MEAS_E1_SM_OBJ, + MEAS_E2_SM_OBJ, + MEAS_X11_SH_OBJ, + MEAS_X12_SH_OBJ, + MEAS_X22_SH_OBJ, + MEAS_E1_SH_OBJ, + MEAS_E2_SH_OBJ, + MEAS_X11_SM_PSF, + MEAS_X12_SM_PSF, + MEAS_X22_SM_PSF, + MEAS_E1_SM_PSF, + MEAS_E2_SM_PSF, + MEAS_X11_SH_PSF, + MEAS_X12_SH_PSF, + MEAS_X22_SH_PSF, + MEAS_E1_SH_PSF, + MEAS_E2_SH_PSF, + MEAS_E_BV, // extinction (mags) + MEAS_E_BV_ERR, + MEAS_DISTANCE_MOD, // distance modulus (mags) + MEAS_DISTANCE_MOD_ERR, + MEAS_M_R, // absolute mag in r-band + MEAS_M_R_ERR, + MEAS_FEH, // metallicity + MEAS_FEH_ERR, + MEAS_URA_GALMODEL, // model pm prediction + MEAS_UDEC_GALMODEL, // model pm prediction + MEAS_RA_GALMODEL, // model pm prediction + MEAS_DEC_GALMODEL, // model pm prediction }; @@ -200,6 +273,19 @@ AVE_EXTID_HI, AVE_EXTID_LO, - AVE_REF_COLOR, + AVE_REF_COLOR_BLUE, + AVE_REF_COLOR_RED, AVE_PHOT, // photometry class of values + AVE_E_BV, // extinction (mags) + AVE_E_BV_ERR, + AVE_DISTANCE_MOD, // distance modulus (mags) + AVE_DISTANCE_MOD_ERR, + AVE_M_R, // absolute mag in r-band + AVE_M_R_ERR, + AVE_FEH, // metallicity + AVE_FEH_ERR, + AVE_URA_GALMODEL, // model pm prediction + AVE_UDEC_GALMODEL, // model pm prediction + AVE_RA_GALMODEL, // model pm prediction + AVE_DEC_GALMODEL, // model pm prediction }; @@ -276,5 +362,6 @@ IMAGE_NLINK_PHOTOM, IMAGE_NLINK_ASTROM, - IMAGE_REF_COLOR + IMAGE_REF_COLOR_BLUE, + IMAGE_REF_COLOR_RED }; @@ -368,6 +455,6 @@ int ParseImageField PROTO((dbField *field, char *fieldName)); -dbValue dbExtractAverages PROTO((Average *average, SecFilt *secfilt, Measure *measure, dbField *field)); -dbValue dbExtractMeasures PROTO((Average *average, SecFilt *secfilt, Measure *measure, dbField *field)); +dbValue dbExtractAverages PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, dbField *field)); +dbValue dbExtractMeasures PROTO((Average *average, SecFilt *secfilt, Measure *measure, Lensing *lensing, StarPar *starpar, dbField *field)); dbValue dbExtractImages PROTO((Image *image, off_t Nimage, off_t N, dbField *field)); @@ -376,4 +463,9 @@ int dbAstroRegionLimits PROTO((dbStack **stack, int *nstack, SkyRegionSelection *selection, int table)); char *strfloat PROTO((float value)); + +int dbFieldNeedMeasure (dbField *fields, int Nfields); +int dbFieldNeedLensobj (dbField *fields, int Nfields); +int dbFieldNeedLensing (dbField *fields, int Nfields); +int dbFieldNeedStarpar (dbField *fields, int Nfields, int isAverage); void FreeImageSelection (void); Index: trunk/Ohana/src/libdvo/include/libdvo_astro.h =================================================================== --- trunk/Ohana/src/libdvo/include/libdvo_astro.h (revision 37729) +++ trunk/Ohana/src/libdvo/include/libdvo_astro.h (revision 37807) @@ -1,16 +1,37 @@ -#ifdef notdef -# include "external.h" -# include "shell.h" -# include "dvomath.h" -# include "convert.h" -# include "display.h" -# include "data.h" -#endif +# ifndef LIBDVO_ASTRO_H +# define LIBDVO_ASTRO_H -# ifndef ASTRO_H -# define ASTRO_H +typedef enum { + COORD_NONE, + COORD_CELESTIAL, + COORD_GALACTIC, // based on Liu et al 2011 + COORD_ECLIPTIC, + COORD_GALACTIC_REID_2004, // an older definition of the galactic plane +} CoordTransformSystem; -void InitAstro (void); -void FreeAstro (void); +typedef enum { + PROJ_NONE, // undefined + PROJ_ZEA, // zenithal + PROJ_ZPL, // zenithal + PROJ_ARC, // zenithal + PROJ_STG, // zenithal + PROJ_SIN, // zenithal + PROJ_TAN, // zenithal + PROJ_TNX, // zenithal + PROJ_DIS, // zenithal (TAN + polyterms) + PROJ_LIN, // cartesian + PROJ_PLY, // cartesian (allow polyterms) + PROJ_WRP, // cartesian (allow polyterms, require mosaic) + PROJ_AIT, // pseudocyl + PROJ_GLS, // pseudocyl + PROJ_PAR, // pseudocyl +} OhanaProjection; + +typedef enum { + PROJ_MODE_NONE, + PROJ_MODE_CARTESIAN, + PROJ_MODE_ZENITHAL, + PROJ_MODE_PSEUDOCYL, +} OhanaProjectionMode; typedef struct { @@ -29,8 +50,136 @@ } CoordTransform; -typedef enum {COORD_NONE, COORD_CELESTIAL, COORD_GALACTIC, COORD_ECLIPTIC} CoordTransformSystem; +typedef struct { + int Nx; + int Ny; + float dX; + float dY; + float **dXv; + float **dYv; + unsigned int tableID; + unsigned int imageID; + int keep; +} AstromOffsetMap; + +typedef struct { + int Nmap; + AstromOffsetMap **map; + int *imageIDtoTableSeq; + int MaxImageID; + int MaxTableID; +} AstromOffsetTable; + +/* Internal version of Coords (see CoordsDisk in libautocode/def/coords-disk.d) */ +typedef struct Coords { + double crval1; // coordinate at reference pixel + double crval2; // coordinate at reference pixel + float crpix1; // coordinate of reference pixel + float crpix2; // coordinate of reference pixel + float cdelt1; // degrees per pixel + float cdelt2; // degrees per pixel + float pc1_1; // rotation matrix + float pc1_2; // rotation matrix + float pc2_1; // rotation matrix + float pc2_2; // rotation matrix + float polyterms[7][2]; // higher order warping terms + char ctype[15]; // coordinate type + char Npolyterms; // order of polynomial + struct Coords *mosaic; // pointer to parent mosaic + AstromOffsetMap *offsetMap; // pointer to offset map transformation +} Coords; + +typedef struct Image { + Coords coords; // astrometric data + e_time tzero; // readout time (row 0) + unsigned int nstar; // number of stars on image + float secz; // airmass (mag) + unsigned short NX; // image width + unsigned short NY; // image height + float apmifit; // aperture correction (mag) + float dapmifit; // apmifit error (mag) + float Mcal; // calibration mag (mag) + float dMcal; // error on Mcal (mag) + short Xm; // image chisq (10*log(value)) + short photcode; // identifier for CCD, + float exptime; // exposure time (seconds) + float sidtime; // sidereal time of exposure + float latitude; // observatory latitude (degrees) + float RAo; // image center (degrees) + float DECo; // image center (degrees) + float Radius; // image radius (degrees) + float refColorBlue; // median astrometry ref color + float refColorRed; // median astrometry ref color + char name[117]; // name of original image + unsigned char detection_limit; // detection limit (10*mag) + unsigned char saturation_limit; // saturation limit (10*mag) + unsigned char cerror; // astrometric error (50*arcsec) + unsigned char fwhm_x; // PSF x width (25*arcsec) + unsigned char fwhm_y; // PSF y width (25*arcsec) + unsigned char trate; // scan rate (100 usec/pixel) + unsigned char ccdnum; // CCD ID number + unsigned int flags; // image quality flags + unsigned int imageID; // internal image ID + unsigned int parentID; // associated ref image + unsigned int externID; // external image ID + unsigned short sourceID; // analysis source ID + short nLinkAstrom; // mean number of matched measurements for astrometry + short nLinkPhotom; // mean number of matched measurements for astrometry + short ubercalDist; // distance to nearest ubercal image + float dXpixSys; // systematic astrometry error in X + float dYpixSys; // systematic astrometry error in Y + float dMagSys; // systematic photometry error + unsigned short nFitAstrom; // number of stars used for astrometry cal + unsigned short nFitPhotom; // number of stars used for photometry cal + unsigned int photom_map_id; // reference to 2D zero point map + unsigned int astrom_map_id; // reference to 2D astrometry map + struct Image *parent; // pointer to parent mosaic (not save to disk) +} Image; CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output); int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform); +/* in coords.c */ +void InitCoords (Coords *coords, char *projection); +void CopyCoords (Coords *tgt, Coords *src); + +int XY_to_LM (double *L, double *M, double x, double y, Coords *coords); +int LM_to_XY (double *x, double *y, double L, double M, Coords *coords); +int RD_to_LM (double *L, double *M, double ra, double dec, Coords *coords); +int LM_to_RD (double *ra, double *dec, double L, double M, Coords *coords); +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords); +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords); +int fXY_to_RD (float *ra, float *dec, double x, double y, Coords *coords); +int fRD_to_XY (float *x, float *y, double ra, double dec, Coords *coords); +int GetCoords (Coords *coords, Header *header); +int PutCoords (Coords *coords, Header *header); +void coords_precess (double *ra, double *dec, double in_epoch, double out_epoch); +OhanaProjection GetProjection (char *ctype); +int SetProjection (char *ctype, OhanaProjection proj); +OhanaProjectionMode GetProjectionMode (OhanaProjection proj); + +/* in AstromOffsetMapIO.c */ +AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE); +int AstromOffsetMapSave (AstromOffsetTable *table, char *filename); +AstromOffsetTable *AstromOffsetMapToTable(AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap); +AstromOffsetMap_Disk_6x6 *AstromOffsetTableToMap(AstromOffsetTable *table, off_t *Nmap); + +/* in AstromOffsetMapOps.c */ +float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir); +int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); + +/* in AstromOffsetMapUtils.c */ +AstromOffsetMap *AstromOffsetMapInit (int Nx, int Ny); +void AstromOffsetMapFree (AstromOffsetMap *map); + +int AstromOffsetTableNewMap (AstromOffsetTable *table, int Nx, int Ny, Image *image); +int AstromOffsetTableMatchChips (Image *images, off_t Nimages, AstromOffsetTable *table); +AstromOffsetTable *AstromOffsetTableInit(); +void AstromOffsetTableFree(AstromOffsetTable *table); +int AstromOffsetTableAddMapFromImage (AstromOffsetTable *table, Image *image); +void AstromOffsetMapPrint (AstromOffsetMap *map, char *filename); +int AstromOffsetMapRepair (AstromOffsetMap *map, int xdir); +AstromOffsetMap *AstromOffsetMapCopy (AstromOffsetMap *map); +void AstromOffsetMapCopyData (AstromOffsetMap *tgt, AstromOffsetMap *src); +void AstromOffsetMapSetOrder (AstromOffsetMap *map, int Nx, int Ny, Image *image); + # endif Index: trunk/Ohana/src/libdvo/include/ps1_sim_defs.h =================================================================== --- trunk/Ohana/src/libdvo/include/ps1_sim_defs.h (revision 37807) +++ trunk/Ohana/src/libdvo/include/ps1_sim_defs.h (revision 37807) @@ -0,0 +1,14 @@ +Image *Image_PS1_SIM_ToInternal (Image_PS1_SIM *in, off_t Nvalues, off_t Nalloc); +Image_PS1_SIM *ImageInternalTo_PS1_SIM (Image *in, off_t Nvalues); + +Average *Average_PS1_SIM_ToInternal (Average_PS1_SIM *in, off_t Nvalues, SecFilt **primary); +Average_PS1_SIM *AverageInternalTo_PS1_SIM (Average *in, off_t Nvalues, SecFilt *primary); + +Measure *Measure_PS1_SIM_ToInternal (Average *ave, Measure_PS1_SIM *in, off_t Nvalues); +Measure_PS1_SIM *MeasureInternalTo_PS1_SIM (Average *ave, Measure *in, off_t Nvalues); + +SecFilt *SecFilt_PS1_SIM_ToInternal (SecFilt_PS1_SIM *in, off_t Nvalues); +SecFilt_PS1_SIM *SecFiltInternalTo_PS1_SIM (SecFilt *in, off_t Nvalues); + +StarPar *StarPar_PS1_SIM_ToInternal (StarPar_PS1_SIM *in, off_t Nvalues); +StarPar_PS1_SIM *StarParInternalTo_PS1_SIM (StarPar *in, off_t Nvalues); Index: trunk/Ohana/src/libdvo/include/ps1_v5_defs.h =================================================================== --- trunk/Ohana/src/libdvo/include/ps1_v5_defs.h (revision 37729) +++ trunk/Ohana/src/libdvo/include/ps1_v5_defs.h (revision 37807) @@ -9,8 +9,12 @@ Lensing *Lensing_PS1_V5_ToInternal (Lensing_PS1_V5 *in, off_t Nvalues); +Lensing *Lensing_PS1_V5_ALT_ToInternal (Lensing_PS1_V5_ALT *in, off_t Nvalues); Lensing_PS1_V5 *LensingInternalTo_PS1_V5 (Lensing *in, off_t Nvalues); Lensobj *Lensobj_PS1_V5_ToInternal (Lensobj_PS1_V5 *in, off_t Nvalues); Lensobj_PS1_V5 *LensobjInternalTo_PS1_V5 (Lensobj *in, off_t Nvalues); +StarPar *StarPar_PS1_V5_ToInternal (StarPar_PS1_V5 *in, off_t Nvalues); +StarPar_PS1_V5 *StarParInternalTo_PS1_V5 (StarPar *in, off_t Nvalues); + PhotCode *PhotCode_PS1_V5_To_Internal (PhotCode_PS1_V5 *in, off_t Nvalues); PhotCode_PS1_V5 *PhotCode_Internal_To_PS1_V5 (PhotCode *in, off_t Nvalues); Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapIO.c (revision 37807) @@ -0,0 +1,207 @@ +# include "dvo.h" + +AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) { + + off_t Nmap; + + Header header; + Matrix matrix; + Header theader; + FTable ftable; + + FILE *f = fopen (filename, "r"); + if (f == NULL) { + if (VERBOSE) fprintf (stderr, "can't find Astrom Offset Map file %s\n", filename); + return (NULL); + } + + /* load in table data */ + ftable.header = &theader; + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map header\n"); + fclose (f); + return (NULL); + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map matrix\n"); + gfits_free_header (&header); + fclose (f); + return (NULL); + } + + // for now, we only have one flavor (6x6) + + if (!gfits_fread_ftable (f, &ftable, "ASTROM_OFFSET_MAP_DISK_6x6")) { + if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map table\n"); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + fclose (f); + return (NULL); + } + AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL); + if (!map_disk) { + fprintf (stderr, "ERROR: failed to read Astrom Offset Map\n"); + exit (2); + } + + // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps. + // Convert the disk array of maps to then internal format in a rich structure: + AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap); + + gfits_free_header (&theader); + gfits_free_matrix (&matrix); + gfits_free_table (&ftable); + + return (table); +} + +int AstromOffsetMapSave (AstromOffsetTable *table, char *filename) { + + off_t Nmap; + + Header header; + Matrix matrix; + Header theader; + FTable ftable; + + FILE *f; + + /* make phu header (no matrix needed) */ + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps. + // Convert the internal format in a rich structure into a disk array of maps. + AstromOffsetMap_Disk_6x6 *map_disk = AstromOffsetTableToMap (table, &Nmap); + + ftable.header = &theader; + gfits_table_set_AstromOffsetMap_Disk_6x6 (&ftable, map_disk, Nmap); + + f = fopen (filename, "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "cannot open %s for output\n", filename); + return (FALSE); + } + + gfits_fwrite_header (f, &header); + gfits_fwrite_matrix (f, &matrix); + gfits_fwrite_Theader (f, &theader); + gfits_fwrite_table (f, &ftable); + fclose (f); + + return (TRUE); +} + +AstromOffsetTable *AstromOffsetMapToTable(AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap) { + + int i, j, k; + + AstromOffsetTable *table = NULL; + ALLOCATE (table, AstromOffsetTable, 1); + + table->Nmap = Nmap; + ALLOCATE (table->map, AstromOffsetMap *, Nmap); + + // find the max value of imageID + int MaxTableID = 0; + int MaxImageID = 0; + for (i = 0; i < Nmap; i++) { + MaxTableID = MAX(map_disk[i].tableID, MaxTableID); + MaxImageID = MAX(map_disk[i].imageID, MaxImageID); + } + table->MaxTableID = MaxTableID; + table->MaxImageID = MaxImageID; + + // generate the index and init values to -1 + ALLOCATE (table->imageIDtoTableSeq, int, MaxImageID + 1); + for (i = 0; i <= MaxImageID; i++) { + table->imageIDtoTableSeq[i] = -1; + } + + // assign the ID values + for (i = 0; i < Nmap; i++) { + int ImageID = map_disk[i].imageID; + myAssert (table->imageIDtoTableSeq[ImageID] == -1, "oops, duplicate image IDs"); + table->imageIDtoTableSeq[ImageID] = i; + } + + // assign the map values (this allocates just the area needed for each image, not the + // full 6x6, saving some memory while doing the analysis) + for (i = 0; i < Nmap; i++) { + ALLOCATE (table->map[i], AstromOffsetMap, 1); + + table->map[i][0].Nx = map_disk[i].Nx; + table->map[i][0].Ny = map_disk[i].Ny; + table->map[i][0].tableID = map_disk[i].tableID; + table->map[i][0].imageID = map_disk[i].imageID; + + // STORE THESE VALUES? + table->map[i][0].dY = map_disk[i].dX; + table->map[i][0].dX = map_disk[i].dY; + + // since this was on disk, we obviously keep it + table->map[i][0].keep = TRUE; + + ALLOCATE (table->map[i][0].dXv, float *, map_disk[i].Nx); + ALLOCATE (table->map[i][0].dYv, float *, map_disk[i].Nx); + + for (j = 0; j < map_disk[i].Nx; j++) { + ALLOCATE (table->map[i][0].dXv[j], float, map_disk[i].Ny); + ALLOCATE (table->map[i][0].dYv[j], float, map_disk[i].Ny); + + for (k = 0; k < map_disk[i].Ny; k++) { + table->map[i][0].dXv[j][k] = map_disk[i].dXv[j][k]; + table->map[i][0].dYv[j][k] = map_disk[i].dYv[j][k]; + } + } + } + return table; +} + +AstromOffsetMap_Disk_6x6 *AstromOffsetTableToMap(AstromOffsetTable *table, off_t *Nmap) { + + int i, j, k; + + AstromOffsetMap_Disk_6x6 *map_disk = NULL; + ALLOCATE (map_disk, AstromOffsetMap_Disk_6x6, table->Nmap); + + // assign the map values (this allocates just the area needed for each image, not the + // full 6x6, saving some memory while doing the analysis) + // some maps in the table should be skipped because their image no longer uses them + int Ndisk = 0; + for (i = 0; i < table->Nmap; i++) { + if (!table->map[i][0].keep) continue; + map_disk[Ndisk].Nx = table->map[i][0].Nx; + map_disk[Ndisk].Ny = table->map[i][0].Ny; + map_disk[Ndisk].tableID = table->map[i][0].tableID; + map_disk[Ndisk].imageID = table->map[i][0].imageID; + + map_disk[Ndisk].dX = table->map[i][0].dX; + map_disk[Ndisk].dY = table->map[i][0].dY; + + + + for (j = 0; j < map_disk[Ndisk].Nx; j++) { + for (k = 0; k < map_disk[Ndisk].Ny; k++) { + map_disk[Ndisk].dXv[j][k] = table->map[i][0].dXv[j][k]; + map_disk[Ndisk].dYv[j][k] = table->map[i][0].dYv[j][k]; + } + for (k = map_disk[Ndisk].Ny; k < 6; k++) { + map_disk[Ndisk].dXv[j][k] = 0.0; + map_disk[Ndisk].dYv[j][k] = 0.0; + } + } + for (j = map_disk[Ndisk].Nx; j < 6; j++) { + for (k = 0; k < 6; k++) { + map_disk[Ndisk].dXv[j][k] = 0.0; + map_disk[Ndisk].dYv[j][k] = 0.0; + } + } + Ndisk ++; + } + *Nmap = Ndisk; + return map_disk; +} + Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapOps.c (revision 37807) @@ -0,0 +1,543 @@ +# include "dvo.h" + +/* AstromOffsetMap functions: + */ + +int dump_map_data (float *x, float *y, float *f, int Npts, char *filename); +int AstromOffsetMapFit_Chisq (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); +int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir); + +float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir) { + + // given x,y find the offset (x-direction) at that location from the map + + float **V = xdir ? map->dXv : map->dYv; + + // if there is no spatial information, the value at the point is the value in the map + if ((map->Nx == 1) && (map->Ny == 1)) return V[0][0]; + + // if there is no spatial information in 1D, the interpolation is 1D + if (map->Nx == 1) { + float ymo = y * map->dY - 0.5; + int ymi = floor(ymo); + ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 (Ny must be > 1) + + float ymf = ymo - ymi; + + float value = ymf * V[0][ymi+1] + (1.0 - ymf) * V[0][ymi]; + return value; + } + + // if there is no spatial information in 1D, the interpolation is 1D + if (map->Ny == 1) { + float xmo = x * map->dX - 0.5; + int xmi = floor(xmo); + xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-2 (Nx must be > 1) + + float xmf = xmo - xmi; + + float value = xmf * V[xmi+1][0] + (1.0 - xmf) * V[xmi][0]; + return value; + } + + // x & y are in Big Image coordinates, convert to fractional map coordinates + + float xmo = x * map->dX - 0.5; + int xmi = floor(xmo); + xmi = MAX(0,MIN(map->Nx - 2, xmi)); // force range of ymi to be 0,Nx-1 + float xmf = xmo - xmi; + + float ymo = y * map->dY - 0.5; + int ymi = floor(ymo); + ymi = MAX(0,MIN(map->Ny - 2, ymi)); // force range of ymi to be 0,Ny-2 + float ymf = ymo - ymi; + + float V00 = V[xmi+0][ymi+0]; + float V01 = V[xmi+0][ymi+1]; + float V10 = V[xmi+1][ymi+0]; + float V11 = V[xmi+1][ymi+1]; + + float Vx0 = V10*xmf + V00*(1.0 - xmf); + float Vx1 = V11*xmf + V01*(1.0 - xmf); + + float value = ymf * Vx1 + (1.0 - ymf) * Vx0; + return value; +} + +int AstromOffsetMapFit (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) { + + int status = AstromOffsetMapFit_Median (map, x, y, f, Npts, xdir); + return status; +} + +// given (x,y),value vector sets, choose the map that minimizes the difference between +// the values and bilinear interpolation of the map +int AstromOffsetMapFit_Chisq (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) { + + int i, j, ix, iy, jx, jy; + + // choose to map direction: + float **V = xdir ? map->dXv : map->dYv; + + // special cases: + if ((map->Nx == 1) && (map->Ny == 1)) { + + // no spatial information, just return (clipped?) median or mean or something + VStatsType stats; + stats.statmode = VSTATS_INNER_MEAN; + + vstats_getstats_f (f, NULL, NULL, Npts, &stats); + V[0][0] = stats.mean; + return TRUE; + } + + // special cases: + if (map->Nx == 1) { + fprintf (stderr, "1 x N not yet coded\n"); + exit (1); + return TRUE; + } + + if (map->Ny == 1) { + fprintf (stderr, "N x 1 not yet coded\n"); + exit (1); + return TRUE; + } + + // this bit is copied nearly intact from psImageMapFit.c + + // set up the redirection table so we can use sA[-1][-1], etc + float SAm[3][3], *SAv[3], **sA; + + for (i = 0; i < 3; i++) { + SAv[i] = SAm[i] + 1; + } + sA = SAv + 1; + + int Nx = map->Nx; + int Ny = map->Ny; + + double **A, **B; + ALLOCATE (A, double *, Nx*Ny); + ALLOCATE (B, double *, Nx*Ny); + for (i = 0; i < Nx*Ny; i++) { + ALLOCATE (A[i], double, Nx*Ny); + memset (A[i], 0, sizeof(double)*Nx*Ny); + ALLOCATE (B[i], double, 1); + memset (B[i], 0, sizeof(double)); + } + + // we are looping over the Nx,Ny image map elements; + // the matrix equation contains Nx*Ny rows and columns + // for (int n = 1; n < Nx - 1; n++) { + // for (int m = 1; m < Ny - 1; m++) { + + if (0) { + FILE *fd = fopen ("stats.dump.txt", "w"); + for (i = 0; i < Npts; i++) { + fprintf (fd, "%d %f %f %f\n", i, x[i], y[i], f[i]); + } + fclose(fd); + } + + // float Total = 0.0; + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + // define & init summing variables + double rx_rx_ry_ry = 0; + double rx_rx_dy_ry = 0; + double dx_rx_ry_ry = 0; + double dx_rx_dy_ry = 0; + double fi_rx_ry = 0; + double rx_rx_py_py = 0; + double rx_rx_qy_py = 0; + double dx_rx_py_py = 0; + double dx_rx_qy_py = 0; + double fi_rx_py = 0; + double px_px_ry_ry = 0; + double px_px_dy_ry = 0; + double qx_px_ry_ry = 0; + double qx_px_dy_ry = 0; + double fi_px_ry = 0; + double px_px_py_py = 0; + double px_px_qy_py = 0; + double qx_px_py_py = 0; + double qx_px_qy_py = 0; + double fi_px_py = 0; + + // generate the sums for the fitting matrix element I,J + // I = n + nX*m + // J = (n + jn) + nX*(m + jm) + for (i = 0; i < Npts; i++) { + + // data value & weight for this point + if (!isfinite(f[i])) continue; + + // if (mask && (mask[i] & maskValue)) continue; + + // base coordinate offset for this point (x,y) relative to this map element (n,m) + // double dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5); + // double dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5); + + double dx = x[i] * map->dX - ix - 0.5; + double dy = y[i] * map->dY - iy - 0.5; + + // XXX do I need to do something different at the edge or not? I want to allow + // x[i] to extend beyond the grid + +# if (0) + // edge cases to include: + bool edgeX = false; + edgeX |= ((n == 1) && (dx < -1.0)); + edgeX |= ((n == Nx - 2) && (dx > +1.0)); + + bool edgeY = false; + edgeY |= ((m == 1) && (dy < -1.0)); + edgeY |= ((m == Ny - 2) && (dy > +1.0)); + + // skip points outside of 2x2 grid centered on n,m: + if (!edgeX && (fabs(dx) > 1.0)) continue; + if (!edgeY && (fabs(dy) > 1.0)) continue; +# endif + + if (fabs(dx) > 1.0) continue; + if (fabs(dy) > 1.0) continue; + + // related offset values + double rx = 1.0 - dx; + double ry = 1.0 - dy; + double px = 1.0 + dx; + double py = 1.0 + dy; + double qx = 0.0 - dx; + double qy = 0.0 - dy; + + // sum the appropriate elements for the different quadrants + int Qx = (dx >= 0) ? 1 : 0; + int Qy = (dy >= 0) ? 1 : 0; + + // XXX why change the selection at the edges? + // if (n == 0) Qx = 1; + // if (n == Nx - 1) Qx = 0; + // if (m == 0) Qy = 1; + // if (m == Ny - 1) Qy = 0; + + float fi = f[i]; + assert (isfinite(fi)); + assert (isfinite(rx)); + assert (isfinite(ry)); + + // points at offset 1,1 + if ((Qx == 1) && (Qy == 1)) { + rx_rx_ry_ry += rx*rx*ry*ry; + rx_rx_dy_ry += rx*rx*dy*ry; + dx_rx_ry_ry += dx*rx*ry*ry; + dx_rx_dy_ry += dx*rx*dy*ry; + fi_rx_ry += fi*rx*ry; + } + // points at offset 1,0 + if ((Qx == 1) && (Qy == 0)) { + rx_rx_py_py += rx*rx*py*py; + rx_rx_qy_py += rx*rx*qy*py; + dx_rx_py_py += dx*rx*py*py; + dx_rx_qy_py += dx*rx*qy*py; + fi_rx_py += fi*rx*py; + } + // points at offset 0,1 + if ((Qx == 0) && (Qy == 1)) { + px_px_ry_ry += px*px*ry*ry; + px_px_dy_ry += px*px*dy*ry; + qx_px_ry_ry += qx*px*ry*ry; + qx_px_dy_ry += qx*px*dy*ry; + fi_px_ry += fi*px*ry; + } + // points at offset 0,0 + if ((Qx == 0) && (Qy == 0)) { + px_px_py_py += px*px*py*py; + px_px_qy_py += px*px*qy*py; + qx_px_py_py += qx*px*py*py; + qx_px_qy_py += qx*px*qy*py; + fi_px_py += fi*px*py; + } + } + + // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm), + // jn,jm = -1 to +1. Convert the sums above into the correct coefficients + sA[-1][-1] = qx_px_qy_py; + sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py; + sA[-1][+1] = qx_px_dy_ry; + sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py; + sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py; + sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry; + sA[+1][-1] = dx_rx_qy_py; + sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py; + sA[+1][+1] = dx_rx_dy_ry; + + // I[ 0][ 0] = index for this n,m element: + int I = ix + Nx * iy; + B[I][0] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py; + + // insert these values into their corresponding locations in A, B + // float Sum = 0.0; + for (jx = -1; jx <= +1; jx++) { + if (ix + jx < 0) continue; + if (ix + jx >= Nx) continue; + for (jy = -1; jy <= +1; jy++) { + if (iy + jy < 0) continue; + if (iy + jy >= Ny) continue; + int J = (ix + jx) + Nx * (iy + jy); + A[J][I] = sA[jx][jy]; + if (abs(A[J][I]) > 1000) { + fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, ix, iy, ix + jx, iy + jy, sA[jx][jy]); + } + // Sum += sA[jn][jm]; + } + } + // fprintf (stderr, "B %d (%d %d) : %f : %f\n", I, n, m, B->data.F32[I], Sum); + // Total += Sum; + } + } + // fprintf (stderr, "Total: %f\n", Total); + + // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0 + int *Empty = NULL; + ALLOCATE (Empty, int, Nx*Ny); + + // find the max value + double MaxPivot = 0.0; + for (i = 0; i < Nx*Ny; i++) { + MaxPivot = MAX(A[i][i], MaxPivot); + } + + // elements of A[i][j] are just the fractional weight of the number of input points to the pixel. + // any pixels which have < 5% contribution are not very well constrained. NAN them for now... + double MinPivot = 0.05 * MaxPivot; + for (i = 0; i < Nx*Ny; i++) { + Empty[i] = 0; + if (fabs(A[i][i]) < MinPivot) { + Empty[i] = 1; + for (j = 0; j < Nx*Ny; j++) { + A[i][j] = 0.0; + A[j][i] = 0.0; + } + A[i][i] = 1.0; + B[i][0] = 0.0; + } + } + + if (1) { + FILE *fd = fopen ("matrix.dat", "w"); + for (i = 0; i < Nx*Ny; i++) { + for (j = 0; j < Nx*Ny; j++) { + fprintf (fd, "%10.4f ", A[i][j]); + } + fprintf (fd, " : %10.4f\n", B[i][0]); + } + fclose (fd); + } + + if (!dgaussjordan(A, B, Nx*Ny, 1)) { + fprintf (stderr, "FAIL\n"); + exit (1); + } + + // set bad values to NaN + for (i = 0; i < Nx*Ny; i++) { + if (Empty[i]) { + B[i][0] = NAN; + A[i][i] = 0; + } + } + + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + int I = ix + Nx * iy; + V[ix][iy] = B[I][0]; + // error[ix][iy] = sqrt(A[I][I]); + } + } + + // XXX free someone? + for (i = 0; i < Nx*Ny; i++) { + free (A[i]); + free (B[i]); + } + free (A); + free (B); + free (Empty); + + return TRUE; +} + +// given (x,y),value vector sets, choose the map that minimizes the difference between +// the values and bilinear interpolation of the map +int AstromOffsetMapFit_Median (AstromOffsetMap *map, float *x, float *y, float *f, int Npts, int xdir) { + + int i, ix, iy; + + // choose to map direction: + float **V = xdir ? map->dXv : map->dYv; + + // measure clipped median in each bin + VStatsType stats; + stats.statmode = VSTATS_INNER_MEAN; + + int Nx = map->Nx; + int Ny = map->Ny; + int Npix = Nx*Ny; + + double **values; + int *Nvalue; + ALLOCATE (Nvalue, int, Npix); + ALLOCATE (values, double *, Npix); + for (i = 0; i < Npix; i++) { + ALLOCATE (values[i], double, Npts); + Nvalue[i] = 0; + } + + // assign the points to the map cells + for (i = 0; i < Npts; i++) { + + // data value & weight for this point + if (!isfinite(f[i])) continue; + + // if (mask && (mask[i] & maskValue)) continue; + + // base coordinate offset for this point (x,y) relative to this map element (n,m) + // double dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5); + // double dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5); + + // bin for this point + ix = MAX(0, MIN(Nx-1, (int)(x[i] * map->dX))); + iy = MAX(0, MIN(Ny-1, (int)(y[i] * map->dY))); + int I = ix + Nx * iy; + + int N = Nvalue[I]; + values[I][N] = f[i]; + Nvalue[I] ++; + } + + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + int I = ix + Nx * iy; + if (Nvalue[I] > 5) { + vstats_getstats (values[I], NULL, NULL, Nvalue[I], &stats); + V[ix][iy] = stats.mean; + } else { + V[ix][iy] = NAN; + } + } + } + + // XXX free someone? + for (i = 0; i < Npix; i++) { + free (values[i]); + } + free (values); + free (Nvalue); + + return TRUE; +} + +// this function repairs an image with NAN pixels (only valid for a small-scale map -- no robust mean) +int AstromOffsetMapRepair (AstromOffsetMap *map, int xdir) { + + // we are going to repair the image by: + // 1) finding NAN pixels + // 2) if any of the neighbors are valid, + // replace with the mean of the neighbors + // 3) otherwise, replace with the image mean + + int Nx = map->Nx; + int Ny = map->Ny; + + int ix, iy; + + float **Vfix = NULL; + ALLOCATE (Vfix, float *, Nx); + for (ix = 0; ix < Nx; ix++) { + ALLOCATE (Vfix[ix], float, Ny); + for (iy = 0; iy < Ny; iy++) { + Vfix[ix][iy] = NAN; + } + } + + // find the global mean + float mean = 0.0; + float npix = 0.0; + + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy]; + if (isnan(value)) continue; + mean += value; + npix += 1.0; + } + } + mean /= npix; + + // find the NAN pixels: + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + + float value = xdir ? map->dXv[ix][iy] : map->dYv[ix][iy]; + if (!isnan(value)) { + Vfix[ix][iy] = value; + continue; + } + + // find mean of all possible neighbors + float meanLocal = 0.0; + float npixLocal = 0.0; + + int jx, jy; + for (jx = -1; jx <= +1; jx++) { + int nx = ix + jx; + if (nx < 0) continue; + if (nx >= Nx) continue; + for (jy = -1; jy <= +1; jy++) { + int ny = iy + jy; + if (ny < 0) continue; + if (ny >= Ny) continue; + float value = xdir ? map->dXv[nx][ny] : map->dYv[nx][ny]; + if (isnan(value)) continue; + + meanLocal += value; + npixLocal += 1.0; + } + } + // if there are no valid (non-NAN) local pixels, use global mean + meanLocal = (npixLocal > 0.0) ? meanLocal / npixLocal : mean; + Vfix[ix][iy] = meanLocal; + } + } + + // replace the bad pixels: + for (ix = 0; ix < Ny; ix++) { + for (iy = 0; iy < Ny; iy++) { + if (xdir) { + map->dXv[ix][iy] = Vfix[ix][iy]; + } else { + map->dYv[ix][iy] = Vfix[ix][iy]; + } + } + } + for (ix = 0; ix < Nx; ix++) { + free (Vfix[ix]); + } + free (Vfix); + + return TRUE; +} + +int dump_map_data (float *x, float *y, float *f, int Npts, char *filename) { + + FILE *fout = fopen (filename, "w"); + + int i; + for (i = 0; i < Npts; i++) { + fprintf (fout, "%d %f %f %f\n", i, x[i], y[i], f[i]); + } + fclose (fout); + return TRUE; +} + Index: trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c =================================================================== --- trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/AstromOffsetMapUtils.c (revision 37807) @@ -0,0 +1,280 @@ +# include + +int AstromOffsetTableMatchChips (Image *images, off_t Nimages, AstromOffsetTable *table) { + + // we have a table of astrometry offset maps. we want to find the chips that match to + // each of the maps so we can assign the map to the image.coords.offsetMap entry + + // we have the lookup table of table->imageIDtoTableSeq to find the match by imageID + + // we are just going to assume that table->imageIDtoTableSeq is correct and complete + off_t i; + for (i = 0; i < Nimages; i++) { + int imageID = images[i].imageID; + if (imageID < 0) continue; + if (imageID > table->MaxImageID) continue; + + int seq = table->imageIDtoTableSeq[imageID]; + if (seq < 0) continue; + if (seq >= table->Nmap) continue; // this one is probably not valid, right? + + images[i].coords.offsetMap = table->map[seq]; + } + return (TRUE); +} + +AstromOffsetMap *AstromOffsetMapInit (int Nx, int Ny) { + + AstromOffsetMap *map = NULL; + ALLOCATE (map, AstromOffsetMap, 1); + + map->dX = NAN; + map->dY = NAN; + map->Nx = Nx; // output map size + map->Ny = Ny; // output map size + + map->keep = TRUE; // output map size + + ALLOCATE (map->dXv, float *, map->Nx); + ALLOCATE (map->dYv, float *, map->Nx); + + int j, k; + for (j = 0; j < map->Nx; j++) { + ALLOCATE (map->dXv[j], float, map->Ny); + ALLOCATE (map->dYv[j], float, map->Ny); + + for (k = 0; k < map->Ny; k++) { + map->dXv[j][k] = 0.0; + map->dYv[j][k] = 0.0; + } + } + + return map; +} + +void AstromOffsetMapFree (AstromOffsetMap *map) { + + if (!map) return; + + int j; + for (j = 0; j < map->Nx; j++) { + if (map->dXv) free (map->dXv[j]); + if (map->dYv) free (map->dYv[j]); + } + + free (map->dXv); + free (map->dYv); + + free (map); + + return; +} + +void AstromOffsetMapSetOrder (AstromOffsetMap *map, int Nx, int Ny, Image *image) { + + int j, k; + + // rather than try to figure out how to resize, and free/allocate, lets just free the old arrays and make new ones + for (j = 0; j < map->Nx; j++) { + if (map->dXv) free (map->dXv[j]); + if (map->dYv) free (map->dYv[j]); + } + free (map->dXv); + free (map->dYv); + + map->Nx = Nx; // output map size + map->Ny = Ny; // output map size + map->dX = Nx / (float) image->NX; + map->dY = Ny / (float) image->NY; + + map->keep = TRUE; // output map size + + ALLOCATE (map->dXv, float *, map->Nx); + ALLOCATE (map->dYv, float *, map->Nx); + + for (j = 0; j < map->Nx; j++) { + ALLOCATE (map->dXv[j], float, map->Ny); + ALLOCATE (map->dYv[j], float, map->Ny); + + for (k = 0; k < map->Ny; k++) { + map->dXv[j][k] = 0.0; + map->dYv[j][k] = 0.0; + } + } + return; +} + +AstromOffsetMap *AstromOffsetMapCopy (AstromOffsetMap *map) { + + if (!map) return NULL; + + AstromOffsetMap *tgt = AstromOffsetMapInit (map->Nx, map->Ny); + + tgt->dX = map->dX; + tgt->dY = map->dY; + tgt->imageID = map->imageID; + tgt->tableID = map->tableID; + + int j, k; + for (j = 0; j < tgt->Nx; j++) { + for (k = 0; k < map->Ny; k++) { + tgt->dXv[j][k] = map->dXv[j][k]; + tgt->dYv[j][k] = map->dYv[j][k]; + } + } + return tgt; +} + +// copy the data from one to another, assuming a pre-allocated structure +void AstromOffsetMapCopyData (AstromOffsetMap *tgt, AstromOffsetMap *src) { + + myAssert (tgt->Nx == src->Nx, "programming error"); + myAssert (tgt->Ny == src->Ny, "programming error"); + + tgt->dX = src->dX; + tgt->dY = src->dY; + tgt->imageID = src->imageID; + tgt->tableID = src->tableID; + + tgt->keep = src->keep; + + int j, k; + for (j = 0; j < tgt->Nx; j++) { + for (k = 0; k < src->Ny; k++) { + tgt->dXv[j][k] = src->dXv[j][k]; + tgt->dYv[j][k] = src->dYv[j][k]; + } + } + return; +} + +void AstromOffsetTableFree (AstromOffsetTable *table) { + + if (!table) return; + + int i, j; + for (i = 0; i < table->Nmap; i++) { + for (j = 0; j < table->map[i][0].Nx; j++) { + free (table->map[i][0].dXv[j]); + free (table->map[i][0].dYv[j]); + } + if (!table->map[i]) continue; + FREE (table->map[i][0].dXv); + FREE (table->map[i][0].dYv); + FREE (table->map[i]); + } + FREE (table->imageIDtoTableSeq); + FREE (table->map); +} + +int AstromOffsetTableNewMap (AstromOffsetTable *table, int Nx, int Ny, Image *image) { + + off_t i; + if (image->imageID > table->MaxImageID) { + int oldMaxID = table->MaxImageID; + table->MaxImageID = image->imageID; + REALLOCATE (table->imageIDtoTableSeq, int, table->MaxImageID + 1); + for (i = oldMaxID + 1; i < table->MaxImageID + 1; i++) { + table->imageIDtoTableSeq[i] = -1; + } + } + int seq = table->imageIDtoTableSeq[image->imageID]; + if (seq != -1) { + // already assigned (but not used) + AstromOffsetMapSetOrder (table->map[seq], Nx, Ny, image); + image[0].coords.offsetMap = table->map[seq]; + return TRUE; + } + + int Nmap = table->Nmap; + + table->Nmap++; + REALLOCATE (table->map, AstromOffsetMap *, table->Nmap); + + myAssert (table->imageIDtoTableSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision"); + table->imageIDtoTableSeq[image->imageID] = Nmap; + + table->map[Nmap] = AstromOffsetMapInit (Nx, Ny); + + table->MaxTableID ++; + table->map[Nmap][0].tableID = table->MaxTableID; + table->map[Nmap][0].imageID = image->imageID; + + table->map[Nmap][0].dX = Nx / (float) image->NX; + table->map[Nmap][0].dY = Ny / (float) image->NY; + + image[0].coords.offsetMap = table->map[Nmap]; + return TRUE; +} + +int AstromOffsetTableAddMapFromImage (AstromOffsetTable *table, Image *image) { + + int Nmap = table->Nmap; + + table->Nmap++; + REALLOCATE (table->map, AstromOffsetMap *, table->Nmap); + + // find the imageID and update the imageIDtoTableSeq allocation if needed + off_t i; + if (image->imageID > table->MaxImageID) { + int oldMaxID = table->MaxImageID; + table->MaxImageID = image->imageID; + REALLOCATE (table->imageIDtoTableSeq, int, table->MaxImageID + 1); + for (i = oldMaxID + 1; i < table->MaxImageID + 1; i++) { + table->imageIDtoTableSeq[i] = -1; + } + } + myAssert (table->imageIDtoTableSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision"); + table->imageIDtoTableSeq[image->imageID] = Nmap; + + table->map[Nmap] = image[0].coords.offsetMap; + return TRUE; +} + +AstromOffsetTable *AstromOffsetTableInit() { + + AstromOffsetTable *table = NULL; + ALLOCATE (table, AstromOffsetTable, 1); + + table->Nmap = 0; + ALLOCATE (table->map, AstromOffsetMap *, 1); + + table->MaxTableID = 0; + table->MaxImageID = 0; + + // generate the index and init values to -1 + ALLOCATE (table->imageIDtoTableSeq, int, 1); + table->imageIDtoTableSeq[0] = -1; + + return table; +} + +void AstromOffsetMapPrint (AstromOffsetMap *map, char *filename) { + + FILE *f = stderr; + if (strcasecmp(filename, "stderr") && strcasecmp(filename, "stdout")) { + f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "failed to open output file %s\n", filename); + } + } + + fprintf (f, "imageID: %d, tableID: %d (dX: %f, dY: %f), keep: %d\n", map->imageID, map->tableID, map->dX, map->dY, map->keep); + int ix, iy; + fprintf (f, "dXv map:\n"); + for (ix = 0; ix < map->Nx; ix++) { + for (iy = 0; iy < map->Ny; iy++) { + fprintf (f, "%9.5f ", map->dXv[ix][iy]); + } + fprintf (f, "\n"); + } + fprintf (f, "dYv map:\n"); + for (ix = 0; ix < map->Nx; ix++) { + for (iy = 0; iy < map->Ny; iy++) { + fprintf (f, "%9.5f ", map->dYv[ix][iy]); + } + fprintf (f, "\n"); + } + if (f != stderr) fclose (f); + return; +} Index: trunk/Ohana/src/libdvo/src/BoundaryTree.c =================================================================== --- trunk/Ohana/src/libdvo/src/BoundaryTree.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/BoundaryTree.c (revision 37807) @@ -371,4 +371,6 @@ int BoundaryTreeCellCoords (BoundaryTree *tree, int *zone, int *band, double ra, double dec) { + ra = ohana_normalize_angle (ra); + // first, find the containing zone Index: trunk/Ohana/src/libdvo/src/HostTable.c =================================================================== --- trunk/Ohana/src/libdvo/src/HostTable.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/HostTable.c (revision 37807) @@ -3,15 +3,20 @@ # include +void InitHost (HostInfo *host) { + host->hostname = NULL; + host->pathname = NULL; + host->results = NULL; + host->stdio[HOST_STDIN] = -1; + host->stdio[HOST_STDOUT] = -1; + host->stdio[HOST_STDERR] = -1; + host->pid = 0; + return; +} + void InitHosts (HostInfo *hosts, int Nhosts, int NHOSTS) { int i; for (i = Nhosts; i < NHOSTS; i++) { - hosts[i].hostname = NULL; - hosts[i].pathname = NULL; - hosts[i].results = NULL; - hosts[i].stdio[HOST_STDIN] = -1; - hosts[i].stdio[HOST_STDOUT] = -1; - hosts[i].stdio[HOST_STDERR] = -1; - hosts[i].pid = 0; + InitHost (&hosts[i]); } return; @@ -36,4 +41,15 @@ } if (table->index) free (table->index); + free (table); + return; +} + +// a HostTableGroup is a pointer set to one or more HostTable entries +void FreeTableGroup (HostTableGroup *table) { + + if (!table) return; + if (table->hosts) { + free (table->hosts); + } free (table); return; @@ -136,4 +152,66 @@ fclose (f); return table; +} + +// split a host table into Ngroups, each with a unique set of hosts +HostTableGroup *HostTableGroups (HostTable *table, int *ngroups) { + + // identify the unique host names and the number of times they each exist + // (Host tables are not very large; I can afford to do this in an inefficient way) + int i, j, k; + + int Nunique = 0; + char **unique_hosts = NULL; // an array of pointers to a set of the names + int *Nunique_hosts = NULL; // an array of pointers to the count of each + ALLOCATE ( unique_hosts, char *, table->Nhosts); + ALLOCATE (Nunique_hosts, int, table->Nhosts); + for (i = 0; i < table->Nhosts; i++) Nunique_hosts[i] = 0; + + for (i = 0; i < table->Nhosts; i++) { + int found = FALSE; + for (j = 0; !found && (j < Nunique); j++) { + if (!strcmp(table->hosts[i].hostname, unique_hosts[j])) { + Nunique_hosts[j] ++; + found = TRUE; + } + } + if (!found) { + unique_hosts[Nunique] = table->hosts[i].hostname; + Nunique_hosts[Nunique] = 1; + Nunique ++; + } + } + + int Ngroups = 0; + for (i = 0; i < Nunique; i++) { + Ngroups = MAX (Nunique_hosts[i], Ngroups); + } + + HostTableGroup *groups = NULL; + int *foundHost = NULL; + + ALLOCATE (foundHost, int, table->Nhosts); + ALLOCATE (groups, HostTableGroup, Ngroups); + + for (i = 0; i < table->Nhosts; i++) foundHost[i] = FALSE; + + // in each group, attempt to add one of each unique host + for (i = 0; i < Ngroups; i++) { + groups[i].Nhosts = 0; + ALLOCATE (groups[i].hosts, HostInfo *, Nunique); + for (j = 0; j < Nunique; j++) { + int found = FALSE; + for (k = 0; !found && (k < table->Nhosts); k++) { + if (foundHost[k]) continue; + if (strcmp(unique_hosts[j], table->hosts[k].hostname)) continue; + groups[i].hosts[groups[i].Nhosts] = &table->hosts[k]; + groups[i].Nhosts ++; + found = TRUE; + foundHost[k] = TRUE; + } + } + } + *ngroups = Ngroups; + return groups; } @@ -255,5 +333,5 @@ // this loop has 2 chunks: (a) check for I/O + (b) check for jobs done - while (1) { + while (Nfound < Nrunning) { // Wait up to 0.5 second for host to provide I/O @@ -284,5 +362,5 @@ // if (table->hosts[i].stdout.Nbuffer > 0x10000) { if (table->hosts[i].stdout.Nbuffer > 0x1000) { - int printHead = VERBOSE || (table->hosts[i].stdout.Nbuffer > 0); + int printHead = VERBOSE && (table->hosts[i].stdout.Nbuffer > 0); if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stdout.Nbuffer); int Nout = write (STDOUT_FILENO, table->hosts[i].stdout.buffer, table->hosts[i].stdout.Nbuffer); @@ -298,5 +376,5 @@ // if (table->hosts[i].stderr.Nbuffer > 0x10000) { if (table->hosts[i].stderr.Nbuffer > 0x1000) { - int printHead = VERBOSE || (table->hosts[i].stderr.Nbuffer > 0); + int printHead = VERBOSE && (table->hosts[i].stderr.Nbuffer > 0); if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v1)\n", table->hosts[i].hostname, table->hosts[i].stderr.Nbuffer); int Nout = write (STDOUT_FILENO, table->hosts[i].stderr.buffer, table->hosts[i].stderr.Nbuffer); @@ -314,5 +392,5 @@ if (!pid) { // fprintf (stderr, "no hosts to harvest\n"); - usleep (500000); + usleep (200000); break; // no outstanding jobs have finished } @@ -367,9 +445,6 @@ } } - if (!found) { - fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); - exit (2); - } - Nfound ++; + if (found) Nfound ++; + // we sometimes harvest children not in the list (eg, closed kapa window). just ignore if (Nfound == Nrunning) goto escape; // we've harvested all jobs } @@ -389,4 +464,180 @@ } +// wait for all children to complete, report output to stdout +int HostTableGroupWaitJobsGetIO (HostTableGroup *table, char *file, int lineno, int VERBOSE) { + + // we have launched table->Nhosts jobs; wait for all of them to complete... + // if one (N) failed to launch, we will get an ECHILD error from the last (N) calls + + // we need to read any data waiting on stderr or stdout from these jobs, or the overfull + // buffers can cause a problem. we alternate between 'select' and 'waitpid' calls with + // timeouts for both + + // add all hosts' sockets to the fd_sets + fd_set rdSet, wtSet; + FD_ZERO (&rdSet); + FD_ZERO (&wtSet); + + // XXX can I set the fd_sets once, since I am not actually closing the fd's? + + int globalStatus = TRUE; + + int i; + int Nmax = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i][0].pid) continue; // any unconnected hosts should be skipped + FD_SET (table->hosts[i][0].stdio[HOST_STDIN], &wtSet); + Nmax = MAX (Nmax, table->hosts[i][0].stdio[HOST_STDIN]); + FD_SET (table->hosts[i][0].stdio[HOST_STDOUT], &rdSet); + Nmax = MAX (Nmax, table->hosts[i][0].stdio[HOST_STDOUT]); + FD_SET (table->hosts[i][0].stdio[HOST_STDERR], &rdSet); + Nmax = MAX (Nmax, table->hosts[i][0].stdio[HOST_STDERR]); + } + Nmax ++; + + // need the list of connected hosts for exit test below + int Nrunning = 0; + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i][0].pid) continue; // any unconnected hosts should be skipped + Nrunning ++; + } + + int Nfound = 0; + + // this loop has 2 chunks: (a) check for I/O + (b) check for jobs done + while (1) { + + // Wait up to 0.5 second for host to provide I/O + // timeout gets mucked: need to reset before each select + struct timeval timeout; + timeout.tv_sec = 10; + timeout.tv_usec = 500000; + + int status = select (Nmax, NULL, &wtSet, NULL, &timeout); + if (status == -1) { + perror("select()"); + exit (2); + } + + // we have some sockets to check, check sockets for all hosts + for (i = 0; (status > 0) && (i < table->Nhosts); i++) { + if (!table->hosts[i][0].pid) continue; // any unconnected hosts should be skipped + + if (FALSE && FD_ISSET (table->hosts[i][0].stdio[HOST_STDIN], &wtSet)) { + // this host is waiting for input : this is an error, so exit + fprintf (stderr, "host %s is waiting for input\n", table->hosts[i][0].hostname); + abort(); + } + + if ((table->hosts[i][0].stdio[HOST_STDOUT] > 0) && FD_ISSET (table->hosts[i][0].stdio[HOST_STDOUT], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i][0].stdout, table->hosts[i][0].stdio[HOST_STDOUT]); + // if (table->hosts[i][0].stdout.Nbuffer > 0x10000) { + if (table->hosts[i][0].stdout.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i][0].stdout.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v1)\n", table->hosts[i][0].hostname, table->hosts[i][0].stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i][0].stdout.buffer, table->hosts[i][0].stdout.Nbuffer); + if (Nout != table->hosts[i][0].stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i][0].stdout); + if (printHead) fprintf (stdout, "\n"); + } + } + + if ((table->hosts[i][0].stdio[HOST_STDERR] > 0) && FD_ISSET (table->hosts[i][0].stdio[HOST_STDERR], &rdSet)) { + // this host has waiting output : read to buffer, and dump if necessary + ReadtoIOBuffer (&table->hosts[i][0].stderr, table->hosts[i][0].stdio[HOST_STDERR]); + // if (table->hosts[i][0].stderr.Nbuffer > 0x10000) { + if (table->hosts[i][0].stderr.Nbuffer > 0x1000) { + int printHead = VERBOSE || (table->hosts[i][0].stderr.Nbuffer > 0); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v1)\n", table->hosts[i][0].hostname, table->hosts[i][0].stderr.Nbuffer); + int Nout = write (STDOUT_FILENO, table->hosts[i][0].stderr.buffer, table->hosts[i][0].stderr.Nbuffer); + if (Nout != table->hosts[i][0].stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&table->hosts[i][0].stderr); + if (printHead) fprintf (stdout, "\n"); + } + } + } + + // now check if any children have finished... + while (TRUE) { + int status = 0; + int pid = waitpid (-1, &status, WNOHANG); + if (!pid) { + // fprintf (stderr, "no hosts to harvest\n"); + usleep (500000); + break; // no outstanding jobs have finished + } + if ((pid == -1) && (errno == ECHILD)) goto escape; // no more jobs on which to wait + if ((pid == -1) && (errno != ECHILD)) { + fprintf (stderr, "programming error (2)? %s %d", file, lineno); + exit (2); + } + + // find the host which has finished + int found = FALSE; + for (i = 0; (i < table->Nhosts) && !found; i++) { + if (table->hosts[i][0].pid != pid) continue; + found = TRUE; + + HostInfo *host = table->hosts[i]; + + // check on the status of this and report any output? + if (VERBOSE) fprintf (stdout, "job finished for %s (%d)\n", host->hostname, pid); + + // read stdout + int printHead; + printHead = VERBOSE || (host->stdout.Nbuffer > 0); + EmptyIOBuffer (&host->stdout, 100, host->stdio[HOST_STDOUT]); + if (printHead) fprintf (stdout, "--- stdout from %s --- (%d bytes, v2)\n", host->hostname, host->stdout.Nbuffer); + int Nout = write (STDOUT_FILENO, host->stdout.buffer, host->stdout.Nbuffer); + if (Nout != host->stdout.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stdout); + if (printHead) fprintf (stdout, "\n"); + + // read stderr + printHead = VERBOSE || (host->stderr.Nbuffer > 0); + EmptyIOBuffer (&host->stderr, 100, host->stdio[HOST_STDERR]); + if (printHead) fprintf (stdout, "--- stderr from %s --- (%d bytes, v2)\n", host->hostname, host->stderr.Nbuffer); + Nout = write (STDOUT_FILENO, host->stderr.buffer, host->stderr.Nbuffer); + if (Nout != host->stderr.Nbuffer) { fprintf (stderr, "(error writing log?)\n"); } + FlushIOBuffer (&host->stderr); + if (printHead) fprintf (stdout, "\n"); + + if (WIFEXITED(status)) { + if (VERBOSE) fprintf (stdout, "normal completion, exit status is %d\n", WEXITSTATUS(status)); + host->status = WEXITSTATUS(status); + if (host->status) { + fprintf (stdout, "job failed on %s\n", host->hostname); + globalStatus = FALSE; + } + } else { + host->status = -1; + fprintf (stdout, "job exited abnormally on %s\n", host->hostname); + globalStatus = FALSE; + continue; + } + } + if (!found) { + fprintf (stderr, "Programming error: failed to matched finished job to known host!\n"); + exit (2); + } + Nfound ++; + if (Nfound == Nrunning) goto escape; // we've harvested all jobs + } + } + +escape: + + // close all opened connections + for (i = 0; i < table->Nhosts; i++) { + if (!table->hosts[i][0].pid) continue; // any unconnected hosts should be skipped + close (table->hosts[i][0].stdio[HOST_STDIN]); + close (table->hosts[i][0].stdio[HOST_STDOUT]); + close (table->hosts[i][0].stdio[HOST_STDERR]); + } + + return globalStatus; +} + int HostTableTestHost (SkyRegion *region, int hostID) { Index: trunk/Ohana/src/libdvo/src/ImageMetadata.c =================================================================== --- trunk/Ohana/src/libdvo/src/ImageMetadata.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/ImageMetadata.c (revision 37807) @@ -156,13 +156,13 @@ // assign the storage arrays for (i = 0; i < Nimage; i++) { - int Nmosaic = FindMosaicForImage (image, Nimage, i); - if (!Nmosaic) continue; - Nmosaic --; + if (!image[i].coords.mosaic) continue; imageID[i] = image[i].imageID; externID[i] = image[i].externID; - crval1[i] = image[Nmosaic].coords.crval1; - crval2[i] = image[Nmosaic].coords.crval2; - - theta[i] = DEG_RAD*atan2(image[Nmosaic].coords.pc1_2, image[Nmosaic].coords.pc1_1); + Coords *mosaic = image[i].coords.mosaic; + + crval1[i] = mosaic->crval1; + crval2[i] = mosaic->crval2; + + theta[i] = DEG_RAD*atan2(mosaic->pc1_2, mosaic->pc1_1); Mcal[i] = image[i].Mcal; Index: trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c =================================================================== --- trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c (revision 37807) @@ -13,18 +13,10 @@ /* field defines a frame with 0,0 at the field center, and 1 arcsec / pixel */ - field.crpix1 = field.crpix2 = 0.0; + InitCoords (&field, "DEC--SIN"); field.cdelt1 = field.cdelt2 = 1.0 / 3600; - field.pc1_1 = field.pc2_2 = 1.0; - field.pc1_2 = field.pc2_1 = 0.0; - field.Npolyterms = 0; - strcpy (field.ctype, "DEC--SIN"); /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */ - mosaic.crpix1 = mosaic.crpix2 = 0.0; + InitCoords (&mosaic, "DEC--SIN"); mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600; - mosaic.pc1_1 = mosaic.pc2_2 = 1.0; - mosaic.pc1_2 = mosaic.pc2_1 = 0.0; - mosaic.Npolyterms = 0; - strcpy (mosaic.ctype, "DEC--SIN"); if ((image = ImageMetadataLoad (filename, &Nimage)) == NULL) return (FALSE); Index: trunk/Ohana/src/libdvo/src/ImageOps.c =================================================================== --- trunk/Ohana/src/libdvo/src/ImageOps.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/ImageOps.c (revision 37807) @@ -55,5 +55,4 @@ if (TimeSelect && ((image[i].tzero < tzero) || (image[i].tzero+image[i].trate*image[i].NY > tzero + trange))) continue; if (selection->useDisplay) { - if (!FindMosaicForImage (image, Nimage, i)) continue; // first check if region center is in image status = RD_to_XY (&X, &Y, Rmid, graph.coords.crval2, &image[i].coords); @@ -89,5 +88,4 @@ } if (selection->useSkyregion) { - if (!FindMosaicForImage (image, Nimage, i)) continue; /* project this image to screen display coords */ x[0] = 0; y[0] = 0; Index: trunk/Ohana/src/libdvo/src/ImageSelection.c =================================================================== --- trunk/Ohana/src/libdvo/src/ImageSelection.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/ImageSelection.c (revision 37807) @@ -22,11 +22,7 @@ if (mosaicMode) { - /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */ - mosaic.crpix1 = mosaic.crpix2 = 0.0; - mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600; - mosaic.pc1_1 = mosaic.pc2_2 = 1.0; - mosaic.pc1_2 = mosaic.pc2_1 = 0.0; - mosaic.Npolyterms = 0; - strcpy (mosaic.ctype, "DEC--SIN"); + /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */ + InitCoords (&mosaic, "DEC--SIN"); + mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600; } @@ -64,5 +60,4 @@ } if (m == -1) return (NULL); - if (!FindMosaicForImage (image, Nimage, m)) return (NULL); return (&image[m]); } @@ -85,5 +80,4 @@ } if (m == -1) return (NULL); - if (!FindMosaicForImage (image, Nimage, m)) return (NULL); return (&image[m]); } Index: trunk/Ohana/src/libdvo/src/LoadImages.c =================================================================== --- trunk/Ohana/src/libdvo/src/LoadImages.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/LoadImages.c (revision 37807) @@ -9,4 +9,5 @@ static Image *imageCache = NULL; static off_t cacheNimage = 0; +AstromOffsetTable *table = NULL; static time_t getLastModified(char *filename); @@ -39,5 +40,4 @@ lastModified = 0; } - gfits_db_init (&db); @@ -77,4 +77,19 @@ } + // assign image->parent and image->coords.mosaic + BuildChipMatch (image, *Nimage); + // MARKTIME("build chip match: %f sec\n", dtime); + + char mapfile[DVO_MAX_PATH]; + snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", catdir); + + if (table) AstromOffsetTableFree(table); + table = AstromOffsetMapLoad (mapfile, FALSE); + + // assign images.coords.offsetMap -> table->map[i] + if (table) { + AstromOffsetTableMatchChips (image, *Nimage, table); + } + return (image); } @@ -92,4 +107,5 @@ if (!dvoUseImageCache && (images != NULL)) { free(images); + AstromOffsetTableFree (table); } else { // defer free until next LoadImages with a different or modified Images.dat Index: trunk/Ohana/src/libdvo/src/RegionHostTable.c =================================================================== --- trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/RegionHostTable.c (revision 37807) @@ -31,4 +31,6 @@ hosts[i].Nneighbors = 0; hosts[i].isNeighbor = FALSE; + + hosts[i].astromTable = NULL; } return; Index: trunk/Ohana/src/libdvo/src/astrom_maps.c =================================================================== --- trunk/Ohana/src/libdvo/src/astrom_maps.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/astrom_maps.c (revision 37807) @@ -0,0 +1,70 @@ +# include + +int MatchChipsToAstromMaps (Image *images, off_t Nimages, AstromOffsetTable *table) { + + // we have a table of astrometry offset maps. we want to find the chips that match to + // each of the maps so we can assign the map to the image.coords.offsetMap entry + + // we have the lookup table of table->IDtoSeq to find the match by imageID + + // we are just going to assume that table->IDtoSeq is correct and complete + for (i = 0; i < Nimages; i++) { + int imageID = images[i].imageID; + if (imageID < 0) continue; + if (imageID > table->maxID) continue; + + int seq = table->IDtoSeq[imageID]; + if (seq < 0) continue; + if (seq >= table->Nmap) continue; // this one is probably not valid, right? + + images[i].coords.offsetMap = &table->map[seq]; + } + + return (TRUE); +} + +int AstromOffsetTableNewMap (AstromOffsetTable *table, int order, Image *image) { + + int Nmap = table->Nmap; + + table->Nmap++; + REALLOCATE (table->map, AstromOffsetMap, table->Nmap); + + int Nx = order; + int Ny = order; + + if (image->imageID > table->MaxID) { + int oldMaxID = table->MaxID; + table->MaxID = image->imageID; + REALLOCATE (table->IDtoSeq, int, table->MaxID + 1); + for (i = oldMaxID + 1; i < table->MaxID + 1; i++) { + table->IDtoSeq[i] = -1; + } + } + myAssert (table->IDtoSeq[image->imageID] == -1, "table IDtoSeq not initiazed or image collision"); + table->IDtoSeq[image->imageID] = Nmap; + + table->map[Nmap].Nx = Nx; + table->map[Nmap].Ny = Ny; + table->map[Nmap].ID = table->maxID; table->maxID ++; + table->map[Nmap].imageID = image->imageID; + + table->map[Nmap].dX = Nx / image->Nx; + table->map[Nmap].dY = Ny / image->Ny; + + ALLOCATE (table->map[Nmap].dXv, float *, Nx); + ALLOCATE (table->map[Nmap].dYv, float *, Nx); + + for (j = 0; Nx; j++) { + ALLOCATE (table->map[i].dXv[j], float, Ny); + ALLOCATE (table->map[i].dYv[j], float, Ny); + + for (k = 0; k < Ny; k++) { + table->map[i].dXv[j][k] = 0.0; + table->map[i].dYv[j][k] = 0.0; + } + } + image[0].coords.offsetMap = &table->map[Nmap]; + return TRUE; +} + Index: trunk/Ohana/src/libdvo/src/coord_systems.c =================================================================== --- trunk/Ohana/src/libdvo/src/coord_systems.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/coord_systems.c (revision 37807) @@ -1,4 +1,3 @@ # include "dvo.h" -# include "libdvo_astro.h" CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output) { @@ -40,6 +39,6 @@ // \delta^p_J = +27d 07' 41.7087" // \theta_J = 122.93192526 deg + // these differ from the above at the several milliarcsec level - // These values were used in the past, but transform to/from B1950.0 @@ -47,4 +46,20 @@ // transform->Xo = 282.25*RAD_DEG; // transform->xo = 33.00; + + break; + case COORD_GALACTIC_REID_2004: + + // J2000 FK5 transformation from Reid et al 2004, ApJ 616, 872 + transform->phi = (27.128336111 - 90.0)*RAD_DEG; // +27d 07' 42.01" (J2000) + transform->Xo = (192.859508333 + 90.0)*RAD_DEG; // 12h 51m 26.282s (J2000) + transform->xo = 32.932; + + // comments from LSD code (lsd/builtins/misc.py, Mario Juric) + // # Appendix of Reid et al. (http://adsabs.harvard.edu/cgi-bin/bib_query?2004ApJ...616..872R) + // # This convention is also used by LAMBDA/WMAP (http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm) + + // _angp = np.radians(192.859508333) # 12h 51m 26.282s (J2000) + // _dngp = np.radians(27.128336111) # +27d 07' 42.01" (J2000) + // _l0 = np.radians(32.932) break; @@ -87,8 +102,33 @@ switch (output) { case COORD_CELESTIAL: - // J2000 transformation from Liu et al 2011 () + // J2000 transformation from Liu et al 2011 (A&A 526, A16) transform->phi = 62.8717488056*RAD_DEG; transform->Xo = 32.9319185700*RAD_DEG; transform->xo = 282.8594812080; + + // These values were used in the past, but transform to/from B1950.0 + // transform->phi = 62.60*RAD_DEG; + // transform->Xo = 33.00*RAD_DEG; + // transform->xo = 282.25; + break; + case COORD_GALACTIC: + transform->phi = 0.0; + transform->Xo = 0.0; + transform->xo = 0.0; + transform->isIdentity = TRUE; + break; + case COORD_ECLIPTIC: + return NULL; + default: + abort(); + } + break; + case COORD_GALACTIC_REID_2004: + switch (output) { + case COORD_CELESTIAL: + // J2000 FK5 transformation from Reid et al 2004, ApJ 616, 872 + transform->phi = (90.0 - 27.128336111)*RAD_DEG; // +27d 07' 42.01" (J2000) + transform->Xo = 32.932*RAD_DEG; + transform->xo = (192.859508333 + 90.0); // 12h 51m 26.282s (J2000) // These values were used in the past, but transform to/from B1950.0 @@ -144,10 +184,16 @@ Y *= RAD_DEG; + double cosY = cos(Y); + double sinY = sin(Y); + + double cosX = cos(X); + double sinX = sin(X); + // recast with constants extracted: - sin_y = cos(Y)*sin(X)*transform->sin_phi_cos_Xo - cos(Y)*cos(X)*transform->sin_phi_sin_Xo + sin(Y)*transform->cos_phi; + sin_y = cosY*sinX*transform->sin_phi_cos_Xo - cosY*cosX*transform->sin_phi_sin_Xo + sinY*transform->cos_phi; cos_y = sqrt (1 - sin_y*sin_y); - sin_x = cos(Y)*sin(X)*transform->cos_phi_cos_Xo - cos(Y)*cos(X)*transform->cos_phi_sin_Xo - sin(Y)*transform->sin_phi; - cos_x = cos(Y)*cos(X)*transform->cos_Xo + cos(Y)*sin(X)*transform->sin_Xo; + sin_x = cosY*sinX*transform->cos_phi_cos_Xo - cosY*cosX*transform->cos_phi_sin_Xo - sinY*transform->sin_phi; + cos_x = cosY*cosX*transform->cos_Xo + cosY*sinX*transform->sin_Xo; // atan2 returns -pi : +pi Index: trunk/Ohana/src/libdvo/src/coordops.c =================================================================== --- trunk/Ohana/src/libdvo/src/coordops.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/coordops.c (revision 37807) @@ -3,20 +3,51 @@ /* note that Coords.ctype carries the DEC (ctype2) value */ -static Coords mosaic; -static int gotMosaic = FALSE; - -void RegisterMosaic (Coords *coords) { - mosaic = *coords; - gotMosaic = TRUE; -} - -int isRegisteredMosaic () { - return (gotMosaic); -} - -int GetMosaicCoords (Coords *coords) { - *coords = mosaic; - return (TRUE); -} +// default coordinates have unit scale and no rotation +void InitCoords (Coords *coords, char *projection) { + coords->crval1 = 0.0; + coords->crval2 = 0.0; + coords->crpix1 = 0.0; + coords->crpix2 = 0.0; + + coords->cdelt1 = 1.0; + coords->cdelt2 = 1.0; + + coords->pc1_1 = 1.0; // EW is flipped relative to sky parity (dRA/dX > 0) + coords->pc2_2 = 1.0; + coords->pc2_1 = 0.0; + coords->pc1_2 = 0.0; + + coords->Npolyterms = 1; + coords->mosaic = NULL; + coords->offsetMap = NULL; + + memset (coords->polyterms, 0, 7*2*sizeof(float)); + + if (projection) { + strcpy (coords->ctype, projection); + } else { + strcpy (coords->ctype, "NONE"); + } +} + +// default coordinates have unit scale and no rotation +void CopyCoords (Coords *tgt, Coords *src) { + tgt->crval1 = src->crval1; + tgt->crval2 = src->crval2; + tgt->crpix1 = src->crpix1; + tgt->crpix2 = src->crpix2; + tgt->cdelt1 = src->cdelt1; + tgt->cdelt2 = src->cdelt2; + tgt->pc1_1 = src->pc1_1; + tgt->pc2_2 = src->pc2_2; + tgt->pc2_1 = src->pc2_1; + tgt->pc1_2 = src->pc1_2; + tgt->Npolyterms = src->Npolyterms; + tgt->mosaic = src->mosaic; // pointer to another structure + tgt->offsetMap = src->offsetMap; // pointer to another structure + + memcpy (tgt->polyterms, src->polyterms, 7*2*sizeof(float)); + strcpy (tgt->ctype, src->ctype); +} int XY_to_LM (double *L, double *M, double x, double y, Coords *coords) { @@ -28,4 +59,17 @@ proj = GetProjection (coords[0].ctype); mode = GetProjectionMode (proj); + + // if we have an image map, it should be applied to the pixel coordinates + // before anything else is done. It is easier to invert the map if it is + // of the form (dX,dY) = f(X,Y), then if (dL,dM) = f(X,Y). + if (coords[0].Npolyterms == -1) { + myAssert (coords[0].offsetMap, "offsetMap is requested by not defined"); + float dX = AstromOffsetMapValue (coords->offsetMap, x, y, TRUE); + float dY = AstromOffsetMapValue (coords->offsetMap, x, y, FALSE); + if (isnan(dX)) dX = 0.0; + if (isnan(dY)) dY = 0.0; + x += dX; + y += dY; + } /** convert pixel coordinates to cartesian system **/ @@ -78,6 +122,9 @@ /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */ if (proj == PROJ_WRP) { - if (!gotMosaic) return (FALSE); - XY_to_RD (ra, dec, L, M, &mosaic); + if (!coords->mosaic) { + myAbort ("missing mosaic element"); + return (FALSE); + } + XY_to_RD (ra, dec, L, M, coords->mosaic); } return (TRUE); @@ -149,5 +196,5 @@ *dec = DEG_RAD*delta; - *ra = ohana_normalize_angle (*ra); + // *ra = ohana_normalize_angle (*ra); return (TRUE); @@ -182,5 +229,5 @@ *dec = delta + coords[0].crval2; - *ra = ohana_normalize_angle (*ra); + // *ra = ohana_normalize_angle (*ra); return (TRUE); @@ -208,6 +255,9 @@ if (mode == PROJ_MODE_CARTESIAN) { if (proj == PROJ_WRP) { - if (!gotMosaic) return (FALSE); - RD_to_XY (L, M, ra, dec, &mosaic); + if (!coords->mosaic) { + myAbort ("missing mosaic element"); + return (FALSE); + } + RD_to_XY (L, M, ra, dec, coords->mosaic); return (TRUE); } @@ -381,4 +431,43 @@ *y = Yo / coords[0].cdelt2 + coords[0].crpix2; + // if we have an image map, it should be applied to the pixel coordinates + // after anything else is done. It is easier to invert the map if it is + // of the form (dX,dY) = f(X,Y), then if (dL,dM) = f(X,Y). + if (coords[0].Npolyterms == -1) { + myAssert (coords[0].offsetMap, "offsetMap is requested by not defined"); + + double xraw = *x; + double yraw = *y; + + double dXo = 0.0; + double dYo = 0.0; + + int i; + for (i = 0; i < 4; i++) { + double dX = AstromOffsetMapValue (coords->offsetMap, xraw, yraw, TRUE); + double dY = AstromOffsetMapValue (coords->offsetMap, xraw, yraw, FALSE); + + if (isnan(dX)) dX = 0.0; + if (isnan(dY)) dY = 0.0; + + dX -= dXo; + dY -= dYo; + + xraw -= dX; + yraw -= dY; + + dXo += dX; + dYo += dY; + } + + *x = xraw; + *y = yraw; + + // I need to iterate since the position at which I have first made the correction is + // not the true position. but if dX,dY is small and the gradient of dX,dY is also + // small, then even a single pass gets us close. + + // XXX test this and use a while (hypot(dX,dY) > XXX) condition + } return (TRUE); } @@ -449,10 +538,6 @@ rotate = 0.0; - coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; - coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; - coords[0].pc1_1 = coords[0].pc2_2 = 1.0; - coords[0].pc2_1 = coords[0].pc1_2 = 0.0; - coords[0].Npolyterms = 1; - strcpy (coords[0].ctype, "NONE"); + + InitCoords (coords, NULL); mode = COORD_TYPE_NONE; @@ -619,5 +704,5 @@ } } - + if (!status) { // fprintf (stderr, "error getting all elements for coordinate mode %s\n", coords[0].ctype); Index: trunk/Ohana/src/libdvo/src/dbExtractAverages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dbExtractAverages.c (revision 37807) @@ -50,5 +50,5 @@ /* return average.field based on the selection */ -dbValue dbExtractAverages (Average *average, SecFilt *secfilt, Measure *measure, dbField *field) { +dbValue dbExtractAverages (Average *average, SecFilt *secfilt, Measure *measure, Lensobj *lensobj, StarPar *starpar, dbField *field) { // off_t i; @@ -161,6 +161,9 @@ value.Flt = average[0].stargal; break; - case AVE_REF_COLOR: + case AVE_REF_COLOR_BLUE: value.Flt = average[0].refColorBlue; + break; + case AVE_REF_COLOR_RED: + value.Flt = average[0].refColorRed; break; @@ -306,7 +309,4 @@ value.Flt = PhotUCdist (field->photcode, average, secfilt); break; - // XX case MAG_OPTION_STACK_DET_ID: - // XX value.Int = PhotStackID (field->photcode, average, secfilt); - // XX break; case MAG_OPTION_FLAGS: { int Nsec = GetPhotcodeNsec (field->photcode->code); @@ -315,8 +315,50 @@ break; } + + // g:X11_SM_OBJ, etc, are only valid for average + case MAG_OPTION_X11_SM_OBJ: { value.Flt = LensValue_X11_sm_obj (field->photcode, lensobj); break; } + case MAG_OPTION_X12_SM_OBJ: { value.Flt = LensValue_X12_sm_obj (field->photcode, lensobj); break; } + case MAG_OPTION_X22_SM_OBJ: { value.Flt = LensValue_X22_sm_obj (field->photcode, lensobj); break; } + case MAG_OPTION_E1_SM_OBJ: { value.Flt = LensValue_E1_sm_obj (field->photcode, lensobj); break; } + case MAG_OPTION_E2_SM_OBJ: { value.Flt = LensValue_E2_sm_obj (field->photcode, lensobj); break; } + + case MAG_OPTION_X11_SH_OBJ: { value.Flt = LensValue_X11_sh_obj (field->photcode, lensobj); break; } + case MAG_OPTION_X12_SH_OBJ: { value.Flt = LensValue_X12_sh_obj (field->photcode, lensobj); break; } + case MAG_OPTION_X22_SH_OBJ: { value.Flt = LensValue_X22_sh_obj (field->photcode, lensobj); break; } + case MAG_OPTION_E1_SH_OBJ: { value.Flt = LensValue_E1_sh_obj (field->photcode, lensobj); break; } + case MAG_OPTION_E2_SH_OBJ: { value.Flt = LensValue_E2_sh_obj (field->photcode, lensobj); break; } + + case MAG_OPTION_X11_SM_PSF: { value.Flt = LensValue_X11_sm_psf (field->photcode, lensobj); break; } + case MAG_OPTION_X12_SM_PSF: { value.Flt = LensValue_X12_sm_psf (field->photcode, lensobj); break; } + case MAG_OPTION_X22_SM_PSF: { value.Flt = LensValue_X22_sm_psf (field->photcode, lensobj); break; } + case MAG_OPTION_E1_SM_PSF: { value.Flt = LensValue_E1_sm_psf (field->photcode, lensobj); break; } + case MAG_OPTION_E2_SM_PSF: { value.Flt = LensValue_E2_sm_psf (field->photcode, lensobj); break; } + + case MAG_OPTION_X11_SH_PSF: { value.Flt = LensValue_X11_sh_psf (field->photcode, lensobj); break; } + case MAG_OPTION_X12_SH_PSF: { value.Flt = LensValue_X12_sh_psf (field->photcode, lensobj); break; } + case MAG_OPTION_X22_SH_PSF: { value.Flt = LensValue_X22_sh_psf (field->photcode, lensobj); break; } + case MAG_OPTION_E1_SH_PSF: { value.Flt = LensValue_E1_sh_psf (field->photcode, lensobj); break; } + case MAG_OPTION_E2_SH_PSF: { value.Flt = LensValue_E2_sh_psf (field->photcode, lensobj); break; } + case MAG_OPTION_NONE: break; - } - break; + + default: + break; + } + break; + + case AVE_E_BV: { value.Flt = starpar->Ebv; break; } + case AVE_E_BV_ERR: { value.Flt = starpar->dEbv; break; } + case AVE_DISTANCE_MOD: { value.Flt = starpar->DistMag; break; } + case AVE_DISTANCE_MOD_ERR: { value.Flt = starpar->dDistMag; break; } + case AVE_M_R: { value.Flt = starpar->M_r; break; } + case AVE_M_R_ERR: { value.Flt = starpar->dM_r; break; } + case AVE_FEH: { value.Flt = starpar->FeH; break; } + case AVE_FEH_ERR: { value.Flt = starpar->dFeH; break; } + case AVE_URA_GALMODEL: { value.Flt = starpar->uRA; break; } + case AVE_UDEC_GALMODEL: { value.Flt = starpar->uDEC; break; } + case AVE_RA_GALMODEL: { value.Flt = starpar->R; break; } + case AVE_DEC_GALMODEL: { value.Flt = starpar->D; break; } case AVE_TYPE: @@ -324,4 +366,6 @@ case AVE_TYPEFRAC: break; + + } return (value); Index: trunk/Ohana/src/libdvo/src/dbExtractImages.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dbExtractImages.c (revision 37807) @@ -59,5 +59,4 @@ time_t t; dbValue value; - off_t Nmosaic; value.Flt = NAN; @@ -68,5 +67,4 @@ case IMAGE_RA: if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -83,5 +81,4 @@ case IMAGE_DEC: if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -99,5 +96,4 @@ if (!haveGalactic) { if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -118,5 +114,4 @@ if (!haveGalactic) { if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -137,5 +132,4 @@ if (!haveEcliptic) { if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -156,5 +150,4 @@ if (!haveEcliptic) { if (!haveCelestial) { - if (!FindMosaicForImage (image, Nimage, N)) return value; if (!strcmp(&image[N].coords.ctype[4], "-DIS")) { x = 0.0; @@ -269,17 +262,14 @@ case IMAGE_FWHM_MEDIAN: - if (!(Nmosaic = FindMosaicForImage (image, Nimage, N))) return value; - Nmosaic --; // XXX kind of a hack: FindMosaicForImage returns 0 or the mosaic seq number + 1 - value.Flt = (image[Nmosaic].fwhm_x + image[Nmosaic].fwhm_y) / 50.0; + if (!image[N].parent) return value; + value.Flt = (image[N].parent->fwhm_x + image[N].parent->fwhm_y) / 50.0; break; case IMAGE_FWHM_MAJ_MEDIAN: - if (!(Nmosaic = FindMosaicForImage (image, Nimage, N))) return value; - Nmosaic --; // XXX kind of a hack: FindMosaicForImage returns 0 or the mosaic seq number + 1 - value.Flt = image[Nmosaic].fwhm_x / 25.0; + if (!image[N].parent) return value; + value.Flt = image[N].parent->fwhm_x / 25.0; break; case IMAGE_FWHM_MIN_MEDIAN: - if (!(Nmosaic = FindMosaicForImage (image, Nimage, N))) return value; - Nmosaic --; // XXX kind of a hack: FindMosaicForImage returns 0 or the mosaic seq number + 1 - value.Flt = image[Nmosaic].fwhm_y / 25.0; + if (!image[N].parent) return value; + value.Flt = image[N].parent->fwhm_y / 25.0; break; @@ -375,6 +365,9 @@ break; - case IMAGE_REF_COLOR: + case IMAGE_REF_COLOR_BLUE: value.Flt = image[N].refColorBlue; + break; + case IMAGE_REF_COLOR_RED: + value.Flt = image[N].refColorRed; break; } Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c (revision 37807) @@ -85,5 +85,5 @@ /* return measure.field based on the selection */ -dbValue dbExtractMeasures (Average *average, SecFilt *secfilt, Measure *measure, dbField *field) { +dbValue dbExtractMeasures (Average *average, SecFilt *secfilt, Measure *measure, Lensing *lensing, StarPar *starpar, dbField *field) { int Nsec; @@ -249,4 +249,33 @@ case MAG_OPTION_UC_DIST: case MAG_OPTION_FLAGS: + case MAG_OPTION_X11_SM_OBJ: + case MAG_OPTION_X12_SM_OBJ: + case MAG_OPTION_X22_SM_OBJ: + case MAG_OPTION_E1_SM_OBJ: + case MAG_OPTION_E2_SM_OBJ: + case MAG_OPTION_X11_SH_OBJ: + case MAG_OPTION_X12_SH_OBJ: + case MAG_OPTION_X22_SH_OBJ: + case MAG_OPTION_E1_SH_OBJ: + case MAG_OPTION_E2_SH_OBJ: + case MAG_OPTION_X11_SM_PSF: + case MAG_OPTION_X12_SM_PSF: + case MAG_OPTION_X22_SM_PSF: + case MAG_OPTION_E1_SM_PSF: + case MAG_OPTION_E2_SM_PSF: + case MAG_OPTION_X11_SH_PSF: + case MAG_OPTION_X12_SH_PSF: + case MAG_OPTION_X22_SH_PSF: + case MAG_OPTION_E1_SH_PSF: + case MAG_OPTION_E2_SH_PSF: + case MAG_OPTION_F_AP_R5: + case MAG_OPTION_F_ERR_AP_R5: + case MAG_OPTION_F_STDEV_AP_R5: + case MAG_OPTION_F_FILL_AP_R5: + case MAG_OPTION_F_AP_R6: + case MAG_OPTION_F_ERR_AP_R6: + case MAG_OPTION_F_STDEV_AP_R6: + case MAG_OPTION_F_FILL_AP_R6: + case MAG_OPTION_NONE: break; @@ -474,6 +503,9 @@ break; - case MEAS_REF_COLOR: + case MEAS_REF_COLOR_BLUE: value.Flt = average[0].refColorBlue; + break; + case MEAS_REF_COLOR_RED: + value.Flt = average[0].refColorRed; break; @@ -528,19 +560,16 @@ # endif break; - case MEAS_XFIX: /* OK */ - value.Flt = measure[0].Xfix; - break; - case MEAS_YFIX: /* OK */ - value.Flt = measure[0].Yfix; - break; - case MEAS_XCCD_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dXccd); - break; - case MEAS_YCCD_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dYccd); - break; - case MEAS_POS_SYS_ERR: /* OK */ - value.Flt = FromShortPixels(measure[0].dRsys); - break; + case MEAS_XFIX: { value.Flt = measure[0].Xfix; break; } + case MEAS_YFIX: { value.Flt = measure[0].Yfix; break; } + case MEAS_XCCD_ERR: { value.Flt = FromShortPixels(measure[0].dXccd); break; } + case MEAS_YCCD_ERR: { value.Flt = FromShortPixels(measure[0].dYccd); break; } + case MEAS_POS_SYS_ERR: { value.Flt = FromShortPixels(measure[0].dRsys); } + + case MEAS_XOFF_KH: { value.Flt = measure[0].XoffKH; break; } + case MEAS_YOFF_KH: { value.Flt = measure[0].YoffKH; break; } + case MEAS_XOFF_DCR: { value.Flt = measure[0].XoffDCR; break; } + case MEAS_YOFF_DCR: { value.Flt = measure[0].YoffDCR; break; } + case MEAS_ROFF_GAL: { value.Flt = measure[0].RoffGAL; break; } + case MEAS_DOFF_GAL: { value.Flt = measure[0].DoffGAL; break; } case MEAS_XFIELD: /* offset relative to exposure center in ra,dec space */ @@ -581,12 +610,5 @@ case MEAS_XMOSAIC: /* offset relative to exposure center in camera coords */ if (!haveMosaicMeas) { - if (REMOTE_CLIENT) { - mosaic = MatchMosaicMetadata (measure[0].imageID); - } else { - // fprintf (stderr, "non-parallel Xmos broken\n"); - // abort(); - // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); - mosaic = MatchMosaicMetadata (measure[0].imageID); - } + mosaic = MatchMosaicMetadata (measure[0].imageID); if (mosaic == NULL) break; double Rm = measure[0].R; @@ -598,12 +620,5 @@ case MEAS_YMOSAIC: /* OK */ if (!haveMosaicMeas) { - if (REMOTE_CLIENT) { - mosaic = MatchMosaicMetadata (measure[0].imageID); - } else { - // fprintf (stderr, "non-parallel Xmos broken\n"); - // abort(); - // mosaic = MatchMosaic (measure[0].t, measure[0].photcode); - mosaic = MatchMosaicMetadata (measure[0].imageID); - } + mosaic = MatchMosaicMetadata (measure[0].imageID); if (mosaic == NULL) break; double Rm = measure[0].R; @@ -680,7 +695,5 @@ break; - case MEAS_MCAL_OFFSET: /* OK */ - value.Flt = measure[0].Mcal; - break; + case MEAS_MCAL_OFFSET: { value.Flt = measure[0].Mcal; break; } case MEAS_FLAT: /* OK */ @@ -717,4 +730,17 @@ } break; + + case MEAS_E_BV: { value.Flt = starpar->Ebv; break; } + case MEAS_E_BV_ERR: { value.Flt = starpar->dEbv; break; } + case MEAS_DISTANCE_MOD: { value.Flt = starpar->DistMag; break; } + case MEAS_DISTANCE_MOD_ERR: { value.Flt = starpar->dDistMag; break; } + case MEAS_M_R: { value.Flt = starpar->M_r; break; } + case MEAS_M_R_ERR: { value.Flt = starpar->dM_r; break; } + case MEAS_FEH: { value.Flt = starpar->FeH; break; } + case MEAS_FEH_ERR: { value.Flt = starpar->dFeH; break; } + case MEAS_URA_GALMODEL: { value.Flt = starpar->uRA; break; } + case MEAS_UDEC_GALMODEL: { value.Flt = starpar->uDEC; break; } + case MEAS_RA_GALMODEL: { value.Flt = starpar->R; break; } + case MEAS_DEC_GALMODEL: { value.Flt = starpar->D; break; } // case MEAS_FLUX_PSF: /* OK */ Index: trunk/Ohana/src/libdvo/src/dbFields.c =================================================================== --- trunk/Ohana/src/libdvo/src/dbFields.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dbFields.c (revision 37807) @@ -28,8 +28,108 @@ } +int dbFieldNeedMeasure (dbField *fields, int Nfields) { + int i; + for (i = 0; i < Nfields; i++) { + if (fields[i].photcode == NULL) continue; // non-measure fields do not have a photcode + if (fields[i].photcode[0].type == PHOT_REF) return TRUE; + if (fields[i].photcode[0].type == PHOT_DEP) return TRUE; + } + return FALSE; +} + +int dbFieldNeedLensing (dbField *fields, int Nfields) { + int i; + for (i = 0; i < Nfields; i++) { + if (fields[i].ID == MEAS_X11_SM_OBJ) return TRUE; + if (fields[i].ID == MEAS_X12_SM_OBJ) return TRUE; + if (fields[i].ID == MEAS_X22_SM_OBJ) return TRUE; + if (fields[i].ID == MEAS_E1_SM_OBJ) return TRUE; + if (fields[i].ID == MEAS_E2_SM_OBJ) return TRUE; + if (fields[i].ID == MEAS_X11_SH_OBJ) return TRUE; + if (fields[i].ID == MEAS_X12_SH_OBJ) return TRUE; + if (fields[i].ID == MEAS_X22_SH_OBJ) return TRUE; + if (fields[i].ID == MEAS_E1_SH_OBJ) return TRUE; + if (fields[i].ID == MEAS_E2_SH_OBJ) return TRUE; + if (fields[i].ID == MEAS_X11_SM_PSF) return TRUE; + if (fields[i].ID == MEAS_X12_SM_PSF) return TRUE; + if (fields[i].ID == MEAS_X22_SM_PSF) return TRUE; + if (fields[i].ID == MEAS_E1_SM_PSF) return TRUE; + if (fields[i].ID == MEAS_E2_SM_PSF) return TRUE; + if (fields[i].ID == MEAS_X11_SH_PSF) return TRUE; + if (fields[i].ID == MEAS_X12_SH_PSF) return TRUE; + if (fields[i].ID == MEAS_X22_SH_PSF) return TRUE; + if (fields[i].ID == MEAS_E1_SH_PSF) return TRUE; + if (fields[i].ID == MEAS_E2_SH_PSF) return TRUE; + } + return FALSE; +} + +int dbFieldNeedLensobj (dbField *fields, int Nfields) { + int i; + for (i = 0; i < Nfields; i++) { + if (fields[i].magOption == MAG_OPTION_NONE) continue; // non-measure fields do not have a photcode + + if (fields[i].magOption == MAG_OPTION_X11_SM_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X12_SM_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X22_SM_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_E1_SM_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_E2_SM_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X11_SH_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X12_SH_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X22_SH_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_E1_SH_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_E2_SH_OBJ) return TRUE; + if (fields[i].magOption == MAG_OPTION_X11_SM_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_X12_SM_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_X22_SM_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_E1_SM_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_E2_SM_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_X11_SH_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_X12_SH_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_X22_SH_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_E1_SH_PSF) return TRUE; + if (fields[i].magOption == MAG_OPTION_E2_SH_PSF) return TRUE; + } + return FALSE; +} + +int dbFieldNeedStarpar (dbField *fields, int Nfields, int isAverage) { + int i; + for (i = 0; i < Nfields; i++) { + if (isAverage) { + if (fields[i].ID == AVE_E_BV) return TRUE; + if (fields[i].ID == AVE_E_BV_ERR) return TRUE; + if (fields[i].ID == AVE_DISTANCE_MOD) return TRUE; + if (fields[i].ID == AVE_DISTANCE_MOD_ERR) return TRUE; + if (fields[i].ID == AVE_M_R) return TRUE; + if (fields[i].ID == AVE_M_R_ERR) return TRUE; + if (fields[i].ID == AVE_FEH) return TRUE; + if (fields[i].ID == AVE_FEH_ERR) return TRUE; + if (fields[i].ID == AVE_URA_GALMODEL) return TRUE; + if (fields[i].ID == AVE_UDEC_GALMODEL) return TRUE; + if (fields[i].ID == AVE_RA_GALMODEL) return TRUE; + if (fields[i].ID == AVE_DEC_GALMODEL) return TRUE; + } else { + if (fields[i].ID == MEAS_E_BV) return TRUE; + if (fields[i].ID == MEAS_E_BV_ERR) return TRUE; + if (fields[i].ID == MEAS_DISTANCE_MOD) return TRUE; + if (fields[i].ID == MEAS_DISTANCE_MOD_ERR) return TRUE; + if (fields[i].ID == MEAS_M_R) return TRUE; + if (fields[i].ID == MEAS_M_R_ERR) return TRUE; + if (fields[i].ID == MEAS_FEH) return TRUE; + if (fields[i].ID == MEAS_FEH_ERR) return TRUE; + if (fields[i].ID == MEAS_URA_GALMODEL) return TRUE; + if (fields[i].ID == MEAS_UDEC_GALMODEL) return TRUE; + if (fields[i].ID == MEAS_RA_GALMODEL) return TRUE; + if (fields[i].ID == MEAS_DEC_GALMODEL) return TRUE; + } + } + return FALSE; +} + dvoMagSourceType GetMagSource (char *string) { - if (!strcasecmp (string, "chip")) return (MAG_SRC_CHP); - if (!strcasecmp (string, "warp")) return (MAG_SRC_WRP); + if (!strcasecmp (string, "chip")) return (MAG_SRC_CHP); + if (!strcasecmp (string, "warp")) return (MAG_SRC_WRP); if (!strcasecmp (string, "stack")) return (MAG_SRC_STK); @@ -66,89 +166,48 @@ dvoMagOptionType GetMagOption (char *string) { - if (!strcasecmp (string, "mag")) return (MAG_OPTION_MAG); - if (!strcasecmp (string, "err")) return (MAG_OPTION_ERR); - if (!strcasecmp (string, "magerr")) return (MAG_OPTION_ERR); - if (!strcasecmp (string, "flux")) return (MAG_OPTION_FLUX); - if (!strcasecmp (string, "fluxerr")) return (MAG_OPTION_FLUX_ERR); - if (!strcasecmp (string, "stdev")) return (MAG_OPTION_STDEV); - if (!strcasecmp (string, "chisq")) return (MAG_OPTION_CHISQ); - if (!strcasecmp (string, "min")) return (MAG_OPTION_MIN); - if (!strcasecmp (string, "max")) return (MAG_OPTION_MAX); - if (!strcasecmp (string, "ncode")) return (MAG_OPTION_NCODE); - if (!strcasecmp (string, "nphot")) return (MAG_OPTION_NPHOT); - if (!strcasecmp (string, "uc_dist")) return (MAG_OPTION_UC_DIST); - if (!strcasecmp (string, "flags")) return (MAG_OPTION_FLAGS); + if (!strcasecmp (string, "mag")) return MAG_OPTION_MAG; + if (!strcasecmp (string, "err")) return MAG_OPTION_ERR; + if (!strcasecmp (string, "magerr")) return MAG_OPTION_ERR; + if (!strcasecmp (string, "flux")) return MAG_OPTION_FLUX; + if (!strcasecmp (string, "fluxerr")) return MAG_OPTION_FLUX_ERR; + if (!strcasecmp (string, "stdev")) return MAG_OPTION_STDEV; + if (!strcasecmp (string, "chisq")) return MAG_OPTION_CHISQ; + if (!strcasecmp (string, "min")) return MAG_OPTION_MIN; + if (!strcasecmp (string, "max")) return MAG_OPTION_MAX; + if (!strcasecmp (string, "ncode")) return MAG_OPTION_NCODE; + if (!strcasecmp (string, "nphot")) return MAG_OPTION_NPHOT; + if (!strcasecmp (string, "uc_dist")) return MAG_OPTION_UC_DIST; + if (!strcasecmp (string, "flags")) return MAG_OPTION_FLAGS; + if (!strcasecmp (string, "X11_SM_OBJ")) return MAG_OPTION_X11_SM_OBJ; + if (!strcasecmp (string, "X12_SM_OBJ")) return MAG_OPTION_X12_SM_OBJ; + if (!strcasecmp (string, "X22_SM_OBJ")) return MAG_OPTION_X22_SM_OBJ; + if (!strcasecmp (string, "E1_SM_OBJ")) return MAG_OPTION_E1_SM_OBJ; + if (!strcasecmp (string, "E2_SM_OBJ")) return MAG_OPTION_E2_SM_OBJ; + if (!strcasecmp (string, "X11_SH_OBJ")) return MAG_OPTION_X11_SH_OBJ; + if (!strcasecmp (string, "X12_SH_OBJ")) return MAG_OPTION_X12_SH_OBJ; + if (!strcasecmp (string, "X22_SH_OBJ")) return MAG_OPTION_X22_SH_OBJ; + if (!strcasecmp (string, "E1_SH_OBJ")) return MAG_OPTION_E1_SH_OBJ; + if (!strcasecmp (string, "E2_SH_OBJ")) return MAG_OPTION_E2_SH_OBJ; + if (!strcasecmp (string, "X11_SM_PSF")) return MAG_OPTION_X11_SM_PSF; + if (!strcasecmp (string, "X12_SM_PSF")) return MAG_OPTION_X12_SM_PSF; + if (!strcasecmp (string, "X22_SM_PSF")) return MAG_OPTION_X22_SM_PSF; + if (!strcasecmp (string, "E1_SM_PSF")) return MAG_OPTION_E1_SM_PSF; + if (!strcasecmp (string, "E2_SM_PSF")) return MAG_OPTION_E2_SM_PSF; + if (!strcasecmp (string, "X11_SH_PSF")) return MAG_OPTION_X11_SH_PSF; + if (!strcasecmp (string, "X12_SH_PSF")) return MAG_OPTION_X12_SH_PSF; + if (!strcasecmp (string, "X22_SH_PSF")) return MAG_OPTION_X22_SH_PSF; + if (!strcasecmp (string, "E1_SH_PSF")) return MAG_OPTION_E1_SH_PSF; + if (!strcasecmp (string, "E2_SH_PSF")) return MAG_OPTION_E2_SH_PSF; + // if (!strcasecmp (string, "F_AP_R5")) return MAG_OPTION_F_AP_R5; + // if (!strcasecmp (string, "F_ERR_AP_R5")) return MAG_OPTION_F_ERR_AP_R5; + // if (!strcasecmp (string, "F_STDEV_AP_R5")) return MAG_OPTION_F_STDEV_AP_R5; + // if (!strcasecmp (string, "F_FILL_AP_R5")) return MAG_OPTION_F_FILL_AP_R5; + // if (!strcasecmp (string, "F_AP_R6")) return MAG_OPTION_F_AP_R6; + // if (!strcasecmp (string, "F_ERR_AP_R6")) return MAG_OPTION_F_ERR_AP_R6; + // if (!strcasecmp (string, "F_STDEV_AP_R6")) return MAG_OPTION_F_STDEV_AP_R6; + // if (!strcasecmp (string, "F_FILL_AP_R6")) return MAG_OPTION_F_FILL_AP_R6; return MAG_OPTION_NONE; } - -# if (0) -int GetMagMode (char *string) { - - // these all imply PSF mags: - if (!strcasecmp (string, "inst")) return (MAG_INST); - if (!strcasecmp (string, "cat")) return (MAG_CAT); - if (!strcasecmp (string, "sys")) return (MAG_SYS); - if (!strcasecmp (string, "rel")) return (MAG_REL); - if (!strcasecmp (string, "cal")) return (MAG_CAL); - if (!strcasecmp (string, "ave")) return (MAG_AVE); - if (!strcasecmp (string, "ref")) return (MAG_REF); - if (!strcasecmp (string, "err")) return (MAG_ERR); - if (!strcasecmp (string, "aveerr")) return (MAG_AVE_ERR); - - // these are explicit PSF mags: - if (!strcasecmp (string, "psf_inst")) return (MAG_INST); - if (!strcasecmp (string, "psf_cat")) return (MAG_CAT); - if (!strcasecmp (string, "psf_sys")) return (MAG_SYS); - if (!strcasecmp (string, "psf_rel")) return (MAG_REL); - if (!strcasecmp (string, "psf_cal")) return (MAG_CAL); - if (!strcasecmp (string, "psf_ave")) return (MAG_AVE); - if (!strcasecmp (string, "psf_ref")) return (MAG_REF); - if (!strcasecmp (string, "psf_err")) return (MAG_ERR); - if (!strcasecmp (string, "psf_aveerr")) return (MAG_AVE_ERR); - - // these are aper mags: - if (!strcasecmp (string, "ap")) return (MAG_APER); - if (!strcasecmp (string, "aper")) return (MAG_APER); - if (!strcasecmp (string, "aperinst")) return (MAG_APER_INST); - if (!strcasecmp (string, "aper_inst")) return (MAG_APER_INST); - if (!strcasecmp (string, "aper_cat")) return (MAG_APER_CAT); - if (!strcasecmp (string, "aper_sys")) return (MAG_APER_SYS); - if (!strcasecmp (string, "aper_rel")) return (MAG_APER_REL); - if (!strcasecmp (string, "aper_cal")) return (MAG_APER_CAL); - if (!strcasecmp (string, "aper_ave")) return (MAG_APER_AVE); - if (!strcasecmp (string, "aper_ref")) return (MAG_APER_REF); - if (!strcasecmp (string, "aper_err")) return (MAG_APER_ERR); - - // these are kron mags: - if (!strcasecmp (string, "kron")) return (MAG_KRON); - if (!strcasecmp (string, "kron_inst")) return (MAG_KRON_INST); - if (!strcasecmp (string, "kron_cat")) return (MAG_KRON_CAT); - if (!strcasecmp (string, "kron_sys")) return (MAG_KRON_SYS); - if (!strcasecmp (string, "kron_rel")) return (MAG_KRON_REL); - if (!strcasecmp (string, "kron_cal")) return (MAG_KRON_CAL); - if (!strcasecmp (string, "kron_ave")) return (MAG_KRON_AVE); - if (!strcasecmp (string, "kron_ref")) return (MAG_KRON_REF); - if (!strcasecmp (string, "kron_err")) return (MAG_KRON_ERR); - if (!strcasecmp (string, "kroninst")) return (MAG_KRON_INST); - if (!strcasecmp (string, "kronerr")) return (MAG_KRON_ERR); - - // other fields: - if (!strcasecmp (string, "photflags")) return (MAG_PHOT_FLAGS); - if (!strcasecmp (string, "flags")) return (MAG_PHOT_FLAGS); - if (!strcasecmp (string, "chisq")) return (MAG_CHISQ); - if (!strcasecmp (string, "ncode")) return (MAG_NCODE); - if (!strcasecmp (string, "nphot")) return (MAG_NPHOT); - if (!strcasecmp (string, "stdev")) return (MAG_STDEV); - if (!strcasecmp (string, "min")) return (MAG_MIN); - if (!strcasecmp (string, "max")) return (MAG_MAX); - if (!strcasecmp (string, "ucdist")) return (MAG_UC_DIST); - if (!strcasecmp (string, "fluxpsf")) return (MAG_FLUX_PSF); - if (!strcasecmp (string, "fluxpsferr")) return (MAG_FLUX_PSF_ERR); - if (!strcasecmp (string, "fluxkron")) return (MAG_FLUX_KRON); - if (!strcasecmp (string, "fluxkronerr")) return (MAG_FLUX_KRON_ERR); - return (MAG_NONE); -} -# endif // field may be of the form mag:psf:inst:wrp. except for the first subword, @@ -285,4 +344,5 @@ fprintf (stderr, "ERROR: unknown mag/photocode argument %s\n", word); free (word); + return FALSE; } @@ -414,4 +474,10 @@ if (!strcasecmp (fieldName, "XFIX")) ESCAPE (MEAS_XFIX, OPIHI_FLT); if (!strcasecmp (fieldName, "YFIX")) ESCAPE (MEAS_YFIX, OPIHI_FLT); + if (!strcasecmp (fieldName, "XOFF_KH")) ESCAPE (MEAS_XOFF_KH, OPIHI_FLT); + if (!strcasecmp (fieldName, "YOFF_KH")) ESCAPE (MEAS_YOFF_KH, OPIHI_FLT); + if (!strcasecmp (fieldName, "XOFF_DCR")) ESCAPE (MEAS_XOFF_DCR, OPIHI_FLT); + if (!strcasecmp (fieldName, "YOFF_DCR")) ESCAPE (MEAS_YOFF_DCR, OPIHI_FLT); + if (!strcasecmp (fieldName, "ROFF_GAL")) ESCAPE (MEAS_ROFF_GAL, OPIHI_FLT); + if (!strcasecmp (fieldName, "DOFF_GAL")) ESCAPE (MEAS_DOFF_GAL, OPIHI_FLT); if (!strcasecmp (fieldName, "POS_SYS_ERR")) ESCAPE (MEAS_POS_SYS_ERR, OPIHI_FLT); if (!strcasecmp (fieldName, "XFIELD")) ESCAPE (MEAS_XFIELD, OPIHI_FLT); @@ -441,5 +507,47 @@ if (!strcasecmp (fieldName, "FLAT")) ESCAPE (MEAS_FLAT, OPIHI_FLT); if (!strcasecmp (fieldName, "CENTER_OFFSET")) ESCAPE (MEAS_CENTER_OFFSET, OPIHI_FLT); - // if (!strcasecmp (fieldName, "FLUX")) ESCAPE (MEAS_FLUX_PSF, OPIHI_FLT); + + // individual lensing measurements are not grouped by photcode: + if (!strcasecmp (fieldName, "X11_SM_OBJ")) ESCAPE (MEAS_X11_SM_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X12_SM_OBJ")) ESCAPE (MEAS_X12_SM_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X22_SM_OBJ")) ESCAPE (MEAS_X22_SM_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "E1_SM_OBJ")) ESCAPE (MEAS_E1_SM_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "E2_SM_OBJ")) ESCAPE (MEAS_E2_SM_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X11_SH_OBJ")) ESCAPE (MEAS_X11_SH_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X12_SH_OBJ")) ESCAPE (MEAS_X12_SH_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X22_SH_OBJ")) ESCAPE (MEAS_X22_SH_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "E1_SH_OBJ")) ESCAPE (MEAS_E1_SH_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "E2_SH_OBJ")) ESCAPE (MEAS_E2_SH_OBJ, OPIHI_FLT); + if (!strcasecmp (fieldName, "X11_SM_PSF")) ESCAPE (MEAS_X11_SM_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "X12_SM_PSF")) ESCAPE (MEAS_X12_SM_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "X22_SM_PSF")) ESCAPE (MEAS_X22_SM_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "E1_SM_PSF")) ESCAPE (MEAS_E1_SM_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "E2_SM_PSF")) ESCAPE (MEAS_E2_SM_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "X11_SH_PSF")) ESCAPE (MEAS_X11_SH_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "X12_SH_PSF")) ESCAPE (MEAS_X12_SH_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "X22_SH_PSF")) ESCAPE (MEAS_X22_SH_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "E1_SH_PSF")) ESCAPE (MEAS_E1_SH_PSF, OPIHI_FLT); + if (!strcasecmp (fieldName, "E2_SH_PSF")) ESCAPE (MEAS_E2_SH_PSF, OPIHI_FLT); + + if (!strcasecmp (fieldName, "E_BV")) ESCAPE (MEAS_E_BV , OPIHI_FLT); + if (!strcasecmp (fieldName, "E_BV_ERR")) ESCAPE (MEAS_E_BV_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "DISTANCE_MOD")) ESCAPE (MEAS_DISTANCE_MOD , OPIHI_FLT); + if (!strcasecmp (fieldName, "DISTANCE_MOD_ERR")) ESCAPE (MEAS_DISTANCE_MOD_ERR, OPIHI_FLT); + if (!strcasecmp (fieldName, "D_M")) ESCAPE (MEAS_DISTANCE_MOD , OPIHI_FLT); + if (!strcasecmp (fieldName, "D_M_ERR")) ESCAPE (MEAS_DISTANCE_MOD_ERR, OPIHI_FLT); + if (!strcasecmp (fieldName, "M_R")) ESCAPE (MEAS_M_R , OPIHI_FLT); + if (!strcasecmp (fieldName, "M_R_ERR")) ESCAPE (MEAS_M_R_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "FEH")) ESCAPE (MEAS_FEH , OPIHI_FLT); + if (!strcasecmp (fieldName, "FEH_ERR")) ESCAPE (MEAS_FEH_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "URA_GALMODEL")) ESCAPE (MEAS_URA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "UDEC_GALMODEL")) ESCAPE (MEAS_UDEC_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "RA_GALMODEL")) ESCAPE (MEAS_RA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "DEC_GALMODEL")) ESCAPE (MEAS_DEC_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "URA_GAL")) ESCAPE (MEAS_URA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "UDEC_GAL")) ESCAPE (MEAS_UDEC_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "RA_GAL")) ESCAPE (MEAS_RA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "DEC_GAL")) ESCAPE (MEAS_DEC_GALMODEL , OPIHI_FLT); + + // if (!strcasecmp (fieldName, "FLUX")) ESCAPE (MEAS_FLUX_PSF, OPIHI_FLT); // if (!strcasecmp (fieldName, "FLUX_ERR")) ESCAPE (MEAS_FLUX_PSF_ERR, OPIHI_FLT); // if (!strcasecmp (fieldName, "FLUX_PSF")) ESCAPE (MEAS_FLUX_PSF, OPIHI_FLT); @@ -447,5 +555,6 @@ // if (!strcasecmp (fieldName, "FLUX_KRON")) ESCAPE (MEAS_FLUX_KRON, OPIHI_FLT); // if (!strcasecmp (fieldName, "FLUX_KRON_ERR")) ESCAPE (MEAS_FLUX_KRON_ERR, OPIHI_FLT); - if (!strcasecmp (fieldName, "REF_COLOR")) ESCAPE (MEAS_REF_COLOR, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_BLUE")) ESCAPE (MEAS_REF_COLOR_BLUE, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_RED")) ESCAPE (MEAS_REF_COLOR_RED, OPIHI_FLT); // for words that don't parse, try a photcode @@ -522,5 +631,20 @@ if (!strcasecmp (fieldName, "EXTID_HI")) ESCAPE (AVE_EXTID_HI, OPIHI_INT); if (!strcasecmp (fieldName, "EXTID_LO")) ESCAPE (AVE_EXTID_LO, OPIHI_INT); - if (!strcasecmp (fieldName, "REF_COLOR")) ESCAPE (AVE_REF_COLOR, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_BLUE")) ESCAPE (AVE_REF_COLOR_BLUE, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_RED")) ESCAPE (AVE_REF_COLOR_RED, OPIHI_FLT); + + if (!strcasecmp (fieldName, "E_BV")) ESCAPE (AVE_E_BV , OPIHI_FLT); + if (!strcasecmp (fieldName, "E_BV_ERR")) ESCAPE (AVE_E_BV_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "DISTANCE_MOD")) ESCAPE (AVE_DISTANCE_MOD , OPIHI_FLT); + if (!strcasecmp (fieldName, "DISTANCE_MOD_ERR")) ESCAPE (AVE_DISTANCE_MOD_ERR, OPIHI_FLT); + if (!strcasecmp (fieldName, "M_R")) ESCAPE (AVE_M_R , OPIHI_FLT); + if (!strcasecmp (fieldName, "M_R_ERR")) ESCAPE (AVE_M_R_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "FEH")) ESCAPE (AVE_FEH , OPIHI_FLT); + if (!strcasecmp (fieldName, "FEH_ERR")) ESCAPE (AVE_FEH_ERR , OPIHI_FLT); + if (!strcasecmp (fieldName, "URA_GALMODEL")) ESCAPE (AVE_URA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "UDEC_GALMODEL")) ESCAPE (AVE_UDEC_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "RA_GALMODEL")) ESCAPE (AVE_RA_GALMODEL , OPIHI_FLT); + if (!strcasecmp (fieldName, "DEC_GALMODEL")) ESCAPE (AVE_DEC_GALMODEL , OPIHI_FLT); + // check for code:mode in photcode name @@ -652,5 +776,6 @@ if (!strcasecmp (fieldName, "NLINK_PHOTOM")) ESCAPE (IMAGE_NLINK_PHOTOM, OPIHI_INT); if (!strcasecmp (fieldName, "NLINK_ASTROM")) ESCAPE (IMAGE_NLINK_ASTROM, OPIHI_INT); - if (!strcasecmp (fieldName, "REF_COLOR")) ESCAPE (IMAGE_REF_COLOR, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_BLUE")) ESCAPE (IMAGE_REF_COLOR_BLUE, OPIHI_FLT); + if (!strcasecmp (fieldName, "REF_COLOR_RED")) ESCAPE (IMAGE_REF_COLOR_RED, OPIHI_FLT); // for words that don't parse, try a photcode Index: trunk/Ohana/src/libdvo/src/dvo_catalog.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_catalog.c (revision 37807) @@ -58,4 +58,5 @@ if (!strcasecmp (catformat, "PS1_V5")) return (DVO_FORMAT_PS1_V5); if (!strcasecmp (catformat, "PS1_REF")) return (DVO_FORMAT_PS1_REF); + if (!strcasecmp (catformat, "PS1_SIM")) return (DVO_FORMAT_PS1_SIM); return (DVO_FORMAT_UNDEF); } @@ -129,4 +130,5 @@ average->Nlensing = 0; average->Nlensobj = 0; + average->Nstarpar = 0; average->Nextend = 0; @@ -135,6 +137,6 @@ average->lensingOffset = -1; average->lensobjOffset = -1; + average->starparOffset = -1; average->extendOffset = -1; - average->paramsOffset = -1; average->refColorBlue = NAN; @@ -144,4 +146,5 @@ average->skycellID = 0; average->projectionID = 0; + average->dummy = 0; @@ -437,4 +440,25 @@ // init all data, or just catalog data +void dvo_starpar_init (StarPar *starpar) { + starpar->R = NAN; + starpar->D = NAN; + starpar->galLat = NAN; + starpar->galLon = NAN; + starpar->Ebv = NAN; + starpar->dEbv = NAN; + starpar->DistMag = NAN; + starpar->dDistMag = NAN; + starpar->M_r = NAN; + starpar->dM_r = NAN; + starpar->FeH = NAN; + starpar->dFeH = NAN; + starpar->uRA = NAN; + starpar->uDEC = NAN; + starpar->averef = -1; + starpar->objID = -1; + starpar->catID = -1; +} + +// init all data, or just catalog data void dvo_catalog_init (Catalog *catalog, int complete) { @@ -462,4 +486,6 @@ catalog[0].lensobj = NULL; + catalog[0].starpar = NULL; + catalog[0].averageT = NULL; catalog[0].measureT = NULL; @@ -477,4 +503,6 @@ catalog[0].Nlensobj = 0; + catalog[0].Nstarpar = 0; + catalog[0].Naves_disk = 0; catalog[0].Nmeas_disk = 0; @@ -485,4 +513,6 @@ catalog[0].Nlensobj_disk = 0; + catalog[0].Nstarpar_disk = 0; + catalog[0].Naves_off = 0; catalog[0].Nmeas_off = 0; @@ -492,4 +522,6 @@ catalog[0].Nlensing_off = 0; catalog[0].Nlensobj_off = 0; + + catalog[0].Nstarpar_off = 0; /* pointers to SPLIT data files */ @@ -499,4 +531,5 @@ catalog[0].lensing_catalog = NULL; catalog[0].lensobj_catalog = NULL; + catalog[0].starpar_catalog = NULL; /* pointers for data manipulation */ @@ -540,5 +573,10 @@ if (catalog[0].f == (FILE *) NULL) return (2); - fflush (catalog[0].f); + + if (fflush (catalog[0].f)) { + fprintf (stderr, "failed to flush file\n"); + return FALSE; + } + // fprintf (stderr, "unlocking: %s\n", catalog[0].filename); @@ -552,12 +590,16 @@ // closes f but does not set back to NULL - fclearlockfile (catalog[0].filename, catalog[0].f, catalog[0].lockmode, &dbstate); + if (!fclearlockfile (catalog[0].filename, catalog[0].f, catalog[0].lockmode, &dbstate)) { + fprintf (stderr, "failed to unlock or close file\n"); + return (0); + } if (catalog[0].catmode == DVO_MODE_SPLIT) { - if (catalog[0].measure_catalog) dvo_catalog_unlock (catalog[0].measure_catalog); - if (catalog[0].missing_catalog) dvo_catalog_unlock (catalog[0].missing_catalog); - if (catalog[0].secfilt_catalog) dvo_catalog_unlock (catalog[0].secfilt_catalog); - if (catalog[0].lensing_catalog) dvo_catalog_unlock (catalog[0].lensing_catalog); - if (catalog[0].lensobj_catalog) dvo_catalog_unlock (catalog[0].lensobj_catalog); + if (catalog[0].measure_catalog) { if (!dvo_catalog_unlock (catalog[0].measure_catalog)) { fprintf (stderr, "failed to unlock measures\n"); return (0); }} + if (catalog[0].missing_catalog) { if (!dvo_catalog_unlock (catalog[0].missing_catalog)) { fprintf (stderr, "failed to unlock missing\n"); return (0); }} + if (catalog[0].secfilt_catalog) { if (!dvo_catalog_unlock (catalog[0].secfilt_catalog)) { fprintf (stderr, "failed to unlock secfilt\n"); return (0); }} + if (catalog[0].lensing_catalog) { if (!dvo_catalog_unlock (catalog[0].lensing_catalog)) { fprintf (stderr, "failed to unlock lensing\n"); return (0); }} + if (catalog[0].lensobj_catalog) { if (!dvo_catalog_unlock (catalog[0].lensobj_catalog)) { fprintf (stderr, "failed to unlock lensobj\n"); return (0); }} + if (catalog[0].starpar_catalog) { if (!dvo_catalog_unlock (catalog[0].starpar_catalog)) { fprintf (stderr, "failed to unlock starpar\n"); return (0); }} } return (1); @@ -597,4 +639,5 @@ } + // NOTE: this only check if we can write a backup; it does not actually write the backup if (!check_file_access (catalog[0].filename, BACKUP, READWRITE, VERBOSE)) { fprintf (stderr, "no permission to access %s\n", catalog[0].filename); @@ -835,4 +878,9 @@ free (catalog[0].lensobj_catalog); } + if (catalog[0].starpar_catalog) { + free (catalog[0].starpar_catalog[0].filename); + dvo_catalog_free (catalog[0].starpar_catalog); + free (catalog[0].starpar_catalog); + } } dvo_catalog_free_data (catalog); @@ -873,4 +921,9 @@ catalog[0].Nlensobj = 0; catalog[0].lensobj = NULL; + } + if (catalog[0].starpar) { + free (catalog[0].starpar); + catalog[0].Nstarpar = 0; + catalog[0].starpar = NULL; } if (catalog[0].nOwn_t) { free (catalog[0].nOwn_t); catalog[0].nOwn_t = NULL; } @@ -985,4 +1038,14 @@ } } + if (catalog[0].starpar_catalog != NULL) { + if (catalog[0].Nstarpar_disk == 0) { + // need to relock (and re-open) file for close elsewhere + status = dvo_catalog_lock (catalog[0].starpar_catalog, lockmode); + } else { + if (!dvo_catalog_backup (catalog[0].starpar_catalog, FALSE)) { + return FALSE; + } + } + } } return TRUE; @@ -1036,4 +1099,9 @@ } } + if ((catalog[0].starpar_catalog != NULL) && (catalog[0].Nstarpar_disk > 0)) { + if (!dvo_catalog_unlink_backup (catalog[0].starpar_catalog, FALSE)) { + return FALSE; + } + } } return TRUE; Index: trunk/Ohana/src/libdvo/src/dvo_catalog_create.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_create.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_create.c (revision 37807) @@ -95,4 +95,16 @@ free (file); + /* define starpar catalog file */ + ALLOCATE (catalog[0].starpar_catalog, Catalog, 1); + dvo_catalog_init (catalog[0].starpar_catalog, TRUE); + + /* create basic data for starpar catalog file */ + gfits_create_header (&catalog[0].starpar_catalog[0].header); + ALLOCATE (catalog[0].starpar_catalog[0].filename, char, length); + sprintf (catalog[0].starpar_catalog[0].filename, "%s/%s.cpz", path, root); + file = filebasename (catalog[0].starpar_catalog[0].filename); + gfits_modify (&catalog[0].header, "STARPAR", "%s", 1, file); + free (file); + free (path); free (root); @@ -117,4 +129,8 @@ } if (dvo_catalog_lock (catalog[0].lensobj_catalog, catalog[0].lockmode) != DVO_CAT_OPEN_EMPTY) { + fprintf (stderr, "error with file lock\n"); + exit (2); + } + if (dvo_catalog_lock (catalog[0].starpar_catalog, catalog[0].lockmode) != DVO_CAT_OPEN_EMPTY) { fprintf (stderr, "error with file lock\n"); exit (2); @@ -146,4 +162,5 @@ ALLOCATE (catalog[0].lensing, Lensing, 1); ALLOCATE (catalog[0].lensobj, Lensobj, 1); + ALLOCATE (catalog[0].starpar, StarPar, 1); /* setup secondary filters to match photcodes: Index: trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_mef.c (revision 37807) @@ -4,9 +4,4 @@ off_t Nbytes; - off_t Naverage; - off_t Nmeasure; - off_t Nmissing; - off_t Nlensing; - off_t Nlensobj; off_t Nitems, Nexpect; int Nsecfilt; @@ -18,12 +13,4 @@ ftable.header = &header; - - /* get the components from the header */ - if (!gfits_scan (&catalog[0].header, "NSTARS", OFF_T_FMT, 1, &Naverage)) return (FALSE); - if (!gfits_scan (&catalog[0].header, "NMEAS", OFF_T_FMT, 1, &Nmeasure)) return (FALSE); - if (!gfits_scan (&catalog[0].header, "NMISS", OFF_T_FMT, 1, &Nmissing)) return (FALSE); - if (!gfits_scan (&catalog[0].header, "NSECFILT", "%d", 1, &Nsecfilt)) Nsecfilt = 0; - if (!gfits_scan (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, &Nlensing)) return (FALSE); - if (!gfits_scan (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, &Nlensobj)) return (FALSE); /* the OBJID is a counter that uniquely defines an average entry and never changes. if @@ -38,26 +25,34 @@ catalog[0].catID = 0; } + + // NSTARS, average, Naves_disk + +# define GET_TABLE_SIZES(HFIELD, CFIELD, DFIELD, REQUIRED) \ + off_t N##CFIELD; \ + if (REQUIRED) { \ + if (!gfits_scan (&catalog[0].header, HFIELD, OFF_T_FMT, 1, &N##CFIELD)) return (FALSE); \ + } else { \ + if (!gfits_scan (&catalog[0].header, HFIELD, OFF_T_FMT, 1, &N##CFIELD)) N##CFIELD = 0; \ + } \ + /* save the current number so we can do partial updates */ \ + catalog[0].DFIELD = N##CFIELD; \ + /* default values, but we will assign these a valid value before we exit (even if empty) */ \ + catalog[0].CFIELD = NULL; - /* save the current number so we can do partial updates */ - catalog[0].Naves_disk = Naverage; - catalog[0].Nmeas_disk = Nmeasure; - catalog[0].Nmiss_disk = Nmissing; - catalog[0].Nsecf_disk = Naverage * Nsecfilt; - catalog[0].Nlensing_disk = Nlensing; - catalog[0].Nlensobj_disk = Nlensobj; + /* get the components and sizes from the header */ + GET_TABLE_SIZES ("NSTARS", average, Naves_disk, TRUE); + GET_TABLE_SIZES ("NMEAS", measure, Nmeas_disk, TRUE); + GET_TABLE_SIZES ("NMISS", missing, Nmiss_disk, FALSE); + GET_TABLE_SIZES ("NLENSING", lensing, Nlensing_disk, FALSE); + GET_TABLE_SIZES ("NLENSOBJ", lensobj, Nlensobj_disk, FALSE); + GET_TABLE_SIZES ("NSTARPAR", starpar, Nstarpar_disk, FALSE); /** Nsecfilt is unusual: it does not list the number of data items in the table instead, the number of items is Nsecfilt * Naverage. **/ + + if (!gfits_scan (&catalog[0].header, "NSECFILT", "%d", 1, &Nsecfilt)) Nsecfilt = 0; catalog[0].Nsecfilt = Nsecfilt; - - /* default values, but we will assign these a valid value before we exit (even if empty) */ - catalog[0].average = NULL; - catalog[0].measure = NULL; - catalog[0].missing = NULL; + catalog[0].Nsecf_disk = Naverage * Nsecfilt; catalog[0].secfilt = NULL; - catalog[0].lensing = NULL; - catalog[0].lensobj = NULL; - - /* validate table mode ?*/ /*** Average Table ***/ @@ -259,4 +254,29 @@ } + /* read StarPar table header */ + if (!gfits_fread_header (catalog[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table starpar header"); + return (FALSE); + } + /* read StarPar table data */ + if (catalog[0].catflags & LOAD_STARPAR) { + if (!gfits_fread_ftable_data (catalog[0].f, &ftable, FALSE)) { + if (VERBOSE) fprintf (stderr, "can't read table starpar data"); + return (FALSE); + } + catalog[0].starpar = FtableToStarPar (&ftable, &catalog[0].Nstarpar, &catalog[0].catformat); + if (Nstarpar != catalog[0].Nstarpar_disk) { + fprintf (stderr, "Warning: mismatch between Nstarpar in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nstarpar, catalog[0].Nstarpar_disk); + } + catalog[0].Nstarpar = catalog[0].Nstarpar_disk; + catalog[0].Nstarpar_off = 0; + } else { + Nbytes = gfits_data_size (&header); + fseeko (catalog[0].f, Nbytes, SEEK_CUR); + ALLOCATE (catalog[0].starpar, StarPar, 1); + catalog[0].Nstarpar = 0; + catalog[0].Nstarpar_off = catalog[0].Nstarpar_disk; + } + return (TRUE); } @@ -309,4 +329,5 @@ gfits_modify (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, catalog[0].Nlensing); gfits_modify (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, catalog[0].Nlensobj); + gfits_modify (&catalog[0].header, "NSTARPAR", OFF_T_FMT, 1, catalog[0].Nstarpar); gfits_modify_alt (&catalog[0].header, "EXTEND", "%t", 1, TRUE); gfits_modify (&catalog[0].header, "OBJID", "%d", 1, catalog[0].objID); @@ -403,4 +424,24 @@ /* write out Lensobj table (convert to FITS table format) */ LensobjToFtable (&ftable, catalog[0].lensobj, catalog[0].Nlensobj, catalog[0].catformat); + if (!gfits_fwrite_Theader (catalog[0].f, &header)) { + fprintf (stderr, "can't write table header"); + goto failure; + } + if (!gfits_fwrite_table (catalog[0].f, &ftable)) { + fprintf (stderr, "can't write table data"); + goto failure; + } + gfits_free_table (&ftable); + gfits_free_header (&header); + + /* free temp storage */ + if (primary != NULL) { + free (primary); + free (secfilt); + } + return (TRUE); + + /* write out StarPar table (convert to FITS table format) */ + StarParToFtable (&ftable, catalog[0].starpar, catalog[0].Nstarpar, catalog[0].catformat); if (!gfits_fwrite_Theader (catalog[0].f, &header)) { fprintf (stderr, "can't write table header"); @@ -446,4 +487,6 @@ lensobj header lensobj table + starpar header + starpar table */ Index: trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_raw.c (revision 37807) @@ -92,4 +92,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -287,4 +288,5 @@ if (catalog[0].catformat == DVO_FORMAT_PS1_V5) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_V5"); if (catalog[0].catformat == DVO_FORMAT_PS1_REF) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_REF"); + if (catalog[0].catformat == DVO_FORMAT_PS1_SIM) gfits_modify (&catalog[0].header, "FORMAT", "%s", 1, "PS1_SIM"); /* rewind file pointers and truncate file */ @@ -391,4 +393,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -443,4 +446,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -498,4 +502,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -550,4 +555,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -606,4 +612,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: @@ -658,4 +665,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); default: Index: trunk/Ohana/src/libdvo/src/dvo_catalog_split.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_catalog_split.c (revision 37807) @@ -75,6 +75,62 @@ Matrix matrix; + // ohana_memcheck_func (TRUE); + + + if (0) { + // XXX a TEST : is the data corrupted or the file? + char name[1024]; + snprintf (name, 1024, "%s.phu.v2", catalog->filename); + FILE *f = fopen (name, "w"); + fwrite (catalog->header.buffer, sizeof(char), catalog->header.datasize, f); + // fprintf (stderr, "save %d bytes in %s (%zx)\n", (int) catalog->header.datasize, name, (size_t) catalog->header.buffer); + // fprintf (stderr, "first bytes are %c%c%c%c%c%c%c%c\n", + // catalog->header.buffer[0], catalog->header.buffer[1], catalog->header.buffer[2], catalog->header.buffer[3], + // catalog->header.buffer[4], catalog->header.buffer[5], catalog->header.buffer[6], catalog->header.buffer[7]); + fclose (f); + } + if (0) { + // XXX a TEST : is the data corrupted or the file? + char name[1024]; + snprintf (name, 1024, "%s.phu", catalog->filename); + FILE *f = fopen (name, "w"); + fwrite (catalog->header.buffer, sizeof(char), catalog->header.datasize, f); + // fprintf (stderr, "save %d bytes in %s (%zx)\n", (int) catalog->header.datasize, name, (size_t) catalog->header.buffer); + // fprintf (stderr, "first bytes are %c%c%c%c%c%c%c%c\n", + // catalog->header.buffer[0], catalog->header.buffer[1], catalog->header.buffer[2], catalog->header.buffer[3], + // catalog->header.buffer[4], catalog->header.buffer[5], catalog->header.buffer[6], catalog->header.buffer[7]); + fclose (f); + } +# if (0) + { + // XXX a TEST : is the data corrupted or the file? + char name[1024]; + snprintf (name, 1024, "%s.phu.v2", catalog->filename); + FILE *f = fopen (name, "w"); + fwrite (catalog->header.buffer, sizeof(char), catalog->header.datasize, f); + fprintf (stderr, "save %d bytes in %s (%zx)\n", (int) catalog->header.datasize, name, (size_t) catalog->header.buffer); + fprintf (stderr, "first bytes are %c%c%c%c%c%c%c%c\n", + catalog->header.buffer[0], catalog->header.buffer[1], catalog->header.buffer[2], catalog->header.buffer[3], + catalog->header.buffer[4], catalog->header.buffer[5], catalog->header.buffer[6], catalog->header.buffer[7]); + fclose (f); + } + { + // XXX a TEST : is the data corrupted or the file? + char name[1024]; + snprintf (name, 1024, "%s.phu.v3", catalog->filename); + FILE *f = fopen (name, "w"); + fwrite (catalog->header.buffer, sizeof(char), catalog->header.datasize, f); + fprintf (stderr, "save %d bytes in %s (%zx)\n", (int) catalog->header.datasize, name, (size_t) catalog->header.buffer); + fprintf (stderr, "first bytes are %c%c%c%c%c%c%c%c\n", + catalog->header.buffer[0], catalog->header.buffer[1], catalog->header.buffer[2], catalog->header.buffer[3], + catalog->header.buffer[4], catalog->header.buffer[5], catalog->header.buffer[6], catalog->header.buffer[7]); + fclose (f); + } +# endif /* rewind file pointers and truncate (file is still open) */ - fseeko (catalog->f, 0LL, SEEK_SET); + if (fseeko (catalog->f, 0LL, SEEK_SET)) { + perror ("dvo_catalog_save_subset: "); + fprintf (stderr, "failed to seek to beginning\n"); + } // write PHU header @@ -134,8 +190,9 @@ if (!gfits_scan (&catalog[0].header, name, "%s", 1, string)) { free (path); - // databases created prior to ~2014.07.01 did not have the LENSING or LENSOBJ paths in their headers. + // databases created prior to ~2014.07.01 did not have the LENSING, LENSOBJ, STARPAR paths in their headers. // in these cases, we do not try to lock or open the relevant file if (!strcmp (name, "LENSING")) return (DVO_CAT_OPEN_EMPTY); if (!strcmp (name, "LENSOBJ")) return (DVO_CAT_OPEN_EMPTY); + if (!strcmp (name, "STARPAR")) return (DVO_CAT_OPEN_EMPTY); return (DVO_CAT_OPEN_FAIL); } @@ -182,4 +239,5 @@ off_t Nlensing; off_t Nlensobj; + off_t Nstarpar; off_t Nitems; int status, Nsecfilt; @@ -196,5 +254,5 @@ /* get the components from the header - these duplicate information in the split files (NAXIS2) */ - // NSTARS, NMEAS, NMISS are required; NLENSING, NLENSOBJ are not (0 if not found) + // NSTARS, NMEAS, NMISS are required; NLENSING, NLENSOBJ, NSTARPAR are not (0 if not found) if (!gfits_scan (&catalog[0].header, "NSTARS", OFF_T_FMT, 1, &Naverage)) return (FALSE); if (!gfits_scan (&catalog[0].header, "NMEAS", OFF_T_FMT, 1, &Nmeasure)) return (FALSE); @@ -203,4 +261,5 @@ if (!gfits_scan (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, &Nlensing)) Nlensing = 0; if (!gfits_scan (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, &Nlensobj)) Nlensobj = 0; + if (!gfits_scan (&catalog[0].header, "NSTARPAR", OFF_T_FMT, 1, &Nstarpar)) Nstarpar = 0; /* the OBJID is a counter that uniquely defines an average entry and never changes. if @@ -223,4 +282,5 @@ catalog[0].Nlensing_disk = Nlensing; catalog[0].Nlensobj_disk = Nlensobj; + catalog[0].Nstarpar_disk = Nstarpar; /** Nsecfilt is unusual: it does not list the number of data items in the table @@ -235,4 +295,5 @@ catalog[0].lensing = NULL; catalog[0].lensobj = NULL; + catalog[0].starpar = NULL; /*** Average Table ***/ @@ -348,4 +409,5 @@ /*** Secfilt Table ***/ status = dvo_catalog_open_subcat (catalog, &catalog[0].secfilt_catalog, ftable.header, "SECFILT", VERBOSE); + // fprintf (stderr, "secfilt header %s : %zx\n", catalog[0].secfilt_catalog[0].filename, (size_t) catalog[0].secfilt_catalog[0].header.buffer); if (status == DVO_CAT_OPEN_FAIL) { return (FALSE); @@ -465,4 +527,44 @@ gfits_free_header (&header); + /*** StarPar Table ***/ + if (!(catalog[0].catflags & SKIP_STARPAR)) { + // unless we specify 'skip', we still need to load the + status = dvo_catalog_open_subcat (catalog, &catalog[0].starpar_catalog, ftable.header, "STARPAR", VERBOSE); + if (status == DVO_CAT_OPEN_FAIL) { + return (FALSE); + } + // the StarPar file need not exist, unless we expect data to exist + if ((status == DVO_CAT_OPEN_EMPTY) && (catalog[0].Nstarpar_disk > 0)) { + return (FALSE); + } + } + if ((status != DVO_CAT_OPEN_EMPTY) && (catalog[0].catflags & LOAD_STARPAR)) { + // only read the StarPar table if the file actually exists + /* read StarPar table data */ + if (!gfits_fread_ftable_data (catalog[0].starpar_catalog[0].f, &ftable, FALSE)) { + if (VERBOSE) fprintf (stderr, "can't read table starpar data\n"); + return (FALSE); + } + /* convert data format to internal : returns number of row read in Nstarpar */ + catalog[0].starpar = FtableToStarPar (&ftable, &Nstarpar, &catalog[0].catformat); + if (Nstarpar != catalog[0].Nstarpar_disk) { + fprintf (stderr, "Warning: mismatch between Nstarpar in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nstarpar, catalog[0].Nstarpar_disk); + } + catalog[0].Nstarpar = catalog[0].Nstarpar_disk; + catalog[0].Nstarpar_off = 0; + } else { + if (catalog[0].starpar_catalog) { + gfits_free_header (&catalog[0].starpar_catalog[0].header); + } else { + ALLOCATE (catalog[0].starpar_catalog, Catalog, 1); + dvo_catalog_init (catalog[0].starpar_catalog, TRUE); + } + gfits_create_header (&catalog[0].starpar_catalog[0].header); + ALLOCATE (catalog[0].starpar, StarPar, 1); + catalog[0].Nstarpar = 0; + catalog[0].Nstarpar_off = catalog[0].Nstarpar_disk; + } + gfits_free_header (&header); + return (TRUE); } @@ -473,5 +575,5 @@ off_t Nbytes; - off_t Naverage, Nexpect, Nitems, Nmeasure, Nmissing, Nlensing, Nlensobj; + off_t Naverage, Nexpect, Nitems, Nmeasure, Nmissing, Nlensing, Nlensobj, Nstarpar; Header header; FTable ftable; @@ -497,5 +599,5 @@ } /* read Average table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (catalog[0].f, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (catalog[0].f, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table average data"); return (FALSE); @@ -527,5 +629,5 @@ } /* read Secfilt table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, &ftable, start*catalog[0].Nsecfilt, catalog[0].Naverage*catalog[0].Nsecfilt)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start*catalog[0].Nsecfilt, catalog[0].Naverage*catalog[0].Nsecfilt)) { if (VERBOSE) fprintf (stderr, "can't read table measure data"); return (FALSE); @@ -563,5 +665,5 @@ } /* read Measure table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table measure data"); return (FALSE); @@ -594,5 +696,5 @@ } /* read Missing table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table missing data"); return (FALSE); @@ -628,5 +730,5 @@ } /* read Lensing table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table lensing data"); return (FALSE); @@ -659,5 +761,5 @@ } /* read Lensobj table data : format is irrelevant here */ - if (!gfits_fread_ftable_range (subcat[0].f, &ftable, start, Nrows)) { + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { if (VERBOSE) fprintf (stderr, "can't read table lensobj data"); return (FALSE); @@ -675,4 +777,35 @@ } + // XXX check the open status of the catalog + if (catalog[0].catflags & LOAD_STARPAR) { + + Catalog *subcat = catalog[0].starpar_catalog; + + /* move pointer past header -- must be already read (load_catalog) */ + Nbytes = subcat[0].header.datasize + gfits_data_size (&subcat[0].header); + fseeko (subcat[0].f, Nbytes, SEEK_SET); + + /* read StarPar table header */ + if (!gfits_fread_header (subcat[0].f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read table starpar header"); + return (FALSE); + } + /* read StarPar table data : format is irrelevant here */ + if (!gfits_fread_ftable_range (subcat[0].f, FALSE, &ftable, start, Nrows)) { + if (VERBOSE) fprintf (stderr, "can't read table starpar data"); + return (FALSE); + } + + /* convert data format to internal : returns number of row read in Nstarpar */ + catalog[0].starpar = FtableToStarPar (&ftable, &Nstarpar, &catalog[0].catformat); + if (Nstarpar != Nrows) { + // XXX this condition denotes the eof has been reached; not an error or a warning + // fprintf (stderr, "Warning: mismatch between Nstarpar in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nstarpar, Nrows); + } + gfits_free_header (&header); + catalog[0].Nstarpar = Nstarpar; + catalog[0].Nstarpar_off = start; + } + return (TRUE); } @@ -686,5 +819,5 @@ SecFilt *primary, *secfilt; int Nsecfilt; - off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new; + off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new, Nstarpar_disk_new; off_t first, start, Nrows; @@ -719,4 +852,5 @@ Nlensing_disk_new = MAX (catalog[0].Nlensing_disk, catalog[0].Nlensing + catalog[0].Nlensing_off); Nlensobj_disk_new = MAX (catalog[0].Nlensobj_disk, catalog[0].Nlensobj + catalog[0].Nlensobj_off); + Nstarpar_disk_new = MAX (catalog[0].Nstarpar_disk, catalog[0].Nstarpar + catalog[0].Nstarpar_off); /* make sure header is consistent with data */ @@ -727,4 +861,5 @@ gfits_modify (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, Nlensing_disk_new); gfits_modify (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, Nlensobj_disk_new); + gfits_modify (&catalog[0].header, "NSTARPAR", OFF_T_FMT, 1, Nstarpar_disk_new); gfits_modify_alt (&catalog[0].header, "EXTEND", "%t", 1, TRUE); gfits_modify (&catalog[0].header, "OBJID", "%d", 1, catalog[0].objID); @@ -845,6 +980,6 @@ } - /*** Lensing Table ***/ - if ((catalog[0].catflags & LOAD_LENSING) && (catalog[0].lensing != NULL)) { + /*** Lensing Table (optional, do not save if not loaded) ***/ + if ((catalog[0].catflags & LOAD_LENSING) && catalog[0].lensing_catalog->f && catalog[0].lensing) { first = 0; // first row in memory to write @@ -872,6 +1007,6 @@ } - /*** Lensobj Table ***/ - if ((catalog[0].catflags & LOAD_LENSOBJ) && (catalog[0].lensobj != NULL)) { + /*** Lensobj Table (optional, do not save if not loaded) ***/ + if ((catalog[0].catflags & LOAD_LENSOBJ) && catalog[0].lensobj_catalog->f && catalog[0].lensobj) { first = 0; // first row in memory to write @@ -893,4 +1028,31 @@ if (!dvo_catalog_save_subcat (catalog[0].lensobj_catalog, &ftable, start, Nrows, catalog[0].Nlensobj_disk, Nlensobj_disk_new)) { fprintf (stderr, "trouble writing Lensobj table\n"); + goto failure; + } + gfits_free_header (&header); + gfits_free_table (&ftable); + } + + /*** StarPar Table (optional, do not save if not loaded) ***/ + if ((catalog[0].catflags & LOAD_STARPAR) && catalog[0].starpar_catalog->f && catalog[0].starpar) { + + first = 0; // first row in memory to write + start = catalog[0].Nstarpar_off; // first disk row to write + Nrows = catalog[0].Nstarpar - first; + + assert (Nrows >= 0); + assert (first >= 0); + assert (first <= catalog[0].Nstarpar); + assert (catalog[0].Nstarpar_disk >= catalog[0].Nstarpar_off); + + // convert to external table format (note that the block above does not damage or free catalog.average) + if (!StarParToFtable (&ftable, &catalog[0].starpar[first], Nrows, catalog[0].catformat)) { + fprintf (stderr, "trouble converting format\n"); + goto failure; + } + + // write out StarPar table + if (!dvo_catalog_save_subcat (catalog[0].starpar_catalog, &ftable, start, Nrows, catalog[0].Nstarpar_disk, Nstarpar_disk_new)) { + fprintf (stderr, "trouble writing StarPar table\n"); goto failure; } @@ -925,5 +1087,5 @@ SecFilt *primary, *secfilt; int Nsecfilt; - off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new; + off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new, Nstarpar_disk_new; off_t first, start, Nrows; @@ -956,4 +1118,8 @@ if (catalog[0].Nlensobj_off > 0) { fprintf (stderr, "ERROR: only partial catalog (Lensobj) was loaded\n"); + goto failure; + } + if (catalog[0].Nstarpar_off > 0) { + fprintf (stderr, "ERROR: only partial catalog (StarPar) was loaded\n"); goto failure; } @@ -984,4 +1150,5 @@ Nlensing_disk_new = catalog[0].Nlensing; Nlensobj_disk_new = catalog[0].Nlensobj; + Nstarpar_disk_new = catalog[0].Nstarpar; /* make sure header is consistent with data */ @@ -992,4 +1159,5 @@ gfits_modify (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, Nlensing_disk_new); gfits_modify (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, Nlensobj_disk_new); + gfits_modify (&catalog[0].header, "NSTARPAR", OFF_T_FMT, 1, Nstarpar_disk_new); gfits_modify_alt (&catalog[0].header, "EXTEND", "%t", 1, TRUE); gfits_modify (&catalog[0].header, "OBJID", "%d", 1, catalog[0].objID); @@ -1158,4 +1326,31 @@ if (!dvo_catalog_save_subcat (catalog[0].lensobj_catalog, &ftable, start, Nrows, catalog[0].Nlensobj_disk, Nlensobj_disk_new)) { fprintf (stderr, "trouble writing Lensobj table\n"); + goto failure; + } + gfits_free_header (&header); + gfits_free_table (&ftable); + } + + /*** StarPar Table ***/ + if ((catalog[0].catflags & LOAD_STARPAR) && (catalog[0].starpar != NULL)) { + + first = 0; // first row in memory to write + start = catalog[0].Nstarpar_off; // first disk row to write + Nrows = catalog[0].Nstarpar - first; + + assert (Nrows >= 0); + assert (first >= 0); + assert (first <= catalog[0].Nstarpar); + assert (catalog[0].Nstarpar_disk >= catalog[0].Nstarpar_off); + + // convert to external table format + if (!StarParToFtable (&ftable, &catalog[0].starpar[first], Nrows, catalog[0].catformat)) { + fprintf (stderr, "trouble converting format\n"); + goto failure; + } + + // write out StarPar table + if (!dvo_catalog_save_subcat (catalog[0].starpar_catalog, &ftable, start, Nrows, catalog[0].Nstarpar_disk, Nstarpar_disk_new)) { + fprintf (stderr, "trouble writing StarPar table\n"); goto failure; } @@ -1189,5 +1384,5 @@ SecFilt *primary, *secfilt; int Nsecfilt; - off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new; + off_t Naves_disk_new, Nmeas_disk_new, Nmiss_disk_new, Nsecf_disk_new, Nlensing_disk_new, Nlensobj_disk_new, Nstarpar_disk_new; off_t first, start, Nrows; @@ -1221,4 +1416,5 @@ Nlensing_disk_new = MAX (catalog[0].Nlensing_disk, catalog[0].Nlensing + catalog[0].Nlensing_off); Nlensobj_disk_new = MAX (catalog[0].Nlensobj_disk, catalog[0].Nlensobj + catalog[0].Nlensobj_off); + Nstarpar_disk_new = MAX (catalog[0].Nstarpar_disk, catalog[0].Nstarpar + catalog[0].Nstarpar_off); /* make sure header is consistent with data */ @@ -1229,4 +1425,5 @@ gfits_modify (&catalog[0].header, "NLENSING", OFF_T_FMT, 1, Nlensing_disk_new); gfits_modify (&catalog[0].header, "NLENSOBJ", OFF_T_FMT, 1, Nlensobj_disk_new); + gfits_modify (&catalog[0].header, "NSTARPAR", OFF_T_FMT, 1, Nstarpar_disk_new); gfits_modify_alt (&catalog[0].header, "EXTEND", "%t", 1, TRUE); gfits_modify (&catalog[0].header, "OBJID", "%d", 1, catalog[0].objID); @@ -1399,4 +1596,32 @@ if (!dvo_catalog_save_subcat (catalog[0].lensobj_catalog, &ftable, start, Nrows, catalog[0].Nlensobj_disk, Nlensobj_disk_new)) { fprintf (stderr, "trouble writing Lensobj table\n"); + goto failure; + } + gfits_free_header (&header); + gfits_free_table (&ftable); + } + + /*** StarPar Table ***/ + if (catalog[0].starpar != NULL) { + + first = catalog[0].Nstarpar_disk - catalog[0].Nstarpar_off; // first row in memory to write + start = catalog[0].Nstarpar_disk; // first disk row to write + Nrows = catalog[0].Nstarpar - first; + + assert (Nrows >= 0); + assert (first >= 0); + assert (first <= catalog[0].Nstarpar); + assert (catalog[0].Nstarpar_disk >= catalog[0].Nstarpar_off); + + // convert to external table format + // XXX does catalog.starpar have averef correctly set up? + if (!StarParToFtable (&ftable, &catalog[0].starpar[first], Nrows, catalog[0].catformat)) { + fprintf (stderr, "trouble converting format\n"); + goto failure; + } + + // write out StarPar table + if (!dvo_catalog_save_subcat (catalog[0].starpar_catalog, &ftable, start, Nrows, catalog[0].Nstarpar_disk, Nstarpar_disk_new)) { + fprintf (stderr, "trouble writing StarPar table\n"); goto failure; } Index: trunk/Ohana/src/libdvo/src/dvo_convert.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert.c (revision 37807) @@ -43,4 +43,5 @@ CONVERT_FORMAT ("DVO_AVERAGE_PS1_V5", PS1_V5); CONVERT_FORMAT ("DVO_AVERAGE_PS1_REF", PS1_REF); + CONVERT_FORMAT ("DVO_AVERAGE_PS1_SIM", PS1_SIM); # undef CONVERT_FORMAT @@ -129,4 +130,5 @@ CONVERT_FORMAT ("DVO_AVERAGE_PS1_V5", PS1_V5, PS1_V5); CONVERT_FORMAT ("DVO_AVERAGE_PS1_REF", PS1_REF, PS1_REF); + CONVERT_FORMAT ("DVO_AVERAGE_PS1_SIM", PS1_SIM, PS1_SIM); # undef CONVERT_FORMAT @@ -165,4 +167,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE @@ -240,4 +243,5 @@ CONVERT_FORMAT ("DVO_MEASURE_PS1_V5", PS1_V5, PS1_V5, TRUE); CONVERT_FORMAT ("DVO_MEASURE_PS1_REF", PS1_REF, PS1_REF, FALSE); + CONVERT_FORMAT ("DVO_MEASURE_PS1_SIM", PS1_SIM, PS1_SIM, TRUE); # undef CONVERT_FORMAT @@ -254,5 +258,5 @@ case DVO_FORMAT_##FORMAT: { \ Measure_##TYPE *tmpMeasure; \ - myAssert (ABS_COORDS || average, "conversion to internal needs average table"); \ + myAssert (ABS_COORDS || average, "conversion from internal needs average table"); \ tmpMeasure = MeasureInternalTo_##TYPE (average, measure, Nmeasure); \ gfits_table_set_Measure_##TYPE (ftable, tmpMeasure, Nmeasure); \ @@ -278,4 +282,5 @@ FORMAT_CASE (PS1_V5, PS1_V5, TRUE); FORMAT_CASE (PS1_REF, PS1_REF, FALSE); + FORMAT_CASE (PS1_SIM, PS1_SIM, TRUE); # undef FORMAT_CASE @@ -350,4 +355,5 @@ CONVERT_FORMAT ("DVO_SECFILT_PS1_V5", PS1_V5, PS1_V5); CONVERT_FORMAT ("DVO_SECFILT_PS1_REF", PS1_REF, PS1_REF); + CONVERT_FORMAT ("DVO_SECFILT_PS1_SIM", PS1_SIM, PS1_SIM); # undef CONVERT_FORMAT @@ -386,4 +392,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE @@ -407,4 +414,16 @@ return (FALSE); } + + if (!strcmp (extname, "DVO_LENSING_PS1_V5") && (ftable[0].header[0].Naxis[0] == 128)) { + Lensing_PS1_V5_ALT *tmpLensing; + tmpLensing = gfits_table_get_Lensing_PS1_V5_ALT (ftable, Nlensing, NULL); + if (!tmpLensing) { + fprintf (stderr, "ERROR: failed to read lensings\n"); + exit (2); + } + lensing = Lensing_PS1_V5_ALT_ToInternal (tmpLensing, *Nlensing); + free (tmpLensing); + *format = DVO_FORMAT_PS1_V5; + return (lensing); } # define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ @@ -585,4 +604,102 @@ default: fprintf (stderr, "table format unknown (lensobj)\n"); + return (FALSE); + } + return (TRUE); +} + +/*** StarPar / FTable conversion functions ***/ + +StarPar *FtableToStarPar (FTable *ftable, off_t *Nstarpar, char *format) { + + StarPar *starpar; + char extname[80]; + + /* convert to the internal format */ + if (!gfits_scan (ftable[0].header, "EXTNAME", "%s", 1, extname)) { + fprintf (stderr, "EXTNAME missing for starpar table\n"); + return (FALSE); + } + +# define CONVERT_FORMAT(NAME, FORMAT, TYPE) \ + if (!strcmp (extname, NAME)) { \ + StarPar_##TYPE *tmpStarPar; \ + tmpStarPar = gfits_table_get_StarPar_##TYPE (ftable, Nstarpar, NULL); \ + if (!tmpStarPar) { \ + fprintf (stderr, "ERROR: failed to read starpar\n"); \ + exit (2); \ + } \ + starpar = StarPar_##TYPE##_ToInternal (tmpStarPar, *Nstarpar); \ + free (tmpStarPar); \ + *format = DVO_FORMAT_##FORMAT; \ + return (starpar); } + + if (!strcmp (extname, "DVO_STARPAR")) { + starpar = gfits_table_get_StarPar (ftable, Nstarpar, NULL); + if (!starpar) { + fprintf (stderr, "ERROR: failed to read starpar\n"); + exit (2); + } + *format = DVO_FORMAT_INTERNAL; + return (starpar); + } + + CONVERT_FORMAT ("DVO_STARPAR_PS1_SIM", PS1_SIM, PS1_SIM); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_REF", PS1_REF, PS1_REF); +// CONVERT_FORMAT ("DVO_STARPAR_ELIXIR", ELIXIR, Elixir); +// CONVERT_FORMAT ("DVO_STARPAR_LONEOS", LONEOS, Loneos); +// CONVERT_FORMAT ("DVO_STARPAR_PANSTARRS_DEV_0", PANSTARRS_DEV_0, Panstarrs_DEV_0); +// CONVERT_FORMAT ("DVO_STARPAR_PANSTARRS_DEV_1", PANSTARRS_DEV_1, Panstarrs_DEV_1); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_DEV_1", PS1_DEV_1, PS1_DEV_1); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_DEV_2", PS1_DEV_2, PS1_DEV_2); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_V1", PS1_V1, PS1_V1); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_V2", PS1_V2, PS1_V2); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_V3", PS1_V3, PS1_V3); +// CONVERT_FORMAT ("DVO_STARPAR_PS1_V4", PS1_V4, PS1_V4); + CONVERT_FORMAT ("DVO_STARPAR_PS1_V5", PS1_V5, PS1_V5); +// CONVERT_FORMAT ("DVO_STELLAR_PARAMS_PS1_V5", PS1_V5, PS1_V5); + // XXX fix in libauto: STELLAR_PARAMS to STARPAR +# undef CONVERT_FORMAT + + fprintf (stderr, "table format unknown: %s\n", extname); + + *Nstarpar = 0; + return (NULL); +} + +// StarParToFtable needs the Average since old formats stored starpar.dR,dD only +int StarParToFtable (FTable *ftable, StarPar *starpar, off_t Nstarpar, char format) { + +# define FORMAT_CASE(FORMAT, TYPE) \ + case DVO_FORMAT_##FORMAT: { \ + StarPar_##TYPE *tmpStarPar; \ + tmpStarPar = StarParInternalTo_##TYPE (starpar, Nstarpar); \ + gfits_table_set_StarPar_##TYPE (ftable, tmpStarPar, Nstarpar); \ + free (tmpStarPar); \ + break; } + + /* convert from the internal format */ + switch (format) { + case DVO_FORMAT_INTERNAL: { + gfits_table_set_StarPar (ftable, starpar, Nstarpar); + break; } + + FORMAT_CASE (PS1_SIM, PS1_SIM); +// FORMAT_CASE (PS1_REF, PS1_REF); +// FORMAT_CASE (ELIXIR, Elixir); +// FORMAT_CASE (LONEOS, Loneos); +// FORMAT_CASE (PANSTARRS_DEV_0, Panstarrs_DEV_0); +// FORMAT_CASE (PANSTARRS_DEV_1, Panstarrs_DEV_1); +// FORMAT_CASE (PS1_DEV_1, PS1_DEV_1); +// FORMAT_CASE (PS1_DEV_2, PS1_DEV_2); +// FORMAT_CASE (PS1_V1, PS1_V1); +// FORMAT_CASE (PS1_V2, PS1_V2); +// FORMAT_CASE (PS1_V3, PS1_V3); +// FORMAT_CASE (PS1_V4, PS1_V4); + FORMAT_CASE (PS1_V5, PS1_V5); +# undef FORMAT_CASE + + default: + fprintf (stderr, "table format unknown (starpar)\n"); return (FALSE); } @@ -643,15 +760,4 @@ return (TRUE); } - if (!strcmp (extname, "DVO_IMAGE")) { - Image *image; - image = gfits_table_get_Image (ftable, &Nimage, NULL); - if (!image) { - fprintf (stderr, "ERROR: failed to read images\n"); - exit (2); - } - *format = DVO_FORMAT_INTERNAL; - return (TRUE); - } - CONVERT_FORMAT ("DVO_IMAGE_ELIXIR", ELIXIR, Elixir); CONVERT_FORMAT ("DVO_IMAGE_LONEOS", LONEOS, Loneos); @@ -666,4 +772,5 @@ CONVERT_FORMAT ("DVO_IMAGE_PS1_V5", PS1_V5, PS1_V5); CONVERT_FORMAT ("DVO_IMAGE_PS1_REF", PS1_REF, PS1_REF); + CONVERT_FORMAT ("DVO_IMAGE_PS1_SIM", PS1_SIM, PS1_SIM); # undef CONVERT_FORMAT @@ -690,8 +797,4 @@ /* convert from the internal format */ switch (format) { - case DVO_FORMAT_INTERNAL: { - gfits_convert_Image ((Image *) ftable[0].buffer, sizeof(Image), Nimage); - break; } - FORMAT_CASE (ELIXIR, Elixir); FORMAT_CASE (LONEOS, Loneos); @@ -706,4 +809,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE @@ -744,10 +848,4 @@ /* convert from the internal format */ switch (format) { - case DVO_FORMAT_INTERNAL: { - for (i = 0; i < Nrow; i++) { - gfits_convert_Image ((Image *) vtable[0].buffer[i], sizeof(Image), 1); - } - return (TRUE); } - FORMAT_CASE (ELIXIR, Elixir); FORMAT_CASE (LONEOS, Loneos); @@ -762,4 +860,5 @@ FORMAT_CASE (PS1_V5, PS1_V5); FORMAT_CASE (PS1_REF, PS1_REF); + FORMAT_CASE (PS1_SIM, PS1_SIM); # undef FORMAT_CASE Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_1.c (revision 37807) @@ -219,5 +219,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -287,5 +290,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_2.c (revision 37807) @@ -214,5 +214,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -279,5 +282,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_3.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_3.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_DEV_3.c (revision 37807) @@ -18,5 +18,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -78,5 +81,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_REF.c (revision 37807) @@ -145,5 +145,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -208,5 +211,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_SIM.c (revision 37807) @@ -0,0 +1,429 @@ +# include + +/* convert PS1_SIM formats to internal formats */ + +Measure *Measure_PS1_SIM_ToInternal (Average *ave, Measure_PS1_SIM *in, off_t Nvalues) { + + off_t i; + Measure *out; + + ALLOCATE_ZERO (out, Measure, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_measure_init (&out[i]); + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].M = in[i].M; + out[i].dM = in[i].dM; + out[i].Mcal = in[i].Mcal; + out[i].dt = in[i].dt; + out[i].airmass = in[i].airmass; + out[i].az = in[i].az; + out[i].Xccd = in[i].Xccd; + out[i].Yccd = in[i].Yccd; + out[i].Xfix = in[i].Xfix; + out[i].Yfix = in[i].Yfix; + out[i].XoffKH = in[i].XoffKH; + out[i].YoffKH = in[i].YoffKH; + out[i].XoffDCR = in[i].XoffDCR; + out[i].YoffDCR = in[i].YoffDCR; + out[i].RoffGAL = in[i].RoffGAL; + out[i].DoffGAL = in[i].DoffGAL; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].detID = in[i].detID; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].imageID = in[i].imageID; + out[i].psfQF = in[i].psfQF; + out[i].psfQFperf = in[i].psfQFperf; + out[i].photcode = in[i].photcode; + out[i].dXccd = in[i].dXccd; + out[i].dYccd = in[i].dYccd; + out[i].dRsys = in[i].dRsys; + out[i].posangle = in[i].posangle; + out[i].pltscale = in[i].pltscale; + out[i].dbFlags = in[i].dbFlags; + out[i].photFlags = in[i].photFlags; + } + return (out); +} + +Measure_PS1_SIM *MeasureInternalTo_PS1_SIM (Average *ave, Measure *in, off_t Nvalues) { + + off_t i; + Measure_PS1_SIM *out; + + ALLOCATE_ZERO (out, Measure_PS1_SIM, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].M = in[i].M; + out[i].dM = in[i].dM; + out[i].Mcal = in[i].Mcal; + out[i].dt = in[i].dt; + out[i].airmass = in[i].airmass; + out[i].az = in[i].az; + out[i].Xccd = in[i].Xccd; + out[i].Yccd = in[i].Yccd; + out[i].Xfix = in[i].Xfix; + out[i].Yfix = in[i].Yfix; + out[i].XoffKH = in[i].XoffKH; + out[i].YoffKH = in[i].YoffKH; + out[i].XoffDCR = in[i].XoffDCR; + out[i].YoffDCR = in[i].YoffDCR; + out[i].RoffGAL = in[i].RoffGAL; + out[i].DoffGAL = in[i].DoffGAL; + out[i].t = in[i].t; + out[i].averef = in[i].averef; + out[i].detID = in[i].detID; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + out[i].imageID = in[i].imageID; + out[i].psfQF = in[i].psfQF; + out[i].psfQFperf = in[i].psfQFperf; + out[i].photcode = in[i].photcode; + out[i].dXccd = in[i].dXccd; + out[i].dYccd = in[i].dYccd; + out[i].dRsys = in[i].dRsys; + out[i].posangle = in[i].posangle; + out[i].pltscale = in[i].pltscale; + out[i].dbFlags = in[i].dbFlags; + out[i].photFlags = in[i].photFlags; + } + return (out); +} + +// 'primary' is needed to conform with the API for Loneos and Elixir, but is not used +Average *Average_PS1_SIM_ToInternal (Average_PS1_SIM *in, off_t Nvalues, SecFilt **primary) { + + off_t i; + Average *out; + + ALLOCATE_ZERO (out, Average, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_average_init (&out[i]); + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].uR = in[i].uR; + out[i].uD = in[i].uD; + out[i].duR = in[i].duR; + out[i].duD = in[i].duD; + out[i].P = in[i].P; + out[i].dP = in[i].dP; + + out[i].ChiSqAve = in[i].ChiSqAve; + out[i].ChiSqPM = in[i].ChiSqPM; + out[i].ChiSqPar = in[i].ChiSqPar; + out[i].Tmean = in[i].Tmean; + out[i].Trange = in[i].Trange; + + out[i].Npos = in[i].Npos; + + out[i].Nmeasure = in[i].Nmeasure; + out[i].Nstarpar = in[i].Nstarpar; + + out[i].measureOffset = in[i].measureOffset; + out[i].starparOffset = in[i].starparOffset; + + out[i].refColorBlue = in[i].refColorBlue; + out[i].refColorRed = in[i].refColorRed; + + out[i].flags = in[i].flags; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + } + return (out); +} + +// 'primary' is needed to conform with the API for Loneos and Elixir, but is not used +Average_PS1_SIM *AverageInternalTo_PS1_SIM (Average *in, off_t Nvalues, SecFilt *primary) { + + off_t i; + Average_PS1_SIM *out; + + ALLOCATE_ZERO (out, Average_PS1_SIM, Nvalues); + + for (i = 0; i < Nvalues; i++) { + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].dR = in[i].dR; + out[i].dD = in[i].dD; + out[i].uR = in[i].uR; + out[i].uD = in[i].uD; + out[i].duR = in[i].duR; + out[i].duD = in[i].duD; + out[i].P = in[i].P; + out[i].dP = in[i].dP; + + out[i].ChiSqAve = in[i].ChiSqAve; + out[i].ChiSqPM = in[i].ChiSqPM; + out[i].ChiSqPar = in[i].ChiSqPar; + out[i].Tmean = in[i].Tmean; + out[i].Trange = in[i].Trange; + + out[i].Npos = in[i].Npos; + + out[i].Nmeasure = in[i].Nmeasure; + out[i].Nstarpar = in[i].Nstarpar; + + out[i].measureOffset = in[i].measureOffset; + out[i].starparOffset = in[i].starparOffset; + + out[i].refColorBlue = in[i].refColorBlue; + out[i].refColorRed = in[i].refColorRed; + + out[i].flags = in[i].flags; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + } + return (out); +} + +SecFilt *SecFilt_PS1_SIM_ToInternal (SecFilt_PS1_SIM *in, off_t Nvalues) { + + off_t i; + SecFilt *out; + + ALLOCATE_ZERO (out, SecFilt, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_secfilt_init (&out[i]); + + out[i].M = in[i].M; + out[i].dM = in[i].dM; + out[i].Ncode = in[i].Ncode; + out[i].Nused = in[i].Nused; + out[i].flags = in[i].flags; + } + return (out); +} + +SecFilt_PS1_SIM *SecFiltInternalTo_PS1_SIM (SecFilt *in, off_t Nvalues) { + + off_t i; + SecFilt_PS1_SIM *out; + + ALLOCATE_ZERO (out, SecFilt_PS1_SIM, Nvalues); + + for (i = 0; i < Nvalues; i++) { + + out[i].M = in[i].M; + out[i].dM = in[i].dM; + out[i].Ncode = in[i].Ncode; + out[i].Nused = in[i].Nused; + out[i].flags = in[i].flags; + } + return (out); +} + +StarPar *StarPar_PS1_SIM_ToInternal (StarPar_PS1_SIM *in, off_t Nvalues) { + + off_t i; + StarPar *out; + + ALLOCATE_ZERO (out, StarPar, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_starpar_init (&out[i]); + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].galLat = in[i].galLat; + out[i].galLon = in[i].galLon; + + out[i].Ebv = in[i].Ebv ; + out[i].dEbv = in[i].dEbv ; + out[i].DistMag = in[i].DistMag ; + out[i].dDistMag = in[i].dDistMag; + out[i].M_r = in[i].M_r ; + out[i].dM_r = in[i].dM_r ; + out[i].FeH = in[i].FeH ; + out[i].dFeH = in[i].dFeH ; + out[i].uRA = in[i].uRA ; + out[i].uDEC = in[i].uDEC ; + + out[i].averef = in[i].averef; + out[i].objID = in[i].objID ; + out[i].catID = in[i].catID ; + } + return (out); +} + +StarPar_PS1_SIM *StarParInternalTo_PS1_SIM (StarPar *in, off_t Nvalues) { + + off_t i; + StarPar_PS1_SIM *out; + + ALLOCATE_ZERO (out, StarPar_PS1_SIM, Nvalues); + + for (i = 0; i < Nvalues; i++) { + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].galLat = in[i].galLat; + out[i].galLon = in[i].galLon; + + out[i].Ebv = in[i].Ebv ; + out[i].dEbv = in[i].dEbv ; + out[i].DistMag = in[i].DistMag ; + out[i].dDistMag = in[i].dDistMag; + out[i].M_r = in[i].M_r ; + out[i].dM_r = in[i].dM_r ; + out[i].FeH = in[i].FeH ; + out[i].dFeH = in[i].dFeH ; + out[i].uRA = in[i].uRA ; + out[i].uDEC = in[i].uDEC ; + + out[i].averef = in[i].averef; + out[i].objID = in[i].objID ; + out[i].catID = in[i].catID ; + } + return (out); +} + +# define RAW_IMAGE_NAME_LEN 117 + +Image *Image_PS1_SIM_ToInternal (Image_PS1_SIM *in, off_t Nvalues, off_t Nalloc) { + + off_t i; + Image *out; + + char *buffer; + ALLOCATE_ZERO (buffer, char, Nalloc); + out = (Image *) buffer; + // ALLOCATE_ZERO (out, Image, Nvalues); + + for (i = 0; i < Nvalues; i++) { + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; + + // RAW_IMAGE_NAME_LEN == DVO_IMAGE_NAME_LEN + strncpy (out[i].name, in[i].name, RAW_IMAGE_NAME_LEN - 1); + out[i].name[RAW_IMAGE_NAME_LEN - 1] = 0; // force termination + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit; + out[i].dapmifit = in[i].dapmifit; + out[i].Mcal = in[i].Mcal; + out[i].dMcal = in[i].dMcal; + out[i].Xm = in[i].Xm; + out[i].photcode = in[i].photcode; + out[i].exptime = in[i].exptime; + out[i].sidtime = in[i].sidtime; + out[i].latitude = in[i].latitude; + + out[i].RAo = in[i].RAo; + out[i].DECo = in[i].DECo; + out[i].Radius = in[i].Radius; + out[i].refColorBlue = in[i].refColorBlue; + out[i].refColorRed = in[i].refColorRed; + + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].ccdnum = in[i].ccdnum; + out[i].flags = in[i].flags; + out[i].imageID = in[i].imageID; + out[i].parentID = in[i].parentID; + out[i].externID = in[i].externID; + out[i].sourceID = in[i].sourceID; + + // as of 2011.02.03, the old Mx,My,..., Mxxxx,Myyyy have been deprecated and replaced + // with the following. (no real databases used those values -- see + // libdvo/doc/dvo-images.txt) + out[i].nLinkAstrom = in[i].nLinkAstrom; + out[i].nLinkPhotom = in[i].nLinkPhotom; + out[i].ubercalDist = in[i].ubercalDist; + out[i].dXpixSys = in[i].dXpixSys; + out[i].dYpixSys = in[i].dYpixSys; + out[i].dMagSys = in[i].dMagSys; + out[i].nFitAstrom = in[i].nFitAstrom; + out[i].nFitPhotom = in[i].nFitPhotom; + out[i].photom_map_id = in[i].photom_map_id; + out[i].astrom_map_id = in[i].astrom_map_id; + } + return (out); +} + +Image_PS1_SIM *ImageInternalTo_PS1_SIM (Image *in, off_t Nvalues) { + + off_t i; + Image_PS1_SIM *out; + + ALLOCATE_ZERO (out, Image_PS1_SIM, Nvalues); + + for (i = 0; i < Nvalues; i++) { + // this is only save because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + + // RAW_IMAGE_NAME_LEN == DVO_IMAGE_NAME_LEN + strncpy (out[i].name, in[i].name, RAW_IMAGE_NAME_LEN - 1); + out[i].name[RAW_IMAGE_NAME_LEN - 1] = 0; // force termination + + out[i].tzero = in[i].tzero; + out[i].nstar = in[i].nstar; + out[i].secz = in[i].secz; + out[i].NX = in[i].NX; + out[i].NY = in[i].NY; + out[i].apmifit = in[i].apmifit; + out[i].dapmifit = in[i].dapmifit; + out[i].Mcal = in[i].Mcal; + out[i].dMcal = in[i].dMcal; + out[i].Xm = in[i].Xm; + out[i].photcode = in[i].photcode; + out[i].exptime = in[i].exptime; + out[i].sidtime = in[i].sidtime; + out[i].latitude = in[i].latitude; + + out[i].RAo = in[i].RAo; + out[i].DECo = in[i].DECo; + out[i].Radius = in[i].Radius; + out[i].refColorBlue = in[i].refColorBlue; + out[i].refColorRed = in[i].refColorRed; + + out[i].detection_limit = in[i].detection_limit; + out[i].saturation_limit = in[i].saturation_limit; + out[i].cerror = in[i].cerror; + out[i].fwhm_x = in[i].fwhm_x; + out[i].fwhm_y = in[i].fwhm_y; + out[i].trate = in[i].trate; + out[i].ccdnum = in[i].ccdnum; + out[i].flags = in[i].flags; + out[i].imageID = in[i].imageID; + out[i].parentID = in[i].parentID; + out[i].externID = in[i].externID; + out[i].sourceID = in[i].sourceID; + + // as of 2011.02.03, the old Mx,My,..., Mxxxx,Myyyy have been deprecated and replaced + // with the following. (no real databases used those values -- see + // libdvo/doc/dvo-images.txt) + out[i].nLinkAstrom = in[i].nLinkAstrom; + out[i].nLinkPhotom = in[i].nLinkPhotom; + out[i].ubercalDist = in[i].ubercalDist; + out[i].dXpixSys = in[i].dXpixSys; + out[i].dYpixSys = in[i].dYpixSys; + out[i].dMagSys = in[i].dMagSys; + out[i].nFitAstrom = in[i].nFitAstrom; + out[i].nFitPhotom = in[i].nFitPhotom; + out[i].photom_map_id = in[i].photom_map_id; + out[i].astrom_map_id = in[i].astrom_map_id; + } + return (out); +} + Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V1.c (revision 37807) @@ -238,5 +238,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -302,5 +305,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V2.c (revision 37807) @@ -249,5 +249,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -313,5 +316,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V3.c (revision 37807) @@ -259,5 +259,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -323,5 +326,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V4.c (revision 37807) @@ -299,5 +299,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only save because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN @@ -363,5 +366,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only save because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN > DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_PS1_V5.c (revision 37807) @@ -191,4 +191,5 @@ out[i].Nlensing = in[i].Nlensing; out[i].Nlensobj = in[i].Nlensobj; + out[i].Nstarpar = in[i].Nstarpar; out[i].Nextend = in[i].Nextend; @@ -197,6 +198,6 @@ out[i].lensingOffset = in[i].lensingOffset; out[i].lensobjOffset = in[i].lensobjOffset; + out[i].starparOffset = in[i].starparOffset; out[i].extendOffset = in[i].extendOffset; - out[i].paramsOffset = in[i].paramsOffset; out[i].refColorBlue = in[i].refColorBlue; @@ -258,4 +259,5 @@ out[i].Nlensing = in[i].Nlensing; out[i].Nlensobj = in[i].Nlensobj; + out[i].Nstarpar = in[i].Nstarpar; out[i].Nextend = in[i].Nextend; @@ -264,6 +266,6 @@ out[i].lensingOffset = in[i].lensingOffset; out[i].lensobjOffset = in[i].lensobjOffset; + out[i].starparOffset = in[i].starparOffset; out[i].extendOffset = in[i].extendOffset; - out[i].paramsOffset = in[i].paramsOffset; out[i].refColorBlue = in[i].refColorBlue; @@ -663,4 +665,69 @@ } +StarPar *StarPar_PS1_V5_ToInternal (StarPar_PS1_V5 *in, off_t Nvalues) { + + off_t i; + StarPar *out; + + ALLOCATE_ZERO (out, StarPar, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_starpar_init (&out[i]); + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].galLat = in[i].galLat; + out[i].galLon = in[i].galLon; + + out[i].Ebv = in[i].Ebv ; + out[i].dEbv = in[i].dEbv ; + out[i].DistMag = in[i].DistMag ; + out[i].dDistMag = in[i].dDistMag; + out[i].M_r = in[i].M_r ; + out[i].dM_r = in[i].dM_r ; + out[i].FeH = in[i].FeH ; + out[i].dFeH = in[i].dFeH ; + out[i].uRA = in[i].uRA ; + out[i].uDEC = in[i].uDEC ; + + out[i].averef = in[i].averef; + out[i].objID = in[i].objID ; + out[i].catID = in[i].catID ; + } + return (out); +} + +StarPar_PS1_V5 *StarParInternalTo_PS1_V5 (StarPar *in, off_t Nvalues) { + + off_t i; + StarPar_PS1_V5 *out; + + ALLOCATE_ZERO (out, StarPar_PS1_V5, Nvalues); + + for (i = 0; i < Nvalues; i++) { + + out[i].R = in[i].R; + out[i].D = in[i].D; + out[i].galLat = in[i].galLat; + out[i].galLon = in[i].galLon; + + out[i].Ebv = in[i].Ebv ; + out[i].dEbv = in[i].dEbv ; + out[i].DistMag = in[i].DistMag ; + out[i].dDistMag = in[i].dDistMag; + out[i].M_r = in[i].M_r ; + out[i].dM_r = in[i].dM_r ; + out[i].FeH = in[i].FeH ; + out[i].dFeH = in[i].dFeH ; + out[i].uRA = in[i].uRA ; + out[i].uDEC = in[i].uDEC ; + + out[i].averef = in[i].averef; + out[i].objID = in[i].objID ; + out[i].catID = in[i].catID ; + } + return (out); +} + # define RAW_IMAGE_NAME_LEN 117 @@ -676,5 +743,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN == DVO_IMAGE_NAME_LEN @@ -741,5 +811,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only save because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN == DVO_IMAGE_NAME_LEN @@ -920,6 +991,6 @@ out[i].lensingOffset = in[i].lensingOffset; out[i].lensobjOffset = in[i].lensobjOffset; + out[i].starparOffset = in[i].starparOffset; out[i].extendOffset = in[i].extendOffset; - out[i].paramsOffset = in[i].paramsOffset; out[i].refColorBlue = in[i].refColorBlue; @@ -1194,2 +1265,55 @@ +Lensing *Lensing_PS1_V5_ALT_ToInternal (Lensing_PS1_V5_ALT *in, off_t Nvalues) { + + off_t i; + Lensing *out; + + ALLOCATE_ZERO (out, Lensing, Nvalues); + + for (i = 0; i < Nvalues; i++) { + dvo_lensing_init (&out[i]); + + out[i].X11_sm_obj = in[i].X11_sm_obj; + out[i].X12_sm_obj = in[i].X12_sm_obj; + out[i].X22_sm_obj = in[i].X22_sm_obj; + out[i].E1_sm_obj = in[i].E1_sm_obj; + out[i].E2_sm_obj = in[i].E2_sm_obj; + + out[i].X11_sh_obj = in[i].X11_sh_obj; + out[i].X12_sh_obj = in[i].X12_sh_obj; + out[i].X22_sh_obj = in[i].X22_sh_obj; + out[i].E1_sh_obj = in[i].E1_sh_obj; + out[i].E2_sh_obj = in[i].E2_sh_obj; + + out[i].X11_sm_psf = in[i].X11_sm_psf; + out[i].X12_sm_psf = in[i].X12_sm_psf; + out[i].X22_sm_psf = in[i].X22_sm_psf; + out[i].E1_sm_psf = in[i].E1_sm_psf; + out[i].E2_sm_psf = in[i].E2_sm_psf; + + out[i].X11_sh_psf = in[i].X11_sh_psf; + out[i].X12_sh_psf = in[i].X12_sh_psf; + out[i].X22_sh_psf = in[i].X22_sh_psf; + out[i].E1_sh_psf = in[i].E1_sh_psf; + out[i].E2_sh_psf = in[i].E2_sh_psf; + + out[i].F_ApR5 = in[i].F_ApR5; + out[i].dF_ApR5 = in[i].dF_ApR5; + out[i].sF_ApR5 = in[i].sF_ApR5; + out[i].fF_ApR5 = in[i].fF_ApR5; + + out[i].F_ApR6 = in[i].F_ApR6; + out[i].dF_ApR6 = in[i].dF_ApR6; + out[i].sF_ApR6 = in[i].sF_ApR6; + out[i].fF_ApR6 = in[i].fF_ApR6; + + out[i].detID = in[i].detID; + out[i].objID = in[i].objID; + out[i].catID = in[i].catID; + + out[i].averef = in[i].averef; + } + return (out); +} + Index: trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_elixir.c (revision 37807) @@ -210,5 +210,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -280,5 +283,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_loneos.c (revision 37807) @@ -188,5 +188,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -259,5 +262,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_0.c (revision 37807) @@ -227,5 +227,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -297,5 +300,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_convert_panstarrs_DEV_1.c (revision 37807) @@ -227,5 +227,8 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); + out[i].coords.mosaic = NULL; + out[i].coords.offsetMap = NULL; // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN @@ -297,5 +300,6 @@ for (i = 0; i < Nvalues; i++) { - memcpy (&out[i].coords, &in[i].coords, sizeof(Coords)); + // this is only safe because the initial 120 bytes in Coords match CoordsDisk + memcpy (&out[i].coords, &in[i].coords, sizeof(CoordsDisk)); // RAW_IMAGE_NAME_LEN < DVO_IMAGE_NAME_LEN Index: trunk/Ohana/src/libdvo/src/dvo_image.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_image.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_image.c (revision 37807) @@ -200,5 +200,5 @@ gfits_create_matrix (&db[0].header, &db[0].matrix); - gfits_table_set_Image (&db[0].ftable, NULL, 0); + gfits_table_set_Image (&db[0].ftable); gfits_modify (&db[0].header, "NIMAGES", "%d", 1, 0); @@ -207,5 +207,5 @@ dvo_image_createID (&db[0].header); - if (db[0].format == DVO_FORMAT_INTERNAL) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "INTERNAL"); + // if (db[0].format == DVO_FORMAT_INTERNAL) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "INTERNAL"); if (db[0].format == DVO_FORMAT_LONEOS) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "LONEOS"); if (db[0].format == DVO_FORMAT_ELIXIR) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "ELIXIR"); @@ -221,4 +221,5 @@ if (db[0].format == DVO_FORMAT_PS1_V5) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_V5"); if (db[0].format == DVO_FORMAT_PS1_REF) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_REF"); + if (db[0].format == DVO_FORMAT_PS1_SIM) gfits_modify (&db[0].header, "FORMAT", "%s", 1, "PS1_SIM"); return; @@ -267,2 +268,121 @@ } +int gfits_table_set_Image (FTable *ftable) { + + Header *header; + + header = ftable[0].header; + + gfits_table_mkheader_Image (header); + + /* create table */ + if (!gfits_create_table (header, ftable)) return (FALSE); + + return (TRUE); +} + +int gfits_table_mkheader_Image (Header *header) { + + /* create table header */ + if (!gfits_create_table_header (header, "BINTABLE", "DVO_IMAGE")) return (FALSE); + + /* define table layout */ + /** TABLE DEFINITION **/ + gfits_define_bintable_column (header, "D", "CRVAL1", "coordinate at reference pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "D", "CRVAL2", "coordinate at reference pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "CRPIX1", "coordinate of reference pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "CRPIX2", "coordinate of reference pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "CDELT1", "degrees per pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "CDELT2", "degrees per pixel", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "PC1_1", "rotation matrix", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "PC1_2", "rotation matrix", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "PC2_1", "rotation matrix", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "PC2_2", "rotation matrix", "", 1.0, 0.0); + gfits_define_bintable_column (header, "14E", "POLYTERMS", "higher order warping terms", "", 1.0, 0.0); + if (sizeof(void *) == 4) { + gfits_define_bintable_column (header, "J", "POINTER1", "", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "POINTER2", "", "", 1.0, 0.0); + } + if (sizeof(void *) == 8) { + gfits_define_bintable_column (header, "K", "POINTER1", "", "", 1.0, 0.0); + gfits_define_bintable_column (header, "K", "POINTER2", "", "", 1.0, 0.0); + } + gfits_define_bintable_column (header, "15A", "CTYPE", "coordinate type", "", 1.0, 0.0); + gfits_define_bintable_column (header, "A", "NPOLYTERMS", "order of polynomial", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "TZERO", "readout time (row 0)", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "NSTAR", "number of stars on image", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "SECZ", "airmass", "mag", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "NX", "image width", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "NY", "image height", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "APMIFIT", "aperture correction", "mag", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "DAPMIFIT", "apmifit error", "mag", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "MCAL", "calibration mag", "mag", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "DMCAL", "error on Mcal", "mag", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "XM", "image chisq", "10*log(value)", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "PHOTCODE", "identifier for CCD,", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "EXPTIME", "exposure time", "seconds", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "ST", "sidereal time of exposure", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "LAT", "observatory latitude", "degrees", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "RA_CENTER", "image center", "degrees", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "DEC_CENTER", "image center", "degrees", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "RADIUS", "image radius", "degrees", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "REF_COLOR_BLUE", "median astrometry ref color", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "REF_COLOR_RED", "median astrometry ref color", "", 1.0, 0.0); + gfits_define_bintable_column (header, "117A", "NAME", "name of original image ", "", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "DETECTION_LIMIT", "detection limit", "10*mag", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "SATURATION_LIMIT", "saturation limit", "10*mag", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "CERROR", "astrometric error", "50*arcsec", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "FWHM_X", "PSF x width", "25*arcsec", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "FWHM_Y", "PSF y width", "25*arcsec", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "TRATE", "scan rate", "100 usec/pixel", 1.0, 0.0); + gfits_define_bintable_column (header, "B", "CCDNUM", "CCD ID number", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "FLAGS", "image quality flags", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "IMAGE_ID", "internal image ID", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "PARENT_ID", "associated ref image", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "EXTERN_ID", "external image ID", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "SOURCE_ID", "analysis source ID", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "NLINK_ASTROM", "mean number of matched measurements for astrometry", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "NLINK_PHOTOM", "mean number of matched measurements for astrometry", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "UBERCAL_DIST", "distance to nearest ubercal image", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "XPIX_SYS_ERR", "systematic astrometry error in X", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "YPIX_SYS_ERR", "systematic astrometry error in Y", "", 1.0, 0.0); + gfits_define_bintable_column (header, "E", "MAG_SYS_ERR", "systematic photometry error", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "N_FIT_ASTROM", "number of stars used for astrometry cal", "", 1.0, 0.0); + gfits_define_bintable_column (header, "I", "N_FIT_PHOTOM", "number of stars used for photometry cal", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "PHOTOM_MAP_ID", "reference to 2D zero point map", "", 1.0, 0.0); + gfits_define_bintable_column (header, "J", "ASTROM_MAP_ID", "reference to 2D astrometry map", "", 1.0, 0.0); + if (sizeof(void *) == 4) { + gfits_define_bintable_column (header, "J", "POINTER3", "", "", 1.0, 0.0); + } + if (sizeof(void *) == 8) { + gfits_define_bintable_column (header, "K", "POINTER3", "", "", 1.0, 0.0); + } + + return (TRUE); +} + +/* return internal structure representation */ +Image *gfits_table_get_Image (FTable *ftable, off_t *Ndata, char *swapped) { + + int Ncols; + Image *data; + + Ncols = ftable[0].header[0].Naxis[0]; + if (Ncols != sizeof(Image)) { + fprintf (stderr, "ERROR: mis-match in table size: width is %d but should be %d bytes\n", Ncols, (int) sizeof(Image)); + return NULL; + } + + *Ndata = ftable[0].header[0].Naxis[1]; + data = (Image *) ftable[0].buffer; + + if (!swapped) { + myAbort ("invalid to call this without suppying 'swapped'"); + } + if (*swapped == FALSE) { + myAbort ("invalid for table to be swapped"); + } + + return (data); +} + Index: trunk/Ohana/src/libdvo/src/dvo_image_raw.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_image_raw.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_image_raw.c (revision 37807) @@ -61,4 +61,5 @@ if (db[0].format == DVO_FORMAT_PS1_V5) ImageSize = sizeof(Image_PS1_V5); if (db[0].format == DVO_FORMAT_PS1_REF) ImageSize = sizeof(Image_PS1_REF); + if (db[0].format == DVO_FORMAT_PS1_SIM) ImageSize = sizeof(Image_PS1_SIM); /* check that filesize makes sense */ @@ -92,4 +93,5 @@ if (db[0].format == DVO_FORMAT_PS1_V5) gfits_table_mkheader_Image_PS1_V5 (&db[0].theader); if (db[0].format == DVO_FORMAT_PS1_REF) gfits_table_mkheader_Image_PS1_REF (&db[0].theader); + if (db[0].format == DVO_FORMAT_PS1_SIM) gfits_table_mkheader_Image_PS1_SIM (&db[0].theader); /* read data from file */ Index: trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c (revision 37807) @@ -1788,2 +1788,36 @@ m2 = (Ns2 == -1) ? average[0].M : secfilt[Ns2].M; */ + +/*** lensing *********************************************************************************/ + +# define LENSFIELD(NAME) \ + float LensValue_##NAME (PhotCode *code, Lensobj *lensobj) { \ + if (code == NULL) return NAN; \ + int Ns = photcodes->hashNsec[code->code]; \ + float value = (Ns == -1) ? NAN : lensobj[Ns].NAME; \ + return (value); \ +} + +LENSFIELD(X11_sm_obj); +LENSFIELD(X12_sm_obj); +LENSFIELD(X22_sm_obj); +LENSFIELD(E1_sm_obj); +LENSFIELD(E2_sm_obj); + +LENSFIELD(X11_sh_obj); +LENSFIELD(X12_sh_obj); +LENSFIELD(X22_sh_obj); +LENSFIELD(E1_sh_obj); +LENSFIELD(E2_sh_obj); + +LENSFIELD(X11_sm_psf); +LENSFIELD(X12_sm_psf); +LENSFIELD(X22_sm_psf); +LENSFIELD(E1_sm_psf); +LENSFIELD(E2_sm_psf); + +LENSFIELD(X11_sh_psf); +LENSFIELD(X12_sh_psf); +LENSFIELD(X22_sh_psf); +LENSFIELD(E1_sh_psf); +LENSFIELD(E2_sh_psf); Index: trunk/Ohana/src/libdvo/src/dvo_tiny_values.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_tiny_values.c (revision 37807) @@ -20,4 +20,6 @@ measureT[0].dM = measure[0].dM; measureT[0].airmass = measure[0].airmass; + measureT[0].RoffGAL = measure[0].RoffGAL; + measureT[0].DoffGAL = measure[0].DoffGAL; measureT[0].Xccd = measure[0].Xccd; measureT[0].Yccd = measure[0].Yccd; Index: trunk/Ohana/src/libdvo/src/dvo_util.c =================================================================== --- trunk/Ohana/src/libdvo/src/dvo_util.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/dvo_util.c (revision 37807) @@ -145,10 +145,5 @@ if ((image->externID == externID) && (image->sourceID == sourceID)) { BuildChipMatch(dvoConfig->images, dvoConfig->nImages); - if (FindMosaicForImage(dvoConfig->images, dvoConfig->nImages, i)) { - return image; - } else { - fprintf(stderr, "FindMosaicForImage failed\n"); - return NULL; - } + return image; } } Index: trunk/Ohana/src/libdvo/src/flatcorr_io.c =================================================================== --- trunk/Ohana/src/libdvo/src/flatcorr_io.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/flatcorr_io.c (revision 37807) @@ -183,4 +183,11 @@ } +// eddie uses cell_x = (int) (X + X_PAD) / (base + X_PAD) +// where base = 600, X_PAD = 8, Y_PAD = 10. denominator is thus (608,610), whic is (4880/8),(4864/8) +// these match Eddie's NCELL_X,Y / CHIP_DX,DY, so I think we are good with this: + +# define X_PAD 8 +# define Y_PAD 10 + float FlatCorrectionOffset (FlatCorrectionTable *flatcorr, int ID, int X, int Y) { @@ -191,6 +198,6 @@ // convert X,Y to Xbin, Ybin: - int Xbin = MAX(MIN(X * flatcorr->image[seq].Nx / flatcorr->image[seq].DX, flatcorr->image[seq].Nx - 1), 0); - int Ybin = MAX(MIN(Y * flatcorr->image[seq].Ny / flatcorr->image[seq].DY, flatcorr->image[seq].Ny - 1), 0); + int Xbin = MAX(MIN((X + X_PAD) * flatcorr->image[seq].Nx / flatcorr->image[seq].DX, flatcorr->image[seq].Nx - 1), 0); + int Ybin = MAX(MIN((Y + Y_PAD) * flatcorr->image[seq].Ny / flatcorr->image[seq].DY, flatcorr->image[seq].Ny - 1), 0); // XXX warn if X,Y are out of range? not super important unless *way* out of range Index: trunk/Ohana/src/libdvo/src/galaxy_model.c =================================================================== --- trunk/Ohana/src/libdvo/src/galaxy_model.c (revision 37807) +++ trunk/Ohana/src/libdvo/src/galaxy_model.c (revision 37807) @@ -0,0 +1,148 @@ +# include + +// velocity / distance / proper-motion scale: +static double iFkap = 0.001 / 4.74047; // (km/sec/kpc) / (arcsec/ year) + +double A_oort = 0.0; // km/sec/kpc +double B_oort = 0.0; // km/sec/kpc +double U_sol = 0.0; // km/sec +double V_sol = 0.0; // km/sec +double W_sol = 0.0; // km/sec + +int InitGalaxyModel (char *version) { + + // galaxy model parameters from Siegfreid Roeser: + if (!strcmp(version, "ROESER")) { + A_oort = +14.50; // km/sec/kpc + B_oort = -13.00; // km/sec/kpc + U_sol = 9.44; // km/sec + V_sol = 11.90; // km/sec + W_sol = 7.20; // km/sec + return TRUE; + } + if (!strcmp(version, "FEAST-HIPPARCOS")) { + // http://arxiv.org/pdf/astro-ph/9706293v1.pdf + A_oort = +14.82; // km/sec/kpc + B_oort = -12.37; // km/sec/kpc + U_sol = 9.32; // km/sec + V_sol = 11.18; // km/sec + W_sol = 7.61; // km/sec + return TRUE; + } + return FALSE; +} + +int GalaxyMotionModel (double *uL_gal, double *uB_gal, double L, double B) { + double Lrad = L * RAD_DEG; + double Brad = B * RAD_DEG; + *uL_gal = (A_oort * cos(2.0*Lrad) + B_oort) * cos(Brad) * iFkap; + *uB_gal = -0.5*A_oort * sin(2.0*Lrad) * sin(Brad*2.0) * iFkap; + return TRUE; +} + +int GalaxyMotionModel_radians (double *uL_gal, double *uB_gal, double Lrad, double Brad) { + *uL_gal = (A_oort * cos(2.0*Lrad) + B_oort) * cos(Brad) * iFkap; + *uB_gal = -0.5*A_oort * sin(2.0*Lrad) * sin(Brad*2.0) * iFkap; + return TRUE; +} + +int SolarMotionModel (double *uL_sol, double *uB_sol, double L, double B, double distance) { + double Lrad = L * RAD_DEG; + double Brad = B * RAD_DEG; + *uL_sol = (U_sol * sin(Lrad) - V_sol * cos(Lrad)) * iFkap / distance; + *uB_sol = ((U_sol * cos(Lrad) + V_sol * sin(Lrad))*sin(Brad) - W_sol*cos(Brad)) * iFkap / distance; + return TRUE; +} + +int SolarMotionModel_radians (double *uL_sol, double *uB_sol, double Lrad, double Brad, double distance) { + *uL_sol = (U_sol * sin(Lrad) - V_sol * cos(Lrad)) * iFkap / distance; + *uB_sol = ((U_sol * cos(Lrad) + V_sol * sin(Lrad))*sin(Brad) - W_sol*cos(Brad)) * iFkap / distance; + return TRUE; +} + +// Use the inverse transform: CELESTIAL->GALACTIC for uL,uB -> uR,uD +int TransformProperMotion (double *uR, double *uD, double uL, double uB, double R, double D, CoordTransform *transform) { + + double Rrad = R*RAD_DEG; + double Drad = D*RAD_DEG; + + // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf + double C1 = + cos(Drad)*transform->cos_phi + + sin(Drad)*cos(Rrad)*transform->sin_phi_sin_Xo - + sin(Drad)*sin(Rrad)*transform->sin_phi_cos_Xo; + + double C2 = + - cos(Rrad)*transform->sin_phi_cos_Xo + - sin(Rrad)*transform->sin_phi_sin_Xo; + + double cosBinv = 1.0 / sqrt(C1*C1 + C2*C2); + + // XXX add errors : I need to be able to choose the stars based on the error distribution + *uR = cosBinv * (C1 * uL - C2 * uB); + *uD = cosBinv * (C1 * uB + C2 * uL); + + return TRUE; +} +int TransformProperMotion_radians (double *uR, double *uD, double uL, double uB, double Rrad, double Drad, CoordTransform *transform) { + + // C1, C2 are from http://arxiv.org/pdf/1306.2945v2.pdf + double C1 = + cos(Drad)*transform->cos_phi + + sin(Drad)*cos(Rrad)*transform->sin_phi_sin_Xo - + sin(Drad)*sin(Rrad)*transform->sin_phi_cos_Xo; + + double C2 = + - cos(Rrad)*transform->sin_phi_cos_Xo + - sin(Rrad)*transform->sin_phi_sin_Xo; + + double cosBinv = 1.0 / sqrt(C1*C1 + C2*C2); + + // XXX add errors : I need to be able to choose the stars based on the error distribution + *uR = cosBinv * (C1 * uL - C2 * uB); + *uD = cosBinv * (C1 * uB + C2 * uL); + + return TRUE; +} + +// I am using the galactic rotation and solar motion model to predict the reflex proper +// motion of stars in the database. + +/***** + +From Siegfreid Roeser: + +Dear Gene, + + please find attached our approach for galactic rotation and + solar motion that we assumed in our approach last fall. +(Hope I understood correctly what you requested at the +telecon). If not, please come back. + + fkap = 4.74047d-0 {km/s/kpc}/{mas/y} + + gl,gb is galactic longitude resp. latitude + + + PARAMETER (distance= 1.d0,Aoort=14.5d0 ,Boort=-13.d0 ) + kpc km/s/kpc km/s/kpc + PARAMETER (Usol = 9.44d0, Vsol =11.90d0, Wsol = 7.20d0) + km/s +C +C Gal. Rotation +C + proper motions in longitude/latitude + + emulgal = (AOORT*dcos(2.d0*gl)+BOORT)*dcos(gb)/fkap + emubgal = -0.5d0*AOORT*dsin(2.d0*gl)*dsin(2.d0*gb)/fkap +c +c Solar Motion +cc + emulsol = Usol*dsin(gl) - Vsol*dcos(gl) + emulsol = emulsol/fkap/distance + emubsol = Usol*dsin(gb)*dcos(gl)+Vsol*dsin(gb)*dsin(gl) + emubsol = emubsol -Wsol*dcos(gb) + emubsol = emubsol/fkap/distance + + +******/ Index: trunk/Ohana/src/libdvo/src/mosaic_astrom.c =================================================================== --- trunk/Ohana/src/libdvo/src/mosaic_astrom.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/mosaic_astrom.c (revision 37807) @@ -1,6 +1,11 @@ # include + +// the block of code below was used to define the chip->mosaic relationship before this +// was assigned to the coords value off_t getDISentry (e_time start, e_time stop, e_time *startMos, off_t *indexMos, off_t Nmosaic); void SortDISindex (e_time *S, off_t *I, off_t N); + +# if (0) /* chip-match table: j = ChipMatch[i], images[j] is DIS for images[i] WRP */ @@ -93,5 +98,6 @@ } -int BuildChipMatch_old (Image *images, off_t Nimages) { +// this is a very old version of BuildChipMatch which does not use bisection +int BuildChipMatch_nosort (Image *images, off_t Nimages) { off_t i, j, NDIS; @@ -140,5 +146,6 @@ } -int BuildChipMatch (Image *images, off_t Nimages) { +// this is an old version of BuildChipMatch which stores the results here in a local static variable +int BuildChipMatch_static (Image *images, off_t Nimages) { off_t i, j, NDIS; @@ -198,8 +205,61 @@ return (TRUE); } - +# endif + +// XXX Note that this links (images[i].coords.mosaic, images[i].parent) will break +// if we reallocate the Image array in the middle of program +int BuildChipMatch (Image *images, off_t Nimages) { + + off_t i, j, NDIS; + + off_t Ndis = 0; + off_t *DISentry = NULL; + e_time *DIStzero = NULL; + + if (DISentry != NULL) free (DISentry); + if (DIStzero != NULL) free (DIStzero); + + // allocate containers for DIS indexing + Ndis = 0; + NDIS = 100; + ALLOCATE (DISentry, off_t, NDIS); + ALLOCATE (DIStzero, e_time, NDIS); + + // find all DIS images, save tzero (& photcode?) + for (i = 0; i < Nimages; i++) { + if (strcmp(&images[i].coords.ctype[4], "-DIS")) continue; + DISentry[Ndis] = i; + DIStzero[Ndis] = images[i].tzero; + Ndis ++; + if (Ndis >= NDIS) { + NDIS += 100; + REALLOCATE (DISentry, off_t, NDIS); + REALLOCATE (DIStzero, e_time, NDIS); + } + } + + // sort the index, start, and stop by the start times: + SortDISindex (DIStzero, DISentry, Ndis); + + /* find all matched WRP images */ + for (i = 0; i < Nimages; i++) { + images[i].parent = NULL; // reset to NULL + images[i].coords.mosaic = NULL; // reset to NULL + if (strcmp(&images[i].coords.ctype[4], "-WRP")) continue; // only define link for WRP coords + + j = getDISentry (images[i].tzero, images[i].tzero + (int) images[i].exptime, DIStzero, DISentry, Ndis); + if (j == -1) { + fprintf (stderr, "WARNING: can't find matching mosaic \n"); + continue; + } + if (j >= Nimages) myAbort("invalid DIS entry"); + images[i].parent = &images[j]; + images[i].coords.mosaic = &images[j].coords; + } + return (TRUE); +} + +// use bisection to find the overlapping mosaic off_t getDISentry (e_time start, e_time stop, e_time *startMos, off_t *indexMos, off_t Nmosaic) { - - // use bisection to find the overlapping mosaic off_t Nlo, Nhi, N; Index: trunk/Ohana/src/libdvo/src/skyregion_gsc.c =================================================================== --- trunk/Ohana/src/libdvo/src/skyregion_gsc.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/skyregion_gsc.c (revision 37807) @@ -204,5 +204,5 @@ for (i = 0; i < Nregions; i++) { strncpy (temp, &buffer[i*48], 48); - temp[49] = 0; + temp[48] = 0; hstgsc_hms_to_deg (&Rmin, &Rmax, &Dmin, &Dmax, &temp[7]); if (Dmax < Dmin) SWAP (Dmin, Dmax); Index: trunk/Ohana/src/libdvo/src/skyregion_ops.c =================================================================== --- trunk/Ohana/src/libdvo/src/skyregion_ops.c (revision 37729) +++ trunk/Ohana/src/libdvo/src/skyregion_ops.c (revision 37807) @@ -16,4 +16,22 @@ if (c) return (c+1); return name; +} + +/* find region which corresponds to the given index */ +SkyList *SkyRegionByIndex (SkyTable *table, int index) { + + SkyList *list; + + ALLOCATE (list, SkyList, 1); + ALLOCATE (list[0].regions, SkyRegion *, 1); + ALLOCATE (list[0].filename, char *, 1); + list[0].Nregions = 0; + list[0].ownElements = FALSE; // this list is only holding a view to the elements + strcpy (list[0].hosts, table[0].hosts); + + list[0].regions[0] = &table[0].regions[index]; + list[0].filename[0] = table[0].filename[index]; + list[0].Nregions = 1; + return (list); } @@ -250,15 +268,13 @@ } - if( Rmax > (Rmin+180.)) { - - double temp = Rmax; - Rmax = Rmin; - Rmin = temp; - } - list = SkyListByBounds (table, depth, Rmin, Rmax, Dmin, Dmax); - return (list); -} - + + // fprintf (stderr, "%s : %f - %f, %f - %f : %d\n", + // image->name, Rmin, Rmax, Dmin, Dmax, (int) list->Nregions); + + return (list); +} + +// a region like -10 10 0 20 will become 350 10 0 20 and return all regions on both sides of 0,360 SkyList *SkyListByBounds (SkyTable *table, int depth, double Rmin, double Rmax, double Dmin, double Dmax) {