aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/photcal/mkobsfile/phagrow.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/photcal/mkobsfile/phagrow.x')
-rw-r--r--noao/digiphot/photcal/mkobsfile/phagrow.x1982
1 files changed, 1982 insertions, 0 deletions
diff --git a/noao/digiphot/photcal/mkobsfile/phagrow.x b/noao/digiphot/photcal/mkobsfile/phagrow.x
new file mode 100644
index 00000000..a8c6dfce
--- /dev/null
+++ b/noao/digiphot/photcal/mkobsfile/phagrow.x
@@ -0,0 +1,1982 @@
+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
generated by cgit (git 2.34.1) at 2025-09-20 08:59:40 -0400