Initial commit

1 files changed, 289 insertions, 0 deletions

diff --git a/noao/onedspec/sensfunc/sfgraph.x b/noao/onedspec/sensfunc/sfgraph.x
new file mode 100644
index 00000000..bb6fb26f
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfgraph.x
@@ -0,0 +1,289 @@
+include	<gset.h>
+include	<error.h>
+include	<math/curfit.h>
+include	"sensfunc.h"
+
+define	NCURVE		50		# Number of vectors in curve
+
+
+# SF_GRAPH -- Graph the desired data on the output graphics device.
+# This entry procedure determines the graph types, sets the device viewports
+# for each graph, and calls a procedure to make the graph.
+
+procedure sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of values in extinction table
+pointer	ecv			# Residual extinction curve.
+
+int	i, image, ngraphs, strlen()
+real	fa[8]
+pointer	sp, id, gio
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	# Clear the graphs, write the title, and set the viewports based on
+	# the number of graphs.
+
+	call smark (sp)
+	call salloc (id, SZ_LINE, TY_CHAR)
+	call sysid (Memc[id], SZ_LINE)
+
+	gio = GP_GIO(gp)
+	call gclear (gio)
+	call gseti (gio, G_FACOLOR, 0)
+	call gseti (gio, G_NMINOR, 0)
+	ngraphs = strlen (GP_GRAPHS(gp,1))
+	switch (ngraphs) {
+	case 1:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .1, .9) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	case 2:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .10, .97, .10, .45) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	case 3:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .10, .50, .1, .45) 
+	    call gseti (gio, G_WCS, 3)
+	    call gsview (gio, .57, .97, .1, .45) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	default:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .50, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .57, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 3)
+	    call gsview (gio, .10, .50, .1, .45) 
+	    call gseti (gio, G_WCS, 4)
+	    call gsview (gio, .57, .97, .1, .45) 
+	    GP_SZMARK(gp) = .01
+	    GP_SZMDEL(gp) = .01
+	}
+
+	# For each graph select the viewport and call a procedure to make
+	# the graph.
+
+	image = 0
+	for (i = 1; GP_GRAPHS(gp,i) != EOS; i = i + 1) {
+	    call gseti (gio, G_WCS, i)
+	    if (i > 1)
+		call gfill (gio, fa, fa[5], 4, GF_SOLID)
+	    switch (GP_GRAPHS(gp,i)) {
+	    case 'a', 's', 'r':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,i))
+		call sf_graph1 (gp, stds, nstds, GP_GRAPHS(gp,i))
+	    case 'e':
+		call sf_egraph (gp, wextn, extn, nextn, ecv)
+	    case 'c':
+		call sf_cgraph (gp, stds, nstds, cv)
+	    case 'i', 'l':
+		if (GP_GRAPHS(gp,i) == 'i')
+		    GP_LOG(gp) = NO
+		else
+		    GP_LOG(gp) = YES
+		image = image + 1
+		iferr (call sf_image (gp, image, stds, nstds, cv, wextn, extn,
+		    nextn, ecv))
+		    call erract (EA_WARN)
+	    }
+	}
+
+	call gseti (gio, G_WCS, 0)
+	call gtext (gio, 0.5, 1., Memc[id], "h=c,v=t,f=b")
+	call gtext (gio, 0.5, 0.97, GP_TITLE(gp), "h=c,v=t,f=b")
+
+	call sfree (sp)
+end
+
+
+# SF_GRAPH1 -- Make graph of sensitivity or residual sensitivity vs wavelength.
+
+procedure sf_graph1 (gp, stds, nstds, graph)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+char	graph			# Graph type
+
+int	i, j, n, mark, mdel, cdel, color
+real	szmark, szmdel, ymin, ymax, y1, y2
+pointer	x, y, w, gio
+
+begin
+	gio = GP_GIO(gp)
+
+	# Autoscale the included data in Y and set wavelength range.
+	j = 0
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    j = j + 1
+	    n = STD_NWAVES(stds[i])
+	    x = STD_X(stds[i])
+	    y = STD_Y(stds[i])
+	    if (j == 1)
+		call alimr (Memr[y], n, ymin, ymax)
+	    else {
+		call alimr (Memr[y], n, y1, y2)
+		ymin = min (ymin, y1)
+		ymax = max (ymax, y2)
+	    }
+	}
+	y1 = 0.05 * (ymax - ymin)
+	ymin = ymin - y1
+	ymax = ymax + y1
+
+	# Draw axes and title based on type of graph.
+	switch (graph) {
+	case 'a':
+	    call gswind (gio, GP_AIRMASS(gp,1), GP_AIRMASS(gp,2), ymin, ymax)
+	    call glabax (gio, "Sensitivity Residuals vs Airmass", "", "")
+	case 's':
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Sensitivity vs Wavelength", "", "")
+	case 'r':
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Sensitivity Residuals vs Wavelength", "", "")
+	}
+
+	# Plot the data with appropriate mark types and sizes.
+	mdel = GP_MDEL(gp)
+	cdel = GP_CDEL(gp)
+	szmdel = GP_SZMDEL(gp)
+	szmark = GP_SZMARK(gp)
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    if (i != nstds-1) {
+		mark = GP_MARK(gp)
+		color = GP_CMARK(gp)
+	    } else {
+		mark = GP_MADD(gp)
+		color = GP_CADD(gp)
+	    }
+	    n = STD_NWAVES(stds[i])
+	    x = STD_X(stds[i]) - 1
+	    y = STD_Y(stds[i]) - 1
+	    w = STD_WTS(stds[i]) - 1
+	    do j = 1, n {
+		if (Memr[w+j] == 0.) {
+		    call gseti (gio, G_PLCOLOR, cdel)
+		    call gmark (gio, Memr[x+j], Memr[y+j], mdel, szmdel, szmdel)
+		} else {
+		    call gseti (gio, G_PLCOLOR, color)
+	            call gmark (gio, Memr[x+j], Memr[y+j], mark, szmark, szmark)
+		}
+	    }
+	}
+end
+
+
+# SF_EGRAPH -- Graph extinction curves with and without residual extinction
+# correction.
+
+procedure sf_egraph (gp, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+pointer	ecv			# Residual extinction curve
+
+int	i, j
+real	xmin, xmax, dx, x, cveval()
+pointer	sp, ext, extnew, gio
+
+begin
+	call smark (sp)
+	call salloc (ext, NCURVE, TY_REAL)
+
+	# Interpolate extinction table to a grid of wavelengths within
+	# the range of the sensitivity data.
+
+	gio = GP_GIO(gp)
+	xmin = GP_WSTART(gp)
+	xmax = GP_WEND(gp)
+	dx = (xmax - xmin) / (NCURVE - 1)
+	x = xmin
+	do i = 0, NCURVE-1 {
+	    call intrp (1, wextn, extn, nextn, x, Memr[ext+i], j)
+	    x = x + dx
+	}
+	call gswind (gio, xmin, xmax, INDEF, INDEF)
+	call gascale (gio, Memr[ext], NCURVE, 2)
+
+	# If there is a residual extinction curve determine a new extinction
+	# curve.
+
+	if (ecv != NULL) {
+	    call salloc (extnew, NCURVE, TY_REAL)
+	    call amovr (Memr[ext], Memr[extnew], NCURVE)
+	    x = xmin
+	    do i = 0, NCURVE-1 {
+	        Memr[extnew+i] = Memr[extnew+i] + cveval (ecv, x)
+		x = x + dx
+	    }
+	    call grscale (gio, Memr[extnew], NCURVE, 2)
+	}
+
+	# Draw the axes and title and extinction curves.
+	call glabax (gio, "Extinction vs Wavelength", "", "")
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+	call gvline (gio, Memr[ext], NCURVE, xmin, xmax)
+	if (ecv != NULL) {
+	    call gseti (gio, G_PLTYPE, 2)
+	    call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)+1)
+	    call gvline (gio, Memr[extnew], NCURVE, xmin, xmax)
+	    call gseti (gio, G_PLTYPE, 1)
+	}
+
+	call sfree (sp)
+end
+
+
+# SF_FITGRAPH -- Overplot the fitted sensitivity curve.
+
+procedure sf_fitgraph (gp, cv)
+
+pointer	gp			# Graphics structure
+pointer	cv			# Sensitivity function curve
+
+int	i, j
+real	x1, x2, y1, y2, cveval()
+pointer	gio
+
+begin
+	gio = GP_GIO(gp)
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+
+	# Only plot on sensitivity curve graph types.
+	for (j = 1; GP_GRAPHS(gp,j) != EOS; j = j + 1) {
+	    if (GP_GRAPHS(gp,j) != 's')
+		next
+	    call gseti (gio, G_WCS, j)
+	    call ggwind (gio, x1, x2, y1, y2)
+	    x2 = (x2 - x1) / NCURVE
+	    y1 = cveval (cv, x1)
+	    call gamove (gio, x1, y1)
+	    do i = 1, NCURVE {
+	        x1 = x1 + x2
+	        y1 = cveval (cv, x1)
+	        call gadraw (gio, x1, y1)
+	    }
+	}
+end