Repatch (from linux) of OSX IRAF

62 files changed, 7299 insertions, 0 deletions

diff --git a/noao/imred/dtoi/README b/noao/imred/dtoi/README
new file mode 100644
index 00000000..e293c6c8
--- /dev/null
+++ b/noao/imred/dtoi/README
@@ -0,0 +1 @@
+DTOI -- Density to intensity transformation package.
diff --git a/noao/imred/dtoi/Revisions b/noao/imred/dtoi/Revisions
new file mode 100644
index 00000000..e6e955ff
--- /dev/null
+++ b/noao/imred/dtoi/Revisions
@@ -0,0 +1,144 @@
+.help revisions Jun88 noao.imred.dtoi
+.nf
+
+imred$dtoi/selftest.x
+	The 'lut' was declared as TY_INT instead of TY_REAL (5/4/13)
+
+imred$dtoi/spolist.x
+	A TY_INT pointer was being used with Memr[] (4/20/13)
+
+=======
+V2.16
+=======
+
+imred$dtoi/hdtoi.x
+	The hd_fogcalc() procedure was being called with too few arguments
+	(7/12/09, MJF)
+
+=======
+V2.14.1
+=======
+
+imred$dtoi/hdicfit/mkpkg
+	Added missing dependencies.  (12/13/01, MJF)
+
+imred$dtoi/mkpkg
+	Added missing dependencies.  (10/11/99, Valdes)
+
+=======
+V2.11.2
+=======
+
+imred$dtoi/doc/hdfit.hlp
+imred$dtoi/doc/hdshift.hlp
+imred$dtoi/doc/hdtoi.hlp
+imred$dtoi/doc/spotlist.hlp
+	Fixed minor formating problems.  (4/22/99, Valdes)
+
+imred$dtoi/hdicggraph.x
+imred$dtoi/hdicebars.x
+	Replace direct access to GTOOLS structure with GTOOLS interface.
+	(12/18/98, Valdes)
+
+=======
+V2.11.1
+=======
+
+imred$dtoi/spotlist.x	7 Sept, 1989 SRo
+	In Suzanne's absence, removed the divide by nfog fog images in
+	hd_fogcalc() after a bug report by Steve Majewski that when multiple
+	fog images are used, the resultant fog values are too low by 1/nfog.
+	The total_pix accumulator is already a total for all images.  Should
+	be verified by Suzanne on her return.
+
+imred$dtoi/hdicfit/hdicgdelte.x  17 April, 1989 ShJ
+	Procedure icg_dl1 was declared as an integer procedure, although
+	it does not return a function values and was not being called
+	as a function.  
+
+imred$dtoi/hdicfit/hdicgcolon.x  27 March, 1989 ShJ
+	Changed only a comment line, to point out that the :ebars command
+	switches the meaning of the error bars between representing 
+	standard deviations or weights.  
+
+imred$dtoi/spotlist.x 11 May, 1988 ShJ
+	Added parameter "option", which selects a mean or median calculation
+	of spot densities and the fog value.  Previously, only a mean 
+	calculation was available.  The choice of "option" is written to
+	the database.  
+
+imred$dtoi/hd_aravr.x  11 May, 1988 ShJ
+	This is a new file, modified from the vops version only in that it
+	puts an upper limit on the number of iterations allowed in the
+	procedure of 10.  This should avoid the problem reported by
+	Steve Majewski of nearly saturated pixels "hanging"; they actually
+	were oscillating endlessly in the mean rejection calculation.
+
+imred$dtoi/hdtoi.x 11 May, 1988 ShJ
+	Added parameter "option", which selects a mean or median calulcation
+	of the fog value, in those cases where the user supplies a fog 
+	image, rather than a fog value.  Previously, only a mean calculation
+	was available.  This newly calculated value of fog is now also
+	written to the database.
+
+imred$dtoi/hdtoi.x 5 April, 1988 ShJ
+	Modified the way in which hdtoi scales intensities to the "ceiling"
+	parameter.  A maximum intensity is calculated, using cveval to
+	evaluate the transformed value of the maximum density above fog.
+	The fog value subtracted is hdtoi.fog, which may or may not
+	equal the value of fog written by spotlist into the database.
+	Previously, the program was incorrectly subtracting the spotlist
+	fog value from the maximum density when calculating the maximum
+	intensity and so the scale factor to represent saturated intensities
+	as "ceiling". 
+
+imred$dtoi/hdicgundel.x	25 March, 1988  ShJ
+	Procedures icg_u1d and icg_undeleted in this source file were
+	declared as functions but not returning a value.  They were not
+	being called as functions, and are now not declared as functions.
+
+imred$dtoi/database.x	14 March, 1988	ShJ
+	Added procedure to put a double precision value:  ddb_putd().
+
+imred$dtoi/spotlist.x	14 March, 1988	ShJ
+	Added code for additional input parameter: "maxad", the maximum
+	allowed input value, i.e., the AD units in a saturated pixel.  From
+	this value and the "scale" parameter, spotlist calculates the
+	maximum permissable density, and writes this to the database as a
+	double precision value. 
+
+imred$dtoi/hdfit.x	14 March, 1988	ShJ
+	Parameter "maxden" was removed, as this value is now calculated
+	precisely by spotlist and can be read from the database.
+
+imred$dtoi/spotlist.par
+imred$dtoi/hdfit.par
+imred$dtoi/doc/spotlist.hlp
+imred$dtoi/doc/hdfit.hlp	14 March, 1988	ShJ
+	Modified help pages and parameter files for spotlist and hdfit to
+	support the above changes to the user interface.
+
+imred$dtoi/hdtoi.x	11 March, 1988	ShJ
+	Bug fixed in hd_aptrans().  The equations for the k75 and k50 
+	transformations were incorrect, having a '*' where a '+' should
+	have been.  Also, the strncmp procedure was only checking the first 
+	character to distinguish between transformation types, in which
+	case k75 and k50 would not be resolved.
+
+imred$dtoi/hdfit.x	10 March, 1988	ShJ
+	Repaired an error in incrementing an offset into the output
+	arrays in hd_rdloge.   This error was seen only when more than
+	two databases contributed input values to be fit.  This was
+	a trivial arithmetic error that had been overlooked due to
+	insufficient testing.  This bug would produce a memory error
+	as values were written into unallocated locations.
+
+imred$dtoi/hdicfit/hdic.com	9 March, 1988	ShJ
+	Reordered the values in the common block, as the largest double
+	precision value was following the integers.  The SUN-4 compiler
+	gave warning of this, then padded the double value.  Also at this
+	time file hdicgundelete.x was renamed to hdicgundel.x to avoid
+	problems associated with filenames longer than 14 characters.  The 
+	hdicfit library was completely rebuilt at this time, to insure the
+	new object module would be retrieved.
+.endhelp
diff --git a/noao/imred/dtoi/database.x b/noao/imred/dtoi/database.x
new file mode 100644
index 00000000..cb501472
--- /dev/null
+++ b/noao/imred/dtoi/database.x
@@ -0,0 +1,611 @@
+include	<time.h>
+include	<ctype.h>
+include	<ctotok.h>
+include	<finfo.h>
+
+# DATABASE.X -- This file contains source for the database utilities used by
+# the DTOI package.  They are based on Frank Valdes's dtext utilities.
+
+# Definition for dbtext structure.
+
+define	DT_LEN		5
+
+define	DT		Memi[$1]	# FIO channel
+define	DT_NRECS	Memi[$1+1]	# Number of records
+define	DT_OFFSETS	Memi[$1+2]	# Pointer to record offsets
+define	DT_NAMES	Memi[$1+3]	# Pointer to name indices
+define	DT_MAP		Memi[$1+4]	# Pointer to record names
+
+define	DT_OFFSET	Meml[DT_OFFSETS($1)+$2-1]
+define	DT_NAMEI	Memi[DT_NAMES($1)+$2-1]
+define	DT_NAME		Memc[DT_MAP($1)+DT_NAMEI($1,$2)]
+
+define	DT_ALLOC	20		# Allocation block size
+
+
+# DDB_MAP -- Map the database.
+
+pointer procedure ddb_map (database, mode)
+
+char	database[ARB]			# Database file
+int	mode				# FIO mode
+
+int	i, nrec
+int	dt_alloc1, dt_alloc2
+pointer	db, str
+
+int	open(), fscan(), strlen()
+bool	streq()
+long	note()
+errchk	delete, open
+
+begin
+	call calloc (db, DT_LEN, TY_STRUCT)
+
+	if (mode == NEW_FILE)
+	    iferr (call delete (database))
+		;
+
+	DT(db) = open (database, mode, TEXT_FILE)
+
+	if (mode != READ_ONLY)
+	    return (db)
+
+	dt_alloc1 = DT_ALLOC
+	dt_alloc2 = DT_ALLOC * SZ_LINE
+	call malloc (DT_OFFSETS(db), dt_alloc1, TY_LONG)
+	call malloc (DT_NAMES(db), dt_alloc1, TY_INT)
+	call malloc (DT_MAP(db), dt_alloc2, TY_CHAR)
+	call malloc (str, SZ_LINE, TY_CHAR)
+
+	nrec = 1
+	DT_NRECS(db) = 0
+	DT_NAMEI(db, nrec) = 0
+
+	while (fscan (DT(db)) != EOF) {
+	    call gargwrd (DT_NAME(db, nrec), SZ_LINE)
+
+	    if (streq (DT_NAME(db, nrec), "begin")) {
+		call gargstr (Memc[str], SZ_LINE)
+		for (i=str; IS_WHITE(Memc[i]); i=i+1)
+		    ;
+		call strcpy (Memc[i], DT_NAME(db, nrec), SZ_LINE)
+
+		for (i = 1; i < nrec; i = i + 1)
+		    if (streq (DT_NAME(db, i), DT_NAME(db, nrec)))
+			break
+
+		if (i < nrec)
+		    DT_OFFSET(db, i) = note (DT(db))
+		else {
+		    DT_NRECS(db) = nrec
+		    DT_OFFSET(db, nrec) = note (DT(db))
+		    DT_NAMEI(db, nrec+1) = DT_NAMEI(db, nrec) +
+		        strlen (DT_NAME(db, nrec)) + 1
+		    nrec = nrec + 1
+		}
+
+		if (nrec == dt_alloc1) {
+		    dt_alloc1 = dt_alloc1 + DT_ALLOC
+		    call realloc (DT_OFFSETS(db), dt_alloc1, TY_LONG)
+		    call realloc (DT_NAMES(db), dt_alloc1, TY_INT)
+	        }
+
+	        if (DT_NAMEI(db, nrec) + SZ_LINE >= dt_alloc2) {
+		    dt_alloc2 = dt_alloc2 + DT_ALLOC * SZ_LINE
+		    call realloc (DT_MAP(db), dt_alloc2, TY_CHAR)
+		}
+	    }
+	}
+
+	call realloc (DT_MAP(db), DT_NAMEI(db, nrec), TY_CHAR)
+	call realloc (DT_OFFSETS(db), DT_NRECS(db), TY_LONG)
+	call realloc (DT_NAMES(db), DT_NRECS(db), TY_INT)
+
+	return (db)
+end
+
+
+# DDB_UNMAP -- Unmap the database.
+
+procedure ddb_unmap (db)
+
+pointer	db				# Database file descriptor
+
+begin
+	call close (DT(db))
+	call mfree (DT_MAP(db), TY_CHAR)
+	call mfree (DT_OFFSETS(db), TY_LONG)
+	call mfree (DT_NAMES(db), TY_INT)
+	call mfree (db, TY_STRUCT)
+end
+
+
+# DDB_PREC -- Write a record to the database.
+
+procedure ddb_prec (db, record)
+
+pointer	db		# Pointer to database
+char	record[ARB]	# Name of record to enter
+
+pointer	sp, entry
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "begin\t%s\n")
+	    call pargstr (record)
+
+	call fprintf (DT(db), Memc[entry])
+
+	call sfree (sp)
+end
+
+
+# DDB_PUTR -- Write a real valued field into the current record.
+
+procedure ddb_putr (db, field, value)
+
+pointer	db		# Pointer to database
+char	field[ARB]	# Format string
+real	value		# Real value to output
+
+pointer	sp, entry
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%g\n")
+	    call pargstr (field)
+	    call pargr (value)
+	
+	call fprintf (DT(db), Memc[entry])
+
+	call sfree (sp)
+end
+
+# DDB_PUTD -- Write a double valued field into the current record.
+
+procedure ddb_putd (db, field, value)
+
+pointer	db		# Pointer to database
+char	field[ARB]	# Format string
+double	value		# Real value to output
+
+pointer	sp, entry
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%g\n")
+	    call pargstr (field)
+	    call pargd (value)
+	
+	call fprintf (DT(db), Memc[entry])
+
+	call sfree (sp)
+end
+
+
+# DDB_PUTI -- Write an integer valued field into the current record.
+
+procedure ddb_puti (db, field, value)
+
+pointer	db		# Pointer to database
+char	field[ARB]	# Format string
+int	value		# Integer value to output
+
+pointer	sp, entry
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%d\n")
+	    call pargstr (field)
+	    call pargi (value)
+	
+	call fprintf (DT(db), Memc[entry])
+
+	call sfree (sp)
+end
+
+
+# DDB_PSTR -- Write a string valued field into the current record.
+
+procedure ddb_pstr (db, field, value)
+
+pointer	db		# Pointer to database
+char	field[ARB]	# Format string
+char	value[ARB]	# String field
+
+pointer	sp, entry
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%s\n")
+	    call pargstr (field)
+	    call pargstr (value)
+	
+	call fprintf (DT(db), Memc[entry])
+
+	call sfree (sp)
+end
+
+
+# DDB_PAR -- Put an array of real values into a field in the current record.
+
+procedure ddb_par (db, field, array, nvals)
+
+pointer	db		# Pointer to database structure
+char	field[ARB]	# Name of field to be added
+real	array[nvals]	# Array of real values 
+int	nvals		# Number of values in array
+
+pointer	sp, entry
+int	i
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%d\n")
+	    call pargstr (field)
+	    call pargi (nvals)
+
+	call fprintf (DT(db), Memc[entry])
+
+	do i = 1, nvals {
+	    call sprintf (Memc[entry], SZ_LINE, "\t\t%g\n")
+		call pargr (array[i])
+	    
+	    call fprintf (DT(db), Memc[entry])
+	}
+
+	call sfree (sp)
+end
+
+
+# DDB_PAD -- Put an array of double values into a field in the current record.
+
+procedure ddb_pad (db, field, array, nvals)
+
+pointer	db		# Pointer to database structure
+char	field[ARB]	# Name of field to be added
+double	array[nvals]	# Array of double values 
+int	nvals		# Number of values in array
+
+pointer	sp, entry
+int	i
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%d\n")
+	    call pargstr (field)
+	    call pargi (nvals)
+
+	call fprintf (DT(db), Memc[entry])
+
+	do i = 1, nvals {
+	    call sprintf (Memc[entry], SZ_LINE, "\t\t%g\n")
+		call pargd (array[i])
+	    
+	    call fprintf (DT(db), Memc[entry])
+	}
+
+	call sfree (sp)
+end
+
+
+# DDB_PAI -- Put an array of integer values into a field in the current 
+# record.
+
+procedure ddb_pai (db, field, array, nvals)
+
+pointer	db		# Pointer to database structure
+char	field[ARB]	# Name of field to be added
+int	array[nvals]	# Array of integer values 
+int	nvals		# Number of values in array
+
+pointer	sp, entry
+int	i
+
+begin
+	call smark (sp)
+	call salloc (entry, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[entry], SZ_LINE, "\t%s\t%d\n")
+	    call pargstr (field)
+	    call pargi (nvals)
+
+	call fprintf (DT(db), Memc[entry])
+
+	do i = 1, nvals {
+	    call sprintf (Memc[entry], SZ_LINE, "\t\t%d\n")
+		call pargi (array[i])
+	    
+	    call fprintf (DT(db), Memc[entry])
+	}
+
+	call sfree (sp)
+end
+
+
+# DDB_LOCATE -- Locate a database record.
+
+int procedure ddb_locate (db, name)
+
+pointer	db				# DTTEXT pointer
+char	name[ARB]			# Record name
+
+int	i
+
+bool	streq()
+pointer	sp, err_string
+
+begin
+	do i = 1, DT_NRECS(db) {
+	    if (streq (name, DT_NAME(db, i)))
+		return (i)
+	}
+
+	# The record was not found
+
+	call smark (sp)
+	call salloc (err_string, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[err_string], SZ_LINE, 
+	    "DDB_LOCATE: Database record '%s' not found")
+	        call pargstr (name)
+
+	call error (21, Memc[err_string])
+	call sfree (sp)
+end
+
+
+# DDB_GSTR -- Get a string field
+
+procedure ddb_gstr (db, record, field, str, maxchar)
+
+pointer	db				# Database file descriptor
+int	record				# Database index
+char	field[ARB]			# Database field
+char	str[maxchar]			# String value
+int	maxchar				# Maximum characters for string
+
+char	name[SZ_LINE]
+
+pointer	sp, esp
+int	i, fscan()
+bool	streq()
+
+begin
+	if ((record < 1) || (record > DT_NRECS(db)))
+	    call error (22, "Database record request out of bounds")
+
+	call seek (DT(db), DT_OFFSET(db, record))
+
+	while (fscan (DT(db)) != EOF) {
+	    call gargwrd (name, SZ_LINE)
+
+	    if (streq (name, "begin"))
+		break
+	    else if (streq (name, field)) {
+		call gargstr (str, maxchar)
+		for (i=1; IS_WHITE(str[i]); i=i+1)
+		    ;
+		if (i>1)
+		    call strcpy (str[i], str, maxchar)
+		return
+	    }
+	}
+
+	call smark (sp)
+	call salloc (esp, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[esp], SZ_LINE, 
+	    "DDB_GSTR: Database field '%s' not found")
+		call pargstr (field)
+
+	call error (23, Memc[esp])
+	call sfree (sp)
+end
+
+
+# DDB_GETI -- Get an integer field.
+
+int procedure ddb_geti (db, record, field)
+
+pointer	db				# DTTEXT pointer
+int	record				# Database index
+char	field[ARB]			# Database field
+
+real	rval
+bool    fp_equalr()
+real	ddb_getr()
+
+errchk	ddb_getr
+
+begin
+	rval = ddb_getr (db, record, field)
+
+	if (fp_equalr (rval, INDEFR))
+	    return (INDEFI)
+	else 
+	    return (int (rval))
+end
+
+
+# DDB_GETR -- Get an real field
+
+real procedure ddb_getr (db, record, field)
+
+pointer	db				# DTTEXT pointer
+int	record				# Database index
+char	field[ARB]			# Database field
+
+pointer	sp, esp
+real	rval
+char	name[SZ_LINE]
+
+int	fscan(), nscan()
+bool	streq()
+
+begin
+	if ((record < 1) || (record > DT_NRECS(db)))
+	    call error (24, "Database record request out of bounds")
+
+	call seek (DT(db), DT_OFFSET(db, record))
+
+	while (fscan (DT(db)) != EOF) {
+	    call gargwrd (name, SZ_LINE)
+
+	    if (streq (name, "begin"))
+		break
+	    else if (streq (name, field)) {
+		call gargr (rval)
+		if (nscan() == 2)
+		   return (rval)
+		else
+		   call error (25, "Error in database field value")
+	    }
+	}
+
+	call smark (sp)
+	call salloc (esp, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[esp], SZ_LINE, 
+	    "DDB_GET: Database field '%s' not found")
+		call pargstr (field)
+
+	call error (26, Memc[esp])
+	call sfree (sp)
+end
+
+
+# DDB_GAR -- Get a real array field
+
+procedure ddb_gar (db, record, field, array, len_array, npts)
+
+pointer	db				# DTTEXT pointer
+int	record				# Database index
+char	field[ARB]			# Database field
+real	array[len_array]		# Array values
+int	len_array			# Length of array
+int	npts				# Number of values in the field
+
+int	i
+double	tmp
+pointer	sp, dubble
+bool	fp_equald()
+
+begin
+	call smark (sp)
+	call salloc (dubble, npts, TY_DOUBLE)
+
+	call ddb_gad (db, record, field, Memd[dubble], len_array, npts)
+	do i = 1, npts {
+	    tmp = Memd[dubble+i-1]
+	    if (fp_equald (tmp, INDEFD))
+		array[i] = INDEFR
+	    else
+		array[i] = real (tmp)
+	}
+
+	call sfree (sp)
+end
+
+
+# DDB_GAD -- Get a double array field
+
+procedure ddb_gad (db, record, field, array, len_array, npts)
+
+pointer	db				# DTTEXT pointer
+int	record				# Database index
+char	field[ARB]			# Database field
+double	array[len_array]		# Array values
+int	len_array			# Length of array
+int	npts				# Number of points in the array
+
+pointer	sp, esp
+char	name[SZ_LINE]
+int	i
+
+int	fscan(), nscan()
+bool	streq()
+
+begin
+	if ((record < 1) || (record > DT_NRECS(db)))
+	    call error (27, "Database record request out of bounds")
+
+	call seek (DT(db), DT_OFFSET(db, record))
+
+	while (fscan (DT(db)) != EOF) {
+	    call gargwrd (name, SZ_LINE)
+
+	    if (streq (name, "begin"))
+		break
+	    else if (streq (name, field)) {
+		call gargi (npts)
+		if (nscan() != 2)
+		    call error (28, "Error in database field value")
+
+		npts = min (npts, len_array)
+		for (i = 1; i <= npts; i = i + 1) {
+		    if (fscan (DT(db)) == EOF)
+		        call error (29, "Error in database field value")
+
+		    call gargd (array[i])
+		    if (nscan() != 1)
+		        call error (30, "Error in database field value")
+		}
+		return
+	    }
+	}
+
+	call smark (sp)
+	call salloc (esp, SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[esp], SZ_LINE, 
+	    "DDB_GA: Database field '%s' not found")
+	        call pargstr (field)
+
+	call error (31, Memc[esp])
+	call sfree (sp)
+end
+
+
+# DDB_PTIME -- Put a time string with a comment
+
+procedure ddb_ptime (db)
+
+pointer	db				# DTTEXT pointer
+
+char	timestr[SZ_TIME]
+long	time, clktime()
+
+begin
+	time = clktime (0)
+	call cnvtime (time, timestr, SZ_TIME)
+	call fprintf (DT(db), "# %s\n")
+	    call pargstr (timestr)
+end
+
+
+# DDB_SCAN -- Scan a line from the database.
+
+int procedure ddb_scan (db)
+
+pointer	db			# DDB pointer
+int	fscan()
+
+begin
+	return (fscan (DT(db)))
+end
diff --git a/noao/imred/dtoi/dematch.par b/noao/imred/dtoi/dematch.par
new file mode 100644
index 00000000..cd7a5cfc
--- /dev/null
+++ b/noao/imred/dtoi/dematch.par
@@ -0,0 +1,8 @@
+# Cl parameters are:
+database,f,a,,,,Database of density information
+wedge,f,a,"",,,Wedge number used for calibration
+filter,f,a,"",,,Filter used for calibration
+emulsion,f,a,"",,,Emulsion used for calibration
+wedgefile,f,h,"noao$lib/hdwedge.dat",,,File containing exposure information
+nskip,i,h,0,,,Number of faint spots skipped
+verbose,b,h,yes,,,Print exposure record from wedgefile
diff --git a/noao/imred/dtoi/dematch.x b/noao/imred/dtoi/dematch.x
new file mode 100644
index 00000000..ff47e90f
--- /dev/null
+++ b/noao/imred/dtoi/dematch.x
@@ -0,0 +1,160 @@
+include	<error.h>
+
+# T_DEMATCH -- matches densities listed in the input database to
+# log exposure values retrieved from a system maintained or user
+# provided file.  The output matches are added to the input database.
+
+procedure t_dematch ()
+
+pointer	filter, emulsion, wedge_db, density_db, db, exp, den, sp
+pointer wedge, expo
+int	nskip, rec, nvalues
+
+pointer	ddb_map()
+bool	clgetb()
+int	ddb_locate(), ddb_geti(), clgeti()
+
+begin
+	call smark (sp)
+	call salloc (filter, SZ_FNAME, TY_CHAR)
+	call salloc (emulsion,  SZ_FNAME, TY_CHAR)
+	call salloc (wedge_db, SZ_FNAME, TY_CHAR)
+	call salloc (density_db, SZ_FNAME, TY_CHAR)
+	call salloc (wedge, SZ_FNAME, TY_CHAR)
+
+	# Get parameters
+	call clgstr ("database", Memc[density_db], SZ_FNAME)
+	call clgstr ("wedge", Memc[wedge], SZ_FNAME)
+	call clgstr ("filter", Memc[filter], SZ_FNAME)
+	call clgstr ("emulsion", Memc[emulsion], SZ_FNAME)
+	call clgstr ("wedgefile", Memc[wedge_db], SZ_FNAME)
+	nskip = clgeti ("nskip")
+
+	# Retrieve exposure information; one wedge per run
+	call hd_rwedge (Memc[wedge_db], exp, Memc[wedge], Memc[filter], 
+	    Memc[emulsion], clgetb("verbose"))
+
+	db = ddb_map (Memc[density_db], READ_ONLY)
+	iferr {
+	    rec = ddb_locate (db, "density")
+	    nvalues = ddb_geti (db, rec, "den_val")
+	} then
+	    call error (0, "Error locating density record in database")
+
+	call salloc (den, nvalues, TY_REAL)
+	iferr (call ddb_gar (db, rec, "den_val", Memr[den], nvalues, nvalues))
+	    call error (0, "Error reading density information")
+
+	# Close db file before reopening for append.
+	call ddb_unmap (db)
+
+	# Add fields for wedge, filter and plate as "calibrate" record
+	db = ddb_map (Memc[density_db], APPEND)
+	call ddb_prec (db, "calibrate")
+	call ddb_pstr (db, "wedge", Memc[wedge])
+	call ddb_pstr (db, "filter", Memc[filter])
+	call ddb_pstr (db, "emulsion", Memc[emulsion])
+
+	# Exposures are returned in increasing order.  Make sure the
+	# exposures are output in the same order as density values.
+
+	call salloc (expo, nvalues, TY_REAL)
+	call amovr (Memr[exp+nskip], Memr[expo], nvalues)
+
+	if (Memr[den] > Memr[den+nvalues-1])
+	    call hd_reorderr (Memr[expo], nvalues)
+
+	# Now add exposure values to database
+	call ddb_ptime (db)
+	call ddb_prec (db, "exposure")
+	call ddb_par (db, "log_exp", Memr[expo], nvalues)
+
+	call ddb_unmap (db)
+
+	call mfree (exp, TY_REAL)
+	call sfree (sp)
+end
+
+
+# HD_RWEDGE -- Read wedge information from database file for a given
+# wedge, filter, emulsion combination.  A pointer to the extracted
+# exposure values is returned as an argument.
+
+procedure hd_rwedge (db_file, exp, wedge, filter, emulsion, verbose)
+
+char	db_file[SZ_FNAME]	# Name of database with exposure information
+pointer	exp			# Pointer to array of exposure values - output
+char	wedge[SZ_FNAME]		# Wedge number
+char	filter[SZ_FNAME]	# Filter used
+char	emulsion[SZ_FNAME]	# Emulsion used
+bool	verbose			# Print record of exposure information?
+
+pointer	db
+char	wfe[SZ_FNAME]
+int	rec, nvalues, stat, i
+
+pointer ddb_map()
+int	ddb_locate(), ddb_geti(), ddb_scan()
+errchk	ddb_map, ddb_scan
+
+begin
+	# Convert strings to upper case for matching in database
+	call strupr (wedge)
+	call strupr (filter)
+	call strupr (emulsion)
+
+	db = ddb_map (db_file, READ_ONLY)
+	iferr (rec = ddb_locate (db, wedge))
+	    call erract (EA_FATAL)
+
+	# Construct field name of filter/emulsion combination in question
+	call sprintf (wfe, SZ_FNAME, "%s/%s")
+	    call pargstr (emulsion)
+	    call pargstr (filter)
+
+	# Retrieve exposure values from database file
+	iferr (nvalues = ddb_geti (db, rec, wfe))
+	    call erract (EA_FATAL)
+
+	call malloc (exp, nvalues, TY_REAL)
+	stat = ddb_scan (db)
+
+	do i = 1, nvalues {
+	    call gargr (Memr[exp+i-1])
+
+	    # Calibration values are stored 8 values per line
+	    if (mod (i, 8) == 0) {
+		if (nvalues > i)
+		    stat = ddb_scan (db)
+	    }
+	}
+	 
+	if (verbose) {
+	    call printf ("\nCalibration log exposure values for %s/%s: \n")
+	        call pargstr (wedge)
+	        call pargstr (wfe)
+	    
+	    do i = 1, nvalues {
+	        call printf ("%8g ")
+		    call pargr (Memr[exp+i-1])
+	        if (mod (i, 8) == 0)
+		    call printf ("\n")
+	    }
+
+	    # Need a newline if the last line didn't have 8 entries.
+	    if (mod (nvalues, 8) != 0)
+		call printf ("\n")
+
+	    call flush (STDOUT)
+	}
+
+	# Exposures are returned sorted in increasing order.  That is,
+	# the exposure for the lightest spot is element one, and the
+	# darkest spot's exposure value is the last array element.
+
+	if (Memr[exp] > Memr[exp+nvalues-1])
+	    # Out of order - flip elements
+	    call hd_reorderr (Memr[exp], nvalues)
+
+	call ddb_unmap (db)
+end
diff --git a/noao/imred/dtoi/doc/dematch.hlp b/noao/imred/dtoi/doc/dematch.hlp
new file mode 100644
index 00000000..7dffad26
--- /dev/null
+++ b/noao/imred/dtoi/doc/dematch.hlp
@@ -0,0 +1,51 @@
+.help dematch Feb87 imred.dtoi
+.ih
+NAME
+dematch -- match density to log exposure values
+.ih
+USAGE
+dematch database 
+.ih
+PARAMETERS
+.ls database
+Database containing density list, probably from \fIspotlist\fR.
+.le
+.ls wedge = "", filter = "", emulsion = ""
+Information used to retrieve log exposure values from \fBwedgefile\fR.
+.le
+.ls wedgefile = "noao$lib/hdwedge.dat"
+Name of file containing wedge intensity information.
+.le
+.ls nskip = 0
+Number of faint spots skipped, used as an offset into the list of
+log exposure values.
+.le
+.ls verbose = yes
+Print the log exposure information to STDOUT as well as to \fBdatabase\fR.
+.le
+.ih
+DESCRIPTION
+Task \fIdematch\fR matches density values to log exposure values.  A database
+of density values is input, as well as information needed to 
+retrieve log exposure values from a reference file.  The two sources of 
+information are matched, and the matching log exposure values are added 
+as a record in the database.
+
+Parameter \fBnskip\fR tells how many faint spots were not
+included in the density \fBdatabase\fR.  This information is
+used to align the density, exposure values.  It doesn't matter if the 
+densities are listed in a monotonically increasing or decreasing
+order, as long as no spots were omitted between the first and last
+measured.
+.ih
+EXAMPLES
+Match densities in db1 to log exposure values for wedge#117
+with a IIIAJ emulsion and a GG385 filter.
+.nf
+
+	cl> dematch db1 wedge=117 filt=gg385 emulsion=IIIAJ
+.fi
+.ih
+SEE ALSO
+spotlist, hdfit, hdtoi
+.endhelp
diff --git a/noao/imred/dtoi/doc/dtoi.ms b/noao/imred/dtoi/doc/dtoi.ms
new file mode 100644
index 00000000..4f999259
--- /dev/null
+++ b/noao/imred/dtoi/doc/dtoi.ms
@@ -0,0 +1,576 @@
+.RP
+.TL
+An Overview of the IRAF DTOI Package
+.AU
+Suzanne Hammond Jacoby
+.AI
+IRAF Group - Central Computer Services
+.K2 "" "" "*"
+February 1987
+.br
+Revised July 1988
+
+.AB
+This document describes the DTOI package, which contains tasks
+for determining and applying a density to intensity transformation to
+photographic data.  The transformation is determined from a set
+of calibration spots with known relative intensities.  A curve is
+interactively fit to the densities and intensities of the calibration
+spots.  The transformation is then applied and a new output image written.
+.AE
+
+.NH
+Introduction
+.PP
+The DTOI package contains tasks for computing and applying a density
+to intensity transformation to photographic data.  These tasks perform the
+standard steps in linearizing data: calculating HD data points from 
+calibration spots, fitting a curve to these points and applying the HD
+curve to the data.  It is also possible to combine related HD curves.
+Communication between the tasks is via text files which the user can 
+inspect or modify.  It is intended
+to be easy for users to introduce data from outside the DTOI package
+into the processing.
+.PP
+There are currently six tasks in the package.  They are:
+
+.ce
+The \fBDTOI\fR Package
+.TS
+center;
+n.
+spotlist \&- Calculate densities and weights of calibration spots.
+dematch \&- Match densities to log exposure values.
+hdfit \&- Fit characteristic curve to density, exposure data.
+hdtoi \&- Apply HD transformation to image data.
+hdshift \&- Align related characteristic curves.
+selftest \&- Test transformation algorithm.
+.TE
+.PP
+The DTOI package does not currently support the self calibration of images, 
+but the addition of this capability is planned.  This would involve
+determining the HD curve from the data itself, by assuming the point spread
+function scales linearly with intensity.
+.PP
+Upon entering the package, your calibration spots and images to be
+transformed should be on the disk in IRAF image format.
+
+.NH 
+Determining the HD Curve Data 
+.PP
+To determine the HD curve, you need two sets of data:  the
+measured photographic densities of a set of calibration spots and
+the log exposure values corresponding to these measurements.  The
+log exposure values must be known a priori.  Tasks \fIspotlist\fR and
+\fIdematch\fR are used to assemble these two data sets.
+.NH 2
+SPOTLIST
+.PP
+The first step is to calculate the density of
+the calibration spots, each of which is a separate IRAF image or image
+section.  The spot density is either the median of the spot pixels or
+the mean of the pixels when pixels more then a user specified number of
+standard deviations away from the mean have been rejected.  The numbers
+in the spot image must be scaled to density; parameter \fBspotlist.scale\fR
+is used such that density = input_value * scale.  Task \fIspotlist\fR also
+calculates the standard deviation of each spot and reports
+the number of good pixels, i.e., the number of pixels not rejected
+when determining the mean density.  
+The final product of this task is a record in the data base containing a
+density for each spot.  The scale factor used is also written to the data
+base; it will be read later in task \fIhdtoi\fR.
+.NH 2
+DEMATCH
+.PP
+Log exposure values must be matched to the measured density values.  These
+log exposure values must be known a priori and will be read from a file.
+Task \fIdematch\fR retrieves the proper exposure information by
+matching the wedge number, emulsion type and filter used.  Once a match
+has been made, the proper log exposure values are written to a record
+in the database.  
+.PP
+A database of log exposure values for the NOAO standard wedges is maintained 
+in a system file; the wedge/emulsion/filter combinations available are listed
+in last section of this document.  This file can be replaced with one specific 
+to any institution; the file name is supplied to task \fIdematch\fR as a 
+parameter.  In this way the wedge file can be personalized to any application 
+and not be lost when the system is updated.
+
+.NH
+Fitting the Curve
+.PP
+The HD curve, or characteristic curve, is a plot of density versus log 
+exposure.  This curve is determined from the data points generated by
+tasks \fIspotlist\fR and \fIdematch\fR.  The objective is to fit
+a curve to these points, such that Log exposure = F(Density).  The
+technique available in this package allows the independent variable of the
+fit to be a transformation of the density (log opacitance, for example).
+The log exposure and density values are
+read from the database.  If multiple entries for a particular record are 
+present in the database, the last one is used.
+.NH 2
+HDFIT
+.PP
+Task \fIhdfit\fR fits a characteristic curve to density and log exposure
+values in preparation for transforming an image from density to intensity.
+Five functional forms of the curve are available:
+.nf
+
+ 	  Power Series
+    	  Linear Spline
+    	  Cubic Spline
+    	  Legendre Polynomials
+    	  Chebyshev Polynomials
+ 	
+.fi
+.LP 
+It is possible to apply a transformation to the 
+independent variable (density above fog) prior to the fit.  The traditional 
+choice is to fit log exposure
+as a function of the log opacitance, rather than density directly.  This is
+sometimes referred to as the Baker, or Seidel, function.  Transforming
+the density has the effect of stretching the low density data points, which
+tend to be relatively oversampled. 
+In the DTOI package, four independendent variables are currently available:
+.nf
+
+ 	Density
+ 	Log Opacitance
+ 	K50 - (Kaiser* Transform with alpha = 0.50)
+ 	K75 - (Kaiser* Transform with alpha = 0.75)
+
+.fi
+.FS
+* Margoshes and Rasberry, Spectrochimica Acta, Vol 24B, p497, (1969)
+.FE
+Any combination of transformation type and fitting function can be used and
+changed interactively.  Two combinations of interest are discussed here.
+
+The default fit is a power series fit where the independent variable is 
+Log Opacitance.  That is:
+.LP
+.EQ
+
+	"Log Exposure = " sum from k=0 to {ncoeff - 1} {A sub k Y sup k}
+
+.EN
+.sp 1
+.EQ
+	"where Y = Log Opacitance = "Log sub 10 (10 sup Density - 1)
+.EN
+.LP
+A fit that is expected to best model a IIIA-J emulsion is a power series
+fit to a K75 transform of the density.  That is,
+.LP
+.EQ
+
+	"Log Exposure = "sum from k=0 to {ncoeff - 1} {A sub k Y sup k}
+
+.EN
+.sp 1
+.EQ
+"where Y = K75 transform = Density + 0.75 " Log sub 10 (1 - 10 sup -Density )
+.EN
+.LP
+Over the expected small dynamic range in variables of the fit, legendre
+and chebyshev functions offer no advantages over a simple power series
+functional form.  The cubic and linear spline fits may follow the data very
+closely, but with typically sparse data sets this is not desirable.  It
+is expected that power series fit will serve satisfactorily in all cases.
+
+.NH 3
+Interactive Curve Fitting
+.PP
+Task \fIhdfit\fR can be run interactively or not.  In interactive mode,
+points in the sample can be edited, added or deleted.  Weighting values
+can be changed as well as the fog value, the type of transformation
+and the fitting function chosen.  To obtain the best fit possible, interactive
+fitting is recommended.  A complete list of the available commands
+is printed here; this list is also available interactively with the 
+keystroke '\fL?\fR'.
+.TS
+center;
+c s s w(3.0i)
+c l s.
+
+	DTOI INTERACTIVE CURVE FITTING OPTIONS 
+
+\fL?\fR	Print options
+\fLa\fR	Add the point at the cursor position to the sample
+\fLc\fR	Print the coordinates and fit of point nearest the cursor
+\fLd\fR	Delete data point nearest the cursor
+\fLf\fR	Fit the data and redraw or overplot
+\fLg\fR	T{
+Redefine graph keys.  Any of the following data types may be along
+either axis:
+T}
+.T&
+l l l.
+ 	\fLx\fR  Independent variable	\fLy\fR  Dependent variable
+ 	\fLf\fR  Fitted value	\fLr\fR  Residual (y - f)
+ 	\fLd\fR  Ratio (y / f)	\fLn\fR  Nonlinear part of y
+ 	\fLu\fR  Density above fog
+
+Graph keys:	 	 
+.T&
+c l s.
+
+\fLh\fR	h = (x,y)  transformed density vs. log exposure
+\fLi\fR	i = (y,x)  log exposure vs. transformed density
+\fLj\fR	j = (x,r)  transformed density vs. residuals
+\fLk\fR	k = (x,d)  transformed density vs. the y(data)/y(fit) ratio
+\fLl\fR	l = (y,u)  log exposure vs. density above fog (HD Curve)
+
+\fLo\fR	Overplot the next graph
+\fLq\fR	T{
+Terminate the interactive curve fitting, updating the database file. 
+T}
+\fLr\fR	Redraw graph
+\fLu\fR	Undelete the deleted point nearest the cursor
+\fLw\fR	Set the graph window.  For help type 'w' followed by '?'.
+\fLx\fR	Change the x value of the point nearest the cursor
+\fLy\fR	Change the y value of the point nearest the cursor
+\fLz\fR	Change the weight of the point nearest the cursor
+
+.T&
+l s s w(3.0i).
+T{
+The parameters are listed or set with the following commands which may be
+abbreviated.  To list the value of a parameter type the command alone.
+T}
+
+.T&
+l l s.
+
+\fL:show \fR[\fIfile\fR]	Show the values of all the parameters
+\fL:vshow \fR[\fIfile\fR]	Show the values of all the parameters verbosely
+\fL:errors \fR[\fIfile\fR]	Print the errors of the fit (default STDOUT)
+\fL:reset \fR	T{
+Return to original conditions of x, y, wts and npts.
+T}
+\fL:ebars \fR[\fIerrors/weights\fR]	T{
+The size of marker type '[hv]ebars' can show either standard deviations or
+relative weights.
+T}
+\fL:function \fR[\fIvalue\fR]	T{
+Fitting function (power, chebyshev, legendre, spline3, or spline1)
+T}
+\fL:transform \fR[\fIvalue\fR]	Set the transform type (none, logo, k50, k75)
+\fL:fog \fR[\fIvalue\fR]	Change the fog level (or ":fog reset")
+\fL:order \fR[\fIvalue\fR]	Fitting function order
+\fL:quit \fR	Terminate HDFIT without updating database
+\fL:/mark \fRstring	T{
+Mark type (point, box, plus, cross, diamond, hline, vline, hebar, vebar, circle)
+T}
+
+.T&
+l s s.
+T{
+Additional commands are available for setting graph formats and manipulating
+the graphics.  Use the following commands for help.
+T}
+
+.T&
+l l s.
+\fL:/help\fR	Print help for graph formatting option
+\fL:.help\fR	Print cursor mode help
+
+.TE
+.PP
+The value of fog can be changed interactively if you have
+reason to override the value written in the database by \fIspotlist\fR.
+You can reset the fog to its original value with the command ":fog reset".
+A common problem with defining the HD curve is that some of 
+the calibration spot densities fall below fog.  This is caused by either
+the low signal to noise at low densities or by making a poor choice of 
+where to scan the fog level.  These points are rejected from the fit
+when a transformation of the density is being made, as the transform cannot
+be evaluated for negative density.  If the fog value or transformation 
+type is interactively changed so this problem no longer exists, 
+the spot densities are restored in the sample.
+
+The parameters of the final fit are written to a database which then
+contains the information
+necessary to reinitialize the curfit package for applying the transformation
+in \fIhdtoi\fR.  
+
+.NH
+Applying the Transform
+.PP
+.NH 2
+HDTOI
+.PP
+Once the HD curve has been defined, it is applied to a density image
+in task \fIhdtoi\fR.
+Here the transformation is applied, as described by the fit parameters
+stored in the database.  If more than one record of fit parameters is
+present, the last one is used.  This means task \fIhdfit\fR can be
+repeated until an acceptable solution is found; the last solution will
+be used by \fIhdtoi\fR.  On output, a new output image is written; the 
+input image is left intact.
+.PP
+The transformation is accomplished by using a look-up table.  All possible
+input values, from the minimum to maximum values found in the image, are
+converted to density using the scale value read from the database, and then 
+to intensity using the fit parameters determined by \fIhdfit\fR.  The input
+value is then the index into the intensity table: 
+intensity = look_up_table (input_value).
+.PP
+A scaling factor can be applied to the final intensities, as typically
+they will be < 1.0.  (The maximum log exposure in the NOAO wedge database
+is 0.0)  By default, a saturated density pixel will be assigned the "ceiling"
+intensity of 30000 and the other pixels are scaled accordingly.  
+The user is responsible for choosing a ceiling value
+that will avoid having significant digits truncated.
+The precision 
+of the transform is unaffected by scaling the
+final intensities, although caution must be used if the output image
+pixel type is an integer.  
+.PP
+The value of fog to be used is entered by the user, and can be either
+a number or a list of file names from which to calculate the fog value.
+The fog value is subtracted from the input image before the transformation
+takes place.
+Again, consider density values below fog.  Two choices are available for
+these densities: the calculated intensity can be equal to the constant 
+value 0.0 or equal -1.0 times the intensity determined for absolute (density).
+
+.NH
+Aligning Related HD curves
+.PP
+Calibration data sets from several plates can be combined once a shift
+particular to each set has been removed.  "Different spot exposures 
+define a series of HD curves which are parallel but mutually separated
+by arbitrary shifts in log exposure, produced by differing lamp intensities
+or exposure times.  Generally, Kodak spectroscopic plates can be
+averaged if [1] they come from the same emulsion batch and box, [2]
+they receive identical hypersensitizing, [3] they are exposed similarly and
+[4] they receive the same development." * 
+.FS
+* "Averaging Photographic Characteristic Curves", John Kormendy, from
+"ESO Workshop on Two Dimensional Photometry", Edited by P. Crane and
+K.Kjar, p 69, (1980), an ESO Publication.
+.FE
+.NH 2
+HDSHIFT
+.PP
+Procedure \fIhdshift\fR calculates and subtracts a zero point shift to 
+bring several related HD curves into alignment.  The individual shifts 
+are calculated by elimination of the first coefficient (Bevington, eqn 9-3):
+.EQ
+
+a0 = y bar - a sub 1 X bar - a sub 2 X bar sup 2 - ~ ...~ - a sub n X bar sup n
+
+.EN
+Here, the averages over y and X refer to individual calibration set averages; 
+the coefficients a1, ... an were previously calculated using data from all 
+calibration sets with task \fIhdfit\fR, and stored in the database.  The
+a0 term is calculated individually for each database; this term represents
+the zero point shift in log exposure and will be different for each database.
+
+On output, the log exposure values in each database have been 
+shifted to the zero point shift of the first database in the list.  The
+log exposure records are now aligned and it would be appropriate
+to run \fIhdfit\fR on the modified database list.
+.NH
+Testing the Transformation Algorithm
+.PP
+A test task is included to see if any numerical errors were introduced
+during the density to intensity transformation.  It also evaluates
+truncation errors produced when an output image with integer pixels,
+rather than reals, is written.
+.NH 2
+SELFTEST
+.PP
+An intensity vector is generated from a density vector in two different
+ways.  The first method uses the density vector and known coefficients
+to compute the intensity.  The second method uses the curfit package
+to generate a look up table of intensities as done in task \fIhdtoi\fR.  The
+residual of the two vectors is plotted; ideally the difference between
+the 'known' and 'calculated' intensity is zero.
+.PP
+Task \fIselftest\fR also plots intensity as a function of density for
+both integer and real output pixels.  The user should investigate the
+plot with the cursor zoom and expand capabilities to determine if
+truncation errors are significant.
+.NH
+The Wedgefile Database
+.PP
+Task \fIdematch\fR reads a database and retrieves log exposure information 
+for certain combinations of wedge number, photographic emulsion and filter.  
+Those combinations included in the NOAO database are listed in the next 
+section, although any calibration data can be included if the values are
+known.  To modify the database, it is recommended that
+you generate a new file rather than add records to the existing file.  This
+way, the modifications will not be lost when a new version of the IRAF
+system is released.  
+
+In the database, the information for each wedge makes up a separate record;
+each record starts with the word \fBbegin\fR.  Each record has a title field 
+and can have multiple emulsion/filter fields.  The number of log exposure
+values must be given, followed by the values written 8 per line.  The order 
+of the exposure data can be either monotonically increasing or decreasing.  
+Here is an example:
+.DS
+begin	115
+	title	MAYALL 4-M PF BEFORE 15APR74 (CHROME) [MP1-MP968]
+	IIIAJ/UG2	16
+		  0.000 -0.160 -0.419 -0.671 -0.872 -1.153 -1.471 -1.765
+		 -2.106 -2.342 -2.614 -2.876 -3.183 -3.555 -3.911 -4.058
+	IIAO/UG2	16
+		  0.000 -0.160 -0.418 -0.670 -0.871 -1.152 -1.468 -1.761
+		 -2.102 -2.338 -2.609 -2.870 -3.176 -3.547 -3.901 -4.047
+
+.DE
+.NH 2
+Contents of the NOAO Wedgefile
+.LP
+The following table lists the wedge/emulsion/filter combinations available in
+the NOAO wedgefile database.  
+.TS
+center;
+l l s s s
+l l l l l.
+
+\fBWedge  24 	CTIO SCHMIDT WESTON TUBE SENSITOMETER.            \fR
+ 	MONO/MONO					
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge  48 	PALOMAR 48-INCH SCHMIDT STEP WEDGE.               \fR
+ 	MONO/MONO			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge  84 	OLD 84-INCH SPOT SENSITOMETER (1967)               \fR
+ 	MONO/MONO			
+
+.TE
+.TS
+l l s s s
+l l l l l.
+\fBWedge 101 	SPOT BOX 4, KEPT IN SCHOENING-S LAB.              \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 115 	MAYALL 4-M PF BEFORE 15APR74 (CHROME) [MP1-MP968] \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4770	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 117 	CTIO 4-METER P.F.                                 \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4770	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 118 	CTIO 4-METER CASSEGRAIN                           \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470	MONO/6900		
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 119	SPOT BOX 5, KEPT AT MAYALL 4-METER.               \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 120 	SPOT BOX 6, KEPT AT 2.1-METER.                    \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 121 	SPOT BOX 8, KEPT IN SCHOENING'S LAB.              \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470        
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 122 	SPOT BOX 7, AVAILABLE AT KPNO NIGHT ASST'S OFFICE \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470C			
+.TE
+.TS
+l l s s s
+l l l l l.
+\fBWedge 123 	MAYALL 4-M P.F. 15APR74 TO 21MAY74 [MP969-MP1051] \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4770	MONO/5200	MONO/5876        
+ 	MONO/6470        
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 129 	MAYALL 4-METER P.F. AFTER 21MAY74 [MP1052-->  ]   \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 130 	MAYALL 4-METER CASS CAMERA.                       \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4760	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 138 	TRAVELLING BOX AFTER 06JAN78.                    \fR
+ 	IIIAJ/UG2	IIAO/UG2	IIIAJ/*5113	IIAO/*5113       
+ 	IIAO/GG385	IIIAJ/CLEAR	IIIAJ/GG385	IIAD/GG495       
+ 	127/GG495	098/RG610	127/RG610	IVN/RG695       
+ 	MONO/4363	MONO/4770	MONO/5200	MONO/5876        
+ 	MONO/6470			
+
+.T&
+l l s s s
+l l l l l.
+\fBWedge 201 	TEN UCLA SPOTS (H. FORD, 10JAN78)                 \fR
+ 	MONO/MONO			
+.TE
diff --git a/noao/imred/dtoi/doc/dtoi.toc b/noao/imred/dtoi/doc/dtoi.toc
new file mode 100644
index 00000000..cd794321
--- /dev/null
+++ b/noao/imred/dtoi/doc/dtoi.toc
@@ -0,0 +1,34 @@
+.LP
+.DS C
+\fBTable of Contents\fR
+.DE
+.sp 3
+1.\h'|0.4i'\fBIntroduction\fP\l'|5.6i.'\0\01
+.sp
+2.\h'|0.4i'\fBDetermining the HD Curve Data \fP\l'|5.6i.'\0\01
+.br
+\h'|0.4i'2.1.\h'|0.9i'SPOTLIST\l'|5.6i.'\0\02
+.br
+\h'|0.4i'2.2.\h'|0.9i'DEMATCH\l'|5.6i.'\0\02
+.sp
+3.\h'|0.4i'\fBFitting the Curve\fP\l'|5.6i.'\0\02
+.br
+\h'|0.4i'3.1.\h'|0.9i'HDFIT\l'|5.6i.'\0\02
+.br
+\h'|0.9i'3.1.1.\h'|1.5i'Interactive Curve Fitting\l'|5.6i.'\0\03
+.sp
+4.\h'|0.4i'\fBApplying the Transform\fP\l'|5.6i.'\0\05
+.br
+\h'|0.4i'4.1.\h'|0.9i'HDTOI\l'|5.6i.'\0\05
+.sp
+5.\h'|0.4i'\fBAligning Related HD curves\fP\l'|5.6i.'\0\06
+.br
+\h'|0.4i'5.1.\h'|0.9i'HDSHIFT\l'|5.6i.'\0\06
+.sp
+6.\h'|0.4i'\fBTesting the Transformation Algorithm\fP\l'|5.6i.'\0\06
+.br
+\h'|0.4i'6.1.\h'|0.9i'SELFTEST\l'|5.6i.'\0\06
+.sp
+7.\h'|0.4i'\fBThe Wedgefile Database\fP\l'|5.6i.'\0\06
+.br
+\h'|0.4i'7.1.\h'|0.9i'Contents of the NOAO Wedgefile\l'|5.6i.'\0\07
diff --git a/noao/imred/dtoi/doc/hdfit.hlp b/noao/imred/dtoi/doc/hdfit.hlp
new file mode 100644
index 00000000..0b55137e
--- /dev/null
+++ b/noao/imred/dtoi/doc/hdfit.hlp
@@ -0,0 +1,79 @@
+.help hdfit Mar88 imred.dtoi
+.ih
+NAME
+hdfit -- fit characteristic curve to density, exposure data
+.ih
+USAGE
+hdfit database 
+.ih
+PARAMETERS
+.ls database
+Database[s] containing the density, log exposure information.
+.le
+.ls function = "power"
+Type of curve to fit; chosen from "power", "legendre", "chebyshev", 
+"spline1" or "spline3".  Abbreviations are permitted.
+.le
+.ls transform = "logopacitance"
+Transformation performed on the density prior to fitting.  Chosen from
+"none", "logopacitance", "k50" or "k75". 
+.le
+.ls weighting = "none"
+Weights can be assigned to the independent variable for fitting a curve.
+Choices are "none", "user" and "calculated".
+.le
+.ls order = 4
+Order of the fit.
+.le
+.ls interactive = yes
+Fit the data interactively?
+.le
+.ls device = "stdgraph"
+Interactive graphics device.
+.le
+.ls cursor = "stdgcur"
+Source of cursor input.
+.le
+.ih
+DESCRIPTION
+Task \fIhdfit\fR is used to fit a curve to density and log exposure
+values in preparation for transforming an image from density to intensity.
+The log exposure and density are read from \fBdatabase\fR.
+More than one database can be input,
+in which case one curve is fit to the combined data and the results
+written to each database in the list.
+
+Weights can be applied to the independent variable of the fit.
+Weights can be changed interactively, and are originally chosen from
+"none", "user" and "calculated".  A weights value can
+be calculated from the standard deviations, read from \fBdatabase\fR,
+as weight = (normalized density) / sdev.  If user weights are to be
+used, they are read from \fBdatabase\fR record "weights" as "wts_vals"
+entries.  
+
+When \fBinteractive\fR = yes, the HD curve is plotted and the cursor
+made available for interactively examining and altering the fit.
+The fitting function, transformation and order can be modified; data
+points can be added, deleted or edited.  Four choices of independent
+variable are available in \fBhdfit\fR by means of the parameter 
+\fBtransform\fR.  No transformation can take place, in which case
+the independent variable is the density.  Other choices are the log
+opacitance or a Kaiser transform with alpha = 0.50 or 0.75.  The
+default choice is to fit log exposure as a function of the log opacitance; 
+this is traditionally known as the Baker-Seidel function.
+.ih
+EXAMPLES
+.nf
+Using the defaults as starting parameters, interactively fit a curve to
+the data points in db1.
+
+	cl> hdfit db1 
+
+A sixth order power series function is fit in batch mode to the db1 data.
+
+	cl> hdfit db1 order=6 interactive-
+.fi
+.ih
+SEE ALSO
+spotlist, dematch, hdtoi
+.endhelp
diff --git a/noao/imred/dtoi/doc/hdshift.hlp b/noao/imred/dtoi/doc/hdshift.hlp
new file mode 100644
index 00000000..aaa59063
--- /dev/null
+++ b/noao/imred/dtoi/doc/hdshift.hlp
@@ -0,0 +1,50 @@
+.help hdshift Feb87 imred.dtoi
+.ih
+NAME
+hdshift - calculate and subtract zero point to align HD curves.
+.ih
+USAGE
+hdshift database
+.ih
+PARAMETERS
+.ls database
+Input list of databases containing density, exposure and fit information.
+.le
+.ih
+DESCRIPTION
+For each file in \fBdatabase\fR, procedure \fBhdshift\fR calculates and 
+subtracts a zero point shift to bring several related HD curves into
+alignment.  The individual shifts are calculated by elimination of the 
+first coefficient (Bevington, eqn 9-3):
+.nf
+                _      _      _               _
+           a0 = y - a1*X - a2*X**2 - ... - an*X**n
+
+.fi
+Here, the averages over y and X refer to individual \fBdatabase\fR averages; 
+the coefficients a1, ... an were previously calculated using data from all 
+\fBdatabase\fRs, in task \fIhdfit\fR, and stored in the database.  The
+a0 term is calculated individually for each database; this term represents
+the zero point shift in log exposure and will be different for each database.
+
+On output, the log exposure values in each \fBdatabase\fR have been 
+shifted to the zero point shift of the first database in the list.  The
+log exposure records are now aligned and it would be appropriate
+to run task \fIhdfit\fR on the modified \fBdatabase\fR list and
+determine the common solution.
+.ih
+EXAMPLES
+.nf
+
+Shift the curves in four databases to a common zero point.  
+
+	cl> hdshift db1,db2,db3,db4
+.fi
+.ih
+SEE ALSO
+hdfit, hdtoi
+.br
+"Averaging Photographic Characteristic Curves", John Kormendy, from
+"ESO Workshop on Two Dimensional Photometry", Edited by P. Crane and
+K.Kjar, p 69, (1980), an ESO Publication.
+.endhelp
diff --git a/noao/imred/dtoi/doc/hdtoi.hlp b/noao/imred/dtoi/doc/hdtoi.hlp
new file mode 100644
index 00000000..bf4355f0
--- /dev/null
+++ b/noao/imred/dtoi/doc/hdtoi.hlp
@@ -0,0 +1,88 @@
+.help hdtoi May88 imred.dtoi
+.ih
+NAME
+hdtoi -- transform images according to hd curve
+.ih
+USAGE
+hdtoi input output database
+.ih
+PARAMETERS
+.ls input
+List of images to be transformed.
+.le
+.ls output
+List of output image names.
+.le
+.ls database
+Name of text database describing HD curve.
+.le
+.ls fog = ""
+Value of fog level, read from database if unspecified.
+.le
+.ls option = "mean"
+Option for calculating fog density when \fBfog\fR is a file list, can be
+either "mean" or "median".
+.le
+.ls sigma = 3.0
+If \fBfog\fR is a file name, and \fBoption\fR = "mean", the mean fog density
+is iteratively calculated using this rejection criteria.
+.le
+.ls floor = 0.0
+Value assigned to levels below fog, can be either 0.0 or -1.0.  
+.le
+.ls ceiling = 30000.
+The final intensities are scaled to this value, such that a saturated
+input density equals \fBceiling\fR on output.
+.le
+.ls datatype = "r"
+Datatype of output image pixels.
+.le
+.ls verbose = yes
+Print log of processing to STDOUT.
+.le
+.ih
+DESCRIPTION
+Task \fIhdtoi\fR transforms one image to another as described by the 
+\fBdatabase\fR.  There is only one HD curve per run; the same 
+transformation is applied to all input images.
+
+The fog value can be obtained in three ways: read from the database, read
+as a floating point number, or calculated from a list of fog images.  If 
+parameter \fBfog\fR is not specified, the fog value is read from 
+\fBdatabase\fR.  If \fBfog\fR is specified, it can be entered
+as either a floating point number or as a list of file names.  If the
+value cannot be read as a number, it is assumed to be a file name.  In that
+case, the density of each file in the fog list is calculated and the 
+average of these values is subtracted from \fBinput\fR before processing.
+The algorithm used to calculate the fog density is selected by the
+\fBoption\fR parameter, and is either a "mean" or "median" calculation.
+The fog density can be the mean value after pixels more than the specified
+number of sigma have been rejected, or the median value of all the fog spot
+pixels.
+
+The fog value is subtracted from the input image before the transformation
+takes place.  It is possible that some density values will fall below
+the fog level; these values are handled in one of two ways.  Values
+below the fog value are set equal to 0.0 when \fBfloor\fR = 0.0.  If 
+\fBfloor\fR = -1.0, the resulting intensity = -1 * intensity (abs (value)).
+
+A scaling factor is applied to the final intensities, as typically
+they will be < 1.0.  The \fBceiling\fR parameter is used to specify what
+value a saturated density is transformed to; all intensities are scaled
+to this upper limit.  The precision of the transformation is unaffected by 
+this parameter, although caution must be used if the output image pixel
+type is an integer.  The user is responsible for choosing
+a \fBceiling\fR that avoids the truncation of significant digits.
+.ih
+EXAMPLES
+Convert three density images to intensity images as described in database db1.
+
+	cl> hdtoi denin* intim1,intim2,intim3 db1
+.ih
+TIME REQUIREMENTS
+Task \fBhdtoi\fR requires 20 cpu seconds to transform a 512 square image, with
+a 12 bit data range, on a VAX 750
+.ih
+SEE ALSO
+spotlist, dematch, hdfit
+.endhelp
diff --git a/noao/imred/dtoi/doc/selftest.hlp b/noao/imred/dtoi/doc/selftest.hlp
new file mode 100644
index 00000000..329c9099
--- /dev/null
+++ b/noao/imred/dtoi/doc/selftest.hlp
@@ -0,0 +1,81 @@
+.help selftest Feb87 imred.dtoi
+.ih
+NAME
+selftest -- test routine to verify \fIdtoi\fR transformation
+.ih
+USAGE
+selftest nbits
+.ih
+PARAMETERS
+.ls nbits = 12
+Dymanic range of data to test.
+.le
+.ls device = "stdgraph" 
+Plotting device for graphical output.
+.le
+.ls verbose = no
+A table of density, intensity values is printed if \fBverbose\fR = yes.
+.le
+.ls ceiling = 30000.
+Maximum intensity to output.
+.le
+.ls max_raw = 0
+The maximum raw data value.  Needed only if \fInbits\fR equals something
+other than 12, 15 or 0.
+.le
+.ls scale = 0.0
+The raw data value to density scale value.  Needed only if \fInbits\fR
+equals something other than 12, 15, or 0.
+.le
+
+.ih
+DESCRIPTION
+Task \fIselftest\fR is a test program for the \fIdtoi\fR package.  Its 
+output can be examined to see if numerical errors are introduced during
+the density to intensity transformation.  It also evaluates truncation
+errors produced when an output image with integer pixels is written.  
+
+Many different PDS setups can be investigated with task \fBselftest\fR.
+Setting parameter \fInbits\fR = 12
+indicates PDS format data, with data range 0 to 3071.  Setting \fInbits\fR = 15 
+indicates FITS format data, with data range 0 to 24575.  The special value of
+\fInbits\fR = 0 means a small test data range from 1 to 144 is investigated.
+If any other value of \fInbits\fR is entered, the user is queried for the
+max raw data values and the raw data to density scaling factor.
+
+An intensity vector is generated from a density vector in two different ways.  
+The first method uses the density vector and known coefficients to compute
+the intensity.  The second method uses the curfit package to generate a
+look up table of intensities as done in task \fBHDTOI\fR.  The residual
+of the two intensity vectors is plotted.  Ideally, the difference between
+the 'known' intensities and 'calculated' intensities is zero.
+
+The second plot output by \fBselftest\fR shows intensity as a function
+of density.  Two lines are overplotted; integer intensity versus density
+and real intensity versus density.  Because truncation errors are most
+pronounced at low density values, the plot covers only the lowest 5%
+of the density range.  The user should investigate the plot with the
+cursor zoom and expand capabilities to determine if truncation errors
+are significant.
+
+In verbose mode, \fBselftest\fR produced a three column table of raw
+data value, density and calculated intensity. 
+
+.ih
+EXAMPLES
+
+.nf
+Run task selftest for 12 bit data with plots appearing on the terminal.
+
+	cl> selftest
+
+.fi
+Run selftest in verbose mode, spooling the output to file 'ditable'.  This
+file is then run through the 'fields' task to extract the density and intensity
+columns which are piped to plot.  The results in a plot of the look up table.
+.nf
+
+	cl> selftest ver+ > ditable
+	cl> fields ditable 2,3 | graph xlab=Density ylab=Intensity
+.fi
+.endhelp
diff --git a/noao/imred/dtoi/doc/splotlist.hlp b/noao/imred/dtoi/doc/splotlist.hlp
new file mode 100644
index 00000000..43b3f223
--- /dev/null
+++ b/noao/imred/dtoi/doc/splotlist.hlp
@@ -0,0 +1,81 @@
+.help spotlist May88 imred.dtoi
+.ih
+NAME
+spotlist -- calculate densities of calibration spots 
+.ih
+USAGE
+spotlist spots fogs database
+.ih
+PARAMETERS
+.ls spots
+List of image files containing the calibration data.
+.le
+.ls fogs
+List of image files containing fog spots.
+.le
+.ls database
+Name for output database.
+.le
+.ls scale = 0.00151			# (4.65 / 3071.)
+The scale factor to convert values in the image files to densities, such
+that scale = density / input_value.
+.le
+.ls maxad = 3071
+The maximum A/D value, that is, the input value corresponding to a
+saturated pixel.
+.le
+.ls option= "mean"
+Option for calculating densities can be either "mean" or "median".
+.le
+.ls sigma = 3.0
+Rejection criteria for iteratively calculating mean density.
+.le
+.ih
+DESCRIPTION
+Task \fIspotlist\fR reads calibration spot images and calculates their
+density and standard deviation.  Three records are entered in the
+database: density, standard deviation and number of unrejected pixels.
+Each record contains as many entries as calibration spots.
+
+All input values are multiplied by the \fBscale\fR parameter to convert
+them to densities.  The value of \fBscale\fR is not critical to the 
+reductions, it is provided so that \fIspotlist\fR output can be in the 
+familiar units of density.  The default value of \fBscale\fR is correct
+for NOAO PDS data written to a PDS format tape.  If a FITS format tape was 
+written, \fBscale\fR = 0.0001893.  These values are appropriate for the PDS 
+with its new 15-bit logarithmic A to D converter.  The value of \fBscale\fR
+used is also entered in the database.
+
+Parameter \fBmaxad\fR is the integer input value that represents a
+saturated pixel.  This value is used by \fIspotlist\fR to accurately
+calculate the density of a saturated pixel, which is then entered in the
+database.  This value of "maxden" will later be used by task \fIhdfit\fR
+to normalize the independent variable vector, and by task \fIhdtoi\fR to
+scale the intensity range precisely to a user specified value.
+
+A fog level is calculated from image \fBfogs\fR, and entered into
+the database file.  If more than one image is given for \fBfogs\fR, 
+a single fog value is calculated from all fog pixels.  The fog level
+is calculated but not subtracted in this procedure.  The fog images to be 
+averaged should be the same size.  An entry for the fog level is made
+in the database.
+
+The \fBspots\fR files are assumed to be ordered such that they are either
+monotonically increasing or decreasing in density, with no omitted spots
+between the first and last measured.  The calculated density can be
+calculated two ways; the algorithm used is selected by the \fBoption\fR
+parameter.  The density is either the mean spot value after pixels more
+than the specified number of sigma from the mean value have been rejected,
+or the median value of all the spot pixels.  
+.ih
+EXAMPLES
+Calculate mean densities of calibration spots which had previously been
+read in from a FITS format tape.  The database "db1" will be created.
+
+.nf
+	cl> spotlist spots* fogspot db1 scale=1.893e-4
+.fi
+.ih
+SEE ALSO
+dematch, hdfit, hdtoi
+.endhelp
diff --git a/noao/imred/dtoi/dtoi.cl b/noao/imred/dtoi/dtoi.cl
new file mode 100644
index 00000000..3dd6c42b
--- /dev/null
+++ b/noao/imred/dtoi/dtoi.cl
@@ -0,0 +1,16 @@
+#{
+# DTOI package -- Density to Intensity Transformation Package.  (Feb87)
+
+set	dtoi		= "iraf$noao/imred/dtoi/"
+
+package	dtoi
+
+task	dematch,
+	hdfit,
+	hdtoi,
+	hdshift,
+	selftest,
+	spotlist	= "dtoi$x_dtoi.e"
+
+
+clbye()
diff --git a/noao/imred/dtoi/dtoi.hd b/noao/imred/dtoi/dtoi.hd
new file mode 100644
index 00000000..af554ac2
--- /dev/null
+++ b/noao/imred/dtoi/dtoi.hd
@@ -0,0 +1,11 @@
+# Help directory for the dtoi package.
+
+$doc		= "iraf$noao/imred/dtoi/doc/"
+
+dematch		hlp=doc$dematch.hlp,		src=dematch.x
+hdfit		hlp=doc$hdfit.hlp,		src=hdfit.x
+hdtoi		hlp=doc$hdtoi.hlp,		src=hdtoi.x
+hdshift		hlp=doc$hdshift.hlp,		src=hdshift.x
+spotlist	hlp=doc$spotlist.hlp,		src=spotlist.x
+selftest	hlp=doc$selftest.hlp,		src=selftest.x
+revisions	sys=Revisions
diff --git a/noao/imred/dtoi/dtoi.men b/noao/imred/dtoi/dtoi.men
new file mode 100644
index 00000000..f22d92eb
--- /dev/null
+++ b/noao/imred/dtoi/dtoi.men
@@ -0,0 +1,6 @@
+	dematch - Match a list of density values to exposure values
+	  hdfit - Fit a curve to density, log exposure values
+	hdshift - Align related HD curves
+          hdtoi - Apply DTOI transformation to density image
+       selftest - Self test program to check DTOI transformation
+       spotlist - Generate a list of calibration spot values
diff --git a/noao/imred/dtoi/dtoi.par b/noao/imred/dtoi/dtoi.par
new file mode 100644
index 00000000..494aa73c
--- /dev/null
+++ b/noao/imred/dtoi/dtoi.par
@@ -0,0 +1,2 @@
+version,s,h,"February 13, 1987"
+mode,s,h,ql
diff --git a/noao/imred/dtoi/hd_aravr.x b/noao/imred/dtoi/hd_aravr.x
new file mode 100644
index 00000000..3376264b
--- /dev/null
+++ b/noao/imred/dtoi/hd_aravr.x
@@ -0,0 +1,50 @@
+define	MAX_ITERATIONS	10
+include	<mach.h>
+
+# HD_ARAV -- Compute the mean and standard deviation of a sample array by
+# iteratively rejecting points further than KSIG from the mean.  If the
+# value of KSIG is given as 0.0, a cutoff value will be automatically
+# calculated from the standard deviation and number of points in the sample.
+# The number of pixels remaining in the sample upon termination is returned
+# as the function value.
+#
+# A max_iterations parameter was added to prevent the rejection scheme
+# from oscillating endlessly for nearly saturated pixels.  This is the
+# only difference between the vops procedure and hd_aravr.  (ShJ 5/88)
+
+int procedure hd_aravr (a, npix, mean, sigma, ksig)
+
+real	a[ARB]			# input data array
+real	mean, sigma, ksig, deviation, lcut, hcut, lgpx
+int	npix, niter, ngpix, old_ngpix, awvgr()
+
+begin
+	lcut = -MAX_REAL				# no rejection to start
+	hcut =  MAX_REAL
+	ngpix = 0
+	niter = 0
+
+	# Iteratively compute mean, sigma and reject outliers until no
+	# more pixels are rejected, or until there are no more pixels,
+	# or until the maximum iterations limit is exceeded.
+
+	repeat {
+	    niter = niter + 1
+	    old_ngpix = ngpix
+	    ngpix = awvgr (a, npix, mean, sigma, lcut, hcut)
+	    if (ngpix <= 1)
+		break
+
+	    if (ksig == 0.0) {				# Chauvenet's relation
+		lgpx = log10 (real(ngpix))
+		deviation = (lgpx * (-0.1042 * lgpx + 1.1695) + .8895) * sigma
+	    } else
+		deviation = sigma * abs(ksig)
+
+	    lcut = mean - deviation			# compute window
+	    hcut = mean + deviation
+
+	} until ((old_ngpix == ngpix) || (niter > MAX_ITERATIONS))
+
+	return (ngpix)
+end
diff --git a/noao/imred/dtoi/hdfit.par b/noao/imred/dtoi/hdfit.par
new file mode 100644
index 00000000..7f42aa2a
--- /dev/null
+++ b/noao/imred/dtoi/hdfit.par
@@ -0,0 +1,9 @@
+# Cl parameters for task hdfit are:
+database,f,a,,,,List of database names
+function,s,h,"power",,,Fitting function
+transform,s,h,"logopacitance",,,Independent variable transformation
+order,i,h,4,,,Initial order (ncoeff) of fit
+interactive,b,h,y,,,Interactive fitting flag
+device,s,h,"stdgraph",,,Name of interactive graphics device
+weighting,s,h,"none",,,Type of weighting
+cursor,*gcur,h,"",,,Source of cursor input
diff --git a/noao/imred/dtoi/hdfit.x b/noao/imred/dtoi/hdfit.x
new file mode 100644
index 00000000..04a82f6e
--- /dev/null
+++ b/noao/imred/dtoi/hdfit.x
@@ -0,0 +1,364 @@
+include	<math/curfit.h>
+include	<imhdr.h>
+include	<fset.h>
+include	<mach.h>
+include	<ctype.h>
+include	<error.h>
+include	<pkg/gtools.h>
+include	<pkg/xtanswer.h>
+include	"hdicfit/hdicfit.h"
+
+# T_HDFIT -- Fit a curve.  This task fits a characteristic
+# curve to density and log exposure data read from an input
+# database.  The interactive curve fitting package is used.
+# The database is updated to contain the values of the fit 
+# necessary to reinitialize the curfit package for performing 
+# the transformation in hdtoi. 
+
+procedure t_hdfit ()
+
+pointer	sp, fcn, device, db, gt, exp, wts, save, trans, weight
+pointer	dbfile, ic, den, errs
+int	db_list, order, interactive, nsave, nvals, wt_type, update
+real	ref_fog, real_fog
+double	fog, maxden
+
+pointer	ddb_map(), gt_init()
+bool	clgetb(), fp_equalr(), fp_equald()
+int	clpopni(), clgeti(), strncmp(), clgfil(), hd_fit()
+real	ic_getr()
+
+begin
+	call smark  (sp)
+	call salloc (fcn, SZ_FNAME, TY_CHAR)
+	call salloc (device, SZ_FNAME, TY_CHAR)
+	call salloc (trans, SZ_FNAME, TY_CHAR)
+	call salloc (weight,  SZ_FNAME, TY_CHAR)
+	call salloc (dbfile, SZ_FNAME, TY_CHAR)
+
+	# Get cl parameters
+	db_list = clpopni ("database")
+	call clgstr ("function", Memc[fcn], SZ_FNAME)
+	call clgstr ("transform", Memc[trans], SZ_FNAME)
+	order = clgeti ("order")
+
+	# Decide which type of weighting the user wants
+	call clgstr ("weighting", Memc[weight], SZ_FNAME)
+	if (strncmp (Memc[weight], "none", 1) == 0)
+	    wt_type = WT_NONE
+	else if (strncmp (Memc[weight], "user", 1) == 0)
+	    wt_type = WT_USER
+	else if (strncmp (Memc[weight], "calc", 1) == 0)
+	    wt_type = WT_CALC
+	else
+	    call error (0, "Unrecognized weighting type")
+
+	# Initialize interactive curve fitting package.
+	gt = gt_init ()
+	call ic_open (ic)
+	call ic_pstr (ic, "function", Memc[fcn])
+	call ic_pstr (ic, "transform", Memc[trans])
+	call ic_puti (ic, "order", order)
+	call ic_pstr (ic, "ylabel", "Log Exposure")
+	call ic_pkey (ic, 5, 'y', 'u')
+
+	if (clgetb ("interactive")) {
+	    interactive = YES
+	    call clgstr ("device", Memc[device], SZ_FNAME)
+	} else {
+	    interactive = NO
+	    call strcpy ("", Memc[device], SZ_FNAME)
+	}
+
+	# Read information from each dlog file; accumulate number of values.
+	# The density (not fog subtracted) is returned.  The density values
+	# are also sorted in increasing order.
+
+	call hd_rdloge (db_list, exp, den, wts, errs, nvals, fog, maxden,
+	    wt_type)
+	if (nvals == 0)
+	    call error (1, "T_HDFIT: No data values in sample")
+
+	call hd_sdloge (Memd[den], Memd[exp], Memd[wts], Memd[errs], nvals)
+	call ic_putr (ic, "fog", real(fog))
+	ref_fog = real (fog)
+	call ic_putr (ic, "rfog", ref_fog)
+
+	# Initialize the dtoi/icgfit interface.
+	if (fp_equald (maxden, 0.0D0))
+	    call error (1, "Saturated pixel density not initialized")
+	call hdic_init (Memd[den], nvals, maxden)
+
+	update = hd_fit (ic, gt, Memd[den], Memd[exp], Memd[wts], Memd[errs],
+	    nvals, save, nsave, Memc[device], interactive)
+	     
+	if (update == YES) {
+	    # Record fit information in (each) database
+	    call ic_gstr  (ic, "function", Memc[fcn], SZ_FNAME)
+	    call ic_gstr  (ic, "transform", Memc[trans], SZ_FNAME)
+	    real_fog = ic_getr (ic, "fog")
+
+	    while (clgfil (db_list, Memc[dbfile], SZ_FNAME) != EOF) {
+	        db = ddb_map (Memc[dbfile], APPEND)
+	        call ddb_ptime (db)
+	        # Add new fog record if it was changed interactively.
+	        if (!fp_equalr (real_fog, ref_fog)) {
+		    call ddb_prec (db, "fog")
+		    call ddb_putr (db, "density", real_fog)
+	        }
+
+	        call ddb_prec (db, "cv")
+	        call ddb_pad  (db, "save", Memd[save], nsave)
+	        call ddb_pstr (db, "function", Memc[fcn])
+	        call ddb_pstr (db, "transformation", Memc[trans])
+
+	        call ddb_unmap (db)
+	    }
+	}
+
+	call ic_closed (ic)
+	call mfree (save, TY_DOUBLE)
+	call mfree (den,  TY_DOUBLE)
+	call mfree (exp,  TY_DOUBLE)
+	call mfree (wts , TY_DOUBLE)
+	call mfree (errs, TY_DOUBLE)
+
+	call gt_free (gt)
+	call clpcls (db_list)
+	call sfree (sp)
+end
+
+
+# HD_RLOGE -- Read log exposure, density and weights from a single dloge
+# database file.  Pointers to the three arrays are returned as arguments.
+# The number of values accumulated is returned also; note that the value
+# of nvals is changed upon return; it should not be given as a constant.
+# If more than one database is being read (as in HDSHIFT applications),
+# the density ABOVE fog is returned, and the reference fog values set = 0.0
+# The maximum density, the density of a saturated pixel, is read from the
+# first databas and returned.
+
+procedure hd_rdloge (db_list, exp, den, wts, errs, nvals, fog, maxden, wt_type)
+
+int	db_list			# File descriptor for data base file
+pointer	exp			# Pointer to exposure array - returned
+pointer	den			# Pointer to density array - returned
+pointer	wts			# Pointer to weights array - returned
+pointer	errs			# Pointer to std deviation array - returned
+int	nvals			# Number of data pairs read - returned
+double	fog			# Value of fog read from database - returned
+double  maxden			# Maximum density, read from db - returned
+int	wt_type			# Type of weighting
+
+pointer	db
+bool	sdevrec
+int	buflen, off, rec
+int	nden, nexp, nwts, nspots, nerrs, nfiles
+char	dloge[SZ_FNAME]
+
+pointer	ddb_map()
+bool	fp_equald()
+int	ddb_locate(), ddb_geti(), clgfil(), imtlen()
+real	ddb_getr()
+errchk	ddb_locate, ddb_gad, malloc, ddb_map, ddb_unmap
+
+begin
+	nvals = 0
+	off = 0
+	buflen = NSPOTS
+	nfiles = imtlen (db_list)
+	maxden = 0.0D0
+
+	# Dynamically allocate memory for arrays; it can be increased later.
+	call malloc (exp, buflen, TY_DOUBLE)
+	call malloc (den, buflen, TY_DOUBLE)
+	call malloc (wts, buflen, TY_DOUBLE)
+	call malloc (errs, buflen, TY_DOUBLE)
+
+	while (clgfil (db_list, dloge, SZ_FNAME) != EOF) {
+	    iferr (db = ddb_map (dloge, READ_ONLY)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Get fog value to be subtracted from density
+	    rec = ddb_locate (db, "fog")
+	    fog = double (ddb_getr (db, rec, "density"))
+
+	    # Get density array
+	    rec = ddb_locate (db, "density")
+	    nden = ddb_geti (db, rec, "den_val")
+
+	    call ddb_gad (db, rec, "den_val", Memd[den+off], nden, nden)
+	    if (nfiles > 1)
+		call asubkd (Memd[den+off], fog, Memd[den+off], nden)
+
+	    # Get log exposure array
+	    rec = ddb_locate (db, "exposure")
+	    nexp = ddb_geti (db, rec, "log_exp")
+	    call ddb_gad (db, rec, "log_exp", Memd[exp+off], nexp, nexp)
+
+	    # Get saturated pixel density if not already set
+	    if (fp_equald (maxden, 0.0D0)) {
+	        iferr (rec = ddb_locate (db, "common")){
+		    ;
+		} else
+	            maxden = double (ddb_getr (db, rec, "maxden"))
+	    }
+
+	    # Get std deviation array
+	    sdevrec = true
+	    iferr {
+	        rec = ddb_locate (db, "standard deviation")
+	        nerrs = ddb_geti (db, rec, "sdev_val")
+	        call ddb_gad (db, rec, "sdev_val", Memd[errs+off], nerrs, nerrs)
+	    } then {
+		call erract (EA_WARN)
+		call eprintf ("Marker type '[ve]bar' can only show weights\n")
+		call amovkd (0.0D0, Memd[errs+off], nden)
+		sdevrec = FALSE
+	    }
+
+	    if (wt_type == WT_CALC) {
+	        if (sdevrec) {
+		    iferr {
+                        nspots = min (nden, nexp, nerrs)
+                        call adivd (Memd[den+off], Memd[errs+off], 
+			    Memd[wts+off], nspots)
+	            } then {
+		        call erract (EA_WARN)
+		        call eprintf ("All weights set to 1.0\n")
+	                call amovkd (double (1.0), Memd[wts+off], nspots)
+	            }
+	        } else {
+		    nspots = min (nden, nexp)	
+		    call eprintf ("No sdev record; All weights set to 1.0\n")
+	            call amovkd (double (1.0), Memd[wts+off], nspots)
+	        }
+	    }
+
+	    if (wt_type == WT_USER) {
+	        # WT_USER: fill "user" weights array
+	        iferr {
+	            rec = ddb_locate (db, "weight")
+	            nwts = ddb_geti (db, rec, "wts_val")
+		    nspots = min (nden, nexp, nwts)
+	            call ddb_gad (db, rec, "wts_val", Memd[wts+off], nwts, nwts)
+	        } then {
+	            # Users weights can't be found. Set weights array to 1.0's.
+		    call erract (EA_WARN)
+		    call eprintf ("All weights set to 1.0\n")
+	            nspots = min (nden, nexp)
+	            call amovkd (double (1.0), Memd[wts+off], nspots)
+		}
+	    }
+
+	    if (wt_type == WT_NONE) {
+	        # WT_NONE: fill "none" weights array
+	        nspots = min (nden, nexp)
+		call amovkd (double (1.0), Memd[wts+off], nspots)
+	    }
+    
+	    # Increment number of counts; reallocate memory if necessary.
+	    nvals = nvals + nspots
+	    off = off + nspots
+
+	    if (nvals > buflen) {
+		buflen = buflen + NSPOTS
+		call realloc (exp, buflen, TY_DOUBLE)
+		call realloc (den, buflen, TY_DOUBLE)
+		call realloc (wts, buflen, TY_DOUBLE)
+		call realloc (errs, buflen, TY_DOUBLE)
+	    }
+
+	    call ddb_unmap (db)
+	}
+
+	if (nfiles > 1)
+	    fog = 0.0D0
+	call clprew (db_list)
+end
+
+
+# HD_SLOGE -- Sort the log exposure, density and weight information in order
+# of increasing exposure value.  The sorting is done is place.  The three
+# data values are assummed matched on input, that is, exposure[i] matches
+# density[i] with weight[i] for all array entries.
+
+procedure hd_sdloge (density, exposure, weights, errors,  nvals)
+
+double	density[nvals]		# Density array
+double	exposure[nvals]		# Exposure array
+double	weights[nvals]		# Weights array
+double	errors[nvals]		# Standard deviation array
+int	nvals			# Number of values in data arrays
+
+int	i, j
+double	temp
+define	swap	{temp=$1;$1=$2;$2=temp}
+
+begin
+	# Bubble sort - inefficient, but sorting is done only once on
+	# an expected small sample size (16 pts typically). 
+
+	for (i = nvals; i > 1; i = i - 1)
+	    for (j = 1; j < i; j = j + 1) 
+		if (density [j] > density [j+1]) {
+
+		    # Out of order; exchange values
+		    swap (exposure[j], exposure[j+1])
+		    swap ( density[j],  density[j+1])
+		    swap ( weights[j],  weights[j+1])
+		    swap (  errors[j],   errors[j+1])
+		}
+end
+
+
+# HD_FIT -- Fit the curve to input density, exposure and weight values.
+# The fit can be performed interactively or not. 
+
+int procedure hd_fit (ic,
+	gt, den, exp, wts, errs, nvals, save, nsave, dev, interact)
+
+pointer	ic
+pointer	gt			# Graphics tools pointer
+double	den[ARB]		# Density values
+double	exp[ARB]		# Exposure values
+double	wts[ARB]		# Weight array
+double	errs[ARB]		# Standard deviation array
+int	nvals			# Number of data pairs to fit
+pointer	save			# ??
+int	nsave			# ??
+char	dev[SZ_FNAME]		# Interactive graphics device
+int	interact		# Flag for interactive graphics
+
+pointer	gp, cv, sp, x, dum
+pointer	gopen()
+int	update, dcvstati()
+errchk	malloc, gopen
+
+begin
+	if (interact == YES) {
+	    gp = gopen (dev, NEW_FILE, STDGRAPH)
+	    call icg_fitd (ic, gp, "cursor", gt, cv, den, exp, wts, errs, nvals)
+	    call gclose (gp)
+	    update = IC_UPDATE(ic)
+
+	} else {
+	    # Do fit non-interactively
+	    call smark (sp)
+	    call salloc (x, nvals, TY_DOUBLE)
+	    call salloc (dum, nvals, TY_INT)
+	    call hdic_transform (ic, den, wts, Memd[x], wts, Memi[dum], nvals)
+	    call ic_fitd (ic, cv, Memd[x], exp, wts, nvals, YES, YES, YES, YES)
+	    call sfree (sp)
+	    update = YES
+	}
+
+	nsave = (dcvstati (cv, CVORDER)) + 7
+	call malloc (save, nsave, TY_DOUBLE)
+	call dcvsave (cv, Memd[save])
+	call dcvfree (cv)
+
+	return (update)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdic.com b/noao/imred/dtoi/hdicfit/hdic.com
new file mode 100644
index 00000000..959df10e
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdic.com
@@ -0,0 +1,6 @@
+# Common block for hdtoi package/ icgfit package interface.
+
+int	nraw
+double	maxden
+pointer	den, big_den
+common	/raw/maxden, den, big_den, nraw
diff --git a/noao/imred/dtoi/hdicfit/hdicadd.x b/noao/imred/dtoi/hdicfit/hdicadd.x
new file mode 100644
index 00000000..1fae152d
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicadd.x
@@ -0,0 +1,47 @@
+include	"hdicfit.h"
+
+# HDIC_ADDPOINT -- Add a density, exposure, weights point into the sample.
+
+procedure hdic_addpoint (ic,
+	nden, nexp, nwts, den, y, wts, userwts, x, wd, sdev, npts)
+
+pointer	ic		# Pointer to ic data structure
+real	nden		# New density value to be added
+real	nexp		# New exposure value to be added
+double	nwts		# New weight value to be added
+pointer	den		# Pointer to existing density array
+pointer	y		# Pointer to existing exposure array
+pointer	wts		# Pointer to existing wts array
+pointer	userwts		# Pointer to existing userwts array
+pointer	x		# Pointer to existing array of ind variables
+pointer	wd		# Pointer to flag array for deletion reasons
+pointer	sdev		# Pointer to standard deviation array
+int	npts		# Number of points, incremented on output
+
+begin
+	npts = npts + 1
+
+	call realloc (den,     npts, TY_DOUBLE)
+	call realloc (y,       npts, TY_DOUBLE)
+	call realloc (wts,     npts, TY_DOUBLE)
+	call realloc (userwts, npts, TY_DOUBLE)
+	call realloc (x,       npts, TY_DOUBLE)
+	call realloc (wd,      npts, TY_INT)
+	call realloc (sdev,    npts, TY_DOUBLE)
+
+	Memd[den+npts-1] = double (nden)
+	Memd[y  +npts-1] = double (nexp)
+	Memd[wts+npts-1] = double (nwts)
+	Memd[userwts+npts-1] = double (nwts)
+	Memi[wd +npts-1] = NDELETE
+	Memd[sdev+npts-1] = ADDED_PT
+
+	# Sort the data and then update the reference vector.
+	call hdic_sort (Memd[den], Memd[y], Memd[wts], Memd[userwts], 
+	    Memi[wd], Memd[sdev],  npts)
+	call hdic_init (Memd[den], npts, Memd[den+npts-1])
+
+	IC_NEWX(ic) = YES
+	IC_NEWY(ic) = YES
+	IC_NEWWTS(ic) = YES
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicclean.x b/noao/imred/dtoi/hdicfit/hdicclean.x
new file mode 100644
index 00000000..6c838123
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicclean.x
@@ -0,0 +1,94 @@
+include <pkg/rg.h>
+include	"hdicfit.h"
+
+# IC_CLEAN -- Replace rejected points by the fitted values.
+
+procedure ic_cleand (ic, cv, x, y, w, npts)
+
+pointer	ic				# ICFIT pointer
+pointer	cv				# Curfit pointer
+double	x[ARB]				# Ordinates
+double	y[ARB]				# Abscissas
+double	w[ARB]				# Weights
+int	npts				# Number of points
+
+int	i, nclean, newreject
+pointer	sp, xclean, yclean, wclean
+double	dcveval()
+
+begin
+	if ((IC_LOW(ic) == 0.) && (IC_HIGH(ic) == 0.))
+	    return
+
+	# If there has been no subsampling and no sample averaging then the
+	# IC_REJPTS(ic) array already contains the rejected points.
+
+	if (npts == IC_NFIT(ic)) {
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+		        y[i] = dcveval (cv, x[i])
+	        }
+	    }
+
+	# If there has been no sample averaging then the rejpts array already
+	# contains indices into the subsampled array.
+
+	} else if (abs(IC_NAVERAGE(ic)) == 1) {
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[IC_YFIT(ic)+i-1] =
+			    dcveval (cv, Memd[IC_XFIT(ic)+i-1])
+		}
+	    }
+	    call rg_unpackd (IC_RG(ic), Memd[IC_YFIT(ic)], y)
+
+	# Because ic_fit rejects points from the fitting data which
+	# has been sample averaged the rejpts array refers to the wrong data.
+	# Do the cleaning using ic_deviant to find the points to reject.
+
+	} else if (RG_NPTS(IC_RG(ic)) == npts) {
+	    call amovki (NO, Memi[IC_REJPTS(ic)], npts)
+	    call ic_deviantd (cv, x, y, w, Memi[IC_REJPTS(ic)], npts,
+		IC_LOW(ic), IC_HIGH(ic), IC_GROW(ic), NO, IC_NREJECT(ic),
+		newreject)
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			y[i] = dcveval (cv, x[i])
+		}
+	    }
+
+	# If there is subsampling then allocate temporary arrays for the
+	# subsample points.
+
+	} else {
+	    call smark (sp)
+
+	    nclean = RG_NPTS(IC_RG(ic))
+	    call salloc (xclean, nclean, TY_DOUBLE)
+	    call salloc (yclean, nclean, TY_DOUBLE)
+	    call salloc (wclean, nclean, TY_DOUBLE)
+
+	    call rg_packd (IC_RG(ic), x, Memd[xclean])
+	    call rg_packd (IC_RG(ic), y, Memd[yclean])
+	    call rg_packd (IC_RG(ic), w, Memd[wclean])
+	    call amovki (NO, Memi[IC_REJPTS(ic)], npts)
+
+	    call ic_deviantd (cv, Memd[xclean], Memd[yclean],
+		Memd[wclean], Memi[IC_REJPTS(ic)], nclean, IC_LOW(ic),
+		IC_HIGH(ic), IC_GROW(ic), NO, IC_NREJECT(ic), newreject)
+
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[yclean+i-1] = dcveval (cv, Memd[xclean+i-1])
+		}
+	    }
+	    call rg_unpackd (IC_RG(ic), Memd[yclean], y)
+
+	    call sfree (sp)
+	}
+end
+
diff --git a/noao/imred/dtoi/hdicfit/hdicdeviant.x b/noao/imred/dtoi/hdicfit/hdicdeviant.x
new file mode 100644
index 00000000..0c6f8f75
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicdeviant.x
@@ -0,0 +1,116 @@
+include	<mach.h>
+include	<math/curfit.h>
+
+# IC_DEVIANT -- Find deviant points with large residuals from the fit
+# and reject from the fit.
+#
+# The sigma of the fit residuals is calculated.  The rejection thresholds
+# are set at +-reject*sigma.  Points outside the rejection threshold are
+# recorded in the reject array.
+
+procedure ic_deviantd (cv, x, y, w, rejpts, npts, low_reject, high_reject,
+	grow, refit, nreject, newreject)
+
+pointer	cv				# Curve descriptor
+double	x[ARB]				# Input ordinates
+double	y[ARB]				# Input data values
+double	w[ARB]				# Weights
+int	rejpts[ARB]			# Points rejected
+int	npts				# Number of input points
+real	low_reject, high_reject		# Rejection thresholds
+real	grow				# Rejection radius
+int	refit				# Refit the curve?
+int	nreject				# Number of points rejected
+int	newreject			# Number of new points rejected
+
+int	i, j, i_min, i_max, ilast
+double	sigma, low_cut, high_cut, residual
+pointer	sp, residuals
+
+begin
+	# If low_reject and high_reject are zero then simply return.
+	if ((low_reject == 0.) && (high_reject == 0.))
+	    return
+
+	# Allocate memory for the residuals.
+	call smark (sp)
+	call salloc (residuals, npts, TY_DOUBLE)
+
+	# Compute the residuals.
+	call dcvvector (cv, x, Memd[residuals], npts)
+	call asubd (y, Memd[residuals], Memd[residuals], npts)
+
+	# Compute the sigma of the residuals.  If there are less than
+	# 5 points return.
+
+	j = 0
+	nreject = 0
+	sigma = 0.
+
+	do i = 1, npts {
+	    if ((w[i] != 0.) && (rejpts[i] == NO)) {
+		sigma = sigma + Memd[residuals+i-1] ** 2
+		j = j + 1
+	    } else if (rejpts[i] == YES)
+		nreject = nreject + 1
+	}
+
+	if (j < 5) {
+	    call sfree (sp)
+	    return
+	} else
+	    sigma = sqrt (sigma / j)
+
+	if (low_reject > 0.)
+	    low_cut = -low_reject * sigma
+	else
+	    low_cut = -MAX_REAL
+	if (high_reject > 0.)
+	    high_cut = high_reject * sigma
+	else
+	    high_cut = MAX_REAL
+
+	# Reject the residuals exceeding the rejection limits.
+	# A for loop is used instead of do because with region growing we
+	# want to modify the loop index.
+
+	newreject = 0
+	for (i = 1; i <= npts; i = i + 1) {
+	    if ((w[i] == 0.) || (rejpts[i] == YES))
+		next
+
+	    residual = Memd[residuals + i - 1]
+	    if ((residual > high_cut) || (residual < low_cut)) {
+		i_min = max (1, int (i - grow))
+		i_max = min (npts, int (i + grow))
+
+		# Reject points from the fit and flag them.
+		do j = i_min, i_max {
+		    if ((abs (x[i] - x[j]) <= grow) && (w[j] != 0.) &&
+		        (rejpts[j] == NO)) {
+			if (refit == YES)
+	        	    call dcvrject (cv, x[j], y[j], w[j])
+			rejpts[j] = YES
+		        newreject = newreject + 1
+			ilast = j
+		    }
+		}
+		i = ilast
+	    }
+	}
+
+	if ((refit == YES) && (newreject > 0)) {
+	    call dcvsolve (cv, i)
+
+	    switch (i) {
+	    case SINGULAR:
+		call eprintf ("ic_reject: Singular solution\n")
+	    case NO_DEG_FREEDOM:
+		call eprintf ("ic_reject: No degrees of freedom\n")
+	    }
+	}
+
+	nreject = nreject + newreject
+
+	call sfree (sp)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicdosetup.x b/noao/imred/dtoi/hdicfit/hdicdosetup.x
new file mode 100644
index 00000000..c9e117ec
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicdosetup.x
@@ -0,0 +1,104 @@
+include	<math/curfit.h>
+include	"hdicfit.h"
+
+# IC_DOSETUP -- Setup the fit.  This is called at the start of each call
+# to ic_fit to update the fitting parameters if necessary.
+
+procedure ic_dosetupd (ic, cv, x, wts, npts, newx, newwts, newfunction, refit)
+
+pointer	ic				# ICFIT pointer
+pointer	cv				# Curfit pointer
+double	x[ARB]				# Ordinates of data
+double	wts[ARB]			# Weights
+int	npts				# Number of points in data
+int	newx				# New x points?
+int	newwts				# New weights?
+int	newfunction			# New function?
+int	refit				# Use cvrefit?
+
+int	ord
+pointer	rg_xrangesd()
+extern	hd_powerd()
+errchk	rg_xrangesd, malloc
+
+begin
+	# Set sample points.
+	if ((newx == YES) || (newwts == YES)) {
+	    if (npts == 0)
+		call error (0, "No data points for fit")
+
+	    call mfree (IC_XFIT(ic), TY_DOUBLE)
+	    call mfree (IC_YFIT(ic), TY_DOUBLE)
+	    call malloc (IC_XFIT(ic), npts, TY_DOUBLE)
+
+	    call mfree (IC_WTSFIT(ic), TY_DOUBLE)
+	    call malloc (IC_WTSFIT(ic), npts, TY_DOUBLE)
+
+	    call mfree (IC_REJPTS(ic), TY_INT)
+	    call malloc (IC_REJPTS(ic), npts, TY_INT)
+	    call amovki (NO, Memi[IC_REJPTS(ic)], npts)
+	    IC_NREJECT(ic) = 0
+
+	    # Set sample points.
+
+	    call rg_free (IC_RG(ic))
+	    IC_RG(ic) = rg_xrangesd (Memc[IC_SAMPLE(ic)], x, npts)
+	    call rg_order (IC_RG(ic))
+	    call rg_merge (IC_RG(ic))
+	    call rg_wtbind (IC_RG(ic), abs (IC_NAVERAGE(ic)), x, wts, npts,
+		Memd[IC_XFIT(ic)], Memd[IC_WTSFIT(ic)], IC_NFIT(ic))
+
+	    if (IC_NFIT(ic) == 0)
+		call error (0, "No sample points for fit")
+
+	    if (IC_NFIT(ic) == npts) {
+	        call rg_free (IC_RG(ic))
+	        call mfree (IC_XFIT(ic), TY_DOUBLE)
+	        call mfree (IC_WTSFIT(ic), TY_DOUBLE)
+	        IC_YFIT(ic) = NULL
+		IC_WTSFIT(ic) = NULL
+	    } else
+	        call malloc (IC_YFIT(ic), IC_NFIT(ic), TY_DOUBLE)
+
+	    refit = NO
+	}
+
+	# Set curve fitting parameters.
+
+	if ((newx == YES) || (newfunction == YES)) {
+	    if (cv != NULL)
+	        call dcvfree (cv)
+
+	    switch (IC_FUNCTION(ic)) {
+	    case LEGENDRE, CHEBYSHEV:
+		ord = min (IC_ORDER(ic), IC_NFIT(ic))
+	        call dcvinit (cv, IC_FUNCTION(ic), ord, double (IC_XMIN(ic)),
+		    double (IC_XMAX(ic)))
+	    case SPLINE1:
+		ord = min (IC_ORDER(ic), IC_NFIT(ic) - 1)
+		if (ord > 0)
+	            call dcvinit (cv, SPLINE1, ord, double (IC_XMIN(ic)),
+			double (IC_XMAX(ic)))
+		else
+	            call dcvinit (cv, LEGENDRE, IC_NFIT(ic),
+			double (IC_XMIN(ic)), double (IC_XMAX(ic)))
+	    case SPLINE3:
+		ord = min (IC_ORDER(ic), IC_NFIT(ic) - 3)
+		if (ord > 0)
+	            call dcvinit (cv, SPLINE3, ord, double (IC_XMIN(ic)),
+			double (IC_XMAX(ic)))
+		else
+	            call dcvinit (cv, LEGENDRE, IC_NFIT(ic),
+			double (IC_XMIN(ic)), double (IC_XMAX(ic)))
+ 	    case USERFNC:
+ 		ord = min (IC_ORDER(ic), IC_NFIT(ic))
+ 		call dcvinit (cv, USERFNC, ord, double (IC_XMIN(ic)),
+ 		    double (IC_XMAX(ic)))
+ 		call dcvuserfnc (cv, hd_powerd)
+	    default:
+		call error (0, "Unknown fitting function")
+	    }
+
+	    refit = NO
+	}
+end
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
diff --git a/noao/imred/dtoi/hdicfit/hdicerrors.x b/noao/imred/dtoi/hdicfit/hdicerrors.x
new file mode 100644
index 00000000..7722e5e2
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicerrors.x
@@ -0,0 +1,143 @@
+include	<math/curfit.h>
+include	"hdicfit.h"
+
+# IC_ERRORS -- Compute and error diagnositic information.
+
+procedure ic_errorsd (ic, file, cv, x, y, wts, npts)
+
+pointer	ic		# ICFIT pointer
+char	file[ARB]	# Output file
+pointer	cv		# Curfit pointer
+double	x[ARB]		# Ordinates
+double	y[ARB]		# Abscissas
+double	wts[ARB]	# Weights
+int	npts		# Number of data points
+
+int	i, n, deleted, ncoeffs, fd
+double	chisqr, rms
+pointer	sp, fit, wts1, coeffs, errors
+
+int	dcvstati(), open()
+double	ic_rmsd()
+errchk	open()
+
+begin
+	# Open the output file.
+	fd = open (file, APPEND, TEXT_FILE)
+
+	# Determine the number of coefficients and allocate memory.
+	ncoeffs = dcvstati (cv, CVNCOEFF)
+	call smark (sp)
+	call salloc (coeffs, ncoeffs, TY_DOUBLE)
+	call salloc (errors, ncoeffs, TY_DOUBLE)
+
+	if (npts == IC_NFIT(ic)) {
+	    # Allocate memory for the fit.
+	    n = npts
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (wts, Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+		do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Get the coefficients and compute the errors.
+	    call dcvvector (cv, x, Memd[fit], n)
+	    call dcvcoeff (cv, Memd[coeffs], ncoeffs)
+	    call dcverrors (cv, y, Memd[wts1], Memd[fit], n, chisqr,
+		Memd[errors])
+	    rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n)
+
+	} else {
+	    # Allocate memory for the fit.
+	    n = IC_NFIT(ic)
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Get the coefficients and compute the errors.
+	    call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n)
+	    rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)],
+		Memd[fit], Memd[wts1], n)
+	    call dcvcoeff (cv, Memd[coeffs], ncoeffs)
+	    call dcverrors (cv, Memd[IC_YFIT(ic)], Memd[wts1], Memd[fit],
+		n, chisqr, Memd[errors])
+	}
+
+	# Print the error analysis.
+	call fprintf (fd, "total points = %d\nsample points = %d\n")
+	    call pargi (npts)
+	    call pargi (n)
+	call fprintf (fd, "nrejected = %d\ndeleted = %d\n")
+	    call pargi (IC_NREJECT(ic))
+	    call pargi (deleted)
+	call fprintf (fd, "RMS = %7.4g\n")
+	    call pargd (rms)
+	call fprintf (fd, "square root of reduced chi square = %7.4g\n")
+	    call pargd (sqrt (chisqr))
+
+	#call fprintf (fd, "\tcoefficent\terror\n")
+	#do i = 1, ncoeffs {
+	#    call fprintf (fd, "\t%10.4e\t%10.4e\n")
+	#	call parg$t (Mem$t[coeffs+i-1])
+	#	call parg$t (Mem$t[errors+i-1])
+	#}
+
+	# Free allocated memory.
+	call sfree (sp)
+	call close (fd)
+end
+
+
+# IC_RMS -- Compute RMS of points which have not been deleted.
+
+double procedure ic_rmsd (x, y, fit, wts, npts)
+
+double	x[ARB]		# Ordinates
+double	y[ARB]		# Abscissas
+double	fit[ARB]	# Fit
+double	wts[ARB]	# Weights
+int	npts		# Number of data points
+
+int	i, n
+double	resid, rms
+
+begin
+	rms = 0.
+	n = 0
+	do i = 1, npts {
+	    if (wts[i] == 0.)
+		next
+	    resid = y[i] - fit[i]
+	    rms = rms + resid * resid
+	    n = n + 1
+	}
+
+	if (n > 0)
+	    rms = sqrt (rms / n)
+
+	return (rms)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicfit.h b/noao/imred/dtoi/hdicfit/hdicfit.h
new file mode 100644
index 00000000..89d9c73d
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicfit.h
@@ -0,0 +1,65 @@
+# Definition file for DTOI task hdfit which uses the hdicfit subdirectory.
+
+define	NSPOTS		64	   # Initially, max number of calibration spots
+define	NVALS_FIT	200	   # Number of points in vector of fitted points
+define	WT_NONE		0	   # No weighting used in fit
+define	WT_USER		1	   # User specifies weighting in fit
+define	WT_CALC		2	   # Weighting in fit calculated from std dev
+define	MIN_DEN		EPSILONR   # Density used for setting curfit minval
+define	HD_NONE		1	   # Ind var is density - no transformation
+define	HD_LOGO		2	   # Ind var is log opacitance
+define	HD_K50		3	   # Ind var is Kaiser transform w/ alpha=0.50
+define	HD_K75		4	   # Ind var is Kaiser transform w/ alpha=0.75
+define	UDELETE		100	   # Point deleted by user flag
+define	PDELETE		101	   # Point deleted by program
+define	NDELETE		102	   # Point not deleted 
+define	ADDED_PT	0.0	   # Indication of added point in sdev array
+
+# The ICFIT data structure - modified for use with the DTOI package.
+
+define	IC_NGKEYS	5		# Number of graph keys
+define	IC_LENSTRUCT	47		# Length of ICFIT structure
+
+# User fitting parameters
+define	IC_FUNCTION	Memi[$1]	# Function type
+define	IC_ORDER	Memi[$1+1]	# Order of function
+define	IC_SAMPLE	Memi[$1+2]	# Pointer to sample string
+define	IC_NAVERAGE	Memi[$1+3]	# Sampling averaging bin
+define	IC_NITERATE	Memi[$1+4]	# Number of rejection interation
+define	IC_TRANSFORM	Memi[$1+5]	# Type of transformation ** DTOI ONLY **
+define	IC_XMIN		Memr[P2R($1+6)]	# Minimum value for curve
+define	IC_XMAX		Memr[P2R($1+7)]	# Maximum value for curve
+define	IC_LOW		Memr[P2R($1+8)]	# Low rejection value
+define	IC_HIGH		Memr[P2R($1+9)]	# Low rejection value
+define	IC_GROW		Memr[P2R($1+10)]# Rejection growing radius
+
+# ICFIT parameters used for fitting
+define	IC_NFIT		Memi[$1+11]	# Number of fit points
+define	IC_NREJECT	Memi[$1+12]	# Number of rejected points
+define	IC_RG		Memi[$1+13]	# Pointer for ranges
+define	IC_XFIT		Memi[$1+14]	# Pointer to ordinates of fit points
+define	IC_YFIT		Memi[$1+15]	# Pointer to abscissas of fit points
+define	IC_WTSFIT	Memi[$1+16]	# Pointer to weights of fit points
+define	IC_REJPTS	Memi[$1+17]	# Pointer to rejected points
+
+# ICFIT parameters used for interactive graphics
+define	IC_NEWX		Memi[$1+18]	# New x fit points?
+define	IC_NEWY		Memi[$1+19]	# New y points?
+define	IC_NEWWTS	Memi[$1+20]	# New weights?
+define	IC_NEWFUNCTION	Memi[$1+21]	# New fitting function?
+define	IC_NEWTRANSFORM Memi[$1+22]	# New transform?  ** DTOI ONLY **
+define	IC_OVERPLOT	Memi[$1+23]	# Overplot next plot?
+define	IC_FITERROR	Memi[$1+24]	# Error in fit
+define	IC_LABELS	Memi[$1+25+$2-1]# Graph axis labels
+define	IC_UNITS	Memi[$1+27+$2-1]# Graph axis units
+
+define	IC_FOG		Memr[P2R($1+29)]# *** DTOI ONLY *** value of fog level
+define	IC_NEWFOG	Memi[$1+30]	# Flag for change in fog
+define	IC_RESET	Memi[$1+31]	# Flag for resetting variables
+define	IC_UPDATE	Memi[$1+32]
+define	IC_EBARS	Memi[$1+33]	# Flag for plotting error bars
+define	IC_RFOG		Memr[P2R($1+34)]# Reference value of fog for resetting
+
+# ICFIT key definitions
+define	IC_GKEY		Memi[$1+35]			# Graph key
+define	IC_AXES		Memi[$1+36+($2-1)*2+$3-1]	# Graph axis codes
diff --git a/noao/imred/dtoi/hdicfit/hdicfit.x b/noao/imred/dtoi/hdicfit/hdicfit.x
new file mode 100644
index 00000000..91e05f9d
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicfit.x
@@ -0,0 +1,80 @@
+include	<math/curfit.h>
+include	"hdicfit.h"
+
+# IC_FIT -- Fit a function.  This is the main fitting task.  It uses
+# flags to define changes since the last fit.  This allows the most
+# efficient use of the curfit and ranges packages.
+
+procedure ic_fitd (ic, cv, x, y, wts, npts, newx, newy, newwts, newfunction)
+
+pointer	ic				# ICFIT pointer
+pointer	cv				# Curfit pointer
+double	x[npts]				# Ordinates
+double	y[npts]				# Data to be fit
+double	wts[npts]			# Weights
+int	npts				# Number of points
+int	newx				# New x points?
+int	newy				# New y points?
+int	newwts				# New weights?
+int	newfunction			# New function?
+
+int	ier, refit
+errchk	ic_dosetupd
+
+begin
+	# Setup the new parameters.
+	call ic_dosetupd (ic, cv, x, wts, npts, newx, newwts, newfunction,
+	    refit)
+
+	if (npts == IC_NFIT(ic)) {
+	    # If not sampling use the data array directly.
+	    if (refit == NO) {
+	        call dcvfit (cv, x, y, wts, npts, WTS_USER, ier)
+	    } else if (newy == YES)
+	        call dcvrefit (cv, x, y, wts, ier)
+
+	} else {
+	    # If sampling first form the sample y values.
+	    if ((newx == YES) || (newy == YES) || (newwts == YES))
+	        call rg_wtbind (IC_RG(ic), IC_NAVERAGE(ic), y, wts, npts,
+		    Memd[IC_YFIT(ic)], Memd[IC_WTSFIT(ic)], IC_NFIT(ic))
+	    if (refit == NO) {
+	        call dcvfit (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)],
+		    Memd[IC_WTSFIT(ic)], IC_NFIT(ic), WTS_USER, ier)
+	    } else if (newy == YES)
+	        call dcvrefit (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)],
+		    Memd[IC_WTSFIT(ic)], ier)
+	}
+
+	# Check for an error in the fit.
+	switch (ier) {
+	case SINGULAR:
+	    call printf ("Singular solution")
+	    call flush (STDOUT)
+	case NO_DEG_FREEDOM:
+	    call printf ("No degrees of freedom")
+	    call flush (STDOUT)
+	    return
+	}
+
+	refit = YES
+
+	# Do pixel rejection if desired.
+	if ((IC_LOW(ic) > 0.) || (IC_HIGH(ic) > 0.)) {
+	    if (npts == IC_NFIT(ic)) {
+	        call ic_rejectd (cv, x, y, wts, Memi[IC_REJPTS(ic)],
+		    IC_NFIT(ic), IC_LOW(ic), IC_HIGH(ic), IC_NITERATE(ic),
+		    IC_GROW(ic), IC_NREJECT(ic))
+	    } else {
+	        call ic_rejectd (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)],
+		    Memd[IC_WTSFIT(ic)], Memi[IC_REJPTS(ic)], IC_NFIT(ic),
+		    IC_LOW(ic), IC_HIGH(ic), IC_NITERATE(ic), IC_GROW(ic),
+		    IC_NREJECT(ic))
+	    }
+
+	    if (IC_NREJECT(ic) > 0)
+		refit = NO
+	} else
+	    IC_NREJECT(ic) = 0
+end
+
diff --git a/noao/imred/dtoi/hdicfit/hdicgaxes.x b/noao/imred/dtoi/hdicfit/hdicgaxes.x
new file mode 100644
index 00000000..4f016c9d
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgaxes.x
@@ -0,0 +1,101 @@
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+# ICG_AXES -- Set axes data.
+# The applications program may set additional axes types.
+
+procedure icg_axesd (ic, gt, cv, axis, x, y, z, npts)
+
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+pointer	cv				# CURFIT pointer
+int	axis				# Output axis
+double	x[npts]				# Independent variable
+double	y[npts]				# Dependent variable
+double	z[npts]				# Output values
+int	npts				# Number of points
+
+int	i, axistype, gtlabel[2], gtunits[2]
+double 	a, b, xmin, xmax
+pointer	label, units
+
+double	dcveval(), icg_dvzd()
+errchk	adivd()
+extern	icg_dvzd()
+
+data	gtlabel/GTXLABEL, GTYLABEL/
+data	gtunits/GTXUNITS, GTYUNITS/
+
+begin
+	axistype = IC_AXES(ic, IC_GKEY(ic), axis)
+	switch (axistype) {
+	case 'x':	# Independent variable
+	    call gt_sets (gt, gtlabel[axis], Memc[IC_LABELS(ic,1)])
+	    call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,1)])
+	    call amovd (x, z, npts)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	case 'y':	# Dependent variable
+	    call gt_sets (gt, gtlabel[axis], Memc[IC_LABELS(ic,2)])
+	    call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)])
+	    call amovd (y, z, npts)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	case 'f':	# Fitted values
+	    call gt_sets (gt, gtlabel[axis], "fit")
+	    call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)])
+	    call dcvvector (cv, x, z, npts)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	case 'r':	# Residuals
+	    call gt_sets (gt, gtlabel[axis], "residuals")
+	    call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)])
+	    call dcvvector (cv, x, z, npts)
+	    call asubd (y, z, z, npts)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	case 'd':	# Ratio
+	    call gt_sets (gt, gtlabel[axis], "ratio")
+	    call gt_sets (gt, gtunits[axis], "")
+	    call dcvvector (cv, x, z, npts)
+#	    iferr (call adiv$t (y, z, z, npts))
+		call advzd (y, z, z, npts, icg_dvzd)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	case 'n':	# Linear component removed
+	    call gt_sets (gt, gtlabel[axis], "non-linear component")
+	    call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)])
+	    xmin = IC_XMIN(ic)
+	    xmax = IC_XMAX(ic)
+	    a = dcveval (cv, double(xmin))
+	    b = (dcveval (cv, double(xmax)) - a) / (xmax - xmin)
+	    do i = 1, npts
+	        z[i] = y[i] - a - b * (x[i] - xmin)
+	    call gt_sets (gt, GTSUBTITLE, "")
+	default:	# User axes types.
+	    call malloc (label, SZ_LINE, TY_CHAR)
+	    call malloc (units, SZ_LINE, TY_CHAR)
+	    if (axis == 1) {
+		call strcpy (Memc[IC_LABELS(ic,1)], Memc[label], SZ_LINE)
+		call strcpy (Memc[IC_UNITS(ic,1)], Memc[units], SZ_LINE)
+	        call amovd (x, z, npts)
+	    } else {
+		call strcpy (Memc[IC_LABELS(ic,2)], Memc[label], SZ_LINE)
+		call strcpy (Memc[IC_UNITS(ic,2)], Memc[units], SZ_LINE)
+	        call amovd (y, z, npts)
+	    }
+	    call icg_uaxesd (ic, axistype, cv, x, y, z, npts, Memc[label], 
+		Memc[units], SZ_LINE)
+	    call gt_sets (gt, gtlabel[axis], Memc[label])
+	    call gt_sets (gt, gtunits[axis], Memc[units])
+	    call gt_sets (gt, GTSUBTITLE, "HD Curve")
+	    call mfree (label, TY_CHAR)
+	    call mfree (units, TY_CHAR)
+	}
+end
+
+
+# ICG_DVZ -- Error action to take on zero division.
+
+double procedure icg_dvzd (x)
+
+double	x			# Numerator
+
+begin
+	return (1.)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgcolon.x b/noao/imred/dtoi/hdicfit/hdicgcolon.x
new file mode 100644
index 00000000..52299df7
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgcolon.x
@@ -0,0 +1,284 @@
+include	<error.h>
+include	<gset.h>
+include	"hdicfit.h"
+
+define	EB_WTS	10
+define	EB_SDEV	11
+
+# List of colon commands.
+define	CMDS "|show|sample|naverage|function|order|low_reject|high_reject|\
+              |niterate|grow|errors|vshow|transform|fog|reset|quit|ebars|"
+
+define	SHOW		1	# Show values of parameters
+define	SAMPLE		2	# Set or show sample ranges
+define	NAVERAGE	3	# Set or show sample averaging or medianing
+define	FUNCTION	4	# Set or show function type
+define	ORDER		5	# Set or show order
+define	LOW_REJECT	6	# Set or show lower rejection factor
+define	HIGH_REJECT	7	# Set or show upper rejection factor
+# newline		8
+define	NITERATE	9	# Set or show rejection iterations
+define	GROW		10	# Set or show rejection growing radius
+define	ERRORS		11	# Show errors of fit
+define	VSHOW		12	# Show verbose information
+define	TRANSFORM	13	# Set or show transformation
+define	FOG		14	# Set or show value of fog
+define	RESET		15	# Reset x, y, wts, npts to original values
+define	QUIT		16	# Terminate without updating database
+define	EBARS		17	# Set error bars to represent weights or
+				#     standard deviations
+
+# ICG_COLON -- Processes colon commands.  
+
+procedure icg_colond (ic, cmdstr, gp, gt, cv, x, y, wts, npts)
+
+pointer	ic				# ICFIT pointer
+char	cmdstr[ARB]			# Command string
+pointer	gp				# GIO pointer
+pointer	gt				# GTOOLS pointer
+pointer	cv				# CURFIT pointer for error listing
+double	x[npts], y[npts], wts[npts]	# Data arrays for error listing
+int	npts				# Number of data points
+
+real	rval
+char	cmd[SZ_LINE]
+int	ncmd, ival, ip, junk
+
+int	nscan(), strdic(), strncmp(), ctor()
+string	funcs "|chebyshev|legendre|spline1|spline3|power|"
+string	tform "|none|logopacitance|k50|k75|"
+
+begin
+	# Use formated scan to parse the command string.
+	# The first word is the command and it may be minimum match
+	# abbreviated with the list of commands.
+
+	call sscan (cmdstr)
+	call gargwrd (cmd, SZ_LINE)
+	ncmd = strdic (cmd, cmd, SZ_LINE, CMDS)
+
+	switch (ncmd) {
+	case SHOW:
+	    # show: Show the values of the fitting parameters.  The terminal
+	    # is cleared and paged using the gtools paging procedures.
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (nscan() == 1) {
+	        call gdeactivate (gp, AW_CLEAR)
+		call ic_show (ic, "STDOUT", gt)
+	        call greactivate (gp, AW_PAUSE)
+	    } else {
+		iferr (call ic_show (ic, cmd, gt))
+		    call erract (EA_WARN)
+	    }
+
+	case SAMPLE:
+	    # sample: List or set the sample points.
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+	        call printf ("sample = %s\n")
+		    call pargstr (Memc[IC_SAMPLE(ic)])
+	    } else {
+		call strcpy (cmd, Memc[IC_SAMPLE(ic)], SZ_LINE)
+		IC_NEWX(ic) = YES
+	    }
+
+	case NAVERAGE:
+	    # naverage: List or set the sample averging.
+
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("naverage = %d\n")
+		    call pargi (IC_NAVERAGE(ic))
+	    } else {
+		IC_NAVERAGE(ic) = ival
+		IC_NEWX(ic) = YES
+	    }
+
+	case FUNCTION:
+	    # function: List or set the fitting function.
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+		call printf ("function = %s\n")
+		    call ic_gstr (ic, "function", cmd, SZ_LINE)
+		    call pargstr (cmd)
+	    } else {
+		if (strdic (cmd, cmd, SZ_LINE, funcs) > 0) {
+		    call ic_pstr (ic, "function", cmd)
+		    IC_NEWFUNCTION(ic) = YES
+		} else
+		    call printf ("Unknown or ambiguous function")
+	    }
+
+	case ORDER:
+	    # order: List or set the function order.
+
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("order = %d\n")
+		    call pargi (IC_ORDER(ic))
+	    } else {
+		IC_ORDER(ic) = ival
+		IC_NEWFUNCTION(ic) = YES
+	    }
+
+	case LOW_REJECT:
+	    # low_reject: List or set the lower rejection threshold.
+
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("low_reject = %g\n")
+		    call pargr (IC_LOW(ic))
+	    } else
+		IC_LOW(ic) = rval
+
+	case HIGH_REJECT:
+	    # high_reject: List or set the high rejection threshold.
+
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("high_reject = %g\n")
+		    call pargr (IC_HIGH(ic))
+	    } else
+		IC_HIGH(ic) = rval
+
+	case NITERATE:
+	    # niterate: List or set the number of rejection iterations.
+
+	    call gargi (ival)
+	    if (nscan() == 1) {
+		call printf ("niterate = %d\n")
+		    call pargi (IC_NITERATE(ic))
+	    } else
+		IC_NITERATE(ic) = ival
+
+	case GROW:
+	    # grow: List or set the rejection growing.
+
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("grow = %g\n")
+		    call pargr (IC_GROW(ic))
+	    } else
+		IC_GROW(ic) = rval
+
+	case ERRORS:
+	    call gargwrd (cmd, SZ_LINE)
+	    if (nscan() == 1) {
+	        call gdeactivate (gp, AW_CLEAR)
+		call ic_show (ic, "STDOUT", gt)
+		call ic_errorsd (ic, "STDOUT", cv, x, y, wts, npts)
+	        call greactivate (gp, AW_PAUSE)
+	    } else {
+		iferr {
+		    call ic_show (ic, cmd, gt)
+		    call ic_errorsd (ic, cmd, cv, x, y, wts, npts)
+		} then
+		    call erract (EA_WARN)
+	    }
+	case VSHOW:
+	    # verbose show:  Show the values of the fitting parameters. 
+	    # The terminal is paged using the gtools paging procedure. 
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+		call gdeactivate (gp, AW_CLEAR)
+		call ic_vshowd (ic, "STDOUT", cv, x, y, wts, npts, gt)
+		call greactivate (gp, AW_PAUSE)
+	    } else {
+		iferr {
+		    call ic_vshowd (ic, cmd, cv, x, y, wts, npts, gt)
+		} then 
+		    call erract (EA_WARN)
+	    }
+	case TRANSFORM:
+	    # transform: List or set the transformation type.  This
+	    # option applies to HDTOI procedures only.
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+		call printf ("transform = %s\n")
+		    call ic_gstr (ic, "transform", cmd, SZ_LINE)
+		    call pargstr (cmd)
+	    } else {
+		ival= strdic (cmd, cmd, SZ_LINE, tform) 
+		if (ival > 0) {
+		    call ic_pstr (ic, "transform", cmd)
+		    IC_NEWTRANSFORM(ic) = YES
+		    IC_NEWX(ic) = YES
+		    switch (IC_TRANSFORM(ic)) {
+		    case HD_NONE:
+	    		call ic_pstr (ic, "xlabel", "Density")
+		    case HD_LOGO:
+	    		call ic_pstr (ic, "xlabel", 
+			    "Log Opacitance: log (10**Den - 1)")
+		    case HD_K50:
+	    		call ic_pstr (ic, "xlabel", 
+		           "Den + 0.50 * Log (1 - (10 ** -Den))")
+		    case HD_K75:
+	    		call ic_pstr (ic, "xlabel", 
+			    "Den + 0.75 * Log (1 - (10 ** -Den))")
+		    }
+		} else
+		    call printf ("Unknown or ambiguous transform")
+	    }
+
+	case FOG:
+	    # fog: DTOI ONLY - change or reset the value of the fog level
+
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+		call printf ("fog = %g\n")
+		    call pargr (IC_FOG(ic))
+	    } else {
+		if (strncmp (cmd, "reset", 1) == 0)
+		    IC_FOG(ic) = IC_RFOG(ic)
+		else {
+		    ip = 1
+		    junk = ctor (cmd, ip, rval)
+		    IC_FOG(ic) = rval
+		}
+		IC_NEWFOG(ic) = YES
+		IC_NEWX(ic) = YES
+	    }
+
+	case RESET:
+	    # Set flag to reset x, y, wts and npts to original values.
+	    IC_RESET(ic) = YES
+	    IC_NEWX(ic) = YES
+	    IC_NEWY(ic) = YES
+	    IC_NEWWTS(ic) = YES
+	    IC_NEWFUNCTION(ic) = YES
+	    IC_NEWTRANSFORM(ic) = YES
+
+	case QUIT:
+	    # Set update flag to know
+	    IC_UPDATE(ic) = NO
+
+	case EBARS:
+	    # [HV]BAR marker can indicate either errors or weights
+	    call gargwrd (cmd, SZ_LINE)
+	    if (cmd[1] == EOS) {
+		if (IC_EBARS(ic) == EB_WTS)
+		    call printf ("ebars = Weights\n")
+		else if (IC_EBARS(ic) == EB_SDEV)
+		    call printf ("ebars = Errors\n")
+	    } else {
+		if (strncmp (cmd, "weights", 1) == 0 || 
+		    strncmp (cmd, "WEIGHTS", 1) == 0)
+		    IC_EBARS(ic) = EB_WTS
+		else if (strncmp (cmd, "errors", 1) == 0 || 
+		    strncmp (cmd, "ERRORS", 1) == 0)
+		    IC_EBARS(ic) = EB_SDEV
+		else
+		    call printf ("Unrecognized value for ebars '%s'\n")
+			call pargstr (cmd)
+	    }
+
+	default:
+	    call eprintf ("Unrecognized command '%s'\n")
+		call pargstr (cmd)
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgdelete.x b/noao/imred/dtoi/hdicfit/hdicgdelete.x
new file mode 100644
index 00000000..82c61f70
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgdelete.x
@@ -0,0 +1,81 @@
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+define	MSIZE		2.		# Mark size
+
+# ICG_DELETE -- Delete data point nearest the cursor.
+# The nearest point to the cursor in NDC coordinates is determined.
+
+procedure icg_deleted (ic, gp, gt, cv, x, y, wts, npts, wx, wy)
+
+pointer	ic					# ICFIT pointer
+pointer	gp					# GIO pointer
+pointer	gt					# GTOOLS pointer
+pointer	cv					# CURFIT pointer
+double	x[npts], y[npts]			# Data points
+double	wts[npts]				# Weight array
+int	npts					# Number of points
+real	wx, wy					# Position to be nearest
+
+int	gt_geti()
+pointer	sp, xout, yout
+
+begin
+	call smark (sp)
+	call salloc (xout, npts, TY_DOUBLE)
+	call salloc (yout, npts, TY_DOUBLE)
+
+	call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts)
+	call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts)
+	if (gt_geti (gt, GTTRANSPOSE) == NO)
+	    call icg_d1d (ic, gp, Memd[xout], Memd[yout], wts, npts, wx, wy)
+	else
+	    call icg_d1d (ic, gp, Memd[yout], Memd[xout], wts, npts, wy, wx)
+
+	call sfree (sp)
+end
+
+
+# ICG_D1D - Do the actual delete.
+
+procedure icg_d1d (ic, gp, x, y, wts, npts, wx, wy)
+
+pointer	ic					# ICFIT pointer
+pointer	gp					# GIO pointer
+double	x[npts], y[npts]			# Data points
+double	wts[npts]				# Weight array
+int	npts					# Number of points
+real	wx, wy					# Position to be nearest
+
+int	i, j
+real	x0, y0, r2, r2min
+
+begin
+	# Transform world cursor coordinates to NDC.
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+
+	# Search for nearest point to a point with non-zero weight.
+	r2min = MAX_REAL
+	do i = 1, npts {
+	    if (wts[i] == 0.)
+		next
+
+	    call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0)
+		
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		j = i
+	    }
+	}
+
+	# Mark the deleted point with a cross and set the weight to zero.
+	if (j != 0) {
+	    call gscur (gp, real (x[j]), real (y[j]))
+	    call gmark (gp, real (x[j]), real (y[j]), GM_CROSS, MSIZE, MSIZE)
+	    wts[j] = 0.
+	    IC_NEWWTS(ic) = YES
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgfit.x b/noao/imred/dtoi/hdicfit/hdicgfit.x
new file mode 100644
index 00000000..b50ed03c
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgfit.x
@@ -0,0 +1,402 @@
+include	<error.h>
+include	<pkg/gtools.h>
+include	<mach.h>
+include	<gset.h>
+include	"hdicfit.h"
+
+define	HELP		"noao$lib/scr/hdicgfit.key"
+define	PROMPT		"hdicfit options"
+define	EB_WTS		10
+define	EB_SDEV		11
+
+# ICG_FIT -- Interactive curve fitting with graphics.  This is the main
+# entry point for the interactive graphics part of the icfit package.
+
+procedure icg_fitd (ic, gp, cursor, gt, cv, density, loge, owts, osdev, npts)
+
+pointer	ic			# ICFIT pointer
+pointer	gp			# GIO pointer
+char	cursor[ARB]		# GIO cursor input
+pointer	gt			# GTOOLS pointer
+pointer	cv			# CURFIT pointer
+double	density[npts]		# Original density, not fog subtracted
+double	loge[npts]		# Original log exposure values
+double	owts[npts]		# Original weights array
+double	osdev[npts]		# Original standard deviation array
+int	npts			# Number of points
+
+real	wx, wy, wwy, oldx, oldy
+int	wcs, key, nptso
+char	cmd[SZ_LINE]
+
+int	i, newgraph, axes[2], linetype
+double	x1, newwt
+pointer	userwts, x, y, wts, den, whydel, sdev, ebw
+
+int	clgcur(), stridxs(), scan(), nscan()
+int	icg_nearestd(), gstati()
+double	dcveval()
+errchk	ic_fitd, malloc
+define	redraw_ 91
+
+begin
+	# Allocate memory for the fit and a copy of the weights.
+	# The weights are copied because they are changed when points are
+	# deleted.  The input x is the untransformed density, and is used
+	# to generate other types of transforms.  Points can be added to
+	# the sample, so the y array and weights array can change as well.
+	# The original number of points is also remembered.
+
+	call malloc (userwts, npts, TY_DOUBLE)
+	call malloc (x,       npts, TY_DOUBLE)
+	call malloc (y,       npts, TY_DOUBLE)
+	call malloc (wts,     npts, TY_DOUBLE)
+	call malloc (den,     npts, TY_DOUBLE)
+	call malloc (whydel,  npts, TY_INT)
+	call malloc (sdev,    npts, TY_DOUBLE)
+	call malloc (ebw,     npts, TY_DOUBLE)
+
+	call amovd (owts,     Memd[userwts], npts)
+	call amovd (owts,     Memd[wts],     npts)
+	call amovd (loge,     Memd[y],       npts)
+	call amovd (density,  Memd[den],     npts)
+	call amovki (NDELETE, Memi[whydel],  npts)
+	call amovd (osdev,    Memd[sdev],    npts)
+	nptso = npts
+
+	# Initialize
+	IC_OVERPLOT(ic) = NO
+	IC_NEWX(ic) = YES
+	IC_NEWY(ic) = YES
+	IC_NEWWTS(ic) = YES
+	IC_NEWFUNCTION(ic) = YES
+	IC_NEWTRANSFORM(ic) = YES
+	IC_UPDATE(ic) = YES
+	IC_EBARS(ic) = EB_SDEV
+
+	# Read cursor commands.
+
+	key = 'f'
+	axes[1] = IC_AXES(ic, IC_GKEY(ic), 1)
+	axes[2] = IC_AXES(ic, IC_GKEY(ic), 2)
+
+	repeat {
+	    switch (key) {
+	    case '?': # Print help text.
+		call gpagefile (gp, HELP, PROMPT)
+
+	    case 'q': # Terminate cursor loop
+		break
+
+	    case ':': # List or set parameters
+		if (stridxs ("/", cmd) == 1)
+	            call gt_colon (cmd, gp, gt, newgraph)
+		else
+		    call icg_colond (ic, cmd, gp, gt, cv, Memd[x], 
+			Memd[y], Memd[wts], npts)
+
+		if (IC_RESET(ic) == YES) {
+		    npts = nptso
+		    call amovd (owts, Memd[userwts], npts)
+		    call amovd (owts, Memd[wts], npts)
+		    call amovd (loge, Memd[y], npts)
+		    call amovd (density, Memd[den], npts)
+		    call amovki (NDELETE, Memi[whydel], npts)
+		    call amovd (osdev, Memd[sdev], npts)
+		    call hdic_init (density, npts, 0.0)
+		}
+
+		# See if user wants to quit without updating
+		if (IC_UPDATE(ic) == NO) {
+		    call mfree (x, TY_DOUBLE)
+		    call mfree (y, TY_DOUBLE)
+		    call mfree (wts, TY_DOUBLE)
+		    call mfree (userwts, TY_DOUBLE)
+		    call mfree (den, TY_DOUBLE)
+		    call mfree (sdev, TY_DOUBLE)
+		    return
+		}
+
+	    case 'a': # Add data points to the sample.  This is only possible 
+		      # from an HD curve plot.
+
+		if ((IC_AXES (ic, IC_GKEY(ic), 1) == 'y') &&
+		    (IC_AXES (ic, IC_GKEY(ic), 2) == 'u')) {
+
+		    # Query for weight after plotting current location
+		    # call gt_plot (gp, gt, wx, wy, 1)
+		    call gmark (gp, wx, wy, GM_CIRCLE, 2.0, 2.0)
+		    newwt = 1.0D0
+		    call printf ("Enter weight of new point (%g): ")
+			call pargd (newwt)
+		    call flush (STDOUT)
+		    if (scan() != EOF) {
+		        call pargd (x1)
+		        if (nscan() == 1) {
+			    if (!IS_INDEFD (x1)) {
+			        newwt = x1
+			    }
+		        }
+		    }
+
+		} else {
+		    call eprintf ("Points can be added only from an HD Curve\n")
+		    next
+		}
+
+		# Add fog into "density above fog" value read from cursor
+		wwy = wy + IC_FOG (ic)
+		if (wwy < 0.0) {
+		    call eprintf (
+			"New density (%g) is below fog and will not be added\n")
+			    call pargr (wwy)
+		    next
+		}
+
+		# Add point into sample
+		call eprintf ("New Point: density above fog = %.4f, log ")
+		    call pargr (wwy)
+		call eprintf ("exposure = %.4f, weight = %.4f\n")
+		    call pargr (wx)
+		    call pargd (newwt)
+
+		call hdic_addpoint (ic, wwy, wx, newwt, den, y, wts, userwts, 
+		    x, whydel, sdev, npts)
+
+		call realloc (ebw, npts, TY_DOUBLE)
+
+		call hdic_transform (ic, Memd[den], Memd[userwts], Memd[x], 
+		    Memd[wts], Memi[whydel], npts)
+
+	    case 'c': # Print the positions of data points.
+		i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, 
+		    wx, wy)
+
+	    	if (i != 0) {
+		    call printf ("den= %7.4g x= %7.4g  exp= %7.4g  fit= %7.4g")
+			call pargd (Memd[den+i-1])
+			call pargd (Memd[x+i-1])
+			call pargd (Memd[y+i-1])
+			call pargd (dcveval (cv, Memd[x+i-1]))
+		}
+
+	    case 'd': # Delete data points.
+		call icg_deleted (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts], 
+		    npts, wx, wy)
+
+	    case 'f': # Fit the function and reset the flags.
+		iferr {
+		    # Copy new transformed vector, if necessary
+		    if (IC_NEWTRANSFORM(ic) == YES || IC_NEWFOG(ic) == YES) 
+			call hdic_transform (ic, Memd[den], Memd[userwts], 
+			    Memd[x], Memd[wts], Memi[whydel], npts)
+
+		    call ic_fitd (ic, cv, Memd[x], Memd[y], Memd[wts], npts, 
+			IC_NEWX(ic), IC_NEWY(ic), IC_NEWWTS(ic), 
+			    IC_NEWFUNCTION(ic))
+
+		    IC_NEWX(ic) = NO
+		    IC_NEWY(ic) = NO
+		    IC_NEWWTS(ic) = NO
+		    IC_NEWFUNCTION(ic) = NO
+		    IC_NEWTRANSFORM(ic) = NO
+		    IC_FITERROR(ic) = NO
+		    IC_NEWFOG(ic) = NO
+		    newgraph = YES
+		} then {
+		    IC_FITERROR(ic) = YES
+		    call erract (EA_WARN)
+		    newgraph = NO
+		}
+
+	    case 'g':	# Set graph axes types.
+		call printf ("Graph key to be defined: ")
+		call flush (STDOUT)
+		if (scan() == EOF)
+		    goto redraw_
+		call gargc (cmd[1])
+
+		switch (cmd[1]) {
+		case '\n':
+		case 'h', 'i', 'j', 'k', 'l':
+		    switch (cmd[1]) {
+		    case 'h':
+		        key = 1
+		    case 'i':
+		        key = 2
+		    case 'j':
+		        key = 3
+		    case 'k':
+		        key = 4
+		    case 'l':
+		        key = 5
+		    }
+
+		    call printf ("Set graph axes types (%c, %c): ")
+		        call pargc (IC_AXES(ic, key, 1))
+		        call pargc (IC_AXES(ic, key, 2))
+		    call flush (STDOUT)
+		    if (scan() == EOF)
+		        goto redraw_
+		    call gargc (cmd[1])
+
+		    switch (cmd[1]) {
+		    case '\n':
+		    default:
+		        call gargc (cmd[2])
+		        call gargc (cmd[2])
+		        if (cmd[2] != '\n') {
+			    IC_AXES(ic, key, 1) = cmd[1]
+			    IC_AXES(ic, key, 2) = cmd[2]
+		        }
+		    }
+		default:
+		    call printf ("Not a graph key")
+		}
+
+	    case 'h':
+		if (IC_GKEY(ic) != 1) {
+		    IC_GKEY(ic) = 1
+		    newgraph = YES
+		}
+
+	    case 'i':
+		if (IC_GKEY(ic) != 2) {
+		    IC_GKEY(ic) = 2
+		    newgraph = YES
+		}
+
+	    case 'j':
+		if (IC_GKEY(ic) != 3) {
+		    IC_GKEY(ic) = 3
+		    newgraph = YES
+		}
+
+	    case 'k':
+		if (IC_GKEY(ic) != 4) {
+		    IC_GKEY(ic) = 4
+		    newgraph = YES
+		}
+
+	    case 'l':
+		if (IC_GKEY(ic) != 5) {
+		    IC_GKEY(ic) = 5
+		    newgraph = YES
+		}
+
+	    case 'o': # Set overplot flag
+		IC_OVERPLOT(ic) = YES
+
+	    case 'r': # Redraw the graph
+		newgraph = YES
+
+	    case 'u': # Undelete data points.
+		call icg_undeleted (ic, gp, gt, cv, Memd[x], Memd[y], 
+		    Memd[wts], Memd[userwts], npts, wx, wy)
+
+	    case 'w':  # Window graph
+		call gt_window (gt, gp, cursor, newgraph)
+
+	    case 'x': # Reset the value of the x point.
+		i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, 
+		    wy)
+
+	    	if (i != 0) {
+		    call printf ("Enter new x (%g): ")
+			call pargd (Memd[x+i-1])
+		    call flush (STDOUT)
+		    if (scan() != EOF) {
+		        call gargd (x1)
+		        if (nscan() == 1) {
+			    if (!IS_INDEF (x1)) {
+				oldx = Memd[x+i-1]
+				oldy = Memd[y+i-1]
+			        Memd[x+i-1] = x1
+				call hd_redraw (gp, oldx, oldy, x1, oldy)
+			        IC_NEWX(ic) = YES
+			    }
+			}
+		    }
+		}
+
+	    case 'y': # Reset the value of the y point.
+		i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, 
+		    wy)
+
+	    	if (i != 0) {
+		    call printf ("Enter new y (%g): ")
+			call pargd (Memd[y+i-1])
+		    call flush (STDOUT)
+		    if (scan() != EOF) {
+		        call gargd (x1)
+		        if (nscan() == 1) {
+			    if (!IS_INDEF (x1)) {
+				oldx = Memd[x+i-1]
+				oldy = Memd[y+i-1]
+			        Memd[y+i-1] = x1
+				call hd_redraw (gp, oldx, oldy, oldx, x1)
+			        IC_NEWY(ic) = YES
+			    }
+			}
+		    }
+		}
+
+	    case 'z': # Reset the weight value of the nearest point
+		i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, 
+		    wy)
+
+	    	if (i != 0) {
+		    call printf ("Enter new weight (%g): ")
+			call pargd (Memd[wts+i-1])
+		    call flush (STDOUT)
+		    if (scan() != EOF) {
+		        call gargd (x1)
+		        if (nscan() == 1) {
+			    if (!IS_INDEF (x1)) {
+			        Memd[wts+i-1] = x1
+			        IC_NEWWTS(ic) = YES
+			    }
+			}
+		    }
+		}
+
+	    default: # Let the user decide on any other keys.
+		call icg_user (ic, gp, gt, cv, wx, wy, wcs, key, cmd)
+	    }
+
+	    # Redraw the graph if necessary.
+redraw_	    if (newgraph == YES) {
+		if (IC_AXES(ic, IC_GKEY(ic), 1) != axes[1]) {
+		    axes[1] = IC_AXES(ic, IC_GKEY(ic), 1)
+		    call gt_setr (gt, GTXMIN, INDEF)
+		    call gt_setr (gt, GTXMAX, INDEF)
+		}
+		if (IC_AXES(ic, IC_GKEY(ic), 2) != axes[2]) {
+		    axes[2] = IC_AXES(ic, IC_GKEY(ic), 2)
+		    call gt_setr (gt, GTYMIN, INDEF)
+		    call gt_setr (gt, GTYMAX, INDEF)
+		}
+
+		# Overplot with a different line type
+		if (IC_OVERPLOT(ic) == YES)
+		    linetype = min ((gstati (gp, G_PLTYPE) + 1), 4)
+		else
+		    linetype = GL_SOLID
+		call gseti (gp, G_PLTYPE, linetype)
+
+		call hdic_ebw (ic, Memd[den], Memd[x], Memd[sdev], Memd[ebw], 
+		    npts)
+
+	    	call icg_graphd (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts], 
+		    Memd[ebw], npts)
+
+		newgraph = NO
+	    }
+	} until (clgcur (cursor, wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+
+	call mfree (x, TY_DOUBLE)
+	call mfree (y, TY_DOUBLE)
+	call mfree (wts, TY_DOUBLE)
+	call mfree (userwts, TY_DOUBLE)
+	call mfree (den, TY_DOUBLE)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicggraph.x b/noao/imred/dtoi/hdicfit/hdicggraph.x
new file mode 100644
index 00000000..dcd9ade3
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicggraph.x
@@ -0,0 +1,329 @@
+include	<mach.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+define	MSIZE		2.		# Mark size
+define	SZ_TPLOT	7		# String length for plot type
+define	EB_WTS		10
+define	EB_SDEV		11
+define	MAX_SZMARKER	4
+
+# ICG_GRAPH -- Graph data and fit.
+
+procedure icg_graphd (ic, gp, gt, cv, x, y, wts, ebw, npts)
+
+pointer	ic				# ICFIT pointer
+pointer	gp				# GIO pointer
+pointer	gt				# GTOOLS pointers
+pointer	cv				# Curfit pointer
+double	x[npts]				# Independent variable
+double	y[npts]				# Dependent variable
+double	wts[npts]			# Weights
+double	ebw[npts]			# Half the error bar width
+int	npts				# Number of points
+
+pointer	xout, yout
+int	nvalues
+errchk	malloc
+
+begin
+	call malloc (xout, npts, TY_DOUBLE)
+	call malloc (yout, npts, TY_DOUBLE)
+
+	call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts)
+	call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts)
+	call icg_paramsd (ic, cv, x, y, wts, npts, gt)
+	
+	call icg_g1d (ic, gp, gt, Memd[xout], Memd[yout], wts, ebw, npts)
+
+	if (npts != IC_NFIT(ic)) {
+	    if ((abs (IC_NAVERAGE(ic)) > 1) || (IC_NREJECT(ic) > 0)) {
+	        call realloc (xout, IC_NFIT(ic), TY_DOUBLE)
+		call realloc (yout, IC_NFIT(ic), TY_DOUBLE)
+		call icg_axesd (ic, gt, cv, 1, Memd[IC_XFIT(ic)],
+		    Memd[IC_YFIT(ic)], Memd[xout], IC_NFIT(ic))
+		call icg_axesd (ic, gt, cv, 2, Memd[IC_XFIT(ic)],
+		    Memd[IC_YFIT(ic)], Memd[yout], IC_NFIT(ic))
+		call icg_g2d (ic, gp, gt, Memd[xout], Memd[yout],
+		    IC_NFIT(ic))
+	    }
+
+	} else if (IC_NREJECT(ic) > 0)
+	    call icg_g2d (ic, gp, gt, Memd[xout], Memd[yout], npts)
+
+	nvalues = NVALS_FIT
+	call icg_gfd (ic, gp, gt, cv, nvalues)
+
+	# Mark the the sample regions.
+	call icg_sampled (ic, gp, gt, x, npts, 1)
+
+	call mfree (xout, TY_DOUBLE)
+	call mfree (yout, TY_DOUBLE)
+end
+
+
+# ICG_G1D -- Plot data vector.
+
+procedure icg_g1d (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[npts]		# Array of independent variables
+double	y[npts]		# Array of dependent variables
+double	wts[npts]	# Array of weights
+double	ebw[npts]	# Error bar half width in WCS (positive density)
+int	npts		# Number of points
+
+pointer	sp, xr, yr, sz, szmk, gt1, gt2
+char	tplot[SZ_TPLOT]
+int	i, xaxis, yaxis, symbols[3], markplot
+real	size, rmin, rmax, big
+bool	fp_equald(), streq(), fp_equalr()
+int	strncmp()
+data	symbols/GM_PLUS, GM_BOX, GM_CIRCLE/
+include	"hdic.com"
+
+begin
+	call smark (sp)
+	call salloc (xr,   npts, TY_REAL)
+	call salloc (yr,   npts, TY_REAL)
+	call salloc (sz,   npts, TY_REAL)
+	call salloc (szmk, npts, TY_REAL)
+
+	call achtdr (x, Memr[xr], npts)
+	call achtdr (y, Memr[yr], npts)
+
+	xaxis = (IC_AXES (ic, IC_GKEY(ic), 1))
+	yaxis = (IC_AXES (ic, IC_GKEY(ic), 2))
+
+	# Set up gtools structure for deleted (gt1) and added (gt2) points
+	call gt_copy (gt, gt1)
+	call gt_setr (gt1, GTXSIZE, MSIZE)
+	call gt_setr (gt1, GTYSIZE, MSIZE)
+	call gt_sets (gt1, GTMARK,  "cross")
+
+	call gt_copy (gt, gt2)
+	call gt_setr (gt2, GTXSIZE, MSIZE)
+	call gt_setr (gt2, GTYSIZE, MSIZE)
+	call gt_sets (gt2, GTMARK,  "circle")
+
+	markplot = NO
+	call gt_gets (gt, GTTYPE, tplot, SZ_TPLOT)
+	if (strncmp (tplot, "mark", 4) == 0)
+	    markplot = YES
+
+	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 (xaxis == 'y' && yaxis == 'u')
+		call gswind (gp, INDEF, INDEF, INDEF, real (maxden))
+
+	    call gt_swind (gp, gt)
+	    call gt_labax (gp, gt)
+	}
+
+	# Calculate size of markers if error bars are being used.  If the
+	# weights are being used as the marker size, they are first scaled
+	# between 1.0 and the maximum marker size.
+
+	call gt_gets (gt, GTMARK, tplot, SZ_TPLOT)
+	if (streq (tplot, "hebar") || streq (tplot, "vebar")) {
+	    if (IC_EBARS(ic) == EB_WTS) {
+		call achtdr (wts, Memr[sz], npts)
+		call alimr  (Memr[sz], npts, rmin, rmax)
+		if (fp_equalr (rmin, rmax))
+		    call amovr (Memr[sz], Memr[szmk], npts)
+		else {
+		    big = real  (MAX_SZMARKER)
+		    call amapr  (Memr[sz], Memr[szmk], npts,rmin,rmax,1.0, big)
+		}
+	    } else {
+		call achtdr (ebw, Memr[sz], npts)
+		call hd_szmk (gp, gt, xaxis, yaxis, tplot, Memr[sz], 
+		    Memr[szmk], npts)
+	    }
+	} else
+	    call amovkr (MSIZE, Memr[szmk], npts)
+
+	do i = 1, npts {
+	    # Check for deleted point
+	    if (fp_equald (wts[i], 0.0D0))
+	        call gt_plot (gp, gt1, Memr[xr+i-1], Memr[yr+i-1], 1)
+
+	    # Check for added point
+	    else if (fp_equald (ebw[i], ADDED_PT))
+	        call gt_plot (gp, gt2, Memr[xr+i-1], Memr[yr+i-1], 1)
+
+	    else {
+	        size = Memr[szmk+i-1]
+	        call gt_setr (gt, GTXSIZE, size)
+	        call gt_setr (gt, GTYSIZE, size)
+                call gt_plot (gp, gt, Memr[xr+i-1], Memr[yr+i-1], 1)
+	    }
+	}
+	
+	IC_OVERPLOT(ic) = NO
+	call sfree (sp)
+end
+
+
+# HD_SZMK -- Calculate size of error bar markers.  This procedure is
+# called when the marker type is hebar or vebar and the marker size is
+# to be the error in the point, not its weight in the fit.
+
+procedure hd_szmk (gp, gt, xaxis, yaxis, mark, insz, outsz, npts)
+
+pointer	gp			# Pointer to gio structure
+pointer	gt			# Pointer to gtools structure
+int	xaxis, yaxis		# Codes for x and y axis types
+char	mark[SZ_TPLOT]		# Type of marker to use
+real	insz[npts]		# Standard deviations in WCS units
+real	outsz[npts]		# Output size array in NDC units
+int	npts			# Number of points in arrays
+
+int	gt_geti()
+char	tplot[SZ_TPLOT]
+real	dx, dy
+bool	streq()
+
+begin
+	# Check validity of axis types
+	if (xaxis != 'x' && xaxis != 'u' && yaxis != 'x' && yaxis != 'u') {
+	    call eprintf ("Choose graph type with axes 'u' or 'x'; ")
+	    call eprintf ("Using marker size 2.0\n")
+	    call amovkr (MSIZE, outsz, npts)
+	    return
+	}
+
+	call gt_gets (gt, GTMARK, tplot, SZ_TPLOT)
+	if (streq (tplot, "hebar")) {
+	    if (yaxis == 'x' || yaxis == 'u') {
+	        call gt_sets (gt, GTMARK, "vebar")
+	        call eprintf ("Marker switched to vebar\n")
+		# call flush (STDOUT)
+	    }
+	} else if (streq (tplot, "vebar")) {
+	    if (xaxis == 'x' || xaxis == 'u') {
+	        call gt_sets (gt, GTMARK, "hebar")
+	        call eprintf ("Marker switched to hebar\n")
+		# call flush (STDOUT)
+	    }
+	}
+
+	# Need to scale standard deviation from density to NDC units
+        call ggscale (gp, 0.0, 0.0, dx, dy)
+	if (gt_geti (gt, GTTRANSPOSE) == NO)
+	    call adivkr (insz, dx, outsz, npts)
+	else
+	    call adivkr (insz, dy, outsz, npts)
+end
+
+
+# ICG_G2D -- Show sample range and rejected data on plot.
+
+procedure icg_g2d (ic, gp, gt, x, y, npts)
+
+pointer	ic			# ICFIT pointer
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+double	x[npts], y[npts]	# Data points
+int	npts			# Number of data points
+
+int	i
+pointer	sp, xr, yr, gt1
+
+begin
+	call smark (sp)
+	call salloc (xr, npts, TY_REAL)
+	call salloc (yr, npts, TY_REAL)
+	call achtdr (x, Memr[xr], npts)
+	call achtdr (y, Memr[yr], npts)
+
+	call gt_copy (gt, gt1)
+	call gt_sets (gt1, GTTYPE, "mark")
+
+	# Mark the sample points.
+
+	if (abs (IC_NAVERAGE(ic)) > 1) {
+	    call gt_sets (gt1, GTMARK, "plus")
+	    call gt_setr (gt1, GTXSIZE, real (-abs (IC_NAVERAGE(ic))))
+	    call gt_setr (gt1, GTYSIZE, 1.)
+	    call gt_plot (gp, gt1, Memr[xr], Memr[yr], npts)
+	}
+
+	# Mark the rejected points.
+
+	if (IC_NREJECT(ic) > 0) {
+	    call gt_sets (gt1, GTMARK, "diamond")
+	    call gt_setr (gt1, GTXSIZE, MSIZE)
+	    call gt_setr (gt1, GTYSIZE, MSIZE)
+	    do i = 1, npts {
+		if (Memi[IC_REJPTS(ic)+i-1] == YES)
+	            call gt_plot (gp, gt1, Memr[xr+i-1], Memr[yr+i-1], 1)
+	    }
+	}
+
+	call gt_free (gt1)
+	call sfree (sp)
+end
+
+
+# ICG_GF[RD] -- Overplot the fit on a graph of the data points.  A vector
+# is constructed and evaluated, then plotted.
+
+procedure icg_gfd (ic, gp, gt, cv, npts)
+
+pointer	ic			# ICFIT pointer
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOL pointer
+pointer	cv			# CURFIT pointer
+int	npts			# Number of points to plot
+
+pointer	sp, xr, yr, x, y, xo, yo, gt1
+int	gstati()
+
+begin
+	call smark (sp)
+
+	if (IC_FITERROR(ic) == YES)
+	    return
+
+	call salloc (xr, npts, TY_REAL)
+	call salloc (yr, npts, TY_REAL)
+	call salloc (x,  npts, TY_DOUBLE)
+	call salloc (y,  npts, TY_DOUBLE)
+	call salloc (xo, npts, TY_DOUBLE)
+	call salloc (yo, npts, TY_DOUBLE)
+
+	# Calculate big vector of independent variable values.  Note value
+	# of npts can change with this call.
+	call hdic_gvec (ic, Memd[x], npts, IC_TRANSFORM(ic))
+
+	# Calculate vector of fit values. 
+	call dcvvector (cv, Memd[x], Memd[y], npts)
+
+	# Convert to user function or transpose axes.  Change type to reals
+	# for plotting.
+	call icg_axesd (ic, gt, cv, 1, Memd[x], Memd[y], Memd[xo], npts)
+	call icg_axesd (ic, gt, cv, 2, Memd[x], Memd[y], Memd[yo], npts)
+	call achtdr (Memd[xo], Memr[xr], npts)
+	call achtdr (Memd[yo], Memr[yr], npts)
+
+	call gt_copy (gt, gt1)
+	call gt_sets (gt1, GTTYPE, "line")
+	call gt_seti (gt1, GTLINE, gstati (gp, G_PLTYPE))
+	call gt_plot (gp, gt1, Memr[xr], Memr[yr], npts)
+	call gt_free (gt1)
+
+	call sfree (sp)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgnearest.x b/noao/imred/dtoi/hdicfit/hdicgnearest.x
new file mode 100644
index 00000000..8d556273
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgnearest.x
@@ -0,0 +1,72 @@
+include	<mach.h>
+include	<pkg/gtools.h>
+
+# ICG_NEAREST -- Find the nearest point to the cursor and return the index.
+# The nearest point to the cursor in NDC coordinates is determined.
+# The cursor is moved to the nearest point selected.
+
+int procedure icg_nearestd (ic, gp, gt, cv, x, y, npts, wx, wy)
+
+pointer	ic					# ICFIT pointer
+pointer	gp					# GIO pointer
+pointer	gt					# GTOOLS pointer
+pointer	cv					# CURFIT pointer
+double	x[npts], y[npts]			# Data points
+int	npts					# Number of points
+real	wx, wy					# Cursor position
+
+int	pt
+pointer	sp, xout, yout
+int	icg_nd(), gt_geti()
+
+begin
+	call smark (sp)
+	call salloc (xout, npts, TY_DOUBLE)
+	call salloc (yout, npts, TY_DOUBLE)
+
+	call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts)
+	call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts)
+
+	if (gt_geti (gt, GTTRANSPOSE) == NO)
+	    pt = icg_nd (gp, Memd[xout], Memd[yout], npts, wx, wy)
+	else
+	    pt = icg_nd (gp, Memd[yout], Memd[xout], npts, wy, wx)
+	call sfree (sp)
+	
+	return (pt)
+end
+
+
+# ICG_ND -- ??
+
+int procedure icg_nd (gp, x, y, npts, wx, wy)
+
+pointer	gp					# GIO pointer
+double	x[npts], y[npts]			# Data points
+int	npts					# Number of points
+real	wx, wy					# Cursor position
+
+int	i, j
+real	x0, y0, r2, r2min
+
+begin
+	# Transform world cursor coordinates to NDC.
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+
+	# Search for nearest point.
+	r2min = MAX_REAL
+	do i = 1, npts {
+	    call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		j = i
+	    }
+	}
+
+	# Move the cursor to the selected point and return the index.
+	if (j != 0)
+	    call gscur (gp, real (x[j]), real (y[j]))
+
+	return (j)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgparams.x b/noao/imred/dtoi/hdicfit/hdicgparams.x
new file mode 100644
index 00000000..e502fba1
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgparams.x
@@ -0,0 +1,94 @@
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+# ICG_PARAMS -- Set parameter string.
+
+procedure icg_paramsd (ic, cv, x, y, wts, npts, gt)
+
+pointer	ic		# ICFIT pointer
+pointer	cv		# Curfit pointer
+double	x[ARB]		# Ordinates
+double	y[ARB]		# Abscissas
+double	wts[ARB]	# Weights
+int	npts		# Number of data points
+pointer	gt		# GTOOLS pointer
+
+double	rms
+int	i, n, deleted
+pointer	sp, fit, wts1, str, params
+double	ic_rmsd()
+
+begin
+	call smark (sp)
+
+	if (npts == IC_NFIT(ic)) {
+	    # Allocate memory for the fit.
+	    n = npts
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (wts, Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+		do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Set the fit and compute the RMS error.
+	    call dcvvector (cv, x, Memd[fit], n)
+	    rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n)
+
+	} else {
+	    # Allocate memory for the fit.
+	    n = IC_NFIT(ic)
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Set the fit and compute the rms error.
+	    call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n)
+	    rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], Memd[fit],
+		Memd[wts1], n)
+	}
+
+	# Print the parameters and errors.
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (params, 2*SZ_LINE, TY_CHAR)
+
+	call sprintf (Memc[str],SZ_LINE, "function=%s, order=%d, transform=%s")
+	    call ic_gstr (ic, "function", Memc[params], 2*SZ_LINE)
+	    call pargstr (Memc[params])
+	    call pargi (IC_ORDER(ic))
+	    call ic_gstr (ic, "transform", Memc[params], SZ_LINE)
+	    call pargstr (Memc[params])
+	call sprintf (Memc[params], 2*SZ_LINE,
+            "%s\nfog=%.5f,  total=%d, deleted=%d, RMS=%7.4g")
+	    call pargstr (Memc[str])
+	    call pargr (IC_FOG(ic))
+	    call pargi (npts)
+	    call pargi (deleted)
+	    call pargd (rms)
+	call gt_sets (gt, GTPARAMS, Memc[params])
+
+	call sfree (sp)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgredraw.x b/noao/imred/dtoi/hdicfit/hdicgredraw.x
new file mode 100644
index 00000000..d42a0740
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgredraw.x
@@ -0,0 +1,22 @@
+include	<gset.h>
+
+define	MSIZE	2.0
+
+# HD_REDRAW -- Redraw a data point.  The old position marker is erased,
+# the new marker drawn, and the cursor moved to the new location.
+
+procedure hd_redraw (gp, oldx, oldy, newx, newy)
+
+pointer	gp		# Pointer to graphics stream
+real	oldx		# Old x coordinate
+real	oldy		# Old y coordinate
+real	newx		# New x coordinate
+real	newy		# New y coordinate
+
+begin
+	call gseti (gp, G_PMLTYPE, GL_CLEAR)
+	call gmark (gp, oldx, oldy, GM_PLUS, MSIZE, MSIZE)
+	call gseti (gp, G_PMLTYPE, GL_SOLID)
+	call gmark (gp, newx, newy, GM_PLUS, MSIZE, MSIZE)
+	call gscur (gp, newx, newy)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgsample.x b/noao/imred/dtoi/hdicfit/hdicgsample.x
new file mode 100644
index 00000000..12eef158
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgsample.x
@@ -0,0 +1,84 @@
+include	<gset.h>
+include	<pkg/rg.h>
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+# ICG_SAMPLE -- Mark sample.
+
+procedure icg_sampled (ic, gp, gt, x, npts, pltype)
+
+pointer	ic				# ICFIT pointer
+pointer	gp				# GIO pointer
+pointer	gt				# GTOOLS pointer
+double	x[npts]				# Ordinates of graph
+int	npts				# Number of data points
+int	pltype				# Plot line type
+
+pointer	rg
+int	i, axis, pltype1
+real	xl, xr, yb, yt, dy
+real	x1, x2, y1, y2, y3
+
+int	gstati(), stridxs(), gt_geti()
+pointer	rg_xrangesd()
+
+begin
+	if (stridxs ("*", Memc[IC_SAMPLE(ic)]) > 0)
+	    return
+
+	# Find axis along which the independent data is plotted.
+	if (IC_AXES(ic,IC_GKEY(ic),1) == 'x')
+	    axis = 1
+	else if (IC_AXES(ic,IC_GKEY(ic),2) == 'x')
+	    axis = 2
+	else
+	    return
+
+	if (gt_geti (gt, GTTRANSPOSE) == YES)
+	    axis = mod (axis, 2) + 1
+
+	pltype1 = gstati (gp, G_PLTYPE)
+	call gseti (gp, G_PLTYPE, pltype)
+	rg = rg_xrangesd (Memc[IC_SAMPLE(ic)], x, npts)
+
+	switch (axis) {
+	case 1:
+	    call ggwind (gp, xl, xr, yb, yt)
+
+	    dy = yt - yb
+	    y1 = yb + dy / 100
+	    y2 = y1 + dy / 20
+	    y3 = (y1 + y2) / 2
+
+	    do i = 1, RG_NRGS(rg) {
+	        x1 = x[RG_X1(rg, i)]
+	        x2 = x[RG_X2(rg, i)]
+	        if ((x1 > xl) && (x1 < xr))
+	            call gline (gp, x1, y1, x1, y2)
+	        if ((x2 > xl) && (x2 < xr))
+	            call gline (gp, x2, y1, x2, y2)
+	        call gline (gp, x1, y3, x2, y3)
+	    }
+	case 2:
+	    call ggwind (gp, yb, yt, xl, xr)
+
+	    dy = yt - yb
+	    y1 = yb + dy / 100
+	    y2 = y1 + dy / 20
+	    y3 = (y1 + y2) / 2
+
+	    do i = 1, RG_NRGS(rg) {
+	        x1 = x[RG_X1(rg, i)]
+	        x2 = x[RG_X2(rg, i)]
+	        if ((x1 > xl) && (x1 < xr))
+	            call gline (gp, y1, x1, y2, x1)
+	        if ((x2 > xl) && (x2 < xr))
+	            call gline (gp, y1, x2, y2, x2)
+	        call gline (gp, y3, x1, y3, x2)
+	    }
+	}
+
+	call gseti (gp, G_PLTYPE, pltype1)
+	call rg_free (rg)
+end
+
diff --git a/noao/imred/dtoi/hdicfit/hdicguaxes.x b/noao/imred/dtoi/hdicfit/hdicguaxes.x
new file mode 100644
index 00000000..f8199c87
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicguaxes.x
@@ -0,0 +1,38 @@
+include	"hdicfit.h"
+
+# ICG_UAXIS -- Set user axis.
+
+procedure icg_uaxesd (ic, key, cv, x, y, z, npts, label, units, maxchars)
+
+pointer	ic				# Pointer to ic structure
+int	key				# Key for axes
+pointer	cv				# CURFIT pointer
+double	x[npts]				# Independent variable
+double	y[npts]				# Dependent variable
+double	z[npts]				# Output values
+int	npts				# Number of points
+char	label[maxchars]			# Axis label
+char	units[maxchars]			# Units for axis
+int	maxchars			# Maximum chars in label
+
+int	offset
+double	fog
+real	ic_getr()
+include	"hdic.com"
+
+begin
+	# Axis type 'u' returns the untransformed independent variable
+	# in the z array.  That is, the original density values after
+	# subtracting the current fog value.  Some density values could be
+	# below fog were excluded from the transformed vector. 
+
+	call strcpy ("Density above fog", label, maxchars)
+	fog = double (ic_getr (ic, "fog"))
+
+	if (npts == nraw)
+	    call asubkd (Memd[den], fog, z, npts)
+	else {
+	    offset = big_den + (NVALS_FIT - npts)
+	    call asubkd (Memd[offset], fog, z, npts)
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgundel.x b/noao/imred/dtoi/hdicfit/hdicgundel.x
new file mode 100644
index 00000000..9a7c68a4
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgundel.x
@@ -0,0 +1,87 @@
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+define	MSIZE		2.		# Mark size
+
+# ICG_UNDELETE -- Undelete data point nearest the cursor.
+# The nearest point to the cursor in NDC coordinates is determined.
+
+procedure icg_undeleted (ic, gp, gt, cv, x, y, wts, userwts, npts, wx, wy)
+
+pointer	ic					# ICFIT pointer
+pointer	gp					# GIO pointer
+pointer	gt					# GTOOLS pointer
+pointer	cv					# CURFIT pointer
+double	x[npts], y[npts]			# Data points
+double	wts[npts], userwts[npts]		# Weight arrays
+int	npts					# Number of points
+real	wx, wy					# Position to be nearest
+
+pointer	sp, xout, yout
+int	gt_geti()
+
+begin
+	call smark (sp)
+	call salloc (xout, npts, TY_DOUBLE)
+	call salloc (yout, npts, TY_DOUBLE)
+
+	call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts)
+	call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts)
+	if (gt_geti (gt, GTTRANSPOSE) == NO) {
+	    call icg_u1d (ic, gp, Memd[xout], Memd[yout], wts, userwts,
+		npts, wx, wy)
+	} else {
+	    call icg_u1d (ic, gp, Memd[yout], Memd[xout], wts, userwts,
+		npts, wy, wx)
+	}
+
+	call sfree (sp)
+end
+
+
+# ICG_U1D -- Do the actual undelete.
+
+procedure icg_u1d (ic, gp, x, y, wts, userwts, npts, wx, wy)
+
+pointer	ic					# ICFIT pointer
+pointer	gp					# GIO pointer
+double	x[npts], y[npts]			# Data points
+double	wts[npts], userwts[npts]		# Weight arrays
+int	npts					# Number of points
+real	wx, wy					# Position to be nearest
+
+int	i, j
+real	x0, y0, r2, r2min
+
+begin
+	# Transform world cursor coordinates to NDC.
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+
+	# Search for nearest point to a point with zero weight.
+	r2min = MAX_REAL
+	do i = 1, npts {
+	    if (wts[i] != 0.)
+		next
+
+	    call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0)
+
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		j = i
+	    }
+	}
+
+	# Unmark the deleted point and reset the weight.
+	if (j != 0) {
+	    call gscur (gp, real (x[j]), real (y[j]))
+	    call gseti (gp, G_PMLTYPE, GL_CLEAR)
+	    call gmark (gp, real (x[j]), real (y[j]), GM_CROSS, MSIZE, MSIZE)
+	    call gseti (gp, G_PMLTYPE, GL_SOLID)
+	    call gmark (gp, real (x[j]), real (y[j]), GM_PLUS, MSIZE, MSIZE)
+	    wts[j] = userwts[j]
+	    IC_NEWWTS(ic) = YES
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicguser.x b/noao/imred/dtoi/hdicfit/hdicguser.x
new file mode 100644
index 00000000..5e070537
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicguser.x
@@ -0,0 +1,17 @@
+# ICG_USER -- User default action
+
+procedure icg_user (ic, gp, gt, cv, wx, wy, wcs, key, cmd)
+
+pointer	ic			# ICFIT pointer
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+pointer	cv			# CURFIT pointer
+real	wx, wy			# Cursor positions
+int	wcs			# GIO WCS
+int	key			# Cursor key
+char	cmd[ARB]		# Cursor command
+
+begin
+	# Ring bell
+	call printf ("\07")
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicgvec.x b/noao/imred/dtoi/hdicfit/hdicgvec.x
new file mode 100644
index 00000000..81909cb8
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicgvec.x
@@ -0,0 +1,74 @@
+include	<mach.h>
+include	"hdicfit.h"
+
+# HDIC_GVEC -- Get a vector of data to be plotted.  The raw density
+# values are stored in common.  From these values, all possible transforms
+# can be generated.
+
+procedure hdic_gvec (ic, xout, npts, transform)
+
+pointer	ic		# Pointer to ic structure
+double	xout[npts]	# Output vector
+int	npts		# Desired npoints to output - possibly changed on output
+int	transform	# Integer code for desired type of transform
+
+pointer	sp, bigdenaf
+int	i, j, npts_desired
+double	fog, dval
+real	ic_getr()
+include	"hdic.com"
+
+begin
+	npts_desired = npts
+	if (npts_desired != NVALS_FIT) 
+	    call error (0, "hdicgvec: nvals != NVALS_FIT")
+
+	call smark (sp)
+	call salloc (bigdenaf, npts, TY_DOUBLE)
+
+	j = 1
+
+	fog = double (ic_getr (ic, "fog"))
+	call asubkd (Memd[big_den], fog, Memd[bigdenaf], npts)
+
+	switch (transform) {
+	case HD_NONE:
+	    call amovd (Memd[bigdenaf], xout, NVALS_FIT)
+
+	case HD_LOGO:
+	    do i = 1, npts_desired {
+		dval = Memd[bigdenaf+i-1]
+		if (dval < 0.0D0)
+		    npts = npts - 1
+		else {
+		    xout[j] = log10 ((10.**dval) - 1.0)
+		    j = j + 1
+		}
+	    }
+
+	case HD_K75:
+	    do i = 1, npts_desired {
+		dval = Memd[bigdenaf+i-1]
+		if (dval < 0.0D0)
+		    npts = npts - 1
+		else {
+		    xout[j] = dval + 0.75 * log10 (1.0 - (10.** (-dval)))
+		    j = j + 1
+		}
+	    }
+
+	case HD_K50:
+	    do i = 1, npts_desired {
+		dval = Memd[bigdenaf+i-1]
+		if (dval < 0.0D0)
+		    npts = npts - 1
+		else {
+		    xout[j] = dval + 0.50 * log10 (1.0 - (10.** (-dval)))
+		    j = j + 1
+		}
+	    }
+
+	}
+
+	call sfree (sp)
+end
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
diff --git a/noao/imred/dtoi/hdicfit/hdicparams.x b/noao/imred/dtoi/hdicfit/hdicparams.x
new file mode 100644
index 00000000..1840dccd
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicparams.x
@@ -0,0 +1,323 @@
+include "hdicfit.h"
+
+define	FUNCTIONS	"|chebyshev|legendre|spline3|spline1|power|"
+define	TRANSFORMS	"|none|logopacitance|k50|k75|"
+
+# IC_OPEN -- Open ICFIT parameter structure.
+
+procedure ic_open (ic)
+
+pointer	ic		# ICFIT pointer
+errchk	malloc
+
+begin
+	# Allocate memory for the package parameter structure.
+	call malloc (ic, IC_LENSTRUCT, TY_STRUCT)
+	call malloc (IC_SAMPLE(ic), SZ_LINE, TY_CHAR)
+	call malloc (IC_LABELS(ic,1), SZ_LINE, TY_CHAR)
+	call malloc (IC_LABELS(ic,2), SZ_LINE, TY_CHAR)
+	call malloc (IC_UNITS(ic,1), SZ_LINE, TY_CHAR)
+	call malloc (IC_UNITS(ic,2), SZ_LINE, TY_CHAR)
+
+	# Set defaults.
+	call ic_pstr (ic, "function", "spline3")
+	call ic_puti (ic, "order", 1)
+	call ic_pstr (ic, "sample", "*")
+	call ic_puti (ic, "naverage", 1)
+	call ic_puti (ic, "niterate", 0)
+	call ic_putr (ic, "low", 0.)
+	call ic_putr (ic, "high", 0.)
+	call ic_putr (ic, "grow", 0.)
+	call ic_pstr (ic, "xlabel", "X")
+	call ic_pstr (ic, "ylabel", "Y")
+	call ic_pstr (ic, "xunits", "")
+	call ic_pstr (ic, "yunits", "")
+	call ic_puti (ic, "key", 1)
+	call ic_pkey (ic, 1, 'x', 'y')
+	call ic_pkey (ic, 2, 'y', 'x')
+	call ic_pkey (ic, 3, 'x', 'r')
+	call ic_pkey (ic, 4, 'x', 'd')
+	call ic_pkey (ic, 5, 'x', 'n')
+
+	# Initialize other parameters
+	IC_RG(ic) = NULL
+	IC_XFIT(ic) = NULL
+	IC_YFIT(ic) = NULL
+	IC_WTSFIT(ic) = NULL
+	IC_REJPTS(ic) = NULL
+end
+
+
+# IC_CLOSER -- Close ICFIT parameter structure.
+
+procedure ic_closer (ic)
+
+pointer	ic		# ICFIT pointer
+
+begin
+	# Free memory for the package parameter structure.
+	call rg_free (IC_RG(ic))
+	call mfree (IC_XFIT(ic), TY_REAL)
+	call mfree (IC_YFIT(ic), TY_REAL)
+	call mfree (IC_WTSFIT(ic), TY_REAL)
+	call mfree (IC_REJPTS(ic), TY_INT)
+	call mfree (IC_SAMPLE(ic), TY_CHAR)
+	call mfree (IC_LABELS(ic,1), TY_CHAR)
+	call mfree (IC_LABELS(ic,2), TY_CHAR)
+	call mfree (IC_UNITS(ic,1), TY_CHAR)
+	call mfree (IC_UNITS(ic,2), TY_CHAR)
+	call mfree (ic, TY_STRUCT)
+end
+
+
+# IC_CLOSED -- Close ICFIT parameter structure.
+
+procedure ic_closed (ic)
+
+pointer	ic		# ICFIT pointer
+
+begin
+	# Free memory for the package parameter structure.
+	call rg_free (IC_RG(ic))
+	call mfree (IC_XFIT(ic), TY_DOUBLE)
+	call mfree (IC_YFIT(ic), TY_DOUBLE)
+	call mfree (IC_WTSFIT(ic), TY_DOUBLE)
+	call mfree (IC_REJPTS(ic), TY_INT)
+	call mfree (IC_SAMPLE(ic), TY_CHAR)
+	call mfree (IC_LABELS(ic,1), TY_CHAR)
+	call mfree (IC_LABELS(ic,2), TY_CHAR)
+	call mfree (IC_UNITS(ic,1), TY_CHAR)
+	call mfree (IC_UNITS(ic,2), TY_CHAR)
+	call mfree (ic, TY_STRUCT)
+end
+
+
+# IC_PSTR -- Put string valued parameters.
+
+procedure ic_pstr (ic, param, str)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be put
+char	str[ARB]		# String value
+
+int	i
+pointer	ptr
+int	strdic()
+bool	streq()
+
+begin
+	if (streq (param, "sample"))
+	    call strcpy (str, Memc[IC_SAMPLE(ic)], SZ_LINE)
+	else if (streq (param, "function")) {
+	    call malloc (ptr, SZ_LINE, TY_CHAR)
+	    i = strdic (str, Memc[ptr], SZ_LINE, FUNCTIONS)
+	    if (i > 0)
+		IC_FUNCTION(ic) = i
+	    call mfree (ptr, TY_CHAR)
+	} else if (streq (param, "transform")) {
+	    call malloc (ptr, SZ_LINE, TY_CHAR)
+	    i = strdic (str, Memc[ptr], SZ_LINE, TRANSFORMS)
+	    if (i > 0)
+		IC_TRANSFORM(ic) = i
+	    call mfree (ptr, TY_CHAR)
+	} else if (streq (param, "xlabel"))
+	    call strcpy (str, Memc[IC_LABELS(ic,1)], SZ_LINE)
+	else if (streq (param, "ylabel"))
+	    call strcpy (str, Memc[IC_LABELS(ic,2)], SZ_LINE)
+	else if (streq (param, "xunits"))
+	    call strcpy (str, Memc[IC_UNITS(ic,1)], SZ_LINE)
+	else if (streq (param, "yunits"))
+	    call strcpy (str, Memc[IC_UNITS(ic,2)], SZ_LINE)
+	else
+	    call error (0, "ICFIT: Unknown parameter")
+end
+
+
+# IC_PUTI -- Put integer valued parameters.
+
+procedure ic_puti (ic, param, ival)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be put
+int	ival			# Integer value
+
+bool	streq()
+
+begin
+	if (streq (param, "naverage"))
+	    IC_NAVERAGE(ic) = ival
+	else if (streq (param, "order"))
+	    IC_ORDER(ic) = ival
+	else if (streq (param, "niterate"))
+	    IC_NITERATE(ic) = ival
+	else if (streq (param, "key"))
+	    IC_GKEY(ic) = ival
+	else if (streq (param, "transform"))
+	    IC_TRANSFORM(ic) = ival
+	else
+	    call error (0, "ICFIT: Unknown parameter")
+end
+
+
+# IC_PKEY -- Put key parameters.
+# Note the key types must be integers not characters.
+
+procedure ic_pkey (ic, key, xaxis, yaxis)
+
+pointer	ic			# ICFIT pointer
+int	key			# Key to be defined
+int	xaxis			# X axis type
+int	yaxis			# Y axis type
+
+begin
+	IC_AXES(ic, key, 1) = xaxis
+	IC_AXES(ic, key, 2) = yaxis
+end
+
+
+# IC_GKEY -- Get key parameters.
+
+procedure ic_gkey (ic, key, xaxis, yaxis)
+
+pointer	ic			# ICFIT pointer
+int	key			# Key to be gotten
+int	xaxis			# X axis type
+int	yaxis			# Y axis type
+
+begin
+	xaxis = IC_AXES(ic, key, 1)
+	yaxis = IC_AXES(ic, key, 2)
+end
+
+
+# IC_PUTR -- Put real valued parameters.
+
+procedure ic_putr (ic, param, rval)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be put
+real	rval			# Real value
+
+bool	streq()
+
+begin
+	if (streq (param, "xmin"))
+	    IC_XMIN(ic) = rval
+	else if (streq (param, "xmax"))
+	    IC_XMAX(ic) = rval
+	else if (streq (param, "low"))
+	    IC_LOW(ic) = rval
+	else if (streq (param, "high"))
+	    IC_HIGH(ic) = rval
+	else if (streq (param, "grow"))
+	    IC_GROW(ic) = rval
+	else if (streq (param, "fog"))
+	    IC_FOG(ic) = rval
+	else if (streq (param, "rfog"))
+	    IC_RFOG(ic) = rval
+	else
+	    call error (0, "ICFIT: Unknown parameter")
+end
+
+
+# IC_GSTR -- Get string valued parameters.
+
+procedure ic_gstr (ic, param, str, maxchars)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be put
+char	str[maxchars]		# String value
+int	maxchars		# Maximum number of characters
+
+bool	streq()
+
+begin
+	if (streq (param, "sample"))
+	    call strcpy (Memc[IC_SAMPLE(ic)], str, maxchars)
+	else if (streq (param, "xlabel"))
+	    call strcpy (Memc[IC_LABELS(ic,1)], str, maxchars)
+	else if (streq (param, "ylabel"))
+	    call strcpy (Memc[IC_LABELS(ic,2)], str, maxchars)
+	else if (streq (param, "xunits"))
+	    call strcpy (Memc[IC_UNITS(ic,1)], str, maxchars)
+	else if (streq (param, "yunits"))
+	    call strcpy (Memc[IC_UNITS(ic,2)], str, maxchars)
+	else if (streq (param, "function")) {
+	    switch (IC_FUNCTION(ic)) {
+	    case 1:
+		call strcpy ("chebyshev", str, maxchars)
+	    case 2:
+		call strcpy ("legendre", str, maxchars)
+	    case 3:
+		call strcpy ("spline3", str, maxchars)
+	    case 4:
+		call strcpy ("spline1", str, maxchars)
+	    case 5:
+		call strcpy ("power", str, maxchars)
+	    }
+	} else if (streq (param, "transform")) {
+	    switch (IC_TRANSFORM(ic)) {
+	    case 1:
+		call strcpy ("none", str, maxchars)
+	    case 2:
+		call strcpy ("logopacitance", str, maxchars)
+	    case 3:
+		call strcpy ("k50", str, maxchars)
+	    case 4:
+		call strcpy ("k75", str, maxchars)
+	    }
+	} else
+	    call error (0, "ICFIT: Unknown parameter")
+end
+
+
+# IC_GETI -- Get integer valued parameters.
+
+int procedure ic_geti (ic, param)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be gotten
+
+bool	streq()
+
+begin
+	if (streq (param, "naverage"))
+	    return (IC_NAVERAGE(ic))
+	else if (streq (param, "order"))
+	    return (IC_ORDER(ic))
+	else if (streq (param, "niterate"))
+	    return (IC_NITERATE(ic))
+	else if (streq (param, "key"))
+	    return (IC_GKEY(ic))
+	else if (streq (param, "transform"))
+	    return (IC_TRANSFORM(ic))
+
+	call error (0, "ICFIT: Unknown parameter")
+end
+
+
+# IC_GETR -- Get real valued parameters.
+
+real procedure ic_getr (ic, param)
+
+pointer	ic			# ICFIT pointer
+char	param[ARB]		# Parameter to be put
+
+bool	streq()
+
+begin
+	if (streq (param, "xmin"))
+	    return (IC_XMIN(ic))
+	else if (streq (param, "xmax"))
+	    return (IC_XMAX(ic))
+	else if (streq (param, "low"))
+	    return (IC_LOW(ic))
+	else if (streq (param, "high"))
+	    return (IC_HIGH(ic))
+	else if (streq (param, "grow"))
+	    return (IC_GROW(ic))
+	else if (streq (param, "fog"))
+	    return (IC_FOG(ic))
+
+	call error (0, "ICFIT: Unknown parameter")
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicreject.x b/noao/imred/dtoi/hdicfit/hdicreject.x
new file mode 100644
index 00000000..82dd093c
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicreject.x
@@ -0,0 +1,39 @@
+# IC_REJECT -- Reject points with large residuals from the fit.
+#
+# The sigma of the fit residuals is calculated.  The rejection thresholds
+# are set at low_reject*sigma and high_reject*sigma.  Points outside the
+# rejection threshold are rejected from the fit and flagged in the rejpts
+# array.  Finally, the remaining points are refit.
+
+procedure ic_rejectd (cv, x, y, w, rejpts, npts, low_reject, high_reject,
+	niterate, grow, nreject)
+
+pointer	cv				# Curve descriptor
+double	x[npts]				# Input ordinates
+double	y[npts]				# Input data values
+double	w[npts]				# Weights
+int	rejpts[npts]			# Points rejected
+int	npts				# Number of input points
+real	low_reject, high_reject		# Rejection threshold
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection radius
+int	nreject				# Number of points rejected
+
+int	i, newreject
+
+begin
+	# Initialize rejection.
+	nreject = 0
+	call amovki (NO, rejpts, npts)
+
+	if (niterate <= 0)
+	    return
+
+	# Find deviant points.
+	do i = 1, niterate {
+	    call ic_deviantd (cv, x, y, w, rejpts, npts, low_reject,
+		high_reject, grow, YES, nreject, newreject)
+	    if (newreject == 0)
+		break
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicshow.x b/noao/imred/dtoi/hdicfit/hdicshow.x
new file mode 100644
index 00000000..521f04d3
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicshow.x
@@ -0,0 +1,52 @@
+include	<pkg/gtools.h>
+include	"hdicfit.h"
+
+# IC_SHOW -- Show the values of the parameters.
+
+procedure ic_show (ic, file, gt)
+
+pointer	ic			# ICFIT pointer
+char	file[ARB]		# Output file
+pointer	gt			# GTOOLS pointer
+
+int	fd
+pointer	str
+int	open()
+long	clktime()
+errchk	open, malloc
+
+begin
+	fd = open (file, APPEND, TEXT_FILE)
+	call malloc (str, SZ_LINE, TY_CHAR)
+
+	call cnvtime (clktime(0), Memc[str], SZ_LINE)
+	call fprintf (fd, "\n# %s\n")
+	    call pargstr (Memc[str])
+
+	call gt_gets (gt, GTTITLE, Memc[str], SZ_LINE)
+	call fprintf (fd, "# %s\n")
+	    call pargstr (Memc[str])
+
+	call gt_gets (gt, GTYUNITS, Memc[str], SZ_LINE)
+	if (Memc[str] != EOS) {
+	    call fprintf (fd, "fit units = %s\n")
+	        call pargstr (Memc[str])
+	}
+
+	call ic_gstr (ic, "function", Memc[str], SZ_LINE)
+	call fprintf (fd, "function = %s\n")
+	    call pargstr (Memc[str])
+
+	call fprintf (fd, "order = %d\n")
+	    call pargi (IC_ORDER(ic))
+
+	call ic_gstr (ic, "transform", Memc[str], SZ_LINE)
+	call fprintf (fd, "transform = %s\n")
+	    call pargstr (Memc[str])
+
+	call fprintf (fd, "fog = %g\n")
+	    call pargr (IC_FOG(ic))
+
+	call mfree (str, TY_CHAR)
+	call close (fd)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicsort.x b/noao/imred/dtoi/hdicfit/hdicsort.x
new file mode 100644
index 00000000..16d6e0a6
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicsort.x
@@ -0,0 +1,38 @@
+# HDIC_SORT -- sort the log exposure, density and weight information in order
+# of increasing density value.  The sorting is done is place.  The four
+# data values are assummed matched on input, that is, exposure[i] matches
+# density[i] with weight[i] (and userwts[i]) for all array entries.
+
+procedure hdic_sort (density, exposure, weights, userwts, whydel, sdev, nvals)
+
+double	density[nvals]		# Density array
+double	exposure[nvals]		# Log exposure array
+double	weights[nvals]		# Weights array
+double	userwts[nvals]		# Reference weights array
+int	whydel[nvals]		# Flag array of reasons for deletion
+double	sdev[nvals]		# Array of standard deviations
+int	nvals			# Number of values to sort
+
+int	i, j
+double	temp
+define	swap	{temp=$1;$1=$2;$2=temp}
+int	itemp
+define	iswap   {itemp=$1;$1=$2;$2=itemp}
+
+begin
+	# Bubble sort - inefficient, but sorting is done infrequently
+	# an expected small sample size (16 pts typically). 
+
+	for (i = nvals; i > 1; i = i - 1)
+	    for (j = 1; j < i; j = j + 1) 
+		if (density [j] > density[j+1]) {
+
+		    # Out of order; exchange values
+		    swap (exposure[j], exposure[j+1])
+		    swap ( density[j],  density[j+1])
+		    swap ( weights[j],  weights[j+1])
+		    swap ( userwts[j],  userwts[j+1])
+		   iswap (  whydel[j],   whydel[j+1])
+		    swap (    sdev[j],     sdev[j+1])
+		}
+end
diff --git a/noao/imred/dtoi/hdicfit/hdictrans.x b/noao/imred/dtoi/hdicfit/hdictrans.x
new file mode 100644
index 00000000..0ee89037
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdictrans.x
@@ -0,0 +1,155 @@
+include	<mach.h>
+include	"hdicfit.h"
+
+# HDIC_TRANSFORM -- Transform density to independent variable of fit.  The
+# desired transform is stored in the ic structure.  A vector of x values
+# is returned, as is a possibly modified weights array.  The minimum and
+# maximum limits of the fit are updated in the ic structure; the labels
+# are set also when IC_NEWTRANSFORM = YES.  The fog value is subtracted
+# from the input density array and the transform performed. 
+
+procedure hdic_transform (ic, density, userwts, xout, wts, whydel, npts)
+
+pointer	ic		# Pointer to ic structure
+double	density[npts]	# Array of original density values
+double	userwts[npts]	# Array of original weights values
+double	xout[npts]	# Transformed density above fog (returned)
+double	wts[npts]	# Input weights array
+int	whydel[npts]	# Reason for deletion array 
+int	npts		# The number of density points - maybe changed on output
+
+int	i
+pointer	denaf, sp
+double	fog, xxmin, xxmax, dval
+bool	fp_equald()
+real	ic_getr()
+include	"hdic.com"
+
+begin
+	# Allocate space for density above fog array
+	call smark (sp)
+	call salloc (denaf, npts, TY_DOUBLE)
+
+	fog = double (ic_getr (ic, "fog"))
+	call asubkd (density, fog, Memd[denaf], npts)
+
+	switch (IC_TRANSFORM(ic)) {
+	case HD_NONE:
+	    do i = 1, npts {
+		xout[i] = Memd[denaf+i-1]
+		# In every case, if the point was deleted by the program, 
+		# restore it.
+		if (whydel[i] == PDELETE) {
+	            wts[i] = userwts[i]
+		    whydel[i] = NDELETE
+	        }
+	    }
+
+	    call ic_pstr (ic, "xlabel", "Density above Fog")
+	    xxmin = MIN_DEN - fog
+	    xxmax = maxden
+	    call ic_putr (ic, "xmin", real (xxmin))
+	    call ic_putr (ic, "xmax", real (xxmax))
+
+	case HD_LOGO:
+	    call ic_pstr (ic, "xlabel", "Log Opacitance: Log (10**Den - 1)")
+	    xxmin = log10 ((10. ** (MIN_DEN)) - 1.0)
+	    xxmax = log10 ((10. ** (maxden)) - 1.0)
+	    call ic_putr (ic, "xmin", real (xxmin))
+	    call ic_putr (ic, "xmax", real (xxmax))
+
+	    do i = 1, npts {
+		dval = Memd[denaf+i-1]
+		if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) {
+		    xout[i] = dval
+		    wts[i] = 0.0D0
+		    whydel[i] = PDELETE
+
+		} else {
+		    xout[i] = log10 ((10. ** (dval)) - 1.0)
+
+		    # If point had been deleted, find out why.  It affects the
+		    # weights value returned.  Only if the point was previously
+		    # deleted by the program, restore it; otherwise, leave it 
+		    # alone.
+
+		    if (fp_equald (wts[i], 0.0D0)) {
+		        if (whydel[i] == PDELETE) {
+		            wts[i] = userwts[i]
+			    whydel[i] = NDELETE
+			}
+		    } else
+			wts[i] = userwts[i]
+	    }
+	}
+
+	case HD_K75:
+	    call ic_pstr (ic, "xlabel", "Den + 0.75 * Log (1 - (10 ** -Den))")
+	    xxmin = MIN_DEN + 0.75 * log10 (1.0 - (10. ** (-MIN_DEN)))
+	    xxmax = maxden + 0.75 * log10 (1.0 - (10. ** (-maxden)))
+	    call ic_putr (ic, "xmin", real (xxmin))
+	    call ic_putr (ic, "xmax", real (xxmax))
+
+	    do i = 1, npts {
+		dval = Memd[denaf+i-1]
+		if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) {
+		    xout[i] = dval
+		    wts[i] = 0.0D0
+		    whydel[i] = PDELETE
+
+		} else {
+		    xout[i] = dval + 0.75 * log10 (1.0 - (10.** (-dval)))
+
+		    # If point had been deleted, find out why.  It affects the
+		    # weights value returned.  Only if the point was previously
+		    # deleted by the program, restore it; otherwise, leave it 
+		    # alone.
+
+		    if (fp_equald (wts[i], 0.0D0)) {
+		        if (wts[i] == PDELETE) {
+		            wts[i] = userwts[i]
+			    whydel[i] = NDELETE
+			}
+		    } else
+			wts[i] = userwts[i]
+	    }
+	}
+
+	case HD_K50:
+	    call ic_pstr (ic, "xlabel", "Den + 0.50 * Log (1 - (10 ** -Den))")
+	    xxmin = MIN_DEN + 0.50 * log10 (1.0 - (10. ** (-MIN_DEN)))
+	    xxmax = maxden + 0.50 * log10 (1.0 - (10. ** (-maxden)))
+	    call ic_putr (ic, "xmin", real (xxmin))
+	    call ic_putr (ic, "xmax", real (xxmax))
+
+	    do i = 1, npts {
+		dval = Memd[denaf+i-1]
+		if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) {
+		    xout[i] = dval
+		    wts[i] = 0.0D0
+		    whydel[i] = PDELETE
+
+		} else {
+		    xout[i] = dval + 0.50 * log10 (1.0 - (10.** (-dval)))
+
+		    # If point had been deleted, find out why.  It affects the
+		    # weights value returned.  Only if the point was previously
+		    # deleted by the program, restore it; otherwise, leave it 
+		    # alone.
+
+		    if (fp_equald (wts[i], 0.0D0)) {
+		        if (wts[i] == PDELETE) {
+		            wts[i] = userwts[i]
+			    whydel[i] = NDELETE
+			}
+		    } else
+			wts[i] = userwts[i]
+		}
+	    }
+
+	default:
+	    call eprintf ("Unrecognizable Transform\n")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/imred/dtoi/hdicfit/hdicvshow.x b/noao/imred/dtoi/hdicfit/hdicvshow.x
new file mode 100644
index 00000000..e62f6522
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicvshow.x
@@ -0,0 +1,155 @@
+include	<math/curfit.h>
+include	"hdicfit.h"
+
+# IC_VSHOW -- Show fit parameters in verbose mode.
+
+procedure ic_vshowd (ic, file, cv, x, y, wts, npts, gt)
+
+pointer	ic		# ICFIT pointer
+char	file[ARB]	# Output file
+pointer	cv		# Curfit pointer
+double	x[ARB]		# Ordinates
+double	y[ARB]		# Abscissas
+double	wts[ARB]	# Weights
+int	npts		# Number of data points
+pointer	gt		# Graphics tools pointer
+
+double	chisqr, rms
+int	i, n, deleted, ncoeffs, fd
+pointer	sp, fit, wts1, coeffs, errors
+
+int	dcvstati(), open()
+double	ic_rmsd()
+errchk	open()
+
+begin
+	# Do the standard ic_show option, then add on the verbose part.
+	call ic_show (ic, file, gt)
+
+	if (npts == 0) {
+	    call eprintf ("Incomplete output - no data points for fit\n")
+	    return
+	}
+
+	# Open the output file.
+	fd = open (file, APPEND, TEXT_FILE)
+
+	# Determine the number of coefficients and allocate memory.
+	ncoeffs = dcvstati (cv, CVNCOEFF)
+
+	call smark (sp)
+	call salloc (coeffs, ncoeffs, TY_DOUBLE)
+	call salloc (errors, ncoeffs, TY_DOUBLE)
+
+	if (npts == IC_NFIT(ic)) {
+	    # Allocate memory for the fit.
+	    n = npts
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (wts, Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+		do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Get the coefficients and compute the errors.
+	    call dcvvector (cv, x, Memd[fit], n)
+	    call dcvcoeff (cv, Memd[coeffs], ncoeffs)
+	    call dcverrors (cv, y, Memd[wts1], Memd[fit], n, chisqr,
+		Memd[errors])
+	    rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n)
+
+	} else {
+	    # Allocate memory for the fit.
+	    n = IC_NFIT(ic)
+	    call salloc (fit, n, TY_DOUBLE)
+	    call salloc (wts1, n, TY_DOUBLE)
+
+	    # Eliminate rejected points and count deleted points.
+	    call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n)
+	    if (IC_NREJECT(ic) > 0) {
+	        do i = 1, npts {
+		    if (Memi[IC_REJPTS(ic)+i-1] == YES)
+			Memd[wts1+i-1] = 0.
+		}
+	    }
+	    deleted = 0
+	    do i = 1, n {
+		if (wts[i] == 0.)
+		    deleted = deleted + 1
+	    }
+
+	    # Get the coefficients and compute the errors.
+	    call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n)
+	    rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)],
+		Memd[fit], Memd[wts1], n)
+	    call dcvcoeff (cv, Memd[coeffs], ncoeffs)
+	    call dcverrors (cv, Memd[IC_YFIT(ic)], Memd[wts1], Memd[fit],
+		n, chisqr, Memd[errors])
+	}
+
+	# Print the error analysis.
+	call fprintf (fd, "total points = %d\n")
+	    call pargi (npts)
+	call fprintf (fd, "deleted = %d\n")
+	    call pargi (deleted)
+	call fprintf (fd, "RMS = %7.4g\n")
+	    call pargd (rms)
+	call fprintf (fd, "square root of reduced chi square = %7.4g\n")
+	    call pargd (sqrt (chisqr))
+
+	call fprintf (fd, "# \t  coefficent\t  error\n")
+	do i = 1, ncoeffs {
+	    call fprintf (fd, "\t%10.4e\t%10.4e\n")
+		call pargd (Memd[coeffs+i-1])
+	 	call pargd (Memd[errors+i-1])
+	}
+
+	# Print x,y pairs and weights
+	call ic_listxywd (fd, cv, x, y, wts, npts)
+
+	call sfree (sp)
+	call close (fd)
+end
+
+
+# IC_LISTXYW -- List data as x,y pairs on output with their weights.  Used 
+# for verbose show procedure.  The untransformed density is also output,
+# regardless of what transformation may have been applied.
+
+procedure ic_listxywd (fd, cv, xvals, yvals, weights, nvalues)
+
+int	fd			# File descriptor of output file
+pointer	cv			# Pointer to curfit structure
+int	nvalues			# Number of data values
+double	xvals[nvalues]		# Array of x data values
+double	yvals[nvalues]		# Array of y data values
+double	weights[nvalues]	# Array of weight values
+
+int	i
+double	dcveval()
+include	"hdic.com"
+
+begin
+	call fprintf (fd,"\n#%15t Density %27t X %39t Yfit %51t LogE %63tWts\n")
+
+	do i = 1, nvalues {
+	    call fprintf (fd, 
+	    "%2d %15t%-12.7f%27t%-12.7f%39t%-12.7f%51t%-12.7f%63t%-12.7f\n")
+		call pargi (i)
+		call pargd (Memd[den+i-1])
+		call pargd (xvals[i])
+		call pargd (dcveval (cv, xvals[i]))
+		call pargd (yvals[i])
+		call pargd (weights[i])
+	}
+end
diff --git a/noao/imred/dtoi/hdicfit/mkpkg b/noao/imred/dtoi/hdicfit/mkpkg
new file mode 100644
index 00000000..5fc5031f
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/mkpkg
@@ -0,0 +1,37 @@
+# HDICFIT package.  *** Modified from icgfit package for DTOI applications ***
+
+update:
+	$update libhdic.a
+
+libhdic.a:
+
+	hdicclean.x	hdicfit.h <pkg/rg.h>
+	hdicdeviant.x	<mach.h> <math/curfit.h>
+	hdicdosetup.x	hdicfit.h <math/curfit.h>
+	hdicerrors.x	hdicfit.h <math/curfit.h>
+	hdicfit.x	hdicfit.h <math/curfit.h>
+	hdicgaxes.x	hdicfit.h <pkg/gtools.h>
+	hdicgcolon.x	hdicfit.h <error.h> <gset.h>
+	hdicgdelete.x	hdicfit.h <gset.h> <mach.h> <pkg/gtools.h>
+	hdicgfit.x	hdicfit.h <error.h> <pkg/gtools.h> <mach.h> <gset.h>
+	hdicggraph.x	hdicfit.h hdic.com <gset.h> <pkg/gtools.h> <mach.h>
+	hdicgnearest.x	<mach.h> <pkg/gtools.h>
+	hdicgparams.x	hdicfit.h <pkg/gtools.h>
+	hdicgsample.x	hdicfit.h <gset.h> <pkg/gtools.h> <pkg/rg.h>
+	hdicguaxes.x	hdic.com hdicfit.h
+	hdicgundel.x	hdicfit.h <gset.h> <mach.h> <pkg/gtools.h>
+	hdicguser.x	
+	hdicparams.x	hdicfit.h
+	hdicreject.x	
+	hdicshow.x	hdicfit.h <pkg/gtools.h>
+	hdicvshow.x	hdicfit.h hdic.com <math/curfit.h>
+
+	hdictrans.x	hdicfit.h <mach.h> hdic.com
+	hdicgvec.x	hdicfit.h <mach.h> hdic.com
+	hdicadd.x	hdicfit.h
+	hdicinit.x 	hdic.com hdicfit.h <mach.h>
+	hdicsort.x
+	userfcn.x	<error.h>
+	hdicgredraw.x	<gset.h>
+	hdicebars.x	hdicfit.h hdic.com <pkg/gtools.h> <gset.h> <mach.h>
+	;
diff --git a/noao/imred/dtoi/hdicfit/userfcn.x b/noao/imred/dtoi/hdicfit/userfcn.x
new file mode 100644
index 00000000..d5dba4ed
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/userfcn.x
@@ -0,0 +1,37 @@
+include	<error.h>
+
+# HD_POWERR -- Construct the basis functions for a power series function.
+# Invoked from curfit as a user function.  Real version.
+
+procedure hd_powerr (x, order, k1, k2, basis)
+
+real	x		# array of data points
+int	order		# order of polynomial, order = 1, constant
+real	k1, k2		# normalizing constants - unused
+real	basis[ARB]	# basis functions
+
+int	i
+
+begin
+	do i = 1, order
+	    iferr (basis[i] = x ** (i-1))
+		call erract (EA_FATAL) 
+end
+
+
+# HD_POWERD -- Double version of above.
+
+procedure hd_powerd (x, order, k1, k2, basis)
+
+double	x		# array of data points
+int	order		# order of polynomial, order = 1, constant
+double	k1, k2		# normalizing constants - unused
+double	basis[ARB]	# basis functions
+
+int	i
+
+begin
+	do i = 1, order
+            iferr (basis[i] = x ** (i-1))
+	        call erract (EA_FATAL)
+end
diff --git a/noao/imred/dtoi/hdshift.par b/noao/imred/dtoi/hdshift.par
new file mode 100644
index 00000000..c1fdf5c9
--- /dev/null
+++ b/noao/imred/dtoi/hdshift.par
@@ -0,0 +1,2 @@
+# Cl parameters for task hdshift are:
+database,s,a,"",,,List of database file
diff --git a/noao/imred/dtoi/hdshift.x b/noao/imred/dtoi/hdshift.x
new file mode 100644
index 00000000..cde846ad
--- /dev/null
+++ b/noao/imred/dtoi/hdshift.x
@@ -0,0 +1,184 @@
+include	<math/curfit.h>
+
+# T_HDSHIFT -- Task hdshift in the dtoi package.  This task is provided
+# to support Kormendy's method of combining related characteristic curves.
+# A zero point shift in log exp unique to each set of spots is calculated and
+# subtracted.  A single curve is fit to the combined, shifted data in
+# a separate task (hdfit).
+
+procedure t_hdshift ()
+
+pointer	sp, de, fun, save, db, cv, ps_coeff, den, exp, cvf
+int	fd, rec, nsave, ncoeff, nvalues, i, nfile, junk
+real	a0, ref_a0
+
+pointer	ddb_map()
+int	clpopni(), ddb_locate(), ddb_geti(), cvstati(), strncmp(), clgfil()
+
+begin
+	# Allocate space on stack for string buffers
+	call smark (sp)
+	call salloc (de, SZ_FNAME, TY_CHAR)
+	call salloc (cvf, SZ_FNAME, TY_CHAR)
+	call salloc (fun, SZ_FNAME, TY_CHAR)
+
+	# Get list of the database names.  The curfit information is retrieved
+	# from the first file in the list, the list is then rewound.
+
+	fd = clpopni ("database")
+	junk = clgfil (fd, Memc[cvf], SZ_FNAME)
+	call clprew (fd)
+
+	# Get coefficients of common fit from cv_file
+	db = ddb_map (Memc[cvf], READ_ONLY)
+	rec = ddb_locate (db, "cv")
+	nsave = ddb_geti (db, rec, "save")
+	call salloc (save, nsave, TY_REAL)
+	call ddb_gar (db, rec, "save", Memr[save], nsave, nsave)
+	call ddb_gstr (db, rec, "function", Memc[fun], SZ_LINE)
+
+	call cvrestore (cv, Memr[save])
+	ncoeff = cvstati (cv, CVNCOEFF)
+	call salloc (ps_coeff, ncoeff, TY_REAL)
+
+	if (strncmp (Memc[fun], "power", 5) == 0)
+	    call cvcoeff (cv, Memr[ps_coeff], ncoeff)
+	else
+	    call cvpower (cv, Memr[ps_coeff], ncoeff)
+
+	do i = 1, ncoeff {
+	    call eprintf ("%d   %.7g\n")
+	        call pargi (i)
+	        call pargr (Memr[ps_coeff+i-1])
+	}
+
+	call ddb_unmap (db)
+
+	nfile = 0
+	while (clgfil (fd, Memc[de], SZ_FNAME) != EOF) {
+	    db = ddb_map (Memc[de], READ_ONLY)
+	    call hds_read (db, den, exp, nvalues)
+	    call hds_calc (den, exp, nvalues, Memr[ps_coeff], ncoeff, a0)
+	    nfile = nfile + 1
+	    if (nfile == 1)
+		ref_a0 = a0
+	    a0 = a0 - ref_a0
+
+	    call printf ("file %s: subtracting zero point a0 = %.7g\n")
+		call pargstr (Memc[de])
+		call pargr (a0)
+
+	    # Write new log exposure information to database
+	    db = ddb_map (Memc[de], APPEND)
+	    call hds_wdb  (db, exp, nvalues, a0)
+	    call mfree (den, TY_REAL)
+	    call mfree (exp, TY_REAL)
+	    call ddb_unmap (db)
+	}
+
+	call clpcls (fd)
+	call sfree (sp)
+end
+
+
+# HDS_READ -- Read the density and exposure values from the database file.
+# The density above fog and log exposure values are returned, as well as
+# the number of data pairs read.
+
+procedure hds_read (db, den, exp, nvalues)
+
+pointer	db			# Pointer to input database file
+pointer	den			# Pointer to density array - returned
+pointer	exp			# Pointer to exposure array - returned
+int	nvalues			# Number of data pairs read - returned
+
+real	fog
+int	nden, nexp, rec
+int	ddb_locate(), ddb_geti()
+real	ddb_getr()
+
+begin
+	# Get fog value to be subtracted from density
+	rec = ddb_locate (db, "fog")
+	fog = ddb_getr (db, rec, "density")
+
+	# Get density array
+	rec = ddb_locate (db, "density")
+	nden = ddb_geti (db, rec, "den_val")
+	call malloc (den, nden, TY_REAL)
+	call ddb_gar (db, rec, "den_val", Memr[den], nden, nden)
+	call asubkr (Memr[den], fog, Memr[den], nden)
+
+	# Get exposure array
+	rec = ddb_locate (db, "exposure")
+	nexp = ddb_geti (db, rec, "log_exp")
+	call malloc (exp, nexp, TY_REAL)
+	call ddb_gar (db, rec, "log_exp", Memr[exp], nexp, nexp)
+
+	nvalues = min (nden, nexp)
+end
+
+
+# HDS_CALC -- Calculate the individual shift, a0.
+
+procedure hds_calc (den, exp, nvalues, ps_coeff, ncoeff, a0)
+
+pointer	den
+pointer	exp
+int	nvalues
+real	ps_coeff[ARB]
+int	ncoeff
+real	a0
+
+int	i
+real	yavg, ycalc, xavg
+
+begin
+	# Calculate average density and log exposure values
+	xavg = 0.0
+	yavg = 0.0
+
+	do i = 1, nvalues {
+	    xavg = xavg + Memr[den+i-1]
+	    yavg = yavg + Memr[exp+i-1]
+	}
+
+	xavg = xavg / real (nvalues)
+	yavg = yavg / real (nvalues)
+
+	ycalc = 0.0
+	do i = 2, ncoeff
+	    ycalc = ycalc + ps_coeff[i] * (xavg ** real (i-1))
+
+	# Subtraction yields the zero point shift in question
+	a0 = yavg - ycalc
+end
+
+
+# HDS_WDB -- Write shifted log exposure values to database.
+
+procedure hds_wdb  (db, exp, nvalues, a0)
+
+pointer	db		# Pointer to database
+pointer	exp		# Pointer to array of exposure values
+int	nvalues		# Number of exposure values in sample
+real	a0		# Shift to be subtracted
+
+pointer	sp, expsub
+
+begin
+	call smark (sp)
+	call salloc (expsub, nvalues, TY_REAL)
+
+	call ddb_ptime (db)
+	call ddb_prec (db, "exposure")
+
+	call eprintf ("a0 = %g\n")
+	    call pargr (a0)
+
+	call asubkr (Memr[exp], a0, Memr[expsub], nvalues)
+	call ddb_par (db, "log_exp", Memr[expsub], nvalues)
+
+	call ddb_putr (db, "A0 shift", a0)
+	call sfree (sp)
+end
diff --git a/noao/imred/dtoi/hdtoi.par b/noao/imred/dtoi/hdtoi.par
new file mode 100644
index 00000000..e26cde83
--- /dev/null
+++ b/noao/imred/dtoi/hdtoi.par
@@ -0,0 +1,11 @@
+# Cl parameters are:
+input,f,a,,,,List of images to be transformed
+output,f,a,,,,List of output image names
+database,f,a,,,,Database containing fit parameters
+fog,s,h,"",,,Value of fog - read from database if unspecified
+sigma,r,h,3.0,,,Rejection criteria for determining mean fog
+floor,r,h,0.0,,,Value assigned to levels below fog
+ceiling,r,h,30000.,,,Scale highest density to this intensity value
+datatype,s,h,r,,,Pixel type of output image
+option,s,h,mean,"mean|median",,Choice of fog algorithm (mean or median)
+verbose,b,h,yes,,,Print log of processing to STDOUT
diff --git a/noao/imred/dtoi/hdtoi.x b/noao/imred/dtoi/hdtoi.x
new file mode 100644
index 00000000..5b550ea8
--- /dev/null
+++ b/noao/imred/dtoi/hdtoi.x
@@ -0,0 +1,407 @@
+include	<imhdr.h>
+include	<mach.h>
+include	<math/curfit.h>
+include	<error.h>
+include	"hdicfit/hdicfit.h"
+
+# T_HDTOI -- transform an image from density to intensity, according
+# to an hd curve described in an input database.  A look up table of
+# all possible values is generated with the curfit package, and
+# then the image is transformed line by line.  A fog value is subtracted
+# from the image prior to transformation, and it can be entered as either
+# a number or a list of fog images from which the fog value is calculated.
+# If a fog value has not been entered by the user, it is read from the database.
+
+procedure t_hdtoi ()
+
+pointer	sp, cv, fog, db, im_in, lut, im_out, imageout, imagein, option
+bool	verbose
+int	minval, maxval, in_list, rec, ip, out_list, fog_list, ngpix
+int	datatype, nluv, updatedb, nfpix
+real	sigma, floor, scale, fog_val, sdev
+
+char	clgetc()
+bool	streq(), clgetb()
+pointer	ddb_map(), immap()
+int	imtopenp(), ddb_locate(), ctor(), imtlen(), imtgetim()
+int	get_data_type(), imtopen()
+real	clgetr(), ddb_getr()
+
+begin
+	call smark (sp)
+	call salloc (cv,       SZ_FNAME, TY_CHAR)
+	call salloc (fog,      SZ_LINE,  TY_CHAR)
+	call salloc (imageout, SZ_FNAME, TY_CHAR)
+	call salloc (imagein,  SZ_FNAME, TY_CHAR)
+
+	# Get cl parameters
+	in_list = imtopenp ("input")
+	out_list = imtopenp ("output")
+	call clgstr ("database", Memc[cv], SZ_FNAME)
+	call clgstr ("fog", Memc[fog], SZ_LINE)
+	sigma = clgetr ("sigma")
+	floor = clgetr ("floor")
+	verbose = clgetb ("verbose")
+	updatedb = NO
+
+	datatype = get_data_type (clgetc ("datatype"))
+	if (datatype == ERR)
+	    call eprintf ("Using input pixel datatype for output\n")
+
+	db = ddb_map (Memc[cv], READ_ONLY)
+	rec = ddb_locate (db, "common")
+	scale = ddb_getr (db, rec, "scale")
+
+	# If not specified by user, get fog value from database.  User can
+	# specify fog as a real number or a list of fog file names.
+
+	if (streq (Memc[fog], "")) {
+	    rec = ddb_locate (db, "fog")
+	    fog_val = ddb_getr (db, rec, "density")
+	} else {
+	    ip = 1
+	    if (ctor (Memc[fog], ip, fog_val) == 0) {
+		if (verbose)
+		    call eprintf ("Calculating fog value ...\n")
+		fog_list = imtopen (Memc[fog])
+		call salloc (option, SZ_FNAME, TY_CHAR)
+		call clgstr ("option", Memc[option], SZ_FNAME)
+	        call hd_fogcalc (fog_list, fog_val, sdev, ngpix, scale, sigma,
+		    Memc[option], nfpix)
+
+		call eprintf ("Fog density = %f, sdev = %f, ngpix = %d\n")
+		    call pargr (fog_val)
+		    call pargr (sdev)
+		    call pargi (ngpix)
+
+		updatedb = YES
+	    }
+	}
+
+	# Generate look up table.  First, the range of input values to
+	# calculate output values for must be determined.  Arguments
+	# minval and maxval are integers because we assume all input
+	# images are short integers.
+
+	call hd_glimits (in_list, minval, maxval)
+	nluv = (maxval - minval) + 1
+	call salloc (lut, nluv, TY_REAL)
+
+	if (verbose)
+	    call eprintf ("Generating look up table ...\n")
+	call hd_wlut (db, Memr[lut], minval, maxval, fog_val, floor)
+
+	# Loop through input images, applying transform
+	if (imtlen (in_list) != imtlen (out_list)) {
+	    call imtclose (in_list)
+	    call imtclose (out_list)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	while ((imtgetim (in_list, Memc[imagein], SZ_FNAME) != EOF) &&
+	    (imtgetim (out_list, Memc[imageout], SZ_FNAME) != EOF)) {
+
+	    iferr (im_in = immap (Memc[imagein], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next 
+	    }
+
+	    iferr (im_out = immap (Memc[imageout], NEW_COPY, im_in)) {
+		call imunmap (im_in)
+		call erract (EA_WARN)
+		next
+	    }
+
+	    if (verbose) {
+	        call eprintf ("Density to Intensity Transform: %s ===> %s\n")
+    		    call pargstr (Memc[imagein])
+    		    call pargstr (Memc[imageout])
+	    }
+
+	    call hd_transform (im_in, im_out, Memr[lut], nluv, minval, datatype)
+
+	    call imunmap (im_in)
+	    call imunmap (im_out)
+	}
+
+	call ddb_unmap (db)
+	call imtclose (in_list)
+	call imtclose (out_list)
+
+	if (updatedb == YES) {
+	    db = ddb_map (Memc[cv], APPEND)
+	    # Write fog information to database as single record
+	    call ddb_prec (db, "fog")
+	    call ddb_putr (db, "density", fog_val)
+	    call ddb_putr (db, "sdev", sdev)
+	    call ddb_puti (db, "ngpix", ngpix)
+	    call ddb_pstr (db, "option", Memc[option])
+	    call ddb_unmap (db)
+	}
+
+	call sfree (sp)
+end
+
+
+# HD_TRANSFORM -- Apply transformation to image.
+
+procedure hd_transform (im, im_out, lu_table, nvals, minval, datatype)
+
+pointer	im		# Input image header pointer
+pointer	im_out		# Transformed image header pointer
+real	lu_table[ARB]	# Array of intensity values
+int	nvals		# Number of values in the lut
+int	minval		# Offset to first value in look up table
+int	datatype	# Pixel type on output
+
+int	j, ncols
+pointer	ptr_in, ptr_out, sp, luti
+pointer	impl2r(), imgl2i(), impl2i()
+
+begin
+	if (datatype == ERR)
+	    IM_PIXTYPE(im_out) = IM_PIXTYPE(im)
+	else
+	    IM_PIXTYPE(im_out) = datatype
+
+	ncols = IM_LEN(im,1)
+
+	switch (datatype) {
+	case TY_REAL, TY_DOUBLE:
+	    # Loop over input image rows.  The look up table is left as
+	    # a real array and a floating point image is written out.
+
+	    do j = 1, IM_LEN(im,2) {
+	        ptr_in = imgl2i (im, j)
+	        ptr_out = impl2r (im_out, j)
+		call asubki (Memi[ptr_in], minval, Memi[ptr_in], ncols)
+		call alutr (Memi[ptr_in], Memr[ptr_out], ncols, lu_table)
+	    }
+
+	default:
+	    # Loop over input image rows. The look up table is truncated
+	    # to type integer.
+
+	    call smark (sp)
+	    call salloc (luti, nvals, TY_INT)
+	    call achtri (lu_table, Memi[luti], nvals)
+
+	    do j = 1, IM_LEN(im,2) {
+	        ptr_in = imgl2i (im, j)
+	        ptr_out = impl2i (im_out, j)
+		call asubki (Memi[ptr_in], minval, Memi[ptr_in], ncols)
+		call aluti (Memi[ptr_in], Memi[ptr_out], ncols, Memi[luti])
+	    }
+
+	    call sfree (sp)
+	}
+end
+
+
+# HD_WLUT -- write look up table, such that intensity = lut [a/d output].
+# An entry is made in the look up table for every possible input value,
+# from minval to maxval.
+
+procedure hd_wlut (db, lut, minval, maxval, fog_val, floor)
+
+pointer	db		# Pointer to database file
+real 	lut[ARB]	# Pointer to look up table, which gets filled here
+int	minval		# Minimum value to transform
+int	maxval		# Maximum value to transform
+real	fog_val		# Fog value to be subtracted from densities
+real	floor		# Value assigned to densities below fog
+
+bool	zerofloor
+pointer	sp, trans, fcn, save, cv, dens, ind_var, value
+int	rec, nsave, i, function, nneg, npos, nvalues
+real	scale, maxcvval, factor, maxexp, maxden, maxdenaf
+
+bool	fp_equalr()
+int	strncmp(), ddb_locate(), ddb_geti(), cvstati()
+real	ddb_getr(), clgetr(), cveval()
+extern	hd_powerr()
+
+begin
+	call smark (sp)
+	call salloc (trans, SZ_FNAME, TY_CHAR)
+	call salloc (fcn, SZ_FNAME, TY_CHAR)
+
+	nvalues = (maxval - minval) + 1
+	call salloc (ind_var, nvalues, TY_REAL)
+	call salloc (dens,    nvalues, TY_REAL)
+	call salloc (value,   nvalues, TY_REAL)
+
+	rec = ddb_locate (db, "common")
+	scale = ddb_getr (db, rec, "scale")
+	maxden = ddb_getr (db, rec, "maxden")
+
+	rec = ddb_locate (db, "cv")
+	nsave = ddb_geti (db, rec, "save")
+	call salloc (save, nsave, TY_REAL)
+	call ddb_gar (db, rec, "save", Memr[save], nsave, nsave)
+	call ddb_gstr (db, rec, "transformation", Memc[trans], SZ_LINE)
+
+	call cvrestore (cv, Memr[save])
+	function = cvstati (cv, CVTYPE)
+
+	if (function == USERFNC) {
+	    # Need to restablish choice of user function
+	    call ddb_gstr (db, rec, "function", Memc[fcn], SZ_FNAME)
+	    if (strncmp (Memc[fcn], "power", 1) == 0) 
+		call cvuserfnc (cv, hd_powerr)
+	    else
+		call error (0, "Unknown user function in database")
+	}
+
+	maxdenaf = maxden - fog_val
+	call hd_aptrans (maxdenaf, maxcvval, 1, Memc[trans])
+	maxcvval = 10.0 ** (cveval (cv, maxcvval))
+	factor = clgetr ("ceiling") / maxcvval
+	maxexp = real (MAX_EXPONENT) - (log10 (factor) + 1.0)
+
+	zerofloor = false
+	if (fp_equalr (0.0, floor))
+	    zerofloor = true
+
+	do i = 1, nvalues
+	    Memr[value+i-1] = real (minval + i - 1)
+
+	# Scale all posible voltage values to density above fog
+	call altmr (Memr[value], Memr[dens], nvalues, scale, -fog_val)
+
+	# Find index of first value greater than MIN_DEN.  Values less than 
+	# this must be handled as the user specified with the floor parameter.
+
+	for (nneg=0; Memr[dens+nneg] < MIN_DEN; nneg=nneg+1)
+	    ;
+	npos = nvalues - nneg
+	
+	# Generate independent variable vector and then lut values.  The
+	# logic is different if there are values below fog.  Evaluating 
+	# the polynomial fit with cvvector yields the log exposure.  This
+	# is then converted to intensity and scaled by a user supplied factor.
+
+	if (nneg > 0) {
+	    if (zerofloor) {
+		call amovkr (0.0, lut, nneg)
+		call hd_aptrans (Memr[dens+nneg],Memr[ind_var],npos,Memc[trans])
+		call cvvector (cv, Memr[ind_var], lut[nneg+1], npos)
+		call argtr (lut[nneg+1], npos, maxexp, maxexp)
+		do i = nneg+1, nvalues
+		    lut[i] = (10. ** lut[i]) * factor
+
+	    } else {
+		call amulkr (Memr[dens], -1.0, Memr[dens], nneg)
+
+		# Care must be taken so that no density of value 0.0 is 
+		# passed to hd_aptrans.  This would cause an overflow.
+
+		do i = 1, nneg {
+		    if (fp_equalr (Memr[dens], 0.0))
+			Memr[dens] = MIN_DEN
+		}
+
+		call hd_aptrans (Memr[dens],Memr[ind_var], nvalues, Memc[trans])
+		call cvvector (cv, Memr[ind_var], lut, nvalues)
+		call argtr (lut, nvalues, maxexp, maxexp)
+	        do i = 1, nvalues
+		    lut[i] = (10.0 ** lut[i]) * factor
+		call amulkr (lut, -1.0, lut, nneg)
+	    }
+		
+	} else {
+	    call hd_aptrans (Memr[dens], Memr[ind_var], nvalues, Memc[trans])
+	    call cvvector (cv, Memr[ind_var], lut, nvalues)
+	    call argtr (lut, nvalues, maxexp, maxexp)
+	    do i = 1, nvalues
+		lut[i] = (10.0 ** lut[i]) * factor
+	}
+
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# HD_APTRANS -- Apply transformation, generating a vector of independent
+# variables from a density vector.  It is assummed all values in the
+# input density vector are valid and will not cause arithmetic errors.
+# No checking for out of bounds values is performed.
+
+procedure hd_aptrans (density, ind_var, nvalues, transform)
+
+real	density[nvalues]	# Density vector - input
+real	ind_var[nvalues]	# Ind variable vector - filled on output
+int	nvalues			# Length of vectors
+char	transform[ARB]		# String containing transformation type
+
+int	i
+int	strncmp()
+
+begin
+	if (strncmp (transform, "logopacitance", 1) == 0) {
+	    do i = 1, nvalues
+	        ind_var[i] = log10 ((10. ** density[i]) - 1.0)
+
+	} else if (strncmp (transform, "k75", 2) == 0) {
+	    do i = 1, nvalues
+	        ind_var[i] = density[i] + .75 * log10(1. - 10. ** (-density[i]))
+
+	} else if (strncmp (transform, "k50", 2) == 0) {
+	    do i = 1, nvalues
+	        ind_var[i] = density[i] + .50 * log10(1. - 10. ** (-density[i]))
+
+	} else if (strncmp (transform, "none", 1) == 0) {
+	    do i = 1, nvalues
+	        ind_var[i] = density[i]
+
+	} else 
+	    call error (0, "Unrecognized transformation in database file")
+end
+
+
+# HD_GLIMITS -- Determinine the range of max and min values for a list
+# of images.
+
+procedure hd_glimits (in_list, minval, maxval)
+
+int	in_list		# File descriptor for list of images
+int	minval		# Smallest pixel value - returned
+int	maxval		# Largest  pixel value - returned
+
+pointer	im
+char	image[SZ_FNAME]
+real	current_min, current_max, min, max
+pointer	immap()
+int	imtgetim()
+errchk	imtgetim, im_minmax, imtrew
+
+begin
+	current_min = MAX_REAL
+	current_max = EPSILONR
+
+	while (imtgetim (in_list, image, SZ_FNAME) != EOF) {
+	    iferr (im = immap (image, READ_ONLY, 0))
+		# Just ignore it, warning will be printed by t_hdtoi
+		next
+
+	    # Update min max values if necessary
+	    if (IM_LIMTIME(im) < IM_MTIME(im))
+	        call im_minmax (im, IM_MIN(im), IM_MAX(im))
+
+	    min = IM_MIN(im)
+	    max = IM_MAX(im)
+
+	    if (min < current_min)
+		current_min = min
+
+	    if (max > current_max)
+		current_max = max
+	    
+	    call imunmap (im)
+	}
+
+	minval = int (current_min)
+	maxval = int (current_max)
+
+	call imtrew (in_list)
+end
diff --git a/noao/imred/dtoi/minmax.x b/noao/imred/dtoi/minmax.x
new file mode 100644
index 00000000..97d113f0
--- /dev/null
+++ b/noao/imred/dtoi/minmax.x
@@ -0,0 +1,73 @@
+include	<imhdr.h>
+
+# IM_MINMAX -- Compute the minimum and maximum pixel values of an image.
+# Works for images of any dimensionality, size, or datatype, although
+# the min and max values can currently only be stored in the image header
+# as real values.
+
+procedure im_minmax (im, min_value, max_value)
+
+pointer	im				# image descriptor
+real	min_value			# minimum pixel value in image (out)
+real	max_value			# maximum pixel value in image (out)
+
+pointer	buf
+bool	first_line
+long	v[IM_MAXDIM]
+short	minval_s, maxval_s
+long	minval_l, maxval_l
+real	minval_r, maxval_r
+int	imgnls(), imgnll(), imgnlr()
+errchk	amovkl, imgnls, imgnll, imgnlr, alims, aliml, alimr
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)		# start vector
+	first_line = true
+	min_value = INDEF
+	max_value = INDEF
+
+	switch (IM_PIXTYPE(im)) {
+	case TY_SHORT:
+	    while (imgnls (im, buf, v) != EOF) {
+		call alims (Mems[buf], IM_LEN(im,1), minval_s, maxval_s)
+		if (first_line) {
+		    min_value = minval_s
+		    max_value = maxval_s
+		    first_line = false
+		} else {
+		    if (minval_s < min_value)
+			min_value = minval_s
+		    if (maxval_s > max_value)
+			max_value = maxval_s
+		}
+	    }
+	case TY_USHORT, TY_INT, TY_LONG:
+	    while (imgnll (im, buf, v) != EOF) {
+		call aliml (Meml[buf], IM_LEN(im,1), minval_l, maxval_l)
+		if (first_line) {
+		    min_value = minval_s
+		    max_value = maxval_s
+		    first_line = false
+		} else {
+		    if (minval_l < min_value)
+			min_value = minval_l
+		    if (maxval_l > max_value)
+			max_value = maxval_l
+		}
+	    }
+	default:
+	    while (imgnlr (im, buf, v) != EOF) {
+		call alimr (Memr[buf], IM_LEN(im,1), minval_r, maxval_r)
+		if (first_line) {
+		    min_value = minval_r
+		    max_value = maxval_r
+		    first_line = false
+		} else {
+		    if (minval_r < min_value)
+			min_value = minval_r
+		    if (maxval_r > max_value)
+			max_value = maxval_r
+		}
+	    }
+	}
+end
diff --git a/noao/imred/dtoi/mkpkg b/noao/imred/dtoi/mkpkg
new file mode 100644
index 00000000..80d5b737
--- /dev/null
+++ b/noao/imred/dtoi/mkpkg
@@ -0,0 +1,40 @@
+# Make the DTOI package.
+
+$call	relink
+$exit
+
+update:
+	$call	relink
+	$call	install
+	;
+
+relink:
+	$call update@hdicfit
+	$update	libpkg.a
+	$call	dtoi
+	;
+
+install:
+	$move	xx_dtoi.e noaobin$x_dtoi.e
+	;
+
+dtoi:
+	$omake  x_dtoi.x
+	$link	x_dtoi.o libpkg.a hdicfit/libhdic.a -lxtools -lcurfit\
+		-o xx_dtoi.e
+	;
+
+libpkg.a:
+	database.x	<ctotok.h> <ctype.h> <finfo.h> <time.h>
+	dematch.x	<error.h>
+	hd_aravr.x	<mach.h>
+	hdfit.x		hdicfit/hdicfit.h <ctype.h> <error.h> <fset.h>\
+			<imhdr.h> <mach.h> <math/curfit.h> <pkg/gtools.h>\
+			<pkg/xtanswer.h>
+	hdshift.x	<math/curfit.h>
+	hdtoi.x		hdicfit/hdicfit.h <error.h> <imhdr.h> <mach.h>\
+			<math/curfit.h>
+	minmax.x	<imhdr.h>
+	selftest.x	<gio.h> <gset.h> <mach.h>
+	spotlist.x	<error.h> <fset.h> <imhdr.h> <mach.h>
+	;
diff --git a/noao/imred/dtoi/selftest.par b/noao/imred/dtoi/selftest.par
new file mode 100644
index 00000000..75445cf7
--- /dev/null
+++ b/noao/imred/dtoi/selftest.par
@@ -0,0 +1,8 @@
+# Parameters for task selftest
+nbits,i,q,12,,,Test data range
+device,s,h,"stdgraph",,,Output device
+verbose,b,h,no,,,"Print density, intensity values?"
+ceiling,r,h,30000.,,,Scale densest point to this intensity
+max_raw,i,q,0,,,Max raw data value
+scale,r,q,,0.0,,Raw value to density scale factor
+mode,s,h,ql
diff --git a/noao/imred/dtoi/selftest.x b/noao/imred/dtoi/selftest.x
new file mode 100644
index 00000000..3a7a4a4a
--- /dev/null
+++ b/noao/imred/dtoi/selftest.x
@@ -0,0 +1,290 @@
+include	<gio.h>
+include	<gset.h>
+include	<mach.h>
+
+define	A0	(-1.74743)
+define	A1	(0.73)
+define	A2	(-0.24)
+define	A3	(0.035)
+define	MAXDEN	6.0
+
+# T_SELFTEST -- a test procedure for the DTOI package.  Two intensity vectors
+# are calculated in different ways and compared.  A plot of the residuals is
+# shown.  A plot showing the extent of truncation errors is also drawn.  Two
+# standard ranges of data values are available: 12 bit, representing PDS
+# format data and 15 bit, representing the FITS data format available on the
+# PDS.  Any other choice results in a small test, ranging from 1 - 144.
+
+procedure t_selftest
+
+bool	verbose
+char	device[SZ_FNAME]
+pointer	sp, intk, intc, raw, den, gp
+int	min_raw, max_raw, nvalues, i, nbits
+real	scale, factor, ceiling
+
+bool	clgetb()
+pointer	gopen()
+int	clgeti()
+real	clgetr()
+
+begin
+	call smark (sp)
+
+	nbits = clgeti ("nbits")
+
+	switch (nbits) {
+	case 12:
+	    min_raw = 1
+	    max_raw = 3072
+	    scale   = 0.00151
+	case 15:
+	    min_raw = 1
+	    max_raw = 24576
+	    scale   = 4.65 / 24575. 
+	case 0:
+	    call eprintf ("Using test data range from 1 - 144\n")
+	    min_raw = 1
+	    max_raw = 144
+	    scale   = 0.0325
+	default:
+	    call eprintf ("Unknown case: nbits = '%d', Please supply values:\n")
+		call pargi (nbits)
+	    min_raw = 1
+	    max_raw = clgeti ("max_raw")
+	    # max density = 6.0.  Density = raw value * scale.
+	    call clputr ("scale.p_maximum", real (MAXDEN / max_raw))
+	    call clputr ("scale.p_default", real (4.65 / (max_raw - 1)))
+	    scale = clgetr ("scale")
+	}
+
+	call clgstr ("device", device, SZ_FNAME)
+	verbose = clgetb ("verbose")
+	ceiling = clgetr ("ceiling")
+
+	gp = gopen (device, NEW_FILE, STDGRAPH)
+
+	nvalues = max_raw - min_raw + 1
+	call salloc (intk, nvalues, TY_REAL)
+	call salloc (intc, nvalues, TY_REAL)
+	call salloc (den,  nvalues, TY_REAL)
+	call salloc (raw,  nvalues, TY_REAL)
+
+	do i = 1, nvalues
+	    Memr[raw+i-1] = min_raw + i - 1
+
+	call amulkr (Memr[raw], scale, Memr[den], nvalues)
+
+	call hd_known (min_raw, max_raw, scale, Memr[intk], nvalues)
+	call hd_calc  (min_raw, max_raw, scale, Memr[intc], nvalues)
+	 
+	if (verbose) {
+	    factor = ceiling / Memr[intc+nvalues-1]
+	    call printf ("# %20tRaw Value %40tDensity %60tIntensity\n\n")
+	    do i = 1, nvalues {
+	        call printf ("%20t%d %40t%g %60t%g\n")
+		    call pargi (i)
+		    call pargr (Memr[den+i-1])
+		    call pargr (Memr[intc+i-1] * factor)
+	    }
+	}
+
+	call hd_plotit (gp, Memr[den], Memr[intk], Memr[intc], nvalues)
+
+	call hd_trunc (gp, Memr[den], nvalues, ceiling, Memr[intc])
+
+	call gclose (gp)
+	call sfree (sp)
+end
+
+
+# HD_KNOWN -- Calculate vector of known intensity values.
+
+procedure hd_known (min_raw, max_raw, scale, intk, nvalues)
+
+int	min_raw			# Minimum raw data value
+int	max_raw			# Maximum raw data value
+real	scale			# Density = raw_value * scale
+real	intk[nvalues]		# Known intensities - filled on return
+int	nvalues			# Number of intensity values requested
+
+int	i
+real	density, logo
+real	exp
+
+begin
+	do i = min_raw, max_raw {
+	    density = max (EPSILONR, i * scale)
+	    logo = log10 ((10. ** density) - 1.0)
+	    exp  = A0 + A1 * logo + A2 * logo ** 2 + A3 * logo ** 3
+	    intk[i] = 10 ** (exp)
+	}
+end
+
+
+# HD_CALC -- Calcuate vector of intensity values as in HDTOI.
+
+procedure hd_calc (min, max, scale, intc, nvalues)
+
+int	min			# Minimum raw data value
+int	max			# Maximum raw data value
+real	scale			# Density = raw_value * scale
+real	intc[nvalues]		# Calculated intensity values - filled on return
+int	nvalues			# Number of intensity values requested
+
+real	cfit[9]
+pointer	sp, lut
+
+begin
+	call smark (sp)
+	call salloc (lut, nvalues, TY_REAL)
+
+	cfit[1] = 5.0
+	cfit[2] = 4.0
+	cfit[3] = -10.0
+	cfit[4] = MAXDEN
+	cfit[5] = 1.
+	cfit[6] = A0
+	cfit[7] = A1
+	cfit[8] = A2
+	cfit[9] = A3
+	call st_wlut (Memr[lut], min, max, scale, cfit)
+	call st_transform (min, max, Memr[lut], nvalues, intc)
+
+	call sfree (sp)
+end
+
+
+# HD_TRUNC -- Investigate truncation errors for real versus int output image.
+
+procedure hd_trunc (gp, density, nvalues, ceiling, intc)
+
+pointer	gp			# Pointer to graphics stream
+real	density[nvalues]	# Density array
+int	nvalues			# Number of density, intensity values
+real	ceiling			# Max intensity to output
+real	intc[nvalues]		# Calculated intensity values
+
+pointer	sp, yint, yreal
+int	npvals
+real	factor
+
+begin
+	call smark (sp)
+
+	# Only the lowest 5% of the data values are plotted
+	npvals = nvalues * 0.05
+
+	call salloc (yint, npvals, TY_INT)
+	call salloc (yreal, npvals, TY_REAL)
+
+	# Scale intensity vector to ceiling 
+	factor = ceiling / intc[nvalues]
+
+	call amulkr (intc, factor, intc, npvals)
+	call achtri (intc, Memi[yint], npvals)
+	call achtir (Memi[yint], Memr[yreal], npvals)
+
+	call gascale (gp, density,  npvals, 1)
+	call gascale (gp, Memr[yreal], npvals, 2)
+	call gsview  (gp, 0.2, 0.9, 0.1, 0.4)
+	call gseti (gp, G_ROUND, YES)
+	call glabax  (gp, 
+	    "Expand to see Truncation Errors\n (real=SOLID, integer=DASHED)", 
+	    "Density (Lowest 5% only)", "Intensity")
+
+	call gseti  (gp, G_PLTYPE, GL_SOLID)
+	call gpline (gp, density, intc, npvals)
+
+	call gseti  (gp, G_PLTYPE, GL_DASHED)
+	call gpline (gp, density, Memr[yreal], npvals)
+
+	call sfree (sp)
+end
+
+
+# HD_PLOTIT -- Plot residuals of calculated versus known itensity.
+
+procedure hd_plotit (gp, density, intk, intc, nvalues)
+
+pointer	gp			# Pointer to graphics stream
+real	density[nvalues]	# Density array
+real	intk[nvalues]		# Array of known intensities
+real	intc[nvalues]		# Array of calculated intensities
+int	nvalues			# Number of density, intensity values
+
+pointer	sp, resid
+
+begin
+	call smark (sp)
+	call salloc (resid, nvalues, TY_REAL)
+
+	call asubr (intk, intc, Memr[resid], nvalues)
+
+	call gascale (gp, density, nvalues, 1)
+	call gascale (gp, Memr[resid], nvalues, 2)
+	call gsview  (gp, 0.2, 0.9, 0.6, 0.9)
+	call gseti   (gp, G_ROUND, YES)
+
+	call glabax  (gp, "Residual Intensity\n (Known - Calculated)", 
+	    "Density", "")
+	call gpline  (gp, density, Memr[resid], nvalues)
+
+	call sfree (sp)
+end
+
+
+# ST_WLUT -- Generate look up table, using technique of HDTOI.
+
+procedure st_wlut (lut, min, max, scale, cfit)
+
+real	lut[ARB]	# Look up table of intensities
+int	min		# Minimum raw data value
+int	max		# Maximum raw data value
+real	scale		# Density = raw_value * scale
+real	cfit[ARB]	# Coefficient array for restoring curfit
+
+pointer	cv, sp, den, indv, kv
+int	nvalues, i
+extern	hd_powerr()
+
+begin
+	call smark (sp)
+	nvalues = max - min + 1
+	call salloc (den, nvalues, TY_REAL)
+	call salloc (indv, nvalues, TY_REAL)
+	call salloc (kv, nvalues, TY_REAL)
+	do i = 1, nvalues
+	    Memr[kv+i-1] = real (i)
+
+	call amulkr (Memr[kv], scale, Memr[den], nvalues)
+
+	call cvrestore (cv, cfit)
+	call cvuserfnc (cv, hd_powerr)
+
+	call hd_aptrans (Memr[den], Memr[indv], nvalues, "logo")
+	call cvvector (cv, Memr[indv], lut, nvalues)
+	do i = 1, nvalues
+	    lut[i] = 10.0 ** lut[i]
+
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# ST_TRANSFORM -- Apply transformation from look up table to input vector.
+
+procedure st_transform (min, max, lut, nvalues, intc)
+
+int	min		# Minimum raw data value
+int	max		# Maximum raw data value
+real	lut[ARB]	# Array of intensity values
+int	nvalues		# Number of density, intensity values
+real	intc[ARB]	# Calculated intensities - returned
+
+int	i
+
+begin
+	do i = 1, nvalues
+	    intc[i] = lut[i]
+end
diff --git a/noao/imred/dtoi/spotlist.par b/noao/imred/dtoi/spotlist.par
new file mode 100644
index 00000000..8134c8de
--- /dev/null
+++ b/noao/imred/dtoi/spotlist.par
@@ -0,0 +1,8 @@
+# The cl parameters for task spotlist in the dtoi package
+spots,f,a,,,,List of image files containing spot data
+fogs,f,a,,,,List of images containing fog spots
+database,f,a,,,,Name of output database file
+scale,r,h,0.00151,,,Input value to density scale factor
+maxad,i,h,3071,,,Integer A/D value of saturated pixel
+sigma,r,h,3.0,,,Rejection criteria when determining mean density
+option,s,h,mean,"mean|median",,Choice of algorithm (mean or median)
diff --git a/noao/imred/dtoi/spotlist.x b/noao/imred/dtoi/spotlist.x
new file mode 100644
index 00000000..8f4ffeee
--- /dev/null
+++ b/noao/imred/dtoi/spotlist.x
@@ -0,0 +1,395 @@
+include	<error.h>
+include	<imhdr.h>
+include	<mach.h>
+include	<fset.h>
+
+# T_SPOTLIST -- calculates the densities and standard deviations of calibration
+# spots read in as IRAF images.  Entries for density, sdev, spot number
+# and number of pixels used in the average are made in the output database.
+# These values are entered for all spots and the fog level.  The
+# fog level is calculated but not subtracted in this procedure.
+
+procedure t_spotlist ()
+
+pointer	db, sp, den, ngpix, sdev, dloge, option, npix
+int	spot_fd, fog_fd, fngpix, nspots, maxad, fnpix
+real	scale, sigma, fsdev, fog
+double	maxden
+
+pointer	ddb_map()
+int	imtopenp(), imtlen(), clgeti(), strncmp()
+real	clgetr()
+
+begin
+	call smark (sp)
+	call salloc (dloge, SZ_FNAME, TY_CHAR)
+	call salloc (option, SZ_FNAME, TY_CHAR)
+
+	# Get parameters and open file name templates
+	spot_fd = imtopenp ("spots")
+	fog_fd  = imtopenp ("fogs")
+	call clgstr ("database", Memc[dloge], SZ_FNAME)
+	scale = clgetr ("scale")
+	maxad = clgeti ("maxad")
+	sigma = clgetr ("sigma")
+	call clgstr ("option", Memc[option], SZ_FNAME)
+
+	maxden = maxad * double (scale)
+
+	# Allocate space for density, standard deviation and ngpix arrays
+	nspots = imtlen (spot_fd)
+	call salloc ( den,  nspots, TY_REAL)
+	call salloc (sdev,  nspots, TY_REAL)
+	call salloc (ngpix, nspots, TY_INT)
+	call salloc ( npix, nspots, TY_INT)
+
+	# Calculate densities depending on algorithm option.  The 
+	# number of saturated pixels per spot is also calculated now.
+
+	if (strncmp (Memc[option], "median", 3) == 0)
+	    call hd_median (spot_fd, Memr[den], Memr[sdev], Memi[ngpix],
+		nspots, scale, Memi[npix])
+	else
+	    call hd_mean (spot_fd, Memr[den], Memr[sdev], Memi[ngpix], 
+		nspots, scale, sigma, Memi[npix])
+
+	# Calculate fog level and count saturated pixels
+	call hd_fogcalc (fog_fd, fog, fsdev, fngpix, scale, sigma, 
+	    Memc[option], fnpix)
+
+	# Now print results to stdout
+	call hd_printit (Memr[den], Memr[sdev], Memi[npix], Memi[ngpix], 
+	    fog, fsdev, fnpix, fngpix, nspots)
+	    
+	# Open output database file and write spot information
+	db = ddb_map (Memc[dloge], APPEND)
+	call hd_wspotdb (db, Memr[den], Memr[sdev], Memi[ngpix], nspots)
+
+	# Write fog information to database as single record
+	call ddb_prec (db, "fog")
+	call ddb_putr (db, "density", fog)
+	call ddb_putr (db, "sdev", fsdev)
+	call ddb_puti (db, "ngpix", fngpix)
+
+	# Scale info gets written to database also (very precisely!)
+	call ddb_prec (db, "common")
+	call ddb_putr (db, "scale", scale)
+	call ddb_putd (db, "maxden", maxden)
+	call ddb_pstr (db, "option", Memc[option])
+
+	call ddb_unmap (db)
+
+	call imtclose (spot_fd)
+	call imtclose (fog_fd)
+	call sfree (sp)
+end
+
+
+# HD_MEAN -- Calculate mean density of calibration spots.
+
+procedure hd_mean (spot_fd, den, sdev, ngpix, nspots, scale, sigma, npix)
+
+int	spot_fd		# File descriptor for list of spots
+real	den[ARB]	# Mean density values - filled on return
+real	sdev[ARB]	# Standard deviation array - filled on return
+int	ngpix[ARB]	# Number of unrejected pixels - filled on return
+int	nspots		# Number of spots in list
+real	scale		# Scale for voltage to density conversion
+real	sigma		# Rejection criteria set by user
+int	npix[ARB]	# Number of pixels per spot
+
+pointer	im, spot, sp, pix
+int	i, junk, ncols, nlines
+pointer	immap(), imgs2r()
+int	imtgetim(), hd_aravr()
+errchk	imgs2r, amulkr, hd_aravr
+
+begin
+	call smark (sp)
+	call salloc (spot, SZ_FNAME, TY_CHAR)
+
+	# Loop over spot rasters.  Calculate density and standard deviation.
+	for (i = 1; i <= nspots; i = i + 1) {
+	    junk = imtgetim (spot_fd, Memc[spot], SZ_FNAME)
+	    iferr (im = immap (Memc[spot], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    ncols = IM_LEN(im,1)
+	    nlines = IM_LEN(im,2)
+	    npix[i] = ncols * nlines
+
+	    # For all pixels in the image, scale the a/d value to density and
+	    # calculate the mean value, using a mean rejection algorithm.
+
+	    pix = imgs2r (im, 1, ncols, 1, nlines)
+	    call amulkr (Memr[pix], scale, Memr[pix], npix[i])
+	    ngpix[i] = hd_aravr (Memr[pix], npix[i], den[i], sdev[i], sigma)
+	    call imunmap (im)
+	}
+
+	call sfree (sp)
+end
+
+
+# HD_MEDIAN -- Calculate median density of calibration spots.
+
+procedure hd_median (spot_fd, den, sdev, ngpix, nspots, scale, npix)
+
+int	spot_fd		# File descriptor for list of spots
+real	den[ARB]	# Mean density values - filled on return
+real	sdev[ARB]	# Standard deviation of input spots
+int	ngpix[ARB]	# Number of pixels not rejected
+int	nspots		# Number of spots in list
+real	scale		# Scale for voltage to density conversion
+int	npix[ARB]	# Number of pixels per spot
+
+pointer	im, spot, sp, pix
+int	i, junk, ncols, nlines
+real	junk_mean
+pointer	immap(), imgs2r()
+int	imtgetim()
+real	amedr()
+errchk	imgs2r, amulkr, amedr
+
+begin
+	call smark (sp)
+	call salloc (spot, SZ_FNAME, TY_CHAR)
+
+	# Loop over spot rasters.  Calculate density and standard deviation.
+	for (i = 1; i <= nspots; i = i + 1) {
+	    junk = imtgetim (spot_fd, Memc[spot], SZ_FNAME)
+	    iferr (im = immap (Memc[spot], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    ncols = IM_LEN(im,1)
+	    nlines = IM_LEN(im,2)
+	    npix[i] = ncols * nlines
+
+	    # For all pixels in the image, scale the a/d value to density and
+	    # calculate the median value.  For the user's information, the
+	    # sigma is also calculated in a call to aavgr.
+
+	    pix = imgs2r (im, 1, ncols, 1, nlines)
+	    call amulkr (Memr[pix], scale, Memr[pix], npix[i])
+	    den[i] = amedr (Memr[pix], npix[i])
+	    ngpix[i] = npix[i]
+	    call aavgr (Memr[pix], npix[i], junk_mean, sdev[i])
+	    call imunmap (im)
+	}
+
+	call sfree (sp)
+end
+
+
+# HD_FOGCALC -- Calculate fog density.
+
+procedure hd_fogcalc (fog_fd, fog, fsdev, fngpix, scale, sigma, option, nfpix)
+
+int	fog_fd		# File descriptor for list of fog spots
+real	fog		# Mean fog value - returned
+real	fsdev		# Standard deviation - returned
+int	fngpix		# Number of pixels used - returned
+real	scale		# Voltage to density scaling factor
+real	sigma		# Rejection criteria - set by user
+char	option[ARB]     # User's choice of mean/median algorithm
+int	nfpix		# Total number of fog pixels
+
+pointer	pix, im, sp, fogfile, ptr
+int	nfog, maxfog, i, junk, ncols, nlines, npix, total_pix, op
+real	junk_mean
+pointer	immap(), imgs2r()
+int	imtlen(), imtgetim(), hd_aravr(), strncmp()
+real	amedr()
+errchk	calloc, imgs2r, aaddr, amulkr, hd_aravr
+
+begin
+	call smark (sp)
+	call salloc (fogfile, SZ_FNAME, TY_CHAR)
+
+	pix = NULL
+	total_pix = 0
+	op = 0
+	nfog = imtlen (fog_fd)
+	maxfog = nfog
+	do i = 1, maxfog {
+	    junk = imtgetim (fog_fd, Memc[fogfile], SZ_FNAME)
+	    iferr (im = immap (Memc[fogfile], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		nfog = nfog - 1
+		next
+	    }
+
+	    ncols = IM_LEN(im,1)
+	    nlines = IM_LEN(im,2)
+	    npix = ncols * nlines
+	    total_pix = total_pix + npix
+
+	    if (pix == NULL)
+		# Initialize space for accumulating pixels
+		call calloc (pix, npix, TY_REAL)
+	    else
+	       # Increase space for accumulating pixels
+	       call realloc (pix, total_pix, TY_REAL)
+
+	    # Build up vector of fog pixels
+	    ptr = imgs2r (im, 1, ncols, 1, nlines)
+	    call amovr (Memr[ptr], Memr[pix+op], npix)
+	    op = op + npix
+
+	    call imunmap (im)
+	}
+
+	# Scale values to density and calculate fog and std deviation
+
+	if (nfog > 0) {
+# 7 Sept 1989, S. Rooke:  in Suzanne's absence, made following bugfix after
+# bug reported by Steve Majewski that fog values are off by 1/n images where
+# multiple fog images are used in a single run.  The total_pix already contains
+# the sum of all pixel values, so the fog pixel values should not be divided
+# by nfog.  This should be verified by Suzanne on her return, and these comments
+# removed.
+#	    call amulkr (Memr[pix], scale / real (nfog), Memr[pix], total_pix)
+	    call amulkr (Memr[pix], scale, Memr[pix], total_pix)
+	    if (strncmp (option, "median", 3) == 0) {
+		fog = amedr (Memr[pix], total_pix)
+		fngpix = total_pix
+		call aavgr (Memr[pix], total_pix, junk_mean, fsdev)
+	    } else
+	        fngpix = hd_aravr (Memr[pix], total_pix, fog, fsdev, sigma)
+	}  else {
+	    fog = 0.0
+	    fsdev = 0.0
+	    fngpix = 0
+	}
+	nfpix = total_pix
+
+	call mfree (pix, TY_REAL)
+	call sfree (sp)
+end
+
+
+# HD_WSPOTDB -- Write spot information to database file.  Values are first
+# sorted in order of increasing density.
+
+procedure hd_wspotdb (db, density, sdev, ngpix, nspots)
+
+pointer	db			# Pointer to database
+real	density[ARB]		# Array of densities
+real	sdev [ARB]		# Array of standard deviations
+int	ngpix [ARB]		# Array of npix used in calculations
+int	nspots			# Number of density spots
+
+begin
+	if (density[1] > density[nspots]) {
+	    # Need to reorder arrays
+	    call hd_reorderr (density, nspots)
+	    call hd_reorderr (   sdev, nspots)
+	    call hd_reorderi (  ngpix, nspots)
+	}
+
+	call ddb_ptime (db)
+
+	# Density record
+	call ddb_prec (db, "density")
+	call ddb_par (db, "den_val", density, nspots)
+
+	# Standard deviation of density is written to a record
+	call ddb_prec (db, "standard deviation")
+	call ddb_par (db, "sdev_val", sdev, nspots)
+
+	# Record for npix_used
+	call ddb_prec (db, "ngpix")
+	call ddb_pai (db, "npix_val", ngpix, nspots)
+end
+
+
+# HD_REORDERR - Flip order of real array in place.
+
+procedure hd_reorderr (array, nvals)
+
+real	array[ARB]		# Real array to be reordered
+int	nvals			# Number of elements in array
+
+pointer	sp, tmp
+int	i
+
+begin
+	call smark (sp)
+	call salloc (tmp, nvals, TY_REAL)
+
+	call amovr (array, Memr[tmp], nvals)
+	do i = 1, nvals
+	    array[i] = Memr[tmp+nvals-i]
+
+	call sfree (sp)
+end
+
+
+# HD_REORDERI -- Flip order of integer array in place.
+
+procedure hd_reorderi (array, nvals)
+
+int	array[ARB]		# Integer array to be ordered
+int	nvals			# Number of elements in array
+
+pointer	sp, tmp
+int	i
+
+begin
+	call smark (sp)
+	call salloc (tmp, nvals, TY_INT)
+
+	call amovi (array, Memi[tmp], nvals)
+	do i = 1, nvals
+	    array[i] = Memi[tmp+nvals-i]
+
+	call sfree (sp)
+end
+
+
+# HD_PRINTIT -- Neatly printout all the accumulated information.
+
+procedure hd_printit (den, sdev, npix, ngpix, fog, fsdev, fnpix, fngpix, nspots)
+
+real	  den[ARB]	# density array 
+real	 sdev[ARB]	# std deviation array 
+int	 npix[ARB]      # npix array
+int	ngpix[ARB]      # ngoodpix array
+real	fog		# fog value
+real	fsdev           # std deviation of fog
+int     fnpix		# npix in fog
+int	fngpix		# ngoodpix in fog
+int	nspots		# number of spots
+int	i
+
+begin
+
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	call printf ("\n                                         # Number of P")
+	call printf ("ixels\n")
+	call printf ("# Spot Number  Density   Std Deviation  Total  Used  Rej")
+	call printf ("ected \n")
+
+	do i = 1, nspots {
+	    call printf ("    %d %17t%.4f %27t%.4f %43t%d %5d %6d \n")
+	        call pargi (i)
+	        call pargr (den[i])
+	        call pargr (sdev[i])
+	        call pargi (npix[i])
+	        call pargi (ngpix[i])
+	        call pargi (npix[i] - ngpix[i])
+	}
+
+	call printf ("   FOG %17t%.4f %27t%.4f %43t%d %5d %6d \n")
+	    call pargr (fog)
+	    call pargr (fsdev)
+	    call pargi (fnpix)
+	    call pargi (fngpix)
+	    call pargi (fnpix - fngpix)
+
+end
diff --git a/noao/imred/dtoi/x_dtoi.x b/noao/imred/dtoi/x_dtoi.x
new file mode 100644
index 00000000..5fd0cc72
--- /dev/null
+++ b/noao/imred/dtoi/x_dtoi.x
@@ -0,0 +1,6 @@
+task spotlist = t_spotlist,
+	hdfit = t_hdfit,
+      hdshift = t_hdshift,
+	hdtoi = t_hdtoi,
+      dematch = t_dematch,
+     selftest = t_selftest