Initial commit

1 files changed, 261 insertions, 0 deletions

diff --git a/noao/imred/src/fibers/fibresponse.cl b/noao/imred/src/fibers/fibresponse.cl
new file mode 100644
index 00000000..2379fcc4
--- /dev/null
+++ b/noao/imred/src/fibers/fibresponse.cl
@@ -0,0 +1,261 @@
+# FIBRESPONSE -- Make a aperture response spectrum using a flat field
+# and a throughput file or image.
+
+procedure fibresponse (flat, throughput, apreference, response)
+
+string	flat			{prompt="Flat field spectrum"}
+string	throughput		{prompt="Throughput file or image"}
+string	apreference		{prompt="Aperture reference spectrum"}
+string	response		{prompt="Response spectrum"}
+
+bool	recenter = no		{prompt="Recenter apertures?"}
+bool	edit = yes		{prompt="Edit/review apertures?"}
+bool	trace = no		{prompt="Trace spectra?"}
+bool	clean = no		{prompt="Detect and replace bad pixels?"}
+bool	fitflat = no		{prompt="Fit and ratio flat field spectrum?"}
+bool	interactive = yes	{prompt="Fit flat field interactively?"}
+string	function = "spline3"	{prompt="Fitting function",
+				 enum="spline3|legendre|chebyshev|spline1"}
+int	order = 20		{prompt="Fitting function order", min=1}
+
+begin
+	string	imtype, mstype
+	file	flat2d, skyflat2d, apref, resp
+	file	temp1, temp2, log1, log2
+	int	n, ap, naxis
+	real	respval
+	struct	err
+
+	imtype = "." // envget ("imtype")
+	i = stridx (",", imtype)
+	if (i > 0)
+	    imtype = substr (imtype, 1, i-1)
+	mstype = ".ms" // imtype
+	n = strlen (imtype)
+
+	flat2d = flat
+	skyflat2d = throughput
+	apref = apreference
+	resp = response
+	temp1 = mktemp ("tmp")
+	temp2 = mktemp ("tmp")
+
+	# Check required input and output.
+ 	if (resp == "" || resp == flat2d || resp == skyflat2d)
+ 	    error (1, "Bad response image name")
+	if (flat2d == "" && skyflat2d == "")
+	    error (1, "No flat field or throughput specified")
+
+	if (flat2d != "") {
+	    i = strlen (flat2d)
+	    if (i > n && substr (flat2d, i-n+1, i) == imtype)
+		flat2d = substr (flat2d, 1, i-n)
+	    if (!access (flat2d // imtype)) {
+		printf ("Flat field spectrum not found - %s%s\n",
+		    flat2d, imtype) | scan (err)
+		error (1, err // "\nCheck settting of imtype")
+	    }
+	}
+	if (skyflat2d != "") {
+	    i = strlen (skyflat2d)
+	    if (i > n && substr (skyflat2d, i-n+1, i) == imtype)
+	        skyflat2d = substr (skyflat2d, 1, i-n)
+	    if (!access (skyflat2d // imtype)) {
+		if (!access (skyflat2d)) {
+		    printf ("Throughput file or image not found - %s%s\n",
+			skyflat2d, imtype) | scan (err)
+		    error (1, err // "\nCheck settting of imtype")
+		}
+
+		if (flat2d == "") {
+		    i = strlen (apref)
+		    if (i > n && substr (apref, i-n+1, i) == imtype)
+			apref = substr (apref, 1, i-n)
+		    if (!access (apref // imtype))
+			error (1, "Aperture reference image required")
+		}
+	    }
+	}
+
+	# Set logging
+	tee.append = yes
+	if (logfile == "")
+	    log1 = "dev$null"
+	else
+	    log1 = logfile
+	if (verbose)
+	    log2 = "STDOUT"
+	else
+	    log2 = "dev$null"
+
+	# Initialize APSCRIPT
+	apscript.references = apref
+	if (recenter)
+	    apscript.ansrecenter = "YES"
+	else
+	    apscript.ansrecenter = "NO"
+	apscript.ansresize = "NO"
+	if (edit)
+	    apscript.ansedit = "yes"
+	else
+	    apscript.ansedit = "NO"
+	if (trace)
+	    apscript.anstrace = "YES"
+	else
+	    apscript.anstrace = "NO"
+	apscript.ansextract = "YES"
+
+	# If using a flat field extract it if necessary and possibly fit it
+	# and ratio the individual apertures by an overall smooth function
+
+	if (flat2d != "") {
+	    if (!access (flat2d // mstype)) {
+		print ("Extract flat field ", flat2d) | tee (log1)
+		if (flat2d != apref)
+		    apscript (flat2d, output=resp, background="none",
+			clean=clean, extras=no)
+		else
+		    apscript (flat2d, output=resp, ansrecenter="NO",
+			ansrecenter="NO", ansresize="NO", ansedit="NO",
+			anstrace="NO", background="none",
+			clean=clean, extras=no)
+	    } else
+		imcopy (flat2d//".ms", resp, verbose=no)
+
+	    if (fitflat) {
+		print ("Fit and ratio flat field ", flat2d) | tee (log1)
+		blkavg (resp, temp1, option="average", b1=1, b2=10000)
+		imcopy (temp1//"[*,1]", temp1, verbose=no)
+		fit1d (temp1, temp1, "fit", axis=1, interactive=interactive,
+		    sample="*", naverage=1, function=function, order=order,
+		    low_reject=0., high_reject=0., niterate=1, grow=0.,
+		    graphics="stdgraph")
+		sarith (resp, "/", temp1, resp, w1=INDEF, w2=INDEF,
+		    apertures="", bands="", beams="", apmodulus=0, reverse=no,
+		    ignoreaps=yes, format="multispec", renumber=no,
+		    offset=0, clobber=yes, merge=no, errval=0, verbose=no)
+		imdelete (temp1, verify=no)
+	    }
+	}
+
+	# If using a throughput image extract it if necessary.
+	# Apply it to the flat field if given and otherwise only
+	# compute the throughput through each aperture.
+
+	if (skyflat2d != "") {
+	    if (access (skyflat2d // imtype)) {
+		if (!access (skyflat2d // mstype)) {
+		    print ("Extract throughput image ", skyflat2d) | tee (log1)
+		    apscript (skyflat2d, output=temp1, background="none",
+			clean=clean, extras=no)
+		    temp2 = temp1
+		} else
+		    temp2 = skyflat2d // ".ms"
+
+		if (flat2d != "") {
+		    print ("Correct flat field to throughput image") |
+			tee (log1)
+		    sarith (temp2, "/", resp, temp1, w1=INDEF, w2=INDEF,
+			apertures="", bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=no, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=no,
+			errval=0, verbose=no)
+		    fit1d (temp1, temp1, type="fit", axis=1,
+			interactive=no, sample="*", naverage=1,
+			function="legendre", order=1, niterate=0)
+		    sarith (resp, "*", temp1, resp, w1=INDEF, w2=INDEF,
+			apertures="", bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=yes, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=no,
+			errval=0, verbose=no)
+		    imdelete (temp1, verify=no)
+		} else {
+		    print ("Compute aperture throughput from image") |
+			tee (log1)
+		    fit1d (temp2, resp, type="fit", axis=1,
+			interactive=no, sample="*", naverage=1,
+			function="legendre", order=1, niterate=0)
+		    if (temp2 == temp1)
+			imdelete (temp2, verify=no)
+		}
+
+	    # If a flat field and throughput file are used scale the average
+	    # flat field in each aperture to those values
+
+	    } else if (flat2d != "") {
+		print ("Correct flat field with throughput file ", skyflat2d) |
+		    tee (log1)
+		fit1d (resp, resp, type="ratio", axis=1,
+		    interactive=no, sample="*", naverage=1,
+		    function="legendre", order=1, niterate=0)
+
+		list = skyflat2d
+		while (fscan (list, ap, respval) != EOF) {
+		    sarith (resp, "*", respval, resp, w1=INDEF, w2=INDEF,
+			apertures=ap, bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=no, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=yes,
+			errval=0., verbose=no, >& "dev$null")
+		}
+		list = ""
+
+	    # If only a throughput file is given create the response from the
+	    # aperture reference and set the aperture response to the specified
+	    # values.
+
+	    } else {
+		print ("Set aperture throughput using ", skyflat2d) | tee (log1)
+		if (!access (apref // mstype)) {
+		    apscript (apref, output=resp, ansrecenter="NO",
+			ansrecenter="NO", ansresize="NO", ansedit="NO",
+			anstrace="NO", background="none",
+			clean=no, extras=no)
+		    sarith (resp, "replace", "0", resp, w1=INDEF,
+			w2=INDEF, apertures="", bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=no, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=yes,
+			errval=0., verbose=no)
+		} else
+		    sarith (apref//".ms", "replace", "0", resp, w1=INDEF,
+			w2=INDEF, apertures="", bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=no, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=yes,
+			errval=0., verbose=no)
+
+		list = skyflat2d
+		while (fscan (list, ap, respval) != EOF) {
+		    sarith (resp, "replace", respval, resp, w1=INDEF, w2=INDEF,
+			apertures=ap, bands="", beams="", apmodulus=0,
+			reverse=no, ignoreaps=no, format="multispec",
+			renumber=no, offset=0, clobber=yes, merge=yes,
+			errval=0., verbose=no)
+		}
+		list = ""
+	    }
+	}
+
+	# The final response is normalized to overall unit mean and the
+	# average aperture response is printed.
+
+	print ("Create the normalized response ", resp) | tee (log1)
+	bscale (resp, resp, bzero="0.", bscale="mean", section="",
+	    step=1, upper=INDEF, lower=INDEF, verbose=yes) | tee (log1, >log2)
+	blkavg (resp, temp1, option="average", b1=10000, b2=1)
+	print ("Average aperture response:") | tee (log1, >log2)
+	naxis = 5
+	#hselect (temp1, "naxis", yes) | scan (naxis)
+	#if (naxis == 1)
+	#    listpixels (temp1) | tee (log1, >log2)
+	#else
+	#    listpixels (temp1//"[1,*]") | tee (log1, >log2)
+
+	hselect (temp1, "naxis", yes, > temp2)
+	list = temp2; ap = fscan (list, naxis)
+	if (naxis == 1)
+	    listpixels (temp1) | tee (log1, >log2)
+	else
+	    listpixels (temp1//"[1,*]") | tee (log1, >log2)
+	list = ""; delete (temp2, verify=no)
+
+	imdelete (temp1, verify=no)
+end