aboutsummaryrefslogtreecommitdiff
path: root/noao/obsutil/src/specfocus/t_specfocus.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/obsutil/src/specfocus/t_specfocus.x')
-rw-r--r--noao/obsutil/src/specfocus/t_specfocus.x762
1 files changed, 762 insertions, 0 deletions
diff --git a/noao/obsutil/src/specfocus/t_specfocus.x b/noao/obsutil/src/specfocus/t_specfocus.x
new file mode 100644
index 00000000..7f176de5
--- /dev/null
+++ b/noao/obsutil/src/specfocus/t_specfocus.x
@@ -0,0 +1,762 @@
+include <error.h>
+include <imhdr.h>
+include <mach.h>
+include <math.h>
+include <math/curfit.h>
+include <math/iminterp.h>
+include "specfocus.h"
+
+
+# T_SPECFOCUS -- Spectral focusing task
+
+procedure t_specfocus ()
+
+int list # List of images
+pointer fvals # List of focus values
+int dispaxis # Default dispersion axis
+int amin # Lower edge of data along slit
+int amax # Upper edge of data along slit
+int nspec # Number of spectra to subdivide width
+int ndisp # Number of dispersion samples
+int lag # Maximum lag
+real level # Level for width
+bool shifts # Measure shifts?
+int log # Log file descriptor
+
+int i, j, k, l, nimages, npix, nprep
+int aaxis, a1, da, na
+int baxis, b1, db, nb
+int c1, dc, nc
+int l1, dl, nl
+pointer sp, image, sys
+pointer rg, sfs, sf, sfd, sfavg, sfbest, im, mw, data, buf1, buf2
+pointer rng_open(), immap(), imgl2r(), mw_openim()
+int clgeti(), imgeti(), imtopenp(), imtlen(), imtgetim(), nowhite()
+int rng_index(), open()
+real rval, clgetr(), imgetr(), asumr()
+bool ms, clgetb(), streq()
+errchk immap, spf_width
+
+int spf_compare()
+extern spf_compare
+
+begin
+ call smark (sp)
+ call salloc (fvals, SZ_LINE, TY_CHAR)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (sys, SZ_FNAME, TY_CHAR)
+
+ # Get task parameters (except log file)
+ list = imtopenp ("images")
+ call clgstr ("focus", Memc[fvals], SZ_LINE)
+ dispaxis = clgeti ("dispaxis")
+ amin = clgeti ("slit1")
+ amax = clgeti ("slit2")
+ nspec = clgeti ("nspectra")
+ ndisp = clgeti ("ndisp")
+ lag = (clgeti ("corwidth") + 1) / 2
+ level = clgetr ("level")
+ shifts = clgetb ("shifts")
+
+ if (level > 1.)
+ level = level / 100.
+ level = max (0.05, min (0.95, level))
+
+ # Initialize focus values
+ if (nowhite (Memc[fvals], Memc[fvals], SZ_LINE) == 0)
+ call strcpy ("1x1", Memc[fvals], SZ_LINE)
+ iferr (rg = rng_open (Memc[fvals], -MAX_REAL, MAX_REAL, 1.))
+ rg = NULL
+
+ # Allocate array for the image focus data structure pointers
+ nimages = imtlen (list)
+ call malloc (sfs, nimages, TY_POINTER)
+
+ # Accumulate the focus data
+ nimages = 0
+ while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+ im = immap (Memc[image], READ_ONLY, 0)
+ mw = mw_openim (im)
+ call mw_gsystem (mw, Memc[sys], SZ_FNAME)
+ ms = streq (Memc[sys], "multispec")
+ call mw_close (mw)
+
+ # Set the focus value
+ if (rg != NULL) {
+ if (rng_index (rg, nimages+1, rval) == EOF)
+ call error (1, "Focus list ended prematurely")
+ } else
+ rval = imgetr (im, Memc[fvals])
+
+ # Set dispersion and cross dispersion axes
+ if (ms) {
+ baxis = 1
+ aaxis = 2
+
+ # Set sampling across the dispersion axis
+ if (IS_INDEFI (amin))
+ i = 1
+ else
+ i = amin
+ if (IS_INDEFI (amax))
+ j = IM_LEN(im,aaxis)
+ else
+ j = amax
+ a1 = max (1, min (i, j))
+ da = min (IM_LEN(im,aaxis), max (i, j)) - a1 + 1
+ if (da < 1)
+ call error (1, "Error in slit limits")
+ na = da
+ da = da / na
+
+ # Set sampling along the dispersion axis
+ npix = IM_LEN(im,baxis)
+ nb = min (ndisp, npix / 100)
+ db = npix / nb
+ b1 = 1 + (npix - nb * db) / 2
+
+ # Set sampling along the columns and lines
+ c1 = b1; dc = db; nc = nb
+ l1 = a1; dl = da; nl = na
+ } else {
+ iferr (baxis = imgeti (im, "dispaxis"))
+ baxis = dispaxis
+ aaxis = 3 - baxis
+
+ # Set sampling across the dispersion axis
+ if (IS_INDEFI (amin))
+ i = 1
+ else
+ i = amin
+ if (IS_INDEFI (amax))
+ j = IM_LEN(im,aaxis)
+ else
+ j = amax
+ a1 = max (1, min (i, j))
+ da = min (IM_LEN(im,aaxis), max (i, j)) - a1 + 1
+ if (da < 1)
+ call error (1, "Error in slit limits")
+ na = min (nspec, da)
+ da = da / na
+
+ # Set sampling along the dispersion axis
+ npix = IM_LEN(im,baxis)
+ nb = min (ndisp, npix / 100)
+ db = npix / nb
+ b1 = 1 + (npix - nb * db) / 2
+
+ # Set sampling along the columns and lines
+ if (baxis == 1) {
+ c1 = b1; dc = db; nc = nb
+ l1 = a1; dl = da; nl = na
+ } else {
+ c1 = a1; dc = da; nc = na
+ l1 = b1; dl = db; nl = nb
+ }
+ }
+
+ # Check for consistency
+ if (nimages > 0) {
+ if (baxis!=SF_AXIS(sf) ||
+ c1!=SF_X1(sf) || dc!=SF_DX(sf) || nc!=SF_NX(sf) ||
+ l1!=SF_Y1(sf) || dl!=SF_DY(sf) || nl!=SF_NY(sf))
+ call error (1, "Input images have different formats")
+ }
+
+ # Allocate the focus data structure for the image
+ call spf_alloc (sf, Memc[image], rval, level, baxis, npix, na,
+ c1, dc, nc, l1, dl, nl)
+
+ # Get the spectrum samples
+ if (baxis == 1) {
+ do i = 1, na {
+ k = a1 + (i - 1) * da
+ l = k + da - 1
+ data = SF_DATA(sf) + (i - 1) * npix
+ do j = k, l
+ call aaddr (Memr[imgl2r(im,j)], Memr[data], Memr[data],
+ npix)
+ }
+ } else {
+ do j = 1, npix {
+ buf1 = imgl2r (im, j)
+ data = SF_DATA(sf) + j - 1
+ do i = 1, na {
+ k = a1 + (i - 1) * da
+ Memr[data] = asumr (Memr[buf1+k-1], da)
+ data = data + npix
+ }
+ }
+ }
+ Memi[sfs+nimages] = sf
+ nimages = nimages + 1
+
+ call imunmap (im)
+ }
+
+ if (nimages == 0)
+ call error (1, "No input data")
+
+ # Sort the structures
+ call qsort (Memi[sfs], nimages, spf_compare)
+
+ # Allocate structure for the best focus
+ call spf_alloc (sfavg, "Best", INDEF, level, baxis, 0, 0, c1, dc, nc,
+ l1, dl, nl)
+
+ # Compute the correlations and profile width and position
+ nprep = db + 2 * lag
+ call malloc (buf1, nprep, TY_REAL)
+ call malloc (buf2, nprep, TY_REAL)
+ l = (na + 1) / 2
+ do k = 1, nimages {
+ sf = Memi[sfs+k-1]
+ do i = 1, nb {
+ if (baxis == 1)
+ sfd = SFD(sf,i,l)
+ else
+ sfd = SFD(sf,l,i)
+ call spf_prep (Memr[SF_SPEC(sfd)], db, Memr[buf1], nprep)
+ call spf_corr (Memr[buf1], Memr[buf1], nprep, lag,
+ SF_ASI(sfd), SF_POS(sfd), SF_WID(sfd), SF_LEVEL(sf))
+ do j = 1, na {
+ if (j != l) {
+ if (baxis == 1)
+ sfd = SFD(sf,i,j)
+ else
+ sfd = SFD(sf,j,i)
+ call spf_prep (Memr[SF_SPEC(sfd)], db, Memr[buf2],
+ nprep)
+ call spf_corr (Memr[buf2], Memr[buf2], nprep, lag,
+ SF_ASI(sfd), SF_POS(sfd), SF_WID(sfd), SF_LEVEL(sf))
+ if (shifts)
+ call spf_corr (Memr[buf2], Memr[buf1], nprep,
+ lag, SF_ASI(sfd), SF_POS(sfd), rval,
+ SF_LEVEL(sf))
+ }
+ }
+ }
+ }
+ call mfree (buf1, TY_REAL)
+ call mfree (buf2, TY_REAL)
+
+ # Set the averages
+ call spf_fitfocus (Memi[sfs], nimages, sfavg, sfbest)
+
+ # Graph the results
+ call spf_graph (sfavg, sfbest, Memi[sfs], nimages, lag)
+
+ # Log the results
+ call spf_log (sfavg, sfbest, Memi[sfs], nimages, shifts, STDOUT)
+ call clgstr ("logfile", Memc[image], SZ_FNAME)
+ ifnoerr (log = open (Memc[image], APPEND, TEXT_FILE)) {
+ call spf_log (sfavg, sfbest, Memi[sfs], nimages, shifts, log)
+ call close (log)
+ }
+
+ # Finish up
+ do i = 1, nimages
+ call spf_free (Memi[sfs+i-1])
+ call spf_free (sfavg)
+ call mfree (sfs, TY_POINTER)
+ call rng_close (rg)
+ call imtclose (list)
+ call sfree (sp)
+end
+
+
+# SPF_ALLOC -- Allocate a focus data structure for an image
+
+procedure spf_alloc (sf, image, focus, level, axis, ndisp, nspec,
+ x1, dx, nx, y1, dy, ny)
+
+pointer sf # Image focus data structure
+char image[ARB] # Image name
+real focus # Focus value
+real level # Level for width
+int axis # Dispersion axis
+int ndisp # Number of pixels (along dispersion)
+int nspec # Number of spectra (across dispersion)
+int x1, dx, nx # X sampling
+int y1, dy, ny # Y sampling
+
+int i, j
+pointer data, sfd
+
+begin
+ call calloc (sf, LEN_SF, TY_STRUCT)
+
+ call strcpy (image, SF_IMAGE(sf), SZ_SFFNAME)
+ SF_FOCUS(sf) = focus
+ SF_WIDTH(sf) = INDEF
+ SF_LEVEL(sf) = level
+ SF_AXIS(sf) = axis
+ SF_X1(sf) = x1
+ SF_DX(sf) = dx
+ SF_NX(sf) = nx
+ SF_Y1(sf) = y1
+ SF_DY(sf) = dy
+ SF_NY(sf) = ny
+ call malloc (SF_SFD(sf), nx*ny, TY_POINTER)
+ SF_NSFD(sf) = nx*ny
+ SF_NPIX(sf) = ndisp
+ if (ndisp > 0)
+ call calloc (SF_DATA(sf), ndisp * nspec, TY_REAL)
+
+ data = SF_DATA(sf)
+ do j = 1, ny {
+ do i = 1, nx {
+ call calloc (sfd, LEN_SFD, TY_STRUCT)
+ SFD(sf,i,j) = sfd
+ SF_X(sfd) = x1 + (i - 0.5) * dx
+ SF_Y(sfd) = y1 + (j - 0.5) * dy
+ if (ndisp > 0) {
+ if (axis == 1)
+ SF_SPEC(sfd) = data + (j-1)*ndisp + (i-1)*dx + x1-1
+ else
+ SF_SPEC(sfd) = data + (i-1)*ndisp + (j-1)*dy + y1-1
+ }
+ call asiinit (SF_ASI(sfd), II_SPLINE3)
+ SF_FOC(sfd) = focus
+ SF_WID(sfd) = INDEF
+ SF_POS(sfd) = INDEF
+ SF_DEL(sfd) = NO
+ }
+ }
+end
+
+
+# SPF_FREE -- Free a focus image data structure
+
+procedure spf_free (sf)
+
+pointer sf # Image focus data structure
+
+int i
+pointer sfd
+
+begin
+ do i = 1, SF_NSFD(sf) {
+ sfd = SFD(sf,i,1)
+ call asifree (SF_ASI(sfd))
+ call mfree (sfd, TY_STRUCT)
+ }
+ call mfree (SF_DATA(sf), TY_REAL)
+ call mfree (SF_SFD(sf), TY_POINTER)
+ call mfree (sf, TY_STRUCT)
+end
+
+
+# SPF_PREP -- Prepare spectra for correlation: fit continuum, subtract, taper
+
+procedure spf_prep (in, nin, out, nout)
+
+real in[nin] # Input spectrum
+int nin # Number of pixels in input spectrum
+real out[nout] # Output spectrum
+int nout # Number of pixels output spectrum (nin+2*lag)
+
+int i, lag
+real cveval()
+pointer sp, x, w, ic, cv
+
+begin
+ call smark (sp)
+ call salloc (x, nin, TY_REAL)
+ call salloc (w, nin, TY_REAL)
+
+ call ic_open (ic)
+ call ic_pstr (ic, "function", "chebyshev")
+ call ic_puti (ic, "order", 3)
+ call ic_putr (ic, "low", 3.)
+ call ic_putr (ic, "high", 1.)
+ call ic_puti (ic, "niterate", 5)
+ call ic_putr (ic, "grow", 1.)
+ call ic_putr (ic, "xmin", 1.)
+ call ic_putr (ic, "xmax", real(nin))
+
+ do i = 1, nin {
+ Memr[x+i-1] = i
+ Memr[w+i-1] = 1
+ }
+ call ic_fit (ic, cv, Memr[x], in, Memr[w], nin, YES, YES, YES, YES)
+
+ lag = (nout - nin) / 2
+ do i = 1-lag, 0
+ out[i+lag] = 0.
+ do i = 1, lag-1
+ out[i+lag] = (1-cos (PI*i/lag))/2 * (in[i] - cveval (cv, real(i)))
+ do i = lag, nin-lag+1
+ out[i+lag] = (in[i] - cveval (cv, real(i)))
+ do i = nin-lag+2, nin
+ out[i+lag] = (1-cos (PI*(nin+1-i)/lag))/2 *
+ (in[i] - cveval (cv, real(i)))
+ do i = nin+1, nin+lag
+ out[i+lag] = 0.
+
+ call cvfree (cv)
+ call ic_closer (ic)
+ call sfree (sp)
+end
+
+
+# SPF_CORR -- Correlate spectra, fit profile, and measure center/width
+
+procedure spf_corr (spec1, spec2, npix, lag, asi, center, width, level)
+
+real spec1[npix] # First spectrum
+real spec2[npix] # Second spectrum
+int npix # Number of pixels in spectra
+int lag # Maximum correlation lag
+pointer asi # Pointer to correlation profile interpolator
+real center # Center of profile
+real width # Width of profile
+real level # Level at which width is determined
+
+int i, j, nprof
+real x, p, pmin, pmax, asieval()
+pointer sp, prof
+
+begin
+ nprof = 2 * lag + 1
+
+ call smark (sp)
+ call salloc (prof, nprof, TY_REAL)
+
+ do j = -lag, lag {
+ p = 0.
+ do i = 1+lag, npix-lag
+ p = p + spec1[i] * spec2[i-j]
+ Memr[prof+j+lag] = p
+ }
+
+ # Fit interpolator
+ call asifit (asi, Memr[prof], nprof)
+
+ # Find the minimum and maximum
+ center = 1.
+ pmin = asieval (asi, 1.)
+ pmax = pmin
+ for (x=1; x<=nprof; x=x+.01) {
+ p = asieval (asi, x)
+ if (p < pmin)
+ pmin = p
+ if (p > pmax) {
+ pmax = p
+ center = x
+ }
+ }
+
+ # Normalize
+ pmax = pmax - pmin
+ do i = 0, nprof-1
+ Memr[prof+i] = (Memr[prof+i] - pmin) / pmax
+
+ call asifit (asi, Memr[prof], nprof)
+
+ # Find the equal flux points
+ for (x=center; x>=1 && asieval (asi,x)>level; x=x-0.01)
+ ;
+ width = x
+ for (x=center; x<=nprof && asieval (asi,x)>level; x=x+0.01)
+ ;
+ width = (x - width - 0.01) / sqrt (2.)
+ center = center - lag - 1
+
+ call sfree (sp)
+end
+
+
+# SPF_FITFOCUS -- Find the best focus at each sample and the average over all
+# samples.
+
+procedure spf_fitfocus (sfs, nimages, sfavg, sfbest)
+
+pointer sfs[nimages] #I Images
+int nimages #I Number of images
+pointer sfavg #U Average image
+pointer sfbest #U Best focus image
+
+int i, j, n, jmin, nims
+pointer sp, x, y, z, sfd
+real focus, fwhm, pos, foc
+bool fp_equalr()
+
+define avg_ 10
+
+begin
+ call smark (sp)
+ call salloc (x, nimages, TY_REAL)
+ call salloc (y, nimages, TY_REAL)
+ call salloc (z, nimages, TY_REAL)
+
+ do i = 1, SF_NSFD(sfavg) {
+ # Collect the focus values
+ nims = 0
+ do j = 1, nimages {
+ sfd = SFD(sfs[j],i,1)
+ if (SF_DEL(sfd) == NO) {
+ Memr[x+nims] = SF_FOC(sfd)
+ Memr[y+nims] = SF_WID(sfd)
+ Memr[z+nims] = SF_POS(sfd)
+ nims = nims + 1
+ }
+ }
+ sfd = SFD(sfavg,i,1)
+
+ # Take the smallest width at each unique focus.
+ if (nims > 0) {
+ call xt_sort3 (Memr[x], Memr[y], Memr[z], nims)
+ n = 0
+ do j = 1, nims-1 {
+ if (fp_equalr (Memr[x+n], Memr[x+j])) {
+ if (Memr[y+n] > Memr[y+j]) {
+ Memr[y+n] = Memr[y+j]
+ Memr[z+n] = Memr[z+j]
+ }
+ } else {
+ n = n + 1
+ Memr[x+n] = Memr[x+j]
+ Memr[y+n] = Memr[y+j]
+ Memr[z+n] = Memr[z+j]
+ }
+ }
+
+ # Find the minimum width
+ jmin = 0
+ do j = 1, n
+ if (Memr[y+j] < Memr[y+jmin])
+ jmin = j
+
+ # Use parabolic interpolation to find the best focus
+ if (jmin == 0 || jmin == n) {
+ focus = Memr[x+jmin]
+ fwhm = Memr[y+jmin]
+ pos = Memr[z+jmin]
+ } else
+ call spf_parab (Memr[x+jmin-1], Memr[y+jmin-1],
+ Memr[z+jmin-1], focus, fwhm, pos)
+
+ SF_FOC(sfd) = focus
+ SF_WID(sfd) = fwhm
+ SF_POS(sfd) = pos
+ SF_DEL(sfd) = NO
+ } else {
+ SF_FOC(sfd) = INDEF
+ SF_WID(sfd) = INDEF
+ SF_POS(sfd) = INDEF
+ SF_DEL(sfd) = YES
+ }
+ }
+
+ call sfree (sp)
+
+avg_
+ # Set the averages over all samples
+ n = 0
+ focus = 0.
+ fwhm = 0.
+ do i = 1, SF_NSFD(sfavg) {
+ sfd = SFD(sfavg,i,1)
+ if (SF_DEL(sfd) == NO) {
+ focus = focus + SF_FOC(sfd)
+ fwhm = fwhm + SF_WID(sfd)
+ n = n + 1
+ }
+ }
+
+ if (n > 0) {
+ SF_FOCUS(sfavg) = focus / n
+ SF_WIDTH(sfavg) = fwhm / n
+ } else {
+ SF_FOCUS(sfavg) = INDEF
+ SF_WIDTH(sfavg) = INDEF
+ }
+
+ do j = 1, nimages {
+ n = 0
+ focus = 0.
+ fwhm = 0.
+ do i = 1, SF_NSFD(sfs[j]) {
+ sfd = SFD(sfs[j],i,1)
+ if (SF_DEL(sfd) == NO) {
+ fwhm = fwhm + SF_WID(sfd)
+ n = n + 1
+ }
+ }
+
+ if (n > 0)
+ SF_WIDTH(sfs[j]) = fwhm / n
+ else
+ SF_WIDTH(sfs[j]) = INDEF
+ }
+
+ # Set the best focus image
+ sfbest = NULL
+ focus = SF_FOCUS(sfavg)
+ if (!IS_INDEF(focus)) {
+ pos = MAX_REAL
+ do j = 1, nimages {
+ foc = SF_FOCUS(sfs[j])
+ if (!IS_INDEF(foc)) {
+ fwhm = abs (focus - foc)
+ if (fwhm < pos) {
+ pos = fwhm
+ sfbest = sfs[j]
+ }
+ }
+ }
+ }
+end
+
+
+# SPF_PARAB -- Find the minimum of a parabolic fit to three points.
+
+procedure spf_parab (x, y, z, xmin, ymin, zmin)
+
+real x[3]
+real y[3]
+real z[3]
+real xmin
+real ymin
+real zmin
+
+real x12, x13, x23, x212, x213, x223, y12, y13, y23, a, b, c
+
+begin
+ x12 = x[1] - x[2]
+ x13 = x[1] - x[3]
+ x23 = x[2] - x[3]
+ x212 = x[1] * x[1] - x[2] * x[2]
+ x213 = x[1] * x[1] - x[3] * x[3]
+ x223 = x[2] * x[2] - x[3] * x[3]
+ y12 = y[1] - y[2]
+ y13 = y[1] - y[3]
+ y23 = y[2] - y[3]
+ c = (y13 - y23 * x13 / x23) / (x213 - x223 * x13 / x23)
+ b = (y23 - c * x223) / x23
+ a = y[3] - b * x[3] - c * x[3] * x[3]
+ xmin = -b / (2 * c)
+ ymin = a + b * xmin + c * xmin * xmin
+
+ if (xmin < x[2])
+ zmin = z[2] + (z[1] - z[2]) / (x[1] - x[2]) * (xmin - x[2])
+ else
+ zmin = z[2] + (z[3] - z[2]) / (x[3] - x[2]) * (xmin - x[2])
+end
+
+
+# SPF_COMPARE -- Compare two structures by focus values
+
+int procedure spf_compare (sf1, sf2)
+
+pointer sf1, sf2 # Structures to be compared.
+
+begin
+ if (SF_FOCUS[sf1] < SF_FOCUS[sf2])
+ return (-1)
+ else if (SF_FOCUS[sf1] > SF_FOCUS[sf2])
+ return (1)
+ else
+ return (0)
+end
+
+
+# SPF_LOG -- Print log of results
+
+procedure spf_log (sfavg, sfbest, sfs, nimages, shifts, log)
+
+pointer sfavg # Average image
+pointer sfbest # Best focus image
+pointer sfs[nimages] # Images
+int nimages # Number of images
+bool shifts # Measure shifts?
+int log # Log file descriptor
+
+int i, j
+pointer sp, str, sfd
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ call sysid (Memc[str], SZ_LINE)
+ call fprintf (log, "SPECFOCUS: %s\n")
+ call pargstr (Memc[str])
+
+ call fprintf (log,
+ " Best average focus at %g with average width of %.2f at %d%% of peak\n\n")
+ call pargr (SF_FOCUS(sfavg))
+ call pargr (SF_WIDTH(sfavg))
+ call pargr (100 * SF_LEVEL(sfavg))
+
+ call fprintf (log, " -- Average Over All Samples\n\n")
+ call fprintf (log, "\t%25wImage Focus Width\n")
+ do i = 1, nimages {
+ call fprintf (log, "\t%30s %5.3g %5.2f\n")
+ call pargstr (SF_IMAGE(sfs[i]))
+ call pargr (SF_FOCUS(sfs[i]))
+ call pargr (SF_WIDTH(sfs[i]))
+ }
+ call fprintf (log, "\n")
+
+ call fprintf (log, " -- Image %s at Focus %g --\n")
+ call pargstr (SF_IMAGE(sfbest))
+ call pargr (SF_FOCUS(sfbest))
+
+ if (SF_NSFD(sfbest) > 1) {
+ call fprintf (log, "\n\n\tWidth at %d%% of Peak:\n")
+ call pargr (100 * SF_LEVEL(sfavg))
+ call fprintf (log, "\n\t%9w Columns\n\t%9w ")
+ do i = 1, SF_NX(sfbest) {
+ call fprintf (log, " %4d-%-4d")
+ call pargi (SF_X1(sfbest)+(i-1)*SF_DX(sfbest))
+ call pargi (SF_X1(sfbest)+i*SF_DX(sfbest)-1)
+ }
+ call fprintf (log, "\n\t Lines +")
+ do i = 1, SF_NX(sfbest)
+ call fprintf (log, "-----------")
+ do j = 1, SF_NY(sfbest) {
+ if (SF_DY(sfbest) > 1.) {
+ call fprintf (log, "\n\t%4d-%-4d |")
+ call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+ call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+ } else {
+ call fprintf (log, "\n\t%9d |")
+ call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+ call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+ }
+ do i = 1, SF_NX(sfbest) {
+ sfd = SFD(sfbest,i,j)
+ call fprintf (log, " %5.2f ")
+ call pargr (SF_WID(sfd))
+ }
+ }
+ if (shifts) {
+ call fprintf (log,
+ "\n\n\tPosition Shifts Relative To Central Sample:\n")
+ call fprintf (log, "\n\t%9w Columns\n\t%9w ")
+ do i = 1, SF_NX(sfbest) {
+ call fprintf (log, " %4d-%-4d")
+ call pargi (SF_X1(sfbest)+(i-1)*SF_DX(sfbest))
+ call pargi (SF_X1(sfbest)+i*SF_DX(sfbest)-1)
+ }
+ call fprintf (log, "\n\t Lines +")
+ do i = 1, SF_NX(sfbest)
+ call fprintf (log, "-----------")
+ do j = 1, SF_NY(sfbest) {
+ call fprintf (log, "\n\t%4d-%-4d |")
+ call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+ call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+ do i = 1, SF_NX(sfbest) {
+ sfd = SFD(sfbest,i,j)
+ call fprintf (log, " %5.2f ")
+ call pargr (SF_POS(sfd))
+ }
+ }
+ }
+ }
+ call fprintf (log, "\n\n")
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:52:55 -0400