Index: branches/eam_branches/ohana.20190329/src/opihi/cmd.data/interpolate_presort.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/cmd.data/interpolate_presort.c (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/cmd.data/interpolate_presort.c (revision 40663) @@ -7,4 +7,11 @@ double x0, dx, dy, y0; Vector *xout, *yout, *xin, *yin; + + int N; + int FillEnds = FALSE; + if ((N = get_argument (argc, argv, "-fill-ends"))) { + FillEnds = TRUE; + remove_argument (N, &argc, argv); + } /** check basic syntax **/ @@ -38,8 +45,15 @@ /* in and out vectors are sorted */ j = 0; + i = 0; + + // fill in the rest with NANs + while (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]) { + yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN; + i++; + } // every output pixel should get a value unless // we are below the lowest in or above the highest in - for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { + for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) { yout[0].elements.Flt[i] = NAN; @@ -66,10 +80,10 @@ // fill in the rest with NANs + int NinLast = xin[0].Nelements - 1; while (i < yout[0].Nelements) { - yout[0].elements.Flt[i] = NAN; + yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN; i++; } return (TRUE); - } Index: branches/eam_branches/ohana.20190329/src/opihi/include/deimos.h =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/include/deimos.h (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/include/deimos.h (revision 40663) @@ -7,5 +7,5 @@ } DeimosResult; -float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row); +float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace, Spline *psf_trace, int Nx, int Ny, int row); void deimos_make_kernel (float stilt, int Nx); void deimos_free_kernel (); @@ -13,8 +13,8 @@ // internal to make_model: -float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row); +float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row); float *deimos_apply_tilt (float *input, int Nx, int Ny); void deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny); void deimos_add_background (opihi_flt *backgnd, float *out, int Nx, int Ny, int row); -float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row); +float *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row); Index: branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/Makefile (revision 40663) @@ -25,4 +25,5 @@ $(SRC)/deimos_fitslit.$(ARCH).o \ $(SRC)/deimos_getobj.$(ARCH).o \ +$(SRC)/deimos_getalt.$(ARCH).o \ $(SRC)/deimos_fitobj.$(ARCH).o \ $(SRC)/deimos_fitalt.$(ARCH).o \ Index: branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/deimos.c (revision 40663) @@ -6,4 +6,5 @@ int deimos_mkalt (int argc, char **argv); int deimos_getobj (int argc, char **argv); +int deimos_getalt (int argc, char **argv); int deimos_fitobj (int argc, char **argv); int deimos_fitalt (int argc, char **argv); @@ -14,4 +15,5 @@ {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"}, {1, "getobj", deimos_getobj, "determine crude object parameters"}, + {1, "getalt", deimos_getalt, "determine crude object parameters"}, {1, "mkobj", deimos_mkobj, "make a full object image"}, {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"}, Index: branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_fitalt.c (revision 40663) @@ -27,5 +27,6 @@ // 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 slit central x pos vs y-coord // * profile : slit window profile (vector) // * PSF : point-spread function vector (flux normalized, x-dir) @@ -39,5 +40,6 @@ int N; - Spline *trace = NULL; + Spline *slit_trace = NULL; + Spline *psf_trace = NULL; Vector *profile = NULL; @@ -106,7 +108,12 @@ // 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; } @@ -244,5 +251,5 @@ // current vector value: // for deimos_make_model, row is the starting point of the subimage - float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0); + float *model_ref = deimos_make_model (objCurr, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0); // in the test functions below, dRow is the central pixel of the subset, e.g., obj[row] => objCurr[dRow] @@ -252,13 +259,13 @@ // delta obj (save dObj) opihi_flt *objtest = deimos_make_test (objCurr, objNoise[row], dRow, Nrow, &dObj); - float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, trace, Nx, Nrow, row_0); + float *model_obj = deimos_make_model (objtest, skyCurr, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0); // delta sky (save dSky) opihi_flt *skytest = deimos_make_test (skyCurr, skyNoise[row], dRow, Nrow, &dSky); - float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, trace, Nx, Nrow, row_0); + float *model_sky = deimos_make_model (objCurr, skytest, bckCurr, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0); // delta bck (save dBck) opihi_flt *bcktest = deimos_make_test (bckCurr, bckNoise[row], dRow, Nrow, &dBck); - float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, trace, Nx, Nrow, row_0); + float *model_bck = deimos_make_model (objCurr, skyCurr, bcktest, psf, profile, slit_trace, psf_trace, Nx, Nrow, row_0); // dObj, dSky, dBck carry the delta used derivatives in and the delta applied to the values out @@ -294,5 +301,5 @@ } - float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, trace, Nx, Ny, 0); + float *model_full = deimos_make_model (objVal, skyVal, bckVal, psf, profile, slit_trace, psf_trace, Nx, Ny, 0); int Npts = 0; @@ -305,5 +312,5 @@ usage: - gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + gprint (GP_ERR, "USAGE: deimos fitobj (buffer) (Xref) -object vector -sky vector -backgnd vector -slit-trace spline -psf-trace spline -profile vector -psf vector -stilt (angle)\n"); return FALSE; } Index: branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getalt.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getalt.c (revision 40663) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getalt.c (revision 40663) @@ -0,0 +1,287 @@ +# include "data.h" + +// temporary storage used in local functions: + +static float *tmpv = NULL; +static float *tmpx = NULL; + +float get_median_flt (float *val, int Nval) { + + if (Nval <= 0) return NAN; + + int Nval2 = Nval / 2; + + fsort (val, Nval); + float value = (Nval % 2) ? val[Nval2] : 0.5*(val[Nval2-1] + val[Nval2]); + + return value; +} + +int deimos_getalt (int argc, char **argv) { + + // determine rough object model parameters based on quick stats + // input : buffer PSF trace profile stilt + // output : object sky background + // not used : LSF -- not needed (measurement per disp pixel) + + // input parameters: + // * buffer : 2D image of slit (full frame or cutout?) + // * 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) + // * stilt : slit tilt response : 2D kernel? + + // output parameters: + // * 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 + + int N; + + Spline *slit_trace = NULL; + Spline *psf_trace = NULL; + Vector *profile = NULL; + + Vector *object = NULL; + Vector *sky = NULL; + Vector *backgnd = NULL; + + Vector *psf = NULL; + + Buffer *buffer = NULL; + + float stilt = 0.0; // angle of the slit + + // Input parameters: + 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; } + if ((N = get_argument (argc, argv, "-profile"))) { + remove_argument (N, &argc, argv); + if ((profile = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-psf"))) { + remove_argument (N, &argc, argv); + if ((psf = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-stilt"))) { + remove_argument (N, &argc, argv); + stilt = atof (argv[N]); + remove_argument (N, &argc, argv); + } else { goto usage; } + + // Output parameters: + if ((N = get_argument (argc, argv, "-object"))) { + remove_argument (N, &argc, argv); + if ((object = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-sky"))) { + remove_argument (N, &argc, argv); + if ((sky = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + if ((N = get_argument (argc, argv, "-backgnd"))) { + remove_argument (N, &argc, argv); + if ((backgnd = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE); + remove_argument (N, &argc, argv); + } else { goto usage; } + + if (argc != 3) goto usage; + + // buffer is the full image, Xref is reference coordinate of the profile in the buffer + if ((buffer = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE); + float Xref = atof(argv[2]); + + // profile and PSF are defined with central reference pixel mapped to Xref (+ trace) + int Nprofile = profile->Nelements; + if (Nprofile % 2 == 0) { + gprint (GP_ERR, "slit profile vector must have an odd number of pixels\n"); + return FALSE; + } + int Xprofile = Nprofile / 2; + + int Npsf = psf->Nelements; + if (Npsf % 2 == 0) { + gprint (GP_ERR, "PSF vector must have an odd number of pixels\n"); + return FALSE; + } + int Xpsf = Npsf / 2; + + // we will generate output vectors (object, sky, background) of length Ny + int Nx = buffer[0].matrix.Naxis[0]; + int Ny = buffer[0].matrix.Naxis[1]; + ResetVector (object, OPIHI_FLT, Ny); + ResetVector (sky, OPIHI_FLT, Ny); + ResetVector (backgnd, OPIHI_FLT, Ny); + + // generate two buffers to store the flux and coordinate for the background, sky, object pixels: + ALLOCATE (tmpv, float, Nx); + ALLOCATE (tmpx, float, Nx); + + float *buffV = (float *) buffer[0].matrix.buffer; + + float dxMax = 0; // maximum distance between sky regions + + // extract sky-subtracted values for the psf fit, using the two sky vectors to predict the sky + float sin_stilt = sin(stilt*RAD_DEG); + float dyMax = dxMax*sin_stilt; // dxMax is based on the profile width, calculated above + + // we are fitting this PSF to the data: + opihi_flt *psfV = psf->elements.Flt; + + // init the regions which are outside the valid region + for (int iy = 0; iy < fabs(dyMax); iy++) { + object->elements.Flt[iy] = 0; + sky->elements.Flt[iy] = 0; + backgnd->elements.Flt[iy] = 0; + } + for (int iy = Ny - fabs(dyMax); iy < Ny; iy++) { + object->elements.Flt[iy] = 0; + sky->elements.Flt[iy] = 0; + backgnd->elements.Flt[iy] = 0; + } + + for (int iy = fabs(dyMax); iy < Ny - fabs(dyMax); iy++) { + + // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center + 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?) + + // XtracePSF is the offset between the buffer ix and the psf elements at this location + float XtracePSF = Xref - Xpsf + dXtraceSLIT + dXtracePSF; + + // XtraceSLIT is the offset between the buffer ix and the slit profile elements at this location + float XtraceSLIT = Xref - Xprofile + dXtraceSLIT; + + // measure the background outside the slit + + // select pixels from the image which are in the background regions (profile < 0.15) + int Nbck = 0; + int Nsky = 0; + for (int ix = 0; ix < Nx; ix ++) { + int islit = ix - XtraceSLIT - 0.5; // slit window coordinate + int ipsf = ix - XtracePSF - 0.5; // psf model coordinate + + if (islit < 0) continue; + if (islit >= Nprofile) continue; + + // coordinate of the PSF center: + int dxpsf = ipsf - Xpsf; + int dypsf = dxpsf * sin_stilt; + + int iY = iy + dypsf; + + if (iy == 6000) { + fprintf (stderr, "bck : %d %d %d : %f : %f\n", ix, islit, ipsf, buffV[ix + iY*Nx], profile->elements.Flt[islit]); + } + + // pixels in the off-slit background: + if (profile->elements.Flt[islit] < 0.15) { + tmpv[Nbck] = buffV[ix + iY*Nx]; + if (!isfinite(tmpv[Nbck])) continue; + Nbck ++; + continue; + } + + // the psf window can be smaller than the slit window + if (ipsf < 0) continue; + if (ipsf >= Npsf) continue; + + // pixels in the on-slit sky: + if ((profile->elements.Flt[islit] > 0.95) && (psfV[ipsf] < 0.15)) { + tmpx[Nsky] = buffV[ix + iY*Nx]; + if (!isfinite(tmpx[Nsky])) continue; + Nsky ++; + } + } + backgnd->elements.Flt[iy] = get_median_flt (tmpv, Nbck); + float So = get_median_flt (tmpx, Nsky); + sky->elements.Flt[iy] = So; + + // now I need to find the limits of the PSF in the buffer at this iy location + int Ntmp = 0; + for (int ix = 0; ix < Nx; ix ++) { + int islit = ix - XtraceSLIT - 0.5; // slit window coordinate + int ipsf = ix - XtracePSF - 0.5; // psf model coordinate + + if (islit < 0) continue; + if (islit >= Nprofile) continue; + + // the psf window can be smaller than the slit window + if (ipsf < 0) continue; + if (ipsf >= Npsf) continue; + + // coordinate of the PSF center: + int dxpsf = ipsf - Xpsf; + int dypsf = dxpsf * sin_stilt; + + int iY = iy + dypsf; + + // filter on the PSF value? + if (profile->elements.Flt[islit] < 0.95) continue; + + if (iy == 6000) { + fprintf (stderr, "%d %d %d : %f %f : %f : %f\n", ix, islit, ipsf, buffV[ix + iY*Nx] - So, psfV[ipsf], (buffV[ix + iY*Nx] - So)*psfV[ipsf], profile->elements.Flt[islit]); + } + + tmpv[Ntmp] = buffV[ix + iY*Nx] - So; + tmpx[Ntmp] = psfV[ipsf]; + Ntmp ++; + } + + // I now have a pair of vectors (tmpx, tmpv) which are aligned to the PSF at this location + // now fit the PSF to tmpv by calculating a normalization + float S1 = 0.0, S2 = 0.0; + for (int ipsf = 0; ipsf < Ntmp; ipsf++) { + if (iy == 6000) { + fprintf (stderr, "%d : %f %f : %f\n", ipsf, tmpv[ipsf], tmpx[ipsf], tmpv[ipsf]*tmpx[ipsf]); + } + 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]) + } + object->elements.Flt[iy] = S1 / S2; + } + + free (tmpv); + free (tmpx); + + return TRUE; + + usage: + gprint (GP_ERR, "USAGE: deimos getobj (buffer) (Xref) -object vector -sky vector -backgnd vector -trace spline -profile vector -psf vector -stilt (angle)\n"); + // can I make default values for trace (0.0), profile (?), psf (fraction of profile?), stilt (0.0) + return FALSE; +} + + /* + // convolve the object flux at this y-coord with PSF and add sky + for (int ix = 0; i < Nx; i++) { + opihi_flt g = s = 0; + for (int n = -Npsf; n <= Npsf; n++) { + if (ix + n < 0) continue; + if (ix + n >= Nx) continue; + opihi_flt value = objV ? objV[ix+n] : 1.0; + s += psfV[n]*value; + g += psfV[n]; + } + out[ix + iy*Nx] = Npsf ? s / g : (objV ? objV[ix+n] : 1.0); + if (skyV) { + out[ix + iy*Nx] += skyV[iy] + } + */ + Index: branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_getobj.c (revision 40663) @@ -3,6 +3,6 @@ // 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,5 +39,6 @@ // 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,5 +52,6 @@ int N; - Spline *trace = NULL; + Spline *slit_trace = NULL; + Spline *psf_trace = NULL; Vector *profile = NULL; @@ -71,7 +73,12 @@ // 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,8 +225,9 @@ // 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,21 +292,36 @@ // 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,13 +329,13 @@ // 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; Index: branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c =================================================================== --- branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c (revision 40662) +++ branches/eam_branches/ohana.20190329/src/opihi/mana/deimos_mkmodel.c (revision 40663) @@ -86,11 +86,11 @@ // generate a subimage of the full model for the range row - row + Ny // obj,sky,bck are subset vectors corresponding to the range row to row+Ny -float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) { - - float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row); +float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace, Spline *psf_trace, int Nx, int Ny, int row) { + + float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row); float *new = deimos_apply_tilt (raw, Nx, Ny); deimos_apply_profile (profile, new, Nx, Ny); deimos_add_background (bck, new, Nx, Ny, row); - float *out = deimos_apply_trace (trace, new, Nx, Ny, row); + float *out = deimos_apply_trace (slit_trace, new, Nx, Ny, row); free (raw); @@ -104,5 +104,5 @@ // *** generate a raw image with no slit tilt and no trace offsets. apply the PSF to the // object flux, add in the sky flux -float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) { +float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) { OHANA_UNUSED_PARAM(row); @@ -126,4 +126,19 @@ opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row] + // 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; + Xref += dxi; + } + + // 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++) { @@ -136,5 +151,14 @@ // 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)) { + value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf; + } + if (n == 0) { + value += objVy*psfV[n]*dxr; + } + } else { + value += objVy * psfV[n]; + } } out[ix + iy*Nx] = value; @@ -219,5 +243,5 @@ // shift pixels in x based on the trace -float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) { +float *deimos_apply_trace (Spline *slit_trace, float *input, int Nx, int Ny, int row) { ALLOCATE_PTR (output, float, Nx*Ny); @@ -227,5 +251,5 @@ // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center // NOTE: spline is evaluated at full dispersion coordinate - float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row); + float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy + row); // extract the integer pixel offset and the fractional offset