diff options
Diffstat (limited to 'noao/digiphot/daophot/psf/dpfitpsf.x')
| -rw-r--r-- | noao/digiphot/daophot/psf/dpfitpsf.x | 1693 |
1 files changed, 1693 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/psf/dpfitpsf.x b/noao/digiphot/daophot/psf/dpfitpsf.x new file mode 100644 index 00000000..2b1c6bbc --- /dev/null +++ b/noao/digiphot/daophot/psf/dpfitpsf.x @@ -0,0 +1,1693 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" +include "../lib/peakdef.h" + +# DP_FITPSF -- Compute the PSF function. + +int procedure dp_fitpsf (dao, im, errmsg, maxch) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +char errmsg[ARB] # string containing error message +int maxch # max characters in errmsg + +int nfunc, func +pointer sp, flist, fstr, psf, psffit +int dp_pstati(), dp_unsatstar(), dp_fitana(), dp_ifitana() +int dp_fitlt(), dp_fctdecode(), dp_strwrd(), strdic() + +begin + # Get some daophot pointers. + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Test to see if there are any psf stars. + if (dp_pstati (dao, PNUM) <= 0) { + call sprintf (errmsg, maxch, "The PSF star list is empty") + return (ERR) + } + + # Test to see if there are any unsaturated stars. At the same time + # make sure that the first star is unsaturated. + if (dp_unsatstar (dao) <= 0) { + call sprintf (errmsg, maxch, + "There are no unsaturated PSF stars") + return (ERR) + } + + # Determine the analytic function. + call smark (sp) + call salloc (flist, SZ_FNAME, TY_CHAR) + call salloc (fstr, SZ_FNAME, TY_CHAR) + call dp_stats (dao, FUNCLIST, Memc[flist], SZ_FNAME) + nfunc = dp_fctdecode (Memc[flist], Memc[fstr], SZ_FNAME) + func = dp_strwrd (1, Memc[fstr], SZ_FNAME, Memc[flist]) + func = strdic (Memc[fstr], Memc[fstr], SZ_LINE, FCTN_FTYPES) + + # Compute the analytic part of the PSF function. + if (nfunc > 1 || func == FCTN_AUTO) { + + # Loop over all the analytic functions. + if (func == FCTN_AUTO) { + call strcpy (FCTN_FTYPES, Memc[flist], SZ_FNAME) + nfunc = FCTN_NFTYPES + } + + # Find the best fitting analytic function. + if (dp_ifitana (dao, im, Memc[flist], nfunc) == ERR) { + call sprintf (errmsg, maxch, + "Analytic function solution failed to converge") + call sfree (sp) + return (ERR) + } else if (DP_VERBOSE(dao) == YES) + call dp_listpars (dao) + + } else { + + # Save the analytic function for the new fit. + call strcpy (Memc[fstr], DP_FUNCTION(dao), SZ_LINE) + + # Initialize the parameters. + call dp_reinit (dao) + + # Fit the analytic part of the function. + if (dp_fitana (dao, im, Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], + Memr[DP_PH(psf)], Memi[DP_PSAT(psf)], DP_PNUM(psf)) == ERR) { + call sprintf (errmsg, maxch, + "Analytic function solution failed to converge") + call sfree (sp) + return (ERR) + } else if (DP_VERBOSE(dao) == YES) { + call printf ("\nFitting function %s norm scatter: %g\n") + call pargstr (DP_FUNCTION(dao)) + call pargr (DP_PSIGANA(psf)) + call dp_listpars (dao) + } + } + call sfree (sp) + + # Compute the look-up table. + if (dp_fitlt (dao, im) == ERR) { + call sprintf (errmsg, maxch, + "Too few stars to compute PSF lookup tables") + return (ERR) + } else if (DP_VERBOSE(dao) == YES) { + call printf ("\nComputed %d lookup table(s)\n") + call pargi (DP_NVLTABLE(psffit)+DP_NFEXTABLE(psffit)) + } + + return (OK) +end + + +# DP_UNSATSTAR -- Make sure there is at least one unsaturated star. + +int procedure dp_unsatstar (dao) + +pointer dao # pointer to the daophot structure + +int i, first_unsat, nstar +pointer psf + +begin + psf = DP_PSF(dao) + first_unsat = 0 + + nstar = 0 + do i = 1, DP_PNUM(psf) { + if (Memi[DP_PSAT(psf)+i-1] == YES) + next + nstar = nstar + 1 + if (first_unsat == 0) + first_unsat = i + } + + if (first_unsat > 1) + call dp_pfswap (dao, 1, first_unsat) + + return (nstar) +end + + +# DP_REINIT -- Reinitialize the psf function parameters. + +procedure dp_reinit (dao) + +pointer dao # pointer to the daophot structure + +pointer psffit +bool streq() + +begin + psffit = DP_PSFFIT(dao) + + # Define the psf function. + if (streq (DP_FUNCTION(dao), "gauss")) + DP_PSFUNCTION(psffit) = FCTN_GAUSS + else if (streq (DP_FUNCTION(dao), "moffat25")) + DP_PSFUNCTION(psffit) = FCTN_MOFFAT25 + else if (streq (DP_FUNCTION(dao), "moffat15")) + DP_PSFUNCTION(psffit) = FCTN_MOFFAT15 + else if (streq (DP_FUNCTION(dao), "penny1")) + DP_PSFUNCTION(psffit) = FCTN_PENNY1 + else if (streq (DP_FUNCTION(dao), "penny2")) + DP_PSFUNCTION(psffit) = FCTN_PENNY2 + else if (streq (DP_FUNCTION(dao), "lorentz")) + DP_PSFUNCTION(psffit) = FCTN_LORENTZ + else + call error (0, "Unknown PSF function\n") + + switch (DP_VARORDER(dao)) { + case -1: + DP_NVLTABLE(psffit) = 0 + case 0: + DP_NVLTABLE(psffit) = 1 + case 1: + DP_NVLTABLE(psffit) = 3 + case 2: + DP_NVLTABLE(psffit) = 6 + } + if (DP_FEXPAND(dao) == NO) + DP_NFEXTABLE(psffit) = 0 + else + DP_NFEXTABLE(psffit) = 5 + + # Set the initial values of the function parameters. + switch (DP_PSFUNCTION(psffit)) { + case FCTN_GAUSS: + DP_PSFNPARS(psffit) = 2 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + case FCTN_MOFFAT25: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + Memr[DP_PSFPARS(psffit)+3] = 2.5 + case FCTN_MOFFAT15: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + Memr[DP_PSFPARS(psffit)+3] = 1.5 + case FCTN_PENNY1: + DP_PSFNPARS(psffit) = 4 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.75 + Memr[DP_PSFPARS(psffit)+3] = 0.0 + case FCTN_PENNY2: + DP_PSFNPARS(psffit) = 5 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.75 + Memr[DP_PSFPARS(psffit)+3] = 0.0 + Memr[DP_PSFPARS(psffit)+4] = 0.0 + case FCTN_LORENTZ: + DP_PSFNPARS(psffit) = 3 + Memr[DP_PSFPARS(psffit)] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+1] = DP_FWHMPSF(dao) / 2.0 + Memr[DP_PSFPARS(psffit)+2] = 0.0 + default: + call error (0, "Unknown PSF function\n") + } +end + + +# DP_IFITANA -- Fit the PSF stars to each of the analytic functions in +# turn to determine which one gives the best fit. + +int procedure dp_ifitana (dao, im, funclist, nfunc) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +char funclist # the list of functions to be fit +int nfunc # number of functions + +int i, psftype, npars +pointer psf, psffit, sp, fstr, func +pointer ixtmp, iytmp, ihtmp, istmp, xtmp, ytmp, htmp, ptmp, stmp +real osig, osum, height, dhdxc, dhdyc, junk, ofactor, factor +int dp_strwrd(), strdic(), dp_fitana() +real dp_profile() + +begin + # Get some pointers. + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Allocate some temporary storage space. + call smark (sp) + call salloc (fstr, SZ_FNAME, TY_CHAR) + call salloc (func, SZ_FNAME, TY_CHAR) + + call salloc (ixtmp, DP_PNUM(psf), TY_REAL) + call salloc (iytmp, DP_PNUM(psf), TY_REAL) + call salloc (ihtmp, DP_PNUM(psf), TY_REAL) + call salloc (istmp, DP_PNUM(psf), TY_INT) + + call salloc (xtmp, DP_PNUM(psf), TY_REAL) + call salloc (ytmp, DP_PNUM(psf), TY_REAL) + call salloc (htmp, DP_PNUM(psf), TY_REAL) + call salloc (stmp, DP_PNUM(psf), TY_INT) + call salloc (ptmp, MAX_NFCTNPARS, TY_REAL) + + # Initialize. + call strcpy (DP_FUNCTION(dao), Memc[func], SZ_FNAME) + npars = 0 + osig = MAX_REAL + call amovr (Memr[DP_PXCEN(psf)], Memr[ixtmp], DP_PNUM(psf)) + call amovr (Memr[DP_PYCEN(psf)], Memr[iytmp], DP_PNUM(psf)) + call amovr (Memr[DP_PH(psf)], Memr[ihtmp], DP_PNUM(psf)) + call amovi (Memi[DP_PSAT(psf)], Memi[istmp], DP_PNUM(psf)) + ofactor = dp_profile (DP_PSFUNCTION(psffit), 0.0, 0.0, + Memr[DP_PSFPARS(psffit)], dhdxc, dhdyc, junk, 0) + + factor = 1 + do i = 1, nfunc { + + # Get the function name and set it. + if (dp_strwrd (i, Memc[fstr], SZ_FNAME, funclist) <= 0) + next + if (strdic (Memc[fstr], Memc[fstr], SZ_FNAME, FCTN_FTYPES) <= 0) + next + call strcpy (Memc[fstr], DP_FUNCTION(dao), SZ_FNAME) + + # Start from the same initial state. + call dp_reinit (dao) + call amovr (Memr[ixtmp], Memr[xtmp], DP_PNUM(psf)) + call amovr (Memr[iytmp], Memr[ytmp], DP_PNUM(psf)) + if (i == 1) + call amovr (Memr[ihtmp], Memr[htmp], DP_PNUM(psf)) + else { + factor = ofactor / dp_profile (DP_PSFUNCTION(psffit), 0.0, 0.0, + Memr[DP_PSFPARS(psffit)], dhdxc, dhdyc, junk, 0) + call amulkr (Memr[ihtmp], factor, Memr[htmp], DP_PNUM(psf)) + } + call amovi (Memi[istmp], Memi[stmp], DP_PNUM(psf)) + + call printf ("Trying function %s norm scatter = ") + call pargstr (Memc[fstr]) + + # Do the fit. + if (dp_fitana (dao, im, Memr[xtmp], Memr[ytmp], Memr[htmp], + Memi[stmp], DP_PNUM(psf)) == ERR) { + call printf ("error\n") + next + } else { + call printf ("%g\n") + call pargr (DP_PSIGANA(psf)) + } + + # Save the better fit. + if (DP_PSIGANA(psf) < osig) { + call strcpy (Memc[fstr], Memc[func], SZ_FNAME) + psftype = DP_PSFUNCTION(psffit) + height = DP_PSFHEIGHT(psffit) + npars = DP_PSFNPARS(psffit) + call amovr (Memr[DP_PSFPARS(psffit)], Memr[ptmp], + MAX_NFCTNPARS) + call amovr (Memr[xtmp], Memr[DP_PXCEN(psf)], DP_PNUM(psf)) + call amovr (Memr[ytmp], Memr[DP_PYCEN(psf)], DP_PNUM(psf)) + call amovr (Memr[htmp], Memr[DP_PH(psf)], DP_PNUM(psf)) + call amovi (Memi[stmp], Memi[DP_PSAT(psf)], DP_PNUM(psf)) + osig = DP_PSIGANA(psf) + osum = DP_PSUMANA(psf) + } + } + + # Restore the best fit parameters. + if (npars > 0) { + call strcpy (Memc[func], DP_FUNCTION(dao), SZ_FNAME) + DP_PSFUNCTION(psffit) = psftype + DP_PSFHEIGHT(psffit) = height + DP_PSFNPARS(psffit) = npars + DP_PSIGANA(psf) = osig + DP_PSUMANA(psf) = osum + call amovr (Memr[ptmp], Memr[DP_PSFPARS(psffit)], + MAX_NFCTNPARS) + call printf ("Best fitting function is %s\n") + call pargstr (DP_FUNCTION(dao)) + } + + # Cleanup. + call sfree (sp) + + if (npars > 0) + return (OK) + else + return (ERR) +end + + +# DP_FITANA -- Fit the analytic part of the psf function + +int procedure dp_fitana (dao, im, pxcen, pycen, ph, pstat, npsfstars) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +real pxcen[ARB] # x coordinates of the psf stars +real pycen[ARB] # y coordinates of the psf stars +real ph[ARB] # heights of the psf stars +int pstat[ARB] # saturation status of psf stars +int npsfstars # the number of psf stars + +int i, niter, istar, mpar, lx, ly, nx, ny, ier +pointer apsel, psf, psffit, data +real fitrad, rsq, oldchi, sumfree +pointer imgs2r() + +begin + # Get the psf fitting structure pointer. + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Define some variables. + oldchi = 0.0 + mpar = 2 + fitrad = DP_FITRAD(dao) + rsq = fitrad ** 2 + + # Get some memory. + call dp_amempsf (dao) + + # Initialize the fit. + call amovkr (0.5, Memr[DP_PCLAMP(psf)], DP_PSFNPARS(psffit)) + call aclrr (Memr[DP_PZ(psf)], DP_PSFNPARS(psffit)) + call aclrr (Memr[DP_POLD(psf)], DP_PSFNPARS(psffit)) + Memr[DP_PCLAMP(psf)] = 2.0 + Memr[DP_PCLAMP(psf)+1] = 2.0 + + call aclrr (Memr[DP_PXOLD(psf)], npsfstars) + call aclrr (Memr[DP_PYOLD(psf)], npsfstars) + call amovkr (1.0, Memr[DP_PXCLAMP(psf)], npsfstars) + call amovkr (1.0, Memr[DP_PYCLAMP(psf)], npsfstars) + + # Iterate. + do niter = 1, MAX_NPSFITER { + + # Initialize the current integration. + call aclrr (Memr[DP_PV(psf)], DP_PSFNPARS(psffit)) + call aclrr (Memr[DP_PC(psf)], DP_PSFNPARS(psffit) * + DP_PSFNPARS(psffit)) + + # Loop over the stars. + DP_PSIGANA(psf) = 0.0 + DP_PSUMANA(psf) = 0.0 + do istar = 1, npsfstars { + + # Test for saturation. + if (pstat[istar] == YES) + next + + # Define the subraster to be read in. + lx = int (pxcen[istar] - fitrad) + 1 + ly = int (pycen[istar] - fitrad) + 1 + nx = (int (pxcen[istar] + fitrad) - lx) + 1 + ny = (int (pycen[istar] + fitrad) - ly) + 1 + + # Is the star off the image? + if (lx > IM_LEN(im,1) || ly > IM_LEN(im,2) || (lx + nx - 1) < + 1 || (ly + ny - 1) < 1) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d is outside the image\n") + call pargi (Memi[DP_APID(apsel)+istar-1]) + } + next + } + + # Is the star too near the edge of the frame? + if (lx < 1 || ly < 1 || (lx + nx - 1) > IM_LEN(im,1) || + (ly + ny - 1) > IM_LEN(im,2)) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Star %d is too near the edge of the image\n") + call pargi (Memi[DP_APID(apsel)+istar-1]) + } + next + } + + # Read in the subraster. + data = imgs2r (im, lx, lx + nx - 1, ly, ly + ny - 1) + + # Fit x, y, and height for the PSF star istar. + call dp_xyhiter (DP_PSFUNCTION(psffit), + Memr[DP_PSFPARS(psffit)], rsq, Memr[data], nx, ny, lx, ly, + pxcen[istar], pycen[istar], + Memr[DP_APMSKY(apsel)+istar-1], ph[istar], + Memr[DP_PXCLAMP(psf)+istar-1], + Memr[DP_PYCLAMP(psf)+istar-1], Memr[DP_PXOLD(psf)+istar-1], + Memr[DP_PYOLD(psf)+istar-1]) + + # Fit the parameters for the entire list of stars + call dp_paccum (DP_PSFUNCTION(psffit), + Memr[DP_PSFPARS(psffit)], DP_PSFNPARS(psffit), mpar, rsq, + Memr[data], nx, ny, lx, ly, pxcen[istar], + pycen[istar], Memr[DP_APMSKY(apsel)+istar-1], + ph[istar], niter, Memr[DP_PC(psf)], + Memr[DP_PV(psf)], Memr[DP_PTMP(psf)], DP_PSIGANA(psf), + DP_PSUMANA(psf)) + } + + # Invert the matrix and compute the new parameters. + call invers (Memr[DP_PC(psf)], DP_PSFNPARS(psffit), mpar, ier) + call mvmul (Memr[DP_PC(psf)], DP_PSFNPARS(psffit), mpar, + Memr[DP_PV(psf)], Memr[DP_PZ(psf)]) + + do i = 1, mpar { + if ((Memr[DP_PZ(psf)+i-1] * Memr[DP_POLD(psf)+i-1]) < 0.0) + Memr[DP_PCLAMP(psf)+i-1] = 0.5 * + Memr[DP_PCLAMP(psf)+i-1] + else + Memr[DP_PCLAMP(psf)+i-1] = 1.1 * + Memr[DP_PCLAMP(psf)+i-1] + } + call amovr (Memr[DP_PZ(psf)], Memr[DP_POLD(psf)], mpar) + call amulr (Memr[DP_PZ(psf)], Memr[DP_PCLAMP(psf)], + Memr[DP_PZ(psf)], mpar) + + Memr[DP_PZ(psf)] = max (-0.1 * Memr[DP_PSFPARS(psffit)], + min (0.1 * Memr[DP_PSFPARS(psffit)], Memr[DP_PZ(psf)])) + Memr[DP_PZ(psf)+1] = max (-0.1 * Memr[DP_PSFPARS(psffit)+1], + min (0.1 * Memr[DP_PSFPARS(psffit)+1], Memr[DP_PZ(psf)+1])) + #if (mpar > 2) + #Memr[DP_PZ(psf)+2] = Memr[DP_PZ(psf)+2] / + #(1.0 + abs (Memr[DP_PZ(psf)+2]) / + #(min (0.1, 1.0 - abs (Memr[DP_PSFPARS(psffit)+2])))) + call aaddr (Memr[DP_PSFPARS(psffit)], Memr[DP_PZ(psf)], + Memr[DP_PSFPARS(psffit)], mpar) + + # Check for convergence. + sumfree = DP_PSUMANA(psf) - real (mpar + 3 * npsfstars) + if (sumfree > 0.0 && DP_PSIGANA(psf) >= 0.0) + DP_PSIGANA(psf) = sqrt (DP_PSIGANA(psf) / sumfree) + else + DP_PSIGANA(psf) = 9.999 + + if (mpar == DP_PSFNPARS(psffit)) { + if (abs (oldchi / DP_PSIGANA(psf) - 1.0) < 1.0e-5) { + DP_PSFHEIGHT(psffit) = ph[1] + if (IS_INDEFR(Memr[DP_PMAG(psf)])) + DP_PSFMAG(psffit) = Memr[DP_APMAG(apsel)] + else + DP_PSFMAG(psffit) = Memr[DP_PMAG(psf)] + DP_PSFX(psffit) = real (IM_LEN(im,1) - 1) / 2.0 + DP_PSFY(psffit) = real (IM_LEN(im,2) - 1) / 2.0 + return (OK) + } else + oldchi = DP_PSIGANA(psf) + } else { + if (abs (oldchi / DP_PSIGANA(psf) - 1.0) < 1.0e-3) { + mpar = mpar + 1 + oldchi = 0.0 + } else + oldchi = DP_PSIGANA(psf) + } + } + + return (ERR) +end + + +# DP_XYHITER -- Increment the initial x, y, and height values for a star. + +procedure dp_xyhiter (psftype, params, rsq, data, nx, ny, lx, ly, x, y, sky, h, + xclamp, yclamp, xold, yold) + +int psftype # analytic point spread function type +real params[ARB] # current function parameter values +real rsq # the fitting radius squared +real data[nx,ARB] # the input image data +int nx, ny # the dimensions of the input image data +int lx, ly # the coordinates of the ll corner of the image data +real x, y # the input/output stellar coordinates +real sky # the input sky value +real h # the input/output height value +real xclamp, yclamp # the input/output clamping factors for x and y +real xold, yold # the input/output x and y correction factors + +int i, j +real dhn, dhd, dxn, dxd, dyn, dyd, dx, dy, wt, dhdxc, dhdyc, junk, p, dp +real prod +real dp_profile() + +begin + dhn = 0.0 + dhd = 0.0 + dxn = 0.0 + dxd = 0.0 + dyn = 0.0 + dyd = 0.0 + + do j = 1, ny { + dy = real ((ly + j) - 1) - y + do i = 1, nx { + dx = real ((lx + i) - 1) - x + wt = (dx ** 2 + dy ** 2) / rsq + #if (wt >= 1.0) + if (wt >= 0.999998) + next + p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc, junk, 0) + dp = data[i,j] - h * p - sky + dhdxc = dhdxc * h + dhdyc = dhdyc * h + wt = 5.0 / (5.0 + (wt / (1.0 - wt))) + prod = wt * p + dhn = dhn + prod * dp + dhd = dhd + prod * p + prod = wt * dhdxc + dxn = dxn + prod * dp + dxd = dxd + prod * dhdxc + prod = wt * dhdyc + dyn = dyn + prod * dp + dyd = dyd + prod * dhdyc + } + } + + h = h + (dhn / dhd) + dxn = dxn / dxd + if ((xold * dxn) < 0.0) + xclamp = 0.5 * xclamp + xold = dxn + x = x + (dxn / (1.0 + (abs(dxn) / xclamp))) + dyn = dyn / dyd + if ((yold * dyn) < 0.0) + yclamp = 0.5 * yclamp + yold = dyn + y = y + (dyn / (1.0 + (abs(dyn) / yclamp))) +end + + +# DP_PACCUM -- Accumulate the data for the parameter fit. + +procedure dp_paccum (psftype, params, npars, mpars, rsq, data, nx, ny, lx, + ly, x, y, sky, h, iter, c, v, temp, chi, sumwt) + +int psftype # analytic point spread function type +real params[ARB] # current function parameter values +int npars # number of function parameters +int mpars # the number of active parameters +real rsq # the fitting radius squared +real data[nx,ARB] # the input image data +int nx, ny # the dimensions of the input image data +int lx, ly # the coordinates of the ll corner of the image data +real x, y # the input/output stellar coordinates +real sky # the input sky value +real h # the input/output height value +int iter # the current iteration +real c[npars,ARB] # accumulation matrix +real v[ARB] # accumulation vector +real temp[ARB] # temporary storage vector +real chi # the chi sum +real sumwt # the number of points sum + +int i, j, k, l +real peak, dx, dy, wt, dhdxc, dhdyc, p, dp +real dp_profile() + +begin + peak = h * dp_profile (psftype, 0.0, 0.0, params, dhdxc, dhdyc, temp, 0) + do j = 1, ny { + dy = real ((ly + j) - 1) - y + do i = 1, nx { + dx = real ((lx + i) - 1) - x + wt = (dx ** 2 + dy ** 2) / rsq + #if (wt >= 1.0) + if (wt >= 0.999998) + next + p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc, temp, 1) + dp = data[i,j] - h * p - sky + do k = 1, mpars + temp[k] = h * temp[k] + chi = chi + (dp / peak) ** 2 + sumwt = sumwt + 1.0 + wt = 5.0 / (5.0 + (wt / (1.0 - wt))) + if (iter >= 4) + wt = wt / (1.0 + abs (20.0 * dp / peak)) + do k = 1, mpars { + v[k] = v[k] + wt * dp * temp[k] + do l = 1, mpars { + c[l,k] = c[l,k] + wt * temp[l] * temp[k] + } + } + } + } +end + + +# DP_FITLT -- Given the analytic function compute the lookup tables. + +int procedure dp_fitlt (dao, im) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image + +int istar, nexp, lx, mx, ly, my, iter, nclean, ndata, fit_saturated, nfit +int nunsat +double volume +pointer apsel, psf, psffit, sp, wimname, wim, data +real datamin, datamax, sumfree, resid, dfdx, dfdy, junk +int dp_dclean(), dp_resana(), dp_ltcompute(), dp_fsaturated() +double dp_pnorm() +pointer immap(), imps2r(), dp_subrast() +real dp_profile(), dp_sweight() +define fitsaturated_ 11 + +begin + # Get some pointers. + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Check to see whether lookup tables are required. + nexp = DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit) + if (nexp <= 0) + return (OK) + + # Return if there are too few stars to fit the lookup tables. + if (DP_PNUM(psf) < nexp) + return (ERR) + + # Determine the number of saturated stars. + nunsat = 0 + do istar = 1, DP_PNUM(psf) { + if (Memi[DP_PSAT(psf)+istar-1] == NO) + next + if ((Memr[DP_PH(psf)+istar-1] * dp_profile (DP_PSFUNCTION(psffit), + 0.0, 0.0, Memr[DP_PSFPARS(psffit)], dfdx, dfdy, junk, 0) + + Memr[DP_APMSKY(apsel)+istar-1]) <= datamax) + next + Memi[DP_PSAT(psf)+istar-1] = YES + Memr[DP_PH(psf)+istar-1] = INDEFR + nunsat = nunsat + 1 + } + nunsat = DP_PNUM(psf) - nunsat + + # Return if there are too few unsaturated psf stars to fit the lookup + # tables. + if (nunsat < nexp) + return (ERR) + + # Allocate memory for computing lookup tables. + call dp_tmempsf (dao) + + # Define some constants. + fit_saturated = DP_SATURATED(dao) + if (IS_INDEFR(DP_MINGDATA(dao))) + datamin = -MAX_REAL + else + datamin = DP_MINGDATA(dao) + if (IS_INDEFR(DP_MAXGDATA(dao))) + datamax = MAX_REAL + else + datamax = DP_MAXGDATA(dao) + sumfree = sqrt (DP_PSUMANA(psf) / (DP_PSUMANA(psf) - (nexp + + 3.0 * DP_PNUM(psf)))) + + # Get the image name. + call smark (sp) + call salloc (wimname, SZ_FNAME, TY_CHAR) + call mktemp ("tmp", Memc[wimname], SZ_FNAME) + + # Open a temporary image to hold the weights. + wim = immap (Memc[wimname], NEW_IMAGE, 0) + IM_NDIM(wim) = 2 + IM_LEN(wim,1) = DP_PNUM(psf) + IM_LEN(wim,2) = DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit) + IM_PIXTYPE(wim) = TY_REAL + + # Compute the constant part of the psf in preparation for normalizing + # the lookup tables. + if (nexp > 1) + call dp_pconst (DP_PSFUNCTION(psffit), Memr[DP_PSFPARS(psffit)], + Memr[DP_PH(psf)], Memr[DP_PCONST(psf)], DP_PSFSIZE(psffit)) + + nfit = 0 + +fitsaturated_ + + # Compute the look-up table. + do iter = 1, DP_NCLEAN(dao) + 1 { + + # Initialize the fitting arrays. + call aclrr (Memr[DP_PV(psf)], nexp) + call aclrr (Memr[DP_PC(psf)], nexp * nexp) + call aclrr (Memr[DP_PSUMN(psf)], DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit)) + call aclrr (Memr[DP_PSUMSQ(psf)], DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit)) + call aclrr (Memr[DP_PSFLUT(psffit)], nexp * DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit)) + + # Loop over the PSF stars. + do istar = 1, DP_PNUM(psf) { + + # Get the weight image. + DP_PSUMW(psf) = imps2r (wim, istar, istar, 1, + DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit)) + call aclrr (Memr[DP_PSUMW(psf)], DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit)) + + # Skip saturated star? + if (IS_INDEFR(Memr[DP_PH(psf)+istar-1])) + next + + # Get the data. + data = dp_subrast (im, Memr[DP_PXCEN(psf)+istar-1], + Memr[DP_PYCEN(psf)+istar-1], DP_PSFRAD(dao), lx, mx, + ly, my) + + # Is the star off the image? + if (lx > IM_LEN(im,1) || ly > IM_LEN(im,2) || mx < 1 || + my < 1) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d is outside the image\n") + call pargi (Memi[DP_APID(apsel)+istar-1]) + } + next + } + + # Clean bad pixels outside the fitting radius but inside + # the psf radius from the subraster. + nclean = dp_dclean (Memr[data], (mx - lx + 1), (my - ly + 1), + lx, ly, Memr[DP_PXCEN(psf)+istar-1], Memr[DP_PYCEN(psf)+ + istar-1], DP_FITRAD(dao), datamin, datamax) + + # Subtract the analytic part of the fit from the data + # and compute the normalized residual of the star. + ndata = dp_resana (im, DP_PSFUNCTION(psffit), + Memr[DP_PSFPARS(psffit)], Memr[data], (mx - lx + 1), + (my - ly + 1), lx, ly, Memr[DP_PXCEN(psf)], + Memr[DP_PYCEN(psf)], Memr[DP_APMSKY(apsel)], + Memr[DP_PH(psf)], DP_PNUM(psf), istar, DP_PSFRAD(dao), + DP_FITRAD(dao), datamax, resid) + + # Compute the proper weight for the star. + if (IS_INDEFR(Memr[DP_PH(psf)+istar-1]) || + Memr[DP_PH(psf)+istar-1] <= 0.0) { + Memr[DP_PWEIGHT(psf)+istar-1] = 0.0 + } else if (Memi[DP_PSAT(psf)+istar-1] == YES) { + Memr[DP_PWEIGHT(psf)+istar-1] = 0.5 * + Memr[DP_PH(psf)+istar-1] / Memr[DP_PH(psf)] + } else { + Memr[DP_PWEIGHT(psf)+istar-1] = (Memr[DP_PH(psf)+ + istar-1] / Memr[DP_PH(psf)]) * dp_sweight (resid, + sumfree, DP_PSIGANA(psf)) + } + + # Compute the expansion vector. + call dp_eaccum (Memr[DP_PXCEN(psf)+istar-1], + Memr[DP_PYCEN(psf)+istar-1], DP_PSFX(psffit), + DP_PSFY(psffit), Memr[DP_PTMP(psf)], DP_NVLTABLE(psffit), + DP_NFEXTABLE(psffit)) + + # Compute the contribution to the lookup table of the + # particular star. + call dp_ltinterp (Memr[data], (mx - lx + 1), (my - ly + 1), + lx, ly, Memr[DP_PXCEN(psf)+istar-1], Memr[DP_PYCEN(psf)+ + istar-1], Memr[DP_APMSKY(apsel)+istar-1], + Memr[DP_PH(psf)+istar-1] / Memr[DP_PH(psf)], + Memr[DP_PWEIGHT(psf)+istar-1], Memr[DP_PSUMN(psf)], + Memr[DP_PSUMW(psf)], Memr[DP_PSUMSQ(psf)], + Memr[DP_PSIGMA(psf)], Memr[DP_POLDLUT(psf)], + Memr[DP_PTMP(psf)], Memr[DP_PSFLUT(psffit)], nexp, + DP_PSFSIZE(psffit), datamax, iter, DP_NCLEAN(dao) + 1) + + call mfree (data, TY_REAL) + } + + call imflush (wim) + + # Compute the lookup table. + if (dp_ltcompute (Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], + DP_PNUM(psf), DP_PSFX(psffit), DP_PSFY(psffit), + Memr[DP_PSUMN(psf)], wim, Memr[DP_PSFLUT(psffit)], + Memr[DP_PC(psf)], Memr[DP_PTMP(psf)], Memr[DP_PV(psf)], + nexp, DP_PSFSIZE(psffit), DP_NVLTABLE(psffit), + DP_NFEXTABLE(psffit)) == ERR) { + call sfree (sp) + return (ERR) + } + + + # Compute the standard deviation arrays for the next pass. + if (iter < (DP_NCLEAN(dao) + 1)) { + if (nunsat <= nexp) + break + call amovr (Memr[DP_PSFLUT(psffit)], Memr[DP_POLDLUT(psf)], + nexp * DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit)) + call dp_stdcompute (Memr[DP_PSUMN(psf)], Memr[DP_PSUMSQ(psf)], + Memr[DP_PSIGMA(psf)], DP_PSFSIZE(psffit), + DP_PSFSIZE(psffit), nexp) + } + + } + + if (nexp > 1) { + + # Accumulate the v vector. + call dp_vaccum (Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], + Memr[DP_PH(psf)], Memr[DP_PWEIGHT(psf)], DP_PNUM(psf), + DP_PSFX(psffit), DP_PSFY(psffit), Memr[DP_PTMP(psf)], + Memr[DP_PV(psf)], nexp, DP_NVLTABLE(psffit), + DP_NFEXTABLE(psffit)) + + # Compute the constant part of the psf. + volume = dp_pnorm (Memr[DP_PCONST(psf)], Memr[DP_PSIGMA(psf)], + DP_PSFSIZE(psffit)) + + # Normalize lookup tables. + call dp_ltnorm (Memr[DP_PCONST(psf)], Memr[DP_PV(psf)], + Memr[DP_PSFLUT(psffit)], Memr[DP_PSIGMA(psf)], nexp, + DP_PSFSIZE(psffit), volume) + + } + + # Make a copy of the psf. + call dp_pcopy (Memr[DP_PSFLUT(psffit)], Memr[DP_POLDLUT(psf)], + DP_PSFSIZE(psffit), DP_PSFSIZE(psffit), nexp) + + # Include the saturated psf stars in the fit. + if (fit_saturated == YES) { + nfit = dp_fsaturated (dao, im, Memi[DP_APID(apsel)], + Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], Memr[DP_PH(psf)], + Memr[DP_APMSKY(apsel)], Memi[DP_PSAT(psf)], DP_PNUM(psf)) + fit_saturated = NO + if (nfit > 0) { + nunsat = nunsat + nfit + if (nexp > 1) + call aaddr (Memr[DP_PCONST(psf)], Memr[DP_POLDLUT(psf)], + Memr[DP_PCONST(psf)], DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit)) + goto fitsaturated_ + } + } + + # Cleanup the temporary images and arrays. + call imunmap (wim) + call imdelete (Memc[wimname]) + call sfree (sp) + + return (OK) +end + + +# DP_DCLEAN -- Clean bad pixels that are outside the fitting radius from +# the data. Note that the star must not be considered to be saturated to +# arrive at this point. + +int procedure dp_dclean (data, nx, ny, lx, ly, x, y, fitrad, datamin, datamax) + +real data[nx,ARB] # the input image data +int nx, ny # the dimensions of the input image data +int lx, ly # the coordinates of the ll image corner +real x, y # the input/output stellar coordinates +real fitrad # the fitting radius +real datamin # the min good data value +real datamax # the max good data value + +bool redo +int i, j, l, k, nclean +real frad2, dy2, dr2, sumd, sumn + +begin + nclean = 0 + repeat { + + redo = false + frad2 = fitrad ** 2 + + do j = 1, ny { + + dy2 = (real ((ly - 1) + j) - y) ** 2 + if (dy2 < frad2) + next + do i = 1, nx { + + if (data[i,j] >= datamin && data[i,j] <= datamax) + next + dr2 = dy2 + (real ((lx - 1) + i) - x) ** 2 + if (dr2 < frad2) + next + + sumd = 0.0 + sumn = 0.0 + do l = max (1, j-1), min (ny, j+2) { + do k = max (1, i-1), min (nx, i+2) { + if (data[k,l] < datamin || data[k,l] > datamax) + next + sumd = sumd + data[k,l] + sumn = sumn + 1.0 + } + } + if (sumn < 2.5) + redo = true + else { + nclean = nclean + 1 + data[i,j] = sumd / sumn + } + } + } + + } until (! redo) + + return (nclean) +end + + +# DP_RESANA -- Compute the residuals from the analytic function. + +int procedure dp_resana (im, psftype, params, data, nx, ny, lx, ly, + x, y, sky, h, nstar, psfstar, psfrad, fitrad, maxgdata, resid) + +pointer im # the input image descriptor +int psftype # analytic point spread function type +real params[ARB] # current function parameter values +real data[nx,ARB] # the input image data +int nx, ny # the dimensions of the input image data +int lx, ly # the coordinates of the ll image corner +real x[ARB] # the input x coords of the psf stars +real y[ARB] # the input y coords of the psf stars +real sky[ARB] # the input sky values of the psf stars +real h[ARB] # the input height values of the psf stars +int nstar # the number of psf stars +int psfstar # the psf star in question +real psfrad # the psf radius +real fitrad # the fitting radius +real maxgdata # the maximum good data value +real resid # standard deviation of fit + +int i, j, istar, rx, ry, x1, x2, y1, y2, nresid +real frad2, dx, dy, dy2, dr2, p, dhdxc, dhdyc, junk +int dp_lsubrast() +real dp_profile() + +begin + frad2 = fitrad ** 2 + rx = lx + nx - 1 + ry = ly + ny - 1 + + resid = 0.0 + nresid = 0 + do istar = 1, nstar { + + # Check for saturation. + if (IS_INDEFR(h[istar])) + next + + # Does the subraster of another PSF star overlap the current + # subraster ?. + if (dp_lsubrast (im, x[istar], y[istar], psfrad, x1, x2, + y1, y2) == ERR) + next + if (x2 < lx || y2 < ly || x1 > rx || y1 > ry) + next + + # Check the limits of overlap. + if (x1 < lx) + x1 = lx + if (x2 > rx) + x2 = rx + if (y1 < ly) + y1 = ly + if (y2 > ry) + y2 = ry + + # Subract off the analytic part of the fits and accumulate + # the residuals for the psf star. + do j = y1 - ly + 1, y2 - ly + 1 { + dy = real ((ly - 1) + j) - y[istar] + dy2 = dy ** 2 + do i = x1 - lx + 1, x2 - lx + 1 { + if (data[i,j] > maxgdata) + next + dx = real ((lx - 1) + i) - x[istar] + p = dp_profile (psftype, dx, dy, params, dhdxc, dhdyc, + junk, 0) + data[i,j] = data[i,j] - h[istar] * p + if (istar != psfstar) + next + dr2 = dy2 + dx ** 2 + if (dr2 >= frad2) + next + resid = resid + (data[i,j] - sky[istar]) ** 2 + nresid = nresid + 1 + } + } + } + + if (nresid <= 0) + resid = 0.0 + else + resid = sqrt (resid / nresid) / (h[psfstar] * dp_profile (psftype, + 0.0, 0.0, params, dhdxc, dhdyc, junk, 0)) + + return (nresid) +end + + +# DP_SWEIGHT -- Compute the weight for the star. + +real procedure dp_sweight (resid, sumfree, sumana) + +real resid # normalized residual wrt analytic fit +real sumfree # number of degrees of freedom +real sumana # number of points contributing to analytic fit + +real weight + +begin + weight = resid * sumfree / sumana + weight = 1.0 / (1.0 + (weight / 2.0) ** 2) + return (weight) +end + + +# DP_EACCUM -- Calcuate the expansion vector for a single PSF star. + +procedure dp_eaccum (x, y, xmid, ymid, junk, nvexp, nfexp) + +real x, y # the stellar coordinates +real xmid, ymid # the psf coordinates +real junk[ARB] # temporary storage vector +int nvexp # the number of variable psf look-up tables +int nfexp # the number of pixel expansion tables + +int j + +begin + # The variable psf terms. + switch (nvexp) { + case 1: + junk[1] = 1.0 + case 3: + junk[1] = 1.0 + junk[2] = ((x - 1.0) / xmid) - 1.0 + junk[3] = ((y - 1.0) / ymid) - 1.0 + case 6: + junk[1] = 1.0 + junk[2] = ((x - 1.0) / xmid) - 1.0 + junk[3] = ((y - 1.0) / ymid) - 1.0 + junk[4] = (1.5 * (junk[2] ** 2)) - 0.5 + junk[5] = junk[2] * junk[3] + junk[6] = (1.5 * (junk[3] ** 2)) - 0.5 + } + + # The fractional pixel expansion terms if any. + if (nfexp > 0) { + j = nvexp + 1 + junk[j] = 2.0 * (x - real (nint(x))) + j = j + 1 + junk[j] = 2.0 * (y - real (nint(y))) + j = j + 1 + junk[j] = (1.5 * (junk[j-2] ** 2)) - 0.5 + j = j + 1 + junk[j] = junk[j-3] * junk[j-2] + j = j + 1 + junk[j] = (1.5 * (junk[j-3] ** 2)) - 0.5 + } + +end + + +# DP_LTINTERP -- Compute the contribution to the lookup table of a single +# PSF star. + +procedure dp_ltinterp (data, nx, ny, lx, ly, x, y, sky, hratio, weight, sumn, + sumw, sumsq, sig, old, temp, psflut, nexp, nxlut, maxgdata, iter, niter) + +real data[nx,ARB] # the input image data +int nx, ny # the dimensions of the input image data +int lx, ly # the coordinates of the ll image corner +real x, y # the input/output stellar coordinates +real sky # sky value for star +real hratio # scale factor for star +real weight # weight for the star +real sumn[nxlut,ARB] # number of points +real sumw[nxlut,ARB] # sum of the weights +real sumsq[nxlut,ARB] # sum of the residuals +real sig[nxlut,ARB] # residuals of previous iteration +real old[nexp,nxlut,ARB] # old lookup table +real temp[ARB] # the single star expansion vector +real psflut[nexp,nxlut,ARB] # the psf lookup table +int nexp, nxlut # the dimensions of the lookup table +real maxgdata # maximum good data value +int iter # the current iteration +int niter # the maximum number of iterations + +bool omit +int i, j, k, kx, ky, ii, jj, middle, jysq, irsq, midsq +real jy, ix, dx, dy, diff, dfdx, dfdy, wt, oldsum +real bicubic() + +begin + middle = (nxlut + 1) / 2 + midsq = middle ** 2 + do j = 1, nxlut { + jysq = (j - middle) ** 2 + if (jysq > midsq) + next + jy = y + real (j - (nxlut + 1) / 2) / 2.0 - real (ly - 1) + ky = int (jy) + if (ky < 2 || (ky + 2) > ny) + next + dy = jy - real (ky) + do i = 1, nxlut { + irsq = jysq + (i - middle) ** 2 + if (irsq > midsq) + next + ix = x + real (i - (nxlut + 1) / 2) / 2.0 - real (lx - 1) + kx = int (ix) + if (kx < 2 || (kx + 2) > nx) + next + omit = false + do jj = ky - 1, ky + 2 { + do ii = kx - 1, kx + 2 { + if (data[ii,jj] <= maxgdata) + next + omit = true + } + } + if (omit) + next + dx = ix - real (kx) + diff = (bicubic (data[kx-1,ky-1], nx, dx, dy, dfdx, dfdy) - + sky) / hratio + if (iter == 1 || sig[i,j] <= 0.0) { + wt = 1.0 + } else { + oldsum = 0.0 + do k = 1, nexp + oldsum = oldsum + old[k,i,j] * temp[k] + if ((iter - 1) <= max (3, (niter - 1) / 2)) + wt = 1.0 / (1.0 + abs (diff - oldsum) / sig[i,j]) + else + wt = 1.0 / (1.0 + ((diff - oldsum) / sig[i,j]) ** 2) + } + wt = wt * weight + sumn[i,j] = sumn[i,j] + 1.0 + sumw[i,j] = wt + sumsq[i,j] = sumsq[i,j] + abs (diff) + psflut[1,i,j] = psflut[1,i,j] + wt * diff + if (nexp <= 1) + next + do k = 2, nexp + psflut[k,i,j] = psflut[k,i,j] + (temp[k] * wt * diff) + } + } +end + + +# DP_LTCOMPUTE -- Compute the lookup table. + +int procedure dp_ltcompute (x, y, npsf, xmid, ymid, sumn, wim, psflut, c, + junk, v, nexp, nxlut, nvexp, nfexp) + +real x[ARB] # array of psf star x coordinates +real y[ARB] # array of psf star y coordinates +int npsf # the number of psf stars +real xmid, ymid # the mid-point of the psf +real sumn[nxlut,ARB] # number of points +pointer wim # pointer to the temporary weight image +real psflut[nexp,nxlut,ARB] # the psf lookup table +real c[nexp,ARB] # the expansion matrix +real junk[ARB] # temporary junk vector +real v[ARB] # temporary vector +int nexp, nxlut # the size of the lookup table +int nvexp, nfexp # size of the expansion look-up tables + +int i, j, k, l, line, istar, ier, middle, midsq, jysq, irsq +pointer sumw +real weight +pointer imgs2r() + +begin + middle = (nxlut + 1) / 2 + midsq = middle ** 2 + do j = 1, nxlut { + jysq = (j - middle) ** 2 + if (jysq > midsq) + next + do i = 1, nxlut { + irsq = (i - middle) ** 2 + jysq + if (irsq > midsq) + next + if (nint (sumn[i,j]) < nexp) + return (ERR) + line = i + (j - 1) * nxlut + sumw = imgs2r (wim, 1, npsf, line, line) + do k = 1, nexp + do l = 1, nexp + c[k,l] = 0.0 + do istar = 1, npsf { + weight = Memr[sumw+istar-1] + if (weight <= 0.0) + next + call dp_eaccum (x[istar], y[istar], xmid, ymid, junk, + nvexp, nfexp) + do k = 1, nexp + do l = 1, nexp + c[k,l] = c[k,l] + weight * junk[k] * junk[l] + } + call invers (c, nexp, nexp, ier) + call mvmul (c, nexp, nexp, psflut[1,i,j], v) + do k = 1, nexp + psflut[k,i,j] = v[k] + } + } + + return (OK) +end + + +# DP_STDCOMPUTE -- Estimate the standard deviation of the fit. + +procedure dp_stdcompute (sumn, sumsq, sigma, nx, ny, nexp) + +real sumn[nx,ARB] # number of point +real sumsq[nx,ARB] # sum of the standard deviations +real sigma[nx,ARB] # output standard deviation array +int nx, ny # size of the arrays +int nexp # number of expansion vectors + +int i, j + +begin + do j = 1, ny { + do i = 1, nx { + if (sumn[i,j] <= nexp) + sigma[i,j] = 0.0 + else + sigma[i,j] = 1.2533 * sumsq[i,j] / sqrt (sumn[i,j] * + (sumn[i,j] - nexp)) + } + } +end + + +# DP_VACCUM -- Accumulate the expansion vector. + +procedure dp_vaccum (x, y, h, weight, nstar, xmid, ymid, junk, avetrm, + nexp, nvexp, nfexp) + +real x[ARB] # the stellar x coordinates +real y[ARB] # the stellar y coordinates +real h[ARB] # the stellar heights +real weight[ARB] # the stellar weights +int nstar # number of stars +real xmid, ymid # the psf coordinates +real junk[ARB] # temporary storage vector +real avetrm[ARB] # the expansion vector +int nexp # the total number of expansion terms +int nvexp # number of variable expansion terms +int nfexp # number of fractional pixel expansion terms + +int k, j + +begin + # Zero the accumulation vector + do j = 1, nexp + avetrm[j] = 0.0 + + do k = 1, nstar { + if (IS_INDEFR(h[k])) + next + + # The variable psf expansion terms. + switch (nvexp) { + case 1: + junk[1] = 1.0 + case 3: + junk[1] = 1.0 + junk[2] = ((x[k] - 1.0) / xmid) - 1.0 + junk[3] = ((y[k] - 1.0) / ymid) - 1.0 + case 6: + junk[1] = 1.0 + junk[2] = ((x[k] - 1.0) / xmid) - 1.0 + junk[3] = ((y[k] - 1.0) / ymid) - 1.0 + junk[4] = (1.5 * (junk[2] ** 2)) - 0.5 + junk[5] = junk[2] * junk[3] + junk[6] = (1.5 * (junk[3] ** 2)) - 0.5 + } + + # The fractional pixel expansion terms if any. + if (nfexp > 0) { + j = nvexp + 1 + junk[j] = 2.0 * (x[k] - real (nint(x[k]))) + j = j + 1 + junk[j] = 2.0 * (y[k] - real (nint(y[k]))) + j = j + 1 + junk[j] = (1.5 * (junk[j-2] ** 2)) - 0.5 + j = j + 1 + junk[j] = junk[j-3] * junk[j-2] + j = j + 1 + junk[j] = (1.5 * (junk[j-3] ** 2)) - 0.5 + } + + # Accumulate the expansion vector. + do j = 1, nexp + avetrm[j] = avetrm[j] + weight[k] * junk[j] + } + + # Average the expansion vector. + do j = nexp, 1, -1 + avetrm[j] = avetrm[j] / avetrm[1] +end + + +# DP_PCONST -- Compute the analytic part of the psf. + +procedure dp_pconst (psftype, params, hpsf1, ana, nxlut) + +int psftype # the type of analytic psf function +real params[ARB] # the analytic parameters array +real hpsf1 # the height of the psf function +real ana[nxlut,ARB] # the computed analytic function +int nxlut # the size of the constant lookup table + +int i, j, middle, midsq, dj2, dr2 +real dx, dy, dfdx, dfdy, junk +real dp_profile() + +begin + # Compute the constant part of the psf. + middle = (nxlut + 1) / 2 + midsq = middle ** 2 + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + dy = real (j - middle) / 2.0 + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + dx = real (i - middle) / 2.0 + ana[i,j] = hpsf1 * dp_profile (psftype, dx, dy, params, + dfdx, dfdy, junk, 0) + } + } +end + + +# DP_PNORM -- Compute the psf normalization parameters. + +double procedure dp_pnorm (ana, pixels, nxlut) + +real ana[nxlut,ARB] # the computed analytic function +real pixels[ARB] # pixel storage array +int nxlut # the size of the constant lookup table + +int i, j, middle, midsq, edgesq, npts, dj2, dr2 +double vol +real median +real pctile() + +begin + # Ensure that the profile which will be added and subtracted + # from the various lookup tables has a median value of zero + # around the rim + + middle = (nxlut + 1) / 2 + midsq = middle ** 2 + edgesq = (middle - 2) ** 2 + npts = 0 + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + if (dr2 < edgesq) + next + npts = npts + 1 + pixels[npts] = ana[i,j] + } + } + median = pctile (pixels, npts, (npts+1) / 2) + + vol = 0.0d0 + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + ana[i,j] = ana[i,j] - median + vol = vol + double (ana[i,j]) + } + } + + return (vol) +end + + +# DP_LTNORM -- Compute the final lookup table. + +procedure dp_ltnorm (ana, v, psflut, pixels, nexp, nxlut, volume) + +real ana[nxlut,ARB] # analytic part of the profile +real v[ARB] # the expansion vector +real psflut[nexp,nxlut,ARB] # the psf lookup table +real pixels[ARB] # scratch array for determining median +int nexp, nxlut # the size of the lookup table +double volume # total flux in the constant psf + +int i, j, k, middle, midsq, edgesq, npts, dj2, dr2 +double sum +real median, dmedian +real pctile() + +begin + # Ensure that the psf which will be added and subtracted from the + # various lookup tables has a median value of zero around the rim. + + middle = (nxlut + 1) / 2 + midsq = middle ** 2 + edgesq = (middle - 2) ** 2 + do k = 2, nexp { + + npts = 0 + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + if (dr2 < edgesq) + next + npts = npts + 1 + pixels[npts] = psflut[k,i,j] + } + } + + median = pctile (pixels, npts, (npts+1) / 2) + dmedian = v[k] * median + + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + psflut[k,i,j] = psflut[k,i,j] - median + psflut[1,i,j] = psflut[1,i,j] + dmedian + } + } + } + + # Determine the net volume of each of the higher order PSF + # tables and force it to zero by subtracting a scaled copy + # of the constant psf. Scale the part that has been subtracted + # off by the mean polynomial term and add it in to the constant + # part of the psf so that at the centroid of the psf stars + # positions the psf remains unchanged. + + do k = 2, nexp { + + sum = 0.0d0 + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + sum = sum + double (psflut[k,i,j]) + } + } + + median = real (sum / volume) + dmedian = v[k] * median + + do j = 1, nxlut { + dj2 = (j - middle) ** 2 + if (dj2 > midsq) + next + do i = 1, nxlut { + dr2 = (i - middle) ** 2 + dj2 + if (dr2 > midsq) + next + psflut[k,i,j] = psflut[k,i,j] - median * ana[i,j] + psflut[1,i,j] = psflut[1,i,j] + dmedian * ana[i,j] + } + } + } +end + + +# DP_FSATURATED -- Fit the saturated stars. + +int procedure dp_fsaturated (dao, im, id, xcen, ycen, h, sky, sat, nstars) + +pointer dao # pointer to the main daophot structure +pointer im # pointer to the input image +int id[ARB] # array of stellar ids +real xcen[ARB] # array of stellar y coordinates +real ycen[ARB] # array of stellar y coordinates +real h[ARB] # array of stellar amplitudes +real sky[ARB] # array of sky values +int sat[ARB] # array of saturation indicators +int nstars # number of stars + +int nfit, nsat, istar, lowx, lowy, nxpix, nypix, recenter, fitsky +int clipexp, ier, niter +pointer psf, psffit, subim, psflut +real x, y, dx, dy, skyval, scale, errmag, chi, sharp, cliprange +int dp_pkfit() +pointer dp_gsubrast() + +begin + # Get some pointers. + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Allocate memory. + call dp_pksetup (dao) + call dp_mempk (dao, 3) + + # Save the default values of some critical parameters. + recenter = DP_RECENTER(dao) + fitsky = DP_FITSKY(dao) + psflut = DP_PSFLUT(psffit) + clipexp = DP_CLIPEXP(dao) + cliprange = DP_CLIPRANGE(dao) + + # Set some fitting parameters. + DP_RECENTER(dao) = YES + DP_FITSKY(dao) = NO + DP_CLIPEXP(dao) = 8 + DP_CLIPRANGE(dao) = 2.5 + DP_PSFLUT(psffit) = DP_POLDLUT(psf) + + nfit = 0 + nsat = 0 + do istar = 1, nstars { + + # Skip saturated stars. + if (sat[istar] == NO) + next + nsat = nsat + 1 + + # Get the data. + subim = dp_gsubrast (im, xcen[istar], ycen[istar], DP_PSFRAD(dao), + lowx, lowy, nxpix, nypix) + if (subim == NULL) + next + + # Set the intial values for the fit parameters. + x = xcen[istar] - lowx + 1.0 + y = ycen[istar] - lowy + 1.0 + dx = (xcen[istar] - 1.0) / DP_PSFX(psffit) - 1.0 + dy = (ycen[istar] - 1.0) / DP_PSFY(psffit) - 1.0 + skyval = sky[istar] + scale = 3.0 + + # Fit the star. + ier = dp_pkfit (dao, Memr[subim], nxpix, nypix, 0.5 * + DP_PSFRAD(dao), x, y, dx, dy, scale, skyval, errmag, chi, + sharp, niter) + + # Compute the fit parameters. + if (ier != PKERR_OK) { + scale = INDEFR + errmag = INDEFR + niter = 0 + chi = INDEFR + sharp = INDEFR + } else { + nfit = nfit + 1 + xcen[istar] = x + lowx - 1.0 + ycen[istar] = y + lowy - 1.0 + h[istar] = scale * h[1] + errmag = 1.085736 * errmag / scale + if (errmag >= 2.0) + errmag = INDEFR + scale = DP_PSFMAG(psffit) - 2.5 * log10 (scale) + } + + if (DP_VERBOSE(dao) == YES) { + if (nsat == 1) + call printf ("\nFit for saturated stars\n") + call printf ( + " %6d %7.2f %7.2f %8.3f %6.3f %3d %7.2f %7.2f\n") + call pargi (id[istar]) + call pargr (xcen[istar]) + call pargr (ycen[istar]) + call pargr (scale) + call pargr (errmag) + call pargi (niter) + call pargr (chi) + call pargr (sharp) + } + } + + if (DP_VERBOSE(dao) == YES && nsat > 0) + call printf ("\n") + + # Restore the default values of some critical parameters. + DP_RECENTER(dao) = recenter + DP_FITSKY(dao) = fitsky + DP_CLIPEXP(dao) = clipexp + DP_CLIPRANGE(dao) = cliprange + DP_PSFLUT(psffit) = psflut + + # Free memory. + call dp_pkclose (dao) + + return (nfit) +end + + +# DP_PCOPY -- Make a copy of the psf in correct storage order. + +procedure dp_pcopy (inlut, outlut, nx, ny, nexp) + +real inlut[nexp,nx,ny] # the input look-up table +real outlut[nx,ny,nexp] # the output look-up table +int nx,ny,nexp # the size of the look-up table + +int k,i,j + +begin + do k = 1, nexp { + do j = 1, ny { + do i = 1, nx { + outlut[i,j,k] = inlut[k,i,j] + } + } + } +end |