Initial commit

1 files changed, 97 insertions, 0 deletions

diff --git a/sys/mwcs/gen/mwctrand.x b/sys/mwcs/gen/mwctrand.x
new file mode 100644
index 00000000..70e575cc
--- /dev/null
+++ b/sys/mwcs/gen/mwctrand.x
@@ -0,0 +1,97 @@
+include	"../mwcs.h"
+
+# MW_CTRAN -- Transform a single N-dimensional point, using the optimized
+# transformation set up by a prior call to MW_SCTRAN.
+
+procedure mw_ctrand (a_ct, p1, p2, ndim)
+
+pointer	a_ct			#I pointer to CTRAN descriptor
+double	p1[ndim]		#I coordinates of point in input system
+double	p2[ndim]		#O coordinates of point in output system
+int	ndim			#I dimensionality of point
+
+int	naxes, i, j
+pointer	ct, fc, ltm, ltv, d_ct
+double	v1[MAX_DIM], v2[MAX_DIM], iv[MAX_DIM], ov[MAX_DIM]
+errchk	zcall3
+
+begin
+	# Get real or double version of descriptor.
+	ct = CT_D(a_ct)
+
+	ltm = CT_LTM(ct)
+	ltv = CT_LTV(ct)
+
+	# Specially optimized cases.
+	if (CT_TYPE(ct) == LNR) {
+	    # Simple linear, nonrotated transformation.
+	    do i = 1, ndim
+		p2[i] = Memd[ltm+(i-1)*(ndim+1)] * p1[i] + Memd[ltv+i-1]
+	    return
+	} else if (CT_TYPE(ct) == LRO) {
+	    # Simple linear, rotated transformation.
+	    call mw_ltrand (p1, p2, Memd[ltm], Memd[ltv], ndim)
+	    return
+	}
+
+	# If we get here the transformation involves a call to one or more
+	# WCS functions.  In this general case, the transformation consists
+	# of zero or more calls to WCS functions to transform the input
+	# world coordinates to the linear input system, followed by a general
+	# linear transformation to the linear output system, followed by zero
+	# or more calls to WCS functions to do the forward transformation
+	# to generate the final output world coordinates.  The WCS function
+	# calls are always evaluated in double precision.
+
+	# Make zero or more WCS function calls for the different axes of the
+	# input system (inverse transform).
+
+	call achtdd (p1, iv, ndim)
+	do j = 1, CT_NCALLI(ct) {
+	    # Get pointer to function call descriptor.
+	    fc = CT_FCI(ct,j)
+	    naxes = FC_NAXES(fc)
+
+	    # Extract the coordinate vector for the function call.
+	    do i = 1, naxes
+		v1[i] = p1[FC_AXIS(fc,i)]
+
+	    # Call the WCS function.
+	    call zcall3 (FC_FCN(fc), fc, v1, v2)
+
+	    # Edit the vector IV, replacing the entries associated with
+	    # the WCS function by the transformed values.
+
+	    do i = 1, naxes
+		iv[FC_AXIS(fc,i)] = v2[i]
+	}
+
+	# Apply the general linear transformation.  We may as well do this in
+	# double since we already have to use double for the function calls.
+
+	d_ct = CT_D(a_ct)
+	call mw_ltrand (iv, ov, Memd[CT_LTM(d_ct)], Memd[CT_LTV(d_ct)], ndim)
+
+	# Make zero or more WCS function calls for the different axes of the
+	# output system (forward transform to final world system).
+
+	call achtdd (ov, p2, ndim)
+	do j = 1, CT_NCALLO(ct) {
+	    # Get pointer to function call descriptor.
+	    fc = CT_FCO(ct,j)
+	    naxes = FC_NAXES(fc)
+
+	    # Extract the coordinate vector for the function call.
+	    do i = 1, naxes
+		v1[i] = ov[FC_AXIS(fc,i)]
+
+	    # Call the WCS function.
+	    call zcall3 (FC_FCN(fc), fc, v1, v2)
+
+	    # Edit the final output vector, replacing the entries for the
+	    # function axes by their transformed values.
+
+	    do i = 1, naxes
+		p2[FC_AXIS(fc,i)] = v2[i]
+	}
+end