1 files changed, 376 insertions, 0 deletions
diff --git a/noao/imred/ccdred/src/icscale.x b/noao/imred/ccdred/src/icscale.x
new file mode 100644
index 00000000..fc4efb2f
--- /dev/null
+++ b/noao/imred/ccdred/src/icscale.x
@@ -0,0 +1,376 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<imset.h>
+include	<error.h>
+include	"icombine.h"
+
+# IC_SCALE -- Get the scale factors for the images.
+# 1. This procedure does CLIO to determine the type of scaling desired.
+# 2. The output header parameters for exposure time and NCOMBINE are set.
+
+procedure ic_scale (in, out, offsets, scales, zeros, wts, nimages)
+
+pointer	in[nimages]		# Input images
+pointer	out[ARB]		# Output images
+int	offsets[nimages,ARB]	# Image offsets
+real	scales[nimages]		# Scale factors
+real	zeros[nimages]		# Zero or sky levels
+real	wts[nimages]		# Weights
+int	nimages			# Number of images
+
+int	stype, ztype, wtype
+int	i, j, k, l, nout
+real	mode, median, mean, exposure, zmean, darktime, dark
+pointer	sp, ncombine, exptime, modes, medians, means
+pointer	section, str, sname, zname, wname, imref
+bool	domode, domedian, domean, dozero, snorm, znorm, wflag
+
+bool	clgetb()
+int	hdmgeti(), strdic(), ic_gscale()
+real	hdmgetr(), asumr(), asumi()
+errchk	ic_gscale, ic_statr
+
+include	"icombine.com"
+
+begin
+	call smark (sp)
+	call salloc (ncombine, nimages, TY_INT)
+	call salloc (exptime, nimages, TY_REAL)
+	call salloc (modes, nimages, TY_REAL)
+	call salloc (medians, nimages, TY_REAL)
+	call salloc (means, nimages, TY_REAL)
+	call salloc (section, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (sname, SZ_FNAME, TY_CHAR)
+	call salloc (zname, SZ_FNAME, TY_CHAR)
+	call salloc (wname, SZ_FNAME, TY_CHAR)
+
+	# Set the defaults.
+	call amovki (1, Memi[ncombine], nimages)
+	call amovkr (0., Memr[exptime], nimages)
+	call amovkr (INDEF, Memr[modes], nimages)
+	call amovkr (INDEF, Memr[medians], nimages)
+	call amovkr (INDEF, Memr[means], nimages)
+	call amovkr (1., scales, nimages)
+	call amovkr (0., zeros, nimages)
+	call amovkr (1., wts, nimages)
+
+	# Get the number of images previously combined and the exposure times.
+	# The default combine number is 1 and the default exposure is 0.
+
+	do i = 1, nimages {
+	    iferr (Memi[ncombine+i-1] = hdmgeti (in[i], "ncombine"))
+		Memi[ncombine+i-1] = 1
+	    iferr (Memr[exptime+i-1] = hdmgetr (in[i], "exptime"))
+		Memr[exptime+i-1] = 0.
+	    if (project) {
+		call amovki (Memi[ncombine], Memi[ncombine], nimages)
+		call amovkr (Memr[exptime], Memr[exptime], nimages)
+		break
+	    }
+	}
+
+	# Set scaling factors.
+
+	stype = ic_gscale ("scale", Memc[sname], STYPES, in, Memr[exptime],
+	    scales, nimages)
+	ztype = ic_gscale ("zero", Memc[zname], ZTYPES, in, Memr[exptime],
+	    zeros, nimages)
+	wtype = ic_gscale ("weight", Memc[wname], WTYPES, in, Memr[exptime],
+	    wts, nimages)
+
+	# Get image statistics only if needed.
+	domode = ((stype==S_MODE)||(ztype==S_MODE)||(wtype==S_MODE))
+	domedian = ((stype==S_MEDIAN)||(ztype==S_MEDIAN)||(wtype==S_MEDIAN))
+	domean = ((stype==S_MEAN)||(ztype==S_MEAN)||(wtype==S_MEAN))
+	if (domode || domedian || domean) {
+	    Memc[section] = EOS
+	    Memc[str] = EOS
+	    call clgstr ("statsec", Memc[section], SZ_FNAME)
+	    call sscan (Memc[section])
+	    call gargwrd (Memc[section], SZ_FNAME)
+	    call  gargwrd (Memc[str], SZ_LINE)
+
+	    i = strdic (Memc[section], Memc[section], SZ_FNAME, S_SECTION)
+	    switch (i) {
+	    case S_INPUT:
+		call strcpy (Memc[str], Memc[section], SZ_FNAME)
+		imref = NULL
+	    case S_OUTPUT:
+		call strcpy (Memc[str], Memc[section], SZ_FNAME)
+		imref = out[1]
+	    case S_OVERLAP:
+		call strcpy ("[", Memc[section], SZ_FNAME)
+		do i = 1, IM_NDIM(out[1]) {
+		    k = offsets[1,i] + 1
+		    l = offsets[1,i] + IM_LEN(in[1],i)
+		    do j = 2, nimages {
+			k = max (k, offsets[j,i]+1)
+			l = min (l, offsets[j,i]+IM_LEN(in[j],i))
+		    }
+		    if (i < IM_NDIM(out[1]))
+			call sprintf (Memc[str], SZ_LINE, "%d:%d,")
+		    else
+			call sprintf (Memc[str], SZ_LINE, "%d:%d]")
+			    call pargi (k)
+			    call pargi (l)
+		    call strcat (Memc[str], Memc[section], SZ_FNAME)
+		}
+		imref = out[1]
+	    default:
+		imref = NULL
+	    }
+
+	    do i = 1, nimages {
+		if (imref != out[1])
+		    imref = in[i]
+		call ic_statr (in[i], imref, Memc[section], offsets,
+		    i, nimages, domode, domedian, domean, mode, median, mean)
+		if (domode) {
+		    Memr[modes+i-1] = mode
+		    if (stype == S_MODE)
+			scales[i] = mode
+		    if (ztype == S_MODE)
+			zeros[i] = mode
+		    if (wtype == S_MODE)
+			wts[i] = mode
+		}
+		if (domedian) {
+		    Memr[medians+i-1] = median
+		    if (stype == S_MEDIAN)
+			scales[i] = median
+		    if (ztype == S_MEDIAN)
+			zeros[i] = median
+		    if (wtype == S_MEDIAN)
+			wts[i] = median
+		}
+		if (domean) {
+		    Memr[means+i-1] = mean
+		    if (stype == S_MEAN)
+			scales[i] = mean
+		    if (ztype == S_MEAN)
+			zeros[i] = mean
+		    if (wtype == S_MEAN)
+			wts[i] = mean
+		}
+	    }
+	}
+
+	do i = 1, nimages
+	    if (scales[i] <= 0.) {
+		call eprintf ("WARNING: Negative scale factors")
+		call eprintf (" -- ignoring scaling\n")
+		call amovkr (1., scales, nimages)
+		break
+	    }
+
+	# Convert to relative factors if needed.
+	snorm = (stype == S_FILE || stype == S_KEYWORD)
+	znorm = (ztype == S_FILE || ztype == S_KEYWORD)
+	wflag = (wtype == S_FILE || wtype == S_KEYWORD)
+	if (snorm)
+	    call arcpr (1., scales, scales, nimages)
+	else {
+	    mean = asumr (scales, nimages) / nimages
+	    call adivkr (scales, mean, scales, nimages)
+	}
+	call adivr (zeros, scales, zeros, nimages)
+	zmean = asumr (zeros, nimages) / nimages
+
+	if (wtype != S_NONE) {
+	    do i = 1, nimages {
+		if (wts[i] <= 0.) {
+		    call eprintf ("WARNING: Negative weights")
+		    call eprintf (" -- using only NCOMBINE weights\n")
+		    do j = 1, nimages
+			wts[j] = Memi[ncombine+j-1]
+		    break
+		}
+		if (ztype == S_NONE || znorm || wflag)
+	            wts[i] = Memi[ncombine+i-1] * wts[i]
+		else {
+		    if (zeros[i] <= 0.) {
+			call eprintf ("WARNING: Negative zero offsets")
+			call eprintf (" -- ignoring zero weight adjustments\n")
+			do j = 1, nimages
+			    wts[j] = Memi[ncombine+j-1] * wts[j]
+			break
+		    }
+		    wts[i] = Memi[ncombine+i-1] * wts[i] * zmean / zeros[i]
+		}
+	    }
+	}
+
+	if (znorm)
+	    call anegr (zeros, zeros, nimages)
+	else {
+	    # Because of finite arithmetic it is possible for the zero offsets
+	    # to be nonzero even when they are all equal.  Just for the sake of
+	    # a nice log set the zero offsets in this case.
+
+	    call asubkr (zeros, zmean, zeros, nimages)
+	    for (i=2; (i<=nimages)&&(zeros[i]==zeros[1]); i=i+1)
+		;
+	    if (i > nimages)
+		call aclrr (zeros, nimages)
+	}
+	mean = asumr (wts, nimages)
+	call adivkr (wts, mean, wts, nimages)
+
+	# Set flags for scaling, zero offsets, sigma scaling, weights.
+	# Sigma scaling may be suppressed if the scales or zeros are
+	# different by a specified tolerance.
+
+	doscale = false
+	dozero = false
+	doscale1 = false
+	dowts = false
+	do i = 2, nimages {
+	    if (snorm || scales[i] != scales[1])
+		doscale = true
+	    if (znorm || zeros[i] != zeros[1])
+		dozero = true
+	    if (wts[i] != wts[1])
+		dowts = true
+	}
+	if (doscale && sigscale != 0.) {
+	    do i = 1, nimages {
+		if (abs (scales[i] - 1) > sigscale) {
+		    doscale1 = true
+		    break
+		}
+	    }
+	    if (!doscale1 && zmean > 0.) {
+		do i = 1, nimages {
+		    if (abs (zeros[i] / zmean) > sigscale) {
+			doscale1 = true
+			break
+		    }
+		}
+	    }
+	}
+		    
+	# Set the output header parameters.
+	nout = asumi (Memi[ncombine], nimages)
+	call hdmputi (out[1], "ncombine", nout)
+	exposure = 0.
+	darktime = 0.
+	mean = 0.
+	do i = 1, nimages {
+	    exposure = exposure + wts[i] * Memr[exptime+i-1] / scales[i]
+	    ifnoerr (dark = hdmgetr (in[i], "darktime"))
+		darktime = darktime + wts[i] * dark / scales[i]
+	    else
+		darktime = darktime + wts[i] * Memr[exptime+i-1] / scales[i]
+	    ifnoerr (mode = hdmgetr (in[i], "ccdmean"))
+		mean = mean + wts[i] * mode / scales[i]
+	}
+	call hdmputr (out[1], "exptime", exposure)
+	call hdmputr (out[1], "darktime", darktime)
+	ifnoerr (mode = hdmgetr (out[1], "ccdmean")) {
+	    call hdmputr (out[1], "ccdmean", mean)
+	    iferr (call imdelf (out[1], "ccdmeant"))
+		;
+	}
+	if (out[2] != NULL) {
+	    call imstats (out[2], IM_IMAGENAME, Memc[str], SZ_FNAME)
+	    call imastr (out[1], "BPM", Memc[str])
+	}
+
+	# Start the log here since much of the info is only available here.
+	if (clgetb ("verbose")) {
+	    i = logfd
+	    logfd = STDOUT
+	    call ic_log (in, out, Memi[ncombine], Memr[exptime], Memc[sname],
+		Memc[zname], Memc[wname], Memr[modes], Memr[medians],
+		Memr[means], scales, zeros, wts, offsets, nimages, dozero,
+		nout, "", exposure)
+
+	    logfd = i
+	}
+	call ic_log (in, out, Memi[ncombine], Memr[exptime], Memc[sname],
+	    Memc[zname], Memc[wname], Memr[modes], Memr[medians], Memr[means],
+	    scales, zeros, wts, offsets, nimages, dozero, nout,
+	    "", exposure)
+
+	doscale = (doscale || dozero)
+
+	call sfree (sp)
+end
+
+
+# IC_GSCALE -- Get scale values as directed by CL parameter
+# The values can be one of those in the dictionary, from a file specified
+# with a @ prefix, or from an image header keyword specified by a ! prefix.
+
+int procedure ic_gscale (param, name, dic, in, exptime, values, nimages)
+
+char	param[ARB]		#I CL parameter name
+char	name[SZ_FNAME]		#O Parameter value
+char	dic[ARB]		#I Dictionary string
+pointer	in[nimages]		#I IMIO pointers
+real	exptime[nimages]	#I Exposure times
+real	values[nimages]		#O Values
+int	nimages			#I Number of images
+
+int	type			#O Type of value
+
+int	fd, i, nowhite(), open(), fscan(), nscan(), strdic()
+real	rval, hdmgetr()
+pointer	errstr
+errchk	open, hdmgetr()
+
+include	"icombine.com"
+
+begin
+	call clgstr (param, name, SZ_FNAME)
+	if (nowhite (name, name, SZ_FNAME) == 0)
+	    type = S_NONE
+	else if (name[1] == '@') {
+	    type = S_FILE
+	    fd = open (name[2], READ_ONLY, TEXT_FILE)
+	    i = 0
+	    while (fscan (fd) != EOF) {
+		call gargr (rval)
+		if (nscan() != 1)
+		    next
+		if (i == nimages) {
+		   call eprintf (
+		       "Warning: Ignoring additional %s values in %s\n")
+		       call pargstr (param)
+		       call pargstr (name[2])
+		   break
+		}
+		i = i + 1
+		values[i] = rval
+	    }
+	    call close (fd)
+	    if (i < nimages) {
+		call salloc (errstr, SZ_LINE, TY_CHAR)
+		call sprintf (Memc[errstr], SZ_FNAME,
+		    "Insufficient %s values in %s")
+		    call pargstr (param)
+		    call pargstr (name[2])
+		call error (1, Memc[errstr])
+	    }
+	} else if (name[1] == '!') {
+	    type = S_KEYWORD
+	    do i = 1, nimages {
+		values[i] = hdmgetr (in[i], name[2])
+		if (project) {
+		    call amovkr (values, values, nimages)
+		    break
+		}
+	    }
+	} else {
+	    type = strdic (name, name, SZ_FNAME, dic)
+	    if (type == 0)
+		call error (1, "Unknown scale, zero, or weight type")
+	    if (type==S_EXPOSURE)
+		do i = 1, nimages
+		    values[i] = max (0.001, exptime[i])
+	}
+
+	return (type)
+end