Changeset 37562

Timestamp:

Nov 6, 2014, 8:17:04 PM (12 years ago)

Author:

eugene

Message:

add frame correction

Location:

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

Files:

• include/relastro.h (modified) (1 diff)
• src/FrameCorrection.c (modified) (2 diffs)
• src/UpdateChips.c (modified) (1 diff)

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

@@ -575 +575 @@
 FrameCorrectionType *FrameCorrectionInit (double scale);
 void FrameCorrectionFree (FrameCorrectionType *frame);
+
 int FrameCorrection (Catalog *catalog, int Ncatalog);
-int FrameCorrectionFit (Catalog *catalog, int Ncatalog, SHterms *dR, SHterms *dD);
+int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dR, SHterms *dD);
 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);

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

@@ -75 +75 @@
   // go from ICRF object to the catalog entry.  
 
-  // Lmax in recipe
-  int Lmax = 20;
-
-  // dR and dD will carry the fit coefficients in RA and DEC
-  SHterms *dRc = SHtermsInit (Lmax);
-  SHterms *dDc = SHtermsInit (Lmax);
-
-  FrameCorrectionFit (catalog, Ncatalog, dRc, dDc);
-
-  double pltscale = 0.5; // degrees per pixel
-  FrameCorrectionType *frame = FrameCorrectionInit (pltscale);
+  // we can have 2 kinds of corrections:
+  // 'frame' is a map in spherical coords of the corrections based on spherical harmonics
+  // 'map' is a correction based on a local projection to a linear coordinate system.  I
+  // am going to use this for the pole, but also for SAS for testing.
+  // We need to be sure only one of the two is applied. For the pole, the boundary is a
+  // fixed line of DEC; for the SAS, the boundary is the projection 'image'
+
+  Coords coords;
+  AstromOffsetMap *map = NULL;
+  FrameCorrectionType *frame = NULL;
+
+# define LOCAL_FRAME 1
+# define SH_FRAME 0
+
+  int POLE = FALSE;
+
+  if (LOCAL_FRAME) {
+    // for the local frame correction, we are going to convert R,D into 
+    // X,Y linear coordinates in an 'image' centered on the field center. 
+    // we can then use an AstromOffsetMap structure to carry the correction values
+
+    // what parameters define a local frame correction?
+    // * Ro, Do -- map/projection center
+    // * scale  -- degrees / pixel
+    // * Nx, Ny -- size of map in pixels
+    // * Rmin,Rmax,Dmin,Dmax
+
+    // example parameters for SAS: 320 - 340, -10 - +10
+    double Ro = 330.0;
+    double Do =   0.0;
+
+    // double Ro =  0.0;
+    // double Do = 90.0;
+    // POLE = TRUE;
+
+    double scale = 2.0; // degrees per patch
+    int Nx = 10;
+    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);
+  } 
+
+  if (SH_FRAME) {
+    // Lmax in recipe
+    int Lmax = 20;
+
+    // dR and dD will carry the fit coefficients in RA and DEC
+    SHterms *dRc = SHtermsInit (Lmax);
+    SHterms *dDc = SHtermsInit (Lmax);
+
+    FrameCorrectionFitSH (catalog, Ncatalog, dRc, dDc);
+
+    double pltscale = 0.5; // degrees per pixel
+    frame = FrameCorrectionInit (pltscale);
   
-  FrameCorrectionFromSH (frame, dRc, dDc);
-
-  // XXX write out an image to represent the correction
+    FrameCorrectionFromSH (frame, dRc, dDc);
+
+    // XXX write out an image to represent the correction
+
+    SHtermsFree (dRc);
+    SHtermsFree (dDc);
+  }
 
   // Now apply the correction to all of the average.R,D values
-
-  int i, j;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      
-      double R = catalog[i].average[j].R;
-      double D = catalog[i].average[j].D;
-
-      // special handling for the polar regions...
-
-      int iD = (D + 89.0) / frame->scale;
-      int iR = R / frame->dR[iD];
-
-      double dR = frame->Roff[iD][iR];
-      double dD = frame->Doff[iD][iR];
-
-      catalog[i].average[j].R -= dR;
-      catalog[i].average[j].D -= dD;
-    }
-  }
+  FrameCorrectionApply (catalog, Ncatalog, frame, map, &coords, POLE);
+
   FrameCorrectionFree(frame);
-  return TRUE;
-}
-
-
-int FrameCorrectionFit (Catalog *catalog, int Ncatalog, SHterms *dRc, SHterms *dDc) {
+  AstromOffsetMapFree (map);
+
+  return TRUE;
+}
+
+int FrameCorrectionFitSH (Catalog *catalog, int Ncatalog, SHterms *dRc, SHterms *dDc) {
 
   myAssert (dRc->lmax == dDc->lmax, "dR and dD must match\n");
@@ -196 +237 @@
   return TRUE;
 }
+
+int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, Coords *coords) {
+
+  double Xave, Yave, Xmeas, Ymeas;
+  float *X, *Y, *dX, *dY;
+
+  int Npts = 0;
+  int NPTS = 100;
+
+  ALLOCATE (X, float, NPTS);
+  ALLOCATE (Y, float, NPTS);
+  ALLOCATE (dX, float, NPTS);
+  ALLOCATE (dY, float, NPTS);
+
+  int Nicrf = ICRFmax();
+
+  int i;
+  for (i = 0; i < Nicrf; i++) {
+
+    int cat, meas, ave;
+    ICRFdata (i, &cat, &meas, &ave);
+
+    Average *average = &catalog[cat].average[ave];
+    Measure *measure = &catalog[cat].measure[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);
+    if (!status) continue;
+    if (Xave < 0.0) continue;
+    if (Xave > map->Nx) continue;
+    if (Yave < 0.0) continue;
+    if (Yave > map->Ny) continue;
+
+    RD_to_XY (&Xmeas, &Ymeas, measure->R, measure->D, coords);
+
+    // record these in arcsec or degree?
+    // correct for cos(D) or not?
+    X[Npts] = Xave;
+    Y[Npts] = Yave;
+    dX[Npts] = Xave - Xmeas;
+    dY[Npts] = Yave - Ymeas;
+    Npts ++;
+    if (Npts == NPTS) {
+      NPTS += 100;
+      REALLOCATE (X, float, NPTS);
+      REALLOCATE (Y, float, NPTS);
+      REALLOCATE (dX, float, NPTS);
+      REALLOCATE (dY, float, NPTS);
+    }
+  }    
+
+  AstromOffsetMapFit (map, X, Y, dX, Npts, TRUE);
+  AstromOffsetMapFit (map, X, Y, dY, Npts, FALSE);
+
+  return TRUE;
+}
+
+int FrameCorrectionApply (Catalog *catalog, int Ncatalog, FrameCorrectionType *frame, AstromOffsetMap *map, Coords *coords, int POLE) {
+
+  double Xave, Yave;
+
+  int i, j;
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      
+      Average *average = &catalog[i].average[j];
+
+      double R = average->R;
+      double D = average->D;
+
+      myAssert (coords, "no local frame coords defined");
+
+      // 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);
+      } else {
+        useLocal = (useLocal && (Xave > 0.0));
+        useLocal = (useLocal && (Xave < map->Nx));
+        useLocal = (useLocal && (Yave > 0.0));
+        useLocal = (useLocal && (Yave < map->Ny));
+      }
+
+      if (useLocal) {
+        myAssert (map, "no local frame map defined");
+        double dX = AstromOffsetMapValue (map, Xave, Yave, TRUE);
+        double dY = AstromOffsetMapValue (map, Xave, Yave, FALSE);
+        double Xmeas = Xave - dX;
+        double Ymeas = Yave - dY;
+        XY_to_RD (&average->R, &average->D, Xmeas, Ymeas, coords);
+      } else {
+        myAssert (frame, "no frame correction defined");
+
+        int iD = (D + 89.0) / frame->scale;
+        int iR = R / frame->dR[iD];
+
+        double dR = frame->Roff[iD][iR];
+        double dD = frame->Doff[iD][iR];
+        average->R -= dR;
+        average->D -= dD;
+      }
+    }
+  }
+  return TRUE;
+}
+

branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/UpdateChips.c

@@ -322 +322 @@
     } 
 
+    // apply the modified R,D back to the measures
+    setImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, raw, Nraw, MODE_MOSAIC);
+
     saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
     threadinfo->mode[i] = 3;