diff options
Diffstat (limited to 'noao/artdata/t_mkechelle.x')
| -rw-r--r-- | noao/artdata/t_mkechelle.x | 1248 |
1 files changed, 1248 insertions, 0 deletions
diff --git a/noao/artdata/t_mkechelle.x b/noao/artdata/t_mkechelle.x new file mode 100644 index 00000000..ca9e5e37 --- /dev/null +++ b/noao/artdata/t_mkechelle.x @@ -0,0 +1,1248 @@ +include <error.h> +include <imhdr.h> +include <math.h> + +define LEN_UA 20000 # Maximum user header +define LEN_COMMENT 70 # Maximum comment length + +define PTYPES "|extracted|gaussian|slit|" +define EXTRACTED 1 # Extracted format +define GAUSS 2 # Gaussian (pexp = 2) +define SLIT 3 # Slit (pexp = 10) + + +# T_MKECHELLE -- Make echelle spectra. +# Extracted or full two dimensional formats may be created. +# The spectrum may consist of a constant continuum, a blackbody continuum, +# and emission and absorption lines of varying widths and strengths. +# The spectral features may come from a line list or be randomly generated. +# A redshift may be applied to the spectrum. The order profiles may +# be either a gaussian or a slit with a specified FWHM. Both the +# spectral features and the profile are subsampled. The echelle format +# includes a blaze function corrected for light losses to other other +# other orders and reflected components. If a focal length is specified +# the wavelength nonlinearity is included. + +procedure t_mkechelle() + +int images # List of echelle spectra to be created +int nc # Number of columns (across dispersion) +int nl # Number of lines (along dispersion) +int norders # Number of orders +int profile # Profile type +real width # Profile width (pixels) +real scattered # Scattered light peak intensity +real xc, yc # Central pixel postion +real pixsize # Pixel size (mm) + +int mc[2] # Central order +real f[2] # Focal length (mm) +real gmm[2] # Grating grooves per mm +real blaze[2] # Blaze angle (degrees) +real t[2] # Angle from blaze angle +real wc[2] # Central wavelength +real disp[2] # Central dispersion + +real z # Redshift +real cont # Continuum at central wavelength +real temp # Blackbody temperture (Kelvin) +int lines # List of files containing lines +int nrandom # Number of spectral lines +real peak # Peak/continuum +real sigma # Sigma of lines (Angstroms) +long seed # Random number seed +real subsample # Subsampling (nxsub param) +real nsigma # Dynamic range of gaussian (dynrange param) + +int ncnew, nlnew, nonew +bool new, flag[2] +int i, j, k, k1, k2, m, m1, m2, dc +long seed1 +real mwc, mw1, mw2, dmw, x, x1, x2, dx, w, p, s, xc1, dx1 +real p1, p2, pcen, fwhm, flux, flux1 +real a[2], b[2], c[2], tb[2], cb[2], tt[2], ctb[2], t2tb[2], xmin[2], xmax[2] +real aplow[2], aphigh[2] +double w1, dw +pointer sp, image, fname, apnum, comment +pointer im, mw, waves, peaks, sigmas, buf, spec, bf1, bf2, asi, data + +long clgetl(), clktime() +int clgeti(), clgwrd(), imtopenp(), imtgetim() +int nowhite(), access(), open(), fscan(), nscan() +real clgetr(), urand(), asigrl() +real ecx2w(), ecxdx(), ecw2x(), ecw2xr, ecdelta() +pointer immap(), mw_open(), impl2r(), imgl2r() +bool clgetb() +errchk open(), ecgrating() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (apnum, SZ_FNAME, TY_CHAR) + call salloc (comment, LEN_COMMENT, TY_CHAR) + + # Get parameters. + if (clgetb ("make")) + images = imtopenp ("images") + ncnew = clgeti ("ncols") + nlnew = clgeti ("nlines") + nonew = clgeti ("norders") + profile = clgwrd ("profile", Memc[comment], LEN_COMMENT, PTYPES) + width = clgetr ("width") + scattered = clgetr ("scattered") + xc = clgetr ("xc") + yc = clgetr ("yc") + pixsize = clgetr ("pixsize") + + f[1] = clgetr ("f") + mc[1] = clgeti ("order") + gmm[1] = clgetr ("gmm") + blaze[1] = clgetr ("blaze") + t[1] = clgetr ("theta") + wc[1] = clgetr ("wavelength") + disp[1] = clgetr ("dispersion") + + f[2] = clgetr ("cf") + mc[2] = clgeti ("corder") + gmm[2] = clgetr ("cgmm") + blaze[2] = clgetr ("cblaze") + t[2] = clgetr ("ctheta") + wc[2] = clgetr ("cwavelength") + disp[2] = clgetr ("cdispersion") + + z = clgetr ("rv") + if (clgetb ("z")) + z = 1 + z + else { + z = z / 299792.5 + z = sqrt ((1 + z) / (1 - z)) + } + cont = clgetr ("continuum") + temp = clgetr ("temperature") + lines = imtopenp ("lines") + peak = clgetr ("peak") + sigma = clgetr ("sigma") + seed = clgetl ("seed") + if (IS_INDEFL(seed)) + seed1 = seed1 + clktime (long(0)) + else + seed1 = seed + subsample = 1. / clgeti ("nxsub") + nsigma = sqrt (2. * log (clgetr ("dynrange"))) + + # Substitute defaults for INDEF center parameters + if (IS_INDEF(xc)) + xc = (ncnew - 1) / 2. + if (IS_INDEF(yc)) + yc = (nlnew - 1) / 2. + + # Derive and check grating parameters. + do i = 1, 2 { + if (mc[i] == 0) { + if (IS_INDEF(wc[i]) || IS_INDEF(disp[i])) + call error (1, "Prism wavelength parameters missing") + next + } + if (!IS_INDEF(pixsize)) { + if (!IS_INDEF(f[i])) + f[i] = f[i] / pixsize + if (!IS_INDEF(disp[i])) + disp[i] = disp[i] * pixsize + } + if (i == 1) + flag[i] = true + else + flag[i] = false + iferr (call ecgrating (flag[i], f[i], gmm[i], blaze[i], t[i], mc[i], + wc[i], disp[i])) { + if (i == 1) + call eprintf ("Echelle grating: ") + else + call eprintf ("Crossdisperser grating: ") + if (!IS_INDEF(mc[i])&&!IS_INDEF(wc[i])&&!IS_INDEF(disp[i])) { + call eprintf ("Using linear dispersion\n") + call erract (EA_WARN) + flag[i] = true + } else { + call eprintf ("\n") + call erract (EA_ERROR) + } + } else + flag[i] = false + } + + # List grating parameters if desired. + if (clgetb ("list")) + call eclist (pixsize, f, gmm, blaze, t, mc, wc, disp) + + # If not making an image return. + if (!clgetb ("make")) { + call sfree (sp) + return + } + + # Loop through images. Line list files may be missing. + Memc[fname] = EOS + while (imtgetim (images, Memc[image], SZ_FNAME) != EOF) { + i = imtgetim (lines, Memc[fname], SZ_FNAME) + + # Map image and check for existing images. + ifnoerr (im = immap (Memc[image], READ_WRITE, LEN_UA)) { + call eprintf ("%s: ") + call pargstr (Memc[image]) + call flush (STDERR) + if (!clgetb ("clobber")) { + call eprintf ("Warning: Image already exists (%s)\n") + call pargstr (Memc[image]) + call imunmap (im) + next + } + new = false + + if (profile == EXTRACTED) { + nl = IM_LEN(im,1) + norders = IM_LEN(im,2) + } else { + nc = IM_LEN(im,1) + nl = IM_LEN(im,2) + } + } else { + iferr (im = immap (Memc[image], NEW_IMAGE, LEN_UA)) { + call erract (EA_WARN) + next + } + new = true + + nc = ncnew + nl = nlnew + norders = nonew + + IM_PIXTYPE(im) = TY_REAL + IM_NDIM(im) = 2 + if (profile == EXTRACTED) { + IM_LEN(im,1) = nl + IM_LEN(im,2) = norders + } else { + IM_LEN(im,1) = nc + IM_LEN(im,2) = nl + } + } + + # Set frequently used constants. + do i = 1, 2 { + if (flag[i]) { + f[i] = INDEF + a[i] = mc[i] * wc[i] + b[i] = mc[i] * disp[i] + c[i] = PI * mc[i] * disp[i] + } else { + b[i] = 1e7 / gmm[i] + a[i] = b[i] * sin (DEGTORAD(t[i])) + c[i] = b[i] * PI * cos (DEGTORAD(t[i])) + ctb[i] = cos (DEGTORAD(blaze[i])) + t2tb[i] = tan (DEGTORAD(2 * blaze[i])) + tt[i] = tan (DEGTORAD(t[i] - blaze[i])) + tb[i] = tan (DEGTORAD(2 * blaze[i] - t[i])) + cb[i] = cos (DEGTORAD(2 * blaze[i] - t[i])) + xmin[i] = -1. / max (tt[i], f[i] / yc) + xmax[i] = -1. / min (tt[i], -f[i] / yc) + } + } + + # Set orders. + m1 = max (1, mc[1] - (norders-1) / 2) + m2 = m1 + norders - 1 + mwc = mc[1] * wc[1] + mw1 = ecx2w (-yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) + mw2 = ecx2w (yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) + dmw = mw2 - mw1 + if (mc[2] == 0) { + disp[2] = disp[2] / wc[2] + dx1 = 1. + } else { + xc1 = xc - ecw2x (mc[2]*wc[1], 2, a, b, f, tb, ctb, t2tb) + dx1 = ecxdx (xc - xc1, 2, f, tb) + } + + # For 2D images adjust orders to exclude those outside image. + if (profile != EXTRACTED) { + for (; m1<mc[1]; m1=m1+1) { + w = mw1 / m1 + if (mc[2] == 0) + x = (w - wc[1]) / (w * disp[2]) + xc + else + x = ecw2x (mc[2]*w, 2, a, b, f, tb, ctb, t2tb) + xc1 + if (x < nc) + break + } + for (; m2>mc[1]; m2=m2-1) { + w = mw2 / m2 + if (mc[2] == 0) + x = (w - wc[1]) / (w * disp[2]) + xc + else + x = ecw2x (mc[2]*w, 2, a, b, f, tb, ctb, t2tb) + xc1 + if (x > 0) + break + } + norders = m2 - m1 + 1 + } + + # Setup header parameters for new images. + if (new) { + call clgstr ("header", Memc[comment], LEN_COMMENT) + iferr (call mkh_header (im, Memc[comment], true, false)) + call erract (EA_WARN) + + call clgstr ("title", IM_TITLE(im), SZ_IMTITLE) + if (profile == EXTRACTED) { + mw = mw_open (NULL, 2) + call mw_newsystem (mw, "multispec", 2) + call mw_swtype (mw, 1, 1, "multispec", "") + call mw_swtype (mw, 2, 1, "multispec", "") + if (IS_INDEF(f[1])) { + call mw_swattrs (mw, 1, "label", "Wavelength") + call mw_swattrs (mw, 1, "units", "Angstroms") + } + call smw_open(mw, NULL, im) + + do m = m1, m2 { + i = m - m1 + 1 + if (IS_INDEF(f[1])) { + w1 = mw1 / m + dw = b[1] / m + dc = 0 + } else { + w1 = 1. + dw = 1. + dc = -1 + } + w = mwc / m + if (mc[2] == 0) + x = (w - wc[1]) / (w * disp[2]) + xc + else + x = ecw2x (mc[2]*w, 2, a, b, f, tb, ctb, t2tb) + xc1 + aplow[1] = 1 + x - width + aphigh[1] = 1 + x + width + aplow[2] = INDEFR + aphigh[2] = INDEFR + call smw_swattrs (mw, i, 1, i, m, dc, w1, dw, nl, + 0D0, aplow, aphigh, "") + } + call smw_saveim (mw, im) + call smw_close (mw) + } else + call imaddi (im, "DISPAXIS", 2) + } + + # Get the line list if given or create random lines. + i = nowhite (Memc[fname], Memc[fname], SZ_FNAME) + if (access (Memc[fname], 0, 0) == YES) { + i = open (Memc[fname], READ_ONLY, TEXT_FILE) + nrandom = 0 + while (fscan (i) != EOF) { + call gargr (w) + call gargr (p) + call gargr (s) + if (nscan() < 1) + next + if (nscan() < 3) + s = sigma + if (nscan() < 2) + p = peak * urand (seed1) + if (nrandom == 0) { + j = 50 + call malloc (waves, j, TY_REAL) + call malloc (peaks, j, TY_REAL) + call malloc (sigmas, j, TY_REAL) + } else if (nrandom == j) { + j = j + 10 + call realloc (waves, j, TY_REAL) + call realloc (peaks, j, TY_REAL) + call realloc (sigmas, j, TY_REAL) + } + Memr[waves+nrandom] = w * z + Memr[peaks+nrandom] = p / z + Memr[sigmas+nrandom] = s * z + nrandom = nrandom + 1 + } + call close (i) + } else { + nrandom = clgeti ("nrandom") + call malloc (waves, nrandom, TY_REAL) + call malloc (peaks, nrandom, TY_REAL) + call malloc (sigmas, nrandom, TY_REAL) + j = max (1, mc[1] - (norders-1) / 2) + do i = 0, nrandom-1 { + w = z * (mw1 + dmw * urand (seed1)) + w = w - dmw * nint ((w - mwc) / dmw) + m = j + norders * urand (seed1) + Memr[waves+i] = w / m + Memr[peaks+i] = peak * urand (seed1) / z + Memr[sigmas+i] = sigma * z + } + if (nrandom > 0 && Memc[fname] != EOS) { + i = open (Memc[fname], NEW_FILE, TEXT_FILE) + do j = 0, nrandom-1 { + call fprintf (i, "%g %g %g\n") + call pargr (Memr[waves+j] / z) + call pargr (Memr[peaks+j] * z) + call pargr (Memr[sigmas+j] / z) + } + call close (i) + } + } + + # Find the absolute response of the gratings at the reference + # blaze peak. + + flux = 1. + w = wc[1] + m = mc[1] + do i = m - 1, 1, -1 { + x = ecw2x (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + do i = m + 1, ARB { + x = ecw2x (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + j = (a[1] + b[1]) / w + do i = j, 1, -1 { + x = ecw2xr (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + if (mc[2] != 0) { + x = ecw2x (mc[2]*w, 2, a, b, f, tb, ctb, t2tb) + p = ecdelta (x, w, 2, f, c, tt) + if (p != 0.) { + p = (sin (p) / p) ** 2 + if (p != 0.) + flux = flux / p + } + } + flux1 = flux + + # Make the 1D spectrum. + call malloc (buf, norders*nl, TY_REAL) + do m = m1, m2 { + spec = buf + (m - m1) * nl + call aclrr (Memr[spec], nl) + + # Make the lines. + do i = 0, nrandom-1 { + w = m * Memr[waves+i] + p = Memr[peaks+i] * subsample + dx = m * Memr[sigmas+i] + x1 = max (-0.499, + ecw2x (w-nsigma*dx, 1, a, b, f, tb, ctb, t2tb)+yc) + x2 = min (nl-0.501, + ecw2x (w+nsigma*dx, 1, a, b, f, tb, ctb, t2tb)+yc) + dx = -0.5 / dx ** 2 + for (x = x1; x <= x2; x = x + subsample) { + s = ecx2w (x-yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) + j = nint (x) + Memr[spec+j] = Memr[spec+j] + p * exp (dx*(s-w)**2) + } + } + + # Initialize blackbody function. + if (temp > 0.) { + w = wc[1] * 1.0e-8 + x = exp (1.4388 / (w * temp)) + p1 = w**5 * (x-1.0) + } + + # Compute blaze peak correction + flux = 1. + w = mc[1] * wc[1] / m + + do i = m - 1, 1, -1 { + x = ecw2x (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + do i = m + 1, ARB { + x = ecw2x (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + j = (a[1] + b[1]) / w + do i = j, 1, -1 { + x = ecw2xr (i*w, 1, a, b, f, tb, ctb, t2tb) + if (IS_INDEF(x)) + break + p = ecdelta (x, w, 1, f, c, tt) + flux = flux + (sin (p) / p) ** 2 + if (abs (p) > 100.) + break + } + + do i = 0, nl-1 { + w = ecx2w (i-yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) / m + + # Scale by continuum. + p = cont + if (temp > 0.) { + x2 = w * 1e-8 + x1 = exp (1.4388 / (x2 * temp)) + p = p * (p1 / x2**5 / (x1-1.0)) + } + if (p > 0.) + Memr[spec+i] = max (0., p * (1. + Memr[spec+i])) + + # Apply blaze functions and pixel size corrections. + x = 1 / flux + p = ecdelta (i-yc, w, 1, f, c, tt) + if (p != 0.) + x = x * (sin (p) / p) ** 2 + + if (mc[2] != 0) { + s = ecw2x (mc[2]*w, 2, a, b, f, tb, ctb, t2tb) + p = ecdelta (s, w, 2, f, c, tt) + if (p != 0.) + x = x * (sin (p) / p) ** 2 + } + + dx = ecxdx (i-yc, 1, f, tb) * mc[1] / m + if (mc[2] != 0) + dx = dx * ecxdx (s, 2, f, tb) + + Memr[spec+i] = Memr[spec+i] * flux1 * x * dx / dx1 + } + } + + # Write 1D spectrum or create 2D spectrum. + if (profile == EXTRACTED) { + do i = 1, norders { + spec = buf + (i - 1) * nl + if (new) + call amovr (Memr[spec], Memr[impl2r(im,i)], nl) + else + call aaddr (Memr[spec], Memr[imgl2r(im,i)], + Memr[impl2r(im,i)], nl) + } + } else { + # Make scattered light response. + if (scattered > 0.) { + call malloc (bf1, nc, TY_REAL) + call malloc (bf2, nl, TY_REAL) + if (mc[2] != 0) { + do i = 0, nc-1 { + s = i - xc1 + w = ecx2w (s, 2, a, b, f, cb, tb, t2tb, + xmin, xmax) / mc[2] + p = ecdelta (s, w, 2, f, c, tt) + if (p == 0.) + Memr[bf1+i] = scattered + else + Memr[bf1+i] = scattered * (sin (p) / p) ** 2 + } + } else + call amovkr (scattered, Memr[bf1], nc) + + s = wc[1] - disp[1] * yc + do i = 0, nl - 1 { + w = ecx2w (i-yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) / + mc[1] + p = ecdelta (i-yc, w, 1, f, c, tt) + if (p == 0.) + Memr[bf2+i] = 1. + else + Memr[bf2+i] = (sin (p) / p) ** 2 + } + } else + bf1 = NULL + + # Make the profile templates stored as an interpolation function + # with the returned center, fwhm, and flux. + + switch (profile) { + case GAUSS: + call mkprof (2., asi, pcen, fwhm, flux) + case SLIT: + call mkprof (10., asi, pcen, fwhm, flux) + } + + # Now expand the 1D spectra into 2D profiles. + dx = fwhm / width + do i = 0, nl-1 { + data = impl2r (im, i+1) + if (new) + call aclrr (Memr[data], nc) + else + call amovr (Memr[imgl2r(im,i+1)], Memr[data], nc) + + if (bf1 != NULL) + do j = 0, nc-1 + Memr[data+j] = Memr[data+j] + + Memr[bf1+j] * Memr[bf2+i] + + w = ecx2w (i-yc, 1, a, b, f, cb, tb, t2tb, xmin, xmax) + do j = 0, norders-1 { + x = w / (m1 + j) + p = Memr[buf+j*nl+i] / flux + if (mc[2] == 0) + x = (x - wc[1]) / (x * disp[2]) + xc + else + x = ecw2x (mc[2]*x, 2, a, b, f, tb, ctb, t2tb) + + xc1 + p1 = max (-0.49, x - pcen / dx) + p2 = min (nc - 0.51, x + pcen / dx) + if (p1 >= p2) + next + + k1 = p1 + 0.5 + k2 = p2 + 0.5 + x1 = (p1 - x) * dx + pcen + 1 + x2 = (min (p2, k1 + 0.5) - x) * dx + pcen + 1 + x1 = max (1., x1) + x2 = max (1., x2) + + m = data + k1 + Memr[m] = Memr[m] + p * asigrl (asi, x1, x2) + do k = k1+1, k2-1 { + x1 = x2 + x2 = x1 + dx + m = m + 1 + Memr[m] = Memr[m] + p * asigrl (asi, x1, x2) + } + x1 = x2 + x2 = (p2 - x) * dx + pcen + 1 + m = m + 1 + Memr[m] = Memr[m] + p * asigrl (asi, x1, x2) + } + } + + call asifree (asi) + } + + call mfree (buf, TY_REAL) + if (bf1 != NULL) { + call mfree (bf1, TY_REAL) + call mfree (bf2, TY_REAL) + } + call mfree (waves, TY_REAL) + call mfree (peaks, TY_REAL) + call mfree (sigmas, TY_REAL) + + # Add comment history of task parameters. + if (clgetb ("comments")) { + call strcpy ("# ", Memc[comment], LEN_COMMENT) + call cnvtime (clktime (0), Memc[comment+2], LEN_COMMENT-2) + call mkh_comment (im, Memc[comment]) + call mkh_comment (im, "begin mkechelle") + call mkh_comment1 (im, "profile", 's') + if (profile != EXTRACTED) { + call mkh_comment1 (im, "width", 'r') + call mkh_comment1 (im, "scattered", 'r') + } + call mkh_comment1 (im, "norders", 'i') + call sprintf (Memc[comment], LEN_COMMENT, "%9txc%24t%g") + call pargr (1+xc) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tyc%24t%g") + call pargr (1+yc) + call mkh_comment (im, Memc[comment]) + call mkh_comment1 (im, "pixsize", 'r') + + call sprintf (Memc[comment], LEN_COMMENT, "%9tf%24t%g") + if (IS_INDEF(pixsize) || IS_INDEF(f[1])) + call pargr (f[1]) + else + call pargr (f[1] * pixsize) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tgmm%24t%g") + call pargr (gmm[1]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tblaze%24t%g") + call pargr (blaze[1]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9ttheta%24t%g") + call pargr (t[1]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9torder%24t%d") + call pargi (mc[1]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9twavelength%24t%g") + call pargr (wc[1]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tdispersion%24t%g") + if (IS_INDEF(pixsize) || IS_INDEF(disp[1])) + call pargr (disp[1]) + else + call pargr (disp[1] / pixsize) + call mkh_comment (im, Memc[comment]) + + call sprintf (Memc[comment], LEN_COMMENT, "%9tcf%24t%g") + if (IS_INDEF(pixsize) || IS_INDEF(f[2])) + call pargr (f[2]) + else + call pargr (f[2] * pixsize) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tcgmm%24t%g") + call pargr (gmm[2]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tcblaze%24t%g") + call pargr (blaze[2]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tctheta%24t%g") + call pargr (t[2]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tcorder%24t%d") + call pargi (mc[2]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tcwavelength%24t%g") + call pargr (wc[2]) + call mkh_comment (im, Memc[comment]) + call sprintf (Memc[comment], LEN_COMMENT, "%9tcdispersion%24t%g") + if (IS_INDEF(pixsize) || IS_INDEF(disp[2])) + call pargr (disp[2]) + else + call pargr (disp[2] / pixsize) + call mkh_comment (im, Memc[comment]) + + call mkh_comment1 (im, "rv", 'r') + call mkh_comment1 (im, "z", 'b') + call mkh_comment1 (im, "continuum", 'r') + call mkh_comment1 (im, "temperature", 'r') + if (nrandom > 0) { + if (Memc[fname] != EOS) + call mkh_comment1 (im, "lines", 's') + call sprintf (Memc[comment], LEN_COMMENT, "%9tnlines%24t%d") + call pargi (nrandom) + call mkh_comment (im, Memc[comment]) + call mkh_comment1 (im, "peak", 'r') + call mkh_comment1 (im, "sigma", 'r') + call mkh_comment1 (im, "seed", 'i') + } + } + + call imunmap (im) + } + + call imtclose (images) + call imtclose (lines) + call sfree (sp) +end + + +# Definitions of INDEF parameter flags. +define F 1B +define G 2B +define B 4B +define T 10B +define M 20B +define W 40B +define D 100B + +# Combinations +define FG 3B +define FB 5B +define FT 11B +define FM 21B +define FW 41B +define GB 6B +define GT 12B +define GW 42B +define GD 102B +define BT 14B +define BM 24B +define BW 44B +define BD 104B +define TM 30B +define TW 50B +define TD 110B +define MW 60B +define MD 120B +define WD 140B + + +# ECGRATING -- Derive and check grating parameters. + +procedure ecgrating (e, f, g, b, t, m, w, d) + +bool e +real f, g, b, t, w, d, x +int m + +int i, flags +define err_ 10 + +begin + if (!IS_INDEF(f)) { + if (f <= 0.) + f = INDEF + } + if (!IS_INDEF(g)) { + if (g <= 0.) + g = INDEF + else + g = g / 1e7 + } + if (!IS_INDEF(b)) { + b = DEGTORAD (b) + if (b == 0. && t == 0.) + t = INDEF + } + if (!IS_INDEF(t)) { + t = DEGTORAD (t) + if (t > PI && !IS_INDEF(b)) + t = t - TWOPI + b + } + if (!IS_INDEFI(m) && m <= 0) + m = INDEFI + if (!IS_INDEF(w) && w <= 0.) + w = INDEF + if (!IS_INDEF(d) && d <= 0.) + d = INDEF + + flags = 0 + if (IS_INDEF(f)) + flags = flags + F + if (IS_INDEF(g)) + flags = flags + G + if (IS_INDEF(b)) + flags = flags + B + if (IS_INDEF(t)) + flags = flags + T + if (IS_INDEFI(m)) + flags = flags + M + if (IS_INDEF(w)) + flags = flags + W + if (IS_INDEF(d)) + flags = flags + D + + switch (flags) { + case 0, F, G, B, T, M, W, D: + switch (flags) { + case F: + f = cos (2 * b - t) / (g * m * d) + case G: + g = (sin (t) + sin (2 * b - t)) / (m * w) + if (g == 0.) + g = INDEF + case B: + if (t > PI) { + x = g * m * w / (2 * cos (t)) + if (abs (x) > 1.) + goto err_ + b = asin (x) + t = t - TWOPI + b + } else { + x = g * m * w - sin (t) + if (abs (x) > 1.) + goto err_ + b = (t + asin (x)) / 2 + } + case T: + x = g * m * w / (2 * sin(b)) + if (abs (x) > 1.) + goto err_ + if (e) + t = b + acos (x) + else + t = b - acos (x) + case M: + m = max (1, nint ((sin(t) + sin(2*b-t)) / (g * w))) + } + if (!IS_INDEF(g)) { + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + } + case FG: + x = (sin (t) + sin (2 * b - t)) / (m * w) + if (x == 0.) + goto err_ + g = x + f = cos (2 * b - t) / (g * m * d) + case FB: + if (t > PI) { + x = g * m * w / (2 * cos (t)) + if (abs (x) > 1.) + goto err_ + b = asin (x) + t = t - TWOPI + b + } else { + x = g * m * w - sin (t) + if (abs (x) > 1.) + goto err_ + b = (t + asin (x)) / 2 + } + f = cos (2 * b - t) / (g * m * d) + case FT: + x = g * m * w / (2 * sin (b)) + if (abs (x) > 1.) + goto err_ + if (e) + t = b + acos (x) + else + t = b - acos (x) + f = cos (2 * b - t) / (g * m * d) + case FM: + m = nint ((sin (t) + sin (2 * b - t)) / (g * w)) + f = cos (2 * b - t) / (g * m * d) + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + case FW: + w = (sin (t) + sin (2 * b - t)) / (g * m) + f = cos (2 * b - t) / (g * m * d) + case GB: + x = f * d / w + if (t > PI) { + b = atan (1 / (2 * x - tan (t))) + t = t - TWOPI + b + } else { + x = (tan (t) - x) / (1 + 2 * x * tan (t)) + b = atan (x + sqrt (1 + x * x)) + } + g = (sin (t) + sin (2 * b - t)) / (m * w) + case GT: + t = b + atan (2 * f * d / w - 1 / tan (b)) + g = (sin (t) + sin (2 * b - t)) / (m * w) + case GW: + g = cos (2 * b - t) / (f * m * d) + if (g == 0.) + g = INDEF + else + w = (sin (t) + sin (2 * b - t)) / (g * m) + case GD: + x = (sin (t) + sin (2 * b - t)) / (m * w) + if (x == 0.) + goto err_ + g = x + d = cos (2 * b - t) / (f * g * m) + case BT: + x = f * g * m * d + if (abs (x) > 1.) + goto err_ + x = acos (x) + x = g * m * w - sin (x) + if (abs (x) > 1.) + goto err_ + t = asin (x) + b = (acos (f * g * m * d) + t) / 2 + case BM: + x = f * d / w + if (t > PI) { + b = atan (1 / (2 * x - tan (t))) + t = t - TWOPI + b + } else { + x = (tan (t) - x) / (1 + 2 * x * tan (t)) + b = atan (x + sqrt (1 + x * x)) + } + m = max (1, nint ((sin(t) + sin(2*b-t)) / (g * w))) + b = (t + asin (g * m * w - sin (t))) / 2 + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + case BW: + b = (t + acos (f * g * m * d)) / 2 + w = (sin (t) + sin (2 * b - t)) / (g * m) + case BD: + if (t > PI) { + x = g * m * w / (2 * cos (t)) + if (abs (x) > 1.) + goto err_ + b = asin (x) + t = t - TWOPI + b + } else { + x = g * m * w - sin (t) + if (abs (x) > 1.) + goto err_ + b = (t + asin (x)) / 2 + } + d = cos (2 * b - t) / (f * g * m) + case TM: + x = f * d / w + x = b + 2 * atan (x - 1 / (2 * tan (b))) + i = max (1, nint ((sin(x) + sin(2*b-x)) / (g * w))) + x = g * i * w / (2 * sin (b)) + if (abs (x) > 1.) + goto err_ + if (e) + t = b + acos (x) + else + t = b - acos (x) + m = i + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + case TW: + x = f * g * m * d + if (abs (x) > 1.) + goto err_ + t = 2 * b - acos (x) + w = (sin (t) + sin (2 * b - t)) / (g * m) + case TD: + x = g * m * w / (2 * sin (b)) + if (abs (x) > 1.) + goto err_ + if (e) + t = b + acos (x) + else + t = b - acos (x) + d = cos (2 * b - t) / (f * g * m) + case MW: + m = max (1, nint (cos (2 * b - t) / (f * g * d))) + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + case MD: + m = max (1, nint ((sin(t) + sin(2*b-t)) / (g * w))) + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + case WD: + w = (sin (t) + sin (2 * b - t)) / (g * m) + d = cos (2 * b - t) / (f * g * m) + } + + if (!IS_INDEF(g)) + g = g * 1e7 + if (!IS_INDEF(b)) + b = RADTODEG (b) + if (!IS_INDEF(t)) + t = RADTODEG (t) + + if (IS_INDEF(f) || IS_INDEF(g) || IS_INDEF(b) || IS_INDEF(t) || + IS_INDEF(m) || IS_INDEF(w) || IS_INDEF(d)) + call error (1, + "Insufficient information to to resolve grating parameters") + + return + +err_ if (!IS_INDEF(g)) + g = g * 1e7 + if (!IS_INDEF(b)) + b = RADTODEG (b) + if (!IS_INDEF(t)) + t = RADTODEG (t) + call error (2, "Impossible combination of grating parameters") +end + + +# ECLIST -- List grating parameters. + +procedure eclist (p, f, g, b, t, m, w, d) + +real p, f[2], g[2], b[2], t[2], w[2], d[2] +int m[2] + +begin + call printf ("Echelle grating parameters:\n") + call printf (" Focal length = %g %s\n") + if (IS_INDEF(p) || IS_INDEF(f[1])) + call pargr (f[1]) + else + call pargr (f[1] * p) + if (IS_INDEF(p)) + call pargstr ("pixels") + else + call pargstr ("mm") + call printf (" Grating = %g grooves/mm\n") + call pargr (g[1]) + call printf (" Blaze angle = %g degrees\n") + call pargr (b[1]) + call printf (" Incidence angle = %g degrees\n") + call pargr (t[1]) + call printf (" Reference order = %d\n") + call pargi (m[1]) + call printf ( + " Blaze wavelength of reference order = %g Angstroms\n") + call pargr (w[1]) + call printf ( + " Blaze dispersion of reference order = %g Angstroms/%s\n") + if (IS_INDEF(p) || IS_INDEF(d[1])) + call pargr (d[1]) + else + call pargr (d[1] / p) + if (IS_INDEF(p)) + call pargstr ("pixels") + else + call pargstr ("mm") + + if (m[2] == 0.) { + call printf ("Crossdisperser prism parameters:\n") + call printf (" Reference wavelength = %g Angstroms/pixel\n") + call pargr (w[2]) + call printf ( + " Dispersion at reference wavelength = %g Angstroms/%s\n") + if (IS_INDEF(p) || IS_INDEF(d[2])) + call pargr (d[2]) + else + call pargr (d[2] / p) + if (IS_INDEF(p)) + call pargstr ("pixels") + else + call pargstr ("mm") + } else { + call printf ("Crossdisperser grating parameters:\n") + call printf (" Focal length = %g %s\n") + if (IS_INDEF(p) || IS_INDEF(f[2])) + call pargr (f[2]) + else + call pargr (f[2] * p) + if (IS_INDEF(p)) + call pargstr ("pixels") + else + call pargstr ("mm") + call printf (" Grating = %g grooves/mm\n") + call pargr (g[2]) + call printf (" Blaze angle = %g degrees\n") + call pargr (b[2]) + call printf (" Incidence angle = %g degrees\n") + call pargr (t[2]) + call printf (" Order = %d\n") + call pargi (m[2]) + call printf ( + " Blaze wavelength = %g Angstroms\n") + call pargr (w[2]) + call printf ( + " Blaze dispersion = %g Angstroms/%s\n") + if (IS_INDEF(p) || IS_INDEF(d[2])) + call pargr (d[2]) + else + call pargr (d[2] / p) + if (IS_INDEF(p)) + call pargstr ("pixels") + else + call pargstr ("mm") + } + call flush (STDOUT) +end + + +# ECX2W -- Given pixel position return wavelength. + +real procedure ecx2w (x, i, a, b, f, cb, tb, t2tb, xmin, xmax) + +real x, a[2], b[2], f[2], cb[2], tb[2], t2tb[2], xmin[2], xmax[2], w +int i + +begin + if (IS_INDEF(f[i])) + return (a[i] + b[i] * x) + w = x / f[1] + w = a[i] + b[i] * cb[i] / sqrt (1 + w * w) * (w + tb[i]) + return (w) +end + + +# ECX2WR -- Given pixel position return wavelength of reflected component. + +real procedure ecx2wr (x, i, a, b, f, cb, tb, t2tb, xmin, xmax) + +real x, a[2], b[2], f[2], cb[2], tb[2], t2tb[2], xmin[2], xmax[2], w +int i + +begin + if (IS_INDEF(f[i])) + return (INDEF) + + w = x / f[i] + if (x <= xmin[i] || x >= xmax[i]) + return (INDEF) + + w = (w - t2tb[i]) / (1 + w * t2tb[i]) + w = a[i] + b[i] * cb[i] / sqrt (1 + w * w) * (w + tb[i]) + return (w) +end + + +# ECXDX -- Given pixel position return pixel size per unit wavelength +# normalized to the central pixel. + +real procedure ecxdx (x, i, f, tb) + +real x, f[2], tb[2], dx +int i + +begin + if (IS_INDEF(f[i])) + return (1.) + + dx = x / f[i] + dx = (1 - dx * tb[i]) / sqrt ((1 + dx * dx) ** 3) + return (dx) +end + + +# ECW2X -- Given wavelength return pixel position. + +real procedure ecw2x (w, i, a, b, f, tb, ctb, t2tb) + +real w, a[2], b[2], f[2], tb[2], ctb[2], t2tb[2], x +int i + +begin + if (IS_INDEF(f[i])) + return ((w - a[i]) / b[i]) + + x = (w - a[i]) / b[i] + if (x >= 1. || x <= -ctb[i]) + return (INDEF) + x = x / sqrt (1 - x * x) + x = f[i] * (x - tb[i]) / (1 + x * tb[i]) + return (x) + +end + + +# ECW2XR -- Given wavelength return pixel position of reflected component. + +real procedure ecw2xr (w, i, a, b, f, tb, ctb, t2tb) + +real w, a[2], b[2], f[2], tb[2], ctb[2], t2tb[2], x +int i + +begin + if (IS_INDEF(f[i])) + return (INDEF) + + x = (w - a[i]) / b[i] + if (x >= 1. || x <= ctb[i]) + return (INDEF) + x = x / sqrt (1 - x * x) + x = (x - tb[i]) / (1 + x * tb[i]) + x = f[i] * (x + t2tb[i]) / (1 - x * t2tb[i]) + return (x) +end + + +# ECDELTA -- Given pixel position and wavelength return blaze function +# phase angle. + +real procedure ecdelta (x, w, i, f, c, tt) + +real x, w, f[2], c[2], tt[2], d +int i + +begin + if (IS_INDEF(f[i])) + return (c[i] / w * x) + + d = x / f[i] + d = 1 / sqrt (1 + d * d) + d = c[i] / w * (d * x / f[i] + tt[i] * (1 - d)) + return (d) +end |