include include include "../lib/daophotdef.h" include "../lib/apseldef.h" include "../lib/nstardef.h" # DP_NSTARFIT -- Fit the stellar group. int procedure dp_nstarfit (dao, im, nin_group, mean_sky, cdimen, iter, converge) pointer dao # pointer to the daophot structure pointer im # pointer to the input image int nin_group # original group size real mean_sky # the mean sky for the group int cdimen # dimensions of the coefficient matrix int iter # the current iteration bool converge # did the fit converge ? bool clip, refit int i, j, k, xpix, ypix, group_size, nterm, ncols, mindex, ifaint, tifaint int ntmin, ixmin, ixmax, iymin, iymax, ifxmin, ifxmax, ifymin, ifymax, flag pointer apsel, psffit, nstar, subim, ypixel, pixel real fitradsq, cutoff, psfradsq, sepcrit, sepmin, perr, peakerr, sky_value real read_noise, mingdata, maxgdata, chigrp, wcrit, xmin, xmax, ymin, ymax real datum, sumres, grpwt, xtemp, ytemp, weight, ds, pred_pixval, sigmasq real relerr, dswt, faint, tfaint, noise bool dp_nstmerge(), dp_ntomit(), dp_ntmin(), dp_ncheckc() pointer imgs2r() real dp_ntskyval(), dp_ntsubtract() begin # Define the daophot pointers. psffit = DP_PSFFIT (dao) apsel = DP_APSEL(dao) nstar = DP_NSTAR (dao) # Set up some daophot constants. At some point these will be computed # when the NSTAR task is started up instead of at the beginning of # each group fit. For the moment it is convenient and not too # costly to compute them here. Also initialize the fit. if (iter == 1) { # Compute the fitting and psf radii. fitradsq = DP_FITRAD (dao) ** 2 psfradsq = DP_PSFRAD(dao) ** 2 cutoff = CUT_FACTOR * fitradsq # Compute the merging parameters. if (IS_INDEFR(DP_MERGERAD(dao))) { sepcrit = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 + Memr[DP_PSFPARS(psffit)+1] ** 2) sepmin = min (1.0, FRACTION_MINSEP * sepcrit) } else { sepcrit = DP_MERGERAD(dao) ** 2 sepmin = min (1.0, FRACTION_MINSEP * sepcrit) } # Compute the noise parameters. read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2 perr = 0.01 * DP_FLATERR(dao) peakerr = 0.01 * DP_PROFERR(dao) / (Memr[DP_PSFPARS(psffit)] * Memr[DP_PSFPARS(psffit)+1]) # Compute the minimum and maximum good pixel values. if (IS_INDEFR(DP_MINGDATA(dao))) mingdata = -MAX_REAL else mingdata = DP_MINGDATA(dao) if (IS_INDEFR(DP_MAXGDATA(dao))) maxgdata = MAX_REAL else maxgdata = DP_MAXGDATA(dao) # Initialize the fit. chigrp = 1.0 clip = false if (DP_RECENTER(dao) == YES) { nterm = 3 * nin_group ntmin = 3 } else { nterm = nin_group ntmin = 1 } if (DP_FITSKY(dao) == YES) { nterm = nterm + 1 ntmin = ntmin + 1 } call aclrr (Memr[DP_NXOLD(nstar)], nterm) call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) } # Start a new iteration. group_size = nin_group if (DP_RECENTER(dao) == YES) nterm = 3 * group_size else nterm = group_size if (DP_FITSKY(dao) == YES) nterm = nterm + 1 # Begin fitting the current group of stars. repeat { # Initialize the convergence criteria. converge = false # Check that there is at least 1 star in the group. if (group_size < 1) return (group_size) # Set up the critical error for star rejection. if (iter >= NITER_MAX) wcrit = WCRIT_MAX else if (iter >= NITER_MED) wcrit = WCRIT_MED else if (iter >= NITER_MIN) wcrit = WCRIT_MIN else wcrit = MAX_REAL # Set the sky fitting derivative. if (DP_FITSKY(dao) == YES) Memr[DP_NX(nstar)+nterm-1] = -1.0 # Initialize arrays. call aclrr (Memr[DP_APCHI(apsel)], group_size) call aclrr (Memr[DP_NSUMWT(nstar)], group_size) call aclrr (Memr[DP_NNUMER(nstar)], group_size) call aclrr (Memr[DP_NDENOM(nstar)], group_size) call amovki (NSTERR_OK, Memi[DP_NIER(nstar)], group_size) # Compute the minimum and maximum x and y values. call alimr (Memr[DP_APXCEN(apsel)], group_size, xmin, xmax) call alimr (Memr[DP_APYCEN(apsel)], group_size, ymin, ymax) # Check to see whether any two stars are within the critical # difference from each other. if ((group_size > 1) && dp_nstmerge (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)], group_size, sepcrit, sepmin, wcrit, i, j, k)) { # Compute the new centroid and brightness. call dp_nstcen (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], i, j, k) # Print out verbose comments. if (DP_VERBOSE (dao) == YES) { call printf ("\tStar %-5d has merged with star %-5d\n") call pargi (Memi[DP_APID(apsel)+k-1]) call pargi (Memi[DP_APID(apsel)+i-1]) } # Recompute the mean sky. if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES)) mean_sky = (mean_sky * group_size - Memr[DP_APMSKY(apsel)+k-1]) / (group_size - 1) # Now remove the k-th star from the group. call dp_remove (k, group_size, Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], NSTERR_MERGE) # After deleting a star, resize the matrix, release all of # the clamps and back up the iteration counter. if (DP_RECENTER(dao) == YES) nterm = 3 * group_size else nterm = group_size if (DP_FITSKY(dao) == YES) nterm = nterm + 1 clip = false call aclrr (Memr[DP_NXOLD(nstar)], nterm) call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) iter = max (1, iter - 1) next } # Now we can proceed with the iteration. If this is the first # iteration read in the subraster. Determine the size of the # subraster we need to extract from the image for this group. if (iter == 1) { ixmin = max (1, int (xmin - DP_PSFRAD(dao) - DP_FITRAD(dao)) + 1) iymin = max (1, int (ymin - DP_PSFRAD(dao) - DP_FITRAD(dao)) + 1) ixmax = min (IM_LEN(im,1), int (xmax + DP_PSFRAD(dao) + DP_FITRAD(dao))) iymax = min (IM_LEN(im,2), int (ymax + DP_PSFRAD(dao) + DP_FITRAD(dao))) subim = imgs2r (im, ixmin, ixmax, iymin, iymax) ncols = ixmax - ixmin + 1 } # Compute the area on the subraster that is off interest to # the current iteration. ifxmin = max (ixmin, int (xmin - DP_FITRAD(dao)) + 1) ifymin = max (iymin, int (ymin - DP_FITRAD(dao)) + 1) ifxmax = min (ixmax, int (xmax + DP_FITRAD(dao))) ifymax = min (iymax, int (ymax + DP_FITRAD(dao))) # Zero the normal matrix and the vector of residuals. call aclrr (Memr[DP_NV(nstar)], nterm) call aclrr (Memr[DP_NC(nstar)], cdimen * cdimen) call aclri (Memi[DP_NNPIX(nstar)], group_size) sumres = 0.0 grpwt = 0.0 # Loop over the pixels in the subraster. ypixel = subim + (ifymin - iymin) * ncols + ifxmin - ixmin - 1 do ypix = ifymin, ifymax { ytemp = ypix pixel = ypixel do xpix = ifxmin, ifxmax { xtemp = xpix pixel = pixel + 1 datum = Memr[pixel] # Reject pixel if not in good data range. if ((datum < mingdata) || (datum > maxgdata)) next # If this pixel is within one fitting radius of at # least one star then include it in the solution. # While figuring this out, compute the squared distance # of this pixel from the centroid of each star in the # group, and sum its contribution to the number # of pixels within one fitting radius of each star. if (dp_ntomit (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)], Memi[DP_NSKIP(nstar)], group_size, xtemp, ytemp, cutoff)) next if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == NO)) { sky_value = dp_ntskyval (Memr[DP_APMSKY(apsel)], Memi[DP_NSKIP(nstar)], group_size) if (IS_INDEFR(sky_value)) sky_value = mean_sky } else sky_value = mean_sky # Subtract the mean sky from the pixel. #ds = datum - mean_sky ds = datum - sky_value # Now loop over the stars and subtract from this pixel # the light contribution from each star within one psf # radius. # # The condition equation for pixel[xpix,ypix] is # # residual = data[xpix,ypix] - sky - sum (scale * # psf[xpix-xc,ypix-yc]) # # The scale, xc, and yc are iterated until # # sum (weight * residual ** 2) # # is minimized. Weight will be a function of 1) the # distance of the pixel from the center of the nearest # star, 2) the model-predicted brightness of the pixel # (taking into consideration the readout noise, the # photons/ADU, and the interpolating error of the PSF, # and, 3) the size of the residual itself. One is # necessary to prevent the non-linear least squares # solution from oscillating: oftimes it will come # to pass that if you include a pixel in the solution # then the predicted shift of the centroid will cause # that pixel to be excluded in the next iteration, and the # new predicted shift of the centroid will cause that # pixel to be included again. This could go on ad # infinitum. The cure is to have the weight of the pixel # go asymptotically to zero as its distance from the # stellar centroid approaches the fitting radius. In # the case like that just described, the solution can # then find a real minimum of the sum of the weighted # squared residuals with that pixel at some low-weight # position just inside the fitting radius. Two is # just sensible weighting. Three is a crude attempt # at making the solution more robust against bad pixels. weight = dp_ntsubtract (dao, im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_NRPIXSQ(nstar)], Memr[DP_APMAG(apsel)], Memi[DP_NSKIP(nstar)], Memr[DP_NX(nstar)], group_size, xtemp, ytemp, ds, psfradsq, fitradsq, DP_RECENTER(dao)) # At this point the vector X contains the first # derivative of the condition equation, for the pixel # under consideration, with respect to each of the # fitting parameters for all of these stars. We now need # to add these values into the normal matrix and the vector # of residuals. # 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 # of the total brightness for the difficulty of flat- # fielding and bias subtraction, plus an estimated error # of the same fraction of the fourth derivative at the # peak of the profile, to account for the difficulty # of accurately interpolating within the point-spread # function. The fourth derivative of the PSF is # proportional to H / sigma ** 4 (sigma is the Gaussian # width parameter for the stellar core); using the geometric # mean of sigma(x) and sigma(y), this becomes H / [sigma(x) # * sigma(y)] ** 2. The ratio of the fitting error to this # quantity is estimated to be approximately 0.027 from a # good-seeing CTIO frame. (This is probably a function of # seeing, sampling etc.) # Get the residual from the PSF fit and the pixel # intensity as predicted by the fit. Pred_pixval = raw data # minus residual = model predicted value of the intensity at # this point. pred_pixval = max (0.0, datum - ds) if ((pred_pixval > maxgdata) && (iter >= 4)) next #sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + #(perr * pred_pixval) ** 2 + (peakerr * #(pred_pixval - mean_sky)) ** 2 sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + (perr * pred_pixval) ** 2 + (peakerr * (pred_pixval - sky_value)) ** 2 if (sigmasq <= 0.0) next relerr = abs (ds) / sqrt (sigmasq) # Add this residual into the weighted sum of the # absolute relative residuals. sumres = sumres + relerr * weight grpwt = grpwt + weight # Add into the accumulating sums of the weighted # absolute relative residuals and of the image sharpness # parameter for each of the stars. Include in the # sharpness index only those pixels within NCORE_SIGMA # sigma of the centroid of the object. This saves time # and floating underflows by excluding pixels # which contribute less than about one part in a # million to the index. call dp_acsharp (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memi[DP_NSKIP(nstar)], Memi[DP_NNPIX(nstar)], Memr[DP_NNUMER(nstar)], Memr[DP_NDENOM(nstar)], Memr[DP_NSUMWT(nstar)], Memr[DP_APCHI(apsel)], group_size, xtemp, ytemp, Memr[DP_PSFPARS(psffit)], Memr[DP_PSFPARS(psffit)+1], ds, sigmasq, relerr, weight) # If clipping is in effect, reduce the weight of a bad # pixel. A pixel having a residual of 2.5 sigma gets # reduced to half weight and one with a rersidual of 5 # sigma gets weight of 1 / 257. weight = weight / sigmasq if (clip) weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) * chigrp)) ** DP_CLIPEXP(dao)) dswt = ds * weight # Work this pixel into the normal matrix. call dp_mataccum (Memr[DP_NX(nstar)], Memi[DP_NSKIP(nstar)], group_size, Memr[DP_NC(nstar)], Memr[DP_NV(nstar)], cdimen, nterm, weight, dswt, DP_RECENTER(dao), DP_FITSKY(dao)) } ypixel = ypixel + ncols } # Make sure that every star in the group has at least MIN_FIT_PIXEL # pixels within one fitting radius. refit = dp_ntmin (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memi[DP_NNPIX(nstar)], Memi[DP_NIER(nstar)], group_size, nterm, DP_RECENTER(dao), DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, DP_VERBOSE(dao)) if (group_size < 1) return (group_size) if (refit) next # Reflect the normal matrix across the diagonal. call dp_mreflect (Memr[DP_NC(nstar)], cdimen, nterm) # Compute the robust estimate of the standard deviation of the # residuals for the group as a whole, and for each star. This # estimate is sqrt (PI/2) * weighted mean absolute relative # residual. Do you like that "absolute relative" stuff? (PBS) # NO! (LED) # # CHI = CHI_NORM * SUM (weight * resid) / (# of pixels) # # This gets corrected for bias by being multiplied by # # SQRT (# of pixels) / (# of pixels - 3) # Then the value is driven towards unity, depending on # exactly how many pixels were involved: if CHIOLD is based # on a total weight of 3, then it is extremely poorly determined # and we just want to keep CHIOLD = 1. The larger GRPWT is, the # better determined CHIOLD is, and the less we want to force it # toward unity. So, just take the weighted average of CHIOLD and # unity, with weights GRPWT - 3 and 1, respectively. if (grpwt > ntmin) { chigrp = CHI_NORM * sumres * sqrt (1.0 / (grpwt * (grpwt - ntmin))) chigrp = ((grpwt - ntmin) * chigrp + ntmin) / grpwt } # CHIOLD has been pulled toward 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 # greated than one, CHI will be nudged toward the group average. # In order to work optimally, of course, this requires that # the # of photons / ADU, the READNOISE and the other noise # contributors are properly specified. call dp_ntchi (Memr[DP_NSUMWT(nstar)], Memr[DP_APCHI(apsel)], group_size, ntmin, chigrp, grpwt) # Invert the matrix. call invers (Memr[DP_NC(nstar)], cdimen, nterm, flag) # Check for a singular matrix. refit = dp_ncheckc (Memr[DP_NC(nstar)], cdimen, nterm, Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size, DP_RECENTER(dao), DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, DP_VERBOSE(dao)) if (group_size < 1) return (group_size) if (refit) next # Solve for position and scale factor increments. call mvmul (Memr[DP_NC(nstar)], cdimen, nterm, Memr[DP_NV(nstar)], Memr[DP_NX(nstar)]) if (iter <= 1) refit = true else refit = false # Fit the sky. if (DP_FITSKY(dao) == YES) { noise = sqrt (abs (mean_sky / DP_PHOTADU(dao)) + read_noise) mean_sky = mean_sky - max (-3.0 * noise, min (Memr[DP_NX(nstar)+nterm-1], 3.0 * noise)) if (abs (Memr[DP_NX(nstar)+nterm-1]) > (1.0e-4 * mean_sky)) refit = true } # 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 pixels per iteration. Any time that the # parameter correction changes sign from one iteration to the # next, the maximum permissible change will be reduced by a factor # of two. These clamps are released any time a star disappears. call dp_ntclamp (Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)], Memr[DP_NSUMWT(nstar)], group_size, Memr[DP_NC(nstar)], cdimen, Memr[DP_NX(nstar)], Memr[DP_NXOLD(nstar)], Memr[DP_NXCLAMP(nstar)], DP_RECENTER(dao), clip, refit) # Check whether the estimated centroid of the any star has # moved so far out of the limits of the picture that it has fewer # than 4 or 5 pixels within one fitting radius. call dp_ntcentroid (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], group_size, ixmin, ixmax, iymin, iymax, fitradsq, DP_FITSKY(dao), DP_GROUPSKY(dao), mean_sky, refit, DP_VERBOSE(dao)) if (group_size < 1) return (group_size) # Update matrix dimensions. if (DP_RECENTER(dao) == YES) nterm = 3 * group_size else nterm = group_size if (DP_FITSKY(dao) == YES) nterm = nterm + 1 # Now check whether any of the stars is too faint (more than 12.5 # magnitudes fainter than the PSF star). If several stars are too # faint, delete the faintest one, and set the brightness of the # other faint ones to 12.5 magnitudes below the PSF star. That way # on the next iteration we will see whether these stars want to # grow or to disappear. faint = 0.0 ifaint = 0 call dp_ntfmag (Memr[DP_APMAG(apsel)], group_size, tfaint, tifaint) # If at least one star is more than 12.5 magnitudes fainter # than the PSF then ifaint is the relative index of the faintest # of them, and faint is the relative brightness of the # faintest of them. # No very faint star was detected. if (tifaint <= 0) { # If the solution has not converged and if the number of # iterations is less than MIN_ITER, perform another iteration # with no questions asked. if ((refit) && (iter < MIN_ITER)) return (group_size) # If the solution doesn't think it has converged, after the # fourth iteration delete the least certain star if it is less # less than a one-sigma detection; after the eighth iteration # delete the least certain star if it is less than a 1.5 sigma # detection; after the twelfth iteration OR if the solution # thinks it has converged, delete the least certain star if it # is less than a two-sigma detection. call dp_fsnoise (Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)], group_size, faint, ifaint) if ((refit) && (iter < DP_MAXITER (dao)) && (faint < wcrit)) return (group_size) } # Reject the appropriate star. if ((tifaint > 0) || (faint >= MIN_FAINT)) { # Either (a) the solution thinks it has not converged # and the faintest star is more uncertain than sqrt(wcrit) # or (b) the solution thinks it has converged and the # faintest star is more uncertain than two-sigma. if (DP_VERBOSE (dao) == YES) { mindex = max (tifaint, ifaint) call printf ( "\tStar %-5d has been deleted because it is too faint\n") call pargi (Memi[DP_APID(apsel)+mindex-1]) } if ((DP_FITSKY(dao) == NO) && (DP_GROUPSKY(dao) == YES) && (group_size > 1)) mean_sky = (mean_sky * group_size - Memr[DP_APMSKY(apsel)+max(tifaint,ifaint)-1]) / (group_size - 1) call dp_remove (max (tifaint, ifaint), group_size, Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memi[DP_NIER(nstar)], NSTERR_FAINT) if (group_size < 1) return (group_size) if (DP_RECENTER(dao) == YES) nterm = 3 * group_size else nterm = group_size if (DP_FITSKY(dao) == YES) nterm = nterm + 1 call aclrr (Memr[DP_NXOLD(nstar)], nterm) call amovkr (1.0, Memr[DP_NXCLAMP(nstar)], nterm) clip = false iter = max (1, iter - 1) next } # Solution has either converged or gone to the maximum number # of iterations. if ((iter < DP_MAXITER(dao)) && (! clip)) { # The first convergence milestone has been reached. Turn on the # clipper, loosen the clamps and keep on going. if (DP_CLIPEXP(dao) > 0) clip = true converge = false call aclrr (Memr[DP_NXOLD(nstar)], nterm) call amaxkr (Memr[DP_NXCLAMP(nstar)], 0.25, Memr[DP_NXCLAMP(nstar)], nterm) return (group_size) } converge = true } until (converge) return (group_size) end # DP_NSTMERGE -- Decide whether two stars in a group should merge. bool procedure dp_nstmerge (xcen, ycen, mag, magerr, group_size, sepcrit, sepmin, wcrit, i, j, k) real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real mag[ARB] # array of magnitudes real magerr[ARB] # array of magnitude errors int group_size # group size real sepcrit # the critical separation squared real sepmin # the minimum separation real wcrit # critical error for rejection int i, j, k # output indices real separation begin do i = 1, group_size { do j = 1, i - 1 { # Compute the separation. separation = (xcen[j] - xcen[i]) ** 2 + (ycen[j] - ycen[i]) ** 2 if (separation > sepcrit) next # Find the fainter of the two stars. k = j if (mag[i] < mag[j]) k = i if ((separation < sepmin) || ((magerr[k] / mag[k]) > wcrit)) return (true) } } return (false) end # DP_NSTCEN -- Recompute the centroid and brightness of the i-th star. procedure dp_nstcen (xcen, ycen, mag, i, j, k) real xcen[ARB] # the x centers real ycen[ARB] # the y centers real mag[ARB] # the magnitudes int i, j, k # array indices begin # Now eliminate the fainter of the two stars. The k-th # star is now the fainter of the two, the i-th the # brighter. if (mag[i] < mag[j]) i = j # Now replace the centroid of the i-th star with the # weighted mean of the most recent estimates of the # centroids of the i-th and the k-th stars, and the # brightness of i-th with the sum of the brightnesses. xcen[i] = xcen[i] * mag[i] + xcen[k] * mag[k] ycen[i] = ycen[i] * mag[i] + ycen[k] * mag[k] mag[i] = mag[i] + mag[k] xcen[i] = xcen[i] / mag[i] ycen[i] = ycen[i] / mag[i] end # DP_NTOMIT -- Check whether a pixel is within one fitting radius of another # star. bool procedure dp_ntomit (xcen, ycen, rpixsq, skip, group_size, fx, fy, cutoff) real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real rpixsq[ARB] # array of distances squared int skip[ARB] # array of skip values int group_size # the group size real fx, fy # pixel position in image real cutoff # radius cuttoff for pixel inclusion bool omit int i begin omit = true do i = 1, group_size { skip[i] = YES rpixsq[i] = (fx - xcen[i]) ** 2 + (fy - ycen[i]) ** 2 if (rpixsq[i] > cutoff) next skip[i] = NO omit = false } return (omit) end # DP_NTSKYVAL -- Compute the average sky value to use for a particular # pixel by averaging the sky values of all stars for which the # pixel is within one fitting radius. real procedure dp_ntskyval (sky, skip, nstar) real sky[ARB] # array of sky values int skip[ARB] # array of skip values int nstar # the number of stars int i, npts real sum begin sum = 0.0 npts = 0 do i = 1, nstar { if (skip[i] == YES) next sum = sum + sky[i] npts = npts + 1 } if (npts <= 0) return (INDEFR) else return (sum / npts) end # DP_NTSUBTRACT -- Procedure to subtract the contribution of a particular # pixel from a particular star. real procedure dp_ntsubtract (dao, im, xcen, ycen, rpixsq, mag, skip, x, group_size, fx, fy, ds, psfradsq, fitradsq, recenter) pointer dao # pointer to the daophot structure pointer im # the input image descriptor real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real rpixsq[ARB] # array of distances squared real mag[ARB] # array of magnitudes int skip[ARB] # array of skip values real x[ARB] # x accumulation vector array int group_size # size of the group real fx, fy # center of pixel in image real ds # pixel value real psfradsq # psf radius squared real fitradsq # fit radius squared int recenter # recenter the coordinates int i, i3, k pointer psffit real weight, dx, dy, deltax, deltay, val, dvdx, dvdy, rsq real dp_usepsf() begin psffit = DP_PSFFIT(dao) weight = 0.0 do i = 1, group_size { if (rpixsq[i] >= psfradsq) next dx = fx - xcen[i] dy = fy - 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 val = 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) ds = ds - mag[i] * val if (skip[i] == YES) next rsq = rpixsq[i] / fitradsq weight = max (weight, 5.0 / (5.0 + rsq / (1.0 - rsq))) if (recenter == YES) { i3 = 3 * i k = i3 - 2 x[k] = -val k = i3 - 1 x[k] = -mag[i] * dvdx x[i3] = -mag[i] * dvdy } else x[i] = -val } return (weight) end # DP_ACSHARP -- Procedure to accumulate sums of the weighted absolute # relative residuals and the image sharpness parameter for each of the # stars. procedure dp_acsharp (xcen, ycen, skip, npix, numer, denom, sumwt, chi, group_size, fx, fy, fwhmx, fwhmy, ds, sigmasq, relerr, weight) real xcen[ARB] # array of object x centers real ycen[ARB] # array of object y centers int skip[ARB] # array of skip values int npix[ARB] # array of numbers of pixels real numer[ARB] # numerator array real denom[ARB] # denominator array real sumwt[ARB] # array of summed weights real chi[ARB] # array of chis int group_size # group size paramter. real fx, fy # position of data in image real fwhmx, fwhmy # gaussian core widths in x and y real ds # the data residual real sigmasq # the sigma squared real relerr # the relative error real weight # the weight int i real rhosq, dfdsig begin do i = 1, group_size { if (skip[i] == YES) next # Accumulate the number of pixels, chi and sum of the weights. npix[i] = npix[i] + 1 chi[i] = chi[i] + relerr * weight sumwt[i] = sumwt[i] + weight # Include in the sharpness index only those pixels # within NCORE_SIGMASQ of the centroid of the # object. (This saves time and floating underflows # by excluding pixels which contribute very little # to the index. rhosq = ((xcen[i] - fx) / fwhmx) ** 2 + ((ycen[i] - fy) / fwhmy) ** 2 if (rhosq > NCORE_SIGMASQ) next rhosq = 0.6931472 * rhosq dfdsig = exp (-rhosq) * (rhosq - 1.0) numer[i] = numer[i] + dfdsig * ds / sigmasq denom[i] = denom[i] + (dfdsig ** 2) / sigmasq } end # DP_MATACCUM -- Procedure to accumulate the data into the matrices. procedure dp_mataccum (x, skip, group_size, c, v, cdimen, nterm, weight, dswt, recenter, fitsky) real x[ARB] # x array int skip[ARB] # skip vector int group_size # size of the group real c[cdimen,ARB] # coefficient matrix real v[ARB] # vector array int cdimen # dimensions of the coefficient matrix int nterm # the number of terms real weight # weight real dswt # data weight int recenter # recenter the coordinates int fitsky # fit the sky value int i, i3, i3m2, k, j, l begin if (fitsky == YES) { c[nterm,nterm] = c[nterm,nterm] + weight v[nterm] = v[nterm] - dswt } do i = 1, group_size { if (skip[i] == YES) next if (recenter == YES) { i3 = i * 3 i3m2 = i3 - 2 do k = i3m2, i3 { if (fitsky == YES) c[nterm,k] = c[nterm,k] - x[k] * weight v[k] = v[k] + x[k] * dswt } do j = 1, i { if (skip[j] == YES) next do k = i3m2, i3 { do l = 3 * j - 2, min (k, 3 * j) c[k,l] = c[k,l] + x[k] * x[l] * weight } } } else { v[i] = v[i] + x[i] * dswt if (fitsky == YES) c[nterm,i] = c[nterm,i] - x[i] * weight do j = 1, i { if (skip[j] == YES) next c[i,j] = c[i,j] + x[i] * x[j] * weight } } } end # DP_NTMIN -- Make sure that every star in the group has at least # MIN_NPIX pixels within one fitting radius. bool procedure dp_ntmin (ids, xcen, ycen, mag, sky, npix, nier, group_size, nterm, recenter, fitsky, groupsky, mean_sky, verbose) int ids[ARB] # array of ids real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real mag[ARB] # array of magnitudes real sky[ARB] # array of sky values int npix[ARB] # array of pixel numbers int nier[ARB] # array of error codes int group_size # size of the group int nterm # number of terms int recenter # recenter the objects int fitsky # fit the sky value int groupsky # use group sky value real mean_sky # the current mean sky value int verbose # verbose flag bool redo int i begin redo = false do i = 1, group_size { if (npix[i] >= MIN_NPIX) next redo = true if (verbose == YES) { call printf ( "\tStar %-5d has been deleted: too few good pixels\n") call pargi (ids[i]) } if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1) call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, nier, NSTERR_NOPIX) if (group_size <= 0) return (redo) if (recenter == YES) nterm = 3 * group_size else nterm = group_size if (fitsky == YES) nterm = nterm + 1 } return (redo) end # DP_MREFLECT -- Reflect the normal matrix around the diagonal. procedure dp_mreflect (c, cdimen, nterm) real c[cdimen,ARB] # coefficient matrix int cdimen # dimension of the c matrix int nterm # number of terms int l, k begin # Reflect the normal matrix across the diagonal. do l = 2, nterm { do k = 1, l - 1 c[k,l] = c[l,k] } end # DP_NTCHI -- Compute the chi value for each star. procedure dp_ntchi (sumwt, chi, group_size, ntmin, chigrp, grpwt) real sumwt[ARB] # sum of the weights real chi[ARB] # the chis:wq int group_size # size of the group int ntmin # minimum number of points real chigrp # the group chi real grpwt # the group weight int i begin do i = 1, group_size { if (sumwt[i] > ntmin) { chi[i] = CHI_NORM * chi[i] / sqrt ((sumwt[i] - ntmin) * sumwt[i]) sumwt[i] = ((sumwt[i] - ntmin) * chi[i] + MIN_SUMWT) / sumwt[i] } else { chi[i] = chigrp sumwt[i] = grpwt } } end # DP_NCHECKC -- Check the inverted matrix for singularity. bool procedure dp_ncheckc (c, cdimen, nterm, ids, xcen, ycen, mag, sky, nier, group_size, recenter, fitsky, groupsky, mean_sky, verbose) real c[cdimen,ARB] # coefficient matrix int cdimen # dimension of the c matrix int nterm # number of terms int ids[ARB] # array of ids real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real mag[ARB] # array of magnitudes real sky[ARB] # array of sky values int nier[ARB] # array of error codes int group_size # size of the group int recenter # recenter the objects int fitsky # fit the sky value int groupsky # use group sky value real mean_sky # the current mean sky value int verbose # verbose flag bool redo int j, starno, i begin redo = false do j = 1, nterm { if (c[j,j] > 0.0) next redo = true if ((j == nterm) && (fitsky == YES)) starno = 0 else if (recenter == YES) starno = (j + 2) / 3 else starno = j if (starno == 0) { if (verbose == YES) { do i = 1, group_size { call printf ( "\tStar %-5d has been deleted: singular matrix\n") call pargi (ids[i]) } } call amovki (NSTERR_SINGULAR, nier, group_size) group_size = 0 } else { if (verbose == YES) { call printf ( "\tStar %-5d has been deleted: singular matrix\n") call pargi (ids[starno]) } if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) mean_sky = (mean_sky * group_size - sky[starno]) / (group_size - 1) call dp_remove (starno, group_size, ids, xcen, ycen, mag, sky, nier, NSTERR_SINGULAR) } if (group_size <= 0) return (redo) if (recenter == YES) nterm = 3 * group_size else nterm = group_size if (fitsky == YES) nterm = nterm + 1 } return (redo) end # DP_NTCLAMP -- Restrict the amount the solution can vary on each iteration. procedure dp_ntclamp (xcen, ycen, mag, magerr, sumwt, group_size, c, cdimen, x, xold, clamp, recenter, clip, redo) real xcen[ARB] # x centers array real ycen[ARB] # y centers array real mag[ARB] # magnitude array real magerr[ARB] # magnitude errors array real sumwt[ARB] # array of weight sums int group_size # size of the group real c[cdimen, ARB] # coefficient matrix int cdimen # dimensions of c real x[ARB] # x vector real xold[ARB] # old x vector real clamp[ARB] # clamp on solution matrix int recenter # recenter the objects bool clip # clip the matrix bool redo # redo the solution int i, l, j, k real df begin do i = 1, group_size { # If any correction has changed sign since the last # iteration, reduce the maximum permissible change by # a factor of two. # 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 (the least-squares thinks that the star should # be brighter) a change of 1 magnitude is a change of a # factor of 2.5; if the brightness correction is positive # (the star should be fainter) a change of 1 magnitude # is a change of 60%. if (recenter == YES) { l = 3 * i k = l - 1 j = l - 2 if ((xold[j] * x[j]) < 0.0) clamp[j] = 0.5 * clamp[j] mag[i] = mag[i] - x[j] / (1.0 + max (x[j] / (MAX_DELTA_FAINTER * mag[i]), -x[j] / (MAX_DELTA_BRIGHTER * mag[i])) / clamp[j]) xold[j] = x[j] if ((xold[k] * x[k]) < 0.0) clamp[k] = 0.5 * clamp[k] if ((xold[l] * x[l]) < 0.0) clamp[l] = 0.5 * clamp[l] xcen[i] = xcen[i] - x[k] / (1.0 + abs(x[k]) / (clamp[k] * MAX_DELTA_PIX)) ycen[i] = ycen[i] - x[l] / (1.0 + abs(x[l]) / (clamp[l] * MAX_DELTA_PIX)) xold[k] = x[k] xold[l] = x[l] magerr[i] =sumwt[i] * sqrt (c[j,j]) } else { if ((xold[i] * x[i]) < 0.0) clamp[i] = 0.5 * clamp[i] mag[i] = mag[i] - x[i] / (1.0 + max (x[i] / (MAX_DELTA_FAINTER * mag[i]), -x[i] / (MAX_DELTA_BRIGHTER * mag[i])) / clamp[i]) xold[i] = x[i] magerr[i] =sumwt[i] * sqrt (c[i,i]) } # There are two milestones in the convergence process: the fits # proceed normally until each star's magnitude changes by less # than its standard error or MAX_NEW_ERRMAG magnitudes, whichever # is greater, and its x and y centroids change by less than 0.02 # pixel. At this point the least squares begins to apply # down-weighting of pixels with large residuals as described # above. The fits then continue until each star's # magnitude changes by less than MAX (MAX_NEW_ERRMAG * std. error, # MAX_NEW_RELBRIGHT2 magnitude), ad its centroids change by # less than 0.002 pixel. if (redo) next if (clip) { if (abs (x[j]) > max (MAX_NEW_ERRMAG * magerr[i], MAX_NEW_RELBRIGHT2 * mag[i])) { redo = true } else if (recenter == YES) { 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 { if (abs (x[j]) > max (magerr[i], MAX_NEW_RELBRIGHT1 * mag[i])) { redo = true } else if (recenter == YES) { df = sumwt[i] ** 2 if (x[k] ** 2 > max (df * c[k,k], MAX_PIXERR1)) redo = true else if (x[l] ** 2 > max (df * c[l,l], MAX_PIXERR1)) redo = true } } } end # DP_NTCENTROID -- Check the new centroids to see if they have moved too # far off the edge of the image. procedure dp_ntcentroid (ids, xcen, ycen, mag, sky, nier, group_size, ixmin, ixmax, iymin, iymax, fitradsq, fitsky, groupsky, mean_sky, redo, verbose) int ids[ARB] # array of ids real xcen[ARB] # array of x centers real ycen[ARB] # array of y centers real mag[ARB] # array of magnitudes real sky[ARB] # array of sky values int nier[ARB] # array of error codes int group_size # size of the group int ixmin,ixmax # subraster x limits int iymin,iymax # subraster y limits real fitradsq # fit radius squared int fitsky # fit the sky value int groupsky # use the group sky value real mean_sky # the mean sky value bool redo # redo fit int verbose # verbose mode int i real dx, dy begin # Check whether the centroid of any star has moved so far outside # the picture that it has fewer than four or five pixels within # one fitting radius. do i = 1, group_size { # If the centroid of the star is outside the picture in x or # y, then DX or DY is its distance from the center of the edge # pixel; otherwise DX and DY are zero. dx = max (ixmin - xcen[i], xcen[i] - ixmax, 0.0) dy = max (iymin - ycen[i], ycen[i] - iymax, 0.0) if ((dx <= MAX_PIX_INCREMENT) && (dy <= MAX_PIX_INCREMENT)) next if (((dx + 1.0) ** 2 + (dy + 1.0) ** 2) < fitradsq) next # Print a warning message about the star. if (verbose == YES) { call printf ( "\tStar %-5d has been deleted: new center too far off image\n") call pargi (i) } # Adjust the sky. if ((fitsky == NO) && (groupsky == YES) && (group_size > 1)) mean_sky = (mean_sky * group_size - sky[i]) / (group_size - 1) # Delete it. call dp_remove (i, group_size, ids, xcen, ycen, mag, sky, nier, NSTERR_OFFIMAGE) if (group_size < 1) break redo = true } end # DP_NTFMAG -- Check for faint stars. procedure dp_ntfmag (mag, group_size, faint, ifaint) real mag[ARB] # array of magnitudes int group_size # size of the group real faint # faintest magnitude int ifaint # index of faintest magnitude int i begin faint = 1.0 ifaint = 0 do i = 1, group_size { if (mag[i] > MIN_REL_BRIGHT) next if (mag[i] <= faint) { faint = mag[i] ifaint = i } mag[i] = MIN_REL_BRIGHT } end # DP_FSNOISE -- Compute the smallest signal to noise ratio. procedure dp_fsnoise (mag, magerr, group_size, faint, ifaint) real mag[ARB] # array of magnitudes real magerr[ARB] # array of magnitude errors int group_size # size of group real faint # faint value int ifaint # faint index int i real weight begin faint = 0.0 ifaint = 0 do i = 1, group_size { weight = magerr[i] / mag[i] if (weight < faint) next faint = weight ifaint = i } end # DP_REMOVE -- Remove the i-th star from the list of stars in the current # group by moving it to the end of the group. procedure dp_remove (i, nstar, ids, xcen, ycen, mag, sky, nier, pier) int i # the star to be removed int nstar # the number of stars in the group int ids[ARB] # the array of star ids real xcen[ARB] # the array of star x positions real ycen[ARB] # the array of star y positions real mag[ARB] # the array of magnitudes real sky[ARB] # the array of sky values. int nier[ARB] # array of error codes int pier # error code for deleted star int ihold, phold real xhold, yhold, shold, mhold begin nier[i] = pier if (i != nstar) { ihold = ids[i] xhold = xcen[i] yhold = ycen[i] shold = sky[i] mhold = mag[i] phold = nier[i] ids[i] = ids[nstar] xcen[i] = xcen[nstar] ycen[i] = ycen[nstar] sky[i] = sky[nstar] mag[i] = mag[nstar] nier[i] = nier[nstar] ids[nstar] = ihold xcen[nstar] = xhold ycen[nstar] = yhold sky[nstar] = shold mag[nstar] = mhold nier[nstar] = phold } nstar = nstar - 1 end