Changeset 32695 for trunk/Ohana/src/relastro

Timestamp:

Nov 17, 2011, 2:17:35 PM (15 years ago)

Author:

eugene

Message:

fix psphotStack 2nd pass (subtract 1st pass sources from the correct image); fix PS color plotting; plug a leak in psphotStack

Location:

trunk

Files:

• . (modified) (1 prop)
• Ohana/src/relastro/include/relastro.h (modified) (1 diff)
• Ohana/src/relastro/src/FitPM.c (modified) (2 diffs)
• Ohana/src/relastro/src/FitPMandPar.c (modified) (2 diffs)
• Ohana/src/relastro/src/ImageOps.c (modified) (1 diff)
• Ohana/src/relastro/src/UpdateObjects.c (modified) (5 diffs)
• Ohana/src/relastro/src/args.c (modified) (2 diffs)

trunk/Ohana/src/relastro/include/relastro.h

@@ -317 +317 @@
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
 int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean);
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB);
 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
-int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB);
 
 Mosaic *getMosaicForImage (off_t N);

trunk/Ohana/src/relastro/src/FitPM.c

@@ -2 +2 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB) {
 
   int i;
@@ -77 +77 @@
     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
   }
   fit[0].Nfit = Npts;

trunk/Ohana/src/relastro/src/FitPMandPar.c

@@ -2 +2 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB) {
 
   int i;
@@ -99 +99 @@
   chisq = 0.0;
   for (i = 0; i < Npts; i++) {
-    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].dp*pR[i];
-    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].dp*pD[i];
+    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
+    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+
   }
   fit[0].Nfit = Npts;

trunk/Ohana/src/relastro/src/ImageOps.c

@@ -913 +913 @@
   if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
   if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
+  if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks?
   
   /* select measurements by photcode, or equiv photcode, if specified */

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -84 +84 @@
       /* calculate the average value of R,D for a single star */
 
+      XVERB = FALSE;
+
       // skip objects which are known to be problematic
       // XXX include this code or not?
@@ -196 +198 @@
       }
       
-      XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
+      // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
+      // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
 
       // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -211 +216 @@
         }         
 
-        FitPM (&fitPM, X, dX, Y, dY, T, N);
-
-        if (XVERB) fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
+        FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
+
+        if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
 
         // project Ro, Do back to RA,DEC
@@ -230 +235 @@
           ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
         }
-        FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N);
+        FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
         XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
         catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
@@ -275 +280 @@
         break;
       }
-      if (XVERB) fprintf (stderr, "A%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
+      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
                           catalog[i].average[j].R, 
                           catalog[i].average[j].D, 

trunk/Ohana/src/relastro/src/args.c

@@ -14 +14 @@
   FIT_MODE = FIT_AVERAGE;
 
+  OBJ_ID_SRC = OBJ_ID_DST = 0;
+  CAT_ID_SRC = CAT_ID_DST = 0;
+
   if ((N = get_argument (argc, argv, "-merge-source"))) {
     if (N > argc - 6) usage_merge_source();
@@ -40 +43 @@
     remove_argument (N, &argc, argv);
     FIT_TARGET = TARGET_OBJECTS;
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj1"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj2"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
   }