Index: /branches/eam_branches/ipp-20140904/Ohana/src/libfits/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libfits/Makefile (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libfits/Makefile (revision 37658) @@ -32,5 +32,5 @@ $(HEAD)/F_scan.$(ARCH).o $(HEAD)/F_print.$(ARCH).o \ $(HEAD)/F_copy_H.$(ARCH).o $(HEAD)/F_delete.$(ARCH).o \ -$(HEAD)/F_find_XH.$(ARCH).o $(HEAD)/F_init_H.$(ARCH).o \ +$(HEAD)/F_find_H.$(ARCH).o $(HEAD)/F_init_H.$(ARCH).o \ $(HEAD)/F_convert_H.$(ARCH).o $(HEAD)/F_read_XH.$(ARCH).o \ $(HEAD)/version.$(ARCH).o Index: /branches/eam_branches/ipp-20140904/Ohana/src/libfits/header/F_find_H.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libfits/header/F_find_H.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libfits/header/F_find_H.c (revision 37658) @@ -5,7 +5,4 @@ int gfits_find_Xheader (FILE *f, Header *header, char *extname) { - int j, found; - off_t Nbytes; - Header *header; char tname[80]; Index: /branches/eam_branches/ipp-20140904/Ohana/src/libfits/include/gfitsio.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libfits/include/gfitsio.h (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libfits/include/gfitsio.h (revision 37658) @@ -119,5 +119,5 @@ int gfits_modify_alt PROTO((Header *header, char *field, char *mode, int N,...)); int gfits_print_alt PROTO((Header *header, char *field, char *mode, int N,...)); -int gfits_find_Xheader PROTO((char *filename, Header *header, char *extname)); +int gfits_find_Xheader PROTO((FILE *f, Header *header, char *extname)); int gfits_read_Xheader PROTO((char *filename, Header *header, int N)); int gfits_read_header PROTO((char *filename, Header *header)); Index: /branches/eam_branches/ipp-20140904/Ohana/src/libfits/table/F_read_T.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/libfits/table/F_read_T.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/libfits/table/F_read_T.c (revision 37658) @@ -19,8 +19,5 @@ int gfits_fread_ftable (FILE *f, FTable *table, char *extname) { - int j, found; - off_t Nbytes; Header *header; - char tname[80]; header = table[0].header; Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/Makefile =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/Makefile (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/Makefile (revision 37658) @@ -77,6 +77,10 @@ $(SRC)/share_mean_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 \ @@ -131,5 +135,4 @@ $(SRC)/plotstuff.$(ARCH).o \ $(SRC)/ICRF.$(ARCH).o \ -$(SRC)/FrameCorrection.$(ARCH).o \ $(SRC)/relastroVisual.$(ARCH).o \ $(SRC)/syncfile.$(ARCH).o \ Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h (revision 37658) @@ -151,4 +151,19 @@ } 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 @@ -203,4 +218,5 @@ int USE_ICRF_LOCAL; int USE_ICRF_SHFIT; +int USE_ICRF_POLE; int RESET; @@ -578,13 +594,42 @@ 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 FrameCorrection (Catalog *catalog, int Ncatalog); -int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dR, SHterms *dD); +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); -int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords); -int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE); +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); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ConfigInit.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ConfigInit.c (revision 37658) @@ -41,4 +41,5 @@ 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: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c (revision 37658) @@ -1,5 +1,3 @@ # include "relastro.h" - -int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection); /* The sequence of calibration is something like this: @@ -17,90 +15,73 @@ */ -FrameCorrectionSet *FrameCorrectionSetInit () { - - FrameCorrectionSet *set = NULL; - ALLOCATE (set, FrameCorrectionType, 1); - - set->frame = NULL; - set->map = NULL; - return set; -} - -void FrameCorrectionSetFree (FrameCorrectionSet *set) { - - if (!set) return; - - FrameCorrectionFree(set->frame); - AstromOffsetMapFree (set->map); - free set; -} - -FrameCorrectionType *FrameCorrectionInit (double scale) { +int FrameCorrectionParallelMaster (RegionHostTable *regionHosts) { + + if (!USE_ICRF_CORRECT) return TRUE; 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 FrameCorrection (Catalog *catalog, int Ncatalog) { + for (i = 0; i < NLOOP; i++) { + 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); + } + return TRUE; +} + +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 @@ -114,9 +95,5 @@ // fixed line of DEC; for the SAS, the boundary is the projection 'image' - Coords coords; - AstromOffsetMap *map = NULL; - FrameCorrectionType *frame = NULL; - - int POLE = FALSE; + FrameCorrectionSet *set = FrameCorrectionSetInit(); if (USE_ICRF_LOCAL) { @@ -134,8 +111,8 @@ double Ro = 330.0; double Do = 0.0; - - // double Ro = 0.0; - // double Do = 90.0; - // POLE = TRUE; + if (USE_ICRF_POLE) { + Ro = 0.0; + Do = 90.0; + } double scale = 2.0; // degrees per patch @@ -143,16 +120,21 @@ int Ny = 10; - map = AstromOffsetMapInit (Nx, Ny); - map->dX = 1.0; // scale from projection (in arcsec) to correction patches - map->dY = 1.0; - - 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 - - FrameCorrectionFitLocal (catalog, Ncatalog, map, &coords); + 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; } @@ -167,34 +149,20 @@ INITTIME; - FrameCorrectionFitSH (catalog, Ncatalog, dRc, dDc); + FrameCorrectionFitSH (icrfobj, dRc, dDc); MARKTIME ("done with FitSH: %f sec\n", dtime); double pltscale = 0.5; // degrees per pixel - frame = FrameCorrectionInit (pltscale); + FrameCorrectionType *frame = FrameCorrectionInit (pltscale); FrameCorrectionFromSH (frame, dRc, dDc); MARKTIME ("done with FromSH: %f sec\n", dtime); - // XXX write out an image to represent the correction - static int version = 0; - char filenameR[1024], filenameD[1024]; - snprintf (filenameR, 1024, "%s/map.dR.%02d.dat", CATDIR, version); - snprintf (filenameD, 1024, "%s/map.dD.%02d.dat", CATDIR, version); - version ++; - - FrameCorrectionWriteImageSH (frame, filenameR, TRUE); - FrameCorrectionWriteImageSH (frame, filenameD, FALSE); - SHtermsFree (dRc); SHtermsFree (dDc); - } - - // Now apply the correction to all of the average.R,D values - FrameCorrectionApply (catalog, Ncatalog, frame, map, &coords, POLE); - - FrameCorrectionFree(frame); - AstromOffsetMapFree (map); - - return TRUE; + + set->frame = frame; + } + + return set; } @@ -205,7 +173,5 @@ // average.R' = average.R - dR(R,d), average.D' = average.D - dD(R,d) -int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit); - -int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dRc, SHterms *dDc) { +int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc) { int i; @@ -213,46 +179,14 @@ myAssert (dRc->lmax == dDc->lmax, "dR and dD must match\n"); - int Nicrf = ICRFmax(); - fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, Nicrf); - - double *Rave, *Dave, *dRoff, *dDoff; - - int Npts = 0; - ALLOCATE (Rave, double, Nicrf); - ALLOCATE (Dave, double, Nicrf); - ALLOCATE (dRoff, double, Nicrf); - ALLOCATE (dDoff, double, Nicrf); + fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, icrfobj->Nicrfobj); 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++) { - - int cat, meas, ave; - 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); - double dD = 3600.0*(average->D - measure->D); - - if (isnan(dR)) continue; - if (isnan(dD)) continue; - - - Rave[Npts] = average->R; - Dave[Npts] = average->D; - dRoff[Npts] = dR; - dDoff[Npts] = dD; - Npts ++; - } - - SHfitWithMask (Rave, Dave, dRoff, mask, Npts, dRc); - SHfitWithMask (Rave, Dave, dDoff, mask, Npts, dDc); + ALLOCATE (mask, int, icrfobj->Nicrfobj); + memset (mask, 0, icrfobj->Nicrfobj*sizeof(int)); + + 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 @@ -266,7 +200,7 @@ version ++; - for (i = 0; i < Nicrf; i++) { + for (i = 0; i < icrfobj->Nicrfobj; i++) { - SHtermsForRD (Ylm, Rave[i], Dave[i]); + SHtermsForRD (Ylm, icrfobj->Rave[i], icrfobj->Dave[i]); double dRfit = 0.0; @@ -278,12 +212,9 @@ 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\n", Rave[i], Dave[i], dRoff[i], dDoff[i], dRfit, dDfit); + fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f\n", + icrfobj->Rave[i], icrfobj->Dave[i], icrfobj->dRoff[i], icrfobj->dDoff[i], dRfit, dDfit); } fclose (f); - free (Rave); - free (Dave); - free (dRoff); - free (dDoff); free (mask); @@ -426,115 +357,17 @@ } -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 -int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection) { - - Header header; - Matrix matrix; - 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); - - 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; - } - } - gfits_write_header (filename, &header); - gfits_write_matrix (filename, &matrix); - return TRUE; -} - -int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords) { +int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords) { int i; double Xave, Yave, Xmeas, Ymeas; float *X, *Y, *dX, *dY; - - int Nicrf = ICRFmax(); - if (!Nicrf) { - fprintf (stderr, "no matched ICRF quasars for reference correction, skipping\n"); - return FALSE; - } - + float *dXf, *dYf, *dPos; double *Rave, *Dave, *Rmeas, *Dmeas; - float *dXf, *dYf, *dPos; + + int Nicrf = icrfobj->Nicrfobj; + + AstromOffsetMap *map = coords->offsetMap; int Npts = 0; - ALLOCATE (Rave, double, Nicrf); - ALLOCATE (Dave, double, Nicrf); - ALLOCATE (Rmeas, double, Nicrf); - ALLOCATE (Dmeas, double, Nicrf); - ALLOCATE (X, float, Nicrf); ALLOCATE (Y, float, Nicrf); @@ -546,4 +379,10 @@ 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); @@ -553,12 +392,6 @@ 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? - // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax) - int status = RD_to_XY (&Xave, &Yave, average->R, average->D, coords); + int status = RD_to_XY (&Xave, &Yave, icrfobj->Rave[i], icrfobj->Dave[i], coords); if (!status) continue; if (Xave < 0.0) continue; @@ -567,5 +400,7 @@ if (Yave > map->Ny) continue; - RD_to_XY (&Xmeas, &Ymeas, measure->R, measure->D, coords); + double Rm = icrfobj->Rave[i] - icrfobj->dRoff[i] / 3600.0; + double Dm = icrfobj->Dave[i] - icrfobj->dDoff[i] / 3600.0; + RD_to_XY (&Xmeas, &Ymeas, Rm, Dm, coords); // record these in arcsec or degree? @@ -575,8 +410,8 @@ dX[Npts] = Xave - Xmeas; dY[Npts] = Yave - Ymeas; - Rave[Npts] = average->R; - Dave[Npts] = average->D; - Rmeas[Npts] = measure->R; - Dmeas[Npts] = measure->D; + Rave[Npts] = icrfobj->Rave[i]; + Dave[Npts] = icrfobj->Dave[i]; + Rmeas[Npts] = Rm; + Dmeas[Npts] = Dm; Npts ++; } @@ -651,46 +486,25 @@ } - // write out the maps as images - { - Header header; - Matrix matrix; - - 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; - - int ix, iy; - for (ix = 0; ix < map->Nx; ix++) { - for (iy = 0; iy < map->Ny; iy++) { - buffer[ix + map->Nx*iy] = map->dXv[ix][iy]; - } - } - - gfits_write_header ("map.dX.fits", &header); - gfits_write_matrix ("map.dX.fits", &matrix); - - for (ix = 0; ix < map->Nx; ix++) { - for (iy = 0; iy < map->Ny; iy++) { - buffer[ix + map->Nx*iy] = map->dYv[ix][iy]; - } - } - gfits_write_header ("map.dY.fits", &header); - gfits_write_matrix ("map.dY.fits", &matrix); - } - - return TRUE; -} - -int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE) { + 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->offsetMap; int i, j; @@ -706,7 +520,7 @@ // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax) - int useLocal = RD_to_XY (&Xave, &Yave, average->R, average->D, coords); - if (POLE) { - useLocal = (average->D > 85.0); + int useLocal = RD_to_XY (&Xave, &Yave, R, D, coords); + if (USE_ICRF_POLE) { + useLocal = (D > 85.0); } else { useLocal = (useLocal && (Xave > 0.0)); @@ -739,7 +553,8 @@ 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", average->R, average->D, dR, dD); + fprintf (stderr, "skip: %10.6f %10.6f : %7.3f %7.3f\n", R, D, dR, dD); continue; } Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionIO.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionIO.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionIO.c (revision 37658) @@ -1,17 +1,17 @@ # include "relastro.h" -typedef struct { - AstromOffsetMap *map; - FrameCorrectionType *frame; -} FrameCorrectionSet; +// 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 *FrameCorrectionLoad(char *filename) { +FrameCorrectionSet *FrameCorrectionSetLoad(char *filename) { Header header; - Header theader; Matrix matrix; - FTable ftable; FILE *f = fopen (filename, "r"); @@ -43,5 +43,5 @@ exit (2); } - if (!gfits_scan (&header, "LOCAL_DR", "%s", 1, LOCAL_DR_)) { + if (!gfits_scan (&header, "LOCAL_DR", "%s", 1, LOCAL_DR)) { fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field LOCAL_DR\n", filename); exit (2); @@ -51,6 +51,8 @@ exit (2); } - - Header Xheader; + gfits_free_header (&header); + + FrameCorrectionSet *set = FrameCorrectionSetInit (); + if (strcmp(SHMAP_DR, "NONE")) { if (!strcmp (SHMAP_DD, "NONE")) { @@ -58,9 +60,9 @@ exit (2); } - if (!gfits_find_Xheader(f, &Xheader, SHMAP_DR)) { + 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, &Xheader)) { + if (!gfits_fread_matrix (f, &matrix, &header)) { if (VERBOSE) fprintf (stderr, "can't read SHMAP_DR extension image\n"); gfits_free_header (&header); @@ -69,9 +71,13 @@ } // convert image to frame - if (!gfits_find_Xheader(f, &Xheader, SHMAP_DD)) { + 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, &Xheader)) { + if (!gfits_fread_matrix (f, &matrix, &header)) { if (VERBOSE) fprintf (stderr, "can't read SHMAP_DD extension image\n"); gfits_free_header (&header); @@ -79,4 +85,7 @@ return NULL; } + set->frame = FrameCorrectionImageToSH (&header, &matrix, set->frame, TRUE); + gfits_free_header (&header); + gfits_free_matrix (&matrix); } else { fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DD but not SHMAP_DR\n", filename); @@ -84,15 +93,15 @@ } - 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, &Xheader, SHMAP_DR)) { - fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DR extension %s\n", filename, SHMAP_DR); - exit (2); - } - if (!gfits_fread_matrix (f, &matrix, &Xheader)) { - if (VERBOSE) fprintf (stderr, "can't read SHMAP_DR extension image\n"); + 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); @@ -100,27 +109,40 @@ } // convert image to frame - if (!gfits_find_Xheader(f, &Xheader, SHMAP_DD)) { - fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DD extension %s\n", filename, SHMAP_DD); - exit (2); - } - if (!gfits_fread_matrix (f, &matrix, &Xheader)) { - if (VERBOSE) fprintf (stderr, "can't read SHMAP_DD extension image\n"); - gfits_free_header (&header); - fclose (f); - return NULL; - } - } else { - fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DD but not SHMAP_DR\n", filename); - exit (2); - } + 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 { + 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 i; + int status; Header header; - Header theader; Matrix matrix; - FTable ftable; + + 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); @@ -129,42 +151,86 @@ gfits_create_matrix (&header, &matrix); - FILE *f = fopen (filename, "w"); - if (!f) { - fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename); - return FALSE; - } - - int status; + // 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 meanpos %s\n", filename); + 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 meanpos %s\n", filename); - - if (set->frame) { - gfits_modify (&header, "SHMAP_DR", "%s", "SHMAP_DR", 1); - gfits_modify (&header, "SHMAP_DD", "%s", "SHMAP_DD", 1); - status = gfits_fwrite_header (f, &Xheader); - 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); - } - if (set->map) { - gfits_modify (&header, "LOCAL_DR", "%s", "LOCAL_DR", 1); - gfits_modify (&header, "LOCAL_DD", "%s", "LOCAL_DD", 1); - status = gfits_fwrite_header (f, &Xheader); - 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); - } - - // convert frame -> image - // convert map -> image + 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); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionUtils.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 37658) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionUtils.c (revision 37658) @@ -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: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ICRFobjIO.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ICRFobjIO.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ICRFobjIO.c (revision 37658) @@ -1,13 +1,14 @@ # include "relastro.h" +// column OUT refers to a variable which is already defined # define GET_COLUMN(OUT,NAME,TYPE) \ - TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \ + 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, off_t *nicrfobj) { +ICRFobj *ICRFobjLoad(char *filename) { - int i, Ncol; + int Ncol; off_t Nrow; Header header; @@ -16,5 +17,4 @@ FTable ftable; - *nicrfobj = 0; ICRFobj *icrfobj = NULL; @@ -53,24 +53,13 @@ fclose (f); - char type[16]; + char type[16]; // used in GET_COLUMN - GET_COLUMN (R_AVE, "R_AVE", double); - GET_COLUMN (D_AVE, "D_AVE", double); - GET_COLUMN (R_MEAS, "R_MEAS", double); - GET_COLUMN (D_MEAS, "D_MEAS", double); + ALLOCATE (icrfobj, ICRFobj, 1); - ALLOCATE (icrfobj, ICRFobj, Nrow); - for (i = 0; i < Nrow; i++) { - icrfobj[i].Rave = Rave[i]; - icrfobj[i].Dave = Dave[i]; - icrfobj[i].Rmeas = Rmeas[i]; - icrfobj[i].Dmeas = Dmeas[i]; - } + 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 (Rave ); - free (Dave ); - free (Rmeas); - free (Dmeas); // free FITS table pieces... @@ -80,5 +69,5 @@ gfits_free_table (&ftable); - *nicrfobj = Nrow; + icrfobj->Nicrfobj = Nrow; return icrfobj; } @@ -91,7 +80,6 @@ } -int ICRFobjSave(char *filename, ICRFobj *icrfobj, off_t Nicrfobj) { +int ICRFobjSave(char *filename, ICRFobj *icrfobj) { - int i; Header header; Header theader; @@ -106,38 +94,17 @@ 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_MEAS", "meas position, ra", "degrees", 1.0, 0.0); - gfits_define_bintable_column (&theader, "D", "D_MEAS", "meas position, dec", "degrees", 1.0, 0.0); + 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); - double *Rave, *Dave, *Rmeas, *Dmeas; - - // create intermediate storage arrays - ALLOCATE (Rave , double, Nicrfobj); - ALLOCATE (Dave , double, Nicrfobj); - ALLOCATE (Rmeas, double, Nicrfobj); - ALLOCATE (Dmeas, double, Nicrfobj); - - // assign the storage arrays - for (i = 0; i < Nicrfobj; i++) { - Rave [i] = icrfobj[i].Rave ; - Dave [i] = icrfobj[i].Dave ; - Rmeas[i] = icrfobj[i].Rmeas; - Dmeas[i] = icrfobj[i].Dmeas; - } - // add the columns to the output array - gfits_set_bintable_column (&theader, &ftable, "R_AVE", Rave , Nicrfobj); - gfits_set_bintable_column (&theader, &ftable, "D_AVE", Dave , Nicrfobj); - gfits_set_bintable_column (&theader, &ftable, "R_MEAS", Rmeas, Nicrfobj); - gfits_set_bintable_column (&theader, &ftable, "D_MEAS", Dmeas, Nicrfobj); - - free (Rave ); - free (Dave ); - free (Rmeas ); - free (Dmeas ); + 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"); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c (revision 37658) @@ -100,4 +100,9 @@ strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS); + 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"); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_catalogs.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_catalogs.c (revision 37658) @@ -171,4 +171,5 @@ 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); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c (revision 37658) @@ -79,5 +79,5 @@ // if GALAXY_MODEL is selected, we want to delay the frame correction until we have // applied the galaxy model a couple of times. - FrameCorrection (catalog, Ncatalog); + FrameCorrectionSerial (catalog, Ncatalog); } UpdateChips (catalog, Ncatalog); // measure.X,Y -> R,D, fit image.coords Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c (revision 37658) @@ -88,9 +88,8 @@ for (i = 0; i < NLOOP; i++) { UpdateObjects (catalog, Ncatalog, (i > 0)); - share_icrf_obj (catalog, Ncatalog, regionHosts, i); - applyFrameCorr (catalog, Ncatalog, regionHosts, i); // read and apply the frame correction + FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i); share_mean_pos (catalog, Ncatalog, regionHosts, i); - // XXX FrameCorrection (catalog, Ncatalog); -- this needs to happen on the master machine? slurp_mean_pos (catalog, Ncatalog, regionHosts, i); + UpdateChips (catalog, Ncatalog); share_image_pos (regionHosts, i); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c (revision 37658) @@ -44,4 +44,8 @@ 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"); Index: /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_icrf_obj.c =================================================================== --- /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_icrf_obj.c (revision 37657) +++ /branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_icrf_obj.c (revision 37658) @@ -7,29 +7,15 @@ if (!USE_ICRF_CORRECT) return TRUE; - int Nicrf = ICRFmax(); + ICRFobj *icrfobj = get_ICRF_data (catalog, Ncatalog); - // we need to write an empty table none are available - ICRFobj *icrfobj = NULL; - ALLOCATE (icrfobj, ICRFobj, Nicrf); - - // select the ICRF QSOS and save the necessary data - 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? - - icrfobj[i].Rave = average->R; - icrfobj[i].Dave = average->D; - icrfobj[i].Rmeas = measure->R; - icrfobj[i].Rmeas = measure->R; - } + // 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, Nicrf); + ICRFobjSave (filename, icrfobj); + free (filename); + free (icrfobj); - free (filename); char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "icrfobj.sync"); @@ -42,9 +28,8 @@ // load mean and specific positions for all known ICRF objects from remote hosts. // this function occurs on the master host in parallel-region processing -int slurp_icrf_obj (RegionHostTable *regionHosts, int nloop) { +ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop) { off_t i; - int Nicrfobj = 0; ICRFobj *icrfobj = NULL; ALLOCATE (icrfobj, ICRFobj, 1); @@ -57,36 +42,86 @@ free (syncfile); - off_t Nsubset = 0; char *filename = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.fits"); - ICRFobj *icrfobjSubset = ICRFobjLoad (filename, &Nsubset); + ICRFobj *icrfobjSubset = ICRFobjLoad (filename); free (filename); // merge_mean_pos reallocs icrfobj and frees the input icrfobjSubset - icrfobj = merge_icrf_obj (icrfobj, &Nicrfobj, icrfobjSubset, Nsubset); + icrfobj = merge_icrf_obj (icrfobj, icrfobjSubset); } - - frame_correction_icrfobj (); - - save_frame_correction (); - - free (icrfobj); - - fprintf (stderr, "DONE grabbing mean object pos from other hosts...\n"); - - return TRUE; + return icrfobj; } -ICRFobj *merge_icrf_obj (ICRFobj *target, int *ntarget, ICRFobj *source, int Nsource) { +ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source) { off_t i; - REALLOCATE (target, ICRFobj, *ntarget + Nsource); - for (i = 0; i < Nsource; i++) { - off_t n = i + *ntarget; - target[n] = source[i]; + REALLOCATE (target, ICRFobj, 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]; } - free (source); + target->Nicrfobj += source->Nicrfobj; - *ntarget += Nsource; + 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); + 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->dRoff[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; +} +