IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37658


Ignore:
Timestamp:
Nov 21, 2014, 4:50:19 PM (12 years ago)
Author:
eugene
Message:

substantial rework to fully support parallel ops (needs testing)

Location:
branches/eam_branches/ipp-20140904/Ohana/src
Files:
1 added
16 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/libfits/Makefile

    r37656 r37658  
    3232$(HEAD)/F_scan.$(ARCH).o                        $(HEAD)/F_print.$(ARCH).o     \
    3333$(HEAD)/F_copy_H.$(ARCH).o                      $(HEAD)/F_delete.$(ARCH).o    \
    34 $(HEAD)/F_find_XH.$(ARCH).o                     $(HEAD)/F_init_H.$(ARCH).o    \
     34$(HEAD)/F_find_H.$(ARCH).o                      $(HEAD)/F_init_H.$(ARCH).o    \
    3535$(HEAD)/F_convert_H.$(ARCH).o                   $(HEAD)/F_read_XH.$(ARCH).o   \
    3636$(HEAD)/version.$(ARCH).o

  • branches/eam_branches/ipp-20140904/Ohana/src/libfits/header/F_find_H.c

    r37656 r37658  
    55int gfits_find_Xheader (FILE *f, Header *header, char *extname) {
    66
    7   int j, found;
    8   off_t Nbytes;
    9   Header *header;
    107  char tname[80];
    118

  • branches/eam_branches/ipp-20140904/Ohana/src/libfits/include/gfitsio.h

    r37656 r37658  
    119119int     gfits_modify_alt               PROTO((Header *header, char *field, char *mode, int N,...));
    120120int     gfits_print_alt                PROTO((Header *header, char *field, char *mode, int N,...));
    121 int     gfits_find_Xheader             PROTO((char *filename, Header *header, char *extname));
     121int     gfits_find_Xheader             PROTO((FILE *f, Header *header, char *extname));
    122122int     gfits_read_Xheader             PROTO((char *filename, Header *header, int N));
    123123int     gfits_read_header              PROTO((char *filename, Header *header));

  • branches/eam_branches/ipp-20140904/Ohana/src/libfits/table/F_read_T.c

    r37656 r37658  
    1919int gfits_fread_ftable (FILE *f, FTable *table, char *extname) {
    2020
    21   int j, found;
    22   off_t Nbytes;
    2321  Header *header;
    24   char tname[80];
    2522
    2623  header = table[0].header;

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/Makefile

    r37545 r37658  
    7777$(SRC)/share_mean_pos.$(ARCH).o \
    7878$(SRC)/share_images_pos.$(ARCH).o \
     79$(SRC)/share_icrf_obj.$(ARCH).o \
    7980$(SRC)/ICRF.$(ARCH).o \
     81$(SRC)/ICRFobjIO.$(ARCH).o \
    8082$(SRC)/FrameCorrection.$(ARCH).o \
     83$(SRC)/FrameCorrectionIO.$(ARCH).o \
     84$(SRC)/FrameCorrectionUtils.$(ARCH).o \
    8185$(SRC)/ImagePosIO.$(ARCH).o \
    8286$(SRC)/ImageTable.$(ARCH).o \
    … …  
    131135$(SRC)/plotstuff.$(ARCH).o           \
    132136$(SRC)/ICRF.$(ARCH).o \
    133 $(SRC)/FrameCorrection.$(ARCH).o \
    134137$(SRC)/relastroVisual.$(ARCH).o \
    135138$(SRC)/syncfile.$(ARCH).o \

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/include/relastro.h

    r37655 r37658  
    151151} FrameCorrectionType;
    152152
     153typedef struct {
     154  // AstromOffsetMap *map;
     155  Coords *coords; // carries a pointer to the AstromOffsetMap
     156  FrameCorrectionType *frame;
     157} FrameCorrectionSet;
     158
     159
     160typedef struct {
     161  double *Rave;
     162  double *Dave;
     163  double *dRoff;
     164  double *dDoff;
     165  int Nicrfobj;
     166} ICRFobj;
     167
    153168/* global variables set in parameter file */
    154169# define DVO_MAX_PATH 1024
    … …  
    203218int    USE_ICRF_LOCAL;
    204219int    USE_ICRF_SHFIT;
     220int    USE_ICRF_POLE;
    205221
    206222int    RESET;
    … …  
    578594int select_catalog_ICRF (Catalog *catalog, int Ncatalog);
    579595
     596void ICRFobjFree (ICRFobj *icrfobj);
     597ICRFobj *get_ICRF_data (Catalog *catalog, int Ncatalog);
     598ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source);
     599ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop);
     600int share_icrf_obj (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     601
     602int ICRFobjSave(char *filename, ICRFobj *icrfobj);
     603ICRFobj *ICRFobjLoad(char *filename);
     604
     605int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit);
     606
    580607FrameCorrectionType *FrameCorrectionInit (double scale);
    581608void FrameCorrectionFree (FrameCorrectionType *frame);
    582609
    583 int FrameCorrection (Catalog *catalog, int Ncatalog);
    584 int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dR, SHterms *dD);
     610int FrameCorrectionParallelMaster (RegionHostTable *regionHosts);
     611int FrameCorrectionParallelSlave (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop);
     612int FrameCorrectionSerial (Catalog *catalog, int Ncatalog);
     613
     614FrameCorrectionSet *FrameCorrectionMeasure (ICRFobj *icrfobj);
     615
    585616int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD);
    586617
    587 int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords);
    588 int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE);
     618Coords *FrameCorrectionImageToMap (Header *header, Matrix *matrix, Coords *coords, int raDirection);
     619int FrameCorrectionMapToImage (Header *header, Matrix *matrix, Coords *coords, int raDirection);
     620int FrameCorrectionSHtoImage (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection);
     621FrameCorrectionType *FrameCorrectionImageToSH (Header *header, Matrix *matrix, FrameCorrectionType *frame, int raDirection);
     622int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD);
     623void FrameCorrectionFree (FrameCorrectionType *frame);
     624FrameCorrectionType *FrameCorrectionInit (double scale);
     625void FrameCorrectionSetFree (FrameCorrectionSet *set);
     626FrameCorrectionSet *FrameCorrectionSetInit ();
     627
     628int FrameCorrectionSetSave(char *filename, FrameCorrectionSet *set);
     629FrameCorrectionSet *FrameCorrectionSetLoad(char *filename);
     630
     631int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc);
     632int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords);
     633int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, Coords *coords);
    589634
    590635int dump_stardata_pts (StarData *raw, int Npts, char *filename);

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ConfigInit.c

    r37655 r37658  
    4141  if (!ScanConfig (config, "USE_ICRF_LOCAL",      "%d", 0, &USE_ICRF_LOCAL))       USE_ICRF_LOCAL = FALSE;
    4242  if (!ScanConfig (config, "USE_ICRF_SHFIT",      "%d", 0, &USE_ICRF_SHFIT))       USE_ICRF_SHFIT = FALSE;
     43  if (!ScanConfig (config, "USE_ICRF_POLE",       "%d", 0, &USE_ICRF_POLE))        USE_ICRF_POLE = FALSE;
    4344
    4445  // force CATDIR to be absolute (so parallel mode will work)

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c

    r37657 r37658  
    11# include "relastro.h"
    2 
    3 int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection);
    42
    53/* The sequence of calibration is something like this:
    … …  
    1715*/
    1816
    19 FrameCorrectionSet *FrameCorrectionSetInit () {
    20 
    21   FrameCorrectionSet *set = NULL;
    22   ALLOCATE (set, FrameCorrectionType, 1);
    23 
    24   set->frame = NULL;
    25   set->map = NULL;
    26   return set;
    27 }
    28 
    29 void FrameCorrectionSetFree (FrameCorrectionSet *set) {
    30 
    31   if (!set) return;
    32 
    33   FrameCorrectionFree(set->frame);
    34   AstromOffsetMapFree (set->map);
    35   free set;
    36 }
    37 
    38 FrameCorrectionType *FrameCorrectionInit (double scale) {
     17int FrameCorrectionParallelMaster (RegionHostTable *regionHosts) {
     18
     19  if (!USE_ICRF_CORRECT) return TRUE;
    3920
    4021  int i;
    41 
    42   FrameCorrectionType *frame = NULL;
    43   ALLOCATE (frame, FrameCorrectionType, 1);
    44  
    45   frame->scale = scale; // degrees / pixel
    46   frame->Ndec = 180.0 / scale;
    47   if (frame->Ndec * scale < 180.0) frame->Ndec ++;
    48 
    49   // binD = (DEC + 89.0) / scale
    50   ALLOCATE (frame->Nra,  int, frame->Ndec);
    51   ALLOCATE (frame->dR,   double, frame->Ndec);
    52   ALLOCATE (frame->Roff, double *, frame->Ndec);
    53   ALLOCATE (frame->Doff, double *, frame->Ndec);
    54 
    55   // I am going to avoid corrections at |D| = 90
    56 
    57   // It may be sleazy, but since I'm using SH, I'm going to just correct the region at |D|
    58   // > 80 with a 2D grid of corrections to X,Y projections
    59 
    60   for (i = 0; i < frame->Ndec; i++) {
    61     double D = i * scale - 89.0;
    62     if (D > 90.0) {
    63       frame->dR[i] = 0.0;
    64       frame->Nra[i] = 0;
    65       frame->Roff[i] = NULL;
    66       frame->Doff[i] = NULL;
    67       continue;
    68     }
    69     frame->dR[i] = scale / cos(RAD_DEG*D);
    70     frame->Nra[i] = MAX (0, 360.0 / frame->dR[i]);
    71     if (frame->Nra[i] * frame->dR[i] < 360.0) frame->Nra[i] ++;
    72 
    73     ALLOCATE (frame->Roff[i], double, frame->Nra[i]);
    74     ALLOCATE (frame->Doff[i], double, frame->Nra[i]);
    75     // ohana_memcheck_func (TRUE);
    76     // fprintf (stderr, "alloc : %d : %lx %lx\n", i, (long unsigned int)(size_t *) frame->Roff[i], (long unsigned int) (size_t *) frame->Doff[i]);
    77   }
    78   return frame;
    79 }
    80 
    81 void FrameCorrectionFree (FrameCorrectionType *frame) {
    82 
    83   if (!frame) return;
    84 
    85   int i;
    86 
    87   for (i = 0; i < frame->Ndec; i++) {
    88     if (frame->Roff[i]) free (frame->Roff[i]);
    89     if (frame->Doff[i]) free (frame->Doff[i]);
    90   }
    91 
    92   if (frame->dR)   free (frame->dR);
    93   if (frame->Nra)  free (frame->Nra);
    94   if (frame->Roff) free (frame->Roff);
    95   if (frame->Doff) free (frame->Doff);
    96 
    97   free (frame);
    98  
    99   return;
    100 }
    101 
    102 int FrameCorrection (Catalog *catalog, int Ncatalog) {
     22  for (i = 0; i < NLOOP; i++) {
     23    ICRFobj *icrfobj = slurp_icrf_obj (regionHosts, i);
     24    FrameCorrectionSet *set = FrameCorrectionMeasure (icrfobj);
     25
     26    char *filename = make_filename (CATDIR, "master", 0, "frame.corr.fits");
     27    FrameCorrectionSetSave (filename, set);
     28    free (filename);
     29
     30    char *syncfile = make_filename (CATDIR, "master", 0, "frame.corr.sync");
     31    update_sync_file (syncfile, i);
     32    free (syncfile);
     33  }
     34  return TRUE;
     35}
     36
     37int FrameCorrectionParallelSlave (Catalog *catalog, int Ncatalog, RegionHostTable *regionHosts, int nloop) {
    10338
    10439  if (!USE_ICRF_CORRECT) return TRUE;
     40
     41  share_icrf_obj (catalog, Ncatalog, regionHosts, nloop);
     42
     43  char *syncfile = make_filename (CATDIR, "master", 0, "frame.corr.sync");
     44  check_sync_file (syncfile, nloop);
     45  free (syncfile);
     46   
     47  char *filename = make_filename (CATDIR, "master", 0, "frame.corr.fits");
     48  FrameCorrectionSet *set = FrameCorrectionSetLoad (filename);
     49  free (filename);
     50
     51  // Now apply the correction to all of the average.R,D values
     52  FrameCorrectionApply (catalog, Ncatalog, set->frame, set->coords);
     53  FrameCorrectionSetFree(set);
     54
     55  return TRUE;
     56}
     57
     58int FrameCorrectionSerial (Catalog *catalog, int Ncatalog) {
     59
     60  if (!USE_ICRF_CORRECT) return TRUE;
     61
     62  ICRFobj *icrfobj = get_ICRF_data (catalog, Ncatalog);
     63  if (!icrfobj->Nicrfobj) {
     64    fprintf (stderr, "no matched ICRF quasars for reference correction, skipping\n");
     65    return FALSE;
     66  }
     67
     68  FrameCorrectionSet *set = FrameCorrectionMeasure (icrfobj);
     69  ICRFobjFree (icrfobj);
     70
     71  // write out an image to represent the correction
     72  static int version = 0;
     73  char filename[1024];
     74  snprintf (filename, 1024, "%s/frame.%03d.corr.fits", CATDIR, version);
     75  FrameCorrectionSetSave (filename, set);
     76  version ++;
     77
     78  // Now apply the correction to all of the average.R,D values
     79  FrameCorrectionApply (catalog, Ncatalog, set->frame, set->coords);
     80  FrameCorrectionSetFree(set);
     81
     82  return TRUE;
     83}
     84
     85FrameCorrectionSet *FrameCorrectionMeasure (ICRFobj *icrfobj) {
    10586
    10687  // we have only a few thousand ICRF reference points.  I need to have a lookup table to
    … …  
    11495  // fixed line of DEC; for the SAS, the boundary is the projection 'image'
    11596
    116   Coords coords;
    117   AstromOffsetMap *map = NULL;
    118   FrameCorrectionType *frame = NULL;
    119 
    120   int POLE = FALSE;
     97  FrameCorrectionSet *set = FrameCorrectionSetInit();
    12198
    12299  if (USE_ICRF_LOCAL) {
    … …  
    134111    double Ro = 330.0;
    135112    double Do =   0.0;
    136 
    137     // double Ro =  0.0;
    138     // double Do = 90.0;
    139     // POLE = TRUE;
     113    if (USE_ICRF_POLE) {
     114      Ro =  0.0;
     115      Do = 90.0;
     116    }
    140117
    141118    double scale = 2.0; // degrees per patch
    … …  
    143120    int Ny = 10;
    144121   
    145     map = AstromOffsetMapInit (Nx, Ny);
    146     map->dX = 1.0; // scale from projection (in arcsec) to correction patches
    147     map->dY = 1.0;
    148  
    149     InitCoords (&coords, "DEC--SIN");
    150     coords.cdelt1 = coords.cdelt2 = scale;
    151     coords.crval1 = Ro; // SAS center
    152     coords.crval2 = Do;
    153     coords.crpix1 = 0.5*Nx; // middle of projection is middle of map
    154     coords.crpix2 = 0.5*Ny; // middle of projection is middle of map
    155 
    156     FrameCorrectionFitLocal (catalog, Ncatalog, map, &coords);
     122    Coords *coords = NULL;
     123    ALLOCATE (coords, Coords, 1);
     124
     125    InitCoords (coords, "DEC--SIN");
     126    coords->cdelt1 = coords->cdelt2 = scale;
     127    coords->crval1 = Ro; // SAS center
     128    coords->crval2 = Do;
     129    coords->crpix1 = 0.5*Nx; // middle of projection is middle of map
     130    coords->crpix2 = 0.5*Ny; // middle of projection is middle of map
     131
     132    // map is a component of coords
     133    coords->offsetMap = AstromOffsetMapInit (Nx, Ny);
     134    coords->offsetMap->dX = 1.0; // scale from projection (in arcsec) to correction patches
     135    coords->offsetMap->dY = 1.0;
     136
     137    FrameCorrectionFitLocal (icrfobj, coords);
     138    set->coords = coords;
    157139  }
    158140
    … …  
    167149    INITTIME;
    168150
    169     FrameCorrectionFitSH (catalog, Ncatalog, dRc, dDc);
     151    FrameCorrectionFitSH (icrfobj, dRc, dDc);
    170152    MARKTIME ("done with FitSH: %f sec\n", dtime);
    171153
    172154    double pltscale = 0.5; // degrees per pixel
    173     frame = FrameCorrectionInit (pltscale);
     155    FrameCorrectionType *frame = FrameCorrectionInit (pltscale);
    174156 
    175157    FrameCorrectionFromSH (frame, dRc, dDc);
    176158    MARKTIME ("done with FromSH: %f sec\n", dtime);
    177159
    178     // XXX write out an image to represent the correction
    179     static int version = 0;
    180     char filenameR[1024], filenameD[1024];
    181     snprintf (filenameR, 1024, "%s/map.dR.%02d.dat", CATDIR, version);
    182     snprintf (filenameD, 1024, "%s/map.dD.%02d.dat", CATDIR, version);
    183     version ++;
    184 
    185     FrameCorrectionWriteImageSH (frame, filenameR, TRUE);
    186     FrameCorrectionWriteImageSH (frame, filenameD, FALSE);
    187 
    188160    SHtermsFree (dRc);
    189161    SHtermsFree (dDc);
    190   }
    191 
    192   // Now apply the correction to all of the average.R,D values
    193   FrameCorrectionApply (catalog, Ncatalog, frame, map, &coords, POLE);
    194 
    195   FrameCorrectionFree(frame);
    196   AstromOffsetMapFree (map);
    197 
    198   return TRUE;
     162
     163    set->frame = frame;
     164  }
     165
     166  return set;
    199167}
    200168
    … …  
    205173// average.R' = average.R - dR(R,d), average.D' = average.D - dD(R,d)
    206174
    207 int SHfitWithMask (double *R, double *D, double *value, int *mask, int Npts, SHterms *fit);
    208 
    209 int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dRc, SHterms *dDc) {
     175int FrameCorrectionFitSH (ICRFobj *icrfobj, SHterms *dRc, SHterms *dDc) {
    210176
    211177  int i;
    … …  
    213179  myAssert (dRc->lmax == dDc->lmax, "dR and dD must match\n");
    214180
    215   int Nicrf = ICRFmax();
    216   fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, Nicrf);
    217 
    218   double *Rave, *Dave, *dRoff, *dDoff;
    219 
    220   int Npts = 0;
    221   ALLOCATE (Rave,  double, Nicrf);
    222   ALLOCATE (Dave,  double, Nicrf);
    223   ALLOCATE (dRoff, double, Nicrf);
    224   ALLOCATE (dDoff, double, Nicrf);
     181  fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, icrfobj->Nicrfobj);
    225182
    226183  int *mask = NULL;
    227   ALLOCATE (mask, int, Nicrf);
    228   memset (mask, 0, Nicrf*sizeof(int));
    229 
    230   // select the ICRF QSOS and save the necessary data
    231   for (i = 0; i < Nicrf; i++) {
    232 
    233     int cat, meas, ave;
    234     ICRFdata (i, &cat, &ave, &meas);
    235 
    236     Average *average = &catalog[cat].average[ave];
    237     MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
    238 
    239     // record these in arcsec or degree?
    240     // correct for cos(D) or not?
    241     double dR = 3600.0*(average->R - measure->R);
    242     double dD = 3600.0*(average->D - measure->D);
    243 
    244     if (isnan(dR)) continue;
    245     if (isnan(dD)) continue;
    246 
    247 
    248     Rave[Npts] = average->R;
    249     Dave[Npts] = average->D;
    250     dRoff[Npts] = dR;
    251     dDoff[Npts] = dD;
    252     Npts ++;
    253   }
    254 
    255   SHfitWithMask (Rave, Dave, dRoff, mask, Npts, dRc);
    256   SHfitWithMask (Rave, Dave, dDoff, mask, Npts, dDc);
     184  ALLOCATE (mask, int, icrfobj->Nicrfobj);
     185  memset (mask, 0, icrfobj->Nicrfobj*sizeof(int));
     186
     187  SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dRoff, mask, icrfobj->Nicrfobj, dRc);
     188  SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dDoff, mask, icrfobj->Nicrfobj, dDc);
     189
     190  // the bit below is for testing
    257191
    258192  // allocate an SHterms structure to hold the Ylm values
    … …  
    266200  version ++;
    267201
    268   for (i = 0; i < Nicrf; i++) {
     202  for (i = 0; i < icrfobj->Nicrfobj; i++) {
    269203   
    270     SHtermsForRD (Ylm, Rave[i], Dave[i]);
     204    SHtermsForRD (Ylm, icrfobj->Rave[i], icrfobj->Dave[i]);
    271205     
    272206    double dRfit = 0.0;
    … …  
    278212      dDfit += dDc->Fr[j]*Ylm->Fr[j] + dDc->Fi[j]*Ylm->Fi[j];
    279213    }
    280     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);
     214    fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f\n",
     215             icrfobj->Rave[i], icrfobj->Dave[i], icrfobj->dRoff[i], icrfobj->dDoff[i], dRfit, dDfit);
    281216  }
    282217  fclose (f);
    283218
    284   free (Rave);
    285   free (Dave);
    286   free (dRoff);
    287   free (dDoff);
    288219  free (mask);
    289220
    … …  
    426357}
    427358
    428 int FrameCorrectionFromSH (FrameCorrectionType *frame, SHterms *dR, SHterms *dD) {
    429 
    430   myAssert (dR->lmax == dD->lmax, "dR and dD must match\n");
    431 
    432   // allocate an SHterms structure to hold the Ylm values
    433   SHterms *SH = SHtermsInit (dR->lmax);
    434 
    435   int i, j, k;
    436   // Now I need to apply the correction coefficients to SH images to generate an image
    437   for (i = 0; i < frame->Ndec; i++) {
    438     double D = i * frame->scale - 89.0;
    439     for (j = 0; j < frame->Nra[i]; j++) {
    440       double R = j * frame->dR[i];
    441      
    442       SHtermsForRD (SH, R, D);
    443      
    444       frame->Roff[i][j] = 0.0;
    445       frame->Doff[i][j] = 0.0;
    446 
    447       for (k = 0; k < SH->Nterms; k++) {
    448         frame->Roff[i][j] += dR->Fr[k]*SH->Fr[k] + dR->Fi[k]*SH->Fi[k];
    449         frame->Doff[i][j] += dD->Fr[k]*SH->Fr[k] + dD->Fi[k]*SH->Fi[k];
    450       }
    451     }     
    452   }
    453   SHtermsFree (SH);
    454   return TRUE;
    455 }
    456 
    457 // write out the maps as images
    458 int FrameCorrectionWriteImageSH (FrameCorrectionType *frame, char *filename, int raDirection) {
    459 
    460   Header header;
    461   Matrix matrix;
    462   Coords coords;
    463    
    464   gfits_init_header (&header);
    465   header.bitpix = -32;
    466   header.Naxes = 2;
    467   header.Naxis[0] = 370.0 / frame->scale;
    468   header.Naxis[1] = 190.0 / frame->scale;
    469 
    470   int Nx = header.Naxis[0];
    471   int Ny = header.Naxis[1];
    472 
    473   gfits_create_header (&header);
    474   gfits_create_matrix (&header, &matrix);
    475 
    476   InitCoords (&coords, "DEC--AIT");
    477   coords.cdelt1 = coords.cdelt2 = frame->scale;
    478   coords.crval1 = 0.0; // SAS center
    479   coords.crval2 = 0.0;
    480   coords.crpix1 = 0.5*Nx; // middle of projection is middle of map
    481   coords.crpix2 = 0.5*Ny; // middle of projection is middle of map
    482    
    483   PutCoords (&coords, &header);
    484 
    485   float *buffer = (float *)matrix.buffer;
    486 
    487   int ix, iy;
    488   for (ix = 0; ix < Nx; ix++) {
    489     for (iy = 0; iy < Ny; iy++) {
    490       double R, D;
    491       int status = XY_to_RD (&R, &D, ix, iy, &coords);
    492       if (!status) continue;
    493 
    494       R = ohana_normalize_angle (R);
    495 
    496       int iD = (D + 89.0) / frame->scale;
    497       if (iD < 0) continue;
    498       if (iD >= frame->Ndec) continue;
    499       int iR = R / frame->dR[iD];
    500       if (iR < 0) continue;
    501       if (iR >= frame->Nra[iD]) continue;
    502 
    503 
    504       float value = NAN;
    505       if (raDirection) {
    506         value = frame->Roff[iD][iR];
    507       } else {
    508         value = frame->Doff[iD][iR];
    509       }
    510       buffer[ix + Nx*iy] = value;
    511     }
    512   }
    513   gfits_write_header (filename, &header);
    514   gfits_write_matrix (filename, &matrix);
    515   return TRUE;
    516 }
    517 
    518 int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords) {
     359int FrameCorrectionFitLocal (ICRFobj *icrfobj, Coords *coords) {
    519360
    520361  int i;
    521362  double Xave, Yave, Xmeas, Ymeas;
    522363  float *X, *Y, *dX, *dY;
    523 
    524   int Nicrf = ICRFmax();
    525   if (!Nicrf) {
    526     fprintf (stderr, "no matched ICRF quasars for reference correction, skipping\n");
    527     return FALSE;
    528   }
    529 
     364  float *dXf, *dYf, *dPos;
    530365  double *Rave, *Dave, *Rmeas, *Dmeas;
    531   float *dXf, *dYf, *dPos;
     366
     367  int Nicrf = icrfobj->Nicrfobj;
     368
     369  AstromOffsetMap *map = coords->offsetMap;
    532370
    533371  int Npts = 0;
    534   ALLOCATE (Rave, double, Nicrf);
    535   ALLOCATE (Dave, double, Nicrf);
    536   ALLOCATE (Rmeas, double, Nicrf);
    537   ALLOCATE (Dmeas, double, Nicrf);
    538 
    539372  ALLOCATE (X,  float, Nicrf);
    540373  ALLOCATE (Y,  float, Nicrf);
    … …  
    546379  ALLOCATE (dPos, float, Nicrf);
    547380
     381  // arrays to store the points actually in the local region
     382  ALLOCATE (Rave, double, Nicrf);
     383  ALLOCATE (Dave, double, Nicrf);
     384  ALLOCATE (Rmeas, double, Nicrf);
     385  ALLOCATE (Dmeas, double, Nicrf);
     386
    548387  int *mask = NULL;
    549388  ALLOCATE (mask, int, Nicrf);
    … …  
    553392  for (i = 0; i < Nicrf; i++) {
    554393
    555     int cat, ave, meas;
    556     ICRFdata (i, &cat, &ave, &meas);
    557 
    558     Average *average = &catalog[cat].average[ave];
    559     MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
    560 
    561394    // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax)
    562     int status = RD_to_XY (&Xave, &Yave, average->R, average->D, coords);
     395    int status = RD_to_XY (&Xave, &Yave, icrfobj->Rave[i], icrfobj->Dave[i], coords);
    563396    if (!status) continue;
    564397    if (Xave < 0.0) continue;
    … …  
    567400    if (Yave > map->Ny) continue;
    568401
    569     RD_to_XY (&Xmeas, &Ymeas, measure->R, measure->D, coords);
     402    double Rm = icrfobj->Rave[i] - icrfobj->dRoff[i] / 3600.0;
     403    double Dm = icrfobj->Dave[i] - icrfobj->dDoff[i] / 3600.0;
     404    RD_to_XY (&Xmeas, &Ymeas, Rm, Dm, coords);
    570405
    571406    // record these in arcsec or degree?
    … …  
    575410    dX[Npts] = Xave - Xmeas;
    576411    dY[Npts] = Yave - Ymeas;
    577     Rave[Npts] = average->R;
    578     Dave[Npts] = average->D;
    579     Rmeas[Npts] = measure->R;
    580     Dmeas[Npts] = measure->D;
     412    Rave[Npts] = icrfobj->Rave[i];
     413    Dave[Npts] = icrfobj->Dave[i];
     414    Rmeas[Npts] = Rm;
     415    Dmeas[Npts] = Dm;
    581416    Npts ++;
    582417  }
    … …  
    651486  }
    652487
    653   // write out the maps as images
    654   {
    655     Header header;
    656     Matrix matrix;
    657    
    658     gfits_init_header (&header);
    659     header.bitpix = -32;
    660     header.Naxes = 2;
    661     header.Naxis[0] = map->Nx;
    662     header.Naxis[1] = map->Ny;
    663 
    664     gfits_create_header (&header);
    665     gfits_create_matrix (&header, &matrix);
    666     PutCoords (coords, &header);
    667 
    668     float *buffer = (float *)matrix.buffer;
    669 
    670     int ix, iy;
    671     for (ix = 0; ix < map->Nx; ix++) {
    672       for (iy = 0; iy < map->Ny; iy++) {
    673         buffer[ix + map->Nx*iy] = map->dXv[ix][iy];
    674       }
    675     }
    676    
    677     gfits_write_header ("map.dX.fits", &header);
    678     gfits_write_matrix ("map.dX.fits", &matrix);
    679 
    680     for (ix = 0; ix < map->Nx; ix++) {
    681       for (iy = 0; iy < map->Ny; iy++) {
    682         buffer[ix + map->Nx*iy] = map->dYv[ix][iy];
    683       }
    684     }
    685     gfits_write_header ("map.dY.fits", &header);
    686     gfits_write_matrix ("map.dY.fits", &matrix);
    687   }
    688 
    689   return TRUE;
    690 }
    691 
    692 int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE) {
     488  free (X);
     489  free (Y);
     490  free (dX);
     491  free (dY);
     492  free (dXf);
     493  free (dYf);
     494  free (dPos);
     495  free (Rave);
     496  free (Dave);
     497  free (Rmeas);
     498  free (Dmeas);
     499  free (mask);
     500
     501  return TRUE;
     502}
     503
     504int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, Coords *coords) {
    693505
    694506  double Xave, Yave;
     507
     508  AstromOffsetMap *map = coords->offsetMap;
    695509
    696510  int i, j;
    … …  
    706520
    707521      // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax)
    708       int useLocal = RD_to_XY (&Xave, &Yave, average->R, average->D, coords);
    709       if (POLE) {
    710         useLocal = (average->D > 85.0);
     522      int useLocal = RD_to_XY (&Xave, &Yave, R, D, coords);
     523      if (USE_ICRF_POLE) {
     524        useLocal = (D > 85.0);
    711525      } else {
    712526        useLocal = (useLocal && (Xave > 0.0));
    … …  
    739553        double dD = frame->Doff[iD][iR];
    740554
     555        // XXX tighten this up??
    741556        // do not apply if the fabs(offset) is > 8 arcsec)
    742557        if ((fabs(dR) > 8.0) || (fabs(dD) > 8.0)) {
    743           fprintf (stderr, "skip: %10.6f %10.6f : %7.3f %7.3f\n", average->R, average->D, dR, dD);
     558          fprintf (stderr, "skip: %10.6f %10.6f : %7.3f %7.3f\n", R, D, dR, dD);
    744559          continue;
    745560        }

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrectionIO.c

    r37657 r37658  
    11# include "relastro.h"
    22
    3 typedef struct {
    4   AstromOffsetMap *map;
    5   FrameCorrectionType *frame;
    6 } FrameCorrectionSet;
     3// STATUS is value expected for success
     4# define CHECK_STATUS(STATUS,MSG,...)                                   \
     5  if (!(STATUS)) {                                                      \
     6    fprintf (stderr, MSG, __VA_ARGS__);                                 \
     7    return FALSE;                                                       \
     8  }
    79
    810// the frame correction consists of N image correction maps
    911// there are local maps and AITOFF maps
    10 FrameCorrectionSet *FrameCorrectionLoad(char *filename) {
     12FrameCorrectionSet *FrameCorrectionSetLoad(char *filename) {
    1113
    1214  Header header;
    13   Header theader;
    1415  Matrix matrix;
    15   FTable ftable;
    1616
    1717  FILE *f = fopen (filename, "r");
    … …  
    4343    exit (2);
    4444  }
    45   if (!gfits_scan (&header, "LOCAL_DR", "%s", 1, LOCAL_DR_)) {
     45  if (!gfits_scan (&header, "LOCAL_DR", "%s", 1, LOCAL_DR)) {
    4646    fprintf (stderr, "ERROR: FrameCorrectionSet %s missing required field LOCAL_DR\n", filename);
    4747    exit (2);
    … …  
    5151    exit (2);
    5252  }
    53 
    54   Header Xheader;
     53  gfits_free_header (&header);
     54   
     55  FrameCorrectionSet *set = FrameCorrectionSetInit ();
     56
    5557  if (strcmp(SHMAP_DR, "NONE")) {
    5658    if (!strcmp (SHMAP_DD, "NONE")) {
    … …  
    5860      exit (2);
    5961    }
    60     if (!gfits_find_Xheader(f, &Xheader, SHMAP_DR)) {
     62    if (!gfits_find_Xheader(f, &header, SHMAP_DR)) {
    6163      fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DR extension %s\n", filename, SHMAP_DR);
    6264      exit (2);
    6365    }
    64     if (!gfits_fread_matrix (f, &matrix, &Xheader)) {
     66    if (!gfits_fread_matrix (f, &matrix, &header)) {
    6567      if (VERBOSE) fprintf (stderr, "can't read SHMAP_DR extension image\n");
    6668      gfits_free_header (&header);
    … …  
    6971    }
    7072    // convert image to frame
    71     if (!gfits_find_Xheader(f, &Xheader, SHMAP_DD)) {
     73    set->frame = FrameCorrectionImageToSH (&header, &matrix, NULL, TRUE);
     74    gfits_free_header (&header);
     75    gfits_free_matrix (&matrix);
     76
     77    if (!gfits_find_Xheader(f, &header, SHMAP_DD)) {
    7278      fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DD extension %s\n", filename, SHMAP_DD);
    7379      exit (2);
    7480    }
    75     if (!gfits_fread_matrix (f, &matrix, &Xheader)) {
     81    if (!gfits_fread_matrix (f, &matrix, &header)) {
    7682      if (VERBOSE) fprintf (stderr, "can't read SHMAP_DD extension image\n");
    7783      gfits_free_header (&header);
    … …  
    7985      return NULL;
    8086    }
     87    set->frame = FrameCorrectionImageToSH (&header, &matrix, set->frame, TRUE);
     88    gfits_free_header (&header);
     89    gfits_free_matrix (&matrix);
    8190  } else {
    8291    fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DD but not SHMAP_DR\n", filename);
    … …  
    8493  }
    8594
    86   if (strcmp(SHMAP_DR, "NONE")) {
    87     if (!strcmp (SHMAP_DD, "NONE")) {
    88       fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DR but not SHMAP_DD\n", filename);
    89       exit (2);
    90     }
    91     if (!gfits_find_Xheader(f, &Xheader, SHMAP_DR)) {
    92       fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DR extension %s\n", filename, SHMAP_DR);
    93       exit (2);
    94     }
    95     if (!gfits_fread_matrix (f, &matrix, &Xheader)) {
    96       if (VERBOSE) fprintf (stderr, "can't read SHMAP_DR extension image\n");
     95  if (strcmp(LOCAL_DR, "NONE")) {
     96    if (!strcmp (LOCAL_DD, "NONE")) {
     97      fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined LOCAL_DR but not LOCAL_DD\n", filename);
     98      exit (2);
     99    }
     100    if (!gfits_find_Xheader(f, &header, LOCAL_DR)) {
     101      fprintf (stderr, "ERROR: FrameCorrectionSet %s missing LOCAL_DR extension %s\n", filename, LOCAL_DR);
     102      exit (2);
     103    }
     104    if (!gfits_fread_matrix (f, &matrix, &header)) {
     105      if (VERBOSE) fprintf (stderr, "can't read LOCAL_DR extension image\n");
    97106      gfits_free_header (&header);
    98107      fclose (f);
    … …  
    100109    }
    101110    // convert image to frame
    102     if (!gfits_find_Xheader(f, &Xheader, SHMAP_DD)) {
    103       fprintf (stderr, "ERROR: FrameCorrectionSet %s missing SHMAP_DD extension %s\n", filename, SHMAP_DD);
    104       exit (2);
    105     }
    106     if (!gfits_fread_matrix (f, &matrix, &Xheader)) {
    107       if (VERBOSE) fprintf (stderr, "can't read SHMAP_DD extension image\n");
    108       gfits_free_header (&header);
    109       fclose (f);
    110       return NULL;
    111     }
    112   } else {
    113     fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined SHMAP_DD but not SHMAP_DR\n", filename);
    114     exit (2);
    115   }
     111    set->coords = FrameCorrectionImageToMap (&header, &matrix, NULL, TRUE);
     112    gfits_free_header (&header);
     113    gfits_free_matrix (&matrix);
     114
     115    if (!gfits_find_Xheader(f, &header, LOCAL_DD)) {
     116      fprintf (stderr, "ERROR: FrameCorrectionSet %s missing LOCAL_DD extension %s\n", filename, LOCAL_DD);
     117      exit (2);
     118    }
     119    if (!gfits_fread_matrix (f, &matrix, &header)) {
     120      if (VERBOSE) fprintf (stderr, "can't read LOCAL_DD extension image\n");
     121      gfits_free_header (&header);
     122      fclose (f);
     123      return NULL;
     124    }
     125    set->coords = FrameCorrectionImageToMap (&header, &matrix, set->coords, FALSE);
     126    gfits_free_header (&header);
     127    gfits_free_matrix (&matrix);
     128  } else {
     129    fprintf (stderr, "ERROR: FrameCorrectionSet %s has defined LOCAL_DD but not LOCAL_DR\n", filename);
     130    exit (2);
     131  }
     132  fclose (f);
     133  return set;
    116134}
    117135
    118136int FrameCorrectionSetSave(char *filename, FrameCorrectionSet *set) {
    119137
    120   int i;
     138  int status;
    121139  Header header;
    122   Header theader;
    123140  Matrix matrix;
    124   FTable ftable;
     141
     142  FILE *f = fopen (filename, "w");
     143  if (!f) {
     144    fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename);
     145    return FALSE;
     146  }
    125147
    126148  gfits_init_header (&header);
    … …  
    129151  gfits_create_matrix (&header, &matrix);
    130152
    131   FILE *f = fopen (filename, "w");
    132   if (!f) {
    133     fprintf (stderr, "ERROR: cannot open meanmag file for output %s\n", filename);
    134     return FALSE;
    135   }
    136 
    137   int status;
     153  // add entries to the PHU
     154  if (set->frame) {
     155    gfits_modify (&header, "SHMAP_DR", "%s", 1, "SHMAP_DR");
     156    gfits_modify (&header, "SHMAP_DD", "%s", 1, "SHMAP_DD");
     157  } else {
     158    gfits_modify (&header, "SHMAP_DR", "%s", 1, "NONE");
     159    gfits_modify (&header, "SHMAP_DD", "%s", 1, "NONE");
     160  }
     161  if (set->coords) {
     162    gfits_modify (&header, "LOCAL_DR", "%s", 1, "LOCAL_DR");
     163    gfits_modify (&header, "LOCAL_DD", "%s", 1, "LOCAL_DD");
     164  } else {
     165    gfits_modify (&header, "LOCAL_DR", "%s", 1, "NONE");
     166    gfits_modify (&header, "LOCAL_DD", "%s", 1, "NONE");
     167  }
     168
     169  // write the PHU header (& empty matrix) to disk
    138170  status = gfits_fwrite_header  (f, &header);
    139   CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename);
     171  CHECK_STATUS (status, "ERROR: cannot write header for frame correction %s\n", filename);
    140172
    141173  status = gfits_fwrite_matrix  (f, &matrix);
    142   CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename);
    143 
    144   if (set->frame) {
    145     gfits_modify (&header, "SHMAP_DR", "%s", "SHMAP_DR", 1);
    146     gfits_modify (&header, "SHMAP_DD", "%s", "SHMAP_DD", 1);
    147   status = gfits_fwrite_header  (f, &Xheader);
    148   CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename);
    149 
    150   status = gfits_fwrite_matrix  (f, &matrix);
    151   CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename);
    152   }
    153   if (set->map) {
    154     gfits_modify (&header, "LOCAL_DR", "%s", "LOCAL_DR", 1);
    155     gfits_modify (&header, "LOCAL_DD", "%s", "LOCAL_DD", 1);
    156   status = gfits_fwrite_header  (f, &Xheader);
    157   CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename);
    158 
    159   status = gfits_fwrite_matrix  (f, &matrix);
    160   CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename);
    161   }
    162 
    163   // convert frame -> image
    164   // convert map -> image
     174  CHECK_STATUS (status, "ERROR: cannot write matrix for correction %s\n", filename);
    165175
    166176  gfits_free_header (&header);
    167177  gfits_free_matrix (&matrix);
    168178
     179  if (set->frame) {
     180    // create and save the RA correction
     181    FrameCorrectionSHtoImage (&header, &matrix, set->frame, TRUE);
     182    gfits_modify (&header, "EXTNAME", "%s", 1, "SHMAP_DR");
     183
     184    status = gfits_fwrite_header  (f, &header);
     185    CHECK_STATUS (status, "ERROR: cannot write header for SHMAP_DR %s\n", filename);
     186
     187    status = gfits_fwrite_matrix  (f, &matrix);
     188    CHECK_STATUS (status, "ERROR: cannot write matrix for SHMAP_DR %s\n", filename);
     189
     190    gfits_free_header (&header);
     191    gfits_free_matrix (&matrix);
     192
     193    // create and save the DEC correction
     194    FrameCorrectionSHtoImage (&header, &matrix, set->frame, FALSE);
     195    gfits_modify (&header, "EXTNAME", "%s", 1, "SHMAP_DD");
     196
     197    status = gfits_fwrite_header  (f, &header);
     198    CHECK_STATUS (status, "ERROR: cannot write header for SHMAP_DD %s\n", filename);
     199
     200    status = gfits_fwrite_matrix  (f, &matrix);
     201    CHECK_STATUS (status, "ERROR: cannot write matrix for SHMAP_DD %s\n", filename);
     202
     203    gfits_free_header (&header);
     204    gfits_free_matrix (&matrix);
     205  }
     206
     207  if (set->coords) {
     208    // create and save the RA correction
     209    FrameCorrectionMapToImage (&header, &matrix, set->coords, TRUE);
     210    gfits_modify (&header, "EXTNAME", "%s", 1, "LOCAL_DR");
     211
     212    status = gfits_fwrite_header  (f, &header);
     213    CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename);
     214
     215    status = gfits_fwrite_matrix  (f, &matrix);
     216    CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename);
     217
     218    gfits_free_header (&header);
     219    gfits_free_matrix (&matrix);
     220
     221    // create and save the DEC correction
     222    FrameCorrectionMapToImage (&header, &matrix, set->coords, FALSE);
     223    gfits_modify (&header, "EXTNAME", "%s", 1, "LOCAL_DD");
     224
     225    status = gfits_fwrite_header  (f, &header);
     226    CHECK_STATUS (status, "ERROR: cannot write header for meanpos %s\n", filename);
     227
     228    status = gfits_fwrite_matrix  (f, &matrix);
     229    CHECK_STATUS (status, "ERROR: cannot write matrix for meanpos %s\n", filename);
     230
     231    gfits_free_header (&header);
     232    gfits_free_matrix (&matrix);
     233  }
     234
    169235  int fd = fileno (f);
    170236

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/ICRFobjIO.c

    r37657 r37658  
    11# include "relastro.h"
    22
     3// column OUT refers to a variable which is already defined
    34# define GET_COLUMN(OUT,NAME,TYPE) \
    4   TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
     5  OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
    56  myAssert (!strcmp(type, #TYPE), "wrong column type");
    67
    78// this is nearly identical to the one in 'uniphot/src' used by setphot / setphot_client.
    89// Here, we use a handful of different columns (if not, we could move to libdvo)
    9 ICRFobj *ICRFobjLoad(char *filename, off_t *nicrfobj) {
     10ICRFobj *ICRFobjLoad(char *filename) {
    1011
    11   int i, Ncol;
     12  int Ncol;
    1213  off_t Nrow;
    1314  Header header;
    … …  
    1617  FTable ftable;
    1718
    18   *nicrfobj = 0;
    1919  ICRFobj *icrfobj = NULL;
    2020
    … …  
    5353  fclose (f);
    5454
    55   char type[16];
     55  char type[16]; // used in GET_COLUMN
    5656
    57   GET_COLUMN (R_AVE,     "R_AVE",          double);
    58   GET_COLUMN (D_AVE,     "D_AVE",          double);
    59   GET_COLUMN (R_MEAS,    "R_MEAS",         double);
    60   GET_COLUMN (D_MEAS,    "D_MEAS",         double);
     57  ALLOCATE (icrfobj, ICRFobj, 1);
    6158
    62   ALLOCATE (icrfobj, ICRFobj, Nrow);
    63   for (i = 0; i < Nrow; i++) {
    64     icrfobj[i].Rave           = Rave[i];
    65     icrfobj[i].Dave           = Dave[i];
    66     icrfobj[i].Rmeas          = Rmeas[i];
    67     icrfobj[i].Dmeas          = Dmeas[i];
    68   }
     59  GET_COLUMN (icrfobj->Rave,     "R_AVE",        double);
     60  GET_COLUMN (icrfobj->Dave,     "D_AVE",        double);
     61  GET_COLUMN (icrfobj->dRoff,    "R_OFF",        double);
     62  GET_COLUMN (icrfobj->dDoff,    "D_OFF",        double);
    6963  fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow);
    70 
    71   free (Rave );
    72   free (Dave );
    73   free (Rmeas);
    74   free (Dmeas);
    7564
    7665  // free FITS table pieces...
    … …  
    8069  gfits_free_table  (&ftable);
    8170
    82   *nicrfobj = Nrow;
     71  icrfobj->Nicrfobj = Nrow;
    8372  return icrfobj;
    8473}
    … …  
    9180  }
    9281
    93 int ICRFobjSave(char *filename, ICRFobj *icrfobj, off_t Nicrfobj) {
     82int ICRFobjSave(char *filename, ICRFobj *icrfobj) {
    9483
    95   int i;
    9684  Header header;
    9785  Header theader;
    … …  
    10694  gfits_create_table_header (&theader, "BINTABLE", "ICRFOBJ");
    10795
    108   gfits_define_bintable_column (&theader, "D", "R_AVE",       "mean position, ra",  "degrees", 1.0, 0.0);
    109   gfits_define_bintable_column (&theader, "D", "D_AVE",       "mean position, dec", "degrees", 1.0, 0.0);
    110   gfits_define_bintable_column (&theader, "D", "R_MEAS",      "meas position, ra",  "degrees", 1.0, 0.0);
    111   gfits_define_bintable_column (&theader, "D", "D_MEAS",      "meas position, dec", "degrees", 1.0, 0.0);
     96  gfits_define_bintable_column (&theader, "D", "R_AVE",      "mean position, ra",  "degrees", 1.0, 0.0);
     97  gfits_define_bintable_column (&theader, "D", "D_AVE",      "mean position, dec", "degrees", 1.0, 0.0);
     98  gfits_define_bintable_column (&theader, "D", "R_OFF",      "position offset, ra",  "arcsec", 1.0, 0.0);
     99  gfits_define_bintable_column (&theader, "D", "D_OFF",      "position offset, dec", "arcsec", 1.0, 0.0);
    112100
    113101  // generate the output array that carries the data
    114102  gfits_create_table (&theader, &ftable);
    115103
    116   double *Rave, *Dave, *Rmeas, *Dmeas;
    117 
    118   // create intermediate storage arrays
    119   ALLOCATE (Rave ,         double,             Nicrfobj);
    120   ALLOCATE (Dave ,         double,             Nicrfobj);
    121   ALLOCATE (Rmeas,         double,             Nicrfobj);
    122   ALLOCATE (Dmeas,         double,             Nicrfobj);
    123 
    124   // assign the storage arrays
    125   for (i = 0; i < Nicrfobj; i++) {
    126     Rave [i]   = icrfobj[i].Rave ;
    127     Dave [i]   = icrfobj[i].Dave ;
    128     Rmeas[i]   = icrfobj[i].Rmeas;
    129     Dmeas[i]   = icrfobj[i].Dmeas;
    130   }
    131 
    132104  // add the columns to the output array
    133   gfits_set_bintable_column (&theader, &ftable, "R_AVE",        Rave ,         Nicrfobj);
    134   gfits_set_bintable_column (&theader, &ftable, "D_AVE",        Dave ,         Nicrfobj);
    135   gfits_set_bintable_column (&theader, &ftable, "R_MEAS",       Rmeas,         Nicrfobj);
    136   gfits_set_bintable_column (&theader, &ftable, "D_MEAS",       Dmeas,         Nicrfobj);
    137 
    138   free (Rave     );
    139   free (Dave     );
    140   free (Rmeas    );
    141   free (Dmeas    );
     105  gfits_set_bintable_column (&theader, &ftable, "R_AVE",       icrfobj->Rave ,         icrfobj->Nicrfobj);
     106  gfits_set_bintable_column (&theader, &ftable, "D_AVE",       icrfobj->Dave ,         icrfobj->Nicrfobj);
     107  gfits_set_bintable_column (&theader, &ftable, "R_OFF",       icrfobj->dRoff,         icrfobj->Nicrfobj);
     108  gfits_set_bintable_column (&theader, &ftable, "D_OFF",       icrfobj->dDoff,         icrfobj->Nicrfobj);
    142109
    143110  FILE *f = fopen (filename, "w");

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/launch_region_hosts.c

    r37038 r37658  
    100100    strextend (command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS);
    101101
     102    strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT);
     103    strextend (command, "-D USE_ICRF_LOCAL %d",   USE_ICRF_LOCAL);
     104    strextend (command, "-D USE_ICRF_SHFIT %d",   USE_ICRF_SHFIT);
     105    strextend (command, "-D USE_ICRF_POLE %d",    USE_ICRF_POLE);
     106
    102107    if (USE_FIXED_PIXCOORDS) strextend (command, "-D USE_FIXED_PIXCOORDS 1");
    103108

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/load_catalogs.c

    r37657 r37658  
    171171    if (FlagOutlier)   strextend (command, "-clip %d", CLIP_THRESH);
    172172    if (ExcludeBogus)  strextend (command, "-exclude-bogus %f", ExcludeBogusRadius);
     173
    173174    if (USE_ICRF_CORRECT) strextend (command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT);
    174175    if (USE_GALAXY_MODEL) strextend (command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL);

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_images.c

    r37655 r37658  
    7979          // if GALAXY_MODEL is selected, we want to delay the frame correction until we have
    8080          // applied the galaxy model a couple of times.
    81           FrameCorrection (catalog, Ncatalog);
     81          FrameCorrectionSerial (catalog, Ncatalog);
    8282        }
    8383        UpdateChips (catalog, Ncatalog);   // measure.X,Y -> R,D, fit image.coords

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_images.c

    r37657 r37658  
    8888      for (i = 0; i < NLOOP; i++) {
    8989        UpdateObjects (catalog, Ncatalog, (i > 0));
    90         share_icrf_obj (catalog, Ncatalog, regionHosts, i);
    91         applyFrameCorr (catalog, Ncatalog, regionHosts, i); // read and apply the frame correction
     90        FrameCorrectionParallelSlave (catalog, Ncatalog, regionHosts, i);
    9291        share_mean_pos (catalog, Ncatalog, regionHosts, i);
    93         // XXX FrameCorrection (catalog, Ncatalog); -- this needs to happen on the master machine?
    9492        slurp_mean_pos (catalog, Ncatalog, regionHosts, i);
     93
    9594        UpdateChips (catalog, Ncatalog);
    9695        share_image_pos (regionHosts, i);

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/relastro_parallel_regions.c

    r36630 r37658  
    4444  if (!launch_region_hosts (regionHosts)) Shutdown ("error launching region hosts");
    4545
     46  // If we are doing FrameCorrection, we have to loop here waiting for ICRF objects from the region
     47  // hosts.  we then measure the frame correction and send back the result
     48  if (!FrameCorrectionParallelMaster (regionHosts)) Shutdown ("error running parallel frame correction");
     49
    4650  // retrieve updated image parameters from the remote hosts (also set Image.mcal)
    4751  if (!slurp_image_pos (regionHosts, -1)) Shutdown ("error loading image updates");

  • branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/share_icrf_obj.c

    r37657 r37658  
    77  if (!USE_ICRF_CORRECT) return TRUE;
    88
    9   int Nicrf = ICRFmax();
     9  ICRFobj *icrfobj = get_ICRF_data (catalog, Ncatalog);
    1010
    11   // we need to write an empty table none are available
    12   ICRFobj *icrfobj = NULL;
    13   ALLOCATE (icrfobj, ICRFobj, Nicrf);
    14 
    15   // select the ICRF QSOS and save the necessary data
    16   for (i = 0; i < Nicrf; i++) {
    17 
    18     int cat, ave, meas;
    19     ICRFdata (i, &cat, &ave, &meas);
    20 
    21     Average *average = &catalog[cat].average[ave];
    22     MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
    23 
    24     icrfobj[i].Rave  = average->R;
    25     icrfobj[i].Dave  = average->D;
    26     icrfobj[i].Rmeas = measure->R;
    27     icrfobj[i].Rmeas = measure->R;
    28   }
     11  // write out the meanmag fits table AND write state in some file
     12  int myHost = regionHosts->index[REGION_HOST_ID];
     13  char *hostname = regionHosts->hosts[myHost].hostname;
    2914
    3015  char *filename = make_filename (CATDIR, hostname, REGION_HOST_ID, "icrfobj.fits");
    31   ICRFobjSave (filename, icrfobj, Nicrf);
     16  ICRFobjSave (filename, icrfobj);
     17  free (filename);
     18
    3219  free (icrfobj);
    33   free (filename);
    3420
    3521  char *syncfile = make_filename (CATDIR, hostname, REGION_HOST_ID, "icrfobj.sync");
    … …  
    4228// load mean and specific positions for all known ICRF objects from remote hosts.
    4329// this function occurs on the master host in parallel-region processing
    44 int slurp_icrf_obj (RegionHostTable *regionHosts, int nloop) {
     30ICRFobj *slurp_icrf_obj (RegionHostTable *regionHosts, int nloop) {
    4531
    4632  off_t i;
    4733
    48   int Nicrfobj = 0;
    4934  ICRFobj *icrfobj = NULL;
    5035  ALLOCATE (icrfobj, ICRFobj, 1);
    … …  
    5742    free (syncfile);
    5843   
    59     off_t Nsubset = 0;
    6044    char *filename = make_filename (CATDIR, regionHosts->hosts[i].hostname, regionHosts->hosts[i].hostID, "icrfobj.fits");
    61     ICRFobj *icrfobjSubset = ICRFobjLoad (filename, &Nsubset);
     45    ICRFobj *icrfobjSubset = ICRFobjLoad (filename);
    6246    free (filename);
    6347
    6448    // merge_mean_pos reallocs icrfobj and frees the input icrfobjSubset
    65     icrfobj = merge_icrf_obj (icrfobj, &Nicrfobj, icrfobjSubset, Nsubset);
     49    icrfobj = merge_icrf_obj (icrfobj, icrfobjSubset);
    6650  }
    67 
    68   frame_correction_icrfobj ();
    69 
    70   save_frame_correction ();
    71 
    72   free (icrfobj);
    73 
    74   fprintf (stderr, "DONE grabbing mean object pos from other hosts...\n");
    75 
    76   return TRUE;
     51  return icrfobj;
    7752}
    7853
    79 ICRFobj *merge_icrf_obj (ICRFobj *target, int *ntarget, ICRFobj *source, int Nsource) {
     54ICRFobj *merge_icrf_obj (ICRFobj *target, ICRFobj *source) {
    8055
    8156  off_t i;
    8257
    83   REALLOCATE (target, ICRFobj, *ntarget + Nsource);
    84   for (i = 0; i < Nsource; i++) {
    85     off_t n = i + *ntarget;
    86     target[n] = source[i];
     58  REALLOCATE (target, ICRFobj, target->Nicrfobj + source->Nicrfobj);
     59  for (i = 0; i < source->Nicrfobj; i++) {
     60    off_t n = i + target->Nicrfobj;
     61    target->Rave[n]  = source->Rave[i];
     62    target->Dave[n]  = source->Dave[i];
     63    target->dRoff[n] = source->dRoff[i];
     64    target->dDoff[n] = source->dDoff[i];
    8765  }
    88   free (source);
     66  target->Nicrfobj += source->Nicrfobj;
    8967
    90   *ntarget += Nsource;
     68  ICRFobjFree (source);
    9169  return (target);
    9270}
     71
     72ICRFobj *get_ICRF_data (Catalog *catalog, int Ncatalog) {
     73
     74  int i;
     75  int Nicrf = ICRFmax();
     76
     77  // we need to write an empty table none are available
     78  ICRFobj *icrfobj = NULL;
     79  ALLOCATE (icrfobj, ICRFobj, 1);
     80
     81  ALLOCATE (icrfobj->Rave,  double, Nicrf);
     82  ALLOCATE (icrfobj->Dave,  double, Nicrf);
     83  ALLOCATE (icrfobj->dRoff, double, Nicrf);
     84  ALLOCATE (icrfobj->dDoff, double, Nicrf);
     85
     86  // select the ICRF QSOS and save the necessary data
     87  int Npts = 0;
     88  for (i = 0; i < Nicrf; i++) {
     89
     90    int cat, ave, meas;
     91    ICRFdata (i, &cat, &ave, &meas);
     92
     93    Average *average = &catalog[cat].average[ave];
     94    MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
     95
     96    // record these in arcsec or degree?
     97    // correct for cos(D) or not?
     98    double dR = 3600.0*(average->R - measure->R);
     99    double dD = 3600.0*(average->D - measure->D);
     100
     101    // dR = 3600*(Rave - Rmeas) -> Rmeas = Rave - dR / 36000
     102
     103    if (isnan(dR)) continue;
     104    if (isnan(dD)) continue;
     105
     106    icrfobj->Rave[Npts]  = average->R;
     107    icrfobj->Dave[Npts]  = average->D;
     108    icrfobj->dRoff[Npts] = dR;
     109    icrfobj->dRoff[Npts] = dD;
     110    Npts ++;
     111  }
     112  icrfobj->Nicrfobj = Npts;
     113
     114  return icrfobj;
     115}
     116
     117void ICRFobjFree (ICRFobj *icrfobj) {
     118
     119  if (!icrfobj) return;
     120
     121  FREE (icrfobj->Rave);
     122  FREE (icrfobj->Dave);
     123  FREE (icrfobj->dRoff);
     124  FREE (icrfobj->dDoff);
     125  return;
     126}
     127
Note: See TracChangeset for help on using the changeset viewer.