Initial commit

1 files changed, 278 insertions, 0 deletions

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