include include 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