Changeset 16798 for branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c

Timestamp:

Mar 4, 2008, 12:18:13 PM (18 years ago)

Author:

eugene

Message:

substantial work to fix FitChips

File:

• branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c (modified) (7 diffs)

branches/eam_branch_20080223/Ohana/src/relastro/src/fitpoly.c

@@ -128 +128 @@
       matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
     }
+
+    // mask the terms not represented by the Coords terms
+    if (ix + iy > fit[0].Norder) {
+      for (j = 0; j < fit[0].Nelems; j++) {
+        matrix[i][j] = 0.0;
+      }
+      vector[i][0] = 0.0;
+      vector[i][1] = 0.0;
+      matrix[i][i] = 1.0;
+    }      
   }
 
@@ -133 +143 @@
     ix = i % fit[0].Nterms;
     iy = i / fit[0].Nterms;
-    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
-             ix, iy, vector[i][0], ix, iy, vector[i][1]);
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }     
 
@@ -142 +152 @@
     ix = i % fit[0].Nterms;
     iy = i / fit[0].Nterms;
-    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
-             ix, iy, vector[i][0], ix, iy, vector[i][1]);
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
   }     
 
@@ -158 +168 @@
 }
 
+void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) {
+  
+  int ix, iy;
+  double xterm, yterm, term;
+  double x, y;
+
+  x = 0.0;
+  y = 0.0;
+
+  xterm = 1;
+  for (ix = 0; ix < fit[0].Nterms; ix++) {
+    yterm = 1;
+    for (iy = 0; iy < fit[0].Nterms; iy++) {
+      term = xterm*yterm;
+      x += fit[0].xfit[ix][iy]*term;
+      y += fit[0].yfit[ix][iy]*term;
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  *x2 = x;
+  *y2 = y;
+}
+
 // Nx, Ny is the number of terms in x and in y
 double **poly2d_dx (double **poly, int Nx, int Ny) {
@@ -209 +243 @@
   for (i = 0; i < Nx; i++) {
     for (j = 0; j < Ny; j++) {
-      out[i][j] = poly[i][j+1];
+      out[i][j] = poly[i][j];
     }
   }
@@ -237 +271 @@
 /* this should only apply to the polynomial, not the projection terms */
 /* compare with psastro supporting code */
-void fit_apply_coords (CoordFit *fit, Coords *coords) {
+CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) {
 
   double c11, c12;
   double c21, c22;
-  double Xo, Yo, R;
+  double Xo, Yo, R1, R2;
   CoordFit *modfit;
 
@@ -256 +290 @@
   modfit = CoordsSetCenter (fit, Xo, Yo);
 
+  /* we do not modify crval1,2: these are kept at the default values */
+
+  // set cdelt1, cdelt2
+  coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]);
+  coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]);
+  R1 = 1 / coords[0].cdelt1;
+  R2 = 1 / coords[0].cdelt2;
+
   // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
-  coords[0].pc1_1 = modfit[0].xfit[1][0];
-  coords[0].pc1_2 = modfit[0].xfit[0][1];
-  coords[0].pc2_1 = modfit[0].yfit[1][0];
-  coords[0].pc2_2 = modfit[0].yfit[0][1];
+  coords[0].pc1_1 = modfit[0].xfit[1][0] * R1;
+  coords[0].pc1_2 = modfit[0].xfit[0][1] * R2;
+  coords[0].pc2_1 = modfit[0].yfit[1][0] * R1;
+  coords[0].pc2_2 = modfit[0].yfit[0][1] * R2;
 
   // set the polyterm elements 
   if (coords->Npolyterms > 1) {
-    coords[0].polyterms[0][0] = modfit[0].xfit[2][0];
-    coords[0].polyterms[1][0] = modfit[0].xfit[1][1];
-    coords[0].polyterms[2][0] = modfit[0].xfit[0][2];
-
-    coords[0].polyterms[0][1] = modfit[0].yfit[2][0];
-    coords[0].polyterms[1][1] = modfit[0].yfit[1][1];
-    coords[0].polyterms[2][1] = modfit[0].yfit[0][2];
+    coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1;
+    coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2;
+    coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2;
+
+    coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1;
+    coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2;
+    coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2;
   }
 
   // I need to validate Norder
   if (coords->Npolyterms > 2) {
-    coords[0].polyterms[3][0] = modfit[0].xfit[3][0];
-    coords[0].polyterms[4][0] = modfit[0].xfit[2][1];
-    coords[0].polyterms[5][0] = modfit[0].xfit[1][2];
-    coords[0].polyterms[6][0] = modfit[0].xfit[0][3];
-
-    coords[0].polyterms[3][1] = modfit[0].yfit[3][0];
-    coords[0].polyterms[4][1] = modfit[0].yfit[2][1];
-    coords[0].polyterms[5][1] = modfit[0].yfit[1][2];
-    coords[0].polyterms[6][1] = modfit[0].yfit[0][3];
-  }
-
-  /* we do not modify crval1,2: these are kept at the default values */
-
-  /* normalize pc11,etc */
-  c11 = coords[0].pc1_1;
-  c21 = coords[0].pc1_1;
-  c12 = coords[0].pc1_1;
-  c22 = coords[0].pc1_1;
-  coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
-  R = 1 / coords[0].cdelt1;
-
-  coords[0].pc1_1  = c11*R;
-  coords[0].pc2_1  = c21*R;
-  coords[0].pc1_2  = c12*R;
-  coords[0].pc2_2  = c22*R;
+    coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2;
+
+    coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2;
+  }
 
   fit_free (modfit);