Changeset 38486

Timestamp:

Jun 17, 2015, 5:33:20 AM (11 years ago)

Author:

eugene

Message:

loadgalphot runs and does not crash (but does not yet work)

Location:

branches/eam_branches/ipp-20150616/Ohana/src/addstar

Files:

• include/loadgalphot.h (modified) (1 diff)
• src/loadgalphot_catalog.c (modified) (1 diff)
• src/loadgalphot_fit2d.c (modified) (5 diffs)
• src/loadgalphot_readstars.c (modified) (8 diffs)

branches/eam_branches/ipp-20150616/Ohana/src/addstar/include/loadgalphot.h

@@ -48 +48 @@
 void fit2d_free (Fit2D *fit);
 int fit2d (Fit2D *fit, float *xval, float *yval, float *zval, float *zfit, char *mask, int Npts);
+int fit2d_reset (Fit2D *fit);

branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_catalog.c

@@ -7 +7 @@
 
   // now we have all of the loaded stars in this catalog
+  dvo_catalog_init (&catalog, TRUE);
   catalog.filename = filename;
   catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data

branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_fit2d.c

@@ -7 +7 @@
 
   Fit2D *fit;
-  ALLOCATE (fit, Fit2D, 1);
+  ALLOCATE_ZERO (fit, Fit2D, 1);
 
   // allocate static arrays
+  fit->order = order;
   fit->nterm = order + 1;
   fit->wterm = fit->nterm*(fit->nterm + 1)/2;
@@ -19 +20 @@
   ALLOCATE (fit->c, double *, fit->wterm);
   for (i = 0; i < fit->wterm; i++) {
-    ALLOCATE (fit->c[i], double, fit->wterm);
-    ALLOCATE (fit->b[i], double, 1);
+    ALLOCATE_ZERO (fit->c[i], double, fit->wterm);
+    ALLOCATE_ZERO (fit->b[i], double, 1);
   }
   for (i = 0; i < fit->mterm; i++) {
-    ALLOCATE (fit->s[i], double, fit->mterm);
-  }
+    ALLOCATE_ZERO (fit->s[i], double, fit->mterm);
+  }
+  ALLOCATE_ZERO (fit->Cii, double, fit->wterm);
   return fit;
+}
+
+int fit2d_reset (Fit2D *fit) {
+
+  int i;
+
+  for (i = 0; i < fit->wterm; i++) {
+    memset (fit->c[i], 0, fit->wterm*sizeof(double));
+    memset (fit->b[i], 0, sizeof(double));
+  }
+  for (i = 0; i < fit->mterm; i++) {
+    memset (fit->s[i], 0, fit->mterm*sizeof(double));
+  }
+  return TRUE;
 }
 
@@ -49 +65 @@
 int fit2d (Fit2D *fit, float *xval, float *yval, float *zval, float *zfit, char *mask, int Npts) {
   
-  int N, i, j, nx, ny, K, k, n, Npt;
+  int N, i, j, nx, ny, K, k, n, Npt, status;
 
   double mean, sigma, maxsigma;
@@ -59 +75 @@
 
   for (N = 0; N < ClipNiter; N++) {
-    /* init registers for current pass */
-    for (i = 0; i < fit->wterm; i++) {
-      memset (fit->c[i], 0, fit->wterm*sizeof(double));
-      memset (fit->b[i], 0, sizeof(double));
-    }
-    for (i = 0; i < fit->mterm; i++) {
-      memset (fit->s[i], 0, fit->mterm*sizeof(double));
-    }
+    fit2d_reset (fit);
 
     // pointers which loop over datapoints
@@ -120 +129 @@
 
     // invert the c,b matrix equation
-    dgaussjordan (fit->c, fit->b, fit->wterm, 1);
+    status = dgaussjordan (fit->c, fit->b, fit->wterm, 1);
+    if (!status) return FALSE;
 
     /* the b[][0] terms are in the following order:

branches/eam_branches/ipp-20150616/Ohana/src/addstar/src/loadgalphot_readstars.c

@@ -93 +93 @@
   Fit2D *fit = fit2d_init (2);
 
+  int Nbad = 0;
   for (i = 0; i < Nrow; i++) {
 
@@ -117 +118 @@
     // I need to find the minimum position (interpolated) in this 2D space
 
-    FitChisqMinimum (fit, &stars[i].galphot, &GAL_CHISQ[Nsample*i], &GAL_FLUX[Nsample*i], &GAL_FLUX_ERR[Nsample*i], Nsample, FR_MAJOR_MIN[i], FR_MAJOR_MAX[i], FR_MAJOR_DEL[i], FR_MINOR_MIN[i], FR_MINOR_MAX[i], FR_MINOR_DEL[i]);
+    if (!FitChisqMinimum (fit, &stars[i].galphot, &GAL_CHISQ[Nsample*i], &GAL_FLUX[Nsample*i], &GAL_FLUX_ERR[Nsample*i], Nsample, FR_MAJOR_MIN[i], FR_MAJOR_MAX[i], FR_MAJOR_DEL[i], FR_MINOR_MIN[i], FR_MINOR_MAX[i], FR_MINOR_DEL[i])) {
+      Nbad ++;
+    }
   }
 
@@ -129 +132 @@
   fit2d_free (fit);
 
+  fprintf (stderr, "Nbad: %d of %d\n", Nbad, (int) Nrow);
+
   *nstars = Nrow;
   return (stars);
@@ -155 +160 @@
   // validate the square dimensions??
 
-  int Nx = (MajorMax - MajorMin) / MajorDel;
-  int Ny = (MinorMax - MinorMin) / MinorDel;
+  int Nx = nearbyint((MajorMax - MajorMin) / MajorDel) + 1;
+  int Ny = nearbyint((MinorMax - MinorMin) / MinorDel) + 1;
   myAssert (Nx*Ny == Npts, "inconsistent grid");
   
@@ -193 +198 @@
     galphot->majorAxisErr = MajorDel;
     galphot->minorAxisErr = MinorDel;
-    return TRUE;
+    return FALSE;
   }
 
@@ -217 +222 @@
 
   // fit the inner 9 points
-  fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts);
+  if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) {
+    // raise a flag?
+    galphot->majorAxis = Xpt[iMin];
+    galphot->minorAxis = Ypt[iMin];
+    galphot->majorAxisErr = MajorDel;
+    galphot->minorAxisErr = MinorDel;
+    return FALSE;
+  }
 
   // exclude any extreme outliers
@@ -230 +242 @@
   
   // re-fit all except the most extreme set
-  fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts);
+  if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) {
+    // raise a flag?
+    galphot->majorAxis = Xpt[iMin];
+    galphot->minorAxis = Ypt[iMin];
+    galphot->majorAxisErr = MajorDel;
+    galphot->minorAxisErr = MinorDel;
+    return FALSE;
+  }
 
   // get X,Y for the chisq min from the fit:
@@ -269 +288 @@
   galphot->magErr = -2.5*log10(fluxErr[iMin]); // correct for exptime?
   galphot->chisq  = chisqMin;
-
   
   return TRUE;