Repatch (from linux) of OSX IRAF

1 files changed, 209 insertions, 0 deletions

diff --git a/noao/onedspec/getextn.x b/noao/onedspec/getextn.x
new file mode 100644
index 00000000..82640152
--- /dev/null
+++ b/noao/onedspec/getextn.x
@@ -0,0 +1,209 @@
+include	<error.h>
+include	<syserr.h>
+
+define	EXTN_LOOKUP	10	# Interp index for de-extinction
+define	DEXTN_LOOKUP	11	# Interp index for differential extn table
+define	TEMP_SPACE	100	# Amount of temporary space to allocate
+
+# GET_EXTN -- Get extinction from calibration file and
+#             any update as indicated from the SENSITIVITY
+#             computation
+
+procedure get_extn (wave_tbl, extn_tbl, nwaves)
+
+pointer	wave_tbl, extn_tbl
+int	nwaves
+
+pointer	waves, extns
+
+begin
+	# Get standard extinction values
+	call ext_load (waves, extns, nwaves)
+
+	# Copy values to external pointer.
+	# Use of salloc is incorrect but this is a hack on old code. FV
+	call salloc (extn_tbl, nwaves, TY_REAL)
+	call salloc (wave_tbl, nwaves, TY_REAL)
+	call amovr (Memr[waves], Memr[wave_tbl], nwaves)
+	call amovr (Memr[extns], Memr[extn_tbl], nwaves)
+	call mfree (waves, TY_REAL)
+	call mfree (extns, TY_REAL)
+end
+
+
+# DE_EXT_SPEC -- Apply extinction correction to a spectrum
+
+procedure de_ext_spec (spec, airm, w0, wpc, wave_tbl, extn_tbl, nwaves, len)
+
+real	spec[ARB], wave_tbl[ARB], extn_tbl[ARB]
+real	airm, w0, wpc
+int	nwaves, len
+
+int	i, ierr
+real	wave, ext
+bool	lin_log
+
+begin
+	# Assume linear dispersion, but possibly in LOG10
+	if (w0 < 5.0 && wpc < 0.05)
+	    lin_log = true
+	else
+	    lin_log = false
+
+	# Initialize interpolator
+	call intrp0 (EXTN_LOOKUP)
+
+	do i = 1, len {
+	    wave = w0 + (i-1) * wpc
+	    if (lin_log)
+		wave = 10.0 ** wave
+
+	    # Table must be in wavelength, not log[]
+	    call intrp (EXTN_LOOKUP, wave_tbl, extn_tbl, nwaves, 
+		wave, ext, ierr)
+
+	    spec[i] = spec[i] * 10.0 ** (0.4 * airm * ext)
+	}
+end
+
+# SUM_SPEC -- Add up counts within a specified region of a spectrum
+#             denoted by a wavelength range.
+#             The summation is active only over those pixels which
+#             are completely within the range specification.
+#             Data referenced outside the spectrum is ignored.
+
+procedure sum_spec (spec, w1, w2, w0, wpc, counts, len)
+
+real	spec[ARB], w1, w2, w0, wpc, counts
+int	len
+
+int	i, pix1, pix2
+
+real	pix_index()
+
+begin
+	# Compute pixel numbers from w1 to w2
+	pix1 = max (int (pix_index (w0, wpc, w1) + 1.0), 1)
+	pix2 = max (int (pix_index (w0, wpc, w2)      ), pix1)
+	pix2 = min (pix2, len)
+
+	counts = 0.0
+
+	do i = pix1, pix2
+	    counts = counts + spec[i]
+
+	# Guarantee that there are no negative counts
+	if (counts < 0.0)
+	    counts = 0.0
+end
+
+# PIX_INDEX -- Returns the pixel index at wavelength for linearly
+#              dispersion corrected spectra
+#
+#              The "Guess" is made that if the start wavelength for the
+#              spectrum is less than 5.0 and the dispersion is less than
+#              0.05, the spectrum has been linearized in LOG10 space.
+#
+# Note that in IRAF, a pixel index effectively refers to the center of a pixel.
+# So a spectrum must actually extend from w0-0.5*wpc to w0+(len+0.5)*wpc
+
+real procedure pix_index (w0, wpc, w)
+
+real	w0, wpc, w
+real	xw
+
+begin
+	# Check for LOG10 dispersion
+
+	if (w0 < 5.0 && wpc < 0.05)
+	    xw = log10 (w)
+	else
+	    xw = w
+
+	pix_index = (xw - w0) / wpc + 1.0
+end
+
+
+define	NALLOC	128	# Allocation block size
+
+# EXT_LOAD -- Read extinction data from database directory.
+
+procedure ext_load (waves, extns, nwaves)
+
+pointer	waves, extns
+int	nwaves
+
+real	wave, extn
+int	fd, nalloc
+pointer	sp, file
+
+int	open(), fscan(), nscan(), errcode()
+
+begin
+	call smark (sp)
+	call salloc (file, SZ_FNAME, TY_CHAR)
+
+	# Get the extinction file and open it.
+	call clgstr ("extinction", Memc[file], SZ_FNAME)
+	iferr (fd = open (Memc[file], READ_ONLY, TEXT_FILE)) {
+	    switch (errcode()) {
+	    case SYS_FNOFNAME:
+		nwaves = 2
+		call malloc (waves, nwaves, TY_REAL)
+		call malloc (extns, nwaves, TY_REAL)
+		Memr[waves] = 1000.
+		Memr[extns] = 0.
+		Memr[waves+1] = 10000.
+		Memr[extns+1] = 0.
+		call eprintf ("No extinction correction applied\n")
+		return
+	    default:
+		call erract (EA_ERROR)
+	    }
+	}
+
+	# Read the extinction data.
+	nalloc = 0
+	nwaves = 0
+	while (fscan (fd) != EOF) {
+	    call gargr (wave)
+	    call gargr (extn)
+	    if (nscan() != 2)
+		next
+
+	    if (nalloc == 0) {
+		nalloc = nalloc + NALLOC
+		call malloc (waves, nalloc, TY_REAL)
+		call malloc (extns, nalloc, TY_REAL)
+	    } else if (nwaves == nalloc) {
+		nalloc = nalloc + NALLOC
+		call realloc (waves, nalloc, TY_REAL)
+		call realloc (extns, nalloc, TY_REAL)
+	    }
+
+	    Memr[waves+nwaves] = wave
+	    Memr[extns+nwaves] = extn
+	    nwaves = nwaves + 1
+	}
+	call close (fd)
+
+	if (nwaves == 0)
+	    call error (1, "No extinction data found")
+
+	call realloc (waves, nwaves, TY_REAL)
+	call realloc (extns, nwaves, TY_REAL)
+
+	call sfree (sp)
+end
+
+
+# EXT_FREE -- Free extinction data arrays.
+
+procedure ext_free (waves, extns)
+
+pointer	waves, extns
+
+begin
+	call mfree (waves, TY_REAL)
+	call mfree (extns, TY_REAL)
+end