aboutsummaryrefslogtreecommitdiff
path: root/sys/mwcs/gen
diff options
context:
space:
mode:
Diffstat (limited to 'sys/mwcs/gen')
-rw-r--r--sys/mwcs/gen/mkpkg29
-rw-r--r--sys/mwcs/gen/mwc1trand.x24
-rw-r--r--sys/mwcs/gen/mwc1tranr.x24
-rw-r--r--sys/mwcs/gen/mwc2trand.x38
-rw-r--r--sys/mwcs/gen/mwc2tranr.x38
-rw-r--r--sys/mwcs/gen/mwctrand.x97
-rw-r--r--sys/mwcs/gen/mwctranr.x97
-rw-r--r--sys/mwcs/gen/mwgctrand.x44
-rw-r--r--sys/mwcs/gen/mwgctranr.x44
-rw-r--r--sys/mwcs/gen/mwltrand.x26
-rw-r--r--sys/mwcs/gen/mwltranr.x26
-rw-r--r--sys/mwcs/gen/mwmmuld.x21
-rw-r--r--sys/mwcs/gen/mwmmulr.x21
-rw-r--r--sys/mwcs/gen/mwv1trand.x32
-rw-r--r--sys/mwcs/gen/mwv1tranr.x32
-rw-r--r--sys/mwcs/gen/mwv2trand.x49
-rw-r--r--sys/mwcs/gen/mwv2tranr.x49
-rw-r--r--sys/mwcs/gen/mwvmuld.x20
-rw-r--r--sys/mwcs/gen/mwvmulr.x20
-rw-r--r--sys/mwcs/gen/mwvtrand.x18
-rw-r--r--sys/mwcs/gen/mwvtranr.x18
21 files changed, 767 insertions, 0 deletions
diff --git a/sys/mwcs/gen/mkpkg b/sys/mwcs/gen/mkpkg
new file mode 100644
index 00000000..bc8fe837
--- /dev/null
+++ b/sys/mwcs/gen/mkpkg
@@ -0,0 +1,29 @@
+# Make the generic portion of MWCS.
+
+$checkout libex.a lib$
+$udate libex.a
+$checkin libex.a lib$
+$exit
+
+libex.a:
+ mwc1trand.x ../mwcs.h
+ mwc1tranr.x ../mwcs.h
+ mwc2trand.x ../mwcs.h
+ mwc2tranr.x ../mwcs.h
+ mwctrand.x ../mwcs.h
+ mwctranr.x ../mwcs.h
+ mwgctrand.x ../mwcs.h
+ mwgctranr.x ../mwcs.h
+ mwltrand.x
+ mwltranr.x
+ mwmmuld.x
+ mwmmulr.x
+ mwv1trand.x ../mwcs.h
+ mwv1tranr.x ../mwcs.h
+ mwv2trand.x ../mwcs.h
+ mwv2tranr.x ../mwcs.h
+ mwvmuld.x
+ mwvmulr.x
+ mwvtrand.x
+ mwvtranr.x
+ ;
diff --git a/sys/mwcs/gen/mwc1trand.x b/sys/mwcs/gen/mwc1trand.x
new file mode 100644
index 00000000..af46e02d
--- /dev/null
+++ b/sys/mwcs/gen/mwc1trand.x
@@ -0,0 +1,24 @@
+include "../mwcs.h"
+
+# MW_C1TRAN -- Optimized 1D coordinate transformation.
+
+double procedure mw_c1trand (a_ct, x)
+
+pointer a_ct #I pointer to CTRAN descriptor
+double x #I coordinates in input system
+
+double y
+pointer ct
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_D(a_ct)
+
+ # Perform the transformation; LNR is a simple linear transformation.
+ if (CT_TYPE(ct) == LNR) {
+ return (Memd[CT_LTM(ct)] * x + Memd[CT_LTV(ct)])
+ } else {
+ call mw_ctrand (a_ct, x, y, 1)
+ return (y)
+ }
+end
diff --git a/sys/mwcs/gen/mwc1tranr.x b/sys/mwcs/gen/mwc1tranr.x
new file mode 100644
index 00000000..06ad0bf7
--- /dev/null
+++ b/sys/mwcs/gen/mwc1tranr.x
@@ -0,0 +1,24 @@
+include "../mwcs.h"
+
+# MW_C1TRAN -- Optimized 1D coordinate transformation.
+
+real procedure mw_c1tranr (a_ct, x)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real x #I coordinates in input system
+
+real y
+pointer ct
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_R(a_ct)
+
+ # Perform the transformation; LNR is a simple linear transformation.
+ if (CT_TYPE(ct) == LNR) {
+ return (Memr[CT_LTM(ct)] * x + Memr[CT_LTV(ct)])
+ } else {
+ call mw_ctranr (a_ct, x, y, 1)
+ return (y)
+ }
+end
diff --git a/sys/mwcs/gen/mwc2trand.x b/sys/mwcs/gen/mwc2trand.x
new file mode 100644
index 00000000..0cb156bd
--- /dev/null
+++ b/sys/mwcs/gen/mwc2trand.x
@@ -0,0 +1,38 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_C2TRAN -- Optimized 2D coordinate transformation.
+
+procedure mw_c2trand (a_ct, x1,y1, x2,y2)
+
+pointer a_ct #I pointer to CTRAN descriptor
+double x1,y1 #I coordinates in input system
+double x2,y2 #O coordinates in output system
+
+pointer ct, ltm, ltv
+double p1[2], p2[2]
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_D(a_ct)
+
+ ltm = CT_LTM(ct)
+ ltv = CT_LTV(ct)
+
+ if (CT_TYPE(ct) == LNR) {
+ # Simple linear, nonrotated transformation.
+ x2 = Memd[ltm ] * x1 + Memd[ltv ]
+ y2 = Memd[ltm+3] * y1 + Memd[ltv+1]
+ } else if (CT_TYPE(ct) == LRO) {
+ # Linear, rotated transformation.
+ p1[1] = x1; p1[2] = y1
+ x2 = Memd[ltm ] * p1[1] + Memd[ltm+1] * p1[2] + Memd[ltv ]
+ y2 = Memd[ltm+2] * p1[1] + Memd[ltm+3] * p1[2] + Memd[ltv+1]
+ } else {
+ # General case involving one or more functional terms.
+ p1[1] = x1; p1[2] = y1
+ call mw_ctrand (a_ct, p1, p2, 2)
+ x2 = p2[1]; y2 = p2[2]
+ }
+end
diff --git a/sys/mwcs/gen/mwc2tranr.x b/sys/mwcs/gen/mwc2tranr.x
new file mode 100644
index 00000000..ef5b5ef7
--- /dev/null
+++ b/sys/mwcs/gen/mwc2tranr.x
@@ -0,0 +1,38 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_C2TRAN -- Optimized 2D coordinate transformation.
+
+procedure mw_c2tranr (a_ct, x1,y1, x2,y2)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real x1,y1 #I coordinates in input system
+real x2,y2 #O coordinates in output system
+
+pointer ct, ltm, ltv
+real p1[2], p2[2]
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_R(a_ct)
+
+ ltm = CT_LTM(ct)
+ ltv = CT_LTV(ct)
+
+ if (CT_TYPE(ct) == LNR) {
+ # Simple linear, nonrotated transformation.
+ x2 = Memr[ltm ] * x1 + Memr[ltv ]
+ y2 = Memr[ltm+3] * y1 + Memr[ltv+1]
+ } else if (CT_TYPE(ct) == LRO) {
+ # Linear, rotated transformation.
+ p1[1] = x1; p1[2] = y1
+ x2 = Memr[ltm ] * p1[1] + Memr[ltm+1] * p1[2] + Memr[ltv ]
+ y2 = Memr[ltm+2] * p1[1] + Memr[ltm+3] * p1[2] + Memr[ltv+1]
+ } else {
+ # General case involving one or more functional terms.
+ p1[1] = x1; p1[2] = y1
+ call mw_ctranr (a_ct, p1, p2, 2)
+ x2 = p2[1]; y2 = p2[2]
+ }
+end
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
diff --git a/sys/mwcs/gen/mwctranr.x b/sys/mwcs/gen/mwctranr.x
new file mode 100644
index 00000000..0574a563
--- /dev/null
+++ b/sys/mwcs/gen/mwctranr.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_ctranr (a_ct, p1, p2, ndim)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real p1[ndim] #I coordinates of point in input system
+real 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_R(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] = Memr[ltm+(i-1)*(ndim+1)] * p1[i] + Memr[ltv+i-1]
+ return
+ } else if (CT_TYPE(ct) == LRO) {
+ # Simple linear, rotated transformation.
+ call mw_ltranr (p1, p2, Memr[ltm], Memr[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 achtrd (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 achtdr (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
diff --git a/sys/mwcs/gen/mwgctrand.x b/sys/mwcs/gen/mwgctrand.x
new file mode 100644
index 00000000..cfdca886
--- /dev/null
+++ b/sys/mwcs/gen/mwgctrand.x
@@ -0,0 +1,44 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_GCTRAN -- Get a coordinate transformation compiled in a previous call
+# to mw_sctran. When the transformation is compiled, it is specified by
+# naming the input and output systems and the axes over which the transform
+# is to be performed. Rather than return this information, which the
+# application already knows, we return the actual transform, i.e., the
+# linear transformation matrix and translation vector comprising the linear
+# portion of the transform, and axis class arrays for the input and output
+# systems defining the axis types. If the axis types are all zero, there
+# are no WCS function calls for any axis in either system, and the
+# transformation is completely linear (hence computable by the application
+# if desired, e.g., with mw_ltr).
+
+int procedure mw_gctrand (a_ct, o_ltm, o_ltv, axtype1, axtype2, maxdim)
+
+pointer a_ct #I pointer to CTRAN descriptor
+double o_ltm[ARB] #O linear tranformation matrix
+double o_ltv[ARB] #O translation matrix
+int axtype1[ARB] #O axis types for input system
+int axtype2[ARB] #O axis types for output system
+int maxdim #I how much stuff to return
+
+pointer ct
+int pdim, ndim, i, j
+
+begin
+ ct = CT_D(a_ct)
+ pdim = CT_NDIM(ct)
+ ndim = min (pdim, maxdim)
+
+ # Output the goods.
+ do j = 1, ndim {
+ axtype1[j] = WCS_AXCLASS(CT_WCSI(ct),j)
+ axtype2[j] = WCS_AXCLASS(CT_WCSO(ct),j)
+ o_ltv[j] = Memd[CT_LTV(ct)+(j-1)]
+ do i = 1, ndim
+ o_ltm[(j-1)*ndim+i] = Memd[CT_LTM(ct)+(j-1)*pdim+(i-1)]
+ }
+
+ return (pdim)
+end
diff --git a/sys/mwcs/gen/mwgctranr.x b/sys/mwcs/gen/mwgctranr.x
new file mode 100644
index 00000000..7825c6df
--- /dev/null
+++ b/sys/mwcs/gen/mwgctranr.x
@@ -0,0 +1,44 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_GCTRAN -- Get a coordinate transformation compiled in a previous call
+# to mw_sctran. When the transformation is compiled, it is specified by
+# naming the input and output systems and the axes over which the transform
+# is to be performed. Rather than return this information, which the
+# application already knows, we return the actual transform, i.e., the
+# linear transformation matrix and translation vector comprising the linear
+# portion of the transform, and axis class arrays for the input and output
+# systems defining the axis types. If the axis types are all zero, there
+# are no WCS function calls for any axis in either system, and the
+# transformation is completely linear (hence computable by the application
+# if desired, e.g., with mw_ltr).
+
+int procedure mw_gctranr (a_ct, o_ltm, o_ltv, axtype1, axtype2, maxdim)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real o_ltm[ARB] #O linear tranformation matrix
+real o_ltv[ARB] #O translation matrix
+int axtype1[ARB] #O axis types for input system
+int axtype2[ARB] #O axis types for output system
+int maxdim #I how much stuff to return
+
+pointer ct
+int pdim, ndim, i, j
+
+begin
+ ct = CT_R(a_ct)
+ pdim = CT_NDIM(ct)
+ ndim = min (pdim, maxdim)
+
+ # Output the goods.
+ do j = 1, ndim {
+ axtype1[j] = WCS_AXCLASS(CT_WCSI(ct),j)
+ axtype2[j] = WCS_AXCLASS(CT_WCSO(ct),j)
+ o_ltv[j] = Memr[CT_LTV(ct)+(j-1)]
+ do i = 1, ndim
+ o_ltm[(j-1)*ndim+i] = Memr[CT_LTM(ct)+(j-1)*pdim+(i-1)]
+ }
+
+ return (pdim)
+end
diff --git a/sys/mwcs/gen/mwltrand.x b/sys/mwcs/gen/mwltrand.x
new file mode 100644
index 00000000..d35670c7
--- /dev/null
+++ b/sys/mwcs/gen/mwltrand.x
@@ -0,0 +1,26 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_LTRAN -- Perform a general N-dimensional linear transformation, i.e.,
+# matrix multiply and translation.
+
+procedure mw_ltrand (p1, p2, ltm, ltv, ndim)
+
+double p1[ndim] #I input point
+double p2[ndim] #O transformed output point
+double ltm[ndim,ndim] #I linear transformation matrix
+double ltv[ndim] #I linear translation vector
+int ndim #I dimension of system
+
+int i, j
+double p3[MAX_DIM]
+
+begin
+ call amovd (p1, p3, ndim)
+ do j = 1, ndim {
+ p2[j] = ltv[j]
+ do i = 1, ndim
+ p2[j] = p2[j] + ltm[i,j] * p3[i]
+ }
+end
diff --git a/sys/mwcs/gen/mwltranr.x b/sys/mwcs/gen/mwltranr.x
new file mode 100644
index 00000000..9cafe4d2
--- /dev/null
+++ b/sys/mwcs/gen/mwltranr.x
@@ -0,0 +1,26 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_LTRAN -- Perform a general N-dimensional linear transformation, i.e.,
+# matrix multiply and translation.
+
+procedure mw_ltranr (p1, p2, ltm, ltv, ndim)
+
+real p1[ndim] #I input point
+real p2[ndim] #O transformed output point
+real ltm[ndim,ndim] #I linear transformation matrix
+real ltv[ndim] #I linear translation vector
+int ndim #I dimension of system
+
+int i, j
+real p3[MAX_DIM]
+
+begin
+ call amovr (p1, p3, ndim)
+ do j = 1, ndim {
+ p2[j] = ltv[j]
+ do i = 1, ndim
+ p2[j] = p2[j] + ltm[i,j] * p3[i]
+ }
+end
diff --git a/sys/mwcs/gen/mwmmuld.x b/sys/mwcs/gen/mwmmuld.x
new file mode 100644
index 00000000..ae35f082
--- /dev/null
+++ b/sys/mwcs/gen/mwmmuld.x
@@ -0,0 +1,21 @@
+# MW_MMUL -- Matrix multiply.
+
+procedure mw_mmuld (a, b, c, ndim)
+
+double a[ndim,ndim] #I left input matrix
+double b[ndim,ndim] #I right input matrix
+double c[ndim,ndim] #O output matrix
+int ndim #I dimensionality of system
+
+int i, j, k
+double v
+
+begin
+ do j = 1, ndim
+ do i = 1, ndim {
+ v = 0
+ do k = 1, ndim
+ v = v + a[k,j] * b[i,k]
+ c[i,j] = v
+ }
+end
diff --git a/sys/mwcs/gen/mwmmulr.x b/sys/mwcs/gen/mwmmulr.x
new file mode 100644
index 00000000..83e14d2c
--- /dev/null
+++ b/sys/mwcs/gen/mwmmulr.x
@@ -0,0 +1,21 @@
+# MW_MMUL -- Matrix multiply.
+
+procedure mw_mmulr (a, b, c, ndim)
+
+real a[ndim,ndim] #I left input matrix
+real b[ndim,ndim] #I right input matrix
+real c[ndim,ndim] #O output matrix
+int ndim #I dimensionality of system
+
+int i, j, k
+real v
+
+begin
+ do j = 1, ndim
+ do i = 1, ndim {
+ v = 0
+ do k = 1, ndim
+ v = v + a[k,j] * b[i,k]
+ c[i,j] = v
+ }
+end
diff --git a/sys/mwcs/gen/mwv1trand.x b/sys/mwcs/gen/mwv1trand.x
new file mode 100644
index 00000000..3c3ac124
--- /dev/null
+++ b/sys/mwcs/gen/mwv1trand.x
@@ -0,0 +1,32 @@
+include "../mwcs.h"
+
+# MW_V1TRAN -- Optimized 1D coordinate transformation for an array of points.
+
+procedure mw_v1trand (a_ct, x1, x2, npts)
+
+pointer a_ct #I pointer to CTRAN descriptor
+double x1[ARB] #I coordinates in input system
+double x2[ARB] #O coordinates in output system
+int npts
+
+int i
+pointer ct
+double scale, offset
+errchk mw_ctrand
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_D(a_ct)
+
+ scale = Memd[CT_LTM(ct)]
+ offset = Memd[CT_LTV(ct)]
+
+ # Perform the transformation; case LNR is a simple linear transform.
+ if (CT_TYPE(ct) == LNR) {
+ do i = 1, npts
+ x2[i] = scale * x1[i] + offset
+ } else {
+ do i = 1, npts
+ call mw_ctrand (a_ct, x1[i], x2[i], 1)
+ }
+end
diff --git a/sys/mwcs/gen/mwv1tranr.x b/sys/mwcs/gen/mwv1tranr.x
new file mode 100644
index 00000000..045f6a33
--- /dev/null
+++ b/sys/mwcs/gen/mwv1tranr.x
@@ -0,0 +1,32 @@
+include "../mwcs.h"
+
+# MW_V1TRAN -- Optimized 1D coordinate transformation for an array of points.
+
+procedure mw_v1tranr (a_ct, x1, x2, npts)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real x1[ARB] #I coordinates in input system
+real x2[ARB] #O coordinates in output system
+int npts
+
+int i
+pointer ct
+real scale, offset
+errchk mw_ctranr
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_R(a_ct)
+
+ scale = Memr[CT_LTM(ct)]
+ offset = Memr[CT_LTV(ct)]
+
+ # Perform the transformation; case LNR is a simple linear transform.
+ if (CT_TYPE(ct) == LNR) {
+ do i = 1, npts
+ x2[i] = scale * x1[i] + offset
+ } else {
+ do i = 1, npts
+ call mw_ctranr (a_ct, x1[i], x2[i], 1)
+ }
+end
diff --git a/sys/mwcs/gen/mwv2trand.x b/sys/mwcs/gen/mwv2trand.x
new file mode 100644
index 00000000..3a1cf329
--- /dev/null
+++ b/sys/mwcs/gen/mwv2trand.x
@@ -0,0 +1,49 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_V2TRAN -- Optimized 2D coordinate transformation for an array of points.
+
+procedure mw_v2trand (a_ct, x1,y1, x2,y2, npts)
+
+pointer a_ct #I pointer to CTRAN descriptor
+double x1[ARB],y1[ARB] #I coordinates in input system
+double x2[ARB],y2[ARB] #O coordinates in output system
+int npts
+
+int i
+pointer ct, ltm, ltv
+double p1[2], p2[2]
+errchk mw_ctrand
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_D(a_ct)
+
+ ltm = CT_LTM(ct)
+ ltv = CT_LTV(ct)
+
+ if (CT_TYPE(ct) == LNR) {
+ # Simple linear, nonrotated transformation.
+ do i = 1, npts {
+ x2[i] = Memd[ltm ] * x1[i] + Memd[ltv ]
+ y2[i] = Memd[ltm+3] * y1[i] + Memd[ltv+1]
+ }
+ } else if (CT_TYPE(ct) == LRO) {
+ # Linear, rotated transformation.
+ do i = 1, npts {
+ p1[1] = x1[i]; p1[2] = y1[i]
+ x2[i] = Memd[ltm ] * p1[1] + Memd[ltm+1] * p2[1] +
+ Memd[ltv ]
+ y2[i] = Memd[ltm+2] * p1[1] + Memd[ltm+3] * p2[1] +
+ Memd[ltv+1]
+ }
+ } else {
+ # General case involving one or more functional terms.
+ do i = 1, npts {
+ p1[1] = x1[i]; p1[2] = y1[i]
+ call mw_ctrand (a_ct, p1, p2, 2)
+ x2[i] = p2[1]; y2[i] = p2[2]
+ }
+ }
+end
diff --git a/sys/mwcs/gen/mwv2tranr.x b/sys/mwcs/gen/mwv2tranr.x
new file mode 100644
index 00000000..dc2fe58f
--- /dev/null
+++ b/sys/mwcs/gen/mwv2tranr.x
@@ -0,0 +1,49 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "../mwcs.h"
+
+# MW_V2TRAN -- Optimized 2D coordinate transformation for an array of points.
+
+procedure mw_v2tranr (a_ct, x1,y1, x2,y2, npts)
+
+pointer a_ct #I pointer to CTRAN descriptor
+real x1[ARB],y1[ARB] #I coordinates in input system
+real x2[ARB],y2[ARB] #O coordinates in output system
+int npts
+
+int i
+pointer ct, ltm, ltv
+real p1[2], p2[2]
+errchk mw_ctranr
+
+begin
+ # Get real or double version of descriptor.
+ ct = CT_R(a_ct)
+
+ ltm = CT_LTM(ct)
+ ltv = CT_LTV(ct)
+
+ if (CT_TYPE(ct) == LNR) {
+ # Simple linear, nonrotated transformation.
+ do i = 1, npts {
+ x2[i] = Memr[ltm ] * x1[i] + Memr[ltv ]
+ y2[i] = Memr[ltm+3] * y1[i] + Memr[ltv+1]
+ }
+ } else if (CT_TYPE(ct) == LRO) {
+ # Linear, rotated transformation.
+ do i = 1, npts {
+ p1[1] = x1[i]; p1[2] = y1[i]
+ x2[i] = Memr[ltm ] * p1[1] + Memr[ltm+1] * p2[1] +
+ Memr[ltv ]
+ y2[i] = Memr[ltm+2] * p1[1] + Memr[ltm+3] * p2[1] +
+ Memr[ltv+1]
+ }
+ } else {
+ # General case involving one or more functional terms.
+ do i = 1, npts {
+ p1[1] = x1[i]; p1[2] = y1[i]
+ call mw_ctranr (a_ct, p1, p2, 2)
+ x2[i] = p2[1]; y2[i] = p2[2]
+ }
+ }
+end
diff --git a/sys/mwcs/gen/mwvmuld.x b/sys/mwcs/gen/mwvmuld.x
new file mode 100644
index 00000000..0af8dfa7
--- /dev/null
+++ b/sys/mwcs/gen/mwvmuld.x
@@ -0,0 +1,20 @@
+# MW_VMUL -- Vector multiply.
+
+procedure mw_vmuld (a, b, c, ndim)
+
+double a[ndim,ndim] #I input matrix
+double b[ndim] #I input vector
+double c[ndim] #O output vector
+int ndim #I system dimension
+
+int i, j
+double v
+
+begin
+ do j = 1, ndim {
+ v = 0
+ do i = 1, ndim
+ v = v + a[i,j] * b[i]
+ c[j] = v
+ }
+end
diff --git a/sys/mwcs/gen/mwvmulr.x b/sys/mwcs/gen/mwvmulr.x
new file mode 100644
index 00000000..54a0776e
--- /dev/null
+++ b/sys/mwcs/gen/mwvmulr.x
@@ -0,0 +1,20 @@
+# MW_VMUL -- Vector multiply.
+
+procedure mw_vmulr (a, b, c, ndim)
+
+real a[ndim,ndim] #I input matrix
+real b[ndim] #I input vector
+real c[ndim] #O output vector
+int ndim #I system dimension
+
+int i, j
+real v
+
+begin
+ do j = 1, ndim {
+ v = 0
+ do i = 1, ndim
+ v = v + a[i,j] * b[i]
+ c[j] = v
+ }
+end
diff --git a/sys/mwcs/gen/mwvtrand.x b/sys/mwcs/gen/mwvtrand.x
new file mode 100644
index 00000000..1a1cb662
--- /dev/null
+++ b/sys/mwcs/gen/mwvtrand.x
@@ -0,0 +1,18 @@
+# MW_VTRAN -- Transform an array of N-dimensional points, expressed as a
+# 2D vector where v[1,i] is point I of vector V.
+
+procedure mw_vtrand (ct, v1, v2, ndim, npts)
+
+pointer ct #I pointer to CTRAN descriptor
+double v1[ndim,npts] #I points to be transformed
+double v2[ndim,npts] #O vector to get the transformed points
+int ndim #I dimensionality of each point
+int npts #I number of points
+
+int i
+errchk mw_ctrand
+
+begin
+ do i = 1, npts
+ call mw_ctrand (ct, v1[1,i], v2[1,i], ndim)
+end
diff --git a/sys/mwcs/gen/mwvtranr.x b/sys/mwcs/gen/mwvtranr.x
new file mode 100644
index 00000000..ca705c8b
--- /dev/null
+++ b/sys/mwcs/gen/mwvtranr.x
@@ -0,0 +1,18 @@
+# MW_VTRAN -- Transform an array of N-dimensional points, expressed as a
+# 2D vector where v[1,i] is point I of vector V.
+
+procedure mw_vtranr (ct, v1, v2, ndim, npts)
+
+pointer ct #I pointer to CTRAN descriptor
+real v1[ndim,npts] #I points to be transformed
+real v2[ndim,npts] #O vector to get the transformed points
+int ndim #I dimensionality of each point
+int npts #I number of points
+
+int i
+errchk mw_ctranr
+
+begin
+ do i = 1, npts
+ call mw_ctranr (ct, v1[1,i], v2[1,i], ndim)
+end