Repatch (from linux) of OSX IRAF

1 files changed, 170 insertions, 0 deletions

diff --git a/noao/imred/quadred/src/quad/gainmeasure.x b/noao/imred/quadred/src/quad/gainmeasure.x
new file mode 100644
index 00000000..592cf7cf
--- /dev/null
+++ b/noao/imred/quadred/src/quad/gainmeasure.x
@@ -0,0 +1,170 @@
+include <imhdr.h>
+include "quadgeom.h"
+
+define  OPT_REFLECT     4
+
+# GAINMEASURE -- Calculate the gain (e/ADU) and RON of a CCD using the CTIO
+# (Bruce Attwood) algorithm.
+# Input are a pair of high signal level exposures (Flat1, Flat2) and a pair of
+# zero exposures (Zero1, Zero2). We then calculate:
+#
+#    epadu = <Flat1> + <Flat2> - (<Zero1> + <Zero2>) / var{Diff_F} - Var{Diff_Z}
+#    RON   = RMS {Diff_Z} * epadu / sqrt(2)
+#
+# Where:
+#
+#    diff_Z = Zero1 - Zero2
+#    diff_F = Flat1 - Flat2
+#
+# The statistics must be calculated for regions free of bad pixels and other
+# defects, and with reasonably uniform illumination. 
+
+procedure t_gainmeasure ()
+
+pointer	flat1, flat2		#TI High level images
+pointer	zero1, zero2		#TI Zero level images
+char	section[SZ_LINE]	#TI Section for calculation
+
+char	buffer[SZ_LINE]
+int	npix, x1, x2, y1, y2, amp
+pointer	sp, f1, f2, z1, z2, fd, zd, qg
+real	f1bar,   f2bar,   z1bar,   z2bar,   fdbar,   zdbar
+real	f1sigma, f2sigma, z1sigma, z2sigma, fdsigma, zdsigma
+real	div, epadu, ron
+bool	headers
+
+pointer	immap(), imgs2r()
+bool	clgetb(), quadsect()
+int	hdmaccf()
+
+begin
+        # Open instrument file
+        call clgstr    ("instrument",  buffer,  SZ_FNAME)
+        call hdmopen   (buffer)
+
+	# Map input images
+	call clgstr ("flat1",   buffer, SZ_LINE)
+	flat1 = immap (buffer,  READ_ONLY, 0)
+
+	call clgstr ("flat2",   buffer, SZ_LINE)
+	flat2 = immap (buffer,  READ_ONLY, 0)
+
+	call clgstr ("zero1",   buffer, SZ_LINE)
+	zero1 = immap (buffer,  READ_ONLY, 0)
+
+	call clgstr ("zero2",   buffer, SZ_LINE)
+	zero2 = immap (buffer,  READ_ONLY, 0)
+
+        # Get section over which measurement is to be made.
+        call clgstr ("section", section, SZ_LINE)
+
+	# See if headers are to be printed
+	headers = clgetb ("print_headers")
+
+        # Set-up quadgeom structure. We blithely assume all images are the same.
+        call quadalloc (qg)
+
+	if (hdmaccf (flat1, "HDR_REV") == NO) {
+	    call quadgeom  (flat1, qg, "", "")
+	} else {
+	    call qghdr2 (flat1, qg)
+	}
+#       call quaddump  (qg)
+
+	if (headers) {
+	    call printf ("#")
+	    do amp = 1, QG_NAMPS (qg) {
+		call printf ("%9wAmp%2s%5w")
+		    call pargstr (Memc[QG_AMPID (qg, amp)])
+	    }
+	    call printf ("\n")
+
+	    call printf ("#")
+	    do amp = 1, QG_NAMPS (qg) {
+		call printf ("%5wGain%4wRON%3w")
+	    }
+	    call printf ("\n")
+	    call printf ("#")
+	    do amp = 1, QG_NAMPS (qg) {
+		call printf ("%3w(e-/ADU)%2w(e-)%2w")
+	    }
+	    call printf ("\n")
+	}
+
+	call printf ("%1w")
+        do amp = 1, QG_NAMPS (qg) {
+ 
+            if (quadsect (qg, section, OPT_REFLECT, amp, x1, x2, y1, y2)) {
+
+		npix = (abs(y2 - y1) + 1) * (abs(x2 - x1) + 1)
+
+		# Allocate working arrays
+		call smark  (sp)
+		call salloc (fd, npix, TY_REAL)
+		call salloc (zd, npix, TY_REAL)
+
+		# Read data 
+		f1 = imgs2r (flat1, x1, x2, y1, y2)
+		f2 = imgs2r (flat2, x1, x2, y1, y2)
+		z1 = imgs2r (zero1, x1, x2, y1, y2)
+		z2 = imgs2r (zero2, x1, x2, y1, y2)
+
+		# Calculate differences
+		call asubr (Memr[f1], Memr[f2], Memr[fd], npix)
+		call asubr (Memr[z1], Memr[z2], Memr[zd], npix)
+
+		# Calculate means and standard deviations
+		call aavgr (Memr[f1], npix, f1bar, f1sigma)
+		call aavgr (Memr[f2], npix, f2bar, f2sigma)
+		call aavgr (Memr[z1], npix, z1bar, z1sigma)
+		call aavgr (Memr[z2], npix, z2bar, z2sigma)
+		call aavgr (Memr[fd], npix, fdbar, fdsigma)
+		call aavgr (Memr[zd], npix, zdbar, zdsigma)
+
+#		call eprintf ("f1bar=%g f1sigma=%g\n")
+#		    call pargr (f1bar)
+#		    call pargr (f1sigma)
+#		call eprintf ("f2bar=%g f2sigma=%g\n")
+#		    call pargr (f2bar)
+#		    call pargr (f2sigma)
+#		call eprintf ("z1bar=%g z1sigma=%g\n")
+#		    call pargr (z1bar)
+#		    call pargr (z1sigma)
+#		call eprintf ("z2bar=%g z2sigma=%g\n")
+#		    call pargr (z2bar)
+#		    call pargr (z2sigma)
+#		call eprintf ("fdbar=%g fdsigma=%g\n")
+#		    call pargr (fdbar)
+#		    call pargr (fdsigma)
+#		call eprintf ("zdbar=%g zdsigma=%g\n")
+#		    call pargr (zdbar)
+#		    call pargr (zdsigma)
+
+		div  = fdsigma**2 - zdsigma**2
+		if (div > 0.0) {
+		    epadu = ((f1bar + f2bar) - (z1bar + z2bar)) / div
+		    ron   = epadu * zdsigma / 1.41421356
+		} else {
+		    epadu = INDEF
+		    ron   = INDEF
+		}
+
+		# Print results
+		call printf ("%3w%6.2f%2w%6.2f%2w")
+		    call pargr (epadu)
+		    call pargr (ron)
+
+	    # Free working arrays
+	    call sfree (sp)
+
+	    }
+	}
+
+	call printf ("\n")
+
+	# Tidy up
+	call imunmap (flat1)
+	call imunmap (flat2)
+	call imunmap (zero1)
+	call imunmap (zero2)
+end