Changeset 41159 for trunk/Ohana/src/opihi/mana/deimos_mkobj.c

Timestamp:

Nov 27, 2019, 11:04:51 AM (7 years ago)

Author:

eugene

Message:

improvments to deimos extraction tools

File:

• trunk/Ohana/src/opihi/mana/deimos_mkobj.c (modified) (7 diffs)

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

@@ -6 +6 @@
 
   // input parameters:
-  // * trace       : spline fit of slit central x pos vs y-coord
-  // * profile     : slit window profile (vector)
-  // * object      : vector of object flux vs y-coord 
-  // * sky         : vector of local sky signal vs y-coord
-  // * background  : vector of extra-slit background flux vs y-coord
-  // * PSF         : point-spread function vector (flux normalized, x-dir)
-  // * LSF         : sline-spread function vector (flux normalized, y-dir)
-  // * stilt       : slit tilt response : 2D kernel? 
+  // * slit_trace_red : spline fit of slit central x pos vs y-coord for red chip
+  // * slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip
+  // * psf_trace      : spline fit of psf central x pos vs y-coord
+  // * profile        : slit window profile (vector)
+  // * object         : vector of object flux vs y-coord 
+  // * sky            : vector of local sky signal vs y-coord
+  // * backgnd        : vector of extra-slit background flux vs y-coord
+  // * PSF            : point-spread function vector (flux normalized, x-dir)
+  // * LSF            : sline-spread function vector (flux normalized, y-dir)
+  // * stilt          : slit tilt response : 2D kernel? 
+
+  // * redlimit       : pixel break between red and blu chips
+  // * Nwave          : number of wavelength pixels in output image
 
   int N;
 
   // if any of these are not defined, they will have assumed identity values
-  Spline *trace   = NULL;
-  Vector *profile = NULL;
-  Vector *object  = NULL;
-  Vector *sky     = NULL;
-  Vector *backgnd = NULL;
-
-  Vector *psf     = NULL;
-  Vector *lsf     = NULL;
+  Spline *slit_trace_red = NULL;
+  Spline *slit_trace_blu = NULL;
+  Spline *psf_trace  = NULL;
+  Vector *profile    = NULL;
+  Vector *object     = NULL;
+  Vector *sky        = NULL;
+  Vector *backgnd    = NULL;
+                     
+  Vector *psf        = NULL;
+  Vector *lsf        = NULL;
   // add tilt later
 
-  Buffer *output  = NULL;
+  Buffer *output     = NULL;
 
   float stilt = 0.0; // angle of the slit
 
-  if ((N = get_argument (argc, argv, "-trace"))) {
-    remove_argument (N, &argc, argv);
-    if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
+  // for a red vs blu spline, we need to specify the split point
+  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
+  // this more generic (need to define the ranges somewhere)
+  int redlimit = 4096;
+  if ((N = get_argument (argc, argv, "-redlimit"))) {
+    remove_argument (N, &argc, argv);
+    redlimit = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
+    remove_argument (N, &argc, argv);
+    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+  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);
   }
@@ -81 +108 @@
   }
 
+  int VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = TRUE;
+  }
+
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
+    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-object vector] [-sky vector] [-backgnd vector] [-slit-trace-red spline] [-slit-trace-blu spline] [-psf-trace spline] [-profile vector] [-psf vector] [-lsf vector] [-stilt angle] [-Nwave Npixel] [-redlimit pixel]\n");
+    gprint (GP_ERR, "  slit_trace_red : spline fit of slit central x pos vs y-coord for red chip\n");
+    gprint (GP_ERR, "  slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip\n");
+    gprint (GP_ERR, "  psf_trace      : spline fit of psf central x pos vs y-coord\n");
+    gprint (GP_ERR, "  profile        : slit window profile (vector)\n");
+    gprint (GP_ERR, "  object         : vector of object flux vs y-coord \n");
+    gprint (GP_ERR, "  sky            : vector of local sky signal vs y-coord\n");
+    gprint (GP_ERR, "  backgnd        : vector of extra-slit background flux vs y-coord\n");
+    gprint (GP_ERR, "  PSF            : point-spread function vector (flux normalized, x-dir)\n");
+    gprint (GP_ERR, "  LSF            : sline-spread function vector (flux normalized, y-dir)\n");
+    gprint (GP_ERR, "  stilt          : slit tilt response : 2D kernel? \n");
+    gprint (GP_ERR, "  redlimit       : pixel break between red and blu chips\n");
+    gprint (GP_ERR, "  Nwave          : number of wavelength pixels in output image\n");
     return FALSE;
   }
@@ -123 +168 @@
   // if we are appying a trace offset spline, we need to generate an output window which
   // is Nx + the full swing of the trace, then window back down
-  if (trace) {
+  if (slit_trace_red && slit_trace_blu) {
     float dXmin = +1000;
     float dXmax = -1000;
     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 dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
+      float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       dXmin = MIN (dx, dXmin);
       dXmax = MAX (dx, dXmax);
@@ -181 +227 @@
     if (psf) Npof -= Npsf;
 
+    // integral and fractional pixel offsets of PSF (none, by default)
+    int dxi = 0;
+    float dxf = 0.0, dxr = 1.0;
+    
+    if (psf_trace) {
+      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
+      dxi = floor(dx);
+      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+      dxr = 1 - dxf;
+      Npof += dxi;
+      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
+      if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
+    }
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
+
     // flux = obj * PSF + sky
     for (int ix = 0; ix < Nx; ix++) {
@@ -191 +254 @@
         // only add in the flux if we are in range of the PSF
         if ((n >= 0) && (n < NpsfFull)) {
-          value += objVy * psfV[n];
+          if (doInterp) {
+            if ((n > 0) && (n < NpsfFull)) {
+              if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
+                value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
+              }
+            }
+            if (n == 0) {
+              if (isfinite(psfV[n])) {
+                value += objVy*psfV[n]*dxr;
+              }
+            }
+          } else {
+            if (isfinite(psfV[n])) {
+              value += objVy*psfV[n];
+            }
+          }
         }
       } else {
@@ -329 +407 @@
 
   // shift pixels in x based on the trace
-  if (trace) {
+  if (slit_trace_red && slit_trace_blu) {
     ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
     float *outTrace = (float *) outTraceBuffer;
@@ -337 +415 @@
 
       // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
-      float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
+      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
+      float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
       
       // extract the integer pixel offset and the fractional offset