Changeset 37604

Timestamp:

Nov 14, 2014, 8:42:00 PM (12 years ago)

Author:

eugene

Message:

testing ICRF and Galaxy Model correctios

Location:

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

Files:

• include/relastro.h (modified) (1 diff)
• src/FitChip.c (modified) (1 diff)
• src/FrameCorrection.c (modified) (3 diffs)
• src/ICRF.c (modified) (1 diff)
• src/ImageOps.c (modified) (1 diff)
• src/UpdateObjects.c (modified) (1 diff)
• src/bcatalog.c (modified) (4 diffs)
• src/initialize.c (modified) (1 diff)
• src/load_catalogs.c (modified) (1 diff)

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

@@ -284 +284 @@
 int           args                PROTO((int argc, char **argv));
 int           args_client         PROTO((int argc, char **argv));
-int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
+int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog, int Nsubcatalog));
 void          clean_images        PROTO((void));
 void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));

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

@@ -72 +72 @@
       return FALSE;
     }
-    // fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
+    // XXX WATCH OUT HERE!!
+    if (!CHIPMAP) order_use = MIN(order_use, 1);
+    fprintf (stderr, "using %d for %s\n", order_use, image[0].name);
 
     // when fitting the map, first fit a linear model (below? change Npolyterms to -1)

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

@@ -180 +180 @@
 
     Average *average = &catalog[cat].average[ave];
-    Measure *measure = &catalog[cat].measure[meas]; // MeasureTiny?
+    MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
 
     // record these in arcsec or degree?
@@ -240 +240 @@
 int FrameCorrectionFitLocal (Catalog *catalog, int Ncatalog, AstromOffsetMap *map, 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;
+  }
+
+  double *Rave, *Dave, *Rmeas, *Dmeas;
+  float *dXf, *dYf, *dPos;
+
   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;
+  ALLOCATE (Rave, double, Nicrf);
+  ALLOCATE (Dave, double, Nicrf);
+  ALLOCATE (Rmeas, double, Nicrf);
+  ALLOCATE (Dmeas, double, Nicrf);
+
+  ALLOCATE (X,  float, Nicrf);
+  ALLOCATE (Y,  float, Nicrf);
+  ALLOCATE (dX, float, Nicrf);
+  ALLOCATE (dY, float, Nicrf);
+  ALLOCATE (dXf, float, Nicrf);
+  ALLOCATE (dYf, float, Nicrf);
+
+  ALLOCATE (dPos, float, Nicrf);
+
+  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, &meas, &ave);
+    int cat, ave, meas;
+    ICRFdata (i, &cat, &ave, &meas);
 
     Average *average = &catalog[cat].average[ave];
-    Measure *measure = &catalog[cat].measure[meas]; // MeasureTiny?
+    MeasureTiny *measure = &catalog[cat].measureT[meas]; // MeasureTiny?
 
     // this local correction is defined for (Rmin < R < Rmax, Dmin < D < Dmax)
@@ -278 +297 @@
     dX[Npts] = Xave - Xmeas;
     dY[Npts] = Yave - Ymeas;
+    Rave[Npts] = average->R;
+    Dave[Npts] = average->D;
+    Rmeas[Npts] = measure->R;
+    Dmeas[Npts] = measure->D;
     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);
+  }
+
+  float *Xfit, *Yfit, *dXfit, *dYfit;
+  ALLOCATE (Xfit, float, Npts);
+  ALLOCATE (Yfit, float, Npts);
+  ALLOCATE (dXfit, float, Npts);
+  ALLOCATE (dYfit, float, Npts);
+
+  int Niter;
+  for (Niter = 0; Niter < 3; Niter ++) {
+    int Nkeep = 0;
+    for (i = 0; i < Npts; i++) {
+      if (mask[i]) continue;
+      
+      Xfit[Nkeep] = X[i];
+      Yfit[Nkeep] = Y[i];
+      dXfit[Nkeep] = dX[i];
+      dYfit[Nkeep] = dY[i];
+      Nkeep ++;
+    }
+
+    AstromOffsetMapFit (map, Xfit, Yfit, dXfit, Nkeep, TRUE);
+    AstromOffsetMapFit (map, Xfit, Yfit, dYfit, Nkeep, FALSE);
+
+    for (i = 0; i < Npts; i++) {
+      dXf[i] = AstromOffsetMapValue (map, X[i], Y[i], TRUE);
+      dYf[i] = AstromOffsetMapValue (map, X[i], Y[i], FALSE);
+    }
+
+    // generate a histogram of the dPos = sqrt((dXf - dX)^2 + (dYf - dY)^2) values
+    Nkeep = 0;
+    for (i = 0; i < Npts; i++) {
+      if (mask[i]) continue;
+      float dP = hypot((dXf[i] - dX[i]), (dYf[i] - dY[i]));
+      if (isnan(dP)) continue; // not all points are in the map
+      dPos[Nkeep] = dP;
+      Nkeep ++;
+    }
+    fsort (dPos, Nkeep);
+    
+    float dPmax = 1.5*dPos[(int)(0.95*Nkeep)];
+
+    // mask points based on dPmax
+    Nkeep = 0;
+    for (i = 0; i < Npts; i++) {
+      float dP = hypot((dXf[i] - dX[i]), (dYf[i] - dY[i]));
+      if (dP > dPmax) mask[i] = 1;
+      if (isnan(dP)) mask[i] = 1;
+      if (mask[i]) Nkeep++;
+    }
+    fprintf (stderr, "keeping %d of %d points\n", Nkeep, Npts);
+  }
+  free (Xfit);
+  free (Yfit);
+  free (dXfit);
+  free (dYfit);
+
+  // *** save the icrf points and fit residuals
+  { 
+    FILE *f = fopen ("map.dat", "w");
+
+    for (i = 0; i < Npts; i++) {
+      fprintf (f, "%3d %12.8f %12.8f  %12.8f %12.8f  %f %f : %f %f : %f %f : %d\n", i, Rave[i], Dave[i], Rmeas[i], Dmeas[i], X[i], Y[i], dX[i], dY[i], dXf[i], dYf[i], mask[i]);
+    }    
+    fclose (f);
+  }
+
+  // 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;

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

@@ -18 +18 @@
 int ICRFsave (int cat, int ave, int meas) {
   
+  fprintf (stderr, "ICRF: %d %d %d\n", cat, ave, meas);
+
   ICRFtoCatalog[Nicrf] = cat;
   ICRFtoAverage[Nicrf] = ave;

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

@@ -810 +810 @@
     // XXX subtract off dR,dD GAL here
     if (USE_GALAXY_MODEL) {
-      ref[i].R += measure[0].RoffGAL;
-      ref[i].D += measure[0].DoffGAL;
+      ref[i].R += measure[0].RoffGAL / 3600.0;
+      ref[i].D += measure[0].DoffGAL / 3600.0;
     }
 

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

@@ -210 +210 @@
 
     if (USE_GALAXY_MODEL) {
-      Ri -= measure[0].RoffGAL;
-      Di -= measure[0].DoffGAL;
+      Ri -= measure[k].RoffGAL / 3600.0;
+      Di -= measure[k].DoffGAL / 3600.0;
     }
 

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

@@ -19 +19 @@
 int LimitDensityCatalog_ByNmeasureGrid (Catalog *subcatalog, Catalog *oldcatalog);
 
-int bcatalog (Catalog *subcatalog, Catalog *catalog) {
+int bcatalog (Catalog *subcatalog, Catalog *catalog, int Nsubcatalog) {
 
   off_t i, j, offset;
@@ -62 +62 @@
 
   int Nicrf = 0;
+
+  FILE *f = fopen ("test.icrf.dat", "a");
 
   /* exclude stars not in range or with too few measurements */
@@ -98 +100 @@
         // we need to save the location of the ICRF QSOs in the database
         if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_ICRF_QSO) {
-          ICRFsave (i, Naverage, Nmeasure);
+          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 ++;
         }
@@ -224 +229 @@
   subcatalog[0].catID    = catalog[0].catID;
   assert (Nsecfilt == catalog[0].Nsecfilt);
+
+  fclose (f);
 
 // limit the total number of stars in the catalog

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

@@ -63 +63 @@
   initstats (STATMODE);
 
+  if (USE_ICRF_CORRECT) ICRFinit();
+
   /* XXX drop irrelevant entries */
   if (SHOW_PARAMS) {

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

@@ -68 +68 @@
       // results are in Average, Secfilt, and MeasureTiny
       // Ncat tracks the actually used catalogs
-      bcatalog (&catalog[Ncat], &tcatalog);
+      bcatalog (&catalog[Ncat], &tcatalog, Ncat);
       dvo_catalog_unlock (&tcatalog);
       dvo_catalog_free (&tcatalog);