Changeset 40622 for trunk

Timestamp:

Feb 8, 2019, 8:20:03 PM (7 years ago)

Author:

eugene

Message:

add deimos fitslit, mkobj, mkslit commands

Location:

trunk/Ohana/src/opihi/mana

Files:

• Makefile (modified) (1 diff)
• deimos.c (modified) (1 diff)
• deimos_commands.c (modified) (2 diffs)
• deimos_fitslit.c (added)
• deimos_mkobj.c (added)
• deimos_mkslit.c (added)
• init.c (modified) (2 diffs)

trunk/Ohana/src/opihi/mana/Makefile

@@ -21 +21 @@
 $(SRC)/init.$(ARCH).o \
 $(SRC)/mana.$(ARCH).o \
+$(SRC)/deimos_mkobj.$(ARCH).o \
+$(SRC)/deimos_mkslit.$(ARCH).o \
+$(SRC)/deimos_fitslit.$(ARCH).o \
 $(SRC)/findrowpeaks.$(ARCH).o 
 
 cmds = \
+$(SRC)/deimos.$(ARCH).o \
 $(SRC)/rawstars.$(ARCH).o \
 $(SRC)/fitcontour.$(ARCH).o \

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

@@ -1 +1 @@
 # include "data.h"
 
-int deimos_list (int argc, char **argv);
+int deimos_mkslit (int argc, char **argv);
+int deimos_fitslit (int argc, char **argv);
+int deimos_mkobj (int argc, char **argv);
 
-static Command book_commands[] = {
-  {1, "list",     book_list,     "list books"},
+static Command deimos_commands[] = {
+  {1, "mkobj", deimos_mkobj, "make a full object image"},
+  {1, "mkslit", deimos_mkslit, "make a slit image"},
+  {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
 };
 
-int book_command (int argc, char **argv) {
+int deimos (int argc, char **argv) {
 
   int i, N, status;
 
   if (argc < 2) {
-    gprint (GP_ERR, "USAGE: book (command)\n");
-    gprint (GP_ERR, "    book list                                  : list books\n");
+    gprint (GP_ERR, "USAGE: deimos (command)\n");
+    gprint (GP_ERR, "    deimos mkslit : make slit image\n");
     return (FALSE);
   }
 
-  N = sizeof (book_commands) / sizeof (Command);
+  N = sizeof (deimos_commands) / sizeof (Command);
 
-  /* find the book sub-command which matches */
+  /* find the deimos sub-command which matches */
   for (i = 0; i < N; i++) {
-    if (!strcmp (book_commands[i].name, argv[1])) {
-      status = (*book_commands[i].func) (argc - 1, argv + 1);
+    if (!strcmp (deimos_commands[i].name, argv[1])) {
+      status = (*deimos_commands[i].func) (argc - 1, argv + 1);
       return (status);
     }
   }
 
-  gprint (GP_ERR, "unknown book command %s\n", argv[1]);
+  gprint (GP_ERR, "unknown deimos command %s\n", argv[1]);
   return (FALSE);
 }

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

@@ -2 +2 @@
 
 int deimos_mkslit (int argc, char **argv) {
-  OHANA_UNUSED_PARAM(argv);
-  if (argc != 5) {
-    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (flux) (buffer)\n");
-    return FALSE;
-  }
 
   // generate model observed flux given a uniformly illuminated slit
@@ -16 +11 @@
 
   Vector *profile = NULL;
-  Vector *flux = NULL;
-  Spline *trace = NULL;
-  Buffer *buff = NULL;
+  Vector *flux    = NULL;
+  Vector *sky     = NULL;
+  Spline *trace   = NULL;
+  Buffer *buff    = NULL;
+
+  // user-specified flux value (otherwise assumed to be 1)
+  if ((N = get_argument (argc, argv, "-flux"))) {
+    remove_argument (N, &argc, argv);
+    if ((flux = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // user-specified sky value (otherwise assumed to be 0)
+  if ((N = get_argument (argc, argv, "-sky"))) {
+    remove_argument (N, &argc, argv);
+    if ((sky = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+    remove_argument (N, &argc, argv);
+  }
+
+  // either flux or fluxbins must be specified to define output size
+  int fluxbins = 0;
+  if ((N = get_argument (argc, argv, "-fluxbins"))) {
+    remove_argument (N, &argc, argv);
+    Nrows = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
+    return FALSE;
+  }
+
+  if (!flux && !fluxbins) {
+    gprint (GP_ERR, "either flux or fluxbins must be specified to define output size\n");
+    return FALSE;
+  }
 
   // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
   if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
   if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
-  if ((flux    = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
-  if ((buff    = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
-  
+  if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
+
   // define the output window
   int Nx = 2*profile[0].Nelements;
-  int Ny = flux[0].Nelements;
+  int Ny = flux ? flux[0].Nelements : fluxbins;
   
   ResetBuffer (buff, Nx, Ny, -32, 0.0, 1.0);
-  
+
+  // the profile is registered to the midpoint of the vector : (int) N / 2
+  // the output window is registered to the midpoint of the window : (int) Nx / 2
+
+  int NxMidProfile = profile->Nelements / 2;
+  int NxMidWindow  = Nx / 2;
+
+  float *out = (float *) buff[0].matrix.buffer;
 
   // loop over the y positions
-  for (iy = 0; iy < flux->Nelements; iy++) {
+  for (int iy = 0; iy < Ny; iy++) {
 
-    // evaluate the trace spline at this y-coord to find the x-coord of the profile center
-    opihi_flt dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, 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);
     
-    // interpolate the profile with offset of dx
-    for (ix = 0; ix < profile->Nelements; ix++) {
-      // interpolate from profile->elements.Flt[ix], interpolate from profile->elements.Flt[ix+1] to ix+dx
+    // extract the integer pixel offset and the fractional offset
+    int dxi = floor(dx); // -1.7 -> -2, -0.5 -> -1, +0.5 -> 0, +1.7 -> 1
+    float dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
+
+    // starting point in output window is :
+    // XXX handle case if profile is wider than window
+    int Sx = (int)(NxMidWindow - NxMidProfile) + dxi;
+
+    // if fractional offset is small, do not interpolate
+    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
+
+    // loop over the relevant pixels in the window:
+    for (int ix = 0; ix < Nx; ix++) {
+
+      // equivalent coord in the profile:
+      int px = ix - Sx;
+      if (px < 0) continue;
+      if (px >= profile->Nelements) continue;
+
+      float vout = NAN;
+
+      if (doInterp) {
+        if ((px > 0) && (px < profile->Nelements - 1)) {
+          // if ((px == 5) && (iy % 50 == 0)) {
+          //   fprintf (stderr, "%d : %d %d : %d, %f\n", iy, Sx, px, dxi, dxf);
+          // }
+          // XXX something is not quite right in my math: I think (1-dx), dx are reversed
+          // but if I use what I expect the interpolation is backwards...
+          vout = profile->elements.Flt[px]*(1 - dxf) + profile->elements.Flt[px-1]*dxf;
+        }
+        if (px == 0) {
+          vout = profile->elements.Flt[px]*dxf;
+        }
+        if (px == profile->Nelements - 1) {
+          vout = profile->elements.Flt[profile->Nelements - 1]*(1.0 - dxf);
+        }
+      } else {
+        vout = profile->elements.Flt[px];
+      }
+
+      if (flux) vout *= flux->elements.Flt[iy];
+      if (sky) vout += sky->elements.Flt[iy];
+
+      out[ix + iy*Nx] = vout;
     }
-
   }
 

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

@@ -1 +1 @@
 # include "opihi.h"
 
+int deimos          PROTO((int, char **));
 int findpeaks       PROTO((int, char **));
 int fitcontour      PROTO((int, char **));
@@ -8 +9 @@
 
 static Command cmds[] = {  
+  {1, "deimos",      deimos,       "deimos multislit spectrograph tools"},
   {1, "findpeaks",   findpeaks,    "find image peaks"},
   {1, "fitcontour",  fitcontour,   "fit ellipse contour"},