diff options
Diffstat (limited to 'noao/digiphot/daophot/psf/dpsubpsf.x')
| -rw-r--r-- | noao/digiphot/daophot/psf/dpsubpsf.x | 183 |
1 files changed, 183 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/psf/dpsubpsf.x b/noao/digiphot/daophot/psf/dpsubpsf.x new file mode 100644 index 00000000..8f15867a --- /dev/null +++ b/noao/digiphot/daophot/psf/dpsubpsf.x @@ -0,0 +1,183 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_SUBPSF -- Add the star at the given position to the PSF if it exists and +# passes the selection criteria. + +int procedure dp_subpsf (dao, im, x, y, idnum, gd, mgd, showplots) + +pointer dao # pointer to daophot structure +pointer im # pointer to image +real x, y # position of proposed PSF star +int idnum # id number of desired star +pointer gd # pointer to the graphics stream +pointer mgd # pointer to the metacode descriptor +bool showplots # show plots? + +real tx, ty +pointer srim +int x1, x2, y1, y2, starnum, saturated +bool star_ok + +real dp_statr() +pointer dp_psubrast() +int dp_locstar(), dp_idstar(), dp_pstati(), dp_issat() + +begin + # Convert coordinates for display. + if (showplots) + call dp_ltov (im, x, y, tx, ty, 1) + else + call dp_wout (dao, im, x, y, tx, ty, 1) + + # Check that the position of the star is within the image. + if (idnum == 0 && (x < 1.0 || x > real (IM_LEN(im,1)) || y < 1.0 || y > + real (IM_LEN(im,2)))) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star at %g,%g is outside the image\n") + call pargr (tx) + call pargr (ty) + } + return (ERR) + } + + # Find the star in the aperture photometry list + if (idnum == 0) + starnum = dp_locstar (dao, im, x, y) + else + starnum = dp_idstar (dao, im, idnum) + if (starnum == 0) { + if (DP_VERBOSE(dao) == YES) { + if (idnum > 0) { + call printf ("Star %d not found in the photometry file\n") + call pargi (idnum) + } else { + call printf ( + "Star at %g,%g not found in the photometry file\n") + call pargr (tx) + call pargr (ty) + } + } + return (ERR) + } else if (starnum < 0 || starnum > dp_pstati (dao, PNUM)) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d not found in the PSF star list\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + # Check to see if the star is saturated. + if (dp_issat (dao, starnum) == YES) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d is saturated\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + # Get the data subraster, check for saturation and bad pixels, + # and compute the min and max data values inside the subraster. + srim = dp_psubrast (dao, im, -MAX_REAL, MAX_REAL, x1, x2, y1, y2, + saturated) + if (srim == NULL) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d error reading data subraster\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + # Now let us do the subtraction. + call dp_spstar (dao, Memr[srim], (x2 - x1 + 1), (y2 - y1 + 1), + x1, y1, starnum, dp_statr (dao, PSFRAD) ** 2) + + # Now let's look at the extracted subraster. + if (showplots) { + call dp_showpsf (dao, im, Memr[srim], (x2 - x1 + 1), (y2 - y1 + 1), + x1, y1, gd, star_ok) + } else + star_ok = true + + if (star_ok) { + if (mgd != NULL) + call dp_plotpsf (dao, im, Memr[srim], (x2 - x1 + 1), + (y2 - y1 + 1), x1, y1, mgd) + if (DP_VERBOSE(dao) == YES) { + call printf ("PSF star %d saved by user\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + call mfree (srim, TY_REAL) + return (ERR) + } + + # Delete the star from the list by moving it to the position + # currently occupied by PNUM and moving everything else up. + call dp_pfreorder (dao, dp_pstati (dao, CUR_PSF), + dp_pstati (dao, PNUM)) + + # Decrement the list of psf stars. + call dp_pseti (dao, PNUM, dp_pstati (dao, PNUM) - 1) + + # Print message. + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d has been deleted from the PSF star list\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + + call mfree (srim, TY_REAL) + return (OK) +end + + +# DP_SPSTAR -- Subtract the fitted star from the data subraster. + +procedure dp_spstar (dao, data, nx, ny, x1, y1, starno, psfradsq) + +pointer dao # pointer to the daophot structure +real data[nx,ARB] # the data subraster +int nx, ny # the dimensions of the data subraster +int x1, y1 # the coordinates of the ll pixel +int starno # the index of the star in question +real psfradsq # the psf radius squared + +int i, j +pointer psf, psffit +real xstar, ystar, scale, deltax, deltay, dx, dy, dysq, rsq, svdx, svdy +real dp_usepsf() + +begin + # Define some daophot pointers + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + # Get the star position and scale factor. + xstar = Memr[DP_PXCEN(psf)+starno-1] + ystar = Memr[DP_PYCEN(psf)+starno-1] + deltax = (xstar - 1.0) / DP_PSFX(psffit) - 1.0 + deltay = (ystar - 1.0) / DP_PSFY(psffit) - 1.0 + xstar = xstar - x1 + 1 + ystar = ystar - y1 + 1 + scale = Memr[DP_PH(psf)+starno-1] / Memr[DP_PH(psf)] + + do j = 1, ny { + dy = real (j) - ystar + dysq = dy ** 2 + do i = 1, nx { + dx = real (i) - xstar + rsq = dx ** 2 + dysq + if (rsq >= psfradsq) + next + data[i,j] = data[i,j] - scale * + dp_usepsf (DP_PSFUNCTION(psffit), dx, dy, + DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)], + Memr[DP_POLDLUT(psf)], DP_PSFSIZE(psffit), + DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), + deltax, deltay, svdx, svdy) + } + } + + call alimr (data, nx * ny, DP_CUR_PSFMIN(psf), DP_CUR_PSFMAX(psf)) +end |