aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/t_tweak.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/onedspec/t_tweak.x')
-rw-r--r--noao/onedspec/t_tweak.x1352
1 files changed, 1352 insertions, 0 deletions
diff --git a/noao/onedspec/t_tweak.x b/noao/onedspec/t_tweak.x
new file mode 100644
index 00000000..e50dccea
--- /dev/null
+++ b/noao/onedspec/t_tweak.x
@@ -0,0 +1,1352 @@
+include <error.h>
+include <gset.h>
+include <imset.h>
+include <imhdr.h>
+include <math.h>
+include <math/iminterp.h>
+include <pkg/gtools.h>
+include <smw.h>
+include <units.h>
+include <pkg/xtanswer.h>
+
+# Tweak data object definitions.
+define TWK_SLEN 999 # Length of sample region string
+define TWK_LEN 580 # Length of data object
+
+define TWK_TYPE Memc[P2C($1)] # Tweak type (maxchars=19)
+define TWK_SH Memi[$1+11] # Spectrum pointer
+define TWK_CAL Memi[$1+12] # Calibration pointer
+define TWK_WAVE Memi[$1+13] # Pointer to wavelengths
+define TWK_SPEC Memi[$1+14] # Pointer to calibrated spectrum
+define TWK_SHIFT Memr[P2R($1+15)] # Shift
+define TWK_DSHIFT Memr[P2R($1+16)] # Shift step
+define TWK_SCALE Memr[P2R($1+17)] # Scaling factor
+define TWK_DSCALE Memr[P2R($1+18)] # Scaling factor step
+define TWK_RG Memi[$1+19] # Range pointer
+define TWK_RMS Memr[P2R($1+20)] # RMS in sample regions
+define TWK_OFFSET Memr[P2R($1+21)] # Offset in graphs
+define TWK_BOX Memi[$1+22] # Boxcar smoothing size
+define TWK_THRESH Memr[P2R($1+23)] # Calibration threshold
+define TWK_SAMPLE Memc[P2C($1+30)] # Sample regions (maxchars=999)
+define TWK_HELP Memc[P2C($1+530)] # Help file (maxchars=99)
+
+# Tweak types.
+define SKYTWEAK 1 # Sky subtraction
+define TELLURIC 2 # Telluric division
+
+# Secondary graph types.
+define GNONE 0 # No graph
+define GCAL 1 # Graph calibration spectrum
+define GDATA 2 # Graph data spectrum
+
+
+# T_SKYTWEAK -- Sky subtract spectra with shift and scale tweaking.
+# The sky calibration spectra are scaled and shifted to best subtract
+# sky features. Automatic and interactive methods are provided.
+
+procedure t_skytweak ()
+
+pointer twk # TWK data object
+
+begin
+ call malloc (twk, TWK_LEN, TY_STRUCT)
+ call strcpy ("SKYTWEAK", TWK_TYPE(twk), 19)
+ call strcpy ("onedspec$doc/skytweak.key", TWK_HELP(twk), 99)
+ call tweak (twk)
+ call mfree (twk, TY_STRUCT)
+end
+
+
+# T_TELLURIC -- Correct spectra for telluric features.
+# The telluric calibration spectra are scaled by raising to a power (Beers law)
+# and shifted to best remove telluric features. Automatic and interactive
+# methods are provided.
+
+procedure t_telluric ()
+
+pointer twk # TWK data object
+
+begin
+ call malloc (twk, TWK_LEN, TY_STRUCT)
+ call strcpy ("TELLURIC", TWK_TYPE(twk), 19)
+ call strcpy ("onedspec$doc/telluric.key", TWK_HELP(twk), 99)
+ call tweak (twk)
+ call mfree (twk, TY_STRUCT)
+end
+
+
+# TWEAK -- Tweak spectra for shift and scale before applying a correction.
+# This procedure implements both sky subtraction and telluric division.
+
+procedure tweak (twk)
+
+pointer twk #I TWK data object
+
+pointer inlist # Input list
+pointer outlist # Output list
+pointer callist # Calibration list
+bool xcorr # Cross correlate for initial shift
+bool tweakrms # Tweak to minimize RMS?
+bool ignoreaps # Ignore aperture numbers?
+int lag # Cross correlation lag
+bool interactive # Interactive?
+
+int i, j, k, n, nout, ncal, answer
+real shift, scale, fcor, ical, mean
+pointer sp, input, output, calname, temp
+pointer in, smw, sh, out, pcal, cal, x, y, data, tmp
+
+int clgeti(), imtgetim(), imtlen()
+bool clgetb(), streq()
+real clgetr(), asieval()
+double shdr_wl(), shdr_lw()
+pointer imtopenp(), immap(), smw_openim(), impl3r(), imgl3r()
+errchk immap, smw_openim, shdr_open, twk_gcal, twk_tweak, impl3r, imgl3r
+
+begin
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (output, SZ_FNAME, TY_CHAR)
+ call salloc (calname, SZ_FNAME, TY_CHAR)
+ call salloc (temp, SZ_LINE, TY_CHAR)
+ call malloc (TWK_WAVE(twk), 1000, TY_DOUBLE)
+ call malloc (TWK_SPEC(twk), 1000, TY_REAL)
+
+ # Get task parameters.
+ inlist = imtopenp ("input")
+ outlist = imtopenp ("output")
+ callist = imtopenp ("cal")
+ ignoreaps = clgetb ("ignoreaps")
+ if (TWK_TYPE(twk) == 'T')
+ TWK_THRESH(twk) = clgetr ("threshold")
+ TWK_SHIFT(twk) = clgetr ("shift")
+ TWK_SCALE(twk) = clgetr ("scale")
+ xcorr = clgetb ("xcorr")
+ tweakrms = clgetb ("tweakrms")
+ interactive = clgetb ("interactive")
+ if (interactive)
+ answer = YES
+ else
+ answer = ALWAYSNO
+ call clgstr ("sample", TWK_SAMPLE(twk), TWK_SLEN)
+ lag = clgeti ("lag")
+ TWK_DSHIFT(twk) = max (0., clgetr ("dshift"))
+ TWK_DSCALE(twk) = max (0., min (0.99, clgetr ("dscale")))
+ TWK_OFFSET(twk) = clgetr ("offset")
+ TWK_BOX(twk) = max (1, clgeti ("smooth"))
+
+ if (imtlen (inlist) != imtlen (callist) && imtlen (callist) != 1) {
+ call imtclose (inlist)
+ call imtclose (outlist)
+ call imtclose (callist)
+ call sfree (sp)
+ call error (1, "Image lists do not match")
+ }
+
+ # Loop over all input images.
+ sh = NULL
+ ncal = 0
+ while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+ if (imtgetim (callist, Memc[calname], SZ_FNAME) != EOF) {
+ if (ncal > 0) {
+ do i = 0, ncal-1 {
+ cal = Memi[pcal+i]
+ call asifree (IM(cal))
+ call smw_close (MW(cal))
+ call shdr_close (cal)
+ }
+ call mfree (pcal, TY_POINTER)
+ ncal = 0
+ }
+ }
+
+ in = NULL; smw = NULL; sh = NULL; out = NULL
+ iferr {
+ # Set output image. Use a temporary image when output=input.
+ if (imtlen (outlist) > 0) {
+ if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF)
+ break
+ } else
+ call strcpy (Memc[input], Memc[output], SZ_FNAME)
+
+ # Map the input image.
+ tmp = immap (Memc[input], READ_ONLY, 0); in = tmp
+ tmp = smw_openim (in); smw = tmp
+ if (smw == SMW_ND)
+ call error (1, "NDSPEC data not supported")
+ call shdr_open (in, smw, 1, 1, INDEFI, SHHDR, sh)
+
+ # Map the output image.
+ if (streq (Memc[input], Memc[output]))
+ call mktemp ("temp", Memc[temp], SZ_LINE)
+ else
+ call strcpy (Memc[output], Memc[temp], SZ_LINE)
+ tmp = immap (Memc[temp], NEW_COPY, in); out = tmp
+ if (IM_PIXTYPE(out) != TY_DOUBLE)
+ IM_PIXTYPE(out) = TY_REAL
+
+ # Determine airmass if needed.
+ if (TWK_TYPE(twk) == 'T') {
+ if (IS_INDEF(AM(sh))) {
+ call printf ("%s: ")
+ call pargstr (Memc[input])
+ call flush (STDOUT)
+ AM(sh) = clgetr ("airmass")
+ }
+ }
+
+ # Calibrate all spectra in the image.
+ # Only the first band is done.
+ do i = 1, IM_LEN(in,2) {
+
+ # Get the spectra.
+ call shdr_open (in, smw, i, 1, INDEFI, SHDATA, sh)
+ call realloc (TWK_WAVE(twk), SN(sh), TY_DOUBLE)
+ x = TWK_WAVE(twk)
+ do k = 1, SN(sh) {
+ Memd[x] = shdr_lw (sh, double(k))
+ x = x + 1
+ }
+ if (ignoreaps)
+ call twk_gcal (twk, Memc[calname], INDEFI,
+ pcal, ncal, cal)
+ else
+ call twk_gcal (twk, Memc[calname], AP(sh),
+ pcal, ncal, cal)
+
+ # Determine the shift and scale.
+ TWK_SH(twk) = sh
+ TWK_CAL(twk) = cal
+ call realloc (SY(cal), SN(sh), TY_REAL)
+ call realloc (TWK_SPEC(twk), SN(sh), TY_REAL)
+ if (answer == NO || answer == YES) {
+ call printf ("%s%s: ")
+ call pargstr (IMNAME(sh))
+ call pargstr (IMSEC(sh))
+ call flush (STDOUT)
+ call xt_clanswer ("answer", answer)
+ }
+ if (answer == YES || answer == ALWAYSYES)
+ interactive = true
+ else
+ interactive = false
+ call twk_tweak (twk, xcorr, tweakrms, interactive, lag)
+ shift = TWK_SHIFT(twk)
+ if (TWK_TYPE(twk) == 'T')
+ scale = TWK_SCALE(twk) * AM(sh) / AM(cal)
+ else
+ scale = TWK_SCALE(twk)
+
+ # Calibrate the output spectrum.
+ nout = 0
+ mean = 0.
+ x = TWK_WAVE(twk)
+ y = SY(sh)
+ n = SN(sh)
+ data = impl3r (out, i, 1)
+ do k = 1, n {
+ ical = shdr_wl (cal, Memd[x]) + shift
+ if (ical < 1. || ical > SN(cal)) {
+ if (ical < 0.5 || ical > SN(cal) + 0.5)
+ nout = nout + 1
+ ical = max (1., min (real(SN(cal)), ical))
+ }
+ if (TWK_TYPE(twk) == 'T') {
+ fcor = max (TWK_THRESH(twk),
+ asieval (IM(cal),ical)) ** scale
+ Memr[data] = Memr[y] / fcor
+ mean = mean + fcor
+ } else {
+ fcor = asieval (IM(cal),ical) * scale
+ Memr[data] = Memr[y] - fcor
+ }
+ x = x + 1
+ y = y + 1
+ data = data + 1
+ }
+ mean = mean / n
+ if (TWK_TYPE(twk) == 'T')
+ call amulkr (Memr[data-n], mean, Memr[data-n], n)
+ do k = n+1, IM_LEN(out,1) {
+ Memr[data] = 0
+ data = data + 1
+ }
+
+ # Log the results.
+ if (i == 1) {
+ call printf ("%s:\n Output: %s - %s\n")
+ call pargstr (TWK_TYPE(twk))
+ call pargstr (Memc[output])
+ call pargstr (IM_TITLE(out))
+ }
+ call printf (" Input: %s%s - %s\n")
+ call pargstr (IMNAME(sh))
+ call pargstr (IMSEC(sh))
+ call pargstr (TITLE(sh))
+ call printf (" Calibration: %s%s - %s\n")
+ call pargstr (IMNAME(cal))
+ call pargstr (IMSEC(cal))
+ call pargstr (TITLE(cal))
+ call printf (" Tweak: shift = %.2f, scale = %.3f")
+ call pargr (shift)
+ call pargr (TWK_SCALE(twk))
+ if (TWK_TYPE(twk) == 'T') {
+ call printf (", normalization = %.4g\n")
+ call pargr (mean)
+ } else
+ call printf ("\n")
+ if (nout > 0) {
+ call printf (
+ " WARNING: %d pixels outside of calibration limits\n")
+ call pargi (nout)
+ }
+ call flush (STDOUT)
+ }
+ do j = 2, IM_LEN(in,3) {
+ do i = 1, IM_LEN(in,2) {
+ y = imgl3r (in, i, j)
+ data = impl3r (out, i, j)
+ call amovr (Memr[y], Memr[data], IM_LEN(out,1))
+ }
+ }
+ } then {
+ call erract (EA_WARN)
+ if (out != NULL) {
+ call imunmap (out)
+ if (!streq (Memc[input], Memc[output]))
+ call imdelete (Memc[output])
+ }
+ }
+
+ # Finish up this image.
+ if (in != NULL)
+ call imunmap (in)
+ if (smw != NULL) {
+ call smw_close (smw)
+ if (sh != NULL)
+ MW(sh) = NULL
+ }
+ if (out != NULL) {
+ call imunmap (out)
+ if (streq (Memc[input], Memc[output])) {
+ call imdelete (Memc[input])
+ call imrename (Memc[temp], Memc[output])
+ }
+ }
+ }
+
+ # Finish up.
+ if (ncal > 0) {
+ do i = 0, ncal-1 {
+ cal = Memi[pcal+i]
+ call asifree (IM(cal))
+ call smw_close (MW(cal))
+ call shdr_close (cal)
+ }
+ call mfree (pcal, TY_POINTER)
+ }
+ if (sh != NULL)
+ call shdr_close (sh)
+ call imtclose (inlist)
+ call imtclose (outlist)
+ call imtclose (callist)
+ call mfree (TWK_SPEC(twk), TY_REAL)
+ call sfree (sp)
+end
+
+
+# TWK_GCAL -- Get calibration data
+# An interpolation function is fit and stored in the image pointer field.
+# For efficiency the calibration data is saved by aperture so that additional
+# calls simply return the data pointer.
+
+procedure twk_gcal (twk, calname, ap, pcal, ncal, cal)
+
+pointer twk # TWK data object
+char calname[ARB] # Calibration image name
+int ap # Aperture
+pointer pcal # Pointer to cal data
+int ncal # Number of active cal data structures
+pointer cal # Calibration data structure
+
+int i, clgwrd()
+pointer sp, str, im, smw, immap(), smw_openim()
+real clgetr()
+errchk immap, smw_openim, shdr_open, asifit
+
+begin
+ # Check for previously saved calibration
+ for (i=0; i<ncal; i=i+1) {
+ cal = Memi[pcal+i]
+ if (AP(cal) == ap)
+ return
+ }
+
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Allocate space for a new data pointer and get the calibration data.
+
+ if (ncal == 0)
+ call malloc (pcal, 10, TY_POINTER)
+ else if (mod (ncal, 10) == 0)
+ call realloc (pcal, ncal+10, TY_POINTER)
+
+ im = immap (calname, READ_ONLY, 0)
+ smw = smw_openim (im)
+ cal = NULL
+ call shdr_open (im, smw, 1, 1, ap, SHDATA, cal)
+ AP(cal) = ap
+ Memi[pcal+ncal] = cal
+ ncal = ncal + 1
+ call imunmap (im)
+
+ call asiinit (im, clgwrd ("interp", Memc[str], SZ_FNAME,II_FUNCTIONS))
+ call asifit (im, Memr[SY(cal)], SN(cal))
+ IM(cal) = im
+
+ # Determine airmass if needed.
+ if (TWK_TYPE(twk) == 'T') {
+ if (IS_INDEF(AM(cal))) {
+ call printf ("%s: ")
+ call pargstr (calname)
+ call flush (STDOUT)
+ AM(cal) = clgetr ("airmass")
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# TWK_TWEAK -- Determine the shift and scale using automatic and interactive
+# methods.
+
+procedure twk_tweak (twk, xcorr, tweakrms, interactive, lag)
+
+pointer twk #I TWK data object
+bool xcorr #I Cross correlate for shift
+bool tweakrms #I Tweak by minimizing RMS?
+bool interactive #I Interactive fitting?
+int lag #I Cross correlation lag
+
+int i, n, nlag
+real ical, asieval()
+double shdr_wl()
+pointer sh, cal, rg, asi, x, y, rg_xrangesd()
+errchk twk_rmsmin, twk_fit
+
+begin
+ sh = TWK_SH(twk)
+ cal = TWK_CAL(twk)
+
+ # Set ranges.
+ rg = rg_xrangesd (TWK_SAMPLE(twk), Memd[TWK_WAVE(twk)], SN(sh))
+ call rg_order (rg)
+ call rg_merge (rg)
+ TWK_RG(twk) = rg
+
+ # Cross correlate for shift.
+ if (xcorr && lag > 0) {
+ n = SN(sh)
+ nlag = n + 2 * lag
+ call malloc (x, nlag, TY_REAL)
+ call malloc (y, nlag, TY_REAL)
+
+ do i = 0, n-1 {
+ ical = max (1D0, min (double(SN(cal)),
+ shdr_wl (cal, Memd[TWK_WAVE(twk)+i])))
+ Memr[y+i] = asieval (IM(cal), ical)
+ }
+
+ call twk_prep (Memr[y], n, Memr[x], nlag)
+ call twk_prep (Memr[SY(sh)], n, Memr[y], nlag)
+ call twk_xcorr (Memr[x], Memr[y], i, rg, lag, asi, TWK_SHIFT(twk),
+ ical, 0.5)
+ call asifree (asi)
+ call mfree (x, TY_REAL)
+ call mfree (y, TY_REAL)
+ }
+
+ # Tweak by minimizing RMS.
+ if (tweakrms)
+ call twk_rmsmin (twk)
+
+ # Do interactive step.
+ if (interactive)
+ call twk_fit (twk)
+
+ call rg_free (TWK_RG(twk))
+end
+
+
+# TWK_PREP -- Prepare spectra for correlation: fit continuum, subtract, taper
+
+procedure twk_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
+
+
+# TWK_XCORR -- Correlate spectra, fit profile, and measure center/width
+
+procedure twk_xcorr (spec1, spec2, npix, rg, lag, asi, center, width, level)
+
+real spec1[npix] # First spectrum
+real spec2[npix] # Second spectrum
+int npix # Number of pixels in spectra
+pointer rg # Ranges
+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, k, n, ishift, nprof, rg_inrange()
+real x, p, pmin, pmax, asieval()
+pointer sp, prof
+
+begin
+ nprof = 2 * lag + 1
+
+ call smark (sp)
+ call salloc (prof, nprof, TY_REAL)
+
+ ishift = nint (center)
+ n = 0
+ do j = -lag, lag {
+ p = 0.
+ do i = 1+lag, npix-lag {
+ if (rg_inrange (rg, i-lag) == NO)
+ next
+ k = i - j - ishift
+ if (k < 1 || k > npix)
+ next
+ p = p + spec1[i] * spec2[k]
+ n = n + 1
+ }
+ Memr[prof+j+lag] = p
+ }
+ if (n < 10 * nprof) {
+ call sfree (sp)
+ return
+ }
+
+ # Fit interpolator
+ call asiinit (asi, II_SPLINE3)
+ 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 + ishift
+
+ call sfree (sp)
+end
+
+
+# TWK_RMSMIN -- Tweak shift and scale to minimize the RMS.
+# This changes the shift and scale parameters but not the step.
+
+procedure twk_rmsmin (twk)
+
+pointer twk #I TWK data object
+
+int i
+real lastshift, lastscale
+errchk twk_ashift, twk_ascale
+
+begin
+ lastshift = INDEFR
+ lastscale = INDEFR
+ do i = 1, 2 {
+ if (TWK_SHIFT(twk) == lastshift && TWK_SCALE(twk) == lastscale)
+ break
+ lastshift = TWK_SHIFT(twk)
+ call twk_ashift (twk)
+
+ if (TWK_SHIFT(twk) == lastshift && TWK_SCALE(twk) == lastscale)
+ break
+ lastscale = TWK_SCALE(twk)
+ call twk_ascale (twk)
+ }
+end
+
+
+# TWK_ASCALE -- Automatically determine scale by minimizing RMS.
+
+procedure twk_ascale (twk)
+
+pointer twk #I TWK data object
+
+int i
+real shift, oscale, dscale, lastscale, scale[3], rms[3]
+errchk twk_spec
+
+begin
+ dscale = TWK_DSCALE(twk)
+ if (dscale == 0.)
+ return
+ oscale = TWK_SCALE(twk)
+ shift = TWK_SHIFT(twk)
+ do i = 1, 3 {
+ scale[i] = (1 - (i - 2) * dscale) * oscale
+ call twk_spec (twk, shift, scale[i])
+ rms[i] = TWK_RMS(twk)
+ lastscale = TWK_SCALE(twk)
+ }
+ while (dscale > 0.01) {
+ if (scale[1] / oscale < 0.5 || scale[3] / oscale > 2.) {
+ TWK_SCALE(twk) = oscale
+ break
+ }
+ if (rms[1] < rms[2]) {
+ scale[3] = scale[2]
+ scale[2] = scale[1]
+ scale[1] = (1 - dscale) * scale[2]
+ rms[3] = rms[2]
+ rms[2] = rms[1]
+ call twk_spec (twk, shift, scale[1])
+ rms[1] = TWK_RMS(twk)
+ lastscale = TWK_SCALE(twk)
+ } else if (rms[3] < rms[2]) {
+ scale[1] = scale[2]
+ scale[2] = scale[3]
+ scale[3] = (1+dscale) * scale[2]
+ rms[1] = rms[2]
+ rms[2] = rms[3]
+ call twk_spec (twk, shift, scale[3])
+ rms[3] = TWK_RMS(twk)
+ lastscale = TWK_SCALE(twk)
+ } else {
+ dscale = dscale / 2
+ scale[1] = (1-dscale) * scale[2]
+ scale[3] = (1+dscale) * scale[2]
+ call twk_spec (twk, shift, scale[1])
+ rms[1] = TWK_RMS(twk)
+ call twk_spec (twk, shift, scale[3])
+ rms[3] = TWK_RMS(twk)
+ lastscale = TWK_SCALE(twk)
+ }
+ if (rms[1] < rms[2])
+ TWK_SCALE(twk) = scale[1]
+ else if (rms[3] < rms[2])
+ TWK_SCALE(twk) = scale[3]
+ else
+ TWK_SCALE(twk) = scale[2]
+ }
+
+ if (TWK_SCALE(twk) != lastscale)
+ call twk_spec (twk, shift, TWK_SCALE(twk))
+end
+
+
+# TWK_ASHIFT -- Automatically determine shift by minimizing RMS.
+
+procedure twk_ashift (twk)
+
+pointer twk #I TWK data object
+
+int i
+real scale, oshift, dshift, lastshift, shift[3], rms[3]
+errchk twk_spec
+
+begin
+ dshift = TWK_DSHIFT(twk)
+ if (dshift == 0.)
+ return
+ oshift = TWK_SHIFT(twk)
+ scale = TWK_SCALE(twk)
+ do i = 1, 3 {
+ shift[i] = oshift + dshift * (i - 2)
+ call twk_spec (twk, shift[i], scale)
+ rms[i] = TWK_RMS(twk)
+ lastshift = TWK_SHIFT(twk)
+ }
+ while (dshift > 0.01) {
+ if (abs (oshift - shift[2]) > 2.) {
+ TWK_SHIFT(twk) = oshift
+ break
+ }
+ if (rms[1] < rms[2]) {
+ shift[3] = shift[2]
+ shift[2] = shift[1]
+ shift[1] = shift[2] - dshift
+ rms[3] = rms[2]
+ rms[2] = rms[1]
+ call twk_spec (twk, shift[1], scale)
+ rms[1] = TWK_RMS(twk)
+ lastshift = TWK_SHIFT(twk)
+ } else if (rms[3] < rms[2]) {
+ shift[1] = shift[2]
+ shift[2] = shift[3]
+ shift[3] = shift[2] + dshift
+ rms[1] = rms[2]
+ rms[2] = rms[3]
+ call twk_spec (twk, shift[3], scale)
+ rms[3] = TWK_RMS(twk)
+ lastshift = TWK_SHIFT(twk)
+ } else {
+ dshift = dshift / 2
+ shift[1] = shift[2] - dshift
+ call twk_spec (twk, shift[1], scale)
+ rms[1] = TWK_RMS(twk)
+ shift[3] = shift[2] + dshift
+ call twk_spec (twk, shift[3], scale)
+ rms[3] = TWK_RMS(twk)
+ lastshift = TWK_SHIFT(twk)
+ }
+ if (rms[1] < rms[2])
+ TWK_SHIFT(twk) = shift[1]
+ else if (rms[3] < rms[2])
+ TWK_SHIFT(twk) = shift[3]
+ else
+ TWK_SHIFT(twk) = shift[2]
+ }
+
+ if (TWK_SHIFT(twk) != lastshift)
+ call twk_spec (twk, TWK_SHIFT(twk), scale)
+end
+
+
+# TWK_SPEC -- Evaluate the calibrated spectrum with the specified shift
+# and scale. Compute the RMS within the sample regions. Apply a
+# smoothing if necessary. The output spectrum and shift and scale
+# used are returned in the TWK data structure.
+
+procedure twk_spec (twk, shift, scale)
+
+pointer twk #I TWK data object
+real shift #I Shift
+real scale #I Scale
+
+char type
+pointer sh, cal, asi, x, y, ycal, z, rg, temp
+int i, j, k, n, ncal, nstat, box, rg_inrange()
+real thresh, amratio, norm, sum1, sum2, xcal, xcal1, zval, asieval()
+double shdr_wl()
+
+begin
+ # Dereference the data structures.
+ type = TWK_TYPE(twk)
+ sh = TWK_SH(twk)
+ cal = TWK_CAL(twk)
+ asi = IM(cal)
+ x = TWK_WAVE(twk)
+ y = SY(sh)
+ ycal = SY(cal)
+ z = TWK_SPEC(twk)
+ n = SN(sh)
+ ncal = SN(cal)
+ rg = TWK_RG(twk)
+ thresh = TWK_THRESH(twk)
+ amratio = AM(sh) / AM(cal)
+
+ # Evaluate the calibrated spectrum and the statistics.
+ norm = 0.
+ sum1 = 0.
+ sum2 = 0.
+ nstat = 0
+ do i = 0, n-1 {
+ # Spectra
+ xcal = shdr_wl (cal, Memd[x+i]) + shift
+ xcal1 = max (1., min (real(ncal), xcal))
+ #Memr[ycal+i] = asieval (asi, xcal1) ** (amratio * scale)
+ #Memr[z+i] = Memr[y+i] / (Memr[ycal+i]
+ Memr[ycal+i] = asieval (asi, xcal1)
+ if (type == 'T') {
+ Memr[ycal+i] = max (thresh, Memr[ycal+i])
+ if (Memr[ycal+i] <= 0.)
+ call error (1,
+ "Calibration spectrum negative or zero (set threshold parameter)")
+ Memr[z+i] = Memr[y+i] / (Memr[ycal+i] ** (amratio * scale))
+ } else
+ Memr[z+i] = Memr[y+i] - (Memr[ycal+i] * scale)
+ norm = norm + Memr[z+i]
+ }
+
+ do i = 3, n-4 {
+ # Statistics
+ if (rg_inrange (rg, i+1) == NO)
+ next
+# if (xcal < 1 || xcal > ncal)
+# next
+# zval = Memr[z+i]
+ zval = Memr[z+i] - (Memr[z+i-3] + Memr[z+i+3]) / 2
+ sum1 = sum1 + zval
+ sum2 = sum2 + zval * zval
+ nstat = nstat + 1
+ }
+
+ # Normalize
+ if (TWK_TYPE(twk) == 'T') {
+ norm = norm / n
+ if (norm > 0.) {
+ call adivkr (Memr[z], norm, Memr[z], n)
+ sum1 = sum1 / norm
+ sum2 = sum2 / norm / norm
+ }
+ }
+
+ # RMS
+ if (nstat == 0)
+ TWK_RMS(twk) = INDEF
+ else
+ TWK_RMS(twk) = sqrt (nstat * sum2 - sum1**2) / nstat
+
+ TWK_SHIFT(twk) = shift
+ TWK_SCALE(twk) = scale
+
+ # Smooth
+ if (TWK_BOX(twk) > 1) {
+ call malloc (temp, n, TY_REAL)
+ box = TWK_BOX(twk)
+ box = min (n, box)
+ i = (1-box) / 2
+ sum1 = 0.
+ for (j=i; j<i+box; j=j+1)
+ sum1 = sum1 + Memr[z+max(0,j)]
+ for (k=0; k<n; k=k+1) {
+ Memr[temp+k] = sum1
+ sum1 = sum1 - Memr[z+max(0,i)] + Memr[z+min(n-1,j)]
+ i = i + 1
+ j = j + 1
+ }
+ call adivkr (Memr[temp], real(box), Memr[z], n)
+ call mfree (temp, TY_REAL)
+ }
+end
+
+
+# TWK_FIT -- Interactive fitting procedure.
+
+procedure twk_fit (twk)
+
+pointer twk #I TWK data object
+
+int i, j, n, newgraph, newdata, key, wcs, pix, clgcur(), gt_geti()
+int graph1, graph2
+real wx, wy, shift[3], scale[3], dy, gt_getr()
+double shdr_wl()
+pointer sp, str, cmd, z[3]
+pointer sh, gp, gt[2], gopen(), gt_init(), rg_xrangesd()
+errchk twk_spec, twk_rmsmin
+
+begin
+ sh = TWK_SH(twk)
+ n = SN(sh)
+
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+ call salloc (z[1], n, TY_REAL)
+ call salloc (z[3], n, TY_REAL)
+ z[2] = TWK_SPEC(twk)
+
+ # Initialize the graphics.
+ gp = gopen ("stdgraph", NEW_FILE+AW_DEFER, STDGRAPH)
+ gt[1] = gt_init ()
+ call sprintf (Memc[str], SZ_LINE,
+ "%s: spectrum = %s%s, calibration = %s%s")
+ call pargstr (TWK_TYPE(twk))
+ call pargstr (IMNAME(sh))
+ call pargstr (IMSEC(sh))
+ call pargstr (IMNAME(TWK_CAL(twk)))
+ call pargstr (IMSEC(TWK_CAL(twk)))
+ call gt_sets (gt[1], GTTITLE, Memc[str])
+ if (UN_LABEL(UN(sh)) != EOS) {
+ call gt_sets (gt[1], GTXLABEL, UN_LABEL(UN(sh)))
+ call gt_sets (gt[1], GTXUNITS, UN_UNITS(UN(sh)))
+ } else
+ call gt_sets (gt[1], GTXLABEL, "Pixels")
+ call gt_sets (gt[1], GTTYPE, "line")
+
+ gt[2] = gt_init ()
+ if (UN_LABEL(UN(sh)) != EOS) {
+ call gt_sets (gt[2], GTXLABEL, UN_LABEL(UN(sh)))
+ call gt_sets (gt[2], GTXUNITS, UN_UNITS(UN(sh)))
+ } else
+ call gt_sets (gt[2], GTXLABEL, "Pixels")
+ call gt_sets (gt[2], GTTYPE, "line")
+
+ # Cursor loop.
+ if (TWK_DSCALE(twk) > 0.)
+ graph1 = 'y'
+ else
+ graph1 = 'x'
+ graph2 = GCAL
+ newdata = YES
+ key = 'r'
+ repeat {
+ switch (key) {
+ case ':':
+ call twk_colon (Memc[cmd], twk, gp, gt, wcs, newdata, newgraph)
+ case '?':
+ call twk_colon ("help", twk, gp, gt, wcs, newdata, newgraph)
+ case 'a':
+ call twk_rmsmin (twk)
+ newdata = YES
+ case 'c':
+ if (graph2 == GCAL)
+ graph2 = GNONE
+ else
+ graph2 = GCAL
+ call gt_setr (gt[2], GTYMIN, INDEF)
+ call gt_setr (gt[2], GTYMAX, INDEF)
+ newgraph = YES
+ case 'd':
+ if (graph2 == GDATA)
+ graph2 = GNONE
+ else
+ graph2 = GDATA
+ call gt_setr (gt[2], GTYMIN, INDEF)
+ call gt_setr (gt[2], GTYMAX, INDEF)
+ newgraph = YES
+ case 'e':
+ switch (graph1) {
+ case 'x':
+ if (TWK_DSHIFT(twk) == 0.)
+ TWK_DSHIFT(twk) = 0.1
+ else
+ TWK_DSHIFT(twk) = 2 * TWK_DSHIFT(twk)
+ case 'y':
+ if (TWK_DSCALE(twk) == 0.)
+ TWK_DSCALE(twk) = 0.1
+ else
+ TWK_DSCALE(twk) = min (0.99, 2 * TWK_DSCALE(twk))
+ }
+ newdata = YES
+ case 'q':
+ break
+ case 'r':
+ newgraph = YES
+ case 's':
+ dy = wx
+ call printf ("s to add sample region or n for new regions:\n")
+ if (clgcur ("cursor",wx,wy,wcs,key,Memc[cmd],SZ_LINE) == EOF)
+ break
+ switch (key) {
+ case 'n':
+ call rg_free (TWK_RG(twk))
+ call sprintf (TWK_SAMPLE(twk), TWK_SLEN, "%g:%g")
+ call pargr (dy)
+ call pargr (wx)
+ TWK_RG(twk) = rg_xrangesd (TWK_SAMPLE(twk),
+ Memd[TWK_WAVE(twk)], SN(sh))
+ newdata = YES
+ case 's':
+ call rg_free (TWK_RG(twk))
+ if (TWK_SAMPLE(twk) == '*') {
+ call sprintf (TWK_SAMPLE(twk), TWK_SLEN, "%g:%g")
+ call pargr (dy)
+ call pargr (wx)
+ } else {
+ call sprintf (Memc[cmd], SZ_LINE, ",%g:%g")
+ call pargr (dy)
+ call pargr (wx)
+ call strcat (Memc[cmd], TWK_SAMPLE(twk), TWK_SLEN)
+ }
+ TWK_RG(twk) = rg_xrangesd (TWK_SAMPLE(twk),
+ Memd[TWK_WAVE(twk)], SN(sh))
+ newdata = YES
+ }
+ case 'w':
+ call gt_window (gt[wcs], gp, "cursor", newgraph)
+ if (wcs == 1) {
+ call gt_setr (gt[2], GTXMIN, gt_getr (gt[1], GTXMIN))
+ call gt_setr (gt[2], GTXMAX, gt_getr (gt[1], GTXMAX))
+ call gt_seti (gt[2], GTXFLIP, gt_geti (gt[1], GTXFLIP))
+ } else {
+ call gt_setr (gt[1], GTXMIN, gt_getr (gt[2], GTXMIN))
+ call gt_setr (gt[1], GTXMAX, gt_getr (gt[2], GTXMAX))
+ call gt_seti (gt[1], GTXFLIP, gt_geti (gt[2], GTXFLIP))
+ }
+ case 'x', 'y':
+ pix = max (1, min (n, nint (shdr_wl (sh, double (wx))))) - 1
+ j = 1
+ dy = abs (wy - Memr[z[j]+pix])
+ do i = 2, 3
+ if (abs (wy - Memr[z[i]+pix]) < dy) {
+ j = i
+ dy = abs (wy - Memr[z[j]+pix])
+ }
+ TWK_SHIFT(twk) = shift[j]
+ TWK_SCALE(twk) = scale[j]
+ if (j == 2 && graph1 == key) {
+ if (key == 'x')
+ TWK_DSHIFT(twk) = TWK_DSHIFT(twk) / 2.
+ else if (key == 'y')
+ TWK_DSCALE(twk) = TWK_DSCALE(twk) / 2.
+ }
+ if (TWK_DSHIFT(twk) == 0.)
+ graph1 = 'y'
+ else if (TWK_DSHIFT(twk) == 0.)
+ graph1 = 'x'
+ else
+ graph1 = key
+ newdata = YES
+ default:
+ call printf ("\007\n")
+ }
+
+ if (newdata == YES) {
+ if (graph1 == 'x') {
+ shift[1] = TWK_SHIFT(twk) - TWK_DSHIFT(twk)
+ shift[2] = TWK_SHIFT(twk)
+ shift[3] = TWK_SHIFT(twk) + TWK_DSHIFT(twk)
+ scale[1] = TWK_SCALE(twk)
+ scale[2] = TWK_SCALE(twk)
+ scale[3] = TWK_SCALE(twk)
+ } else if (graph1 == 'y') {
+ shift[1] = TWK_SHIFT(twk)
+ shift[2] = TWK_SHIFT(twk)
+ shift[3] = TWK_SHIFT(twk)
+ scale[1] = TWK_SCALE(twk) * (1 - TWK_DSCALE(twk))
+ scale[2] = TWK_SCALE(twk)
+ scale[3] = TWK_SCALE(twk) * (1 + TWK_DSCALE(twk))
+ }
+ iferr {
+ TWK_SPEC(twk) = z[1]
+ call twk_spec (twk, shift[1], scale[1])
+ call asubkr (Memr[z[1]], TWK_OFFSET(twk), Memr[z[1]], n)
+ TWK_SPEC(twk) = z[3]
+ call twk_spec (twk, shift[3], scale[3])
+ call aaddkr (Memr[z[3]], TWK_OFFSET(twk), Memr[z[3]], n)
+ TWK_SPEC(twk) = z[2]
+ call twk_spec (twk, shift[2], scale[2])
+ newdata = NO
+ } then {
+ TWK_SPEC(twk) = z[2]
+ call gt_free (gt[1])
+ call gt_free (gt[2])
+ call gclose (gp)
+ call sfree (sp)
+ call erract (EA_ERROR)
+ }
+
+ call sprintf (Memc[str], SZ_LINE, "scale = %5g")
+ call pargr (TWK_SCALE(twk))
+ if (graph1 == 'y') {
+ call sprintf (Memc[cmd], SZ_LINE, " +/- %6g")
+ call pargr (TWK_DSCALE(twk))
+ call strcat (Memc[cmd], Memc[str], SZ_LINE)
+ }
+ call sprintf (Memc[cmd], SZ_LINE, ", shift = %.2f")
+ call pargr (TWK_SHIFT(twk))
+ call strcat (Memc[cmd], Memc[str], SZ_LINE)
+ if (graph1 == 'x') {
+ call sprintf (Memc[cmd], SZ_LINE, " +/- %.2f")
+ call pargr (TWK_DSHIFT(twk))
+ call strcat (Memc[cmd], Memc[str], SZ_LINE)
+ }
+ call sprintf (Memc[cmd], SZ_LINE, ", offset = %3g")
+ call pargr (TWK_OFFSET(twk))
+ call strcat (Memc[cmd], Memc[str], SZ_LINE)
+ call sprintf (Memc[cmd], SZ_LINE, ", rms = %.3g")
+ call pargr (TWK_RMS(twk))
+ call strcat (Memc[cmd], Memc[str], SZ_LINE)
+ call gt_sets (gt[1], GTCOMMENTS, Memc[str])
+
+ newgraph = YES
+ }
+
+ if (newgraph == YES) {
+ call twk_graph (twk, gp, gt, graph1, graph2, Memr[SX(sh)],
+ Memr[z[1]], Memr[z[2]], Memr[z[3]], SN(sh))
+ newgraph = NO
+ }
+ } until (clgcur ("cursor", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+
+ call gt_free (gt[1])
+ call gt_free (gt[2])
+ call gclose (gp)
+ call sfree (sp)
+end
+
+
+# TWK_GRAPH -- Make the interactive graph.
+
+procedure twk_graph (twk, gp, gt, graph1, graph2, x, y1, y2, y3, npts)
+
+pointer twk #I TWK data object
+pointer gp #I GIO pointer
+pointer gt[2] #I GTOOLS pointer
+int graph1 #I Type for graph 1
+int graph2 #I Type for graph 2
+real x[npts] #I X values
+real y1[npts] #I Y values
+real y2[npts] #I Y values
+real y3[npts] #I Y values
+int npts #I Number of values
+
+real xmin, xmax, ymin, ymax, xmin1, xmax1, ymin1, ymax1
+
+begin
+ call gclear (gp)
+ call gseti (gp, G_WCS, 1)
+ if (graph2 != GNONE) {
+ call gsview (gp, 0.1, 0.9, 0.4, 0.9)
+ call gseti (gp, G_XLABELTICKS, NO)
+ call gt_seti (gt[1], GTDRAWXLABELS, NO)
+ }
+ call gt_ascale (gp, gt[1], x, y1, npts)
+ call ggwind (gp, xmin, xmax, ymin, ymax)
+ call gt_ascale (gp, gt[1], x, y2, npts)
+ call ggwind (gp, xmin1, xmax1, ymin1, ymax1)
+ xmin = min (xmin, xmin1)
+ xmax = max (xmax, xmax1)
+ ymin = min (ymin, ymin1)
+ ymax = max (ymax, ymax1)
+ call gt_ascale (gp, gt[1], x, y3, npts)
+ call ggwind (gp, xmin1, xmax1, ymin1, ymax1)
+ xmin = min (xmin, xmin1)
+ xmax = max (xmax, xmax1)
+ ymin = min (ymin, ymin1)
+ ymax = max (ymax, ymax1)
+ call gswind (gp, xmin, xmax, ymin, ymax)
+ call gt_swind (gp, gt[1])
+ call gt_labax (gp, gt[1])
+
+ call gt_plot (gp, gt[1], x, y1, npts)
+ call gt_plot (gp, gt[1], x, y2, npts)
+ call gt_plot (gp, gt[1], x, y3, npts)
+ call rg_gxmarkr (gp, TWK_SAMPLE(twk), x, npts, 1)
+
+ switch (graph2) {
+ case GCAL:
+ call gseti (gp, G_WCS, 2)
+ call gseti (gp, G_YNMAJOR, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gt_seti (gt[2], GTDRAWXLABELS, YES)
+ call gt_seti (gt[2], GTDRAWTITLE, NO)
+ call gt_ascale (gp, gt[2], x, Memr[SY(TWK_CAL(twk))], npts)
+ call gsview (gp, 0.1, 0.9, 0.1, 0.4)
+ call gswind (gp, xmin, xmax, INDEF, INDEF)
+ call gt_swind (gp, gt[2])
+ call gt_labax (gp, gt[2])
+ call gt_plot (gp, gt[2], x, Memr[SY(TWK_CAL(twk))], npts)
+ case GDATA:
+ call gseti (gp, G_WCS, 2)
+ call gseti (gp, G_YNMAJOR, 3)
+ call gseti (gp, G_XLABELTICKS, YES)
+ call gt_seti (gt[2], GTDRAWXLABELS, YES)
+ call gt_seti (gt[2], GTDRAWTITLE, NO)
+ call gt_ascale (gp, gt[2], x, Memr[SY(TWK_SH(twk))], npts)
+ call gsview (gp, 0.1, 0.9, 0.1, 0.4)
+ call gswind (gp, xmin, xmax, INDEF, INDEF)
+ call gt_swind (gp, gt[2])
+ call gt_labax (gp, gt[2])
+ call gt_plot (gp, gt[2], x, Memr[SY(TWK_SH(twk))], npts)
+ }
+end
+
+
+# List of colon commands.
+define CMDS "|help|shift|scale|dshift|dscale|offset|smooth|sample|"
+define HELP 1 # Print help
+define SHIFT 2 # Shift
+define SCALE 3 # Scale factor
+define DSHIFT 4 # Shift intervale
+define DSCALE 5 # Scale factor interval
+define OFFSET 6 # Offset
+define SMOOTH 7 # Boxcar smoothing
+define SAMPLE 8 # Sample
+
+# TWK_COLON -- Act on colon commands.
+
+procedure twk_colon (command, twk, gp, gt, wcs, newdata, newgraph)
+
+char command[ARB] #I Colon command
+pointer twk #I TWK data object
+pointer gp #I GIO
+pointer gt[2] #I GTOOLS
+int wcs #I WCS
+int newdata #O New data flag
+int newgraph #O New graph flag
+
+int ncmd, ival, gt_geti(), strdic(), nscan()
+real rval, gt_getr()
+pointer sp, cmd, rg, rg_xrangesd()
+
+begin
+ # Check for GTOOLS command.
+ if (command[1] == '/') {
+ call gt_colon (command, gp, gt[wcs], newgraph)
+ if (wcs == 1) {
+ call gt_setr (gt[2], GTXMIN, gt_getr (gt[1], GTXMIN))
+ call gt_setr (gt[2], GTXMAX, gt_getr (gt[1], GTXMAX))
+ call gt_seti (gt[2], GTXFLIP, gt_geti (gt[1], GTXFLIP))
+ } else {
+ call gt_setr (gt[1], GTXMIN, gt_getr (gt[2], GTXMIN))
+ call gt_setr (gt[1], GTXMAX, gt_getr (gt[2], GTXMAX))
+ call gt_seti (gt[1], GTXFLIP, gt_geti (gt[2], GTXFLIP))
+ }
+ return
+ }
+
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ # Scan the command string.
+ call sscan (command)
+ call gargwrd (Memc[cmd], SZ_LINE)
+ ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CMDS)
+
+ # Execute command.
+ switch (ncmd) {
+ case HELP:
+ call gpagefile (gp, TWK_HELP(twk), TWK_TYPE(twk))
+ case SHIFT:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("shift %g\n")
+ call pargr (TWK_SHIFT(twk))
+ } else {
+ TWK_SHIFT(twk) = rval
+ newdata = YES
+ }
+ case SCALE:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("scale %g\n")
+ call pargr (TWK_SCALE(twk))
+ } else {
+ TWK_SCALE(twk) = rval
+ newdata = YES
+ }
+ case DSHIFT:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("dshift %g\n")
+ call pargr (TWK_DSHIFT(twk))
+ } else {
+ TWK_DSHIFT(twk) = rval
+ newdata = YES
+ }
+ case DSCALE:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("dscale %g\n")
+ call pargr (TWK_DSCALE(twk))
+ } else {
+ if (rval < 0. || rval >= 1.)
+ call printf ("dscale must be between zero and one\007\n")
+ else {
+ TWK_DSCALE(twk) = rval
+ newdata = YES
+ }
+ }
+ case OFFSET:
+ call gargr (rval)
+ if (nscan() == 1) {
+ call printf ("offset %g\n")
+ call pargr (TWK_OFFSET(twk))
+ } else if (rval != TWK_OFFSET(twk)) {
+ TWK_OFFSET(twk) = rval
+ call gt_setr (gt[1], GTYMIN, INDEF)
+ call gt_setr (gt[1], GTYMAX, INDEF)
+ newdata = YES
+ }
+ case SMOOTH:
+ call gargi (ival)
+ if (nscan() == 1) {
+ call printf ("smooth %d\n")
+ call pargi (TWK_BOX(twk))
+ } else {
+ ival = max (1, ival)
+ if (ival != TWK_BOX(twk)) {
+ TWK_BOX(twk) = max (1, ival)
+ newdata = YES
+ }
+ }
+ case SAMPLE:
+ call gargstr (Memc[cmd], SZ_LINE)
+ if (Memc[cmd] == EOS) {
+ call printf ("sample %s\n")
+ call pargstr (TWK_SAMPLE(twk))
+ } else {
+ ifnoerr (rg = rg_xrangesd (Memc[cmd+1], Memd[TWK_WAVE(twk)],
+ SN(TWK_SH(twk)))) {
+ call rg_free (TWK_RG(twk))
+ call strcpy (Memc[cmd+1], TWK_SAMPLE(twk), TWK_SLEN)
+ TWK_RG(twk) = rg
+ newdata = YES
+ } else
+ call erract (EA_WARN)
+ }
+ default:
+ call printf ("\007\n")
+ }
+
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 10:36:22 -0400