From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- noao/onedspec/splot/anshdr.x | 84 ++++ noao/onedspec/splot/autoexp.x | 79 ++++ noao/onedspec/splot/avgsnr.x | 72 ++++ noao/onedspec/splot/conflam.x | 28 ++ noao/onedspec/splot/confnu.x | 28 ++ noao/onedspec/splot/deblend.x | 627 ++++++++++++++++++++++++++++++ noao/onedspec/splot/eqwidth.x | 109 ++++++ noao/onedspec/splot/eqwidthcp.x | 240 ++++++++++++ noao/onedspec/splot/fixx.x | 27 ++ noao/onedspec/splot/flatten.x | 110 ++++++ noao/onedspec/splot/fudgept.x | 38 ++ noao/onedspec/splot/fudgex.x | 46 +++ noao/onedspec/splot/getimage.x | 159 ++++++++ noao/onedspec/splot/gfit.x | 391 +++++++++++++++++++ noao/onedspec/splot/mkpkg | 38 ++ noao/onedspec/splot/mktitle.x | 41 ++ noao/onedspec/splot/plotstd.x | 70 ++++ noao/onedspec/splot/replot.x | 27 ++ noao/onedspec/splot/smooth.x | 54 +++ noao/onedspec/splot/spdeblend.x | 819 +++++++++++++++++++++++++++++++++++++++ noao/onedspec/splot/splabel.x | 112 ++++++ noao/onedspec/splot/splot.key | 116 ++++++ noao/onedspec/splot/splot.log | 8 + noao/onedspec/splot/splot.x | 605 +++++++++++++++++++++++++++++ noao/onedspec/splot/splotcolon.x | 263 +++++++++++++ noao/onedspec/splot/splotfun.x | 127 ++++++ noao/onedspec/splot/stshelp.key | 7 + noao/onedspec/splot/stshelp.x | 34 ++ noao/onedspec/splot/sumflux.x | 165 ++++++++ noao/onedspec/splot/usercoord.x | 94 +++++ noao/onedspec/splot/voigt.x | 71 ++++ noao/onedspec/splot/wrspect.x | 397 +++++++++++++++++++ 32 files changed, 5086 insertions(+) create mode 100644 noao/onedspec/splot/anshdr.x create mode 100644 noao/onedspec/splot/autoexp.x create mode 100644 noao/onedspec/splot/avgsnr.x create mode 100644 noao/onedspec/splot/conflam.x create mode 100644 noao/onedspec/splot/confnu.x create mode 100644 noao/onedspec/splot/deblend.x create mode 100644 noao/onedspec/splot/eqwidth.x create mode 100644 noao/onedspec/splot/eqwidthcp.x create mode 100644 noao/onedspec/splot/fixx.x create mode 100644 noao/onedspec/splot/flatten.x create mode 100644 noao/onedspec/splot/fudgept.x create mode 100644 noao/onedspec/splot/fudgex.x create mode 100644 noao/onedspec/splot/getimage.x create mode 100644 noao/onedspec/splot/gfit.x create mode 100644 noao/onedspec/splot/mkpkg create mode 100644 noao/onedspec/splot/mktitle.x create mode 100644 noao/onedspec/splot/plotstd.x create mode 100644 noao/onedspec/splot/replot.x create mode 100644 noao/onedspec/splot/smooth.x create mode 100644 noao/onedspec/splot/spdeblend.x create mode 100644 noao/onedspec/splot/splabel.x create mode 100644 noao/onedspec/splot/splot.key create mode 100644 noao/onedspec/splot/splot.log create mode 100644 noao/onedspec/splot/splot.x create mode 100644 noao/onedspec/splot/splotcolon.x create mode 100644 noao/onedspec/splot/splotfun.x create mode 100644 noao/onedspec/splot/stshelp.key create mode 100644 noao/onedspec/splot/stshelp.x create mode 100644 noao/onedspec/splot/sumflux.x create mode 100644 noao/onedspec/splot/usercoord.x create mode 100644 noao/onedspec/splot/voigt.x create mode 100644 noao/onedspec/splot/wrspect.x (limited to 'noao/onedspec/splot') 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 +include +include + +# 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 +include +include + +# 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 +include + +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 +include + +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 +include + +# 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 + + +# 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 + +# 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 + +# 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 +include +include +include + +# 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 +include +include + +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 + autoexp.x + avgsnr.x + conflam.x + confnu.x + deblend.x + eqwidth.x + eqwidthcp.x + fixx.x + flatten.x + fudgept.x + fudgex.x + getimage.x + gfit.x + mktitle.x + plotstd.x + replot.x + smooth.x + spdeblend.x + splabel.x + splot.x \ + + splotcolon.x + splotfun.x + stshelp.x + sumflux.x + usercoord.x + voigt.x + wrspect.x + ; 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 +include +include + +# 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 +include +include + +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 + +# 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 +include + +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 +include + +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 - 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. + +:# - Add comment to log file +:dispaxis - Change summing parameter for 2D images +:log - Enable logging to save_file +:nolog - Disable logging to save_file +:nsum - Change summing parameter for 2D images +:show - Show full output of deblending and equiv. width measurments +:units - Change coordinate units (see below) + +:label