From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/onedspec/splot/wrspect.x | 397 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 397 insertions(+) create mode 100644 noao/onedspec/splot/wrspect.x (limited to 'noao/onedspec/splot/wrspect.x') diff --git a/noao/onedspec/splot/wrspect.x b/noao/onedspec/splot/wrspect.x new file mode 100644 index 00000000..b744a180 --- /dev/null +++ b/noao/onedspec/splot/wrspect.x @@ -0,0 +1,397 @@ +include +include +include +include +include +include + +# SP_WRSPECT -- Write spectrum to the same image or another image. + +procedure sp_wrspect (sh1) + +pointer sh1 # Spectrum pointer to be written + +bool overwrite +pointer sp, str +int nowhite(), errcode() +bool clgetb(), xt_imnameeq() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Initially set overwrite to false in order to warn the user. + overwrite = false + + # Get new image name. + call clgstr ("new_image", Memc[str], SZ_LINE) + if (nowhite (Memc[str], Memc[str], SZ_LINE) == 0) { + call sfree (sp) + return + } + + # Check for overwriting the current file. + if (xt_imnameeq (IMNAME(sh1), Memc[str])) { + overwrite = clgetb ("overwrite") + if (!overwrite) { + call sfree (sp) + return + } + } + + # Write spectrum. + iferr (call wrspect (sh1, Memc[str], overwrite)) { + switch (errcode()) { + case SYS_IKICLOB: + call erract (EA_WARN) + # Try again if overwrite is requested. + if (!overwrite) + overwrite = clgetb ("overwrite") + if (overwrite) { + iferr (call wrspect (sh1, Memc[str], overwrite)) + call erract (EA_WARN) + } + default: + call erract (EA_WARN) + } + } + call sfree (sp) +end + + +# WRSPECT -- Write spectrum to the same image or another image. +# +# If overwriting reopen the image READ_WRITE. If this is not possible it is +# an error which may be trapped by the calling routine if desired. +# +# If writing to another image determine if the image exists. If not make a +# NEW_COPY of the image and copy all spectra and associated data. NDSPEC +# format spectra, i.e. 2D or 3D images, are copied to a 1D spectrum. +# +# If the image exists check the overwrite parameter. If overwriting, open the +# image READ_WRITE and return an error if this is not possible. If the +# output image has only one spectrum delete the image and create a NEW_COPY +# of the current spectrum image. Otherwise we will be replacing only the +# current spectrum so copy all spectra from the current image. +# +# When the input and output images are not the same open the output WCS and +# select the spectrum of the same aperture to replace. It is an error if the +# output spectrum does not contain a spectrum of the same aperture. It is +# also an error if the output spectrum is an NDSPEC image. + +procedure wrspect (sh1, output, overwrite) + +pointer sh1 # Spectrum pointer to be written +char output[ARB] # Output spectrum filename +bool overwrite # Overwrite existing spectrum? + +bool delim +char errstr[SZ_LINE] +int i, j, np1, np2, dtype[2], nw[2], err +real r[2] +double w1[2], dw[2], z[2] +pointer coeff, im, in, out, mw1, mw2, sh2, outbuf, ptr + +int imaccf(), errget() +bool xt_imnameeq(), fp_equald() +pointer immap(), smw_openim(), imgl3r(), impl3r(), imps3r() +errchk immap, imgl3r, impl3r, imps3r, imdelf, shdr_open, wrspect1 +errchk smw_openim, smw_gwattrs, smw_swattrs, smw_saveim + +begin + in = IM(sh1) + mw1 = MW(sh1) + out = NULL + mw2 = NULL + sh2 = NULL + ptr = NULL + delim = false + + iferr { + # Open and initialize the output image. + if (xt_imnameeq (IMNAME(sh1), output)) { + if (!overwrite) { + call sprintf (errstr, SZ_LINE, "No overwrite set (%s)") + call pargstr (output) + call error (1, errstr) + } + + call imunmap (in) + iferr (im = immap (IMNAME(sh1), READ_WRITE, 0)) { + in = immap (IMNAME(sh1), READ_ONLY, 0) + call erract (EA_ERROR) + } + in = im + IM(sh1) = in + out = in + mw2 = MW(sh1) + sh2 = sh1 + + } else { + iferr (im = immap (output, NEW_COPY, in)) { + if (!overwrite) + call erract (EA_ERROR) + im = immap (output, READ_WRITE, 0); out = im + + if (IM_LEN(out,2) == 1) { + call imunmap (out) + call imdelete (output) + im = immap (output, NEW_COPY, in); out = im + if (IM_PIXTYPE(out) != TY_DOUBLE) + IM_PIXTYPE(out) = TY_REAL + do j = 1, IM_LEN(out,3) + do i = 1, IM_LEN(out,2) + call amovr (Memr[imgl3r(in,i,j)], + Memr[impl3r(out,i,j)], IM_LEN(out,1)) + } + + im = smw_openim (out); mw2 = im + switch (SMW_FORMAT(mw1)) { + case SMW_ND: + if (SMW_FORMAT(mw2) != SMW_ND) + call error (1, "Incompatible spectral formats") + if (IM_NDIM(in) != IM_NDIM(out)) + call error (2, "Incompatible dimensions") + do i = 1, IM_NDIM(in) + if (IM_LEN(in,i) != IM_LEN(out,i)) + call error (2, "Incompatible dimensions") + coeff = NULL + call smw_gwattrs (mw1, 1, 1, i, i, + dtype[1], w1[1], dw[1], nw[1], z, r, r, coeff) + call smw_gwattrs (mw2, 1, 1, i, i, + dtype[2], w1[2], dw[2], nw[2], z, r, r, coeff) + call mfree (coeff, TY_CHAR) + if (dtype[1]!=dtype[2] || !fp_equald (w1[1],w1[2]) || + !fp_equald (dw[1],dw[2])) + call error (3, + "Incompatible dispersion coordinates") + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + case SMW_ES, SMW_MS: + if (SMW_FORMAT(mw2) == SMW_ND) + call error (1, "Incompatible spectral formats") + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + } + + } else { + delim = true + out = im + IM_PIXTYPE(out) = TY_REAL + im = smw_openim (out); mw2 = im + call shdr_open (out, mw2, APINDEX(sh1), LINDEX(sh1,2), + AP(sh1), SHHDR, ptr) + sh2 = ptr + + do j = 1, IM_LEN(out,3) + do i = 1, IM_LEN(out,2) + call amovr (Memr[imgl3r(in,i,j)], + Memr[impl3r(out,i,j)], IM_LEN(out,1)) + } + } + + # Check, set, and update the WCS information. Note that + # wrspect1 may change the smw pointers. + + call wrspect1 (sh1, sh2) + mw1 = MW(sh1) + mw2 = MW(sh2) + call smw_saveim (mw2, out) + + # Update spectrum calibration parameters. + if (EC(sh1) == ECYES) + call imaddi (out, "EX-FLAG", EC(sh1)) + else if (imaccf (out, "EX-FLAG") == YES) + call imdelf (out, "EX-FLAG") + if (FC(sh1) == FCYES) + call imaddi (out, "CA-FLAG", FC(sh1)) + else if (imaccf (out, "CA-FLAG") == YES) + call imdelf (out, "CA-FLAG") + if (RC(sh1) != EOS) + call imastr (out, "DEREDDEN", RC(sh1)) + else if (imaccf (out, "DEREDDEN") == YES) + call imdelf (out, "DEREDDEN") + + # Copy the spectrum. + i = max (1, LINDEX(sh2,1)) + j = max (1, LINDEX(sh2,2)) + np1 = NP1(sh1) + np2 = NP2(sh1) + switch (SMW_FORMAT(mw1)) { + case SMW_ND: + switch (SMW_LAXIS(mw1,1)) { + case 1: + outbuf = imps3r (out, np1, np2, i, i, j, j) + case 2: + outbuf = imps3r (out, i, i, np1, np2, j, j) + case 3: + outbuf = imps3r (out, i, i, j, j, np1, np2) + } + call amovr (Memr[SY(sh1)], Memr[outbuf], SN(sh1)) + case SMW_ES, SMW_MS: + outbuf = impl3r (out, i, j) + call amovr (Memr[SY(sh1)], Memr[outbuf+np1-1], SN(sh1)) + if (np1 > 1) + call amovkr (Memr[outbuf+np1-1], Memr[outbuf], np1-1) + if (np2 < IM_LEN(out,1)) + call amovkr (Memr[outbuf+np2-1], Memr[outbuf+np2], + IM_LEN(out,1)-np2) + } + + # Close output image if not the same as the input image. + if (out != in) { + call shdr_close (sh2) + call smw_close (mw2) + call imunmap (out) + } + } then { + err = errget (errstr, SZ_LINE) + if (out != in) { + if (sh2 != NULL) + call shdr_close (sh2) + if (mw2 != NULL) + call smw_close (mw2) + if (out != NULL) { + call imunmap (out) + if (delim) + iferr (call imdelete (output)) + ; + } + } + call error (err, errstr) + } + +end + + +# WRSPECT1 -- Set output WCS attributes. +# This requires checking compatibility of the WCS with other spectra +# in the image. + +procedure wrspect1 (sh1, sh2) + +pointer sh1 # Input +pointer sh2 # Output + +int i, j, beam, dtype, nw +double w1, wb, dw, z, a, b, p1, p2, p3, shdr_lw() +real aplow[2], aphigh[2] +pointer in, out, smw1, smw2, mw, smw_sctran() +pointer sp, key, str, ltm, ltv, coeff +bool fp_equald(), strne() +errchk mw_glterm, smw_gwattrs, smw_swattrs, smw_sctran + +begin + in = IM(sh1) + out = IM(sh2) + smw1 = MW(sh1) + smw2 = MW(sh2) + mw = SMW_MW(smw2,0) + + # The output format must not be NDSPEC and there must be a + # matching aperture in the output image. + + if (AP(sh2) != AP(sh1) || LINDEX(sh2,1) != LINDEX(sh1,1)) + call error (6, "Matching aperture not found in output image") + + call smark (sp) + call salloc (key, SZ_LINE, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (ltm, 3*3, TY_DOUBLE) + call salloc (ltv, 3, TY_DOUBLE) + call malloc (coeff, SZ_LINE, TY_CHAR) + + # Check dispersion function compatibility. + # Nonlinear functions can't be copied to a different physical + # coordinate system though the linear dispersion can be + # adjusted. + + i = SMW_PDIM(smw2) + j = SMW_PAXIS(smw2,1) + call mw_gltermd (mw, Memd[ltm], Memd[ltv], SMW_PDIM(smw2)) + a = Memd[ltv+(j-1)] + b = Memd[ltm+(i+1)*(j-1)] + if (DC(sh1) == DCFUNC) { + i = SMW_PDIM(smw1) + j = SMW_PAXIS(smw1,1) + call mw_gltermd (SMW_MW(smw1,0), Memd[ltm], Memd[ltv], i) + Memd[ltv] = Memd[ltv+(j-1)] + Memd[ltm] = Memd[ltm+(i+1)*(j-1)] + if (!fp_equald (a, Memd[ltv]) || !fp_equald (b ,Memd[ltm])) { + call error (7, + "Physical basis for nonlinear dispersion functions don't match") + } + } + + call smw_gwattrs (smw1, LINDEX(sh1,1), LINDEX(sh1,2), + AP(sh1), beam, dtype, w1, dw, nw, z, aplow, aphigh, coeff) + + w1 = shdr_lw (sh1, 1D0) + wb = shdr_lw (sh1, double(SN(sh1))) + iferr { + call un_ctrand (UN(sh1), MWUN(sh1), w1, w1, 1) + call un_ctrand (UN(sh1), MWUN(sh1), wb, wb, 1) + } then + ; + + p1 = (NP1(sh1) - a) / b + p2 = (NP2(sh1) - a) / b + p3 = (IM_LEN(out,1) - a) / b + nw = nint (min (max (p1 ,p3), max (p1, p2))) + NP1(sh1) - 1 + if (dtype == DCLOG) { + if (p1 != p2) + dw = (log10(wb*(1+z)) - log10(w1*(1+z))) / (p2 - p1) + w1 = log10 (w1*(1+z)) - (p1 - 1) * dw + w1 = 10. ** w1 + dw = (w1 * 10D0 ** ((nw-1)*dw) - w1) / (nw - 1) + } else { + if (p1 != p2) + dw = (wb - w1) / (p2 - p1) * (1 + z) + w1 = w1 * (1 + z) - (p1 - 1) * dw + } + + # Note that this may change the smw pointer. + call smw_swattrs (smw2, LINDEX(sh2,1), 1, AP(sh2), beam, dtype, + w1, dw, nw, z, aplow, aphigh, Memc[coeff]) + if (smw2 != MW(sh2)) { + switch (SMW_FORMAT(smw2)) { + case SMW_ND, SMW_ES: + i = 2 ** (SMW_PAXIS(smw2,1) - 1) + case SMW_MS: + i = 3B + } + CTLW1(sh2) = smw_sctran (smw2, "logical", "world", i) + CTWL1(sh2) = smw_sctran (smw2, "world", "logical", i) + CTLW(sh2) = CTLW1(sh2) + CTWL(sh2) = CTWL1(sh2) + MW(sh2) = smw2 + mw = SMW_MW(smw2,0) + } + + # Copy title + call smw_sapid (smw2, LINDEX(sh2,1), 1, TITLE(sh1)) + if (Memc[SID(sh1,1)] != EOS) { + call sprintf (Memc[key], SZ_LINE, "BANDID%d") + call pargi (LINDEX(sh1,2)) + iferr (call imgstr (out, Memc[key], Memc[str], SZ_LINE)) + call imastr (out, Memc[key], Memc[SID(sh1,1)]) + else { + if (strne (Memc[SID(sh1,1)], Memc[str])) + call eprintf ( + "Warning: Input and output types (BANDID) differ\n") + } + } + + # Copy label and units + if (UN_LABEL(MWUN(sh1)) != EOS) + call mw_swattrs (mw, 1, "label", UN_LABEL(MWUN(sh1))) + if (UN_UNITS(MWUN(sh1)) != EOS) + call mw_swattrs (mw, 1, "units", UN_UNITS(MWUN(sh1))) + if (UN_USER(UN(sh1)) != EOF) + call mw_swattrs (mw, 1, "units_display", UN_USER(UN(sh1))) + + call mfree (coeff, TY_CHAR) + call sfree (sp) +end -- cgit