Repatch (from linux) of OSX IRAF

10 files changed, 2684 insertions, 0 deletions

diff --git a/noao/digiphot/photcal/evaluate/README b/noao/digiphot/photcal/evaluate/README
new file mode 100644
index 00000000..879dfcee
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/README
@@ -0,0 +1,3 @@
+This subdirectory contains the code for evaluatin the fit found by FITCOEFFS
+either in the forward sense using the EVALUATE task or in the inverse sense
+using the INVEVAL task.
diff --git a/noao/digiphot/photcal/evaluate/invert.com b/noao/digiphot/photcal/evaluate/invert.com
new file mode 100644
index 00000000..96af60d9
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/invert.com
@@ -0,0 +1,15 @@
+pointer	py		# pointer to the reference equation values
+pointer	pyfit		# pointer to the fit equation values
+pointer	pa		# pointer to the parameters to be fit
+pointer	pdelta		# pointer to the parameter increments
+pointer	pda		# pointer to the temporary parameter increments
+pointer	palpha		# pointer to the accumulated matrix
+pointer	pbeta		# pointer to the accumulated right side vector
+pointer	pik		# pointer to working index array for matrix inversion
+pointer	pjk		# pointer to working index array for matrix inversion
+
+pointer	pyerr		# pointer to the error equation values
+pointer	pafit		# pointer to the temporary fit array
+
+common /invertcom / py, pyfit, pa, pdelta, pda, palpha, pbeta, pik, pjk
+common /invertcom / pyerr, pafit
diff --git a/noao/digiphot/photcal/evaluate/mkpkg b/noao/digiphot/photcal/evaluate/mkpkg
new file mode 100644
index 00000000..667d92e3
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/mkpkg
@@ -0,0 +1,18 @@
+#  The MKPKG file for the INVEVALUATE and PHOTPROC task.
+
+$checkout	libpkg.a ".."
+$update		libpkg.a
+$checkin	libpkg.a ".."
+$exit
+
+libpkg.a:
+	t_invertfit.x   "../lib/io.h" "../lib/parser.h" <math/nlfit.h>  \
+			<error.h>     "../lib/preval.h"
+	t_evalfit.x	"../lib/io.h" "../lib/parser.h" <math/nlfit.h>  \
+			<error.h>
+	phinvert.x	<mach.h> "invert.com" "../lib/parser.h"
+	pherrors.x	"../lib/parser.h" "invert.com"
+	phprint.x	"../lib/parser.h"
+	phcheck.x	"../lib/parser.h"  "../lib/preval.h"
+	phminv.f
+	;
diff --git a/noao/digiphot/photcal/evaluate/phcheck.x b/noao/digiphot/photcal/evaluate/phcheck.x
new file mode 100644
index 00000000..2f729930
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/phcheck.x
@@ -0,0 +1,240 @@
+include "../lib/parser.h"
+include "../lib/preval.h"
+
+# PH_SETEQN -- Determine whether a PHOTCAL expression includes any set
+# equations which contain references to catalog variables. Return a YES
+# or NO status.
+
+int procedure ph_seteqn (code)
+
+pointer	code			# the photcal equation code
+
+int	ip, ins
+pointer	setcode
+int	pr_geti(), pr_gsym(), ph_cateqn()
+pointer	pr_gsymp()
+
+begin
+	# Return NO if there are no defined set equations.
+	if (pr_geti (NSETEQS) <= 0)
+	    return (NO)
+
+	# Initialize.
+	ip = 0  
+	ins = Memi[code+ip]
+
+	# Loop over the equation parsing instructions.
+	while (ins != PEV_EOC) {
+
+	    switch (ins) {
+	    case PEV_SETEQ:
+		ip = ip + 1
+		setcode = pr_gsymp (pr_gsym (Memi[code+ip], PTY_SETEQ),
+		    PSEQRPNEQ)
+		if (ph_cateqn (setcode) == YES)
+		    return (YES)
+	    case PEV_NUMBER, PEV_CATVAR, PEV_OBSVAR, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		# do nothing
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[code+ip]
+	}
+
+	return (NO)
+end
+
+
+# PH_CATEQN -- Given a PHOTCAL set equation determine whether there are any
+# references in it to catalog variables. Return a YES or NO status.
+
+int procedure ph_cateqn (code)
+
+pointer	code		# pointer to the equation code
+
+int	ip, ins
+int	pr_geti()
+
+begin
+	# Return NO if there are no catalog variables.
+	if (pr_geti (NCATVARS) <= 0)
+	    return (NO)
+
+	# Initialize.
+	ip = 0  
+	ins = Memi[code+ip]
+
+	# Loop over the equation parsing instructions.
+	while (ins != PEV_EOC) {
+
+	    switch (ins) {
+	    case PEV_CATVAR:
+		return (YES)
+	    case  PEV_NUMBER, PEV_OBSVAR, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_SETEQ, PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		# do nothing
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[code+ip]
+	}
+
+	return (NO)
+end
+
+
+# PH_SETVAR -- Check to see whether the fit expression portion of a PHOTCAL
+# transformation equation contains references to set equations which include
+# previously unreferenced catalog variables. Return the number of previously
+# unreferenced set equations.
+
+int procedure ph_setvar (eqn, sym, cmap, omap, tuservars, uservars, usererrs,
+	nstdvars, nobsvars, eqset, maxnset, userset, nset)
+
+int	eqn			# the equation number
+pointer	sym			# pointer to the equation symbol
+pointer	cmap			# pointer to catalog variable column map
+pointer	omap			# pointer to observations variable column map
+int	tuservars[nstdvars,ARB]	# list of active catalog variables per equation
+int	uservars[ARB]		# list of active catalog variables
+int	usererrs[ARB]		# list of active observational error columns
+int	nstdvars		# total number of catalog variables
+int	nobsvars		# total number of observations variables
+int	eqset[maxnset,ARB]	# set equation table
+int	maxnset			# maximum number of set equations
+int	userset[ARB]		# the list of active set equations
+int	nset			# number of set equations
+
+bool	new
+int	i, ip, ins, nvar
+pointer	fitcode, setsym
+int	pr_gsym(), ph_catvar()
+pointer	pr_gsymp()
+
+begin
+	# Get the fit expression symbols.
+	fitcode = pr_gsymp (sym, PTEQRPNFIT)
+
+	# Initialize.
+	nvar = 0
+	ip = 0  
+	ins = Memi[fitcode+ip]
+
+	# Loop over the fit expression.
+	while (ins != PEV_EOC) {
+
+	    switch (ins) {
+	    case PEV_SETEQ:
+		ip = ip + 1
+		setsym = pr_gsym (Memi[fitcode+ip], PTY_SETEQ)
+		if (ph_catvar (eqn, setsym, cmap, omap, tuservars, uservars,
+		    usererrs, nstdvars, nobsvars) == YES) {
+		    new = true
+		    do i = 1, nset + nvar {
+			if (setsym != userset[i])
+			    next
+			eqset[i,eqn] = setsym
+			new = false
+			break
+		    }
+		    if (new) {
+			nvar = nvar + 1
+			eqset[nset+nvar,eqn] = setsym
+			userset[nset+nvar] = setsym 
+		    }
+		}
+	    case PEV_NUMBER, PEV_CATVAR, PEV_OBSVAR, PEV_PARAM:
+		ip = ip + 1
+	    case PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+	    default:
+		# do nothing
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[fitcode+ip]
+	}
+
+	return (nvar)
+end
+
+
+# PH_CATVAR -- Determine whether any catalog variables in the set
+# equation have not been previously referenced in the transformation
+# equations.
+
+int procedure ph_catvar (eqn, setsym, cmap, omap, tuservars, uservars,
+	usererrs, nstdvars, nobsvars)
+
+int	eqn			# the transformation equation number
+int	setsym			# the set equation symbol
+pointer	cmap			# pointer to catalog variable column map
+pointer	omap			# pointer to observations variable column map
+int	tuservars[nstdvars,ARB]	# active variables as a function of eqn
+int	uservars[ARB]		# currently active catalog variables
+int	usererrs[ARB]		# currently active observations varaiables errs
+int	nstdvars		# the total number of catalog variables
+int	nobsvars		# the total number of observations variables
+
+int	setcode, ip, ins, stat, vcol, ecol
+pointer	catsym, obsym
+int	pr_gsym(), pr_gsymi(), pr_findmap1()
+pointer	pr_gsymp()
+
+begin
+	# Get the set equation code.
+	setcode = pr_gsymp (setsym, PSEQRPNEQ)
+
+	# Initialize.
+	ip = 0  
+	ins = Memi[setcode+ip]
+	stat = NO
+
+	while (ins != PEV_EOC) {
+
+	    switch (ins) {
+	    case PEV_CATVAR:
+		ip = ip + 1
+		catsym = pr_gsym (Memi[setcode+ip] - nobsvars, PTY_CATVAR)
+		vcol = pr_gsymi (catsym, PINPCOL)
+		vcol = pr_findmap1 (cmap, vcol)
+		if (uservars[vcol] <= 0) {
+		    tuservars[vcol,eqn] = -(vcol + nobsvars)
+		    uservars[vcol] = -(vcol + nobsvars)
+		    stat = YES
+		}
+
+	    case PEV_OBSVAR:
+		ip = ip + 1
+		obsym = pr_gsym (Memi[setcode+ip], PTY_OBSVAR)
+		vcol = pr_gsymi (obsym, PINPCOL)
+		ecol = pr_gsymi (obsym, PINPERRCOL)
+		vcol = pr_findmap1 (omap, vcol)
+		if (! IS_INDEFI(ecol))
+		    ecol = pr_findmap1 (omap, ecol)
+		if (! IS_INDEFI(ecol))
+		    usererrs[vcol] = ecol
+
+	    case PEV_NUMBER, PEV_PARAM:
+		ip = ip + 1
+
+	    case PEV_SETEQ, PEV_EXTEQ, PEV_TRNEQ:
+		ip = ip + 1
+
+	    default:
+		# do nothing
+	    }
+
+	    ip = ip + 1
+	    ins = Memi[setcode+ip]
+	}
+
+	return (stat)
+end
diff --git a/noao/digiphot/photcal/evaluate/pherrors.x b/noao/digiphot/photcal/evaluate/pherrors.x
new file mode 100644
index 00000000..77a95391
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/pherrors.x
@@ -0,0 +1,278 @@
+include "../lib/parser.h"
+
+define	DET_TOL		1.0e-20
+
+# PH_IEREQN -- Compute error estimates for the photometric indices based
+# on the value of any user-supplied errors equations. Use the chi-fitting
+# method developed by Bevington.
+
+int procedure ph_iereqn (params, a, nobs, deltaa, aindex, errors, nstd,
+	seteqn, setvals, serrors, nset, nustd, eqindex, neq)
+
+pointer	params[ARB]	# input array of pointers to the fitted parameters
+real	a[ARB]		# array of observed and catalog variables
+int	nobs		# the number of observed variables
+real	deltaa[ARB]	# array of increments for the catalog variables
+int	aindex[ARB]	# list of active catalog variables
+real	errors[ARB]	# output array of error estimates
+int	nstd		# number of catalog variables to be fit
+int	seteqn[ARB]	# array of set equation codes
+real	setvals[ARB]	# the set equation values
+real	serrors[ARB]	# output array of set equation errors
+int	nset		# the number of set equations
+int	nustd		# the number of active catalog variables
+int	eqindex[ARB]	# array of equations indices
+int	neq		# number of equations
+
+int	i, j, sym, nerrors
+pointer	errcode
+real	rms, det
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval(), ph_accum(), ph_incrms()
+
+include "invert.com"
+
+begin
+	# Evaluate the reference and error equations.
+	nerrors = 0
+	do i = 1, neq {
+
+	    sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	    Memr[py+i-1] = pr_eval (pr_gsymp (sym, PTEQRPNREF), a,
+	        Memr[params[eqindex[i]]])
+	    if (IS_INDEFR(Memr[py+i-1]))
+		return (0)
+
+	    errcode = pr_gsymp (sym, PTEQRPNERROR)
+	    if (errcode == NULL)
+		Memr[pyerr+i-1] = 0.0 
+	    else {
+	        Memr[pyerr+i-1] = pr_eval (errcode, a,
+		    Memr[params[eqindex[i]]])
+	        if (IS_INDEFR(Memr[pyerr+i-1]))
+		    Memr[pyerr+i-1] = 0.0 
+	        else
+		    nerrors = nerrors + 1
+	    }
+	}
+
+	# Return if there are no error equations.
+	if (nerrors <= 0)
+	    return (0)
+
+	# Compute each equations contribution to the to the total error.
+	call aclrr (errors, nstd)
+	call aclrr (serrors, nset)
+	do i = 1, neq {
+
+	    # Add in the appropriate error term.
+	    if (Memr[pyerr+i-1] <= 0.0)
+		next
+	    call amovr (a[nobs+1], Memr[pafit], nstd)
+	    Memr[py+i-1] = Memr[py+i-1] + Memr[pyerr+i-1]
+
+	    # Accumulate the matrices and vectors, do the inversion, and
+	    # compute the new parameter increments.
+	    rms = ph_accum (params, Memr[py], Memr[pyfit], eqindex, neq, a,
+	        nobs, deltaa, aindex, Memr[pda], nstd, Memd[palpha],
+		Memr[pbeta], Memi[pik], Memi[pjk], nustd, det)
+
+	    # Compute the contribution to the sum of the squares of the errors.
+	    #if ((! IS_INDEFR(rms)) && (abs (det) >= DET_TOL)) {
+	    if (! IS_INDEFR(rms)) {
+	        rms = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, neq,
+		    a, nobs, Memr[pda], aindex, nstd, rms)
+		do j = 1, nstd {
+		    if (aindex[j] == 0)
+			next
+		    errors[j] = errors[j] + (a[nobs+j] - Memr[pafit+j-1]) ** 2
+		}
+		do j = 1, nset {
+		    serrors[j] = serrors[j] + (pr_eval (seteqn[j], a,
+		        Memr[params[1]]) - setvals[j]) ** 2
+		}
+	    }
+
+	    # Reset the error term.
+	    Memr[py+i-1] = Memr[py+i-1] - Memr[pyerr+i-1]
+	    call amovr (Memr[pafit], a[nobs+1], nstd)
+	}
+
+	# Compute the errors themselves.
+	do i = 1, nstd {
+	    if (errors[i] <= 0.0)
+		errors[i] = 0.0
+	    else
+	        errors[i] = sqrt (errors[i])
+	}
+	do i = 1, nset {
+	    if (serrors[i] <= 0.0)
+		serrors[i] = 0.0
+	    else
+	        serrors[i] = sqrt (serrors[i])
+	}
+
+
+	return (nerrors)
+end
+
+
+# PH_IERVAL -- Compute error estimates for the photometric indices based
+# on the value of any user-supplied errors equations.
+
+int procedure ph_ierval (params, a, eindex, nobs, deltaa, aindex, errors,
+	nstd, seteqn, setvals, serrors, nset, nustd, eqindex, neq)
+
+pointer	params[ARB]	# input array of pointers to the fitted parameters
+real	a[ARB]		# array of observed and catalog variables
+int	eindex[ARB]	# indices of observed and catalog variables with errors
+int	nobs		# the number of observed variables
+real	deltaa[ARB]	# array of increments for the catalog variables
+int	aindex[ARB]	# list of active catalog variables
+real	errors[ARB]	# output array of error estimates
+int	nstd		# number of catalog variables to be fit
+int	seteqn[ARB]	# array of set equation codes
+real	setvals[ARB]	# the set equation values
+real	serrors[ARB]	# output array of set equation errors
+int	nset		# the number of set equations
+int	nustd		# the number of active catalog variables
+int	eqindex[ARB]	# array of equation indices
+int	neq		# number of equations
+
+int	i, j, k, sym, nerrors
+real	atemp, rms, det
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval(), ph_accum(), ph_incrms()
+
+include "invert.com"
+
+begin
+	# Initialize.
+	nerrors = 0
+	call aclrr (errors, nstd)
+	call aclrr (serrors, nset)
+
+	# Loop over the observational variables.
+	do j = 1, nobs {
+
+	    # Use only variables with errors.
+	    if (eindex[j] <= 0)
+		next
+	    nerrors = nerrors + 1
+	    if (IS_INDEFR(a[eindex[j]]) || (a[eindex[j]] <= 0.0))
+		next
+
+	    atemp = a[j]
+	    a[j] = a[j] + a[eindex[j]]
+	    call amovr (a[nobs+1], Memr[pafit], nstd)
+
+	    # Evaluate the reference equations.
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        Memr[py+i-1] = pr_eval (pr_gsymp (sym, PTEQRPNREF), a,
+	            Memr[params[eqindex[i]]])
+	        if (IS_INDEFR(Memr[py+i-1]))
+		    return (0)
+	    }
+
+	    # Accumulate the matrices and vectors, do the inversion, and
+	    # compute the new parameter increments.
+	    #call eprintf ("errors before accum\n")
+	    rms = ph_accum (params, Memr[py], Memr[pyfit], eqindex, neq, a,
+	        nobs, deltaa, aindex, Memr[pda], nstd, Memd[palpha],
+		Memr[pbeta], Memi[pik], Memi[pjk], nustd, det)
+	    #call eprintf ("errors after accum\n")
+
+	    # Compute the contribution to the sum of the squares of the errors.
+	    #if ((! IS_INDEFR(rms)) && (abs (det) >= DET_TOL)) {
+	    if (! IS_INDEFR(rms)) {
+	        #call eprintf ("errors before accum\n")
+	        rms = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, neq,
+		    a, nobs, Memr[pda], aindex, nstd, rms)
+	        #call eprintf ("errors before accum\n")
+
+		do k = 1, nstd {
+		    if (aindex[k] <= 0)
+			next
+		    errors[k] = errors[k] + (a[nobs+k] - Memr[pafit+k-1]) ** 2
+		}
+
+	        a[j] = atemp
+		do k = 1, nset {
+		    serrors[k] = serrors[k] + (pr_eval (seteqn[k], a,
+		        Memr[params[1]]) - setvals[k]) ** 2
+		}
+	        call amovr (Memr[pafit], a[nobs+1], nstd)
+
+	    } else {
+
+	        a[j] = atemp
+	        call amovr (Memr[pafit], a[nobs+1], nstd)
+	    }
+	}
+
+	# Compute the errors themselves.
+	do i = 1, nstd {
+	    if (errors[i] <= 0.0)
+		errors[i] = 0.0
+	    else
+	        errors[i] = sqrt (errors[i])
+	}
+	do i = 1, nset {
+	    if (serrors[i] <= 0.0)
+		serrors[i] = 0.0
+	    else
+	        serrors[i] = sqrt (serrors[i])
+	}
+
+	return (nerrors)
+end
+
+
+# PH_ERVAL -- Compute the error in an equation by summing the effects of errors
+# in the observational variables.
+
+real procedure ph_erval (symfit, fitval, params, a, aindex, eindex, nobsvars)
+
+pointer	symfit			# pointer to the fit equation
+real	fitval			# the best fit value
+real	params[ARB]		# the fitted equation parameters
+real	a[ARB]			# the observed and catalog variable values
+int	aindex[ARB]		# the list of active variables
+int	eindex[ARB]		# the list 
+int	nobsvars		# the number of observed variables
+
+int	i, nerrors
+real	errval, afit
+real	pr_eval()
+
+begin
+	# Initialize.
+	nerrors = 0
+	errval = 0.0
+
+	# Loop over the observed variables.
+	do i = 1, nobsvars {
+
+	    # Skip observed variables with no errors.
+	    if (eindex[i] <= 0)
+		next
+	    nerrors = nerrors + 1
+	    if ((IS_INDEFR(a[eindex[i]])) || (a[eindex[i]] <= 0.0))
+		next
+
+	    # Evaluate the contribution of each observed variable to the
+	    # total error.
+	    afit = a[i]
+	    a[i] = a[i] + a[eindex[i]]
+	    errval = errval + (pr_eval (symfit, a, params) - fitval) ** 2
+	    a[i] = afit
+	}
+
+	if (nerrors <= 0)
+	    return (INDEFR)
+	else
+	    return (sqrt (errval))
+end
diff --git a/noao/digiphot/photcal/evaluate/phinvert.x b/noao/digiphot/photcal/evaluate/phinvert.x
new file mode 100644
index 00000000..9d7ee14b
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/phinvert.x
@@ -0,0 +1,619 @@
+include <mach.h>
+include "../lib/parser.h"
+
+# PH_IVINIT -- Preallocate the space required by the inversion code.
+
+procedure ph_ivinit (nstd, nustd, neq)
+
+int	nstd		# number of catalog variables
+int	nustd		# number of catalog variables to be fit
+int	neq		# number of equations
+
+include "invert.com"
+
+begin
+	call malloc (py, neq, TY_REAL)
+	call malloc (pyfit, neq, TY_REAL)
+	call malloc (pa, nstd, TY_REAL)
+	call malloc (pdelta, nstd, TY_REAL)
+	call malloc (pda, nstd, TY_REAL)
+
+	call malloc (palpha, nustd * nustd, TY_DOUBLE)
+	call malloc (pbeta, nustd, TY_REAL)
+	call malloc (pik, nustd, TY_INT)
+	call malloc (pjk, nustd, TY_INT)
+
+	call malloc (pyerr, neq, TY_REAL)
+	call malloc (pafit, nstd, TY_REAL)
+end
+
+
+# PH_IVFREE -- Free the space required by the inversion code.
+
+procedure ph_ivfree ()
+
+include "invert.com"
+
+begin
+        call mfree (py, TY_REAL)
+        call mfree (pyfit, TY_REAL)
+        call mfree (pa, TY_REAL)
+        call mfree (pdelta, TY_REAL)
+        call mfree (pda, TY_REAL)
+
+        call mfree (palpha, TY_DOUBLE)
+        call mfree (pbeta, TY_REAL)
+        call mfree (pik, TY_INT)
+        call mfree (pjk, TY_INT)
+
+        call mfree (pyerr, TY_REAL)
+        call mfree (pafit, TY_REAL)
+end
+
+
+# PH_OBJCHECK -- Check that the equations for this particular star are
+# invertable.
+
+int procedure ph_objcheck (params, a, vartable, nstdvars, nreq, eqset,
+	maxnset, vindex, nvar, eqindex, neq)
+
+pointer	params[ARB]		# array of pointers to the fitted parameters
+real	a[ARB]			# array of observed and catalog variables
+int	vartable[nstdvars,ARB]	# table of variables as a function of equation
+int	nstdvars		# the total number of catalog variables
+int	nreq			# the total number of reference equations
+int	eqset[maxnset,ARB]	# set equation table
+int	maxnset			# maximum number of set equations
+int	vindex[ARB]		# output index of variables
+int	nvar			# number of variables used in the equations
+int	eqindex[ARB]		# output index of equations
+int	neq			# number of equations to be reduced
+
+int	i, j, sym, ncat, nset
+real	rval
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval()
+
+begin
+	# Initialize
+	call aclri (vindex, nstdvars)
+	call aclri (eqindex, nreq)
+
+	# Evalute the reference equations.
+	neq = 0
+	do i = 1, nreq {
+	    sym = pr_gsym (i, PTY_TRNEQ)
+	    rval = pr_eval (pr_gsymp (sym, PTEQRPNREF), a, Memr[params[i]])
+	    if (IS_INDEFR(rval))
+		next
+	    neq = neq + 1
+	    eqindex[neq] = i
+	}
+
+	# If there is no data return.
+	if (neq <= 0)
+	    return (ERR)
+
+	# Determine which variables are used by the reduced set of equations.
+	do i = 1, neq {
+	    do j = 1, nstdvars {
+		if (vartable[j,eqindex[i]] == 0)
+		    next
+		vindex[j] = vartable[j,eqindex[i]]
+	    }
+	}
+
+	# Deterine which set equations are used by the reduced set of equations
+	nset = 0
+	do j = 1, maxnset {
+	    do i = 1, neq {
+		if (eqset[j,eqindex[i]] == 0)
+		    next
+		nset = nset + 1
+		break
+	    }
+	}
+
+	# Count the number of variables.
+	nvar = 0
+	ncat = 0
+	do j = 1, nstdvars {
+	    if (vindex[j] == 0)
+		next
+	    nvar = nvar + 1
+	    if (vindex[j] > 0)
+		ncat = ncat + 1
+	}
+
+	if ((ncat + nset) > neq)
+	    return (ERR)
+
+	return (OK)
+end
+
+
+define	MAX_NITER1	10		# maximum number of iterations
+define	MAX_NITER2	50		# maximum number of trials
+define	DET_TOL		1.0E-20		# minimum value of the determinant
+#define	DET_TOL		0.0		# minimum value of the determinant
+define	MIN_DELTA	0.01		# minimum absolute value of
+					# parameter increments
+
+# PH_INVERT -- Invert the transformation to compute the standard indices.
+
+int procedure ph_invert (params, a, nobs, deltaa, aindex, nstd,
+	nustd, eqindex, nueq)
+
+pointer	params[ARB]	# input array of pointers to the fitted parameters
+real	a[ARB]		# array of observed and catalog variables
+int	nobs		# the number of observed variables
+real	deltaa[ARB]	# array of increments for the catalog variables
+int	aindex[ARB]	# index of active catalog variables
+int	nstd		# total number of catalog variables
+int	nustd		# number of catalog variables to be fit
+int	eqindex[ARB]	# the equation index
+int	nueq		# total number of equations used
+
+int	i, sym, niter
+real	stdev1, stdev2, det, rms
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval(), ph_accum(), ph_incrms()
+
+include "invert.com"
+
+begin
+	# Evalute the reference equations.
+	do i = 1, nueq {
+	    sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	    Memr[py+i-1] = pr_eval (pr_gsymp (sym, PTEQRPNREF), a,
+	        Memr[params[eqindex[i]]])
+	    if (IS_INDEFR(Memr[py+i-1]))
+		return (ERR)
+	}
+
+	# Initialize the parameter increments. This will be incremented
+	# each time through the fitting loop.
+	call amovr (deltaa, Memr[pdelta], nstd)
+
+	# Accumulate the matrices and vectors, do the inversion and compute
+	# the new parameter increments. The fit will terminate when the
+	# determinant of the inversion matrix becomes < 1.0e-20, if the
+	# standard deviation of the fit begins to increase, if the rms of
+	# the fit < EPSILONR, or the maximum number of iterations is
+	# exceeded, whichever comes first. 
+
+	niter = 0
+
+	repeat {
+
+	    # Compute the curvature currection. Return INDEFR if there are
+	    # bad data in the fit. Terminate the fit if the determinant of
+	    # the curvature matrix is too close to zero.
+
+	    stdev1 = ph_accum (params, Memr[py], Memr[pyfit], eqindex, nueq,
+	        a, nobs, Memr[pdelta], aindex, Memr[pda], nstd,
+		Memd[palpha], Memr[pbeta], Memi[pik], Memi[pjk], nustd, det)
+
+	    #call eprintf ("acc: niter=%d det=%g stdev1=%g\n")
+		#call pargi (niter+1)
+		#call pargr (det)
+		#call pargr (stdev1)
+
+	    # Return if there is INDEF data in the fit.
+	    if (IS_INDEFR(stdev1))
+	        return (ERR)
+
+	    # Check the size of the determinant but force at least one fit.
+	    if ((abs (det) < DET_TOL) && (niter > 0))
+		break
+
+	    # Find the new parameter values.
+	    call amovr (a, Memr[pa], nstd)
+	    stdev2 = ph_incrms (params, Memr[py], Memr[pyfit], eqindex, nueq,
+	        a, nobs, Memr[pda], aindex, nstd, stdev1)
+
+	    #call eprintf ("inc: niter=%d det=%g stdev2=%g\n")
+		#call pargi (niter+1)
+		#call pargr (det)
+		#call pargr (stdev2)
+
+	    # Check the new values.
+	    if (IS_INDEFR(stdev2) || (stdev2 >= stdev1)) {
+		if (niter == 0)
+	            return (ERR)
+		else {
+		    call amovr (Memr[pa], a, nstd)
+		    return (OK)
+		}
+	    }
+
+	    # User the new deltas and increment the fit counters.
+	    call anegr (Memr[pda], Memr[pdelta], nstd)
+	    call ph_deltamin (Memr[pdelta], MIN_DELTA, Memr[pdelta], nstd)
+	    niter = niter + 1
+	    rms = sqrt (stdev2)
+
+	} until ((niter == MAX_NITER1) || (rms <= EPSILONR))
+
+	return (OK)
+end
+
+
+# PH_ACCUM -- Accumulate the matrix of second derivatives and vector
+# of first derivatives  required for parabolic expansion of the rms
+# non-linear least squares fitting technique. This code is a modification
+# of the Bevington CHIFIT subroutine where the reduced chi-squared has
+# been replaced by the rms.  The original CHIFIT cannot be used to fit the
+# case when there are n data points and n unknowns since the number of
+# degrees of freedom is  zero and hence so is the reduced chi-squared.
+
+real procedure ph_accum (params, y, yfit, eqindex, neq, a, nobs, deltaa,
+	aindex, da, nstd, alpha, beta, ik, jk, nustd, det)
+
+pointer	params[ARB]	    # input array of ptrs to the fitted parameters
+real	y[ARB]		    # array of reference equation values
+real	yfit[ARB]	    # array of fitted reference equation values
+int	eqindex[ARB]	    # the equation indices
+int	neq		    # number of equations to be inverted
+real	a[ARB]		    # array of observed and catalog variables
+int	nobs		    # the number of observed variables
+real	deltaa[ARB]	    # array of increments for the catalog variables
+int	aindex[ARB]	    # index of active catalog variables
+real	da[ARB]		    # new array of parameter increments
+int	nstd		    # number of catalog variables
+double	alpha[nustd,ARB]    # the matrix of second derivatives
+real	beta[ARB]	    # the vector of first derivatives
+int	ik[ARB]		    # working array for the matrix inversion
+int	jk[ARB]		    # working array for the matrix inversion
+int	nustd		    # The number of catalog variables to be fit
+real	det		    # determinant of inverted matrix
+
+int	i, j, k, nj, nk, sym
+real	aj, ak, rms1, rms2, rms3
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval()
+
+begin
+	# Compute the initial rms, checking for INDEF valued variables.
+	rms1 = 0.0
+	do i = 1, neq {
+	    sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	    yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+	        Memr[params[eqindex[i]]])
+	    if (IS_INDEFR(yfit[i]))
+		return (INDEFR)
+	    rms1 = rms1 + (y[i] - yfit[i]) ** 2
+	}
+	rms1 = rms1 / neq
+
+	nj = 0
+	do j = 1, nstd {
+
+	    # Check the status of the parameter.
+	    if (aindex[j] == 0)
+		next
+	    nj = nj + 1
+
+	    # Increment each parameter.
+	    aj = a[j+nobs]
+	    a[j+nobs] = aj + deltaa[j]
+
+	    # Compute a new rms.
+	    rms2 = 0.0
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		    Memr[params[eqindex[i]]])
+	        rms2 = rms2 + (y[i] - yfit[i]) ** 2
+	    }
+	    rms2 = rms2 / neq
+
+	    # Begin accumulating the diagonal elements .
+	    alpha[nj,nj] = rms2 - 2.0 * rms1
+	    beta[nj] = -rms2
+
+	    # Accumulate the non-diagonal elements.
+	    nk = 0
+	    do k = 1, nstd {
+
+		if (aindex[k] == 0)
+		    next
+		nk = nk + 1
+
+		if ((nk - nj) == 0)
+		    next
+	        else if ((nk - nj) < 0) {
+		    alpha[nk,nj] = (alpha[nk,nj] - rms2) / 2.0
+		    alpha[nj,nk] = alpha[nk,nj]
+		    next
+	        }
+
+		alpha[nj,nk] = rms1 - rms2
+		ak = a[k+nobs]
+		a[k+nobs] = ak + deltaa[k]
+
+	        rms3 = 0.0
+	        do i = 1, neq {
+	            sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	            yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		        Memr[params[eqindex[i]]])
+	            rms3 = rms3 + (y[i] - yfit[i]) ** 2
+	        }
+	        rms3 = rms3 / neq
+
+		alpha[nj,nk] = alpha[nj,nk] + rms3
+		a[k+nobs] = ak
+	    }
+
+	    # Continue accumulating the diagonal elements.
+	    a[j+nobs] = aj - deltaa[j]
+	    rms3 = 0.0
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		    Memr[params[eqindex[i]]])
+	        rms3 = rms3 + (y[i] - yfit[i]) ** 2
+	    }
+	    rms3 = rms3 / neq
+
+	    a[j+nobs] = aj
+	    alpha[nj,nj] = (alpha[nj,nj] + rms3) / 2.0
+	    beta[nj] = (beta[nj] + rms3) / 4.0
+	}
+
+	# Eliminate any curvature from the matrix of second derivatives.
+	do j = 1, nj {
+	    if (alpha[j,j] > 0.0)
+		next
+	    if (alpha[j,j] < 0.0)
+		alpha[j,j] = -alpha[j,j]
+	    else
+		alpha[j,j] = 0.01
+	    do k = 1, nk {
+		if ((k - j) == 0)
+		    next
+		alpha[j,k] = 0.0
+		alpha[k,j] = 0.0
+	    }
+	}
+
+	# Invert the matrix.
+	call phminv (alpha, ik, jk, nj, det)
+
+	# Increment the parameters.
+	nj = 0
+	do j = 1, nstd {
+	    da[j] = 0.0
+	    if (aindex[j] == 0)
+		next
+	    nj = nj + 1
+	    nk = 0
+	    do k = 1, nstd {
+	        if (aindex[k] == 0)
+		    next
+		nk = nk + 1
+		da[j] = da[j] + beta[nk] * alpha[nj,nk]
+	    }
+	    da[j] = 0.2 * da[j] * deltaa[j]
+	}
+
+	# If the determinate is too small increment the parameters by
+	# deltas and try again.
+	#if (abs (det) < DET_TOL) {
+	    #call eprintf ("using approx\n")
+	    #do j = 1, nstd {
+	        #if (aindex[j] == 0) {
+		    #da[j] = 0.0
+		    #next
+		#}
+		#call eprintf ("i=%d dain=%g ")
+		    #call pargi (j)
+		    #call pargr (da[j])
+		#if (da[j] > 0.0)
+		    #da[j] = abs (deltaa[j])
+		#else if (da[j] < 0.0)
+		    #da[j] = -abs (deltaa[j])
+		#else
+		    #da[j] = 0.0
+		#call eprintf ("daout=%g\n")
+		    #call pargr (da[j])
+	    #}
+	#}
+
+	return (rms1)
+end
+
+
+# PH_INCRMS -- Increment the parameters until three values of the rms are found
+# which bracket the best fitting data point and fit a parabola to the three
+# different rms points.
+
+real procedure ph_incrms (params, y, yfit, eqindex, neq, a, nobs, da, aindex,
+	nstd, rms1)
+
+pointer	params[ARB]	# input array of ptrs to the fitted parameters
+real	y[ARB]		# the array of reference equation values
+real	yfit[ARB]	# the array of fitted reference equation values
+int	eqindex[ARB]	# list of active equations
+int	neq		# the number of equations
+real	a[ARB]		# array of observed and fitted variables
+int	nobs		# number of observed variables
+real	da[ARB]		# the parameter increments
+int	aindex[ARB]	# the index of active catalog variables
+int	nstd		# number of catalog variables
+real	rms1		# the first rms point
+
+int	j, i, sym, niter
+real	orms2, rms2, orms3, rms3, delta, rms
+int	pr_gsym()
+pointer	pr_gsymp()
+real	pr_eval()
+
+begin
+	# Adjust the parameters.
+	rms = rms1
+	do j = 1, nstd {
+	    if (aindex[j] == 0)
+		next
+	    a[j+nobs] = a[j+nobs] + da[j]
+	}
+
+	# Alter the parameter increments until the rms starts to decrease.
+	orms2 = MAX_REAL
+	niter = 0
+	repeat {
+
+	    rms2 = 0.0
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		    Memr[params[eqindex[i]]])
+	        rms2 = rms2 + (y[i] - yfit[i]) ** 2
+	    }
+	    rms2 = rms2 / neq
+	    #call eprintf ("    niter=%d rms1=%g rms2=%g\n")
+		#call pargi (niter)
+		#call pargr (rms1)
+		#call pargr (rms2)
+	    if (rms2 <= 0.0)
+		break
+	    if ((rms1 - rms2) >= 0.0)  
+		break
+
+	    # If rms2 does not decrease and does not change from one iteration
+	    # to the next we are probably near the precision limits of the
+	    # computer or the computed curvature has the wrong sign. In that
+	    # case quit.
+
+	    if (orms2 < MAX_REAL) {
+		if (abs ((rms2 - orms2) / orms2) < EPSILONR)
+		    return (rms)
+	    }
+
+	    orms2 = rms2
+	    niter = niter + 1
+	    if (niter >= MAX_NITER2)
+		return (rms)
+
+	    do j = 1, nstd {
+	        if (aindex[j] == 0)
+		    next
+		#da[j] = da[j] / 2.0
+		a[j+nobs] = a[j+nobs] - da[j]
+	    }
+
+
+	}
+
+	# Alter the parameter increments until the rms starts to increase.
+	orms3 = MAX_REAL
+	niter = 0
+	repeat {
+
+	    do j = 1, nstd {
+		if (aindex[j] == 0)
+		    next
+	        a[j+nobs] = a[j+nobs] + da[j]
+	    }
+
+	    rms3 = 0.0
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		    Memr[params[eqindex[i]]])
+	        rms3 = rms3 + (y[i] - yfit[i]) ** 2
+	    }
+	    rms3 = rms3 / neq
+	    #call eprintf ("    niter=%d rms1=%g rms2=%g rms3=%g\n")
+		#call pargi (niter)
+		#call pargr (rms1)
+		#call pargr (rms2)
+		#call pargr (rms3)
+	    if ((rms3 - rms2) >= 0.0)
+		break
+
+	    if (orms3 < MAX_REAL) {
+		if (abs ((rms3 - orms3) / orms3) < EPSILONR)
+		    return (rms)
+	    }
+
+	    niter = niter + 1
+	    if (niter >= MAX_NITER2)
+		return (rms)
+
+
+	    orms3 = rms3
+	    rms1 = rms2
+	    rms2 = rms3
+
+	}
+
+	# Fit a parabola to the three values of the rms that bracket the fit.
+	if (rms3 <= rms2)
+	    delta = 1.0
+	else
+	    delta = 1.0 / (1.0 + (rms1 - rms2) / (rms3 - rms2)) + 0.5
+	do j = 1, nstd {
+	    if (aindex[j] == 0)
+		next
+	    a[nobs+j] = a[nobs+j] - delta * da[j]
+	}
+
+	rms = 0.0
+	do i = 1, neq {
+	    sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	    yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+	        Memr[params[eqindex[i]]])
+	    rms = rms + (y[i] - yfit[i]) ** 2
+	}
+	rms = rms / neq
+
+	if ((rms2 - rms) < 0.0) {
+	    do j = 1, nstd {
+	        if (aindex[j] == 0)
+		    next
+	        a[nobs+j] = a[nobs+j] + (delta - 1.) * da[j]
+	    }
+	    do i = 1, neq {
+	        sym = pr_gsym (eqindex[i], PTY_TRNEQ)
+	        yfit[i] = pr_eval (pr_gsymp (sym, PTEQRPNFIT), a,
+		    Memr[params[eqindex[i]]])
+	    }
+	    rms = rms2
+	}
+
+	#call eprintf ("    incr: rms1=%g rms2=%g rms3=%g rms=%g\n")
+	    #call pargr (rms1)
+	    #call pargr (rms2)
+	    #call pargr (rms3)
+	    #call pargr (rms)
+
+	return (rms)
+end
+
+
+# PH_DELTAMIN -- Check to make sure that the absolute value of the deltaa
+# is always greater than or equal to min_delta.
+
+procedure ph_deltamin (a, min_delta, b, npix)
+
+real	a[ARB]			# input vector
+real	min_delta		# minimum permitted absolute value 
+real	b[ARB]			# output vector
+int	npix			# number of points
+
+int	i
+
+begin
+	do i = 1, npix {
+	    if (abs(a[i]) >= min_delta)
+		b[i] = a[i]
+	    else if (a[i] < 0.0)
+		b[i] = -min_delta
+	    else
+		b[i] = min_delta
+	}
+end
diff --git a/noao/digiphot/photcal/evaluate/phminv.f b/noao/digiphot/photcal/evaluate/phminv.f
new file mode 100644
index 00000000..44d7430e
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/phminv.f
@@ -0,0 +1,114 @@
+C SUBROUTINE PHMINV.F 
+C
+C SOURCE								       
+C   BEVINGTON, PAGES 302-303.
+C
+C PURPOSE
+C   INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT
+C
+C USAGE 
+C   CALL PHMINV (ARRAY, IK, JK, NORDER, DET)
+C
+C DESCRIPTION OF PARAMETERS
+C   ARRAY  - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE
+C   IK     - WORK ARRAY REQUIRED BY PHMINV
+C   JK     - WORK ARRAY REQUIRED BY PHMINV
+C   NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT)
+C   DET    - DETERMINANT OF INPUT MATRIX
+C
+C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
+C   NONE
+C
+C COMMENT
+C   DIMENSION STATEMENT VALID FOR NORDER UP TO 10
+C
+	SUBROUTINE PHMINV (ARRAY,IK,JK,NORDER,DET)
+	DOUBLE PRECISION ARRAY,AMAX,SAVE
+	DIMENSION ARRAY(NORDER,NORDER),IK(NORDER),JK(NORDER)
+C
+10	DET=1.
+11	DO 100 K=1,NORDER
+C
+C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
+C
+	AMAX=0. 
+21	CONTINUE
+	IK(K) = K
+	JK(K) = K
+	DO 30 I=K,NORDER
+	DO 30 J=K,NORDER
+23	IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30
+24	AMAX=ARRAY(I,J) 
+	IK(K)=I 
+	JK(K)=J 
+30	CONTINUE
+C
+C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
+C
+41	I=IK(K) 
+	IF (I-K) 21,51,43
+43	DO 50 J=1,NORDER
+	SAVE=ARRAY(K,J) 
+	ARRAY(K,J)=ARRAY(I,J)
+50	ARRAY(I,J)=-SAVE
+51	J=JK(K) 
+	IF (J-K) 21,600,53
+53	DO 60 I=1,NORDER
+	SAVE=ARRAY(I,K) 
+	ARRAY(I,K)=ARRAY(I,J)
+60	ARRAY(I,J)=-SAVE
+
+C
+C CHECK FOR SINGULAR MATRIX
+C
+600  	IF (AMAX) 61,610,61
+610	DET=0.
+	DO 700 I=1,NORDER
+	IF (I-K) 630,700,630
+630	ARRAY(I,K)=0.
+700	CONTINUE
+	DO 800 J=1,NORDER
+	IF (J-K) 730,800,730
+730	ARRAY(K,J)=0.
+800	CONTINUE
+	ARRAY(K,K)=0.
+	GOTO 100
+C
+C ACCUMULATE ELEMENTS OF INVERSE MATRIX 
+C
+61	DO 70 I=1,NORDER
+	IF (I-K) 63,70,63
+63	ARRAY(I,K)=-ARRAY(I,K)/AMAX
+70	CONTINUE
+71	DO 80 I=1,NORDER
+	DO 80 J=1,NORDER
+	IF (I-K) 74,80,74
+74	IF (J-K) 75,80,75
+75	ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
+80	CONTINUE
+81	DO 90 J=1,NORDER
+	IF (J-K) 83,90,83
+83	ARRAY(K,J)=ARRAY(K,J)/AMAX
+90	CONTINUE
+	ARRAY(K,K)=1./AMAX
+100	DET=DET*AMAX
+C
+C RESTORE ORDERING OF MATRIX
+C
+101	DO 130 L=1,NORDER
+	K=NORDER-L+1
+	J=IK(K) 
+	IF (J-K) 111,111,105
+105	DO 110 I=1,NORDER
+	SAVE=ARRAY(I,K) 
+	ARRAY(I,K)=-ARRAY(I,J)
+110	ARRAY(I,J)=SAVE 
+111	I=JK(K) 
+	IF (I-K) 130,130,113
+113	DO 120 J=1,NORDER
+	SAVE=ARRAY(K,J) 
+	ARRAY(K,J)=-ARRAY(I,J)
+120	ARRAY(I,J)=SAVE 
+130	CONTINUE
+140	RETURN
+	END
diff --git a/noao/digiphot/photcal/evaluate/phprint.x b/noao/digiphot/photcal/evaluate/phprint.x
new file mode 100644
index 00000000..94150792
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/phprint.x
@@ -0,0 +1,110 @@
+include "../lib/parser.h"
+
+# PH_MKPLIST -- Construct the list of variables to be printed.
+
+int procedure ph_mkplist (plist, cmap, omap, nobsvars, psym, pcols,
+	max_len_plist)
+
+int	plist		# pointer to the list of variables
+int	cmap		# catalog column to data column mapping
+int	omap		# observations column to data column mapping
+int	nobsvars	# number of observations variables
+int	psym[ARB]	# the output array of variable symbols
+int	pcols[ARB]	# offset in the data array of the symbol
+int	max_len_plist	# the maximum length of the variables list
+
+int	len_plist, sym, col
+pointer	sp, pname
+int	fntgfnb(), pr_getsym(), pr_gsymi(), pr_findmap1()
+
+begin
+	call smark (sp)
+	call salloc (pname, SZ_FNAME, TY_CHAR)
+
+	len_plist = 0
+	while (fntgfnb (plist, Memc[pname], SZ_FNAME) != EOF) {
+	    if (len_plist >= max_len_plist)
+		break
+	    sym = pr_getsym (Memc[pname])
+	    if (IS_INDEFI(sym))
+		next
+	    switch (pr_gsymi (sym, PSYMTYPE)) {
+	    case PTY_CATVAR:
+		psym[len_plist+1] = sym
+		col  = pr_gsymi (sym, PINPCOL)
+		pcols[len_plist+1] = pr_findmap1 (cmap, col) + nobsvars 
+	    case PTY_OBSVAR:
+		psym[len_plist+1] = sym
+		col  = pr_gsymi (sym, PINPCOL)
+		pcols[len_plist+1] = pr_findmap1 (omap, col)
+	    case PTY_SETEQ:
+		psym[len_plist+1] = sym
+		pcols[len_plist+1] = 0 
+	    default:
+		next
+	    }
+	    len_plist = len_plist + 1
+	}
+
+	call sfree (sp)
+	return (len_plist)
+end
+
+
+# PH_OFORMAT -- Construct the output formatstr string.
+
+procedure ph_oformatstr (getid, ncols, formatstr, maxch)
+
+int	getid		# output the object id
+int	ncols		# number of columns in the output text file
+char	formatstr[ARB]	# the output formatstr string
+int	maxch		# maximum number of characters in the formatstr string
+
+int	i, fcol
+
+begin
+	if (getid == YES) {
+	    call strcpy ("%-10s ", formatstr, maxch)
+	    fcol = 2
+	} else {
+	    formatstr[1] = EOS
+	    fcol = 1
+	}
+
+	do i = fcol, ncols {
+	    if (i == ncols)
+		call strcat ("%-7.3f\n", formatstr, maxch)
+	    else
+		call strcat ("%-7.3f ", formatstr, maxch)
+	}
+end
+
+
+# PH_OFIELDS -- Count the number of fields in the formatstr string.
+
+int procedure ph_ofields (formatstr)
+
+char	formatstr[ARB]		# the input formatstr string
+
+char	percent
+int	ip, findex, nfields
+int	stridx()
+data	percent /'%'/
+
+begin
+	nfields = 0
+
+	ip = 1
+	while (formatstr[ip] != EOS) {
+	    findex = stridx (percent, formatstr[ip])
+	    if (findex == 0)
+		break
+	    ip = findex + ip
+	    if (formatstr[ip] == percent)
+		ip = ip + 1
+	    else
+	        nfields = nfields + 1
+	}
+
+	return (nfields)
+end
diff --git a/noao/digiphot/photcal/evaluate/t_evalfit.x b/noao/digiphot/photcal/evaluate/t_evalfit.x
new file mode 100644
index 00000000..daf407c2
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/t_evalfit.x
@@ -0,0 +1,495 @@
+include	<error.h>
+include <math/nlfit.h>
+include	"../lib/io.h"
+include	"../lib/parser.h"
+
+# Define the pointer Mem
+define	MEMP	Memi
+
+# T_EVALFIT - Main photometry processing task. This task will convert 
+# intrumental photometric indices into standard indices, by using the 
+# same configuration table and coefficients found by FITCOEFFS.
+
+procedure t_evalfit()
+
+pointer	observations		# pointer to the observations files list
+pointer	catalogs		# pointer to the catalogs files list
+pointer	config			# pointer to configuration file name
+pointer paramfile		# pointer to fitted parameters file name
+pointer calibfile		# pointer to the output file name
+int	type			# type of output to be processed
+int	etype			# algorithm for computing errors
+pointer	print			# pointer to the output variables list
+pointer	formatstr		# pointer to the output format string
+pointer	catdir			# pointer to the standard catalogs directory
+
+int	i, j, getid, matchid, vcol, ecol, pindex, dummy, stat
+int	obslist, stdlist, plist, ofd, ifd, sym, symvar, ncols
+int	ntrneqs, nparams, nstd, nobsvars, nvars, len_plist
+pointer	sp, input, starname, dummyname, stdtable, cmap, omap
+pointer	vars, uservars, usererrs, fsym, rsym, esym, params, errors, psym, pcols
+real	ref, fit, errval, resid, chisqr, rms, pval
+
+bool	clgetb()
+int	clgwrd(), open(), io_gcoeffs(), io_gobs(), pr_parse(), pr_geti()
+int	pr_gsym(), pr_gsymi(), pr_findmap1(), pr_gvari(), fntopnb()
+int	fntgfnb(), fntlenb(), ph_mkplist(), ph_header(), ph_ofields()
+pointer	pr_gsymp(), pr_xgetname()
+real	pr_eval(), ph_erval()
+errchk	io_gcoeffs()
+
+begin
+	# Allocate space for file names and character strings.
+	call smark (sp)
+	call salloc (observations, SZ_LINE, TY_CHAR)
+	call salloc (catalogs, SZ_LINE, TY_CHAR)
+	call salloc (config, SZ_FNAME, TY_CHAR)
+	call salloc (paramfile, SZ_FNAME, TY_CHAR)
+	call salloc (calibfile, SZ_FNAME, TY_CHAR)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (starname, SZ_LINE, TY_CHAR)
+	call salloc (dummyname, SZ_LINE, TY_CHAR)
+	call salloc (print, SZ_LINE, TY_CHAR)
+	call salloc (formatstr, SZ_LINE, TY_CHAR)
+	call salloc (catdir, SZ_FNAME, TY_CHAR)
+
+	# Get the input and output file lists and  names.
+	call clgstr ("observations", Memc[observations], SZ_LINE)
+	call clgstr ("config", Memc[config], SZ_FNAME)
+	call clgstr ("parameters", Memc[paramfile], SZ_FNAME)
+	call clgstr ("calib", Memc[calibfile], SZ_FNAME)
+	call clgstr ("catalogs", Memc[catalogs], SZ_LINE)
+	call clgstr ("catdir", Memc[catdir], SZ_FNAME)
+
+	# Get output type flags.
+	etype = clgwrd ("errors", Memc[dummyname], SZ_LINE, ERR_OPTIONS)
+	type = clgwrd ("objects", Memc[dummyname], SZ_LINE, TYPE_STRING)
+	call clgstr ("print", Memc[print], SZ_LINE)
+	call clgstr ("format", Memc[formatstr], SZ_LINE)
+
+	# Open the output file.
+	iferr {
+	    if (clgetb ("append"))
+	        ofd = open (Memc[calibfile], APPEND, TEXT_FILE)
+	    else
+	        ofd = open (Memc[calibfile], NEW_FILE, TEXT_FILE)
+	} then {
+	    call erract (EA_WARN)
+	    call sfree (sp)
+	    return
+	}
+
+	# Parse the configuration table.
+	if (pr_parse (Memc[config]) == ERR) {
+	    call eprintf ("Error: Cannot parse the configuration file\n")
+	    call close (ofd)
+	    call sfree (sp)
+	    return
+	}
+
+	# Get number of variables and allocate memory for their values.
+	# Program stars do not have standard indices (last values in
+	# the array), and they will remain undefined, and this will give
+	# undefined equation values for equations that use them.
+
+	nobsvars = pr_geti (NOBSVARS)
+	nvars = nobsvars + pr_geti (NCATVARS)
+	call salloc (vars, nvars, TY_REAL)
+	call salloc (uservars, nvars, TY_INT)
+	call aclri (Memi[uservars], nvars)
+	call salloc (usererrs, nvars, TY_INT)
+	call aclri (Memi[usererrs], nvars)
+
+	# Map observational and catalog columns. This has to be done here in
+	# order to avoid mapping columns every time a new program
+	# star is read.
+
+	call pr_catmap (cmap, dummy)
+	call pr_obsmap (omap, dummy)
+
+	# Get all the parameter values for all the equations. This has
+	# to be done outside the main loop to minimize disk i/o.
+	# Otherwise it would be necessary to fetch them for each
+	# star observation.
+
+	ntrneqs = pr_geti (NTRNEQS)
+	call salloc (fsym, ntrneqs, TY_POINTER)
+	call salloc (rsym, ntrneqs, TY_POINTER)
+	call salloc (esym, ntrneqs, TY_POINTER)
+	call salloc (params, ntrneqs, TY_POINTER)
+	call salloc (errors, ntrneqs, TY_POINTER)
+	do i = 1, ntrneqs {
+
+	    # Get equation symbol and number of parameters for it.
+	    sym = pr_gsym (i, PTY_TRNEQ)
+	    nparams = pr_gsymi (sym, PTEQNPAR)
+
+	    # Get the fit, reference and error equation pointers.
+	    Memi[fsym+i-1] = pr_gsymp (sym, PTEQRPNFIT)
+	    Memi[rsym+i-1] = pr_gsymp (sym, PTEQRPNREF)
+	    Memi[esym+i-1] = pr_gsymp (sym, PTEQRPNERROR)
+
+	    # Determine which catalog and observational variables were
+	    # actually used in the reference equations and which have
+	    # defined error columns.
+
+	    do j = 1, pr_gsymi (sym, PTEQNRCAT) + pr_gsymi(sym, PTEQNROBS) {
+		symvar = pr_gvari (sym, j, PTEQREFVAR)
+		vcol = pr_gsymi (symvar, PINPCOL)
+		ecol = pr_gsymi (symvar, PINPERRCOL)
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) {
+		    vcol = pr_findmap1 (cmap, vcol) + nobsvars
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (cmap, ecol) + nobsvars
+		} else {
+		    vcol = pr_findmap1 (omap, vcol)
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (omap, ecol)
+		}
+		Memi[uservars+vcol-1] = vcol
+		if (! IS_INDEFI(ecol))
+		    Memi[usererrs+vcol-1] = ecol
+	    }
+
+	    # Determine which catalog and observational variables were
+	    # actually used in the fit equations and which have
+	    # defined error columns.
+
+	    do j = 1, pr_gsymi (sym, PTEQNFCAT) + pr_gsymi(sym, PTEQNFOBS) {
+		symvar = pr_gvari (sym, j, PTEQFITVAR)
+		vcol = pr_gsymi (symvar, PINPCOL)
+		ecol = pr_gsymi (symvar, PINPERRCOL)
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) {
+		    vcol = pr_findmap1 (cmap, vcol) + nobsvars
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (cmap, ecol) + nobsvars
+		} else {
+		    vcol = pr_findmap1 (omap, vcol)
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (omap, ecol)
+		}
+		Memi[uservars+vcol-1] = vcol
+		if (! IS_INDEFI(ecol))
+		    Memi[usererrs+vcol-1] = ecol
+	    }
+
+	    # Allocate space for parameter values and errors,
+	    # for the current equation, and read them from the
+	    # parameters file.
+
+	    call salloc (MEMP[params+i-1], nparams, TY_REAL)
+	    call salloc (MEMP[errors+i-1], nparams, TY_REAL)
+	    iferr {
+	        if (io_gcoeffs (Memc[paramfile], sym, stat, chisqr, rms,
+		    Memr[MEMP[params+i-1]], Memr[MEMP[errors+i-1]],
+		    nparams) != nparams) {
+		    call eprintf ("Warning: Error reading parameters for ")
+		    call eprintf ("equation %s from %s\n")
+			call pargstr (Memc[pr_xgetname (sym)])
+			call pargstr (Memc[paramfile])
+		    call amovkr (INDEFR, Memr[MEMP[params+i-1]], nparams)
+		    call amovkr (INDEFR, Memr[MEMP[errors+i-1]], nparams)
+		}
+	    } then {
+		call erract (EA_WARN)
+		call close (ofd)
+		call pr_unmap (cmap)
+		call pr_unmap (omap)
+		call pr_free()
+		call sfree (sp)
+		return
+	    }
+
+	    # Issue a warning if some of the equations in the system did
+	    # not converge but proceed with the evaluation.
+	    # didn't converge.
+
+	    if (stat != DONE) {
+		call eprintf (
+		    "Warning: The solution for equation %s did not converge")
+		    call pargstr (Memc[pr_xgetname (sym)])
+	    }
+	}
+
+	# Decide whether to use id and/or catalog matching or not.
+	if (pr_geti (MINCOL) <= 1) {
+	    getid = NO
+	    matchid = NO
+	    stdtable = NULL
+	}  else if (Memc[catalogs] == EOS) {
+	    getid = YES
+	    matchid = NO
+	    stdtable = NULL
+	} else {
+	    getid = YES
+	    matchid = YES
+	}
+
+	# Get the list of optional variables to be printed. These variables
+	# may include any of the catalog or observations variables or the
+	# set equations. Variables which match any of these catagories
+	# will be ommited from the list to be printed.
+
+	plist = fntopnb (Memc[print], NO)
+	len_plist = fntlenb (plist)
+	if (len_plist > 0) {
+	    call salloc (psym, len_plist, TY_INT)
+	    call salloc (pcols, len_plist, TY_INT)
+	    len_plist = ph_mkplist (plist, cmap, omap, nobsvars, Memi[psym],
+	        Memi[pcols], len_plist)
+	}
+	call fntclsb (plist)
+
+	# Print the output header.
+	ncols = ph_header (ofd, Memc[observations], Memc[catalogs],
+	    Memc[config], Memc[paramfile], type, getid, Memi[psym],
+	    Memi[pcols], len_plist, etype, matchid)
+
+	# Set the formatstr string.
+	if (Memc[formatstr] == EOS) {
+	    call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE)
+	} else if (ph_ofields (Memc[formatstr]) != ncols) {
+	    call eprintf ("Warning: The number of formatstr string fields ")
+	    call eprintf ("does not match the number of output columns\n")
+	    call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE)
+	}
+
+	# If catalog matching is enabled, read in the catalog data.
+	if (matchid == YES) {
+	    stdlist = fntopnb (Memc[catalogs], NO)
+	    call io_gcatdat (Memc[catdir], stdlist, stdtable, nstd, dummy)
+	    call fntclsb (stdlist)
+	} 
+
+	# Loop over the observation files.
+	obslist = fntopnb (Memc[observations], NO)
+	while (fntgfnb (obslist, Memc[input], SZ_FNAME) != EOF) {
+
+	    # Open the input file.
+	    iferr (ifd = open (Memc[input], READ_ONLY, TEXT_FILE)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Get observations and names for program stars. The call
+	    # to io_getline_init() is necessary since it's not called
+	    # inside io_gobs(), and it should be called at least once
+	    # per input file.
+
+	    call io_getline_init()
+	    while (io_gobs (ifd, stdtable, omap, type, Memr[vars], nvars,
+		   getid, Memc[starname], Memc[dummyname], SZ_LINE) != EOF) {
+
+		# Print the star name.
+		call fprintf (ofd, Memc[formatstr])
+		if (getid == YES)
+		    call pargstr (Memc[starname])
+
+		# Print the extra variables.
+		do i = 1, len_plist {
+		    pindex = Memi[pcols+i-1]
+		    if (pindex > 0)
+			pval = Memr[vars+pindex-1]
+		    else
+			pval = pr_eval (Memi[psym+i-1], Memr[vars],
+			    Memr[MEMP[params]])
+		    call pargr (pval)
+		}
+
+		# Loop over all the transformatstrion equations.
+		do i = 1, pr_geti (NTRNEQS) {
+
+		    # Get the equation symbol.
+		    sym = pr_gsym (i, PTY_TRNEQ)
+
+		    # Compute the fit. 
+		    fit = pr_eval (Memi[fsym+i-1], Memr[vars],
+		        Memr[MEMP[params+i-1]])
+		    call pargr (fit)
+
+		    # Compute the error.
+		    if (etype != ERR_UNDEFINED) {
+		        if (IS_INDEFR(fit))
+			    errval = INDEFR
+			switch (etype) {
+		        case ERR_OBSERRORS:
+		            if (IS_INDEFR(fit))
+			        errval = INDEFR
+			    else
+			        errval = ph_erval (Memi[fsym+i-1], fit,
+			            Memr[MEMP[params+i-1]], Memr[vars],
+				    Memi[uservars], Memi[usererrs], nobsvars)
+		        case ERR_EQUATIONS:
+			    if (Memi[esym+i-1] == NULL)
+			        errval = INDEFR
+			    else
+			        errval = pr_eval (Memi[esym+i-1], Memr[vars],
+			            Memr[MEMP[params+i-1]])
+		        default:
+			    errval = INDEFR
+			}
+			call pargr (errval)
+		    }
+
+		    # Compute the residual.
+		    if (getid == NO || (matchid == YES &&
+		        type != TYPE_PROGRAM)) {
+		        if (IS_INDEFR (fit))
+			    resid = INDEFR
+		        else {
+		            ref = pr_eval (Memi[rsym+i-1], Memr[vars],
+		                Memr[MEMP[params+i-1]])
+			    if (IS_INDEFR(ref))
+			        resid = INDEFR
+			    else
+		                resid = ref - fit
+		        }
+		        call pargr (resid)
+		    }
+		}
+	    }
+
+	    # Close the input file.
+	    call close (ifd)
+	}
+
+	# Free memory.
+	call sfree (sp)
+	call pr_free()
+	call pr_unmap (cmap)
+	call pr_unmap (omap)
+
+	# Close all files.
+	if (stdtable != NULL)
+	    call stclose (stdtable)
+	call close (ofd)
+	call fntclsb (obslist)
+end
+
+
+# PH_HEADER - Print the output file header.
+
+int procedure ph_header (fd, observations, catalogs, config, paramfile, type,
+	getid, psym, pcols, len_plist, etype, matchid)
+
+int	fd			# output file descriptor
+char	observations[ARB]	# observation file list
+char	catalogs[ARB]		# catalog file list
+char	config[ARB]		# configuration file name
+char	paramfile[ARB]		# fitted parameters file
+int	type			# type of object to output
+int	getid			# output the object id
+int	psym[ARB]		# list of additional variables to be output
+int	pcols[ARB]		# columns numbers of additional variables
+int	len_plist		# length of symbol list
+int	etype			# type of error output
+int	matchid			# is it possible to match ids
+
+int	i, olist, slist, ncols
+pointer sp, time, name
+int	fntopnb(), fntgfnb(), pr_geti(), pr_gsym()
+long	clktime()
+pointer	pr_gsymc(), pr_gsymp(), pr_xgetname()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (time, SZ_LINE, TY_CHAR)
+	call salloc (name, SZ_FNAME, TY_CHAR)
+
+	# Add the time stamp.
+	call cnvtime (clktime (0), Memc[time], SZ_LINE)
+	call fprintf (fd, "\n# %s\n")
+	    call pargstr (Memc[time])
+
+	# Add the observation file names.
+	olist = fntopnb (observations, NO)
+	call fprintf (fd, "# List of observations files:\n") 
+	while (fntgfnb (olist, Memc[name], SZ_FNAME) != EOF) {
+	    call fprintf (fd, "#\t\t%s\n")
+		call pargstr (Memc[name])
+	}
+	call fntclsb (olist)
+
+	# Add the catalog file names.
+	if (matchid == YES) {
+	    slist = fntopnb (catalogs, NO)
+	    call fprintf (fd, "# List of catalog files:\n") 
+	    while (fntgfnb (slist, Memc[name], SZ_FNAME) != EOF) {
+		call fprintf (fd, "#\t\t%s\n")
+		    call pargstr (Memc[name])
+	    }
+	    call fntclsb (slist)
+	}
+
+	# Add the configuration file name.
+	call fprintf (fd, "# Config:\t%s\n")
+	    call pargstr (config)
+
+	# Add the parameters file name.
+	call fprintf (fd, "# Parameters:\t%s\n")
+	    call pargstr (paramfile)
+
+	# Add the output options.
+	call fprintf (fd, "#\n")
+	if (matchid == YES) {
+	    if (type == TYPE_ALL)
+		call fprintf (fd,
+		    "# Computed indices for program and standard objects\n")
+	    else if (type == TYPE_PROGRAM)
+		call fprintf (fd,
+		    "# Computed indices for program objects only\n")
+	    else if (type == TYPE_STANDARDS)
+		call fprintf (fd,
+		    "# Computed indices for standard objects only\n")
+	} else
+	    call fprintf (fd,
+	        "# Computed indices for program and standard objects\n")
+
+	# Print the optional id header.
+	call fprintf (fd, "#\n")
+	call fprintf (fd, "# Columns:\n")
+	if (getid == YES) {
+	    call fprintf (fd, "#\t1\tobject id\n")
+	    ncols = 1
+	} else
+	    ncols = 0
+
+	# Print headers for the variables in the print list and set any set
+	# equation offsets to pointers.
+	do i = 1, len_plist {
+	    ncols = ncols + 1
+	    call fprintf (fd, "#\t%d\t%s\n")
+		call pargi (ncols)
+		call pargstr (Memc[pr_xgetname(psym[i])])
+	    if (pcols[i] <= 0)
+		psym[i] = pr_gsymp (psym[i], PSEQRPNEQ)
+	}
+
+
+	# Print the equation headers.
+	do i = 1, pr_geti (NTRNEQS) {
+	    ncols = ncols + 1
+	    call fprintf (fd, "#\t%d\t%s\n")
+		call pargi (ncols)
+	        call pargstr (Memc[pr_gsymc(pr_gsym (i,PTY_TRNEQ),PTEQREF)])
+	    if (etype != ERR_UNDEFINED) {
+	        ncols = ncols + 1
+	        call fprintf (fd, "#\t%d\terror(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_gsymc(pr_gsym (i,PTY_TRNEQ),PTEQREF)])
+	    }
+	    if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) {
+	        ncols = ncols + 1
+	        call fprintf (fd, "#\t%d\tresid(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_gsymc (pr_gsym(i,PTY_TRNEQ),PTEQREF)])
+	    }
+	}
+
+	call fprintf (fd,"\n\n")
+
+	call sfree (sp)
+
+	return (ncols)
+end
diff --git a/noao/digiphot/photcal/evaluate/t_invertfit.x b/noao/digiphot/photcal/evaluate/t_invertfit.x
new file mode 100644
index 00000000..532d7f08
--- /dev/null
+++ b/noao/digiphot/photcal/evaluate/t_invertfit.x
@@ -0,0 +1,792 @@
+include	<error.h>
+include <math/nlfit.h>
+include	"../lib/io.h"
+include	"../lib/parser.h"
+include "../lib/preval.h"
+
+# Define the pointer Mem
+define	MEMP	Memi
+
+# T_INVERTFIT - INVERTFIT converts intrumental photometric indices into
+# standard indices by using the configuration file, the coefficients
+# determined by FITPARAMS and inverting the transformations.
+
+procedure t_invertfit ()
+
+pointer	observations		# pointer to the observations file list
+pointer	catalogs		# pointer to the catalogs file list
+pointer	config			# pointer to configuration file name
+pointer paramfile		# pointer to fitted parameters file name
+pointer calibfile		# pointer to the output file name
+int	type			# type of output to be processed
+int	etype			# algorithm for computing the errors
+pointer	print			# pointer to output variables list
+pointer	formatstr		# pointer to the output format string
+pointer	catdir			# pointer to the standard star directory
+
+int	i, j, vcol, ecol, pindex, dummy, stat, getid, matchid
+int	obslist, stdlist, plist, ofd, ifd, sym, symvar, ncols, nstd, nset
+int	ntrneqs, nparams, nstdvars, nustdvars, nobsvars, nvars, nueq, maxnset
+int	len_plist, refcode
+pointer	sp, input, starname, dummyname, stdtable, omap, cmap
+pointer	vars, eqvartable, uservars, usererrs, userset, eqset, tvars
+pointer	dtvars, avtvars, ervars, svars, servars, params, errors, psym, pcols
+pointer	varindex, eqindex	
+real	resid, chisqr, rms, pval
+
+bool	clgetb()
+int	fntopnb(), fntgfnb(), clgwrd(), open(), io_gcoeffs(), io_gobs()
+int	pr_parse(), pr_geti(), pr_gsym(), pr_gsymi(), pr_gvari(), ph_objcheck()
+int	ph_invert(), pr_findmap1(), ph_iereqn(), ph_ierval(), ph_seteqn()
+int	ph_mkplist(), fntlenb(), ph_ofields(), ph_iheader(), ph_setvar()
+pointer	pr_xgetname(), pr_gsymp()
+real	pr_eval()
+errchk	io_gcoeffs()
+
+begin
+	# Allocate space for file names and character strings.
+
+	call smark (sp)
+	call salloc (observations, SZ_LINE, TY_CHAR)
+	call salloc (catalogs, SZ_LINE, TY_CHAR)
+	call salloc (config, SZ_FNAME, TY_CHAR)
+	call salloc (paramfile, SZ_FNAME, TY_CHAR)
+	call salloc (calibfile, SZ_FNAME, TY_CHAR)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (starname, SZ_LINE, TY_CHAR)
+	call salloc (dummyname, SZ_LINE, TY_CHAR)
+	call salloc (print, SZ_LINE, TY_CHAR)
+	call salloc (formatstr, SZ_LINE, TY_CHAR)
+	call salloc (catdir, SZ_FNAME, TY_CHAR)
+
+	# Get the observations list, the catalog list, the configuration
+	# file, the parameters file, and the output file names.
+
+	call clgstr ("observations", Memc[observations], SZ_LINE)
+	call clgstr ("config", Memc[config], SZ_FNAME)
+	call clgstr ("parameters", Memc[paramfile], SZ_FNAME)
+	call clgstr ("calib", Memc[calibfile], SZ_FNAME)
+	call clgstr ("catalogs", Memc[catalogs], SZ_LINE)
+	call clgstr ("catdir", Memc[catdir], SZ_LINE)
+
+	# Get the output type flags.
+
+	etype = clgwrd ("errors", Memc[dummyname], SZ_LINE, ERR_OPTIONS)
+	type = clgwrd ("objects", Memc[dummyname], SZ_LINE, TYPE_STRING)
+	call clgstr ("print", Memc[print], SZ_LINE)
+	call clgstr ("format", Memc[formatstr], SZ_LINE)
+
+	# Open the output file.
+
+	iferr {
+	    if (clgetb ("append"))
+	        ofd = open (Memc[calibfile], APPEND, TEXT_FILE)
+	    else
+	        ofd = open (Memc[calibfile], NEW_FILE, TEXT_FILE)
+	} then {
+	    call erract (EA_WARN)
+	    call sfree (sp)
+	    return
+	}
+
+	# Parse the configuration file.
+
+	if (pr_parse (Memc[config]) == ERR) {
+	    call eprintf ("Error: Cannot parse the configuration file\n")
+	    call close (ofd)
+	    call sfree (sp)
+	    return
+	}
+
+	# Fetch the total number of observed and catalog variables and
+	# the total number of equations in the configuration file.
+
+	nobsvars = pr_geti (NOBSVARS)
+	nstdvars = pr_geti (NCATVARS)
+	nvars = nobsvars + nstdvars
+	ntrneqs = pr_geti (NTRNEQS)
+
+	# Map observations file and catalog file columns.
+
+	call pr_catmap (cmap, dummy)
+	call pr_obsmap (omap, dummy)
+
+	# Check which catalog and observations variables are actually
+	# used in the equations to be inverted and determine whether the
+	# system of equations can actually be inverted. This step is the
+	# first pass through the system of equations. A second pass is
+	# necessary to deal with the set equations. Use this first pass
+	# to fetch the fitted parameters for the equations.
+
+	call salloc (eqvartable, nstdvars * ntrneqs, TY_INT)
+	call aclri (Memi[eqvartable], nstdvars * ntrneqs)
+	call salloc (uservars, nstdvars, TY_INT)
+	call aclri (Memi[uservars], nstdvars)
+	call salloc (usererrs, nobsvars, TY_INT)
+	call aclri (Memi[usererrs], nobsvars)
+	call salloc (params, ntrneqs, TY_POINTER)
+	call salloc (errors, ntrneqs, TY_POINTER)
+
+	do i = 1, ntrneqs {
+
+	    # Get the equation symbol.
+	    sym = pr_gsym (i, PTY_TRNEQ)
+
+	    # Get the reference equation symbol.
+	    refcode = pr_gsymp (sym, PTEQRPNREF)
+
+	    # Quit if there are catalog variables in the function expression
+	    # or if there are set equations containing references to the
+	    # catalog variables in the reference equation.
+
+	    if ((pr_gsymi (sym, PTEQNRCAT) > 0) || (ph_seteqn (refcode) ==
+	        YES)) {
+		call eprintf ("Error: Cannot invert equations with catalog ")
+		call eprintf ("variables in the function expression\n")
+		call pr_unmap (cmap)
+		call pr_unmap (omap)
+	        call close (ofd)
+	        call sfree (sp)
+	        return
+	    }
+
+	    # Determine which catalog and observational variables were
+	    # actually used in the reference equations and which have
+	    # defined error columns.
+
+	    do j = 1, pr_gsymi (sym, PTEQNRCAT) + pr_gsymi(sym, PTEQNROBS) {
+		symvar = pr_gvari (sym, j, PTEQREFVAR)
+		vcol = pr_gsymi (symvar, PINPCOL)
+		ecol = pr_gsymi (symvar, PINPERRCOL)
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) {
+		    vcol = pr_findmap1 (cmap, vcol)
+		    Memi[eqvartable+(i-1)*nstdvars+vcol-1] = vcol + nobsvars
+		    Memi[uservars+vcol-1] = vcol + nobsvars
+		} else {
+		    vcol = pr_findmap1 (omap, vcol)
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (omap, ecol)
+		    if (! IS_INDEFI(ecol))
+		        Memi[usererrs+vcol-1] = ecol
+		}
+	    }
+
+	    # Determine which catalog and observational variables were
+	    # actually used in the fit equations and which have
+	    # defined error columns.
+
+	    do j = 1, pr_gsymi (sym, PTEQNFCAT) + pr_gsymi(sym, PTEQNFOBS) {
+		symvar = pr_gvari (sym, j, PTEQFITVAR)
+		vcol = pr_gsymi (symvar, PINPCOL)
+		ecol = pr_gsymi (symvar, PINPERRCOL)
+		if (pr_gsymi (symvar, PSYMTYPE) == PTY_CATVAR) {
+		    vcol = pr_findmap1 (cmap, vcol)
+		    Memi[eqvartable+(i-1)*nstdvars+vcol-1] = vcol + nobsvars
+		    Memi[uservars+vcol-1] = vcol + nobsvars
+		} else {
+		    vcol = pr_findmap1 (omap, vcol)
+		    if (! IS_INDEFI(ecol))
+		        ecol = pr_findmap1 (omap, ecol)
+		    if (! IS_INDEFI(ecol))
+		        Memi[usererrs+vcol-1] = ecol
+		}
+	    }
+
+	    # Get the number of parameters for the equation.
+	    # Allocate space for parameter values and errors, for the
+	    # current equation, and read them from the coefficient file.
+
+	    nparams = pr_gsymi (sym, PTEQNPAR)
+	    call salloc (MEMP[params+i-1], nparams, TY_REAL)
+	    call salloc (MEMP[errors+i-1], nparams, TY_REAL)
+	    iferr {
+	        if (io_gcoeffs (Memc[paramfile], sym, stat, chisqr, rms,
+	            Memr[MEMP[params+i-1]], Memr[MEMP[errors+i-1]], nparams) !=
+		    nparams) {
+		    call eprintf ("Warning: Error reading parameters for ")
+		    call eprintf ("equation %s from %s\n")
+			call pargstr (Memc[pr_xgetname(sym)])
+			call pargstr (Memc[paramfile])
+		    call amovkr (INDEFR, Memr[MEMP[params+i-1]], nparams)
+		    call amovkr (INDEFR, Memr[MEMP[errors+i-1]], nparams)
+		}
+	    } then {
+		call erract (EA_WARN)
+		call close (ofd)
+		call pr_unmap (cmap)
+		call pr_unmap (omap)
+		call pr_free()
+		call sfree (sp)
+		return
+	    }
+
+	    # Issue a warning if any of the equations in the system to be
+	    # inverted did not converge but proceed with the fit.
+
+	    if (stat != DONE) {
+		call eprintf (
+		    "Warning: The solution for equation %s did not converge")
+		    call pargstr (Memc[pr_xgetname(sym)])
+	    }
+
+	}
+
+	# Count the number of catalog variables used in the transformation
+	# equations.
+
+	nustdvars = 0 
+	do i = 1, nstdvars {
+	    if (Memi[uservars+i-1] != 0)
+		nustdvars = nustdvars + 1
+	}
+
+	# If the number of equations is less than the number of referenced
+	# catalog variables it is not possible to invert the transformation
+	# equations.
+
+	if (nustdvars > ntrneqs) {
+	    call eprintf ("Error: The number of equations is less than ")
+	    call eprintf ("the number of unknowns\n")
+	    call close (ofd)
+	    call pr_unmap (cmap)
+	    call pr_unmap (omap)
+	    call pr_free()
+	    call sfree (sp)
+	    return
+	}
+
+	# Loop over the transformation equations a second time searching for
+	# references to the set equations in the fit expressions. If a set
+	# equation reference is found, check to see whether it contains any
+	# reference to a catalog variable which is not referenced elsewhere
+	# in the transformation equations.  Recompute the number of catalog
+	# variables used in the fit equations.
+
+	maxnset = max (1, pr_geti (NSETEQS))
+	call salloc (eqset, maxnset * ntrneqs, TY_INT)
+	call aclri (Memi[eqset], maxnset * ntrneqs)
+	nset = 0
+	if (pr_geti (NSETEQS) > 0) {
+
+	    # Find the number of independent set equations.
+	    call salloc (userset, pr_geti (NSETEQS), TY_INT)
+	    call amovki (NULL, Memi[userset], pr_geti (NSETEQS))
+	    do i = 1, ntrneqs {
+	        sym = pr_gsym (i, PTY_TRNEQ)
+	        nset = nset + ph_setvar (i, sym, cmap, omap, Memi[eqvartable],
+		    Memi[uservars], Memi[usererrs], nstdvars, nobsvars,
+		    Memi[eqset], pr_geti (NSETEQS), Memi[userset], nset)
+	    }
+
+	    # Is the system invertable?
+	    if ((nustdvars + nset) > ntrneqs) {
+	        call eprintf ("Error: The number of equations is less than ")
+	        call eprintf ("the number of unknowns\n")
+	        call close (ofd)
+	        call pr_unmap (cmap)
+	        call pr_unmap (omap)
+	        call pr_free()
+	        call sfree (sp)
+	        return
+	    }
+
+	    # Recompute the number of independent catalog variables.
+	    nustdvars = 0 
+	    do i = 1, nstdvars {
+	        if (Memi[uservars+i-1] != 0)
+		    nustdvars = nustdvars + 1
+	    }
+
+	}
+
+
+	# Decide whether to use id/catalog matching or not. If id matching
+	# is enabled fetch the catalog data.
+
+	if (pr_geti (MINCOL) <= 1) {
+	    getid = NO
+	    matchid = NO
+	    stdtable = NULL
+	} else if (Memc[catalogs] == EOS) {
+	    getid = YES
+	    matchid = NO
+	    stdtable = NULL
+	} else {
+	    getid = YES
+	    matchid = YES
+	    stdtable = NULL
+	}
+
+	# Get the list of optional variables to be printed. These variables
+	# may include any of the catalog or observations variables or the
+	# set equations. Variables which do not match any of these categories
+	# will be ommitted from the list to be printed.
+
+	plist = fntopnb (Memc[print], NO)
+	len_plist = fntlenb (plist)
+	if (len_plist > 0) {
+	    call salloc (psym, len_plist, TY_INT)
+	    call salloc (pcols, len_plist, TY_INT)
+	    len_plist = ph_mkplist (plist, cmap, omap, nobsvars, Memi[psym],
+		Memi[pcols], len_plist)
+
+	}
+	call fntclsb (plist)
+
+	# Print the output header.
+	ncols = ph_iheader (ofd, Memc[observations], Memc[catalogs],
+	    Memc[config], Memc[paramfile], type, getid, Memi[psym],
+	    Memi[pcols], len_plist, Memi[uservars], nstdvars,
+	    Memi[userset], nset, etype, matchid)
+
+	# Set the format string.
+	if (Memc[formatstr] == EOS) {
+	    call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE)
+	} else if (ph_ofields (Memc[formatstr]) != ncols) {
+	    call eprintf ("Warning: The number of format string fields ")
+	    call eprintf ("does not match the number of output columns\n")
+	    call ph_oformatstr (getid, ncols, Memc[formatstr], SZ_LINE)
+	}
+
+	# Read in the catalog data.
+	if (matchid == YES) {
+	    stdlist = fntopnb (Memc[catalogs], NO)
+	    call io_gcatdat (Memc[catdir], stdlist, stdtable, nstd, dummy)
+	    call fntclsb (stdlist)
+	}
+
+	# Compute the initial values for the catalog variables and initial
+	# values for the delta-variable estimates. If a catalog is present
+	# the initial value of each catalog variable is the average value of
+	# that variable in the catalog. The standard deviation of the catalog
+	# variables is used as the initial value for the catalog variable
+	# increments.
+
+	call salloc (avtvars, nstdvars, TY_REAL)
+	call salloc (dtvars, nstdvars, TY_REAL)
+	call ph_avstdvars (stdtable, Memi[uservars], Memr[avtvars],
+	    Memr[dtvars], nstdvars)
+
+	# Allocated space for the catalog variables and the fitted variables
+	# and their errors.
+	call salloc (vars, nvars, TY_REAL)
+	call salloc (tvars, nvars, TY_REAL)
+	if (nset > 0)
+	    call salloc (svars, nset, TY_REAL)
+	call salloc (ervars, nstdvars, TY_REAL)
+	if (nset > 0)
+	    call salloc (servars, nset, TY_REAL)
+	call salloc (varindex, nstdvars, TY_INT)
+	call salloc (eqindex, ntrneqs, TY_INT)
+
+	# Initialize the fit inversion code.
+	call ph_ivinit (nstdvars, nustdvars, ntrneqs)
+
+	# Loop over the observations files.
+	obslist = fntopnb (Memc[observations], NO)
+	while (fntgfnb (obslist, Memc[input], SZ_FNAME) != EOF) {
+
+	    # Open the input file.
+	    iferr (ifd = open (Memc[input], READ_ONLY, TEXT_FILE)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Get observations and names for program stars. The call
+	    # to io_getline_init() is necessary since it's not called
+	    # inside io_gobs(), and it should be called at least once
+	    # per input file.
+
+	    call io_getline_init()
+	    while (io_gobs (ifd, stdtable, omap, type, Memr[vars], nvars,
+		   getid, Memc[starname], Memc[dummyname], SZ_LINE) != EOF) {
+
+		# Print star name.
+		call fprintf (ofd, Memc[formatstr])
+		if (getid == YES)
+		    call pargstr (Memc[starname])
+		#call eprintf ("name=%s\n")
+		    #call pargstr (Memc[starname])
+
+		# Print the extra variables.
+		do i = 1, len_plist {
+		    pindex = Memi[pcols+i-1]
+		    if (pindex > 0)
+			pval = Memr[vars+pindex-1]
+		    else
+			pval = pr_eval (Memi[psym+i-1], Memr[vars],
+			    Memr[MEMP[params]])
+		    call pargr (pval)
+		}
+
+		# Initialize the fit.
+		call amovr (Memr[vars], Memr[tvars], nobsvars)
+		call amovr (Memr[avtvars],  Memr[tvars+nobsvars], nstdvars)
+
+		# Invert the transformations and compute the errors.
+		if (ph_objcheck (MEMP[params], Memr[tvars], Memi[eqvartable],
+		    nstdvars, ntrneqs, Memi[eqset], maxnset,
+		    Memi[varindex], nustdvars, Memi[eqindex],
+		    nueq) == ERR) {
+		    call amovkr (INDEFR, Memr[tvars+nobsvars], nstdvars)
+		    if (nset > 0)
+			call amovkr (INDEFR, Memr[svars], nset)
+		    call amovkr (INDEFR, Memr[ervars], nstdvars)
+		    if (nset > 0)
+			call amovkr (INDEFR, Memr[servars], nset)
+
+		} else if (ph_invert (MEMP[params], Memr[tvars], nobsvars,
+		    Memr[dtvars], Memi[varindex], nstdvars,
+		    nustdvars, Memi[eqindex], nueq) == ERR) {
+
+		    call amovkr (INDEFR, Memr[tvars+nobsvars], nstdvars)
+		    if (nset > 0)
+			call amovkr (INDEFR, Memr[svars], nset)
+		    call amovkr (INDEFR, Memr[ervars], nstdvars)
+		    if (nset > 0)
+			call amovkr (INDEFR, Memr[servars], nset)
+
+		} else {
+
+		    # Set any unfitted variables to INDEF.
+		    do i = nobsvars + 1, nvars {
+			if (Memi[varindex+i-nobsvars-1] == 0)
+			    Memr[tvars+i-1] = INDEFR
+		    }
+
+		    # Evaluate the set equations.
+		    if (nset > 0) {
+		        do i = 1, nset
+		            Memr[svars+i-1] = pr_eval (Memi[userset+i-1],
+			        Memr[tvars], Memr[MEMP[params+i-1]])
+		    }
+
+		    # Evaluate the errors.
+		    switch (etype) {
+		    case ERR_UNDEFINED:
+			call amovkr (INDEFR, Memr[ervars], nstdvars)
+			if (nset > 0)
+			    call amovkr (INDEFR, Memr[servars], nset)
+		    case ERR_EQUATIONS:
+		        if (ph_iereqn (MEMP[params], Memr[tvars], nobsvars,
+		            Memr[dtvars], Memi[varindex], Memr[ervars],
+			    nstdvars, Memi[userset], Memr[svars],
+			    Memr[servars], nset, nustdvars, Memi[eqindex],
+			    nueq) <= 0) {
+			    call amovkr (INDEFR, Memr[ervars], nstdvars)
+			    if (nset > 0)
+			        call amovkr (INDEFR, Memr[servars], nset)
+			}
+		    case ERR_OBSERRORS:
+		        if (ph_ierval (MEMP[params], Memr[tvars],
+			    Memi[usererrs], nobsvars, Memr[dtvars],
+			    Memi[varindex], Memr[ervars], nstdvars,
+			    Memi[userset], Memr[svars], Memr[servars], nset,
+			    nustdvars, Memi[eqindex], nueq) <= 0) {
+			    call amovkr (INDEFR, Memr[ervars], nstdvars)
+			    if (nset > 0)
+			        call amovkr (INDEFR, Memr[servars], nset)
+			}
+		    default:
+			call amovkr (INDEFR, Memr[ervars], nstdvars)
+		    }
+		}
+
+		# Write out the standard indices, errors and residuals.
+		do i = nobsvars + 1, nvars {
+
+		    # Skip catalog variables not used in the equation.
+		    if (Memi[uservars+i-nobsvars-1] <= 0)
+			next
+
+		    call pargr (Memr[tvars+i-1])
+		    if (etype != ERR_UNDEFINED) {
+			if (IS_INDEFR(Memr[tvars+i-1]))
+			    call pargr (INDEFR)
+			else
+			    call pargr (Memr[ervars+i-nobsvars-1])
+		    }
+
+		    # Compute the residual of the fit.
+		    if (getid == NO || (matchid == YES &&
+		        type != TYPE_PROGRAM)) {
+		        if (IS_INDEFR (Memr[vars+i-1]) ||
+		            IS_INDEFR (Memr[tvars+i-1]))
+			    resid = INDEFR
+		        else
+		            resid = Memr[vars+i-1] - Memr[tvars+i-1]
+		        call pargr (resid)
+		    }
+		}
+
+		# Write out the set equations.
+		do i = 1, nset {
+
+		    # Write out the set equation.
+		    call pargr (Memr[svars+i-1])
+
+		    # Evaluate the error. This is tricky at the moment.
+		    if (etype != ERR_UNDEFINED) {
+			if (IS_INDEFR(Memr[svars+i-1]))
+			    call pargr (INDEFR)
+			else
+			    call pargr (Memr[servars+i-1])
+		    }
+
+		    # Compute the residual of the fit.
+		    if (getid == NO || (matchid == YES &&
+		        type != TYPE_PROGRAM)) {
+		        resid = pr_eval (Memi[userset+i-1], Memr[vars],
+		            Memr[MEMP[params+i-1]])
+		        if (IS_INDEFR (Memr[svars+i-1]) || IS_INDEFR (resid))
+			    resid = INDEFR
+		        else
+		            resid = resid - Memr[svars+i-1]
+		        call pargr (resid)
+		    }
+		}
+	    }
+
+	    # Close the input file.
+	    call close (ifd)
+	}
+
+	# Free memory.
+	call sfree (sp)
+	call ph_ivfree()
+	call pr_free()
+	call pr_unmap (cmap)
+	call pr_unmap (omap)
+
+	# Close all the files.
+	if (stdtable != NULL)
+	    call stclose (stdtable)
+	call close (ofd)
+	call fntclsb (obslist)
+end
+
+
+# PH_IHEADER - Print the output file header.
+
+int procedure ph_iheader (fd, observations, catalogs, config, paramfile, type,
+	getid, psym, pcols, len_plist, catvars, ncatvars, userset, nset,
+	etype, matchid)
+
+int	fd			# output file descriptor
+char	observations[ARB]	# the list of observations files
+char	catalogs[ARB]		# the list of catalog files
+char	config[ARB]		# the configuration file name
+char	paramfile[ARB]		# the fitted parameters file
+int	type			# the type of object to output
+int	getid			# output the object id ?
+int	psym[ARB]		# list of additional symbols to be printed
+int	pcols[ARB]		# column numbers of additional output variables
+int	len_plist		# length of symbol list
+int	catvars[ARB]		# the list of active catalog variables
+int	ncatvars		# number of catalog variables
+int	userset[ARB]		# list of set equations included in fit
+int	nset			# number of set equations
+int	etype			# type of error output
+int	matchid			# is it possible to match ids ?
+
+int	i, olist, slist, ncols
+pointer sp, time, name, sym
+int	fntopnb(), fntgfnb(), pr_gsym()
+long	clktime()
+pointer	pr_xgetname(), pr_gsymp()
+
+begin
+	# Allocate some working space.
+	call smark (sp)
+	call salloc (time, SZ_LINE, TY_CHAR)
+	call salloc (name, SZ_FNAME, TY_CHAR)
+
+	# Add the time stamp.
+	call cnvtime (clktime (0), Memc[time], SZ_LINE)
+	call fprintf (fd, "\n# %s\n")
+	    call pargstr (Memc[time])
+
+	# Add the observation file names.
+	olist = fntopnb (observations, NO)
+	call fprintf (fd, "# List of observations files:\n") 
+	while (fntgfnb (olist, Memc[name], SZ_FNAME) != EOF) {
+	    call fprintf (fd, "#\t\t%s\n")
+		call pargstr (Memc[name])
+	}
+	call fntclsb (olist)
+
+	# Add the catalog file names.
+	if (matchid == YES) {
+	    slist = fntopnb (catalogs, NO)
+	    call fprintf (fd, "# Number of catalog files:\n") 
+	    while (fntgfnb (slist, Memc[name], SZ_FNAME) != EOF) {
+		call fprintf (fd, "#\t\t%s\n")
+		    call pargstr (Memc[name])
+	    }
+	    call fntclsb (slist)
+	}
+
+	# Add the configuration file name.
+	call fprintf (fd, "# Config:\t%s\n")
+	    call pargstr (config)
+
+	# Add the parameters file name.
+	call fprintf (fd, "# Parameters:\t%s\n")
+	    call pargstr (paramfile)
+
+	# Write the output options.
+	call fprintf (fd, "#\n")
+	if (matchid == YES) {
+	    if (type == TYPE_ALL)
+		call fprintf (fd,
+		    "# Computed indices for program and standard objects\n")
+	    else if (type == TYPE_PROGRAM)
+		call fprintf (fd,
+		    "# Computed indices for program objects only\n")
+	    else if (type == TYPE_STANDARDS)
+		call fprintf (fd,
+		    "# Computed indices for standard objects only\n")
+	} else
+	    call fprintf (fd,
+	        "# Computed indices for program and standard objects\n")
+
+	# Print the optional id header
+	call fprintf (fd, "#\n")
+	call fprintf (fd, "# Columns: \n")
+	if (getid == YES) {
+	    call fprintf (fd, "#\t1\tobject id\n")
+	    ncols = 1
+	} else
+	    ncols = 0
+
+	# Print headers for the variables in the print list and set any set
+	# equation offsets to pointers.
+	do i = 1, len_plist {
+	    ncols = ncols + 1
+	    call fprintf (fd, "#\t%d\t%s\n")
+		call pargi (ncols)
+		call pargstr (Memc[pr_xgetname(psym[i])])
+	    if (pcols[i] <= 0)
+		psym[i] = pr_gsymp (psym[i], PSEQRPNEQ)
+	}
+
+	# Print the catalog variables headers.
+	do i = 1, ncatvars {
+	    if (catvars[i] <= 0)
+		next
+	    sym = pr_gsym (i, PTY_CATVAR)
+	    ncols = ncols + 1
+	    call fprintf (fd, "#\t%d\t%s\n")
+		call pargi (ncols)
+	        call pargstr (Memc[pr_xgetname(sym)])
+	    if (etype != ERR_UNDEFINED) {
+		ncols = ncols + 1
+		call fprintf (fd, "#\t%d\terror(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_xgetname(sym)])
+	    }
+	    if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) {
+		ncols = ncols + 1
+		call fprintf (fd, "#\t%d\tresid(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_xgetname(sym)])
+	    }
+	}
+
+	# Print the set equation variables headers.
+	do i = 1, nset {
+	    ncols = ncols + 1
+	    call fprintf (fd, "#\t%d\t%s\n")
+		call pargi (ncols)
+		call pargstr (Memc[pr_xgetname(userset[i])])
+	    if (etype != ERR_UNDEFINED) {
+		ncols = ncols + 1
+		call fprintf (fd, "#\t%d\terror(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_xgetname(userset[i])])
+	    }
+	    if (getid == NO || (matchid == YES && type != TYPE_PROGRAM)) {
+		ncols = ncols + 1
+		call fprintf (fd, "#\t%d\tresid(%s)\n")
+		    call pargi (ncols)
+	            call pargstr (Memc[pr_xgetname(userset[i])])
+	    }
+	    userset[i] = pr_gsymp (userset[i], PSEQRPNEQ)
+	}
+
+	call fprintf (fd, "\n\n")
+
+	call sfree (sp)
+
+	return (ncols)
+end
+
+
+# PH_AVSTDVARS -- Compute the mean and sigma of the variables in the
+# standard table to use as initial guesses for the math routines.
+# If the standard catalog is undefined or empty the mean values of the
+# variables are set to 0.0 and the sigmas are set to 1.0. If there is only
+# one catalog value the sigma is also set to 1.0.
+
+procedure ph_avstdvars (stable, index, avg, sigma, npts)
+
+pointer	stable			# pointer to the symbol table
+int	index[ARB]		# index of active catalog variables
+real	avg[ARB]		# the output array of averages
+real	sigma[ARB]		# the output array of standard deviations
+int	npts			# number of points
+
+int	i, ndata
+pointer	sym, sp, n
+real	data
+pointer	sthead(), stnext()
+
+begin
+	# Initialize and return if the standard table is undefined.
+	call amovkr (0.0, avg, npts)
+	if (stable == NULL) {
+	    call amovkr (0.1, sigma, npts)
+	    return
+	}
+
+	# Allocate some temporary space and intialize the accumulators.
+	call smark (sp)
+	call salloc (n, npts, TY_INT)
+	call aclri (Memi[n], npts)
+	call aclrr (sigma, npts)
+
+	# Read the symbol table and accumlate the sums.
+	sym = sthead (stable)
+	while (sym != NULL) {
+	    do i = 1, npts {
+		if (index[i] == 0)
+		    next
+		data = Memr[P2R(sym+i-1)]
+		if (IS_INDEFR(data))
+		    next
+		avg[i] = avg[i] + data
+		sigma[i] = sigma[i] + data ** 2
+	        Memi[n+i-1] = Memi[n+i-1] + 1
+	    }
+	    sym = stnext (stable, sym)
+	}
+
+	# Compute the averages.
+	do i = 1, npts {
+	    ndata = Memi[n+i-1]
+	    if (ndata == 0)
+		sigma[i] = 0.1
+	    else if (ndata == 1)
+		sigma[i] = 0.1
+	    else {
+		avg[i] = avg[i] / ndata
+		sigma[i] = (sigma[i] - avg[i] * avg[i] * ndata) / (ndata - 1)
+		if (sigma[i] <= 0.0)
+		    sigma[i] = 0.1
+		else
+		    sigma[i] = min (0.1, sqrt (sigma[i]))
+	    }
+	}
+
+	call sfree (sp)
+end