Changeset 40663 for branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c

Timestamp:

Mar 31, 2019, 3:08:13 PM (7 years ago)

Author:

eugene

Message:

deimos updates to add psf trace

File:

• branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c (modified) (7 diffs)

branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c

@@ -3 +3 @@
 // temporary storage used in local functions:
 
-float *tmpv = NULL;
-float *tmpx = NULL;
+static float *tmpv = NULL;
+static float *tmpx = NULL;
 
 // pass in a pointer to the row (&buffV[iy*Nx])
@@ -39 +39 @@
   // input parameters:
   // * buffer      : 2D image of slit (full frame or cutout?)
-  // * trace       : spline fit of slit central x pos vs y-coord
+  // * slit_trace  : spline fit of slit central x pos vs y-coord
+  // * psf_trace   : spline fit of psf central x pos vs y-coord
   // * profile     : slit window profile (vector)
   // * PSF         : point-spread function vector (flux normalized, x-dir)
@@ -51 +52 @@
   int N;
 
-  Spline *trace   = NULL;
+  Spline *slit_trace = NULL;
+  Spline *psf_trace  = NULL;
   Vector *profile = NULL;
 
@@ -71 +73 @@
 
   // Input parameters:
-  if ((N = get_argument (argc, argv, "-trace"))) {
-    remove_argument (N, &argc, argv);
-    if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
+  if ((N = get_argument (argc, argv, "-slit-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  } else { goto usage; }
+  if ((N = get_argument (argc, argv, "-psf-trace"))) {
+    remove_argument (N, &argc, argv);
+    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
     remove_argument (N, &argc, argv);
   } else { goto usage; }
@@ -218 +225 @@
 
   // loop over the y positions
+  // examine the slit window profile to measure background and slit sky
   for (int iy = 0; iy < Ny; iy++) {
 
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
-    float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+    float dXtrace = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
 
     // Xref is the reference coordinate of the profile on the image
@@ -284 +292 @@
     
     // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
-    float dXtrace = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+    float dXtraceSLIT = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
+    float dXtracePSF  = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
       
     // Xref is the reference coordinate of the profile on the image
     // Xpsf is the center of the PSF.  
     // dXtrace follows the x-dir variations of the slit profile (or PSF?)
-    float Xtrace = Xref - Xpsf + dXtrace;
-    
+    float XtracePSF = Xref - Xpsf + dXtraceSLIT + dXtracePSF;
+    float XtraceSLIT = Xref - Xprofile + dXtraceSLIT;
+
     // 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
+    int Ntmp = 0;
     for (int ipsf = 0; ipsf < Npsf; ipsf ++) {
 
-      int ix = ipsf + Xtrace + 0.5; // image coordinate (nearest pixel center)
+      int ix = ipsf + XtracePSF + 0.5; // image coordinate (nearest pixel center)
+      int islit = ix - XtraceSLIT - 0.5; // image coordinate (nearest pixel center)
+
       float So = S_lf + (ix - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
 
-      tmpv[ipsf] = buffV[ix + iy*Nx] - So;
-      tmpx[ipsf] = ipsf;
+      if (islit < 0) continue;
+      if (islit >= profile->Nelements) continue;
+
+      if (iy == 6000) {
+        fprintf (stderr, "%d : %f %f : %f : %f\n", ipsf, buffV[ix + iy*Nx] - So, psfV[ipsf], (buffV[ix + iy*Nx] - So)*psfV[ipsf], profile->elements.Flt[islit]);
+      }
+
+      if (profile->elements.Flt[islit] < 0.95) continue;
+
+      tmpv[Ntmp] = buffV[ix + iy*Nx] - So;
+      tmpx[Ntmp] = psfV[ipsf];
+      Ntmp ++;
     }
 
@@ -306 +329 @@
     // now fit the PSF to tmpv by calculating a normalization
     float S1 = 0.0, S2 = 0.0;
-    for (int ipsf = 0; ipsf < Npsf; ipsf++) {
+    for (int ipsf = 0; ipsf < Ntmp; ipsf++) {
       if (iy == 6000) {
-        fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], psfV[ipsf], tmpv[ipsf]*psfV[ipsf]);
+        fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], tmpx[ipsf], tmpv[ipsf]*tmpx[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])
+      S1 += tmpv[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
+      S2 += tmpx[ipsf]*tmpx[ipsf]; // if we have an errorbar, include / SQ(sigv[i])
     }
 
-    int Xo = Xref - dXtrace; // center of PSF in image x-coordinates
+    int Xo = Xref - dXtraceSLIT - dXtracePSF; // center of PSF in image x-coordinates
     float So = S_lf + (Xo - X_lf)*(S_rt - S_lf)/(X_rt - X_lf); // interpolated sky under PSF peak
     object->elements.Flt[iy] = S1 / S2;