diff options
Diffstat (limited to 'noao/onedspec/splot')
32 files changed, 5086 insertions, 0 deletions
diff --git a/noao/onedspec/splot/anshdr.x b/noao/onedspec/splot/anshdr.x new file mode 100644 index 00000000..e314454b --- /dev/null +++ b/noao/onedspec/splot/anshdr.x @@ -0,0 +1,84 @@ +include <time.h> +include <fset.h> +include <smw.h> + +# ANS_HDR -- Add answer header in answer file + +procedure ans_hdr (sh, newimage, key, fname1, fname2, fd1, fd2) + +pointer sh +int newimage +int key +char fname1[SZ_FNAME] +char fname2[SZ_FNAME] +int fd1, fd2 + +pointer sp, time +long clktime() +int key1, open() +errchk open +data key1/0/ + +begin + # Check for valid file name + if (fd1 == NULL && fname1[1] != EOS) { + fd1 = open (fname1, APPEND, TEXT_FILE) + call fseti (fd1, F_FLUSHNL, YES) + } + if (fd2 == NULL && fname2[1] != EOS) { + fd2 = open (fname2, APPEND, TEXT_FILE) + call fseti (fd2, F_FLUSHNL, YES) + } + + # Print image name. + if (newimage == YES) { + call smark (sp) + call salloc (time, SZ_DATE, TY_CHAR) + call cnvdate (clktime(0), Memc[time], SZ_DATE) + + if (fd1 != NULL) { + call fprintf (fd1, "\n%s [%s%s]: %s\n") + call pargstr (Memc[time]) + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargstr (TITLE(sh)) + } + if (fd2 != NULL) { + call fprintf (fd2, "\n%s [%s%s]: %s\n") + call pargstr (Memc[time]) + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargstr (TITLE(sh)) + } + call sfree (sp) + } + + # Print key dependent header. + if (key != key1) { + if (key != 'm') { + if (fd1 != NULL) { + call fprintf (fd1, "%10s%10s%10s%10s%10s%10s%10s\n") + call pargstr ("center") + call pargstr ("cont") + call pargstr ("flux") + call pargstr ("eqw") + call pargstr ("core") + call pargstr ("gfwhm") + call pargstr ("lfwhm") + call flush (fd1) + } + if (fd2 != NULL) { + call fprintf (fd2, "%10s%10s%10s%10s%10s%10s%10s\n") + call pargstr ("center") + call pargstr ("cont") + call pargstr ("flux") + call pargstr ("eqw") + call pargstr ("core") + call pargstr ("gfwhm") + call pargstr ("lfwhm") + call flush (fd2) + } + } + key1 = key + } +end diff --git a/noao/onedspec/splot/autoexp.x b/noao/onedspec/splot/autoexp.x new file mode 100644 index 00000000..c36f6a8b --- /dev/null +++ b/noao/onedspec/splot/autoexp.x @@ -0,0 +1,79 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> + +# AUTO_EXP -- Auto expand around the marked region + +procedure auto_exp (gp, gt, key, wx1, x, y, n) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +int key # Key +real wx1 # Cursor position +real x[n] # Pixel coordinates +real y[n] # Pixel data for Y scaling +int n # Number of pixels + +char cmd[1] +int i, wcs +real x1, x2, y1, y2, wx2, wy, dx, xmin, xmax, ymin, ymax + +int clgcur() + +begin + # Get the current window. + call ggwind (gp, x1, x2, y1, y2) + + # Compute the new window in x. + dx = x2 - x1 + switch (key) { + case 'a': # Expand + call printf ("again:\n") + i = clgcur ("cursor", wx2, wy, wcs, key, cmd, SZ_LINE) + x1 = wx1 + x2 = wx2 + case ',': # Shift left + x1 = x1 - 0.85 * dx + x2 = x2 - 0.85 * dx + case '.': # Shift right + x1 = x1 + 0.85 * dx + x2 = x2 + 0.85 * dx + case 'z': # Zoom x axis + x1 = x1 + 0.25 * dx + x2 = x2 - 0.25 * dx + } + + if (x1 == x2) { + # Autoscale. + x1 = INDEF + x2 = INDEF + ymin = INDEF + ymax = INDEF + } else { + # Determine the y limits for pixels between x1 and x2. + xmin = min (x1, x2) + xmax = max (x1, x2) + ymin = MAX_REAL + ymax = -MAX_REAL + do i = 1, n { + if (x[i] < xmin || x[i] > xmax) + next + ymin = min (y[i], ymin) + ymax = max (y[i], ymax) + } + if (ymin > ymax) { + ymin = y1 + ymax = y2 + } else if (y1 > y2) { + y1 = ymin + ymin = ymax + ymax = y1 + } + } + + call gt_setr (gt, GTXMIN, x1) + call gt_setr (gt, GTXMAX, x2) + call gt_setr (gt, GTYMIN, ymin) + call gt_setr (gt, GTYMAX, ymax) + call replot (gp, gt, x, y, n, YES) +end diff --git a/noao/onedspec/splot/avgsnr.x b/noao/onedspec/splot/avgsnr.x new file mode 100644 index 00000000..a4ad9ceb --- /dev/null +++ b/noao/onedspec/splot/avgsnr.x @@ -0,0 +1,72 @@ +# AVGSNR -- Compute average value and signal-to-noise in region + +procedure avgsnr (sh, wx1, wy1, y, n, fd1, fd2) + +pointer sh +real wx1, wy1 +real y[n] +int n +int fd1, fd2 + +char command[SZ_FNAME] +real wx2, wy2 +real avg, snr, rms +int i, i1, i2, nsum +int wc, key + +int clgcur() + +begin + # Get second position + call printf ("m again:") + call flush (STDOUT) + i = clgcur ("cursor", wx2, wy2, wc, key, command, SZ_FNAME) + + # Fix pixel indices + call fixx (sh, wx1, wx2, wy1, wy2, i1, i2) + if (i1 == i2) { + call printf ("Cannot determine SNR - move cursor") + return + } + + # Compute avg, rms, snr + nsum = i2 - i1 + 1 + avg = 0. + rms = 0. + snr = 0. + + if (nsum > 0) { + do i = i1, i2 + avg = avg + y[i] + avg = avg / nsum + } + + if (nsum > 1) { + call alimr (y[i1], nsum, wy1, wy2) + wy1 = wy2 - wy1 + if (wy1 > 0.) { + do i = i1, i2 + rms = rms + ((y[i] - avg) / wy1) ** 2 + rms = wy1 * sqrt (rms / (nsum-1)) + snr = avg / rms + } + } + + # Print out + call printf ("avg: %10.4g rms: %10.4g snr: %8.2f\n") + call pargr (avg) + call pargr (rms) + call pargr (snr) + if (fd1 != NULL) { + call fprintf (fd1, "avg: %10.4g rms: %10.4g snr: %8.2f\n") + call pargr (avg) + call pargr (rms) + call pargr (snr) + } + if (fd2 != NULL) { + call fprintf (fd2, "avg: %10.4g rms: %10.4g snr: %8.2f\n") + call pargr (avg) + call pargr (rms) + call pargr (snr) + } +end diff --git a/noao/onedspec/splot/conflam.x b/noao/onedspec/splot/conflam.x new file mode 100644 index 00000000..c322d566 --- /dev/null +++ b/noao/onedspec/splot/conflam.x @@ -0,0 +1,28 @@ +include <error.h> +include <smw.h> + +define VLIGHT 2.997925e18 + +# CONFLAM -- Convert to FLAMBDA from FNU + +procedure conflam (sh) + +pointer sh # SHDR pointer + +int i +real lambda +pointer ang, un_open() +errchk un_open, un_ctranr + +begin + ang = un_open ("angstroms") + iferr { + do i = 0, SN(sh)-1 { + call un_ctranr (UN(sh), ang, Memr[SX(sh)+i], lambda, 1) + Memr[SY(sh)+i] = Memr[SY(sh)+i] * VLIGHT / lambda**2 + } + } then + call erract (EA_WARN) + + call un_close (ang) +end diff --git a/noao/onedspec/splot/confnu.x b/noao/onedspec/splot/confnu.x new file mode 100644 index 00000000..228cea6f --- /dev/null +++ b/noao/onedspec/splot/confnu.x @@ -0,0 +1,28 @@ +include <error.h> +include <smw.h> + +define VLIGHT 2.997925e18 + +# CONFNU -- Convert to FNU from FLAMBDA + +procedure confnu (sh) + +pointer sh # SHDR pointer + +int i +real lambda +pointer ang, un_open() +errchk un_open, un_ctranr + +begin + ang = un_open ("angstroms") + iferr { + do i = 0, SN(sh)-1 { + call un_ctranr (UN(sh), ang, Memr[SX(sh)+i], lambda, 1) + Memr[SY(sh)+i] = Memr[SY(sh)+i] * lambda**2 / VLIGHT + } + } then + call erract (EA_WARN) + + call un_close (ang) +end diff --git a/noao/onedspec/splot/deblend.x b/noao/onedspec/splot/deblend.x new file mode 100644 index 00000000..d43a9d52 --- /dev/null +++ b/noao/onedspec/splot/deblend.x @@ -0,0 +1,627 @@ +include <math.h> +include <mach.h> + +# Profile types. +define GAUSS 1 # Gaussian profile +define LORENTZ 2 # Lorentzian profile +define VOIGT 3 # Voigt profile + +# Elements of fit array. +define BKG 1 # Background +define POS 2 # Position +define INT 3 # Intensity +define GAU 4 # Gaussian FWHM +define LOR 5 # Lorentzian FWHM + +# Type of constraints. +define FIXED 1 # Fixed parameter +define SINGLE 2 # Fit a single value for all lines +define INDEP 3 # Fit independent values for all lines + + +# DOFIT -- Fit line profiles. This is an interface to DOFIT1 +# which puts parameters into the required form and vice-versa. +# It also implements a constrained approach to the solution. + +procedure dofit (fit, x, y, s, npts, dx, nsub, y1, dy, + xp, yp, gp, lp, tp, np, chisq) + +int fit[5] # Fit constraints +real x[npts] # X data +real y[npts] # Y data +real s[npts] # Sigma data +int npts # Number of points +real dx # Pixel size +int nsub # Number of subpixels +real y1 # Continuum offset +real dy # Continuum slope +real xp[np] # Profile positions +real yp[np] # Profile intensities +real gp[np] # Profile Gaussian FWHM +real lp[np] # Profile Lorentzian FWHM +int tp[np] # Profile type +int np # Number of profiles +real chisq # Chi squared + +int i, j, fit1[5] +pointer sp, a, b +errchk dofit1 + +begin + call smark (sp) + call salloc (a, 8 + 5 * np, TY_REAL) + + # Convert positions and widths relative to first component. + Memr[a] = dx + Memr[a+1] = nsub + Memr[a+2] = y1 + Memr[a+3] = dy + Memr[a+4] = yp[1] + Memr[a+5] = xp[1] + Memr[a+6] = gp[1] + Memr[a+7] = lp[1] + do i = 1, np { + b = a + 5 * i + 3 + Memr[b] = yp[i] / Memr[a+4] + Memr[b+1] = xp[i] - Memr[a+5] + switch (tp[i]) { + case GAUSS: + if (Memr[a+6] == 0.) + Memr[a+6] = gp[i] + Memr[b+2] = gp[i] / Memr[a+6] + case LORENTZ: + if (Memr[a+7] == 0.) + Memr[a+7] = lp[i] + Memr[b+3] = lp[i] / Memr[a+7] + case VOIGT: + if (Memr[a+6] == 0.) + Memr[a+6] = gp[i] + Memr[b+2] = gp[i] / Memr[a+6] + if (Memr[a+7] == 0.) + Memr[a+7] = lp[i] + Memr[b+3] = lp[i] / Memr[a+7] + } + Memr[b+4] = tp[i] + } + + # Do fit. + fit1[INT] = fit[INT] + do i = 1, fit[BKG] { + fit1[BKG] = i + fit1[GAU] = min (SINGLE, fit[GAU]) + fit1[LOR] = min (SINGLE, fit[LOR]) + do j = FIXED, fit[POS] { + fit1[POS] = j + if (np > 1 || j != INDEP) + call dofit1 (fit1, x, y, s, npts, Memr[a], np, chisq) + } + if (np > 1 && (fit[GAU] == INDEP || fit[LOR] == INDEP)) { + fit1[GAU] = fit[GAU] + fit1[LOR] = fit[LOR] + call dofit1 (fit1, x, y, s, npts, Memr[a], np, chisq) + } + } + + y1 = Memr[a+2] + dy = Memr[a+3] + do i = 1, np { + b = a + 5 * i + 3 + yp[i] = Memr[b] * Memr[a+4] + xp[i] = Memr[b+1] + Memr[a+5] + switch (tp[i]) { + case GAUSS: + gp[i] = abs (Memr[b+2] * Memr[a+6]) + case LORENTZ: + lp[i] = abs (Memr[b+3] * Memr[a+7]) + case VOIGT: + gp[i] = abs (Memr[b+2] * Memr[a+6]) + lp[i] = abs (Memr[b+3] * Memr[a+7]) + } + } + + call sfree (sp) +end + + +# DOREFIT -- Refit line profiles. This assumes the input is very close +# to the final solution and minimizes the number of calls to the +# fitting routines. This is intended for efficient use in the +# in computing bootstrap error estimates. + +procedure dorefit (fit, x, y, s, npts, dx, nsub, y1, dy, + xp, yp, gp, lp, tp, np, chisq) + +int fit[5] # Fit constraints +real x[npts] # X data +real y[npts] # Y data +real s[npts] # Sigma data +int npts # Number of points +real dx # Pixel size +int nsub # Number of subpixels +real y1 # Continuum offset +real dy # Continuum slope +real xp[np] # Profile positions +real yp[np] # Profile intensities +real gp[np] # Profile Gaussian FWHM +real lp[np] # Profile Lorentzian FWHM +int tp[np] # Profile type +int np # Number of profiles +real chisq # Chi squared + +int i +pointer sp, a, b +errchk dofit1 + +begin + call smark (sp) + call salloc (a, 8 + 5 * np, TY_REAL) + + # Convert positions and widths relative to first component. + Memr[a] = dx + Memr[a+1] = nsub + Memr[a+2] = y1 + Memr[a+3] = dy + Memr[a+4] = yp[1] + Memr[a+5] = xp[1] + Memr[a+6] = gp[1] + Memr[a+7] = lp[1] + do i = 1, np { + b = a + 5 * i + 3 + Memr[b] = yp[i] / Memr[a+4] + Memr[b+1] = xp[i] - Memr[a+5] + switch (tp[i]) { + case GAUSS: + if (Memr[a+6] == 0.) + Memr[a+6] = gp[i] + Memr[b+2] = gp[i] / Memr[a+6] + case LORENTZ: + if (Memr[a+7] == 0.) + Memr[a+7] = lp[i] + Memr[b+3] = lp[i] / Memr[a+7] + case VOIGT: + if (Memr[a+6] == 0.) + Memr[a+6] = gp[i] + Memr[b+2] = gp[i] / Memr[a+6] + if (Memr[a+7] == 0.) + Memr[a+7] = lp[i] + Memr[b+3] = lp[i] / Memr[a+7] + } + Memr[b+4] = tp[i] + } + + # Do fit. + call dofit1 (fit, x, y, s, npts, Memr[a], np, chisq) + + y1 = Memr[a+2] + dy = Memr[a+3] + do i = 1, np { + b = a + 5 * i + 3 + yp[i] = Memr[b] * Memr[a+4] + xp[i] = Memr[b+1] + Memr[a+5] + switch (tp[i]) { + case GAUSS: + gp[i] = abs (Memr[b+2] * Memr[a+6]) + case LORENTZ: + lp[i] = abs (Memr[b+3] * Memr[a+7]) + case VOIGT: + gp[i] = abs (Memr[b+2] * Memr[a+6]) + lp[i] = abs (Memr[b+3] * Memr[a+7]) + } + } + + call sfree (sp) +end + + +# MODEL -- Compute model. + +real procedure model (x, dx, nsub, xp, yp, gp, lp, tp, np) + +real x # X value to be evaluated +real dx # Pixel width +int nsub # Number of subpixels +real xp[np] # Profile positions +real yp[np] # Profile intensities +real gp[np] # Profile Gaussian FWHM +real lp[np] # Profile Lorentzian FWHM +int tp[np] # Profile type +int np # Number of profiles + +int i, j +real delta, x1, y, arg1, arg2, v, v0, u + +begin + delta = dx / nsub + x1 = x - (dx + delta) / 2 + y = 0. + do j = 1, nsub { + x1 = x1 + delta + do i = 1, np { + switch (tp[i]) { + case GAUSS: + arg1 = 1.66511 * abs ((x1 - xp[i]) / gp[i]) + if (arg1 < 5.) + y = y + yp[i] * exp (-arg1**2) + case LORENTZ: + arg2 = abs ((x1 - xp[i]) / (lp[i] / 2)) + y = y + yp[i] / (1 + arg2**2) + case VOIGT: + arg1 = 1.66511 * (x1 - xp[i]) / gp[i] + arg2 = 0.832555 * lp[i] / gp[i] + call voigt (0., arg2, v0, u) + call voigt (arg1, arg2, v, u) + y = y + yp[i] * v / v0 + } + } + } + y = y / nsub + return (y) +end + + +# DERIVS -- Compute model and derivatives for MR_SOLVE procedure. +# This could be optimized more for the Voigt profile by reversing +# the do loops since v0 need only be computed once per line. + +procedure derivs (x, a, y, dyda, na) + +real x # X value to be evaluated +real a[na] # Parameters +real y # Function value +real dyda[na] # Derivatives +int na # Number of parameters + +int i, j, nsub +real dx, dx1, delta, x1, wg, wl, arg1, arg2, I0, dI, c, u, v, v0 + +begin + dx = a[1] + nsub = a[2] + delta = dx / nsub + dx1 = .1 * delta + x1 = x - (dx + delta) / 2 + y = 0. + do i = 1, na + dyda[i] = 0. + do j = 1, nsub { + x1 = x1 + delta + y = y + a[3] + a[4] * x1 + dyda[3] = dyda[3] + 1. + dyda[4] = dyda[4] + x1 + do i = 9, na, 5 { + switch (a[i+4]) { + case GAUSS: + I0 = a[5] * a[i] + wg = a[7] * a[i+2] + arg1 = 1.66511 * (x1 - a[6] - a[i+1]) / wg + if (abs (arg1) < 5.) { + dI = exp (-arg1**2) + c = I0 * dI * arg1 + y = y + I0 * dI + dyda[5] = dyda[5] + a[i] * dI + dyda[6] = dyda[6] + c / wg + dyda[7] = dyda[7] + c * arg1 / a[7] + dyda[i] = dyda[i] + a[5] * dI + dyda[i+1] = dyda[i+1] + c / wg + dyda[i+2] = dyda[i+2] + c * arg1 / a[i+2] + } + case LORENTZ: + I0 = a[5] * a[i] + wl = (a[8] * a[i+3] / 2) + arg2 = (x1 - a[6] - a[i+1]) / wl + dI = 1 / (1 + arg2**2) + c = 2 * I0 * dI * dI * arg2 + y = y + I0 * dI + dyda[5] = dyda[5] + a[i] * dI + dyda[6] = dyda[6] + c / wl + dyda[8] = dyda[8] + c * arg2 / a[8] + dyda[i] = dyda[i] + a[5] * dI + dyda[i+1] = dyda[i+1] + c / wl + dyda[i+3] = dyda[i+3] + c * arg2 / a[i+3] + case VOIGT: + a[7] = max (dx1, abs(a[7])) + a[8] = max (dx1, abs(a[8])) + a[i+2] = max (1E-6, abs(a[i+2])) + a[i+3] = max (1E-6, abs(a[i+3])) + + I0 = a[5] * a[i] + wg = a[7] * a[i+2] + wl = a[8] * a[i+3] + arg1 = 1.66511 * (x1 - a[6] - a[i+1]) / wg + arg2 = 0.832555 * wl / wg + call voigt (0., arg2, v0, u) + call voigt (arg1, arg2, v, u) + v = v / v0; u = u / v0 + dI = (1 - v) / (v0 * SQRTOFPI) + c = 2 * I0 * arg2 + y = y + I0 * v + dyda[5] = dyda[5] + a[i] * v + dyda[6] = dyda[6] + 2 * c * (arg1 * v - arg2 * u) / wl + dyda[7] = dyda[7] + + c * (dI + arg1 * (arg1 / arg2 * v - 2 * u)) / a[7] + dyda[8] = dyda[8] + c * (arg1 * u - dI) / a[8] + dyda[i] = dyda[i] + a[5] * v + dyda[i+1] = dyda[i+1] + 2 * c * (arg1 * v - arg2 * u) / wl + dyda[i+2] = dyda[i+2] + + c * (dI + arg1 * (arg1 / arg2 * v - 2 * u)) / a[i+2] + dyda[i+3] = dyda[i+3] + c * (arg1 * u - dI) / a[i+3] + } + } + } + y = y / nsub + do i = 1, na + dyda[i] = dyda[i] / nsub +end + + +# DOFIT1 -- Perform nonlinear iterative fit for the specified parameters. +# This uses the Levenberg-Marquardt method from NUMERICAL RECIPES. + +procedure dofit1 (fit, x, y, s, npts, a, nlines, chisq) + +int fit[5] # Fit constraints +real x[npts] # X data +real y[npts] # Y data +real s[npts] # Sigma data +int npts # Number of points +real a[ARB] # Fitting parameters +int nlines # Number of lines +real chisq # Chi squared + +int i, np, nfit +real mr, chi2 +pointer sp, flags, ptr +errchk mr_solve + +begin + # Number of terms is 5 for each line plus common background, center, + # intensity and widths. Also the pixel size and number of subpixels. + + np = 5 * nlines + 8 + + call smark (sp) + call salloc (flags, np, TY_INT) + ptr = flags + + # Background. + switch (fit[BKG]) { + case SINGLE: + Memi[ptr] = 3 + Memi[ptr+1] = 4 + ptr = ptr + 2 + } + + # Peaks. + switch (fit[INT]) { + case SINGLE: + Memi[ptr] = 5 + ptr = ptr + 1 + case INDEP: + do i = 1, nlines { + Memi[ptr] = 5 * i + 4 + ptr = ptr + 1 + } + } + + # Positions. + switch (fit[POS]) { + case SINGLE: + Memi[ptr] = 6 + ptr = ptr + 1 + case INDEP: + do i = 1, nlines { + Memi[ptr] = 5 * i + 5 + ptr = ptr + 1 + } + } + + # Gaussian FWHM. + switch (fit[GAU]) { + case SINGLE: + Memi[ptr] = 7 + ptr = ptr + 1 + case INDEP: + do i = 1, nlines { + Memi[ptr] = 5 * i + 6 + ptr = ptr + 1 + } + } + + # Lorentzian FWHM. + switch (fit[LOR]) { + case SINGLE: + Memi[ptr] = 8 + ptr = ptr + 1 + case INDEP: + do i = 1, nlines { + Memi[ptr] = 5 * i + 7 + ptr = ptr + 1 + } + } + + nfit = ptr - flags + mr = -1. + i = 0 + chi2 = MAX_REAL + repeat { + call mr_solve (x, y, s, npts, a, Memi[flags], np, nfit, mr, chisq) + if (chi2 - chisq > 0.0001) + i = 0 + else + i = i + 1 + chi2 = chisq + } until (i == 5) + + mr = 0. + call mr_solve (x, y, s, npts, a, Memi[flags], np, nfit, mr, chisq) + + call sfree (sp) +end + + +# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization. +# +# Use the Levenberg-Marquardt method to minimize the chi squared of a set +# of paraemters. The parameters being fit are indexed by the flag array. +# To initialize the Marquardt parameter, MR, is less than zero. After that +# the parameter is adjusted as needed. To finish set the parameter to zero +# to free memory. This procedure requires a subroutine, DERIVS, which +# takes the derivatives of the function being fit with respect to the +# parameters. There is no limitation on the number of parameters or +# data points. For a description of the method see NUMERICAL RECIPES +# by Press, Flannery, Teukolsky, and Vetterling, p523. + +procedure mr_solve (x, y, s, npts, params, flags, np, nfit, mr, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +real s[npts] # Sigma data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +int nfit # Number of parameters to fit +real mr # MR parameter +real chisq # Chi square of fit + +int i +real chisq1 +pointer new, a1, a2, delta1, delta2 + +errchk mr_invert + +begin + # Allocate memory and initialize. + if (mr < 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + + call malloc (new, np, TY_REAL) + call malloc (a1, nfit*nfit, TY_REAL) + call malloc (a2, nfit*nfit, TY_REAL) + call malloc (delta1, nfit, TY_REAL) + call malloc (delta2, nfit, TY_REAL) + + call amovr (params, Memr[new], np) + call mr_eval (x, y, s, npts, Memr[new], flags, np, Memr[a2], + Memr[delta2], nfit, chisq) + mr = 0.001 + } + + # Restore last good fit and apply the Marquardt parameter. + call amovr (Memr[a2], Memr[a1], nfit * nfit) + call amovr (Memr[delta2], Memr[delta1], nfit) + do i = 1, nfit + Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr) + + # Matrix solution. + call mr_invert (Memr[a1], Memr[delta1], nfit) + + # Compute the new values and curvature matrix. + do i = 1, nfit + Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1] + call mr_eval (x, y, s, npts, Memr[new], flags, np, Memr[a1], + Memr[delta1], nfit, chisq1) + + # Check if chisq has improved. + if (chisq1 < chisq) { + mr = max (EPSILONR, 0.1 * mr) + chisq = chisq1 + call amovr (Memr[a1], Memr[a2], nfit * nfit) + call amovr (Memr[delta1], Memr[delta2], nfit) + call amovr (Memr[new], params, np) + } else + mr = 10. * mr + + if (mr == 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + } +end + + +# MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS. + +procedure mr_eval (x, y, s, npts, params, flags, np, a, delta, nfit, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +real s[npts] # Sigma data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +real a[nfit,nfit] # Curvature matrix +real delta[nfit] # Delta array +int nfit # Number of parameters to fit +real chisq # Chi square of fit + +int i, j, k +real ymod, dy, dydpj, dydpk, sig2i +pointer sp, dydp + +begin + call smark (sp) + call salloc (dydp, np, TY_REAL) + + do j = 1, nfit { + do k = 1, j + a[j,k] = 0. + delta[j] = 0. + } + + chisq = 0. + do i = 1, npts { + call derivs (x[i], params, ymod, Memr[dydp], np) + if (IS_INDEF(ymod)) + next + sig2i = 1. / (s[i] * s[i]) + dy = y[i] - ymod + do j = 1, nfit { + dydpj = Memr[dydp+flags[j]-1] * sig2i + delta[j] = delta[j] + dy * dydpj + do k = 1, j { + dydpk = Memr[dydp+flags[k]-1] + a[j,k] = a[j,k] + dydpj * dydpk + } + } + chisq = chisq + dy * dy * sig2i + } + + do j = 2, nfit + do k = 1, j-1 + a[k,j] = a[j,k] + + call sfree (sp) +end + + +# MR_INVERT -- Solve a set of linear equations using Householder transforms. + +procedure mr_invert (a, b, n) + +real a[n,n] # Input matrix and returned inverse +real b[n] # Input RHS vector and returned solution +int n # Dimension of input matrices + +int krank +real rnorm +pointer sp, h, g, ip + +begin + call smark (sp) + call salloc (h, n, TY_REAL) + call salloc (g, n, TY_REAL) + call salloc (ip, n, TY_INT) + + call hfti (a, n, n, n, b, n, 1, 1E-10, krank, rnorm, + Memr[h], Memr[g], Memi[ip]) + + call sfree (sp) +end diff --git a/noao/onedspec/splot/eqwidth.x b/noao/onedspec/splot/eqwidth.x new file mode 100644 index 00000000..0594041a --- /dev/null +++ b/noao/onedspec/splot/eqwidth.x @@ -0,0 +1,109 @@ +# EQWIDTH -- Compute equivalent width, flux and center + +procedure eqwidth (sh, gfd, wx1, wy1, x, y, n, fd1, fd2) + +pointer sh +int gfd +real wx1, wy1 +real x[ARB] +real y[ARB] +int n +int fd1, fd2 + +char command[SZ_FNAME] +real wx2, wy2, sigma0, invgain +real flux_diff, rsum[2], esum[2], sum[2], cont, ctr[2] +int i, wc, key +pointer sp, s + +real clgetr() +int clgcur() +double shdr_wl() + +begin + # Get second position + call printf ("e again:") + i = clgcur ("cursor", wx2, wy2, wc, key, command, SZ_FNAME) + + if (wx1 == wx2) { + call printf ("Cannot get EQW - move cursor") + return + } + + # Set noise. + sigma0 = clgetr ("sigma0") + invgain = clgetr ("invgain") + if (IS_INDEF(sigma0) || IS_INDEF(invgain) || sigma0<0. || invgain<0.) { + sigma0 = INDEF + invgain = INDEF + } + call smark (sp) + call salloc (s, n, TY_REAL) + if (IS_INDEF(invgain)) + call amovkr (INDEF, Memr[s], n) + else { + do i = 1, n { + if (y[i] > 0) + Memr[s+i-1] = sqrt (sigma0 ** 2 + invgain * y[i]) + else + Memr[s+i-1] = sqrt (sigma0 ** 2) + } + } + + # Derive the needed values + call sumflux (sh, x, y, Memr[s], n, wx1, wx2, wy1, wy2, + sum, rsum, esum, ctr) + + # Compute difference in flux between ramp and spectrum + flux_diff = sum[1] - rsum[1] + + # Compute eq. width of feature using ramp midpoint as + # continuum + cont = 0.5 * (wy1 + wy2) + + # Print on status line - save in answer buffer + call printf ( + "center = %9.7g, eqw = %9.4f, continuum = %9.7g flux = %9.6g\n") + call pargr (ctr[1]) + call pargr (esum[1]) + call pargr (cont) + call pargr (flux_diff) + + if (fd1 != NULL) { + call fprintf (fd1, " %9.7g %9.7g %9.6g %9.4g\n") + call pargr (ctr[1]) + call pargr (cont) + call pargr (flux_diff) + call pargr (esum[1]) + } + if (fd2 != NULL) { + call fprintf (fd2, " %9.7g %9.7g %9.6g %9.4g\n") + call pargr (ctr[1]) + call pargr (cont) + call pargr (flux_diff) + call pargr (esum[1]) + } + if (!IS_INDEF(sigma0)) { + if (fd1 != NULL) { + call fprintf (fd1, + " (%7.5g) %9w (%7.4g) (%7.2g)\n") + call pargr (ctr[2]) + call pargr (sum[2]) + call pargr (esum[2]) + } + if (fd2 != NULL) { + call fprintf (fd2, + " (%7.5g) %9w (%7.4g) (%7.2g)\n") + call pargr (ctr[2]) + call pargr (sum[2]) + call pargr (esum[2]) + } + } + + # Draw cursor position + i = max (1, min (n, nint (shdr_wl (sh, double(ctr[1]))))) + call gline (gfd, wx1, wy1, wx2, wy2) + call gline (gfd, ctr[1], cont, ctr[1], y[i]) + + call sfree (sp) +end diff --git a/noao/onedspec/splot/eqwidthcp.x b/noao/onedspec/splot/eqwidthcp.x new file mode 100644 index 00000000..4fc4fd3d --- /dev/null +++ b/noao/onedspec/splot/eqwidthcp.x @@ -0,0 +1,240 @@ +include <gset.h> + + +# EQWIDTH_CP -- Equivalent width following algorithm provided by +# Caty Pilachowski. This assumes a Gaussian line profile +# and fits to the specified core level, the width at the +# specified flux level above the core, and the specified +# continuum. The line position is found by searching +# near the vertical cursor for the nearest minimum. + +define LEFT 1 # Fit to left edge +define RIGHT 2 # Fit to right edge +define BOTH 3 # Fit to both edges + +procedure eqwidth_cp (sh, gfd, center, cont, ylevel, y, n, key, fd1, fd2, + xg, yg, sg, lg, pg, ng) + +pointer sh +int gfd +real center, cont, ylevel +real y[n] +int n +int key +int fd1, fd2 +pointer xg, yg, sg, lg, pg # Pointers to fit parameters +int ng # Number of components + +int i, i1, i2, isrch, icore, edge +double xleft, xright, rcore, rinter, yl, gfwhm, lfwhm, flux, eqw, w, w1, w2 +double xpara[3], ypara[3], coefs[3], xcore, ycore +double shdr_lw(), shdr_wl() + +# Initialize reasonable values +# isrch -- nr of pixels on either side of cursor to search for min + +data isrch /3/ + +begin + # Check continuum. + if (cont <= 0.) { + call eprintf ("Continuum cannot be less than zero.\n") + return + } + + # Determine which edges of the line to use. + switch (key) { + case 'a', 'l': + edge = LEFT + case 'b', 'r': + edge = RIGHT + default: + edge = BOTH + } + + # Search for local minimum or maximum + icore = max (1, min (n, nint (shdr_wl (sh, double(center))))) + i1 = max (1, icore-isrch) + i2 = min (n, icore+isrch) + + # If half lines is selected, restrict the search + if (edge == LEFT) + i2 = max (i2-2, icore+1) + if (edge == RIGHT) + i1 = min (i1+2, icore-1) + + # Search for core. + # Someday it may be desirable to use parabolic interpolation + # to locate an estimated minimum or maximum for the region + do i = i1, i2 { + if (abs (y[i] - cont) > abs (y[icore] - cont)) + icore = i + } + + # Fit parabola to three points around minimum pixel + xpara[1] = icore - 1 + xpara[2] = icore + xpara[3] = icore + 1 + ypara[1] = y[icore-1] + ypara[2] = y[icore] + ypara[3] = y[icore+1] + + call para (xpara, ypara, coefs) + + # Compute pixel value at minimum + xcore = -coefs[2] / 2.0 / coefs[3] + ycore = coefs[1] + coefs[2] * xcore + coefs[3] * xcore**2 + + # Locate left and right line edges. If the ylevel is INDEF then use + # the half flux point. + if (IS_INDEF (ylevel)) + yl = (cont + ycore) / 2. + else + yl = ylevel + + rcore = abs (ycore - cont) + rinter = abs (yl - cont) + + if (rcore <= rinter) { + call eprintf ( + "Y cursor must be between the continuum and the line core\n") + return + } + + # Bound flux level of interest + if ((edge == LEFT) || (edge == BOTH)) { + for (i=icore; i >= 1; i=i-1) + if (abs (y[i] - cont) < rinter) + break + + if (i < 1) { + call eprintf ("Can't find left edge of line\n") + return + } + + xleft = float (i) + (yl - y[i]) / (y[i+1] - y[i]) + if (edge == LEFT) + xright = xcore + (xcore - xleft) + } + + # Now bound the right side + if ((edge == RIGHT) || (edge == BOTH)) { + for (i=icore; i <= n; i=i+1) + if (abs (y[i] - cont) < rinter) + break + + if (i > n) { + call eprintf ("Can't find right edge of line\n") + return + } + + xright = float (i) - (yl - y[i]) / (y[i-1] - y[i]) + if (edge == RIGHT) + xleft = xcore - (xright - xcore) + } + + # Compute in wavelength + w = shdr_lw (sh, double(xcore)) + w1 = shdr_lw (sh, double(xleft)) + w2 = shdr_lw (sh, double(xright)) + + # Apply Gaussian model + gfwhm = 1.665109 * abs (w2 - w1) / 2. / sqrt (log (rcore/rinter)) + lfwhm = 0. + rcore = ycore - cont + flux = 1.064467 * rcore * gfwhm + eqw = -flux / cont + + call printf ( + "center = %9.7g, eqw = %9.4g, gfwhm = %9.4g\n") + call pargd (w) + call pargd (eqw) + call pargd (gfwhm) + + if (fd1 != NULL) { + call fprintf (fd1, " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargd (w) + call pargr (cont) + call pargd (flux) + call pargd (eqw) + call pargd (ycore - cont) + call pargd (gfwhm) + call pargd (lfwhm) + } + if (fd2 != NULL) { + call fprintf (fd2, " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargd (w) + call pargr (cont) + call pargd (flux) + call pargd (eqw) + call pargd (ycore - cont) + call pargd (gfwhm) + call pargd (lfwhm) + } + + # Mark line computed + call gline (gfd, real(w), cont, real(w), real(ycore)) + call gline (gfd, real(w1), real(yl), real(w2), real(yl)) + + w1 = w - 2 * gfwhm + w2 = cont + rcore * exp (-(1.665109*(w1-w)/gfwhm)**2) + call gseti (gfd, G_PLTYPE, 2) + call gseti (gfd, G_PLCOLOR, 2) + call gamove (gfd, real(w1), real(w2)) + for (; w1 <= w+2*gfwhm; w1=w1+0.05*gfwhm) { + w2 = cont + rcore * exp (-(1.665109*(w1-w)/gfwhm)**2) + call gadraw (gfd, real(w1), real(w2)) + } + call gseti (gfd, G_PLTYPE, 1) + call gseti (gfd, G_PLCOLOR, 1) + + # Save fit parameters + if (ng == 0) { + call malloc (xg, 1, TY_REAL) + call malloc (yg, 1, TY_REAL) + call malloc (sg, 1, TY_REAL) + call malloc (lg, 1, TY_REAL) + call malloc (pg, 1, TY_INT) + } else if (ng != 1) { + call realloc (xg, 1, TY_REAL) + call realloc (yg, 1, TY_REAL) + call realloc (sg, 1, TY_REAL) + call realloc (lg, 1, TY_REAL) + call realloc (pg, 1, TY_INT) + } + Memr[xg] = w + Memr[yg] = rcore + Memr[sg] = gfwhm + Memr[lg] = lfwhm + Memi[pg] = 1 + ng = 1 +end + +# PARA -- Fit a parabola to three points + +procedure para (x, y, c) + +double x[3], y[3], c[3] +double x12, x13, x23, x213, x223, y13, y23 + +begin + x12 = x[1] - x[2] + x13 = x[1] - x[3] + x23 = x[2] - x[3] + + if (x12 == 0. || x13 == 0. || x23 == 0.) + call error (1, "X points are not distinct") + + # Compute relative to an origin at x[3] + x213 = x13 * x13 + x223 = x23 * x23 + y13 = y[1] - y[3] + y23 = y[2] - y[3] + c[3] = (y13 - y23 * x13 / x23) / (x213 - x223 * x13 / x23) + c[2] = (y23 - c[3] * x223) / x23 + c[1] = y[3] + + # Compute relative to an origin at 0. + c[1] = c[1] - x[3] * (c[2] - c[3] * x[3]) + c[2] = c[2] - 2 * c[3] * x[3] +end diff --git a/noao/onedspec/splot/fixx.x b/noao/onedspec/splot/fixx.x new file mode 100644 index 00000000..65bd4e38 --- /dev/null +++ b/noao/onedspec/splot/fixx.x @@ -0,0 +1,27 @@ +include <smw.h> + +# FIXX - Adjust so that pixel indices are increasing. + +procedure fixx (sh, x1, x2, y1, y2, i1, i2) + +pointer sh +real x1, x2, y1, y2 +int i1, i2 + +double z, z1, z2, shdr_wl(), shdr_lw() + +begin + z1 = x1 + z2 = x2 + z1 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z1))) + z2 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z2))) + if (z1 > z2) { + z = y1; y1 = y2; y2 = z + z = z1; z1 = z2; z2 = z + } + + x1 = shdr_lw (sh, z1) + x2 = shdr_lw (sh, z2) + i1 = nint (z1) + i2 = nint (z2) +end diff --git a/noao/onedspec/splot/flatten.x b/noao/onedspec/splot/flatten.x new file mode 100644 index 00000000..aa038d27 --- /dev/null +++ b/noao/onedspec/splot/flatten.x @@ -0,0 +1,110 @@ +include <pkg/gtools.h> + +# FLATTEN -- Flatten a spectrum and normalize to 1.0 +# Use ICFIT for fitting the spectrum + +procedure flatten (gp, gt, x, y, n) + +pointer gp, gt +real x[n] +real y[n] +int n + +bool b +real wx, z +int i, key +pointer sp, str, w, gt2, ic, cv + +bool clgetb() +real clgetr(), ic_getr(), cveval() +int clgeti(), ic_geti(), btoi(), clgcur() +errchk icg_fit + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (w, n, TY_REAL) + + key = '?' + repeat { + switch (key) { + case '/', '-', 'f', 'c', 'n': + call ic_open (ic) + call clgstr ("function", Memc[str], SZ_FNAME) + call ic_pstr (ic, "function", Memc[str]) + call ic_puti (ic, "order", clgeti ("order")) + call ic_putr (ic, "low", clgetr ("low_reject")) + call ic_putr (ic, "high", clgetr ("high_reject")) + call ic_puti (ic, "niterate", clgeti ("niterate")) + call ic_putr (ic, "grow", clgetr ("grow")) + call ic_puti (ic, "markrej", btoi (clgetb ("markrej"))) + switch (key) { + case '/': + call ic_puti (ic, "key", 4) + case '-': + call ic_puti (ic, "key", 3) + case 'f', 'n', 'c': + call ic_puti (ic, "key", 1) + } + + call ic_putr (ic, "xmin", min (x[1], x[n])) + call ic_putr (ic, "xmax", max (x[1], x[n])) + + call gt_copy (gt, gt2) + call gt_gets (gt2, GTXLABEL, Memc[str], SZ_FNAME) + call ic_pstr (ic, "xlabel", Memc[str]) + call gt_gets (gt2, GTYLABEL, Memc[str], SZ_FNAME) + call ic_pstr (ic, "ylabel", Memc[str]) + call gt_gets (gt2, GTXUNITS, Memc[str], SZ_FNAME) + call ic_pstr (ic, "xunits", Memc[str]) + call gt_gets (gt2, GTYUNITS, Memc[str], SZ_FNAME) + call ic_pstr (ic, "yunits", Memc[str]) + + call amovkr (1., Memr[w], n) + call icg_fit (ic, gp, "cursor", gt2, cv, x, y, Memr[w], n) + + switch (key) { + case '/': + do i = 1, n { + z = cveval (cv, x[i]) + if (abs (z) < 1e-30) + y[i] = 1. + else + y[i] = y[i] / z + } + case '-': + do i = 1, n + y[i] = y[i] - cveval (cv, x[i]) + case 'f': + do i = 1, n + y[i] = cveval (cv, x[i]) + case 'c': + call ic_clean (ic, cv, x, y, Memr[w], n) + case 'n': + ; + } + + call ic_gstr (ic, "function", Memc[str], SZ_FNAME) + call clpstr ("function", Memc[str]) + call clputi ("order", ic_geti (ic, "order")) + call clputr ("low_reject", ic_getr (ic, "low")) + call clputr ("high_reject", ic_getr (ic, "high")) + call clputi ("niterate", ic_geti (ic, "niterate")) + call clputr ("grow", ic_getr (ic, "grow")) + b = (ic_geti (ic, "markrej") == YES) + call clputb ("markrej", b) + + call cv_free (cv) + call gt_free (gt2) + call ic_closer (ic) + break + case 'q': + break + default: + call printf ( + "/=normalize, -=subtract, f=fit, c=clean, n=nop, q=quit") + } + } until (clgcur ("cursor", wx, z, i, key, Memc[str], SZ_FNAME) == EOF) + + call sfree (sp) +end diff --git a/noao/onedspec/splot/fudgept.x b/noao/onedspec/splot/fudgept.x new file mode 100644 index 00000000..c2aa3740 --- /dev/null +++ b/noao/onedspec/splot/fudgept.x @@ -0,0 +1,38 @@ +# FUDGEPT -- Fudge a point + +procedure fudgept (sh, gfd, x, y, n, wx, wy) + +pointer sh +int gfd +real x[n] +real y[n] +int n +real wx, wy + +int i1, nplot, istart +double shdr_wl() + +begin + # Get pixel number + i1 = max (1, min (n, nint (shdr_wl (sh, double(wx))))) + + # Replace with Y-value + if (i1 > 0 && i1 <= n) + y[i1] = wy + else + return + + # Plot region around new point + if (i1 > 1 && i1 < n) { + nplot = 3 + istart = i1 - 1 + } else if (i1 == 1) { + nplot = 2 + istart = i1 + } else if (i1 == n) { + nplot = 2 + istart = n - 1 + } + + call gpline (gfd, x[istart], y[istart], nplot) +end diff --git a/noao/onedspec/splot/fudgex.x b/noao/onedspec/splot/fudgex.x new file mode 100644 index 00000000..f1612b31 --- /dev/null +++ b/noao/onedspec/splot/fudgex.x @@ -0,0 +1,46 @@ +# FUDGEX -- Fudge an extended region marked by the cursor + +procedure fudgex (sh, gfd, x, y, n, wx1, wy1, xydraw) + +pointer sh +int gfd +real x[n] +real y[n] +int n +real wx1, wy1 +int xydraw + +char command[SZ_FNAME] +int i, i1, i2, wc, key +real slope +real wx2, wy2 + +int clgcur() +bool fp_equalr() + +begin + # Get second point + call printf ("x again:") + call flush (STDOUT) + i = clgcur ("cursor", wx2, wy2, wc, key, command, SZ_FNAME) + + # Fix order + call fixx (sh, wx1, wx2, wy1, wy2, i1, i2) + + if (xydraw == NO) { + wy1 = y[i1] + wy2 = y[i2] + } + if (fp_equalr (wx1, wx2)) + slope = 0. + else + slope = (wy2-wy1) / (wx2-wx1) + + # Replace pixels + do i = i1, i2 + y[i] = wy1 + (x[i] - wx1) * slope + + # Plot replaced pixels + i = i2 - i1 + 1 + call gpline (gfd, x[i1], y[i1], i) +end diff --git a/noao/onedspec/splot/getimage.x b/noao/onedspec/splot/getimage.x new file mode 100644 index 00000000..671f81de --- /dev/null +++ b/noao/onedspec/splot/getimage.x @@ -0,0 +1,159 @@ +include <error.h> +include <imhdr.h> +include <pkg/gtools.h> +include <smw.h> + +# GETIMAGE -- Read new image pixels. + +procedure getimage (image, nline, nband, nap, wave_scl, w0, wpc, units, + im, mw, sh, gt) + +char image[ARB] +int nline, nband, nap +bool wave_scl +double w0, wpc +real a, b +char units[ARB] +pointer sp, imsect, im, mw, sh, gt + +int da, n, sec[3,3], clgeti() +real gt_getr() +double shdr_lw() +pointer immap(), smw_openim() +errchk immap, shdr_open, shdr_system, un_changer + +begin + call smark (sp) + call salloc (imsect, SZ_FNAME, TY_CHAR) + + # Map the image if necessary. Don't allow image sections but + # determine requested spectrum from any explicit specification. + + da = 0 + if (im == NULL) { + call imgsection (image, Memc[imsect], SZ_FNAME) + call imgimage (image, image, SZ_FNAME) + im = immap (image, READ_ONLY, 0) + mw = smw_openim (im) + n = IM_NDIM(im) + if (Memc[imsect] != EOS) { + call amovki (1, sec[1,1], n) + call amovi (IM_LEN(im,1), sec[1,2], n) + call amovki (1, sec[1,3], n) + call id_section (Memc[imsect], sec[1,1], sec[1,2], sec[1,3], n) + switch (SMW_FORMAT(mw)) { + case SMW_ND: + if (n == 1) + da = 1 + if (n == 2) { + if (abs (sec[1,2]-sec[1,1]) == 0) { + nline = sec[1,1] + da = 2 + } else if (abs (sec[2,2]-sec[2,1]) == 0) { + nline = sec[2,1] + da = 1 + } + } else { + if (abs (sec[1,2]-sec[1,1]) == 0) { + nline = sec[1,1] + if (abs (sec[2,2]-sec[2,1]) == 0) { + nband = sec[2,1] + if (abs (sec[3,2]-sec[3,1]) > 0) + da = 3 + } else if (abs (sec[3,2]-sec[3,1]) == 0) { + nband = sec[3,1] + da = 2 + } + } else if (abs (sec[2,2]-sec[2,1]) == 0) { + nline = sec[2,1] + if (abs (sec[3,2]-sec[3,1]) == 0) { + nband = sec[3,1] + da = 1 + } + } + } + if (da > 0) { + call smw_daxis (mw, im, da, INDEFI, INDEFI) + call smw_saxis (mw, NULL, im) + } + default: + da = 1 + if (n > 1 && abs (sec[2,2]-sec[2,1]) == 0) + nline = sec[2,1] + if (n > 2 && abs (sec[3,2]-sec[3,1]) == 0) + nband = sec[3,1] + } + } + } + + # Get header info. + switch (SMW_FORMAT(mw)) { + case SMW_ND: + nap = INDEFI + n = SMW_LLEN(mw,2) + if (n > 1) { + if (nline == 0) + nline = max (1, min (n, clgeti ("line"))) + } else + nline = 0 + n = SMW_LLEN(mw,3) + if (n > 1) { + if (nband == 0) + nband = max (1, min (n, clgeti ("band"))) + } else + nband = 0 + default: + n = SMW_NSPEC(mw) + if (n > 1) { + if (nline == 0) { + nline = clgeti ("line") + nap = nline + } + } else { + nline = 0 + nap = INDEFI + } + n = SMW_NBANDS(mw) + if (n > 1) { + if (nband == 0) + nband = max (1, min (n, clgeti ("band"))) + } else + nband = 0 + } + + call shdr_open (im, mw, nline, nband, nap, SHDATA, sh) + nap = AP(sh) + nline = LINDEX(sh,1) + + if (DC(sh) == DCNO && !IS_INDEFD(w0)) + call usercoord (sh, 'l', 1D0, w0, 2D0, w0+wpc) + + # Cancel wavelength coordinates if not desired or set units. + if (!wave_scl) + call shdr_system (sh, "physical") + else { + iferr (call shdr_units (sh, units)) + ; + } + + if (da > 0) { + a = gt_getr (gt, GTXMIN) + b = gt_getr (gt, GTXMAX) + if (IS_INDEF(a) && IS_INDEF(b)) { + if (!wave_scl) { + call gt_setr (gt, GTXMIN, real(sec[da,1])) + call gt_setr (gt, GTXMAX, real(sec[da,2])) + } else { + a = shdr_lw (sh, double(sec[da,1])) + b = shdr_lw (sh, double(sec[da,2])) + call gt_setr (gt, GTXMIN, a) + call gt_setr (gt, GTXMAX, b) + } + } + } + + # Make a title. + call mktitle (sh, gt) + + call sfree (sp) +end diff --git a/noao/onedspec/splot/gfit.x b/noao/onedspec/splot/gfit.x new file mode 100644 index 00000000..2e60d8c4 --- /dev/null +++ b/noao/onedspec/splot/gfit.x @@ -0,0 +1,391 @@ +include <error.h> +include <mach.h> +include <gset.h> + +define NSUB 3 # Number of pixel subsamples +define MC_N 50 # Monte-Carlo samples +define MC_P 10 # Percent done interval (percent) +define MC_SIG 68 # Sigma sample point (percent) + +# GFIT -- Fit Gaussian + +procedure gfit (sh, gfd, wx1, wy1, wcs, pix, n, fd1, fd2, xg, yg, sg, lg, pg,ng) + +pointer sh # SHDR pointer +pointer gfd # GIO file descriptor +real wx1, wy1 # Cursor position +real wcs[n] # Spectrum data +real pix[n] # Spectrum data +int n # Number of points +int fd1, fd2 # Output file descriptors +pointer xg, yg, sg, lg, pg # Pointers to fit parameters +int ng # Number of components + +int fit[5], nsub, mc_p, mc_sig, mc_n +int i, j, i1, npts, nlines, wc, key +long seed +real w, dw, wyc, wx, wy, wx2, wy2, v, u +real slope, peak, flux, cont, gfwhm, lfwhm, eqw, scale, sscale, chisq +real sigma0, invgain, wyc1, slope1, flux1, cont1, eqw1 +bool fitit +pointer xg1, yg1, sg1, lg1 +pointer sp, cmd, x, y, s, z, ym, conte, xge, yge, sge, lge, fluxe, eqwe + +int clgeti(), clgcur() +real clgetr(), model(), gasdev(), asumr() +errchk dofit, dorefit + +define done_ 99 + +begin + call smark (sp) + call salloc (cmd, SZ_FNAME, TY_CHAR) + + # Input cursor is first continuum point now get second continuum point. + call printf ("k again:") + if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) { + call sfree (sp) + return + } + + # Set pixel indices and determine number of points to fit. + call fixx (sh, wx1, wx2, wy1, wy2, i1, j) + npts = j - i1 + 1 + if (npts < 3) { + call eprintf ("At least 3 points are required\n") + call sfree (sp) + return + } + + # Allocate space for the points to be fit. + call salloc (x, npts, TY_REAL) + call salloc (y, npts, TY_REAL) + call salloc (s, npts, TY_REAL) + call salloc (z, npts, TY_REAL) + + # Scale the data. + mc_n = clgeti ("nerrsample") + sigma0 = clgetr ("sigma0") + invgain = clgetr ("invgain") + if (IS_INDEF(sigma0) || IS_INDEF(invgain) || sigma0<0. || + invgain<0. || (sigma0 == 0. && invgain == 0.)) { + sigma0 = INDEF + invgain = INDEF + } + scale = 0. + do i = 1, npts { + Memr[x+i-1] = wcs[i1+i-1] + Memr[y+i-1] = pix[i1+i-1] + if (Memr[y+i-1] <= 0.) + if (!IS_INDEF(invgain) && invgain != 0.) { + sigma0 = INDEF + invgain = INDEF + call eprintf ( + "WARNING: Cannot compute errors with non-zero gain") + call eprintf ( + " and negative pixel values.\n") + } + scale = max (scale, abs (Memr[y+i-1])) + } + if (IS_INDEF(sigma0)) { + call amovkr (1., Memr[s], npts) + sscale = 1. + } else { + do i = 1, npts + Memr[s+i-1] = sqrt (sigma0 ** 2 + invgain * Memr[y+i-1]) + sscale = asumr (Memr[s], npts) / npts + } + call adivkr (Memr[y], scale, Memr[y], npts) + call adivkr (Memr[s], sscale, Memr[s], npts) + + # Allocate memory. + nlines = 1 + if (ng == 0) { + call malloc (xg, nlines, TY_REAL) + call malloc (yg, nlines, TY_REAL) + call malloc (sg, nlines, TY_REAL) + call malloc (lg, nlines, TY_REAL) + call malloc (pg, nlines, TY_INT) + } else if (ng != nlines) { + call realloc (xg, nlines, TY_REAL) + call realloc (yg, nlines, TY_REAL) + call realloc (sg, nlines, TY_REAL) + call realloc (lg, nlines, TY_REAL) + call realloc (pg, nlines, TY_INT) + } + ng = nlines + + # Do fit. + fit[1] = 1 + fit[2] = 2 + fit[3] = 2 + fit[4] = 2 + fit[5] = 2 + + # Setup initial estimates. + slope = (wy2-wy1) / (wx2-wx1) / scale + wyc = wy1 / scale - slope * wx1 + wx = 0 + do i = 0, npts-1 { + w = Memr[x+i] + wy = Memr[y+i] - wyc - slope * w + if (abs (wy) > wx) { + wx = abs (wy) + j = i + Memr[xg] = w + Memr[yg] = wy + } + } + + if (j > 0 && j < npts-1) { + w = Memr[x+j-1] + wy = min (0.99, max (0.01, abs (Memr[y+j-1] - wyc - slope*w) / wx)) + gfwhm = 2.355 * sqrt (-0.5 * (w-Memr[xg])**2 / log (wy)) + w = Memr[x+j+1] + wy = min (0.99, max (0.01, abs (Memr[y+j+1] - wyc - slope*w) / wx)) + gfwhm = (gfwhm + 2.355 * sqrt (-0.5*(w-Memr[xg])**2/log (wy))) / 2 + } else + gfwhm = 0.3 * abs (Memr[x+npts-1] - Memr[x]) + + switch (key) { + case 'l': + Memr[sg] = 0. + Memr[lg] = gfwhm + Memi[pg] = 2 + case 'v': + Memr[sg] = 0.5 * gfwhm + Memr[lg] = 0.5 * gfwhm + Memi[pg] = 3 + default: + Memr[sg] = gfwhm + Memr[lg] = 0. + Memi[pg] = 1 + } + + nsub = NSUB + dw = (wcs[n] - wcs[1]) / (n - 1) + iferr (call dofit (fit, Memr[x], Memr[y], Memr[s], npts, dw, nsub, + wyc, slope, Memr[xg], Memr[yg], Memr[sg], Memr[lg], Memi[pg], + ng, chisq)) { + fitit = false + goto done_ + } + + # Compute Monte-Carlo errors. + if (mc_n > 9 && !IS_INDEF(sigma0)) { + mc_p = nint (mc_n * MC_P / 100.) + mc_sig = nint (mc_n * MC_SIG / 100.) + + call salloc (ym, npts, TY_REAL) + call salloc (xg1, ng, TY_REAL) + call salloc (yg1, ng, TY_REAL) + call salloc (sg1, ng, TY_REAL) + call salloc (lg1, ng, TY_REAL) + call salloc (conte, mc_n*ng, TY_REAL) + call salloc (xge, mc_n*ng, TY_REAL) + call salloc (yge, mc_n*ng, TY_REAL) + call salloc (sge, mc_n*ng, TY_REAL) + call salloc (lge, mc_n*ng, TY_REAL) + call salloc (fluxe, mc_n*ng, TY_REAL) + call salloc (eqwe, mc_n*ng, TY_REAL) + do i = 1, npts { + w = Memr[x+i-1] + Memr[ym+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg], + Memr[sg], Memr[lg], Memi[pg], ng) + wyc + slope * w + } + seed = 1 + do i = 0, mc_n-1 { + if (i > 0 && mod (i, mc_p) == 0) { + call printf ("%2d ") + call pargi (100 * i / mc_n) + call flush (STDOUT) + } + do j = 1, npts + Memr[y+j-1] = Memr[ym+j-1] + + sscale / scale * Memr[s+j-1] * gasdev (seed) + wyc1 = wyc + slope1 = slope + call amovr (Memr[xg], Memr[xg1], ng) + call amovr (Memr[yg], Memr[yg1], ng) + call amovr (Memr[sg], Memr[sg1], ng) + call amovr (Memr[lg], Memr[lg1], ng) + call dorefit (fit, Memr[x], Memr[y], Memr[s], npts, + dw, nsub, wyc1, slope1, + Memr[xg1], Memr[yg1], Memr[sg1], Memr[lg1], Memi[pg], ng, + chisq) + + do j = 0, ng-1 { + cont = wyc + slope * Memr[xg+j] + cont1 = wyc1 + slope1 * Memr[xg+j] + switch (Memi[pg+j]) { + case 1: + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] + case 2: + flux = 1.570795 * Memr[yg+j] * Memr[lg+j] + flux1 = 1.570795 * Memr[yg1+j] * Memr[lg1+j] + case 3: + call voigt (0., 0.832555*Memr[lg+j]/Memr[sg+j], v, u) + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] / v + call voigt (0., 0.832555*Memr[lg1+j]/Memr[sg1+j], v, u) + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] / v + } + if (cont > 0. && cont1 > 0.) { + eqw = -flux / cont + eqw1 = -flux1 / cont1 + } else { + eqw = 0. + eqw1 = 0. + } + Memr[conte+j*mc_n+i] = abs (cont1 - cont) + Memr[xge+j*mc_n+i] = abs (Memr[xg1+j] - Memr[xg+j]) + Memr[yge+j*mc_n+i] = abs (Memr[yg1+j] - Memr[yg+j]) + Memr[sge+j*mc_n+i] = abs (Memr[sg1+j] - Memr[sg+j]) + Memr[lge+j*mc_n+i] = abs (Memr[lg1+j] - Memr[lg+j]) + Memr[fluxe+j*mc_n+i] = abs (flux1 - flux) + Memr[eqwe+j*mc_n+i] = abs (eqw1 - eqw) + } + } + do j = 0, ng-1 { + call asrtr (Memr[conte+j*mc_n], Memr[conte+j*mc_n], mc_n) + call asrtr (Memr[xge+j*mc_n], Memr[xge+j*mc_n], mc_n) + call asrtr (Memr[yge+j*mc_n], Memr[yge+j*mc_n], mc_n) + call asrtr (Memr[sge+j*mc_n], Memr[sge+j*mc_n], mc_n) + call asrtr (Memr[lge+j*mc_n], Memr[lge+j*mc_n], mc_n) + call asrtr (Memr[fluxe+j*mc_n], Memr[fluxe+j*mc_n], mc_n) + call asrtr (Memr[eqwe+j*mc_n], Memr[eqwe+j*mc_n], mc_n) + } + call amulkr (Memr[conte], scale, Memr[conte], mc_n*ng) + call amulkr (Memr[yge], scale, Memr[yge], mc_n*ng) + call amulkr (Memr[fluxe], scale, Memr[fluxe], mc_n*ng) + } + + call amulkr (Memr[yg], scale, Memr[yg], ng) + wyc = (wyc + slope * wx1) * scale + slope = slope * scale + + # Compute model spectrum with continuum and plot. + fitit = true + do i = 1, npts { + w = wcs[i1+i-1] + Memr[z+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg], + Memr[sg], Memr[lg], Memi[pg], ng) + wyc + slope * (w - wx1) + } + + call gseti (gfd, G_PLTYPE, 2) + call gseti (gfd, G_PLCOLOR, 2) + call gpline (gfd, wcs[i1], Memr[z], npts) + call gseti (gfd, G_PLTYPE, 3) + call gseti (gfd, G_PLCOLOR, 3) + call gline (gfd, wx1, wyc, wx2, wyc + slope * (wx2 - wx1)) + call gseti (gfd, G_PLTYPE, 1) + call gseti (gfd, G_PLCOLOR, 1) + call gflush (gfd) + +done_ + # Log computed values + if (fitit) { + do i = 1, nlines { + w = Memr[xg+i-1] + cont = wyc + slope * (w - wx1) + peak = Memr[yg+i-1] + gfwhm = Memr[sg+i-1] + lfwhm = Memr[lg+i-1] + switch (Memi[pg+i-1]) { + case 1: + flux = 1.064467 * peak * gfwhm + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, gfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (flux) + call pargr (eqw) + call pargr (gfwhm) + case 2: + flux = 1.570795 * peak * lfwhm + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, lfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (flux) + call pargr (eqw) + call pargr (lfwhm) + case 3: + call voigt (0., 0.832555*lfwhm/gfwhm, v, u) + flux = 1.064467 * peak * gfwhm / v + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, eqw = %6.4g, gfwhm = %6.4g, lfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (eqw) + call pargr (gfwhm) + call pargr (lfwhm) + } + if (fd1 != NULL) { + call fprintf (fd1, + " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargr (w) + call pargr (cont) + call pargr (flux) + call pargr (eqw) + call pargr (peak) + call pargr (gfwhm) + call pargr (lfwhm) + } + if (fd2 != NULL) { + call fprintf (fd2, + " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargr (w) + call pargr (cont) + call pargr (flux) + call pargr (eqw) + call pargr (peak) + call pargr (gfwhm) + call pargr (lfwhm) + } + if (mc_n > 9 && !IS_INDEF(sigma0)) { + if (fd1 != NULL) { + call fprintf (fd1, + " (%7.5g) (%7w) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n") + call pargr (Memr[xge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[yge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[sge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[lge+(i-1)*mc_n+mc_sig]) + } + if (fd2 != NULL) { + call fprintf (fd2, + " (%7.5g) (%7w) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n") + call pargr (Memr[xge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[yge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[sge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[lge+(i-1)*mc_n+mc_sig]) + } + } + } + } else { + call mfree (xg, TY_REAL) + call mfree (yg, TY_REAL) + call mfree (sg, TY_REAL) + call mfree (lg, TY_REAL) + call mfree (pg, TY_INT) + ng = 0 + } + + call sfree (sp) +end diff --git a/noao/onedspec/splot/mkpkg b/noao/onedspec/splot/mkpkg new file mode 100644 index 00000000..43deb993 --- /dev/null +++ b/noao/onedspec/splot/mkpkg @@ -0,0 +1,38 @@ +# SPLOT task. + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + anshdr.x <smw.h> <time.h> + autoexp.x <gset.h> <mach.h> <pkg/gtools.h> + avgsnr.x + conflam.x <error.h> <smw.h> + confnu.x <error.h> <smw.h> + deblend.x <mach.h> <math.h> + eqwidth.x + eqwidthcp.x <gset.h> + fixx.x <smw.h> + flatten.x <pkg/gtools.h> + fudgept.x + fudgex.x + getimage.x <error.h> <imhdr.h> <pkg/gtools.h> <smw.h> + gfit.x <error.h> <gset.h> <mach.h> + mktitle.x <pkg/gtools.h> <smw.h> <units.h> + plotstd.x <error.h> <gset.h> <smw.h> + replot.x <gset.h> + smooth.x + spdeblend.x <error.h> <gset.h> + splabel.x <gset.h> <smw.h> + splot.x <error.h> <gset.h> <imhdr.h> <pkg/gtools.h> <units.h>\ + <smw.h> + splotcolon.x <ctype.h> <error.h> <pkg/gtools.h> <smw.h> <units.h> + splotfun.x <mach.h> <smw.h> <error.h> + stshelp.x <error.h> + sumflux.x + usercoord.x <error.h> <smw.h> <units.h> + voigt.x + wrspect.x <error.h> <imhdr.h> <imio.h> <smw.h> <units.h> + ; diff --git a/noao/onedspec/splot/mktitle.x b/noao/onedspec/splot/mktitle.x new file mode 100644 index 00000000..554599bf --- /dev/null +++ b/noao/onedspec/splot/mktitle.x @@ -0,0 +1,41 @@ +include <pkg/gtools.h> +include <smw.h> +include <units.h> + +# MKTITLE -- Make a spectrum title (IIDS style) + +procedure mktitle (sh, gt) + +pointer sh, gt + +pointer sp, str + +begin + # Do nothing if the GTOOLS pointer is undefined. + if (gt == NULL) + return + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call sprintf (Memc[str], SZ_LINE, + "[%s%s]: %s %.2s ap:%d beam:%d") + call pargstr (IMNAME(sh)) + call pargstr (IMSEC(sh)) + call pargstr (TITLE(sh)) + call pargr (IT(sh)) + call pargi (AP(sh)) + call pargi (BEAM(sh)) + + # Set GTOOLS labels. + call gt_sets (gt, GTTITLE, Memc[str]) + if (UN_LABEL(UN(sh)) != EOS) { + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + } else { + call gt_sets (gt, GTXLABEL, LABEL(sh)) + call gt_sets (gt, GTXUNITS, UNITS(sh)) + } + + call sfree (sp) +end diff --git a/noao/onedspec/splot/plotstd.x b/noao/onedspec/splot/plotstd.x new file mode 100644 index 00000000..dab1554d --- /dev/null +++ b/noao/onedspec/splot/plotstd.x @@ -0,0 +1,70 @@ +include <error.h> +include <gset.h> +include <smw.h> + +define VLIGHT 2.997925e18 + +# PLOT_STD -- Plot the flux values for a standard star on current screen + +procedure plot_std (sh, gfd, fnu) + +pointer sh +int gfd +bool fnu + +pointer waves, bands, mags +int i, nwaves +real w1, w2 +real fnuzero, clgetr() +double shdr_lw() + +begin + # Get calibration data. + iferr (call getcalib (waves, bands, mags, nwaves)) { + call erract (EA_WARN) + return + } + + # Convert to fnu or flambda + fnuzero = clgetr ("fnuzero") + do i = 1, nwaves { + Memr[mags+i-1] = fnuzero * 10.0**(-0.4 * Memr[mags+i-1]) + if (!fnu) + Memr[mags+i-1] = Memr[mags+i-1] * VLIGHT / Memr[waves+i-1]**2 + } + + # Overplot boxes on current plot + w1 = shdr_lw (sh, double(1)) + w2 = shdr_lw (sh, double(SN(sh))) + + do i = 1, nwaves + if (Memr[waves+i-1] > w1 && Memr[waves+i-1] < w2) + call plbox2 (gfd, Memr[waves+i-1]-Memr[bands+i-1]/2, + Memr[mags+i-1], Memr[waves+i-1]+Memr[bands+i-1]/2, .015) + + call freecalib (waves, bands, mags) +end + +# PLBOX2 -- Plot a box of given height and width + +procedure plbox2 (gfd, x1, y1, x2, ndcy) + +int gfd +real x1, x2, y1, ndcy + +real ya1, ya2 +real wx1, wx2, wy1, wy2 + +begin + # Get current WCS range + call ggwind (gfd, wx1, wx2, wy1, wy2) + + # Adjust vertical spacing + ya1 = y1 - ndcy * (wy2 - wy1) + ya2 = y1 + ndcy * (wy2 - wy1) + + call gline (gfd, x1, ya1, x2, ya1) + call gline (gfd, x2, ya1, x2, ya2) + call gline (gfd, x2, ya2, x1, ya2) + call gline (gfd, x1, ya2, x1, ya1) +end diff --git a/noao/onedspec/splot/replot.x b/noao/onedspec/splot/replot.x new file mode 100644 index 00000000..9157846a --- /dev/null +++ b/noao/onedspec/splot/replot.x @@ -0,0 +1,27 @@ +include <gset.h> + +# REPLOT -- Replot the current array + +procedure replot (gfd, gt, x, y, npts, clear) + +pointer gfd +pointer gt +real x[ARB] +real y[ARB] +int npts +int clear + +int wc, gstati() + +begin + if (clear == YES) { + wc = gstati (gfd, G_WCS) + call gclear (gfd) + call gseti (gfd, G_WCS, wc) + call gt_ascale (gfd, gt, x, y, npts) + call gt_swind (gfd, gt) + call gt_labax (gfd, gt) + } + + call gt_plot (gfd, gt, x, y, npts) +end diff --git a/noao/onedspec/splot/smooth.x b/noao/onedspec/splot/smooth.x new file mode 100644 index 00000000..1418fc4f --- /dev/null +++ b/noao/onedspec/splot/smooth.x @@ -0,0 +1,54 @@ +# SMOOTH -- Box smooth the array + +procedure smooth (y, n) + +real y[ARB] +int n + +int i, j, boxsize, halfbox, del +int nsum +real sum +pointer sp, smy + +int clgeti() + +begin + call smark (sp) + call salloc (smy, n, TY_REAL) + + # Get boxsize + boxsize = clgeti ("boxsize") + if (mod (boxsize, 2) == 0) { + boxsize = boxsize + 1 + call eprintf ("WARNING: Using a box size of %d") + call pargi (boxsize) + } + + halfbox = boxsize/2 + + # This is not efficiently coded, but easy to code + # A running box mean would be faster + do i = 1, n { + sum = 0.0 + nsum = 0 + + if (i > halfbox && i < (n-halfbox)) + del = halfbox + else + if (i <= halfbox) + del = i/2 + else + del = (n - i + 1)/2 + + do j = i-del, i+del { + nsum = nsum + 1 + sum = sum + y[j] + } + + Memr[smy+i-1] = sum / nsum + } + + # Replace pixels back + call amovr (Memr[smy], y, n) + call sfree (sp) +end diff --git a/noao/onedspec/splot/spdeblend.x b/noao/onedspec/splot/spdeblend.x new file mode 100644 index 00000000..a07cd52d --- /dev/null +++ b/noao/onedspec/splot/spdeblend.x @@ -0,0 +1,819 @@ +include <error.h> +include <gset.h> + +define NSUB 3 # Number of pixel subsamples +define MC_N 50 # Monte-Carlo samples +define MC_P 10 # Percent done interval (percent) +define MC_SIG 68 # Sigma sample point (percent) + +# Profile types. +define PTYPES "|gaussian|lorentzian|voigt|" +define GAUSS 1 # Gaussian profile +define LORENTZ 2 # Lorentzian profile +define VOIGT 3 # Voigt profile + + +# SP_DEBLEND -- Deblend lines in a spectral region. + +procedure sp_deblend (sh, gfd, wx1, wy1, wcs, pix, n, fd1, fd2, + xg, yg, sg, lg, pg, ng) + +pointer sh # SHDR pointer +pointer gfd # GIO file descriptor +real wx1, wy1 # Cursor position +real wcs[n] # Coordinates +real pix[n] # Spectrum data +int n # Number of points +int fd1, fd2 # Output file descriptors +pointer xg, yg, sg, lg, pg # Pointers to fit parameters +int ng # Number of components + +int fit[5], nsub, mc_p, mc_sig, mc_n +int i, j, i1, npts, nlines, maxlines, wc, key, type, ifit +long seed +real w, dw, wyc, wx, wy, wx2, wy2, u, v +real slope, peak, flux, cont, gfwhm, lfwhm, eqw, scale, sscale, chisq, rms +real sigma0, invgain, wyc1, slope1, flux1, cont1, eqw1 +bool fitit, fitg, fitl +pointer xg1, yg1, sg1, lg1 +pointer sp, cmd, x, y, s, z, waves, types, gfwhms, lfwhms, peaks, ym +pointer conte, xge, yge, sge, lge, fluxe, eqwe + +int clgeti(), clgcur(), open(), fscan(), nscan(), strdic() +real clgetr(), model(), gasdev(), asumr() +double shdr_wl() +errchk dofit, dorefit + +define fitp_ 95 +define fitg_ 96 +define fitl_ 97 +define fitb_ 98 +define done_ 99 + +begin + call smark (sp) + call salloc (cmd, SZ_FNAME, TY_CHAR) + + # Input cursor is first continuum point now get second continuum point. + call printf ("d again:") + if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) { + call sfree (sp) + return + } + + # Set pixel indices and determine number of points to fit. + call fixx (sh, wx1, wx2, wy1, wy2, i1, j) + + npts = j - i1 + 1 + if (npts < 3) { + call eprintf ("At least 3 points are required\n") + call sfree (sp) + return + } + + # Allocate space for the points to be fit. + call salloc (x, npts, TY_REAL) + call salloc (y, npts, TY_REAL) + call salloc (s, npts, TY_REAL) + call salloc (z, npts, TY_REAL) + + # Scale the data. + sigma0 = clgetr ("sigma0") + invgain = clgetr ("invgain") + mc_n = clgeti ("nerrsample") + if (IS_INDEF(sigma0) || IS_INDEF(invgain) || sigma0<0. || + invgain<0. || (sigma0 == 0. && invgain == 0.)) { + sigma0 = INDEF + invgain = INDEF + } + scale = 0. + do i = 1, npts { + Memr[x+i-1] = wcs[i1+i-1] + Memr[y+i-1] = pix[i1+i-1] + if (Memr[y+i-1] <= 0.) + if (invgain != 0.) { + sigma0 = INDEF + invgain = INDEF + call eprintf ( + "WARNING: Cannot compute errors with non-zero gain") + call eprintf ( + " and negative pixel values.\n") + } + scale = max (scale, abs (Memr[y+i-1])) + } + if (IS_INDEF(sigma0)) { + call amovkr (1., Memr[s], npts) + sscale = 1. + } else { + do i = 1, npts + Memr[s+i-1] = sqrt (sigma0 ** 2 + invgain * Memr[y+i-1]) + sscale = asumr (Memr[s], npts) / npts + } + call adivkr (Memr[y], scale, Memr[y], npts) + call adivkr (Memr[s], sscale, Memr[s], npts) + + # Select the lines to be fit. If no lines return. + fitit = false + fitg = false + fitl = false + maxlines = 5 + call malloc (waves, maxlines, TY_REAL) + call malloc (peaks, maxlines, TY_REAL) + call malloc (gfwhms, maxlines, TY_REAL) + call malloc (lfwhms, maxlines, TY_REAL) + call malloc (types, maxlines, TY_INT) + nlines = 0 + call printf ( + "Lines ('f'ile, 'g'aussian, 'l'orentzian, 'v'oigt, 't'ype, 'q'uit:") + while (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME) != EOF) { + switch (key) { + case 'f': + call clgstr ("linelist", Memc[cmd], SZ_FNAME) + call printf ( + "Lines ('f'ile, 'g'aussian, 'l'orentzian, 'v'oigt, 't'ype, 'q'uit:") + iferr (j = open (Memc[cmd], READ_ONLY, TEXT_FILE)) { + call erract (EA_WARN) + next + } + while (fscan (j) != EOF) { + call gargr (wx) + if (nscan() < 1) + next + if (wx < min (wcs[1], wcs[n]) || wx > max (wcs[1], wcs[n])) + next + call gargr (peak) + call gargwrd (Memc[cmd], SZ_FNAME) + call gargr (gfwhm) + call gargr (lfwhm) + type = strdic (Memc[cmd], Memc[cmd], SZ_FNAME, PTYPES) + if (type == 0) + type = GAUSS + switch (nscan()) { + case 0: + next + case 1: + peak = INDEF + type = GAUSS + gfwhm = INDEF + lfwhm = INDEF + case 2: + type = GAUSS + gfwhm = INDEF + lfwhm = INDEF + case 3: + gfwhm = INDEF + lfwhm = INDEF + case 4: + switch (type) { + case GAUSS: + lfwhm = INDEF + case LORENTZ: + lfwhm = gfwhm + gfwhm = INDEF + case VOIGT: + lfwhm = INDEF + } + } + for (i = 0; i < nlines && wx != Memr[waves+i]; i = i + 1) + ; + if (i == nlines) { + if (nlines == maxlines) { + maxlines = maxlines + 5 + call realloc (waves, maxlines, TY_REAL) + call realloc (peaks, maxlines, TY_REAL) + call realloc (gfwhms, maxlines, TY_REAL) + call realloc (lfwhms, maxlines, TY_REAL) + call realloc (types, maxlines, TY_INT) + } + Memr[waves+i] = wx + Memr[peaks+i] = peak + Memr[gfwhms+i] = gfwhm + Memr[lfwhms+i] = lfwhm + Memi[types+i] = type + switch (type) { + case GAUSS: + fitg = true + case LORENTZ: + fitl = true + case VOIGT: + fitg = true + fitl = true + } + nlines = nlines + 1 + call gmark (gfd, wx, wy, GM_VLINE, 3., 3.) + } + } + call close (j) + next + case 'g': + type = GAUSS + peak = INDEF + gfwhm = INDEF + lfwhm = INDEF + case 'l': + type = LORENTZ + peak = INDEF + gfwhm = INDEF + lfwhm = INDEF + case 'v': + type = VOIGT + peak = INDEF + gfwhm = INDEF + lfwhm = INDEF + case 't': + type = GAUSS + wx = clgetr ("wavelength") + peak = INDEF + gfwhm = INDEF + lfwhm = INDEF + call printf ( + "Lines ('f'ile, 'g'aussian, 'l'orentzian, 'v'oigt, 't'ype, 'q'uit:") + case 'q': + call printf ("\n") + break + case 'I': + call fatal (0, "Interrupt") + default: + call printf ( + "Lines ('f'ile, 'g'aussian, 'l'orentzian, 'v'oigt, 't'ype, 'q'uit:\007") + next + } + for (i = 0; i < nlines && wx != Memr[waves+i]; i = i + 1) + ; + if (i == nlines) { + if (nlines == maxlines) { + maxlines = maxlines + 5 + call realloc (waves, maxlines, TY_REAL) + call realloc (peaks, maxlines, TY_REAL) + call realloc (gfwhms, maxlines, TY_REAL) + call realloc (lfwhms, maxlines, TY_REAL) + call realloc (types, maxlines, TY_INT) + } + Memr[waves+i] = wx + Memr[peaks+i] = peak + Memr[gfwhms+i] = gfwhm + Memr[lfwhms+i] = lfwhm + Memi[types+i] = type + switch (type) { + case GAUSS: + fitg = true + case LORENTZ: + fitl = true + case VOIGT: + fitg = true + fitl = true + } + nlines = nlines + 1 + call gmark (gfd, wx, wy, GM_VLINE, 3., 3.) + } + } + if (nlines == 0) + goto done_ + + # Allocate memory. + if (ng == 0) { + call malloc (xg, nlines, TY_REAL) + call malloc (yg, nlines, TY_REAL) + call malloc (sg, nlines, TY_REAL) + call malloc (lg, nlines, TY_REAL) + call malloc (pg, nlines, TY_INT) + } else if (ng != nlines) { + call realloc (xg, nlines, TY_REAL) + call realloc (yg, nlines, TY_REAL) + call realloc (sg, nlines, TY_REAL) + call realloc (lg, nlines, TY_REAL) + call realloc (pg, nlines, TY_INT) + } + ng = nlines + + # Do fits. + ifit = 0 + fit[3] = 3 + repeat { +fitp_ call printf ("Fit positions (fixed, single, all, quit) ") + if (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME) == EOF) + break + switch (key) { + case 'f': + fit[2] = 1 + case 's': + fit[2] = 2 + case 'a': + fit[2] = 3 + case 'q': + break + default: + goto fitp_ + } + if (fitg) { +fitg_ call printf ("Fit Gaussian widths (fixed, single, all, quit) ") + if (clgcur("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME)==EOF) + break + switch (key) { + case 'f': + fit[4] = 1 + case 's': + fit[4] = 2 + case 'a': + fit[4] = 3 + case 'q': + break + default: + goto fitg_ + } + } + if (fitl) { +fitl_ call printf ( + "Fit Lorentzian widths (fixed, single, all, quit) ") + if (clgcur("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME)==EOF) + break + switch (key) { + case 'f': + fit[5] = 1 + case 's': + fit[5] = 2 + case 'a': + fit[5] = 3 + case 'q': + break + default: + goto fitl_ + } + } +fitb_ call printf ("Fit background (no, yes, quit) ") + if (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME) == EOF) + break + switch (key) { + case 'n': + fit[1] = 1 + case 'y': + fit[1] = 2 + case 'q': + break + default: + goto fitb_ + } + call printf ("Fitting...") + call flush (STDOUT) + + # Setup initial estimates. + if (ifit == 0) { + slope = (wy2-wy1) / (wx2-wx1) / scale + wyc = wy1 / scale - slope * wx1 + eqw = abs (Memr[x+npts-1] - Memr[x]) / nlines + do i = 0, nlines-1 { + w = Memr[waves+i] + peak = Memr[peaks+i] + gfwhm = Memr[gfwhms+i] + lfwhm = Memr[lfwhms+i] + type = Memi[types+i] + j = max (1, min (n, nint (shdr_wl (sh, double(w))))) + Memr[xg+i] = w + if (IS_INDEF(peak)) + Memr[yg+i] = pix[j] / scale - wyc - slope * w + else + Memr[yg+i] = peak / scale + Memr[sg+i] = 0. + Memr[lg+i] = 0. + Memi[pg+i] = type + switch (type) { + case GAUSS: + if (IS_INDEF(gfwhm)) + Memr[sg+i] = 0.3 * eqw + else + Memr[sg+i] = gfwhm + case LORENTZ: + if (IS_INDEF(lfwhm)) + Memr[lg+i] = 0.3 * eqw + else + Memr[lg+i] = lfwhm + case VOIGT: + if (IS_INDEF(Memr[gfwhms+i])) + Memr[sg+i] = 0.1 * eqw + else + Memr[sg+i] = gfwhm + if (IS_INDEF(Memr[lfwhms+i])) + Memr[lg+i] = 0.1 * eqw + else + Memr[lg+i] = lfwhm + } + } + } else { + call adivkr (Memr[yg], scale, Memr[yg], ng) + slope = slope / scale + wyc = wyc / scale - slope * wx1 + } + + nsub = NSUB + dw = (wcs[n] - wcs[1]) / (n - 1) + iferr (call dofit (fit, Memr[x], Memr[y], Memr[s], npts, + dw, nsub, wyc, slope, + Memr[xg], Memr[yg], Memr[sg], Memr[lg], Memi[pg], ng, chisq)) { + call erract (EA_WARN) + next + } + ifit = ifit + 1 + + # Compute Monte-Carlo errors. + if (mc_n > 9 && !IS_INDEF(sigma0)) { + mc_p = nint (mc_n * MC_P / 100.) + mc_sig = nint (mc_n * MC_SIG / 100.) + + call salloc (ym, npts, TY_REAL) + call salloc (xg1, ng, TY_REAL) + call salloc (yg1, ng, TY_REAL) + call salloc (sg1, ng, TY_REAL) + call salloc (lg1, ng, TY_REAL) + call salloc (conte, mc_n*ng, TY_REAL) + call salloc (xge, mc_n*ng, TY_REAL) + call salloc (yge, mc_n*ng, TY_REAL) + call salloc (sge, mc_n*ng, TY_REAL) + call salloc (lge, mc_n*ng, TY_REAL) + call salloc (fluxe, mc_n*ng, TY_REAL) + call salloc (eqwe, mc_n*ng, TY_REAL) + do i = 1, npts { + w = Memr[x+i-1] + Memr[ym+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg], + Memr[sg], Memr[lg], Memi[pg], ng) + wyc + slope * w + } + seed = 1 + do i = 0, mc_n-1 { + if (i > 0 && mod (i, mc_p) == 0) { + call printf ("%2d ") + call pargi (100 * i / mc_n) + call flush (STDOUT) + } + do j = 1, npts + Memr[y+j-1] = Memr[ym+j-1] + + sscale / scale * Memr[s+j-1] * gasdev (seed) + wyc1 = wyc + slope1 = slope + call amovr (Memr[xg], Memr[xg1], ng) + call amovr (Memr[yg], Memr[yg1], ng) + call amovr (Memr[sg], Memr[sg1], ng) + call amovr (Memr[lg], Memr[lg1], ng) + call dorefit (fit, Memr[x], Memr[y], Memr[s], npts, + dw, nsub, wyc1, slope1, + Memr[xg1], Memr[yg1], Memr[sg1], Memr[lg1], Memi[pg], + ng, chisq) + + do j = 0, ng-1 { + cont = wyc + slope * Memr[xg+j] + cont1 = wyc1 + slope1 * Memr[xg+j] + switch (Memi[pg+j]) { + case GAUSS: + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] + case LORENTZ: + flux = 1.570795 * Memr[yg+j] * Memr[lg+j] + flux1 = 1.570795 * Memr[yg1+j] * Memr[lg1+j] + case VOIGT: + call voigt (0., 0.832555*Memr[lg+j]/Memr[sg+j], + v, u) + flux = 1.064467 * Memr[yg+j] * Memr[sg+j] / v + call voigt (0., 0.832555*Memr[lg1+j]/Memr[sg1+j], + v, u) + flux1 = 1.064467 * Memr[yg1+j] * Memr[sg1+j] / v + } + if (cont > 0. && cont1 > 0.) { + eqw = -flux / cont + eqw1 = -flux1 / cont1 + } else { + eqw = 0. + eqw1 = 0. + } + Memr[conte+j*mc_n+i] = abs (cont1 - cont) + Memr[xge+j*mc_n+i] = abs (Memr[xg1+j] - Memr[xg+j]) + Memr[yge+j*mc_n+i] = abs (Memr[yg1+j] - Memr[yg+j]) + Memr[sge+j*mc_n+i] = abs (Memr[sg1+j] - Memr[sg+j]) + Memr[lge+j*mc_n+i] = abs (Memr[lg1+j] - Memr[lg+j]) + Memr[fluxe+j*mc_n+i] = abs (flux1 - flux) + Memr[eqwe+j*mc_n+i] = abs (eqw1 - eqw) + } + } + do j = 0, ng-1 { + call asrtr (Memr[conte+j*mc_n], Memr[conte+j*mc_n], mc_n) + call asrtr (Memr[xge+j*mc_n], Memr[xge+j*mc_n], mc_n) + call asrtr (Memr[yge+j*mc_n], Memr[yge+j*mc_n], mc_n) + call asrtr (Memr[sge+j*mc_n], Memr[sge+j*mc_n], mc_n) + call asrtr (Memr[lge+j*mc_n], Memr[lge+j*mc_n], mc_n) + call asrtr (Memr[fluxe+j*mc_n], Memr[fluxe+j*mc_n], mc_n) + call asrtr (Memr[eqwe+j*mc_n], Memr[eqwe+j*mc_n], mc_n) + } + call amulkr (Memr[conte], scale, Memr[conte], mc_n*ng) + call amulkr (Memr[yge], scale, Memr[yge], mc_n*ng) + call amulkr (Memr[fluxe], scale, Memr[fluxe], mc_n*ng) + } + + call amulkr (Memr[yg], scale, Memr[yg], ng) + wyc = (wyc + slope * wx1) * scale + slope = slope * scale + + fitit = true + + # Compute model spectrum with continuum and plot. + call printf ("Overplot (total, components, both, none) ") + if (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], SZ_FNAME) == EOF) + break + + rms = 0. + do i = 1, npts { + w = Memr[x+i-1] + Memr[z+i-1] = model (w, dw, nsub, Memr[xg], Memr[yg], + Memr[sg], Memr[lg], Memi[pg], ng) + Memr[z+i-1] = Memr[z+i-1] + wyc + slope * (w - wx1) + rms = rms + (Memr[z+i-1] / scale - Memr[y+i-1]) ** 2 + } + + # Total. + if (key == 't' || key == 'b') { + call gseti (gfd, G_PLTYPE, 2) + call gseti (gfd, G_PLCOLOR, 2) + call gpline (gfd, Memr[x], Memr[z], npts) + call gseti (gfd, G_PLTYPE, 1) + call gflush (gfd) + } + + # Components. + if (key == 'c' || key == 'b') { + call gseti (gfd, G_PLTYPE, 3) + call gseti (gfd, G_PLCOLOR, 5) + do j = 0, ng-1 { + do i = 1, npts { + w = Memr[x+i-1] + Memr[z+i-1] = model (w, dw, nsub, Memr[xg+j], Memr[yg+j], + Memr[sg+j], Memr[lg+j], Memi[pg+j], 1) + Memr[z+i-1] = Memr[z+i-1] + wyc + slope * (w - wx1) + } + call gpline (gfd, Memr[x], Memr[z], npts) + } + call gseti (gfd, G_PLTYPE, 1) + call gflush (gfd) + } + + if (key != 'n') { + call gseti (gfd, G_PLTYPE, 4) + call gseti (gfd, G_PLCOLOR, 3) + call gline (gfd, wx1, wyc, wx2, wyc + slope * (wx2 - wx1)) + call gseti (gfd, G_PLTYPE, 1) + call gflush (gfd) + } + + + # Print computed values on status line. + i = 1 + key = '' + repeat { + switch (key) { + case '-': + i = i - 1 + if (i < 1) + i = nlines + case '+': + i = i + 1 + if (i > nlines) + i = 1 + case 'q': + call printf ("\n") + break + } + + if (key == 'r') { + call printf ("\nrms = %8.4g") + call pargr (scale * sqrt (chisq / npts)) + } else { + w = Memr[xg+i-1] + cont = wyc + slope * (w - wx1) + peak = Memr[yg+i-1] + gfwhm = Memr[sg+i-1] + lfwhm = Memr[lg+i-1] + switch (Memi[pg+i-1]) { + case GAUSS: + flux = 1.064467 * peak * gfwhm + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, gfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (flux) + call pargr (eqw) + call pargr (gfwhm) + case LORENTZ: + flux = 1.570795 * peak * lfwhm + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, flux = %8.4g, eqw = %6.4g, lfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (flux) + call pargr (eqw) + call pargr (lfwhm) + case VOIGT: + call voigt (0., 0.832555*lfwhm/gfwhm, v, u) + flux = 1.064467 * peak * gfwhm / v + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + call printf ( + "\n%d: center = %8.6g, eqw = %6.4g, gfwhm = %6.4g, lfwhm = %6.4g") + call pargi (i) + call pargr (w) + call pargr (eqw) + call pargr (gfwhm) + call pargr (lfwhm) + } + } + + call printf (" (+,-,r,q):") + call flush (STDOUT) + } until (clgcur ("cursor", + wx, wy, wc, key, Memc[cmd], SZ_FNAME) == EOF) + } + +done_ + call printf ("Deblending done\n") + # Log computed values + if (fitit) { + do i = 1, nlines { + w = Memr[xg+i-1] + cont = wyc + slope * (w - wx1) + peak = Memr[yg+i-1] + gfwhm = Memr[sg+i-1] + lfwhm = Memr[lg+i-1] + switch (Memi[pg+i-1]) { + case GAUSS: + flux = 1.064467 * peak * gfwhm + case LORENTZ: + flux = 1.570795 * peak * lfwhm + case VOIGT: + call voigt (0., 0.832555*lfwhm/gfwhm, v, u) + flux = 1.064467 * peak * gfwhm / v + } + + if (cont > 0.) + eqw = -flux / cont + else + eqw = INDEF + if (fd1 != NULL) { + call fprintf (fd1, + " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargr (w) + call pargr (cont) + call pargr (flux) + call pargr (eqw) + call pargr (peak) + call pargr (gfwhm) + call pargr (lfwhm) + } + if (fd2 != NULL) { + call fprintf (fd2, + " %9.7g %9.7g %9.6g %9.4g %9.6g %9.4g %9.4g\n") + call pargr (w) + call pargr (cont) + call pargr (flux) + call pargr (eqw) + call pargr (peak) + call pargr (gfwhm) + call pargr (lfwhm) + } + if (mc_n > 9 && !IS_INDEF(sigma0)) { + if (fd1 != NULL) { + call fprintf (fd1, + " (%7.5g) (%7.5g) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n") + call pargr (Memr[xge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[conte+(i-1)*mc_n+mc_sig]) + call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[yge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[sge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[lge+(i-1)*mc_n+mc_sig]) + } + if (fd2 != NULL) { + call fprintf (fd2, + " (%7.5g) (%7.5g) (%7.4g) (%7.4g) (%7.4g) (%7.4g) (%7.4g)\n") + call pargr (Memr[xge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[conte+(i-1)*mc_n+mc_sig]) + call pargr (Memr[fluxe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[eqwe+(i-1)*mc_n+mc_sig]) + call pargr (Memr[yge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[sge+(i-1)*mc_n+mc_sig]) + call pargr (Memr[lge+(i-1)*mc_n+mc_sig]) + } + } + } + } else { + call mfree (xg, TY_REAL) + call mfree (yg, TY_REAL) + call mfree (sg, TY_REAL) + call mfree (lg, TY_REAL) + call mfree (pg, TY_INT) + ng = 0 + } + + call mfree (waves, TY_REAL) + call mfree (peaks, TY_REAL) + call mfree (gfwhms, TY_REAL) + call mfree (lfwhms, TY_REAL) + call mfree (types, TY_INT) + call sfree (sp) +end + + +# SUBBLEND -- Subtract last fit. + +procedure subblend (sh, gfd, x, y, n, wx1, wy1, xg, yg, sg, lg, pg, ng) + +pointer sh # SHDR pointer +pointer gfd # GIO file descriptor +real wx1, wy1 # Cursor position +real x[n] # Spectrum data +real y[n] # Spectrum data +int n # Number of points +pointer xg, yg, sg, lg, pg # Pointers to fit parameters +int ng # Number of components + +int i, j, i1, npts, wc, key, nsub +real wx2, wy2, dw +pointer sp, cmd + +int clgcur() +real model() + +begin + if (ng == 0) + return + + call smark (sp) + call salloc (cmd, SZ_FNAME, TY_CHAR) + + # Determine fit range + call printf ("- again:") + call flush (STDOUT) + if (clgcur ("cursor", wx2, wy2, wc, key, Memc[cmd], SZ_FNAME) == EOF) { + call sfree (sp) + return + } + + call fixx (sh, wx1, wx2, wy1, wy2, i1, j) + npts = j - i1 + 1 + + dw = (x[n] - x[1]) / (n - 1) + nsub = NSUB + do i = 1, npts { + y[i1+i-1] = y[i1+i-1] - + model (x[i1+i-1], dw, nsub, Memr[xg], Memr[yg], Memr[sg], + Memr[lg], Memi[pg], ng) + } + + # Plot subtracted curve + call gpline (gfd, x[i1], y[i1], npts) + call gflush (gfd) + + call mfree (xg, TY_REAL) + call mfree (yg, TY_REAL) + call mfree (sg, TY_REAL) + call mfree (lg, TY_REAL) + call mfree (pg, TY_INT) + ng = 0 + call sfree (sp) +end + + +# GASDEV -- Return a normally distributed deviate with zero mean and unit +# variance. The method computes two deviates simultaneously. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +real procedure gasdev (seed) + +long seed # Seed for random numbers + +real v1, v2, r, fac, urand() +int iset +data iset/0/ + +begin + if (iset == 0) { + repeat { + v1 = 2 * urand (seed) - 1. + v2 = 2 * urand (seed) - 1. + r = v1 ** 2 + v2 ** 2 + } until ((r > 0) && (r < 1)) + fac = sqrt (-2. * log (r) / r) + + iset = 1 + return (v1 * fac) + } else { + iset = 0 + return (v2 * fac) + } +end diff --git a/noao/onedspec/splot/splabel.x b/noao/onedspec/splot/splabel.x new file mode 100644 index 00000000..fcba2584 --- /dev/null +++ b/noao/onedspec/splot/splabel.x @@ -0,0 +1,112 @@ +include <gset.h> +include <smw.h> + +define OPTIONS "|label|mabove|mbelow|" +define LABEL 1 # Label at cursor position +define MABOVE 2 # Tick mark plus label above spectrum +define MBELOW 3 # Tick mark plus label below spectrum + + +# SPLABEL -- Add a label. + +procedure splabel (option, sh, gp, x, y, label, format) + +char option[ARB] #I Label option +pointer sh #I Spectrum object +pointer gp #I Graphics object +real x, y #I Label position +char label[ARB] #I Label +char format[ARB] #I Format + +int op, pix, color, markcolor, strdic(), gstati() +real mx, my, x1, x2, y1, y2 +pointer sp, fmt, lab +double shdr_wl() + +define TICK .03 # Tick size in NDC +define GAP .02 # Gap size in NDC + +begin + call smark (sp) + call salloc (fmt, SZ_LINE, TY_CHAR) + call salloc (lab, SZ_LINE, TY_CHAR) + + op = strdic (option, Memc[lab], SZ_LINE, OPTIONS) + if (op == 0) { + call sfree (sp) + return + } + call ggwind (gp, x1, x2, y1, y2) + if ((x < min (x1, x2)) || (x > max (x1, x2))) { + call sfree (sp) + return + } + + # Set label position and draw tick mark. + switch (op) { + case LABEL: + call gctran (gp, x, y, mx, my, 1, 0) + call gctran (gp, mx, my, x1, y2, 0, 1) + markcolor = gstati (gp, G_TICKLABELCOLOR) + if (format[1] == EOS) + call strcpy ("h=c;v=c;s=1.0", Memc[fmt], SZ_LINE) + else + call strcpy (format, Memc[fmt], SZ_LINE) + + case MABOVE: + pix = max (2, min (SN(sh)-3, int (shdr_wl (sh, double (x))))) + y1 = max (Memr[SY(sh)+pix-2], Memr[SY(sh)+pix-1], + Memr[SY(sh)+pix], Memr[SY(sh)+pix+1]) + call gctran (gp, x, y1, mx, my, 1, 0) + call gctran (gp, mx, my + GAP, x1, y1, 0, 1) + call gctran (gp, mx, my + GAP + TICK, x1, y2, 0, 1) + + color = gstati (gp, G_PLCOLOR) + markcolor = gstati (gp, G_TICKLABELCOLOR) + call gseti (gp, G_PLCOLOR, markcolor) + call gline (gp, x1, y1, x1, y2) + call gseti (gp, G_PLCOLOR, color) + + call gctran (gp, mx, my + TICK + 2 * GAP, x1, y2, 0, 1) + if (format[1] == EOS) + call strcpy ("u=180;h=c;v=b;s=0.5", Memc[fmt], SZ_LINE) + else + call strcpy (format, Memc[fmt], SZ_LINE) + + case MBELOW: + pix = max (2, min (SN(sh)-3, int (shdr_wl (sh, double (x))))) + y1 = min (Memr[SY(sh)+pix-2], Memr[SY(sh)+pix-1], + Memr[SY(sh)+pix], Memr[SY(sh)+pix+1]) + call gctran (gp, x, y1, mx, my, 1, 0) + call gctran (gp, mx, my - GAP, x1, y1, 0, 1) + call gctran (gp, mx, my - GAP - TICK, x1, y2, 0, 1) + + color = gstati (gp, G_PLCOLOR) + markcolor = gstati (gp, G_TICKLABELCOLOR) + call gseti (gp, G_PLCOLOR, markcolor) + call gline (gp, x1, y1, x1, y2) + call gseti (gp, G_PLCOLOR, color) + + call gctran (gp, mx, my - TICK - 2 * GAP, x1, y2, 0, 1) + if (format[1] == EOS) + call strcpy ("u=0;h=c;v=t;s=0.5", Memc[fmt], SZ_LINE) + else + call strcpy (format, Memc[fmt], SZ_LINE) + } + + # Draw the label. + if (label[1] != EOS) { + color = gstati (gp, G_TXCOLOR) + call gseti (gp, G_TXCOLOR, markcolor) + if (label[1] == '%') { + call sprintf (Memc[lab], SZ_LINE, label) + call pargr (x) + call gtext (gp, x1, y2, Memc[lab], Memc[fmt]) + } else + call gtext (gp, x1, y2, label, Memc[fmt]) + call gseti (gp, G_TXCOLOR, color) + } + + call gflush (gp) + call sfree (sp) +end diff --git a/noao/onedspec/splot/splot.key b/noao/onedspec/splot/splot.key new file mode 100644 index 00000000..b78c722d --- /dev/null +++ b/noao/onedspec/splot/splot.key @@ -0,0 +1,116 @@ +? - This display r - Redraw the current window +/ - Cycle thru short help on stat line s - Smooth (boxcar) +a - Autoexpand between cursors t - Fit continuum(*) +b - Toggle base plot level to 0.0 u - Adjust coordinate scale(*) +c - Clear and redraw full spectrum v - Velocity scale (toggle) +d - Deblend lines using profile models w - Window the graph +e - Equiv. width, integ flux, center x - Connects 2 cursor positions +f - Arithmetic functions: log, sqrt... y - Plot std star flux from calib file +g - Get new image and plot z - Expand x range by factor of 2 +h - Equivalent widths(*) ) - Go to next spectrum in image +i - Write current image as new image ( - Go to previous spectrum in image +j - Fudge a point to Y-cursor value # - Select new line/aperture +k - Profile fit to single line(*) % - Select new band +l - Convert to F-lambda $ - Toggle wavelength/pixel scale +m - Mean, RMS, snr in marked region - - Subtract deblended fit +n - Convert to F-nu , - Down slide spectrum +o - Toggle overplot of following plot . - Up slide spectrum +p - Convert to wavelength scale I - Interrupt task immediately +q - Quit and exit <space> - Cursor position and flux + +(*) For 'h' key: Measure equivalent widths + a - Left side for width at 1/2 flux l - Left side for continuum = 1 + b - Right side for width at 1/2 flux r - Right side for continuum = 1 + c - Both sides for width at 1/2 flux k - Both sides for continuum = 1 + +(*) For 'k' key: Second key may be used to select profile type + g - Gaussian, l - Lorentzian, v - Voigt, all others - Gaussiank + +(*) For 't' key: Fit the continuum with ICFIT and apply to spectrum + / = normalize by the continuum fit + - = subtract the continuum fit (residuals) + f = replace spectrum by the continuum fit + c = clean spectrum of rejected points + n = do the fitting but leave the spectrum unchanged + q = quit without fitting or modifying spectrum + +(*) For 'u' key: Adjust the coordinate scale by marking features + d = apply doppler correction to bring marked feature to specified coord. + l = set linear (in wavelength) coordinates based on two marked features + z = apply zero point shift to bring marked feature to + specified coordinate + +The colon commands do not allow abbreviations. + +:# <comment> - Add comment to log file +:dispaxis <val> - Change summing parameter for 2D images +:log - Enable logging to save_file +:nolog - Disable logging to save_file +:nsum <val> - Change summing parameter for 2D images +:show - Show full output of deblending and equiv. width measurments +:units <value> - Change coordinate units (see below) + +:label <label> <format> - Add label at cursor position +:mabove <label> <format> - Add tick mark and label above spectrum +:mbelow <label> <format> - Add tick mark and label below spectrum + The label must be quoted if it contains blanks. A label beginning + with % (i.e. %.2f) is treated as a format for the x cursor position. + The optional format is a gtext string (see help on "cursors"). + The labels are not remembered between redraws. + +:auto [yes|no] - Enable/disable autodraw option +:zero [yes|no] - Enable/disable zero baseline option +:xydraw [yes|no] - Enable/disable xydraw option +:hist [yes|no] - Enable/disable histogram line type option +:nosysid [yes|no] - Enable/disable system ID option +:wreset [yes|no] - Enable/disable window reset for new spectra option +:flip [yes|no] - Enable/disable dispersion coordinate flip +:overplot [yes|no]- Enable/disable permanent overplot mode + +:/help Get help on GTOOLS options +:.help Get help on cursor mode options + + + UNITS + +The units are specified by strings having a unit type from the list below +along with the possible preceding modifiers, "inverse", to select the +inverse of the unit and "log" to select logarithmic units. For example "log +angstroms" to plot the logarithm of wavelength in Angstroms and "inv +microns" to plot inverse microns. The various identifiers may be +abbreviated as words but the syntax is not sophisticated enough to +recognized standard scientific abbreviations except as noted below. + + angstroms - Wavelength in Angstroms + nanometers - Wavelength in nanometers + millimicrons - Wavelength in millimicrons + microns - Wavelength in microns + millimeters - Wavelength in millimeters + centimeter - Wavelength in centimeters + meters - Wavelength in meters + hertz - Frequency in hertz (cycles per second) + kilohertz - Frequency in kilohertz + megahertz - Frequency in megahertz + gigahertz - Frequency in gigahertz + m/s - Velocity in meters per second + km/s - Velocity in kilometers per second + ev - Energy in electron volts + kev - Energy in kilo electron volts + mev - Energy in mega electron volts + z - Redshift + + nm - Wavelength in nanometers + mm - Wavelength in millimeters + cm - Wavelength in centimeters + m - Wavelength in meters + Hz - Frequency in hertz (cycles per second) + KHz - Frequency in kilohertz + MHz - Frequency in megahertz + GHz - Frequency in gigahertz + wn - Wave number (inverse centimeters) + +The velocity and redshift units require a trailing value and unit defining the +velocity zero point. For example to plot velocity relative to +a wavelength of 1 micron the unit string would be: + + km/s 1 micron diff --git a/noao/onedspec/splot/splot.log b/noao/onedspec/splot/splot.log new file mode 100644 index 00000000..1efd91c0 --- /dev/null +++ b/noao/onedspec/splot/splot.log @@ -0,0 +1,8 @@ + +Feb 6 15:44 [tofu$s2ndR136002]: NDR136002[1/1] + 4103.416 5.470E-14 7.082E-14 -1.296 + 4060.36 5.637E-14 -6.87E-14 1.218 + +Feb 6 15:45 [tofu$s2ndR136002]: NDR136002[1/1] + 964.983 5.606E-14 1.030E-13 -1.84 + 907.3833 5.425E-14 -8.26E-14 1.525 diff --git a/noao/onedspec/splot/splot.x b/noao/onedspec/splot/splot.x new file mode 100644 index 00000000..4b676660 --- /dev/null +++ b/noao/onedspec/splot/splot.x @@ -0,0 +1,605 @@ +include <error.h> +include <imhdr.h> +include <gset.h> +include <pkg/gtools.h> +include <smw.h> +include <units.h> + +define KEY "noao$onedspec/splot/splot.key" +define HELP "noao$onedspec/splot/stshelp.key" +define PROMPT "splot options" + +define OPTIONS ",auto,zero,xydraw,histogram,nosysid,wreset,flip,overplot," +define NOPTIONS 8 +define AUTO 1 # Option number for auto graph +define ZERO 2 # Option number of zero y minimum +define XYDRAW 3 # Draw connection X,Y pairs +define HIST 4 # Draw histogram style lines +define NOSYSID 5 # Don't include system id +define WRESET 6 # Reset window for each new spectrum +define FLIP 7 # Flip spectra +define OVERPLOT 8 # Overplot toggle + + +# SPLOT -- Plot an image line and play with it - Most appropriate for spectra + +procedure splot () + +int list +int i, j, npts, nline, nband, nap +int wc, key, keyu +real wx, wy +double w1, u1, w2, u2, w0, wpc +real avg_pix, sigma_pix + +int fd1, fd2, ng, hline, hlines +int newgraph, newimage, overplot, options[NOPTIONS] +pointer sp, image, units, units1, units2, units3, cmd, save1, save2 +pointer gp, gt, im, mw, x, y, sh, xg, yg, sg, lg, pg, hptr +bool wave_scl, fnu + +pointer gopen(), gt_init() +int clgcur(), imtopen(), imtgetim(), imaccess(), gt_geti(), nowhite() +real clgetr(), gt_getr() +double clgetd(), shdr_wl() +bool streq(), fp_equald() +errchk getimage, fun_do, ans_hdr, un_changer + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (units, SZ_FNAME, TY_CHAR) + call salloc (units1, SZ_FNAME, TY_CHAR) + call salloc (units2, SZ_FNAME, TY_CHAR) + call salloc (units3, SZ_FNAME, TY_CHAR) + call salloc (cmd, SZ_FNAME, TY_CHAR) + call salloc (save1, SZ_FNAME, TY_CHAR) + call salloc (save2, SZ_FNAME, TY_CHAR) + + # Get task parameters. + + call clgstr ("images", Memc[image], SZ_FNAME) + list = imtopen (Memc[image]) + call clgstr ("save_file", Memc[save1], SZ_FNAME) + call clgstr ("options", Memc[save2], SZ_FNAME) + call xt_gids (Memc[save2], OPTIONS, options, NOPTIONS) + call clgstr ("units", Memc[units], SZ_FNAME) + call mktemp ("tmp$splot", Memc[save2], SZ_FNAME) + + # Allocate space for User area + x = NULL + y = NULL + + # Initialize + im = NULL + sh = NULL + fd1 = NULL + fd2 = NULL + hptr = NULL + ng = 0 + hline = 1 + nline = 0 + nband = 0 + nap = 0 + if (nowhite (Memc[units], Memc[units1], SZ_FNAME) == 0) + call strcpy ("display", Memc[units], SZ_FNAME) + call strcpy (Memc[units], Memc[units1], SZ_FNAME) + call strcpy (Memc[units], Memc[units2], SZ_FNAME) + w0 = INDEFD + wpc = INDEFD + + call clgstr ("graphics", Memc[cmd], SZ_FNAME) + gp = gopen (Memc[cmd], NEW_FILE+AW_DEFER, STDGRAPH) + call gseti (gp, G_WCS, 1) +# call gseti (gp, G_YNMINOR, 0) + + gt = gt_init() + call gt_setr (gt, GTXMIN, clgetr ("xmin")) + call gt_setr (gt, GTXMAX, clgetr ("xmax")) + call gt_setr (gt, GTYMIN, clgetr ("ymin")) + call gt_setr (gt, GTYMAX, clgetr ("ymax")) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + if (options[HIST] == YES) + call gt_sets (gt, GTTYPE, "histogram") + else + call gt_sets (gt, GTTYPE, "line") + call gt_seti (gt, GTXFLIP, options[FLIP]) + if (options[NOSYSID] == YES) + call gt_seti (gt, GTSYSID, NO) + + while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) { + + # Initialize to plot a wavelength scale + wave_scl = true + + # Open image and get pixels + if (imaccess (Memc[image], READ_ONLY) == NO) { + call eprintf ("Can't get image %s\n") + call pargstr (Memc[image]) + next + } + call getimage (Memc[image], nline, nband, nap, wave_scl, + w0, wpc, Memc[units], im, mw, sh, gt) + x = SX(sh) + y = SY(sh) + npts = SN(sh) + newimage = YES + overplot = options[OVERPLOT] + + # Enter cursor loop with 'r' redraw. + key = 'r' + repeat { + switch (key) { + case ':': + if (Memc[cmd] == '/') + call gt_colon (Memc[cmd], gp, gt, newgraph) + else { + call splot_colon (Memc[cmd], options, gp, gt, sh, + wx, wy, Memc[units], Memc[save1], Memc[save2], + fd1, fd2, newgraph) + overplot = options[OVERPLOT] + if (sh == NULL) { + call getimage (Memc[image], nline, nband, nap, + wave_scl, w0, wpc, Memc[units], im, mw, sh, gt) + x = SX(sh) + y = SY(sh) + npts = SN(sh) + newgraph = YES + newimage = YES + } + } + + case 'a': # Autoexpand + call auto_exp (gp, gt, key, wx, Memr[x], Memr[y], npts) + + case 'b': # Toggle base to 0.0 + if (options[ZERO] == NO) { + call gt_setr (gt, GTYMIN, 0.) + options[ZERO] = YES + } else { + call gt_setr (gt, GTYMIN, INDEF) + options[ZERO] = NO + } + newgraph = options[AUTO] + overplot = NO + + case 'c': + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + newgraph = YES + overplot = NO + + case 'd': # De-blend a group of lines + call ans_hdr (sh, newimage, key, Memc[save1], Memc[save2], + fd1, fd2) + call sp_deblend (sh, gp, wx, wy, Memr[x], Memr[y], npts, + fd1, fd2, xg, yg, sg, lg, pg, ng) + newimage = NO + + case 'k': # Fit gaussian + call ans_hdr (sh, newimage, key, Memc[save1], Memc[save2], + fd1, fd2) + call gfit (sh, gp, wx, wy, Memr[x], Memr[y], npts, + fd1, fd2, xg, yg, sg, lg, pg, ng) + newimage = NO + + case 'e': # Equivalent width + call ans_hdr (sh, newimage, key, Memc[save1], Memc[save2], + fd1, fd2) + call eqwidth (sh, gp, wx, wy, Memr[x], Memr[y], npts, + fd1, fd2) + newimage = NO + + case 'v': + iferr { + if (UN_CLASS(UN(sh)) == UN_VEL) { + call strcpy (Memc[units1], Memc[units], SZ_FNAME) + call strcpy (Memc[units2], Memc[units3], SZ_FNAME) + } else { + call strcpy (Memc[units], Memc[units1], SZ_FNAME) + call strcpy (UNITS(sh), Memc[units2], SZ_FNAME) + call un_changer (UN(sh), "angstroms", wx, 1, NO) + call sprintf (Memc[units], SZ_FNAME, + "km/s %g angstroms") + call pargr (wx) + call strcpy (Memc[units], Memc[units3], SZ_FNAME) + } + wx = gt_getr (gt, GTXMIN) + if (!IS_INDEF(wx)) { + call un_changer (UN(sh), Memc[units3], wx, 1, NO) + call gt_setr (gt, GTXMIN, wx) + } + wx = gt_getr (gt, GTXMAX) + if (!IS_INDEF(wx)) { + call un_changer (UN(sh), Memc[units3], wx, 1, NO) + call gt_setr (gt, GTXMAX, wx) + } + call un_changer (UN(sh), Memc[units3], Memr[x], npts, + YES) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + newgraph = YES + overplot = NO + } then + call erract (EA_WARN) + + case 'h': # Equivalent widths -- C. Pilachowski style + call ans_hdr (sh, newimage, key, Memc[save1], Memc[save2], + fd1, fd2) + repeat { + switch (key) { + case 'a', 'b', 'c': # Continuum at cursor width at 1/2 + call eqwidth_cp (sh, gp, wx, wy, INDEF, + Memr[y], npts, key, fd1, fd2, xg, yg, sg, + lg, pg, ng) + break + case 'l', 'r', 'k': # Continuum at 1 + call eqwidth_cp (sh, gp, wx, 1., wy, + Memr[y], npts, key, fd1, fd2, xg, yg, sg, + lg, pg, ng) + break + default: + call printf ( + "Set cursor and type a, b, c, l, r, or k:") + } + } until (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], + SZ_FNAME) == EOF) + newimage = NO + + case 'o': # Set overplot + overplot = YES + + case 'g', '#', '%', '(', ')': # Get new image to plot + i = nline + j = nband + + switch (key) { + case '(': + if (IM_LEN(im,2) > 1) { + nline = max (1, min (IM_LEN(im,2), nline-1)) + nap = INDEFI + } else if (IM_LEN(im,3) > 1) { + nband = max (1, min (IM_LEN(im,3), nband-1)) + } + case ')': + if (IM_LEN(im,2) > 1) { + nline = max (1, min (IM_LEN(im,2), nline+1)) + nap = INDEFI + } else if (IM_LEN(im,3) > 1) { + nband = max (1, min (IM_LEN(im,3), nband+1)) + } + case '#': + nline = 0 + case '%': + nband = 0 + default: + call clgstr ("next_image", Memc[cmd], SZ_FNAME) + if (streq (Memc[image], Memc[cmd])) { + call shdr_close (sh) + } else if (imaccess (Memc[cmd], READ_ONLY) == YES) { + call shdr_close (sh) + call smw_close (mw) + call imunmap (im) + newimage = YES + } else { + call eprintf ("Can't get %s\n") + call pargstr (Memc[cmd]) + next + } + call strcpy (Memc[cmd], Memc[image], SZ_FNAME) + nline = 0 + nband = 0 + } + + call getimage (Memc[image], nline, nband, nap, wave_scl, + w0, wpc, Memc[units], im, mw, sh, gt) + x = SX(sh) + y = SY(sh) + npts = SN(sh) + + if (options[WRESET] == YES && overplot == NO) { + call gt_setr (gt, GTXMIN, clgetr ("xmin")) + call gt_setr (gt, GTXMAX, clgetr ("xmax")) + call gt_setr (gt, GTYMIN, clgetr ("ymin")) + call gt_setr (gt, GTYMAX, clgetr ("ymax")) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + } + + if (nline != i || nband != j) + newimage = YES + newgraph = YES + + case 'w': # Window the graph + call gt_window (gt, gp, "cursor", newgraph) + if (newgraph == YES) { + newgraph = options[AUTO] + overplot = NO + } + + case 'l': # Convert to f-lambda - issue warning if not a + # calibrated image + if (FC(sh) == FCNO) + call eprintf ( + "Warning: (>flam) spectrum not calibrated\n") + + call conflam (sh) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + newgraph = options[AUTO] + overplot = NO + + case 'f': # Function operators + call fun_help () + while (clgcur ("cursor", wx, wy, wc, key, Memc[cmd], + SZ_FNAME) != EOF) { + switch (key) { + case '?': + call fun_help () + case 'q': + break + case 'I': + call fatal (0, "Interrupt") + default: + iferr { + call fun_do (key, sh, Memr[y], npts, w0, wpc) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + if (options[AUTO] == YES) + call replot (gp, gt, Memr[x], Memr[y], + npts, YES) + overplot = NO + call fun_help () + } then { + call erract (EA_WARN) + call tsleep (2) + call fun_help () + } + } + } + call printf ("\n") + + case 'm': # Signal-to-noise + call ans_hdr (sh, newimage, key, Memc[save1], Memc[save2], + fd1, fd2) + call avgsnr (sh, wx, wy, Memr[y], npts, fd1, fd2) + newimage = NO + + case 'n': # Convert to f-nu + if (FC(sh) == FCNO) + call eprintf ( + "Warning: (>fnu) spectrum not calibrated\n") + + call confnu (sh) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + newgraph = options[AUTO] + overplot = NO + + case 'q': + if (options[WRESET] == YES) { + call gt_setr (gt, GTXMIN, clgetr ("xmin")) + call gt_setr (gt, GTXMAX, clgetr ("xmax")) + call gt_setr (gt, GTYMIN, clgetr ("ymin")) + call gt_setr (gt, GTYMAX, clgetr ("ymax")) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + } + + if (nline != i || nband != j) + newimage = YES + newgraph = YES + break + + case 'r': # Replot + newgraph = YES + overplot = NO + + case 's': # Smooth + call smooth (Memr[y], npts) + newgraph = options[AUTO] + + case 't': # FlaTTen spectrum + call flatten (gp, gt, Memr[x], Memr[y], npts) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + newgraph = options[AUTO] + overplot = NO + + case 'p', 'u': # Set user coordinates + if (!wave_scl) { + call shdr_system (sh, "world") + wave_scl = true + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + } + switch (key) { + case 'p': + keyu = 'l' + w1 = Memr[x] + u1 = clgetd ("wstart") + w2 = Memr[x+npts-1] + u2 = clgetd ("wend") + if (IS_INDEFD(u1)) { + u1 = clgetd ("dw") + u1 = u2 - (npts - 1) * u1 + } else if (IS_INDEFD(u2)) { + u2 = clgetd ("dw") + u2 = u1 + (npts - 1) * u2 + } + case 'u': + call printf ( + "Set cursor and select correction: d(oppler), z(eropoint), l(inear)\n") + call flush (STDOUT) + i = clgcur ("cursor", wx, wy, wc, keyu, Memc[cmd], + SZ_FNAME) + w1 = wx + u1 = clgetd ("wavelength") + if (keyu == 'l') { + repeat { + call printf ("Set cursor to second position:") + call flush (STDOUT) + i = clgcur ("cursor", wx, wy, wc, key, + Memc[cmd], SZ_FNAME) + w2 = wx + if (!fp_equald (w1, w2)) { + u2 = clgetd ("wavelength") + break + } + call printf ("Cursor not moved: ") + } + } + } + call usercoord (sh, keyu, w1, u1, w2, u2) + w0 = Memr[x] + wpc = Memr[x+1] - w0 + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + newgraph = options[AUTO] + overplot = NO + + case 'i': # Write image spectrum out + call sp_wrspect (sh) + im = IM(sh) + mw = MW(sh) + + case 'j': # Fudge (fix) a data point + call fudgept (sh, gp, Memr[x], Memr[y], npts, wx, wy) + + case 'x': # Fudge eXtended over a line + call fudgex (sh, gp, Memr[x], Memr[y], npts, wx, wy, + options[XYDRAW]) + + case 'y': # Over plot standard star data + # Estimate data is fnu or flambda: cutoff around dexp[-20] + fnu = false + call aavgr (Memr[y], npts, avg_pix, sigma_pix) + if (log10 (avg_pix) < -19.5) + fnu = true + call plot_std (sh, gp, fnu) + call printf ("\n") + + case 'z': # Zoom x region to larger range + call auto_exp (gp, gt, key, wx, Memr[x], Memr[y], npts) + + case '-': # Subtract deblended fit + call subblend (sh, gp, Memr[x], Memr[y], npts, wx, wy, + xg, yg, sg, lg, pg, ng) + + case '.': # Slide upward + call auto_exp (gp, gt, key, wx, Memr[x], Memr[y], npts) + + case ',': # Slide downward + call auto_exp (gp, gt, key, wx, Memr[x], Memr[y], npts) + + case '$': # Toggle wavelength scale + if (wave_scl) { + call shdr_system (sh, "physical") + wave_scl = false + call gt_sets (gt, GTXLABEL, "Pixel") + call gt_sets (gt, GTXUNITS, "") + } else { + call shdr_system (sh, "world") + wave_scl = true + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + } + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + if (options[ZERO] == YES) + call gt_setr (gt, GTYMIN, 0.) + newgraph = options[AUTO] + overplot = NO + + case '/': # Help on status line + call sts_help (hline, hlines, HELP, hptr) + hline = mod (hline, hlines) + 1 + + case '?': # Help screen + call gpagefile (gp, KEY, PROMPT) + + case 'I': # Interrupt + call fatal (0, "Interrupt") + + default: # Default - print cursor info + i = max (1, min (npts, nint (shdr_wl (sh, double(wx))))) + call printf ("x,y,z(x): %10.3f %10.4g %10.4g\n") + call pargr (wx) + call pargr (wy) + call pargr (Memr[y+i-1]) + } + + if (newgraph == YES) { + if (overplot == YES) { + call printf ("Overplotting: %s") + call pargstr (Memc[image]) + if (nline > 0) { + if (nband > 0) { + call printf ("(%d,%d)") + call pargi (nline) + call pargi (nband) + } else { + call printf ("(%d)") + call pargi (nline) + } + } + call flush (STDOUT) + i = gt_geti (gt, GTLINE) + j = gt_geti (gt, GTCOLOR) + if (options[OVERPLOT] == NO) { + call gt_seti (gt, GTLINE, i+1) + call gt_seti (gt, GTCOLOR, j+1) + } + call replot (gp, gt, Memr[x], Memr[y], npts, NO) + call gt_seti (gt, GTLINE, i) + call gt_seti (gt, GTCOLOR, j) + } else + call replot (gp, gt, Memr[x], Memr[y], npts, YES) + newgraph = NO + overplot = options[OVERPLOT] + } + } until (clgcur ("cursor",wx,wy,wc,key,Memc[cmd],SZ_FNAME) == EOF) + if (im != ERR) { + call shdr_close (sh) + call smw_close (mw) + call imunmap (im) + } + } + + call gclose (gp) + if (fd1 != NULL) + call close (fd1) + if (fd2 != NULL) { + call close (fd2) + call delete (Memc[save2]) + } + if (hptr != NULL) + call mfree (hptr, TY_CHAR) + if (ng > 0) { + call mfree (xg, TY_REAL) + call mfree (yg, TY_REAL) + call mfree (sg, TY_REAL) + call mfree (lg, TY_REAL) + call mfree (pg, TY_INT) + } + call smw_daxis (NULL, NULL, 0, 0, 0) + call gt_free (gt) + call imtclose (list) +end diff --git a/noao/onedspec/splot/splotcolon.x b/noao/onedspec/splot/splotcolon.x new file mode 100644 index 00000000..e68bbecc --- /dev/null +++ b/noao/onedspec/splot/splotcolon.x @@ -0,0 +1,263 @@ +include <error.h> +include <pkg/gtools.h> +include <smw.h> +include <units.h> +include <ctype.h> + +# List of colon commands. +define CMDS "|show|nolog|log|dispaxis|nsum|#|units|auto|zero\ + |xydraw|histogram|nosysid|wreset|flip|overplot\ + |label|mabove|mbelow|" +define SHOW 1 # Show logged data +define NOLOG 2 # Turn off logging +define LOG 3 # Turn on logging +define DA 4 # Dispersion axis +define NS 5 # Summing parameter +define COMMENT 6 # Comment +define UNITS 7 # Units +define AUTO 8 # Option auto graph +define ZERO 9 # Option for zero y minimum +define XYDRAW 10 # Draw connection X,Y pairs +define HIST 11 # Draw histogram style lines +define NOSYSID 12 # Don't include system id +define WRESET 13 # Reset window for each new spectrum +define FLIP 14 # Flip the dispersion coordinates +define OVERPLOT 15 # Toggle overplot +define LABEL 16 # Label spectrum +define MABOVE 17 # Tick mark plus label above spectrum +define MBELOW 18 # Tick mark plus label below spectrum + +define OPOFF 7 # Offset in options array + +# SPLOT_COLON -- Respond to colon command. + +procedure splot_colon (command, options, gp, gt, sh, x, y, units, + fname1, fname2, fd1, fd2, newgraph) + +char command[ARB] # Colon command +int options[ARB] # Options +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer sh # SHIO pointer +real x, y # Coordinate +char units[SZ_FNAME] # Units string +char fname1[SZ_FNAME] # Log file +char fname2[SZ_FNAME] # Temporary log file +int fd1, fd2 # Log file descriptors +int newgraph # New graph? + +bool bval +char cmd[SZ_LINE] +real xval, gt_getr() +int ncmd, ival, access(), nscan(), strdic(), btoi(), gt_geti() +pointer sp, str, smw +errchk un_changer + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Scan the command string and get the first word. + call sscan (command) + call gargwrd (cmd, SZ_LINE) + ncmd = strdic (cmd, cmd, SZ_LINE, CMDS) + + smw = MW(sh) + + switch (ncmd) { + case SHOW: + if (fd2 != NULL) { + call close (fd2) + fd2 = NULL + } + if (access (fname2, 0, 0) == YES) + call gpagefile (gp, fname2, "splot data") + else + call printf ("No measurements\n") + case NOLOG: + call printf ("Logging to %s disabled") + call pargstr (fname1) + fname1[1] = EOS + if (fd1 != NULL) { + call close (fd1) + fd1 = NULL + } + case LOG: + call clgstr ("save_file", fname1, SZ_FNAME) + call printf ("Logging to %s enabled") + call pargstr (fname1) + case DA: + if (SMW_FORMAT(smw) == SMW_ND) { + call gargi (ival) + if (nscan() == 2) { + if (ival < 1) { + call printf ("Bad value for dispaxis (%d)\n") + call pargi (ival) + } else if (ival != SMW_PAXIS(smw,1)) { + call smw_daxis (smw, IM(sh), ival, INDEFI, INDEFI) + call smw_saxes (smw, NULL, IM(sh)) + call shdr_close (sh) + } + } else { + call printf ("dispaxis %d\n") + call pargi (SMW_PAXIS(smw,1)) + } + } else + call printf ("Image is not two dimensional\n") + case NS: + if (SMW_FORMAT(smw) == SMW_ND) { + call gargi (ival) + call gargi (ncmd) + if (nscan() == 1) { + call printf ("nsum %d %d\n") + call pargi (SMW_NSUM(smw,1)) + call pargi (SMW_NSUM(smw,2)) + } else { + if (nscan() == 2) + ncmd = INDEFI + if ((!IS_INDEFI(ival) && ival != SMW_NSUM(smw,1)) || + (!IS_INDEFI(ncmd) && ncmd != SMW_NSUM(smw,2))) { + call smw_daxis (smw, IM(sh), INDEFI, ival, ncmd) + call smw_saxes (smw, NULL, IM(sh)) + call shdr_close (sh) + } + } + } else + call printf ("Invalid image format\n") + case COMMENT: + call ans_hdr (sh, NO, 'm', fname1, fname2, fd1, fd2) + call gargstr (cmd, SZ_LINE) + if (fd1 != NULL) { + call fprintf (fd1, "%s\n") + call pargstr (command) + } + if (fd2 != NULL) { + call fprintf (fd2, "%s\n") + call pargstr (command) + } + case UNITS: + call gargstr (cmd, SZ_LINE) + for (ival=1; IS_WHITE(cmd[ival]); ival=ival+1) + ; + iferr { + xval = gt_getr (gt, GTXMIN) + if (!IS_INDEF(xval)) { + call un_changer (UN(sh), cmd[ival], xval, 1, NO) + call gt_setr (gt, GTXMIN, xval) + } + xval = gt_getr (gt, GTXMAX) + if (!IS_INDEF(xval)) { + call un_changer (UN(sh), cmd[ival], xval, 1, NO) + call gt_setr (gt, GTXMAX, xval) + } + call un_changer (UN(sh), cmd[ival], Memr[SX(sh)], SN(sh), YES) + call strcpy (cmd[ival], units, SZ_FNAME) + call gt_sets (gt, GTXLABEL, UN_LABEL(UN(sh))) + call gt_sets (gt, GTXUNITS, UN_UNITS(UN(sh))) + newgraph = YES + } then + call erract (EA_WARN) + case AUTO: + call gargb (bval) + if (nscan() == 2) + options[AUTO-OPOFF] = btoi (bval) + else { + call printf ("auto %b\n") + call pargi (options[AUTO-OPOFF]) + } + case ZERO: + call gargb (bval) + if (nscan() == 2) { + options[ZERO-OPOFF] = btoi (bval) + if (bval) + call gt_setr (gt, GTYMIN, 0.) + newgraph = options[AUTO-OPOFF] + } else { + call printf ("zero %b\n") + call pargi (options[ZERO-OPOFF]) + } + case XYDRAW: + call gargb (bval) + if (nscan() == 2) + options[XYDRAW-OPOFF] = btoi (bval) + else { + call printf ("xydraw %b\n") + call pargi (options[XYDRAW-OPOFF]) + } + case HIST: + call gargb (bval) + if (nscan() == 2) { + options[HIST-OPOFF] = btoi (bval) + if (bval) + call gt_sets (gt, GTTYPE, "histogram") + else + call gt_sets (gt, GTTYPE, "line") + newgraph = options[AUTO-OPOFF] + } else { + call printf ("hist %b\n") + call pargi (options[HIST-OPOFF]) + } + case NOSYSID: + call gargb (bval) + if (nscan() == 2) { + options[NOSYSID-OPOFF] = btoi (bval) + if (bval) + call gt_seti (gt, GTSYSID, NO) + else + call gt_seti (gt, GTSYSID, YES) + newgraph = options[AUTO-OPOFF] + } else { + call printf ("nosysid %b\n") + call pargi (options[NOSYSID-OPOFF]) + } + case WRESET: + call gargb (bval) + if (nscan() == 2) + options[WRESET-OPOFF] = btoi (bval) + else { + call printf ("wreset %b\n") + call pargi (options[WRESET-OPOFF]) + } + case FLIP: + call gargb (bval) + if (nscan() == 2) { + options[FLIP-OPOFF] = btoi (bval) + call gt_seti (gt, GTXFLIP, options[FLIP-OPOFF]) + } else { + options[FLIP-OPOFF] = gt_geti (gt, GTXFLIP) + call printf ("flip %b\n") + call pargi (options[FLIP-OPOFF]) + } + case OVERPLOT: + call gargb (bval) + if (nscan() == 2) { + options[OVERPLOT-OPOFF] = btoi (bval) + } else { + call printf ("overplot %b\n") + call pargi (options[OVERPLOT-OPOFF]) + } + case LABEL, MABOVE, MBELOW: + call gargwrd (cmd, SZ_LINE) + for (ival=1; IS_WHITE(cmd[ival]); ival=ival+1) + ; + call strcpy (cmd[ival], cmd, SZ_LINE) + call gargwrd (Memc[str], SZ_LINE) + for (ival=1; IS_WHITE(Memc[str+ival-1]); ival=ival+1) + ; + call strcpy (Memc[str+ival-1], Memc[str], SZ_LINE) + + switch (ncmd) { + case LABEL: + call splabel ("label", sh, gp, x, y, cmd, Memc[str]) + case MABOVE: + call splabel ("mabove", sh, gp, x, y, cmd, Memc[str]) + case MBELOW: + call splabel ("mbelow", sh, gp, x, y, cmd, Memc[str]) + } + + default: + call printf ("Unrecognized or ambiguous command\007") + } + + call sfree (sp) +end diff --git a/noao/onedspec/splot/splotfun.x b/noao/onedspec/splot/splotfun.x new file mode 100644 index 00000000..4c94350f --- /dev/null +++ b/noao/onedspec/splot/splotfun.x @@ -0,0 +1,127 @@ +include <error.h> +include <mach.h> +include <smw.h> + +# Function Mode for STEK + +# FUN_DO -- Branch and execute proper function + +procedure fun_do (key, sh1, y, n, w0, wpc) + +int key +pointer sh1 +real y[n] +int n +double w0, wpc + +char spec2[SZ_FNAME] +int i, nline, nband, nap, strlen() +real const, clgetr() +pointer im, mw, sh2 +bool wave_scl +errchk getimage, shdr_rebin + +begin + switch (key) { + case 'a': # Absolute value + do i = 1, n + y[i] = abs (y[i]) + case 'd': # Dexp (base 10) + const = log10 (MAX_REAL) + do i = 1, n + if (abs (y[i]) < const) + y[i] = 10.0 ** y[i] + else if (y[i] >= const) + y[i] = MAX_REAL + else + y[i] = 0.0 + case 'e': # Exp base e + const = log (MAX_REAL) + do i = 1, n + if (abs (y[i]) < const) + y[i] = exp (y[i]) + else if (y[i] >= const) + y[i] = MAX_REAL + else + y[i] = 0.0 + case 'i': # Inverse + do i = 1, n + if (y[i] != 0.0) + y[i] = 1.0/y[i] + else + y[i] = 0.0 + case 'l': # Log10 + do i = 1, n + if (y[i] > 0.0) + y[i] = log10 (y[i]) + else + y[i] = -0.5 + case 'm': # Multiply by constant + const = clgetr ("constant") + call amulkr (y, const, y, n) + case 'n': # Log base e + do i = 1, n + if (y[i] > 0.0) + y[i] = log (y[i]) + else + y[i] = -0.5 + case 'p': # Add constant + const = clgetr ("constant") + call aaddkr (y, const, y, n) + case 's': # Square root + do i = 1, n + if (y[i] >= 0.0) + y[i] = sqrt (y[i]) + else + y[i] = 0.0 + + case '+', '-', '*', '/': # Binary operations + call printf ("Second spectrum ") + call clgstr ("spec2", spec2, SZ_FNAME) + if (strlen (spec2) == 0) + return + + wave_scl = true + nline = 0 + nband = 0 + nap = 0 + im = NULL + mw = NULL + sh2 = NULL + call getimage (spec2, nline, nband, nap, wave_scl, w0, wpc, + "angstroms", im, mw, sh2, NULL) + call shdr_rebin (sh2, sh1) + switch (key) { + case '+': + call aaddr (y, Memr[SY(sh2)], y, n) + case '-': + call asubr (y, Memr[SY(sh2)], y, n) + case '*': + call amulr (y, Memr[SY(sh2)], y, n) + case '/': + do i = 1, n + if (Memr[SY(sh2)+i-1] == 0.0) + y[i] = 0.0 + else + y[i] = y[i] / Memr[SY(sh2)+i-1] + } + call shdr_close (sh2) + call smw_close (mw) + call imunmap (im) + + # Redraw + case 'r': + ; + default: + call error (0, "Unknown function") + } +end + +# FUN_HELP + +procedure fun_help () + +begin + call printf ("q=quit l,n=log10,e d,e=d,exp s=sqrt a=abs i=1/s") + call printf (" p=+k m=*k +,-,*,/=2spec ops\n") +end diff --git a/noao/onedspec/splot/stshelp.key b/noao/onedspec/splot/stshelp.key new file mode 100644 index 00000000..fe351182 --- /dev/null +++ b/noao/onedspec/splot/stshelp.key @@ -0,0 +1,7 @@ +a=expand b=zero base ,=left .=right w=window z=zoom +c=redraw full scale o=overplot r=redraw current scale +d=deblend e=eq. width f=functions h=1 sided eqw i=write sp j=fix pix +k=gauss fit l=flambda m=mean/snr n=fnu p=>wavelth q=quit s=smooth +t=curfit u=set wave v=velocity scale x=fix line y=plot std +/=status help ?=help -=subtr fit $=wavelength/pixel +g=new spectrum #=aperture %=band (=previous )=next diff --git a/noao/onedspec/splot/stshelp.x b/noao/onedspec/splot/stshelp.x new file mode 100644 index 00000000..f34de38a --- /dev/null +++ b/noao/onedspec/splot/stshelp.x @@ -0,0 +1,34 @@ +include <error.h> + + +# STS_HELP -- Issue a help line + +procedure sts_help (line, nlines, fname, ptr) + +int line # Line to print +int nlines # Number of lines of help +char fname[ARB] # Help file +pointer ptr # Cache help + +int fd, open(), getline() + +begin + if (ptr == NULL) { + iferr (fd = open (fname, READ_ONLY, TEXT_FILE)) { + call erract (EA_WARN) + return + } + nlines = 0 + call malloc (ptr, SZ_LINE, TY_CHAR) + while (getline (fd, Memc[ptr+nlines*SZ_LINE]) != EOF) { + nlines = nlines + 1 + call realloc (ptr, (nlines+1)*SZ_LINE, TY_CHAR) + } + call close (fd) + } + + if (line >= 1 && line <= nlines) { + call putline (STDOUT, Memc[ptr+(line-1)*SZ_LINE]) + call flush (STDOUT) + } +end diff --git a/noao/onedspec/splot/sumflux.x b/noao/onedspec/splot/sumflux.x new file mode 100644 index 00000000..36f6ca3b --- /dev/null +++ b/noao/onedspec/splot/sumflux.x @@ -0,0 +1,165 @@ +# SUMFLUX -- Sum up the flux in a specified bandpass + +procedure sumflux (sh, x, y, s, n, eqx1, eqx2, eqy1, eqy2, + sum, rsum, esum, ctr) + +pointer sh +real x[n], y[n], s[n] +int n +real eqx1, eqx2, eqy1, eqy2 +real sum[2], rsum[2], esum[2], ctr[2] + +real slope, csum[2], sum2[2], rampval, scale, delta, wpc +real w1, w2 +int i, i1, i2 +bool fp_equalr() + +begin + call fixx (sh, eqx1, eqx2, eqy1, eqy2, i1, i2) + slope = (eqy2-eqy1) / (eqx2-eqx1) + + sum[1] = 0.0 + rsum[1] = 0.0 + esum[1] = 0.0 + csum[1] = 0.0 + sum2[1] = 0.0 + scale = 0.0 + + for (i=i1+1; i <= i2-1; i = i+1) + scale = max (scale, y[i]) + if (scale <= 0.) + scale = 1. + + for (i=i1+1; i <= i2-1; i = i+1) { + rampval = eqy1 + slope * (x[i] - eqx1) + sum[1] = sum[1] + y[i] + rsum[1] = rsum[1] + rampval + if (!IS_INDEF(esum[1])) { + if (fp_equalr (0., rampval/scale)) + esum[1] = INDEF + else + esum[1] = esum[1] + (1. - y[i] / rampval) + } + } + + for (i=i1+1; i <= i2-1; i = i+1) { + rampval = eqy1 + slope * (x[i] - eqx1) + delta = (y[i] - rampval) / scale + csum[1] = csum[1] + abs(delta)**1.5 * x[i] + sum2[1] = sum2[1] + abs(delta)**1.5 + } + + # end points + if (eqx1 < x[i1]) { + if (i1 > 1) + w1 = (x[i1] - eqx1) / (x[i1] - x[i1-1]) + else + w1 = (x[i1] - eqx1) / (x[i1+1] - x[i1]) + } else { + if (i1 < n) + w1 = (x[i1] - eqx1) / (x[i1+1] - x[i1]) + else + w1 = (x[i1] - eqx1) / (x[i1] - x[i1-1]) + } + if (eqx2 < x[i2]) { + if (i2 > 1) + w2 = (x[i2] - eqx2) / (x[i2] - x[i2-1]) + else + w2 = (x[i2] - eqx2) / (x[i2+1] - x[i2]) + } else { + if (i2 < n) + w2 = (x[i2] - eqx2) / (x[i2+1] - x[i2]) + else + w2 = (x[i2] - eqx2) / (x[i2] - x[i2-1]) + } + w2 = 1.0 - w2 + + sum[1] = sum[1] + w1 * y[i1] + w2 * y[i2] + rsum[1] = rsum[1] + w1 * eqy1 + w2 * eqy2 + if (!IS_INDEF(esum[1])) { + if (fp_equalr (0., eqy1/scale)|| fp_equalr (0., eqy2/scale)) + esum[1] = INDEF + else + esum[1] = esum[1] + w1 * (1. - y[i1] / eqy1) + + w2 * (1. - y[i2] / eqy2) + } + + delta = (y[i1] - eqy1) / scale + csum[1] = csum[1] + w1 * abs(delta)**1.5 * eqx1 + sum2[1] = sum2[1] + w1 * abs(delta)**1.5 + + delta = (y[i2] - eqy2) / scale + csum[1] = csum[1] + w2 * abs(delta)**1.5 * eqx2 + sum2[1] = sum2[1] + w2 * abs(delta)**1.5 + + if (sum2[1] != 0.0) + ctr[1] = csum[1] / sum2[1] + else + ctr[1] = 0.0 + + # Correct for angstroms/channel + if (i1 != i2) + wpc = abs ((x[i2] - x[i1]) / (i2 - i1)) + else if (i1 < n) + wpc = abs (x[i1+1] - x[i1]) + else + wpc = abs (x[i1-1] - x[i1]) + sum[1] = sum[1] * wpc + if (!IS_INDEF(esum[1])) + esum[1] = esum[1] * wpc + rsum[1] = rsum[1] * wpc + + # Errors (Note there are no errors in the ramp values). + if (!IS_INDEF(s[1])) { + sum[2] = 0.0 + rsum[2] = 0.0 + esum[2] = 0.0 + csum[2] = 0.0 + sum2[2] = 0.0 + for (i=i1+1; i <= i2-1; i = i+1) { + rampval = eqy1 + slope * (x[i] - eqx1) + sum[2] = sum[2] + s[i]**2 + if (!IS_INDEF(esum[1])) { + if (fp_equalr (0., rampval/scale)) + esum[2] = INDEF + else + esum[2] = esum[2] + (s[i] / rampval) ** 2 + } + } + + for (i=i1+1; i <= i2-1; i = i+1) { + rampval = eqy1 + slope * (x[i] - eqx1) + delta = (y[i] - rampval) / scale + csum[2] = csum[2] + abs(delta)*((x[i]-ctr[1])*s[i]) ** 2 + } + + # endpoints + sum[2] = sum[2] + (w1 * s[i1])**2 + (w2 * s[i2])**2 + if (!IS_INDEF(esum[1])) { + if (fp_equalr (0., eqy1/scale)|| fp_equalr (0., eqy2/scale)) + esum[2] = INDEF + else + esum[2] = esum[2] + (w1 * s[i1] / eqy1) ** 2 + + (w2 * s[i2] / eqy2) ** 2 + } + + delta = (y[i1] - eqy1) / scale + csum[2] = csum[2] + abs(delta)*(w1*(eqx1-ctr[1])*s[i1]) ** 2 + + delta = (y[i2] - eqy2) / scale + csum[2] = csum[2] + abs(delta)*(w2*(eqx2-ctr[1])*s[i2]) ** 2 + + if (sum2[1] != 0.0) + ctr[2] = 1.5 / scale * sqrt (csum[2]) / sum2[1] + else + ctr[2] = 0.0 + + sum[2] = sqrt (sum[2]) + esum[2] = sqrt (esum[2]) + + # Correct for angstroms/channel + sum[2] = sum[2] * wpc + if (!IS_INDEF(esum[1])) + esum[2] = esum[2] * wpc + } +end diff --git a/noao/onedspec/splot/usercoord.x b/noao/onedspec/splot/usercoord.x new file mode 100644 index 00000000..2a9b3584 --- /dev/null +++ b/noao/onedspec/splot/usercoord.x @@ -0,0 +1,94 @@ +include <error.h> +include <smw.h> +include <units.h> + +# USERCOORD -- Set user coordinates + +procedure usercoord (sh, key, w1, u1, w2, u2) + +pointer sh +int key +double w1, u1, w2, u2 + +int i, format, ap, beam, dtype, nw +double shift, wa, wb, ua, ub, w0, dw, z, smw_c1trand() +real aplow[2], aphigh[2] +pointer coeff, smw, mw, ct, smw_sctran() +errchk smw_sctran + +begin + coeff = NULL + smw = MW(sh) + mw = SMW_MW(smw,0) + format = SMW_FORMAT(smw) + + iferr { + call un_ctrand (UN(sh), MWUN(sh), w1, wa, 1) + call un_ctrand (UN(sh), MWUN(sh), u1, ua, 1) + + call smw_gwattrs (MW(sh), APINDEX(sh), LINDEX(sh,2), + ap, beam, dtype, w0, dw, nw, z, aplow, aphigh, coeff) + + switch (key) { + case 'd': + wa = wa * (1 + z) + switch (UN_CLASS(MWUN(sh))) { + case UN_WAVE: + z = (wa - ua) / ua + case UN_FREQ, UN_ENERGY: + z = (ua - wa) / wa + default: + call error (1, "Inappropriate coordinate units") + } + case 'z': + shift = ua - wa + w0 = w0 + shift + if (dtype == 2) + call sshift1 (shift, coeff) + case 'l': + call un_ctrand (UN(sh), MWUN(sh), w2, wb, 1) + call un_ctrand (UN(sh), MWUN(sh), u2, ub, 1) + + switch (format) { + case SMW_ND: + i = 2 ** (SMW_PAXIS(smw,1) - 1) + ct = smw_sctran (smw, "world", "physical", i) + wa = smw_c1trand (ct, wa) + wb = smw_c1trand (ct, wb) + case SMW_ES, SMW_MS: + ct = smw_sctran (smw, "world", "physical", 3) + call smw_c2trand (ct, wa, double (ap), wa, shift) + call smw_c2trand (ct, wb, double (ap), wb, shift) + } + call smw_ctfree (ct) + + dw = (ub - ua) / (wb - wa) + w0 = ua - (wa - 1) * dw + dtype = 0 + if (UNITS(sh) == EOS) { + call mw_swattrs (mw, SMW_PAXIS(smw,1), + "label", "Wavelength") + call mw_swattrs (mw, SMW_PAXIS(smw,1), + "units", "angstroms") + } + default: + call error (1, "Unknown correction") + } + + call smw_swattrs (smw, LINDEX(sh,1), 1, ap, beam, dtype, w0, + dw, nw, z, aplow, aphigh, Memc[coeff]) + if (smw != MW(sh)) { + CTLW1(sh) = NULL + CTWL1(sh) = NULL + MW(sh) = smw + } + + DC(sh) = dtype + call shdr_system (sh, "world") + if (UN_CLASS(UN(sh)) == UN_UNKNOWN) + call un_copy (MWUN(sh), UN(sh)) + } then + call erract (EA_WARN) + + call mfree (coeff, TY_CHAR) +end diff --git a/noao/onedspec/splot/voigt.x b/noao/onedspec/splot/voigt.x new file mode 100644 index 00000000..08a44c78 --- /dev/null +++ b/noao/onedspec/splot/voigt.x @@ -0,0 +1,71 @@ +# VOIGT -- Compute the real (Voigt function) and imaginary parts of the +# complex function w(z)=exp(-z**2)*erfc(-i*z) in the upper half-plane +# z=x+iy. The maximum relative error of the real part is 2E-6 and the +# imaginary part is 5E-6. +# +# From: Humlicek, J. Quant. Spectrosc. Radiat. Transfer, V21, p309, 1979. + +procedure voigt (xarg, yarg, wr, wi) + +real xarg #I Real part of argument +real yarg #I Imaginary part of argument +real wr #O Real part of function +real wi #O Imaginary part of function + +int i +real x, y, y1, y2, y3, d, d1, d2, d3, d4, r, r2 +real t[6], c[6], s[6] + +data t/.314240376,.947788391,1.59768264,2.27950708,3.02063703,3.8897249/ +data c/1.01172805,-.75197147,1.2557727e-2,1.00220082e-2,-2.42068135e-4, + 5.00848061e-7/ +data s/1.393237,.231152406,-.155351466,6.21836624e-3,9.19082986e-5, + -6.27525958e-7/ + +begin + x = xarg + y = abs (yarg) + wr = 0. + wi = 0. + y1 = y + 1.5 + y2 = y1 * y1 + + # Region II + if (y < 0.85 && abs(x) > 18.1*y+1.65) { + if (abs(x) < 12) + wr = exp (-x * x) + y3 = y + 3 + do i = 1, 6 { + r = x - t[i] + r2 = r * r + d = 1 / (r2 + y2) + d1 = y1 * d + d2 = r * d + wr = wr + y * (c[i] * (r * d2 - 1.5 * d1) + s[i] * y3 * d2) / + (r2 + 2.25) + r = x + t[i] + r2 = r * r + d = 1 / (r2 + y2) + d3 = y1 * d + d4 = r * d + wr = wr + y * (c[i] * (r * d4 - 1.5 * d3) - s[i] * y3 * d4) / + (r2 + 2.25) + wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3) + } + + # Region I + } else { + do i = 1, 6 { + r = x - t[i] + d = 1 / (r * r + y2) + d1 = y1 * d + d2 = r * d + r = x + t[i] + d = 1 / (r * r + y2) + d3 = y1 * d + d4 = r * d + wr = wr + c[i] * (d1 + d3) - s[i] * (d2 - d4) + wi = wi + c[i] * (d2 + d4) + s[i] * (d1 - d3) + } + } +end diff --git a/noao/onedspec/splot/wrspect.x b/noao/onedspec/splot/wrspect.x new file mode 100644 index 00000000..b744a180 --- /dev/null +++ b/noao/onedspec/splot/wrspect.x @@ -0,0 +1,397 @@ +include <error.h> +include <syserr.h> +include <imhdr.h> +include <imio.h> +include <smw.h> +include <units.h> + +# SP_WRSPECT -- Write spectrum to the same image or another image. + +procedure sp_wrspect (sh1) + +pointer sh1 # Spectrum pointer to be written + +bool overwrite +pointer sp, str +int nowhite(), errcode() +bool clgetb(), xt_imnameeq() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Initially set overwrite to false in order to warn the user. + overwrite = false + + # Get new image name. + call clgstr ("new_image", Memc[str], SZ_LINE) + if (nowhite (Memc[str], Memc[str], SZ_LINE) == 0) { + call sfree (sp) + return + } + + # Check for overwriting the current file. + if (xt_imnameeq (IMNAME(sh1), Memc[str])) { + overwrite = clgetb ("overwrite") + if (!overwrite) { + call sfree (sp) + return + } + } + + # Write spectrum. + iferr (call wrspect (sh1, Memc[str], overwrite)) { + switch (errcode()) { + case SYS_IKICLOB: + call erract (EA_WARN) + # Try again if overwrite is requested. + if (!overwrite) + overwrite = clgetb ("overwrite") + if (overwrite) { + iferr (call wrspect (sh1, Memc[str], overwrite)) + call erract (EA_WARN) + } + default: + call erract (EA_WARN) + } + } + call sfree (sp) +end + + +# WRSPECT -- Write spectrum to the same image or another image. +# +# If overwriting reopen the image READ_WRITE. If this is not possible it is +# an error which may be trapped by the calling routine if desired. +# +# If writing to another image determine if the image exists. If not make a +# NEW_COPY of the image and copy all spectra and associated data. NDSPEC +# format spectra, i.e. 2D or 3D images, are copied to a 1D spectrum. +# +# If the image exists check the overwrite parameter. If overwriting, open the +# image READ_WRITE and return an error if this is not possible. If the +# output image has only one spectrum delete the image and create a NEW_COPY +# of the current spectrum image. Otherwise we will be replacing only the +# current spectrum so copy all spectra from the current image. +# +# When the input and output images are not the same open the output WCS and +# select the spectrum of the same aperture to replace. It is an error if the +# output spectrum does not contain a spectrum of the same aperture. It is +# also an error if the output spectrum is an NDSPEC image. + +procedure wrspect (sh1, output, overwrite) + +pointer sh1 # Spectrum pointer to be written +char output[ARB] # Output spectrum filename +bool overwrite # Overwrite existing spectrum? + +bool delim +char errstr[SZ_LINE] +int i, j, np1, np2, dtype[2], nw[2], err +real r[2] +double w1[2], dw[2], z[2] +pointer coeff, im, in, out, mw1, mw2, sh2, outbuf, ptr + +int imaccf(), errget() +bool xt_imnameeq(), fp_equald() +pointer immap(), smw_openim(), imgl3r(), impl3r(), imps3r() +errchk immap, imgl3r, impl3r, imps3r, imdelf, shdr_open, wrspect1 +errchk smw_openim, smw_gwattrs, smw_swattrs, smw_saveim + +begin + in = IM(sh1) + mw1 = MW(sh1) + out = NULL + mw2 = NULL + sh2 = NULL + ptr = NULL + delim = false + + iferr { + # Open and initialize the output image. + if (xt_imnameeq (IMNAME(sh1), output)) { + if (!overwrite) { + call sprintf (errstr, SZ_LINE, "No overwrite set (%s)") + call pargstr (output) + call error (1, errstr) + } + + call imunmap (in) + iferr (im = immap (IMNAME(sh1), READ_WRITE, 0)) { + in = immap (IMNAME(sh1), READ_ONLY, 0) + call erract (EA_ERROR) + } + in = im + IM(sh1) = in + out = in + mw2 = MW(sh1) + sh2 = sh1 + + } else { + iferr (im = immap (output, NEW_COPY, in)) { + if (!overwrite) + call erract (EA_ERROR) + im = immap (output, READ_WRITE, 0); out = im + + if (IM_LEN(out,2) == 1) { + call imunmap (out) + call imdelete (output) + im = immap (output, NEW_COPY, in); out = im + if (IM_PIXTYPE(out) != TY_DOUBLE) + IM_PIXTYPE(out) = TY_REAL + do j = 1, IM_LEN(out,3) + do i = 1, IM_LEN(out,2) + call amovr (Memr[imgl3r(in,i,j)], + Memr[impl3r(out,i,j)], IM_LEN(out,1)) + } + + im = smw_openim (out); mw2 = im + switch (SMW_FORMAT(mw1)) { + case SMW_ND: + if (SMW_FORMAT(mw2) != SMW_ND) + call error (1, "Incompatible spectral formats") + if (IM_NDIM(in) != IM_NDIM(out)) + call error (2, "Incompatible dimensions") + do i = 1, IM_NDIM(in) + if (IM_LEN(in,i) != IM_LEN(out,i)) + call error (2, "Incompatible dimensions") + coeff = NULL + call smw_gwattrs (mw1, 1, 1, i, i, + dtype[1], w1[1], dw[1], nw[1], z, r, r, coeff) + call smw_gwattrs (mw2, 1, 1, i, i, + dtype[2], w1[2], dw[2], nw[2], z, r, r, coeff) + call mfree (coeff, TY_CHAR) + if (dtype[1]!=dtype[2] || !fp_equald (w1[1],w1[2]) || + !fp_equald (dw[1],dw[2])) + call error (3, + "Incompatible dispersion coordinates") + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + case SMW_ES, SMW_MS: + if (SMW_FORMAT(mw2) == SMW_ND) + call error (1, "Incompatible spectral formats") + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + } + + } else { + delim = true + out = im + IM_PIXTYPE(out) = TY_REAL + im = smw_openim (out); mw2 = im + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + + do j = 1, IM_LEN(out,3) + do i = 1, IM_LEN(out,2) + call amovr (Memr[imgl3r(in,i,j)], + Memr[impl3r(out,i,j)], IM_LEN(out,1)) + } + } + + # Check, set, and update the WCS information. Note that + # wrspect1 may change the smw pointers. + + call wrspect1 (sh1, sh2) + mw1 = MW(sh1) + mw2 = MW(sh2) + call smw_saveim (mw2, out) + + # Update spectrum calibration parameters. + if (EC(sh1) == ECYES) + call imaddi (out, "EX-FLAG", EC(sh1)) + else if (imaccf (out, "EX-FLAG") == YES) + call imdelf (out, "EX-FLAG") + if (FC(sh1) == FCYES) + call imaddi (out, "CA-FLAG", FC(sh1)) + else if (imaccf (out, "CA-FLAG") == YES) + call imdelf (out, "CA-FLAG") + if (RC(sh1) != EOS) + call imastr (out, "DEREDDEN", RC(sh1)) + else if (imaccf (out, "DEREDDEN") == YES) + call imdelf (out, "DEREDDEN") + + # Copy the spectrum. + i = max (1, LINDEX(sh2,1)) + j = max (1, LINDEX(sh2,2)) + np1 = NP1(sh1) + np2 = NP2(sh1) + switch (SMW_FORMAT(mw1)) { + case SMW_ND: + switch (SMW_LAXIS(mw1,1)) { + case 1: + outbuf = imps3r (out, np1, np2, i, i, j, j) + case 2: + outbuf = imps3r (out, i, i, np1, np2, j, j) + case 3: + outbuf = imps3r (out, i, i, j, j, np1, np2) + } + call amovr (Memr[SY(sh1)], Memr[outbuf], SN(sh1)) + case SMW_ES, SMW_MS: + outbuf = impl3r (out, i, j) + call amovr (Memr[SY(sh1)], Memr[outbuf+np1-1], SN(sh1)) + if (np1 > 1) + call amovkr (Memr[outbuf+np1-1], Memr[outbuf], np1-1) + if (np2 < IM_LEN(out,1)) + call amovkr (Memr[outbuf+np2-1], Memr[outbuf+np2], + IM_LEN(out,1)-np2) + } + + # Close output image if not the same as the input image. + if (out != in) { + call shdr_close (sh2) + call smw_close (mw2) + call imunmap (out) + } + } then { + err = errget (errstr, SZ_LINE) + if (out != in) { + if (sh2 != NULL) + call shdr_close (sh2) + if (mw2 != NULL) + call smw_close (mw2) + if (out != NULL) { + call imunmap (out) + if (delim) + iferr (call imdelete (output)) + ; + } + } + call error (err, errstr) + } + +end + + +# WRSPECT1 -- Set output WCS attributes. +# This requires checking compatibility of the WCS with other spectra +# in the image. + +procedure wrspect1 (sh1, sh2) + +pointer sh1 # Input +pointer sh2 # Output + +int i, j, beam, dtype, nw +double w1, wb, dw, z, a, b, p1, p2, p3, shdr_lw() +real aplow[2], aphigh[2] +pointer in, out, smw1, smw2, mw, smw_sctran() +pointer sp, key, str, ltm, ltv, coeff +bool fp_equald(), strne() +errchk mw_glterm, smw_gwattrs, smw_swattrs, smw_sctran + +begin + in = IM(sh1) + out = IM(sh2) + smw1 = MW(sh1) + smw2 = MW(sh2) + mw = SMW_MW(smw2,0) + + # The output format must not be NDSPEC and there must be a + # matching aperture in the output image. + + if (AP(sh2) != AP(sh1) || LINDEX(sh2,1) != LINDEX(sh1,1)) + call error (6, "Matching aperture not found in output image") + + call smark (sp) + call salloc (key, SZ_LINE, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (ltm, 3*3, TY_DOUBLE) + call salloc (ltv, 3, TY_DOUBLE) + call malloc (coeff, SZ_LINE, TY_CHAR) + + # Check dispersion function compatibility. + # Nonlinear functions can't be copied to a different physical + # coordinate system though the linear dispersion can be + # adjusted. + + i = SMW_PDIM(smw2) + j = SMW_PAXIS(smw2,1) + call mw_gltermd (mw, Memd[ltm], Memd[ltv], SMW_PDIM(smw2)) + a = Memd[ltv+(j-1)] + b = Memd[ltm+(i+1)*(j-1)] + if (DC(sh1) == DCFUNC) { + i = SMW_PDIM(smw1) + j = SMW_PAXIS(smw1,1) + call mw_gltermd (SMW_MW(smw1,0), Memd[ltm], Memd[ltv], i) + Memd[ltv] = Memd[ltv+(j-1)] + Memd[ltm] = Memd[ltm+(i+1)*(j-1)] + if (!fp_equald (a, Memd[ltv]) || !fp_equald (b ,Memd[ltm])) { + call error (7, + "Physical basis for nonlinear dispersion functions don't match") + } + } + + call smw_gwattrs (smw1, LINDEX(sh1,1), LINDEX(sh1,2), + AP(sh1), beam, dtype, w1, dw, nw, z, aplow, aphigh, coeff) + + w1 = shdr_lw (sh1, 1D0) + wb = shdr_lw (sh1, double(SN(sh1))) + iferr { + call un_ctrand (UN(sh1), MWUN(sh1), w1, w1, 1) + call un_ctrand (UN(sh1), MWUN(sh1), wb, wb, 1) + } then + ; + + p1 = (NP1(sh1) - a) / b + p2 = (NP2(sh1) - a) / b + p3 = (IM_LEN(out,1) - a) / b + nw = nint (min (max (p1 ,p3), max (p1, p2))) + NP1(sh1) - 1 + if (dtype == DCLOG) { + if (p1 != p2) + dw = (log10(wb*(1+z)) - log10(w1*(1+z))) / (p2 - p1) + w1 = log10 (w1*(1+z)) - (p1 - 1) * dw + w1 = 10. ** w1 + dw = (w1 * 10D0 ** ((nw-1)*dw) - w1) / (nw - 1) + } else { + if (p1 != p2) + dw = (wb - w1) / (p2 - p1) * (1 + z) + w1 = w1 * (1 + z) - (p1 - 1) * dw + } + + # Note that this may change the smw pointer. + call smw_swattrs (smw2, LINDEX(sh2,1), 1, AP(sh2), beam, dtype, + w1, dw, nw, z, aplow, aphigh, Memc[coeff]) + if (smw2 != MW(sh2)) { + switch (SMW_FORMAT(smw2)) { + case SMW_ND, SMW_ES: + i = 2 ** (SMW_PAXIS(smw2,1) - 1) + case SMW_MS: + i = 3B + } + CTLW1(sh2) = smw_sctran (smw2, "logical", "world", i) + CTWL1(sh2) = smw_sctran (smw2, "world", "logical", i) + CTLW(sh2) = CTLW1(sh2) + CTWL(sh2) = CTWL1(sh2) + MW(sh2) = smw2 + mw = SMW_MW(smw2,0) + } + + # Copy title + call smw_sapid (smw2, LINDEX(sh2,1), 1, TITLE(sh1)) + if (Memc[SID(sh1,1)] != EOS) { + call sprintf (Memc[key], SZ_LINE, "BANDID%d") + call pargi (LINDEX(sh1,2)) + iferr (call imgstr (out, Memc[key], Memc[str], SZ_LINE)) + call imastr (out, Memc[key], Memc[SID(sh1,1)]) + else { + if (strne (Memc[SID(sh1,1)], Memc[str])) + call eprintf ( + "Warning: Input and output types (BANDID) differ\n") + } + } + + # Copy label and units + if (UN_LABEL(MWUN(sh1)) != EOS) + call mw_swattrs (mw, 1, "label", UN_LABEL(MWUN(sh1))) + if (UN_UNITS(MWUN(sh1)) != EOS) + call mw_swattrs (mw, 1, "units", UN_UNITS(MWUN(sh1))) + if (UN_USER(UN(sh1)) != EOF) + call mw_swattrs (mw, 1, "units_display", UN_USER(UN(sh1))) + + call mfree (coeff, TY_CHAR) + call sfree (sp) +end |