Changeset 40659

Timestamp:

Mar 29, 2019, 4:23:48 PM (7 years ago)

Author:

eugene

Message:

fix some fitting details (slit tilt interpolation)

Location:

trunk/Ohana/src/opihi/mana

Files:

• deimos_fitslit.c (modified) (1 diff)
• deimos_getobj.c (modified) (5 diffs)
• deimos_mkmodel.c (modified) (1 diff)
• deimos_mkobj.c (modified) (3 diffs)
• deimos_mkslit.c (modified) (1 diff)

trunk/Ohana/src/opihi/mana/deimos_fitslit.c

@@ -32 +32 @@
   int Ny = wind[0].matrix.Naxis[1];
   
-  if (Nx != slit[0].matrix.Naxis[0]) return FALSE;
-  if (Ny != slit[0].matrix.Naxis[1]) return FALSE;
+  if (Nx != slit[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
+  if (Ny != slit[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
 
   ResetVector (flux, OPIHI_FLT, Ny);

trunk/Ohana/src/opihi/mana/deimos_getobj.c

@@ -63 +63 @@
 
   float stilt = 0.0; // angle of the slit
+
+  int SaveSegments = FALSE;
+  if ((N = get_argument (argc, argv, "-save-segments"))) {
+    remove_argument (N, &argc, argv);
+    SaveSegments = TRUE;
+  } 
 
   // Input parameters:
@@ -194 +200 @@
   float dxMax = 0; // maximum distance between sky regions
 
+  Vector *skyLsV = NULL, *skyLeV = NULL, *skyLfV = NULL, *skyRsV = NULL, *skyReV = NULL, *skyRfV = NULL;
+  Vector *bckLsV = NULL, *bckLeV = NULL, *bckLfV = NULL, *bckRsV = NULL, *bckReV = NULL, *bckRfV = NULL;
+  if (SaveSegments) {
+    skyLsV = SelectVector ("skyLs", ANYVECTOR, TRUE);  ResetVector (skyLsV, OPIHI_FLT, Ny);
+    skyLeV = SelectVector ("skyLe", ANYVECTOR, TRUE);  ResetVector (skyLeV, OPIHI_FLT, Ny);
+    skyLfV = SelectVector ("skyLf", ANYVECTOR, TRUE);  ResetVector (skyLfV, OPIHI_FLT, Ny);
+    skyRsV = SelectVector ("skyRs", ANYVECTOR, TRUE);  ResetVector (skyRsV, OPIHI_FLT, Ny);
+    skyReV = SelectVector ("skyRe", ANYVECTOR, TRUE);  ResetVector (skyReV, OPIHI_FLT, Ny);
+    skyRfV = SelectVector ("skyRf", ANYVECTOR, TRUE);  ResetVector (skyRfV, OPIHI_FLT, Ny);
+    bckLsV = SelectVector ("bckLs", ANYVECTOR, TRUE);  ResetVector (bckLsV, OPIHI_FLT, Ny);
+    bckLeV = SelectVector ("bckLe", ANYVECTOR, TRUE);  ResetVector (bckLeV, OPIHI_FLT, Ny);
+    bckLfV = SelectVector ("bckLf", ANYVECTOR, TRUE);  ResetVector (bckLfV, OPIHI_FLT, Ny);
+    bckRsV = SelectVector ("bckRs", ANYVECTOR, TRUE);  ResetVector (bckRsV, OPIHI_FLT, Ny);
+    bckReV = SelectVector ("bckRe", ANYVECTOR, TRUE);  ResetVector (bckReV, OPIHI_FLT, Ny);
+    bckRfV = SelectVector ("bckRf", ANYVECTOR, TRUE);  ResetVector (bckRfV, OPIHI_FLT, Ny);
+  }
+
   // loop over the y positions
   for (int iy = 0; iy < Ny; iy++) {
@@ -211 +234 @@
     dxMax = MAX (dxMax, skyX_rt[iy] - skyX_lf[iy]); // this is a constant value unless we ignore some of the pixels
     // assert skyX_rt > skyX_lf ?
+
+    if (SaveSegments) {
+      skyLsV->elements.Flt[iy] = S_s_l + Xtrace;
+      skyLeV->elements.Flt[iy] = S_e_l + Xtrace;
+      skyLfV->elements.Flt[iy] = skyF_lf[iy];
+      skyRsV->elements.Flt[iy] = S_s_r + Xtrace;
+      skyReV->elements.Flt[iy] = S_e_r + Xtrace;
+      skyRfV->elements.Flt[iy] = skyF_rt[iy];
+      bckLsV->elements.Flt[iy] = B_s_l + Xtrace;
+      bckLeV->elements.Flt[iy] = B_e_l + Xtrace;
+      bckLfV->elements.Flt[iy] = bckF_lf[iy];
+      bckRsV->elements.Flt[iy] = B_s_r + Xtrace;
+      bckReV->elements.Flt[iy] = B_e_r + Xtrace;
+      bckRfV->elements.Flt[iy] = bckF_rt[iy];
+    }
   }
 
@@ -253 +291 @@
     float Xtrace = Xref - Xpsf + dXtrace;
     
+    // include the PSF trace and restrict the fit below to portions where the profile is > 0.95
+
     // now I need to find the limits of the PSF in the buffer at this iy location
     for (int ipsf = 0; ipsf < Npsf; ipsf ++) {
@@ -267 +307 @@
     float S1 = 0.0, S2 = 0.0;
     for (int ipsf = 0; ipsf < Npsf; ipsf++) {
+      if (iy == 6000) {
+        fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], psfV[ipsf], tmpv[ipsf]*psfV[ipsf]);
+      }
       S1 += tmpv[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
       S2 += psfV[ipsf]*psfV[ipsf]; // if we have an errorbar, include / SQ(sigv[i])

trunk/Ohana/src/opihi/mana/deimos_mkmodel.c

@@ -164 +164 @@
         opihi_flt g = 0.0, s = 0.0;
         for (int n = -Nkernel_y2; n <= Nkernel_y2; n++) {
-          if (iy + n < 0) continue; // bottom edge of full image
-          if (iy + n >= Ny) continue; // top edge of full image
+          if (iy - n < 0) continue; // bottom edge of full image
+          if (iy - n >= Ny) continue; // top edge of full image
           int pix = ix + (n + Nkernel_y2)*Nx;
           myAssert (pix < Nkernel_pix, "oops (apply_tilt)");
-          s += kernel[pix]*input[ix + (iy + n)*Nx];
+          s += kernel[pix]*input[ix + (iy - n)*Nx];
           g += kernel[pix]; // optimization: save g[ix], all should be the same
         }

trunk/Ohana/src/opihi/mana/deimos_mkobj.c

@@ -208 +208 @@
 
   float sin_stilt = sin(stilt*RAD_DEG);
-  int Nyk = ceil(Nx * fabs(sin_stilt));
+  int Nyk = ceil(Nx * fabs(sin_stilt)) + 1; // one extra pixel buffer
   if (Nyk < 3) Nyk = 3; // minimum of 3 pixels (or kernel assumption below fails)
   Nyk += NlsfFull;
@@ -260 +260 @@
       myAssert (Nkpix_i < Nkpix, "oops");
       kernel[Nkpix_i] = dyf;
-      // fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]);
+      fprintf (stderr, "%d, %d: %f %f\n", ix, (dyi + Nyk2), kernel[ix + (dyi + Nyk2 + 0)*Nx], kernel[ix + (dyi + Nyk2 + 1)*Nx]);
     }
   }
@@ -280 +280 @@
         opihi_flt g = 0.0, s = 0.0;
         for (int n = -Nyk2; n <= Nyk2; n++) {
-          if (iy + n < 0) continue; // bottom edge of full image
-          if (iy + n >= Ny) continue; // top edge of full image
+          if (iy - n < 0) continue; // bottom edge of full image
+          if (iy - n >= Ny) continue; // top edge of full image
           int Nkpix_i = ix + (n + Nyk2)*Nx;
           myAssert (Nkpix_i < Nkpix, "oops");
-          s += kernel[Nkpix_i]*s1[ix + (iy + n)*Nx];
+          s += kernel[Nkpix_i]*s1[ix + (iy - n)*Nx];
           g += kernel[Nkpix_i];
         }

trunk/Ohana/src/opihi/mana/deimos_mkslit.c

@@ -56 +56 @@
 
   // define the output window
-  int Nx = 2*profile[0].Nelements;
+  int Nx = 2*profile[0].Nelements + 1;
   int Ny = flux ? flux[0].Nelements : fluxbins;