1 files changed, 97 insertions, 0 deletions
diff --git a/pkg/images/lib/rgxymatch.x b/pkg/images/lib/rgxymatch.x
new file mode 100644
index 00000000..3b2b45cb
--- /dev/null
+++ b/pkg/images/lib/rgxymatch.x
@@ -0,0 +1,97 @@
+include <mwset.h>
+
+# RG_RXYL -- Compute the grid of logical coordinates.
+
+procedure rg_rxyl (xl, yl, nx, ny, x1, x2, y1, y2)
+
+double	xl[ARB]			#O array of output x coordinates
+double	yl[ARB]			#O array of output y coordinates
+int	nx			#I the size of the grid in x
+int	ny			#I the size of the grid in y
+double	x1			#I the lower limit of the grid in x
+double	x2			#I the upper limit of the grid in x
+double	y1			#I the lower limit of the grid in y
+double	y2			#I the upper limit of the grid in y
+
+double	xstep, ystep, x, y
+int	i, j, npts
+
+begin
+	if (nx == 1)
+	    xstep = 0.0d0
+	else
+	    xstep = (x2 - x1) / (nx - 1)
+	if (ny == 1)
+	    ystep = 0.0d0
+	else
+	    ystep = (y2 - y1) / (ny - 1)
+	npts = 0
+
+	y = y1
+	do j = 1, ny {
+	    x = x1
+	    do i = 1, nx {
+		npts = npts + 1
+		xl[npts] = x
+		yl[npts] = y
+		x = x + xstep
+	    }
+	    y = y + ystep
+	}
+end
+
+
+# RG_XYTOXY -- Compute the world coordinate list give the wcs descriptor
+# and the logical coordinates.
+
+pointer procedure rg_xytoxy (mw, xl, yl, xw, yw, npts, inwcs, outwcs, ax1, ax2)
+
+pointer	mw			#I the wcs descriptor
+double	xl[ARB]			#I the input logical x coordinate
+double	yl[ARB]			#I the input logical y coordinate
+double	xw[ARB]			#O the output world x coordinate
+double	yw[ARB]			#O the output world y coordinate
+int	npts			#I the number of coordinates.
+char	inwcs[ARB]		#I the input wcs
+char	outwcs[ARB]		#I the output wcs
+int	ax1			#I the logical x axis
+int	ax2			#I the logical y axis
+
+int	i, axbits
+pointer	ct
+double	mw_c1trand()
+int	mw_stati()
+pointer	mw_sctran()
+errchk	mw_sctran()
+
+begin
+	# Compile the transformation.
+	if (mw == NULL) {
+	    ct = NULL
+	} else if (mw_stati (mw, MW_NDIM) >= 2) {
+	    axbits = 2 ** (ax1 - 1) + 2 ** (ax2 - 1)
+	    iferr (ct = mw_sctran (mw, inwcs, outwcs, axbits))
+	        ct = NULL
+	} else {
+	    axbits = 2 ** (ax1 - 1)
+	    iferr (ct = mw_sctran (mw, inwcs, outwcs, axbits))
+	        ct = NULL
+	}
+
+	# Compute the world coordinates.
+	if (ct == NULL) {
+	    call amovd (xl, xw, npts)
+	    call amovd (yl, yw, npts)
+	} else if (mw_stati (mw, MW_NDIM) == 2) {
+	    do i = 1, npts
+		call mw_c2trand (ct, xl[i], yl[i], xw[i], yw[i])
+	} else {
+	    do i = 1, npts {
+		xw[i] = mw_c1trand (ct, xl[i])
+		yw[i] = yl[i]
+	    }
+	}
+
+	return (ct)
+end
+