Initial commit

1 files changed, 217 insertions, 0 deletions

diff --git a/noao/imred/dtoi/hdicfit/hdicebars.x b/noao/imred/dtoi/hdicfit/hdicebars.x
new file mode 100644
index 00000000..335c161d
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicebars.x
@@ -0,0 +1,217 @@
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+define	MSIZE	2.0
+
+# HDIC_EBARS -- Plot data points with their error bars as markers.  The
+# independent variable is fog subtracted transformed density.
+
+procedure hdic_ebars (ic, gp, gt, x, y, wts, ebw, npts)
+
+pointer	ic		# Pointer to ic structure
+pointer	gp		# Pointer to graphics stream
+pointer	gt		# Pointer to gtools structure
+double	x[ARB]		# Array of independent variables
+double	y[ARB]		# Array of dependent variables
+double	wts[ARB]	# Array of weights
+double	ebw[ARB]	# Error bar half width (positive number)
+int	npts		# Number of points
+
+pointer	sp, xr, yr, sz
+int	orig_mark, i, xaxis, yaxis, mark, added_mark
+real	size, dx, dy, szmk
+int	gt_geti()
+bool	fp_equald()
+include	"hdic.com"
+
+begin
+	call smark (sp)
+
+	xaxis = (IC_AXES (ic, IC_GKEY(ic), 1))
+	yaxis = (IC_AXES (ic, IC_GKEY(ic), 2))
+
+	if (xaxis == 'x' || xaxis == 'u')
+	    orig_mark = GM_HEBAR
+	else if (yaxis == 'x' || yaxis == 'u')
+	    orig_mark = GM_VEBAR
+	else {
+	    call eprintf ("Choose graph type with axes 'u' or 'x'\n")
+	    call sfree (sp)
+	    return
+	}
+
+	added_mark = GM_CIRCLE
+
+	call salloc (xr, npts, TY_REAL)
+	call salloc (yr, npts, TY_REAL)
+	call salloc (sz, npts, TY_REAL)
+
+	call achtdr (x, Memr[xr], npts)
+	call achtdr (y, Memr[yr], npts)
+	call achtdr (ebw, Memr[sz], npts)
+
+	if (IC_OVERPLOT(ic) == NO) {
+	    # Start a new plot
+	    call gclear (gp)
+
+	    # Set the graph scale and axes
+	    call gascale (gp, Memr[xr], npts, 1)
+	    call gascale (gp, Memr[yr], npts, 2)
+
+	    # If plotting HD curve, set wy2 to maxden, which may have
+	    # been updated if a new endpoint was added.
+
+	    if ((IC_AXES (ic, IC_GKEY(ic), 1) == 'y') &&
+		(IC_AXES (ic, IC_GKEY(ic), 2) == 'u')) 
+		call gswind (gp, INDEF, INDEF, INDEF, real (maxden))
+
+	    call gt_swind (gp, gt)
+	    call gt_labax (gp, gt)
+	}
+
+	call ggscale (gp, 0.0, 0.0, dx, dy)
+
+	do i = 1, npts {
+	    size = Memr[sz+i-1]		# Sizes are WCS units; transform them
+	    if (gt_geti (gt, GTTRANSPOSE) == NO) {
+	        size = size / dx
+		mark = orig_mark
+		szmk = size
+		# Check for added point
+		if (fp_equald (ebw[i], ADDED_PT)) {
+		    szmk = MSIZE
+		    mark = added_mark
+		}
+
+		# Check for deleted point
+		if (fp_equald (wts[i], 0.0D0)) {
+		    szmk = MSIZE
+		    mark = mark + GM_CROSS
+		}
+
+	        call gmark (gp, Memr[xr+i-1], Memr[yr+i-1], mark, szmk, szmk)
+
+	    } else {
+		size = size / dy
+		szmk = size
+		mark = orig_mark
+
+		# Check for added point
+		if (fp_equald (ebw[i], ADDED_PT)) {
+		    szmk = MSIZE
+		    mark = added_mark
+		}
+
+		# Check for deleted point
+		if (fp_equald (wts[i], 0.0D0)) {
+		    szmk = MSIZE
+		    mark = mark + GM_CROSS
+		}
+
+		call gmark (gp, Memr[yr+i-1], Memr[xr+i-1], mark, szmk, szmk)
+	    }
+	}
+	
+	IC_OVERPLOT(ic) = NO
+	call sfree (sp)
+end
+
+
+# HDIC_EBW -- Calculate error bar width for plotting points.  Width is
+# returned in NDC units.
+
+procedure hdic_ebw (ic, density, indv, sdev, ebw, npts)
+
+pointer	ic		# Pointer to ic structure
+double	density[ARB]	# Untransformed density NOT fog subtracted
+double	indv[ARB]	# Transformed density above fog
+double	sdev[ARB]	# Array of standard deviation values, density units
+double	ebw[ARB]	# Error bar half width (positive numbers)
+int	npts		# Number of data points
+
+double	fog
+pointer	sp, denaf, outwe
+int	xaxis, yaxis, i
+bool	fp_equald()
+real	ic_getr()
+
+begin
+	xaxis = (IC_AXES (ic, IC_GKEY(ic), 1))
+	yaxis = (IC_AXES (ic, IC_GKEY(ic), 2))
+
+	if (xaxis == 'u' || yaxis == 'u') {
+	    call amovd (sdev, ebw, npts) 
+	    return
+	}
+
+	call smark (sp)
+	call salloc (denaf, npts, TY_DOUBLE)
+	call salloc (outwe, npts, TY_DOUBLE)
+
+	fog = double (ic_getr (ic, "fog"))
+
+	call asubkd (density, fog, Memd[denaf], npts)
+	call aaddd  (Memd[denaf], sdev, Memd[denaf], npts)
+
+	call hdic_ebtran (Memd[denaf], Memd[outwe], npts, IC_TRANSFORM(ic))
+
+	# Subtract transformed values to get errors.  Then check for
+	# added points, which are flagged with ebw=ADDED_PT.
+
+	call asubd (Memd[outwe], indv, ebw, npts)
+	do i = 1, npts {
+	    if (fp_equald (sdev[i], ADDED_PT))
+		ebw[i] = ADDED_PT
+	}
+
+	call sfree (sp)
+end
+
+
+# HDIC_EBTRAN -- Apply transformation, generating a vector of independent
+# variables from a density vector.  The independent variable vector is the
+# standard deviation values and the output values are the errors.  This
+# routine checks for ADDED_PT valued standard deviation, which indicates an
+# added point.
+
+procedure hdic_ebtran (density, ind_var, npts, transform)
+
+double	density[npts]		# Density vector - input
+double	ind_var[npts]		# Ind variable vector - filled on output
+int	npts			# Length of data vectors
+int	transform		# Integer code for transform type
+
+int	i
+bool	fp_equald()
+
+begin
+	switch (transform) {
+	case HD_LOGO:
+	    do i = 1, npts {
+		if (fp_equald (density[i], 0.0))
+		    ind_var[i] = 0.0
+		else
+		    ind_var[i] = log10 ((10. ** density[i]) - 1.0)
+	    }
+	case HD_K75:
+	    do i = 1, npts {
+		if (fp_equald (density[i], 0.0))
+		    ind_var[i] = 0.0
+		else
+		    ind_var[i] = density[i] + 0.75*log10(1.- 10.**(-density[i]))
+	    }
+	case HD_K50:
+	    do i = 1, npts {
+		if (fp_equald (density[i], 0.0))
+		    ind_var[i] = 0.0
+		else
+		    ind_var[i] = density[i] + 0.50*log10(1.- 10.**(-density[i]))
+	    }
+	case HD_NONE:
+	    call amovd (density, ind_var, npts)
+	default:
+	    call error (0, "Unrecognized transformation in HDIC_EBTRAN")
+	}
+end