aboutsummaryrefslogtreecommitdiff
path: root/noao/digiphot/daophot/allstar/dpalphot.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/digiphot/daophot/allstar/dpalphot.x')
-rw-r--r--noao/digiphot/daophot/allstar/dpalphot.x1438
1 files changed, 1438 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/allstar/dpalphot.x b/noao/digiphot/daophot/allstar/dpalphot.x
new file mode 100644
index 00000000..74741958
--- /dev/null
+++ b/noao/digiphot/daophot/allstar/dpalphot.x
@@ -0,0 +1,1438 @@
+include <mach.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/allstardef.h"
+
+# DP_ALPHOT -- Perform one iteration on each group of stars.
+
+int procedure dp_alphot (dao, im, ntot, istar, niter, sepcrit, sepmin, wcrit,
+ clip, clampmax, version)
+
+pointer dao # pointer to the daophot structure
+pointer im # pointer to the input image
+int ntot # the total number of stars left unfit
+int istar # first star in the current group
+int niter # the current iteration
+real sepcrit # critical separation for merging
+real sepmin # minimum separation for merging
+real wcrit # the critical error for merging
+bool clip # clip the data
+real clampmax # maximum clamping factor
+int version # version number
+
+bool regroup
+int cdimen, lstar, nstar, fstar, nterm, ixmin, ixmax, iymin, iymax
+int fixmin, fixmax, fiymin, fiymax, nxpix, nypix, flag, i, j, k, ifaint
+pointer apsel, psffit, allstar, data, subt, weights
+real radius, totradius, mean_sky, pererr, peakerr, faint
+real xmin, xmax, ymin, ymax, sumres, grpwt, chigrp
+
+bool dp_alredo(), dp_checkc(), dp_almerge(), dp_alclamp(), dp_alfaint()
+int dp_laststar(), dp_delfaintest()
+pointer dp_gwt(), dp_gst(), dp_gdc()
+real dp_almsky(), asumr()
+
+begin
+ # Get some daophot pointers.
+ psffit = DP_PSFFIT (dao)
+ apsel = DP_APSEL(dao)
+ allstar = DP_ALLSTAR (dao)
+
+ # Define some constants. At some point these should be stored
+ # in the allstar structure at task initialization. When the final
+ # rewrite gets done this will occur.
+
+ radius = DP_FITRAD(dao)
+ totradius = DP_FITRAD(dao) + DP_PSFRAD(dao) + 1
+ if (DP_RECENTER(dao) == YES)
+ cdimen = 3 * DP_MAXGROUP(dao) + 1
+ else
+ cdimen = DP_MAXGROUP(dao) + 1
+ pererr = 0.01 * DP_FLATERR(dao)
+ peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] *
+ Memr[DP_PSFPARS(psffit)+1])
+ regroup = false
+
+ repeat {
+
+ # Given the last star arrary as computed by dp_regroup, the first
+ # star in the current group, and the total number of stars,
+ # find the last star in the current group and compute the
+ # total number of stars in the group.
+
+ lstar = dp_laststar (Memi[DP_ALAST(allstar)], istar, ntot)
+ nstar = lstar - istar + 1
+
+ # Don't compute the subraster limits if no regroup has been
+ # performed.
+ if (! regroup) {
+ call alimr (Memr[DP_APXCEN(apsel)+istar-1], nstar, xmin, xmax)
+ call alimr (Memr[DP_APYCEN(apsel)+istar-1], nstar, ymin, ymax)
+ }
+
+ # Determine whether the group is too large to be fit.
+
+ if (nstar > DP_MAXGROUP (dao)) {
+
+ if (DP_VERBOSE(dao) == YES) {
+ call printf (
+ "Group too large: %5d Fitting Radius: %5.2f Regrouping\n")
+ call pargi (nstar)
+ call pargr (radius)
+ }
+
+ # If size of group too small go to next group.
+ radius = RADIUS_FRACTION * radius
+ if ((DP_RECENTER(dao) == YES) && (niter >= 2)) {
+ if (radius < DENSE_RADIUS1) {
+ fstar = dp_delfaintest (Memr[DP_APMAG(apsel)], istar,
+ lstar)
+ Memr[DP_APMAG(apsel)+fstar-1] = INDEFR
+ Memi[DP_ASKIP(allstar)+fstar-1] = YES
+ Memi[DP_AIER(allstar)+fstar-1] = ALLERR_BIGGROUP
+ if (DP_VERBOSE(dao) == YES) {
+ call printf (
+ "REJECTING: Star ID: %d Group too dense to reduce\n")
+ call pargi (Memi[DP_APID(apsel)+fstar-1])
+ }
+ return (lstar)
+ }
+ } else {
+ if (radius < DENSE_RADIUS2) {
+ fstar = dp_delfaintest (Memr[DP_APMAG(apsel)], istar,
+ lstar)
+ Memr[DP_APMAG(apsel)+fstar-1] = INDEFR
+ Memi[DP_ASKIP(allstar)+fstar-1] = YES
+ Memi[DP_AIER(allstar)+fstar-1] = ALLERR_BIGGROUP
+ if (DP_VERBOSE(dao) == YES) {
+ call printf (
+ "REJECTING: Star ID: %d Group too dense to reduce\n")
+ call pargi (Memi[DP_APID(apsel)+fstar-1])
+ }
+ return (lstar)
+ }
+ }
+
+ # Regroup the stars.
+ call dp_regroup (Memi[DP_APID(apsel)+istar-1],
+ Memr[DP_APXCEN(apsel)+istar-1], Memr[DP_APYCEN(apsel)+
+ istar-1], Memr[DP_APMAG(apsel)+istar-1],
+ Memr[DP_APMSKY(apsel)+istar-1],
+ Memr[DP_ASUMWT(allstar)+istar-1],
+ Memr[DP_AXOLD(allstar)+ istar-1],
+ Memr[DP_AYOLD(allstar)+istar-1],
+ Memr[DP_AXCLAMP(allstar)+istar-1],
+ Memr[DP_AYCLAMP(allstar)+istar-1], nstar, radius,
+ Memi[DP_ALAST(allstar)+istar-1])
+ regroup = true
+ next
+ }
+
+ # Compute the mean sky for the group. If the sky is undefined
+ # reject the group.
+ mean_sky = dp_almsky (Memr[DP_APMSKY(apsel)+istar-1], nstar)
+ if (IS_INDEFR(mean_sky)) {
+ call amovkr (INDEFR, Memr[DP_APMAG(apsel)+istar-1], nstar)
+ call amovki (YES, Memi[DP_ASKIP(allstar)+istar-1], nstar)
+ call amovki (ALLERR_INDEFSKY, Memi[DP_AIER(allstar)+istar-1],
+ nstar)
+ if (DP_VERBOSE(dao) == YES) {
+ do i = istar, lstar {
+ call printf (
+ "REJECTING: Star ID: %d Group sky value is undefined\n")
+ call pargi (Memi[DP_APID(apsel)+i-1])
+ }
+ }
+ return (lstar)
+ }
+
+ # Re-compute the number of terms in the fitting equation.
+ if ((DP_RECENTER(dao) == YES) && (niter >= 2))
+ nterm = 3 * nstar
+ else
+ nterm = nstar
+
+ # Zero the fitting arrays.
+ chigrp = asumr (Memr[DP_ASUMWT(allstar)+istar-1], nstar) / nstar
+ call aclrr (Memr[DP_APCHI(apsel)+istar-1], nstar)
+ call aclrr (Memr[DP_ASUMWT(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_ANUMER(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_ADENOM(allstar)+istar-1], nstar)
+ call aclri (Memi[DP_ANPIX(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_AV(allstar)], nterm)
+ call aclrr (Memr[DP_AC(allstar)], cdimen * cdimen)
+
+ # Compute the subraster limits.
+
+ ixmin = min (IM_LEN (im,1), max (1, int (xmin - totradius) + 1))
+ ixmax = min (IM_LEN (im,1), max (1,int (xmax + totradius)))
+ iymin = min (IM_LEN (im,2), max (1, int (ymin - totradius) + 1))
+ iymax = min (IM_LEN (im,2), max (1, int (ymax + totradius)))
+
+ # Get pointer to the required weight, scratch image and
+ # subtracted image pixels. Need to modify this so writing
+ # is only enabled if iter >= MIN_ITER.
+
+ subt = dp_gst (dao, im, iymin, iymax, READ_ONLY, NO)
+ weights = dp_gwt (dao, im, iymin, iymax, READ_WRITE, NO)
+ data = dp_gdc (dao, im, iymin, iymax, READ_WRITE, NO)
+
+ # Compute the fitting limits in the subraster.
+
+ fixmin = min (ixmax, max (ixmin, int (xmin - DP_FITRAD(dao)) + 1))
+ fixmax = min (ixmax, max (ixmin, int (xmax + DP_FITRAD(dao))))
+ fiymin = min (iymax, max (iymin, int (ymin - DP_FITRAD(dao)) + 1))
+ fiymax = min (iymax, max (iymin, int (ymax + DP_FITRAD(dao))))
+ nypix = fiymax - fiymin + 1
+ nxpix = fixmax - fixmin + 1
+
+ # Now we deal with the pixels one by one.
+ sumres = 0.0
+ grpwt = 0.0
+
+ # Accumulate the data into the matrix.
+ call dp_alaccum (dao, im, Memr[data], DP_DNX(allstar),
+ DP_DNY(allstar), DP_DXOFF(allstar), DP_DYOFF(allstar),
+ Memr[subt], DP_SNX(allstar), DP_SNY(allstar), DP_SXOFF(allstar),
+ DP_SYOFF(allstar), Memr[weights], DP_WNX(allstar),
+ DP_WNY(allstar), DP_WXOFF(allstar), DP_WYOFF(allstar), nxpix,
+ nypix, fixmin, fiymin, mean_sky, istar, lstar, niter, clip,
+ pererr, peakerr, sumres, grpwt, chigrp, cdimen, nterm)
+
+ # Reflect the normal matrix across the diagonal.
+
+ call dp_mreflect (Memr[DP_AC(allstar)], cdimen, nterm)
+
+ # Compute the estimate of the standard deviation of the residuals
+ # for the group as a whole, and for each star. This estimate starts
+ # out as SQRT (PI / 2) * [sum [weight * ABS (residual / sigma)] /
+ # sum [weight] # and then gets corrected for bias by SQRT (# of
+ # pixels / [# of pixel - # degrees of freedom]).
+ #
+ # But now we drive the value toward unity, depending upon exactly
+ # how many pixels were involved. If CHI is based on exactly a total
+ # weight of 3, then it is extremely poorly determined, and we just
+ # want to keep CHI = 1. The larger the GRPWT is, the better
+ # determined CHI is, and the less we want to force it toward unity.
+ # So, just take the weighted average of CHI and unity, with weights
+ # GRPWT = 3, and 1, respectively.
+
+ if (grpwt > nterm) {
+ chigrp = CHI_NORM * sumres / sqrt (grpwt * (grpwt - nterm))
+ chigrp = ((grpwt - 3.0) * chigrp + 3.0) / grpwt
+ } else {
+ chigrp = 1.0
+ }
+
+ # CHIGRP has been pulled towards its expected value of unity to
+ # keep the statistics of a small number of pixels from completely
+ # dominating the error analysis. Similarly, the photometric errors
+ # for the individual stars will be pulled toward unity now. Later
+ # on, if the number of stars in the group is greater than one,
+ # CHI(I) will nudged toward the group average. In order to work
+ # optimally this requires that the gain and read noise and any
+ # other parameters which represent the noise sources have been
+ # properly specified.
+ #
+ # At the same time, make sure that every star in the group contains
+ # at least MIN_FIT_PIXEL valid pixels if re-centroiding is being
+ # performed, 1 valid pixel if not. If any star in the group fails
+ # to meet this criterion, mark that star for deletion and skip
+ # ahead to the next group.
+
+ if (dp_alredo (Memi[DP_APID(apsel)], Memr[DP_APMAG(apsel)],
+ Memr[DP_APCHI(apsel)], Memr[DP_ASUMWT(allstar)],
+ Memi[DP_ANPIX(allstar)], Memi[DP_ASKIP(allstar)],
+ Memi[DP_AIER(allstar)], istar, lstar, niter, chigrp,
+ DP_RECENTER(dao), DP_VERBOSE(dao)))
+ return (lstar)
+
+ # Invert the matrix.
+ if (version == 1)
+ call invers (Memr[DP_AC(allstar)], cdimen, nterm, flag)
+ else
+ call invers2 (Memr[DP_AC(allstar)], cdimen, nterm, flag)
+
+ if (dp_checkc (Memr[DP_AC(allstar)], cdimen, nterm,
+ Memi[DP_APID(apsel)], Memr[DP_APMAG(apsel)],
+ Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)], istar, niter,
+ DP_RECENTER(dao), DP_VERBOSE(dao)))
+ return (lstar)
+
+ # Compute the coefficients.
+ call mvmul (Memr[DP_AC(allstar)], cdimen, nterm,
+ Memr[DP_AV(allstar)], Memr[DP_AX(allstar)])
+
+ # In the beginning, the brightness of each star will be permitted
+ # to change by no more than 2 magnitudes per iteration, and the
+ # x,y coordinates of each centroid will be permitted to change by
+ # no more than 0.4 pixel per iteration. Any time that the
+ # parameter correction changes sign from one iteration to the next
+ # the maximum permissable change will be reduced by a factor of
+ # two. These clamps are released any time a star in the group
+ # disappears.
+
+ if (dp_alclamp (dao, im, Memi[DP_APID(apsel)],
+ Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+ Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+ Memr[DP_ASUMWT(allstar)], Memi[DP_ASKIP(allstar)],
+ Memr[DP_AXOLD(allstar)], Memr[DP_AXCLAMP(allstar)],
+ Memr[DP_AYOLD(allstar)], Memr[DP_AYCLAMP(allstar)], istar,
+ lstar, Memr[DP_AC(allstar)], Memr[DP_AX(allstar)], cdimen,
+ niter, clampmax, pererr, peakerr)) {
+
+ # Flush the new data to the output buffers. Actually
+ # only data which fits this criterion need be written.
+ # However this makes sequential i/o management tricky.
+ # Leave this alone for the moment. Note only the cache-
+ # state is affected.
+ }
+
+ # If there is more than one star, check to see whether any two
+ # stars have merged.
+
+ if (nstar > 1) {
+
+ if (dp_almerge (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+ Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+ Memi[DP_ASKIP(allstar)], istar, lstar, sepcrit, sepmin,
+ wcrit, j, k)) {
+
+ call dp_alcentroid (Memr[DP_APXCEN(apsel)],
+ Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], j, k)
+
+ # Remove the k-th star from this group.
+ if (DP_VERBOSE(dao) == YES) {
+ call printf (
+ "REJECTING: Star ID: %d has merged with star ID: %d\n")
+ call pargi (Memi[DP_APID(apsel)+j-1])
+ call pargi (Memi[DP_APID(apsel)+k-1])
+ }
+ Memr[DP_APMAG(apsel)+j-1] = INDEFR
+ Memi[DP_ASKIP(allstar)+j-1] = YES
+ Memi[DP_AIER(allstar)+j-1] = ALLERR_MERGE
+
+ # Loosen the clamps of every star in the group
+ call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+
+ }
+ }
+
+ # If the number of iterations completed is less than or equal
+ # to 3, perform another iteration no questions asked.
+
+ if (niter <= (MIN_ITER - 1))
+ return (lstar)
+
+ # If not star has been removed from the group in the previous
+ # iteration, check to see if any of the stars is too faint (more
+ # than 12.5 magnitudes fainter thanthe PSF star). If several stars
+ # are too faint, delete the faintest one, and set the
+ # brightnesses of the other faint ones exactly to 12.5
+ # magnitudes below the PSF star. That way on the next iteration
+ # we will see whether these stars want to grow or disappear.
+
+ if (dp_alfaint (Memr[DP_APMAG(apsel)], Memi[DP_ASKIP(allstar)],
+ istar, lstar, faint, ifaint))
+ return (lstar)
+
+ # If at least one star is more than 12.5 magnitudes fainter
+ # than the PSF, then IFAINT is the index of the faintest of them
+ # and faint is the relative brightness of the faintest of them.
+
+ if (ifaint > 0) {
+
+ if (DP_VERBOSE(dao) == YES) {
+ call printf ("REJECTING: Star ID: %d is too faint\n")
+ call pargi (Memi[DP_APID(apsel)+ifaint-1])
+ }
+ Memr[DP_APMAG(apsel)+ifaint-1] = INDEFR
+ Memi[DP_ASKIP(allstar)+ifaint-1] = YES
+ Memi[DP_AIER(allstar)+ifaint-1] = ALLERR_FAINT
+
+ call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+
+ # If no star is more than 12.5 magnitudes fainter than the PSF
+ # then after the 50th iteration delete the least certain star i
+ # if it is less than a one sigma detection; after the 100th
+ # iteration delete the least certain star if it is less than a
+ # two sigma detection; after the 150th delete the least certain
+ # star if it is less than a 3 sigma detection.
+
+ } else if (niter >= 5) {
+
+ call dp_almfaint (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
+ istar, lstar, faint, ifaint)
+
+ if (faint < wcrit) {
+
+ if (DP_VERBOSE(dao) == YES) {
+ call printf ("REJECTING: Star ID: %d is too faint\n")
+ call pargi (Memi[DP_APID(apsel)+ifaint-1])
+ }
+ Memr[DP_APMAG(apsel)+ifaint-1] = INDEFR
+ Memi[DP_ASKIP(allstar)+ifaint-1] = YES
+ Memi[DP_AIER(allstar)+ifaint-1] = ALLERR_FAINT
+
+ # Loosen the clamps of every star in the group.
+ call aclrr (Memr[DP_AXOLD(allstar)+istar-1], nstar)
+ call aclrr (Memr[DP_AYOLD(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AXCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AXCLAMP(allstar)+istar-1], nstar)
+ call amaxkr (Memr[DP_AYCLAMP(allstar)+istar-1], 0.5 *
+ clampmax, Memr[DP_AYCLAMP(allstar)+istar-1], nstar)
+ }
+ }
+
+ break
+ }
+
+ return (lstar)
+end
+
+
+# DP_LASTSTAR -- Find the last star in the current group.
+
+int procedure dp_laststar (last, istar, ntot)
+
+int last[ARB] # the grouping information array
+int istar # index of the first star in the group
+int ntot # total number of stars
+
+int lstar
+
+begin
+ do lstar = istar, ntot, 1 {
+ if (last[lstar] == YES)
+ break
+ }
+
+ return (lstar)
+end
+
+
+# DP_DELFAINTEST -- Delete the faintest star in the group.
+
+int procedure dp_delfaintest (mag, istar, lstar)
+
+real mag[ARB] # the array of magnitudes
+int istar # the first star in the group
+int lstar # the last star in the group
+
+int i, starno
+real faint
+
+begin
+ starno = 0
+ faint = MAX_REAL
+ do i = istar, lstar {
+ if (mag[i] >= faint)
+ next
+ faint = mag[i]
+ starno = i
+ }
+
+ return (starno)
+end
+
+
+# DP_ALMSKY -- Determine the mean sky value for the current group of stars.
+
+real procedure dp_almsky (sky, nstar)
+
+real sky[ARB] # array of sky values
+int nstar # number of stars in the group
+
+int i, nsky
+real sky_sum
+
+begin
+ sky_sum = 0.0
+ nsky = 0
+
+ do i = 1, nstar {
+ if (IS_INDEFR(sky[i]))
+ next
+ sky_sum = sky_sum + sky[i]
+ nsky = nsky + 1
+ }
+
+ if (nsky <= 0)
+ return (INDEFR)
+ else
+ return (sky_sum / nsky)
+end
+
+
+# DP_ALACCUM -- Accumulate the data into the matrix.
+
+procedure dp_alaccum (dao, im, data, dnx, dny, dxoff, dyoff, subt, snx, sny,
+ sxoff, syoff, weights, wnx, wny, wxoff, wyoff, nxpix, nypix, ixmin,
+ iymin, mean_sky, istar, lstar, niter, clip, pererr, peakerr, sumres,
+ grpwt, chigrp, cdimen, nterm)
+
+pointer dao # pointer to the daophot structure
+pointer im # the input image descriptor
+real data[dnx,dny] # the subtracted data array
+int dnx, dny # dimenions of the data array
+int dxoff, dyoff # lower left corner of the data array
+real subt[snx,sny] # the scratch array
+int snx, sny # dimensions of the scratch array
+int sxoff, syoff # lower left corner of the scratch array
+real weights[wnx,wny] # the weight array
+int wnx, wny # dimensions of the weight array
+int wxoff, wyoff # lower left corner of the weight array
+int nxpix, nypix # the dimensions of the area of interest
+int ixmin,iymin # lower left corner of area of interest
+real mean_sky # the group sky value
+int istar # the index of the first star in current group
+int lstar # the index of the last star in current group
+int niter # the current interation
+bool clip # clip the errors ?
+real pererr # flat fielding error factor
+real peakerr # profile error factor
+real sumres # sum of the residuals
+real grpwt # the group weight
+real chigrp # the group chi value
+int cdimen # maximum number of maxtrix dimensions
+int nterm # number of terms in the matrix to fit
+
+real fitradsq, maxgdata, fix, fiy, d, delta, sigmasq, relerr, wt, dwt
+pointer psffit, apsel, allstar
+int dix, diy, sxdiff, sydiff, wxdiff, wydiff
+real sky_value, dp_alskyval()
+bool dp_alomit()
+
+begin
+ # Set up some pointers.
+ apsel = DP_APSEL(dao)
+ psffit = DP_PSFFIT(dao)
+ allstar = DP_ALLSTAR(dao)
+
+ # These constants need to be stored more permanently in the
+ # allstar structure at some point. They should all be defined
+ # once and for all at task startup. Leave for next phase
+ # of code cleanup.
+
+ fitradsq = DP_FITRAD(dao) ** 2
+ if (IS_INDEFR(DP_MAXGDATA(dao)))
+ maxgdata = MAX_REAL
+ else
+ maxgdata = DP_MAXGDATA(dao)
+
+ # Compute the array offsets.
+ sxdiff = dxoff - sxoff
+ sydiff = dyoff - syoff
+ wxdiff = dxoff - wxoff
+ wydiff = dyoff - wyoff
+
+ do diy = iymin - dyoff + 1, iymin - dyoff + nypix, 1 {
+ fiy = real (diy + dyoff - 1)
+ do dix = ixmin - dxoff + 1, ixmin - dxoff + nxpix, 1 {
+
+ # Skip data with negative weights.
+ if (weights[dix+wxdiff,diy+wydiff] < 0.0)
+ next
+ fix = real (dix + dxoff - 1)
+
+ # If this current pixel is within one fitting radius of
+ # at least one star in the current group, include it in the
+ # calculation. Otherwise skip it. While figuring this out,
+ # compute the squared distance of this pixel from the
+ # centroid of each star in the group.
+
+ if (dp_alomit (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
+ Memr[DP_ARPIXSQ(allstar)], istar, lstar, fix, fiy,
+ fitradsq))
+ next
+
+ call dp_alsetskip (Memr[DP_ARPIXSQ(allstar)],
+ Memi[DP_ASKIP(allstar)], istar, lstar, fitradsq)
+
+ if (DP_GROUPSKY(dao) == NO) {
+ sky_value = dp_alskyval (Memr[DP_APMSKY(apsel)],
+ Memi[DP_ASKIP(allstar)], istar, lstar)
+ if (IS_INDEFR(sky_value))
+ sky_value = mean_sky
+ } else
+ sky_value = mean_sky
+
+ # The expected random error in the pixel is the quadratic
+ # sum of the Poisson statistics, plus the readout noise,
+ # plus an estimated error of 0.75% of the total brightness
+ # for the difficulty of flat-fielding and bias-correcting
+ # the chip, plus an estimated error of some fraction of the
+ # of the fourth derivative at the peak of the profile, to
+ # account for the difficulty of accurately interpolating
+ # within the PSF. The fourth derivative of the PSF is
+ # proportional to H/hwhm ** 4 (hwhm is the width of the
+ # stellar core); using the geometric mean of hwhmx and
+ # hwhmy, this becomes H/(hwhmx * hwhmy) ** 2. The ratio of
+ # the fitting error to this quantity is estimated from a
+ # good-seeing CTIO frame to be approimxately 0.027.
+
+ #d = subt[dix+sxdiff,diy+sydiff] - mean_sky
+ d = subt[dix+sxdiff,diy+sydiff] - sky_value
+ delta = max (0.0, data[dix,diy] - d)
+ if ((delta > maxgdata) && (niter >= MIN_ITER))
+ next
+
+ # Dpos = raw data - residual
+ # = model-predicted brightness in this pixel consisting
+ # of sky plus all stellar profiles, which is
+ # presumably non-negative
+ #
+ # The four noise sources in the model are
+ # readout noise
+ # poisson noise
+ # flat-field errors
+ # profile errors
+ #
+ # Numerically the squares of these quantities are
+ # ronoise = sigma[i,j]
+ # poisson noise = delta / gain
+ # flat-field error = constant * delta ** 2
+ # profile error = constant * sum of profile ** 2
+
+ #sigmasq = weights[dix+wxdiff,diy+wydiff] + delta /
+ #DP_PHOTADU(dao) + (pererr * delta) ** 2 +
+ #(peakerr * (delta - mean_sky)) ** 2
+ sigmasq = weights[dix+wxdiff,diy+wydiff] + delta /
+ DP_PHOTADU(dao) + (pererr * delta) ** 2 +
+ (peakerr * (delta - sky_value)) ** 2
+ if (sigmasq <= 0.0)
+ next
+ relerr = abs (d) / sqrt (sigmasq)
+
+ if (clip && (relerr > MAX_RELERR))
+ next
+
+ # Now include this pixel in the fitting equation for the
+ # group.
+
+ wt = 0.0
+ call dp_alxaccum (dao, im, Memr[DP_APXCEN(apsel)],
+ Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
+ Memr[DP_ARPIXSQ(allstar)], Memi[DP_ASKIP(allstar)],
+ Memr[DP_AX(allstar)], Memr[DP_ANUMER(allstar)],
+ Memr[DP_ADENOM(allstar)], istar, lstar, fix, fiy, d,
+ wt, sigmasq, nterm, fitradsq)
+
+ # At this point the vector X contains the first
+ # derivative of the condition equation for pixel (I,J)
+ # with respect to each of the fitting parameters for
+ # all of the stars. Now these derivatives will be added
+ # into the normal matrix and the vector of residuals.
+ # Add this residual into the weighted sum of the absolute
+ # relative residuals
+
+ dwt = wt * relerr
+ sumres = sumres + dwt
+ grpwt = grpwt + wt
+
+ # SUMRES is the weighted sum for all the pixels in the group.
+ # Now also add the weigghted value of the residuals into the
+ # accumulating sum for each of the stars.
+
+ call dp_alchi (Memi[DP_ASKIP(allstar)], Memr[DP_APCHI(apsel)],
+ Memr[DP_ASUMWT(allstar)], Memi[DP_ANPIX(allstar)], istar,
+ lstar, wt, dwt)
+
+ # Up until now, WT represents only the radial weighting
+ # profile. Now figure in the anticipated standard error
+ # of the pixel.
+
+ wt = wt / sigmasq
+
+ # After the third iteration, reduce the weight of a bad pixel
+ # Note that for the first iteration, only the stellar magnitude
+ # is being solved for, which is a problem in LINEAR least
+ # squares, and so should be solved exactly the first time.
+ # After that, the star is given two iterations to adjust it's
+ # centroid before the clipping is turned on. After that a
+ # pixel having a residual of DP_CLIPRANGE times sigma gets
+ # reduced to half weight; a pixel having a residual of n
+ # sigma gets a weight of 1 / [1 + (n/DP_CLIPRANGE) **
+ # DP_CLIPEXP].
+
+ if (clip)
+ wt = wt / (1.0 + (relerr / (chigrp * DP_CLIPRANGE (dao))) **
+ DP_CLIPEXP (dao))
+
+ # Now work this pixel into the normal matrix
+ dwt = d * wt
+ call dp_alcaccum (Memr[DP_AX(allstar)], Memr[DP_AC(allstar)],
+ Memr[DP_AV(allstar)], Memi[DP_ASKIP(allstar)], cdimen,
+ nterm, istar, lstar, wt, dwt)
+ }
+ }
+end
+
+
+# DP_ALOMIT - Omit pixels which are too far from the center of any star
+# in the group.
+
+bool procedure dp_alomit (xcen, ycen, rpixsq, istar, lstar, fix, fiy,
+ fitradsq)
+
+real xcen[ARB] # array of x centers
+real ycen[ARB] # array of y centers
+real rpixsq[ARB] # array of pixel differences
+int istar, lstar # first and last star
+real fix, fiy # current star position
+real fitradsq # fit radius squared
+
+bool omit
+int i
+
+begin
+ omit = true
+ do i = istar, lstar {
+ rpixsq[i] = (fix - xcen[i]) ** 2 + (fiy - ycen[i]) ** 2
+ if (rpixsq[i] < fitradsq)
+ omit = false
+ }
+
+ return (omit)
+end
+
+
+# DP_ALSETSKIP -- Initialize the skip array.
+
+procedure dp_alsetskip (rpixsq, skip, istar, lstar, fitradsq)
+
+real rpixsq[ARB] # pixel distance squared
+int skip[ARB] # skip array
+int istar # the index of the first star
+int lstar # the index of the last star
+real fitradsq # the fitting radius squared
+
+int i
+
+begin
+ do i = istar, lstar {
+ if (rpixsq[i] < fitradsq)
+ skip[i] = NO
+ else
+ skip[i] = YES
+ }
+end
+
+
+# DP_ALSKYVAL -- Compute the sky value to be subtracted from a given pixel
+# by averaging the sky values of all the the stars for which the pixel in
+# question is within one fitting radius.
+
+real procedure dp_alskyval (sky, skip, istar, lstar)
+
+real sky[ARB] # the array of sky values
+int skip[ARB] # the array of skip values
+int istar # the index of the first star
+int lstar # the index of the last star
+
+int i, npts
+real sum
+
+begin
+ sum = 0.0
+ npts = 0
+ do i = istar, lstar {
+ if (skip[i] == YES)
+ next
+ sum = sum + sky[i]
+ npts = npts + 1
+ }
+
+ if (npts <= 0)
+ return (INDEFR)
+ else
+ return (sum / npts)
+end
+
+
+# DP_ALXACCUM - Accumulate the x vector.
+
+procedure dp_alxaccum (dao, im, xcen, ycen, mag, rpixsq, skip, x, numer1,
+ denom1, istar, lstar, fix, fiy, resid, wt, sigmasq, nterm, fitradsq)
+
+pointer dao # pointer to the daopot structure
+pointer im # pointer to the input image
+real xcen[ARB] # array of x centers
+real ycen[ARB] # array of y centers
+real mag[ARB] # array of relative brightnesses
+real rpixsq[ARB] # array of pixel distances squared
+int skip[ARB] # the include / exclude array
+real x[ARB] # the x vector to be accumulated
+real numer1[ARB] # first numerator array
+real denom1[ARB] # first denominator array
+int istar # index of the first star in group
+int lstar # index of the last star in group
+real fix # the x position of the current pixel
+real fiy # the y position of the current pixel
+real resid # the input data residual
+real wt # the output weight value
+real sigmasq # the sigma squared value
+int nterm # number of terms in the matrix
+real fitradsq # fitting radius squared
+
+real psfsigsqx, psfsigsqy, dx, dy, deltax, deltay, value, radsq, rhosq
+real dvdx, dvdy, dfdsig
+pointer psffit
+int nstar, i, i3, k
+real dp_usepsf()
+
+begin
+ psffit = DP_PSFFIT(dao)
+
+ nstar = lstar - istar + 1
+ psfsigsqx = Memr[DP_PSFPARS(psffit)]
+ psfsigsqy = Memr[DP_PSFPARS(psffit)+1]
+
+ do i = istar, lstar {
+
+ if (skip[i] == YES)
+ next
+ radsq = rpixsq[i] / fitradsq
+ if (radsq >= MAX_RSQ)
+ next
+ wt = max (wt, 5.0 / (5.0 + radsq / (1.0 - radsq)))
+
+ # The condition equation for pixel[I,J] is the following.
+ #
+ # data[I,J] - sum [scale * psf[I,J]] - mean_sky = residual
+ #
+ # Then we will adjust the scale factors, x centers and ycenters
+ # such that sum [weight * residual ** 2] is minimized.
+ # WT will be a function (1) of the distance of this pixel from
+ # the center of the nearest star, (2) of the model predicted
+ # brightness of the pixel (taking into consideration the
+ # readout noise, the photons/ADU, and the interpolation error
+ # of the PSF), and (3) of the size of the residual itself. (1) is
+ # necessary to prevent the non-linear least-squares from
+ # fit from oscillating. (2) is simply sensible weighting and (3)
+ # is a crude attempt at making the solution more robust
+ # against bad pixels.
+
+ dx = fix - xcen[i]
+ dy = fiy - ycen[i]
+ call dp_wpsf (dao, im, xcen[i], ycen[i], deltax, deltay, 1)
+ deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+ deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+ value = dp_usepsf (DP_PSFUNCTION(psffit), dx, dy,
+ DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+ Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit),
+ DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), deltax, deltay,
+ dvdx, dvdy)
+
+ if (nterm > nstar) {
+ i3 = (i - istar + 1) * 3
+ k = i3 - 2
+ x[k] = -value
+ k = i3 - 1
+ x[k] = -mag[i] * dvdx
+ x[i3] = -mag[i] * dvdy
+ } else {
+ k = i - istar + 1
+ x[k] = -value
+ }
+
+ rhosq = (dx / psfsigsqx) ** 2 + (dy / psfsigsqy) ** 2
+ if (rhosq > MAX_RHOSQ)
+ next
+ rhosq = 0.6931472 * rhosq
+ dfdsig = exp (-rhosq) * (rhosq - 1.)
+ numer1[i] = numer1[i] + dfdsig * resid / sigmasq
+ denom1[i] = denom1[i] + dfdsig ** 2 / sigmasq
+ }
+end
+
+
+# DP_ALCHI -- Compute the weights and chis.
+
+procedure dp_alchi (skip, chi, sumwt, npix, istar, lstar, wt, dwt)
+
+int skip[ARB] # include / exclude array for the stars
+real chi[ARB] # array of chi values
+real sumwt[ARB] # array of weight sums
+int npix[ARB] # array of pixel counts
+int istar # index of first star
+int lstar # index of last star
+real wt # input weight
+real dwt # weighted residual
+
+int i
+
+begin
+ do i = istar, lstar {
+ if (skip[i] == YES)
+ next
+ chi[i] = chi[i] + dwt
+ sumwt[i] = sumwt[i] + wt
+ npix[i] = npix[i] + 1
+ }
+end
+
+
+# DP_ALCACCUM -- Accumulate the main matrix.
+
+procedure dp_alcaccum (x, c, v, skip, cdimen, nterm, istar, lstar, wt, dwt)
+
+
+real x[ARB] # x array
+real c[cdimen,ARB] # maxtrix array
+real v[ARB] # v vector
+int skip[ARB] # include / exclude array
+int cdimen # maximum size of c matrix
+int nterm # number of terms to be fit
+int istar # index of the first star
+int lstar # index of the last star
+real wt # input weight
+real dwt # input residual weight
+
+int nstar, i, j, k, l, i3, i3m2, j3
+real xkwt
+
+begin
+ nstar = lstar - istar + 1
+ do i = istar, lstar {
+ if (skip[i] == YES)
+ next
+ if (nterm > nstar) {
+ i3 = (i - istar + 1) * 3
+ i3m2 = i3 - 2
+ do k = i3m2, i3
+ v[k] = v[k] + x[k] * dwt
+ do j = istar, i {
+ if (skip[j] == YES)
+ next
+ j3 = (j - istar + 1) * 3
+ do k = i3m2, i3 {
+ xkwt = x[k] * wt
+ do l = j3 - 2, min (k, j3)
+ c[k,l] = c[k,l] + x[l] * xkwt
+ }
+ }
+ } else {
+ k = i - istar + 1
+ v[k] = v[k] + x[k] * dwt
+ xkwt = x[k] * wt
+ do j = istar, i {
+ if (skip[j] == YES)
+ next
+ l = j - istar + 1
+ c[k,l] = c[k,l] + x[l] * xkwt
+ }
+ }
+ }
+end
+
+
+# DP_ALREDO -- Check to see that there are enough good pixels to fit the
+# star and compute the individual chi values.
+
+bool procedure dp_alredo (id, mag, chi, sumwt, npix, skip, aier, istar, lstar,
+ niter, chigrp, center, verbose)
+
+int id[ARB] # the array of star ids
+real mag[ARB] # the array of relative brightnesses
+real chi[ARB] # the array of chi values
+real sumwt[ARB] # the array of weight values
+int npix[ARB] # the array of pixel counts
+int skip[ARB] # the include / exclude array
+int aier[ARB] # the array of error codes
+int istar # index of the first star
+int lstar # index of the last star
+int niter # the current iteration
+real chigrp # chi value of the group
+int center # recenter the values ?
+int verbose # verbose mode ?
+
+bool redo
+int i
+
+begin
+ redo = false
+ do i = istar, lstar {
+ if (center == YES && (niter >= 2)) {
+ if (npix[i] < MIN_NPIX) {
+ redo = true
+ skip[i] = YES
+ if (verbose == YES) {
+ call printf (
+ "REJECTING: Star ID: %d has too few good pixels\n")
+ call pargi (id[i])
+ }
+ aier[i] = ALLERR_NOPIX
+ mag[i] = INDEFR
+ } else {
+ skip[i] = NO
+ if (sumwt[i] > MIN_SUMWT) {
+ chi[i] = CHI_NORM * chi[i] / sqrt (sumwt[i] *
+ (sumwt[i] - MIN_SUMWT))
+ sumwt[i] = ((sumwt[i] - MIN_SUMWT) * chi[i] +
+ MIN_SUMWT) / sumwt[i]
+ } else
+ chi[i] = chigrp
+ }
+ } else {
+ if (npix[i] < MIN_NPIX) {
+ redo = true
+ skip[i] = YES
+ if (verbose == YES) {
+ call printf (
+ "REJECTING: Star ID: %d has too few good pixels\n")
+ call pargi (id[i])
+ }
+ aier[i] = ALLERR_NOPIX
+ mag[i] = INDEFR
+ } else {
+ skip[i] = NO
+ if (sumwt[i] > 1.0) {
+ chi[i] = CHI_NORM * chi[i] / sqrt (sumwt[i] *
+ (sumwt[i] - 1.0))
+ sumwt[i] = ((sumwt[i] - MIN_SUMWT) * chi[i] +
+ MIN_SUMWT) / sumwt[i]
+ } else
+ chi[i] = chigrp
+ }
+ }
+ }
+
+ return (redo)
+end
+
+
+# DP_CHECKC - Check the c matrix for valid diagonal elements.
+
+bool procedure dp_checkc (c, cdimen, nterm, id, mag, skip, aier, istar, niter,
+ recenter, verbose)
+
+real c[cdimen,ARB] # the c matrix
+int cdimen # maximum size of the c matrix
+int nterm # number of terms to be fit
+int id[ARB] # the array of ids
+real mag[ARB] # the array of relative brightnesses
+int skip[ARB] # the include / exclude array
+int aier[ARB] # the array of error codes
+int istar # index of the first star
+int niter # the iteration number
+int recenter # recenter ?
+int verbose # verbose mode ?
+
+bool redo
+int j, starno
+
+begin
+ redo = false
+ do j = 1, nterm {
+ if (c[j,j] > 0.0)
+ next
+ if ((recenter == YES) && (niter >= 2))
+ starno = (j + 2) / 3
+ else
+ starno = j
+ starno = starno + istar - 1
+ skip[starno] = YES
+ if (verbose == YES) {
+ call printf (
+ "REJECTING: Star ID: %d generated a fitting error\n")
+ call pargi (id[starno])
+ }
+ aier[starno] = ALLERR_SINGULAR
+ mag[starno] = INDEFR
+ redo = true
+ break
+ }
+
+ return (redo)
+end
+
+
+# DP_ALCLAMP -- Get the answers.
+
+bool procedure dp_alclamp (dao, im, id, xcen, ycen, mag, magerr, sumwt,
+ skip, xold, xclamp, yold, yclamp, istar, lstar, c, x,
+ cdimen, niter, clampmax, pererr, peakerr)
+
+pointer dao # pointer to the daophot structure
+pointer im # the input image descriptor
+int id[ARB] # array of star ids
+real xcen[ARB] # array of star x centers
+real ycen[ARB] # array of star y centers
+real mag[ARB] # array of relative brightnesses
+real magerr[ARB] # array of relative brightness errors
+real sumwt[ARB] # array of pixel distance squared values
+int skip[ARB] # include / exclude array
+real xold[ARB] # xold array
+real xclamp[ARB] # xclamp array
+real yold[ARB] # yold array
+real yclamp[ARB] # yclamp array
+int istar # the index of the first star
+int lstar # the index of the last star
+real c[cdimen,ARB] # c matrix
+real x[ARB] # x vector
+int cdimen # the maximum size of c matrix
+int niter # the current iteration
+real clampmax # the maximum clamp value
+real pererr # flat field error
+real peakerr # the profile error
+
+bool redo, bufwrite
+int i, l, k, j, lx, mx, ly, my, nx, ny
+pointer psffit, allstar
+real dwt, psfrad, psfradsq, df
+
+begin
+ # Get some daophot pointers.
+ allstar = DP_ALLSTAR(dao)
+ psffit = DP_PSFFIT(dao)
+
+ # Get some constants.
+ if (DP_PSFSIZE(psffit) == 0)
+ psfrad = DP_PSFRAD(dao)
+ else
+ psfrad = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0
+ psfradsq = psfrad * psfrad
+
+ # Initialize.
+ bufwrite = false
+
+ do i = istar, lstar {
+
+ if ((DP_RECENTER(dao) == YES) && (niter >= 2)) {
+
+ l = 3 * (i - istar + 1)
+ k = l - 1
+ j = l - 2
+
+ # Note that the sign of the correction is such that it must be
+ # SUBTRACTED from the current value of the parameter to get the
+ # improved parameter value. This being the case, if the
+ # correction to the brightness is negative a change of one
+ # magnitude is a change of factor of 2.5; if the brightness
+ # correction is positive a change of one magnitude is a change
+ # of 60%.
+
+ dwt = xold[i] * x[k]
+ if (dwt < 0.0)
+ xclamp[i] = max (MIN_XYCLAMP, MIN_XYCLAMP_FRACTION *
+ xclamp[i])
+ else
+ xclamp[i] = min (clampmax, MAX_XYCLAMP_FRACTION * xclamp[i])
+ xcen[i] = xcen[i] - x[k] / (1.0 + abs (x[k] / xclamp[i]))
+ xold[i] = x[k]
+
+ dwt = yold[i] * x[l]
+ if (dwt < 0.0)
+ yclamp[i] = max (MIN_XYCLAMP, MIN_XYCLAMP_FRACTION *
+ yclamp[i])
+ else
+ yclamp[i] = min (clampmax, MAX_XYCLAMP_FRACTION * yclamp[i])
+ ycen[i] = ycen[i] - x[l] / (1.0 + abs (x[l] / yclamp[i]))
+ yold[i] = x[l]
+
+ mag[i] = mag[i] - x[j] / (1.0 + max (x[j] / (MAX_DELTA_FAINTER *
+ mag[i]), -x[j] / (MAX_DELTA_BRIGHTER * mag[i])))
+ magerr[i] = sumwt[i] * sqrt (c[j,j])
+
+ if (niter >= MIN_ITER) {
+ redo = false
+ if (abs(x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+ MAX_NEW_RELBRIGHT2 * mag[i]))
+ redo = true
+ else {
+ df = (MAX_NEW_ERRMAG * sumwt[i]) ** 2
+ if ((x[k] ** 2) > max (df * c[k,k], MAX_PIXERR2))
+ redo = true
+ else if ((x[l] ** 2) > max (df * c[l,l], MAX_PIXERR2))
+ redo = true
+ }
+ } else
+ redo = true
+
+ } else {
+ j = i - istar + 1
+ mag[i] = mag[i] - x[j] / (1.0 + 1.2 * abs (x[j] / mag[i]))
+ magerr[i] = sumwt[i] * sqrt (c[j,j])
+ if (niter >= 2) {
+ redo = false
+ if (abs(x[j]) > max (MAX_NEW_ERRMAG * magerr[i],
+ MAX_NEW_RELBRIGHT2 * mag[i]))
+ redo = true
+ } else
+ redo = true
+ }
+
+ if (mag[i] < 2.0 * magerr[i])
+ redo = true
+ if (niter >= DP_MAXITER(dao))
+ redo = false
+ if (redo && (niter < DP_MAXITER(dao)))
+ next
+
+ # If this star converged, write out the results for it,
+ # flag it for deletion from the star list and subtract it
+ # from the reference copy of the image.
+
+ call dp_glim (xcen[i], ycen[i], psfrad, DP_DLX(allstar),
+ DP_DMX(allstar), DP_DLY(allstar), DP_DMY(allstar),
+ lx, mx, ly, my)
+ nx = mx - lx + 1
+ ny = my - ly + 1
+
+ call dp_swstar (dao, im, Memr[DP_DBUF(allstar)], DP_DNX(allstar),
+ DP_DNY(allstar), DP_DXOFF(allstar), DP_DYOFF(allstar),
+ Memr[DP_WBUF(allstar)], DP_WNX(allstar), DP_WNY(allstar),
+ DP_WXOFF(allstar), DP_WYOFF(allstar), xcen[i], ycen[i],
+ mag[i], lx, ly, nx, ny, psfradsq, DP_PHOTADU(dao), pererr,
+ peakerr)
+
+ skip[i] = YES
+
+ if (DP_VERBOSE(dao) == YES) {
+ call dp_wout (dao, im, xcen[i], ycen[i], xcen[i], ycen[i], 1)
+ call printf (
+ "FITTING: ID: %5d XCEN: %8.2f YCEN: %8.2f MAG: %8.2f\n")
+ call pargi (id[i])
+ call pargr (xcen[i])
+ call pargr (ycen[i])
+ if (mag[i] <= 0.0)
+ call pargr (INDEFR)
+ else
+ call pargr (DP_PSFMAG(psffit) - 2.5 * log10 (mag[i]))
+
+ }
+ bufwrite = true
+ }
+
+ return (bufwrite)
+end
+
+
+# DP_SWSTAR -- Subtract the fitted star for the data and adjust the
+# weights.
+
+procedure dp_swstar (dao, im, data, dnx, dny, dxoff, dyoff, weights, wnx, wny,
+ wxoff, wyoff, xcen, ycen, mag, lx, ly, nx, ny, psfradsq, gain,
+ pererr, peakerr)
+
+pointer dao # pointer to the daophot structure
+pointer im # pointer to the input image
+real data[dnx,dny] # the data array
+int dnx, dny # dimensions of the data array
+int dxoff, dyoff # lower left corner of data array
+real weights[wnx,wny] # the weights array
+int wnx, wny # dimensions of the weights array
+int wxoff, wyoff # lower left corner of weights array
+real xcen, ycen # the position of the star
+real mag # relative brightness of the star
+int lx, ly # lower left corner region of interest
+int nx, ny # size of region of interest
+real psfradsq # the psf radius squared
+real gain # the gain in photons per adu
+real pererr # the flat field error factor
+real peakerr # the peak error factor
+
+real deltax, deltay, dx, dy, dysq, diff, dvdx, dvdy
+pointer psffit
+int di, dj, wxdiff, wydiff
+real dp_usepsf()
+
+begin
+ psffit = DP_PSFFIT(dao)
+
+ wxdiff = dxoff - wxoff
+ wydiff = dyoff - wyoff
+ call dp_wpsf (dao, im, xcen, ycen, deltax, deltay, 1)
+ deltax = (deltax - 1.0) / DP_PSFX(psffit) - 1.0
+ deltay = (deltay - 1.0) / DP_PSFY(psffit) - 1.0
+
+ do dj = ly - dyoff + 1, ly - dyoff + ny {
+ dy = (dj + dyoff - 1) - ycen
+ dysq = dy * dy
+ do di = lx - dxoff + 1, lx - dxoff + nx {
+ if (weights[di+wxdiff,dj+wydiff] <= -MAX_REAL)
+ next
+ dx = (di + dxoff - 1) - xcen
+ if ((dx * dx + dysq) >= psfradsq)
+ next
+ diff = mag * dp_usepsf (DP_PSFUNCTION(psffit), dx, dy,
+ DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
+ Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit),
+ DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), deltax, deltay,
+ dvdx, dvdy)
+ data[di,dj] = data[di,dj] - diff
+ if (diff > 0.0) {
+ diff = diff / gain + pererr *
+ 2.0 * max (0.0, data[di,dj]) * diff +
+ (peakerr * diff) ** 2
+ if (weights[di+wxdiff,dj+wydiff] >= 0.0)
+ weights[di+wxdiff,dj+wydiff] = weights[di+wxdiff,
+ dj+wydiff] + diff
+ else
+ weights[di+wxdiff,dj+wydiff] = - (abs (weights[di+
+ wxdiff,dj+wydiff]) + diff)
+ }
+ }
+ }
+end
+
+
+# DP_ALMERGE -- Determine whether or not two stars should merge.
+
+bool procedure dp_almerge (xcen, ycen, mag, magerr, skip, istar, lstar,
+ sepcrit, sepmin, wcrit, k, l)
+
+real xcen[ARB] # array of x centers
+real ycen[ARB] # array of y centers
+real mag[ARB] # array of relative brightnesses
+real magerr[ARB] # array of relative brightness errors
+int skip[ARB] # the include / exclude array
+int istar # index to the first star
+int lstar # index to the last star
+real sepcrit # the critical separation
+real sepmin # the minimum separation
+real wcrit # the critical error
+int k, l # indices of stars to be merged
+
+bool merge
+int i, j
+real rsq, sep
+
+begin
+ # Initialize.
+ rsq = sepcrit
+ merge = false
+ k = 0
+ l = 0
+
+ # Find the closest two stars within the critical separation.
+ do i = istar + 1, lstar {
+
+ if (skip[i] == YES)
+ next
+
+ do j = istar, i - 1 {
+ if (skip[j] == YES)
+ next
+ sep = (xcen[i] - xcen[j]) ** 2 + (ycen[i] - ycen[j]) ** 2
+ if (sep >= rsq)
+ next
+
+ # Two stars are overlapping. Identify the fainter of the two.
+ rsq = sep
+ if (mag[i] < mag[j]) {
+ k = i
+ l = j
+ } else {
+ k = j
+ l = i
+ }
+
+ merge = true
+ }
+ }
+
+ # No stars overlap.
+ if (! merge)
+ return (merge)
+
+ # The stars are not close enough.
+ if ((rsq > sepmin) && (mag[k] > wcrit * magerr[k]))
+ merge = false
+
+ return (merge)
+end
+
+
+# DP_ALCENTROID -- Merge two stars by adding their relative brightnesses
+# and replacing their individual positions with the relative brightness
+# weighted mean position.
+
+procedure dp_alcentroid (xcen, ycen, mag, k, l)
+
+real xcen[ARB] # array of x centers
+real ycen[ARB] # array of y centers
+real mag[ARB] # array of magnitudes
+int k, l # input indices, k is fainter
+
+begin
+ xcen[l] = xcen[l] * mag[l] + xcen[k] * mag[k]
+ ycen[l] = ycen[l] * mag[l] + ycen[k] * mag[k]
+ mag[l] = mag[l] + mag[k]
+ xcen[l] = xcen[l] / mag[l]
+ ycen[l] = ycen[l] / mag[l]
+end
+
+
+# DP_ALFAINT -- Find all the stars with realtive brightnesss less than
+# a given number and set their magnitudes to a given minimum relative
+# brightness.
+
+bool procedure dp_alfaint (mag, skip, istar, lstar, faint, ifaint)
+
+real mag[ARB] # array of magnitudes
+int skip[ARB] # skip array
+int istar, lstar # first and last stars
+real faint # faintest star
+int ifaint # index of faintest star
+
+int i
+
+begin
+ # Initialize
+ faint = MIN_REL_BRIGHT
+ ifaint = 0
+
+ do i = istar, lstar {
+ if (skip[i] == YES)
+ return (true)
+ if (mag[i] < faint) {
+ faint = mag[i]
+ ifaint = i
+ }
+ if (mag[i] < MIN_REL_BRIGHT)
+ mag[i] = MIN_REL_BRIGHT
+ }
+
+ return (false)
+end
+
+
+# DP_ALMFAINT -- Find the star with the greatest error in its relative
+# brightness.
+
+procedure dp_almfaint (mag, magerr, istar, lstar, faint, ifaint)
+
+real mag[ARB] # array of relative brightnesses
+real magerr[ARB] # array of relative brightness errors
+int istar # index of first star
+int lstar # index of last star
+real faint # computed faintness limit
+int ifaint # index of faintest star
+
+int i
+real wt
+
+begin
+ faint = MAX_REAL
+ ifaint = 0
+ do i = istar, lstar {
+ wt = mag[i] / magerr[i]
+ if (wt < faint) {
+ faint = wt
+ ifaint = i
+ }
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 12:06:24 -0400