1 files changed, 50 insertions, 0 deletions
diff --git a/sys/qpoe/qpiolwcs.x b/sys/qpoe/qpiolwcs.x
new file mode 100644
index 00000000..76c1650f
--- /dev/null
+++ b/sys/qpoe/qpiolwcs.x
@@ -0,0 +1,50 @@
+include	"qpio.h"
+
+# QPIO_LOADWCS -- Load the WCS, if any, from the QPOE file associated with the
+# given QPIO descriptor, into an open MWCS descriptor.  This is equivalent to
+# QP_LOADWCS except that the Lterm is updated to reflect the current blocking
+# factor and rect (if any) used for rasterization.  In the resultant WCS, the
+# logical coordinate system gives the pixel coordinates of the sampled rect.
+
+pointer	procedure qpio_loadwcs (io)
+
+pointer	io				#I QPIO descriptor
+
+pointer	qp, mw
+int	ndim, i, j
+double	ltv_1[NDIM], ltv_2[NDIM], ltm[NDIM,NDIM]
+pointer	qp_loadwcs()
+errchk	qp_loadwcs
+
+begin
+	qp = IO_QP(io)
+	mw = qp_loadwcs (qp)
+	ndim = NDIM
+
+	# Formalize the transformation.
+	ltv_1[1] = IO_VSDEF(io,1) - 1
+	ltv_1[2] = IO_VSDEF(io,2) - 1
+
+	# L(i) :==  LTM=(1 / block) * P(i)  +  Vx
+	# At pixel {P(i) :==  (block + 1) / 2}  L(i) is 1.0.
+	# Solve for Vx :==  1.0 - (1 / block) * ((block + 1) / 2)
+	#	       -->  0.5 - 1 / (2 * block)
+
+	ltv_2[1] = 0.5d0 - 1.0d0 / double (max (1.0, IO_XBLOCK(io))) / 2.0d0
+	ltv_2[2] = 0.5d0 - 1.0d0 / double (max (1.0, IO_YBLOCK(io))) / 2.0d0
+
+	do j = 1, ndim
+	    do i = 1, ndim
+		if (i == j) {
+		    if (i == 1)
+			ltm[i,j] = 1.0D0 / max (1.0, IO_XBLOCK(io))
+		    else
+			ltm[i,j] = 1.0D0 / max (1.0, IO_YBLOCK(io))
+		} else
+		    ltm[i,j] = 0
+
+	# Apply the transformation.
+	call mw_translated (mw, ltv_1, ltm, ltv_2, ndim)
+
+	return (mw)
+end