1 files changed, 60 insertions, 0 deletions
diff --git a/noao/imred/dtoi/hdicfit/hdicinit.x b/noao/imred/dtoi/hdicfit/hdicinit.x
new file mode 100644
index 00000000..d9730051
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicinit.x
@@ -0,0 +1,60 @@
+include	<mach.h>
+include	"hdicfit.h"
+
+# HDIC_INIT -- Initialize hdfit/icgfit interface.
+
+procedure hdic_init (density, nvalues, dmax)
+
+double	density[nvalues]	# Reference density above fog values
+int	nvalues			# Number of values in sample
+double	dmax			# Maximum possible density
+
+int	i
+pointer	sp, base
+double	xxmax, xxmin, delta_den
+bool	fp_equald()
+include	"hdic.com"
+errchk	malloc
+
+begin
+	call smark (sp)
+
+	if (den == NULL || big_den == NULL) {
+	    call malloc (den, nvalues, TY_DOUBLE)
+	    call malloc (big_den, NVALS_FIT, TY_DOUBLE)
+	} else if (nvalues != nraw) 
+	    call realloc (den, nvalues, TY_DOUBLE)
+
+	nraw = nvalues
+	call salloc (base, NVALS_FIT, TY_DOUBLE)
+
+	# Copy density array to pointer location
+	call amovd (density, Memd[den], nraw)
+
+	# Calculate big vector of density values.  The points are spaced 
+	# linear in log space, to yield adequate spacing at low density values.
+
+	call alimd (density, nraw, xxmin, xxmax)
+
+	# Put user value for maximum density in common block if it is valid.
+	if (! fp_equald (dmax, 0.0D0))
+	    maxden = dmax
+
+	# Make big_den go all the way to maxden, not just xxmax.  Make sure
+	# the value of xxmin won't cause the log function to blow up.
+
+	if (xxmin > 0.0D0)
+	    ;
+	else
+	    xxmin = 2.0 * EPSILOND
+
+	delta_den = (log10 (maxden) - log10 (xxmin)) / double (NVALS_FIT - 1)
+
+	do i = 1, NVALS_FIT
+	    Memd[big_den+i-1] = log10 (xxmin) + double (i-1) * delta_den
+
+	call amovkd (10.0D0, Memd[base], NVALS_FIT)
+	call aexpd (Memd[base], Memd[big_den], Memd[big_den], NVALS_FIT)
+
+	call sfree (sp)
+end