1 files changed, 89 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/sigma.gx b/noao/imred/ccdred/src/sigma.gx
new file mode 100644
index 00000000..8b59f1f6
--- /dev/null
+++ b/noao/imred/ccdred/src/sigma.gx
@@ -0,0 +1,89 @@
+$for (sr)
+# SIGMA -- Compute sigma line from image lines with rejection.
+
+procedure sigma$t (data, nimages, mean, sigma, npts)
+
+pointer	data[nimages]		# Data vectors
+int	nimages			# Number of data vectors
+$if (datatype == sil)
+real	mean[npts]		# Mean vector
+real	sigma[npts]		# Sigma vector (returned)
+$else
+PIXEL	mean[npts]		# Mean vector
+PIXEL	sigma[npts]		# Sigma vector (returned)
+$endif
+int	npts			# Number of points in each vector
+
+$if (datatype == sil)
+real	val, sig, pixval
+$else
+PIXEL	val, sig, pixval
+$endif
+int	i, j, n, n1
+
+begin
+	n = nimages - 1
+	do i = 1, npts {
+	    val = mean[i]
+	    sig = 0.
+	    n1 = n
+	    do j = 1, nimages {
+		pixval = Mem$t[data[j]+i-1]
+		if (IS_INDEF (pixval))
+		    n1 = n1 - 1
+		else
+		    sig = sig + (pixval - val) ** 2
+	    }
+	    if (n1 > 0)
+	        sigma[i] = sqrt (sig / n1)
+	    else
+	        sigma[i] = 0.
+	}
+end
+
+
+# WTSIGMA -- Compute scaled and weighted sigma line from image lines with
+# rejection.
+
+procedure wtsigma$t (data, scales, zeros, wts, nimages, mean, sigma, npts)
+
+pointer	data[nimages]		# Data vectors
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero levels
+real	wts[nimages]		# Weights
+int	nimages			# Number of data vectors
+$if (datatype == sil)
+real	mean[npts]		# Mean vector
+real	sigma[npts]		# Sigma vector (returned)
+real	val, sig, pixval
+$else
+PIXEL	mean[npts]		# Mean vector
+PIXEL	sigma[npts]		# Sigma vector (returned)
+PIXEL	val, sig, pixval
+$endif
+int	npts			# Number of points in each vector
+
+int	i, j, n
+real	sumwts
+
+begin
+	do i = 1, npts {
+	    val = mean[i]
+	    n = 0
+	    sig = 0.
+	    sumwts = 0.
+	    do j = 1, nimages {
+		pixval = Mem$t[data[j]+i-1]
+		if (!IS_INDEF (pixval)) {
+		    n = n + 1
+		    sig = sig + wts[j]*(pixval/scales[j]-zeros[j]-val) ** 2
+		    sumwts = sumwts + wts[j]
+		}
+	    }
+	    if (n > 1)
+	        sigma[i] = sqrt (sig / sumwts * n / (n - 1))
+	    else
+	        sigma[i] = 0.
+	}
+end
+$endfor