Initial commit

1 files changed, 315 insertions, 0 deletions

diff --git a/noao/twodspec/longslit/response.x b/noao/twodspec/longslit/response.x
new file mode 100644
index 00000000..dd61ecc4
--- /dev/null
+++ b/noao/twodspec/longslit/response.x
@@ -0,0 +1,315 @@
+include	<imhdr.h>
+include	<pkg/gtools.h>
+include <pkg/xtanswer.h>
+
+# T_RESPONSE -- Determine the response function for 2D spectra.
+#
+# A calibration image is divided by a normalization spectrum to form
+# a response image.  The normalization spectrum is derived by averaging
+# the normalization image across dispersion.  The normalization spectrum
+# is then smoothed by curve fitting.  The smoothed normalization
+# spectrum is divided into the calibration image to form the response
+# function image.  The curve fitting may be performed interactively
+# using the icfit package.  A response function is determined for each
+# input image.  Image sections in the calibration image may be used to determine
+# the response for only part of an image such as with multiple slits.
+
+# CL callable task.
+#
+# The images are given by image templates.  The number of images must
+# in each list must match.  Image sections are allowed in the calibration
+# image.
+
+procedure t_response ()
+
+int	list1				# List of calibration images
+int	list2				# List of normalization images
+int	list3				# List of response images
+real	threshold			# Response threshold
+int	naverage			# Sample averaging size
+int	order				# Order of curve fitting function
+real	low_reject, high_reject		# Rejection thresholds
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection growing radius
+int	interactive			# Interactive?
+
+pointer	cal, norm, resp, ic, gt
+pointer	sp, image1, image2, image3, history
+
+int	clgeti(), imtopen(), imtgetim(), imtlen(), gt_init(), ic_geti()
+bool	clgetb()
+real	clgetr(), ic_getr()
+pointer	immap()
+
+errchk	immap, ls_immap
+
+begin
+	call smark (sp)
+	call salloc (image1, SZ_LINE, TY_CHAR)
+	call salloc (image2, SZ_LINE, TY_CHAR)
+	call salloc (image3, SZ_LINE, TY_CHAR)
+	call salloc (history, SZ_LINE, TY_CHAR)
+
+	# Get the calibration, normalization, and response image lists and
+	# check that the they match.
+
+	call clgstr ("calibration", Memc[image1], SZ_LINE)
+	call clgstr ("normalization", Memc[image2], SZ_LINE)
+	call clgstr ("response", Memc[image3], SZ_LINE)
+
+	list1 = imtopen (Memc[image1])
+	list2 = imtopen (Memc[image2])
+	list3 = imtopen (Memc[image3])
+	if ((imtlen(list1)!=imtlen(list3)) || (imtlen(list2)!=imtlen(list3))) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call imtclose (list3)
+	    call error (0,  "Image lists do not match")
+	}
+
+	# Get remaining parameters and initialize the curve fitting package.
+
+	threshold = clgetr ("threshold")
+	call clgstr ("sample", Memc[image1], SZ_LINE)
+	naverage = clgeti ("naverage")
+	call clgstr ("function", Memc[image2], SZ_LINE)
+	order = clgeti ("order")
+	low_reject = clgetr ("low_reject")
+	high_reject = clgetr ("high_reject")
+	niterate = clgeti ("niterate")
+	grow = clgetr ("grow")
+	if (clgetb ("interactive"))
+	    interactive = YES
+	else
+	    interactive = ALWAYSNO
+
+	# Set the ICFIT pointer structure.
+	call ic_open (ic)
+	call ic_pstr (ic, "sample", Memc[image1])
+	call ic_puti (ic, "naverage", naverage)
+	call ic_pstr (ic, "function", Memc[image2])
+	call ic_puti (ic, "order", order)
+	call ic_putr (ic, "low", low_reject)
+	call ic_putr (ic, "high", high_reject)
+	call ic_puti (ic, "niterate", niterate)
+	call ic_putr (ic, "grow", grow)
+	call ic_pstr (ic, "ylabel", "")
+
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Create the response image for each calibration image.
+
+	while ((imtgetim (list1, Memc[image1], SZ_LINE) != EOF) &&
+	    (imtgetim (list2, Memc[image2], SZ_LINE) != EOF) &&
+	    (imtgetim (list3, Memc[image3], SZ_LINE) != EOF)) {
+
+	    # Map the images.  If the response image does not exist it
+	    # is created and initialized to unit response everywhere.
+	    # If the calibration image is an image section then the response
+	    # image is opened as a section also.
+
+	    call ls_immap (Memc[image1], Memc[image3], cal, resp)
+	    norm = immap (Memc[image2], READ_ONLY, 0)
+
+	    # Determine whether the normalization spectrum is to be fit
+	    # interactively and if so set the graphics title.
+
+	    call sprintf (Memc[image2], SZ_LINE,
+		"Fit the normalization spectrum for %s interactively")
+	        call pargstr (Memc[image1])
+	    call xt_answer (Memc[image2], interactive)
+
+	    if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	        call sprintf (Memc[image2], SZ_LINE,
+		    "Fit the normalization spectrum for %s\n%s")
+	            call pargstr (Memc[image1])
+	            call pargstr (IM_TITLE(cal))
+	        call gt_sets (gt, GTTITLE, Memc[image2])
+	    }
+
+	    # Make the response.
+	    call re_make (cal, norm, resp, ic, gt, threshold, interactive)
+
+	    # Document the fit.
+	    call ic_gstr (ic, "sample", Memc[history], SZ_LINE)
+	    call clpstr ("sample", Memc[history])
+	    naverage = ic_geti (ic, "naverage")
+	    call clputi ("naverage", naverage)
+	    call ic_gstr (ic, "function", Memc[history], SZ_LINE)
+	    call clpstr ("function", Memc[history])
+	    order = ic_geti (ic, "order")
+	    call clputi ("order", order)
+	    low_reject = ic_getr (ic, "low")
+	    call clputr ("low_reject", low_reject)
+	    high_reject = ic_getr (ic, "high")
+	    call clputr ("high_reject", high_reject)
+	    niterate = ic_geti (ic, "niterate")
+	    call clputi ("niterate", niterate)
+	    grow = ic_getr (ic, "grow")
+	    call clputr ("grow", grow)
+
+	    call imaddr (resp, "ccdmean", 1.)
+	    call sprintf (Memc[history], SZ_LINE,
+		"Response determined from %s.")
+		call pargstr (Memc[image2])
+	    call xt_phistory (resp, Memc[history])
+	    call imunmap (cal)
+	    call imunmap (norm)
+	    call imunmap (resp)
+	}
+
+	# Finish up.
+
+	call ic_closer (ic)
+	call imtclose (list1)
+	call imtclose (list2)
+	call imtclose (list3)
+	call gt_free (gt)
+	call sfree (sp)
+end
+
+
+# RE_MAKE -- Given the calibration image determine the response.
+
+procedure re_make (cal, norm, resp, ic, gt, threshold, interactive)
+
+pointer	cal			# Calibration IMIO pointer
+pointer	norm			# Normalization IMIO pointer
+pointer	resp			# Response IMIO pointer
+pointer	ic			# ICFIT pointer
+pointer	gt			# GTOOLS pointer
+real	threshold		# Response threshold
+int	interactive		# Interactive?
+
+char	graphics[SZ_FNAME]	# Graphics output device
+int	laxis, paxis, npts
+pointer	cv, gp, sp, wavelengths, spectrum, wts
+
+pointer	gopen()
+errchk	get_daxis
+
+begin
+	# Determine the dispersion axis and set the axis labels.
+	call get_daxis (cal, laxis, paxis)
+
+	switch (laxis) {
+	case 1:
+	    call ic_pstr (ic, "xlabel", "Column")
+	case 2:
+	    call ic_pstr (ic, "xlabel", "Line")
+	}
+
+	# Get the normalization spectrum.
+
+	call ls_aimavg (norm, laxis, 1, IM_LEN(norm, 1), 1, IM_LEN(norm, 2),
+	    wavelengths, spectrum, npts)
+
+	# Allocate memory for the fit.
+
+	call smark (sp)
+	call salloc (wts, npts, TY_REAL)
+	call amovkr (1., Memr[wts], npts)
+
+	# Smooth the normalization spectrum.
+
+	call ic_putr (ic, "xmin", Memr[wavelengths])
+	call ic_putr (ic, "xmax", Memr[wavelengths+npts-1])
+
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    call clgstr ("graphics", graphics, SZ_FNAME)
+	    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+	    call icg_fit (ic, gp, "cursor", gt, cv, Memr[wavelengths],
+		Memr[spectrum], Memr[wts], npts)
+	    call gclose (gp)
+	} else {
+	    call ic_fit (ic, cv, Memr[wavelengths], Memr[spectrum], Memr[wts],
+		npts, YES, YES, YES, YES)
+	}
+
+	call cvvector (cv, Memr[wavelengths], Memr[spectrum], npts)
+	call cvfree (cv)
+
+	# Compute the response image by normalizing the calibration
+	# image by the normalization spectrum.
+
+	call re_normalize (cal, resp, laxis, threshold, Memr[spectrum], npts)
+
+	# Free allocated memory.
+
+	call sfree (sp)
+	call mfree (wavelengths, TY_REAL)
+	call mfree (spectrum, TY_REAL)
+end
+	
+
+# RE_NORMALIZE -- Divide each calibration image pixel by the normalization
+# spectrum at that pixel.
+
+procedure re_normalize (cal, resp, axis, threshold, spectrum, npts)
+
+pointer	cal			# Calibration IMIO pointer
+pointer	resp			# Response IMIO pointer
+int	axis			# Dispersion axis
+real	threshold		# Normalization treshold
+real	spectrum[npts]		# Pointer to normalization spectrum
+int	npts			# Number of points in spectrum
+
+int	i, j, ncols, nlines
+real	norm
+pointer	datain, dataout
+
+pointer	imgl2r(), impl2r()
+
+begin
+	ncols = IM_LEN (cal, 1)
+	nlines = IM_LEN (cal, 2)
+
+	# Compute the response image.
+	if (IS_INDEF (threshold)) {
+	    do i = 1, nlines {
+		datain = imgl2r (cal, i)
+		dataout = impl2r (resp, i)
+
+		switch (axis) {
+		case 1:
+		    call adivr (Memr[datain], spectrum, Memr[dataout], ncols)
+		case 2:
+		    call adivkr (Memr[datain], spectrum[i], Memr[dataout],
+			ncols)
+		}
+	    }
+	} else {
+	    do i = 1, nlines {
+		datain = imgl2r (cal, i)
+		dataout = impl2r (resp, i)
+
+		switch (axis) {
+		case 1:
+		    do j = 1, ncols {
+			norm = spectrum[j]
+			if (norm < threshold || Memr[datain] < threshold)
+			    Memr[dataout] = 1.
+			else
+			    Memr[dataout] = Memr[datain] / norm
+			datain = datain + 1
+			dataout = dataout + 1
+		    }
+		case 2:
+		    norm = spectrum[i]
+		    if (norm < threshold)
+			call amovkr (1., Memr[dataout], ncols)
+		    else {
+			do j = 1, ncols {
+			    if (Memr[datain] < threshold)
+				Memr[dataout] = 1.
+			    else
+				Memr[dataout] = Memr[datain] / norm
+			    datain = datain + 1
+			    dataout = dataout + 1
+			}
+		    }
+		}
+	    }
+	}
+end