1 files changed, 102 insertions, 0 deletions
diff --git a/noao/nproto/ace/noisemodel.x b/noao/nproto/ace/noisemodel.x
new file mode 100644
index 00000000..0503c4a1
--- /dev/null
+++ b/noao/nproto/ace/noisemodel.x
@@ -0,0 +1,102 @@
+# NOISEMODEL -- Compute noise model.
+#
+# var = (var(sky) + (image-sky)/gain) / sqrt (exposure)
+#
+# What is actually returned is the square root of the variance.
+# The variance of the sky and the effective gain are for a unit
+# exposure in the exposure map.
+
+procedure noisemodel (image, sky, sig, gain, exp, sigma, npix)
+
+real	image[npix]		#I Image
+real	sky[npix]		#I Sky
+real	sig[npix]		#I Sky sigma
+real	gain[npix]		#I Gain
+real	exp[npix]		#I Exposure
+real	sigma[npix]		#O Sigma
+int	npix			#I Number of pixels
+
+int	i
+real	e, elast, sqrte
+
+begin
+	if (IS_INDEFR(exp[1])) {
+	    if (IS_INDEFR(gain[1]))
+		call amovr (sig, sigma, npix)
+	    else {
+		do i = 1, npix
+		    sigma[i] = sqrt (sig[i] * sig[1] +
+			(image[i] - sky[i]) / gain[i])
+	    }
+	} else if (IS_INDEFR(gain[1])) {
+	    elast = INDEFR
+	    do i = 1, npix {
+		e = exp[i]
+		if (e == 0.) {
+		    sigma[i] = sig[i]
+		    next
+		}
+		if (e != elast) {
+		    sqrte = sqrt (e)
+		    elast = e
+		}
+		sigma[i] = sig[i] / sqrte
+	    }
+	} else {
+	    do i = 1, npix {
+		e = exp[i]
+		if (e == 0.) {
+		    sigma[i] = sqrt (sig[i] * sig[i] +
+			(image[i] - sky[i]) / gain[i])
+		    next
+		}
+		sigma[i] = sqrt ((sig[i] * sig[i] +
+		    (image[i] - sky[i]) / gain[i]) / e)
+	    }
+	}
+end
+
+
+# EXPSIGMA -- Apply exposure map to correct sky sigma.
+# Assume the exposure map has region of contiguous constant values so
+# that the number of square roots can be minimized.  An exposure map
+# value of zero leaves the sigma unchanged.
+
+procedure expsigma (sigma, expmap, npix, mode)
+
+real	sigma[npix]		#U Sigma values
+real	expmap[npix]		#I Exposure map values
+int	npix			#I Number of pixels
+int	mode			#I 0=divide, 1=multiply
+
+int	i
+real	exp, lastexp, scale
+
+begin
+	switch (mode) {
+	case 0:
+	    lastexp = INDEFR
+	    do i = 1, npix {
+		exp = expmap[i]
+		if (exp == 0.)
+		    next
+		if (exp != lastexp) {
+		    scale = sqrt (exp)
+		    lastexp = exp
+		}
+		sigma[i] = sigma[i] / scale
+	    }
+	case 1:
+	    lastexp = INDEFR
+	    do i = 1, npix {
+		exp = expmap[i]
+		if (exp == 0.)
+		    next
+		if (exp != lastexp) {
+		    scale = sqrt (exp)
+		    lastexp = exp
+		}
+		sigma[i] = sigma[i] * scale
+	    }
+	}
+end