1 files changed, 104 insertions, 0 deletions
diff --git a/noao/onedspec/sensfunc/sfcgraph.x b/noao/onedspec/sensfunc/sfcgraph.x
new file mode 100644
index 00000000..2f689d47
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfcgraph.x
@@ -0,0 +1,104 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+# SF_CGRAPH -- Graph of composite points and errors
+
+procedure sf_cgraph (gp, stds, nstds, cv)
+
+pointer	gp		# Graphics structure
+pointer	stds[nstds]	# Standard star data
+int	nstds		# Number of standard stars
+pointer	cv		# Sensitivity function curve
+
+int	i, j, n, nwaves
+real	w, s, ymin, ymax, cveval()
+double	sum, sum2
+pointer	sp, waves, sens, errors, xp, yp, zp, gio
+
+begin
+	nwaves = 0
+	do i = 1, nstds-2
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE)
+	        nwaves = nwaves + STD_NWAVES(stds[i])
+
+	call smark (sp)
+	call salloc (waves, nwaves, TY_REAL)
+	call salloc (sens, nwaves, TY_REAL)
+	call salloc (errors, nwaves, TY_REAL)
+	    
+	nwaves = 0
+	do i = 1, nstds-2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    zp = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[zp] != 0.) {
+		    Memr[waves+nwaves] = Memr[xp]
+		    Memr[sens+nwaves] = Memr[yp]
+		    nwaves = nwaves + 1
+		}
+		xp = xp + 1
+		yp = yp + 1
+		zp = zp + 1
+	    }
+	}
+	call xt_sort2 (Memr[waves], Memr[sens], nwaves)
+
+	n = 0
+	sum = 0.
+	sum2 = 0.
+	ymin = 0.
+	ymax = 0.
+	j = 0
+	do i = 1, nwaves {
+	    w = Memr[waves+i-1]
+	    s = Memr[sens+i-1]
+	    n = n + 1
+	    sum = sum + s
+	    sum2 = sum2 + s * s
+
+	    if ((i < nwaves) && (w == Memr[waves+i]))
+		next
+
+	    if (n > 2) {
+	        sum = sum / n
+	        sum2 = sum2 / n - sum * sum
+	        if (sum2 > 0)
+	            sum2 = sqrt (sum2 / n)
+	        else
+		    sum2 = 0.
+	        sum = sum - cveval (cv, w)
+
+		Memr[waves+j] = w
+		Memr[sens+j] = sum
+		Memr[errors+j] = sum2
+		j = j + 1
+
+		if (sum + sum2 > ymax)
+		    ymax = sum + sum2
+		if (sum - sum2 < ymin)
+		    ymin = sum - sum2
+	    }
+	    n = 0
+	    sum = 0.
+	    sum2 = 0.
+	}
+	nwaves = j
+
+	if (j == 0) {
+	    call printf ("No wavelength overlap for composite points")
+	} else {
+	    gio = GP_GIO(gp)
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Composite Points vs Wavelength", "", "")
+	    call gseti (gio, G_PLCOLOR, GP_CMARK(gp))
+	    do i = 0, nwaves-1
+	        call gmark (gio, Memr[waves+i], Memr[sens+i], GM_VEBAR,
+		    1., -Memr[errors+i])
+	}
+
+	call sfree (sp)
+end