diff options
Diffstat (limited to 'noao/digiphot/apphot/fitpsf/apsfradgauss.x')
| -rw-r--r-- | noao/digiphot/apphot/fitpsf/apsfradgauss.x | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitpsf/apsfradgauss.x b/noao/digiphot/apphot/fitpsf/apsfradgauss.x new file mode 100644 index 00000000..0b4a93c6 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfradgauss.x @@ -0,0 +1,183 @@ +include <math/nlfit.h> +include "../lib/noise.h" +include "../lib/fitpsf.h" + +define NPARAMETERS 5 +define TOL 0.001 + +# APSFRADGAUSS -- Fit a radial Gaussian function to the data. + +int procedure apsfradgauss (ctrpix, nx, ny, emission, fwhmpsf, datamin, + datamax, noise, gain, sigma, maxiter, k2, nreject, par, perr, npar) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +int emission # emission or absorption version +real fwhmpsf # full width half max of the psf +real datamin # minimum good data value +real datamax # maximum good data value +int noise # noise model to be used +real gain # the gain in the data +real sigma # sigma of constant noise term +int maxiter # maximum number of iterations +real k2 # k-sigma rejection criterion +int nreject # maximum number of rejection cycles +real par[ARB] # parameters +real perr[ARB] # errors in parameters +int npar # number of parameters + +extern gaussr, dgaussr +int i, j, npts, list, imin, imax, fier +pointer sp, x, w, zfit, nl, ptr +real sumw, dummy, chisqr, locut, hicut +int locpr(), apreject() +real asumr(), apwssqr() + +begin + # Initialize. + npts = nx * ny + if (npts < NPARAMETERS) + return (AP_NPSF_TOO_SMALL) + + call smark (sp) + call salloc (x, 2 * npts, TY_REAL) + call salloc (w, npts, TY_REAL) + call salloc (zfit, npts, TY_REAL) + call salloc (list, NPARAMETERS, TY_INT) + + # Define the active parameters. + do i = 1, NPARAMETERS + Memi[list+i-1] = i + + # Set variables array. + ptr = x + do j = 1, ny { + do i = 1, nx { + Memr[ptr] = i + Memr[ptr+1] = j + ptr = ptr + 2 + } + } + + # Define the weight array. + switch (noise) { + case AP_NCONSTANT: + call amovkr (1.0, Memr[w], npts) + case AP_NPOISSON: + call amaxkr (ctrpix, 0.0, Memr[w], npts) + if (gain > 0.0) + call adivkr (Memr[w], gain, Memr[w], npts) + if (! IS_INDEFR(sigma)) + call aaddkr (Memr[w], sigma ** 2, Memr[w], npts) + call apreciprocal (Memr[w], Memr[w], npts, 1.0) + default: + call amovkr (1.0, Memr[w], npts) + } + + # Make an initial guess at the parameters. + if (emission == YES) + call ap_wlimr (ctrpix, Memr[w], npts, datamin, datamax, + par[5], par[1], imin, imax) + else + call ap_wlimr (ctrpix, Memr[w], npts, datamin, datamax, + par[1], par[5], imax, imin) + par[1] = par[1] - par[5] + if (mod (imax, nx) == 0) + imin = imax / nx + else + imin = imax / nx + 1 + par[3] = imin + imin = imax - (imin - 1) * nx + par[2] = imin + par[4] = (fwhmpsf ** 2 / 4.0) + + # Fit the function and the errors. + call nlinitr (nl, locpr (gaussr), locpr (dgaussr), par, perr, + NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, fier) + + # Perform the rejection cycle. + if ((nreject > 0) && (k2 > 0.0)) { + do i = 1, nreject { + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call asubr (ctrpix, Memr[zfit], Memr[zfit], npts) + chisqr = apwssqr (Memr[zfit], Memr[w], npts) + sumw = asumr (Memr[w], npts) + if (sumw <= 0.0) + break + else if (chisqr <= 0.0) + break + else + chisqr = sqrt (chisqr / sumw) + locut = - k2 * chisqr + hicut = k2 * chisqr + if (apreject (Memr[zfit], Memr[w], npts, locut, hicut) == 0) + break + call nlpgetr (nl, par, npar) + call nlfreer (nl) + call nlinitr (nl, locpr (gaussr), locpr (dgaussr), par, perr, + NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, + fier) + } + } + + # Get the parameters. + call nlpgetr (nl, par, npar) + par[4] = sqrt (abs(par[4])) + + # Get the errors. + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call nlerrorsr (nl, ctrpix, Memr[zfit], Memr[w], npts, dummy, + chisqr, perr) + perr[4] = sqrt (abs(perr[4])) + + # Compute the mean errors in the parameters. + dummy = 0.0 + do i = 1, npts { + if (Memr[w+i-1] > 0.0) + dummy = dummy + 1.0 + } + dummy = sqrt (dummy) + if (dummy > 0.0) + call adivkr (perr, dummy, perr, npar) + + call nlfreer (nl) + + call sfree (sp) + + # Return the appropriate error code. + if (fier == NO_DEG_FREEDOM) { + return (AP_NPSF_TOO_SMALL) + } else if (fier == SINGULAR) { + return (AP_PSF_SINGULAR) + } else if (fier == NOT_DONE) { + return (AP_PSF_NOCONVERGE) + } else { + return (AP_OK) + } +end + + +# APREJECT -- Reject points outside of the specified intensity limits by +# setting their weights to zero. + +int procedure apreject (pix, w, npts, locut, hicut) + +real pix[ARB] # data +real w[ARB] # weights +int npts # number of data points +real locut, hicut # data limits + +int i, nreject + +begin + nreject = 0 + do i = 1, npts { + if ((pix[i] < locut || pix[i] > hicut) && w[i] > 0.0) { + w[i] = 0.0 + nreject = nreject + 1 + } + } + return (nreject) +end |