path: root/noao/digiphot/daophot/psf/dpfitpsf.x

include <mach.h>
include <imhdr.h>
include "../lib/daophotdef.h"
include "../lib/apseldef.h"
include "../lib/psfdef.h"
include "../lib/peakdef.h"

# DP_FITPSF -- Compute the PSF function.

int procedure dp_fitpsf (dao, im, errmsg, maxch)

pointer	dao				# pointer to the daophot structure
pointer	im				# pointer to the input image
char	errmsg[ARB]			# string containing error message
int	maxch				# max characters in errmsg

int	nfunc, func
pointer	sp, flist, fstr, psf, psffit
int	dp_pstati(), dp_unsatstar(), dp_fitana(), dp_ifitana()
int	dp_fitlt(), dp_fctdecode(), dp_strwrd(), strdic()

begin
	# Get some daophot pointers.
	psf = DP_PSF(dao)
	psffit = DP_PSFFIT(dao)

	# Test to see if there are any psf stars.
	if (dp_pstati (dao, PNUM) <= 0) {
	    call sprintf (errmsg, maxch, "The PSF star list is empty")
	    return (ERR)
	}

	# Test to see if there are any unsaturated stars. At the same time
	# make sure that the first star is unsaturated.
	if (dp_unsatstar (dao) <= 0) {
	    call sprintf (errmsg, maxch,
	        "There are no unsaturated PSF stars")
	    return (ERR)
	}

	# Determine the analytic function.
	call smark (sp)
	call salloc (flist, SZ_FNAME, TY_CHAR)
	call salloc (fstr, SZ_FNAME, TY_CHAR)
	call dp_stats (dao, FUNCLIST, Memc[flist], SZ_FNAME)
	nfunc = dp_fctdecode (Memc[flist], Memc[fstr], SZ_FNAME)
	func = dp_strwrd (1, Memc[fstr], SZ_FNAME, Memc[flist])
	func = strdic (Memc[fstr], Memc[fstr], SZ_LINE, FCTN_FTYPES)

	# Compute the analytic part of the PSF function.
	if (nfunc > 1 || func == FCTN_AUTO) {

	    # Loop over all the analytic functions.
	    if (func == FCTN_AUTO) {
		call strcpy (FCTN_FTYPES, Memc[flist], SZ_FNAME)
		nfunc = FCTN_NFTYPES
	    }

	    # Find the best fitting analytic function.
	    if (dp_ifitana (dao, im, Memc[flist], nfunc) == ERR) {
	        call sprintf (errmsg, maxch,
	            "Analytic function solution failed to converge")
		call sfree (sp)
	        return (ERR)
	    } else if (DP_VERBOSE(dao) == YES)
	        call dp_listpars (dao)

	} else {

	    # Save the analytic function for the new fit.
	    call strcpy (Memc[fstr], DP_FUNCTION(dao), SZ_LINE)

	    # Initialize the parameters.
	    call dp_reinit (dao)

	    # Fit the analytic part of the function.
	    if (dp_fitana (dao, im, Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)],
	        Memr[DP_PH(psf)], Memi[DP_PSAT(psf)], DP_PNUM(psf)) == ERR) {
	        call sprintf (errmsg, maxch,
	            "Analytic function solution failed to converge")
		call sfree (sp)
	        return (ERR)
	    } else if (DP_VERBOSE(dao) == YES) {
	        call printf ("\nFitting function %s    norm scatter: %g\n")
		    call pargstr (DP_FUNCTION(dao))
		    call pargr (DP_PSIGANA(psf))
	        call dp_listpars (dao)
	    }
	}
	call sfree (sp)

	# Compute the look-up table.
	if (dp_fitlt (dao, im) == ERR) {
	    call sprintf (errmsg, maxch,
	        "Too few stars to compute PSF lookup tables")
	    return (ERR)
	} else if (DP_VERBOSE(dao) == YES) {
	    call printf ("\nComputed %d lookup table(s)\n")
	        call pargi (DP_NVLTABLE(psffit)+DP_NFEXTABLE(psffit))
	}

	return (OK)
end


# DP_UNSATSTAR -- Make sure there is at least one unsaturated star.

int procedure dp_unsatstar (dao)

pointer	dao			# pointer to the daophot structure

int	i, first_unsat, nstar
pointer	psf

begin
	psf = DP_PSF(dao)
	first_unsat = 0

	nstar = 0
	do i = 1, DP_PNUM(psf) {
	    if (Memi[DP_PSAT(psf)+i-1] == YES)
		next
	    nstar = nstar + 1
	    if (first_unsat == 0)
		first_unsat = i
	}

	if (first_unsat > 1)
	    call dp_pfswap (dao, 1, first_unsat)

	return (nstar)
end


# DP_REINIT -- Reinitialize the psf function parameters.

procedure dp_reinit (dao)

pointer	dao			# pointer to the daophot structure

pointer	psffit
bool	streq()

begin
	psffit = DP_PSFFIT(dao)

	# Define the psf function.
	if (streq (DP_FUNCTION(dao), "gauss"))
	    DP_PSFUNCTION(psffit) = FCTN_GAUSS
	else if (streq (DP_FUNCTION(dao), "moffat25"))
	    DP_PSFUNCTION(psffit) = FCTN_MOFFAT25
	else if (streq (DP_FUNCTION(dao), "moffat15"))
	    DP_PSFUNCTION(psffit) = FCTN_MOFFAT15
	else if (streq (DP_FUNCTION(dao), "penny1"))
	    DP_PSFUNCTION(psffit) = FCTN_PENNY1
	else if (streq (DP_FUNCTION(dao), "penny2"))
	    DP_PSFUNCTION(psffit) = FCTN_PENNY2
	else if (streq (DP_FUNCTION(dao), "lorentz"))
	    DP_PSFUNCTION(psffit) = FCTN_LORENTZ
	else
	    call error (0, "Unknown PSF function\n")

	switch (DP_VARORDER(dao)) {
	case -1:
	    DP_NVLTABLE(psffit) = 0
	case 0:
	    DP_NVLTABLE(psffit) = 1
	case 1:
	    DP_NVLTABLE(psffit) = 3
	case 2:
	    DP_NVLTABLE(psffit) = 6
	}
	if (DP_FEXPAND(dao) == NO)
	    DP_NFEXTABLE(psffit) = 0
	else
	    DP_NFEXTABLE(psffit) = 5

	# Set the initial values of the function parameters.
	switch (DP_PSFUNCTION(psffit)) {
	case FCTN_GAUSS:
	    DP_PSFNPARS(psffit) = 2
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	case FCTN_MOFFAT25:
	    DP_PSFNPARS(psffit) = 3
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+2] = 0.0
	    Memr[DP_PSFPARS(psffit)+3] = 2.5
	case FCTN_MOFFAT15:
	    DP_PSFNPARS(psffit) = 3
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+2] = 0.0
	    Memr[DP_PSFPARS(psffit)+3] = 1.5
	case FCTN_PENNY1:
	    DP_PSFNPARS(psffit) = 4
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+2] = 0.75
	    Memr[DP_PSFPARS(psffit)+3] = 0.0
	case FCTN_PENNY2:
	    DP_PSFNPARS(psffit) = 5
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+2] = 0.75
	    Memr[DP_PSFPARS(psffit)+3] = 0.0
	    Memr[DP_PSFPARS(psffit)+4] = 0.0
	case FCTN_LORENTZ:
	    DP_PSFNPARS(psffit) = 3
	    Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0
	    Memr[DP_PSFPARS(psffit)+2] = 0.0
	default:
	    call error (0, "Unknown PSF function\n")
	}
end


# DP_IFITANA -- Fit the PSF stars to each of the analytic functions in
# turn to determine which one gives the best fit.

int procedure dp_ifitana (dao, im, funclist, nfunc)

pointer	dao				# pointer to the daophot structure
pointer	im				# pointer to the input image
char	funclist			# the list of functions to be fit
int	nfunc				# number of functions

int	i, psftype, npars
pointer	psf, psffit, sp, fstr, func
pointer	ixtmp, iytmp, ihtmp, istmp, xtmp, ytmp, htmp, ptmp, stmp
real	osig, osum, height, dhdxc, dhdyc, junk, ofactor, factor
int	dp_strwrd(), strdic(), dp_fitana()
real	dp_profile()

begin
	# Get some pointers.
	psf = DP_PSF(dao)
	psffit = DP_PSFFIT(dao)

	# Allocate some temporary storage space.
	call smark (sp)
	call salloc (fstr, SZ_FNAME, TY_CHAR)
	call salloc (func, SZ_FNAME, TY_CHAR)

	call salloc (ixtmp, DP_PNUM(psf), TY_REAL)
	call salloc (iytmp, DP_PNUM(psf), TY_REAL)
	call salloc (ihtmp, DP_PNUM(psf), TY_REAL)
	call salloc (istmp, DP_PNUM(psf), TY_INT)

	call salloc (xtmp, DP_PNUM(psf), TY_REAL)
	call salloc (ytmp, DP_PNUM(psf), TY_REAL)
	call salloc (htmp, DP_PNUM(psf), TY_REAL)
	call salloc (stmp, DP_PNUM(psf), TY_INT)
	call salloc (ptmp, MAX_NFCTNPARS, TY_REAL)

	# Initialize.
	call strcpy (DP_FUNCTION(dao), Memc[func], SZ_FNAME)
	npars = 0
	osig = MAX_REAL
	call amovr (Memr[DP_PXCEN(psf)], Memr[ixtmp], DP_PNUM(psf))
	call amovr (Memr[DP_PYCEN(psf)], Memr[iytmp], DP_PNUM(psf))
	call amovr (Memr[DP_PH(psf)], Memr[ihtmp], DP_PNUM(psf))
	call amovi (Memi[DP_PSAT(psf)], Memi[istmp], DP_PNUM(psf))
	ofactor = dp_profile (DP_PSFUNCTION(psffit), 0.0, 0.0,
	    Memr[DP_PSFPARS(psffit)], dhdxc, dhdyc, junk, 0)

	factor = 1
	do i = 1, nfunc {

	    # Get the function name and set it.
	    if (dp_strwrd (i, Memc[fstr], SZ_FNAME, funclist) <= 0)
		next
	    if (strdic (Memc[fstr], Memc[fstr], SZ_FNAME, FCTN_FTYPES) <= 0)
		next
	    call strcpy (Memc[fstr], DP_FUNCTION(dao), SZ_FNAME)

	    # Start from the same initial state.
	    call dp_reinit (dao)
	    call amovr (Memr[ixtmp], Memr[xtmp], DP_PNUM(psf))
	    call amovr (Memr[iytmp], Memr[ytmp], DP_PNUM(psf))
	    if (i == 1)
	        call amovr (Memr[ihtmp], Memr[htmp], DP_PNUM(psf))
	    else {
		factor = ofactor / dp_profile (DP_PSFUNCTION(psffit), 0.0, 0.0,
	    	    Memr[DP_PSFPARS(psffit)], dhdxc, dhdyc, junk, 0)
		call amulkr (Memr[ihtmp], factor, Memr[htmp], DP_PNUM(psf))
	    }
	    call amovi (Memi[istmp], Memi[stmp], DP_PNUM(psf))

	    call printf ("Trying function %s norm scatter = ")
		call pargstr (Memc[fstr])

	    # Do the fit.
	    if (dp_fitana (dao, im, Memr[xtmp], Memr[ytmp], Memr[htmp],
	        Memi[stmp], DP_PNUM(psf)) == ERR) {
		call printf ("error\n")
		next
	    } else {
		call printf ("%g\n")
		    call pargr (DP_PSIGANA(psf))
	    }

	    # Save the better fit.
	    if (DP_PSIGANA(psf) < osig) {
		call strcpy (Memc[fstr], Memc[func], SZ_FNAME)
		psftype = DP_PSFUNCTION(psffit)
		height = DP_PSFHEIGHT(psffit)
		npars = DP_PSFNPARS(psffit)
		call amovr (Memr[DP_PSFPARS(psffit)], Memr[ptmp],
		    MAX_NFCTNPARS)
	        call amovr (Memr[xtmp], Memr[DP_PXCEN(psf)], DP_PNUM(psf))
	        call amovr (Memr[ytmp], Memr[DP_PYCEN(psf)], DP_PNUM(psf))
	        call amovr (Memr[htmp], Memr[DP_PH(psf)], DP_PNUM(psf))
	        call amovi (Memi[stmp], Memi[DP_PSAT(psf)], DP_PNUM(psf))
		osig = DP_PSIGANA(psf)
		osum = DP_PSUMANA(psf)
	    }
	}

	# Restore the best fit parameters.
	if (npars > 0) {
	    call strcpy (Memc[func], DP_FUNCTION(dao), SZ_FNAME)
	    DP_PSFUNCTION(psffit) = psftype
	    DP_PSFHEIGHT(psffit) = height
	    DP_PSFNPARS(psffit) = npars
	    DP_PSIGANA(psf) = osig
	    DP_PSUMANA(psf) = osum
	    call amovr (Memr[ptmp], Memr[DP_PSFPARS(psffit)],
	        MAX_NFCTNPARS)
	    call printf ("Best fitting function is %s\n")
		call pargstr (DP_FUNCTION(dao))
	}

	# Cleanup.
	call sfree (sp)

	if (npars > 0)
	    return (OK)
	else
	    return (ERR)
end


# DP_FITANA -- Fit the analytic part of the psf function

int procedure dp_fitana (dao, im, pxcen, pycen, ph, pstat, npsfstars)

pointer	dao				# pointer to the daophot structure
pointer	im				# pointer to the input image
real	pxcen[ARB]			# x coordinates of the psf stars
real	pycen[ARB]			# y coordinates of the psf stars
real	ph[ARB]				# heights of the psf stars
int	pstat[ARB]			# saturation status of psf stars
int	npsfstars			# the number of psf stars

int	i, niter, istar, mpar, lx, ly, nx, ny, ier
pointer	apsel, psf, psffit, data
real	fitrad, rsq, oldchi, sumfree
pointer	imgs2r()

begin
	# Get the psf fitting structure pointer.
	apsel = DP_APSEL(dao)
	psf = DP_PSF(dao)
	psffit = DP_PSFFIT(dao)

	# Define some variables.
	oldchi = 0.0
	mpar = 2
	fitrad = DP_FITRAD(dao)
	rsq = fitrad ** 2

	# Get some memory.
	call dp_amempsf (dao)

	# Initialize the fit.
	call amovkr (0.5, Memr[DP_PCLAMP(psf)], DP_PSFNPARS(psffit))
	call aclrr (Memr[DP_PZ(psf)], DP_PSFNPARS(psffit))
	call aclrr (Memr[DP_POLD(psf)], DP_PSFNPARS(psffit))
	Memr[DP_PCLAMP(psf)] = 2.0
	Memr[DP_PCLAMP(psf)+1] = 2.0

	call aclrr (Memr[DP_PXOLD(psf)], npsfstars)
	call aclrr (Memr[DP_PYOLD(psf)], npsfstars)
	call amovkr (1.0, Memr[DP_PXCLAMP(psf)], npsfstars)
	call amovkr (1.0, Memr[DP_PYCLAMP(psf)], npsfstars)

	# Iterate.
	do niter = 1, MAX_NPSFITER {

	    # Initialize the current integration.
	    call aclrr (Memr[DP_PV(psf)], DP_PSFNPARS(psffit)) 
	    call aclrr (Memr[DP_PC(psf)], DP_PSFNPARS(psffit) *
	        DP_PSFNPARS(psffit)) 

	    # Loop over the stars.
	    DP_PSIGANA(psf) = 0.0
	    DP_PSUMANA(psf) = 0.0
	    do istar = 1, npsfstars {

		# Test for saturation.
		if (pstat[istar] == YES)
		    next

		# Define the subraster to be read in.
		lx = int (pxcen[istar] - fitrad) + 1
		ly = int (pycen[istar] - fitrad) + 1
		nx = (int (pxcen[istar] + fitrad) - lx) + 1
		ny = (int (pycen[istar] + fitrad) - ly) + 1

		# Is the star off the image?
		if (lx > IM_LEN(im,1) || ly > IM_LEN(im,2) || (lx + nx - 1) <
		    1 || (ly + ny - 1) < 1) {
		    if (DP_VERBOSE(dao) == YES) {
		        call printf ("Star %d is outside the image\n")
			    call pargi (Memi[DP_APID(apsel)+istar-1])
		    }
		    next
		}

		# Is the star too near the edge of the frame?
		if (lx < 1 || ly < 1 || (lx + nx - 1) > IM_LEN(im,1) ||
		    (ly + ny - 1) > IM_LEN(im,2)) {
		    if (DP_VERBOSE(dao) == YES) {
		        call printf (
			    "Star %d is too near the edge of the image\n")
			    call pargi (Memi[DP_APID(apsel)+istar-1])
		    }
		    next
		}

		# Read in the subraster.
		data = imgs2r (im, lx, lx + nx - 1, ly, ly + ny - 1)

		# Fit x, y, and height for the PSF star istar.
		call dp_xyhiter (DP_PSFUNCTION(psffit),
		    Memr[DP_PSFPARS(psffit)], rsq, Memr[data], nx, ny, lx, ly,
		    pxcen[istar], pycen[istar],
		    Memr[DP_APMSKY(apsel)+istar-1], ph[istar],
		    Memr[DP_PXCLAMP(psf)+istar-1],
		    Memr[DP_PYCLAMP(psf)+istar-1], Memr[DP_PXOLD(psf)+istar-1],
		    Memr[DP_PYOLD(psf)+istar-1])

		# Fit the parameters for the entire list of stars
		call dp_paccum (DP_PSFUNCTION(psffit),
		    Memr[DP_PSFPARS(psffit)], DP_PSFNPARS(psffit), mpar, rsq,
		    Memr[data], nx, ny, lx, ly, pxcen[istar],
		    pycen[istar], Memr[DP_APMSKY(apsel)+istar-1],
		    ph[istar], niter, Memr[DP_PC(psf)],
		    Memr[DP_PV(psf)], Memr[DP_PTMP(psf)], DP_PSIGANA(psf),
		    DP_PSUMANA(psf))
	    }

	    # Invert the matrix and compute the new parameters.
	    call invers (Memr[DP_PC(psf)], DP_PSFNPARS(psffit), mpar, ier)
	    call mvmul (Memr[DP_PC(psf)], DP_PSFNPARS(psffit), mpar,
	        Memr[DP_PV(psf)], Memr[DP_PZ(psf)])

	    do  i = 1, mpar {
	        if ((Memr[DP_PZ(psf)+i-1] * Memr[DP_POLD(psf)+i-1]) < 0.0)
		    Memr[DP_PCLAMP(psf)+i-1] = 0.5 *
		        Memr[DP_PCLAMP(psf)+i-1]
		else
		    Memr[DP_PCLAMP(psf)+i-1] = 1.1 *
		        Memr[DP_PCLAMP(psf)+i-1]
	    }
	    call amovr (Memr[DP_PZ(psf)], Memr[DP_POLD(psf)], mpar)
	    call amulr (Memr[DP_PZ(psf)], Memr[DP_PCLAMP(psf)],
		Memr[DP_PZ(psf)], mpar)

	    Memr[DP_PZ(psf)] = max (-0.1 * Memr[DP_PSFPARS(psffit)],
		min (0.1 * Memr[DP_PSFPARS(psffit)], Memr[DP_PZ(psf)]))
	    Memr[DP_PZ(psf)+1] = max (-0.1 * Memr[DP_PSFPARS(psffit)+1],
		min (0.1 * Memr[DP_PSFPARS(psffit)+1], Memr[DP_PZ(psf)+1]))
	    #if (mpar > 2)
	        #Memr[DP_PZ(psf)+2] = Memr[DP_PZ(psf)+2] /
		    #(1.0 + abs (Memr[DP_PZ(psf)+2]) /
		    #(min (0.1, 1.0 - abs (Memr[DP_PSFPARS(psffit)+2]))))
	    call aaddr (Memr[DP_PSFPARS(psffit)], Memr[DP_PZ(psf)],
		Memr[DP_PSFPARS(psffit)], mpar)

	    # Check for convergence.
	    sumfree = DP_PSUMANA(psf) - real (mpar + 3 * npsfstars)
	    if (sumfree > 0.0 && DP_PSIGANA(psf) >= 0.0)
		DP_PSIGANA(psf) = sqrt (DP_PSIGANA(psf) / sumfree)
	    else
		DP_PSIGANA(psf) = 9.999

	    if (mpar == DP_PSFNPARS(psffit)) {
		if (abs (oldchi / DP_PSIGANA(psf) - 1.0) < 1.0e-5) {
		    DP_PSFHEIGHT(psffit) = ph[1]
		    if (IS_INDEFR(Memr[DP_PMAG(psf)]))
		        DP_PSFMAG(psffit)  = Memr[DP_APMAG(apsel)]
		    else
		        DP_PSFMAG(psffit)  = Memr[DP_PMAG(psf)]
		    DP_PSFX(psffit) = real (IM_LEN(im,1) - 1) / 2.0
		    DP_PSFY(psffit) = real (IM_LEN(im,2) - 1) / 2.0
		    return (OK)
		} else
	    	    oldchi = DP_PSIGANA(psf)
	    } else {
		if (abs (oldchi / DP_PSIGANA(psf) - 1.0) < 1.0e-3) {
		    mpar = mpar + 1
		    oldchi = 0.0
		} else
	    	    oldchi = DP_PSIGANA(psf)
	    }
	}

	return (ERR)
end


# DP_XYHITER -- Increment the initial x, y, and height values for a star.

procedure dp_xyhiter (psftype, params, rsq, data, nx, ny, lx, ly, x, y, sky, h,
	xclamp, yclamp, xold, yold)

int	psftype		# analytic point spread function type
real	params[ARB]	# current function parameter values
real	rsq		# the fitting radius squared
real	data[nx,ARB]	# the input image data
int	nx, ny		# the dimensions of the input image data
int	lx, ly		# the coordinates of the ll corner of the image data
real	x, y		# the input/output stellar coordinates
real	sky		# the input sky value
real	h		# the input/output height value
real	xclamp, yclamp	# the input/output clamping factors for x and y
real	xold, yold	# the input/output x and y correction factors

int	i, j
real	dhn, dhd, dxn, dxd, dyn, dyd, dx, dy, wt, dhdxc, dhdyc, junk, p, dp
real	prod
real	dp_profile()

begin
	dhn = 0.0
	dhd = 0.0
	dxn = 0.0
	dxd = 0.0
	dyn = 0.0
	dyd = 0.0

	do j = 1, ny {
	    dy = real ((ly + j) - 1) - y
	    do i = 1, nx {
	        dx = real ((lx + i) - 1) - x
		wt = (dx ** 2 + dy ** 2) / rsq
		#if (wt >= 1.0)
		if (wt >= 0.999998)
		    next
		p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc, junk, 0)
		dp = data[i,j] - h * p - sky
		dhdxc = dhdxc * h
		dhdyc = dhdyc * h
		wt = 5.0 / (5.0 + (wt / (1.0 - wt)))
		prod = wt * p
		dhn = dhn + prod * dp
		dhd = dhd + prod * p
		prod = wt * dhdxc
		dxn = dxn + prod * dp
		dxd = dxd + prod * dhdxc
		prod = wt * dhdyc
		dyn = dyn + prod * dp
		dyd = dyd + prod * dhdyc
	    }
	}

	h = h + (dhn / dhd)
	dxn = dxn / dxd
	if ((xold * dxn) < 0.0)
	    xclamp = 0.5 * xclamp
	xold = dxn
	x = x + (dxn / (1.0 + (abs(dxn) / xclamp)))
	dyn = dyn / dyd
	if ((yold * dyn) < 0.0)
	    yclamp = 0.5 * yclamp
	yold = dyn
	y = y + (dyn / (1.0 + (abs(dyn) / yclamp)))
end


# DP_PACCUM -- Accumulate the data for the parameter fit.

procedure dp_paccum (psftype, params, npars, mpars, rsq, data, nx, ny, lx,
	ly, x, y, sky, h, iter, c, v, temp, chi, sumwt)

int	psftype		# analytic point spread function type
real	params[ARB]	# current function parameter values
int	npars		# number of function parameters
int	mpars		# the number of active parameters
real	rsq		# the fitting radius squared
real	data[nx,ARB]	# the input image data
int	nx, ny		# the dimensions of the input image data
int	lx, ly		# the coordinates of the ll corner of the image data
real	x, y		# the input/output stellar coordinates
real	sky		# the input sky value
real	h		# the input/output height value
int	iter		# the current iteration
real	c[npars,ARB]	# accumulation matrix
real	v[ARB]		# accumulation vector
real	temp[ARB]	# temporary storage vector
real	chi		# the chi sum
real	sumwt		# the number of points sum

int	i, j, k, l
real	peak, dx, dy, wt, dhdxc, dhdyc, p, dp
real	dp_profile()

begin
	peak = h * dp_profile (psftype, 0.0, 0.0, params, dhdxc, dhdyc, temp, 0)
	do j = 1, ny {
	    dy = real ((ly + j) - 1) - y
	    do i = 1, nx {
	        dx = real ((lx + i) - 1) - x
		wt = (dx ** 2 + dy ** 2) / rsq
		#if (wt >= 1.0)
		if (wt >= 0.999998)
		    next
		p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc, temp, 1) 
		dp = data[i,j] - h * p - sky
		do k = 1, mpars
		    temp[k] = h * temp[k]
		chi = chi + (dp / peak) ** 2
		sumwt = sumwt + 1.0
		wt = 5.0 / (5.0 + (wt / (1.0 - wt)))
		if (iter >= 4)
		    wt = wt / (1.0 + abs (20.0 * dp / peak))
		do k = 1, mpars {
		    v[k] = v[k] + wt * dp * temp[k]
		    do l = 1, mpars {
			c[l,k] = c[l,k] + wt * temp[l] * temp[k]
		    }
		}
	    }
	}
end


# DP_FITLT -- Given the analytic function compute the lookup tables.

int procedure dp_fitlt (dao, im)

pointer	dao		# pointer to the daophot structure
pointer	im		# pointer to the input image

int	istar, nexp, lx, mx, ly, my, iter, nclean, ndata, fit_saturated, nfit
int	nunsat
double	volume
pointer	apsel, psf, psffit, sp, wimname, wim, data
real	datamin, datamax, sumfree, resid, dfdx, dfdy, junk
int	dp_dclean(), dp_resana(), dp_ltcompute(), dp_fsaturated()
double	dp_pnorm()
pointer	immap(), imps2r(), dp_subrast()
real	dp_profile(), dp_sweight()
define	fitsaturated_	11

begin
	# Get some pointers.
	apsel = DP_APSEL(dao)
	psf = DP_PSF(dao)
	psffit = DP_PSFFIT(dao)

	# Check to see whether lookup tables are required.
	nexp = DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit)
	if (nexp <= 0)
	    return (OK)

	# Return if there are too few stars to fit the lookup tables.
	if (DP_PNUM(psf) < nexp)
	    return (ERR)

	# Determine the number of saturated stars.
	nunsat = 0
	do istar = 1, DP_PNUM(psf) {
	    if (Memi[DP_PSAT(psf)+istar-1] == NO)
		next
	    if ((Memr[DP_PH(psf)+istar-1] * dp_profile (DP_PSFUNCTION(psffit),
	        0.0, 0.0, Memr[DP_PSFPARS(psffit)], dfdx, dfdy, junk, 0) +
		Memr[DP_APMSKY(apsel)+istar-1]) <= datamax)
		next
	    Memi[DP_PSAT(psf)+istar-1] = YES
	    Memr[DP_PH(psf)+istar-1] = INDEFR
	    nunsat = nunsat + 1
	}
	nunsat = DP_PNUM(psf) - nunsat

	# Return if there are too few unsaturated psf stars to fit the lookup
	# tables.
	if (nunsat < nexp)
	    return (ERR)

	# Allocate memory for computing lookup tables.
	call dp_tmempsf (dao)

	# Define some constants.
	fit_saturated = DP_SATURATED(dao)
	if (IS_INDEFR(DP_MINGDATA(dao)))
	    datamin = -MAX_REAL
	else
	    datamin = DP_MINGDATA(dao)
	if (IS_INDEFR(DP_MAXGDATA(dao)))
	    datamax = MAX_REAL
	else
	    datamax = DP_MAXGDATA(dao)
	sumfree = sqrt (DP_PSUMANA(psf) / (DP_PSUMANA(psf) - (nexp +
	    3.0 * DP_PNUM(psf))))

	# Get the image name.
	call smark (sp)
	call salloc (wimname, SZ_FNAME, TY_CHAR)
	call mktemp ("tmp", Memc[wimname], SZ_FNAME)

	# Open a temporary image to hold the weights.
	wim = immap (Memc[wimname], NEW_IMAGE, 0)
	IM_NDIM(wim) = 2
	IM_LEN(wim,1) = DP_PNUM(psf)
	IM_LEN(wim,2) = DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit) 
	IM_PIXTYPE(wim) = TY_REAL

	# Compute the constant part of the psf in preparation for normalizing
	# the lookup tables.
	if (nexp > 1)
	    call dp_pconst (DP_PSFUNCTION(psffit), Memr[DP_PSFPARS(psffit)],
	        Memr[DP_PH(psf)], Memr[DP_PCONST(psf)], DP_PSFSIZE(psffit))

	nfit = 0

fitsaturated_

	# Compute the look-up table.
	do iter = 1, DP_NCLEAN(dao) + 1 {

	    # Initialize the fitting arrays.
	    call aclrr (Memr[DP_PV(psf)], nexp)
	    call aclrr (Memr[DP_PC(psf)], nexp * nexp)
	    call aclrr (Memr[DP_PSUMN(psf)], DP_PSFSIZE(psffit) *
	        DP_PSFSIZE(psffit))
	    call aclrr (Memr[DP_PSUMSQ(psf)], DP_PSFSIZE(psffit) *
	        DP_PSFSIZE(psffit))
	    call aclrr (Memr[DP_PSFLUT(psffit)], nexp * DP_PSFSIZE(psffit) *
	        DP_PSFSIZE(psffit))

	    # Loop over the PSF stars.
	    do istar = 1, DP_PNUM(psf) {

		# Get the weight image.
		DP_PSUMW(psf) = imps2r (wim, istar, istar, 1,
		    DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit)) 
	        call aclrr (Memr[DP_PSUMW(psf)], DP_PSFSIZE(psffit) *
	            DP_PSFSIZE(psffit))

		# Skip saturated star?
		if (IS_INDEFR(Memr[DP_PH(psf)+istar-1]))
		    next

	        # Get the data.
	        data = dp_subrast (im, Memr[DP_PXCEN(psf)+istar-1],
		    Memr[DP_PYCEN(psf)+istar-1], DP_PSFRAD(dao), lx, mx,
		    ly, my)

	        # Is the star off the image?
	        if (lx > IM_LEN(im,1) || ly > IM_LEN(im,2) || mx < 1 ||
		    my < 1) {
		    if (DP_VERBOSE(dao) == YES) {
		        call printf ("Star %d is outside the image\n")
		            call pargi (Memi[DP_APID(apsel)+istar-1])
		    }
		    next
	        }

		# Clean bad pixels outside the fitting radius but inside
		# the psf radius from the subraster.
		nclean = dp_dclean (Memr[data], (mx - lx + 1), (my - ly + 1),
		    lx, ly, Memr[DP_PXCEN(psf)+istar-1], Memr[DP_PYCEN(psf)+
		    istar-1], DP_FITRAD(dao), datamin, datamax)

	        # Subtract the analytic part of the fit from the data
		# and compute the normalized residual of the star.
	        ndata = dp_resana (im, DP_PSFUNCTION(psffit),
		    Memr[DP_PSFPARS(psffit)], Memr[data], (mx - lx + 1),
		    (my - ly + 1), lx, ly, Memr[DP_PXCEN(psf)],
		    Memr[DP_PYCEN(psf)], Memr[DP_APMSKY(apsel)],
		    Memr[DP_PH(psf)], DP_PNUM(psf), istar, DP_PSFRAD(dao),
		    DP_FITRAD(dao), datamax, resid)

		# Compute the proper weight for the star.
		if (IS_INDEFR(Memr[DP_PH(psf)+istar-1]) ||
		    Memr[DP_PH(psf)+istar-1] <= 0.0) {
		    Memr[DP_PWEIGHT(psf)+istar-1] = 0.0
		} else if (Memi[DP_PSAT(psf)+istar-1] == YES) {
		    Memr[DP_PWEIGHT(psf)+istar-1] = 0.5 *
		        Memr[DP_PH(psf)+istar-1] / Memr[DP_PH(psf)]
		} else {
		    Memr[DP_PWEIGHT(psf)+istar-1] = (Memr[DP_PH(psf)+
			istar-1] / Memr[DP_PH(psf)]) * dp_sweight (resid,
			sumfree, DP_PSIGANA(psf))
		}

	        # Compute the expansion vector.
	        call dp_eaccum (Memr[DP_PXCEN(psf)+istar-1],
		    Memr[DP_PYCEN(psf)+istar-1], DP_PSFX(psffit),
		    DP_PSFY(psffit), Memr[DP_PTMP(psf)], DP_NVLTABLE(psffit),
		    DP_NFEXTABLE(psffit))

	        # Compute the contribution to the lookup table of the
		# particular star.
	        call dp_ltinterp (Memr[data], (mx - lx + 1), (my - ly + 1),
	            lx, ly, Memr[DP_PXCEN(psf)+istar-1], Memr[DP_PYCEN(psf)+
		    istar-1], Memr[DP_APMSKY(apsel)+istar-1],
		    Memr[DP_PH(psf)+istar-1] / Memr[DP_PH(psf)],
		    Memr[DP_PWEIGHT(psf)+istar-1], Memr[DP_PSUMN(psf)],
		    Memr[DP_PSUMW(psf)], Memr[DP_PSUMSQ(psf)],
		    Memr[DP_PSIGMA(psf)], Memr[DP_POLDLUT(psf)],
		    Memr[DP_PTMP(psf)], Memr[DP_PSFLUT(psffit)], nexp,
		    DP_PSFSIZE(psffit), datamax, iter, DP_NCLEAN(dao) + 1)

	        call mfree (data, TY_REAL)
	    }

	    call imflush (wim)

	    # Compute the lookup table.
	    if (dp_ltcompute (Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)],
	        DP_PNUM(psf), DP_PSFX(psffit), DP_PSFY(psffit),
		Memr[DP_PSUMN(psf)], wim, Memr[DP_PSFLUT(psffit)],
	        Memr[DP_PC(psf)], Memr[DP_PTMP(psf)], Memr[DP_PV(psf)],
		nexp, DP_PSFSIZE(psffit), DP_NVLTABLE(psffit),
		DP_NFEXTABLE(psffit)) == ERR) {
		call sfree (sp)
		return (ERR)
	    }


	    # Compute the standard deviation arrays for the next pass.
	    if (iter < (DP_NCLEAN(dao) + 1)) {
	        if (nunsat <= nexp)
		    break
		call amovr (Memr[DP_PSFLUT(psffit)], Memr[DP_POLDLUT(psf)],
		    nexp * DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit))
		call dp_stdcompute (Memr[DP_PSUMN(psf)], Memr[DP_PSUMSQ(psf)],
		    Memr[DP_PSIGMA(psf)], DP_PSFSIZE(psffit),
		    DP_PSFSIZE(psffit), nexp)
	    }

	}

	if (nexp > 1) {

	    # Accumulate the v vector.
	    call dp_vaccum (Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], 
	        Memr[DP_PH(psf)], Memr[DP_PWEIGHT(psf)], DP_PNUM(psf),
	        DP_PSFX(psffit), DP_PSFY(psffit), Memr[DP_PTMP(psf)],
		Memr[DP_PV(psf)], nexp, DP_NVLTABLE(psffit),
		DP_NFEXTABLE(psffit))

	    # Compute the constant part of the psf.
	    volume = dp_pnorm (Memr[DP_PCONST(psf)], Memr[DP_PSIGMA(psf)],
	        DP_PSFSIZE(psffit))

	    # Normalize lookup tables.
	    call dp_ltnorm (Memr[DP_PCONST(psf)], Memr[DP_PV(psf)],
	        Memr[DP_PSFLUT(psffit)], Memr[DP_PSIGMA(psf)], nexp,
		DP_PSFSIZE(psffit), volume)

	}

	# Make a copy of the psf.
	call dp_pcopy (Memr[DP_PSFLUT(psffit)], Memr[DP_POLDLUT(psf)],
	    DP_PSFSIZE(psffit), DP_PSFSIZE(psffit), nexp)

	# Include the saturated psf stars in the fit.
	if (fit_saturated == YES) {
	    nfit = dp_fsaturated (dao, im, Memi[DP_APID(apsel)],
	        Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], Memr[DP_PH(psf)],
		Memr[DP_APMSKY(apsel)], Memi[DP_PSAT(psf)], DP_PNUM(psf))
	    fit_saturated = NO
	    if (nfit > 0) {
		nunsat = nunsat + nfit
		if (nexp > 1)
		    call aaddr (Memr[DP_PCONST(psf)], Memr[DP_POLDLUT(psf)],
		        Memr[DP_PCONST(psf)], DP_PSFSIZE(psffit) *
		        DP_PSFSIZE(psffit))
	        goto fitsaturated_
	    }
	}

	# Cleanup the temporary images and arrays.
	call imunmap (wim)
	call imdelete (Memc[wimname])
	call sfree (sp)

	return (OK)
end


# DP_DCLEAN -- Clean bad pixels that are outside the fitting radius from
# the data. Note that the star must not be considered to be saturated to
# arrive at this point.

int procedure dp_dclean (data, nx, ny, lx, ly, x, y, fitrad, datamin, datamax)

real	data[nx,ARB]		# the input image data
int	nx, ny			# the dimensions of the input image data
int	lx, ly			# the coordinates of the ll image corner
real	x, y			# the input/output stellar coordinates
real	fitrad			# the fitting radius
real	datamin			# the min good data value
real	datamax			# the max good data value

bool	redo
int	i, j, l, k, nclean
real	frad2, dy2, dr2, sumd, sumn

begin
	nclean = 0
	repeat {

	    redo = false
	    frad2 = fitrad ** 2

	    do j = 1, ny {

	        dy2 = (real ((ly - 1) + j) - y) ** 2
	        if (dy2 < frad2)
		    next
	        do i = 1, nx {

		    if (data[i,j] >= datamin && data[i,j] <= datamax)
		        next
	            dr2 = dy2 + (real ((lx - 1) + i) - x) ** 2
	            if (dr2 < frad2)
		        next

		    sumd = 0.0
		    sumn = 0.0
		    do l = max (1, j-1), min (ny, j+2) {
		        do k = max (1, i-1), min (nx, i+2) {
			    if (data[k,l] < datamin || data[k,l] > datamax)
			        next
			    sumd = sumd + data[k,l]
			    sumn = sumn + 1.0
		        }
		    }
		    if (sumn < 2.5)
		        redo = true
		    else {
			nclean = nclean + 1
		        data[i,j] = sumd / sumn 
		    }
	        }
	    }

	} until (! redo)

	return (nclean)
end


# DP_RESANA -- Compute the residuals from the analytic function.

int procedure dp_resana (im, psftype, params, data, nx, ny, lx, ly,
	x, y, sky, h, nstar, psfstar, psfrad, fitrad, maxgdata, resid)

pointer	im			# the input image descriptor
int	psftype			# analytic point spread function type
real	params[ARB]		# current function parameter values
real	data[nx,ARB]		# the input image data
int	nx, ny			# the dimensions of the input image data
int	lx, ly			# the coordinates of the ll image corner
real	x[ARB]			# the input x coords of the psf stars
real	y[ARB]			# the input y coords of the psf stars
real	sky[ARB]		# the input sky values of the psf stars
real	h[ARB]			# the input height values of the psf stars
int	nstar			# the number of psf stars
int	psfstar			# the psf star in question
real	psfrad			# the psf radius
real	fitrad			# the fitting radius
real	maxgdata		# the maximum good data value
real	resid			# standard deviation of fit

int	i, j, istar, rx, ry, x1, x2, y1, y2, nresid
real	frad2, dx, dy, dy2, dr2, p, dhdxc, dhdyc, junk
int	dp_lsubrast()
real	dp_profile()

begin
	frad2 = fitrad ** 2
	rx = lx + nx - 1
	ry = ly + ny - 1

	resid = 0.0
	nresid = 0
	do istar = 1, nstar {

	    # Check for saturation.
	    if (IS_INDEFR(h[istar]))
		next

	    # Does the subraster of another PSF star overlap the current
	    # subraster ?.
	    if (dp_lsubrast (im, x[istar], y[istar], psfrad, x1, x2,
	        y1, y2) == ERR)
		next
	    if (x2 < lx || y2 < ly || x1 > rx || y1 > ry)
		next

	    # Check the limits of overlap.
	    if (x1 < lx)
		x1 = lx
	    if (x2 > rx)
		x2 = rx
	    if (y1 < ly)
		y1 = ly
	    if (y2 > ry)
		y2 = ry

	    # Subract off the analytic part of the fits and accumulate
	    # the residuals for the psf star.
	    do j = y1 - ly + 1, y2 - ly + 1 {
	        dy = real ((ly - 1) + j) - y[istar]
	        dy2 = dy ** 2
	        do i = x1 - lx + 1, x2 - lx + 1 {
		    if (data[i,j] > maxgdata)
			next
	            dx = real ((lx - 1) + i) - x[istar]
		    p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc,
		        junk, 0)
		    data[i,j] = data[i,j] - h[istar] * p
		    if (istar != psfstar)
		        next
		    dr2 = dy2 + dx ** 2
		    if (dr2 >= frad2)
		        next
		    resid = resid + (data[i,j] - sky[istar]) ** 2
		    nresid = nresid + 1
	        }
	    }
	}

	if (nresid <= 0)
	    resid = 0.0
	else
	    resid = sqrt (resid / nresid) / (h[psfstar] * dp_profile (psftype,
	        0.0, 0.0, params, dhdxc, dhdyc, junk, 0))

	return (nresid)
end


# DP_SWEIGHT -- Compute the weight for the star.

real procedure dp_sweight (resid, sumfree, sumana)

real	resid		# normalized residual wrt analytic fit
real	sumfree		# number of degrees of freedom
real	sumana		# number of points contributing to analytic fit

real	weight

begin
	weight = resid * sumfree / sumana
	weight = 1.0 / (1.0 + (weight / 2.0) ** 2)
	return (weight)
end


# DP_EACCUM -- Calcuate the expansion vector for a single PSF star.

procedure dp_eaccum (x, y, xmid, ymid, junk, nvexp, nfexp)

real	x, y			# the stellar coordinates
real	xmid, ymid		# the psf coordinates
real	junk[ARB]		# temporary storage vector
int	nvexp			# the number of variable psf look-up tables
int	nfexp			# the number of pixel expansion tables

int	j

begin
	# The variable psf terms.
	switch (nvexp) {
	case 1:
	    junk[1] = 1.0
	case 3:
	    junk[1] = 1.0
	    junk[2] = ((x - 1.0) / xmid) - 1.0
	    junk[3] = ((y - 1.0) / ymid) - 1.0
	case 6:
	    junk[1] = 1.0
	    junk[2] = ((x - 1.0) / xmid) - 1.0
	    junk[3] = ((y - 1.0) / ymid) - 1.0
	    junk[4] = (1.5 * (junk[2] ** 2)) - 0.5
	    junk[5] = junk[2] * junk[3]
	    junk[6] = (1.5 * (junk[3] ** 2)) - 0.5
	}

	# The fractional pixel expansion terms if any.
	if (nfexp > 0) {
	    j = nvexp + 1
	    junk[j] = 2.0 * (x - real (nint(x)))
	    j = j + 1
	    junk[j] = 2.0 * (y - real (nint(y)))
	    j = j + 1
	    junk[j] = (1.5 * (junk[j-2] ** 2)) - 0.5
	    j = j + 1
	    junk[j] = junk[j-3] * junk[j-2]
	    j = j + 1
	    junk[j] = (1.5 * (junk[j-3] ** 2)) - 0.5
	}

end


# DP_LTINTERP -- Compute the contribution to the lookup table of a single
# PSF star.

procedure dp_ltinterp (data, nx, ny, lx, ly, x, y, sky, hratio, weight, sumn,
	sumw, sumsq, sig, old, temp, psflut, nexp, nxlut, maxgdata, iter, niter)

real	data[nx,ARB]		# the input image data
int	nx, ny			# the dimensions of the input image data
int	lx, ly			# the coordinates of the ll image corner
real	x, y			# the input/output stellar coordinates
real	sky			# sky value for star
real	hratio			# scale factor for star
real	weight			# weight for the star
real	sumn[nxlut,ARB]		# number of points
real	sumw[nxlut,ARB]		# sum of the weights
real	sumsq[nxlut,ARB]	# sum of the residuals
real	sig[nxlut,ARB]		# residuals of previous iteration
real	old[nexp,nxlut,ARB]	# old lookup table
real	temp[ARB]		# the single star expansion vector
real	psflut[nexp,nxlut,ARB]	# the psf lookup table
int	nexp, nxlut		# the dimensions of the lookup table
real	maxgdata		# maximum good data value
int	iter			# the current iteration
int	niter			# the maximum number of iterations

bool	omit
int	i, j, k, kx, ky, ii, jj, middle, jysq, irsq, midsq
real	jy, ix, dx, dy, diff, dfdx, dfdy, wt, oldsum
real	bicubic()

begin
	middle = (nxlut + 1) / 2
	midsq = middle ** 2
	do j = 1, nxlut {
	    jysq = (j - middle) ** 2
	    if (jysq > midsq)
		next
	    jy = y + real (j - (nxlut + 1) / 2) / 2.0 - real (ly - 1)
	    ky = int (jy)
	    if (ky < 2 || (ky + 2) > ny) 
		next
	    dy = jy - real (ky)
	    do i = 1, nxlut {
	        irsq = jysq + (i - middle) ** 2
		if (irsq > midsq)
		    next
		ix = x + real (i - (nxlut + 1) / 2) / 2.0 - real (lx - 1)
		kx = int (ix)
	        if (kx < 2 || (kx + 2) > nx) 
		    next
		omit = false
		do jj = ky - 1, ky + 2 {
		    do ii = kx - 1, kx + 2 {
			if (data[ii,jj] <= maxgdata)
			    next
			omit = true
		    }
		}
		if (omit)
		    next
		dx = ix - real (kx)
		diff = (bicubic (data[kx-1,ky-1], nx, dx, dy, dfdx, dfdy) -
		    sky) / hratio
		if (iter == 1 || sig[i,j] <= 0.0) {
		    wt = 1.0
		} else {
		    oldsum = 0.0
		    do k = 1, nexp
		        oldsum = oldsum + old[k,i,j]  * temp[k]
		    if ((iter - 1) <= max (3, (niter  - 1) / 2))
		        wt = 1.0 / (1.0 + abs (diff - oldsum) / sig[i,j])
		    else
		        wt = 1.0 / (1.0 + ((diff - oldsum) / sig[i,j]) ** 2)
		}
		wt = wt * weight
		sumn[i,j] = sumn[i,j] + 1.0
		sumw[i,j] = wt
		sumsq[i,j] = sumsq[i,j] + abs (diff)
		psflut[1,i,j] = psflut[1,i,j] + wt * diff
		if (nexp <= 1)
		    next
		do k = 2, nexp
		    psflut[k,i,j] = psflut[k,i,j] + (temp[k] * wt * diff)
	    }
	}
end


# DP_LTCOMPUTE -- Compute the lookup table.

int procedure dp_ltcompute (x, y, npsf, xmid, ymid, sumn, wim, psflut, c,
        junk, v, nexp, nxlut, nvexp, nfexp)

real	x[ARB]			# array of psf star x coordinates
real	y[ARB]			# array of psf star y coordinates
int	npsf			# the number of psf stars
real	xmid, ymid		# the mid-point of the psf
real	sumn[nxlut,ARB]		# number of points
pointer	wim			# pointer to the temporary weight image
real	psflut[nexp,nxlut,ARB]	# the psf lookup table
real	c[nexp,ARB]		# the expansion matrix
real	junk[ARB]		# temporary junk vector
real	v[ARB]			# temporary vector
int	nexp, nxlut		# the size of the lookup table
int	nvexp, nfexp		# size of the expansion look-up tables

int	i, j, k, l, line, istar, ier, middle, midsq, jysq, irsq
pointer	sumw
real	weight
pointer	imgs2r()

begin
	middle = (nxlut + 1) / 2
	midsq = middle ** 2 
	do j = 1, nxlut {
	    jysq = (j - middle) ** 2
	    if (jysq > midsq)
		next
	    do i = 1, nxlut {
		irsq = (i - middle) ** 2 + jysq 
		if (irsq > midsq)
		    next
		if (nint (sumn[i,j]) < nexp)
		    return (ERR)
		line = i + (j - 1) * nxlut
		sumw = imgs2r (wim, 1, npsf, line, line)
		do k = 1, nexp
		    do l = 1, nexp
			c[k,l] = 0.0
		do istar = 1, npsf {
		    weight = Memr[sumw+istar-1]
		    if (weight <= 0.0)
			next
		    call dp_eaccum (x[istar], y[istar], xmid, ymid, junk,
		        nvexp, nfexp)
		    do k = 1, nexp
		        do l = 1, nexp
			    c[k,l] = c[k,l] + weight * junk[k] * junk[l]
		}
		call invers (c, nexp, nexp, ier)
		call mvmul (c, nexp, nexp, psflut[1,i,j], v)
		do k = 1, nexp
		    psflut[k,i,j] = v[k]
	    }
	}

	return (OK)
end


# DP_STDCOMPUTE -- Estimate the standard deviation of the fit.

procedure dp_stdcompute (sumn, sumsq, sigma, nx, ny, nexp)

real	sumn[nx,ARB]		# number of point
real	sumsq[nx,ARB]		# sum of the standard deviations
real	sigma[nx,ARB]		# output standard deviation array
int	nx, ny			# size of the arrays
int	nexp			# number of expansion vectors

int	i, j

begin
	do j = 1, ny {
	    do i = 1, nx {
		if (sumn[i,j] <= nexp)
		    sigma[i,j] = 0.0
		else
		    sigma[i,j] = 1.2533 * sumsq[i,j] / sqrt (sumn[i,j] *
		        (sumn[i,j] - nexp))
	    }
	}
end


# DP_VACCUM -- Accumulate the expansion vector.

procedure dp_vaccum (x, y, h, weight, nstar, xmid, ymid, junk, avetrm,
	nexp, nvexp, nfexp)

real	x[ARB]			# the stellar x coordinates
real	y[ARB]			# the stellar y coordinates
real	h[ARB]			# the stellar heights
real	weight[ARB]		# the stellar weights
int	nstar			# number of stars
real	xmid, ymid		# the psf coordinates
real	junk[ARB]		# temporary storage vector
real	avetrm[ARB]		# the expansion vector
int	nexp			# the total number of expansion terms
int	nvexp			# number of variable expansion terms
int	nfexp			# number of fractional pixel expansion terms

int	k, j

begin
	# Zero the accumulation vector
	do j = 1, nexp
	    avetrm[j] = 0.0

	do k = 1, nstar {
	    if (IS_INDEFR(h[k]))
		next

	    # The variable psf expansion terms.
	    switch (nvexp) {
	    case 1:
	        junk[1] = 1.0
	    case 3:
	        junk[1] = 1.0
	        junk[2] = ((x[k] - 1.0) / xmid) - 1.0
	        junk[3] = ((y[k] - 1.0) / ymid) - 1.0
	    case 6:
	        junk[1] = 1.0
	        junk[2] = ((x[k] - 1.0) / xmid) - 1.0
	        junk[3] = ((y[k] - 1.0) / ymid) - 1.0
	        junk[4] = (1.5 * (junk[2] ** 2)) - 0.5
	        junk[5] = junk[2] * junk[3]
	        junk[6] = (1.5 * (junk[3] ** 2)) - 0.5
	    }

	    # The fractional pixel expansion terms if any.
	    if (nfexp > 0) {
	        j = nvexp + 1
	        junk[j] = 2.0 * (x[k] - real (nint(x[k])))
	        j = j + 1
	        junk[j] = 2.0 * (y[k] - real (nint(y[k])))
	        j = j + 1
	        junk[j] = (1.5 * (junk[j-2] ** 2)) - 0.5
	        j = j + 1
	        junk[j] = junk[j-3] * junk[j-2]
	        j = j + 1
	        junk[j] = (1.5 * (junk[j-3] ** 2)) - 0.5
	    }

	    # Accumulate the expansion vector.
	    do j = 1, nexp
	        avetrm[j] = avetrm[j] + weight[k] * junk[j]
	}

	# Average the expansion vector.
	do j = nexp, 1, -1
	    avetrm[j] = avetrm[j] / avetrm[1]
end


# DP_PCONST -- Compute the analytic part of the psf.

procedure dp_pconst (psftype, params, hpsf1, ana, nxlut)

int	psftype			# the type of analytic psf function
real	params[ARB]		# the analytic parameters array
real	hpsf1			# the height of the psf function
real	ana[nxlut,ARB]		# the computed analytic function
int	nxlut			# the size of the constant lookup table

int	i, j, middle, midsq, dj2, dr2
real	dx, dy, dfdx, dfdy, junk
real	dp_profile()

begin
	# Compute the constant part of the psf.
	middle = (nxlut + 1) / 2
	midsq = middle ** 2
	do j = 1, nxlut {
	    dj2 = (j - middle) ** 2
	    if (dj2 > midsq)
		next
	    dy = real (j - middle) / 2.0
	    do i = 1, nxlut {
		dr2 = (i - middle) ** 2 + dj2
		if (dr2 > midsq)
		    next
	        dx = real (i - middle) / 2.0
		ana[i,j] = hpsf1 * dp_profile (psftype, dx, dy, params,
		    dfdx, dfdy, junk, 0)
	    }
	}
end


# DP_PNORM -- Compute the psf normalization parameters.

double procedure dp_pnorm (ana, pixels, nxlut)

real	ana[nxlut,ARB]		# the computed analytic function
real	pixels[ARB]		# pixel storage array
int	nxlut			# the size of the constant lookup table

int	i, j, middle, midsq, edgesq, npts, dj2, dr2
double	vol
real	median
real	pctile()

begin
	# Ensure that the profile which will be added and subtracted
	# from the various lookup tables has a median value of zero 
	# around the rim

	middle = (nxlut + 1) / 2
	midsq = middle ** 2
	edgesq = (middle - 2) ** 2
	npts = 0
	do j = 1, nxlut {
	    dj2 = (j - middle) ** 2
	    if (dj2 > midsq)
		next
	    do i = 1, nxlut {
		dr2 = (i - middle) ** 2 + dj2
	        if (dr2 > midsq)
		    next
		if (dr2 < edgesq)
		    next
		npts = npts + 1
		pixels[npts] = ana[i,j]
	    }
	}
	median = pctile (pixels, npts, (npts+1) / 2)

	vol = 0.0d0
	do j = 1, nxlut {
	    dj2 = (j - middle) ** 2
	    if (dj2 > midsq)
		next
	    do i = 1, nxlut {
		dr2 = (i - middle) ** 2 + dj2
	        if (dr2 > midsq)
		    next
		ana[i,j] = ana[i,j] - median
		vol = vol + double (ana[i,j])
	    }
	}

	return (vol)
end


# DP_LTNORM -- Compute the final lookup table.

procedure dp_ltnorm (ana, v, psflut, pixels, nexp, nxlut, volume)

real	ana[nxlut,ARB]		# analytic part of the profile
real	v[ARB]			# the expansion vector
real	psflut[nexp,nxlut,ARB]	# the psf lookup table
real	pixels[ARB]		# scratch array for determining median
int	nexp, nxlut		# the size of the lookup table
double	volume			# total flux in the constant psf

int	i, j, k, middle, midsq, edgesq, npts, dj2, dr2
double	sum
real	median, dmedian
real	pctile()

begin
	# Ensure that the psf which will be added and subtracted from the
	# various lookup tables has a median value of zero around the rim.

	middle = (nxlut + 1) / 2
	midsq = middle ** 2
	edgesq = (middle - 2) ** 2
	do k = 2, nexp {

	    npts = 0
	    do j = 1, nxlut {
	        dj2 = (j - middle) ** 2
		if (dj2 > midsq)
		    next
	        do i = 1, nxlut {
		    dr2 = (i - middle) ** 2 + dj2
		    if (dr2 > midsq)
			next
		    if (dr2 < edgesq)
		        next
		    npts = npts + 1
		    pixels[npts] = psflut[k,i,j]
	        }
	    }

	    median = pctile (pixels, npts, (npts+1) / 2)
	    dmedian = v[k] * median

	    do j = 1, nxlut {
	        dj2 = (j - middle) ** 2
		if (dj2 > midsq)
		    next
		do i = 1, nxlut {
		    dr2 = (i - middle) ** 2 + dj2
		    if (dr2 > midsq)
			next
		    psflut[k,i,j] = psflut[k,i,j] - median
		    psflut[1,i,j] = psflut[1,i,j] + dmedian
		}
	    }
	}

	# Determine the net volume of each of the higher order PSF 
	# tables and force it to zero by subtracting a scaled copy
	# of the constant psf. Scale the part that has been subtracted
	# off by the mean polynomial term and add it in to the constant
	# part of the psf so that at the centroid of the psf stars
	# positions the psf remains unchanged.

	do k = 2, nexp {

	    sum = 0.0d0
	    do j = 1, nxlut {
	        dj2 = (j - middle) ** 2
		if (dj2 > midsq)
		    next
		do i = 1, nxlut {
		    dr2 = (i - middle) ** 2 + dj2
		    if (dr2 > midsq)
			next
		    sum = sum + double (psflut[k,i,j])
		}
	    }

	    median = real (sum / volume)
	    dmedian = v[k] * median

	    do j = 1, nxlut {
	        dj2 = (j - middle) ** 2
		if (dj2 > midsq)
		    next
		do i = 1, nxlut {
		    dr2 = (i - middle) ** 2 + dj2
		    if (dr2 > midsq)
			next
		    psflut[k,i,j] = psflut[k,i,j] - median * ana[i,j]
		    psflut[1,i,j] = psflut[1,i,j] + dmedian * ana[i,j]
		}
	    }
	}
end


# DP_FSATURATED -- Fit the saturated stars.

int procedure dp_fsaturated (dao, im, id, xcen, ycen, h, sky, sat, nstars)

pointer	dao			# pointer to the main daophot structure
pointer	im			# pointer to the input image
int	id[ARB]			# array of stellar ids
real	xcen[ARB]		# array of stellar y coordinates
real	ycen[ARB]		# array of stellar y coordinates
real	h[ARB]			# array of stellar amplitudes
real	sky[ARB]		# array of sky values
int	sat[ARB]		# array of saturation indicators
int	nstars			# number of stars

int	nfit, nsat, istar, lowx, lowy, nxpix, nypix, recenter, fitsky
int	clipexp, ier, niter
pointer	psf, psffit, subim, psflut
real	x, y, dx, dy, skyval, scale, errmag, chi, sharp, cliprange
int	dp_pkfit()
pointer dp_gsubrast()

begin
	# Get some pointers.
	psf = DP_PSF(dao)
	psffit = DP_PSFFIT(dao)

	# Allocate memory.
	call dp_pksetup (dao)
	call dp_mempk (dao, 3)

	# Save the default values of some critical parameters.
	recenter = DP_RECENTER(dao)
	fitsky = DP_FITSKY(dao)
	psflut = DP_PSFLUT(psffit)
	clipexp = DP_CLIPEXP(dao)
	cliprange = DP_CLIPRANGE(dao)

	# Set some fitting parameters.
	DP_RECENTER(dao) = YES
	DP_FITSKY(dao) = NO
	DP_CLIPEXP(dao) = 8
	DP_CLIPRANGE(dao) = 2.5
	DP_PSFLUT(psffit) = DP_POLDLUT(psf)

	nfit = 0
	nsat = 0
	do istar = 1, nstars {

	    # Skip saturated stars.
	    if (sat[istar] == NO)
		next
	    nsat = nsat + 1

	    # Get the data.
	    subim = dp_gsubrast (im, xcen[istar], ycen[istar], DP_PSFRAD(dao),
	        lowx, lowy, nxpix, nypix)
	    if (subim == NULL)
		next

	    # Set the intial values for the fit parameters.
	    x = xcen[istar] - lowx + 1.0
	    y = ycen[istar] - lowy + 1.0
	    dx = (xcen[istar] - 1.0) / DP_PSFX(psffit) - 1.0
	    dy = (ycen[istar] - 1.0) / DP_PSFY(psffit) - 1.0
	    skyval = sky[istar]
	    scale = 3.0 

	    # Fit the star.
	    ier = dp_pkfit (dao, Memr[subim], nxpix, nypix, 0.5 *
	        DP_PSFRAD(dao), x, y, dx, dy, scale, skyval, errmag, chi,
		sharp, niter)

	    # Compute the fit parameters.
	    if (ier != PKERR_OK) {
		scale = INDEFR
		errmag = INDEFR
		niter = 0
		chi = INDEFR
		sharp = INDEFR
	    } else {
		nfit = nfit + 1
	        xcen[istar] = x + lowx - 1.0
	        ycen[istar] = y + lowy - 1.0
	        h[istar] = scale * h[1]
		errmag = 1.085736 * errmag / scale
		if (errmag >= 2.0)
		    errmag = INDEFR
		scale = DP_PSFMAG(psffit) - 2.5 * log10 (scale)
	    }

	    if (DP_VERBOSE(dao) == YES) {
		if (nsat == 1)
	            call printf ("\nFit for saturated stars\n")
		call printf (
		    "    %6d  %7.2f  %7.2f  %8.3f  %6.3f  %3d  %7.2f  %7.2f\n")
		    call pargi (id[istar])
		    call pargr (xcen[istar])
		    call pargr (ycen[istar])
		    call pargr (scale)
		    call pargr (errmag)
		    call pargi (niter)
		    call pargr (chi)
		    call pargr (sharp)
	    }
	}

	if (DP_VERBOSE(dao) == YES && nsat > 0)
	    call printf ("\n")

	# Restore the default values of some critical parameters.
	DP_RECENTER(dao) = recenter
	DP_FITSKY(dao) = fitsky
	DP_CLIPEXP(dao) = clipexp
	DP_CLIPRANGE(dao) = cliprange
	DP_PSFLUT(psffit) = psflut

	# Free memory.
	call dp_pkclose (dao)

	return (nfit)
end


# DP_PCOPY -- Make a copy of the psf in correct storage order.

procedure dp_pcopy (inlut, outlut, nx, ny, nexp)

real	inlut[nexp,nx,ny]		# the input look-up table
real	outlut[nx,ny,nexp]		# the output look-up table
int	nx,ny,nexp			# the size of the look-up table

int	k,i,j

begin
	do k = 1, nexp {
	    do j = 1, ny {
		do i = 1, nx {
		    outlut[i,j,k] = inlut[k,i,j]
		}
	    }
	}
end