1 files changed, 85 insertions, 0 deletions
diff --git a/noao/imred/vtel/lstsq.x b/noao/imred/vtel/lstsq.x
new file mode 100644
index 00000000..9091fb48
--- /dev/null
+++ b/noao/imred/vtel/lstsq.x
@@ -0,0 +1,85 @@
+include	<mach.h>
+
+# LSTSQ -- Do a least squares fit to the data contained in the zz array.
+# Algorithm is from Jack Harvey. (Yes, it's a black box...)
+
+procedure lstsq (zz, mz, fno)
+
+real	zz[mz, mz]
+int	mz
+real	fno
+
+int	n, m, m1, i, j, k, l, l1
+real	fn, pp
+
+begin
+	n = mz - 2
+	m = n + 1
+	m1 = m + 1
+	fn = n
+
+	do i = 1, m {
+	    l = i + 1
+	    do k = 1, i-1 {
+		zz[i,l] = zz[i,l] - zz[k,l]**2
+	    }
+
+	    if (i == m)
+		break
+	    if (zz[i,l] >= 0.0)
+	        zz[i,l] = zz[i,l]**.5
+	    else {
+		call eprintf ("square root of negitive number in lstsq\n")
+	        zz[i,l] = 0.0
+	    }
+	    l1 = l + 1
+
+	    do j = l1, m1 {
+		do k = 1, i-1 {
+		    zz[i,j] = zz[i,j] - zz[k,l] * zz[k,j]
+		}
+		if (zz[i,l] >= EPSILONR)
+		    zz[i,j] = zz[i,j] / zz[i,l]
+		else
+		    call eprintf ("divide by zero in lstsq\n")
+	    }
+
+	    if (zz[i,l] >= EPSILONR)
+	        zz[i,i] = 1. / zz[i,l]
+	    else
+		call eprintf ("divide by zero in lstsq\n")
+	    do j = 1, i-1 {
+		pp = 0.
+		l1 = i - 1
+		do k = j, l1 {
+		    pp = pp + zz[k,l] * zz[k,j]
+		}
+	        zz[i,j] = -zz[i,i] * pp
+	    }
+	}
+
+	if ((fno - fn) >= EPSILONR)
+	    if ((zz[m,m1] / (fno - fn)) >= 0.0)
+	        zz[m1,m1] = .6745 * (zz[m,m1] / (fno - fn))**.5
+	    else {
+		call eprintf ("square root of negitive number in lstsq\n")
+	        zz[m1,m1] = 0.0
+	    }
+	else
+	    call eprintf ("divide by zero in lstsq\n")
+
+	do i = 1, n {
+	    zz[m,i] = 0.
+	    pp = 0.
+	    do j = i, n {
+		zz[m,i] = zz[m,i] + zz[j,i] * zz[j,m1]
+		pp = pp + zz[j,i] * zz[j,i]
+	    }
+	    if (pp >= 0.0)
+	        zz[m1,i] = zz[m1,m1] * pp**.5
+	    else {
+		call eprintf ("square root of negitive number in lstsq\n")
+	        zz[m1,i] = 0.0
+	    }
+	}
+end