Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 37261) +++ trunk/Ohana/src/relastro/Makefile (revision 37807) @@ -60,5 +60,4 @@ $(SRC)/CoordOps.$(ARCH).o \ $(SRC)/extra.$(ARCH).o \ -$(SRC)/FixProblemImages.$(ARCH).o \ $(SRC)/StarMaps.$(ARCH).o \ $(SRC)/high_speed_catalogs.$(ARCH).o \ @@ -75,6 +74,14 @@ $(SRC)/relastro_parallel_regions.$(ARCH).o \ $(SRC)/MeanPosIO.$(ARCH).o \ +$(SRC)/MeasPosIO.$(ARCH).o \ $(SRC)/share_mean_pos.$(ARCH).o \ +$(SRC)/share_meas_pos.$(ARCH).o \ $(SRC)/share_images_pos.$(ARCH).o \ +$(SRC)/share_icrf_obj.$(ARCH).o \ +$(SRC)/ICRF.$(ARCH).o \ +$(SRC)/ICRFobjIO.$(ARCH).o \ +$(SRC)/FrameCorrection.$(ARCH).o \ +$(SRC)/FrameCorrectionIO.$(ARCH).o \ +$(SRC)/FrameCorrectionUtils.$(ARCH).o \ $(SRC)/ImagePosIO.$(ARCH).o \ $(SRC)/ImageTable.$(ARCH).o \ @@ -82,4 +89,5 @@ $(SRC)/indexCatalogs.$(ARCH).o \ $(SRC)/syncfile.$(ARCH).o \ +$(SRC)/client_logger.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o @@ -117,5 +125,4 @@ $(SRC)/relastro_objects.$(ARCH).o \ $(SRC)/save_catalogs.$(ARCH).o \ -$(SRC)/CoordOps.$(ARCH).o \ $(SRC)/extra.$(ARCH).o \ $(SRC)/high_speed_catalogs.$(ARCH).o \ @@ -128,16 +135,8 @@ $(SRC)/plot_scatter.$(ARCH).o \ $(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/ICRF.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o \ $(SRC)/syncfile.$(ARCH).o \ $(SRC)/BrightCatalog.$(ARCH).o - -#$(SRC)/FitChip.$(ARCH).o \ -#$(SRC)/FitMosaic.$(ARCH).o \ -#$(SRC)/UpdateChips.$(ARCH).o \ -#$(SRC)/UpdateMosaic.$(ARCH).o \ -#$(SRC)/UpdateSimple.$(ARCH).o \ -#$(SRC)/write_coords.$(ARCH).o \ -#$(SRC)/FixProblemImages.$(ARCH).o \ -#$(SRC)/StarMaps.$(ARCH).o \ $(RELASTRO_CLIENT): $(INC)/relastro.h Index: trunk/Ohana/src/relastro/doc/pv2.txt =================================================================== --- trunk/Ohana/src/relastro/doc/pv2.txt (revision 37261) +++ trunk/Ohana/src/relastro/doc/pv2.txt (revision 37807) @@ -1,2 +1,20 @@ + +2014.10.28 + +Going over the sequence of steps in relastro. I'm concerned that the +iteration is not modifying the measure.R,D values, preventing any +actual iterative improvements. + +* UpdateObjects : takes measure.R,D to calculate average.R,D and + motion terms + +* ImageOps : + * getImageRaw : takes measure.X,Y, applies image.coords to get (L,M), + (P,Q), etc, calculates + * getImageRef : takes average.R,D, applies image.coords to get + (P,Q), (L,M), etc + * FitChip : fit raw.X,Y to ref.L,M, plus residual map (dX,dY) + + *** modified raw.L,M values are then dropped and not used to update measure.R.D 2014.07.03 Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 37261) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 37807) @@ -34,9 +34,17 @@ typedef struct { - float R; - float D; + double R; + double D; unsigned int objID; unsigned int catID; } MeanPos; + +typedef struct { + double R; + double D; + unsigned int objID; + unsigned int catID; + unsigned int imageID; +} MeasPos; typedef struct { @@ -47,4 +55,6 @@ int nFitAstrom; int flags; + float refColorBlue; + float refColorRed; } ImagePos; @@ -142,4 +152,28 @@ } FitStats; +typedef struct { + double scale; + double **Roff; + double **Doff; + double *dR; + int *Nra; + int Ndec; +} FrameCorrectionType; + +typedef struct { + // AstromOffsetMap *map; + Coords *coords; // carries a pointer to the AstromOffsetMap + FrameCorrectionType *frame; +} FrameCorrectionSet; + + +typedef struct { + double *Rave; + double *Dave; + double *dRoff; + double *dDoff; + int Nicrfobj; +} ICRFobj; + /* global variables set in parameter file */ # define DVO_MAX_PATH 1024 @@ -172,4 +206,5 @@ double MIN_ERROR; +int IMFIT_TOO_FEW; // need more than this number of stars to fit an image int IMFIT_CLIP_NITER; // number of clipping iterations to perform in FitChip double IMFIT_CLIP_NSIGMA; // number of sigma to clip in FitChip @@ -188,6 +223,13 @@ int USE_FIXED_PIXCOORDS; +int USE_GALAXY_MODEL; +int USE_ICRF_CORRECT; + +int USE_ICRF_LOCAL; +int USE_ICRF_SHFIT; +int USE_ICRF_POLE; int RESET; +int RESET_IMAGES; int RESET_BAD_IMAGES; int NLOOP; @@ -204,4 +246,5 @@ int PLOTDELAY; int CHIPORDER; +int CHIPMAP; int MaxDensityUse; @@ -242,4 +285,5 @@ int CLIP_THRESH; int USE_BASIC_CHECK; +int APPLY_OFFSETS; // in parallel-regions mode, do not launch UpdateObjectOffsets unless -apply-offsets is called int ExcludeBogus; @@ -397,5 +441,5 @@ void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords); void initObjectData (Catalog *catalog, int Ncatalog); -int UpdateObjects (Catalog *catalog, int Ncatalog); +int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop); int UpdateSimple (Catalog *catalog, int Ncatalog); int UpdateChips (Catalog *catalog, int Ncatalog); @@ -417,5 +461,7 @@ StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode); +int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode); off_t getImageByID (off_t ID); +int updateImageRaw (Catalog *catalog, int Ncatalog, off_t im); Mosaic *getmosaics (off_t *N); @@ -447,15 +493,8 @@ int relastroGetVisual(void); -int FixProblemImages (SkyList *skylist); int *getCatlist (int *N, off_t im); -void initCoords (void); -void getOffsets (double *dPos, off_t *nPos, off_t N); -void saveOffsets (double dPos, off_t nPos, off_t N); -void setBadCoords (off_t N); -int badCoords (off_t N); -Coords *getCoords (off_t N); -int saveCoords (Coords *coords, off_t N); -void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im); +void SaveCoords (Coords *tgt, Coords *src); +void RestoreCoords (Coords *tgt, Coords *src, Image *image); int high_speed_catalogs (SkyTable *sky, SkyList *skylist, int hostID, char *hostpath); @@ -468,4 +507,5 @@ int checkStarMap(int N); int createStarMap (Catalog *catalog, int Ncatalog); +int printStarMap(int N, char *filename); int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit); @@ -519,4 +559,12 @@ MeanPos *MeanPosLoad(char *filename, off_t *nmeanpos); +int share_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int slurp_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int set_meas_pos (MeasPos *measpos, Average *average, MeasureTiny *measure); +MeasPos *merge_meas_pos (MeasPos *target, int *ntarget, MeasPos *source, int Nsource); + +int MeasPosSave(char *filename, MeasPos *measpos, off_t Nmeaspos); +MeasPos *MeasPosLoad(char *filename, off_t *nmeaspos); + int indexCatalogs (Catalog *catalog, int Ncatalog); int catID_and_objID_to_seq (int catID, int objID, int *catSeq, off_t *objSeq); @@ -528,5 +576,5 @@ int share_image_pos (RegionHostTable *regionHosts, int nloop); -int slurp_image_pos (RegionHostTable *regionHosts, int nloop); +int slurp_image_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); ImagePos *merge_image_pos (ImagePos *target, int *ntarget, ImagePos *source, int Nsource); int set_image_pos (ImagePos *image_pos, Image *image); @@ -547,2 +595,66 @@ int isGPC1stack (int photcode); int isGPC1warp (int photcode); + +int save_astrom_table (); +AstromOffsetTable *get_astrom_table (); +void put_astrom_table (AstromOffsetTable *myTable); + +int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts); + +// ICRF QSOs : these must be marked in the database (flag on average, flag on measure) + +void ICRFinit (); +int ICRFsave (int cat, int ave, int meas); +int ICRFdata (int n, int *cat, int *ave, int *meas); +int ICRFmax (); +int select_catalog_ICRF (Catalog *catalog, int Ncatalog); + +void ICRFobjFree (ICRFobj *icrfobj); +ICRFobj *get_ICRF_data (Catalog *catalog, int Ncatalog); +ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source); +ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop); +int share_icrf_obj (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); + +int ICRFobjSave(char *filename, ICRFobj *icrfobj); +ICRFobj *ICRFobjLoad(char *filename); + +int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit); + +FrameCorrectionType *FrameCorrectionInit (double scale); +void FrameCorrectionFree (FrameCorrectionType *frame); + +int FrameCorrectionParallelMaster (RegionHostTable *regionHosts); +int FrameCorrectionParallelSlave (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop); +int FrameCorrectionSerial (Catalog *catalog, int Ncatalog); + +FrameCorrectionSet *FrameCorrectionMeasure (ICRFobj *icrfobj); + +int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD); + +Coords *FrameCorrectionImageToMap (Header *header, Matrix *matrix, Coords *coords, int raDirection); +int FrameCorrectionMapToImage (Header *header, Matrix *matrix, Coords *coords, int raDirection); +int FrameCorrectionSHtoImage (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection); +FrameCorrectionType *FrameCorrectionImageToSH (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection); +int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD); +void FrameCorrectionFree (FrameCorrectionType *frame); +FrameCorrectionType *FrameCorrectionInit (double scale); +void FrameCorrectionSetFree (FrameCorrectionSet *set); +FrameCorrectionSet *FrameCorrectionSetInit (); + +int FrameCorrectionSetSave(char *filename, FrameCorrectionSet *set); +FrameCorrectionSet *FrameCorrectionSetLoad(char *filename); + +int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc); +int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords); +int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, Coords *coords); + +int dump_stardata_pts (StarData *raw, int Npts, char *filename); +void printNcatTotal (); + +int HarvestRegionHosts (RegionHostTable *regionHosts); + +void lockUpdateChips (); +void unlockUpdateChips (); + +int client_logger_init (char *dirname); +int client_logger_message (char *format,...); Index: trunk/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 37261) +++ trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 37807) @@ -43,11 +43,8 @@ if (!gfits_fread_header (f, &header)) { if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); - fclose (f); goto escape; } if (!gfits_fread_matrix (f, &matrix, &header)) { if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); - gfits_free_header (&header); - fclose (f); goto escape; } @@ -70,4 +67,6 @@ GET_COLUMN(R, "RA", double); GET_COLUMN(D, "DEC", double); + GET_COLUMN(RoffGAL, "ROFF_GAL", float); + GET_COLUMN(DoffGAL, "DOFF_GAL", float); GET_COLUMN(M, "MAG_SYS", float); GET_COLUMN(Mcal, "MAG_CAL", float); @@ -94,4 +93,6 @@ measure[i].R = R[i]; measure[i].D = D[i]; + measure[i].RoffGAL = RoffGAL[i]; + measure[i].DoffGAL = DoffGAL[i]; measure[i].M = M[i]; measure[i].Mcal = Mcal[i]; @@ -115,4 +116,6 @@ free (R ); free (D ); + free (RoffGAL ); + free (DoffGAL ); free (M ); free (Mcal ); @@ -310,5 +313,5 @@ if (catalog) free (catalog); - fclose (f); + if (f) fclose (f); return NULL; } @@ -343,6 +346,8 @@ gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); - gfits_define_bintable_column (&theader, "E", "RA_OFF", "ra offset", "arcsec", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "DEC_OFF", "dec offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "RA", "ra", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "dec", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "ROFF_GAL", "ra offset", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "DOFF_GAL", "dec offset", "arcsec", 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_SYS", "magnitude (sys)", NULL, 1.0, 0.0); gfits_define_bintable_column (&theader, "E", "MAG_CAL", "magnitude (cal)", NULL, 1.0, 0.0); @@ -368,4 +373,6 @@ double *R ; ALLOCATE (R , double, catalog->Nmeasure); double *D ; ALLOCATE (D , double, catalog->Nmeasure); + float *RoffGAL ; ALLOCATE (RoffGAL , float, catalog->Nmeasure); + float *DoffGAL ; ALLOCATE (DoffGAL , float, catalog->Nmeasure); float *M ; ALLOCATE (M , float, catalog->Nmeasure); float *Mcal ; ALLOCATE (Mcal , float, catalog->Nmeasure); @@ -390,4 +397,6 @@ R[i] = measure[i].R ; D[i] = measure[i].D ; + RoffGAL[i] = measure[i].RoffGAL ; + DoffGAL[i] = measure[i].DoffGAL ; M[i] = measure[i].M ; Mcal[i] = measure[i].Mcal ; @@ -411,4 +420,6 @@ gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "ROFF_GAL", RoffGAL, catalog->Nmeasure); + gfits_set_bintable_column (&theader, &ftable, "DOFF_GAL", DoffGAL, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_SYS", M, catalog->Nmeasure); gfits_set_bintable_column (&theader, &ftable, "MAG_CAL", Mcal, catalog->Nmeasure); @@ -430,4 +441,6 @@ free (R ); free (D ); + free (RoffGAL ); + free (DoffGAL ); free (M ); free (Mcal ); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 37261) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 37807) @@ -25,4 +25,5 @@ if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NSIGMA", "%lf", 0, &IMFIT_CLIP_NSIGMA)) IMFIT_CLIP_NSIGMA = 3.0; if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01; + if (!ScanConfig (config, "RELASTRO_IMFIT_TOO_FEW", "%d", 0, &IMFIT_TOO_FEW)) IMFIT_TOO_FEW = 5; if (!ScanConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN)) PM_DT_MIN = 0.25; @@ -36,4 +37,9 @@ if (!ScanConfig (config, "USE_FIXED_PIXCOORDS", "%d", 0, &USE_FIXED_PIXCOORDS)) USE_FIXED_PIXCOORDS = FALSE; + if (!ScanConfig (config, "USE_GALAXY_MODEL", "%d", 0, &USE_GALAXY_MODEL)) USE_GALAXY_MODEL = FALSE; + if (!ScanConfig (config, "USE_ICRF_CORRECT", "%d", 0, &USE_ICRF_CORRECT)) USE_ICRF_CORRECT = FALSE; + if (!ScanConfig (config, "USE_ICRF_LOCAL", "%d", 0, &USE_ICRF_LOCAL)) USE_ICRF_LOCAL = FALSE; + if (!ScanConfig (config, "USE_ICRF_SHFIT", "%d", 0, &USE_ICRF_SHFIT)) USE_ICRF_SHFIT = FALSE; + if (!ScanConfig (config, "USE_ICRF_POLE", "%d", 0, &USE_ICRF_POLE)) USE_ICRF_POLE = FALSE; // force CATDIR to be absolute (so parallel mode will work) Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 37261) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 37807) @@ -14,5 +14,5 @@ for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) { - + // measure the scatter for the unmarked (good) measurements with dM < limit // SIGMA_LIM here is redundant with ImageOps.c:915 @@ -49,15 +49,25 @@ // use that count to set order_use below + // 20140925 : the same logic below is used to constrain the dimensions of the image map + // figures to assess the fitting process: // x vs dx, x vs dy, y vs dx, y vs dy : // residual vector field // dx vs mag, dy vs mag + if (relastroGetVisual()) { + fprintf (stderr, "fitting chip %s, %d matched stars\n", image->name, Nmatch); + } relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch); + + // for polynomial fits, I used (5,30,60) = (5,6,7) stars per term // set the maximum order for the polynomial (based on number of stars kept above) int order_use = 0; - if (Nkeep > 5) order_use = 1; // 4 stars per polynomial term (per dimension) - if (Nkeep > 30) order_use = 2; // 5 stars per polynomial term (per dimension) - if (Nkeep > 60) order_use = 3; // 6 stars per polynomial term (per dimension) + if (Nkeep > 5) order_use = 1; // 5 stars per cell + if (Nkeep > 24) order_use = 2; // 6 stars per cell + if (Nkeep > 63) order_use = 3; // 7 stars per cell + if (Nkeep > 128) order_use = 4; // 8 stars per cell + if (Nkeep > 225) order_use = 5; // 9 stars per cell + if (Nkeep > 360) order_use = 6; // 10 stars per cell if (order_use < 1) { if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nkeep); @@ -65,5 +75,14 @@ return FALSE; } - image[0].coords.Npolyterms = order_use; + if (CHIPMAP) { + order_use = MIN(MIN(order_use, CHIPMAP), 6); // can only go up to 6th order map (can be user limited) + } else { + order_use = MIN(order_use, 3); // can only go up to 3rd order for polynomials + } + + if (VERBOSE2) fprintf (stderr, "using %d for %s\n", order_use, image[0].name); + + // when fitting the map, first fit a linear model (below? change Npolyterms to -1) + image[0].coords.Npolyterms = CHIPMAP ? 1 : order_use; if (fit) fit_free (fit); @@ -91,8 +110,36 @@ } + // apply fit to get the fitted X,Y coordinates. we need these to fit the residual map below for (i = 0; i < Nmatch; i++) { + // we have not yet fitted the map, so Npolyterms needs to be 1 here: + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[0].coords); + } + + if (CHIPMAP) { + if (image[0].coords.offsetMap == NULL) { + // allocate a new table and assign to this image + // need to lock threads here because 'table' is a static with + // modifiable elements + lockUpdateChips (); + AstromOffsetTable *table = get_astrom_table (); + AstromOffsetTableNewMap(table, order_use, order_use, image); + unlockUpdateChips (); + } + // XXX Need to check / update the order of the map here. if we already have a map, it might have been + // allocated with too small arrays + AstromOffsetMapSetOrder (image[0].coords.offsetMap, order_use, order_use, image); + image[0].coords.offsetMap->keep = TRUE; // if we are trying or re-trying this map, we should keep it unless it fails + fit_map (image[0].coords.offsetMap, raw, ref, Nmatch); + image[0].coords.Npolyterms = -1; + } + + if (!strcmp(image[0].name, "o5745g0403o.356777.cm.982487.smf[XY17.hdr]")) { + fprintf (stderr, "test image\n"); + } + + for (i = 0; i < Nmatch; i++) { + // if we have fitted the map above, Npolyterms needs to be -1 here: XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords); } - } @@ -223,19 +270,19 @@ L,M -> X,Y (polynomial transformation : WRP) -Some details about this function: - -- the initial value of the clipping radius is set based on the distribution of the dR - values for the bright sources. - -- NITER clipping passes are performed - -- the per-star astrometric errors are included in the fit. The photcode table controls - whether only the reported positional errors are used, or if the magnitudes errors are - scaled to determine the error, and if there is a systematic floor for all sources. - -- input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog - -- outlying detections are clipped in the iterative passes, but the clipped detections are - not recorded + Some details about this function: + + - the initial value of the clipping radius is set based on the distribution of the dR + values for the bright sources. + + - NITER clipping passes are performed + + - the per-star astrometric errors are included in the fit. The photcode table controls + whether only the reported positional errors are used, or if the magnitudes errors are + scaled to determine the error, and if there is a systematic floor for all sources. + + - input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog + + - outlying detections are clipped in the iterative passes, but the clipped detections are + not recorded */ @@ -243,13 +290,72 @@ /* example using fit_apply() : - f = fopen ("test3.dat", "w"); - - // apply new coords to raw (X,Y -> L,M) - for (i = 0; i < Nmatch; i++) { - fprintf (f, "%f %f %f %f ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M); - fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2); - fprintf (f, "%f %f\n", L1, M1); + f = fopen ("test3.dat", "w"); + + // apply new coords to raw (X,Y -> L,M) + for (i = 0; i < Nmatch; i++) { + fprintf (f, "%f %f %f %f ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M); + fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2); + fprintf (f, "%f %f\n", L1, M1); + } + fclose (f); +*/ + +int fit_map (AstromOffsetMap *map, StarData *raw, StarData *ref, int Npts) { + + // we are actually fitting the residual after the linear fit has been taken off + + // fit the linear terms as above + // calculate dX (raw.X - ref.X) and dY + + int i, N; + + float *x, *y, *dX, *dY; + ALLOCATE (x, float, Npts); + ALLOCATE (y, float, Npts); + ALLOCATE (dX, float, Npts); + ALLOCATE (dY, float, Npts); + + N = 0; + for (i = 0; i < Npts; i++) { + if (raw[i].mask) continue; + x[N] = raw[i].X; + y[N] = raw[i].Y; + dX[N] = ref[i].X - raw[i].X; + dY[N] = ref[i].Y - raw[i].Y; + N++; + } + + // in coordsops.c:XY_to_LM, the map is defined to carry dX,dY so that: + // (L,M) = f(X',Y') : (X',Y') = (X,Y) + (dX,dY) + + AstromOffsetMapFit (map, x, y, dX, N, TRUE); + AstromOffsetMapFit (map, x, y, dY, N, FALSE); + + AstromOffsetMapRepair (map, TRUE); + AstromOffsetMapRepair (map, FALSE); + + free (x); + free (y); + free (dX); + free (dY); + + return TRUE; +} + +int dump_stardata_pts (StarData *raw, int Npts, char *filename) { + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "failed to open file %s\n", filename); + return FALSE; + } + + int i; + for (i = 0; i < Npts; i++) { + fprintf (f, "%4d %10.6f %10.6f %10.6f %10.6f %9.2f %9.2f %9.2f %9.2f %9.4f %3d\n", + i, raw[i].R, raw[i].D, raw[i].P, raw[i].Q, raw[i].L, raw[i].M, raw[i].X, raw[i].Y, raw[i].dPos, raw[i].mask); } fclose (f); -*/ - + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/FrameCorrection.c =================================================================== --- trunk/Ohana/src/relastro/src/FrameCorrection.c (revision 37807) +++ trunk/Ohana/src/relastro/src/FrameCorrection.c (revision 37807) @@ -0,0 +1,617 @@ +# include "relastro.h" + +/* The sequence of calibration is something like this: + + * apply the current image.coords to measure.X,Y to get measure.R,D + * use the measure.R,D values to find the mean position average.R,D + * fit the measure.X,Y positions to average.R,D to get image.coords + + At some point in this sequence, we want to compare our current average.R,D (for specific + objects) to the ICRF values for those objects to get a correction set in R,D. + + Apply that correction to get a new set of average.R,D values. The impact of the + adjustment will then be to modify the image.coords to compensate for the offset. + +*/ + +int FrameCorrectionParallelMaster (RegionHostTable *regionHosts) { + + if (!USE_ICRF_CORRECT) return TRUE; + + int i; + for (i = 0; i < NLOOP; i++) { + if ((i > 1) || !USE_GALAXY_MODEL) { + // if GALAXY_MODEL is selected, we want to delay the frame correction until we have + // applied the galaxy model a couple of times. + + ICRFobj *icrfobj = slurp_icrf_obj (regionHosts, i); + FrameCorrectionSet *set = FrameCorrectionMeasure (icrfobj); + + char *filename = make_filename (CATDIR, "master", 0, "frame.corr.fits"); + FrameCorrectionSetSave (filename, set); + free (filename); + + char *syncfile = make_filename (CATDIR, "master", 0, "frame.corr.sync"); + update_sync_file (syncfile, i); + free (syncfile); + } + } + + int status = HarvestRegionHosts (regionHosts); + return status; +} + +int FrameCorrectionParallelSlave (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + if (!USE_ICRF_CORRECT) return TRUE; + + share_icrf_obj (catalog, Ncatalog, regionHosts, nloop); + + char *syncfile = make_filename (CATDIR, "master", 0, "frame.corr.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + char *filename = make_filename (CATDIR, "master", 0, "frame.corr.fits"); + FrameCorrectionSet *set = FrameCorrectionSetLoad (filename); + free (filename); + + // Now apply the correction to all of the average.R,D values + FrameCorrectionApply (catalog, Ncatalog, set->frame, set->coords); + FrameCorrectionSetFree(set); + + return TRUE; +} + +int FrameCorrectionSerial (Catalog *catalog, int Ncatalog) { + + if (!USE_ICRF_CORRECT) return TRUE; + + ICRFobj *icrfobj = get_ICRF_data (catalog, Ncatalog); + if (!icrfobj->Nicrfobj) { + fprintf (stderr, "no matched ICRF quasars for reference correction, skipping\n"); + return FALSE; + } + + FrameCorrectionSet *set = FrameCorrectionMeasure (icrfobj); + ICRFobjFree (icrfobj); + + // write out an image to represent the correction + static int version = 0; + char filename[1024]; + snprintf (filename, 1024, "%s/frame.%03d.corr.fits", CATDIR, version); + FrameCorrectionSetSave (filename, set); + version ++; + + // Now apply the correction to all of the average.R,D values + FrameCorrectionApply (catalog, Ncatalog, set->frame, set->coords); + FrameCorrectionSetFree(set); + + return TRUE; +} + +FrameCorrectionSet *FrameCorrectionMeasure (ICRFobj *icrfobj) { + + // we have only a few thousand ICRF reference points. I need to have a lookup table to + // go from ICRF object to the catalog entry. + + // we can have 2 kinds of corrections: + // 'frame' is a map in spherical coords of the corrections based on spherical harmonics + // 'map' is a correction based on a local projection to a linear coordinate system. I + // am going to use this for the pole, but also for SAS for testing. + // We need to be sure only one of the two is applied. For the pole, the boundary is a + // fixed line of DEC; for the SAS, the boundary is the projection 'image' + + FrameCorrectionSet *set = FrameCorrectionSetInit(); + if (!USE_ICRF_LOCAL && !USE_ICRF_SHFIT) { + fprintf (stderr, "ERROR: no ICRF correction mode selected\n"); + exit (0); + } + + if (USE_ICRF_LOCAL) { + // for the local frame correction, we are going to convert R,D into + // X,Y linear coordinates in an 'image' centered on the field center. + // we can then use an AstromOffsetMap structure to carry the correction values + + // what parameters define a local frame correction? + // * Ro, Do -- map/projection center + // * scale -- degrees / pixel + // * Nx, Ny -- size of map in pixels + // * Rmin,Rmax,Dmin,Dmax + + // example parameters for SAS: 320 - 340, -10 - +10 + double Ro = 330.0; + double Do = 0.0; + double scale = 2.0; // degrees per patch + int Nx = 10; + int Ny = 10; + if (USE_ICRF_POLE) { + Ro = 0.0; + Do = 90.0; + scale = 5.0; // degrees per patch + Nx = 4; + Ny = 4; + } + + + Coords *coords = NULL; + ALLOCATE (coords, Coords, 1); + + InitCoords (coords, "DEC--SIN"); + coords->cdelt1 = coords->cdelt2 = scale; + coords->crval1 = Ro; // SAS center + coords->crval2 = Do; + coords->crpix1 = 0.5*Nx; // middle of projection is middle of map + coords->crpix2 = 0.5*Ny; // middle of projection is middle of map + + // map is a component of coords + coords->offsetMap = AstromOffsetMapInit (Nx, Ny); + coords->offsetMap->dX = 1.0; // scale from projection (in arcsec) to correction patches + coords->offsetMap->dY = 1.0; + + FrameCorrectionFitLocal (icrfobj, coords); + set->coords = coords; + } + + if (USE_ICRF_SHFIT) { + // Lmax in recipe + int Lmax = 5; + + // dR and dD will carry the fit coefficients in RA and DEC + SHterms *dRc = SHtermsInit (Lmax); + SHterms *dDc = SHtermsInit (Lmax); + + INITTIME; + + FrameCorrectionFitSH (icrfobj, dRc, dDc); + MARKTIME ("done with FitSH: %f sec\n", dtime); + + double pltscale = 0.5; // degrees per pixel + FrameCorrectionType *frame = FrameCorrectionInit (pltscale); + + FrameCorrectionFromSH (frame, dRc, dDc); + MARKTIME ("done with FromSH: %f sec\n", dtime); + + SHtermsFree (dRc); + SHtermsFree (dDc); + + set->frame = frame; + } + + return set; +} + +// I need to generate a collection of values dR,dD(R,D), where dR = average.R - ICRF.R, etc. +// I will then fit the set of dR,dD values to a set of spherical harmonics (up to Lmax). +// I will then generate correction images dR,dD(R,D) +// I will then correct all values average.R,D by the interpolated values from the images +// average.R' = average.R - dR(R,d), average.D' = average.D - dD(R,d) + +int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc) { + + int i; + + myAssert (dRc->lmax == dDc->lmax, "dR and dD must match\n"); + + fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, icrfobj->Nicrfobj); + + int *mask = NULL; + ALLOCATE (mask, int, icrfobj->Nicrfobj); + memset (mask, 0, icrfobj->Nicrfobj*sizeof(int)); + + double *dPosObs, *dPosSrt = NULL; + ALLOCATE (dPosObs, double, icrfobj->Nicrfobj); + ALLOCATE (dPosSrt, double, icrfobj->Nicrfobj); + + int Niter; + for (Niter = 0; Niter < 3; Niter ++) { + + SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dRoff, mask, icrfobj->Nicrfobj, dRc); + SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dDoff, mask, icrfobj->Nicrfobj, dDc); + + // the bit below is for testing + + // allocate an SHterms structure to hold the Ylm values + SHterms *Ylm = SHtermsInit (dRc->lmax); + + static int version = 0; + char filename[1024]; + snprintf (filename, 1024, "%s/sh.%02d.dat", CATDIR, version); + FILE *f = fopen (filename, "w"); + myAssert (f, "oops"); + version ++; + + int Nkeep = 0; + for (i = 0; i < icrfobj->Nicrfobj; i++) { + + SHtermsForRD (Ylm, icrfobj->Rave[i], icrfobj->Dave[i]); + + double dRfit = 0.0; + double dDfit = 0.0; + + int j; + for (j = 0; j < Ylm->Nterms; j++) { + dRfit += dRc->Fr[j]*Ylm->Fr[j] + dRc->Fi[j]*Ylm->Fi[j]; + dDfit += dDc->Fr[j]*Ylm->Fr[j] + dDc->Fi[j]*Ylm->Fi[j]; + + } + fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f : %d\n", + icrfobj->Rave[i], icrfobj->Dave[i], icrfobj->dRoff[i], icrfobj->dDoff[i], dRfit, dDfit, mask[i]); + dPosObs[i] = hypot(icrfobj->dRoff[i] - dRfit, icrfobj->dDoff[i] - dDfit); + if (mask[i]) continue; + if (isnan(dPosObs[i])) continue; + dPosSrt[Nkeep] = dPosObs[i]; + Nkeep ++; + } + fclose (f); + + dsort (dPosSrt, Nkeep); + double dPmax = 1.5*dPosSrt[(int)(0.95*Nkeep)]; + + Nkeep = 0; + for (i = 0; i < icrfobj->Nicrfobj; i++) { + if (dPosObs[i] > dPmax) mask[i] = 1; + if (isnan(dPosObs[i])) mask[i] = 1; + if (!mask[i]) Nkeep++; + } + fprintf (stderr, "keeping %d of %d points\n", Nkeep, icrfobj->Nicrfobj); + } + + free (mask); + free (dPosObs); + free (dPosSrt); + + return TRUE; +} + +int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit) { + + int i, j, k; + + // allocate an SHterms structure to hold the Ylm values + SHterms *Ylm = SHtermsInit (fit->lmax); + + // we only fit the linearly independent terms: Re(m >= 0), Im(m > 0) + int NtermRE = 0; + int NtermIM = 0; + for (i = 0; i < Ylm->Nterms; i++) { + if (Ylm->m[i] >= 0) NtermRE ++; + if (Ylm->m[i] > 0) NtermIM ++; + } + int *Nre = NULL; + int *Nim = NULL; + ALLOCATE (Nre, int, NtermRE); + ALLOCATE (Nim, int, NtermIM); + + NtermRE = NtermIM = 0; + for (i = 0; i < Ylm->Nterms; i++) { + if (Ylm->m[i] >= 0) { + Nre[NtermRE] = i; + NtermRE ++; + } + if (Ylm->m[i] > 0) { + Nim[NtermIM] = i; + NtermIM ++; + } + } + + double **Are, **bre, **Aim, **bim; + ALLOCATE (Are, double *, NtermRE); + ALLOCATE (bre, double *, NtermRE); + ALLOCATE (Aim, double *, NtermIM); + ALLOCATE (bim, double *, NtermIM); + for (i = 0; i < NtermRE; i++) { + ALLOCATE_ZERO (Are[i], double, NtermRE); + ALLOCATE_ZERO (bre[i], double, 1); + } + for (i = 0; i < NtermIM; i++) { + ALLOCATE_ZERO (Aim[i], double, NtermIM); + ALLOCATE_ZERO (bim[i], double, 1); + } + + // measure the dot product \sum(F_i * Ylm_i) and the cross terms (\sum(Y_lm * Y_jk)) + int Nfit = 0; + for (i = 0; i < Npts; i++) { + if (mask[i]) continue; + Nfit ++; + + // set the values of Ylm(R[i],D[i]) + SHtermsForRD (Ylm, R[i], D[i]); + + double Fv = value[i]; + + for (j = 0; j < NtermRE; j++) { + int jre = Nre[j]; + bre[j][0] += Fv * Ylm->Fr[jre]; + } + for (j = 0; j < NtermIM; j++) { + int jim = Nim[j]; + bim[j][0] += Fv * Ylm->Fi[jim]; + } + + for (j = 0; j < NtermRE; j++) { + int jre = Nre[j]; + for (k = j; k < NtermRE; k++) { + int kre = Nre[k]; + Are[j][k] += Ylm->Fr[jre] * Ylm->Fr[kre]; + } + } + + for (j = 0; j < NtermIM; j++) { + int jim = Nim[j]; + for (k = j; k < NtermIM; k++) { + int kim = Nim[k]; + Aim[j][k] += Ylm->Fi[jim] * Ylm->Fi[kim]; + } + } + } + + for (j = 1; j < NtermRE; j++) { + for (k = 0; k < j; k++) { + Are[j][k] = Are[k][j]; + } + } + for (j = 1; j < NtermIM; j++) { + for (k = 0; k < j; k++) { + Aim[j][k] = Aim[k][j]; + } + } + + if (!dgaussjordan (Are, bre, NtermRE, 1)) { + gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n"); + return FALSE; + } + if (!dgaussjordan (Aim, bim, NtermIM, 1)) { + gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n"); + return FALSE; + } + + for (j = 0; j < fit->Nterms; j++) { + fit->Fr[j] = 0.0; + fit->Fi[i] = 0.0; + } + for (j = 0; j < NtermRE; j++) { + int jre = Nre[j]; + // fit->Fr[jre] = bre[j][0] * fit->Nterms / (float) Nfit; // XXX EAM : why is this factor needed? + fit->Fr[jre] = bre[j][0]; + } + for (j = 0; j < NtermIM; j++) { + int jim = Nim[j]; + // fit->Fi[jim] = bim[j][0] * fit->Nterms / (float) Nfit; // XXX EAM : why is this factor needed? + fit->Fi[jim] = bim[j][0]; + } + + SHtermsFree (Ylm); + + for (i = 0; i < NtermRE; i++) { + free (Are[i]); + free (bre[i]); + } + for (i = 0; i < NtermIM; i++) { + free (Aim[i]); + free (bim[i]); + } + free (Are); + free (bre); + free (Aim); + free (bim); + + return TRUE; +} + +int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords) { + + int i; + double Xave, Yave, Xmeas, Ymeas; + float *X, *Y, *dX, *dY; + float *dXf, *dYf, *dPos; + double *Rave, *Dave, *Rmeas, *Dmeas; + + int Nicrf = icrfobj->Nicrfobj; + + AstromOffsetMap *map = coords->offsetMap; + + int Npts = 0; + ALLOCATE (X, float, Nicrf); + ALLOCATE (Y, float, Nicrf); + ALLOCATE (dX, float, Nicrf); + ALLOCATE (dY, float, Nicrf); + ALLOCATE (dXf, float, Nicrf); + ALLOCATE (dYf, float, Nicrf); + + ALLOCATE (dPos, float, Nicrf); + + // arrays to store the points actually in the local region + ALLOCATE (Rave, double, Nicrf); + ALLOCATE (Dave, double, Nicrf); + ALLOCATE (Rmeas, double, Nicrf); + ALLOCATE (Dmeas, double, Nicrf); + + int *mask = NULL; + ALLOCATE (mask, int, Nicrf); + memset (mask, 0, Nicrf*sizeof(int)); + + // select the ICRF QSOS and save the necessary data + for (i = 0; i < Nicrf; i++) { + + // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax) + int status = RD_to_XY (&Xave, &Yave, icrfobj->Rave[i], icrfobj->Dave[i], coords); + if (!status) continue; + if (Xave < 0.0) continue; + if (Xave > map->Nx) continue; + if (Yave < 0.0) continue; + if (Yave > map->Ny) continue; + + double Rm = icrfobj->Rave[i] - icrfobj->dRoff[i] / 3600.0 / cos(icrfobj->Rave[i]*RAD_DEG); + double Dm = icrfobj->Dave[i] - icrfobj->dDoff[i] / 3600.0; + RD_to_XY (&Xmeas, &Ymeas, Rm, Dm, coords); + + // record these in arcsec or degree? + // correct for cos(D) or not? + X[Npts] = Xave; + Y[Npts] = Yave; + dX[Npts] = Xave - Xmeas; + dY[Npts] = Yave - Ymeas; + Rave[Npts] = icrfobj->Rave[i]; + Dave[Npts] = icrfobj->Dave[i]; + Rmeas[Npts] = Rm; + Dmeas[Npts] = Dm; + Npts ++; + } + + float *Xfit, *Yfit, *dXfit, *dYfit; + ALLOCATE (Xfit, float, Npts); + ALLOCATE (Yfit, float, Npts); + ALLOCATE (dXfit, float, Npts); + ALLOCATE (dYfit, float, Npts); + + int Niter; + for (Niter = 0; Niter < 3; Niter ++) { + int Nkeep = 0; + for (i = 0; i < Npts; i++) { + if (mask[i]) continue; + + Xfit[Nkeep] = X[i]; + Yfit[Nkeep] = Y[i]; + dXfit[Nkeep] = dX[i]; + dYfit[Nkeep] = dY[i]; + Nkeep ++; + } + + AstromOffsetMapFit (map, Xfit, Yfit, dXfit, Nkeep, TRUE); + AstromOffsetMapFit (map, Xfit, Yfit, dYfit, Nkeep, FALSE); + + for (i = 0; i < Npts; i++) { + dXf[i] = AstromOffsetMapValue (map, X[i], Y[i], TRUE); + dYf[i] = AstromOffsetMapValue (map, X[i], Y[i], FALSE); + } + + // generate a histogram of the dPos = sqrt((dXf - dX)^2 + (dYf - dY)^2) values + Nkeep = 0; + for (i = 0; i < Npts; i++) { + if (mask[i]) continue; + float dP = hypot((dXf[i] - dX[i]), (dYf[i] - dY[i])); + if (isnan(dP)) continue; // not all points are in the map + dPos[Nkeep] = dP; + Nkeep ++; + } + fsort (dPos, Nkeep); + + float dPmax = 1.5*dPos[(int)(0.95*Nkeep)]; + + // mask points based on dPmax + Nkeep = 0; + for (i = 0; i < Npts; i++) { + float dP = hypot((dXf[i] - dX[i]), (dYf[i] - dY[i])); + if (dP > dPmax) mask[i] = 1; + if (isnan(dP)) mask[i] = 1; + if (!mask[i]) Nkeep++; + } + fprintf (stderr, "keeping %d of %d points\n", Nkeep, Npts); + } + free (Xfit); + free (Yfit); + free (dXfit); + free (dYfit); + + // *** save the icrf points and fit residuals + { + static int version = 0; + char filename[1024]; + snprintf (filename, 1024, "%s/map.%02d.dat", CATDIR, version); + FILE *f = fopen (filename, "w"); + version ++; + + for (i = 0; i < Npts; i++) { + fprintf (f, "%3d %12.8f %12.8f %12.8f %12.8f %f %f : %f %f : %f %f : %d\n", i, Rave[i], Dave[i], Rmeas[i], Dmeas[i], X[i], Y[i], dX[i], dY[i], dXf[i], dYf[i], mask[i]); + } + fclose (f); + } + + free (X); + free (Y); + free (dX); + free (dY); + free (dXf); + free (dYf); + free (dPos); + free (Rave); + free (Dave); + free (Rmeas); + free (Dmeas); + free (mask); + + return TRUE; +} + +int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, Coords *coords) { + + double Xave, Yave; + + AstromOffsetMap *map = coords ? coords->offsetMap : NULL; + + int i, j; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + Average *average = &catalog[i].average[j]; + + double R = average->R; + double D = average->D; + + int inLocalRegion = FALSE; + if (USE_ICRF_LOCAL) { + myAssert (coords, "no local frame coords defined"); + + // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax) + inLocalRegion = RD_to_XY (&Xave, &Yave, R, D, coords); + if (USE_ICRF_POLE) { + inLocalRegion = (D > 85.0); + } else { + inLocalRegion = (inLocalRegion && (Xave > 0.0)); + inLocalRegion = (inLocalRegion && (Xave < map->Nx)); + inLocalRegion = (inLocalRegion && (Yave > 0.0)); + inLocalRegion = (inLocalRegion && (Yave < map->Ny)); + } + } + + if (inLocalRegion) { + myAssert (map, "no local frame map defined"); + double dX = AstromOffsetMapValue (map, Xave, Yave, TRUE); + double dY = AstromOffsetMapValue (map, Xave, Yave, FALSE); + if (!isnan(dX) && !isnan(dY)) { + double Xmeas = Xave - dX; + double Ymeas = Yave - dY; + XY_to_RD (&average->R, &average->D, Xmeas, Ymeas, coords); + } + continue; + } + + if (USE_ICRF_SHFIT) { + myAssert (frame, "no frame correction defined"); + + int iD = (D + 89.0) / frame->scale; + if (iD < 0) continue; + if (iD >= frame->Ndec) continue; + + int iR = R / frame->dR[iD]; + if (iR < 0) continue; + if (iR >= frame->Nra[iD]) continue; + + double dR = frame->Roff[iD][iR]; + double dD = frame->Doff[iD][iR]; + + // XXX tighten this up?? + // do not apply if the fabs(offset) is > 8 arcsec) + if ((fabs(dR) > 8.0) || (fabs(dD) > 8.0)) { + fprintf (stderr, "skip: %10.6f %10.6f : %7.3f %7.3f\n", R, D, dR, dD); + continue; + } + + average->R -= dR / 3600.0 / cos(average->R*RAD_DEG); + average->D -= dD / 3600.0; + } + } + } + return TRUE; +} + Index: trunk/Ohana/src/relastro/src/FrameCorrectionIO.c =================================================================== --- trunk/Ohana/src/relastro/src/FrameCorrectionIO.c (revision 37807) +++ trunk/Ohana/src/relastro/src/FrameCorrectionIO.c (revision 37807) @@ -0,0 +1,251 @@ +# include "relastro.h" + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +// the frame correction consists of N image correction maps +// there are local maps and AITOFF maps +FrameCorrectionSet *FrameCorrectionSetLoad(char *filename) { + + Header header; + Matrix matrix; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + // FrameCorrectionSet consists: + // PHU header listing the contents: + // SHMAP_dR : EXTNAME or NONE + // SHMAP_dD : EXTNAME or NONE + // LOCAL_dR : EXTNAME or NONE + // LOCAL_dD : EXTNAME or NONE + + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + + char SHMAP_DR[80], SHMAP_DD[80], LOCAL_DR[80], LOCAL_DD[80]; + if (!gfits_scan (&header, "SHMAP_DR", "%s", 1, SHMAP_DR)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field SHMAP_DR\n", filename); + exit (2); + } + if (!gfits_scan (&header, "SHMAP_DD", "%s", 1, SHMAP_DD)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field SHMAP_DD\n", filename); + exit (2); + } + if (!gfits_scan (&header, "LOCAL_DR", "%s", 1, LOCAL_DR)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field LOCAL_DR\n", filename); + exit (2); + } + if (!gfits_scan (&header, "LOCAL_DD", "%s", 1, LOCAL_DD)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field LOCAL_DD\n", filename); + exit (2); + } + gfits_free_header (&header); + + FrameCorrectionSet *set = FrameCorrectionSetInit (); + + if (strcmp(SHMAP_DR, "NONE")) { + if (!strcmp (SHMAP_DD, "NONE")) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DR but not SHMAP_DD\n", filename); + exit (2); + } + if (!gfits_find_Xheader(f, &header, SHMAP_DR)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DR extension %s\n", filename, SHMAP_DR); + exit (2); + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read SHMAP_DR extension image\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + // convert image to frame + set->frame = FrameCorrectionImageToSH (&header, &matrix, NULL, TRUE); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + if (!gfits_find_Xheader(f, &header, SHMAP_DD)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DD extension %s\n", filename, SHMAP_DD); + exit (2); + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read SHMAP_DD extension image\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + set->frame = FrameCorrectionImageToSH (&header, &matrix, set->frame, TRUE); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } else { + if (strcmp (SHMAP_DD, "NONE")) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DD but not SHMAP_DR\n", filename); + exit (2); + } + } + + if (strcmp(LOCAL_DR, "NONE")) { + if (!strcmp (LOCAL_DD, "NONE")) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined LOCAL_DR but not LOCAL_DD\n", filename); + exit (2); + } + if (!gfits_find_Xheader(f, &header, LOCAL_DR)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing LOCAL_DR extension %s\n", filename, LOCAL_DR); + exit (2); + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read LOCAL_DR extension image\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + // convert image to frame + set->coords = FrameCorrectionImageToMap (&header, &matrix, NULL, TRUE); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + if (!gfits_find_Xheader(f, &header, LOCAL_DD)) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s missing LOCAL_DD extension %s\n", filename, LOCAL_DD); + exit (2); + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read LOCAL_DD extension image\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + set->coords = FrameCorrectionImageToMap (&header, &matrix, set->coords, FALSE); + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } else { + if (strcmp (LOCAL_DD, "NONE")) { + fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined LOCAL_DD but not LOCAL_DR\n", filename); + exit (2); + } + } + fclose (f); + return set; +} + +int FrameCorrectionSetSave(char *filename, FrameCorrectionSet *set) { + + int status; + Header header; + Matrix matrix; + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + // add entries to the PHU + if (set->frame) { + gfits_modify (&header, "SHMAP_DR", "%s", 1, "SHMAP_DR"); + gfits_modify (&header, "SHMAP_DD", "%s", 1, "SHMAP_DD"); + } else { + gfits_modify (&header, "SHMAP_DR", "%s", 1, "NONE"); + gfits_modify (&header, "SHMAP_DD", "%s", 1, "NONE"); + } + if (set->coords) { + gfits_modify (&header, "LOCAL_DR", "%s", 1, "LOCAL_DR"); + gfits_modify (&header, "LOCAL_DD", "%s", 1, "LOCAL_DD"); + } else { + gfits_modify (&header, "LOCAL_DR", "%s", 1, "NONE"); + gfits_modify (&header, "LOCAL_DD", "%s", 1, "NONE"); + } + + // write the PHU header (& empty matrix) to disk + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for frame correction %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for correction %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + if (set->frame) { + // create and save the RA correction + FrameCorrectionSHtoImage (&header, &matrix, set->frame, TRUE); + gfits_modify (&header, "EXTNAME", "%s", 1, "SHMAP_DR"); + + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for SHMAP_DR %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for SHMAP_DR %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + // create and save the DEC correction + FrameCorrectionSHtoImage (&header, &matrix, set->frame, FALSE); + gfits_modify (&header, "EXTNAME", "%s", 1, "SHMAP_DD"); + + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for SHMAP_DD %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for SHMAP_DD %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + + if (set->coords) { + // create and save the RA correction + FrameCorrectionMapToImage (&header, &matrix, set->coords, TRUE); + gfits_modify (&header, "EXTNAME", "%s", 1, "LOCAL_DR"); + + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + + // create and save the DEC correction + FrameCorrectionMapToImage (&header, &matrix, set->coords, FALSE); + gfits_modify (&header, "EXTNAME", "%s", 1, "LOCAL_DD"); + + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + } + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file meanpos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file meanpos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing meanpos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c =================================================================== --- trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 37807) +++ trunk/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 37807) @@ -0,0 +1,277 @@ +# include "relastro.h" + +// functions related to structure management + +FrameCorrectionSet *FrameCorrectionSetInit () { + + FrameCorrectionSet *set = NULL; + ALLOCATE (set, FrameCorrectionSet, 1); + + set->frame = NULL; + set->coords = NULL; + return set; +} + +void FrameCorrectionSetFree (FrameCorrectionSet *set) { + + if (!set) return; + + FrameCorrectionFree(set->frame); + if (set->coords) AstromOffsetMapFree (set->coords->offsetMap); + free (set->coords); + + free (set); + return; +} + +FrameCorrectionType *FrameCorrectionInit (double scale) { + + int i; + + FrameCorrectionType *frame = NULL; + ALLOCATE (frame, FrameCorrectionType, 1); + + frame->scale = scale; // degrees / pixel + frame->Ndec = 180.0 / scale; + if (frame->Ndec * scale < 180.0) frame->Ndec ++; + + // binD = (DEC + 89.0) / scale + ALLOCATE (frame->Nra, int, frame->Ndec); + ALLOCATE (frame->dR, double, frame->Ndec); + ALLOCATE (frame->Roff, double *, frame->Ndec); + ALLOCATE (frame->Doff, double *, frame->Ndec); + + // I am going to avoid corrections at |D| = 90 + + // It may be sleazy, but since I'm using SH, I'm going to just correct the region at |D| + // > 80 with a 2D grid of corrections to X,Y projections + + for (i = 0; i < frame->Ndec; i++) { + double D = i * scale - 89.0; + if (D > 90.0) { + frame->dR[i] = 0.0; + frame->Nra[i] = 0; + frame->Roff[i] = NULL; + frame->Doff[i] = NULL; + continue; + } + frame->dR[i] = scale / cos(RAD_DEG*D); + frame->Nra[i] = MAX (0, 360.0 / frame->dR[i]); + if (frame->Nra[i] * frame->dR[i] < 360.0) frame->Nra[i] ++; + + ALLOCATE (frame->Roff[i], double, frame->Nra[i]); + ALLOCATE (frame->Doff[i], double, frame->Nra[i]); + // ohana_memcheck_func (TRUE); + // fprintf (stderr, "alloc : %d : %lx %lx\n", i, (long unsigned int)(size_t *) frame->Roff[i], (long unsigned int) (size_t *) frame->Doff[i]); + } + return frame; +} + +void FrameCorrectionFree (FrameCorrectionType *frame) { + + if (!frame) return; + + int i; + + for (i = 0; i < frame->Ndec; i++) { + if (frame->Roff[i]) free (frame->Roff[i]); + if (frame->Doff[i]) free (frame->Doff[i]); + } + + if (frame->dR) free (frame->dR); + if (frame->Nra) free (frame->Nra); + if (frame->Roff) free (frame->Roff); + if (frame->Doff) free (frame->Doff); + + free (frame); + + return; +} + +int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD) { + + myAssert (dR->lmax == dD->lmax, "dR and dD must match\n"); + + // allocate an SHterms structure to hold the Ylm values + SHterms *SH = SHtermsInit (dR->lmax); + + int i, j, k; + // Now I need to apply the correction coefficients to SH images to generate an image + for (i = 0; i < frame->Ndec; i++) { + double D = i * frame->scale - 89.0; + for (j = 0; j < frame->Nra[i]; j++) { + double R = j * frame->dR[i]; + + SHtermsForRD (SH, R, D); + + frame->Roff[i][j] = 0.0; + frame->Doff[i][j] = 0.0; + + for (k = 0; k < SH->Nterms; k++) { + frame->Roff[i][j] += dR->Fr[k]*SH->Fr[k] + dR->Fi[k]*SH->Fi[k]; + frame->Doff[i][j] += dD->Fr[k]*SH->Fr[k] + dD->Fi[k]*SH->Fi[k]; + } + } + } + SHtermsFree (SH); + return TRUE; +} + +// write out the maps as images +FrameCorrectionType *FrameCorrectionImageToSH (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection) { + + int i, j; + + if (!frame) { + double pltscale; + if (!gfits_scan (header, "PLTSCALE", "%lf", 1, &pltscale)) { + fprintf (stderr, "failed to find PLTSCALE in header\n"); + exit (2); + } + frame = FrameCorrectionInit (pltscale); + } + + Coords coords; + GetCoords (&coords, header); + + int Nx = header->Naxis[0]; + // int Ny = header->Naxis[1]; + + float *buffer = (float *)matrix->buffer; + double **value = (raDirection) ? frame->Roff : frame->Doff; + + // Now I need to apply the correction coefficients to SH images to generate an image + for (i = 0; i < frame->Ndec; i++) { + double D = i * frame->scale - 89.0; + for (j = 0; j < frame->Nra[i]; j++) { + double R = j * frame->dR[i]; + + value[i][j] = 0.0; + + double X, Y; + if (!RD_to_XY (&X, &Y, R, D, &coords)) continue; + + int ix = X; + int iy = Y; + value[i][j] = buffer[ix + Nx*iy]; + } + } + return frame; +} + +int FrameCorrectionSHtoImage (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection) { + + Coords coords; + + gfits_init_header (header); + header->bitpix = -32; + header->Naxes = 2; + header->Naxis[0] = 370.0 / frame->scale; + header->Naxis[1] = 190.0 / frame->scale; + + int Nx = header->Naxis[0]; + int Ny = header->Naxis[1]; + + gfits_create_header (header); + gfits_create_matrix (header, matrix); + + gfits_modify (header, "PLTSCALE", "%lf", 1, frame->scale); + + InitCoords (&coords, "DEC--AIT"); + coords.cdelt1 = coords.cdelt2 = frame->scale; + coords.crval1 = 0.0; // SAS center + coords.crval2 = 0.0; + coords.crpix1 = 0.5*Nx; // middle of projection is middle of map + coords.crpix2 = 0.5*Ny; // middle of projection is middle of map + + PutCoords (&coords, header); + + float *buffer = (float *)matrix->buffer; + + int ix, iy; + for (ix = 0; ix < Nx; ix++) { + for (iy = 0; iy < Ny; iy++) { + double R, D; + int status = XY_to_RD (&R, &D, ix, iy, &coords); + if (!status) continue; + + R = ohana_normalize_angle (R); + + int iD = (D + 89.0) / frame->scale; + if (iD < 0) continue; + if (iD >= frame->Ndec) continue; + int iR = R / frame->dR[iD]; + if (iR < 0) continue; + if (iR >= frame->Nra[iD]) continue; + + + float value = NAN; + if (raDirection) { + value = frame->Roff[iD][iR]; + } else { + value = frame->Doff[iD][iR]; + } + buffer[ix + Nx*iy] = value; + } + } + return TRUE; +} + +// write out the AstromOffsetMap linked by the Coords structure +int FrameCorrectionMapToImage (Header *header, Matrix *matrix, Coords *coords, int raDirection) { + // write out the maps as images + + AstromOffsetMap *map = coords->offsetMap; + if (!map) { + fprintf (stderr, "ERROR: missing AstromOffsetMap link in coords\n"); + return FALSE; + } + + gfits_init_header (header); + header->bitpix = -32; + header->Naxes = 2; + header->Naxis[0] = map->Nx; + header->Naxis[1] = map->Ny; + + gfits_create_header (header); + gfits_create_matrix (header, matrix); + PutCoords (coords, header); + + float *buffer = (float *)matrix->buffer; + + float **value = (raDirection) ? map->dXv : map->dYv; + + int ix, iy; + for (ix = 0; ix < map->Nx; ix++) { + for (iy = 0; iy < map->Ny; iy++) { + buffer[ix + map->Nx*iy] = value[ix][iy]; + } + } + return TRUE; +} + +Coords *FrameCorrectionImageToMap (Header *header, Matrix *matrix, Coords *coords, int raDirection) { + // write out the maps as images + + if (!coords) { + ALLOCATE (coords, Coords, 1); + GetCoords (coords, header); + coords->offsetMap = AstromOffsetMapInit (header->Naxis[0], header->Naxis[1]); + coords->offsetMap->dX = 1.0; // scale from projection (in arcsec) to correction patches + coords->offsetMap->dY = 1.0; + } + AstromOffsetMap *map = coords->offsetMap; + + float *buffer = (float *) matrix->buffer; + + float **value = (raDirection) ? map->dXv : map->dYv; + + int ix, iy; + for (ix = 0; ix < map->Nx; ix++) { + for (iy = 0; iy < map->Ny; iy++) { + value[ix][iy] = buffer[ix + map->Nx*iy]; + } + } + return coords; +} + Index: trunk/Ohana/src/relastro/src/ICRF.c =================================================================== --- trunk/Ohana/src/relastro/src/ICRF.c (revision 37807) +++ trunk/Ohana/src/relastro/src/ICRF.c (revision 37807) @@ -0,0 +1,92 @@ +# include "relastro.h" + +int Nicrf = 0; +int NICRF = 0; +int *ICRFtoCatalog = NULL; +int *ICRFtoAverage = NULL; +int *ICRFtoMeasure = NULL; + +void ICRFinit () { + + Nicrf = 0; + NICRF = 500; + ALLOCATE (ICRFtoCatalog, int, NICRF); + ALLOCATE (ICRFtoAverage, int, NICRF); + ALLOCATE (ICRFtoMeasure, int, NICRF); +} + +int ICRFsave (int cat, int ave, int meas) { + + // fprintf (stderr, "ICRF: %d %d %d\n", cat, ave, meas); + + ICRFtoCatalog[Nicrf] = cat; + ICRFtoAverage[Nicrf] = ave; + ICRFtoMeasure[Nicrf] = meas; + + Nicrf ++; + if (Nicrf == NICRF) { + NICRF += 500; + + REALLOCATE (ICRFtoCatalog, int, NICRF); + REALLOCATE (ICRFtoAverage, int, NICRF); + REALLOCATE (ICRFtoMeasure, int, NICRF); + } + return TRUE; +} + +int ICRFmax () { + return Nicrf; +} + +int ICRFdata (int n, int *cat, int *ave, int *meas) { + + myAssert (n >= 0 && n < Nicrf, "out of range"); + + *cat = ICRFtoCatalog[n]; + *ave = ICRFtoAverage[n]; + *meas = ICRFtoMeasure[n]; + return TRUE; +} + +int select_catalog_ICRF (Catalog *catalog, int Ncatalog) { + + if (!USE_ICRF_CORRECT) return TRUE; + + int N = 0; + + char filename[1024]; + snprintf (filename, 1024, "%s/test.icrf.dat", CATDIR); + FILE *f = fopen (filename, "w"); + + int c, i, j; + for (c = 0; c < Ncatalog; c++) { + + for (i = 0; i < catalog[c].Naverage; i++) { + if (!(catalog[c].average[i].flags & ID_ICRF_QSO)) continue; + + // only save a single value + int savedICRF = FALSE; + for (j = 0; j < catalog[c].average[i].Nmeasure; j++) { + + int offset = catalog[c].average[i].measureOffset + j; + + if (!(catalog[c].measureT[offset].dbFlags & ID_MEAS_ICRF_QSO)) continue; + + fprintf (f, "%d %d %d : %f %f : %f %f\n", c, i, offset, + catalog[c].average[i].R, catalog[c].average[i].D, + catalog[c].measureT[offset].R, catalog[c].measureT[offset].D); + + if (savedICRF) continue; + + ICRFsave (c, i, offset); + N ++; + savedICRF = TRUE; + } + } + } + fclose (f); + + fprintf (stderr, "added %d ICRF QSO\n", N); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/ICRFobjIO.c =================================================================== --- trunk/Ohana/src/relastro/src/ICRFobjIO.c (revision 37807) +++ trunk/Ohana/src/relastro/src/ICRFobjIO.c (revision 37807) @@ -0,0 +1,146 @@ +# include "relastro.h" + +// column OUT refers to a variable which is already defined +# define GET_COLUMN(OUT,NAME,TYPE) \ + OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +ICRFobj *ICRFobjLoad(char *filename) { + + int Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + ICRFobj *icrfobj = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + char type[16]; // used in GET_COLUMN + + ALLOCATE (icrfobj, ICRFobj, 1); + + GET_COLUMN (icrfobj->Rave, "R_AVE", double); + GET_COLUMN (icrfobj->Dave, "D_AVE", double); + GET_COLUMN (icrfobj->dRoff, "R_OFF", double); + GET_COLUMN (icrfobj->dDoff, "D_OFF", double); + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + // free FITS table pieces... + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + icrfobj->Nicrfobj = Nrow; + return icrfobj; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int ICRFobjSave(char *filename, ICRFobj *icrfobj) { + + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "ICRFOBJ"); + + gfits_define_bintable_column (&theader, "D", "R_AVE", "mean position, ra", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "D_AVE", "mean position, dec", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "R_OFF", "position offset, ra", "arcsec", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "D_OFF", "position offset, dec", "arcsec", 1.0, 0.0); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "R_AVE", icrfobj->Rave , icrfobj->Nicrfobj); + gfits_set_bintable_column (&theader, &ftable, "D_AVE", icrfobj->Dave , icrfobj->Nicrfobj); + gfits_set_bintable_column (&theader, &ftable, "R_OFF", icrfobj->dRoff, icrfobj->Nicrfobj); + gfits_set_bintable_column (&theader, &ftable, "D_OFF", icrfobj->dDoff, icrfobj->Nicrfobj); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for icrfobj %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for icrfobj %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for icrfobj %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for icrfobj %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file icrfobj %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot fsync file icrfobj %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing icrfobj file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 37261) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 37807) @@ -28,5 +28,14 @@ static off_t *imageIDs; // list of all image IDs static off_t *imageIdx; // list of index for image IDs -# else + +// as an alternative, we generate imageSeq, which directly maps imageID -> seq +static off_t *imageSeq; // list of index for image IDs + +// MAX_ID requires 512M to store the image index +# define MAX_ID 0x8000000 +static off_t minImageID = MAX_ID; +static off_t maxImageID = 0; + +# else /* not using IMAGE_ID */ static unsigned int *start; static unsigned int *stop; @@ -84,6 +93,21 @@ imageIdx[i] = i; imageIDs[i] = image[i].imageID; + minImageID = MIN(minImageID, image[i].imageID); + maxImageID = MAX(maxImageID, image[i].imageID); + myAssert (image[i].imageID < MAX_ID, "image IDs too large for index memory"); } llsortpair (imageIDs, imageIdx, Nimage); + + ALLOCATE (imageSeq, off_t, maxImageID + 1); + for (i = 0; i < maxImageID + 1; i++) { + imageSeq[i] = -1; // not yet assigned + } + + for (i = 0; i < Nimage; i++) { + off_t N = image[i].imageID; + myAssert (imageSeq[N] == -1, "previously assigned"); + imageSeq[N] = i; + } + # else ALLOCATE (start, unsigned, Nimage); @@ -98,9 +122,16 @@ off_t getImageByID (off_t ID) { +# if USE_IMAGE_ID + + if (imageSeq) { + if (ID < minImageID) return (-1); + if (ID > maxImageID) return (-1); + off_t N = imageSeq[ID]; + return N; + } // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs // use bisection to find the specified image ID -# if USE_IMAGE_ID off_t Nlo, Nhi, N; @@ -207,18 +238,19 @@ if (VERBOSE2 || (MATCHCAT && (N_onImage[i] < 20))) { name = GetPhotcodeNamebyCode (image[i].photcode); - fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s)", i, image[i].imageID, image[i].name, N_onImage[i], image[i].nstar, + int showExample = (Nfew < 30); + if (showExample) fprintf (stderr, "image "OFF_T_FMT" (%d, %s) has "OFF_T_FMT" of %d measures (%s, %s) ", i, image[i].imageID, image[i].name, N_onImage[i], image[i].nstar, ohana_sec_to_date(image[i].tzero), name); if (N_onImage[i] < 20) { - fprintf (stderr, "*"); + if (showExample) fprintf (stderr, "*"); Nfew ++; } if (N_onImage[i] < 15) { - fprintf (stderr, "*** warning ***"); + if (showExample) fprintf (stderr, "*** warning ***"); Nbad ++; } - fprintf (stderr, "\n"); - } - } - fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad); + if (showExample) fprintf (stderr, "\n"); + } + } + if (MATCHCAT) fprintf (stderr, OFF_T_FMT" total images, %d with < 20 measurements, %d with < 15\n", Nimage, Nfew, Nbad); } @@ -245,4 +277,6 @@ } + if (!ID) return; // detection not associated with an image + idx = getImageByID (ID); if (idx == -1) { @@ -250,4 +284,5 @@ return; } + measure->myDet = TRUE; // I 'own' this detection (I am calibrating its image) // index for (catalog, measure) -> image @@ -334,15 +369,4 @@ # endif -/* -Coords *getCoords (off_t meas, int cat) { - - off_t i; - - i = MeasureToImage[cat][meas]; - if (i == -1) return (NULL); - return (&image[i].coords); -} -*/ - float getColorBlue (off_t meas, int cat) { @@ -433,7 +457,4 @@ Coords *moscoords, *imcoords; - // check if this image is bad and should be skipped - if (badCoords(im)) return; - // WRP images need to have an associated mosaic moscoords = NULL; @@ -491,7 +512,6 @@ int TESTPT = FALSE; - - TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->imageID == CAT_ID_SRC) && (measure->detID == OBJ_ID_SRC); - TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->imageID == CAT_ID_DST) && (measure->detID == OBJ_ID_DST); + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (measure->catID == CAT_ID_SRC) && (measure->objID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (measure->catID == CAT_ID_DST) && (measure->objID == OBJ_ID_DST); if (TESTPT) { fprintf (stderr, "got test det\n"); @@ -597,5 +617,4 @@ int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori; if (VERBOSE2 && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); - saveOffsets (dPos, nPos, im); return; @@ -604,86 +623,4 @@ void printNcatTotal () { fprintf (stderr, "NcatTotal: %d\n", NcatTotal); -} - -// return StarData values for detections in the specified image, converting coordinates from the -// chip positions: X,Y -> L,M -> P,Q -> R,D -void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im) { - - off_t i, m, c, n; - double X, Y, L, M, P, Q, R, D; - - Mosaic *mosaic; - Coords *moscoords, *imcoords, *oldcoords; - - fprintf (stderr, "fix resetImageRaw wrt MeasureTiny\n"); - abort(); - - // check if this image is bad and should be skipped - if (!badCoords(im)) { - fprintf (stderr, "ERROR: inconsistent result?"); - abort(); - } - - // replace the current coords with the old coords: - oldcoords = getCoords (im); - memcpy (&image[im].coords, oldcoords, sizeof(Coords)); - - // WRP images need to have an associated mosaic - moscoords = NULL; - if (!strcmp(&image[im].coords.ctype[4], "-WRP")) { - mosaic = getMosaicForImage (im); - if (mosaic == NULL) return; // if we cannot find the associated image, skip it - moscoords = &mosaic[0].coords; - } - imcoords = &image[im].coords; - - for (i = 0; i < N_onImage[im]; i++) { - m = ImageToMeasure[im][i]; - c = ImageToCatalog[im][i]; - - // XXX unclear if this should use Measure or MeasureTiny; it is only called by FixProblemImages, which has not been updated - Measure *measure = &catalog[c].measure[m]; - - X = measure[0].Xccd; - Y = measure[0].Yccd; - if (USE_FIXED_PIXCOORDS) { - if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) { - float dX = measure[0].Xfix - measure[0].Xccd; - float dY = measure[0].Yfix - measure[0].Yccd; - if (hypot(dX,dY) < 2.0) { - X = measure[0].Xfix; - Y = measure[0].Yfix; - } - } - } - n = measure[0].averef; - - if (moscoords == NULL) { - // this is a Simple image (not a mosaic) - // note that for a Simple image, L,M = P,Q - XY_to_LM (&L, &M, X, Y, imcoords); - LM_to_RD (&R, &D, L, M, imcoords); - } else { - XY_to_LM (&L, &M, X, Y, imcoords); - XY_to_LM (&P, &Q, L, M, moscoords); - LM_to_RD (&R, &D, P, Q, moscoords); - } - measure[0].R = R; - measure[0].D = D; - - float dRoff = dvoOffsetR(measure, &catalog[c].average[n]); - - if (dRoff > +180.0*3600.0) { - // average on high end of boundary, move star up - measure[0].R += 360.0; - dRoff -= 360.0*3600.0; - } - if (dRoff < -180.0*3600.0) { - measure[0].R -= 360.0; - dRoff += 360.0*3600.0; - } - } - - return; } @@ -736,4 +673,11 @@ n = measure[0].averef; + + int TESTPT = FALSE; + TESTPT |= CAT_ID_SRC && OBJ_ID_SRC && (catalog[c].average[n].catID == CAT_ID_SRC) && (catalog[c].average[n].objID == OBJ_ID_SRC); + TESTPT |= CAT_ID_DST && OBJ_ID_DST && (catalog[c].average[n].catID == CAT_ID_DST) && (catalog[c].average[n].objID == OBJ_ID_DST); + if (TESTPT) { + fprintf (stderr, "got test det\n"); + } // an object with only one detection provides no information about the image calibration @@ -808,4 +752,12 @@ ref[i].D = catalog[c].average[n].D; + // if we are correcting for the Galaxy Motion Model, we assume the mean R,D is at the J2000 epoch position + if (USE_GALAXY_MODEL) { + myAssert (!isnan(measure[0].RoffGAL), "oops"); + myAssert (!isnan(measure[0].DoffGAL), "oops"); + ref[i].R += measure[0].RoffGAL / 3600.0; + ref[i].D += measure[0].DoffGAL / 3600.0; + } + ref[i].Mag = measure[0].M; ref[i].dMag = measure[0].dM; @@ -842,4 +794,95 @@ *Nstars = N_onImage[im]; return (ref); +} + +// return StarData values for detections in the specified image, converting coordinates from the +// chip positions: X,Y -> L,M -> P,Q -> R,D. This function is used by the image fitting steps, for +// which the detections have already been filtered when they were loaded (bcatalog) +int setImageRaw (Catalog *catalog, int Ncatalog, off_t im, StarData *raw, off_t Nraw, CoordMode mode) { + + off_t i, m, c; + + Coords *moscoords; + + moscoords = NULL; + if (mode == MODE_MOSAIC) { + moscoords = image[im].coords.mosaic; + myAssert (moscoords, "coords.mosaic not defined for image %s (%d)", image[im].name, (int) im); + } + + myAssert (Nraw == N_onImage[im], "impossible!"); + + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + // XXX should I use the raw coords or just measure.X,Y -> R,D? + + switch (mode) { + case MODE_SIMPLE: + /* note that for a Simple image, L,M = P,Q */ + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); + raw[i].P = raw[i].L; + raw[i].Q = raw[i].M; + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords); + break; + case MODE_MOSAIC: + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); + XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords); + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords); + break; + default: + fprintf (stderr, "error: invalid mode in getImageRaw"); + abort (); + } + + MeasureTiny *measure = &catalog[c].measureT[m]; + measure->R = raw[i].R; + measure->D = raw[i].D; + if (catalog[c].measure) { + catalog[c].measure[m].R = raw[i].R; + catalog[c].measure[m].D = raw[i].D; + } + } + return TRUE; +} + +// return StarData values for detections in the specified image, converting coordinates from the +// chip positions: X,Y -> L,M -> P,Q -> R,D. This function is used by the image fitting steps, for +// which the detections have already been filtered when they were loaded (bcatalog) +int updateImageRaw (Catalog *catalog, int Ncatalog, off_t im) { + + off_t i, m, c; + + for (i = 0; i < N_onImage[im]; i++) { + m = ImageToMeasure[im][i]; + c = ImageToCatalog[im][i]; + + MeasureTiny *measure = &catalog[c].measureT[m]; + + /* apply the current image transformation or use the current value of R+dR, D+dD? */ + double X = measure[0].Xccd; + double Y = measure[0].Yccd; + if (USE_FIXED_PIXCOORDS) { + if (isfinite(measure[0].Xfix) && isfinite(measure[0].Yfix)) { + float dX = measure[0].Xfix - measure[0].Xccd; + float dY = measure[0].Yfix - measure[0].Yccd; + if (hypot(dX,dY) < 2.0) { + X = measure[0].Xfix; + Y = measure[0].Yfix; + } + } + } + + double R, D; + XY_to_RD (&R, &D, X, Y, &image[im].coords); + measure->R = R; + measure->D = D; + if (catalog[c].measure) { + catalog[c].measure[m].R = R; + catalog[c].measure[m].D = D; + } + } + return TRUE; } Index: trunk/Ohana/src/relastro/src/ImagePosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 37261) +++ trunk/Ohana/src/relastro/src/ImagePosIO.c (revision 37807) @@ -70,7 +70,9 @@ GET_COLUMN (polyterms, "POLYTERMS", float); // verify that we got 14 columns? GET_COLUMN (ctype, "CTYPE", char); // verify that we got 15 columns? - GET_COLUMN (Npolyterms, "NPOLYTERMS", char); + GET_COLUMN (Npolyterms, "NPOLYTERMS", byte); GET_COLUMN (dXpixSys, "XPIX_SYS_ERR", float); GET_COLUMN (dYpixSys, "YPIX_SYS_ERR", float); + GET_COLUMN (refColorBlue, "REF_COLOR_BLUE", float); + GET_COLUMN (refColorRed, "REF_COLOR_RED", float); GET_COLUMN (imageID, "ID", int); GET_COLUMN (nFitAstrom, "NFIT", int); @@ -85,20 +87,23 @@ ALLOCATE (image_pos, ImagePos, Nrow); for (i = 0; i < Nrow; i++) { - image_pos[i].coords.crval1 = crval1 [i]; - image_pos[i].coords.crval2 = crval2 [i]; - image_pos[i].coords.crpix1 = crpix1 [i]; - image_pos[i].coords.crpix2 = crpix2 [i]; - image_pos[i].coords.cdelt1 = cdelt1 [i]; - image_pos[i].coords.cdelt2 = cdelt2 [i]; - image_pos[i].coords.pc1_1 = pc1_1 [i]; - image_pos[i].coords.pc1_2 = pc1_2 [i]; - image_pos[i].coords.pc2_1 = pc2_1 [i]; - image_pos[i].coords.pc2_2 = pc2_2 [i]; - image_pos[i].coords.Npolyterms = Npolyterms[i]; - image_pos[i].dXpixSys = dXpixSys [i]; - image_pos[i].dYpixSys = dYpixSys [i]; - image_pos[i].imageID = imageID [i]; - image_pos[i].nFitAstrom = nFitAstrom[i]; - image_pos[i].flags = flags [i]; + InitCoords(&image_pos[i].coords, NULL); + image_pos[i].coords.crval1 = crval1 [i]; + image_pos[i].coords.crval2 = crval2 [i]; + image_pos[i].coords.crpix1 = crpix1 [i]; + image_pos[i].coords.crpix2 = crpix2 [i]; + image_pos[i].coords.cdelt1 = cdelt1 [i]; + image_pos[i].coords.cdelt2 = cdelt2 [i]; + image_pos[i].coords.pc1_1 = pc1_1 [i]; + image_pos[i].coords.pc1_2 = pc1_2 [i]; + image_pos[i].coords.pc2_1 = pc2_1 [i]; + image_pos[i].coords.pc2_2 = pc2_2 [i]; + image_pos[i].coords.Npolyterms = Npolyterms [i]; + image_pos[i].dXpixSys = dXpixSys [i]; + image_pos[i].dYpixSys = dYpixSys [i]; + image_pos[i].refColorBlue = refColorBlue[i]; + image_pos[i].refColorRed = refColorRed [i]; + image_pos[i].imageID = imageID [i]; + image_pos[i].nFitAstrom = nFitAstrom [i]; + image_pos[i].flags = flags [i]; // polyterms and ctype are a bit different @@ -108,22 +113,24 @@ fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow); - free (crval1 ); - free (crval2 ); - free (crpix1 ); - free (crpix2 ); - free (cdelt1 ); - free (cdelt2 ); - free (pc1_1 ); - free (pc1_2 ); - free (pc2_1 ); - free (pc2_2 ); - free (polyterms ); - free (ctype ); - free (Npolyterms); - free (dXpixSys ); - free (dYpixSys ); - free (imageID ); - free (nFitAstrom); - free (flags ); + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms ); + free (dXpixSys ); + free (dYpixSys ); + free (refColorBlue); + free (refColorRed ); + free (imageID ); + free (nFitAstrom ); + free (flags ); *nimage_pos = Nrow; @@ -153,44 +160,48 @@ gfits_create_table_header (&theader, "BINTABLE", "IMAGE_POS"); - gfits_define_bintable_column (&theader, "D", "CRVAL1", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "CRVAL2", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "CRPIX1", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "CRPIX2", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "CDELT1", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "CDELT2", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "PC1_1", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "PC1_2", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "PC2_1", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "PC2_2", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "14E", "POLYTERMS", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "15A", "CTYPE", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "L", "NPOLYTERMS", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "XPIX_SYS_ERR", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "E", "YPIX_SYS_ERR", "word", "unit", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); - gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "CRVAL1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "CRVAL2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CRPIX2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "CDELT2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC1_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_1", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "PC2_2", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "14E", "POLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "15A", "CTYPE", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "L", "NPOLYTERMS", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "XPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "YPIX_SYS_ERR", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "REF_COLOR_BLUE", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "E", "REF_COLOR_RED", "word", "unit", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "ID", "image ID", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "NFIT", "number of fitted stars", "unitless", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "FLAGS", "analysis flags", "unitless", 1.0, 0.0); // generate the output array that carries the data gfits_create_table (&theader, &ftable); - double *crval1 ; - double *crval2 ; - float *crpix1 ; - float *crpix2 ; - float *cdelt1 ; - float *cdelt2 ; - float *pc1_1 ; - float *pc1_2 ; - float *pc2_1 ; - float *pc2_2 ; - float *polyterms ; - char *ctype ; - char *Npolyterms; - float *dXpixSys ; - float *dYpixSys ; - int *imageID ; - int *nFitAstrom; - int *flags ; + double *crval1 ; + double *crval2 ; + float *crpix1 ; + float *crpix2 ; + float *cdelt1 ; + float *cdelt2 ; + float *pc1_1 ; + float *pc1_2 ; + float *pc2_1 ; + float *pc2_2 ; + float *polyterms ; + char *ctype ; + char *Npolyterms ; + float *dXpixSys ; + float *dYpixSys ; + float *refColorBlue; + float *refColorRed ; + int *imageID ; + int *nFitAstrom ; + int *flags ; // create intermediate storage arrays @@ -210,4 +221,6 @@ ALLOCATE (dXpixSys , float , Nimage_pos); ALLOCATE (dYpixSys , float , Nimage_pos); + ALLOCATE (refColorBlue, float , Nimage_pos); + ALLOCATE (refColorRed , float , Nimage_pos); ALLOCATE (imageID , int , Nimage_pos); ALLOCATE (nFitAstrom, int , Nimage_pos); @@ -216,20 +229,22 @@ // assign the storage arrays for (i = 0; i < Nimage_pos; i++) { - crval1 [i] = image_pos[i].coords.crval1 ; - crval2 [i] = image_pos[i].coords.crval2 ; - crpix1 [i] = image_pos[i].coords.crpix1 ; - crpix2 [i] = image_pos[i].coords.crpix2 ; - cdelt1 [i] = image_pos[i].coords.cdelt1 ; - cdelt2 [i] = image_pos[i].coords.cdelt2 ; - pc1_1 [i] = image_pos[i].coords.pc1_1 ; - pc1_2 [i] = image_pos[i].coords.pc1_2 ; - pc2_1 [i] = image_pos[i].coords.pc2_1 ; - pc2_2 [i] = image_pos[i].coords.pc2_2 ; - Npolyterms[i] = image_pos[i].coords.Npolyterms; - dXpixSys [i] = image_pos[i].dXpixSys ; - dYpixSys [i] = image_pos[i].dYpixSys ; - imageID [i] = image_pos[i].imageID ; - nFitAstrom[i] = image_pos[i].nFitAstrom ; - flags [i] = image_pos[i].flags ; + crval1 [i] = image_pos[i].coords.crval1 ; + crval2 [i] = image_pos[i].coords.crval2 ; + crpix1 [i] = image_pos[i].coords.crpix1 ; + crpix2 [i] = image_pos[i].coords.crpix2 ; + cdelt1 [i] = image_pos[i].coords.cdelt1 ; + cdelt2 [i] = image_pos[i].coords.cdelt2 ; + pc1_1 [i] = image_pos[i].coords.pc1_1 ; + pc1_2 [i] = image_pos[i].coords.pc1_2 ; + pc2_1 [i] = image_pos[i].coords.pc2_1 ; + pc2_2 [i] = image_pos[i].coords.pc2_2 ; + Npolyterms [i] = image_pos[i].coords.Npolyterms; + dXpixSys [i] = image_pos[i].dXpixSys ; + dYpixSys [i] = image_pos[i].dYpixSys ; + refColorBlue[i] = image_pos[i].refColorBlue ; + refColorRed [i] = image_pos[i].refColorRed ; + imageID [i] = image_pos[i].imageID ; + nFitAstrom [i] = image_pos[i].nFitAstrom ; + flags [i] = image_pos[i].flags ; // polyterms and ctype are a bit different @@ -239,41 +254,45 @@ // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "CRVAL1", crval1 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CRVAL2", crval2 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CRPIX1", crpix1 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CRPIX2", crpix2 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CDELT1", cdelt1 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CDELT2", cdelt2 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "PC1_1", pc1_1 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "PC1_2", pc1_2 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "PC2_1", pc2_1 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "PC2_2", pc2_2 , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "POLYTERMS", polyterms , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "CTYPE", ctype , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS", Npolyterms, Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR",dXpixSys , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR",dYpixSys , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitAstrom, Nimage_pos); - gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_pos); - - free (crval1 ); - free (crval2 ); - free (crpix1 ); - free (crpix2 ); - free (cdelt1 ); - free (cdelt2 ); - free (pc1_1 ); - free (pc1_2 ); - free (pc2_1 ); - free (pc2_2 ); - free (polyterms ); - free (ctype ); - free (Npolyterms); - free (dXpixSys ); - free (dYpixSys ); - free (imageID ); - free (nFitAstrom); - free (flags ); + gfits_set_bintable_column (&theader, &ftable, "CRVAL1", crval1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRVAL2", crval2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX1", crpix1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CRPIX2", crpix2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT1", cdelt1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CDELT2", cdelt2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_1", pc1_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC1_2", pc1_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_1", pc2_1 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "PC2_2", pc2_2 , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "POLYTERMS", polyterms , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "CTYPE", ctype , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NPOLYTERMS", Npolyterms, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "XPIX_SYS_ERR", dXpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "YPIX_SYS_ERR", dYpixSys , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "REF_COLOR_BLUE", refColorBlue, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "REF_COLOR_RED", refColorRed , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "ID", imageID , Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "NFIT", nFitAstrom, Nimage_pos); + gfits_set_bintable_column (&theader, &ftable, "FLAGS", flags , Nimage_pos); + + free (crval1 ); + free (crval2 ); + free (crpix1 ); + free (crpix2 ); + free (cdelt1 ); + free (cdelt2 ); + free (pc1_1 ); + free (pc1_2 ); + free (pc2_1 ); + free (pc2_2 ); + free (polyterms ); + free (ctype ); + free (Npolyterms ); + free (dXpixSys ); + free (dYpixSys ); + free (refColorBlue); + free (refColorRed ); + free (imageID ); + free (nFitAstrom ); + free (flags ); FILE *f = fopen (filename, "w"); Index: trunk/Ohana/src/relastro/src/MeasPosIO.c =================================================================== --- trunk/Ohana/src/relastro/src/MeasPosIO.c (revision 37807) +++ trunk/Ohana/src/relastro/src/MeasPosIO.c (revision 37807) @@ -0,0 +1,191 @@ +# include "relastro.h" + +# define GET_COLUMN(OUT,NAME,TYPE) \ + TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + myAssert (!strcmp(type, #TYPE), "wrong column type"); + +// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client. +// Here, we use a handful of different columns (if not, we could move to libdvo) +MeasPos *MeasPosLoad(char *filename, off_t *nmeaspos) { + + int i, Ncol; + off_t Nrow; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + *nmeaspos = 0; + MeasPos *measpos = NULL; + + FILE *f = fopen (filename, "r"); + if (!f) { + fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename); + return NULL; + } + + /* load in PHU segment (ignore) */ + if (!gfits_fread_header (f, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset header\n"); + fclose (f); + return NULL; + } + if (!gfits_fread_matrix (f, &matrix, &header)) { + if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n"); + gfits_free_header (&header); + fclose (f); + return NULL; + } + + ftable.header = &theader; + + // load data for this header + if (!gfits_load_header (f, &theader)) { + fclose (f); + return NULL; + } + + // read the fits table bytes + if (!gfits_fread_ftable_data (f, &ftable, FALSE)) { + fclose (f); + return (NULL); + } + fclose (f); + + // a bit annoying : we read the entire block of data, then extract the columns, then set the image structure values. + // this means I need 3 copies in memory at some point. ugh. + + char type[16]; + + GET_COLUMN (R, "RA", double); + GET_COLUMN (D, "DEC", double); + GET_COLUMN (objID, "OBJ_ID", int); + GET_COLUMN (catID, "CAT_ID", int); + GET_COLUMN (imageID, "IMAGE_ID", int); + + ALLOCATE (measpos, MeasPos, Nrow); + for (i = 0; i < Nrow; i++) { + measpos[i].R = R [i]; + measpos[i].D = D [i]; + measpos[i].objID = objID[i]; + measpos[i].catID = catID[i]; + measpos[i].imageID = imageID[i]; + } + fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow); + + free (R ); + free (D ); + free (objID); + free (catID); + free (imageID); + + // free FITS table pieces... + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + *nmeaspos = Nrow; + return measpos; +} + +// STATUS is value expected for success +# define CHECK_STATUS(STATUS,MSG,...) \ + if (!(STATUS)) { \ + fprintf (stderr, MSG, __VA_ARGS__); \ + return FALSE; \ + } + +int MeasPosSave(char *filename, MeasPos *measpos, off_t Nmeaspos) { + + int i; + Header header; + Header theader; + Matrix matrix; + FTable ftable; + + gfits_init_header (&header); + header.extend = TRUE; + gfits_create_header (&header); + gfits_create_matrix (&header, &matrix); + + gfits_create_table_header (&theader, "BINTABLE", "MEASPOS"); + + gfits_define_bintable_column (&theader, "D", "RA", "meas position, ra", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "D", "DEC", "meas position, dec", "degrees", 1.0, 0.0); + gfits_define_bintable_column (&theader, "J", "OBJ_ID", "object ID", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "CAT_ID", "catalog ID", NULL, 1.0, FT_BZERO_INT32); + gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, FT_BZERO_INT32); + + // generate the output array that carries the data + gfits_create_table (&theader, &ftable); + + double *R, *D; + unsigned int *objID, *catID, *imageID; + + // create intermediate storage arrays + ALLOCATE (R, double, Nmeaspos); + ALLOCATE (D, double, Nmeaspos); + ALLOCATE (objID, unsigned int, Nmeaspos); + ALLOCATE (catID, unsigned int, Nmeaspos); + ALLOCATE (imageID, unsigned int, Nmeaspos); + + // assign the storage arrays + for (i = 0; i < Nmeaspos; i++) { + R [i] = measpos[i].R ; + D [i] = measpos[i].D ; + objID[i] = measpos[i].objID; + catID[i] = measpos[i].catID; + imageID[i] = measpos[i].imageID; + } + + // add the columns to the output array + gfits_set_bintable_column (&theader, &ftable, "RA", R, Nmeaspos); + gfits_set_bintable_column (&theader, &ftable, "DEC", D, Nmeaspos); + gfits_set_bintable_column (&theader, &ftable, "OBJ_ID", objID, Nmeaspos); + gfits_set_bintable_column (&theader, &ftable, "CAT_ID", catID, Nmeaspos); + gfits_set_bintable_column (&theader, &ftable, "IMAGE_ID", imageID, Nmeaspos); + + free (R ); + free (D ); + free (objID); + free (catID); + free (imageID); + + FILE *f = fopen (filename, "w"); + if (!f) { + fprintf (stderr, "ERROR: cannot open measmag file for output %s\n", filename); + return FALSE; + } + + int status; + status = gfits_fwrite_header (f, &header); + CHECK_STATUS (status, "ERROR: cannot write header for measpos %s\n", filename); + + status = gfits_fwrite_matrix (f, &matrix); + CHECK_STATUS (status, "ERROR: cannot write matrix for measpos %s\n", filename); + + status = gfits_fwrite_Theader (f, &theader); + CHECK_STATUS (status, "ERROR: cannot write table header for measpos %s\n", filename); + + status = gfits_fwrite_table (f, &ftable); + CHECK_STATUS (status, "ERROR: cannot write table data for measpos %s\n", filename); + + gfits_free_header (&header); + gfits_free_matrix (&matrix); + gfits_free_header (&theader); + gfits_free_table (&ftable); + + int fd = fileno (f); + + status = fflush (f); + CHECK_STATUS (!status, "ERROR: cannot flush file measpos %s\n", filename); + + status = fsync (fd); + CHECK_STATUS (!status, "ERROR: cannot flush file measpos %s\n", filename); + + status = fclose (f); + CHECK_STATUS (!status, "ERROR: problem closing measpos file %s\n", filename); + + return TRUE; +} Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 37261) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 37807) @@ -176,2 +176,20 @@ } +int printStarMap(int N, char *filename) { + + off_t i, Nimages; + Image *images = getimages(&Nimages, NULL); + double L, M; + + FILE *f = fopen (filename, "w"); + + for (i = 0; i < starmap[N].Npoints; i++) { + // set the transformed coordinates + XY_to_LM (&L, &M, starmap[N].points[i].X, starmap[N].points[i].Y, &images[N].coords); + + fprintf (f, "%d %7.2f %f : %7.2f %7.2f : %7.2f %7.2f : %7.2f %7.2f\n", (int) i, starmap[N].points[i].X, starmap[N].points[i].Y, starmap[N].points[i].L, starmap[N].points[i].M, L, M, starmap[N].points[i].L - L, starmap[N].points[i].M - M); + } + fclose (f); + return (TRUE); +} + Index: trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateChips.c (revision 37261) +++ trunk/Ohana/src/relastro/src/UpdateChips.c (revision 37807) @@ -36,5 +36,4 @@ off_t i, Nraw, Nref, nFitAstr; StarData *raw, *ref; - Coords *oldCoords; float dXpixSys, dYpixSys; double *Ro, *Do; @@ -74,4 +73,10 @@ Nskip ++; mode[i] = 0; + continue; + } + if (Nraw <= IMFIT_TOO_FEW) { + Nskip ++; + mode[i] = 0; + free (raw); continue; } @@ -90,5 +95,7 @@ // save these in case of failure - saveCoords (&image[i].coords, i); + Coords oldCoords; + SaveCoords (&oldCoords, &image[i].coords); + dXpixSys = image[i].dXpixSys; dYpixSys = image[i].dYpixSys; @@ -100,36 +107,47 @@ if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); - // restore status quo ante - oldCoords = getCoords (i); - memcpy (&image[i].coords, oldCoords, sizeof(Coords)); - image[i].dXpixSys = dXpixSys; - image[i].dYpixSys = dYpixSys; - image[i].nFitAstrom = nFitAstr; - - saveCenter (image, &Ro[i], &Do[i], i); - mode[i] = 1; - NoldFit ++; - free (raw); - free (ref); - continue; + if (1) { + // restore status quo ante (replace truMap with tmpMap) + RestoreCoords (&image[i].coords, &oldCoords, &image[i]); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &Ro[i], &Do[i], i); + mode[i] = 1; + NoldFit ++; + free (raw); + free (ref); + continue; + } } if (!checkStarMap (i)) { if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); - // restore status quo ante - oldCoords = getCoords (i); - memcpy (&image[i].coords, oldCoords, sizeof(Coords)); - image[i].dXpixSys = dXpixSys; - image[i].dYpixSys = dYpixSys; - image[i].nFitAstrom = nFitAstr; - - saveCenter (image, &Ro[i], &Do[i], i); - mode[i] = 2; - image[i].flags |= ID_IMAGE_ASTROM_POOR; - NoldFit ++; - free (raw); - free (ref); - continue; + + if (1) { + // restore status quo ante (replace truMap with tmpMap) + RestoreCoords (&image[i].coords, &oldCoords, &image[i]); + image[i].dXpixSys = dXpixSys; + image[i].dYpixSys = dYpixSys; + image[i].nFitAstrom = nFitAstr; + + saveCenter (image, &Ro[i], &Do[i], i); + mode[i] = 2; + image[i].flags |= ID_IMAGE_ASTROM_POOR; + NoldFit ++; + free (raw); + free (ref); + continue; + } } + + // Apply the modified coords back to the measure.R,D. Note that raw.R,D, ref.L,M, etc + // are all automatically updated in this block because they are re-generated from + // image.coords on each pass. + setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC); + if (USE_GALAXY_MODEL) { + image[i].flags |= ID_IMAGE_ASTROM_GMM; + } saveCenter (image, &Ro[i], &Do[i], i); @@ -235,5 +253,4 @@ off_t Nraw, Nref, nFitAstr; StarData *raw, *ref; - Coords *oldCoords; float dXpixSys, dYpixSys; @@ -260,4 +277,10 @@ threadinfo->Nskip ++; threadinfo->mode[i] = 0; + continue; + } + if (Nraw <= IMFIT_TOO_FEW) { + threadinfo->Nskip ++; + threadinfo->mode[i] = 0; + free (raw); continue; } @@ -276,5 +299,6 @@ // save these in case of failure - saveCoords (&image[i].coords, i); + Coords oldCoords; + SaveCoords (&oldCoords, &image[i].coords); dXpixSys = image[i].dXpixSys; dYpixSys = image[i].dYpixSys; @@ -286,7 +310,6 @@ if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name, i, Nraw, image[i].nstar); - // restore status quo ante - oldCoords = getCoords (i); - memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + // restore status quo ante (replace truMap with tmpMap) + RestoreCoords (&image[i].coords, &oldCoords, &image[i]); image[i].dXpixSys = dXpixSys; image[i].dYpixSys = dYpixSys; @@ -303,7 +326,7 @@ if (!checkStarMap (i)) { if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name, i, Nraw); - // restore status quo ante - oldCoords = getCoords (i); - memcpy (&image[i].coords, oldCoords, sizeof(Coords)); + + // restore status quo ante (replace truMap with tmpMap) + RestoreCoords (&image[i].coords, &oldCoords, &image[i]); image[i].dXpixSys = dXpixSys; image[i].dYpixSys = dYpixSys; @@ -319,4 +342,7 @@ } + // apply the modified R,D back to the measures + setImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, raw, Nraw, MODE_MOSAIC); + saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i); threadinfo->mode[i] = 3; @@ -333,11 +359,18 @@ static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER; +void lockUpdateChips () { + pthread_mutex_lock (&UpdateChips_mutex); +} +void unlockUpdateChips () { + pthread_mutex_unlock (&UpdateChips_mutex); +} + // we have an array of chips (image, Nimage). we need to hand out images one at a time to // the worker threads as they need off_t getNextImageForThread () { - pthread_mutex_lock (&UpdateChips_mutex); + lockUpdateChips (); if (nextImage >= Nimage) { - pthread_mutex_unlock (&UpdateChips_mutex); + unlockUpdateChips (); return (-1); } @@ -345,5 +378,5 @@ nextImage ++; - pthread_mutex_unlock (&UpdateChips_mutex); + unlockUpdateChips (); return (thisImage); } @@ -494,3 +527,48 @@ } -// XXX if (!FindMosaicForImage (image, Nimage, i)) { } +// save the coords and offset map in src at tgt +void SaveCoords (Coords *tgt, Coords *src) { + myAssert (tgt, "oops"); + myAssert (src, "oops"); + CopyCoords (tgt, src); // src retains a pointer to the AstromOffsetMap table, to be saved + tgt->offsetMap = AstromOffsetMapCopy(src->offsetMap); // oldCoords now saves the old values in a new structure +} + +// restore the coords and offset map at tgt from src +void RestoreCoords (Coords *tgt, Coords *src, Image *image) { + + AstromOffsetMap *truMap = tgt->offsetMap; + AstromOffsetMap *tmpMap = src->offsetMap; + + // copy the coords structure data first, then fix the offsetMap pointers + // this modifies tgt.coords, but we need to keep the original pointer + CopyCoords (tgt, src); + tgt->offsetMap = truMap; + + // we want to keep the old solution, which is (maybe) a pointer to a stand-alone + // map. the image has a bad solution, but it is a pointer to the map in the I/O + // table. those values will be saved on exit. + + if (!tmpMap && !truMap) return; + + if (tmpMap && truMap) { + AstromOffsetMapSetOrder (truMap, tmpMap->Nx, tmpMap->Ny, image); + AstromOffsetMapCopyData (truMap, tmpMap); + return; + } + + if (tmpMap && !truMap) { + // is this even possible? + lockUpdateChips (); + AstromOffsetTable *table = get_astrom_table (); + AstromOffsetTableNewMap(table, tmpMap->Nx, tmpMap->Ny, image); // registers the map with the image + unlockUpdateChips (); + AstromOffsetMapCopyData (image->coords.offsetMap, tmpMap); + return; + } + if (!tmpMap && truMap) { + truMap->keep = FALSE; + return; + } + myAbort ("oops"); +} Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 37261) +++ trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 37807) @@ -1,8 +1,142 @@ # include "relastro.h" +void dump_measures (Average *average, Measure *measure); // in ImageOps.c -void printNcatTotal (); +// this function is called by UpdateObjectOffsets or relastro_images.c for reset +int UpdateMeasures (Catalog *catalog, int Ncatalog) { + + off_t i, j, Nimage; + Image *image; + double DPOS_MAX_ASEC; + + int badImage = + ID_IMAGE_ASTROM_NOCAL | + ID_IMAGE_ASTROM_POOR | + ID_IMAGE_ASTROM_FAIL | + ID_IMAGE_ASTROM_SKIP | + ID_IMAGE_ASTROM_FEW; + + if (RESET_BAD_IMAGES) badImage = 0; + + image = getimages (&Nimage, NULL); + + // track measurements which are far from their original position (poor at pole) + int NoffRAave = 0; int NoffRAori = 0; + int NoffDECave = 0; int NoffDECori = 0; + + for (i = 0; i < Ncatalog; i++) { + myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error"); + for (j = 0; j < catalog[i].Nmeasure; j++) { + MeasureTiny *measureT = &catalog[i].measureT[j]; + Measure *measureB = catalog[i].measure ? &catalog[i].measure[j] : NULL; + + off_t im = getImageByID (measureT->imageID); + if (im < 0) continue; // detections without imageIDs are not associated with a known image (eg, refs) + + // check that we have the right image: + myAssert (measureT->photcode == image[im].photcode, "image photcode mismatch"); + myAssert (measureT->t >= image[im].tzero, "image time mismatch (1)"); + myAssert (measureT->t <= image[im].tzero + image[im].NX*image[im].trate/5000.0, "image time mismatch (2)"); + // I'm allowing the trate*NX to be a factor of 2x too small + + // skip measurements on images that have failed solutions (divergent or otherwise) + // XXX apply this or not?? + if (FALSE && (image[im].flags & badImage)) continue; + + Coords *moscoords = image[im].coords.mosaic; + Coords *imcoords = &image[im].coords; + + if (moscoords) { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords->cdelt1, moscoords->cdelt2); + } else { + DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords->cdelt1, imcoords->cdelt2); + } + + double X = measureT->Xccd; + double Y = measureT->Yccd; + if (USE_FIXED_PIXCOORDS) { + if (isfinite(measureT->Xfix) && isfinite(measureT->Yfix)) { + float dX = measureT->Xfix - measureT->Xccd; + float dY = measureT->Yfix - measureT->Yccd; + if (hypot(dX,dY) < 2.0) { + X = measureT->Xfix; + Y = measureT->Yfix; + } + } + } + + double R, D; + XY_to_RD (&R, &D, X, Y, imcoords); + + int n = measureT[0].averef; + Average *average = &catalog[i].average[n]; + + double dR = 3600.0*(average[0].R - R); + double dD = 3600.0*(average[0].D - D); + + // make sure detection is on the same side of the 0,360 boundary as the average + // this will give some funny results withing ~1 arcsec of the pol + if (dR > +180.0*3600.0) { + // average on high end of boundary, move star up + R += 360.0; + dR = 3600.0*(average[0].R - R); + } + if (dR < -180.0*3600.0) { + // average on low end of boundary, move star down + R -= 360.0; + dR = 3600.0*(average[0].R - R); + } + + float csdec = cos(average[0].D * RAD_DEG); + + // complain if the new location is far from the average location + // NOTE: This should never happen, or our StarMap tests are not working + if (fabs(dR*csdec) > 3.0*ADDSTAR_RADIUS) { + NoffRAave ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, dR*csdec, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + if (fabs(dD) > 3.0*ADDSTAR_RADIUS) { + NoffDECave ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", average[0].R, average[0].D, dR, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + + // complain if the new location is far from the old location + if (fabs(csdec*(measureT[0].R - R)) > DPOS_MAX_ASEC) { + NoffRAori ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement moves far from original location (R): %f %f (%f %f %f %f)\n", average[0].R, average[0].D, measureT[0].R, dR, csdec*(measureT[0].R - R), dD); + dump_measures (&average[0], catalog[i].measure); + } + } + if (fabs(measureT[0].D - D) > DPOS_MAX_ASEC) { + NoffDECori ++; + if (VERBOSE2 && catalog[i].measure) { + fprintf (stderr, "measurement moves far from original location (D): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, measureT[0].D, dD); + dump_measures (&average[0], catalog[i].measure); + } + } + measureT->R = R; + measureT->D = D; + if (measureB) { + measureB->R = R; + measureB->D = D; + } + } + } + // printNcatTotal (); + + int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori; + if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori); + + return (TRUE); +} // this function operates on Measure, not MeasureTiny -int UpdateMeasures (Catalog *catalog, int Ncatalog) { +int UpdateMeasuresOld (Catalog *catalog, int Ncatalog) { off_t i, Nimage; @@ -32,5 +166,5 @@ } - printNcatTotal (); + // printNcatTotal (); return (TRUE); Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 37261) +++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 37807) @@ -59,5 +59,6 @@ // uses MeasureTiny values, but will also update Measure values when populated - UpdateObjects (&catalog, 1); + // if Galaxy Motion Model correction is active, do not apply at this last stage + UpdateObjects (&catalog, 1, FALSE); free_tiny_values(&catalog); @@ -74,7 +75,10 @@ save_catalogs (&catalog, 1); } - + printNcatTotal(); return (TRUE); } + +int UpdateObjectOffsets_parallel_group (HostTableGroup *group, SkyList *sky); +int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky); // CATDIR is supplied globally @@ -90,4 +94,148 @@ exit (1); } + +# if (0) + + UpdateObjectOffsets_parallel_table (table, sky); + +# else + + int Ngroups; + HostTableGroup *groups = HostTableGroups (table, &Ngroups); + // split the table into Ngroups, each with a unique set of machines (this avoids overloading the remote host machines) + + int i; + for (i = 0; i < Ngroups; i++) { + // update only a group of unique machines at a time + UpdateObjectOffsets_parallel_group (&groups[i], sky); + } +# endif + + return TRUE; +} + +// CATDIR is supplied globally +# define DEBUG 1 +int UpdateObjectOffsets_parallel_group (HostTableGroup *group, SkyList *sky) { + + // launch the relastro_client jobs to the parallel hosts + + int i, j; + for (i = 0; i < group->Nhosts; i++) { + + if (sky->Nregions < group->Nhosts) { + // do any of the regions want this host? + int wantThisHost = FALSE; + for (j = 0; j < sky->Nregions; j++) { + if (HostTableTestHost (sky->regions[j], group->hosts[i][0].hostID)) { + wantThisHost = TRUE; + break; + } + } + if (!wantThisHost) { + // fprintf (stderr, "skip host %s\n", group->hosts[i][0].hostname); + continue; + } + } + + // ensure that the paths are absolute path names + char *tmppath = abspath (group->hosts[i][0].pathname, DVO_MAX_PATH); + free (group->hosts[i][0].pathname); + group->hosts[i][0].pathname = tmppath; + + // options / arguments that can affect relastro_client -load: + // VERBOSE, VERBOSE2 + // RESET (-reset) + // TimeSelect -time + // (note that psfQF is applied rigidly at 0.85, as is the galaxy test) + // ImagSelect, ImagMin, ImagMax + // MaxDensityUse, MaxDensityValue + + char command[1024]; + snprintf (command, 1024, "relastro_client -update-offsets -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s -minerror %f", + group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE, MIN_ERROR); + + if (FIT_MODE == FIT_PM_ONLY) strextend (command, "-pm"); + if (FIT_MODE == FIT_PAR_ONLY) strextend (command, "-par"); + if (FIT_MODE == FIT_PM_AND_PAR) strextend (command, "-pmpar"); + + if (VERBOSE) strextend (command, "-v"); + if (VERBOSE2) strextend (command, "-vv"); + if (RESET) strextend (command, "-reset"); + if (UPDATE) strextend (command, "-update"); + + if (RESET_BAD_IMAGES) strextend (command, "-reset-bad-images"); + + if (ImagSelect) strextend (command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (command, "-max-density %f", MaxDensityValue); + + if (USE_BASIC_CHECK) strextend (command, "-basic-image-search"); + if (FlagOutlier) strextend (command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); + + if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); + + if (PHOTCODE_KEEP_LIST) strextend (command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PhotFlagSelect) strextend (command, "+photflags"); + if (PhotFlagBad) strextend (command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { + strextend (command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + } + if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { + strextend (command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); + } + + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); + if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); + + if (TimeSelect) { + char *tstart = ohana_sec_to_date (TSTART); + char *tstop = ohana_sec_to_date (TSTOP); + strextend (command, "-time %s %s", tstart, tstop); + free (tstart); + free (tstop); + } + fprintf (stderr, "command: %s\n", command); + + if (PARALLEL_MANUAL) continue; + + if (PARALLEL_SERIAL) { + int status = system (command); + if (status) { + fprintf (stderr, "ERROR running relastro_client\n"); + exit (2); + } + } else { + // launch the job on the remote machine (no handshake) + int errorInfo = 0; + int pid = rconnect ("ssh", group->hosts[i][0].hostname, command, group->hosts[i][0].stdio, &errorInfo, FALSE); + if (!pid) { + if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", group->hosts[i][0].hostname, errorInfo); + exit (1); + } + group->hosts[i][0].pid = pid; // save for future reference + } + } + + if (PARALLEL_MANUAL) { + fprintf (stderr, "run the relastro_client commands above. when these are done, hit return\n"); + getchar(); + } + if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) { + HostTableGroupWaitJobsGetIO (group, __FILE__, __LINE__, VERBOSE); + } + + return TRUE; +} + +// CATDIR is supplied globally +# define DEBUG 1 +int UpdateObjectOffsets_parallel_table (HostTable *table, SkyList *sky) { + + // launch the relastro_client jobs to the parallel hosts int i, j; Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 37261) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 37807) @@ -2,5 +2,5 @@ # define PAR_TOOFEW 5 -int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m); +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset); int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats); @@ -69,13 +69,6 @@ /* project coordinates to a plane centered on the object with units of arcsec */ - coords.crval1 = 0; - coords.crval2 = 0; - coords.crpix1 = 0; - coords.crpix2 = 0; + InitCoords (&coords, "DEC--SIN"); coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; - coords.pc1_1 = coords.pc2_2 = 1.0; - coords.pc1_2 = coords.pc2_1 = 0.0; - coords.Npolyterms = 1; - strcpy (coords.ctype, "DEC--SIN"); // use J2000 as a reference time @@ -106,5 +99,5 @@ // This function operates on both Measure and MeasureTiny. In the big stages, this should // be called with just MeasureTiny set and Measure == NULL -int UpdateObjects (Catalog *catalog, int Ncatalog) { +int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) { initObjectData (catalog, Ncatalog); @@ -137,8 +130,8 @@ UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack); - UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m); - } - if (VERBOSE) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsChips.Nave, fitStatsChips.Npm, fitStatsChips.Npar, fitStatsChips.Nskip, fitStatsChips.Noffset); - if (VERBOSE) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsStack.Nave, fitStatsStack.Npm, fitStatsStack.Npar, fitStatsStack.Nskip, fitStatsStack.Noffset); + UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop); + } + if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsChips.Nave, fitStatsChips.Npm, fitStatsChips.Npar, fitStatsChips.Nskip, fitStatsChips.Noffset); + if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n", i, fitStatsStack.Nave, fitStatsStack.Npm, fitStatsStack.Npar, fitStatsStack.Nskip, fitStatsStack.Noffset); sumFitStats (&fitStatsChips, &sumStatsChips); sumFitStats (&fitStatsStack, &sumStatsStack); @@ -146,6 +139,6 @@ freeObjectData (); - if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar), sumStatsChips.Nave, sumStatsChips.Npm, sumStatsChips.Npar, sumStatsChips.Nskip, sumStatsChips.Noffset); - if (VERBOSE) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar), sumStatsStack.Nave, sumStatsStack.Npm, sumStatsStack.Npar, sumStatsStack.Nskip, sumStatsStack.Noffset); + if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar), sumStatsChips.Nave, sumStatsChips.Npm, sumStatsChips.Npar, sumStatsChips.Nskip, sumStatsChips.Noffset); + if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n", (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar), sumStatsStack.Nave, sumStatsStack.Npm, sumStatsStack.Npar, sumStatsStack.Nskip, sumStatsStack.Noffset); return (TRUE); } @@ -153,5 +146,5 @@ // This function operates on both Measure and MeasureTiny. In the big stages, this should // be called with just MeasureTiny set and Measure == NULL -int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m) { +int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) { int setRefColor = areImagesMatched(); @@ -215,4 +208,14 @@ double Ri = getMeanR (&measure[k], average, secfilt); double Di = getMeanD (&measure[k], average, secfilt); + + // if we are correcting for the Galaxy Motion Model, only should apply it here + // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and + // (b) if we think the measure R,D is already at the image epoch position + if (USE_GALAXY_MODEL && applyGalaxyOffset) { + Ri -= measure[k].RoffGAL / 3600.0; + Di -= measure[k].DoffGAL / 3600.0; + } + + // XXX add in dR,dD GAL here // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius @@ -329,4 +332,14 @@ // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq + // project all of the R,D coordinates to a plane centered on this coordinate. set + // the times to be relative to Tmean (this is required for parallax as well) + for (k = 0; k < N; k++) { + RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + T[k] -= Tmean; + if (XVERB) { + fprintf (stderr, OFF_T_FMT" %f %f %f %f %f +/- %f %f\n", k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); + } + } + // *** first fit for the proper motion (skip fit if Trange or Npts is too small) *** if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { @@ -340,14 +353,4 @@ } - // project all of the R,D coordinates to a plane centered on this coordinate. set - // the times to be relative to Tmean (this is required for parallax as well) - for (k = 0; k < N; k++) { - RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); - T[k] -= Tmean; - if (XVERB) { - fprintf (stderr, OFF_T_FMT" %f %f %f %f %f +/- %f %f\n", k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]); - } - } - FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB); @@ -357,4 +360,6 @@ XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords); if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p); + if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + // XXX : does this make sense at 0,360 boundary? fitPM.p = fitPM.dp = 0.0; @@ -409,4 +414,6 @@ fitStats->Npar ++; + if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + // XXX a hard-wired hack... if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) { @@ -419,12 +426,13 @@ // ALWAYS fit the average model StatType statsR, statsD; - liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) - liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) - - fitAve.Ro = statsR.mean; - fitAve.dRo = 3600.0*statsR.sigma; - - fitAve.Do = statsD.mean; - fitAve.dDo = 3600.0*statsD.sigma; + liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here) + liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here) + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords); + if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do); + + fitAve.dRo = statsR.sigma; + fitAve.dDo = statsD.sigma; fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN; @@ -466,11 +474,14 @@ break; } - if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n", - average[0].R, - average[0].D, - fit.Ro, fit.Do, - 3600*(average[0].R - fit.Ro), - 3600*(average[0].D - fit.Do), - fit.uR, fit.uD, fit.p, fit.dp); + if (XVERB) { + fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n", + average[0].R, + average[0].D, + fit.Ro, fit.Do, + 3600*(average[0].R - fit.Ro), + 3600*(average[0].D - fit.Do), + fit.uR, fit.uD, fit.p, fit.dp); + fprintf (stderr, "-----\n"); + } // make sure that the fit succeeded @@ -494,4 +505,10 @@ coords.crval1 = average[0].R; coords.crval2 = average[0].D; + if (isnan(coords.crval1)) { + return (FALSE); + } + if (isnan(coords.crval2)) { + return (FALSE); + } double dXoff, dYoff; @@ -505,28 +522,4 @@ return FALSE; // XXX ?? } - - average[0].R = fit.Ro; // RA in degrees - average[0].D = fit.Do; // DEC in degrees - average[0].dR = fit.dRo; // RA scatter in arcsec - average[0].dD = fit.dDo; // DEC scatter in arcsec - - average[0].uR = fit.uR; // RA proper motion in arcsec/year - average[0].uD = fit.uD; // DEC proper motion in arcsec/year - average[0].duR = fit.duR; // RA proper motion error in arcsec/year - average[0].duD = fit.duD; // DEC proper motion error in arcsec/year - - average[0].P = fit.p; // parallax in arcsec - average[0].dP = fit.dp; // parallax error in arcsec - - average[0].ChiSqAve = fitAve.chisq; - average[0].ChiSqPM = fitPM.chisq; - average[0].ChiSqPar = fitPAR.chisq; - - average[0].Tmean = (Tmean * 86400 * 365.25) + T2000; - average[0].Trange = (Trange * 86400 * 365.25); - average[0].Npos = fit.Nfit; - - // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) - average[0].flags &= ~ID_STAR_NO_ASTROM; if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n", @@ -540,4 +533,28 @@ fitAve.chisq, fitPM.chisq, fitPAR.chisq); + average[0].R = fit.Ro; // RA in degrees + average[0].D = fit.Do; // DEC in degrees + average[0].dR = fit.dRo; // RA scatter in arcsec + average[0].dD = fit.dDo; // DEC scatter in arcsec + + average[0].uR = fit.uR; // RA proper motion in arcsec/year + average[0].uD = fit.uD; // DEC proper motion in arcsec/year + average[0].duR = fit.duR; // RA proper motion error in arcsec/year + average[0].duD = fit.duD; // DEC proper motion error in arcsec/year + + average[0].P = fit.p; // parallax in arcsec + average[0].dP = fit.dp; // parallax error in arcsec + + average[0].ChiSqAve = fitAve.chisq; + average[0].ChiSqPM = fitPM.chisq; + average[0].ChiSqPar = fitPAR.chisq; + + average[0].Tmean = (Tmean * 86400 * 365.25) + T2000; + average[0].Trange = (Trange * 86400 * 365.25); + average[0].Npos = fit.Nfit; + + // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM) + average[0].flags &= ~ID_STAR_NO_ASTROM; + return (TRUE); } Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 37261) +++ trunk/Ohana/src/relastro/src/args.c (revision 37807) @@ -265,4 +265,10 @@ remove_argument (N, &argc, argv); MaxDensityUse = TRUE; + } + + APPLY_OFFSETS = FALSE; + if ((N = get_argument (argc, argv, "-apply-offsets"))) { + remove_argument (N, &argc, argv); + APPLY_OFFSETS = TRUE; } @@ -372,4 +378,12 @@ } + CHIPMAP = 0; + if ((N = get_argument (argc, argv, "-chipmap"))) { + remove_argument (N, &argc, argv); + CHIPMAP = atoi(argv[N]); + remove_argument (N, &argc, argv); + + } + SAVEPLOT = FALSE; PLOTSTUFF = FALSE; @@ -399,9 +413,13 @@ RESET = TRUE; } - RESET_BAD_IMAGES = FALSE; if ((N = get_argument (argc, argv, "-reset-bad-images"))) { remove_argument (N, &argc, argv); RESET_BAD_IMAGES = TRUE; + } + RESET_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-reset-images"))) { + remove_argument (N, &argc, argv); + RESET_IMAGES = TRUE; } @@ -868,4 +886,5 @@ fprintf (stderr, " -high-speed (code[,code,code]) (code[,code,code]) (radius) (output catdir)\n"); fprintf (stderr, " -hpm (radius) (output catdir)\n"); + fprintf (stderr, " -chipmap (MaxOrder)\n"); fprintf (stderr, " -merge-source (objID) (catID) into (objID) (catID)\n\n"); @@ -880,4 +899,5 @@ fprintf (stderr, " -statmode (mode)\n"); fprintf (stderr, " -reset\n"); + fprintf (stderr, " -reset-images\n"); fprintf (stderr, " -nloop (N) : number of image-fit iterations\n"); fprintf (stderr, " -update : apply new fit to database\n"); Index: trunk/Ohana/src/relastro/src/assign_images.c =================================================================== --- trunk/Ohana/src/relastro/src/assign_images.c (revision 37261) +++ trunk/Ohana/src/relastro/src/assign_images.c (revision 37807) @@ -23,4 +23,16 @@ BuildChipMatch (image, Nimage); MARKTIME("build chip match for %d images: %f sec\n", (int) Nimage, dtime); + + char mapfile[DVO_MAX_PATH]; + snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR); + AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE); + + // assign images.coords.offsetMap -> table->map[i] + if (table) { + AstromOffsetTableMatchChips (image, Nimage, table); + } else { + table = AstromOffsetTableInit (); + } + put_astrom_table (table); initMosaics (image, Nimage); @@ -89,10 +101,4 @@ if (!strcmp(&image[j].coords.ctype[4], "-DIS")) continue; - // match the image to its corresponding ASTROMETRIC mosaic (not the same as the Mosaic above) - if (!FindMosaicForImage (image, Nimage, j)) { - if (VERBOSE2) fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); - continue; - } - // Exclude images with crazy astrometry // XXX NOTE : this is gpc1-specific @@ -150,13 +156,21 @@ // this is a bit memory expensive : I am making a complete copy of the image table here - off_t Nsubset = regionHosts->hosts[i].Nimage; - regionHosts->hosts[i].image[Nsubset] = image[j]; - regionHosts->hosts[i].imseq[Nsubset] = j; - - regionHosts->hosts[i].Nimage ++; - if (regionHosts->hosts[i].Nimage == regionHosts->hosts[i].NIMAGE) { - regionHosts->hosts[i].NIMAGE += D_NIMAGE; - REALLOCATE (regionHosts->hosts[i].image, Image, regionHosts->hosts[i].NIMAGE); - REALLOCATE (regionHosts->hosts[i].imseq, off_t, regionHosts->hosts[i].NIMAGE); + off_t Nsubset = host->Nimage; + host->image[Nsubset] = image[j]; + host->imseq[Nsubset] = j; + + host->Nimage ++; + if (host->Nimage == host->NIMAGE) { + host->NIMAGE += D_NIMAGE; + REALLOCATE (host->image, Image, host->NIMAGE); + REALLOCATE (host->imseq, off_t, host->NIMAGE); + } + + // save the astrometry maps, where they exist + if (image[j].coords.offsetMap) { + if (!host->astromTable) { + host->astromTable = AstromOffsetTableInit(); + } + AstromOffsetTableAddMapFromImage(host->astromTable, &image[j]); } } @@ -199,64 +213,4 @@ } -int calculate_host_image_bounds (RegionHostTable *regionHosts) { - - int i, n; - off_t j; - - for (i = 0; i < regionHosts->Nhosts; i++) { - - RegionHostInfo *host = ®ionHosts->hosts[i]; - - // XXX clear the chip match? - - BuildChipMatch (host->image, host->Nimage); - - double Rmin = 360.0; - double Rmax = 0.0; - double Dmin = +90.0; - double Dmax = -90.0; - - for (j = 0; j < host->Nimage; j++) { - - Image *image = &host->image[j]; - - if (!FindMosaicForImage (host->image, host->Nimage, j)) { - fprintf (stderr, "missing astrometry? programming error?\n"); - abort (); - } - - // define image corners - for (n = 0; n < 4; n++) { - double Xc, Yc, Rc, Dc; - Xc = Xf[n]*image->NX; - Yc = Yf[n]*image->NY; - XY_to_RD (&Rc, &Dc, Xc, Yc, &image->coords); - Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0); - - Rmin = MIN (Rmin, Rc); - Dmin = MIN (Dmin, Rc); - - Rmax = MAX (Rmax, Rc); - Dmax = MAX (Dmax, Rc); - } - - } - - host->RminCat = Rmin; - host->DminCat = Dmin; - host->RmaxCat = Rmax; - host->DmaxCat = Dmax; - - // regionHosts needs to have the full outer boundary - // (so reload_catalogs covers the correct region) - regionHosts->Rmin = MIN(Rmin, regionHosts->Rmin); - regionHosts->Rmax = MAX(Rmax, regionHosts->Rmax); - regionHosts->Dmin = MIN(Dmin, regionHosts->Dmin); - regionHosts->Dmax = MAX(Dmax, regionHosts->Dmax); - } - return TRUE; -} - - // XXX add a search tree to speed this up? int find_host_for_coords (RegionHostTable *regionHosts, double Rc, double Dc) { Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 37261) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 37807) @@ -27,13 +27,6 @@ /* for outlier rejection, project coordinates to a plane centered on the object with units of arcsec */ - coords.crval1 = 0; - coords.crval2 = 0; - coords.crpix1 = 0; - coords.crpix2 = 0; + InitCoords (&coords, "DEC--SIN"); coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; - coords.pc1_1 = coords.pc2_2 = 1.0; - coords.pc1_2 = coords.pc2_1 = 0.0; - coords.Npolyterms = 1; - strcpy (coords.ctype, "DEC--SIN"); // XXX in the future, use catalog[0].Nsecfilt only? allow catalogs to have variable Nsecfilt? @@ -72,5 +65,4 @@ if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) { continue; - // XXX can't (easily) count this for a given chip } @@ -84,4 +76,43 @@ if (RESET) { subcatalog[0].average[Naverage].flags &= ~averageBits; + } + + // if we want to do static correction to ICRF, keep the ICRF QSOs (regardless of other data quality) + if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) { + Nm = 0; + for (j = 0; j < catalog[0].average[i].Nmeasure; j++) { + + offset = catalog[0].average[i].measureOffset + j; + + CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL; + subcatalog[0].measureT[Nmeasure].averef = Naverage; + if (RESET) { + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + } + + Nmeasure ++; + Nm ++; + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + } + } + subcatalog[0].average[Naverage].Nmeasure = Nm; + Naverage ++; + if (Naverage == NAVERAGE) { + NAVERAGE += 50; + REALLOCATE (subcatalog[0].average, Average, NAVERAGE); + REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt); + } + continue; + } + + // if we want to correct to the galay model, we need to use objects with distances determined by Green & Schlafly + if (USE_GALAXY_MODEL && !catalog[0].average[i].Nstarpar) { + // we may want to add some starpar quality filters... + continue; } @@ -116,21 +147,21 @@ // exclude bogus if (ExcludeBogus) { - double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); - double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); - coords.crval1 = catalog[0].average[i].R; - coords.crval2 = catalog[0].average[i].D; - double Xi, Yi; - RD_to_XY (&Xi, &Yi, Ri, Di, &coords); - double radius = hypot(Xi, Yi); - if (radius > ExcludeBogusRadius) { - NskipBogus ++; - if (VERBOSE2) { - FILE *foutput = fbogus ? fbogus : stderr; - char *date = ohana_sec_to_date(catalog[0].measure[offset].t); - fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date); - free (date); - } - continue; + double Ri = getMeanR_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); + double Di = getMeanD_Big (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]); + coords.crval1 = catalog[0].average[i].R; + coords.crval2 = catalog[0].average[i].D; + double Xi, Yi; + RD_to_XY (&Xi, &Yi, Ri, Di, &coords); + double radius = hypot(Xi, Yi); + if (radius > ExcludeBogusRadius) { + NskipBogus ++; + if (VERBOSE2) { + FILE *foutput = fbogus ? fbogus : stderr; + char *date = ohana_sec_to_date(catalog[0].measure[offset].t); + fprintf (foutput, "exclude bogus: %10.6f %10.6f : %10.6f %10.6f : %6.2f %6.2f : %6.2f : %5d %s\n", catalog[0].average[i].R, catalog[0].average[i].D, Ri, Di, Xi, Yi, radius, catalog[0].measure[offset].photcode, date); + free (date); } + continue; + } } @@ -157,6 +188,6 @@ subcatalog[0].measureT[Nmeasure].averef = Naverage; if (RESET) { - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; - subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM; + subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA; } @@ -164,7 +195,12 @@ Nm ++; if (Nmeasure == NMEASURE) { - NMEASURE += 1000; - REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); - } + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE); + } + } + // we have filtered this object. do not keep it if we now have too few meas + if (Nm <= SRC_MEAS_TOOFEW) { + Nmeasure -= Nm; + continue; } subcatalog[0].average[Naverage].Nmeasure = Nm; @@ -186,5 +222,5 @@ assert (Nsecfilt == catalog[0].Nsecfilt); - // limit the total number of stars in the catalog +// limit the total number of stars in the catalog if (MaxDensityUse) { LimitDensityCatalog_ByNmeasureGrid (subcatalog, catalog); @@ -217,7 +253,7 @@ # define SWAPFUNC(A,B){ int tmpI; off_t tmpT; \ - tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \ - tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \ -} + tmpI = X[A]; X[A] = X[B]; X[B] = tmpI; \ + tmpT = Y[A]; Y[A] = Y[B]; Y[B] = tmpT; \ + } # define COMPARE(A,B)(X[A] > X[B]) @@ -231,5 +267,5 @@ /* generate a grid in a locally-projected coordinate system, try to select average entries from each grid cell in decending Nmeasure order. - */ +*/ int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog) { @@ -269,13 +305,8 @@ /* project coordinates to a plane centered on the object with units of arcsec */ Coords coords; + InitCoords (&coords, "DEC--SIN"); coords.crval1 = Rc; coords.crval2 = Dc; - coords.crpix1 = 0; - coords.crpix2 = 0; coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; - coords.pc1_1 = coords.pc2_2 = 1.0; - coords.pc1_2 = coords.pc2_1 = 0.0; - coords.Npolyterms = 1; - strcpy (coords.ctype, "DEC--SIN"); // convert all average R,D values to X,Y: Index: trunk/Ohana/src/relastro/src/client_logger.c =================================================================== --- trunk/Ohana/src/relastro/src/client_logger.c (revision 37807) +++ trunk/Ohana/src/relastro/src/client_logger.c (revision 37807) @@ -0,0 +1,37 @@ +# include "relastro.h" + +// I'm getting unlogged errors and failures. I need a log ouput that the clients can +// write independent of the master + +static FILE *logfile = NULL; +int client_logger_init (char *dirname) { + + char filename[DVO_MAX_PATH]; + + snprintf (filename, DVO_MAX_PATH, "%s/log.rlas.XXXXXX", dirname); + + int fd = mkstemp (filename); + if (fd == -1) { + fprintf (stderr, "failed to open client logger %s, exiting\n", filename); + exit (50); + } + + logfile = fdopen (fd, "w"); + if (!logfile) { + fprintf (stderr, "failed to fdopen client logger, exiting\n"); + exit (51); + } + return TRUE; +} + +int client_logger_message (char *format,...) { + + va_list argp; + + va_start (argp, format); + vfprintf (logfile, format, argp); + va_end (argp); + + fflush (logfile); + return TRUE; +} Index: trunk/Ohana/src/relastro/src/fitpoly.c =================================================================== --- trunk/Ohana/src/relastro/src/fitpoly.c (revision 37261) +++ trunk/Ohana/src/relastro/src/fitpoly.c (revision 37807) @@ -78,4 +78,39 @@ } +int array_print (double **matrix, int Nx, int Ny) { + + int i, j; + for (i = 0; i < Nx; i++) { + for (j = 0; j < Ny; j++) { + fprintf (stderr, "%10.3e ", matrix[i][j]); + } + fprintf (stderr, "\n"); + } + return TRUE; +} + +int fit_print (CoordFit *fit) { + + fprintf (stderr, "Npts: %d, Norder: %d, Nterms: %d, Nsums: %d, Nelems: %d\n", + fit[0].Npts, fit[0].Norder, fit[0].Nterms, fit[0].Nsums, fit[0].Nelems); + + fprintf (stderr, "sums: \n"); + array_print (fit[0].sum, fit[0].Nsums, fit[0].Nsums); + + fprintf (stderr, "xsums: \n"); + array_print (fit[0].xsum, fit[0].Nterms, fit[0].Nterms); + + fprintf (stderr, "ysums: \n"); + array_print (fit[0].ysum, fit[0].Nterms, fit[0].Nterms); + + fprintf (stderr, "xfits: \n"); + array_print (fit[0].xfit, fit[0].Nterms, fit[0].Nterms); + + fprintf (stderr, "yfits: \n"); + array_print (fit[0].yfit, fit[0].Nterms, fit[0].Nterms); + + return TRUE; +} + // XXX use implicit masks as below or explicit masks (with function to set?) // XXX eg, add a global mask to this file and @@ -134,4 +169,5 @@ for (j = 0; j < fit[0].Nelems; j++) { matrix[i][j] = 0.0; + matrix[j][i] = 0.0; } vector[i][0] = 0.0; @@ -148,5 +184,5 @@ } - if (!dgaussjordan (matrix, vector, fit[0].Nelems, 2)) { + if (!dgaussjordan_pivot (matrix, vector, fit[0].Nelems, 2, 1e-16)) { array_free (matrix, fit[0].Nelems); array_free (vector, fit[0].Nelems); @@ -277,4 +313,5 @@ /* this should only apply to the polynomial, not the projection terms */ /* compare with psastro supporting code */ +// this code will work for linear (Npolyterm == 0 or 1) and linear + map fits int fit_apply_coords (CoordFit *fit, Coords *coords, int keepRef) { @@ -352,4 +389,17 @@ /* keep the order and type from initial values */ + if (isnan(coords[0].crval1)) { + return FALSE; + } + if (isnan(coords[0].crval2)) { + return FALSE; + } + if (isnan(coords[0].crpix1)) { + return FALSE; + } + if (isnan(coords[0].crpix2)) { + return FALSE; + } + return (TRUE); } Index: trunk/Ohana/src/relastro/src/high_speed_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 37261) +++ trunk/Ohana/src/relastro/src/high_speed_objects.c (revision 37807) @@ -149,4 +149,5 @@ // define a local projection + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -155,11 +156,5 @@ tcoords.crval2 = 90.0; } - tcoords.crpix1 = 0; - tcoords.crpix2 = 0; tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - tcoords.Npolyterms = 1; - strcpy (tcoords.ctype, "DEC--ARC"); /* build spatial index (RA sort) referencing input array sequence */ @@ -271,5 +266,5 @@ populate_tiny_values (&catalogOut, DVO_TV_MEASURE); - UpdateObjects (&catalogOut, 1); + UpdateObjects (&catalogOut, 1, 0); free_tiny_values(&catalogOut); Index: trunk/Ohana/src/relastro/src/high_speed_utils.c =================================================================== --- trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 37261) +++ trunk/Ohana/src/relastro/src/high_speed_utils.c (revision 37807) @@ -141,5 +141,5 @@ off_t n; for (n = 0; n < NfieldsA; n++) { - valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsA[n]); + valuesA[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsA[n]); } return dbBooleanCond(stackA, NstackA, valuesA); @@ -155,5 +155,5 @@ off_t n; for (n = 0; n < NfieldsB; n++) { - valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], &fieldsB[n]); + valuesB[n] = dbExtractAverages (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], &catalog[0].measure[m], NULL, NULL, &fieldsB[n]); } return dbBooleanCond(stackB, NstackB, valuesB); Index: trunk/Ohana/src/relastro/src/hpm_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/hpm_objects.c (revision 37261) +++ trunk/Ohana/src/relastro/src/hpm_objects.c (revision 37807) @@ -126,4 +126,5 @@ // define a local projection + InitCoords (&tcoords, "DEC--ARC"); tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax); if (region[0].Dmax < 90) { @@ -132,10 +133,5 @@ tcoords.crval2 = 90.0; } - tcoords.crpix1 = 0; - tcoords.crpix2 = 0; tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - tcoords.Npolyterms = 1; strcpy (tcoords.ctype, "DEC--ARC"); @@ -256,5 +252,5 @@ // distinguishing a bad fit from a fit with 1 or 2 bad points FIT_MODE = FIT_PM_ONLY; - UpdateObjects (&testcat, 1); + UpdateObjects (&testcat, 1, 0); free_tiny_values(&testcat); @@ -286,5 +282,5 @@ populate_tiny_values (&catalogOut, DVO_TV_MEASURE); - UpdateObjects (&catalogOut, 1); + UpdateObjects (&catalogOut, 1, 0); free_tiny_values(&catalogOut); Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 37261) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 37807) @@ -9,5 +9,5 @@ if (DCR_BLUE_COLOR_POS) fprintf (stderr, "DCR_BLUE_COLOR_POS: %s - %s\n", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); - if (DCR_BLUE_COLOR_POS) fprintf (stderr, "DCR_BLUE_COLOR_POS: %s - %s\n", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + if (DCR_RED_COLOR_POS) fprintf (stderr, "DCR_RED_COLOR_POS: %s - %s\n", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); if (PHOTCODE_KEEP_LIST) fprintf (stderr, "PHOTCODE_KEEP_LIST: %s\n", PHOTCODE_KEEP_LIST); @@ -63,4 +63,12 @@ initstats (STATMODE); + if (USE_ICRF_CORRECT) { + if (!USE_ICRF_LOCAL && !USE_ICRF_SHFIT) { + fprintf (stderr, "no ICRF correction method chosen\n"); + exit (2); + } + } + if (USE_ICRF_CORRECT) ICRFinit(); + /* XXX drop irrelevant entries */ if (SHOW_PARAMS) { Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 37261) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 37807) @@ -6,25 +6,52 @@ int i; - // clear the I/O files - for (i = 0; i < regionHosts->Nhosts; i++) { - char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.sync"); - if (truncate (syncfile, 0)) fprintf (stderr, "trouble clearing syncfile %s\n", syncfile); - free (syncfile); + // do not remote the sync and fits files if we do a manual run -- user must clear if needed + if (!PARALLEL_REGIONS_MANUAL) { + // clear the I/O files + for (i = 0; i < regionHosts->Nhosts; i++) { + char *meansync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.sync"); + if (truncate (meansync, 0)) fprintf (stderr, "trouble clearing meansync %s\n", meansync); + free (meansync); - char *fitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanmags.fits"); - if (truncate (fitsfile, 0)) fprintf (stderr, "trouble clearing fitsfile %s\n", fitsfile); - free (fitsfile); + char *meanfits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "meanpos.fits"); + if (truncate (meanfits, 0)) fprintf (stderr, "trouble clearing meanfits %s\n", meanfits); + free (meanfits); - char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.sync"); - if (truncate (imsyncfile, 0)) fprintf (stderr, "trouble clearing imsyncfile %s\n", imsyncfile); - free (imsyncfile); + char *meassync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.sync"); + if (truncate (meassync, 0)) fprintf (stderr, "trouble clearing meassync %s\n", meassync); + free (meassync); - char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagemags.fits"); - if (truncate (imfitsfile, 0)) fprintf (stderr, "trouble clearing imfitsfile %s\n", imfitsfile); - free (imfitsfile); + char *measfits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.fits"); + if (truncate (measfits, 0)) fprintf (stderr, "trouble clearing measfits %s\n", measfits); + free (measfits); - char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); - if (truncate (loopsyncfile, 0)) fprintf (stderr, "trouble clearing loopsyncfile %s\n", loopsyncfile); - free (loopsyncfile); + char *icrfsync = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.sync"); + if (truncate (icrfsync, 0)) fprintf (stderr, "trouble clearing icrfsync %s\n", icrfsync); + free (icrfsync); + + char *icrffits = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.fits"); + if (truncate (icrffits, 0)) fprintf (stderr, "trouble clearing icrffits %s\n", icrffits); + free (icrffits); + + char *imsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.sync"); + if (truncate (imsyncfile, 0)) fprintf (stderr, "trouble clearing imsyncfile %s\n", imsyncfile); + free (imsyncfile); + + char *imfitsfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "imagepos.fits"); + if (truncate (imfitsfile, 0)) fprintf (stderr, "trouble clearing imfitsfile %s\n", imfitsfile); + free (imfitsfile); + + char *loopsyncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "loop.sync"); + if (truncate (loopsyncfile, 0)) fprintf (stderr, "trouble clearing loopsyncfile %s\n", loopsyncfile); + free (loopsyncfile); + } + + char *framesync = make_filename (CATDIR, "master", 0, "frame.corr.sync"); + if (truncate (framesync, 0)) fprintf (stderr, "trouble clearing framesync %s\n", framesync); + free (framesync); + + char *framefits = make_filename (CATDIR, "master", 0, "frame.corr.fits"); + if (truncate (framefits, 0)) fprintf (stderr, "trouble clearing framefits %s\n", framefits); + free (framefits); } @@ -37,7 +64,16 @@ char filename[1024]; snprintf (filename, 1024, "%s/Image.%d.fits", CATDIR, host->hostID); + if (unlink (filename)) fprintf (stderr, "trouble clearing image %s\n", filename); // write the image subset for this host ImageTableSave (filename, host->image, host->Nimage); + + if (host->astromTable) { + char mapname[1024]; + snprintf (mapname, 1024, "%s/AstroMap.%d.fits", CATDIR, host->hostID); + + // write the image subset for this host + AstromOffsetMapSave (host->astromTable, mapname); + } char command[1024]; @@ -92,4 +128,7 @@ if (PhotFlagPoor) strextend (command, "+photflagpoor %d", PhotFlagPoor); + if (CHIPMAP) strextend (command, "-chipmap %d", CHIPMAP); + if (RESET_IMAGES) strextend (command, "-reset-images"); + if (MinBadQF > 0.0) strextend (command, "-min-bad-psfqf %f", MinBadQF); if (MaxMeanOffset != 10.0) strextend (command, "-max-mean-offset %f", MaxMeanOffset); @@ -99,4 +138,11 @@ strextend (command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX); strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS); + + strextend (command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); + + strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); + strextend (command, "-D USE_ICRF_LOCAL %d", USE_ICRF_LOCAL); + strextend (command, "-D USE_ICRF_SHFIT %d", USE_ICRF_SHFIT); + strextend (command, "-D USE_ICRF_POLE %d", USE_ICRF_POLE); if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1"); @@ -133,16 +179,20 @@ } + if (USE_ICRF_CORRECT) return TRUE; + + int status = HarvestRegionHosts (regionHosts); + return status; +} + +int HarvestRegionHosts (RegionHostTable *regionHosts) { + + int status = TRUE; + if (PARALLEL_REGIONS_MANUAL) { fprintf (stderr, "run the relastro_client commands above. when these are done, hit return\n"); getchar(); } else { - RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE); + status = RegionHostTableWaitJobsGetIO (regionHosts, __FILE__, __LINE__, VERBOSE); } - - int status = TRUE; - for (i = 0; i < regionHosts->Nhosts; i++) { - status = status && (regionHosts->hosts[i].status == 0); - } - return status; } Index: trunk/Ohana/src/relastro/src/liststats.c =================================================================== --- trunk/Ohana/src/relastro/src/liststats.c (revision 37261) +++ trunk/Ohana/src/relastro/src/liststats.c (revision 37807) @@ -252,5 +252,5 @@ X2 += M / dM; dS += M; - if (XVERB) fprintf (stderr, "chisq contrib : %f : %f : %f : %f : %f\n", Mo, 3600.0*(value[i] - Mo), 3600.0*dvalue[i], (value[i] - Mo) / dvalue[i], X2); + if (XVERB) fprintf (stderr, "chisq contrib : %f : %f : %f : %f : %f\n", Mo, (value[i] - Mo), dvalue[i], (value[i] - Mo) / dvalue[i], X2); } X2 = X2 / (Nm - 1); Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 37261) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 37807) @@ -6,5 +6,5 @@ Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath, char *syncfile) { - int i, j, Nstar; + int i, j; // int k, m; Catalog *catalog, *pcatalog, tcatalog; @@ -88,8 +88,11 @@ bcatalog_show_skips(); - - Nstar = 0; + fprintf (stderr, "included %d ICRF QSOs\n", ICRFmax()); + + int Nstar = 0; + int Nmeas = 0; for (i = 0; i < Ncat; i++) { Nstar += catalog[i].Naverage; + Nmeas += catalog[i].Nmeasure; if ((sizeof(IDX_T) == 8) && (catalog[i].Nmeasure > 0xffffffff)) { fprintf (stderr, "ERROR: using small-sized IDX_T on data with more than 2G detections per table will cause errors\n"); @@ -101,4 +104,5 @@ fprintf (stderr, "warning: insufficient stars %d\n", Nstar); } + fprintf (stderr, "using %d stars, %d measurements to calibrate\n", Nstar, Nmeas); // if we are running with parallel_images but not a parallel database, we need to @@ -170,4 +174,7 @@ if (FlagOutlier) strextend (command, "-clip %d", CLIP_THRESH); if (ExcludeBogus) strextend (command, "-exclude-bogus %f", ExcludeBogusRadius); + + if (USE_ICRF_CORRECT) strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); + if (USE_GALAXY_MODEL) strextend (command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 37261) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 37807) @@ -1,3 +1,9 @@ # include "relastro.h" + +static AstromOffsetTable *table = NULL; + +/* AstromOffsetMap is the map of Astrometry Offsets for non-polynomial astrometric + * corrections for each chip. + */ int load_images (FITS_DB *db, SkyList *skylist, int UseFullOverlap) { @@ -17,7 +23,27 @@ MARKTIME(" convert image table: %f sec\n", dtime); + // 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); + table = AstromOffsetMapLoad (mapfile, VERBOSE); + + // assign images.coords.offsetMap -> table->map[i] + if (table) { + AstromOffsetTableMatchChips (image, Nimage, table); + } else { + table = AstromOffsetTableInit (); + } + // select the images which overlap the selected sky regions subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset, UseFullOverlap); MARKTIME(" select images: %f sec\n", dtime); + + if (Nsubset == Nimage) { + free (subset); + subset = image; + } initImages (subset, LineNumber, Nsubset); @@ -51,2 +77,20 @@ return TRUE; } + +int save_astrom_table () { + + char mapfile[DVO_MAX_PATH]; + snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.fits", CATDIR); + AstromOffsetMapSave (table, mapfile); + + return TRUE; +} + +AstromOffsetTable *get_astrom_table () { + return table; +} + +void put_astrom_table (AstromOffsetTable *myTable) { + table = myTable; +} + Index: trunk/Ohana/src/relastro/src/markObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/markObjects.c (revision 37261) +++ trunk/Ohana/src/relastro/src/markObjects.c (revision 37807) @@ -14,6 +14,16 @@ int m = catalog[i].average[j].measureOffset; for (n = 0; n < catalog[i].average[j].Nmeasure; n++) { - if (!catalog[i].measureT[m+n].myDet) continue; - nOwn ++; + // this is a detection I own (it is on one of my images) + if (catalog[i].measureT[m+n].myDet) { + nOwn ++; + continue; + } + // 2MASS detections are counted as 'owned' (I completely control this object is nOwn_t == Nmeasure; no sharing needed) + if ((catalog[i].measureT[m+n].photcode == 2011) || + (catalog[i].measureT[m+n].photcode == 2012) || + (catalog[i].measureT[m+n].photcode == 2013)) { + nOwn ++; + continue; + } } catalog[i].nOwn_t[j] = nOwn; Index: trunk/Ohana/src/relastro/src/relastroVisual.c =================================================================== --- trunk/Ohana/src/relastro/src/relastroVisual.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastroVisual.c (revision 37807) @@ -24,4 +24,8 @@ /** Spawn a kapa window */ static int initWindow(int *kapid) { + if (NTHREADS) { + fprintf (stderr, "-visual and -threads are incompatible, quitting\n"); + exit (2); + } if (*kapid == -1) { *kapid = KapaOpenNamedSocket("kapa", "relastro_plots"); @@ -348,5 +352,5 @@ int i; - for(i = 0; i < numObj; i++) { + for (i = 0; i < numObj; i++) { if (raw[i].mask) continue; // XXX Index: trunk/Ohana/src/relastro/src/relastro_client.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_client.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastro_client.c (revision 37807) @@ -79,4 +79,6 @@ load_images (&db, skylist, FALSE); + // I can free the database after I have loaded the images... + // // load the image subset table from the specified location // off_t Nimage; Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 37807) @@ -23,9 +23,5 @@ MARKTIME("load images: %f sec\n", dtime); - initCoords(); - /* load catalog data from region files : subselect high-quality measurements */ - // XXX pass in the image table - // XXX who carries the image grid? // photcodesKeep is used here to allow measurements from the images being calibrated @@ -33,4 +29,7 @@ catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, NULL); MARKTIME("load catalog data: %f sec\n", dtime); + + // find ICRF QSOs for reference downstream (only if USE_ICRF_CORRECT) + select_catalog_ICRF (catalog, Ncatalog); if (photcodesReset) { @@ -65,5 +64,5 @@ case TARGET_SIMPLE: for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, i); UpdateSimple (catalog, Ncatalog); } @@ -71,7 +70,13 @@ case TARGET_CHIPS: + if (RESET_IMAGES) UpdateMeasures (catalog, Ncatalog); for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); - UpdateChips (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, (i > 0)); // calculate ,; if (i > 0), apply Galaxy Motion Model (if desired) + if ((i > 1) || !USE_GALAXY_MODEL) { + // if GALAXY_MODEL is selected, we want to delay the frame correction until we have + // applied the galaxy model a couple of times. + FrameCorrectionSerial (catalog, Ncatalog); + } + UpdateChips (catalog, Ncatalog); // measure.X,Y -> R,D, fit image.coords MARKTIME("update chips: %f sec\n", dtime); } @@ -82,5 +87,5 @@ case TARGET_MOSAICS: for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, i); UpdateMosaic (catalog, Ncatalog); } @@ -102,4 +107,6 @@ if (PARALLEL) { // save the updated image parameters + // need to also save the image map table... + save_astrom_table (); dvo_image_update (&db, VERBOSE); dvo_image_unlock (&db); @@ -109,9 +116,7 @@ UpdateObjectOffsets (skylist, 0, NULL); - // iterate over catalogs to make detection coordinates consistant - // FixProblemImages (skylist); - if (!PARALLEL) { // save the updated image parameters + save_astrom_table (); dvo_image_update (&db, VERBOSE); dvo_image_unlock (&db); Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 37807) @@ -77,5 +77,5 @@ populate_tiny_values(&catalog, DVO_TV_MEASURE); - UpdateObjects (&catalog, 1); + UpdateObjects (&catalog, 1, 0); free_tiny_values(&catalog); Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37807) @@ -17,8 +17,11 @@ INITTIME; + client_logger_init (CATDIR); + // load the RegionTable (UserRegion should not be used at this level) RegionHostTable *regionHosts = RegionHostTableLoad (CATDIR, REGION_FILE); int myHost = regionHosts->index[REGION_HOST_ID]; RegionHostFindNeighbors (regionHosts, myHost); + client_logger_message ("started parallel images on %s\n", regionHosts->hosts[myHost].hostname); // load the subset images belonging to this host @@ -29,4 +32,20 @@ exit (2); } + + // assign image->parent and image->coords.mosaic + BuildChipMatch (image, Nimage); + + char mapfile[DVO_MAX_PATH]; + snprintf (mapfile, DVO_MAX_PATH, "%s/AstroMap.%d.fits", CATDIR, REGION_HOST_ID); + AstromOffsetTable *table = AstromOffsetMapLoad (mapfile, VERBOSE); + + // assign images.coords.offsetMap -> table->map[i] + if (table) { + AstromOffsetTableMatchChips (image, Nimage, table); + } else { + table = AstromOffsetTableInit (); + } + put_astrom_table (table); + client_logger_message ("loaded images on %s\n", regionHosts->hosts[myHost].hostname); // once we have read this table, we should remove it for repeat runs @@ -42,12 +61,15 @@ SkyList *skylist = SkyListByBounds (sky, -1, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); - initCoords(); - /* load catalog data from region files (hostID is 0 since we are not a client */ char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync"); catalog = load_catalogs (skylist, &Ncatalog, TRUE, 0, NULL, syncfile); - MARKTIME("-- load catalog data: %f sec\n", dtime); + MARKTIME("-- load catalog data, host %d: %f sec\n", REGION_HOST_ID, dtime); free (syncfile); + client_logger_message ("loaded catalog data %s\n", regionHosts->hosts[myHost].hostname); + + // find ICRF QSOs for reference downstream (only if USE_ICRF_CORRECT) + select_catalog_ICRF (catalog, Ncatalog); + if (photcodesReset) { photcodesKeep = photcodesReset; @@ -60,11 +82,13 @@ /* match measurements with images, mosaics */ initImageBins (catalog, Ncatalog, TRUE); - MARKTIME("-- make image bins: %f sec\n", dtime); + MARKTIME("-- make image bins, host %d: %f sec\n", REGION_HOST_ID, dtime); findImages (catalog, Ncatalog, TRUE); - MARKTIME("-- set up image indexes: %f sec\n", dtime); + MARKTIME("-- set up image indexes, host %d: %f sec\n", REGION_HOST_ID, dtime); // set test points based on the starmap createStarMap (catalog, Ncatalog); + + client_logger_message ("set up image indexes on %s\n", regionHosts->hosts[myHost].hostname); markObjects (catalog, Ncatalog); @@ -72,26 +96,52 @@ SAVEPLOT = FALSE; + client_logger_message ("starting the loops: %s\n", regionHosts->hosts[myHost].hostname); + /* major modes */ switch (FIT_TARGET) { case TARGET_SIMPLE: for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, i); share_mean_pos (catalog, Ncatalog, regionHosts, i); slurp_mean_pos (catalog, Ncatalog, regionHosts, i); UpdateSimple (catalog, Ncatalog); share_image_pos (regionHosts, i); - slurp_image_pos (regionHosts, i); + slurp_image_pos (catalog, Ncatalog, regionHosts, i); } break; case TARGET_CHIPS: + if (RESET_IMAGES) { + UpdateMeasures (catalog, Ncatalog); + MARKTIME("reset measures, host %d: %f sec\n", REGION_HOST_ID, dtime); + client_logger_message ("reset measures: %s\n", regionHosts->hosts[myHost].hostname); + } for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, (i > 0)); + MARKTIME("update objects loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime); + client_logger_message ("update objects loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); + if ((i > 1) || !USE_GALAXY_MODEL) { + // if GALAXY_MODEL is selected, we want to delay the frame correction until we have + // applied the galaxy model a couple of times. + FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i); + MARKTIME("frame correction loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime); + client_logger_message ("frame correction loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); + } share_mean_pos (catalog, Ncatalog, regionHosts, i); slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + client_logger_message ("slurp mags loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); + UpdateChips (catalog, Ncatalog); + MARKTIME("update chips loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime); + client_logger_message ("update chips loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); + share_image_pos (regionHosts, i); - slurp_image_pos (regionHosts, i); - MARKTIME("update chips: %f sec\n", dtime); + slurp_image_pos (catalog, Ncatalog, regionHosts, i); + client_logger_message ("slurp image pos loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); + + share_meas_pos (catalog, Ncatalog, regionHosts, i); + slurp_meas_pos (catalog, Ncatalog, regionHosts, i); + MARKTIME("done exchange loop %d, host %d: %f sec\n", i, REGION_HOST_ID, dtime); + client_logger_message ("slurp meas loop %d: %s\n", i, regionHosts->hosts[myHost].hostname); } // create summary plots of the process @@ -101,10 +151,10 @@ case TARGET_MOSAICS: for (i = 0; i < NLOOP; i++) { - UpdateObjects (catalog, Ncatalog); + UpdateObjects (catalog, Ncatalog, i); share_mean_pos (catalog, Ncatalog, regionHosts, i); slurp_mean_pos (catalog, Ncatalog, regionHosts, i); UpdateMosaic (catalog, Ncatalog); share_image_pos (regionHosts, i); - slurp_image_pos (regionHosts, i); + slurp_image_pos (catalog, Ncatalog, regionHosts, i); } break; @@ -128,4 +178,5 @@ share_image_pos (regionHosts, -1); + client_logger_message ("share image pos loop %d: %s\n", -1, regionHosts->hosts[myHost].hostname); exit (0); Index: trunk/Ohana/src/relastro/src/relastro_parallel_regions.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37261) +++ trunk/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37807) @@ -44,6 +44,10 @@ if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts"); + // If we are doing FrameCorrection, we have to loop here waiting for ICRF objects from the region + // hosts. we then measure the frame correction and send back the result + if (!FrameCorrectionParallelMaster (regionHosts)) Shutdown ("error running parallel frame correction"); + // retrieve updated image parameters from the remote hosts (also set Image.mcal) - if (!slurp_image_pos (regionHosts, -1)) Shutdown ("error loading image updates"); + if (!slurp_image_pos (NULL, 0, regionHosts, -1)) Shutdown ("error loading image updates"); if (!UPDATE) { @@ -59,5 +63,7 @@ if (PARALLEL) { - // save the updated image parameters + // save the updated image parameters. if we are NOT PARALLEL, we have to wait + // until after UpdateObjectOffsets, or the image parameters are swapped and/or freed + save_astrom_table (); dvo_image_save (&db, VERBOSE); dvo_image_unlock (&db); @@ -65,8 +71,11 @@ // iterate over catalogs to make detection coordinates consistant - UpdateObjectOffsets (skylist, 0, NULL); + if (APPLY_OFFSETS) { + UpdateObjectOffsets (skylist, 0, NULL); + } if (!PARALLEL) { // save the changes to the image parameters + save_astrom_table (); dvo_image_save (&db, VERBOSE); dvo_image_unlock (&db); Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 37261) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 37807) @@ -23,5 +23,4 @@ int InRange, found; double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; - Coords tcoords; SkyRegionCoords *skycoords; @@ -46,9 +45,7 @@ // the comparison is made in the catalog local projection. below we set crval1,2 - tcoords.crpix1 = tcoords.crpix2 = 0.0; + Coords tcoords; + InitCoords (&tcoords, "DEC--TAN"); tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; - tcoords.pc1_1 = tcoords.pc2_2 = 1.0; - tcoords.pc1_2 = tcoords.pc2_1 = 0.0; - strcpy (tcoords.ctype, "DEC--TAN"); ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); @@ -104,6 +101,4 @@ if (VERBOSE) fprintf (stderr, "finding images\n"); - BuildChipMatch (timage, Ntimage); - MARKTIME("build chip match: %f sec\n", dtime); nimage = 0; @@ -148,10 +143,4 @@ } - // this adds 1.3 sec for 3M images - if (!FindMosaicForImage (timage, Ntimage, i)) { - fprintf (stderr, "cannot find mosaic for "OFF_T_FMT"\n", i); - continue; - } - /* define image corners - note the DIS images (mosaic phu) are special */ if (!strcmp(&timage[i].coords.ctype[4], "-DIS")) { Index: trunk/Ohana/src/relastro/src/share_icrf_obj.c =================================================================== --- trunk/Ohana/src/relastro/src/share_icrf_obj.c (revision 37807) +++ trunk/Ohana/src/relastro/src/share_icrf_obj.c (revision 37807) @@ -0,0 +1,136 @@ +# include "relastro.h" + +// we are sharing mean and specific positions for all known ICRF objects which I own +// this function occurs on the remote hosts in parallel-region processing +int share_icrf_obj (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + if (!USE_ICRF_CORRECT) return TRUE; + + ICRFobj *icrfobj = get_ICRF_data (catalog, Ncatalog); + + // write out the meanmag fits table AND write state in some file + int myHost = regionHosts->index[REGION_HOST_ID]; + char *hostname = regionHosts->hosts[myHost].hostname; + + char *filename = make_filename (CATDIR, hostname, REGION_HOST_ID, "icrfobj.fits"); + ICRFobjSave (filename, icrfobj); + free (filename); + + free (icrfobj); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "icrfobj.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +// load mean and specific positions for all known ICRF objects from remote hosts. +// this function occurs on the master host in parallel-region processing +ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop) { + + off_t i; + + ICRFobj *icrfobj = NULL; + ALLOCATE (icrfobj, ICRFobj, 1); + ALLOCATE (icrfobj->Rave, double, 1); + ALLOCATE (icrfobj->Dave, double, 1); + ALLOCATE (icrfobj->dRoff, double, 1); + ALLOCATE (icrfobj->dDoff, double, 1); + icrfobj->Nicrfobj = 0; + + fprintf (stderr, "grabbing icrf objects from remote hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + char *filename = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.fits"); + ICRFobj *icrfobjSubset = ICRFobjLoad (filename); + free (filename); + + // merge_mean_pos reallocs icrfobj and frees the input icrfobjSubset + icrfobj = merge_icrf_obj (icrfobj, icrfobjSubset); + } + return icrfobj; +} + +ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source) { + + off_t i; + + REALLOCATE (target->Rave, double, target->Nicrfobj + source->Nicrfobj); + REALLOCATE (target->Dave, double, target->Nicrfobj + source->Nicrfobj); + REALLOCATE (target->dRoff, double, target->Nicrfobj + source->Nicrfobj); + REALLOCATE (target->dDoff, double, target->Nicrfobj + source->Nicrfobj); + + for (i = 0; i < source->Nicrfobj; i++) { + off_t n = i + target->Nicrfobj; + target->Rave[n] = source->Rave[i]; + target->Dave[n] = source->Dave[i]; + target->dRoff[n] = source->dRoff[i]; + target->dDoff[n] = source->dDoff[i]; + } + target->Nicrfobj += source->Nicrfobj; + + ICRFobjFree (source); + return (target); +} + +ICRFobj *get_ICRF_data (Catalog *catalog, int Ncatalog) { + + int i; + int Nicrf = ICRFmax(); + + // we need to write an empty table none are available + ICRFobj *icrfobj = NULL; + ALLOCATE (icrfobj, ICRFobj, 1); + + ALLOCATE (icrfobj->Rave, double, Nicrf); + ALLOCATE (icrfobj->Dave, double, Nicrf); + ALLOCATE (icrfobj->dRoff, double, Nicrf); + ALLOCATE (icrfobj->dDoff, double, Nicrf); + + // select the ICRF QSOS and save the necessary data + int Npts = 0; + for (i = 0; i < Nicrf; i++) { + + int cat, ave, meas; + ICRFdata (i, &cat, &ave, &meas); + + Average *average = &catalog[cat].average[ave]; + MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny? + + // record these in arcsec or degree? + // correct for cos(D) or not? + double dR = 3600.0*(average->R - measure->R)*cos(average->R*RAD_DEG); + double dD = 3600.0*(average->D - measure->D); + + // dR = 3600*(Rave - Rmeas) -> Rmeas = Rave - dR / 36000 + + if (isnan(dR)) continue; + if (isnan(dD)) continue; + + icrfobj->Rave[Npts] = average->R; + icrfobj->Dave[Npts] = average->D; + icrfobj->dRoff[Npts] = dR; + icrfobj->dDoff[Npts] = dD; + Npts ++; + } + icrfobj->Nicrfobj = Npts; + + return icrfobj; +} + +void ICRFobjFree (ICRFobj *icrfobj) { + + if (!icrfobj) return; + + FREE (icrfobj->Rave); + FREE (icrfobj->Dave); + FREE (icrfobj->dRoff); + FREE (icrfobj->dDoff); + return; +} + Index: trunk/Ohana/src/relastro/src/share_images_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_images_pos.c (revision 37261) +++ trunk/Ohana/src/relastro/src/share_images_pos.c (revision 37807) @@ -34,4 +34,14 @@ free (iposfile); + // save the per-chip residual table + AstromOffsetTable *table = get_astrom_table (); + if (table) { + char mapname[1024]; + snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, REGION_HOST_ID); + + // write the image subset for this host + AstromOffsetMapSave (table, mapname); + } + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "imagepos.sync"); update_sync_file (syncfile, nloop); @@ -41,7 +51,7 @@ } -int slurp_image_pos (RegionHostTable *regionHosts, int nloop) { +int slurp_image_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { - off_t Nimage, i; + off_t Nimage, i, j; Image *images = getimages (&Nimage, NULL); @@ -74,10 +84,55 @@ continue; } + Coords *moscoords = images[seq].coords.mosaic; images[seq].coords = image_pos[i].coords ; - images[seq].dXpixSys = image_pos[i].dXpixSys ; - images[seq].dYpixSys = image_pos[i].dYpixSys ; - images[seq].imageID = image_pos[i].imageID ; - images[seq].nFitAstrom = image_pos[i].nFitAstrom; - images[seq].flags = image_pos[i].flags ; + images[seq].coords.mosaic = moscoords; + + images[seq].dXpixSys = image_pos[i].dXpixSys ; + images[seq].dYpixSys = image_pos[i].dYpixSys ; + images[seq].refColorBlue = image_pos[i].refColorBlue; + images[seq].refColorRed = image_pos[i].refColorRed; + images[seq].imageID = image_pos[i].imageID ; + images[seq].nFitAstrom = image_pos[i].nFitAstrom; + images[seq].flags = image_pos[i].flags ; + } + + // load the astrometry offset maps, if they exist, and apply + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char mapname[1024]; + snprintf (mapname, 1024, "%s/AstroMapUpdate.%d.fits", CATDIR, regionHosts->hosts[i].hostID); + + AstromOffsetTable *table = AstromOffsetMapLoad (mapname, VERBOSE); + + // apply table entries here to existing images + for (j = 0; j < table->Nmap; j++) { + off_t seq = getImageByID (table->map[j][0].imageID); + // I do not necessarily own all images listed in the table. skip if not found + if (seq < 0) continue; + + AstromOffsetMap *oldMap = images[seq].coords.offsetMap; + AstromOffsetMap *newMap = table->map[j]; + + if (oldMap) { + AstromOffsetMapSetOrder (oldMap, newMap->Nx, newMap->Ny, &images[seq]); + AstromOffsetMapCopyData (oldMap, newMap); + } else { + lockUpdateChips (); + AstromOffsetTable *FullTable = get_astrom_table (); + AstromOffsetTableNewMap(FullTable, newMap->Nx, newMap->Ny, &images[seq]); // registers the map with the image + unlockUpdateChips (); + AstromOffsetMapCopyData (images[seq].coords.offsetMap, newMap); // oldMap is now not valid + } + } + AstromOffsetTableFree(table); + } + + // apply the modified image parameters to my detections + for (i = 0; FALSE && (Ncatalog > 0) && (i < Nimage_pos); i++) { + off_t seq = getImageByID (image_pos[i].imageID); + if (seq < 0) continue; + updateImageRaw (catalog, Ncatalog, seq); } free (image_pos); @@ -93,4 +148,6 @@ image_pos->dXpixSys = image->dXpixSys; image_pos->dYpixSys = image->dYpixSys; + image_pos->refColorBlue = image->refColorBlue; + image_pos->refColorRed = image->refColorRed; image_pos->imageID = image->imageID; image_pos->nFitAstrom = image->nFitAstrom; Index: trunk/Ohana/src/relastro/src/share_meas_pos.c =================================================================== --- trunk/Ohana/src/relastro/src/share_meas_pos.c (revision 37807) +++ trunk/Ohana/src/relastro/src/share_meas_pos.c (revision 37807) @@ -0,0 +1,133 @@ +# include "relastro.h" +// we are sharing meas positions for all objects which (a) I own and (b) which have unowned detections + +# define D_NMEASPOS 10000 +int share_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + int i; + off_t j, k; + + off_t Nmeaspos = 0; + off_t NMEASPOS = D_NMEASPOS; + + MeasPos *measpos = NULL; + ALLOCATE (measpos, MeasPos, NMEASPOS); + + int myHost = regionHosts->index[REGION_HOST_ID]; + + // XXX skip some catalogs based on UserPatch? + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + // does this object have missing detections (does someone else need it?) + if (catalog[i].nOwn_t[j] == catalog[i].average[j].Nmeasure) continue; + + off_t m = catalog[i].average[j].measureOffset; + for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) { + if (!catalog[i].measureT[m].myDet) continue; + + set_meas_pos (&measpos[Nmeaspos], &catalog[i].average[j], &catalog[i].measureT[m]); + Nmeaspos ++; + CHECK_REALLOCATE (measpos, MeasPos, NMEASPOS, Nmeaspos, D_NMEASPOS); + } + } + } + + // write out the measmag fits table AND write state in some file + char *hostname = regionHosts->hosts[myHost].hostname; + + char *posfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "measpos.fits"); + MeasPosSave (posfile, measpos, Nmeaspos); + free (measpos); + free (posfile); + + char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "measpos.sync"); + update_sync_file (syncfile, nloop); + free (syncfile); + + return TRUE; +} + +int slurp_meas_pos (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) { + + off_t i, k; + + int Nmeaspos = 0; + MeasPos *measpos = NULL; + ALLOCATE (measpos, MeasPos, 1); + + fprintf (stderr, "grabbing meas object pos from other hosts...\n"); + + for (i = 0; i < regionHosts->Nhosts; i++) { + if (regionHosts->hosts[i].hostID == REGION_HOST_ID) continue; + if (REGION_HOST_ID && !regionHosts->hosts[i].isNeighbor) continue; + + char *syncfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.sync"); + check_sync_file (syncfile, nloop); + free (syncfile); + + off_t Nsubset = 0; + char *posfile = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "measpos.fits"); + MeasPos *measposSubset = MeasPosLoad (posfile, &Nsubset); + free (posfile); + + // merge_meas_pos reallocs measpos and frees the input measposSubset + measpos = merge_meas_pos (measpos, &Nmeaspos, measposSubset, Nsubset); + } + + for (i = 0; i < Nmeaspos; i++) { + int objID = measpos[i].objID; + int catID = measpos[i].catID; + + // set the meas mag + int catSeq; + off_t objSeq; + if (!catID_and_objID_to_seq (catID, objID, &catSeq, &objSeq)) { + // XXX what should I do if this does not match? + continue; + } + off_t m = catalog[catSeq].average[objSeq].measureOffset; + for (k = 0; k < catalog[catSeq].average[objSeq].Nmeasure; k++, m++) { + if (catalog[catSeq].measureT[m].imageID != measpos[i].imageID) continue; + catalog[catSeq].measureT[m].R = measpos[i].R; + catalog[catSeq].measureT[m].D = measpos[i].D; + if (catalog[catSeq].measure) { + catalog[catSeq].measure[m].R = measpos[i].R; + catalog[catSeq].measure[m].D = measpos[i].D; + } + } + } + free (measpos); + + fprintf (stderr, "DONE grabbing meas object pos from other hosts...\n"); + + return TRUE; +} + +int set_meas_pos (MeasPos *measpos, Average *average, MeasureTiny *measure) { + + measpos->R = measure->R; + measpos->D = measure->D; + measpos->objID = average->objID; + measpos->catID = average->catID; + measpos->imageID = measure->imageID; + + return TRUE; +} + +MeasPos *merge_meas_pos (MeasPos *target, int *ntarget, MeasPos *source, int Nsource) { + + off_t i; + + REALLOCATE (target, MeasPos, *ntarget + Nsource); + for (i = 0; i < Nsource; i++) { + off_t n = i + *ntarget; + target[n] = source[i]; + } + + free (source); + + *ntarget += Nsource; + return (target); +} + Index: trunk/Ohana/src/relastro/src/testparallax.c =================================================================== --- trunk/Ohana/src/relastro/src/testparallax.c (revision 37261) +++ trunk/Ohana/src/relastro/src/testparallax.c (revision 37807) @@ -49,13 +49,6 @@ /* project coordinates to a plane centered on the object with units of arcsec */ - coords.crval1 = 0; - coords.crval2 = 0; - coords.crpix1 = 0; - coords.crpix2 = 0; + InitCoords (&coords, "DEC--SIN"); coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; - coords.pc1_1 = coords.pc2_2 = 1.0; - coords.pc1_2 = coords.pc2_1 = 0.0; - coords.Npolyterms = 1; - strcpy (coords.ctype, "DEC--SIN"); // use one point as a local reference