Changeset 37623

Timestamp:

Nov 18, 2014, 9:30:22 AM (12 years ago)

Author:

eugene

Message:

fix up frame correction (Ylm)

Location:

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

Files:

• FrameCorrection.c (modified) (9 diffs)
• bcatalog.c (modified) (2 diffs)

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

@@ -40 +40 @@
   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] = 360.0 / frame->dR[i];
+    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;
@@ -61 +70 @@
   }
 
+  if (frame->dR)   free (frame->dR);
   if (frame->Nra)  free (frame->Nra);
   if (frame->Roff) free (frame->Roff);
@@ -138 +148 @@
     SHterms *dDc = SHtermsInit (Lmax);
 
+    INITTIME;
+
     FrameCorrectionFitSH (catalog, Ncatalog, dRc, dDc);
+    MARKTIME ("done with FitSH: %f sec\n", dtime);
 
     double pltscale = 0.5; // degrees per pixel
@@ -144 +157 @@
   
     FrameCorrectionFromSH (frame, dRc, dDc);
+    MARKTIME ("done with FromSH: %f sec\n", dtime);
 
     // XXX write out an image to represent the correction
@@ -172 +186 @@
   fprintf (stderr, "start Fit %d SH l-modes for %d QSOs\n", dRc->lmax, Nicrf);
 
+  FILE *f = fopen ("sh.dat", "w");
+  myAssert (f, "oops");
+
   int i, j;
   for (i = 0; i < Nicrf; i++) {
 
     int cat, meas, ave;
-    ICRFdata (i, &cat, &meas, &ave);
+    ICRFdata (i, &cat, &ave, &meas);
 
     // I need to generate a collection of values dR,dD(R,D), where dR = average.R - ICRF.R, etc.  
@@ -189 +206 @@
     // record these in arcsec or degree?
     // correct for cos(D) or not?
-    double dR = average->R - measure->R;
-    double dD = average->D - measure->D;
-
+    double dR = 3600.0*(average->R - measure->R);
+    double dD = 3600.0*(average->D - measure->D);
+
+    fprintf (f, "%12.8f %12.8f %7.3f %7.3f\n", average->R, average->D, dR, dD);
     if (isnan(dR)) continue;
     if (isnan(dD)) continue;
@@ -205 +223 @@
       dDc->Fi[j] += dD * SH->Fi[j];
     }
-  }    
+  }
+  fclose (f);
 
   for (j = 0; j < SH->Nterms; j++) {
-    dRc->Fr[j] *= 4.0*M_PI/Nicrf;
-    dRc->Fi[j] *= 4.0*M_PI/Nicrf;
-    dDc->Fr[j] *= 4.0*M_PI/Nicrf;
-    dDc->Fi[j] *= 4.0*M_PI/Nicrf;
+    dRc->Fr[j] *= 4.0*M_PI/(float)Nicrf;
+    dRc->Fi[j] *= 4.0*M_PI/(float)Nicrf;
+    dDc->Fr[j] *= 4.0*M_PI/(float)Nicrf;
+    dDc->Fi[j] *= 4.0*M_PI/(float)Nicrf;
   }
   SHtermsFree (SH);
@@ -283 +302 @@
       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;
@@ -507 +533 @@
 
         int iD = (D + 89.0) / frame->scale;
+        if (iD < 0) continue;
+        if (iD >= frame->Ndec) continue;
+
         int iR = R / frame->dR[iD];
+        if (iR < 0) continue;
+        if (iR >= frame->Nra[iD]) continue;
 
         double dR = frame->Roff[iD][iR];
         double dD = frame->Doff[iD][iR];
-        average->R -= dR;
-        average->D -= dD;
+        average->R -= dR / 3600.0;
+        average->D -= dD / 3600.0;
       }
     }

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

@@ -85 +85 @@
     if (USE_ICRF_CORRECT && (subcatalog[0].average[Naverage].flags & ID_ICRF_QSO)) {
       Nm = 0;
+      int savedICRF = FALSE;
       for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
           
@@ -100 +101 @@
         // we need to save the location of the ICRF QSOs in the database
         if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_ICRF_QSO) {
-          ICRFsave (Nsubcatalog, Naverage, Nmeasure);
-          fprintf (f, "%d %d %d : %f %f : %f %f\n", (int) Nsubcatalog, (int) Naverage, (int) Nmeasure,
-                   catalog[0].average[i].R, catalog[0].average[i].D,
-                   catalog[0].measure[offset].R, catalog[0].measure[offset].D);
-          Nicrf ++;
+          if (!savedICRF) {
+            ICRFsave (Nsubcatalog, Naverage, Nmeasure);
+            fprintf (f, "%d %d %d : %f %f : %f %f\n", (int) Nsubcatalog, (int) Naverage, (int) Nmeasure,
+                     catalog[0].average[i].R, catalog[0].average[i].D,
+                     catalog[0].measure[offset].R, catalog[0].measure[offset].D);
+            Nicrf ++;
+            savedICRF = TRUE;
+          }
         }