Initial commit

1 files changed, 266 insertions, 0 deletions

diff --git a/noao/onedspec/sensfunc/sfstds.x b/noao/onedspec/sensfunc/sfstds.x
new file mode 100644
index 00000000..07219729
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfstds.x
@@ -0,0 +1,266 @@
+include "sensfunc.h"
+
+
+# SF_STDS -- Get the standard observations for the specified apertures.
+# If ignoring aperture set all apertures to 1.
+# This routine knows the output of the task STANDARD.
+
+procedure sf_stds (standards, aps, ignoreaps, stds, nstds)
+
+char	standards	# Standard star data file
+pointer	aps		# Aperture list
+bool	ignoreaps	# Ignore apertures?
+pointer	stds		# Pointer to standard observations (returned)
+int	nstds		# Number of standard observations (returned)
+
+int	i, j, fd, beam, npts, nwaves, nalloc
+real	exptime, airmass, wstart, wend
+real	wavelength, flux, dwave, count
+pointer	sp, image, title, std
+pointer	waves, fluxes, dwaves, counts, sens, fit, wts, iwts, x, y
+
+bool	rng_elementi()
+int	open(), fscan(), nscan(), stridxs()
+errchk	open, malloc, realloc
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_STDIMAGE, TY_CHAR)
+	call salloc (title, SZ_STDTITLE, TY_CHAR)
+
+	# Open the standard observation data file.
+	fd = open (standards, READ_ONLY, TEXT_FILE)
+
+	# Read the standard observations and create a structure for each one.
+	# The beginning of a new star is found by a line whose first word
+	# begins with the character '['.  Otherwise the line is interpreted
+	# as a data line.  All unrecognized formats are skipped.
+
+	nwaves = 0
+	nstds = 0
+	while (fscan (fd) != EOF) {
+	    call gargwrd (Memc[image], SZ_STDIMAGE)
+	    if (Memc[image] == '[') {
+	        call gargi (beam)
+	        call gargi (npts)
+		call gargr (exptime)
+	        call gargr (airmass)
+	        call gargr (wstart)
+	        call gargr (wend)
+	        call gargstr (Memc[title], SZ_STDTITLE)
+	        if (nscan() < 7)
+		    next
+		if (!rng_elementi (aps, beam))
+		    next
+		if (IS_INDEF (exptime) || exptime <= 0.) {
+		    call eprintf (
+		"%s: Warning - exposure time missing or zero, using 1 second\n")
+			call pargstr (Memc[image])
+		    exptime = 1.
+		}
+
+		# For the first one create the pointer to the array of
+		# structures.  For the following stars increase the size
+		# of the pointer array and finish up the previous standard
+		# star.
+
+	        if (nstds == 0) {
+		    nstds = nstds + 1
+		    call calloc (stds, nstds, TY_INT)
+		    call calloc (std, LEN_STD, TY_STRUCT)
+		    Memi[stds+nstds-1] = std
+	        } else {
+	    	    if (nwaves > 0) {
+			call realloc (waves, nwaves, TY_REAL)
+			call realloc (fluxes, nwaves, TY_REAL)
+			call realloc (dwaves, nwaves, TY_REAL)
+			call realloc (counts, nwaves, TY_REAL)
+			call realloc (wts, nwaves, TY_REAL)
+			call malloc (sens, nwaves, TY_REAL)
+			call malloc (fit, nwaves, TY_REAL)
+			call malloc (iwts, nwaves, TY_REAL)
+			call malloc (x, nwaves, TY_REAL)
+			call malloc (y, nwaves, TY_REAL)
+			call amovr (Memr[wts], Memr[iwts], nwaves)
+			STD_NWAVES(std) = nwaves
+			STD_WAVES(std) = waves
+			STD_FLUXES(std) = fluxes
+			STD_DWAVES(std) = dwaves
+			STD_COUNTS(std) = counts
+			STD_SENS(std) = sens
+			STD_FIT(std) = fit
+			STD_WTS(std) = wts
+			STD_IWTS(std) = iwts
+			STD_X(std) = x
+			STD_Y(std) = y
+
+			nstds = nstds + 1
+			call realloc (stds, nstds, TY_INT)
+			call calloc (std, LEN_STD, TY_STRUCT)
+			Memi[stds+nstds-1] = std
+		    }
+	        }
+
+		# Start a new standard star.
+	        std = Memi[stds+nstds-1]
+		if (ignoreaps)
+	            STD_BEAM(std) = 1
+		else
+	            STD_BEAM(std) = beam
+		STD_NPTS(std) = npts
+	        STD_EXPTIME(std) = exptime
+	        STD_AIRMASS(std) = airmass
+		STD_WSTART(std) = wstart
+		STD_WEND(std) = wend
+		STD_SHIFT(std) = 0.
+		STD_NWAVES(std) = 0
+
+		# Decode the image and sky strings.
+	        call strcpy (Memc[title], STD_TITLE(std), SZ_STDTITLE)
+		i = stridxs ("]", Memc[image])
+		if (Memc[image+i] == ']')
+		    i = i + 1
+		Memc[image+i-1] = EOS
+	        call strcpy (Memc[image+1], STD_IMAGE(std), SZ_STDIMAGE)
+		if (Memc[image+i] == '-') {
+		    i = i + 2
+		    j = stridxs ("]", Memc[image+i]) + i
+		    Memc[image+j-1] = EOS
+	            call strcpy (Memc[image+i], STD_SKY(std), SZ_STDIMAGE)
+		} else
+		    STD_SKY(std) = EOS
+	        nwaves = 0
+
+	    # Interprete the line as standard star wavelength point.
+	    } else if (nstds > 0) {
+		call reset_scan()
+		call gargr (wavelength)
+		call gargr (flux)
+		call gargr (dwave)
+		call gargr (count)
+		if (nscan() < 3)
+		    next
+		if (wavelength < min (wstart, wend) ||
+		    wavelength > max (wstart, wend) ||
+		    flux<=0. || dwave<=0. || count<=0.)
+		    next
+		if (!rng_elementi (aps, beam))
+		    next
+		nwaves = nwaves + 1
+
+		# Allocate in blocks to minimize the number of reallocations.
+		if (nwaves == 1) {
+		    nalloc = 100
+		    call malloc (waves, nalloc, TY_REAL)
+		    call malloc (fluxes, nalloc, TY_REAL)
+		    call malloc (dwaves, nalloc, TY_REAL)
+		    call malloc (counts, nalloc, TY_REAL)
+		    call malloc (wts, nalloc, TY_REAL)
+		} else if (nwaves > nalloc) {
+		    nalloc = nalloc + 100
+		    call realloc (waves, nalloc, TY_REAL)
+		    call realloc (fluxes, nalloc, TY_REAL)
+		    call realloc (dwaves, nalloc, TY_REAL)
+		    call realloc (counts, nalloc, TY_REAL)
+		    call realloc (wts, nalloc, TY_REAL)
+		}
+
+		# Record the data and compute the sensitivity.
+		Memr[waves+nwaves-1] = wavelength
+		Memr[fluxes+nwaves-1] = flux
+		Memr[dwaves+nwaves-1] = dwave
+		Memr[counts+nwaves-1] = count
+		Memr[wts+nwaves-1] = 1.
+	    }
+	}
+
+	# Finish up the last standard star and close the file.
+	if (nstds > 0) {
+	    STD_NWAVES(std) = nwaves
+    	    if (nwaves > 0) {
+		call realloc (waves, nwaves, TY_REAL)
+		call realloc (fluxes, nwaves, TY_REAL)
+		call realloc (dwaves, nwaves, TY_REAL)
+		call realloc (counts, nwaves, TY_REAL)
+		call realloc (wts, nwaves, TY_REAL)
+		call malloc (sens, nwaves, TY_REAL)
+		call malloc (fit, nwaves, TY_REAL)
+		call malloc (iwts, nwaves, TY_REAL)
+		call malloc (x, nwaves, TY_REAL)
+		call malloc (y, nwaves, TY_REAL)
+		call amovr (Memr[wts], Memr[iwts], nwaves)
+		STD_WAVES(std) = waves
+		STD_FLUXES(std) = fluxes
+		STD_DWAVES(std) = dwaves
+		STD_COUNTS(std) = counts
+		STD_SENS(std) = sens
+		STD_FIT(std) = fit
+		STD_WTS(std) = wts
+		STD_IWTS(std) = iwts
+		STD_X(std) = x
+		STD_Y(std) = y
+	    }
+	}
+	call close (fd)
+	call sfree (sp)
+
+	# Add standard stars for any added points and composite points.
+	nstds = nstds + 2
+	call realloc (stds, nstds, TY_INT)
+	call calloc (std, LEN_STD, TY_STRUCT)
+	Memi[stds+nstds-2] = std
+	call strcpy ("Added", STD_IMAGE(std), SZ_STDIMAGE)
+	STD_BEAM(std) = STD_BEAM(Memi[stds])
+	STD_NPTS(std) = STD_NPTS(Memi[stds])
+	STD_EXPTIME(std) = 1.
+	STD_AIRMASS(std) = 1.
+	STD_WSTART(std) = STD_WSTART(Memi[stds])
+	STD_WEND(std) = STD_WEND(Memi[stds])
+	STD_SHIFT(std) = 0.
+	STD_NWAVES(std) = 0
+	call calloc (std, LEN_STD, TY_STRUCT)
+	Memi[stds+nstds-1] = std
+	call strcpy ("Composite", STD_IMAGE(std), SZ_STDIMAGE)
+	STD_BEAM(std) = STD_BEAM(Memi[stds])
+	STD_NPTS(std) = STD_NPTS(Memi[stds])
+	STD_EXPTIME(std) = 1.
+	STD_AIRMASS(std) = 1.
+	STD_WSTART(std) = STD_WSTART(Memi[stds])
+	STD_WEND(std) = STD_WEND(Memi[stds])
+	STD_SHIFT(std) = 0.
+	STD_NWAVES(std) = 0
+end
+
+
+# SF_FREE -- Free the standard observations and aperture array.
+
+procedure sf_free (stds, nstds, apertures, napertures)
+
+pointer	stds		# Pointer to standard observations
+int	nstds		# Number of standard observations
+pointer	apertures	# Pointer to apertures array
+int	napertures	# Number of apertures
+
+int	i
+pointer	std
+
+begin
+	do i = 1, nstds {
+	    std = Memi[stds+i-1]
+	    if (STD_NWAVES(std) > 0) {
+		call mfree (STD_WAVES(std), TY_REAL)
+		call mfree (STD_FLUXES(std), TY_REAL)
+		call mfree (STD_DWAVES(std), TY_REAL)
+		call mfree (STD_COUNTS(std), TY_REAL)
+		call mfree (STD_SENS(std), TY_REAL)
+		call mfree (STD_FIT(std), TY_REAL)
+		call mfree (STD_WTS(std), TY_REAL)
+		call mfree (STD_IWTS(std), TY_REAL)
+		call mfree (STD_X(std), TY_REAL)
+		call mfree (STD_Y(std), TY_REAL)
+	    }
+	    call mfree (std, TY_STRUCT)
+	}
+	call mfree (stds, TY_INT)
+	call mfree (apertures, TY_INT)
+end