1 files changed, 129 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/imtab/itbwcs.x b/pkg/utilities/nttools/imtab/itbwcs.x
new file mode 100644
index 00000000..4ecea420
--- /dev/null
+++ b/pkg/utilities/nttools/imtab/itbwcs.x
@@ -0,0 +1,129 @@
+include <imhdr.h>
+include <mwset.h>
+include "imtab.h"
+
+# This file contains three routines:  itb_wcs_init and the single and
+# double precision routines itb_ctranr & itb_ctrand.
+#
+# Phil Hodge, 30-Sep-1993  Subroutines created.
+
+# itb_wcs_init -- open wcs, etc
+# This routine gets the wcs, turns axis mapping off, and initializes
+# the transformation for physical or world coordinates.  The dimension
+# of the original image and the mapping from logical to physical axis 
+# numbers are returned for use by itb_ctrand and itb_ctranr.
+
+procedure itb_wcs_init (im, wcs_type, mw, ct, ax, wcsdim)
+
+pointer im			# i: imhdr pointer for image
+int	wcs_type		# i: wcs type
+pointer mw, ct			# o: mwcs pointers
+int	ax[IM_MAXDIM]		# o: ax[i] is physical axis for logical axis i
+int	wcsdim			# o: dimension of physical image coord system
+#--
+int	axno[IM_MAXDIM]		# axno[j] is logical axis for physical axis j
+int	axval[IM_MAXDIM]	# axval[j] is value if axno[j] is zero
+int	ndim			# number of "logical" axes
+int	i, j
+pointer	mw_openim(), mw_sctran()
+int	mw_stati()
+
+begin
+	if (wcs_type == IMTAB_NO_WCS || wcs_type == IMTAB_LOGICAL) {
+	    mw = NULL
+	    ct = NULL
+	    wcsdim = IM_NDIM(im)
+	    return
+	}
+
+	# Get the wcs.
+	mw = mw_openim (im)
+
+	# Set up the transformation.
+	call mw_seti (mw, MW_USEAXMAP, NO)
+	if (wcs_type == IMTAB_PHYSICAL)
+	    ct = mw_sctran (mw, "logical", "physical", 0)
+	else if (wcs_type == IMTAB_WORLD)
+	    ct = mw_sctran (mw, "logical", "world", 0)
+	wcsdim = mw_stati (mw, MW_NPHYSDIM)
+	ndim = IM_NDIM(im)
+
+	# Get the logical axis number corresponding to each physical axis.
+	call mw_gaxmap (mw, axno, axval, wcsdim)
+
+	# Invert axno:  get the physical axis number for each logical axis.
+	do i = 1, ndim				# initialize
+	    ax[i] = 0
+	do j = 1, wcsdim {
+	    do i = 1, ndim {
+		if (axno[j] == i) {
+		    ax[i] = j
+		    break
+		}
+	    }
+	}
+
+	# Check to be sure each axis was found.
+	do i = 1, ndim {
+	    if (ax[i] < 1)
+		call error (1, "itb_mwcs_init:  an axis was not found")
+	}
+end
+
+# itb_ctran -- translate coordinates with axis mapping = NO
+# This routine translates "logical" coordinates to "physical" or "world".
+# Axis mapping must have been turned off, and the mapping from logical
+# to physical axes is given by the array AX:  if I is a logical axis
+# number, AX[I] is the corresponding physical axis number.  Each element
+# of the array IN must have been initialized to one by the calling routine.
+# Separate single and double precision versions are included.
+
+procedure itb_ctrand (im, ct, ax, in, incoords, outcoords, wcsdim)
+
+pointer im			# i: imhdr pointer
+pointer ct			# i: coordinate transformation pointer
+int	ax[wcsdim]		# i: "logical" to "physical" mapping
+double	in[IM_MAXDIM]		# io: copy of incoords but includes axis mapping
+double	incoords[wcsdim]	# i: input "logical" coordinates
+double	outcoords[wcsdim]	# o: output coordinates
+int	wcsdim			# i: length of incoords & outcoords arrays
+#--
+int	i
+
+begin
+	if (ct == NULL) {
+	    call amovd (incoords, outcoords, wcsdim)
+	    return
+	}
+
+	# Take account of axis mapping; i is the logical axis number.
+	do i = 1, IM_NDIM(im)
+	    in[ax[i]] = incoords[i]
+
+	call mw_ctrand (ct, in, outcoords, wcsdim)
+end
+
+procedure itb_ctranr (im, ct, ax, in, incoords, outcoords, wcsdim)
+
+pointer im			# i: imhdr pointer
+pointer ct			# i: coordinate transformation pointer
+int	ax[wcsdim]		# i: "logical" to "physical" mapping
+real	in[IM_MAXDIM]		# io: copy of incoords but includes axis mapping
+real	incoords[wcsdim]	# i: input "logical" coordinates
+real	outcoords[wcsdim]	# o: output coordinates
+int	wcsdim			# i: length of incoords & outcoords arrays
+#--
+int	i
+
+begin
+	if (ct == NULL) {
+	    call amovr (incoords, outcoords, wcsdim)
+	    return
+	}
+
+	# Take account of axis mapping; i is the logical axis number.
+	do i = 1, IM_NDIM(im)
+	    in[ax[i]] = incoords[i]
+
+	call mw_ctranr (ct, in, outcoords, wcsdim)
+end