path: root/noao/digiphot/photcal/mkobsfile/phagrow.x

include <time.h>
include <mach.h>
include <fset.h>
include "../lib/apfile.h"


# PH_AGROW - Compute the curve of growth using the data from the input
# photometry files, optional airmasses from the airmass file, and the
# initial values of the parameters.

procedure ph_agrow (apcat, magfd, logfd, mgd, imtable, id, x, y, nap, rap,
	mag, merr, naperts, params, lterms, smallap, largeap)

int	apcat			# the aperture correction file descriptor
int	magfd			# the output magnitude file descriptor
int	logfd			# the log file descriptor
pointer	mgd			# pointer to the plot metacode file
pointer	imtable			# pointer to the image symbol table
int	id[ARB]			# the star ids
real	x[ARB]			# the star x coordinates
real	y[ARB]			# the star y coordinates
int	nap[ARB]		# the array of aperture numbers
real	rap[naperts,ARB]	# input array of aperture radii
real	mag[naperts,ARB]	# input array of magnitudes / differences
real	merr[naperts,ARB]	# input array of magnitude errors
int	naperts			# the number of input apertures
double	params[ARB]		# the initial values of the parameters
int	lterms			# the number of terms to be fit
int	smallap			# the small aperture number
int	largeap			# the large aperture number

int	nimages
pointer	agr
int	stnsymbols()
int	ph_amfit()

begin
	if (logfd != NULL)
	    call ph_logtitle (logfd, apcat)

	nimages = stnsymbols (imtable, 0)
	if (nimages <= 0) {
	    call printf ("Error: There is no input data to fit\n")
	    if (logfd != NULL)
	        call fprintf (logfd, "There is no input data to fit\n")
	    return
	}

	# Allocate memory required for fitting.
	call ph_ameminit (agr, lterms, naperts)

	# Fit the seeing radii and the cog model parameters.
	if (ph_amfit (imtable, agr, nap, rap, mag, merr, naperts, params,
	    lterms) == ERR) {
	    call printf ("Error: The cog model fit did not converge\n")
	    if (logfd != NULL) {
		call ph_ishow (logfd, params, lterms)
	        call ph_pshow (logfd, agr, lterms)
	        call ph_rshow (logfd, imtable)
	        call fprintf (logfd,
		    "Error: The cog model fit did not converge\n")
	    }
	} else {
	    if (logfd != NULL) {
		call ph_ishow (logfd, params, lterms)
	        call ph_pshow (logfd, agr, lterms)
	        call ph_rshow (logfd, imtable)
	    }
	    call ph_aeval (imtable, agr, apcat, magfd, logfd, mgd, id, x, y,
	        nap, rap, mag, merr, naperts, smallap, largeap)
	    if (logfd != NULL)
		call ph_tfshow (logfd, agr, naperts)
	}

	# Free memory required for fitting.
	call ph_amemfree (agr)
end


# PH_LOGTITLE -- Print the title for the new log file entry

procedure ph_logtitle (logfd, apcat)

int	logfd			# the log file descriptor
int	apcat			# the aperture correction file descriptor

pointer	sp, afname, date
long	clktime()

begin
	# Allocate working space.
	call smark (sp)
	call salloc (afname, SZ_FNAME, TY_CHAR)
	call salloc (date, SZ_TIME, TY_CHAR)

	# Get the file name and the times.
	call fstats (apcat, F_FILENAME, Memc[afname], SZ_FNAME)
	call cnvtime (clktime (long(0)), Memc[date], SZ_TIME)

	call fprintf (logfd, "NEW LOGFILE ENTRY AT: %s\n")
	    call pargstr (Memc[date])
	call fprintf (logfd, "APFILE: %s\n\n")
	    call pargstr (Memc[afname])

	call sfree (sp)
end


# PH_AMEMINIT -- Allocate memory for doing the fit.

procedure ph_ameminit (agr, lterms, naperts)

pointer	agr			# pointer to the fitting structure
int	lterms			# the number of terms to fit
int	naperts			# the number of apertures

begin
	call malloc (agr, LEN_AGRSTRUCT, TY_STRUCT)

	call malloc (AGR_DM(agr), naperts, TY_DOUBLE)
	call malloc (AGR_DDDR(agr), naperts, TY_DOUBLE)
	call malloc (AGR_T(agr), lterms * naperts, TY_DOUBLE)
	call malloc (AGR_U(agr), lterms * lterms, TY_DOUBLE)
	call malloc (AGR_V(agr), lterms, TY_DOUBLE)
	call malloc (AGR_POLD(agr), lterms, TY_DOUBLE)
	call malloc (AGR_DP(agr), lterms, TY_DOUBLE)
	call malloc (AGR_PARAMS(agr), MAX_MTERMS, TY_DOUBLE)
	call malloc (AGR_PERRORS(agr), MAX_MTERMS, TY_DOUBLE)
	call malloc (AGR_PCLAMPS(agr), MAX_MTERMS, TY_DOUBLE)

	call malloc (AGR_RBAR(agr), naperts, TY_REAL)
	call malloc (AGR_W(agr), naperts, TY_REAL)
	call malloc (AGR_THEO(agr), naperts, TY_REAL)
	call malloc (AGR_ADOPT(agr), naperts, TY_REAL)
	call malloc (AGR_WOBS(agr), naperts, TY_REAL)
	call malloc (AGR_OBS(agr), naperts, TY_REAL)
	call malloc (AGR_WADO(agr), naperts, TY_REAL)

	call malloc (AGR_CUM(agr), naperts + 1, TY_REAL)
	call malloc (AGR_TCUM(agr), naperts + 1, TY_REAL)
	call malloc (AGR_WCUM(agr), naperts + 1, TY_REAL)
	call malloc (AGR_MAGS(agr), naperts, TY_REAL)
	call malloc (AGR_CMAGS(agr), naperts, TY_REAL)
	call malloc (AGR_TMAGS(agr), naperts, TY_REAL)
	call malloc (AGR_WMAGS(agr), naperts, TY_REAL)

	call calloc (AGR_AVE(agr), naperts, TY_REAL)
	call calloc (AGR_RESID(agr), naperts, TY_REAL)
	call calloc (AGR_RESSQ(agr), naperts, TY_REAL)
	call calloc (AGR_WR(agr), naperts, TY_REAL)
	call calloc (AGR_TAVE(agr), naperts, TY_REAL)
	call calloc (AGR_TRESID(agr), naperts, TY_REAL)
	call calloc (AGR_TRESSQ(agr), naperts, TY_REAL)
	call calloc (AGR_TWR(agr), naperts, TY_REAL)
end


# PH_AMFIT -- Fit the seeing disk widths and the model parameters.

int procedure ph_amfit (imtable, agr, nap, rap, mag, merr, naperts, params,
	lterms)

pointer	imtable			# pointer to the image symbol table
pointer	agr			# pointer to the fitting structure
int	nap[ARB]		# the array of aperture numbers
real	rap[naperts,ARB]	# input array of aperture radii
real	mag[naperts,ARB]	# input array of magnitudes / differences
real	merr[naperts,ARB]	# input array of magnitude errors
int	naperts			# the number of input apertures
double	params[ARB]		# the parameters to be fit
int	lterms			# number of terms to fit

double	pold, sumd, sumw, sumn, sumr, sumy, old
int	niter, i, j, k, nimages, ipow, nterms, iflag, ntimes, ngood
pointer	sp, sym, symbol
real	gain, rold, rclamp, dr
int	stnsymbols()
pointer	sthead(), stnext()

begin
	nimages = stnsymbols (imtable, 0)
	if (nimages <= 0)
	    return (ERR)

	# Allocate temporary space
	call smark (sp)
	call salloc (sym, nimages, TY_POINTER)

	# Order the symbol table.
	symbol = sthead (imtable)
	do k = nimages, 1, -1 {
	    Memi[sym+k-1] = symbol
	    symbol = stnext (imtable, symbol)
	}

	# Initialize some variables
	ipow = 1
	nterms = 0
	gain = 0.0
	pold = 0.0d0
	call aclrd (Memd[AGR_POLD(agr)], lterms)
	call ph_apinit (params, Memd[AGR_PARAMS(agr)], Memd[AGR_PERRORS(agr)],
	    Memd[AGR_PCLAMPS(agr)])

	# Compute an improved value for ro and for the stellar profile
	# model parameters.

	for (niter = 1; niter <= AGR_ITMAX; niter = niter + 1) {

	    sumd = 0.0d0
	    sumw = 0.0d0
	    sumn = 0.0d0

	    # Compute the number of terms to be fit this iteration.
	    call ph_anterms (niter, lterms, gain, nterms)

	    # Zero the accumulation matrix and vector
	    call aclrd (Memd[AGR_U(agr)], lterms * lterms)
	    call aclrd (Memd[AGR_V(agr)], lterms)

	    do k = 1, nimages {

		# Get the current symbol.
		symbol = Memi[sym+k-1]

		# Check to see if there is any data.
		if (IMT_NENTRIES(symbol) <= 0) {
		    IMT_RO(symbol) = INDEFR
		    next
		}

		# Check to see if there is any good data.
		ngood = 0
		do i = 1, IMT_NENTRIES(symbol) {
		    do j = 2, nap[i]
			ngood = ngood + 1
		}
		if (ngood <= 0) {
		    IMT_RO(symbol) = INDEFR
		    next
		}

		# Get better estimage of ro.
		if (niter == 1) {
		    if (k > 1 && ! IS_INDEFR(IMT_RO(Memi[sym+k-2]))) {
		        IMT_RO(symbol) = IMT_RO(Memi[sym+k-2])
		    } else
			IMT_RO(symbol) = 0.5 * rap[1,1]
		}

		# Set the model derivative vector
		call ph_adddr (rap[1,IMT_OFFSET(symbol)], Memd[AGR_DDDR(agr)],
		    naperts, Memd[AGR_PARAMS(agr)], IMT_RO(symbol),
		    IMT_NXAIRMASS(symbol))

		# Get the new estimate of ro
		rold = 0.0
		rclamp = IMT_RO(symbol) / 4.0
		ntimes = 0
		repeat {
		    
		    sumr = 0.0d0
		    sumy = 0.0d0

		    call ph_adddrdm (rap[1,IMT_OFFSET(symbol)], 
		        Memd[AGR_DM(agr)], Memd[AGR_DDDR(agr)], naperts,
			Memd[AGR_PARAMS(agr)], IMT_RO(symbol),
			IMT_NXAIRMASS(symbol), niter)

		    call ph_awsum (nap[IMT_OFFSET(symbol)],
		        mag[1,IMT_OFFSET(symbol)],
		        merr[1,IMT_OFFSET(symbol)], Memd[AGR_DM(agr)],
			Memr[AGR_W(agr)], Memd[AGR_DDDR(agr)], naperts,
			IMT_NENTRIES(symbol), ipow, sumr, sumy)

		    if (sumy > 0.0d0) {
			dr = sumr / sumy
			if (dr * rold < 0.0)
			    rclamp = 0.5 * rclamp
			IMT_RO(symbol) = IMT_RO(symbol) - dr / (1.0 + abs (dr /
			    rclamp))
			if (IMT_RO(symbol) <= EPSILONR) {
			    IMT_RO(symbol) = INDEFR
			    call sfree (sp)
			    return (ERR)
			}
			rold = dr
			if (abs (dr / IMT_RO(symbol)) <= 3.0e-4)
			    break
		    } else {
			IMT_RO(symbol) = INDEFR
		        break
		    }

		    ntimes = ntimes + 1
		    if (ntimes >= 100) {
			IMT_RO(symbol) = INDEFR
			call sfree (sp)
			return (ERR)
		    }
		}

		call ph_adp (rap[1,IMT_OFFSET(symbol)], Memd[AGR_DM(agr)],
		    Memd[AGR_T(agr)], naperts, lterms, nterms,
		    Memd[AGR_PARAMS(agr)], IMT_RO(symbol),
		    IMT_NXAIRMASS(symbol))

		call ph_aaccum (nap[IMT_OFFSET(symbol)],
		    mag[1,IMT_OFFSET(symbol)], merr[1,IMT_OFFSET(symbol)],
		    Memd[AGR_DM(agr)], Memr[AGR_W(agr)], naperts,
		    IMT_NENTRIES(symbol), Memd[AGR_T(agr)], Memd[AGR_U(agr)],
		    Memd[AGR_V(agr)], lterms, nterms, ipow, sumd, sumw, sumn)
	    }

	    # Invert the matrix.
	    call dinver (Memd[AGR_U(agr)], lterms, nterms, iflag)
	    if (iflag != 0) {
		call sfree (sp)
		return (ERR)
	    }

	    # Get new values for the parameters.
	    call dmvmul (Memd[AGR_U(agr)], lterms, nterms, Memd[AGR_V(agr)],
		Memd[AGR_DP(agr)])
	    call ph_apsolve (Memd[AGR_POLD(agr)], Memd[AGR_PARAMS(agr)],
	        Memd[AGR_DP(agr)], Memd[AGR_PCLAMPS(agr)], nterms, ipow)

	    # Test the fit.
	    #if (sumn > 6.0)
	    if (sumn > (nterms + 1.0))
	        #sumn = sqrt (sumd / (sumn - 6.0))
	        sumn = sqrt (sumd / (sumn - (nterms + 1.0)))
	    else
		sumn = 0.0
	    sumd = sqrt (sumd / sumw)
	    gain = 2.0 * abs (old - sumd) / (old + sumd)
	    old = sumd
	    if ((nterms < lterms) || (gain > 0.001)) {
		next
	    } else if (ipow <= 2) {
		ipow = 3
		call amovkd (0.1d0, Memd[AGR_PCLAMPS(agr)], 4)
		next
	    } else
		break
	}

	if (niter > AGR_ITMAX) {
	    call sfree (sp)
	    return (ERR)
	}

	# Compute the errors.
	do k = 1, nterms {
	    if (Memd[AGR_U(agr)+(k-1)*lterms+k-1] > 0.0)
	        Memd[AGR_PERRORS(agr)+k-1] =
	            sumn * sqrt (Memd[AGR_U(agr)+(k-1)*lterms+k-1])
	    else
	        Memd[AGR_PERRORS(agr)+k-1] = 0.0d0
	}


	call sfree (sp)
	return (OK)
end


# PH_ISHOW -- Show the initial values of the parameters.

procedure ph_ishow (fd, params, mterms)

int	fd			# pointer to the output file descriptor
double	params[ARB]		# the parameter values
int	mterms			# the number of terms to be fit

int	i

begin
	# Print out the best fit parameters.
	call fprintf (fd,
	    "\nThe initial cog model parameter values for nparams %d\n")
	    call pargi (mterms)
	do i = 1, MAX_MTERMS {
	    call fprintf (fd, "\t%10s: %15.7g\n")
		switch (i) {
		case 1:
		    call pargstr ("swings")
		case 2:
		    call pargstr ("pwings")
		case 3:
		    call pargstr ("pgauss")
		case 4:
		    call pargstr ("rgescale")
		case 5:
		    call pargstr ("xwings")
		}
		call pargd (params[i])
	}
	call fprintf (fd, "\n")
end


# PH_PSHOW -- Show the results of the parameter fitting.

procedure ph_pshow (fd, agr, mterms)

int	fd			# the output file descriptor
pointer	agr			# the pointer to the fitting descriptor
int	mterms			# the number of parameters to fit

begin
	call ph_pwrite (fd, Memd[AGR_PARAMS(agr)], Memd[AGR_PERRORS(agr)],
	    MAX_MTERMS, mterms)
end


# PH_RSHOW -- Show the computed seeing parameters as a function of image.

procedure ph_rshow (fd, imtable)

int	fd			# the output file descriptor
pointer	imtable			# pointer to the symbol table

int	i, nimages
pointer	sp, sym, symbol
int	stnsymbols()
pointer	sthead(), stnext()

begin
	nimages = stnsymbols (imtable, 0)
	if (nimages <= 0)
	    return

	# Allocate temporary space
	call smark (sp)
	call salloc (sym, nimages, TY_POINTER)

	# Order the symbol table.
	symbol = sthead (imtable)
	do i = nimages, 1, -1 {
	    Memi[sym+i-1] = symbol
	    symbol = stnext (imtable, symbol)
	}

	call fprintf (fd,
	    "\nThe seeing radius and assumed airmass for each image\n")
	do i = 1, nimages {
	    symbol = Memi[sym+i-1]
	    call fprintf (fd, "%30s  %8.4f  %8.3f\n")
		call pargstr (IMT_IMNAME(symbol))
		call pargr (IMT_RO(symbol))
		call pargr (IMT_XAIRMASS(symbol))
	}
	call fprintf (fd, "\n")

	call sfree (sp)
end


# PH_AEVAL -- Evaluate the fit for all the images.

procedure ph_aeval (imtable, agr, apcat, magfd, logfd, mgd, id, x, y, nap,
        rap, mag, merr, naperts, smallap, largeap)

pointer	imtable			# pointer to the image symbol table
pointer	agr			# pointer to the fitting structure
int	apcat			# the aperture correction file descriptor
int	magfd			# the best magnitudes file descriptor
int	logfd			# the log file descriptor
pointer	mgd			# pointer to the plot metacode file
int	id[ARB]			# the star ids
real	x[ARB]			# the star x coordinates
real	y[ARB]			# the star y coordinates
int	nap[ARB]		# the array of aperture numbers
real	rap[naperts,ARB]	# input array of aperture radii
real	mag[naperts,ARB]	# input array of magnitudes / differences
real	merr[naperts,ARB]	# input array of magnitude errors
int	naperts			# the number of input apertures
int	smallap			# the small aperture number
int	largeap			# the large aperture number

int	k, nimages
pointer	sp, sym, symbol
int	stnsymbols()
pointer	sthead(), stnext()

begin
	nimages = stnsymbols (imtable, 0)
	if (nimages <= 0)
	    return

	# Initialize some vectors.
	call aclrr (Memr[AGR_TAVE(agr)], naperts)
	call aclrr (Memr[AGR_TRESID(agr)], naperts)
	call aclrr (Memr[AGR_TRESSQ(agr)], naperts)
	call aclrr (Memr[AGR_TWR(agr)], naperts)

	# Allocate temporary space
	call smark (sp)
	call salloc (sym, nimages, TY_POINTER)

	# Order the symbol table.
	symbol = sthead (imtable)
	do k = nimages, 1, -1 {
	    Memi[sym+k-1] = symbol
	    symbol = stnext (imtable, symbol)
	}

	# Write the banner for the apfile.
	call fprintf (apcat,
        "# The aperture correction from apertures %d to %d in magnitudes\n\n")
	    call pargi (smallap)
	    call pargi (largeap)

	# Write the banner for the magfile.
	if (magfd != NULL) {
	    call fprintf (magfd, "# Magnitudes corrected to aperture %d\n\n")
	        call pargi (largeap)
	    call fprintf (magfd,
"#%19tImage%30tFilter%39tExptime%47tAirmass%62tOtime%70tXcenter%80tYcenter\
%93tMag%101tMerr%107tRadius\n\n")
	}

	# Compute the aperture corrections.
	do k = 1, nimages {

	    symbol = Memi[sym+k-1]

	    # Initialize some vectors.
	    call aclrr (Memr[AGR_AVE(agr)], naperts)
	    call aclrr (Memr[AGR_RESID(agr)], naperts)
	    call aclrr (Memr[AGR_RESSQ(agr)], naperts)
	    call aclrr (Memr[AGR_WR(agr)], naperts)

	    call aclrr (Memr[AGR_ADOPT(agr)], naperts)
	    call aclrr (Memr[AGR_WADO(agr)], naperts)
	    call aclrr (Memr[AGR_WOBS(agr)], naperts)

	    call ph_rinit (rap[1,IMT_OFFSET(symbol)], Memr[AGR_RBAR(agr)],
	        naperts)
	    call ph_tinit (rap[1,IMT_OFFSET(symbol)], Memr[AGR_THEO(agr)],
	        naperts, Memd[AGR_PARAMS(agr)], IMT_RO(symbol),
		IMT_NXAIRMASS(symbol))

	    # Accumulate the difference data.
	    call ph_taccum (nap[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
		merr[1,IMT_OFFSET(symbol)], Memr[AGR_THEO(agr)],
		Memr[AGR_W(agr)], Memr[AGR_ADOPT(agr)],
		Memr[AGR_WOBS(agr)], Memr[AGR_AVE(agr)],
		Memr[AGR_RESID(agr)], Memr[AGR_RESSQ(agr)],
		Memr[AGR_WR(agr)], Memr[AGR_OBS(agr)], Memr[AGR_WADO(agr)],
		naperts, IMT_NENTRIES(symbol))

	    # Compute the cumulative differences.
	    call ph_tcum (rap[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)],
		Memr[AGR_WADO(agr)], Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)],
		Memr[AGR_WCUM(agr)], naperts, Memd[AGR_PARAMS(agr)],
		IMT_RO(symbol), IMT_NXAIRMASS(symbol))

	    # Write the aperture photometry file.
	    call ph_wapcat (apcat, IMT_IMNAME(symbol), IMT_RO(symbol),
	        Memr[AGR_ADOPT(agr)], Memr[AGR_WADO(agr)], naperts, smallap,
		largeap)

	    if (logfd != NULL) {
		call ph_twrite (logfd, IMT_IMNAME(symbol), IMT_RO(symbol),
		    IMT_XAIRMASS(symbol), Memr[AGR_RBAR(agr)],
		    Memr[AGR_THEO(agr)], Memr[AGR_OBS(agr)],
		    Memr[AGR_WOBS(agr)], Memr[AGR_ADOPT(agr)],
		    Memr[AGR_WADO(agr)], rap[1,IMT_OFFSET(symbol)],
		    Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)],
		    Memr[AGR_WCUM(agr)], naperts)
		call fprintf (logfd, "\n")
		call ph_tmags (logfd, id[IMT_OFFSET(symbol)],
		    x[IMT_OFFSET(symbol)], y[IMT_OFFSET(symbol)],
		    nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)],
		    mag[1,IMT_OFFSET(symbol)], merr[1,IMT_OFFSET(symbol)],
		    Memr[AGR_ADOPT(agr)], Memr[AGR_W(agr)],
		    Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)],
		    Memr[AGR_WCUM(agr)], Memr[AGR_MAGS(agr)],
		    Memr[AGR_CMAGS(agr)], Memr[AGR_TMAGS(agr)],
		    Memr[AGR_WMAGS(agr)], naperts, IMT_NENTRIES(symbol))
		call ph_papcor (logfd, IMT_IMNAME(symbol), IMT_RO(symbol),
		    rap[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)],
		    Memr[AGR_WADO(agr)], naperts, smallap, largeap)
		call fprintf (logfd, "\n")
	    }

	    if (magfd != NULL) {
		call ph_wmags (magfd, symbol, x[IMT_OFFSET(symbol)],
		    y[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
		    rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
		    merr[1,IMT_OFFSET(symbol)], Memr[AGR_ADOPT(agr)],
		    Memr[AGR_WADO(agr)], Memr[AGR_W(agr)], Memr[AGR_CUM(agr)],
		    Memr[AGR_WCUM(agr)], Memr[AGR_MAGS(agr)],
		    Memr[AGR_CMAGS(agr)], Memr[AGR_WMAGS(agr)], naperts,
		    IMT_NENTRIES(symbol), largeap)
	    }

	    if (mgd != NULL) {
		call ph_gimfit (mgd, agr, IMT_IMNAME(symbol), IMT_RO(symbol),
	            IMT_XAIRMASS(symbol), nap[IMT_OFFSET(symbol)],
		    nap[IMT_OFFSET(symbol)], rap[1,IMT_OFFSET(symbol)],
		    mag[1,IMT_OFFSET(symbol)], naperts,
		    IMT_NENTRIES(symbol), YES)
	        if (! IS_INDEFR(IMT_RO(symbol))) {
		    call ph_gaimres (mgd, agr, IMT_IMNAME(symbol),
		        IMT_RO(symbol), IMT_XAIRMASS(symbol),
			nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
			rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
			naperts, IMT_NENTRIES(symbol), YES)
		    call ph_gbimres (mgd, agr, IMT_IMNAME(symbol),
		        IMT_RO(symbol), IMT_XAIRMASS(symbol),
			nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
			mag[1,IMT_OFFSET(symbol)], naperts,
			IMT_NENTRIES(symbol), YES)
		    call ph_gaximres (mgd, agr, IMT_IMNAME(symbol),
		        IMT_RO(symbol), IMT_XAIRMASS(symbol),
		        nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
		        x[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
			naperts, IMT_NENTRIES(symbol), YES)
		    call ph_gayimres (mgd, agr, IMT_IMNAME(symbol),
		        IMT_RO(symbol), IMT_XAIRMASS(symbol),
			nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
			y[IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
			naperts, IMT_NENTRIES(symbol), YES)
		    call ph_gacum (mgd, agr, IMT_IMNAME(symbol),
		        IMT_RO(symbol), IMT_XAIRMASS(symbol),
			nap[IMT_OFFSET(symbol)], nap[IMT_OFFSET(symbol)],
			rap[1,IMT_OFFSET(symbol)], mag[1,IMT_OFFSET(symbol)],
			naperts, IMT_NENTRIES(symbol), smallap, largeap, YES)
		}
	    }

	    if (! IS_INDEFR(IMT_RO(symbol))) {
	        call aaddr (Memr[AGR_AVE(agr)], Memr[AGR_TAVE(agr)],
	            Memr[AGR_TAVE(agr)], naperts)
	        call aaddr (Memr[AGR_RESID(agr)], Memr[AGR_TRESID(agr)],
	            Memr[AGR_TRESID(agr)], naperts)
	        call aaddr (Memr[AGR_RESSQ(agr)], Memr[AGR_TRESSQ(agr)],
	            Memr[AGR_TRESSQ(agr)], naperts)
	        call aaddr (Memr[AGR_WR(agr)], Memr[AGR_TWR(agr)],
	            Memr[AGR_TWR(agr)], naperts)
	    }
	}

	call sfree (sp)
end


# PH_A1EVAL -- Evaluate the fit for all the images.

procedure ph_a1eval (agr, image, r0, xairmass, nap, rap, mag, merr,
	naperts, npts)

pointer	agr			# pointer to the fitting structure
char	image[ARB]		# the image name
real	r0			# the seeing radius
real	xairmass		# the airmass value
int	nap[ARB]		# the number of apertures array 
real	rap[naperts,ARB]	# the list of aperture radii
real	mag[naperts,ARB]	# the magnitude difference array
real	merr[naperts,ARB]	# the magnitude error array
int	naperts			# the number of apertures
int	npts			# the number of points

begin
	# Initialize some vectors.
	call aclrr (Memr[AGR_AVE(agr)], naperts)
	call aclrr (Memr[AGR_RESID(agr)], naperts)
	call aclrr (Memr[AGR_RESSQ(agr)], naperts)
	call aclrr (Memr[AGR_WR(agr)], naperts)
	call aclrr (Memr[AGR_ADOPT(agr)], naperts)
	call aclrr (Memr[AGR_WOBS(agr)], naperts)

	# Compute the theoretical curve
	call ph_rinit (rap, Memr[AGR_RBAR(agr)], naperts)
	call ph_tinit (rap, Memr[AGR_THEO(agr)], naperts,
	    Memd[AGR_PARAMS(agr)], r0, xairmass)

	# Accumulate the difference data.
	call ph_taccum (nap, mag, merr, Memr[AGR_THEO(agr)], Memr[AGR_W(agr)],
	    Memr[AGR_ADOPT(agr)], Memr[AGR_WOBS(agr)], Memr[AGR_AVE(agr)],
	    Memr[AGR_RESID(agr)], Memr[AGR_RESSQ(agr)], Memr[AGR_WR(agr)],
	    Memr[AGR_OBS(agr)], Memr[AGR_WADO(agr)], naperts, npts)

	# Compute the cumulative differences.
	call ph_tcum (rap, Memr[AGR_ADOPT(agr)], Memr[AGR_WADO(agr)],
	    Memr[AGR_CUM(agr)], Memr[AGR_TCUM(agr)], Memr[AGR_WCUM(agr)],
	    naperts, Memd[AGR_PARAMS(agr)], r0, xairmass)
end


# PH_TFSHOW -- Compute and print the summary statistics.

procedure ph_tfshow (fd, agr, naperts)

int	fd			# the output file descriptor
pointer	agr			# the pointer to the fit structure
int	naperts			# the number of apertures

begin
	call ph_rwrite (fd, Memr[AGR_TAVE(agr)], Memr[AGR_TWR(agr)],
	    Memr[AGR_TRESID(agr)], Memr[AGR_TRESSQ(agr)], naperts)
end


# PH_AMEMFREE -- Free the memory used for doing the fit

procedure ph_amemfree (agr)

pointer	agr			# pointer to the fitting structure

begin
	call mfree (AGR_DM(agr), TY_DOUBLE)
	call mfree (AGR_DDDR(agr), TY_DOUBLE)
	call mfree (AGR_T(agr), TY_DOUBLE)
	call mfree (AGR_U(agr), TY_DOUBLE)
	call mfree (AGR_V(agr), TY_DOUBLE)
	call mfree (AGR_POLD(agr), TY_DOUBLE)
	call mfree (AGR_DP(agr), TY_DOUBLE)
	call mfree (AGR_PARAMS(agr), TY_DOUBLE)
	call mfree (AGR_PERRORS(agr), TY_DOUBLE)
	call mfree (AGR_PCLAMPS(agr), TY_DOUBLE)

	call mfree (AGR_RBAR(agr), TY_REAL)
	call mfree (AGR_W(agr), TY_REAL)
	call mfree (AGR_THEO(agr), TY_REAL)
	call mfree (AGR_WR(agr), TY_REAL)
	call mfree (AGR_ADOPT(agr), TY_REAL)
	call mfree (AGR_WOBS(agr), TY_REAL)
	call mfree (AGR_OBS(agr), TY_REAL)
	call mfree (AGR_WADO(agr), TY_REAL)

	call mfree (AGR_CUM(agr), TY_REAL)
	call mfree (AGR_TCUM(agr), TY_REAL)
	call mfree (AGR_WCUM(agr), TY_REAL)
	call mfree (AGR_MAGS(agr), TY_REAL)
	call mfree (AGR_CMAGS(agr), TY_REAL)
	call mfree (AGR_TMAGS(agr), TY_REAL)
	call mfree (AGR_WMAGS(agr), TY_REAL)

	call mfree (AGR_AVE(agr), TY_REAL)
	call mfree (AGR_RESID(agr), TY_REAL)
	call mfree (AGR_RESSQ(agr), TY_REAL)
	call mfree (AGR_TAVE(agr), TY_REAL)
	call mfree (AGR_TRESID(agr), TY_REAL)
	call mfree (AGR_TRESSQ(agr), TY_REAL)
	call mfree (AGR_TWR(agr), TY_REAL)

	call mfree (agr, TY_STRUCT)
end


# PH_AGETP -- Fetch the initial values of the parameters for the cog model.

procedure ph_agetp (params)

double	params[ARB]

real	clgetr()

begin
	params[1] = clgetr ("swings")
	params[2] = clgetr ("pwings")
	params[3] = clgetr ("pgauss")
	params[4] = clgetr ("rgescale")
	params[5] = clgetr ("xwings")
end


# PH_APINIT -- Set the initial values for the curve of growth model parameters.

procedure ph_apinit (iparams, oparams, perrors, pclamps) 

double	iparams[ARB]		# the input parameters array
double	oparams[ARB]		# the output parameters array
double	perrors[ARB]		# the parameter errors array
double	pclamps[ARB]		# array of parameter clamps

begin
	call amovd (iparams, oparams, MAX_MTERMS)
	call aclrd (perrors, MAX_MTERMS)
	call amovkd (0.2d0, pclamps, MAX_MTERMS)
	pclamps[4] = 0.5d0
end


# PH_ANTERMS -- Compute the number of terms to use in the fit.

procedure ph_anterms (niter, lterms, gain, nterms)

int	niter			# the current iteration
int	lterms			# the maximum number of terms to be fit
real	gain			# the current value of the gain
int	nterms			# the current number of terms to be fit

int	n

begin
	n = 1
	if (nterms >= 1 && gain <= 0.04) {
	    n = 2
	    if (nterms >= 2 && gain <= 0.016) {
	        n = 3
		if (nterms >= 3 && gain <= 0.006) {
	            n = 4
		    if (nterms >= 4 && gain <= 0.0025)
	    	        n = 5
		}
	    }
	}

	nterms = min (n, lterms)
end


# PH_ADDDR -- Compute the derivative wrt ro vector.

procedure ph_adddr (rap, dddr, naperts, params, r0, xairmass)

real	rap[ARB]		# array of aperture radii
double  dddr[ARB]		# the output derivatives array
int	naperts			# the number of apertures
double  params[ARB]		# the input parameter array
real	r0			# the current value of r0
real	xairmass		# the current value of the airmass

int	i
double	ph_dmag()

begin
	do i = 2, naperts {
	    dddr[i] = 500.0d0 * (ph_dmag (rap[i-1], rap[i], xairmass,
	        r0 - 0.001, params) - ph_dmag (rap[i-1], rap[i], xairmass,
		r0 + 0.001, params))
	}
end


# PH_ADDDRDM -- Compute the derivative wrt ro vector.

procedure ph_adddrdm (rap, dm, dddr, naperts, params, r0, xairmass, niter)

real	rap[ARB]		# array of aperture radii
double	dm[ARB]			# the output model differences array
double  dddr[ARB]		# the output derivatives array
int	naperts			# the number of apertures
double  params[ARB]		# the input parameter array
real	r0			# the current value of r0
real	xairmass		# the current value of the airmass
int	niter			# the current iteratiom

int	i
double	ph_dmag()

begin
	do i = 2, naperts {
	    dm[i] = ph_dmag (rap[i-1], rap[i], xairmass, r0, params)
	    #if (niter == 1)
	        dddr[i] = 500.0d0 * (ph_dmag (rap[i-1], rap[i], xairmass,
	            r0 - 0.001, params) - ph_dmag (rap[i-1], rap[i], xairmass,
		    r0 + 0.001, params))
	}
end


# PH_AWSUM -- Accumulate the weighted sums necessary to fit r0.

procedure ph_awsum (nap, mag, merr, dm, w, dddr, naperts, npts, ipow,
	sumr, sumy)

int	nap[ARB]			# array of aperture numbers
real	mag[naperts,ARB]		# array of magnitude difference
real	merr[naperts,ARB]		# array of magnitude errors
double	dm[ARB]				# array of model differences
real	w[ARB]				# array of weights
double	dddr[ARB]			# array of model derivatives
int	naperts				# number of apertures
int	npts				# number of data points
int	ipow				# weighting power factor
double	sumr, sumy			# the output sums

int	i, j
real	diff, wt

begin
	do i = 1, npts {
	    w[1] = 1.0d0
	    if (nap[i] <= 1)
		next
	    do j = 2, nap[i] {
		diff = mag[j,i] - dm[j]
		w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) **
		    ipow)
		w[j] = min (w[j], w[j-1])
		wt = w[j] / merr[j,i] ** 2
		sumr = sumr + wt * diff * dddr[j]
		sumy = sumy + wt * dddr[j] ** 2
	    }
	}
end


# PH_ADP -- Compute the parameter derivative vector with the new
# value of r0.

procedure ph_adp (rap, dm, t, naperts, lterms, nterms, params, r0, xairmass)

real	rap[ARB]		# array of aperture radii
double	dm[ARB]			# the new model differences
double	t[lterms,ARB]		# the parameter derivatives matrix
int	naperts			# the number of apertures
int	lterms			# the size of the derivatives matrix
int	nterms			# the number of terms to be fit
double	params[ARB]		# the current model parameter values
real	r0			# the current r0 value
real	xairmass		# the current airmass

int	i, j
double	ph_dmag()

begin
	do j = 2, naperts {
	    dm[j] = ph_dmag (rap[j-1], rap[j], xairmass, r0, params)
	    do i = 1, nterms {
		params[i] = params[i] - 0.001
		t[i,j] = ph_dmag (rap[j-1], rap[j], xairmass, r0, params)
		params[i] = params[i] + 0.002
		t[i,j] = 500.0d0 * (t[i,j] - ph_dmag (rap[j-1], rap[j],
		    xairmass, r0, params))
		params[i] = params[i] - 0.001
	    }
	}
end


# PH_AACCUM -- Accumulate the matrix and vector required to solve for the
# parameter increments

procedure ph_aaccum (nap, mag, merr, dm, w, naperts, npts, t, u, v, lterms,
	nterms, ipow, sumd, sumw, sumn)

int	nap[ARB]			# array of aperture numbers
real	mag[naperts,ARB]		# array of magnitude difference
real	merr[naperts,ARB]		# array of magnitude errors
double	dm[ARB]				# array of model differences
real	w[ARB]				# array of weights
int	naperts				# number of apertures
int	npts				# number of data points
double	t[lterms,ARB]			# the array of parameter derivatives
double  u[lterms,ARB]			# the matrix to be accumulated
double	v[ARB]				# the vector to be accumulated
int	lterms				# the maximum number of terms to fit
int	nterms				# the current number of terms to fit
int	ipow				# power factor for the weights
double	sumd				# sum of the differences
double	sumw				# sum of the scatter
double	sumn				# sum of the weights

int	i, j, l, m
real	diff, wt

begin
	do i = 1, npts {
	    w[1] = 1.0
	    if (nap[i] <= 1)
		next
	    do j = 2, nap[i] {
		diff = mag[j,i] - dm[j]
		w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** ipow)
		w[j] = min (w[j], w[j-1])
		wt = w[j] / merr[j,i] ** 2
		sumd = sumd + wt * diff ** 2
		sumw = sumw + wt
		sumn = sumn + w[j]
		do l = 1, nterms {
		    v[l] = v[l] + wt * diff * t[l,j]
		    do m = 1, nterms
			u[l,m] = u[l,m] + wt * t[l,j] * t[m,j]
		}
	    }
	}
end


# PH_APSOLVE -- Solve for new values of the parameters.

procedure ph_apsolve (pold, params, dp, pclamps, nterms, ipow)

double	pold[ARB]		# the previous parameter imcrement values
double	params[ARB]		# the current values of the parameters
double	dp[ARB]			# the parameter increment values
double	pclamps[ARB]		# the parameter clamp values
int	nterms			# the number of terms that were fit
int	ipow			# the current weighting factor

int	i

begin
	do i = 1, nterms {
	    if (dp[i] * pold[i] < 0.0d0)
		pclamps[i] = 0.5 * pclamps[i]
	    pold[i] = dp[i]
	}

	params[1] = params[1] - dp[1] / ( 1.0d0 + abs (dp[1] / (pclamps[1] *
	    (params[1] - 1.0d0))))
	if (nterms >= 2) {
	    params[2] = params[2] - dp[2] / (1.0d0 + abs (dp[2] / (pclamps[2] *
		params[2] * (1.0d0 - params[2]))))
	    if (nterms >= 3) {
		params[3] = params[3] - dp[3] / (1.0d0 + abs (dp[3] /
		    (pclamps[3] * params[3] * (1.0d0 - params[3]))))
		if (ipow == 1) {
		    ipow = 2
		    pclamps[1] = 0.1d0
		    pclamps[2] = 0.1d0
		}
		if (nterms >= 4) {
		    params[4] = params[4] - dp[4] / (1.0d0 + abs (dp[4] /
			(pclamps[4] * params[4])))
		    if (nterms >= 5)
			params[5] = params[5] - dp[5] / (1.0d0 + abs (dp[5] /
			    (pclamps[5] * params[2])))
		}
	    }
	}
end

# PH_PWRITE -- Write out the theoetical model parameters.

procedure ph_pwrite (fd, params, perrors, max_nterms, nterms)

int	fd			# the output file descriptor
double	params[ARB]		# the current model parameter values
double	perrors[ARB]		# the current parameter errors
int	max_nterms		# the number of terms
int	nterms			# the number of terms to fit

int	i

begin
	# Print out the best fit parameters.
	call fprintf (fd,
	"\nThe computed cog model parameters and their errors for nparams %d\n")
	    call pargi (nterms)
	do i = 1, max_nterms {
	    call fprintf (fd, "\t%10s: %15.7g +/- %15.7g\n")
		switch (i) {
		case 1:
		    call pargstr ("swings")
		case 2:
		    call pargstr ("pwings")
		case 3:
		    call pargstr ("pgauss")
		case 4:
		    call pargstr ("rgescale")
		case 5:
		    call pargstr ("xwings")
		}
		call pargd (params[i])
		call pargd (perrors[i])
	}
	call fprintf (fd, "\n")
end

# PH_RINIT -- Initialize the rbar vector.

procedure ph_rinit (rap, rbar, naperts)

real	rap[ARB]		# the array of aperture radii
real	rbar[ARB]		# the mean radius estimates
int	naperts			# the number of aperture radii

int	j

begin
	do j = 2, naperts 
	    rbar[j] = 0.5 * (rap[j-1] + rap[j])
end

# PH_TINIT -- Initialize the rbar vector and compute the theoretical estimates.

procedure ph_tinit (rap, theo, naperts, params, r0, airmass)

real	rap[ARB]		# the array of aperture radii
real	theo[ARB]		# the current model estimates
int	naperts			# the number of aperture radii
double	params[ARB]		# the current parameter estimates
real	r0			# the seeing radius
real	airmass			# the airmass value

int	j
double	ph_dmag()

begin
	if (IS_INDEFR(r0))
	    call amovkr (INDEFR, theo[2], naperts - 1)
	else {
	    do j = 2, naperts 
	        theo[j] = ph_dmag (rap[j-1], rap[j], airmass, r0, params)
	}
end


# PH_TACCUM -- Accumulate the data.

procedure ph_taccum (nap, mag, merr, theo, w, adopt, wobs, ave, resid, ressq,
	wr, obs, wado, naperts, npts)

int	nap[ARB]			# array of aperture numbers
real	mag[naperts,ARB]		# the array of magnitude differences
real	merr[naperts,ARB]		# the array of magnitude errors
real	theo[ARB]			# the current model values
real	w[ARB]				# the working weight array
real	adopt[ARB]			# the adopted weight differences
real	wobs[ARB]			# the sum of the weights
real	ave[ARB]			# the average model values
real	resid[ARB]			# the residuals
real	ressq[ARB]			# the residuals squared
real	wr[ARB]				# the sum of the weights
real	obs[ARB]			# the observations
real	wado[ARB]			# the error in the observations
int	naperts				# the number of apertures
int	npts				# the number of npts

int	i, j
real	diff, wt, scale

begin
	do i = 1, npts {
	    w[1] = 1.0
	    if (nap[i] <= 1)
		next
	    do j = 2, nap[i] {
		if (IS_INDEFR(theo[j])) {
		    wt = 1.0 / merr[j,i] ** 2
		    ave[1] = INDEFR
		    ave[j] = INDEFR
		    resid[1] = INDEFR
		    resid[j] = INDEFR
		    ressq[1] = INDEFR
		    ressq[j] = INDEFR
		    wr[1] = INDEFR
		    wr[j] = INDEFR
		} else {
		    diff = mag[j,i] - theo[j]
		    w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])))
		    w[j] = min (w[j], w[j-1])
		    wt = w[j] / merr[j,i] ** 2
		    ave[1] = ave[1] + wt * theo[j]
		    ave[j] = ave[j] + wt * theo[j]
		    resid[1] = resid[1] + wt * diff
		    resid[j] = resid[j] + wt * diff
		    ressq[1] = ressq[1] + wt * diff ** 2
		    ressq[j] = ressq[j] + wt * diff ** 2
		    wr[1] = wr[1] + wt
		    wr[j] = wr[j] + wt
		}
		adopt[j] = adopt[j] + wt * mag[j,i]
		wobs[j] = wobs[j] + wt
	    }
	}

	do j = 2, naperts {
	    if (wobs[j] <= 0.0)
		next
	    obs[j] = adopt[j] / wobs[j]
	    adopt[j] = 0.0
	    wobs[j] = 0.0
	}

	do i = 1, npts {
	    w[1] = 1.0
	    if (nap[i] <= 1)
		next
	    do j = 2, nap[i] {
		diff = mag[j,i] - obs[j]
		w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** 2)
		w[j] = min (w[j], w[j-1])
		wt = w[j] / merr[j,i] ** 2
		adopt[j] = adopt[j] + wt * mag[j,i]
		wobs[j] = wobs[j] + wt
	    }
	}

	do j = 2, naperts {
	    if (wobs[j] <= 0.0)
		next
	    obs[j] = adopt[j] / wobs[j]
	    adopt[j] = 0.0
	    wobs[j] = 0.0
	}

	do i = 1, npts {
	    w[1] = 1.0
	    if (nap[i] <= 1)
		next
	    do j = 2, nap[i] {
		diff = mag[j,i] - obs[j]
		w[j] = 1.0 / (1.0 + abs (diff / (2.0 * merr[j,i])) ** 3)
		w[j] = min (w[j], w[j-1])
		wt = w[j] / merr[j,i] ** 2
		adopt[j] = adopt[j] + wt * mag[j,i]
		wobs[j] = wobs[j] + wt
	    }
	}

	scale = 0.1
	do j = 2, naperts {
	    if (wobs[j] > 0.0) {
		obs[j] = adopt[j] / wobs[j]
		if (IS_INDEFR(theo[j]))
		    wt = 0.0
		else {
		    wt = 1.0 / (scale * theo[j]) ** 2
		    adopt[j] = adopt[j] + wt * theo[j]
		}
		wt = 1.0 / (wobs[j] + wt)
		adopt[j] = wt * adopt[j]
		#wado[j] = sqrt (wt)
		wado[j] = sqrt (wt + (scale * adopt[j]) ** 2)
		wobs[j] = sqrt (1.0 / wobs[j])
	    } else {
		adopt[j] = theo[j]
		if (IS_INDEFR(theo[j])) {
		    wado[j] = INDEFR
		    obs[j] = INDEFR
		    wobs[j] = INDEFR
		} else
		    wado[j] = 2.0 * scale * abs (theo[j])
	    }
	}
end


# PH_TCUM -- Compute the cumulative differences.

procedure ph_tcum (rap, adopt, wado, cum, tcum, wcum, naperts, params,
	r0, airmass)

real	rap[ARB]		# the list of aperture radii
real	adopt[ARB]		# the adopted differences
real	wado[ARB]		# the errors in the adopted differences
real	cum[ARB]		# the accumulated differences
real	tcum[ARB]		# the total accumulated differences
real	wcum[ARB]		# the accumulated difference errors
int	naperts			# the number of aperture radii
double	params[ARB]		# the current parameter values
real	r0			# the seeing radius
real	airmass			# the airmass

int	i
double	ph_dmag()

begin
	if (IS_INDEFR(r0)) {
	    call amovkr (INDEFR, cum, naperts + 1)
	    call amovkr (INDEFR, tcum, naperts + 1)
	    call amovkr (INDEFR, wcum, naperts + 1)
	} else {
	    cum[naperts+1] = 0.0
	    tcum[naperts+1] = ph_dmag (rap[naperts], 2.0 * rap[naperts],
	        airmass, r0, params)
	    wcum[naperts+1] = 0.0
	    do i = naperts, 2, -1 {
	        cum[i] = adopt[i] + cum[i+1]
	        tcum[i] = adopt[i] + tcum[i+1]
	        wcum[i] = wado[i] ** 2 + wcum[i+1]
	    }
	}
end


# PH_WAPCAT -- Write the image entry to the aperture correction catalog.

procedure ph_wapcat (apcat, image, r0, adopt, wado, naperts, smallap, largeap)

int	apcat			# the aperture correction catalog descriptor
char	image[ARB]		# the image name
real	r0			# the seeing radius
real	adopt[ARB]		# the adopted difference values
real	wado[ARB]		# the adopted difference errors
int	naperts			# the number of apertures
int	smallap			# the small aperture number
int	largeap			# the large aperture number

int	i
real	cum, wcum

begin
	if (IS_INDEFR(r0)) {
	    cum = INDEFR
	    wcum = INDEFR
	} else {
	    cum = 0.0
	    wcum = 0.0
	    do i = largeap, smallap + 1, -1 {
	        cum = cum + adopt[i]
	        wcum = wcum + wado[i] ** 2
	    }
	    wcum = sqrt (wcum)
	}
	call fprintf (apcat, "%s %g %g\n")
	    call pargstr (image)
	    call pargr (cum)
	    call pargr (wcum)
end


# PH_PAPCOR -- Print the aperture correction on the standard output.

procedure ph_papcor (fd, image, r0, rap, adopt, wado, naperts, smallap, largeap)

int	fd			# the output file descriptor
char	image[ARB]		# the image name
real	r0			# the seeing radius
real	rap[ARB]		# the aperture radii
real	adopt[ARB]		# the adopted difference values
real	wado[ARB]		# the adopted difference errors
int	naperts			# the number of apertures
int	smallap			# the small aperture number
int	largeap			# the large aperture number

int	i
real	cum, wcum

begin
	if (IS_INDEFR(r0)) {
	    cum = INDEFR
	    wcum = INDEFR
	}  else {
	    cum = 0.0
	    wcum = 0.0
	    do i = largeap, smallap + 1, -1 {
	        cum = cum + adopt[i]
	        wcum = wcum + wado[i] ** 2
	    }
	    wcum = sqrt (wcum)
	}

	call fprintf (fd,
	    "Image: %s rin=%.2f rout=%.2f apercor=%.3f +/- %.4f\n")
	    call pargstr (image)
	    call pargr (rap[smallap])
	    call pargr (rap[largeap])
	    call pargr (cum)
	    call pargr (wcum)
end


define	NPERLINE	7

# PH_TWRITE -- Write the results to the output file.

procedure ph_twrite (fd, image, r0, xairmass, rbar, theo, obs, wobs, adopt,
	wado, rap, cum, tcum, wcum, naperts)

int	fd			# the file descriptor
char	image			# the iamge name
real	r0			# the seeing disk
real	xairmass		# the assumed airmass
real	rbar[ARB]		# the list of mean aperture radii
real	theo[ARB]		# the theoretical model differences
real	obs[ARB]		# the observed differences
real	wobs[ARB]		# the observed difference errors
real	adopt[ARB]		# the adopted differences
real	wado[ARB]		# the adopted difference errors
real	rap[ARB]		# the list of aperture radii
real	cum[ARB]		# the cumulative differences
real	tcum[ARB]		# the total cumulative differences
real	wcum[ARB]		# the errors in the cumulative differences
int	naperts			# the number of aperture

int	j
real	sqwcum

begin
	# Print the title.
	call fprintf (fd, "\nThe cog for image %s from radius %.2f to %.2f\n")
	    call pargstr (image)
	    call pargr (rbar[2])
	    call pargr (rbar[naperts])

	# Print the mean aperture radius.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, "  radius  %9.4f")
		    call pargr (rbar[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (rbar[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (rbar[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (rbar[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the theoretical model.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, "   model  %9.4f")
		    call pargr (theo[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (theo[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (theo[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (theo[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the observed cog.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, "observed  %9.4f")
		    call pargr (obs[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (obs[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (obs[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (obs[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the errors in the observed cog.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, "   sigma  %9.4f")
		    call pargr (wobs[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (wobs[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (wobs[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (wobs[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the adopted cog.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, " adopted  %9.4f")
		    call pargr (adopt[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (adopt[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (adopt[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (adopt[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the errors in the adopted cog.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, "   sigma  %9.4f")
		    call pargr (wado[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (wado[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (wado[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (wado[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the title.
	call fprintf (fd,
	    "\nThe aperture correction for image %s from radius %.2f to %.2f\n")
	    call pargstr (image)
	    call pargr (rap[1])
	    call pargr (rap[naperts])

	# Print the aperture radii.
	do j = 1, naperts {
	    if (j == 1) {
		call fprintf (fd, "  radius  %9.4f")
		    call pargr (rap[j])
	    } else if (mod (j, NPERLINE) == 0) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (rap[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "          %9.4f")
		    call pargr (rap[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (rap[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 0)
	    call fprintf (fd, "\n")

	# Print the aperture correction.
	do j = 2, naperts {
	    if (j == 2) {
		call fprintf (fd, " apercor  %9.4f")
		    call pargr (cum[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (cum[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (cum[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (cum[j])
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")


	# Print the error in the aperture correction.
	do j = 2, naperts {
	    sqwcum = wcum[j]
	    if (j == 2) {
		call fprintf (fd, "   sigma  %9.4f")
		    call pargr (sqwcum)
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (sqwcum)
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (sqwcum)
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (sqwcum)
	    }
	}
	if (mod (naperts, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	# Print the title.
	call fprintf (fd,
	    "\nThe aperture correction for image %s from radius %.2f to %.2f\n")
	    call pargstr (image)
	    call pargr (rap[1])
	    call pargr (2.0 * rap[naperts])

	# Print the total aperture correction.
	do j = 2, naperts + 1 {
	    if (j == 2) {
		call fprintf (fd, "tapercor  %9.4f")
		    call pargr (tcum[j])
	    } else if (mod (j, NPERLINE) == 1) {
		call fprintf (fd, "%9.4f\n")
		    call pargr (tcum[j])
	    } else if (mod (j, NPERLINE) == 2) {
		call fprintf (fd, "          %9.4f")
		    call pargr (tcum[j])
	    } else {
		call fprintf (fd, "%9.4f")
		    call pargr (tcum[j])
	    }
	}
	if (mod (naperts + 1, NPERLINE) != 1)
	    call fprintf (fd, "\n")

	call fprintf (fd, "\n")
end


# PH_TMAGS -- Compute the correction to the last computed magnitude for
# each star and the total magnitude as a function of aperture.

procedure ph_tmags (logfd, id, x, y, nap, rap, mag, merr, adopt, w, cum,
	tcum, wcum, tmpmags, obs, tobs, wobs, naperts, npts)

int	logfd			# the output file descriptor
int	id[ARB]			# stellar id numbers
real	x[ARB]			# stellar x coordinates
real	y[ARB]			# stellar y coordinates
int	nap[ARB]		# array of aperture numbers
real	rap[naperts,ARB]	# the list of aperture radii
real	mag[naperts,ARB]	# the input magnitude difference array
real	merr[naperts,ARB]	# the input magnitude error array
real	adopt[ARB]		# the adopted difference values
real	w[ARB]			# the working weight array
real	cum[ARB]		# the accumulated difference array
real	tcum[ARB]		# the total accumulated difference array
real	wcum[ARB]		# the errors in the accumulated diff array
real	tmpmags[ARB]		# temporary magnitude array
real	obs[ARB]		# the accumulated magnitude array
real	tobs[ARB]		# the total accumulated magnitude array
real	wobs[ARB]		# the observations array
int	naperts			# number of apertures
int	npts			# number of points

int	i, j
real	diff

begin
	# Print the logfile banner.
	if (logfd != NULL) {
	    call fprintf (logfd, "\nThe observed, and corrected to radii %.2f ")
		call pargr (rap[naperts,1])
	    call fprintf (logfd, "and %.2f, magnitudes\n")
		call pargr (2.0 * rap[naperts,1])
	}

	do i = 1, npts {

	    # Compute the observed, correctged and estimated total 
	    # magnitudes.
	    tmpmags[1] = mag[1,i]
	    do j = 1, nap[i] {
		if (j == 1)
		    w[j] = 1.0
		else {
		    diff = mag[j,i] - adopt[j]
		    tmpmags[j] = tmpmags[j-1] + mag[j,i]
		    w[j] = 1.0 / (1.0 + (diff / (2.0 * merr[j,i])) ** 2)
		    w[j] = min (w[j], w[j-1])
		}
		if (IS_INDEFR(cum[j+1])) {
		    obs[j]  = INDEFR
		    tobs[j]  = INDEFR
		    wobs[j] = INDEFR
		} else {
		    obs[j]  = tmpmags[j] + cum[j+1]
		    tobs[j]  = tmpmags[j] + tcum[j+1]
		    wobs[j] = sqrt (wcum[j+1] + merr[j,i] ** 2 / w[j])
		}
	    }
	    do j = nap[i] + 1, naperts {
		obs[j]  = INDEFR
		tobs[j]  = INDEFR
		wobs[j] = INDEFR
	    }

	    # Write out the results for the star to the log file.

	    # Print the banner of the star.
	    call fprintf (logfd, "Star: %d  x: %.3f y: %.3f\n")
	        call pargi (id[i])
	        call pargr (x[i])
	        call pargr (y[i])

	    # Print the aperture radii.
	    do j = 1, naperts {
	        if (j == 1) {
		    call fprintf (logfd, "  radius  %9.4f")
		        call pargr (rap[j,i])
	        } else if (mod (j, NPERLINE) == 0) {
		    call fprintf (logfd, "%9.4f\n")
		        call pargr (rap[j,i])
	        } else if (mod (j, NPERLINE) == 1) {
		    call fprintf (logfd, "          %9.4f")
		        call pargr (rap[j,i])
	        } else {
		    call fprintf (logfd, "%9.4f")
		        call pargr (rap[j,i])
	        }
	    }
	    if (mod (naperts, NPERLINE) != 0)
	        call fprintf (logfd, "\n")

	    # Print the observed magnitudes.
	    do j = 1, naperts {
	        if (j == 1) {
		    call fprintf (logfd, "    mags  %9.4f")
	                call pargr (tmpmags[j])
	        } else if (mod (j, NPERLINE) == 0) {
		    call fprintf (logfd, "%9.4f\n")
		        call pargr (tmpmags[j])
	        } else if (mod (j, NPERLINE) == 1) {
		    call fprintf (logfd, "          %9.4f")
		        call pargr (tmpmags[j])
	        } else {
		    call fprintf (logfd, "%9.4f")
		        call pargr (tmpmags[j])
	        }
	    }
	    if (mod (naperts, NPERLINE) != 0)
	        call fprintf (logfd, "\n")

	    # Print the corrected magnitudes.
	    do j = 1, naperts {
	        if (j == 1) {
		    call fprintf (logfd, "   cmags  %9.4f")
		        call pargr (obs[j])
	        } else if (mod (j, NPERLINE) == 0) {
		    call fprintf (logfd, "%9.4f\n")
		        call pargr (obs[j])
	        } else if (mod (j, NPERLINE) == 1) {
		    call fprintf (logfd, "          %9.4f")
		        call pargr (obs[j])
	        } else {
		    call fprintf (logfd, "%9.4f")
		        call pargr (obs[j])
	        }
	    }
	    if (mod (naperts, NPERLINE) != 0)
	        call fprintf (logfd, "\n")

	    # Print the estimated total corrected magnitudes.
	    do j = 1, naperts {
	        if (j == 1) {
		    call fprintf (logfd, "  tcmags  %9.4f")
		        call pargr (tobs[j])
	        } else if (mod (j, NPERLINE) == 0) {
		    call fprintf (logfd, "%9.4f\n")
		        call pargr (tobs[j])
	        } else if (mod (j, NPERLINE) == 1) {
		    call fprintf (logfd, "          %9.4f")
		        call pargr (tobs[j])
	        } else {
		    call fprintf (logfd, "%9.4f")
		        call pargr (tobs[j])
	        }
	    }
	    if (mod (naperts, NPERLINE) != 0)
	        call fprintf (logfd, "\n")

	    # Print the errors in the total magnitudes
	    do j = 1, naperts {
	        if (j == 1) {
		    call fprintf (logfd, "   sigma  %9.4f")
		        call pargr (wobs[j])
	        } else if (mod (j, NPERLINE) == 0) {
		    call fprintf (logfd, "%9.4f\n")
		        call pargr (wobs[j])
	        } else if (mod (j, NPERLINE) == 1) {
		    call fprintf (logfd, "          %9.4f")
		        call pargr (wobs[j])
	        } else {
		    call fprintf (logfd, "%9.4f")
		        call pargr (wobs[j])
	        }
	    }
	    if (mod (naperts, NPERLINE) != 0)
	        call fprintf (logfd, "\n")
	}

	call fprintf (logfd, "\n")
end


# PH_WMAGS -- Compute the correction to the last computed magnitude for
# each star and the total magnitude as a function of aperture.

procedure ph_wmags (magfd, imsymbol, x, y, nap, rap, mag, merr, adopt, wado,
    w, cum, wcum, tmpmags, obs, wobs, naperts, npts, largeap)

int	magfd			# the best magnitudes file descriptor
pointer	imsymbol		# the image symbol
real	x[ARB]			# stellar x coordinates
real	y[ARB]			# stellar y coordinates
int	nap[ARB]		# array of aperture numbers
real	rap[naperts,ARB]	# the input aperture radii
real	mag[naperts,ARB]	# the input magnitude difference array
real	merr[naperts,ARB]	# the input magnitude error array
real	adopt[ARB]		# the adopted difference values
real	wado[ARB]		# the adopted difference errors
real	w[ARB]			# the working weight array
real	cum[ARB]		# the accumulated difference array
real	wcum[ARB]		# the errors in the accumulated diff array
real	tmpmags[ARB]		# temporary magnitude array
real	obs[ARB]		# the accumulated magnitude array
real	wobs[ARB]		# the observations array
int	naperts			# number of apertures
int	npts			# number of points
int	largeap			# the largest aperture

int	i, j, jfinal
real	diff, sigmin, sigcor
real	assqr()

begin
	do i = 1, npts {

	    # Compute the observed, correctged and estimated total magnitudes.
	    tmpmags[1] = mag[1,i]
	    sigmin = MAX_REAL
	    #jfinal = min (nap[i], largeap)
	    jfinal = largeap
	    if (largeap < nap[i])
		sigcor = assqr (wado[largeap+1], nap[i] - largeap)
	    else
		sigcor = 0.0
	    do j = 1, min (nap[i], largeap) {
		if (j == 1)
		    w[j] = 1.0
		else {
		    diff = mag[j,i] - adopt[j]
		    tmpmags[j] = tmpmags[j-1] + mag[j,i]
		    w[j] = 1.0 / (1.0 + (diff / (2.0 * merr[j,i])) ** 2)
		    w[j] = min (w[j], w[j-1])
		}
		if (IS_INDEFR(cum[j+1])) {
		    obs[j] = INDEFR
		    wobs[j] = INDEFR
		} else {
		    obs[j]  = tmpmags[j] + cum[j+1] - cum[largeap+1]
		    wobs[j] = sqrt (wcum[j+1] - sigcor + merr[j,i] ** 2 / w[j])
		    if (wobs[j] < sigmin) {
		        jfinal = j
		        sigmin = wobs[j]
		    }
		}
	    }

	    # Write out the best magnitude.
	    if (magfd != NULL) {
		call fprintf (magfd,
"%23.23s  %10.10s  %8.2f  %6.3f  %11.1h  %8.2f  %8.2f  %7.3f  %7.3f  %6.3f\n")
		    call pargstr (IMT_IMNAME(imsymbol))
		    call pargstr (IMT_IFILTER(imsymbol))
		    call pargr (IMT_ITIME(imsymbol))
		    call pargr (IMT_XAIRMASS(imsymbol))
		    call pargr (IMT_OTIME(imsymbol))
		    call pargr (x[i])
		    call pargr (y[i])
		    call pargr (obs[jfinal])
		    call pargr (wobs[jfinal])
		    call pargr (rap[jfinal,i])
	    }
	}
end


# PH_ACHI -- Compute the chi statistic.

real procedure ph_achi (wr, resid, ressq, naperts)

real	wr[ARB]			# the weights
real	resid[ARB]		# the residuals
real	ressq[ARB]		# the residuals
int	naperts			# the number of apertures

double	sumwr, sumres, sumressq
real	chi
real	asumr()

begin
	sumwr = asumr (wr, naperts)
	sumres = asumr (resid, naperts)
	sumressq = asumr (ressq, naperts)
	if (sumwr <= 0.0d0)
	    chi = 0.0
	else
	    chi = sumressq / sumwr - (sumres / sumwr) ** 2
	if (chi <= 0.0 || chi > MAX_REAL)
	    return (INDEFR)
	else
	    return (chi)
end


# PH_RWRITE -- Write out the fit statistics for each aperture.

procedure ph_rwrite (fd, ave, wr, resid, ressq, naperts)

int	fd			# the output file descriptor
real	ave[ARB]		# the average values
real	wr[ARB]			# the weights
real	resid[ARB]		# the residuals
real	ressq[ARB]		# the residuals
int	naperts			# the number of apertures

int	j
real	rtmp

begin
	call fprintf (fd, "\nAverage model cog, residual, and rms residual ")
	call fprintf (fd, "for each aperture all images\n")
	do j = 1, naperts {
	    if (wr[j] <= 0.0) {
		ave[j] = INDEFR
		resid[j] = INDEFR
		ressq[j] = INDEFR
	    } else {
	        ave[j] = ave[j] / wr[j]
	        resid[j] = resid[j] / wr[j]
	        rtmp = ressq[j] / wr[j] - resid[j] ** 2
		if (rtmp <= 0.0)
		    ressq[j] = 0.0
		else
	            ressq[j] = sqrt (rtmp)
	    }
	    call fprintf (fd, "\t%3d  %9.5f  %9.5f  %9.5f\n")
		call pargi (j)
		call pargr (ave[j])
		call pargr (resid[j])
		call pargr (ressq[j])
	}
	call fprintf (fd, "\n")
end


# PH_DMAG -- Compute the integral of the magnitude between two radii
# assuming that the stellar profile can be approximated by a circular
# Moffat function.

double procedure ph_dmag (r1, r2, x, r0, params)

real	r1			# beginning radius for the integration
real	r2			# ending radius for the integration
real	x			# the airmass value
real	r0			# the hwhm of the psf
double	params[ARB]		# the model parameter array

double	bpex, pm1, x1, x2, x1sq, x2sq, d1, d2, i1, i2, dmag 

begin
	bpex = params[2] + params[5] * x
	pm1 = params[1] - 1.0d0
	x1 = (r1 / r0)
	x2 = (r2 / r0)
	x1sq = x1 ** 2
	x2sq = x2 ** 2
	x1 = x1 / params[4]
	x2 = x2 / params[4]
	d1 = 1.0d0 + r1 ** 2
	d2 = 1.0d0 + r2 ** 2
	i1 = bpex * (1.0d0 - 1.0d0 / d1 ** pm1) + (1.0d0 - bpex) *
	    (params[3] * (1.0d0 - exp (-0.5d0 * x1sq)) + (1.0d0 - params[3]) *
	    (1.0d0 - (1.0d0 + x1) * exp (-x1)))
	i2 = bpex * (1.0d0 - 1.0d0 / d2 ** pm1) + (1.0d0 - bpex) *
	    (params[3] * (1.0d0 - exp (-0.5d0 * x2sq)) + (1.0d0 - params[3]) *
	    (1.0d0 - (1.0d0 + x2) * exp (-x2)))
	dmag = -2.5d0 * log10 (i2 / i1)

	return (dmag)
end


# DMVMUL -- Multply a matrix (left-hand side) by a one dimensional vector
# (right-hand side) and return the resultant vector.

procedure dmvmul (matrix, maxdim, dim, vector, result)

double	matrix [maxdim, maxdim]	# input matrix
int	maxdim			# maximum size of input matrix
int	dim			# dimension of matrix and vectors
double	vector[maxdim]		# input vector
double	result[maxdim]		# iutput vector

double	sum
int	i, j

begin
	do i = 1, dim {
	    sum = 0.0d0
	    do j = 1, dim
		sum = sum + (matrix[j,i] * vector[j])
	    result[i] = sum
	}
end