Initial commit

9 files changed, 1718 insertions, 0 deletions

diff --git a/noao/digiphot/photcal/fitparams/README b/noao/digiphot/photcal/fitparams/README
new file mode 100644
index 00000000..26c7617b
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/README
@@ -0,0 +1,2 @@
+This subdirectory contains the code for the task FITCOEFFS which computes
+the coefficients of a user-defined fitting function.
diff --git a/noao/digiphot/photcal/fitparams/fteval.com b/noao/digiphot/photcal/fitparams/fteval.com
new file mode 100644
index 00000000..bba04823
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/fteval.com
@@ -0,0 +1,8 @@
+# Fitting function and derivatives evaluation common
+
+pointer	eqcode		# equation code
+pointer	dercode		# equation derivative codes
+pointer	xpcode		# x plotting equation code
+pointer	ypcode		# y plotting equation code
+
+common	/ftevalcom/ eqcode, dercode, xpcode, ypcode
diff --git a/noao/digiphot/photcal/fitparams/fteval.x b/noao/digiphot/photcal/fitparams/fteval.x
new file mode 100644
index 00000000..cde0ff6a
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/fteval.x
@@ -0,0 +1,221 @@
+.help fteval
+Fit evaluation procedures.
+
+These procedures are called by the (I)NLFIT package to compute the fitting
+function, its derivatives at different points, and the plotting equations
+for all points.
+Pointers to the previous equation codes are stored in the evaluator common
+to speed up computations, since they are accessed for every data point.
+A buffer is allocated to store the derivative code pointers, since the
+number of derivatives is variable for each equation.
+.sp
+Code pointers are copied into the evaluator common by ft_evinit(). This
+procedure also allocates space for the derivative code pointers, before
+copying them. Space is freed by ft_evfree(). The fitting function is
+evaluated for a single point with ft_func(), derivatives are evaluated
+for a single point with ft_dfunc(), and the plot equations are evaluated
+for all points with ft_plot(). Plot equations are computed for each
+axis separately (two calls).
+
+Entry points:
+
+.nf
+
+    ft_evinit (sym, npars)			        Initialize fit
+    ft_evfree ()				        Free space allocated
+    ft_func (x, nvars, p, npars, z)		        Fitting function called
+    ft_dfunc (x, nvars, p, dp, npars, z, der)	        Derivatives of fit func.
+    ft_plot (number, p, npars, x, y, z, npts, nvars)    Plotting functions
+.fi
+.endhelp
+
+include	<error.h>
+include	"../lib/parser.h"
+
+# Pointer Mem
+define	MEMP	Memi
+
+
+# FT_EVINIT - Initialize fit by copying fitting function, its derivative,
+# and plotting function codes into the fit common.
+
+procedure ft_evinit (sym, npars)
+
+pointer	sym			# equation symbol
+int	npars			# number of parameters
+
+int	np			# parameter counter
+#bool	clgetb()
+pointer	pr_gsymp(), pr_gderp()
+
+include	"fteval.com"
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_evinit: (sym=%d) (npars=%d)\n")
+		#call pargi (sym)
+		#call pargi (npars)
+	#}
+
+	# Set fitting function code.
+	eqcode = pr_gsymp (sym, PTEQRPNFIT)
+	if (eqcode == NULL)
+	    call error (0, "ft_evinit: Null function equation code")
+
+	# Allocate space for derivative code pointers, and copy them
+	# from the symbol table. This is to avoid an unnecessary 
+	# overhead during the derivative evaluations.
+
+	call malloc (dercode, npars, TY_POINTER)
+	do np = 1, npars {
+	    MEMP[dercode + np - 1] = pr_gderp (sym, np, PTEQRPNDER)
+	    #if (MEMP[dercode + np - 1] == NULL)
+		#call error (0, "ft_evinit: Null derivative equation code")
+	}
+
+	# Set plotting equation codes. They could be null.
+	xpcode = pr_gsymp (sym, PTEQRPNXPLOT)
+	ypcode = pr_gsymp (sym, PTEQRPNYPLOT)
+
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_evinit: (eqcode=%d) ")
+		#call pargi (eqcode)
+	    #do np = 1, npars {
+	        #call eprintf (" (dercode%d=%d)")
+		    #call pargi (np)
+		    #call pargi (MEMP[dercode + np - 1])
+	    #}
+	    #call eprintf (" (xpcode=%d) (ypcode=%d)\n")
+		#call pargi (xpcode)
+		#call pargi (ypcode)
+	#}
+end
+
+
+# FT_EVFREE - Free space used in the fit common.
+
+procedure ft_evfree ()
+
+include	"fteval.com"
+
+#bool	clgetb()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf (
+	       #"ft_evfree: (eqcode=%d) (dercode=%d) (xpcode=%d) (ypcode=%d)\n")
+		#call pargi (eqcode)
+		#call pargi (dercode)
+		#call pargi (xpcode)
+		#call pargi (ypcode)
+	#}
+
+	# Clear code pointers.
+	eqcode = NULL
+	xpcode = NULL
+	ypcode = NULL
+
+	# Free derivative buffer.
+	call mfree (dercode, TY_POINTER)
+end
+
+
+# FT_FUNC - Evaluate fitting function. This function must conform in
+# number and type of arguments to the requirements of the nlfit/inlfit
+# packages even though some arguments are not used in this case.
+
+procedure ft_func (x, nvars, p, npars, z)
+
+real	x[ARB]			# independent values
+int	nvars			# number of variables
+real	p[ARB]			# parameter values
+int	npars			# number of parameters
+real	z			# function value (output)
+
+include	"fteval.com"
+
+real	pr_eval()
+
+begin
+	# Evaluate the function value.
+	z = pr_eval (eqcode, x, p)
+end
+
+
+# FT_DFUNC - Evaluate fitting function, and derivatives of the fitting
+# function with respect to all the parameters.  This function must conform in
+# number and type of arguments to the requirements of the nlfit/inlfit
+# packages even though some arguments are not used in this case.
+
+procedure ft_dfunc (x, nvars, p, dp, npars, z, der)
+
+real	x[ARB]			# independent values
+int	nvars			# number of variables
+real	p[ARB]			# parameter values
+real	dp[ARB]			# parameter derivatives
+int	npars			# number of parameters
+real	z			# function value (output)
+real	der[ARB]		# derivative values (output)
+
+int	i
+pointer	code
+real	pi, zplus, zminus
+real	pr_eval()
+
+include	"fteval.com"
+
+begin
+	# Evaluate the function value.
+	z = pr_eval (eqcode, x, p)
+
+	# Evaluate each one of the derivatives.
+	do i = 1, npars {
+	    code = MEMP[dercode+i-1]
+	    if (code != NULL)
+	        der[i] = pr_eval (code, x, p)
+	    else {
+		pi = p[i]
+		p[i] = pi + 0.5 * dp[i] 
+		zplus = pr_eval (eqcode, x, p)
+		p[i] = pi - 0.5 * dp[i]
+		zminus = pr_eval (eqcode, x, p)
+		der[i] = (zplus - zminus) / dp[i]
+		p[i] = pi
+	    }
+	}
+end
+
+
+# FT_PLOT - Evaluate plot function(s) for all points.
+
+procedure ft_plot (number, p, npars, x, y, z, npts, nvars)
+
+int	number			# plot number
+real	p[npars]		# parameter values
+int	npars			# number of parameters (not used)
+real	x[ARB]			# independent values
+real	y[npts]			# dependent values (not used)
+real	z[npts]			# function values (output)
+int	npts			# number of points
+int	nvars			# number of variables
+
+int	i
+pointer	code
+real	pr_eval()
+
+include	"fteval.com"
+
+begin
+	# Determine plot equation to evaluate.
+	if (number == 1)
+	    code = xpcode
+	else
+	    code = ypcode
+
+	# Iterate over input points.
+	do i = 1, npts
+	    z[i] = pr_eval (code, x[(i-1)*nvars+1], p)
+end
diff --git a/noao/digiphot/photcal/fitparams/ftindef.x b/noao/digiphot/photcal/fitparams/ftindef.x
new file mode 100644
index 00000000..943b82fd
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/ftindef.x
@@ -0,0 +1,184 @@
+include	<error.h>
+include	"../lib/fitparams.h"
+include	"../lib/parser.h"
+
+
+# FT_INDEF - Set zero weight for all undefined input data.
+
+procedure ft_indef (sym, otable, rtable, wtable)
+
+int	sym			# equation symbol
+pointer	otable			# 2d observation table (modified)
+pointer	rtable			# 1d reference table (modified)
+pointer	wtable			# 1d weight table (modified)
+
+#bool	clgetb()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf (
+	        #"ft_indef: (sym=%d) (otable=%d) (rtable=%d) (wtable=%d)\n")
+		#call pargi (sym)
+		#call pargi (otable)
+		#call pargi (rtable)
+		#call pargi (wtable)
+	#}
+
+	# Check for INDEF values in reference table.
+	call ft_indefr (rtable, wtable)
+
+	# Check for INDEF values in fitting equation.
+	call ft_indefo (sym, otable, wtable)
+end
+
+
+# FT_INDEFR - Check reference table for INDEF values. If an INDEF value is
+# found, its corresponding weight (in the weight table) is set to zero, and
+# the INDEF value (in the refence table) replaced by a more suitable one.
+# The latter is because the INLFIT package does not handle INDEF values at
+# all, and it's better to feed it with reasonable values to avoid an overflow
+# or underflow condition.
+
+procedure ft_indefr (rtable, wtable)
+
+pointer	rtable			# reference table
+pointer	wtable			# weight table (modified)
+
+int	npts			# number of points
+int	n
+real	rval, replace
+
+#bool	clgetb()
+int	mct_nrows()
+real	mct_getr()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_indefr: (rtable=%d) (wtable=%d)\n")
+		#call pargi (rtable)
+		#call pargi (wtable)
+	#}
+
+	# Get the number of points.
+	npts = mct_nrows (rtable)
+
+	# Initialize replace value to first non-INDEF value,
+	# if any. Otherwise set it t zero.
+	replace = 0.0
+	do n = 1, npts {
+	    rval = mct_getr (rtable, n, 1)
+	    if (!IS_INDEFR (rval)) {
+		replace = rval
+		break
+	    }
+	}
+
+	# Loop over all data in the table.
+	do n = 1, npts {
+
+	    # Replace values if is INDEF. Otherwise just
+	    # update the replace value.
+	    rval = mct_getr (rtable, n, 1)
+	    if (IS_INDEFR (rval)) {
+		call mct_putr (wtable, n, 1, 0.0)
+		call mct_putr (rtable, n, 1, replace)
+	    } else
+		replace = rval
+	}
+
+	# Debug ?
+	#call dg_dweights ("from ft_indefr", wtable)
+	#call dg_dref ("from ft_indefr", rtable)
+end
+
+
+# FT_INDEFO - Check fitting equation for INDEF values. If an INDEF value is
+# found, its corresponding weight (in the weight table) is set to zero.
+# Undefined values in the table are set to more suitable values, so there
+# won't be problems when plotting data.
+
+procedure ft_indefo (sym, otable, wtable)
+
+int	sym			# equation symbol
+pointer	otable			# observation table (modified)
+pointer	wtable			# weight table (modified)
+
+int	i, n
+int	npts			# number of points
+int	nvars			# number of variables
+real	rval
+pointer	code			# fitting equation code
+pointer	parval			# parameter values
+pointer	replace			# replace values
+pointer	sp
+
+#bool	clgetb()
+int	mct_nrows(), mct_maxcol()
+real	mct_getr()
+real	pr_eval()
+pointer	mct_getrow()
+pointer	pr_gsymp()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_indef: (sym=%d) (otable=%d) (wtable=%d)\n")
+		#call pargi (sym)
+		#call pargi (otable)
+		#call pargi (wtable)
+	#}
+
+	# Get the number of variables and points.
+	npts = mct_nrows (otable)
+	nvars = mct_maxcol (otable)
+
+	# Allocate space for replace values.
+	call smark (sp)
+	call salloc (replace, nvars, TY_REAL)
+
+	# Initialize replace values to first non-undefined
+	# value, if any. Otherwise set it t zero.
+	call aclrr (Memr[replace], nvars)
+	do i = 1, nvars {
+	    do n = 1, npts {
+		rval = mct_getr (otable, n, i)
+		if (!IS_INDEFR (rval)) {
+		    Memr[replace + i - 1] = rval
+		    break
+		}
+	    }
+	}
+
+	# Get the parameter values, and equation code.
+	parval = pr_gsymp (sym, PTEQSPARVAL)
+	code = pr_gsymp (sym, PTEQRPNFIT)
+
+	# Iterate over all the observations.
+	do n = 1, npts {
+
+	    # Evaluate fitting equation.
+	    rval = pr_eval (code, Memr[mct_getrow (otable, n)], Memr[parval])
+
+	    # Substitute weight.
+	    if (IS_INDEFR (rval))
+		call mct_putr (wtable, n, 1, 0.0)
+
+	    # Substitude undefined variable values.
+	    do i = 1, nvars {
+		rval = mct_getr (otable, n, i)
+		if (IS_INDEFR (rval))
+		    call mct_putr (otable, n, i, Memr[replace + i - 1])
+		else
+		    Memr[replace + i - 1] = rval
+	    }
+	}
+
+	# Free memory.
+	call sfree (sp)
+
+	# Debug ?
+	#call dg_dweights ("from ft_indefo", wtable)
+	#call dg_dcatobs ("from ft_indefo", otable)
+end
diff --git a/noao/digiphot/photcal/fitparams/ftref.x b/noao/digiphot/photcal/fitparams/ftref.x
new file mode 100644
index 00000000..0fab41c4
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/ftref.x
@@ -0,0 +1,50 @@
+include	<error.h>
+include	"../lib/parser.h"
+
+
+# FT_RFEVAL - Evaluate reference equation (left hand side) for all observations
+# and store them into a one column table, called the reference table (rtable).
+
+procedure ft_rfeval (code, otable, rtable)
+
+pointer	code			# equation code
+pointer	otable			# 2d observation table
+pointer	rtable			# 1d reference table (output)
+
+int	n
+real	rval
+real	dummy[1]
+
+#bool	clgetb()
+int	mct_nrows()
+real	pr_eval()
+pointer	mct_getrow()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_rfeval (code=%d) (otable=%d)\n")
+		#call pargi (code)
+		#call pargi (otable)
+	#}
+
+	# Allocate space for reference table.
+	# call mct_alloc (rtable, mct_nrows (otable), 1, TY_REAL)
+
+	# Loop over all data in the table.
+	do n = 1, mct_nrows (otable) {
+
+	    # Evaluate the equation.
+	    iferr (rval = pr_eval (code, Memr[mct_getrow (otable, n)], dummy)) {
+		call eprintf ("ft_ref (%d)\n")
+		    call pargi (n)
+		call erract (EA_ERROR)
+	    }
+
+	    # Put data into reference table.
+	    call mct_putr (rtable, n, 1, rval)
+	}
+
+	# Debug ?
+	#call dg_dref ("from ft_rfeval", rtable)
+end
diff --git a/noao/digiphot/photcal/fitparams/fttrneq.x b/noao/digiphot/photcal/fitparams/fttrneq.x
new file mode 100644
index 00000000..6cfb092a
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/fttrneq.x
@@ -0,0 +1,541 @@
+include	<error.h>
+include	<pkg/xtanswer.h>
+include	<pkg/gtools.h>
+include	<math/nlfit.h>
+include	<pkg/inlfit.h>
+include	"../lib/parser.h"
+include	"../lib/fitparams.h"
+
+# Interactive commands
+
+define	CMD_OPTIONS	"|first|last|prev|next|again|go|quit|"
+define	CMD_FIRST	1
+define	CMD_LAST	2
+define	CMD_PREV	3
+define	CMD_NEXT	4
+define	CMD_AGAIN	5
+define	CMD_GO		6
+define	CMD_QUIT	7
+define	CMD_MIN		CMD_FIRST
+define	CMD_MAX		CMD_QUIT
+
+# Prompt to save fit results into output file and equation structure
+
+define	SAVE_PROMPT	"Do you want to save fit results ?"
+define	QUIT_PROMPT	"Program will quit. Do you really want to do so ?"
+
+# Program labels
+
+define	retry		9999
+
+
+# FT_TRNEQS - Fit all the transformation equations.
+
+procedure ft_trneqs (output, logfile, graphics, otable, ntable, wtflag,
+	addscatter, tol, itmax, interactive, high, low, niterate, grow,
+	log_fit, log_results)
+
+char	output[ARB]		# output file name
+char	logfile[ARB]		# the log file name
+char	graphics[ARB]		# graphics output device
+pointer	otable			# standard observation table
+pointer	ntable			# standard name table
+int	wtflag			# the type of weighting
+int	addscatter		# add a scatter term to the weights
+real	tol			# fit tolerance
+int	itmax			# max number of iterations
+bool	interactive		# interactive fit ?
+real	high, low		# rejection factors
+int	niterate		# number of rejection iterations
+real	grow			# rejection growing radius
+bool	log_fit			# log the fit statistics
+bool	log_results		# log the fit results
+
+int	cmd, neq, maxeq, nobs, sym
+pointer	rtable, wtable, totable, gp, gt
+
+#bool	clgetb()
+int	mct_nrows(), pr_gsym(), pr_geti(), ft_trnanswer()
+pointer	gopen(), gt_init(), pr_gsymp()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_trneqs (output=%s) (graphics=%s) (ot=%d) "
+		#call pargstr (output)
+		#call pargstr (graphics)
+		#call pargi (otable)
+	    #call eprintf ("(tol=%g) (maxiter=%d) (int=%b)\n")
+		#call pargr (tol)
+		#call pargi (itmax)
+		#call pargb (interactive)
+	    #call eprintf ("ft_trneqs (gp=%d) (gt=%d) (high=%g) (lo=%g) ")
+		#call pargi (gp)
+		#call pargi (gt)
+		#call pargr (high)
+		#call pargr (low)
+	    #call eprintf ("(niter=%d) (grow=%g)\n")
+		#call pargi (niterate)
+		#call pargr (grow)
+	#}
+
+	# Get number of observations. 
+	nobs = mct_nrows (otable)
+
+	# Allocate space for reference and weight tables. These tables
+	# store the values of the reference equation for all observations,
+	# and the weight for each of them, respectively.
+
+	totable = NULL
+	call mct_alloc (rtable, nobs, 1, TY_REAL)
+	call mct_alloc (wtable, nobs, 1, TY_REAL)
+
+	# Open graphics.
+
+	if (interactive) {
+	    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+	    gt = gt_init()
+	} else {
+	    gp = NULL
+	    gt = NULL
+	}
+
+	# Loop over transformation equations, until the user
+	# sends a quit command.
+
+	neq = 1
+	cmd = CMD_NEXT
+	maxeq = pr_geti (NTRNEQS)
+	repeat {
+
+	    # Get the equation symbol.
+
+	    sym = pr_gsym (neq, PTY_TRNEQ)
+
+	    # Evaluate the reference equation for catalog
+	    # observations and store them into a table.
+
+	    call ft_rfeval (pr_gsymp (sym, PTEQRPNREF), otable, rtable)
+
+	    # Evaluate weights for the current equation. If an error
+	    # condition is raised (due to a negative weight value) the
+	    # program continues with the next equation in the list.
+
+	    switch (wtflag) {
+	    case FWT_UNIFORM:
+		call ft_wtuniform (otable, wtable)
+
+	    case FWT_PHOTOMETRIC:
+		iferr (call ft_wtphoterrs (sym, otable, wtable)) {
+		    call erract (EA_WARN)
+		    next
+		}
+	    case FWT_EQUATIONS:
+	        iferr (call ft_wteqns (sym, otable, wtable)) {
+		    call erract (EA_WARN)
+		    next
+	        }
+	    default:
+		call ft_wtuniform (otable, wtable)
+	    }
+
+	    # Save original observations, and check for undefined values.
+
+	    call mct_copy (otable, totable)
+	    call ft_indef (sym, totable, rtable, wtable)
+
+	    # Fit transformation equation. 
+
+	    call ft_trneq (output, logfile, sym, totable, wtable, rtable,
+	        ntable, wtflag, addscatter, tol, itmax, interactive, gp, gt,
+		high, low, niterate, grow, log_fit, log_results)
+
+	    # Reset the reference table.
+
+	    call mct_reset (rtable)
+
+retry
+	    # Prompt the user for command.
+	    if (interactive)
+		call ft_trncmd (cmd)
+	    else
+		cmd = CMD_NEXT
+
+	    # Branch on command.
+	    switch (cmd) {
+	    case CMD_FIRST:
+		neq = 1
+	    case CMD_LAST:
+		neq = maxeq
+	    case CMD_PREV:
+		if (neq > 1)
+		    neq = neq - 1
+		else
+		    goto retry
+	    case CMD_NEXT:
+		neq = neq + 1
+	    case CMD_AGAIN:
+		;
+	    case CMD_GO:
+		interactive = false
+		neq = neq + 1
+	    case CMD_QUIT:
+		if (ft_trnanswer (QUIT_PROMPT, YES) == YES)
+		    break
+	    default:
+		call error (0, "fttrneqs: Unknown command")
+	    }
+
+	    # Prompt to quit.
+	    if (neq > maxeq) {
+		if (interactive && cmd != CMD_GO) {
+		    if (ft_trnanswer (QUIT_PROMPT, YES) == NO) {
+			neq = maxeq
+			goto retry
+		    } else 
+			break
+		} else
+		    break
+	    }
+	}
+
+	# Close graphics.
+	if (gp != NULL)
+	    call gclose (gp)
+	if (gt != NULL)
+	    call gt_free (gt)
+
+	# Free tables.
+	call mct_free (rtable)
+	call mct_free (wtable)
+	call mct_free (totable)
+end
+
+
+# FT_TRNEQ - Fit single transformation equation.
+
+procedure ft_trneq (output, logfile, sym, otable, wtable, rtable, ntable,
+	wtflag, addscatter, tol, itmax, interactive, gp, gt high, low,
+	niterate, grow, log_fit, log_results)
+
+char	output[ARB]		# output file name
+char	logfile[ARB]		# log file name
+int	sym			# transformation equation symbol
+pointer	otable			# standard observation table
+pointer	wtable			# weight table
+pointer	rtable			# reference equation table
+pointer	ntable			# names table
+int	wtflag			# type of weighting
+int	addscatter		# add a scatter term to the weight equation
+real	tol			# fit tolerance
+int	itmax			# max number of iterations
+bool	interactive		# interactive fit ?
+pointer	gp, gt			# GIO and GTOOLS descriptors
+real	high, low		# rejection factors
+int	niterate		# number of rejection iterations
+real	grow			# rejection growing radius
+bool	log_fit			# log the fit statistics
+bool	log_results		# log the fit results
+
+int	nobs, nvars, nparams, nfparams, len_name
+int	i, stat, answer, psym, nlwtflag
+pointer	sp, params, dparams, name, aux, plist, nl, in
+real	chisqr, variance, scatter, rms
+
+#bool	clgetb()
+bool	streq()
+int	locpr(), mct_nrows(), mct_maxcol(), pr_geti(), pr_gsym(), pr_gsymi()
+int	pr_gpari(), ft_trnanswer()
+pointer	pr_xgetname(), pr_gsymp(), pr_gsymc(), mct_getbuf(), in_getp()
+real	pr_gsymr()
+extern	ft_func(), ft_dfunc(), ft_plot()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode")) {
+	    #call eprintf ("ft_trneq (output=%s) (sym=%d) (ot=%d) (wt=%d) "
+		#call pargstr (output)
+		#call pargi (sym)
+		#call pargi (otable)
+		#call pargi (wtable)
+	    #call eprintf ("(rt=%d) (tol=%g) (maxiter=%d)\n")
+		#call pargi (rtable)
+		#call pargr (tol)
+		#call pargi (itmax)
+	    #call eprintf ("ft_trneq (int=%b) (gp=%d) (gt=%d) (high=%g) "
+		#call pargb (interactive)
+		#call pargi (gp)
+		#call pargi (gt)
+		#call pargr (high)
+	    #call eprintf ("(low=%g) (niter=%d) (grow=%g)\n")
+		#call pargr (low)
+		#call pargi (niterate)
+		#call pargr (grow)
+	#}
+
+	# Get number of observations, and variables for the observational data.
+
+	nobs = mct_nrows (otable)
+	nvars = mct_maxcol (otable)
+	len_name = mct_maxcol (ntable)
+
+	# Get number of parameters for the current equation.
+
+	nparams = pr_gsymi (sym, PTEQNPAR)
+
+	# Allocate stack space.
+	call smark (sp)
+	call salloc (aux, SZ_LINE, TY_CHAR)
+	call salloc (name, SZ_LINE, TY_CHAR)
+	call salloc (params, nparams, TY_REAL)
+	call salloc (dparams, nparams, TY_REAL)
+	call salloc (plist, nparams, TY_INT)
+
+	# Initialize parameters, errors, and parameter list.
+
+	call aclrr (Memr[params], nparams)
+	call aclrr (Memr[dparams], nparams)
+	do i = 1, nparams {
+	    psym = pr_gpari (sym, i, PTEQPAR)
+	    Memr[dparams+i-1] = pr_gsymr (psym, PFITDELTA)
+	}
+	call amovi (Memi[pr_gsymp (sym, PTEQSPLIST)], Memi[plist], nparams)
+
+	# Initialize the fit evaluation. This is necessary in order to
+	# set the equations called by the INLFIT procedures.
+
+	call ft_evinit (sym, nparams)
+
+	# Get number of fitting parameters for the current equation.
+	nfparams = pr_gsymi (sym, PTEQNFPAR)
+
+	# Initialize INLFIT.
+	call in_initr (in, locpr (ft_func), locpr (ft_dfunc),
+	    Memr[pr_gsymp (sym, PTEQSPARVAL)], Memr[dparams],
+	    nparams, Memi[plist], nfparams)
+
+	# Set INLFIT fitting parameters.
+	call in_putr (in, INLTOLERANCE, tol)
+	call in_puti (in, INLMAXITER, itmax)
+	call in_putr (in, INLHIGH, high)
+	call in_putr (in, INLLOW, low)
+	call in_puti (in, INLNREJECT, niterate)
+	call in_putr (in, INLGROW, grow)
+
+	# Put in the reference and fit equation names.
+	call sprintf (Memc[aux], SZ_LINE, "|%s|%s|")
+	    call pargstr (Memc[pr_gsymc (sym, PTEQREF)])
+	    call pargstr (Memc[pr_gsymc (sym, PTEQFIT)])
+	call in_pstr (in, INLFLABELS, Memc[aux])
+
+	# Put in the parameter names.
+	call strcpy ("|", Memc[aux], SZ_LINE)
+	do i = 1, nparams {
+	    call strcat (Memc[pr_xgetname (pr_gpari (sym, i, PTEQPAR))],
+			 Memc[aux], SZ_LINE)
+	    call strcat ("|", Memc[aux], SZ_LINE)
+	}
+	call in_pstr (in, INLPLABELS, Memc[aux])
+
+	# Put in the variable names.
+	call strcpy ("|", Memc[aux], SZ_LINE)
+	do i = 1, pr_geti (NOBSVARS) {
+	    call pr_vtran (Memc[pr_xgetname (pr_gsym (i, PTY_OBSVAR))],
+		Memc[name], SZ_LINE)
+	    call strcat (Memc[name], Memc[aux], SZ_LINE)
+	    call strcat ("|", Memc[aux], SZ_LINE)
+	}
+	do i = 1, pr_geti (NCATVARS) {
+	    call pr_vtran (Memc[pr_xgetname (pr_gsym (i, PTY_CATVAR))],
+		Memc[name], SZ_LINE)
+	    call strcat (Memc[name], Memc[aux], SZ_LINE)
+	    call strcat ("|", Memc[aux], SZ_LINE)
+	}
+	call in_pstr (in, INLVLABELS, Memc[aux])
+
+	# Put plot equations and redefine graph keys, but only if
+	# both plot equations are defined. Otherwise leave the defaults.
+
+	if (pr_gsymp (sym, PTEQRPNXPLOT) != NULL &&
+	    pr_gsymp (sym, PTEQRPNYPLOT) != NULL) {
+
+	    # Put in the plot equation names.
+	    call sprintf (Memc[aux], SZ_LINE, "|%s|%s|")
+		call pargstr (Memc[pr_gsymc (sym, PTEQXPLOT)])
+		call pargstr (Memc[pr_gsymc (sym, PTEQYPLOT)])
+	    call in_pstr (in, INLUSERLABELS, Memc[aux])
+
+	    # Put in the plot equation.
+	    call in_puti (in, INLUAXES, locpr (ft_plot))
+
+	    # Redefine the last graph key.
+	    call in_puti (in, INLGKEY, INLNGKEYS)
+	    call in_pkey (in, INLNGKEYS, 1, KEY_UAXIS, 1)
+	    call in_pkey (in, INLNGKEYS, 2, KEY_UAXIS, 2)
+	}
+
+	# Call the appropiate version of fitting routine according to the
+	# value of the interactive flag.
+
+	if (wtflag == FWT_UNIFORM)
+	    nlwtflag = WTS_USER
+	else if (addscatter == YES)
+	    nlwtflag = WTS_SCATTER
+	else
+	    nlwtflag = WTS_USER
+
+	if (interactive) {
+	    call gt_setr (gt, GTXMIN, INDEFR)
+	    call gt_setr (gt, GTXMAX, INDEFR)
+	    call gt_setr (gt, GTYMIN, INDEFR)
+	    call gt_setr (gt, GTYMAX, INDEFR)
+	    call ing_fitr (in, gp, "cursor", gt, nl, Memr[mct_getbuf (otable)],
+	        Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)],
+		Memc[mct_getbuf(ntable)], nobs, nvars, len_name, nlwtflag,
+		stat)
+	} else {
+	    call in_fitr (in, nl, Memr[mct_getbuf (otable)],
+	        Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)],
+		nobs, nvars, nlwtflag, stat)
+	}
+
+	# Decide whether to prompt the user for saving answer.
+	if (interactive)
+	    answer = ft_trnanswer (SAVE_PROMPT, YES)
+	else
+	    answer = YES
+
+	# Write the parameter values into output file.
+	if (answer == YES || answer == ALWAYSYES) {
+
+	    # Get parameter values, and parameter list.
+	    call nlpgetr (nl, Memr[params], nparams)
+	    call amovi (Memi[in_getp (in, INLPLIST)], Memi[plist], nparams)
+
+	    # Get the parameter errors.
+	    call in_errorsr (in, nl, Memr[mct_getbuf (otable)],
+	        Memr[mct_getbuf (rtable)], Memr[mct_getbuf (wtable)], nobs,
+		nvars, variance, chisqr, scatter, rms, Memr[dparams])
+
+	    # Write parameters and errors into output file.
+	    iferr (call io_pcoeffs (output, sym, stat, wtflag, variance, chisqr,
+	        scatter, rms, Memr[params], Memr[dparams], Memi[plist],
+		nparams))
+	        call erract (EA_WARN)
+
+	    # Log the fit and results.
+	    if (log_fit || log_results) {
+		if (interactive && streq (logfile, "STDOUT"))
+		    call gdeactivate (gp, 0)
+		if (log_fit) {
+		    iferr (call io_title (logfile, "#EQUATION:", sym))
+		        ;
+		    iferr (call ing_showr (in, logfile))
+			;
+		    iferr (call ing_errorsr (in, logfile, nl,
+		        Memr[mct_getbuf(otable)], Memr[mct_getbuf(rtable)],
+		        Memr[mct_getbuf(wtable)], nobs, nvars))
+			;
+		}
+		if (log_results) {
+		    iferr (call io_title (logfile, "#RESULTS:", sym))
+		        ;
+		    iferr (call ing_resultsr (in, logfile, nl,
+		        Memr[mct_getbuf(otable)], Memr[mct_getbuf(rtable)],
+		        Memr[mct_getbuf(wtable)], Memc[mct_getbuf(ntable)],
+		        nobs, nvars, len_name))
+			    ;
+		}
+		if (interactive && streq (logfile, "STDOUT"))
+		    call greactivate (gp, 0)
+	    }
+
+
+	    # Update fitted parameter values into the equation symbol
+	    # substructure if the fit was succesfull.
+	    if (stat == DONE)
+		call amovr (Memr[params], 
+		    Memr[P2R(pr_gsymp (sym, PTEQSPARVAL))], nparams)
+	}
+
+	# Debug ?
+	#if (clgetb ("debug.nlfit"))
+	    #call dg_inldump (in, nl)
+
+
+	# Free inlfit and nlfit descriptors.
+	call in_freer (in)
+	call nlfreer (nl)
+
+	# Free fit evaluation.
+	call ft_evfree ()
+
+	# Free stack space.
+	call sfree (sp)
+end
+
+
+# FT_TRNCMD -- Prompt the user for a command.
+
+procedure ft_trncmd (cmd)
+
+int	cmd			# command code
+
+char	command[SZ_LINE]
+int	scan(), strlen()
+int	strdic(), io_strwrd()
+
+begin
+	# Get current command string.
+	if (io_strwrd (cmd, command, SZ_LINE, CMD_OPTIONS) == 0)
+	    call error (0, "ft_trneqs: Unknown command")
+
+	# Keep prompting the user until there is a valid command.
+	repeat {
+
+	    # Print current setting.
+	    call printf ("Command (%s) (%s) : ")
+		call pargstr (CMD_OPTIONS)
+		call pargstr (command)
+	    call flush (STDOUT)
+
+	    # Get new command.
+	    if (scan () == EOF)
+		command[1] = EOS
+	    else
+	        call gargstr (command, SZ_LINE)
+	    if (strlen (command) == 0) {
+		if (io_strwrd (cmd, command, SZ_LINE, CMD_OPTIONS) <= 0)
+		    ;
+	    }
+
+	    # Search for command in the dictionary.
+	    cmd = strdic (command, command, SZ_LINE, CMD_OPTIONS)
+
+	    # Test command.
+	    if (cmd >= CMD_MIN && cmd <= CMD_MAX)
+		break
+	    else
+		call printf ("\007")
+	}
+end
+
+
+# FT_TRNANSWER -- Get a YES/NO answer from the user, and return it
+# as the procedure value.
+
+int procedure ft_trnanswer (prompt, default)
+
+char	prompt			# prompt to the user
+int	default			# default answer
+
+int	answer
+
+begin
+	answer = default
+	call xt_answer (prompt, answer)
+	if (answer == YES || answer == ALWAYSYES)
+	    return (YES)
+	else
+	    return (NO)
+end
diff --git a/noao/digiphot/photcal/fitparams/ftweights.x b/noao/digiphot/photcal/fitparams/ftweights.x
new file mode 100644
index 00000000..67482867
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/ftweights.x
@@ -0,0 +1,491 @@
+include	"../lib/parser.h"
+include "../lib/preval.h"
+
+# FT_WTUNIFORM -- Set all the weights to 1.0.
+
+procedure ft_wtuniform (otable, wtable)
+
+pointer	otable				# observation table
+pointer	wtable				# weight table (output)
+
+int	nobs
+int	mct_nrows()
+
+begin
+	# Get number of observations.
+	nobs = mct_nrows (otable)
+
+	# Clear weight table with ones (uniform weighting), and enter
+	# data at the last observation set weight table counters.
+
+	call mct_clearr (wtable, 1.0)
+	call mct_putr (wtable, nobs, 1, 1.0)
+end
+
+
+# FT_WTEQNS -- Compute the value of the weights equation for each point. 
+
+procedure ft_wteqns (sym, otable, wtable)
+
+int	sym				# equation symbol
+pointer	otable				# observation table
+pointer	wtable				# weight table (output)
+
+int	n, nobs
+real	wtval, minval, maxval
+pointer	varptr, parptr
+pointer	wtcode, mincode, maxcode
+
+#bool	clgetb()
+int	mct_nrows()
+pointer	pr_gsymp(), mct_getrow()
+real	pr_eval()
+
+begin
+	# Debug ?
+	#if (clgetb ("debug.fitcode"))
+	    #call eprintf ("ft_wteval (sym=%d) (ot=%d) (wt=%d)\n") {
+		#call pargi (sym)
+		#call pargi (otable)
+		#call pargi (wtable)
+	#}
+
+	# Get number of observations.
+	nobs = mct_nrows (otable)
+
+	# Clear weight table with ones (uniform weighting), and enter
+	# data at the last observation set weight table counters.
+
+	call mct_clearr (wtable, 1.0)
+	call mct_putr (wtable, nobs, 1, 1.0)
+
+	# Compute weigths using the weight equation, if there is one defined.
+	wtcode = pr_gsymp (sym, PTEQRPNWEIGHT)
+	if (wtcode != NULL) {
+
+	    # Get pointer to parameter values for the current equation.
+	    parptr = pr_gsymp (sym, PTEQSPARVAL)
+
+	    # Get minimum and maximum equation codes.
+	    mincode = pr_gsymp (sym, PTEQRPNWTSMIN)
+	    maxcode = pr_gsymp (sym, PTEQRPNWTSMAX)
+
+	    # Iterate over all observations.
+	    do n = 1, nobs {
+
+		# Get variable values for the current observation.
+		varptr = mct_getrow (otable, n)
+
+		# Compute weight value.
+		wtval = pr_eval (wtcode, Memr[varptr], Memr[parptr])
+		if (IS_INDEFR(wtval) || wtval < 0.0) {
+
+		    call mct_putr (wtable, n, 1, 0.0)
+
+		} else {
+
+		    # Check value against minimum.
+		    if (mincode != NULL) {
+		        minval = pr_eval (mincode, Memr[varptr], Memr[parptr])
+		        if (! IS_INDEFR(minval) && minval >= 0.0)
+		            wtval = min (wtval, minval)
+		    }
+
+		    # Check value against maximum.
+		    if (maxcode != NULL) {
+		        maxval = pr_eval (maxcode, Memr[varptr], Memr[parptr])
+		        if (! IS_INDEFR(maxval) && maxval >= 0.0)
+		            wtval = max (wtval, maxval)
+		    }
+
+		    # Enter value into weight table.
+		    call mct_putr (wtable, n, 1, wtval)
+		}
+	    }
+	}
+
+	# Debug ?
+	#call dg_dweigths ("from ft_wteval", wtable)
+end
+
+
+# FT_WTPHOTERRS -- Set all the weights to the correct statistical value assuming
+# that all the errors are in the photometric indices, that the errors are
+# independent, that the equations are linear in the photometric indices,
+# and that at least one of the observed or catalog variables in the equation
+# has a measured error.
+
+procedure ft_wtphoterrs (sym, otable, wtable)
+
+int	sym				# equation symbol
+pointer	otable				# observation table
+pointer	wtable				# weight table (output)
+
+double	errval
+int	n, i, nobs, nobsvars, ncatvars, nrvars, nfvars, nerrors, symvar, icol
+int	itemp
+pointer	omap, cmap, sp, rcount, rerrcol, fcount, ferrcol, refcode
+real	val
+int	mct_nrows(), pr_geti(), pr_gsymi(), pr_gvari(), pr_findmap1()
+int	ft_seval()
+pointer	pr_gsymp()
+#int	mct_ncols()
+real	mct_getr()
+errchk	pr_gsymi()
+
+begin
+	# Get number of observations.
+	nobs = mct_nrows (otable)
+
+	# Clear weight table with ones (uniform weighting), and enter
+	# data at the last observation set weight table counters.
+
+	call mct_clearr (wtable, 1.0)
+	call mct_putr (wtable, nobs, 1, 1.0)
+
+	# Map the column numbers for the observed and catalog variables.
+	call pr_obsmap (omap, nobsvars)
+	call pr_catmap (cmap, ncatvars)
+
+	# Get the number of reference and fit equation variables.
+	nrvars = pr_gsymi (sym, PTEQNRCAT) + pr_gsymi (sym, PTEQNROBS) 
+	nfvars = pr_gsymi (sym, PTEQNFCAT) + pr_gsymi (sym, PTEQNFOBS) 
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (rerrcol, nobsvars + ncatvars, TY_INT)
+	call amovi (NULL, Memi[rerrcol], nobsvars + ncatvars)
+	call salloc (rcount, nobsvars + ncatvars, TY_INT)
+	call aclri (Memi[rcount], nobsvars + ncatvars)
+	call salloc (ferrcol, nobsvars + ncatvars, TY_INT)
+	call amovi (NULL, Memi[ferrcol], nobsvars + ncatvars)
+	call salloc (fcount, nobsvars + ncatvars, TY_INT)
+	call aclri (Memi[fcount], nobsvars + ncatvars)
+
+	# Initialize.
+	nerrors = 0
+
+	# Compute the positions of the reference equation variable errors in
+	# the observations table.
+
+	do i = 1, nrvars {
+	    Memi[rcount+i-1] = pr_gvari (sym, i, PTEQREFCNT)
+	    symvar = pr_gvari (sym, i, PTEQREFVAR)
+	    icol = pr_gsymi (symvar, PINPERRCOL)
+	    if (IS_INDEFI(icol))
+		Memi[rerrcol+i-1] = NULL
+	    else {
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR)
+		    Memi[rerrcol+i-1] = pr_findmap1 (cmap, icol) + nobsvars
+		else
+		    Memi[rerrcol+i-1] = pr_findmap1 (omap, icol)
+		nerrors = nerrors + 1
+	    }
+	}
+
+	# Evaluate the contribution of the set equations to the errors
+	# in the reference equation.
+
+	if (pr_geti (NSETEQS) > 0) {
+	    refcode = pr_gsymp (sym, PTEQRPNREF)
+	    itemp = ft_seval (refcode, cmap, omap, nobsvars,
+	        Memi[rerrcol], Memi[rcount], nrvars)
+	    nrvars = nrvars + itemp 
+	    nerrors = nerrors + itemp
+	}
+
+	# Compute the positions of the fit equation variable errors in the
+	# observations table.
+
+	do i = 1, nfvars {
+	    Memi[fcount+i-1] = pr_gvari (sym, i, PTEQFITCNT)
+	    symvar = pr_gvari (sym, i, PTEQFITVAR)
+	    icol = pr_gsymi (symvar, PINPERRCOL)
+	    if (IS_INDEFI(icol))
+		Memi[ferrcol+i-1] = NULL
+	    else {
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR)
+		    Memi[ferrcol+i-1] = pr_findmap1 (cmap, icol) + nobsvars
+		else
+		    Memi[ferrcol+i-1] = pr_findmap1 (omap, icol)
+		nerrors = nerrors + 1
+	    }
+	}
+
+	# Evaluate the contribution of the set equations to the errors
+	# in the fit equation.
+
+	if (pr_geti (NSETEQS) > 0) {
+	    refcode = pr_gsymp (sym, PTEQRPNFIT)
+	    itemp = ft_seval (refcode, cmap, omap, nobsvars,
+	        Memi[ferrcol], Memi[fcount], nfvars)
+	    nfvars = nfvars + itemp
+	    nerrors = nerrors + itemp
+	}
+
+	# Loop over the table rows.
+	if (nerrors > 0) {
+	    do n = 1, nobs {
+
+	        errval = 0.0d0
+
+	        # Add contributions from the reference equation variables.
+	        do i = 1, nrvars {
+		    if (Memi[rerrcol+i-1] == NULL)
+		        next
+		    val = mct_getr (otable, n, Memi[rerrcol+i-1])
+		    if (IS_INDEFR(val) || val < 0.0)
+			next
+		    errval = errval + real (Memi[rcount+i-1]) * (val * val)
+	        }
+
+	        # Add contributions from the fitting equation variables.
+	        do i = 1, nfvars {
+		    if (Memi[ferrcol+i-1] == NULL)
+		        next
+		    val = mct_getr (otable, n, Memi[ferrcol+i-1])
+		    if (IS_INDEFR(val) || val < 0.0)
+			next
+		    errval = errval + real (Memi[fcount+i-1]) * (val * val)
+	        }
+
+	        # Check for negative and zero error values and enter value
+		# into weight table.
+	        if (errval <= 0.0d0)
+	            call mct_putr (wtable, n, 1, 0.0)
+		else
+	            call mct_putr (wtable, n, 1, real (1.0d0 / errval))
+	    }
+	}
+
+	call pr_unmap (cmap)
+	call pr_unmap (omap)
+
+	call sfree (sp)
+end
+
+
+# FT_SEVAL -- Locate any catalog and observations variables that are
+# defined in up to two levels of set equations.
+
+int procedure ft_seval (code, cmap, omap, nobsvars, errcol, ecount, nerrors)
+
+pointer	code		# the equation code
+pointer	cmap		# pointer to the catalog variable column map
+pointer	omap		# pointer to the observations variable column map
+int	nobsvars	# the number of observations variables
+int	errcol[ARB]	# the input/output error columns
+int	ecount[ARB]	# the input/output error column counts
+int	nerrors		# the previous number of errors
+
+int	ip, ins, newerr
+pointer	setcode
+int	pr_gsym(), ft_scoeval()
+pointer	pr_gsymp()
+
+begin
+	ip = 0
+	ins = Memi[code+ip]
+
+	newerr = 0
+	while (ins != PEV_EOC) {
+	    switch (ins) {
+	    case PEV_SETEQ:
+		ip = ip + 1
+		setcode = pr_gsymp (pr_gsym (Memi[code+ip], PTY_SETEQ),
+		    PSEQRPNEQ)
+		newerr = newerr + ft_scoeval (setcode, cmap, omap, nobsvars,
+		    errcol, ecount, nerrors + newerr)
+	    case PEV_NUMBER, PEV_CATVAR, PEV_OBSVAR, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		;
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[code+ip]
+	}
+
+	return (newerr)
+end
+
+
+# FT_SCOEVAL -- Locate any catalog and observations variables as well as any
+# set equations that are defined in a previous set equation.
+
+int procedure ft_scoeval (code, cmap, omap, nobsvars, errcol, ecount, nerrors)
+
+pointer	code		# the equation code
+pointer	cmap		# pointer to the catalog variable column map
+pointer	omap		# pointer to the observations variable column map
+int	nobsvars	# the number of observations variables
+int	errcol[ARB]	# the input/output error columns
+int	ecount[ARB]	# the input/output error column counts
+int	nerrors		# the previous number of errors
+
+bool	new
+int	i, ip, ins, newerr, ecol
+pointer	sym, setcode
+int	pr_gsym(), pr_gsymi(), pr_findmap1(), ft_coeval()
+pointer	pr_gsymp()
+
+begin
+	ip = 0
+	ins = Memi[code+ip]
+
+	newerr = 0
+	while (ins != PEV_EOC) {
+	    switch (ins) {
+	    case PEV_SETEQ:
+		ip = ip + 1
+		setcode = pr_gsymp (pr_gsym (Memi[code+ip], PTY_SETEQ),
+		    PSEQRPNEQ)
+		newerr = newerr + ft_coeval (setcode, cmap, omap, nobsvars,
+		    errcol, ecount, nerrors + newerr)
+	    case PEV_CATVAR:
+		ip = ip + 1
+		sym = pr_gsym (Memi[code+ip] - nobsvars, PTY_CATVAR)
+		ecol = pr_gsymi (sym, PINPERRCOL)
+		if (! IS_INDEFI(ecol))
+		    ecol = pr_findmap1 (cmap, ecol) + nobsvars
+		else
+		    ecol = NULL
+		if (ecol != NULL) {
+		    new = true
+		    do i = 1, nerrors + newerr {
+			if (ecol != errcol[i])
+			    next
+			new = false
+			ecount[i] = ecount[i] + 1 
+		    }
+		    if (new) {
+			newerr = newerr + 1
+			errcol[nerrors+newerr] = ecol
+			ecount[nerrors+newerr] = 1 
+		    }
+		}
+	    case PEV_OBSVAR:
+		ip = ip + 1
+		sym = pr_gsym (Memi[code+ip], PTY_OBSVAR)
+		ecol = pr_gsymi (sym, PINPERRCOL)
+		if (! IS_INDEFI(ecol))
+		    ecol = pr_findmap1 (omap, ecol)
+		else
+		    ecol = NULL
+		if (ecol != NULL) {
+		    new = true
+		    do i = 1, nerrors + newerr {
+			if (ecol != errcol[i])
+			    next
+			new = false
+			ecount[i] = ecount[i] + 1 
+		    }
+		    if (new) {
+			newerr = newerr + 1
+			errcol[nerrors+newerr] = ecol
+			ecount[nerrors+newerr] = 1 
+		    }
+		}
+	    case PEV_NUMBER, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		;
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[code+ip]
+	}
+
+	return (newerr)
+end
+
+
+# FT_COEVAL -- Locate any catalog and observations variables that are defined i
+# a previous set equation.
+
+int procedure ft_coeval (code, cmap, omap, nobsvars, errcol, ecount, nerrors)
+
+pointer	code		# the equation code
+pointer	cmap		# pointer to the catalog variable column map
+pointer	omap		# pointer to the observations variable column map
+int	nobsvars	# the number of observations variables
+int	errcol[ARB]	# the input/output error columns
+int	ecount[ARB]	# the input/output error column counts
+int	nerrors		# the previous number of errors
+
+bool	new
+int	i, ip, ins, newerr, ecol
+pointer	sym
+int	pr_gsym(), pr_gsymi(), pr_findmap1()
+
+begin
+	ip = 0
+	ins = Memi[code+ip]
+
+	newerr = 0
+	while (ins != PEV_EOC) {
+	    switch (ins) {
+	    case PEV_SETEQ:
+		ip = ip + 1
+	    case PEV_CATVAR:
+		ip = ip + 1
+		sym = pr_gsym (Memi[code+ip] - nobsvars, PTY_CATVAR)
+		ecol = pr_gsymi (sym, PINPERRCOL)
+		if (! IS_INDEFI(ecol))
+		    ecol = pr_findmap1 (cmap, ecol) + nobsvars
+		else
+		    ecol = NULL
+		if (ecol != NULL) {
+		    new = true
+		    do i = 1, nerrors + newerr {
+			if (ecol != errcol[i])
+			    next
+			new = false
+			ecount[i] = ecount[i] + 1 
+		    }
+		    if (new) {
+			newerr = newerr + 1
+			errcol[nerrors+newerr] = ecol
+			ecount[nerrors+newerr] = 1 
+		    }
+		}
+	    case PEV_OBSVAR:
+		ip = ip + 1
+		sym = pr_gsym (Memi[code+ip], PTY_OBSVAR)
+		ecol = pr_gsymi (sym, PINPERRCOL)
+		if (! IS_INDEFI(ecol))
+		    ecol = pr_findmap1 (omap, ecol)
+		else
+		    ecol = NULL
+		if (ecol != NULL) {
+		    new = true
+		    do i = 1, nerrors + newerr {
+			if (ecol != errcol[i])
+			    next
+			new = false
+			ecount[i] = ecount[i] + 1 
+		    }
+		    if (new) {
+			newerr = newerr + 1
+			errcol[nerrors+newerr] = ecol
+			ecount[nerrors+newerr] = 1 
+		    }
+		}
+	    case PEV_NUMBER, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		;
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[code+ip]
+	}
+
+	return (newerr)
+end
diff --git a/noao/digiphot/photcal/fitparams/mkpkg b/noao/digiphot/photcal/fitparams/mkpkg
new file mode 100644
index 00000000..c4681932
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/mkpkg
@@ -0,0 +1,16 @@
+# The MKPKG file for the FITPARAMS task.
+
+$checkout	libpkg.a ../
+$update		libpkg.a
+$checkin	libpkg.a ../
+$exit
+
+libpkg.a:
+	fteval.x	"../lib/parser.h" "fteval.com" <error.h>
+	ftindef.x	"../lib/parser.h" "../lib/fitparams.h" <error.h>
+	ftref.x		"../lib/parser.h" <error.h>
+	fttrneq.x	<pkg/inlfit.h> <math/nlfit.h> <error.h> <pkg/gtools.h> \
+			<pkg/xtanswer.h> "../lib/parser.h" "../lib/fitparams.h" 
+	ftweights.x	"../lib/preval.h" "../lib/parser.h"
+	t_fitparams.x	"../lib/parser.h" "../lib/fitparams.h"
+	;
diff --git a/noao/digiphot/photcal/fitparams/t_fitparams.x b/noao/digiphot/photcal/fitparams/t_fitparams.x
new file mode 100644
index 00000000..d306124d
--- /dev/null
+++ b/noao/digiphot/photcal/fitparams/t_fitparams.x
@@ -0,0 +1,205 @@
+include	"../lib/parser.h"
+include	"../lib/fitparams.h"
+
+
+# T_FITPARAMS - Main fitting task. This task will determine the value of the
+# fitting parameters, either interactively or non-interactively, by using
+# the INLFIT package.
+
+procedure t_fitparams()
+
+pointer	observations		# list of observations files
+pointer	catalogs		# list of standard catalogs
+int	stdlist			# file list of standards
+pointer	config			# pointer to configuration file name
+pointer	output			# pointer to output file name
+pointer	logfile			# pointer to the output log file
+int	wtflag			# weighting type
+int	addscatter		# compute an additional scatter term
+real	tol			# fit tolerance
+int	itmax			# max number of iterations
+int	niterate		# number of rejection iterations
+real	high, low		# rejection thresholds
+real	grow			# rejection growing radius
+bool	interactive		# interactive fit ?
+pointer	catdir			# the standard catalogs directory
+pointer graphics		# pointer to the graphics device name
+
+int	obslist, nstd, nobs, nstdvars, getid
+pointer	sp, dir, str, otable, ntable, stable
+
+bool	clgetb()
+int	fntopnb(), fntlenb(), clgeti(), clgwrd(), btoi(), pr_parse(), pr_geti()
+int	fnldir(), access()
+real	clgetr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (observations, SZ_FNAME, TY_CHAR)
+	call salloc (catalogs, SZ_FNAME, TY_CHAR)
+	call salloc (config, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (logfile, SZ_FNAME, TY_CHAR)
+	call salloc (catdir, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+	call salloc (dir, SZ_PATHNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Open the observations files list.
+	call clgstr ("observations", Memc[observations], SZ_FNAME)
+	obslist = fntopnb (Memc[observations], NO)
+	if (Memc[observations] == EOS || fntlenb (obslist) <= 0) {
+	    call eprintf ("ERROR: The observations files list is empty\n")
+	    call fntclsb (obslist)
+	    call sfree (sp)
+	    return
+	}
+
+	# Get the list of catalog files.
+	call clgstr ("catdir", Memc[catdir], SZ_FNAME)
+	call clgstr ("catalogs", Memc[catalogs], SZ_FNAME)
+
+	# Get the configuration file name. 
+	call clgstr ("config", Memc[config], SZ_FNAME)
+	if (access (Memc[config], READ_ONLY, TEXT_FILE) == NO) {
+	    call eprintf (
+	        "ERROR: Cannot open the configuration file for reading\n")
+	    call fntclsb (obslist)
+	    call sfree (sp)
+	    return
+	}
+
+	# Get the output parameters database nam and check that the user
+	# can create or append to the database.
+
+	call clgstr ("parameters", Memc[output], SZ_FNAME)
+	if (Memc[output] == EOS || access (Memc[output], 0,
+	    DIRECTORY_FILE) == YES) {
+	    call eprintf ("ERROR: The parameters file is undefined\n")
+	    call fntclsb (obslist)
+	    call sfree (sp)
+	    return
+	} else if (access (Memc[output], 0, 0) == NO) {
+	    if (fnldir (Memc[output], Memc[dir], SZ_LINE) == 0)
+		call strcpy (".", Memc[dir], SZ_LINE)
+	    if (access (Memc[dir], APPEND, DIRECTORY_FILE) == NO) {
+	    	call eprintf ("ERROR: Cannot open directory %s for writing\n")
+		    call pargstr (Memc[dir])
+	        call fntclsb (obslist)
+	        call sfree (sp)
+		return
+	    }
+	} else if (access (Memc[output], APPEND, TEXT_FILE) == NO) {
+	    call eprintf (
+	        "ERROR: Cannot open existing parameters file %s for writing\n")
+		call pargstr (Memc[output])
+	    call fntclsb (obslist)
+	    call sfree (sp)
+	    return
+	}
+
+	# Get the log file.
+	call clgstr ("logfile", Memc[logfile], SZ_FNAME)
+	if (Memc[logfile] != EOS)
+	    call io_logtime (Memc[logfile])
+
+	# Get weighting parameters.
+	wtflag = clgwrd ("weighting", Memc[str], SZ_LINE, FWT_OPTIONS)
+	addscatter = btoi (clgetb ("addscatter"))
+
+	# Get the fitting parameters.
+	tol = clgetr ("tolerance")
+	itmax = clgeti ("maxiter")
+
+	# Get the rejection parameters.
+	low = clgetr ("low_reject")
+	high = clgetr ("high_reject")
+	niterate = clgeti ("nreject")
+	grow = clgetr ("grow")
+
+	# Get the graphics parameters.
+	interactive = clgetb ("interactive")
+	call clgstr ("graphics", Memc[graphics], SZ_LINE)
+
+	# Parse the configuration table.
+	if (pr_parse (Memc[config]) == ERR) {
+	    call eprintf ("ERROR: Cannot parse the configuration file\n")
+	    call fntclsb (obslist)
+	    call sfree (sp)
+	    return
+	}
+
+	# Read standard data catalog if the catalog section was in the
+	# configuration file, the catalog section is not empty and the
+	# catalogs files list. This can be tested by checking the value
+	# of the minimum input column, the number of catalog variables
+	# and the contents of the standard catalogs list respectively.
+
+	stable = NULL
+	ntable = NULL
+	nstd = 0
+
+	nstdvars = pr_geti (NCATVARS)
+	if (pr_geti (MINCOL) == 1) {
+	    getid = NO
+	} else if (nstdvars == 0) {
+	    getid = YES
+	} else if (Memc[catalogs] == EOS) {
+	    getid = YES
+	    call eprintf ("WARNING: Cannot load catalog variables from ")
+	    call eprintf ("the empty catalog files list\n")
+	} else { 
+	    getid = YES
+	    stdlist = fntopnb (Memc[catalogs], NO)
+	    call io_gcatdat (Memc[catdir], stdlist, stable, nstd, nstdvars)
+	    call fntclsb (stdlist)
+	} 
+
+	# Quit if there is no data.
+	if (stable != NULL && nstd <= 0) {
+	    call eprintf ("ERROR: No data was read from the catalog files")
+	    call stclose (stable)
+	    call fntclsb (obslist)
+	    call pr_free ()
+	    call sfree (sp)
+	    return
+	} 
+
+	# Read in the observations.
+	if (clgetb ("log_unmatched"))
+	    call io_gcatobs (obslist, stable, nstdvars, getid, Memc[logfile],
+	        otable, ntable, nobs)
+	else
+	    call io_gcatobs (obslist, stable, nstdvars, getid, "", otable,
+	        ntable, nobs)
+
+	# Free standard data table since it's no longer needed for the
+	# parameter fitting.
+
+	if (stable != NULL)
+	    call stclose (stable)
+
+	# Process all transformation equations if there are enough observations.
+	if (nobs > 0) {
+	    call ft_trneqs (Memc[output], Memc[logfile], Memc[graphics],
+	        otable, ntable, wtflag, addscatter, tol, itmax, interactive,
+		high, low, niterate, grow, clgetb ("log_fit"),
+		clgetb ("log_results"))
+	} else if (nstd > 0) {
+	    call eprintf ("ERROR: No observations could be matched with ")
+	    call eprintf ("the catalog entries\n")
+	} else {
+	    call eprintf ("ERROR: No observations could be read from ")
+	    call eprintf ("the observations files\n")
+	}
+
+	# Free all space.
+	call pr_free()
+	call mct_free (otable)
+	if (ntable != NULL)
+	    call mct_free (ntable)
+	call fntclsb (obslist)
+
+	call sfree (sp)
+end