Changeset 37675

Timestamp:

Nov 26, 2014, 6:23:50 AM (12 years ago)

Author:

eugene

Message:

add outlier clipping to Ylm fit; l = 5 (still hard-wired); Ylm used to fit dR*cos(dec), not dR

File:

• branches/eam_branches/ipp-20140904/Ohana/src/relastro/src/FrameCorrection.c (modified) (5 diffs)

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

@@ -117 +117 @@
     double Ro = 330.0;
     double Do =   0.0;
+    double scale = 2.0; // degrees per patch
+    int Nx = 10;
+    int Ny = 10;
     if (USE_ICRF_POLE) {
        Ro =  0.0;
        Do = 90.0;
-    }
-
-    double scale = 2.0; // degrees per patch
-    int Nx = 10;
-    int Ny = 10;
+       scale = 5.0; // degrees per patch
+       Nx = 4;
+       Ny = 4;
+    }
+
     
     Coords *coords = NULL;
@@ -147 +150 @@
   if (USE_ICRF_SHFIT) {
     // Lmax in recipe
-    int Lmax = 20;
+    int Lmax = 5;
 
     // dR and dD will carry the fit coefficients in RA and DEC
@@ -191 +194 @@
   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 
-  SHterms *Ylm = SHtermsInit (dRc->lmax);
-
-  static int version = 0;
-  char filename[1024];
-  snprintf (filename, 1024, "%s/sh.%02d.dat", CATDIR, version);
-  FILE *f = fopen (filename, "w");
-  myAssert (f, "oops");
-  version ++;
-
-  for (i = 0; i < icrfobj->Nicrfobj; i++) {
+  double *dPosObs, *dPosSrt = NULL;
+  ALLOCATE (dPosObs, double, icrfobj->Nicrfobj);
+  ALLOCATE (dPosSrt, double, icrfobj->Nicrfobj);
+
+  int Niter;
+  for (Niter = 0; Niter < 3; Niter ++) {
+
+    SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dRoff, mask, icrfobj->Nicrfobj, dRc);
+    SHfitWithMask (icrfobj->Rave, icrfobj->Dave, icrfobj->dDoff, mask, icrfobj->Nicrfobj, dDc);
+
+    // the bit below is for testing
+
+    // allocate an SHterms structure to hold the Ylm values 
+    SHterms *Ylm = SHtermsInit (dRc->lmax);
+
+    static int version = 0;
+    char filename[1024];
+    snprintf (filename, 1024, "%s/sh.%02d.dat", CATDIR, version);
+    FILE *f = fopen (filename, "w");
+    myAssert (f, "oops");
+    version ++;
+
+    int Nkeep = 0;
+    for (i = 0; i < icrfobj->Nicrfobj; i++) {
     
-    SHtermsForRD (Ylm, icrfobj->Rave[i], icrfobj->Dave[i]);
+      SHtermsForRD (Ylm, icrfobj->Rave[i], icrfobj->Dave[i]);
       
-    double dRfit = 0.0;
-    double dDfit = 0.0;
-
-    int j;
-    for (j = 0; j < Ylm->Nterms; j++) {
-      dRfit += dRc->Fr[j]*Ylm->Fr[j] + dRc->Fi[j]*Ylm->Fi[j];
-      dDfit += dDc->Fr[j]*Ylm->Fr[j] + dDc->Fi[j]*Ylm->Fi[j];
-    }
-    fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f\n", 
-             icrfobj->Rave[i], icrfobj->Dave[i], icrfobj->dRoff[i], icrfobj->dDoff[i], dRfit, dDfit);
-  }
-  fclose (f);
+      double dRfit = 0.0;
+      double dDfit = 0.0;
+
+      int j;
+      for (j = 0; j < Ylm->Nterms; j++) {
+        dRfit += dRc->Fr[j]*Ylm->Fr[j] + dRc->Fi[j]*Ylm->Fi[j];
+        dDfit += dDc->Fr[j]*Ylm->Fr[j] + dDc->Fi[j]*Ylm->Fi[j];
+        
+      }
+      fprintf (f, "%12.8f %12.8f %7.3f %7.3f : %7.3f %7.3f : %d\n", 
+               icrfobj->Rave[i], icrfobj->Dave[i], icrfobj->dRoff[i], icrfobj->dDoff[i], dRfit, dDfit, mask[i]);
+      dPosObs[i] = hypot(icrfobj->dRoff[i] - dRfit, icrfobj->dDoff[i] - dDfit);
+      if (mask[i]) continue;
+      if (isnan(dPosObs[i])) continue;
+      dPosSrt[Nkeep] = dPosObs[i];
+      Nkeep ++;
+    }
+    fclose (f);
+    
+    dsort (dPosSrt, Nkeep);
+    double dPmax = 1.5*dPosSrt[(int)(0.95*Nkeep)];
+
+    Nkeep = 0;
+    for (i = 0; i < icrfobj->Nicrfobj; i++) {
+      if (dPosObs[i] > dPmax) mask[i] = 1;
+      if (isnan(dPosObs[i])) mask[i] = 1;
+      if (!mask[i]) Nkeep++;
+    }
+    fprintf (stderr, "keeping %d of %d points\n", Nkeep, icrfobj->Nicrfobj);
+  }
 
   free (mask);
+  free (dPosObs);
+  free (dPosSrt);
 
   return TRUE;
@@ -406 +437 @@
     if (Yave > map->Ny) continue;
 
-    double Rm = icrfobj->Rave[i] - icrfobj->dRoff[i] / 3600.0;
+    double Rm = icrfobj->Rave[i] - icrfobj->dRoff[i] / 3600.0 / cos(icrfobj->Rave[i]*RAD_DEG);
     double Dm = icrfobj->Dave[i] - icrfobj->dDoff[i] / 3600.0;
     RD_to_XY (&Xmeas, &Ymeas, Rm, Dm, coords);
@@ -572 +603 @@
         }
 
-        average->R -= dR / 3600.0;
+        average->R -= dR / 3600.0 / cos(average->R*RAD_DEG);
         average->D -= dD / 3600.0;
       }