diff options
Diffstat (limited to 'noao/onedspec/smw/smwsaxes.x')
-rw-r--r-- | noao/onedspec/smw/smwsaxes.x | 247 |
1 files changed, 247 insertions, 0 deletions
diff --git a/noao/onedspec/smw/smwsaxes.x b/noao/onedspec/smw/smwsaxes.x new file mode 100644 index 00000000..f5e31b63 --- /dev/null +++ b/noao/onedspec/smw/smwsaxes.x @@ -0,0 +1,247 @@ +include <imhdr.h> +include <mwset.h> +include <smw.h> + + +# SMW_SAXES -- Set axes parameters based on previously set dispersion axis. +# If the dispersion axis has been excluded for NDSPEC allow another axis to +# be chosen with a warning. For EQUISPEC and MULTISPEC require the dispersion +# to be 1 and also to be present. + +procedure smw_saxes (smw, smw1, im) + +pointer smw #I SMW pointer +pointer smw1 #I Template SMW pointer +pointer im #I Template IMIO pointer + +int i, pdim, ldim, paxis, laxis, nw, dtype, nspec +real smw_c1tranr() +double w1, dw +pointer sp, str, axno, axval, r, w, cd, mw, ct, smw_sctran() +int mw_stati(), imgeti() +bool streq(), fp_equald() +errchk smw_sctran + +begin + # If a template SMW pointer is specified just copy the axes parameters. + if (smw1 != NULL) { + call strcpy (Memc[SMW_APID(smw1)], Memc[SMW_APID(smw)], SZ_LINE) + SMW_NSPEC(smw) = SMW_NSPEC(smw1) + SMW_NBANDS(smw) = SMW_NBANDS(smw1) + SMW_TRANS(smw) = SMW_TRANS(smw1) + call amovi (SMW_PAXIS(smw1,1), SMW_PAXIS(smw,1), 3) + SMW_LDIM(smw) = SMW_LDIM(smw1) + call amovi (SMW_LAXIS(smw1,1), SMW_LAXIS(smw,1), 3) + call amovi (SMW_LLEN(smw1,1), SMW_LLEN(smw,1), 3) + call amovi (SMW_NSUM(smw1,1), SMW_NSUM(smw,1), 2) + + mw = SMW_MW(smw,0) + SMW_PDIM(smw) = mw_stati (mw, MW_NDIM) + if (SMW_PDIM(smw) > SMW_PDIM(smw1)) + do i = SMW_PDIM(smw1)+1, SMW_PDIM(smw) + SMW_PAXIS(smw,i) = i + + return + } + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (axno, 3, TY_INT) + call salloc (axval, 3, TY_INT) + call aclri (Memi[axno], 3) + + # Determine the dimensions. + mw = SMW_MW(smw,0) + pdim = mw_stati (mw, MW_NDIM) + ldim = IM_NDIM(im) + call mw_gaxmap (mw, Memi[axno], Memi[axval], pdim) + + # Set the physical dispersion axis. + switch (SMW_FORMAT(smw)) { + case SMW_ND: + call salloc (r, pdim, TY_DOUBLE) + call salloc (w, pdim, TY_DOUBLE) + call salloc (cd, pdim*pdim, TY_DOUBLE) + + # Check for a transposed or rotated 2D image. + SMW_TRANS(smw) = NO + if (pdim == 2) { + call mw_gltermd (mw, Memd[cd], Memd[w], pdim) + if (Memd[cd] == 0D0 && Memd[cd+3] == 0D0) { + Memd[cd] = Memd[cd+1] + Memd[cd+1] = 0. + Memd[cd+3] = Memd[cd+2] + Memd[cd+2] = 0. + call mw_sltermd (mw, Memd[cd], Memd[w], pdim) + paxis = SMW_PAXIS(smw,1) + if (paxis == 1) + SMW_PAXIS(smw,1) = 2 + else + SMW_PAXIS(smw,1) = 1 + SMW_TRANS(smw) = YES + } else if (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0) { + Memd[w] = 0 + Memd[w+1] = 0 + Memd[cd] = 1 + Memd[cd+1] = 0 + Memd[cd+2] = 0 + Memd[cd+3] = 1 + call mw_sltermd (mw, Memd[cd], Memd[w], pdim) + } + } + + # If the dispersion axis is of length 1 or has been excluded find + # the first longer axis and print a warning. + + paxis = SMW_PAXIS(smw,1) + i = max (1, min (pdim, paxis)) + laxis = max (1, Memi[axno+i-1]) + if (IM_LEN(im,laxis) == 1) + do laxis = 1, ldim + if (IM_LEN(im,laxis) != 1) + break + + # Determine the number of spectra. + nspec = 1 + do i = 1, ldim + if (i != laxis) + nspec = nspec * IM_LEN(im,i) + SMW_NSPEC(smw) = nspec + SMW_NBANDS(smw) = 1 + + i = paxis + do paxis = 1, pdim + if (Memi[axno+paxis-1] == laxis) + break + + if (i != paxis && nspec > 1) { + call eprintf ( + "WARNING: Dispersion axis %d not found. Using axis %d.\n") + call pargi (i) + call pargi (paxis) + } + + # Set the dispersion system. + call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + if (SMW_TRANS(smw) == YES) { + Memd[cd] = Memd[cd+1] + Memd[cd+1] = 0. + Memd[cd+3] = Memd[cd+2] + Memd[cd+2] = 0. + } + if (pdim == 2 && (Memd[cd+1] != 0D0 || Memd[cd+2] != 0D0)) { + iferr (dtype = imgeti (im, "DC-FLAG")) + dtype = DCNO + if (dtype != DCNO) { + call sfree (sp) + call error (1, + "Rotated, dispersion calibrated spectra are not allowed") + } + Memd[r] = 0 + Memd[r+1] = 0 + Memd[w] = 0 + Memd[w+1] = 0 + Memd[cd] = 1 + Memd[cd+1] = 0 + Memd[cd+2] = 0 + Memd[cd+3] = 1 + } + do i = 0, pdim-1 { + dw = Memd[cd+i*(pdim+1)] + if (dw == 0D0) + Memd[cd+i*(pdim+1)] = 1D0 + } + call mw_swtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) + + dw = Memd[cd+(paxis-1)*(pdim+1)] + w1 = Memd[w+paxis-1] - (Memd[r+paxis-1] - 1) * dw + nw = IM_LEN(im,laxis) + + i = 2 ** (paxis - 1) + ct = smw_sctran (smw, "logical", "physical", i) + nw = max (smw_c1tranr (ct, 0.5), smw_c1tranr (ct, nw+0.5)) + call smw_ctfree (ct) + + iferr (dtype = imgeti (im, "DC-FLAG")) { + iferr (call mw_gwattrs (mw,paxis,"axtype",Memc[str],SZ_LINE)) + Memc[str] = EOS + if (streq (Memc[str], "ra") || streq (Memc[str], "dec")) + dtype = DCNO + else if (fp_equald (1D0, w1) || fp_equald (1D0, dw)) + dtype = DCNO + else + dtype = DCLINEAR + } + if (dtype==DCLOG) { + if (abs(w1)>20. || abs(w1+(nw-1)*dw)>20.) + dtype = DCLINEAR + else { + w1 = 10D0 ** w1 + dw = w1 * (10D0 ** ((nw-1)*dw) - 1) / (nw - 1) + } + } + + if (dtype != DCNO) { + + + iferr (call mw_gwattrs (mw,paxis,"label",Memc[str],SZ_LINE)) { + iferr (call mw_gwattrs(mw,paxis,"units",Memc[str],SZ_LINE)) { + call mw_swattrs (mw, paxis, "units", "angstroms") + call mw_swattrs (mw, paxis, "label", "Wavelength") + } + } + } + + SMW_DTYPE(smw) = INDEFI + call smw_swattrs (smw, 1, 1, INDEFI, INDEFI, + dtype, w1, dw, nw, 0D0, INDEFR, INDEFR, "") + case SMW_ES, SMW_MS: + paxis = 1 + i = Memi[axno+1] + if (i == 0) + SMW_NSPEC(smw) = 1 + else + SMW_NSPEC(smw) = IM_LEN(im,i) + i = Memi[axno+2] + if (i == 0) + SMW_NBANDS(smw) = 1 + else + SMW_NBANDS(smw) = IM_LEN(im,i) + } + + # Check and set the physical and logical dispersion axes. + laxis = Memi[axno+paxis-1] + if (laxis == 0) { + if (Memi[axval+paxis-1] == 0) + laxis = paxis + else + call error (1, "No dispersion axis") + } + + SMW_PDIM(smw) = pdim + SMW_LDIM(smw) = ldim + SMW_PAXIS(smw,1) = paxis + SMW_LAXIS(smw,1) = laxis + SMW_LLEN(smw,1) = IM_LEN(im,laxis) + SMW_LLEN(smw,2) = 1 + SMW_LLEN(smw,3) = 1 + + # Set the spatial axes. + i = 2 + do laxis = 1, ldim { + if (laxis != SMW_LAXIS(smw,1)) { + do paxis = 1, pdim + if (Memi[axno+paxis-1] == laxis) + break + SMW_PAXIS(smw,i) = paxis + SMW_LAXIS(smw,i) = laxis + SMW_LLEN(smw,i) = IM_LEN(im,laxis) + i = i + 1 + } + } + + # Set the default title. + call smw_sapid (smw, 0, 1, IM_TITLE(im)) + + call sfree (sp) +end |