Initial commit

1 files changed, 198 insertions, 0 deletions

diff --git a/sys/mwcs/mwloadim.x b/sys/mwcs/mwloadim.x
new file mode 100644
index 00000000..231f5b9a
--- /dev/null
+++ b/sys/mwcs/mwloadim.x
@@ -0,0 +1,198 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<error.h>
+include	<syserr.h>
+include	<imhdr.h>
+include	<imio.h>
+include	"mwcs.h"
+include	"imwcs.h"
+
+# MW_LOADIM -- Load a MWCS object saved in an image header in FITS format
+# into an MWCS descriptor.  Note that the MWCS descriptor is allocated
+# if the input is NULL.  This is to allow the WCS cards to be read to
+# determine the WCS dimensionality.
+
+procedure mw_loadim (mw, im)
+
+pointer	mw			#U pointer to MWCS descriptor
+pointer	im			#I pointer to image header
+
+bool	have_wcs
+int	ndim, i, j, ea_type
+int	axno[MAX_DIM], axval[MAX_DIM]
+double	maxval
+pointer	sp, sysname, iw, ct, wp, cp, bufp, ip
+
+int	mw_allocd(), mw_refstr(), ctoi(), envgeti()
+pointer	iw_rfits(), iw_findcard(), iw_gbigfits(), mw_open()
+errchk	iw_rfits, mw_allocd, mw_newsystem, mw_swtype, iw_enterwcs, mw_saxmap
+errchk	mw_open
+string	s_physical "physical"
+define	axerr_ 91
+define	axinit_ 92
+
+begin
+	call smark (sp)
+	call salloc (sysname, SZ_FNAME, TY_CHAR)
+
+	# Read the FITS image header into an IMWCS descriptor.
+	iw = iw_rfits (mw, im, RF_REFERENCE)
+	if (mw == NULL) {
+	    ndim = max (IW_NDIM(iw), IM_NPHYSDIM(im))
+	    mw = mw_open (NULL, ndim)
+	}
+	ndim = IW_NDIM(iw)
+
+	# Initialize the MWCS descriptor from the IMWCS descriptor.
+	# Free any storage associated with the old descriptor.
+	# Start with any still allocated CTRAN descriptors.
+
+	do i = 1, MAX_CTRAN {
+	    ct = MI_CTRAN(mw,i)
+	    if (ct != NULL)
+		iferr (call mw_ctfree (ct))
+		    call erract (EA_WARN)
+	}
+
+	# Free the old string and data buffers.
+	if (MI_SBUF(mw) != NULL)
+	    call mfree (MI_SBUF(mw), TY_CHAR)
+	if (MI_DBUF(mw) != NULL)
+	    call mfree (MI_DBUF(mw), TY_DOUBLE)
+
+	# Initialize the new descriptor.
+	call aclri (Memi[mw], LEN_MWCS)
+
+	MI_MAGIC(mw) = MWCS_MAGIC
+	MI_REFIM(mw) = im
+	MI_NDIM(mw) = ndim
+	MI_LTV(mw) = mw_allocd (mw, ndim)
+	MI_LTM(mw) = mw_allocd (mw, ndim * ndim)
+
+	# Set the Lterm.  Set axes with no LTM scales to unit scales.
+	# Issue a warning by default but use "wcs_matrix_err" to allow
+	# setting other error actions.
+
+	call amovd (IW_LTV(iw,1), D(mw,MI_LTV(mw)), ndim)
+	if (iw_findcard (iw, TY_LTM, ERR, 0) != NULL) {
+	    do i = 1, ndim {
+		maxval = 0D0
+		do j = 1, ndim {
+		    D(mw,MI_LTM(mw)+(j-1)*ndim+(i-1)) = IW_LTM(iw,i,j)
+		    maxval = max (maxval, abs (IW_LTM(iw,i,j)))
+		}
+		if (maxval == 0D0) {
+		    iferr (ea_type = envgeti ("wcs_matrix_err"))
+			ea_type = EA_WARN
+		    iferr {
+			switch (ea_type) {
+			case EA_FATAL, EA_ERROR:
+			    call sprintf (Memc[sysname], SZ_FNAME, 
+				"LTM keywords for axis %d undefined")
+				call pargi (i)
+			    call error (SYS_MWMISSAX, Memc[sysname])
+			case EA_WARN:
+			    IW_LTM(iw,i,i) = 1D0
+			    D(mw,MI_LTM(mw)+(i-1)*ndim+(i-1)) = IW_LTM(iw,i,i)
+			    call sprintf (Memc[sysname], SZ_FNAME, 
+				"setting LTM%d_%d to %.4g")
+				call pargi (i)
+				call pargi (i)
+				call pargd (IW_LTM(iw,i,i))
+			    call error (SYS_MWMISSAX, Memc[sysname])
+			default:
+			    IW_LTM(iw,i,i) = 1D0
+			    D(mw,MI_LTM(mw)+(i-1)*ndim+(i-1)) = IW_LTM(iw,i,i)
+			}
+		    } then
+			call erract (ea_type)
+		}
+	    }
+	} else
+	    call mw_mkidmd (D(mw,MI_LTM(mw)), ndim)
+
+	# Set up the builtin world systems "physical" and "logical".
+	# Both are linear systems.  The physical system is a unitary
+	# transformation (since world systems are defined relative to
+	# the physical system), and the logical system has the Lterm
+	# for its linear term.  No wcs attributes other than wtype are
+	# defined.
+
+	# Create the physical system.
+	call mw_newsystem (mw, s_physical, ndim)
+	do i = 1, ndim
+	    call mw_swtype (mw, i, 1, "linear", "")
+
+	# Create the logical system.
+	call mw_newsystem (mw, "logical", ndim)
+	do i = 1, ndim
+	    call mw_swtype (mw, i, 1, "linear", "")
+
+	# Set W and CD for the logical system to point to the Lterm.
+	wp = MI_WCS(mw)
+	WCS_W(wp) = MI_LTV(mw)
+	WCS_CD(wp) = MI_LTM(mw)
+
+	# Did the image header specify a WCS?
+	have_wcs = false
+	do i = 1, IW_NCARDS(iw) {
+	    cp = IW_CARD(iw,i)
+	    switch (C_TYPE(cp)) {
+	    case TY_CTYPE, TY_CRPIX, TY_CRVAL, TY_CD, TY_CDELT:
+		have_wcs = true
+		break
+	    }
+	}
+
+	# Enter the saved WCS.  We make up a system name for now, and patch
+	# it up later once the real name has been recalled along with the
+	# attributes.
+
+	if (have_wcs) {
+	    call mw_newsystem (mw, "image", ndim)
+	    call iw_enterwcs (mw, iw, ndim)
+	    ifnoerr {
+		call mw_gwattrs (mw, 0, "system", Memc[sysname], SZ_FNAME)
+	    } then
+		WCS_SYSTEM(MI_WCS(mw)) = mw_refstr (mw, Memc[sysname])
+	}
+
+	# Restore the saved WCS axis map if any.
+	if (iw_findcard (iw, TY_WAXMAP, ERR, 0) != NULL) {
+	    bufp = iw_gbigfits (iw, TY_WAXMAP, ERR)
+
+	    ip = bufp
+	    do i = 1, ndim {
+		if (ctoi (Memc, ip, axno[i]) <= 0)
+		    goto axerr_
+		if (ctoi (Memc, ip, axval[i]) <= 0) {
+axerr_		    call eprintf ("Image %s: cannot decode WAXMAP\n")
+			call pargstr (IM_NAME(IW_IM(iw)))
+		    goto axinit_
+		}
+	    }
+
+	    call mfree (bufp, TY_CHAR)
+	    call mw_saxmap (mw, axno, axval, ndim)
+
+	} else {
+axinit_	    do i = 1, ndim {
+		MI_AXNO(mw,i) = i
+		MI_AXVAL(mw,i) = 0
+	    }
+	    MI_USEAXMAP(mw) = NO
+	    MI_NLOGDIM(mw) = ndim
+	}
+
+	# Apply the section transform, if the image was opened with an image
+	# section.  This edits the axis map restored above, if any, and must
+	# be done after restoring the original WCS axis map.
+
+	call iw_setaxmap (mw, im)
+
+	# Set the default world system.
+	call mw_sdefwcs (mw)
+
+	call iw_cfits (iw)
+	call sfree (sp)
+end