diff options
Diffstat (limited to 'noao/digiphot/daophot/psf')
34 files changed, 8216 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/psf/README b/noao/digiphot/daophot/psf/README new file mode 100644 index 00000000..4d390489 --- /dev/null +++ b/noao/digiphot/daophot/psf/README @@ -0,0 +1,2 @@ +This directory contains the routines for making up the PSF to be used by +the fitting routines in the PSF. diff --git a/noao/digiphot/daophot/psf/dpaddstar.x b/noao/digiphot/daophot/psf/dpaddstar.x new file mode 100644 index 00000000..c723a2b3 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpaddstar.x @@ -0,0 +1,188 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_ADDSTAR -- Add the star at the given position to the PSF if it exists and +# passes the selection criteria. + +int procedure dp_addstar (dao, im, x, y, mag, idnum, gd, mgd, showplots) + +pointer dao # pointer to daophot structure +pointer im # pointer to image +real x, y # position of proposed PSF star +real mag # mag 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, logood, higood +pointer srim +int x1, x2, y1, y2, starnum, saturated +bool star_ok + +real dp_statr(), dp_pstatr() +pointer dp_psubrast() +int dp_locstar(), dp_idstar(), dp_stati(), dp_pstati() + +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) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Star %d is too near the edge of the image\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } else if (starnum <= dp_pstati (dao, PNUM)) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d is already a PSF star\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + # Check for INDEF valued sky. + if (IS_INDEFR (dp_pstatr (dao, CUR_PSFSKY))) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d has an undefined sky value\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + + logood = dp_statr (dao, MINGDATA) + if (IS_INDEFR(logood)) + logood = -MAX_REAL + higood = dp_statr (dao, MAXGDATA) + if (IS_INDEFR(higood)) + higood = MAX_REAL + + # 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, logood, higood, x1, x2, y1, y2, saturated) + if (srim == NULL) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Star %d has low bad pixels inside fitrad\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + return (ERR) + } + + # Check for saturation. + if (saturated == YES && dp_stati (dao, SATURATED) == NO) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Star %d has high bad pixels inside fitrad\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + call mfree (srim, TY_REAL) + return (ERR) + } + + # 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 if (saturated == YES) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Warning: Star %d contains high bad pixels inside fitrad\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + star_ok = true + } else if (dp_pstatr (dao, CUR_PSFMIN) < logood || dp_pstatr (dao, + CUR_PSFMAX) > higood) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Warning: Star %d contains bad pixels outside fitrad\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + star_ok = true + } else + star_ok = true + + # The star is rejected by the user. + if (! star_ok) { + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d rejected by user\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + call mfree (srim, TY_REAL) + return (ERR) + } + + # Save the plot in the metacode file. + if (mgd != NULL) + call dp_plotpsf (dao, im, Memr[srim], (x2 - x1 + 1), (y2 - y1 + 1), + x1, y1, mgd) + + # Add the star to the PSF star list by swapping its position with the + # position of the star currently in PNUM + 1. + call dp_aplswap (dao, dp_pstati (dao, CUR_PSF), dp_pstati (dao, + PNUM) + 1) + + # Increment the number of psf stars. + call dp_pseti (dao, PNUM, dp_pstati (dao, PNUM) + 1) + + # Reallocate the fitting array space. + call dp_lmempsf (dao) + + # Enter the new initial values. + call dp_xyhpsf (dao, dp_pstati (dao, PNUM), mag, saturated) + + # Print message. + if (DP_VERBOSE(dao) == YES) { + call printf ("Star %d has been added to the PSF star list\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + call dp_ltov (im, dp_pstatr (dao, CUR_PSFX), + dp_pstatr(dao, CUR_PSFY), tx, ty, 1) + call printf ( + "\tX: %7.2f Y: %7.2f Mag: %7.3f Dmin: %g Dmax: %g\n") + call pargr (tx) + call pargr (ty) + call pargr (dp_pstatr (dao, CUR_PSFMAG)) + call pargr (dp_pstatr (dao, CUR_PSFMIN)) + call pargr (dp_pstatr (dao, CUR_PSFMAX)) + } + + call mfree (srim, TY_REAL) + return (OK) +end diff --git a/noao/digiphot/daophot/psf/dpcontpsf.x b/noao/digiphot/daophot/psf/dpcontpsf.x new file mode 100644 index 00000000..27b6bdfc --- /dev/null +++ b/noao/digiphot/daophot/psf/dpcontpsf.x @@ -0,0 +1,449 @@ +include <error.h> +include <mach.h> +include <gset.h> +include <config.h> +include <xwhen.h> +include <fset.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +define DUMMY 6 +define XCEN 0.5 +define YCEN 0.52 +define EDGE1 0.1 +define EDGE2 0.93 +define SZ_LABEL 10 +define SZ_FMT 20 + +# DP_CONTPSF -- Draw a contour plot of a data subraster containing. The +# data floor and ceiling are set by the user, but the contour interval +# is chosen by the routine. + +procedure dp_contpsf (dao, subras, ncols, nlines, title, gp) + +pointer dao # pointer to DAOPHOT structure +real subras[ncols,nlines] # data subraster +int ncols, nlines # dimesnions of subraster +char title[ARB] # title string +pointer gp # pointer to graphics descriptor + +bool perimeter +char system_id[SZ_LINE], label[SZ_LINE] +int epa, status, old_onint, tcojmp[LEN_JUMPBUF] +int wkid, nset, ncontours, dashpat, nhi +pointer sp, temp, psf +real interval, floor, ceiling, zero, finc, ybot +real vx1, vx2, vy1, vy2, wx1, wx2, wy1, wy2 +real first_col, last_col, first_row, last_row + +bool fp_equalr() +extern dp_conint() +common /tcocom/ tcojmp + +int first +int isizel, isizem, isizep, nrep, ncrt, ilab, nulbll, ioffd +int ioffm, isolid, nla, nlm +real xlt, ybt, side, ext, hold[5] +common /conflg/ first +common /conre4/ isizel, isizem , isizep, nrep, ncrt, ilab, nulbll, + ioffd, ext, ioffm, isolid, nla, nlm, xlt, ybt, side +common /noaolb/ hold + +begin + # Get the pointer to the DAOPHOT PSF fitting substructure. + psf = DP_PSF (dao) + + # First of all, intialize conrec's block data before altering any + # parameters in common. + first = 1 + call conbd + + # Set local variables. + zero = 0.0 + floor = DP_CFLOOR (psf) + ceiling = DP_CCEILING (psf) + nhi = -1 + dashpat = 528 + + # Suppress the contour labelling by setting the common + # parameter "ilab" to zero. + ilab = 0 + + # The floor and ceiling are in absolute units, but the zero shift is + # applied first, so correct the numbers for the zero shift. Zero is + # a special number for the floor and ceiling, so do not change value + # if set to zero. + + if (abs (floor) > EPSILON) + floor = floor - zero + if (abs (ceiling) > EPSILON) + ceiling = ceiling - zero + + # User can specify either the number of contours or the contour + # interval, or let conrec pick a nice number. Set ncontours to 0 + # and encode the FINC param expected by conrec. + + ncontours = 0 + if (ncontours <= 0) { + interval = 0 + if (interval <= 0) + finc = 0 + else + finc = interval + } else + finc = - abs (ncontours) + + # Make a copy of the data and do the contouring on this. + call smark (sp) + call salloc (temp, ncols * nlines, TY_REAL) + call amovr (subras, Memr[temp], nlines * ncols) + + first_col = 1.0 + last_col = real (ncols) + first_row = 1.0 + last_row = real (nlines) + + # Apply the zero point shift. + if (abs (zero) > EPSILON) + call asubkr (Memr[temp], zero, Memr[temp], ncols * nlines) + + # Open device and make contour plot. + call gopks (STDERR) + wkid = 1 + call gclear (gp) + call gopwk (wkid, DUMMY, gp) + call gacwk (wkid) + + # The viewport can be set by the user. If not, the viewport is + # assumed to be centered on the device. In either case, the + # viewport to window mapping is established in dp_map_viewport + # and conrec's automatic mapping scheme is avoided by setting nset=1. + vx1 = 0.10 + vx2 = 0.90 + vy1 = 0.10 + vy2 = 0.90 + call dp_map_viewport (gp, ncols, nlines, vx1, vx2, vy1, vy2, false) + nset = 1 + + perimeter = TRUE + if (perimeter) + # Suppress conrec's plot label generation. + ioffm = 1 + else { + # Draw plain old conrec perimeter, set ioffm = 0 to enable label. + ioffm = 0 + call perim (ncols - 1, 1, nlines - 1, 1) + } + + # Install interrupt exception handler. + call zlocpr (dp_conint, epa) + call xwhen (X_INT, epa, old_onint) + + # Make the contour plot. If an interrupt occurs ZSVJMP is reeentered + # with an error status. + call zsvjmp (tcojmp, status) + if (status == OK) { + call conrec (Memr[temp], ncols, ncols, nlines, floor, ceiling, + finc, nset, nhi, -dashpat) + } else { + call gcancel (gp) + call fseti (STDOUT, F_CANCEL, OK) + } + + # Now find window and output text string title. The window is + # set to the full image coordinates for labelling. + call gswind (gp, first_col, last_col, first_row, last_row) + if (perimeter) + call dp_cperimeter (gp) + + call ggview (gp, wx1, wx2, wy1, wy2) + call gseti (gp, G_WCS, 0) + ybot = min (wy2 + .06, 0.99) + call gtext (gp, (wx1 + wx2) / 2.0, ybot, title, "h=c;v=t;f=b;s=.7") + + # Add system id banner to plot. + call gseti (gp, G_CLIP, NO) + call sysid (system_id, SZ_LINE) + ybot = max (wy1 - 0.08, 0.01) + call gtext (gp, (wx1+wx2)/2.0, ybot, system_id, "h=c;v=b;s=.5") + + if (perimeter) { + if (fp_equalr (hold(5), 1.0)) { + call sprintf (label, SZ_LINE, + "contoured from %g to %g, interval = %g") + call pargr (hold(1)) + call pargr (hold(2)) + call pargr (hold(3)) + } else { + call sprintf (label, SZ_LINE, + "contoured from %g to %g, interval = %g, labels scaled by %g") + call pargr (hold(1)) + call pargr (hold(2)) + call pargr (hold(3)) + call pargr (hold(5)) + } + ybot = max (wy1 - 0.06, .03) + call gtext (gp, (wx1 + wx2) / 2.0, ybot, label, "h=c;v=b;s=.6") + } + + call gswind (gp, first_col, last_col, first_row, last_row) + call gdawk (wkid) + call gclks () + call sfree (sp) +end + + +# DP_CONINT -- Interrupt handler for the task contour. Branches back to +# ZSVJMP in the main routine to permit shutdown without an error message. + +procedure dp_conint (vex, next_handler) + +int vex # virtual exception +int next_handler # not used + +int tcojmp[LEN_JUMPBUF] +common /tcocom/ tcojmp + +begin + call xer_reset() + call zdojmp (tcojmp, vex) +end + + +# DP_CPERIMETER -- draw and annotate the axes drawn around the perimeter +# of the image pixels. The viewport and window have been set by +# the calling procedure. Plotting is done in window coordinates. +# This procedure is called by both crtpict and the ncar plotting routines +# contour and hafton. + +procedure dp_cperimeter (gp) + +pointer gp # graphics descriptor +real xs, xe, ys, ye # WCS coordinates of pixel window + +char label[SZ_LABEL], fmt1[SZ_FMT], fmt2[SZ_FMT], fmt3[SZ_FMT], fmt4[SZ_FMT] +int i, first_col, last_col, first_tick, last_tick, bias +int nchar, first_row, last_row, cnt_step, cnt_label +real dist, kk, col, row, dx, dy, sz_char, cw, xsz, label_pos +real xdist, ydist, xspace, yspace, k[3] +bool ggetb() +int itoc() +real ggetr() +data k/1.0,2.0,3.0/ +errchk ggwind, gseti, gctran, gline, gtext, itoc + +begin + # First, get window coordinates and turn off clipping. + call ggwind (gp, xs, xe, ys, ye) + call gseti (gp, G_CLIP, NO) + + # A readable character width seems to be about 1.mm. A readable + # perimeter seperation seems to be about .80mm. If the physical + # size of the output device is contained in the graphcap file, the + # NDC sizes of these measurements can be determined. If not, + # the separation between perimeter axes equals one quarter character + # width or one quarter percent of frame, which ever is larger, and + # the character size is set to 0.40. + + cw = max (ggetr (gp, "cw"), 0.01) + if (ggetb (gp, "xs")) { + xsz = ggetr (gp, "xs") + dist = .80 / (xsz * 1000.) + sz_char = dist / cw + } else { + # Get character width and calculate perimeter separation. + dist = cw * 0.25 + sz_char = 0.40 + } + + # Convert distance to user coordinates. + call ggscale (gp, xs, ys, dx, dy) + xdist = dist * dx + ydist = dist * dy + + # Generate four possible format strings for gtext. + call sprintf (fmt1, SZ_LINE, "h=c;v=t;s=%.2f") + call pargr (sz_char) + call sprintf (fmt2, SZ_LINE, "h=c;v=b;s=%.2f") + call pargr (sz_char) + call sprintf (fmt3, SZ_LINE, "h=r;v=c;s=%.2f") + call pargr (sz_char) + call sprintf (fmt4, SZ_LINE, "h=l;v=c;s=%.2f") + call pargr (sz_char) + + # Draw inner and outer perimeter + kk = k[1] + do i = 1, 2 { + xspace = kk * xdist + yspace = kk * ydist + call gline (gp, xs - xspace, ys - yspace, xe + xspace, ys - yspace) + call gline (gp, xe + xspace, ys - yspace, xe + xspace, ye + yspace) + call gline (gp, xe + xspace, ye + yspace, xs - xspace, ye + yspace) + call gline (gp, xs - xspace, ye + yspace, xs - xspace, ys - yspace) + kk = k[2] + } + + # Now draw x axis tick marks, along both the bottom and top of + # the picture. First find the endpoint integer pixels. + + first_col = int (xs) + last_col = int (xe) + + # Determine increments of ticks and tick labels for x axis. + cnt_step = 1 + cnt_label = 10 + if (last_col - first_col > 256) { + cnt_step = 10 + cnt_label = 100 + } else if (last_col - first_col < 26) { + cnt_step = 1 + cnt_label = 1 + } + + first_tick = first_col + bias = mod (first_tick, cnt_step) + last_tick = last_col + bias + + do i = first_tick, last_tick, cnt_step { + col = real (i - bias) + call gline (gp, col, ys - k[1] * ydist, col, ys - k[2] * ydist) + call gline (gp, col, ye + k[1] * ydist, col, ye + k[2] * ydist) + + if (mod ((i - bias), cnt_label) == 0) { + + # Label tick mark; calculate number of characters needed. + nchar = 3 + if (int (col) == 0) + nchar = 1 + if (int (col) >= 1000) + nchar = 4 + if (itoc (int(col), label, nchar) <= 0) + label[1] = EOS + + # Position label slightly below outer perimeter. Seperation + # is twenty percent of a character width, in WCS. + label_pos = ys - (k[2] * ydist + (cw * 0.20 * dy)) + call gtext (gp, col, label_pos, label, fmt1) + + # Position label slightly above outer perimeter. + label_pos = ye + (k[2] * ydist + (cw * 0.20 * dy)) + call gtext (gp, col, label_pos, label, fmt2) + } + } + + # Label the y axis tick marks along the left and right sides of the + # picture. First find the integer pixel endpoints. + + first_row = int (ys) + last_row = int (ye) + + # Determine increments of ticks and tick labels for y axis. + cnt_step = 1 + cnt_label = 10 + if (last_row - first_row > 256) { + cnt_step = 10 + cnt_label = 100 + } else if (last_row - first_row < 26) { + cnt_step = 1 + cnt_label = 1 + } + + first_tick = first_row + bias = mod (first_tick, cnt_step) + last_tick = last_row + bias + + do i = first_tick, last_tick, cnt_step { + row = real (i - bias) + call gline (gp, xs - k[1] * xdist, row, xs - k[2] * xdist, row) + call gline (gp, xe + k[1] * xdist, row, xe + k[2] * xdist, row) + + if (mod ((i - bias), cnt_label) == 0) { + + # Label tick mark; calculate number of characters needed + nchar = 3 + if (int (row) == 0) + nchar = 1 + else if (int (row) >= 1000) + nchar = 4 + if (itoc (int(row), label, nchar) <= 0) + label[1] = EOS + + # Position label slightly to the left of outer perimeter. + # Separation twenty percent of a character width, in WCS. + label_pos = xs - (k[2] * xdist + (cw * 0.20 * dx)) + call gtext (gp, label_pos, row, label, fmt3) + + # Position label slightly to the right of outer perimeter. + label_pos = xe + (k[2] * xdist + (cw * 0.20 * dx)) + call gtext (gp, label_pos, row, label, fmt4) + } + } +end + + +# DP_MAP_VIEWPORT -- Set device viewport for contour and hafton plots. If not +# specified by user, a default viewport centered on the device is used. + +procedure dp_map_viewport (gp, ncols, nlines, ux1, ux2, uy1, uy2, fill) + +pointer gp # graphics pointer +int ncols # number of image cols +int nlines # number of image lines +real ux1, ux2, uy1, uy2 # NDC coordinates of requested viewort +bool fill # fill viewport (vs enforce unity aspect ratio?) + +real ncolsr, nlinesr, ratio, aspect_ratio +real x1, x2, y1, y2, ext, xdis, ydis +data ext /0.25/ +bool fp_equalr() +real ggetr() + +begin + ncolsr = real (ncols) + nlinesr = real (nlines) + if (fp_equalr (ux1, 0.0) && fp_equalr (ux2, 0.0) && + fp_equalr (uy1, 0.0) && fp_equalr (uy2, 0.0)) { + + x1 = EDGE1 + x2 = EDGE2 + y1 = EDGE1 + y2 = EDGE2 + + # Calculate optimum viewport, as in NCAR's conrec, hafton. + ratio = min (ncolsr, nlinesr) / max (ncolsr, nlinesr) + if (ratio >= ext) { + if (ncols > nlines) + y2 = (y2 - y1) * nlinesr / ncolsr + y1 + else + x2 = (x2 - x1) * ncolsr / nlinesr + x1 + } + + xdis = x2 - x1 + ydis = y2 - y1 + + # So far, the viewport has been calculated so that equal numbers of + # image pixels map to equal distances in NDC space, regardless of + # the aspect ratio of the device. If the parameter "fill" has been + # set to no, the user wants to compensate for a non-unity aspect + # ratio and make equal numbers of image pixels map to into the same + # physical distance on the device, not the same NDC distance. + + if (! fill) { + aspect_ratio = ggetr (gp, "ar") + if (fp_equalr (aspect_ratio, 0.0)) + aspect_ratio = 1.0 + xdis = xdis * aspect_ratio + } + + ux1 = XCEN - (xdis / 2.0) + ux2 = XCEN + (xdis / 2.0) + uy1 = YCEN - (ydis / 2.0) + uy2 = YCEN + (ydis / 2.0) + } + + # Set window and viewport for WCS 1. + call gseti (gp, G_WCS, 1) + call gsview (gp, ux1, ux2, uy1, uy2) + call gswind (gp, 1.0, ncolsr, 1.0, nlinesr) + call set (ux1, ux2, uy1, uy2, 1.0, ncolsr, 1.0, nlinesr, 1) +end diff --git a/noao/digiphot/daophot/psf/dpdelstar.x b/noao/digiphot/daophot/psf/dpdelstar.x new file mode 100644 index 00000000..22ba579f --- /dev/null +++ b/noao/digiphot/daophot/psf/dpdelstar.x @@ -0,0 +1,112 @@ +include <imhdr.h> +include <mach.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_DELSTAR -- Delete the star at the given position from the list of PSF +# stars if it exists and passes the selection criteria. + +int procedure dp_delstar (dao, im, x, y, idnum, gd, 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 +bool showplots # show plots? + +real tx, ty +pointer srim +int starnum, saturated, x1, x2, y1, y2 +bool star_ok + +pointer dp_psubrast() +int dp_locstar(), dp_idstar(), dp_pstati() + +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) + } + + # Get the data subraster, check for bad pixels and compute the min + # and max. + if (showplots) { + 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)) + } + star_ok = false + } else { + call dp_showpsf (dao, im, Memr[srim], (x2 - x1 + 1), + (y2 - y1 + 1), x1, y1, gd, star_ok) + call mfree (srim, TY_REAL) + } + } else + star_ok = true + + if (star_ok) { + if (DP_VERBOSE(dao) == YES) { + call printf ("PSF star %d saved by user\n") + call pargi (dp_pstati (dao, CUR_PSFID)) + } + 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 number 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)) + } + + return (OK) +end 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 diff --git a/noao/digiphot/daophot/psf/dpispstars.x b/noao/digiphot/daophot/psf/dpispstars.x new file mode 100644 index 00000000..85a4c669 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpispstars.x @@ -0,0 +1,329 @@ +include <fset.h> +include <ctype.h> +include <gset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +define HELPFILE "daophot$psf/mkpsflist.key" + +# DP_IPFSTARS -- Construct a stellar PSF from one or more stars in an image +# frame. + +int procedure dp_ipfstars (dao, im, maxnpsf, lolimit, radius, gd, mgd, id, + mkstars, interactive, showplots) + +pointer dao # pointer to DAOPHOT structure +pointer im # pointer to IRAF image +int maxnpsf # maximum number of psf stars +real lolimit # lower magnitude limit for stars +real radius # the confusion radius for good psf stars +pointer gd # pointer to graphics descriptor +pointer mgd # pointer to the metacode file +pointer id # pointer to image display stream +bool mkstars # marked the selected and deleted psf stars +bool interactive # interactive mode +bool showplots # show the plots + +int key, nxstar, wcs, idnum, istar, ip +pointer apsel, sp, cmd, str +real wx, wy + +real dp_pstatr() +int clgcur(), dp_pqverify(), dp_locstar(), dp_idstar(), dp_nxstar() +int dp_pstati(), ctoi(), dp_addstar(), dp_delstar() + +begin + # Get some pointers + apsel = DP_APSEL(dao) + + # Allocate some working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize some variables. + key = 'a' + nxstar = 0 + call dp_pseti (dao, PNUM, 0) + + # Begin to build the PSF. + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Correct for tv coordinates if necessary. + call dp_vtol (im, wx, wy, wx, wy, 1) + + switch (key) { + + # Quit the interactive cursor loop. + case 'q': + if (interactive) { + if (dp_pqverify() == YES) + break + } else + break + + # Print the help page. + case '?': + if (id == NULL) + call pagefile (HELPFILE, "") + else + call gpagefile (id, HELPFILE, "") + + # Add next candidate star to the PSF. + case 'n': + if (dp_pstati (dao, PNUM) >= maxnpsf) { + call printf ("Number of psf stars is >= to maxnpsf\n") + next + } + idnum = dp_nxstar (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMSKY(apsel)], + Memr[DP_APMAG(apsel)], DP_APNUM(apsel), nxstar, lolimit, + radius) + if (idnum <= 0) { + call printf ("No more good psf stars in photometry list\n") + next + } + if (dp_addstar (dao, im, wx, wy, INDEFR, idnum, gd, mgd, + showplots) == ERR) + next + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), dp_pstatr (dao, + CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # Add the star nearest the cursor position to the PSF. + case 'a': + if (dp_pstati (dao, PNUM) >= maxnpsf) { + call printf ("Number of psf stars is >= to maxnpsf\n") + next + } + if (dp_addstar (dao, im, wx, wy, INDEFR, 0, gd, mgd, + showplots) == ERR) + next + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), dp_pstatr (dao, + CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # Delete the star nearest the cursor position from the PSF. + case 'd': + if (dp_delstar (dao, im, wx, wy, 0, gd, showplots) == ERR) + next + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), dp_pstatr (dao, + CUR_PSFY), GM_CROSS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # List all the current psf stars. + case 'l': + if (interactive) + call dp_listpsf (dao, im) + + # Locate the star in the aperture photometry file and print out + # the photometry. + case 'p': + if (interactive) + istar = dp_locstar (dao, im, wx, wy) + else + next + if (istar > 0) + call dp_pshow (dao, im, istar) + else if (istar == 0) + call printf ("Star not found in the photometry file\n") + else + call printf ( + "Star off or too near the edge of the image\n") + + # Colon command mode. + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + switch (Memc[cmd+ip-1]) { + + case 'p': + if (! interactive) + next + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + istar = dp_locstar (dao, im, wx, wy) + else + istar = dp_idstar (dao, im, idnum) + if (istar > 0) + call dp_pshow (dao, im, istar) + else if (istar == 0) + call printf ( + "Star not found in the photometry file\n") + else + call printf ( + "Star is off or too near the edge of the image.\n") + + case 'a': + if (dp_pstati (dao, PNUM) >= maxnpsf) { + call printf ( + "Number of psf stars is >= to maxnpsf\n") + next + } + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + idnum = 0 + if (dp_addstar (dao, im, wx, wy, INDEFR, idnum, gd, mgd, + showplots) == ERR) + next + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), + dp_pstatr (dao, CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + case 'd': + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + idnum = 0 + if (dp_delstar (dao, im, wx, wy, idnum, gd, showplots) == + ERR) + next + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), + dp_pstatr (dao, CUR_PSFY), GM_CROSS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + default: + call printf ("Unknown keystroke command\n") + } + + default: + call printf ("Unknown keystroke command\n") + } + } + + # Free up memory. + call sfree (sp) + + return (dp_pstati (dao, PNUM)) +end + + +define QUERY "[Hit return to continue, q to quit]" + +# DP_PQVERIFY -- Print a message in the status line asking the user if they +# really want to quit, returning YES if they really want to quit, NO otherwise. + +int procedure dp_pqverify () + +int ch +pointer tty +int getci() +pointer ttyodes() + +begin + # Open terminal and print query. + tty = ttyodes ("terminal") + call ttyclearln (STDOUT, tty) + call ttyso (STDOUT, tty, YES) + call printf (QUERY) + call flush (STDOUT) + + # Get character. + call fseti (STDIN, F_RAW, YES) + if (getci (STDIN, ch) == EOF) + ; + + # Reset and close terminal. + call fseti (STDIN, F_RAW, NO) + call ttyso (STDOUT, tty, NO) + call ttyclearln (STDOUT, tty) + call printf ("\n") + call flush (STDOUT) + call ttycdes (tty) + + # Return YES for the quit command, otherwise NO. + if (ch == 'q') { + return (YES) + } else { + return (NO) + } +end + + + +# DP_NXSTAR -- Select the next psf star form the photometry list. + +int procedure dp_nxstar (ids, xcen, ycen, sky, mag, nstar, nxstar, lolimit, + radius) + +int ids[ARB] # array of star ids +real xcen[ARB] # array of x coordinates +real ycen[ARB] # array of y coordinates +real sky[ARB] # array of sky values +real mag[ARB] # array of magnitudes +int nstar # the number of stars +int nxstar # the current star +real lolimit # lower data limit +real radius # minimum separation + +bool omit +int istar, jstar +real radsq, dy2, dr2 + +begin + radsq = radius * radius + + # Check the first star. + if ((mag[1] > lolimit) && (nxstar == 0)) { + nxstar = 1 + return (ids[1]) + } + + # Loop over the candidate psf stars. + do istar = nxstar + 1, nstar { + + # Test that the candidate psf stars are not saturated and + # are sufficiently far from the edge of the frame. + + if (mag[istar] <= lolimit) + next + + # Text that there are no brighter stars with a distance squared + # of radsq. + omit = false + do jstar = 1, istar - 1 { + dy2 = abs (ycen[jstar] - ycen[istar]) + if (dy2 >= radius) + next + dr2 = (xcen[jstar] - xcen[istar]) ** 2 + dy2 ** 2 + if (dr2 >= radsq) + next + omit = true + break + } + + if (omit) + next + + nxstar = istar + return (ids[istar]) + } + + return (0) +end diff --git a/noao/digiphot/daophot/psf/dplocstar.x b/noao/digiphot/daophot/psf/dplocstar.x new file mode 100644 index 00000000..9925f44b --- /dev/null +++ b/noao/digiphot/daophot/psf/dplocstar.x @@ -0,0 +1,109 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +# DP_LOCSTAR -- Given an x,y position locate the star in the aperture +# photometry list which is within a critical radius of the input position. + +int procedure dp_locstar (dao, im, x, y) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +real x, y # input position + +int i +pointer psf, apsel +real crit_rad, fitrad, rad, ax, ay + +begin + # Set up some constants. + psf = DP_PSF(dao) + apsel = DP_APSEL(dao) + crit_rad = DP_MATCHRAD(dao) * DP_MATCHRAD(dao) + fitrad = DP_FITRAD(dao) + + # Search the list. + do i = 1, DP_APNUM(apsel) { + + # Compute distance of star from the cursur position. + ax = Memr[DP_APXCEN(apsel)+i-1] + ay = Memr[DP_APYCEN(apsel)+i-1] + if (! IS_INDEFR(ax) && ! IS_INDEFR(ay)) + rad = (ax - x) * (ax - x) + (ay - y) * (ay - y) + else + rad = MAX_REAL + + # Found the star. + if (rad <= crit_rad) { + + DP_CUR_PSF(psf) = i + DP_CUR_PSFID(psf) = Memi[DP_APID(apsel)+i-1] + DP_CUR_PSFX(psf) = ax + DP_CUR_PSFY(psf) = ay + DP_CUR_PSFSKY(psf) = Memr[DP_APMSKY(apsel)+i-1] + DP_CUR_PSFMAG(psf) = Memr[DP_APMAG(apsel)+i-1] + + # Is the star too close to the edge ? + if (int (ax - fitrad) < 0 || int (ax + fitrad) > IM_LEN(im,1) || + int (ay - fitrad) < 0 || int (ay + fitrad) > IM_LEN(im,2)) + return (-i) + else + return (i) + } + } + + return (0) +end + + +# DP_IDSTAR -- Given an id number locate the star in the aperture +# photometry list. + +int procedure dp_idstar (dao, im, idnum) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +int idnum # id number from photometry list + +int i +pointer psf, apsel +real fitrad, x, y + +begin + # Set up some constants. + psf = DP_PSF(dao) + apsel = DP_APSEL(dao) + fitrad = DP_FITRAD(dao) + + # Search the list. + do i = 1, DP_APNUM (apsel) { + + # Test the id number. + if (idnum != Memi[DP_APID(apsel)+i-1]) + next + + # Found the star. + x = Memr[DP_APXCEN(apsel)+i-1] + y = Memr[DP_APYCEN(apsel)+i-1] + DP_CUR_PSF(psf) = i + DP_CUR_PSFID(psf) = Memi[DP_APID(apsel)+i-1] + DP_CUR_PSFX(psf) = Memr[DP_APXCEN(apsel)+i-1] + DP_CUR_PSFY(psf) = Memr[DP_APYCEN(apsel)+i-1] + DP_CUR_PSFSKY(psf) = Memr[DP_APMSKY(apsel)+i-1] + DP_CUR_PSFMAG(psf) = Memr[DP_APMAG(apsel)+i-1] + + # Is the star too close to the edge ? + if (IS_INDEFR(x) || IS_INDEFR(y)) + return (0) + else if (int (x - fitrad) < 0 || int (x + fitrad) > + IM_LEN(im,1) || int (y - fitrad) < 0 || int (y + fitrad) > + IM_LEN(im,2)) + return (-i) + else + return (i) + } + + return (0) +end diff --git a/noao/digiphot/daophot/psf/dpmempsf.x b/noao/digiphot/daophot/psf/dpmempsf.x new file mode 100644 index 00000000..d198c778 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpmempsf.x @@ -0,0 +1,217 @@ +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + + +# DP_PSFSETUP -- Procedure to set up the PSF parameters. + +procedure dp_psfsetup (dp) + +pointer dp # pointer to daophot structure + +pointer psf + +begin + # Allocate memory for the psf fitting structure. + call malloc (DP_PSF(dp), LEN_PSFSTRUCT, TY_STRUCT) + psf = DP_PSF(dp) + + # Set the size of the data box t be extracted. + DP_LENUSERAREA(psf) = MIN_LENUSERAREA + + # Initialize the psf fitting structure. + DP_PC(psf) = NULL + DP_PV(psf) = NULL + DP_PTMP(psf) = NULL + DP_PZ(psf) = NULL + DP_PCLAMP(psf) = NULL + DP_POLD(psf) = NULL + DP_PSIGANA(psf) = 0.0 + DP_PSUMANA(psf) = 0.0 + + # Initialize the parameters for the current PSF star. + DP_CUR_PSF(psf) = 0 + DP_CUR_PSFID(psf) = 0 + DP_CUR_PSFX(psf) = INDEFR + DP_CUR_PSFY(psf) = INDEFR + DP_CUR_PSFSKY(psf) = INDEFR + DP_CUR_PSFMAG(psf) = INDEFR + + # Initialize the PSF star list arrays. + DP_PNUM (psf) = 0 + DP_PXCEN(psf) = NULL + DP_PYCEN(psf) = NULL + DP_PMAG(psf) = NULL + DP_PH(psf) = NULL + DP_PWEIGHT(psf) = NULL + DP_PSAT(psf) = NULL + DP_PXCLAMP(psf) = NULL + DP_PYCLAMP(psf) = NULL + DP_PXOLD(psf) = NULL + DP_PYOLD(psf) = NULL + + # Initialize the look-up table arrays. + DP_PSUMN(psf) = NULL + DP_PSUMW(psf) = NULL + DP_PSUMSQ(psf) = NULL + DP_PSIGMA(psf) = NULL + DP_PCONST(psf) = NULL + DP_POLDLUT(psf) = NULL + + # Allocate space for and initialize the plot sub-structure. + DP_PLOTTYPE(psf) = PSF_MESHPLOT + DP_MANGV(psf) = 30. + DP_MANGH(psf) = -30. + DP_MFLOOR(psf) = 0. + DP_MCEILING(psf) = 0. + DP_CFLOOR(psf) = 0.0 + DP_CCEILING(psf) = 0.0 +end + + +# DP_LMEMPSF -- Allocate memory required for fitting the list of psf stars. + +procedure dp_lmempsf (dao) + +pointer dao # pointer to daophot structure + +pointer psf + +begin + psf = DP_PSF(dao) + if (DP_PNUM(psf) > 0) { + call realloc (DP_PXCEN(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PYCEN(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PMAG(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PH(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PWEIGHT(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PSAT(psf), DP_PNUM(psf), TY_INT) + call realloc (DP_PXCLAMP(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PYCLAMP(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PXOLD(psf), DP_PNUM(psf), TY_REAL) + call realloc (DP_PYOLD(psf), DP_PNUM(psf), TY_REAL) + } +end + + +# DP_AMEMPSF -- Allocate the memory required for fitting the analytic +# PSF. + +procedure dp_amempsf (dao) + +pointer dao # pointer to daophot structure + +int npars +pointer psf, psffit + +begin + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + npars = DP_PSFNPARS(psffit) + + call realloc (DP_PC(psf), npars * npars, TY_REAL) + call realloc (DP_PV(psf), npars, TY_REAL) + call realloc (DP_PZ(psf), npars, TY_REAL) + call realloc (DP_PCLAMP(psf), npars, TY_REAL) + call realloc (DP_POLD(psf), npars, TY_REAL) + call realloc (DP_PTMP(psf), MAX_NFCTNPARS, TY_REAL) +end + + +# DP_TMEMPSF -- Allocate the memory required for fitting the lookup tables. + +procedure dp_tmempsf (dao) + +pointer dao # pointer to daophot structure + +int nexp +pointer psf, psffit + +begin + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + nexp = DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit) + DP_PSFSIZE(psffit) = 2 * (nint (2.0 * DP_PSFRAD(dao)) + 1) + 1 + DP_PSFRAD(dao) = (real (DP_PSFSIZE(psffit) - 1) / 2.0 - 1.0) / 2.0 + DP_SPSFRAD(dao) = DP_PSFRAD(dao) * DP_SCALE(dao) + DP_RPSFRAD(dao) = DP_PSFRAD(dao) * DP_SCALE(dao) + + call realloc (DP_PC(psf), nexp * nexp, TY_REAL) + call realloc (DP_PV(psf), nexp, TY_REAL) + call realloc (DP_PTMP(psf), nexp, TY_REAL) + + call realloc (DP_PSUMN(psf), DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit), + TY_REAL) + call realloc (DP_PSUMSQ(psf), DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit), + TY_REAL) + call realloc (DP_PSIGMA(psf), DP_PSFSIZE(psffit) * DP_PSFSIZE(psffit), + TY_REAL) + call realloc (DP_POLDLUT(psf), nexp * DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit), TY_REAL) + if (nexp > 1) + call realloc (DP_PCONST(psf), DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit), TY_REAL) + + call realloc (DP_PSFLUT(psffit), nexp * DP_PSFSIZE(psffit) * + DP_PSFSIZE(psffit), TY_REAL) +end + + +# DP_PSFCLOSE -- Procedure to set up the PSF parameters. + +procedure dp_psfclose (dp) + +pointer dp # pointer to daophot structure + +pointer psf + +begin + psf = DP_PSF(dp) + + if (DP_PC(psf) != NULL) + call mfree (DP_PC(psf), TY_REAL) + if (DP_PV(psf) != NULL) + call mfree (DP_PV(psf), TY_REAL) + if (DP_PTMP(psf) != NULL) + call mfree (DP_PTMP(psf), TY_REAL) + if (DP_PZ(psf) != NULL) + call mfree (DP_PZ(psf), TY_REAL) + if (DP_PCLAMP(psf) != NULL) + call mfree (DP_PCLAMP(psf), TY_REAL) + if (DP_POLD(psf) != NULL) + call mfree (DP_POLD(psf), TY_REAL) + + if (DP_PXCEN(psf) != NULL) + call mfree (DP_PXCEN(psf), TY_REAL) + if (DP_PYCEN(psf) != NULL) + call mfree (DP_PYCEN(psf), TY_REAL) + if (DP_PMAG(psf) != NULL) + call mfree (DP_PMAG(psf), TY_REAL) + if (DP_PH(psf) != NULL) + call mfree (DP_PH(psf), TY_REAL) + if (DP_PWEIGHT(psf) != NULL) + call mfree (DP_PWEIGHT(psf), TY_REAL) + if (DP_PSAT(psf) != NULL) + call mfree (DP_PSAT(psf), TY_INT) + if (DP_PXCLAMP(psf) != NULL) + call mfree (DP_PXCLAMP(psf), TY_REAL) + if (DP_PYCLAMP(psf) != NULL) + call mfree (DP_PYCLAMP(psf), TY_REAL) + if (DP_PXOLD(psf) != NULL) + call mfree (DP_PXOLD(psf), TY_REAL) + if (DP_PYOLD(psf) != NULL) + call mfree (DP_PYOLD(psf), TY_REAL) + + if (DP_PSUMN(psf) != NULL) + call mfree (DP_PSUMN(psf), TY_REAL) + if (DP_PSUMSQ(psf) != NULL) + call mfree (DP_PSUMSQ(psf), TY_REAL) + if (DP_PSIGMA(psf) != NULL) + call mfree (DP_PSIGMA(psf), TY_REAL) + if (DP_PCONST(psf) != NULL) + call mfree (DP_PCONST(psf), TY_REAL) + if (DP_POLDLUT(psf) != NULL) + call mfree (DP_POLDLUT(psf), TY_REAL) + + call mfree (psf, TY_STRUCT) +end diff --git a/noao/digiphot/daophot/psf/dpmkpsf.x b/noao/digiphot/daophot/psf/dpmkpsf.x new file mode 100644 index 00000000..2af7dbed --- /dev/null +++ b/noao/digiphot/daophot/psf/dpmkpsf.x @@ -0,0 +1,361 @@ +include <gset.h> +include <ctype.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +define HELPFILE "daophot$psf/mkpsf.key" + +# DP_MKPSF -- Construct a stellar PSF from one or more stars in an image frame. + +procedure dp_mkpsf (dao, im, psfim, opst, psfgr, gd, mgd, id, mkstars, + interactive, showplots) + +pointer dao # pointer to the main daophot structure +pointer im # pointer to the input image +pointer psfim # pointer to the output psfimage +int opst # the psf star list file descriptor +int psfgr # the psf group file descriptor +pointer gd # pointer to graphics descriptor +pointer mgd # pointer to the metacode file +pointer id # pointer to image display stream +bool mkstars # mark the added and deleted psf stars +bool interactive # interactive mode +bool showplots # show the plots + +real wx, wy +pointer apsel, sp, cmd, str +int i, key, wcs, idnum, istar, ip, npsf, verbose +bool psf_new, psf_computed, psf_written + +real dp_pstatr() +int clgcur(), dp_qverify(), dp_locstar(), dp_idstar(), dp_stati() +int dp_pstati(), ctoi(), dp_addstar(), dp_delstar(), dp_subpsf() +int dp_fitpsf() +bool dp_updatepsf() + +begin + # Get some pointers + apsel = DP_APSEL(dao) + + # Allocate some working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize some variables. + key = 'a' + if (dp_pstati (dao, PNUM) > 0) + psf_new = false + else + psf_new = true + psf_computed = false + psf_written = false + + # Begin to build the PSF. + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Convert coordinates if necessary. + call dp_vtol (im, wx, wy, wx, wy, 1) + + switch (key) { + + # Quit the interactive cursor loop. + case 'q': + + if (interactive) { + if (dp_qverify (dao, im, psfim, opst, psfgr, psf_new, + psf_computed, psf_written) == NO) + next + } else if (dp_updatepsf (dao, im, psfim, opst, psfgr, psf_new, + psf_computed, psf_written)) { + psf_computed = true + psf_written = true + } + + # Delete any empty psf and group files lying around. + if (! psf_written) + call dp_rmpsf (dao, psfim, opst, psfgr) + + # Delete the PSF stars. + call dp_pseti (dao, PNUM, 0) + + break + + # Print the help page. + case '?': + if (id == NULL) + call pagefile (HELPFILE, "") + else + call gpagefile (id, HELPFILE, "") + + # Add the star nearest the cursor position to the PSF. + case 'a': + if (dp_addstar (dao, im, wx, wy, INDEFR, 0, gd, mgd, + showplots) == ERR) + next + psf_new = false + psf_computed = false + psf_written = false + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), dp_pstatr (dao, + CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # Subtract the star nearest the cursor position from the PSF + # using the current best fit. + case 's': + if (psf_new) { + call printf ("The PSF star list is empty\n") + next + } + if (! psf_computed) { + call printf ("The PSF is not uptodate\n") + next + } + if (dp_subpsf (dao, im, wx, wy, 0, gd, mgd, showplots) == ERR) + next + if (dp_pstati (dao, PNUM) <= 0) + psf_new = true + psf_computed = false + psf_written = false + + # Delete the star nearest the cursor position from the PSF. + case 'd': + if (dp_delstar (dao, im, wx, wy, 0, gd, showplots) == ERR) + next + if (dp_pstati (dao, PNUM) <= 0) + psf_new = true + psf_computed = false + psf_written = false + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), dp_pstatr (dao, + CUR_PSFY), GM_CROSS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # List all the current psf stars. + case 'l': + if (interactive) + call dp_listpsf (dao, im) + + # Fit the PSF. + case 'f': + + # Reread the stars is necessary. + if (dp_pstati (dao, PNUM) > 0 && (psf_new || psf_computed)) { + npsf = dp_pstati (dao, PNUM) + call dp_pseti (dao, PNUM, 0) + call dp_reinit (dao) + verbose = dp_stati (dao, VERBOSE) + call dp_seti (dao, VERBOSE, NO) + do i = 1, npsf + if (dp_addstar (dao, im, wx, wy, INDEFR, + Memi[DP_APID(apsel)+i-1], gd, mgd, + false) == OK) + psf_new = false + call dp_seti (dao, VERBOSE, verbose) + } + + # Fit the psf. + if (psf_new) { + call printf ("The PSF star list is empty\n") + } else if (dp_fitpsf (dao, im, Memc[str], SZ_LINE) == OK) { + psf_computed = true + } else { + call printf ("%s\n") + call pargstr (Memc[str]) + } + + # Review the fit. + case 'r': + if (psf_new) { + call printf ("The PSF star list is empty\n") + } else if (! psf_computed) { + call printf ("The PSF is not uptodate\n") + } else { + i = 1 + repeat { + if (dp_subpsf (dao, im, wx, wy, + Memi[DP_APID(apsel)+i-1], gd, mgd, + showplots) == OK) { + if (dp_pstati (dao, PNUM) <= 0) + psf_new = true + psf_computed = false + psf_written = false + } else + i = i + 1 + } until (i > dp_pstati(dao, PNUM)) + } + + # Rebuild the PSF from scratch. + case 'z': + + # Print message. + if (interactive) { + call dp_stats (dao, OUTPHOTFILE, Memc[str], SZ_FNAME) + call printf ( + "PSF star list, image, and output files deleted\n") + } + + # Delete the previous psf an dgroup files if any. + call dp_rmpsf (dao, psfim, opst, psfgr) + + # Delete the PSF stars. + call dp_pseti (dao, PNUM, 0) + + # Reopen the psf image and group files. + call dp_oppsf (dao, psfim, opst, psfgr) + + # Reset the reduction flags. + psf_new = true + psf_computed = false + psf_written = false + + # Write out the PSF. + case 'w': + if (dp_updatepsf (dao, im, psfim, opst, psfgr, psf_new, + psf_computed, psf_written)) { + psf_computed = true + psf_written = true + } + + # Locate the star in the aperture photometry file and print out + # the photometry. + case 'p': + if (interactive) + istar = dp_locstar (dao, im, wx, wy) + else + next + if (istar > 0) + call dp_pshow (dao, im, istar) + else if (istar == 0) + call printf ("Star not found in the photometry file.\n") + else + call printf ( + "Star off or too near the edge of the image.\n") + + # Command mode + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + switch (Memc[cmd+ip-1]) { + + case 'p': + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call dp_pcolon (dao, psfim, opst, psfgr, Memc[cmd], + psf_new, psf_computed, psf_written) + next + } + if (! interactive) + next + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + istar = dp_locstar (dao, im, wx, wy) + else + istar = dp_idstar (dao, im, idnum) + if (istar > 0) + call dp_pshow (dao, im, istar) + else if (istar == 0) + call printf ( + "Star not found in the photometry file\n") + else + call printf ( + "Star is off or too near the edge of the image.\n") + + case 'a': + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call dp_pcolon (dao, psfim, opst, psfgr, Memc[cmd], + psf_new, psf_computed, psf_written) + next + } + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + idnum = 0 + if (dp_addstar (dao, im, wx, wy, INDEFR, idnum, gd, mgd, + showplots) == ERR) + next + psf_new = false + psf_computed = false + psf_written = false + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), + dp_pstatr (dao, CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + case 's': + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call dp_pcolon (dao, psfim, opst, psfgr, Memc[cmd], + psf_new, psf_computed, psf_written) + next + } + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + idnum = 0 + if (! psf_computed) { + call printf ("The PSF has not been fit\n") + next + } + if (! psf_computed) { + call printf ("Warning: The PSF is not uptodate\n") + next + } + if (dp_subpsf (dao, im, wx, wy, idnum, gd, mgd, + showplots) == ERR) + next + if (dp_pstati (dao, PNUM) <= 0) + psf_new = true + psf_computed = false + psf_written = false + + case 'd': + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call dp_pcolon (dao, psfim, opst, psfgr, Memc[cmd], + psf_new, psf_computed, psf_written) + next + } + ip = ip + 1 + if (ctoi (Memc[cmd], ip, idnum) <= 0) + idnum = 0 + if (dp_delstar (dao, im, wx, wy, idnum, gd, + showplots) == ERR) + next + if (dp_pstati (dao, PNUM) <= 0) + psf_new = true + psf_computed = false + psf_written = false + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr (dao, CUR_PSFX), + dp_pstatr (dao, CUR_PSFY), GM_CROSS, -5.0, -5.0) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + default: + call dp_pcolon (dao, psfim, opst, psfgr, Memc[cmd], + psf_new, psf_computed, psf_written) + } + + default: + call printf ("Unknown keystroke command\n") + } + } + + # Free up memory. + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dppcolon.x b/noao/digiphot/daophot/psf/dppcolon.x new file mode 100644 index 00000000..9934f244 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppcolon.x @@ -0,0 +1,271 @@ +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_PCOLON -- Process the PSF fitting colon commands. + +procedure dp_pcolon (dao, psfim, opst, psfgr, cmdstr, psf_new, psf_computed, + psf_written) + +pointer dao # pointer to the daophot structure +pointer psfim # pointer to the output psf image +int opst # the output psf star list file descriptor +int psfgr # the output psf group file descriptor +char cmdstr[ARB] # the input command string +bool psf_new # is the psf star list defined ? +bool psf_computed # has the psf been fit ? +bool psf_written # has the psf been updated ? + +bool bval +int ncmd, ival +pointer sp, cmd, str1, str2, str3 +real rval + +bool itob() +int strdic(), nscan(), btoi(), dp_stati(), open(), dp_fctdecode() +int dp_pstati(), strlen() +pointer immap(), tbtopn() +real dp_statr() +errchk immap(), open(), tbtopn() + +begin + # Allocate working space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + call salloc (str3, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PSF_CMDS) + switch (ncmd) { + + case PSFCMD_PSFIMAGE: + call gargwrd (Memc[str1], SZ_LINE) + call gargwrd (Memc[str2], SZ_LINE) + call gargwrd (Memc[str3], SZ_LINE) + if (nscan() == 1) { + call dp_stats (dao, PSFIMAGE, Memc[str1], SZ_LINE) + call dp_stats (dao, OUTREJFILE, Memc[str2], SZ_LINE) + call dp_stats (dao, OUTPHOTFILE, Memc[str3], SZ_LINE) + call printf ("psfim = %s opstfile = %s grpfile = %s\n") + call pargstr (Memc[str1]) + call pargstr (Memc[str2]) + call pargstr (Memc[str3]) + } else if (nscan() == 4) { + if (psfim != NULL) + call imunmap (psfim) + psfim = NULL + if (opst != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (opst) + else + call tbtclo (opst) + } + opst = NULL + if (psfgr != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (psfgr) + else + call tbtclo (psfgr) + } + psfgr = NULL + iferr { + psfim = immap (Memc[str1], NEW_IMAGE, dp_pstati (dao, + LENUSERAREA)) + if (dp_stati (dao, TEXT) == YES) + opst = open (Memc[str2], NEW_FILE, TEXT_FILE) + else + opst = tbtopn (Memc[str2], NEW_FILE, 0) + if (dp_stati (dao, TEXT) == YES) + psfgr = open (Memc[str3], NEW_FILE, TEXT_FILE) + else + psfgr = tbtopn (Memc[str3], NEW_FILE, 0) + } then { + if (psfim != NULL) + call imunmap (psfim) + psfim = NULL + if (opst != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (opst) + else + call tbtclo (opst) + } + opst = NULL + if (psfgr != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (psfgr) + else + call tbtclo (psfgr) + } + psfgr = NULL + call dp_sets (dao, PSFIMAGE, "") + call dp_sets (dao, OUTREJFILE, "") + call dp_sets (dao, OUTPHOTFILE, "") + } else { + call dp_sets (dao, PSFIMAGE, Memc[str1]) + call dp_sets (dao, OUTREJFILE, Memc[str2]) + call dp_sets (dao, OUTPHOTFILE, Memc[str2]) + } + psf_written = false + } + + case PSFCMD_FUNCTION: + call gargwrd (Memc[str1], SZ_LINE) + if (nscan() == 1) { + call dp_stats (dao, FUNCLIST, Memc[str1], SZ_LINE) + call strcpy (Memc[str1+1], Memc[str2], + strlen (Memc[str1]) - 2) + call printf ("function = %s\n") + call pargstr (Memc[str2]) + } else if (dp_fctdecode (Memc[str1], Memc[str2], SZ_FNAME) > 0) { + call dp_sets (dao, FUNCLIST, Memc[str2]) + psf_computed = false + psf_written = false + } + + case PSFCMD_VARORDER: + call gargi (ival) + if (nscan() == 1) { + call printf ("varorder = %d\n") + call pargi (dp_stati (dao, VARORDER)) + } else if (ival >= -1 && ival <= 2) { + call dp_seti (dao, VARORDER, ival) + psf_computed = false + psf_written = false + } + + case PSFCMD_FEXPAND: + call gargb (bval) + if (nscan() == 1) { + call printf ("fexpand = %b\n") + call pargb (itob (dp_stati (dao, FEXPAND))) + } else { + call dp_seti (dao, FEXPAND, btoi (bval)) + psf_computed = false + psf_written = false + } + + case PSFCMD_PSFRAD: + call gargr (rval) + if (nscan() == 1) { + call printf ("psfrad = %g scale units\n") + call pargr (dp_statr (dao, SPSFRAD)) + } else { + call dp_setr (dao, SPSFRAD, rval) + call dp_setr (dao, RPSFRAD, rval) + call dp_setr (dao, PSFRAD, rval / dp_statr (dao, SCALE)) + psf_computed = false + psf_written = false + } + + case PSFCMD_FITRAD: + call gargr (rval) + if (nscan() == 1) { + call printf ("fitrad = %g scale units\n") + call pargr (dp_statr (dao, SFITRAD)) + } else { + call dp_setr (dao, SFITRAD, rval) + call dp_setr (dao, FITRAD, rval / dp_statr (dao, SCALE)) + psf_new = true + psf_computed = false + psf_written = false + } + + case PSFCMD_MATCHRAD: + call gargr (rval) + if (nscan() == 1) { + call printf ("matchrad = %g scale units\n") + call pargr (dp_statr (dao, SMATCHRAD)) + } else { + call dp_setr (dao, SMATCHRAD, rval) + call dp_setr (dao, MATCHRAD, rval / dp_statr (dao, SCALE)) + } + + case PSFCMD_NCLEAN: + call gargi (ival) + if (nscan() == 1) { + call printf ("nclean = %d\n") + call pargi (dp_stati (dao, NCLEAN)) + } else if (ival >= 0) { + call dp_seti (dao, NCLEAN, ival) + psf_computed = false + psf_written = false + } + + case PSFCMD_SATURATED: + call gargb (bval) + if (nscan() == 1) { + call printf ("saturated = %b\n") + call pargb (itob (dp_stati (dao, SATURATED))) + } else { + call dp_seti (dao, SATURATED, btoi (bval)) + psf_computed = false + psf_written = false + } + + case PSFCMD_SCALE: + call gargr (rval) + if (nscan() == 1) { + call printf ("scale = %g units per pixel\n") + call pargr (dp_statr (dao, SCALE)) + } else { + call dp_setr (dao, FWHMPSF, dp_statr (dao, SFWHMPSF) / rval) + call dp_setr (dao, PSFRAD, dp_statr (dao, SPSFRAD) / rval) + call dp_setr (dao, FITRAD, dp_statr (dao, SFITRAD) / rval) + call dp_setr (dao, MATCHRAD, dp_statr (dao, SMATCHRAD) / rval) + call dp_setr (dao, SCALE, rval) + psf_new = true + psf_computed = false + psf_written = false + } + + case PSFCMD_FWHMPSF: + call gargr (rval) + if (nscan() == 1) { + call printf ("fwhmpsf = %g scale units\n") + call pargr (dp_statr (dao, SFWHMPSF)) + } else { + call dp_setr (dao, SFWHMPSF, rval) + call dp_setr (dao, FWHMPSF, rval / dp_statr (dao, SCALE)) + psf_computed = false + psf_written = false + } + + case PSFCMD_DATAMIN: + call gargr (rval) + if (nscan() == 1) { + call printf ("datamin = %g counts\n") + call pargr (dp_statr (dao, MINGDATA)) + } else { + call dp_setr (dao, MINGDATA, rval) + psf_new = true + psf_computed = false + psf_written = false + } + + case PSFCMD_DATAMAX: + call gargr (rval) + if (nscan() == 1) { + call printf ("datamax = %g counts\n") + call pargr (dp_statr (dao, MAXGDATA)) + } else { + call dp_setr (dao, MAXGDATA, rval) + psf_new = true + psf_computed = false + psf_written = false + } + + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dppconfirm.x b/noao/digiphot/daophot/psf/dppconfirm.x new file mode 100644 index 00000000..092cfec3 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppconfirm.x @@ -0,0 +1,26 @@ +# DP_PCONFIRM -- Procedure to confirm the critical psf parameters. + +procedure dp_pconfirm (dao) + +pointer dao # pointer to the daophot structure + +begin + call printf ("\n") + + # Verify the functional form of the psf. + call dp_vfunction(dao) + call dp_vvarorder (dao) + #call dp_vfexpand (dao) + + # Confirm the psf radius. + call dp_vpsfrad (dao) + + # Confirm the fitting radius. + call dp_vfitrad (dao) + + # Confirm the minimum and maximum good data values. + call dp_vdatamin (dao) + call dp_vdatamax (dao) + + call printf ("\n") +end diff --git a/noao/digiphot/daophot/psf/dpplotpsf.x b/noao/digiphot/daophot/psf/dpplotpsf.x new file mode 100644 index 00000000..4b80c50f --- /dev/null +++ b/noao/digiphot/daophot/psf/dpplotpsf.x @@ -0,0 +1,49 @@ +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_PLOTPSF -- Plot the psf using the default plot type. + +procedure dp_plotpsf (dao, im, subrast, ncols, nlines, x1, y1, gd) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +real subrast[ncols,nlines] # image subraster +int ncols, nlines # dimensions of the subraster +int x1, y1 # coordinates of the lower left corner +pointer gd # pointer to the graphics stream + +real tx, ty +pointer sp, title +real dp_pstatr() +int dp_pstati() + +begin + # Return if the graphics stream is undefined. + if (gd == NULL) + return + + # Comvert the coordinates if necessary. + call dp_wout (dao, im, dp_pstatr(dao, CUR_PSFX), dp_pstatr(dao, + CUR_PSFY), tx, ty, 1) + + # Construct the title. + call smark (sp) + call salloc (title, SZ_LINE, TY_CHAR) + call sprintf (Memc[title], SZ_LINE, "Star: %d X: %g Y: %g Mag: %g") + call pargi (dp_pstati (dao, CUR_PSFID)) + call pargr (tx) + call pargr (ty) + call pargr (dp_pstatr (dao, CUR_PSFMAG)) + + # Initialize plot. + if (dp_pstati (dao, PLOTTYPE) == PSF_MESHPLOT) + call dp_surfpsf (dao, subrast, ncols, nlines, Memc[title], gd) + else if (dp_pstati (dao, PLOTTYPE) == PSF_CONTOURPLOT) + call dp_contpsf (dao, subrast, ncols, nlines, Memc[title], gd) + else if (dp_pstati (dao, PLOTTYPE) == PSF_RADIALPLOT) + call dp_radpsf (dao, subrast, ncols, nlines, x1, y1, + Memc[title], gd) + + call gdeactivate (gd, 0) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dppset.x b/noao/digiphot/daophot/psf/dppset.x new file mode 100644 index 00000000..d98b22bd --- /dev/null +++ b/noao/digiphot/daophot/psf/dppset.x @@ -0,0 +1,81 @@ +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_PSETS -- Set a psf fitting string parameter. +# +#procedure dp_psets (dao, param, str) + +#pointer dao # pointer to daophot structure +#int param # parameter +#char str[ARB] # string value + +#begin +# switch (param) { +# default: +# call error (0, "DP_PSETS: Unknown psf fitting string parameter") +# } +#end + + +# DP_PSETI -- Set a daophot psf fitting integer parameter. + +procedure dp_pseti (dao, param, ival) + +pointer dao # pointer to daophot structure +int param # parameter +int ival # integer value + +pointer psf + +begin + psf = DP_PSF(dao) + + switch (param) { + case CUR_PSF: + DP_CUR_PSF(psf) = ival + case CUR_PSFID: + DP_CUR_PSFID(psf) = ival + case PNUM: + DP_PNUM(psf) = ival + case PLOTTYPE: + DP_PLOTTYPE(psf) = ival + case LENUSERAREA: + DP_LENUSERAREA(psf) = ival + default: + call error (0, "DP_PSETI: Unknown integer psf fitting parameter") + } +end + + +# DP_PSETR -- Set a real psf fitting parameter. + +procedure dp_psetr (dao, param, rval) + +pointer dao # pointer to daophot structure +int param # parameter +real rval # real value + +pointer psf + +begin + psf = DP_PSF(dao) + + switch (param) { + case CUR_PSFX: + DP_CUR_PSFX(psf) = rval + case CUR_PSFY: + DP_CUR_PSFY(psf) = rval + case CUR_PSFSKY: + DP_CUR_PSFSKY(psf) = rval + case CUR_PSFMAG: + DP_CUR_PSFMAG(psf) = rval + case CUR_PSFMIN: + DP_CUR_PSFMIN(psf) = rval + case CUR_PSFMAX: + DP_CUR_PSFMAX(psf) = rval + case CUR_PSFGMAX: + DP_CUR_PSFGMAX(psf) = rval + default: + call error (0, "DP_SETR: Unknown real psf fitting parameter") + } +end diff --git a/noao/digiphot/daophot/psf/dppsfutil.x b/noao/digiphot/daophot/psf/dppsfutil.x new file mode 100644 index 00000000..17df56f6 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppsfutil.x @@ -0,0 +1,381 @@ +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +# DP_PFSWAP -- Swap two stars in the PSF star list. + +procedure dp_pfswap (dao, star1, star2) + +pointer dao # pointer to the daophot structure +int star1 # index of first star to be swapped +int star2 # index of second star to be swapped + +pointer apsel, psf + +begin + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + call dp_10swap (star1, star2, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], Memr[DP_PXCEN(psf)], + Memr[DP_PYCEN(psf)], Memr[DP_PH(psf)], Memr[DP_PMAG(psf)], + Memi[DP_PSAT(psf)]) +end + + +# DP_PFREORDER -- Move a star in the PSF star list to the end of the +# list. + +procedure dp_pfreorder (dao, star, nstars) + +pointer dao # pointer to the daophot structure +int star # index of star to be deleted +int nstars # index of second star to be swapped + +pointer apsel, psf + +begin + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + call dp_10reorder (star, Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], Memr[DP_PXCEN(psf)], Memr[DP_PYCEN(psf)], + Memr[DP_PH(psf)], Memr[DP_PMAG(psf)], Memi[DP_PSAT(psf)], nstars) +end + + +# DP_APLSWAP -- Swap two stars in the daophot photometry substructure. + +procedure dp_aplswap (dao, star1, star2) + +pointer dao # pointer to the daophot structure +int star1 # index of first star to be swapped +int star2 # index of second star to be swapped + +pointer apsel + +begin + apsel = DP_APSEL(dao) + call dp_5swap (star1, star2, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)]) +end + + +# DP_XYHPSF -- Compute the initial x, y, and height of the current PSF star. + +procedure dp_xyhpsf (dao, star, mag, saturated) + +pointer dao # pointer to the daophot structure +int star # star for which x, y, h is to be computed +real mag # magnitude of proposed psf star +int saturated # is the star saturated + +pointer apsel, psf, psffit +real dhdxc, dhdyc, junk +real dp_profile() + +begin + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + psffit = DP_PSFFIT(dao) + + Memr[DP_PXCEN(psf)+star-1] = Memr[DP_APXCEN(apsel)+star-1] + Memr[DP_PYCEN(psf)+star-1] = Memr[DP_APYCEN(apsel)+star-1] + if (saturated == YES) + Memr[DP_PH(psf)+star-1] = INDEFR + else + Memr[DP_PH(psf)+star-1] = (DP_CUR_PSFGMAX(psf) - + Memr[DP_APMSKY(apsel)+star-1]) / + dp_profile (DP_PSFUNCTION(psffit), + 0.0, 0.0, Memr[DP_PSFPARS(psffit)], dhdxc, dhdyc, junk, 0) + Memr[DP_PMAG(psf)+star-1] = mag + Memi[DP_PSAT(psf)+star-1] = saturated +end + + +# DP_LISTPSF -- List the PSF stars. + +procedure dp_listpsf (dao, im) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor + +real x, y +pointer apsel, psf +int i + +begin + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + + call printf ("\nCurrent PSF star list\n") + do i = 1, DP_PNUM(psf) { + call dp_ltov (im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], x, y, 1) + call printf ( + " Star: %4d X: %7.2f Y: %7.2f Mag: %7.2f Sky: %10.1f\n") + call pargi (Memi[DP_APID(apsel)+i-1]) + call pargr (x) + call pargr (y) + call pargr (Memr[DP_APMAG(apsel)+i-1]) + call pargr (Memr[DP_APMSKY(apsel)+i-1]) + } + call printf ("\n") +end + + +# DP_LISTPARS -- List the analytic PSF parameters. + +procedure dp_listpars (dao) + +pointer dao # pointer to the daophot structure + +pointer psffit + +begin + psffit = DP_PSFFIT(dao) + call printf ("\nAnalytic PSF fit \n") + call printf ( + " Function: %s X: %g Y: %g Height: %g Psfmag: %g\n") + call pargstr (DP_FUNCTION(dao)) + call pargr (DP_PSFX(psffit)) + call pargr (DP_PSFY(psffit)) + call pargr (DP_PSFHEIGHT(psffit)) + call pargr (DP_PSFMAG(psffit)) + + switch (DP_PSFUNCTION(psffit)) { + case FCTN_GAUSS: + call printf (" Par1: %g Par2: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + case FCTN_MOFFAT15: + call printf ( + " Par1: %g Par2: %g XYterm: %g Moffat: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + call pargr (Memr[DP_PSFPARS(psffit)+2]) + call pargr (1.5) + case FCTN_PENNY1: + call printf (" Par1: %g Par2: %g Par3: %g Par4: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + call pargr (Memr[DP_PSFPARS(psffit)+2]) + call pargr (Memr[DP_PSFPARS(psffit)+3]) + case FCTN_MOFFAT25: + call printf ( + " Par1: %g Par2: %g XYterm: %g Moffat: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + call pargr (Memr[DP_PSFPARS(psffit)+2]) + call pargr (2.5) + case FCTN_PENNY2: + call printf ( + " Par1: %g Par2: %g Par3: %g Par4: %g Par5: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + call pargr (Memr[DP_PSFPARS(psffit)+2]) + call pargr (Memr[DP_PSFPARS(psffit)+3]) + call pargr (Memr[DP_PSFPARS(psffit)+4]) + case FCTN_LORENTZ: + call printf (" Par1: %g Par2: %g Par3: %g\n") + call pargr (Memr[DP_PSFPARS(psffit)]) + call pargr (Memr[DP_PSFPARS(psffit)+1]) + call pargr (Memr[DP_PSFPARS(psffit)+2]) + } +end + + +# DP_PSHOW -- Print photometry for the given star + +procedure dp_pshow (dao, im, istar) + +pointer dao # pointer to the main daophot descriptor +pointer im # the input image descriptor +int istar # star to be printed + +real x, y +pointer apsel + +begin + apsel = DP_APSEL(dao) + call dp_ltov (im, Memr[DP_APXCEN(apsel)+istar-1], + Memr[DP_APYCEN(apsel)+istar-1], x, y, 1) + call printf ( + "Star: %4d X: %7.2f Y: %7.2f Mag: %7.2f Sky: %10.1f\n") + call pargi (Memi[DP_APID(apsel)+istar-1]) + call pargr (x) + call pargr (y) + call pargr (Memr[DP_APMAG(apsel)+istar-1]) + call pargr (Memr[DP_APMSKY(apsel)+istar-1]) +end + + +# DP_10REORDER -- Move a PSF star to the end of the list. + +procedure dp_10reorder (star, id, x, y, mag, sky, xfit, yfit, hfit, pmag, + sat, nstars) + +int star # star to be moved to the end of the list +int id[ARB] # the ids of the stars +real x[ARB] # the x positions of the stars +real y[ARB] # the y positions of the stars +real mag[ARB] # the magnitudes of the stars +real sky[ARB] # the sky values of the stars +real xfit[ARB] # the current x fit array +real yfit[ARB] # the current y fit array +real hfit[ARB] # the current height of the stars +real pmag[ARB] # the psf star list magnitude +int sat[ARB] # are the star saturated +int nstars # number of stars in the list + +int i, ihold, sfhold +real xhold, yhold, mhold, shold, xfhold, yfhold, hfhold, mfhold + +begin + ihold = id[star] + xhold = x[star] + yhold = y[star] + mhold = mag[star] + shold = sky[star] + xfhold = xfit[star] + yfhold = yfit[star] + hfhold = hfit[star] + mfhold = pmag[star] + sfhold = sat[star] + + do i = star + 1, nstars { + id[i-1] = id[i] + x[i-1] = x[i] + y[i-1] = y[i] + mag[i-1] = mag[i] + sky[i-1] = sky[i] + xfit[i-1] = xfit[i] + yfit[i-1] = yfit[i] + hfit[i-1] = hfit[i] + pmag[i-1] = pmag[i] + sat[i-1] = sat[i] + } + + id[nstars] = ihold + x[nstars] = xhold + y[nstars] = yhold + mag[nstars] = mhold + sky[nstars] = shold + xfit[nstars] = xfhold + yfit[nstars] = yfhold + hfit[nstars] = hfhold + pmag[nstars] = mfhold + sat[nstars] = sfhold +end + + +# DP_5SWAP -- Exchange the position of two stars in the APPHOT photometry +# results. + +procedure dp_5swap (star1, star2, id, x, y, mag, sky) + +int star1, star2 # the indices of the two stars to exchange +int id[ARB] # the ids of the stars +real x[ARB] # the x postions of the stars +real y[ARB] # the y positions of the stars +real mag[ARB] # the magnitudes of the stars +real sky[ARB] # the sky values of the stars + +int ihold +real xhold, yhold, mhold, shold + +begin + ihold = id[star1] + xhold = x[star1] + yhold = y[star1] + mhold = mag[star1] + shold = sky[star1] + + id[star1] = id[star2] + x[star1] = x[star2] + y[star1] = y[star2] + mag[star1] = mag[star2] + sky[star1] = sky[star2] + + id[star2] = ihold + x[star2] = xhold + y[star2] = yhold + mag[star2] = mhold + sky[star2] = shold +end + + +# DP_10SWAP -- Exchange the position of two stars in the APPHOT photometry +# and PSF fitting results. + +procedure dp_10swap (star1, star2, id, x, y, mag, sky, xfit, yfit, hfit, pmag, + sat) + +int star1, star2 # the indices of the two stars to exchange +int id[ARB] # the ids of the stars +real x[ARB] # the x postions of the stars +real y[ARB] # the y positions of the stars +real mag[ARB] # the magnitudes of the stars +real sky[ARB] # the sky values of the stars +real xfit[ARB] # the current x fit array +real yfit[ARB] # the current y fit array +real hfit[ARB] # the current height of the stars +real pmag[ARB] # the psf star list magnitudes +int sat[ARB] # are the star saturated + +int ihold, sfhold +real xhold, yhold, mhold, shold, xfhold, yfhold, hfhold, mfhold + +begin + ihold = id[star1] + xhold = x[star1] + yhold = y[star1] + mhold = mag[star1] + shold = sky[star1] + xfhold = xfit[star1] + yfhold = yfit[star1] + hfhold = hfit[star1] + mfhold = pmag[star1] + sfhold = sat[star1] + + id[star1] = id[star2] + x[star1] = x[star2] + y[star1] = y[star2] + mag[star1] = mag[star2] + sky[star1] = sky[star2] + xfit[star1] = xfit[star2] + yfit[star1] = yfit[star2] + hfit[star1] = hfit[star2] + pmag[star1] = pmag[star2] + sat[star1] = sat[star2] + + id[star2] = ihold + x[star2] = xhold + y[star2] = yhold + mag[star2] = mhold + sky[star2] = shold + xfit[star2] = xfhold + yfit[star2] = yfhold + hfit[star2] = hfhold + pmag[star2] = mfhold + sat[star2] = sfhold +end + + +# DP_ISSAT -- Is the candidate PSF star saturated ? + +int procedure dp_issat (dao, starno) + +pointer dao # pointer to the daophot structure +int starno # the star index number + +pointer psf + +begin + psf = DP_PSF(dao) + if (Memi[DP_PSAT(psf)+starno-1] == YES) + return (YES) + else + return (NO) +end diff --git a/noao/digiphot/daophot/psf/dppstat.x b/noao/digiphot/daophot/psf/dppstat.x new file mode 100644 index 00000000..4b3f903f --- /dev/null +++ b/noao/digiphot/daophot/psf/dppstat.x @@ -0,0 +1,80 @@ +include "../lib/daophotdef.h" +include "../lib/psfdef.h" +# +## DP_PSTATS -- Fetch a psf fitting string parameter. +# +#procedure dp_pstats (dao, param, str, maxch) +# +#pointer dao # pointer to daophot structure +#int param # parameter +#char str[ARB] # string value +#int maxch # maximum number of characters +# +#begin +# switch (param) { +# default: +# call error (0, "DP_PSTATS: Unknown psf fitting string parameter") +# } +#end +# + +# DP_PSTATI -- Fetch a psf fitting integer parameter. + +int procedure dp_pstati (dao, param) + +pointer dao # pointer to daophot structure +int param # parameter + +pointer psf + +begin + psf = DP_PSF(dao) + + switch (param) { + case CUR_PSF: + return (DP_CUR_PSF(psf)) + case CUR_PSFID: + return (DP_CUR_PSFID(psf)) + case PNUM: + return (DP_PNUM(psf)) + case PLOTTYPE: + return (DP_PLOTTYPE(psf)) + case LENUSERAREA: + return (DP_LENUSERAREA(psf)) + default: + call error (0, "DP_PSTATI: Unknown psf fitting integer parameter") + } +end + + +# DP_PSTATR -- Fetch a psf fitting real parameter. + +real procedure dp_pstatr (dao, param) + +pointer dao # pointer to daophot structure +int param # parameter + +pointer psf + +begin + psf = DP_PSF(dao) + + switch (param) { + case CUR_PSFX: + return (DP_CUR_PSFX(psf)) + case CUR_PSFY: + return (DP_CUR_PSFY(psf)) + case CUR_PSFSKY: + return (DP_CUR_PSFSKY(psf)) + case CUR_PSFMAG: + return (DP_CUR_PSFMAG(psf)) + case CUR_PSFMIN: + return (DP_CUR_PSFMIN(psf)) + case CUR_PSFMAX: + return (DP_CUR_PSFMAX(psf)) + case CUR_PSFGMAX: + return (DP_CUR_PSFGMAX(psf)) + default: + call error (0, "DP_PSTATR: Unknown psf fitting parameter") + } +end diff --git a/noao/digiphot/daophot/psf/dppsubrast.x b/noao/digiphot/daophot/psf/dppsubrast.x new file mode 100644 index 00000000..1fafb408 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppsubrast.x @@ -0,0 +1,172 @@ +include <mach.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +# DP_PSUBRAST -- Fetch the prospective PSF star data and check it for bad +# pixel values. + +pointer procedure dp_psubrast (dao, im, lowbad, highbad, x1, x2, y1, y2, + saturated) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +real lowbad, highbad # minimum and maximum good data values +int x1, x2 # output x limits of the extracted subraster +int y1, y2 # output y limits of the extracted subraster +int saturated # is the star saturated ? + +pointer psf, buf +real psfrad, fitrad +int dp_chksr() +pointer dp_subrast() + +begin + # Initialize. + psf = DP_PSF(dao) + buf = NULL + psfrad = DP_PSFRAD(dao) + fitrad = DP_FITRAD(dao) + + # Get the data. + buf = dp_subrast (im, DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), psfrad, + x1, x2, y1, y2) + if (buf == NULL) + return (NULL) + + # Check for bad pixels in subraster, compute the min and max. + if (dp_chksr (DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), Memr[buf], + x2 - x1 + 1, y2 - y1 + 1, x1, y1, psfrad, fitrad, + lowbad, highbad, saturated, DP_CUR_PSFMIN(psf), + DP_CUR_PSFMAX(psf), DP_CUR_PSFGMAX(psf)) == ERR) { + call mfree (buf, TY_REAL) + return (NULL) + } + + return (buf) +end + + +# DL_LSUBRAST -- Give a valid PSF star compute the limits of the data to +# be extracted around it. + +int procedure dp_lsubrast (im, xcen, ycen, radius, x1, x2, y1, y2) + +pointer im # input image descriptor +real xcen, ycen # center of subraster +real radius # radius of the box +int x1, y1, x2, y2 # boundaries of subraster + +begin + # Calculate start position of extraction box. + x1 = int (xcen - radius) - 2 + x2 = int (xcen + radius) + 3 + y1 = int (ycen - radius) - 2 + y2 = int (ycen + radius) + 3 + if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1) + return (ERR) + + x1 = max (1, x1) + x2 = min (IM_LEN(im,1), x2) + y1 = max (1, y1) + y2 = min (IM_LEN(im,2), y2) + return (OK) +end + + +# DP_SUBRAST -- Given a valid PSF star extract the data around it. + +pointer procedure dp_subrast (im, xcen, ycen, radius, x1, x2, y1, y2) + +pointer im # input image descriptor +real xcen, ycen # center of subraster +real radius # radius of the box +int x1, y1, x2, y2 # boundaries of subraster + +int j, ncols +pointer buf, ptr, imbuf +pointer imgs2r() + +begin + # Calculate start position of extraction box. + x1 = int (xcen - radius) - 2 + x2 = int (xcen + radius) + 3 + y1 = int (ycen - radius) - 2 + y2 = int (ycen + radius) + 3 + if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1) + return (NULL) + + x1 = max (1, x1) + x2 = min (IM_LEN(im,1), x2) + y1 = max (1, y1) + y2 = min (IM_LEN(im,2), y2) + call malloc (buf, (x2 - x1 + 1) * (y2 - y1 + 1), TY_REAL) + + ptr = buf + ncols = x2 - x1 + 1 + do j = y1, y2 { + imbuf = imgs2r (im, x1, x2, j, j) + call amovr (Memr[imbuf], Memr[ptr], ncols) + ptr = ptr + ncols + } + + return (buf) +end + + +# DP_CHKSR -- Check the input subraster for bad pixels. + +int procedure dp_chksr (x, y, sr, xdim, ydim, x1, y1, psfrad, fitrad, lowbad, + highbad, saturated, dmin, dmax, gmax) + +real x, y # position of the star +real sr[xdim,ydim] # the data subraster +int xdim, ydim # the dimensions of the subraster +int x1, y1 # the lower left hand coordinates of the array +real psfrad # the psf radius +real fitrad # the fitting radius +real lowbad, highbad # the good data limits +int saturated # is the star saturated +real dmin, dmax # output data limits +real gmax # maximum good data limit + +int i,j +real pradsq, fradsq, dy2, r2 + +begin + pradsq = psfrad * psfrad + fradsq = fitrad * fitrad + dmin = MAX_REAL + dmax = -MAX_REAL + gmax = -MAX_REAL + saturated = NO + + # Loop through the pixels looking for bad values. + do j = 1, ydim { + dy2 = (y - (y1 + j -1)) ** 2 + if (dy2 > pradsq) + next + do i = 1, xdim { + r2 = (x - (x1 + i - 1)) ** 2 + dy2 + if (r2 > pradsq) + next + if (sr[i,j] < dmin) + dmin = sr[i,j] + if (sr[i,j] > dmax) + dmax = sr[i,j] + if (sr[i,j] < lowbad || sr[i,j] > highbad) { + if (r2 <= fradsq) { + if (sr[i,j] < lowbad) + return (ERR) + else if (saturated == NO) + saturated = YES + } + } else { + if (sr[i,j] > gmax) + gmax = sr[i,j] + } + } + } + + return (OK) +end diff --git a/noao/digiphot/daophot/psf/dpptconfirm.x b/noao/digiphot/daophot/psf/dpptconfirm.x new file mode 100644 index 00000000..21047440 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpptconfirm.x @@ -0,0 +1,21 @@ +# DP_PTCONFIRM -- Confirm the critical PSTSELECT parameters. + +procedure dp_ptconfirm (dao) + +pointer dao # pointer to the daophot structure + +begin + call printf ("\n") + + # Confirm the psf radius. + call dp_vpsfrad (dao) + + # Confirm the fitting radius. + call dp_vfitrad (dao) + + # Confirm the data minimum and maximum values. + call dp_vdatamin (dao) + call dp_vdatamax (dao) + + call printf ("\n") +end diff --git a/noao/digiphot/daophot/psf/dppwrtgrp.x b/noao/digiphot/daophot/psf/dppwrtgrp.x new file mode 100644 index 00000000..e0aad167 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppwrtgrp.x @@ -0,0 +1,642 @@ +include <time.h> +include <tbset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + + +# DP_WNEISTARS -- Identify the neighbour stars of the psf stars and write them +# out in groups. + +procedure dp_wneistars (dao, im, psfgr) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int psfgr # the output group file descriptor + +bool newfile +int top_star, psf_star, nei1, nei2 +pointer psf +int dp_neistars() + +begin + psf = DP_PSF(dao) + newfile = true + top_star = DP_PNUM(psf) + 1 + do psf_star = 1, DP_PNUM(psf) { + if (dp_neistars (dao, psf_star, top_star, nei1, nei2) <= 0) + ; + call dp_pwrtgrp (dao, im, psfgr, psf_star, nei1, nei2, newfile) + if (newfile) + newfile = false + } +end + + +# DP_NEISTARS -- Identify the neighbours and friends of the neighbours for +# an individual PSF star. + +int procedure dp_neistars (dao, psf_star, top_star, nei1, nei2) + +pointer dao # pointer to the daophot structure +int psf_star # the psf star in question +int top_star # pointer to the current top_star +int nei1, nei2 # pointer to to the list of neighbours + +int j, nei_star +pointer apsel +real rsq1, rsq2, rsq + +begin + # Define some pointers + apsel = DP_APSEL(dao) + + # These are thhe values used by daophot ii. I have decided to keep + # the old values because I have kept the old grouping algorithm. + #rsq1 = (1.5 * DP_PSFRAD(dao) + 2.0 * DP_FITRAD(dao) + 1.0) ** 2 + #rsq2 = (2.0 * DP_FITRAD(dao) + 1.0) ** 2 + + rsq1 = (DP_PSFRAD(dao) + 2.0 * DP_FITRAD(dao) + 1.0) ** 2 + rsq2 = (2.0 * DP_FITRAD(dao)) ** 2 + + # Find the neighbour stars for a given psf star. This step is the + # same as the daophot ii step although I am using a smaller critical + # radius. + + nei1 = top_star + for (j = top_star; j <= DP_APNUM(apsel); j = j + 1) { + rsq = (Memr[DP_APXCEN(apsel)+j-1] - + Memr[DP_APSEL(apsel)+psf_star-1]) ** 2 + + (Memr[DP_APYCEN(apsel)+j-1] - + Memr[DP_APYCEN(apsel)+psf_star-1]) ** 2 + if (rsq > rsq1) + next + call dp_5swap (j, top_star, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)]) + top_star = top_star +1 + } + nei2 = top_star - 1 + + # Find the friends of the neighbor stars. I do this on a per psf + # star basis. I do not find friends of all the neighbor stars + # only the neighbour stars for a particular psf star. This is + # because I found the daophot ii algorithm could produce too + # may odd stars. + + do nei_star = nei1, nei2 { + for (j = top_star; j <= DP_APNUM(apsel); j = j + 1) { + rsq = (Memr[DP_APXCEN(apsel)+j-1] - + Memr[DP_APSEL(apsel)+nei_star-1]) ** 2 + + (Memr[DP_APYCEN(apsel)+j-1] - + Memr[DP_APYCEN(apsel)+nei_star-1]) ** 2 + if (rsq > rsq2) + next + call dp_5swap (j, top_star, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)]) + top_star = top_star +1 + } + } + nei2 = top_star - 1 + + return (nei2 - nei1 + 1) +end + + +# DP_PWRTGRP -- Add a group to the PSF output group file. + +procedure dp_pwrtgrp (dao, im, psfgr, psf_star, nei1_star, nei2_star, new_table) + +pointer dao # the pointer to the daophot structure +pointer im # the input image descriptor +int psfgr # the group file descriptor +int psf_star # the psf star index +int nei1_star, nei2_star # the first and last neighbour star indices +bool new_table # should a new table be created + +int group + +begin + # Chexk to see if the PSF group file is open. + if (psfgr == NULL) + return + + # Create the table. + if (new_table) { + + # Initialize. + group = 1 + + # Make a new group. + if (DP_TEXT(dao) == YES) + call dp_pxnewgrp (dao, im, psfgr, psf_star, nei1_star, + nei2_star, group) + else + call dp_ptnewgrp (dao, im, psfgr, psf_star, nei1_star, + nei2_star, group) + + } else { + + # Increment. + group = group + 1 + + # Add to the file. + if (DP_TEXT(dao) == YES) + call dp_pxaddgrp (dao, im, psfgr, psf_star, nei1_star, + nei2_star, group) + else + call dp_ptaddgrp (dao, im, psfgr, psf_star, nei1_star, + nei2_star, group) + } + +end + + +define NCOLUMN 5 + +# DP_WPLIST -- Write the list of psf stars to the output psf star list. + +procedure dp_wplist (dao, im, opst) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int opst # the output psf star list descriptor + +pointer apsel, psf, sp, ocolpoint, tx, ty +int dp_stati() +bool itob() + +begin + # Get some daophot pointers. + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + + # Allocate some working memory. + call smark (sp) + call salloc (ocolpoint, NCOLUMN, TY_POINTER) + call salloc (tx, DP_PNUM(psf), TY_REAL) + call salloc (ty, DP_PNUM(psf), TY_REAL) + + # Initialize the output file. + if (dp_stati (dao, TEXT) == YES) { + call dp_pxgrppars (dao, opst) + call dp_xpbanner (opst) + } else { + call dp_tpdefcol (opst, Memi[ocolpoint]) + call dp_ptgrppars (dao, opst) + } + + # Write out the stars. + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[tx], Memr[ty], DP_PNUM(psf)) + call dp_wpstars (opst, Memi[ocolpoint], itob (dp_stati (dao, TEXT)), + Memi[DP_APID(apsel)], Memr[tx], Memr[ty], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], DP_PNUM(psf)) + + # Free memory. + call sfree (sp) +end + + +define PGR_NAMESTR "#N%4tID%10tGROUP%16tXCENTER%26tYCENTER%36tMAG%48t\ +MSKY%80t\\\n" +define PGR_UNITSTR "#U%4t##%10t##%16tpixels%26tpixels%36tmagnitudes%48t\ +counts%80t\\\n" +define PGR_FORMATSTR "#F%4t%%-9d%10t%%-6d%16t%%-10.3f%26t%%-10.3f%36t\ +%%-12.3f%48t%%-15.7g%80t \n" +define PGR_DATASTR "%-9d%10t%-6d%16t%-10.3f%26t%-10.3f%36t%-12.3f%48t\ +%-15.7g%80t \n" + +# DP_PXNEWGRP -- Create a new PSF output group text file and write the +# first group to it. + +procedure dp_pxnewgrp (dao, im, tp, psf_star, nei1_star, nei2_star, group) + +pointer dao # pointer to the daophot structure. +pointer im # pointer to the input image +int tp # the output file descriptor +int psf_star # the psf star index +int nei1_star, nei2_star # the first and last neighbour star indices +int group # current group + +real tx, ty +pointer apsel +int i + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Define some pointers. + apsel = DP_APSEL(dao) + + # Write out the header parameters. + call dp_pxgrppars (dao, tp) + + # Set up the column definitions. + call fprintf (tp, "#\n") + call fprintf (tp, PGR_NAMESTR) + call fprintf (tp, PGR_UNITSTR) + call fprintf (tp, PGR_FORMATSTR) + call fprintf (tp, "#\n") + + # Write out the psf star. + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+psf_star-1], + Memr[DP_APYCEN(apsel)+psf_star-1], tx, ty, 1) + call fprintf (tp, PGR_DATASTR) + call pargi (Memi[DP_APID(apsel)+psf_star-1]) + call pargi (group) + call pargr (tx) + call pargr (ty) + call pargr (Memr[DP_APMAG(apsel)+psf_star-1]) + call pargr (Memr[DP_APMSKY(apsel)+psf_star-1]) + + # Write out the neighbour stars. + do i = nei1_star, nei2_star { + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], tx, ty, 1) + call fprintf (tp, PGR_DATASTR) + call pargi (Memi[DP_APID(apsel)+i-1]) + call pargi (group) + call pargr (tx) + call pargr (ty) + call pargr (Memr[DP_APMAG(apsel)+i-1]) + call pargr (Memr[DP_APMSKY(apsel)+i-1]) + } +end + + +# DP_PTNEWGRP -- Create a new PSF output group ST table and add to a new +# group to it. + +procedure dp_ptnewgrp (dao, im, tp, psf_star, nei1_star, nei2_star, group) + +pointer dao # pointer to the daophot structure. +pointer im # the input image descriptor +pointer tp # pointer to table +int psf_star # the psf star index +int nei1_star, nei2_star # the first and last psf star indices +int group # current group + +real tx, ty +pointer sp, colnames, colunits, colformat, coldtype, collen, colpoint +pointer apsel +int i, j, row + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Define some pointers. + apsel = DP_APSEL(dao) + + # Allocate space for table definition. + call smark (sp) + call salloc (colpoint, PSF_NOUTCOLS, TY_INT) + call salloc (colnames, PSF_NOUTCOLS * (SZ_COLNAME + 1), TY_CHAR) + call salloc (colunits, PSF_NOUTCOLS * (SZ_COLUNITS + 1), TY_CHAR) + call salloc (colformat, PSF_NOUTCOLS * (SZ_COLFMT + 1), TY_CHAR) + call salloc (coldtype, PSF_NOUTCOLS, TY_INT) + call salloc (collen, PSF_NOUTCOLS, TY_INT) + + # Set up the column definitions. + call strcpy (ID, Memc[colnames], SZ_COLNAME) + call strcpy (GROUP, Memc[colnames+SZ_COLNAME+1], SZ_COLNAME) + call strcpy (XCENTER, Memc[colnames+2*SZ_COLNAME+2], SZ_COLNAME) + call strcpy (YCENTER, Memc[colnames+3*SZ_COLNAME+3], SZ_COLNAME) + call strcpy (MAG, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME) + call strcpy (SKY, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME) + + # Set up the format definitions. + call strcpy ("%6d", Memc[colformat], SZ_COLFMT) + call strcpy ("%6d", Memc[colformat+SZ_COLFMT+1], SZ_COLFMT) + call strcpy ("%10.2f", Memc[colformat+2*SZ_COLFMT+2], SZ_COLFMT) + call strcpy ("%10.2f", Memc[colformat+3*SZ_COLFMT+3], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT) + call strcpy ("%15.7g", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT) + + # Define the column units. + call strcpy ("NUMBER", Memc[colunits], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+SZ_COLUNITS+1], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+2*SZ_COLUNITS+2], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+3*SZ_COLUNITS+3], SZ_COLUNITS) + call strcpy ("MAGNITUDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS) + call strcpy ("COUNTS", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS) + + # Define the datatypes of the columns. + Memi[coldtype] = TY_INT + Memi[coldtype+1] = TY_INT + Memi[coldtype+2] = TY_REAL + Memi[coldtype+3] = TY_REAL + Memi[coldtype+4] = TY_REAL + Memi[coldtype+5] = TY_REAL + + # Initialize the column length parameter. + do i = 1, PSF_NOUTCOLS { + j = i - 1 + Memi[collen+j] = 1 + } + + # Define the table. + call tbcdef (tp, Memi[colpoint], Memc[colnames], Memc[colunits], + Memc[colformat], Memi[coldtype], Memi[collen], PSF_NOUTCOLS) + + # Create the table. + call tbtcre (tp) + + # Write out the psf star. + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+psf_star-1], + Memr[DP_APYCEN(apsel)+psf_star-1], tx, ty, 1) + call tbrpti (tp, Memi[colpoint], Memi[DP_APID(apsel)+psf_star-1], 1, 1) + call tbrpti (tp, Memi[colpoint+1], group, 1, 1) + call tbrptr (tp, Memi[colpoint+2], tx, 1, 1) + call tbrptr (tp, Memi[colpoint+3], ty, 1, 1) + call tbrptr (tp, Memi[colpoint+4], Memr[DP_APMAG(apsel)+psf_star-1], + 1, 1) + call tbrptr (tp, Memi[colpoint+5], Memr[DP_APMSKY(apsel)+psf_star-1], + 1, 1) + + # Now write out the group. + row = 2 + do i = nei1_star, nei2_star { + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], tx, ty, 1) + call tbrpti (tp, Memi[colpoint], Memi[DP_APID(apsel)+i-1], 1, row) + call tbrpti (tp, Memi[colpoint+1], group, 1, row) + call tbrptr (tp, Memi[colpoint+2], tx, 1, row) + call tbrptr (tp, Memi[colpoint+3], ty, 1, row) + call tbrptr (tp, Memi[colpoint+4], Memr[DP_APMAG(apsel)+i-1], + 1, row) + call tbrptr (tp, Memi[colpoint+5], Memr[DP_APMSKY(apsel)+i-1], + 1, row) + row = row + 1 + } + + # Add the header parameters to the table. + call dp_ptgrppars (dao, tp) + + call sfree (sp) +end + + +# DP_PTGRPPARS -- Add parameters to the header of the output PSF group ST +# table. + +procedure dp_ptgrppars (dao, tp) + +pointer dao # pointer to the daophot structure +pointer tp # pointer to the table + +pointer sp, outstr, date, time +int envfind() + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Allocate workin space. + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + call salloc (date, SZ_DATE, TY_CHAR) + call salloc (time, SZ_DATE, TY_CHAR) + + # Write the task and date identifiers. + if (envfind ("version", Memc[outstr], SZ_LINE) <= 0) + call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE) + call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE) + call tbhadt (tp, "IRAF", Memc[outstr]) + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call tbhadt (tp, "USER", Memc[outstr]) + call gethost (Memc[outstr], SZ_LINE) + call tbhadt (tp, "HOST", Memc[outstr]) + call dp_date (Memc[date], Memc[time], SZ_DATE) + call tbhadt (tp, "DATE", Memc[date]) + call tbhadt (tp, "TIME", Memc[time]) + + # Define the package and task. + call tbhadt (tp, "PACKAGE", "daophot") + call tbhadt (tp, "TASK", "psf") + + # Define the input and output files. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "IMAGE", Memc[outstr]) + + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PHOTFILE", Memc[outstr]) + if (DP_COORDS(dao) == EOS) + call tbhadt (tp, "PSTFILE", "\"\"") + else { + call dp_froot (DP_COORDS(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PSTFILE", Memc[outstr]) + } + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PSFIMAGE", DP_PSFIMAGE(dao)) + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "OPSTFILE", Memc[outstr]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "GRPSFILE", Memc[outstr]) + + # Data dependent parameters. + call tbhadr (tp, "SCALE", DP_SCALE(dao)) + + # Observing parameters. + call tbhadt (tp, "OTIME", DP_OTIME(dao)) + call tbhadt (tp, "IFILTER", DP_IFILTER(dao)) + call tbhadr (tp, "XAIRMASS", DP_XAIRMASS(dao)) + + # Grouping parameters. + call tbhadr (tp, "PSFRAD", DP_SPSFRAD(dao)) + call tbhadr (tp, "FITRAD", DP_SFITRAD(dao)) + + call sfree(sp) +end + + +# DP_PXGRPPARS -- Add parameters to the header of the output PSF group text +# file. + +procedure dp_pxgrppars (dao, tp) + +pointer dao # pointer to the daophot structure +int tp # the output file descriptor + +pointer sp, outstr, date, time +int envfind() + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Allocate workin space. + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + call salloc (date, SZ_DATE, TY_CHAR) + call salloc (time, SZ_DATE, TY_CHAR) + + # Write the task and date identifiers. + if (envfind ("version", Memc[outstr], SZ_LINE) <= 0) + call strcpy ("NOAO/IRAF", Memc[outstr], SZ_LINE) + call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE) + call dp_sparam (tp, "IRAF", Memc[outstr], "version", "") + if (envfind ("userid", Memc[outstr], SZ_LINE) > 0) + call dp_sparam (tp, "USER", Memc[outstr], "name", "") + call gethost (Memc[outstr], SZ_LINE) + call dp_sparam (tp, "HOST", Memc[outstr], "computer", "") + call dp_date (Memc[date], Memc[time], SZ_DATE) + call dp_sparam (tp, "DATE", Memc[date], "yyyy-mm-dd", "") + call dp_sparam (tp, "TIME", Memc[time], "hh:mm:ss", "") + + # Define the package and task. + call dp_sparam (tp, "PACKAGE", "daophot", "name", "") + call dp_sparam (tp, "TASK", "psf", "name", "") + + # Define the input and output files. + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "IMAGE", Memc[outstr], "imagename", "") + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "PHOTFILE", Memc[outstr], "filename", "") + if (DP_COORDS(dao) == EOS) + call dp_sparam (tp, "PSTFILE", "\"\"", "filename", "") + else { + call dp_froot (DP_COORDS(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "PSTFILE", Memc[outstr], "filename", "") + } + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "PSFIMAGE", Memc[outstr], "imagename", "") + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "GRPSFILE", Memc[outstr], "filename", "") + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "OPSTFILE", Memc[outstr], "filename", "") + + # Define the data dependent parameters. + call dp_rparam (tp, "SCALE", DP_SCALE(dao), "units/pix", "") + + # Observing parameters. + call dp_sparam (tp, "OTIME", DP_OTIME(dao), "timeunit", "") + call dp_sparam (tp, "IFILTER", DP_IFILTER(dao), "filter", "") + call dp_rparam (tp, "XAIRMASS", DP_XAIRMASS(dao), "number", "") + + # Grouping parameters. + call dp_rparam (tp, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", "") + call dp_rparam (tp, "FITRAD", DP_SFITRAD(dao), "scaleunit", "") + + call sfree(sp) +end + + +# DP_PXADDGRP -- Add a new group to the existing PSF output group text file. + +procedure dp_pxaddgrp (dao, im, tp, psf_star, nei1_star, nei2_star, group) + +pointer dao # pointer to daophot structure +pointer im # the input image descriptor +int tp # the output file descriptor +int psf_star # the psf star index +int nei1_star, nei2_star # the first and last neighbour star indices +int group # current group + +real tx, ty +pointer apsel +int i + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Get some daophot pointers. + apsel = DP_APSEL(dao) + + # Write out the psf star. + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+psf_star-1], + Memr[DP_APYCEN(apsel)+psf_star-1], tx, ty, 1) + call fprintf (tp, PGR_DATASTR) + call pargi (Memi[DP_APID(apsel)+psf_star-1]) + call pargi (group) + call pargr (tx) + call pargr (ty) + call pargr (Memr[DP_APMAG(apsel)+psf_star-1]) + call pargr (Memr[DP_APMSKY(apsel)+psf_star-1]) + + # Now write out the group. + do i = nei1_star, nei2_star { + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], tx, ty, 1) + call fprintf (tp, PGR_DATASTR) + call pargi (Memi[DP_APID(apsel)+i-1]) + call pargi (group) + call pargr (tx) + call pargr (ty) + call pargr (Memr[DP_APMAG(apsel)+i-1]) + call pargr (Memr[DP_APMSKY(apsel)+i-1]) + } +end + + +# DP_PTADDGRP -- Add a new group to the existing PSF output group ST table. + +procedure dp_ptaddgrp (dao, im, tp, psf_star, nei1_star, nei2_star, group) + +pointer dao # pointer to daophot structure +pointer im # the input image descriptor +pointer tp # pointer to output table +int psf_star # the psf star index +int nei1_star, nei2_star # the first and last neighbor star indices +int group # current group + +real tx, ty +pointer apsel, sp, colpoint +int i, nrows +int tbpsta() + +begin + # Check to see if the PSF group file is open. + if (tp == NULL) + return + + # Allocate space for the column pointers. + call smark (sp) + call salloc (colpoint, PSF_NOUTCOLS, TY_INT) + + # Get some daophot pointers. + apsel = DP_APSEL(dao) + + # Find the number of rows in the table and add on at the end. + nrows = tbpsta (tp, TBL_NROWS) + + # Write out the psf star + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+psf_star-1], + Memr[DP_APYCEN(apsel)+psf_star-1], tx, ty, 1) + nrows = nrows + 1 + call tbrpti (tp, Memi[colpoint], Memi[DP_APID(apsel)+psf_star-1], + 1, nrows) + call tbrpti (tp, Memi[colpoint+1], group, 1, nrows) + call tbrptr (tp, Memi[colpoint+2], tx, 1, nrows) + call tbrptr (tp, Memi[colpoint+3], ty, 1, nrows) + call tbrptr (tp, Memi[colpoint+4], Memr[DP_APMAG(apsel)+psf_star-1], + 1, nrows) + call tbrptr (tp, Memi[colpoint+5], Memr[DP_APMSKY(apsel)+psf_star-1], + 1, nrows) + + # Now write out the group. + do i = nei1_star, nei2_star { + nrows = nrows + 1 + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], tx, ty, 1) + call tbrpti (tp, Memi[colpoint], Memi[DP_APID(apsel)+i-1], 1, nrows) + call tbrpti (tp, Memi[colpoint+1], group, 1, nrows) + call tbrptr (tp, Memi[colpoint+2], tx, 1, nrows) + call tbrptr (tp, Memi[colpoint+3], ty, 1, nrows) + call tbrptr (tp, Memi[colpoint+4], Memr[DP_APMAG(apsel)+i-1], 1, + nrows) + call tbrptr (tp, Memi[colpoint+5], Memr[DP_APMSKY(apsel)+i-1], 1, + nrows) + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dppwselmer.x b/noao/digiphot/daophot/psf/dppwselmer.x new file mode 100644 index 00000000..24684b67 --- /dev/null +++ b/noao/digiphot/daophot/psf/dppwselmer.x @@ -0,0 +1,220 @@ +include <tbset.h> +include "../lib/apseldef.h" + +define NCOLUMN 5 + +define PS_DATA1STR "%-9d%10t%-10.3f%20t%-10.3f%30t%-12.3f%42t%-15.7g%80t \n" + +# DP_XPSELMER -- Write the output photometry record to a text file. + +procedure dp_xpselmer (tpout, id, x, y, mag, sky) + +pointer tpout # pointer to the output table +int id # id of the star +real x, y # position of the star +real mag # magnitude of the star +real sky # value of sky + +begin + call fprintf (tpout, PS_DATA1STR) + call pargi (id) + call pargr (x) + call pargr (y) + call pargr (mag) + call pargr (sky) +end + + +# DP_TPSELMER -- Write out the PSF stars into an ST Table. + +procedure dp_tpselmer (tp_out, id, x, y, mag, sky, colpoint, row) + +int tp_out # the output table descriptor +int id # the object id +real x # the object x coordinate +real y # the object y coordinate +real mag # the object mangitude +real sky # the object sky value +int colpoint[ARB] # the column pointers +int row # current table row + +begin + # Write out the data. + call tbrpti (tp_out, colpoint[1], id, 1, row) + call tbrptr (tp_out, colpoint[2], x, 1, row) + call tbrptr (tp_out, colpoint[3], y, 1, row) + call tbrptr (tp_out, colpoint[4], mag, 1, row) + call tbrptr (tp_out, colpoint[5], sky, 1, row) +end + + +# DP_XPSELPARS -- Add various parameters to the header of the photometry table. + +procedure dp_xpselpars (tp, image, maxnpsf, scale, psfrad, fitrad) + +pointer tp # pointer to the table +char image[ARB] # input image name +int maxnpsf # maximum number of psfstars +real scale # the image scale +real psfrad # the psf radius +real fitrad # the fitting radius + +pointer sp, str + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Add the image name nad maxnpsf parameters. + call dp_imroot (image, Memc[str], SZ_FNAME) + call dp_sparam (tp, "IMAGE", Memc[str], "imagename", "") + call dp_iparam (tp, "MAXNPSF", maxnpsf, "number", "") + call dp_rparam (tp, "NEWSCALE", scale, "units", "") + call dp_rparam (tp, "PSFRAD", psfrad, "scaleunit", "") + call dp_rparam (tp, "FITRAD", fitrad, "scaleunit", "") + + call sfree (sp) +end + + +define PS_NAME1STR "#N%4tID%10tXCENTER%20tYCENTER%30tMAG%42tMSKY%80t\\\n" +define PS_UNIT1STR "#U%4t##%10tpixels%20tpixels%30tmagnitudes%42tcounts\ +%80t\\\n" +define PS_FORMAT1STR "#F%4t%%-9d%10t%%-10.3f%20t%%-10.3f%30t%%-12.3f%42t\ +%%-15.7g%80t \n" + +# DP_XPBANNER -- Create a new text file banner. + +procedure dp_xpbanner (tp) + +pointer tp # pointer to the output file + +begin + # Print out the banner file. + call fprintf (tp, "#\n") + call fprintf (tp, PS_NAME1STR) + call fprintf (tp, PS_UNIT1STR) + call fprintf (tp, PS_FORMAT1STR) + call fprintf (tp, "#\n") +end + + +# DP_TPDEFCOL -- Define the columns for the output table + +procedure dp_tpdefcol (tp, colpoint) + +pointer tp # pointer to the output table +int colpoint[ARB] # array of column pointers + +int i +pointer sp, colnames, colunits, colformat, col_dtype, col_len + +begin + # Allocate space for the table definition. + call smark (sp) + call salloc (colnames, NCOLUMN * (SZ_COLNAME + 1), TY_CHAR) + call salloc (colunits, NCOLUMN * (SZ_COLUNITS + 1), TY_CHAR) + call salloc (colformat, NCOLUMN * (SZ_COLFMT + 1), TY_CHAR) + call salloc (col_dtype, NCOLUMN, TY_INT) + call salloc (col_len, NCOLUMN, TY_INT) + + # Set up the column definitions. + call strcpy (ID, Memc[colnames], SZ_COLNAME) + call strcpy (XCENTER, Memc[colnames+SZ_COLNAME+1], SZ_COLNAME) + call strcpy (YCENTER, Memc[colnames+2*SZ_COLNAME+2], SZ_COLNAME) + call strcpy (MAG, Memc[colnames+3*SZ_COLNAME+3], SZ_COLNAME) + call strcpy (SKY, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME) + + # Define the column formats. + call strcpy ("%6d", Memc[colformat], SZ_COLFMT) + call strcpy ("10.3f", Memc[colformat+SZ_COLFMT+1], SZ_COLFMT) + call strcpy ("10.3f", Memc[colformat+2*SZ_COLFMT+2], SZ_COLFMT) + call strcpy ("12.3f", Memc[colformat+3*SZ_COLFMT+3], SZ_COLFMT) + call strcpy ("15.7g", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT) + + # Define the column units. + call strcpy ("NUMBER", Memc[colunits], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+SZ_COLUNITS+1], SZ_COLUNITS) + call strcpy ("PIXELS", Memc[colunits+2*SZ_COLUNITS+2], SZ_COLUNITS) + call strcpy ("MAGNITUDES", Memc[colunits+3*SZ_COLUNITS+3], SZ_COLUNITS) + call strcpy ("ADU", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS) + + # Define the column data types. + Memi[col_dtype] = TY_INT + Memi[col_dtype+1] = TY_REAL + Memi[col_dtype+2] = TY_REAL + Memi[col_dtype+3] = TY_REAL + Memi[col_dtype+4] = TY_REAL + + # Define the column lengths. + do i = 1, NCOLUMN + Memi[col_len+i-1] = 1 + + # Define and create the table. + call tbcdef (tp, colpoint, Memc[colnames], Memc[colunits], + Memc[colformat], Memi[col_dtype], Memi[col_len], NCOLUMN) + call tbtcre (tp) + + call sfree (sp) +end + + +# DP_TPSELPARS -- Add various parameters to the header of the photometry table. + +procedure dp_tpselpars (tp, image, maxnpsf, scale, psfrad, fitrad) + +pointer tp # pointer to the table +char image[ARB] # the input image name +int maxnpsf # maximum number of psf stars +real scale # the image scale +real psfrad # the psf radius +real fitrad # the fitting radius + +pointer sp, str + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Add the min_group and max_group parameters. + call dp_imroot (image, Memc[str], SZ_FNAME) + call tbhadt (tp, "IMAGE", Memc[str]) + call tbhadi (tp, "MAXNPSF", maxnpsf) + call tbhadr (tp, "SCALE", scale) + call tbhadr (tp, "PSFRAD", psfrad) + call tbhadr (tp, "FITRAD", fitrad) + + call sfree (sp) +end + + +# DP_WPSTARS -- Write the psf stars to the output file. + +procedure dp_wpstars (tp_out, colpoint, text_file, ids, xcen, ycen, mag, + sky, npsf) + +int tp_out # the output file descriptor +int colpoint[ARB] # array of column pointers +bool text_file # is the output file a text file +int ids[ARB] # array of star ids +real xcen[ARB] # array of x coordinates +real ycen[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real sky[ARB] # array of sky values +int npsf # the number of stars + +int istar, row + +begin + row = 0 + do istar = 1, npsf { + if (text_file) + call dp_xpselmer (tp_out, ids[istar], xcen[istar], ycen[istar], + mag[istar], sky[istar]) + else { + row = row + 1 + call dp_tpselmer (tp_out, ids[istar], xcen[istar], ycen[istar], + mag[istar], sky[istar], colpoint, row) + } + } +end diff --git a/noao/digiphot/daophot/psf/dpqverify.x b/noao/digiphot/daophot/psf/dpqverify.x new file mode 100644 index 00000000..6fff239f --- /dev/null +++ b/noao/digiphot/daophot/psf/dpqverify.x @@ -0,0 +1,68 @@ +include <fset.h> + +define QUERY "[Hit return to continue, q to quit, w to force update of PSF]" + +# DP_QVERIFY -- Print a message in the status line asking the user if they +# really want to quit, returning YES if they really want to quit, NO otherwise. + +int procedure dp_qverify (dao, im, psfim, opst, psfgr, psf_new, psf_computed, + psf_written) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +pointer psfim # the output psf image descriptor +int opst # the output psf star list file descriptor +int psfgr # the output psf group file descriptor +bool psf_new # has the psf star list been defined ? +bool psf_computed # is the psf fit defined ? +bool psf_written # has the psf been saved ? + +int ch +pointer tty +bool dp_updatepsf() +int getci() +pointer ttyodes() + +begin + # Print status warning message. + if (psf_new) + call printf ("Warning: The PSF star list is undefined.\n") + else if (! psf_computed) + call printf ("Warning: The PSF fit is not current.\n") + else if (! psf_written) + call printf ("Warning: The PSF has not been saved.\n") + + # Open terminal and print query. + tty = ttyodes ("terminal") + call ttyclearln (STDOUT, tty) + call ttyso (STDOUT, tty, YES) + call printf (QUERY) + call flush (STDOUT) + + # Get character. + call fseti (STDIN, F_RAW, YES) + if (getci (STDIN, ch) == EOF) + ; + + # Reset and close terminal. + call fseti (STDIN, F_RAW, NO) + call ttyso (STDOUT, tty, NO) + call ttyclearln (STDOUT, tty) + call printf ("\n") + call flush (STDOUT) + call ttycdes (tty) + + # Return YES for the quit command, otherwise NO. + if (ch == 'q') { + return (YES) + } else if (ch == 'w') { + if (dp_updatepsf (dao, im, psfim, opst, psfgr, psf_new, + psf_computed, psf_written)) { + psf_computed = true + psf_written = true + } + return (NO) + } else { + return (NO) + } +end diff --git a/noao/digiphot/daophot/psf/dpradpsf.x b/noao/digiphot/daophot/psf/dpradpsf.x new file mode 100644 index 00000000..7e649e7e --- /dev/null +++ b/noao/digiphot/daophot/psf/dpradpsf.x @@ -0,0 +1,75 @@ +include <gset.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +define FRACTION 0.10 + +# DP_RADPSF -- Draw a radial profile plot of a data subraster containing a +# candidate psf star. + +procedure dp_radpsf (dao, subras, ncols, nlines, x1, y1, title, gp) + +pointer dao # pointer to DAOPHOT structure +real subras[ncols,nlines] # data subraster +int ncols, nlines # dimesnions of subraster +int x1, y1 # coordinates of left hand corner +char title[ARB] # title string +pointer gp # pointer to graphics descriptor + +int npts +pointer psf, sp, radius, intensity, str +real ymin, ymax, r1, r2, i1, i2 +int dp_rivectors() + +begin + # Get the pointer to the DAOPHOT PSF fitting substructure. + psf = DP_PSF (dao) + + # Allocate temporary space. + call smark (sp) + call salloc (radius, ncols * nlines, TY_REAL) + call salloc (intensity, ncols * nlines, TY_REAL) + call salloc (str, SZ_LINE, TY_CHAR) + + # Compute the radial profile. + npts = dp_rivectors (subras, ncols, nlines, x1, y1, DP_CUR_PSFX(psf), + DP_CUR_PSFY(psf), DP_PSFRAD(dao), Memr[radius], Memr[intensity]) + call alimr (Memr[intensity], npts, ymin, ymax) + + # Make the plot. + call gclear (gp) + r1 = -FRACTION * DP_PSFRAD(dao) + r2 = DP_PSFRAD(dao) + FRACTION * DP_PSFRAD(dao) + i1 = ymin - FRACTION * (ymax - ymin) + i2 = ymax + FRACTION * (ymax - ymin) + call gswind (gp, r1, r2, i1, i2) + call glabax (gp, title, "Radius (pixels)", "Intensity (counts)") + call gpmark (gp, Memr[radius], Memr[intensity], npts, GM_PLUS, 1.0, + 1.0) + + # Mark the zero radius line. + call gamove (gp, 0.0, i1) + call gadraw (gp, 0.0, i2) + + # Mark the sky level. + call gamove (gp, r1, DP_CUR_PSFSKY(psf)) + call gadraw (gp, r2, DP_CUR_PSFSKY(psf)) + + # Mark the half-width at half-maximum. + call gamove (gp, DP_FWHMPSF(dao) / 2.0, i1) + call gadraw (gp, DP_FWHMPSF(dao) / 2.0, i2) + call sprintf (Memc[str], SZ_LINE, "Half-width half-maximum = %0.2f") + call pargr (DP_FWHMPSF(dao) / 2.0) + call gtext (gp, DP_FWHMPSF(dao) / 2.0, i2, Memc[str], + "q=h;u=180;v=t;p=r") + + # Mark the fitting radius. + call gamove (gp, DP_FITRAD(dao), i1) + call gadraw (gp, DP_FITRAD(dao), i2) + call sprintf (Memc[str], SZ_LINE, "Fitting radius = %0.2f") + call pargr (DP_FITRAD(dao)) + call gtext (gp, DP_FITRAD(dao), i2, Memc[str], "q=h;u=180;v=t;p=r") + + call gflush (gp) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dprmpsf.x b/noao/digiphot/daophot/psf/dprmpsf.x new file mode 100644 index 00000000..dd76c4d4 --- /dev/null +++ b/noao/digiphot/daophot/psf/dprmpsf.x @@ -0,0 +1,156 @@ +include <imhdr.h> +include <fset.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + + +# DP_OPPSF -- Open the current psf image and the psf group file. + +procedure dp_oppsf (dao, psfim, opst, psfgr) + +pointer dao # pointer to the daophot structure +pointer psfim # pointer to the psf image +int opst # the psf star list descriptor +int psfgr # the psf group file descriptor + +pointer sp, str +int open(), dp_stati(), dp_pstati() +pointer immap(), tbtopn() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Reopen the PSF star and group files. + call dp_stats (dao, OUTREJFILE, Memc[str], SZ_FNAME) + if (dp_stati (dao, TEXT) == YES) + opst = open (Memc[str], NEW_FILE, TEXT_FILE) + else + opst = tbtopn (Memc[str], NEW_FILE, 0) + call dp_stats (dao, OUTPHOTFILE, Memc[str], SZ_FNAME) + if (dp_stati (dao, TEXT) == YES) + psfgr = open (Memc[str], NEW_FILE, TEXT_FILE) + else + psfgr = tbtopn (Memc[str], NEW_FILE, 0) + + # Reopen the psf image. + call dp_stats (dao, PSFIMAGE, Memc[str], SZ_FNAME) + psfim = immap (Memc[str], NEW_IMAGE, dp_pstati (dao, LENUSERAREA)) + + call sfree (sp) +end + + +# DP_UPDATEPSF -- Update the psf on disk. + +bool procedure dp_updatepsf (dao, im, psfim, opst, psfgr, psf_new, psf_current, + psf_written) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +pointer psfim # pointer to the psf image +int opst # the psf star list file descriptor +int psfgr # the psf group file descriptor +bool psf_new # is the psf star list defined ? +bool psf_current # is the psf fit uptodate +bool psf_written # has the psf been saved on disk + +bool update +pointer sp, str +int dp_fitpsf() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + update = false + if (psfim == NULL) { + call printf ("Warning: The PSF image is undefined\n") + } else if (psfgr == NULL) { + call printf ("Warning: The PSF group file is undefined\n") + } else if (psf_new) { + call printf ("Warning: The PSF star list is undefined\n") + } else if (! psf_current) { + if (dp_fitpsf (dao, im, Memc[str], SZ_LINE) == OK) { + call dp_writepsf (dao, im, psfim) + call dp_wplist (dao, im, opst) + call dp_wneistars (dao, im, psfgr) + update = true + } else { + call printf ("%s\n") + call pargstr (Memc[str]) + } + } else if (! psf_written) { + call dp_writepsf (dao, im, psfim) + call dp_wplist (dao, im, opst) + call dp_wneistars (dao, im, psfgr) + update = true + } + + if (DP_VERBOSE(dao) == YES && update) { + call printf ("\nWriting PSF image %s\n") + call pargstr (IM_HDRFILE (psfim)) + call dp_stats (dao, OUTREJFILE, Memc[str], SZ_FNAME) + call printf ("Writing output PSF star list %s\n") + call pargstr (Memc[str]) + call dp_stats (dao, OUTPHOTFILE, Memc[str], SZ_FNAME) + call printf ("Writing output PSF star group file %s\n") + call pargstr (Memc[str]) + } + + call sfree (sp) + + return (update) +end + + +# DP_RMPSF -- Close and delete the psf image and the psf group file. + +procedure dp_rmpsf (dao, psfim, opst, psfgr) + +pointer dao # pointer to the daophot structure +pointer psfim # pointer to the psf image +int opst # the psf star list file descriptor +int psfgr # the psf group file descriptor + +pointer sp, temp +int dp_stati(), access() + +begin + call smark (sp) + call salloc (temp, SZ_FNAME, TY_CHAR) + + # Delete the psf star file. The access check is necessary because + # an empty tables file is automatically deleted by the system. + if (dp_stati (dao, TEXT) == YES) { + call fstats (opst, F_FILENAME, Memc[temp], SZ_FNAME) + call close (opst) + } else { + call tbtnam (opst, Memc[temp], SZ_FNAME) + call tbtclo (opst) + } + if (access (Memc[temp], 0, 0) == YES) + call delete (Memc[temp]) + opst = NULL + + # Delete the neighbours file. The access check is necessary because + # an empty tables file is automatically deleted by the system. + if (dp_stati (dao, TEXT) == YES) { + call fstats (psfgr, F_FILENAME, Memc[temp], SZ_FNAME) + call close (psfgr) + } else { + call tbtnam (psfgr, Memc[temp], SZ_FNAME) + call tbtclo (psfgr) + } + if (access (Memc[temp], 0, 0) == YES) + call delete (Memc[temp]) + psfgr = NULL + + # Delete the PSF image. + call strcpy (IM_HDRFILE(psfim), Memc[temp], SZ_FNAME) + call imunmap (psfim) + call imdelete (Memc[temp]) + psfim = NULL + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dprstars.x b/noao/digiphot/daophot/psf/dprstars.x new file mode 100644 index 00000000..21ca67d1 --- /dev/null +++ b/noao/digiphot/daophot/psf/dprstars.x @@ -0,0 +1,156 @@ +include <gset.h> +include <tbset.h> +include "../lib/apseldef.h" +include "../lib/psfdef.h" + + +# DP_RPSTARS -- Read in the IDS and x and y positions of the PSF stars. + +procedure dp_rpstars (dao, im, pst, text_file, gd, mgd, id, mkstars, + matchbyid, showplots) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int pst # the psf star list file descriptor +bool text_file # text or table file ? +pointer gd # the graphics descriptor +pointer mgd # the plot file descriptor +pointer id # the display device descriptor +bool mkstars # mark the stars added to the psf +bool matchbyid # match psf stars by id or position +bool showplots # show the psf star plots + +real x, y, mag, rjunk +pointer sp, fields, indices, key +int i, nrow, idno +real dp_pstatr() +int tbpsta(), dp_apsel(), dp_addstar() + +begin + call smark (sp) + call salloc (fields, SZ_LINE, TY_CHAR) + call salloc (indices, PSF_NINCOLS, TY_INT) + + if (text_file) { + call pt_kyinit (key) + Memi[indices] = DP_PAPID + Memi[indices+1] = DP_PAPXCEN + Memi[indices+2] = DP_PAPYCEN + Memi[indices+3] = DP_PAPMAG1 + call dp_gappsf (Memi[indices], Memc[fields], PSF_NINCOLS) + } else { + call dp_tptinit (pst, Memi[indices]) + nrow = tbpsta (pst, TBL_NROWS) + } + + i = 1 + repeat { + + # Read the next star. + + if (text_file) { + if (dp_apsel (key, pst, Memc[fields], Memi[indices], idno, + x, y, rjunk, mag) == EOF) + break + } else { + if (i > nrow) + break + call dp_tptread (pst, Memi[indices], idno, x, y, mag, i) + } + + call dp_win (dao, im, x, y, x, y, 1) + + # Add it to the PSF star list. + if (idno > 0) { + if (matchbyid) { + if (dp_addstar (dao, im, x, y, mag, idno, gd, mgd, + showplots) == OK) { + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr(dao, CUR_PSFX), + dp_pstatr(dao, CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + } + } else { + if (dp_addstar (dao, im, x, y, INDEFR, 0, gd, mgd, + showplots) == OK) { + if (mkstars && id != NULL) { + call gmark (id, dp_pstatr(dao, CUR_PSFX), + dp_pstatr(dao, CUR_PSFY), GM_PLUS, -5.0, -5.0) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + } + } + } + + i = i + 1 + } + + if (text_file) + call pt_kyfree (key) + call sfree (sp) +end + + +# DP_TPTINIT -- Set up the input psf star list ST table column pointers. + +procedure dp_tptinit (pst, column) + +int pst # the psf star list file descriptor +int column[ARB] # array of column pointers + +begin + call tbcfnd (pst, ID, column[1], 1) + if (column[1] == NULL) + call tbcfnd (pst, "ID", column[1], 1) + if (column[1] == NULL) + call error (0, "Error reading ID column from PSF star file\n") + + call tbcfnd (pst, XCENTER, column[2], 1) + if (column[2] == NULL) + call tbcfnd (pst, "XCENTER", column[2], 1) + if (column[2] == NULL) + call error (0, "Error reading XCENTER column from PSF star file\n") + + call tbcfnd (pst, YCENTER, column[3], 1) + if (column[3] == NULL) + call tbcfnd (pst, "YCENTER", column[3], 1) + if (column[3] == NULL) + call error (0, "Error reading YCENTER column from PSF star file\n") + + call tbcfnd (pst, MAG, column[4], 1) + if (column[4] == NULL) + call tbcfnd (pst, APMAG, column[4], 1) + if (column[4] == NULL) + call error (0, "Error reading MAG column from PSF star file\n") +end + + +# DP_TPTREAD -- Read the id from an ST table. + +procedure dp_tptread (pst, column, idno, x, y, mag, rowno) + +pointer pst # pointer to the ST table +int column[ARB] # array of column pointers +int idno # the output id number +real x, y # the output x and y position +real mag # the output magnitude +int rowno # the row number + +bool nullflag + +begin + call tbrgti (pst, column[1], idno, nullflag, 1, rowno) + if (nullflag) + idno = 0 + call tbrgtr (pst, column[2], x, nullflag, 1, rowno) + call tbrgtr (pst, column[3], y, nullflag, 1, rowno) + call tbrgtr (pst, column[4], mag, nullflag, 1, rowno) +end diff --git a/noao/digiphot/daophot/psf/dpshowpsf.x b/noao/digiphot/daophot/psf/dpshowpsf.x new file mode 100644 index 00000000..ec5a7766 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpshowpsf.x @@ -0,0 +1,287 @@ +include <mach.h> +include <ctype.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +define HELPFILE "daophot$psf/showpsf.key" + +# DP_SHOWPSF -- Interactively make surface and/or contour plots of the +# data subraster around the current PSF star. + +procedure dp_showpsf (dao, im, subrast, ncols, nlines, x1, y1, gd, star_ok) + +pointer dao # pointer to DAOPHOT structure +pointer im # the input image descriptor +real subrast[ncols,nlines] # image subraster +int ncols, nlines # dimensions of the subraster +int x1, y1 # coordinates of left hand corner +pointer gd # pointer to the graphics stream +bool star_ok # true if PSF star ok + +real hibad, wx, wy, rval +pointer sp, cmd, title, psf +int wcs, key, ip +bool do_replot + +int clgcur(), ctor() + +begin + # Return if the graphics stream is undefined. + if (gd == NULL) + return + + # Allocate working space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + + # Initialize various daophot pointers. + psf = DP_PSF (dao) + + # Initialize the contour plot parameters. + DP_CCEILING (psf) = 0.0 + DP_CFLOOR (psf) = 0.0 + + # Define the hibad parameter. + hibad = DP_MAXGDATA(dao) + if (IS_INDEFR(hibad)) + hibad = MAX_REAL + + # Create the plot title. + call dp_ltov (im, DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), wx, wy, 1) + call sprintf (Memc[title], SZ_LINE, "Star: %d X: %g Y: %g Mag: %g\n") + call pargi (DP_CUR_PSFID(psf)) + call pargr (wx) + call pargr (wy) + call pargr (DP_CUR_PSFMAG(psf)) + + # Initialize plot. + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) + call dp_surfpsf (dao, subrast, ncols, nlines, Memc[title], gd) + else if (DP_PLOTTYPE(psf) == PSF_CONTOURPLOT) + call dp_contpsf (dao, subrast, ncols, nlines, Memc[title], gd) + else if (DP_PLOTTYPE(psf) == PSF_RADIALPLOT) + call dp_radpsf (dao, subrast, ncols, nlines, x1, y1, Memc[title], + gd) + + if (DP_CUR_PSFMAX(psf) > hibad) + call printf ("Warning: Star is probably saturated\n") + do_replot = false + + while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + switch (key) { + + # Print the help page + case '?': + if (gd == NULL) + call pagefile (HELPFILE, "") + else + call gpagefile (gd, HELPFILE, "") + + # Print out the photometry for the star. + case 'p': + call printf ("Star: %d Mag: %7.2f Coords: %7.2f %7.2f\n") + call pargi (DP_CUR_PSFID(psf)) + call pargr (DP_CUR_PSFMAG(psf)) + call pargr (DP_CUR_PSFX(psf)) + call pargr (DP_CUR_PSFY(psf)) + + # Print out the plot parameters. + case 't': + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) { + call printf ("Surface plot: angv = %6.1f angh = %6.1f ") + call pargr (DP_MANGV(psf)) + call pargr (DP_MANGH(psf)) + call printf (" Minimum: %7.1f Maximum: %7.1f\n") + call pargr (DP_CUR_PSFMIN(psf)) + call pargr (DP_CUR_PSFMAX(psf)) + } else if (DP_PLOTTYPE(psf) == PSF_CONTOURPLOT) { + call printf ("Contour plot: floor = %6.1f ceil = %6.1f ") + call pargr (DP_CFLOOR(psf)) + call pargr (DP_CCEILING(psf)) + call printf (" Minimum: %7.1f Maximum: %7.1f\n") + call pargr (DP_CUR_PSFMIN(psf)) + call pargr (DP_CUR_PSFMAX(psf)) + } else if (DP_PLOTTYPE(psf) == PSF_RADIALPLOT) { + call printf ( + "Profile plot: Minimum: %7.1f Maximum: %7.1f\n") + call pargr (DP_CUR_PSFMIN(psf)) + call pargr (DP_CUR_PSFMAX(psf)) + } else + call printf ("Unknown plot type.\n") + + # Accept star and quit cursor loop. + case 'a': + star_ok = true + break + + # Star rejected + case 'd': + star_ok = false + break + + # Increase vertical viewing angle of mesh plot by 15 degrees. + case 'n': + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) { + DP_MANGV (psf) = DP_MANGV (psf) + 15. + do_replot = true + } else + call printf ("This key only active for mesh plots.\n") + + # Decrease vertical viewing angle of mesh plot by 15 degrees. + case 's': + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) { + DP_MANGV (psf) = DP_MANGV (psf) - 15. + do_replot = true + } else + call printf ("This key only active for mesh plots.\n") + + # Decrease horizontal viewing angle of mesh plot by 15 degrees. + case 'w': + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) { + DP_MANGH (psf) = DP_MANGH (psf) - 15. + do_replot = true + } else + call printf ("This key only active for mesh plots.\n") + + # Increase horizontal viewing angle of mesh plot by 15 degrees. + case 'e': + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) { + DP_MANGH (psf) = DP_MANGH (psf) + 15. + do_replot = true + } else + call printf ("This key only active for mesh plots.\n") + + # Plot the default mesh plot. + case 'm': + DP_PLOTTYPE(psf) = PSF_MESHPLOT + DP_MANGV (psf) = 30. + DP_MANGH (psf) = -30. + DP_MFLOOR(psf) = 0.0 + DP_MCEILING(psf) = 0.0 + do_replot = true + + # Plot the default contour plots. + case 'c': + DP_PLOTTYPE(psf) = PSF_CONTOURPLOT + DP_CFLOOR(psf) = 0.0 + DP_CCEILING(psf) = 0.0 + do_replot = true + + # Plot the radial profile plots. + case 'r': + DP_PLOTTYPE(psf) = PSF_RADIALPLOT + do_replot = true + + # Command mode. + case ':': + for (ip=1; IS_WHITE (Memc[cmd+ip-1]); ip=ip+1) + ; + + switch (Memc[cmd+ip-1]) { + + # Set surface plot and angles + case 'm': + DP_PLOTTYPE(psf) = PSF_MESHPLOT + ip = ip + 1 + if (ctor (Memc[cmd], ip, DP_MANGV(psf)) <= 0) + DP_MANGV(psf) = 30. + if (ctor (Memc[cmd], ip, DP_MANGH(psf)) <= 0) + DP_MANGH(psf) = -30. + do_replot = true + + case 'c': + # Set surface contour and levels + DP_PLOTTYPE(psf) = PSF_CONTOURPLOT + ip = ip + 1 + if (ctor (Memc[cmd], ip, DP_CFLOOR(psf)) <= 0) + DP_CFLOOR(psf) = 0.0 + if (ctor (Memc[cmd], ip, DP_CCEILING(psf)) <= 0) + DP_CCEILING(psf) = 0.0 + do_replot = true + + case 'r': + # PLot radial profile. + DP_PLOTTYPE(psf) = PSF_RADIALPLOT + do_replot = true + + # Set vertical angle + case 'v': + ip = ip + 1 + if (ctor (Memc[cmd], ip, rval) <= 0) { + call printf ( + "Surface plot vertical viewing angle: %g\n") + call pargr (DP_MANGV(psf)) + } else { + DP_MANGV(psf) = rval + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) + do_replot = true + } + + # Set horizontal angle + case 'h': + ip = ip + 1 + if (ctor (Memc[cmd], ip, rval) <= 0) { + call printf ( + "Surface plot horizontal viewing angle: %g\n") + call pargr (DP_MANGH(psf)) + } else { + DP_MANGH(psf) = rval + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) + do_replot = true + } + + # Set the floor value. + case 'l': + ip = ip + 1 + if (ctor (Memc[cmd], ip, rval) <= 0) { + call printf ( + "Contour plot floor value: %g\n") + call pargr (DP_CFLOOR(psf)) + } else { + DP_CFLOOR(psf) = rval + if (DP_PLOTTYPE(psf) == PSF_CONTOURPLOT) + do_replot = true + } + + # Set the ceiling value. + case 'u': + ip = ip + 1 + if (ctor (Memc[cmd], ip, rval) <= 0) { + call printf ( + "Contour plot ceiling value: %g\n") + call pargr (DP_CCEILING(psf)) + } else { + DP_CCEILING(psf) = rval + if (DP_PLOTTYPE(psf) == PSF_CONTOURPLOT) + do_replot = true + } + + default: + call printf ("Unknown keystroke or cursor command.\007\n") + } + + default: + call printf ("Unknown keystroke or cursor command.\007\n") + } + + # Replot the data. + if (do_replot) { + if (DP_PLOTTYPE(psf) == PSF_MESHPLOT) + call dp_surfpsf (dao, subrast, ncols, nlines, Memc[title], + gd) + else if (DP_PLOTTYPE(psf) == PSF_CONTOURPLOT) + call dp_contpsf (dao, subrast, ncols, nlines, Memc[title], + gd) + else if (DP_PLOTTYPE(psf) == PSF_RADIALPLOT) + call dp_radpsf (dao, subrast, ncols, nlines, x1, y1, + Memc[title], gd) + do_replot = false + } + } + + call gdeactivate (gd, 0) + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/dpspstars.x b/noao/digiphot/daophot/psf/dpspstars.x new file mode 100644 index 00000000..69b14933 --- /dev/null +++ b/noao/digiphot/daophot/psf/dpspstars.x @@ -0,0 +1,194 @@ +include <mach.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +define NCOLUMN 5 + +# DP_GPSTARS -- Select the psf stars. + +procedure dp_gpstars (dao, im, tp_in, tp_out, text_file, maxnpsf, gd, mgd, id, + mkstars, interactive, use_cmdfile) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +int tp_in # the input file descriptor +int tp_out # the output file descriptor +bool text_file # text or table file +int maxnpsf # the maximum number of psf stars +pointer gd # pointer to the graphics stream +pointer mgd # pointer to the plot metacode file +pointer id # pointer to the image display stream +bool mkstars # mark the accepted and deleted stars +bool interactive # interactive mode +bool use_cmdfile # cursor command file mode + +int ier, npsf +pointer apsel, sp, ocolpoint, index +real radius +int dp_pfstars(), dp_ipfstars() + +begin + # Get some daophot pointers. + apsel = DP_APSEL(dao) + + # Define the radius for determining proximity. + radius = DP_PSFRAD(dao) + DP_FITRAD(dao) + 2.0 + + # Allocate some working memory. + call smark (sp) + call salloc (ocolpoint, NCOLUMN, TY_POINTER) + call salloc (index, DP_APNUM(apsel), TY_INT) + + # Initialize the output file. + if (text_file) { + call seek (tp_in, BOF) + call dp_apheader (tp_in, tp_out) + call dp_xpselpars (tp_out, DP_INIMAGE(dao), maxnpsf, DP_SCALE(dao), + DP_SPSFRAD(dao), DP_SFITRAD(dao)) + call dp_xpbanner (tp_out) + } else { + call dp_tpdefcol (tp_out, Memi[ocolpoint]) + call tbhcal (tp_in, tp_out) + call dp_tpselpars (tp_out, DP_INIMAGE(dao), maxnpsf, DP_SCALE(dao), + DP_SPSFRAD(dao), DP_SFITRAD(dao)) + } + + # Change all the INDEFS to a large negative number. + call dp_chindef (Memr[DP_APMAG(apsel)], Memr[DP_APMAG(apsel)], + DP_APNUM(apsel), -MAX_REAL) + + # Sort the stars. + call quick (Memr[DP_APMAG(apsel)], DP_APNUM(apsel), Memi[index], ier) + call dp_irectify (Memi[DP_APID(apsel)], Memi[index], DP_APNUM(apsel)) + call dp_rectify (Memr[DP_APXCEN(apsel)], Memi[index], DP_APNUM(apsel)) + call dp_rectify (Memr[DP_APYCEN(apsel)], Memi[index], DP_APNUM(apsel)) + call dp_rectify (Memr[DP_APMSKY(apsel)], Memi[index], DP_APNUM(apsel)) + + # Select the stars and write them to the output file. + if (use_cmdfile) + npsf = dp_ipfstars (dao, im, maxnpsf, -MAX_REAL, radius, gd, + mgd, id, mkstars, false, false) + else if (interactive) + npsf = dp_ipfstars (dao, im, maxnpsf, -MAX_REAL, radius, gd, + mgd, id, mkstars, true, true) + else + npsf = dp_pfstars (dao, im, Memi[DP_APID(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APMAG(apsel)], Memr[DP_APMSKY(apsel)], + DP_APNUM(apsel), maxnpsf, -MAX_REAL, radius, mgd) + + if (DP_VERBOSE(dao) == YES) { + call printf ("\nTotal of %d PSF stars selected\n") + call pargi (npsf) + } + + # Write out the stars. + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], + Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)], npsf) + call dp_wpstars (tp_out, Memi[ocolpoint], text_file, + Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)], + Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)], + Memr[DP_APMSKY(apsel)], npsf) + + # Free memory. + call sfree (sp) +end + + +# DP_CHINDEF -- Change all the INDEFS to a legal number which is a lower limit. + +procedure dp_chindef (a, b, npix, lolimit) + +real a[ARB] # the input array. +real b[ARB] # the output array. +int npix # number of points +real lolimit # the lower limit + +int i + +begin + do i = 1, npix { + if (IS_INDEFR(a[i])) + b[i] = lolimit + else + b[i] = a[i] + } +end + + +# DP_PFSTARS -- Select the psf stars. + +int procedure dp_pfstars (dao, im, ids, xcen, ycen, mag, sky, nstar, maxnpsf, + lolimit, radius, mgd) + +pointer dao # pointer to the daophot structure +pointer im # pointer to the input image +int ids[ARB] # array of star ids +real xcen[ARB] # array of x coordinates +real ycen[ARB] # array of y coordinates +real mag[ARB] # array of magnitudes +real sky[ARB] # array of sky values +int nstar # the number of stars +int maxnpsf # the maximum number of psf stars +real lolimit # lower data limit +real radius # minimum separation +pointer mgd # pointer to the plot metacode file + +bool omit +int istar, jstar, npsf +real radsq, dy2, dr2 +int dp_addstar() + +begin + # Initialize the list. + call dp_pseti (dao, PNUM, 0) + + # Set the radius. + radsq = radius * radius + + # Get the first star. + if ((mag[1] > lolimit) && (dp_addstar (dao, im, xcen[1], ycen[1], + INDEFR, ids[1], NULL, mgd, false)) == OK) { + npsf = 1 + } else + npsf = 0 + + # Loop over the candidate psf stars. + do istar = 2, nstar { + + # Test for the maximum number of psf stars. + if (npsf >= maxnpsf) + break + + # Test that the candidate psf stars are not saturated and + # are sufficiently far from the edge of the frame. + + if (mag[istar] <= lolimit) + next + + # Text that there are no brighter stars with a distance squared + # of radsq. + omit = false + do jstar = 1, istar - 1 { + dy2 = abs (ycen[jstar] - ycen[istar]) + if (dy2 >= radius) + next + dr2 = (xcen[jstar] - xcen[istar]) ** 2 + dy2 ** 2 + if (dr2 >= radsq) + next + omit = true + break + } + + if (omit) + next + + if (dp_addstar (dao, im, xcen[istar], ycen[istar], INDEFR, + ids[istar], NULL, mgd, false) == ERR) + next + npsf = npsf + 1 + } + + return (npsf) +end 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 diff --git a/noao/digiphot/daophot/psf/dpsurfpsf.x b/noao/digiphot/daophot/psf/dpsurfpsf.x new file mode 100644 index 00000000..de4e9d6e --- /dev/null +++ b/noao/digiphot/daophot/psf/dpsurfpsf.x @@ -0,0 +1,437 @@ +include <fset.h> +include <mach.h> +include <error.h> +include <gset.h> +include <config.h> +include <xwhen.h> +include "../lib/daophotdef.h" +include "../lib/psfdef.h" + +define DUMMY 6 + +# DP_SURFPSF -- Draw a perspective view of a subraster with a given altitude +# and azimuth of the viewing angle. Floor and ceiling constraints may be +# applied to the data before plotting. + +procedure dp_surfpsf (dao, subras, ncols, nlines, title, gd) + +pointer dao # pointer to DAOPHOT structure +real subras[ncols,nlines] # pointer to subraster +int ncols, nlines # dimensions of subraster +char title[ARB] # title string +pointer gd # pointer to graphics stream + +char sysidstr[SZ_LINE] +int first, wkid, epa, status, old_onint, tsujmp[LEN_JUMPBUF] +pointer sp, temp, work, psf +real angh, angv, imcols, imlines, floor, ceiling, vpx1, vpx2, vpy1, vpy2 + +extern dp_sonint() +common /tsucom/ tsujmp +common /noaovp/ vpx1, vpx2, vpy1, vpy2 +common /frstfg/ first + +begin + # Get the psf fitting substructure pointer. + psf = DP_PSF(dao) + + # Initialize surface common blocks before changing any parameters. + first = 1 + call srfabd () + + # Set local variables. + angh = DP_MANGH (psf) + angv = DP_MANGV (psf) + floor = DP_MFLOOR (psf) + ceiling = DP_MCEILING (psf) + floor = min (floor, ceiling) + ceiling = max (floor, ceiling) + + # Allow room for axes and labels. + vpx1 = 0.10 + vpx2 = 0.90 + vpy1 = 0.10 + vpy2 = 0.90 + + # Make a copy of the subraster so we can subtract the zero level. + imcols = real (ncols) + imlines = real (nlines) + call smark (sp) + call salloc (temp, ncols * nlines, TY_REAL) + call amovr (subras, Memr[temp], nlines * ncols) + + # Allocate the working storage needed by EZSRFC. + call malloc (work, 2 * (2 * ncols * nlines + ncols + nlines), TY_REAL) + + # Take off floor and ceiling if enabled (nonzero). + call dp_slimits (Memr[temp], ncols, nlines, floor, ceiling) + + # Set up the titles and the viewport. + call gopks (STDERR) + wkid = 1 + call gclear (gd) + call gopwk (wkid, DUMMY, gd) + call gacwk (wkid) + call gtext (gd, 0.5, .96, title, "s=0.8;f=b;h=c") + call sysid (sysidstr, SZ_LINE) + call gtext (gd, 0.5, .04, sysidstr, "h=c;v=b;s=.5") + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + + # Install interrupt exception handler. + call zlocpr (dp_sonint, epa) + call xwhen (X_INT, epa, old_onint) + + # Plot the surface. + call zsvjmp (tsujmp, status) + if (status == OK) + call ezsrfc (Memr[temp], ncols, nlines, angh, angv, Memr[work]) + else { + call gcancel (gd) + call fseti (STDOUT, F_CANCEL, OK) + } + + # Draw the perimeter. + call gswind (gd, 1.0, imcols, 1.0, imlines) + call gseti (gd, G_CLIP, NO) + call dp_sperimeter (gd, Memr[temp], ncols, nlines, angh, angv) + + # Clean up. + call gdawk (wkid) + call gclks () + call mfree (work, TY_REAL) + call sfree (sp) +end + + +# DP_SONINT -- Interrupt handler for the task surface. Branches back to ZSVJMP +# in the main routine to permit shutdown without an error message. + +procedure dp_sonint (vex, next_handler) + +int vex # virtual exception +int next_handler # not used + +int tsujmp[LEN_JUMPBUF] +common /tsucom/ tsujmp + +begin + call xer_reset() + call zdojmp (tsujmp, vex) +end + + +# DP_SLIMITS -- Apply the floor and ceiling constraints to the subraster. +# If either value is exactly zero, it is not applied. + +procedure dp_slimits (ras, m, n, floor, ceiling) + +real ras[m,n] +int m, n +real floor, ceiling +int i + +begin + do i = 1, n { + if (floor != 0) + call amaxkr (ras[1,i], floor, ras[1,i], m) + if (ceiling != 0) + call aminkr (ras[1,i], ceiling, ras[1,i], m) + } +end + + +define SZ_TLABEL 10 + +# DP_SPERIMETER -- draw and label axes around the surface plot. + +procedure dp_sperimeter (gp, z, ncols, nlines, angh, angv) + +pointer gp # Graphics pointer +int ncols # Number of image columns +int nlines # Number of image lines +real z[ncols, nlines] # Array of intensity values +real angh # Angle of horizontal inclination +real angv # Angle of vertical inclination + +pointer sp, x_val, y_val, kvec +char tlabel[SZ_TLABEL] +real xmin, ymin, delta, fact1, flo, hi, xcen, ycen +real x1_perim, x2_perim, y1_perim, y2_perim, z1, z2 +real wc1, wc2, wl1, wl2, del +int i, j +int itoc() +data fact1 /2.0/ +real vpx1, vpx2, vpy1, vpy2 +common /noaovp/ vpx1, vpx2, vpy1, vpy2 + +begin + call smark (sp) + call salloc (x_val, ncols + 2, TY_REAL) + call salloc (y_val, nlines + 2, TY_REAL) + call salloc (kvec, max (ncols, nlines) + 2, TY_REAL) + + # Get window coordinates set up in calling procedure. + call ggwind (gp, wc1, wc2, wl1, wl2) + + # Set up window, viewport for output. The coordinates returned + # from trn32s are in the range [1-1024]. + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + + # Find range of z for determining perspective. + flo = MAX_REAL + hi = -flo + do j = 1, nlines { + call alimr (z[1,j], ncols, z1, z2) + flo = min (flo, z1) + hi = max (hi, z2) + } + + # Set up linear endpoints and spacing as used in surface. + + delta = (hi-flo) / (max (ncols,nlines) -1.) * fact1 + xmin = -(real (ncols/2) * delta + real (mod (ncols+1, 2)) * delta) + ymin = -(real (nlines/2) * delta + real (mod (nlines+1, 2)) * delta) + del = 2.0 * delta + + # The perimeter is separated from the surface plot by the + # width of delta. + + x1_perim = xmin - delta + y1_perim = ymin - delta + x2_perim = xmin + (real (ncols) * delta) + y2_perim = ymin + (real (nlines) * delta) + + # Set up linear arrays over full perimeter range + do i = 1, ncols + 2 + Memr[x_val+i-1] = x1_perim + (i-1) * delta + do i = 1, nlines + 2 + Memr[y_val+i-1] = y1_perim + (i-1) * delta + + # Draw and label axes and tick marks. + # It is important that frame has not been called after calling srface. + # First to draw the perimeter. Which axes get drawn depends on the + # values of angh and angv. Get angles in the range [-180, 180]. + + if (angh > 180.) + angh = angh - 360. + else if (angh < -180.) + angh = angh + 360. + + if (angv > 180.) + angv = angv - 360. + else if (angv < -180.) + angv = angv + 360. + + # Calculate positions for the axis labels + xcen = 0.5 * (x1_perim + x2_perim) + ycen = 0.5 * (y1_perim + y2_perim) + + if (angh >= 0) { + if (angv >= 0) { + # Case 1: xy rotation positive, looking down from above mid z. + + # First draw x axis. + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call dp_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call dp_label_axis (xcen, y2_perim+del, flo, "X-AXIS", -1, -2) + call dp_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + call dp_label_axis (xmin, y2_perim+del, flo, "1", -1, -2) + if (itoc (int (wc2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, -1, -2) + + # Now draw y axis. + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call dp_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call dp_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call dp_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + call dp_label_axis (x2_perim+del, ymin, flo, "1", 2, -1) + if (itoc (int (wl2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + + } else { + # Case 2: xy rotation positive, looking up from below mid z. + + # First draw x axis. + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call dp_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call dp_label_axis (xcen, y1_perim-del, flo, "X-AXIS", -1, 2) + call dp_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + call dp_label_axis (xmin, y1_perim-del, flo, "1", -1, 2) + if (itoc (int (wc2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, -1, 2) + + # Now draw y axis. + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call dp_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call dp_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call dp_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + call dp_label_axis (x1_perim-del, ymin, flo, "1", 2, 1) + if (itoc (int (wl2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + if (angh < 0) { + if (angv > 0) { + + # Case 3: xy rotation negative, looking down from above mid z + # (default). First draw x axis. + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call dp_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call dp_label_axis (xcen, y1_perim-del, flo, "X-AXIS", 1, 2) + call dp_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + call dp_label_axis (xmin, y1_perim-del, flo, "1", 1, 2) + if (itoc (int (wc2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, 1, 2) + + # Now draw y axis. + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call dp_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call dp_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call dp_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + call dp_label_axis (x2_perim+del, ymin, flo, "1", 2, -1) + if (itoc (int (wl2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + } else { + + # Case 4: xy rotation negative, looking up from below mid z. + # First draw x axis. + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call dp_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call dp_label_axis (xcen, y2_perim+del, flo, "X-AXIS", 1, -2) + call dp_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + call dp_label_axis (xmin, y2_perim+del, flo, "1", 1, -2) + if (itoc (int (wc2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, 1, -2) + + # Now draw y axis. + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call dp_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call dp_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call dp_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + call dp_label_axis (x1_perim-del, ymin, flo, "1", 2, 1) + if (itoc (int (wl2), tlabel, SZ_TLABEL) <= 0) + tlabel[1] = EOS + call dp_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + # Flush plotit buffer before returning. + call plotit (0, 0, 2) + call sfree (sp) +end + + +# DP_DRAW_AXIS -- Draw the axes around the surface plot. + +procedure dp_draw_axis (xvals, yvals, zval, nvals) + +int nvals +real xvals[nvals] +real yvals[nvals] +real zval +pointer sp, xt, yt +int i +real dum + +begin + call smark (sp) + call salloc (xt, nvals, TY_REAL) + call salloc (yt, nvals, TY_REAL) + + do i = 1, nvals + call trn32s (xvals[i], yvals[i], zval, Memr[xt+i-1], Memr[yt+i-1], + dum, 1) + call gpl (nvals, Memr[xt], Memr[yt]) + call sfree (sp) +end + + +define CSIZE 24 + + +# DP_LABEL_AXIS -- Draw the axes labels. + +procedure dp_label_axis (xval, yval, zval, sppstr, path, up) + +real xval +real yval +real zval +char sppstr[SZ_LINE] +int path +int up + +int nchars +int strlen() +% character*64 fstr + +begin + nchars = strlen (sppstr) +% call f77pak (sppstr, fstr, 64) + call pwrzs (xval, yval, zval, fstr, nchars, CSIZE, path, up, 0) +end + + +# DP_DRAW_TICKS -- Draw the x tick marks. + +procedure dp_draw_ticksx (x, y1, y2, zval, nvals) + +int nvals +real x[nvals] +real y1, y2 +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x[i], y1, zval, tkx[1], tky[1], dum, 1) + call trn32s (x[i], y2, zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end + + +# DP_DRAW_TICKSY -- Draw the y tick marks. + +procedure dp_draw_ticksy (x1, x2, y, zval, nvals) + +int nvals +real x1, x2 +real y[nvals] +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x1, y[i], zval, tkx[1], tky[1], dum, 1) + call trn32s (x2, y[i], zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end diff --git a/noao/digiphot/daophot/psf/dpwritepsf.x b/noao/digiphot/daophot/psf/dpwritepsf.x new file mode 100644 index 00000000..239e4faf --- /dev/null +++ b/noao/digiphot/daophot/psf/dpwritepsf.x @@ -0,0 +1,270 @@ +include <time.h> +include <imhdr.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +# DP_WRITEPSF -- Write out the PSF into an IRAF image. + +procedure dp_writepsf (dao, im, psfim) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +pointer psfim # pointer to the output psf image + +begin + # Check that the psfimage is open. + if (psfim == NULL) + return + + # Write out the id and fitting parameters. + call dp_widpars (dao, psfim) + + # Write out the psf function definition parameters. + call dp_wfuncpars (dao, psfim) + + # Write out the list of PSF stars. + call dp_wstars (dao, im, psfim) + + # Write out the lookup table. + call dp_wlt (dao, psfim) +end + + +# DP_WIDPARS -- Add the id and fitting parameters to the PSF image header + +procedure dp_widpars (dao, psfim) + +pointer dao # pointer to the daophot structure +pointer psfim # the psf image descriptor + +pointer sp, outstr, date, time +bool itob() +int envfind() + +begin + # Allocate working space. + call smark (sp) + call salloc (outstr, SZ_LINE, TY_CHAR) + call salloc (date, SZ_DATE, TY_CHAR) + call salloc (time, SZ_DATE, TY_CHAR) + + # Record IRAF version, user, host, date, time, package and task. + if (envfind ("version", Memc[outstr], SZ_LINE) <= 0) + call strcpy ("IRAF", Memc[outstr], SZ_LINE) + call dp_rmwhite (Memc[outstr], Memc[outstr], SZ_LINE) + call imastr (psfim, "IRAF", Memc[outstr]) + call gethost (Memc[outstr], SZ_LINE) + call imastr (psfim, "HOST", Memc[outstr]) + if (envfind ("userid", Memc[outstr], SZ_LINE) <= 0) + Memc[outstr] = EOS + call imastr (psfim, "USER", Memc[outstr]) + call dp_date (Memc[date], Memc[time], SZ_DATE) + call imastr (psfim, "DATE", Memc[date]) + call imastr (psfim, "TIME", Memc[time]) + + # Write out the package, task, and input/output file names. + call imastr (psfim, "PACKAGE", "daophot") + call imastr (psfim, "TASK", "psf") + call dp_imroot (DP_INIMAGE(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "IMAGE", Memc[outstr]) + call dp_froot (DP_INPHOTFILE(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "PHOTFILE", Memc[outstr]) + call dp_froot (DP_COORDS(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "PSTFILE", Memc[outstr]) + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "PSFIMAGE", Memc[outstr]) + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "OPSTFILE", Memc[outstr]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call imastr (psfim, "GRPSFILE", Memc[outstr]) + + # Add information about fitting parameters. + call imaddr (psfim, "SCALE", DP_SCALE(dao)) + call imaddr (psfim, "PSFRAD", DP_SPSFRAD (dao)) + call imaddr (psfim, "FITRAD", DP_SFITRAD(dao)) + call imaddr (psfim, "DATAMIN", DP_MINGDATA(dao)) + call imaddr (psfim, "DATAMAX", DP_MAXGDATA(dao)) + call imaddi (psfim, "NCLEAN", DP_NCLEAN(dao)) + call imaddb (psfim, "USESAT", itob (DP_SATURATED(dao))) + + # Define the image title. + call sprintf (IM_TITLE(psfim), SZ_IMTITLE, "PSF for image: %s") + call pargstr (DP_INIMAGE(dao)) + + call sfree (sp) +end + + +# DP_WFUNCPARS -- Write out the the parameters of the PSF function +# to the PSF image. + +procedure dp_wfuncpars (dao, psfim) + +pointer dao # pointer to the daophot structure +pointer psfim # image descriptor + +int i +pointer sp, str, psffit +bool itob() + +begin + psffit = DP_PSFFIT(dao) + + call imastr (psfim, "FUNCTION", DP_FUNCTION(dao)) + call imaddr (psfim, "PSFX", DP_PSFX(psffit)) + call imaddr (psfim, "PSFY", DP_PSFY(psffit)) + call imaddr (psfim, "PSFHEIGHT", DP_PSFHEIGHT(psffit)) + call imaddr (psfim, "PSFMAG", DP_PSFMAG (psffit)) + + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + switch (DP_PSFUNCTION(psffit)) { + case FCTN_MOFFAT25: + call imaddi (psfim, "NPARS", DP_PSFNPARS(psffit)+1) + do i = 1, DP_PSFNPARS(psffit) { + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (i) + call imaddr (psfim, Memc[str], Memr[DP_PSFPARS(psffit)+i-1]) + } + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (DP_PSFNPARS(psffit)+1) + call imaddr (psfim, Memc[str], 2.5) + case FCTN_MOFFAT15: + call imaddi (psfim, "NPARS", DP_PSFNPARS(psffit)+1) + do i = 1, DP_PSFNPARS(psffit) { + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (i) + call imaddr (psfim, Memc[str], Memr[DP_PSFPARS(psffit)+i-1]) + } + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (DP_PSFNPARS(psffit)+1) + call imaddr (psfim, Memc[str], 1.5) + default: + call imaddi (psfim, "NPARS", DP_PSFNPARS(psffit)) + do i = 1, DP_PSFNPARS(psffit) { + call sprintf (Memc[str], SZ_FNAME, "PAR%d") + call pargi (i) + call imaddr (psfim, Memc[str], Memr[DP_PSFPARS(psffit)+i-1]) + } + } + + call imaddi (psfim, "VARORDER", DP_VARORDER(dao)) + call imaddb (psfim, "FEXPAND", itob (DP_FEXPAND(dao))) + + call sfree (sp) +end + + +# DP_WSTARS -- Write out the PSF star list to the PSF image. + +procedure dp_wstars (dao, im, psfim) + +pointer dao # pointer to the daophot descriptor +pointer im # the input image descriptor +pointer psfim # the psfimage descriptor + +real tx, ty +pointer apsel, psf, sp, str +int i + +begin + apsel = DP_APSEL(dao) + psf = DP_PSF(dao) + + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Write out the number of PSF stars. + call sprintf (Memc[str], SZ_FNAME, "NPSFSTAR") + call imaddi (psfim, Memc[str], DP_PNUM(psf)) + + # Write out the ids of all the PSF stars. + do i = 1, DP_PNUM(psf) { + + call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+i-1], + Memr[DP_APYCEN(apsel)+i-1], tx, ty, 1) + call sprintf (Memc[str], SZ_FNAME, "ID%d") + call pargi (i) + call imaddi (psfim, Memc[str], Memi[DP_APID(apsel)+i-1]) + + call sprintf (Memc[str], SZ_FNAME, "X%d") + call pargi (i) + call imaddr (psfim, Memc[str], tx) + + call sprintf (Memc[str], SZ_FNAME, "Y%d") + call pargi (i) + call imaddr (psfim, Memc[str], ty) + + call sprintf (Memc[str], SZ_FNAME, "MAG%d") + call pargi (i) + call imaddr (psfim, Memc[str], Memr[DP_APMAG(apsel)+i-1]) + } + + call sfree (sp) +end + + +# DP_WLT -- Write out the PSF lookup table to the output PSF image. + +procedure dp_wlt (dao, psfim) + +pointer dao # pointer to DAO Structure +pointer psfim # image descriptor + +int nexp +pointer psffit + +begin + psffit = DP_PSFFIT(dao) + nexp = DP_NVLTABLE(psffit) + DP_NFEXTABLE(psffit) + + IM_PIXTYPE(psfim) = TY_REAL + if (nexp == 0) { + IM_NDIM(psfim) = 0 + } else if (nexp == 1) { + IM_NDIM(psfim) = 2 + IM_LEN(psfim,1) = DP_PSFSIZE(psffit) + IM_LEN(psfim,2) = DP_PSFSIZE(psffit) + } else { + IM_NDIM(psfim) = 3 + IM_LEN(psfim,1) = DP_PSFSIZE(psffit) + IM_LEN(psfim,2) = DP_PSFSIZE(psffit) + IM_LEN(psfim,3) = nexp + } + + if (nexp > 0) + call dp_wltim (psfim, Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_PSFSIZE(psffit), nexp) +end + + +# DP_WLTIM -- Write the lookup tables into the image pixels. + +procedure dp_wltim (psfim, psflut, nxlut, nylut, nexp) + +pointer psfim # image descriptor +real psflut[nexp,nxlut,ARB] # the psf lookup table +int nxlut, nylut,nexp # the dimensions of the psf look-up table + +int i, j, k +pointer sp, v, buf +int impnlr() + +begin + call smark (sp) + call salloc (v, IM_MAXDIM, TY_LONG) + + call amovkl (long(1), Meml[v], IM_MAXDIM) + do k = 1, nexp { + do j = 1, nylut { + if (impnlr (psfim, buf, Meml[v]) == EOF) + ; + do i = 1, nxlut + Memr[buf+i-1] = psflut[k,i,j] + } + } + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/mkpkg b/noao/digiphot/daophot/psf/mkpkg new file mode 100644 index 00000000..5dd38047 --- /dev/null +++ b/noao/digiphot/daophot/psf/mkpkg @@ -0,0 +1,70 @@ +# PSF task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + dpaddstar.x <imhdr.h> <mach.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpcontpsf.x <error.h> <mach.h> \ + <gset.h> <config.h> \ + <xwhen.h> <fset.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpdelstar.x <imhdr.h> ../lib/daophotdef.h \ + ../lib/psfdef.h <mach.h> + dpfitpsf.x <mach.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/psfdef.h ../lib/peakdef.h + dplocstar.x <mach.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/psfdef.h + dpispstars.x <fset.h> <ctype.h> \ + <gset.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/psfdef.h + dppcolon.x ../lib/daophotdef.h ../lib/psfdef.h + dpqverify.x <fset.h> + dpmempsf.x ../lib/daophotdef.h ../lib/psfdef.h + dpmkpsf.x <imhdr.h> <ctype.h> \ + <gset.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/psfdef.h + dppconfirm.x + dppsfutil.x ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/psfdef.h + dpplotpsf.x ../lib/daophotdef.h ../lib/psfdef.h + dppset.x ../lib/daophotdef.h ../lib/psfdef.h + dppstat.x ../lib/daophotdef.h ../lib/psfdef.h + dppsubrast.x <mach.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpptconfirm.x + dppwrtgrp.x <time.h> <tbset.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/psfdef.h + dppwselmer.x <tbset.h> ../lib/apseldef.h + dpradpsf.x <gset.h> ../lib/daophotdef.h \ + ../lib/psfdef.h + dprmpsf.x <fset.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dprstars.x <gset.h> <tbset.h> \ + ../lib/apseldef.h ../lib/psfdef.h + dpshowpsf.x <mach.h> <ctype.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpspstars.x <mach.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/psfdef.h + dpsubpsf.x <mach.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpsurfpsf.x <error.h> <mach.h> \ + <gset.h> <config.h> \ + <xwhen.h> <fset.h> \ + ../lib/daophotdef.h ../lib/psfdef.h + dpwritepsf.x <time.h> <imhdr.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/psfdef.h + t_psf.x <fset.h> <gset.h> \ + <imhdr.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/psfdef.h + t_pstselect.x <fset.h> <gset.h> \ + <imhdr.h> ../lib/daophotdef.h \ + ../lib/apseldef.h ../lib/psfdef.h + ; diff --git a/noao/digiphot/daophot/psf/mkpsf.key b/noao/digiphot/daophot/psf/mkpsf.key new file mode 100644 index 00000000..ed511025 --- /dev/null +++ b/noao/digiphot/daophot/psf/mkpsf.key @@ -0,0 +1,40 @@ + Keystroke Commands + +? Print help +p Print photometry for star nearest the cursor +l List the current psf stars +a Add star nearest cursor to psf star list +f Fit the psf +r Review the fit for all the psf stars +s Subtract fitted psf from psf star nearest cursor +d Delete psf star nearest cursor from psf star list +w Write the psf to the psf image +z Rebuild the psf from scratch +q Quit task + + Colon Commands + +:p [n] Print photometry for star n +:a [n] Add star n to psf star list +:d [n] Delete star n from psf star list +:s [n] Subtract fitted psf from psf star n + + Colon Parameter Editing Commands + +# Data dependent parameters which affect the psf computation + +:scale [value] Show/set the image scale (units / pixel) +:fwhmpsf [value] Show/set the fwhm of psf (scale units) +:datamin [value] Show/set the minimum good data value (counts) +:datamax [value] Show/set the maximum good data value (counts) +:matchrad [value] Show/set matching radius (scale units) + +# Psf computation parameters + +:psfimage [name,name] Show/set the psf image and groupfile +:function [string] Show/set the analytic psf function +:varorder [integer] Show/set order of psf function variability +:nclean [integer] Show/set number of cleaning iterations +:saturated [y/n] Show/set the use saturated star flag +:psfrad [value] Show/set the psf radius (scale units) +:fitrad [value] Show/set the fitting radius (scale units) diff --git a/noao/digiphot/daophot/psf/mkpsflist.key b/noao/digiphot/daophot/psf/mkpsflist.key new file mode 100644 index 00000000..0ed43c8a --- /dev/null +++ b/noao/digiphot/daophot/psf/mkpsflist.key @@ -0,0 +1,15 @@ + Keystroke Commands + +? Print help +p Print photometry for star nearest the cursor +l List the current psf stars +n Select the next good candidate psf star from the list +a Add star nearest cursor to psf star list +d Delete psf star nearest cursor from psf star list +q Quit task + + Colon Commands + +:p [n] Print photometry for star n +:a [n] Add star n to psf star list +:d [n] Delete star n from psf star list diff --git a/noao/digiphot/daophot/psf/showpsf.key b/noao/digiphot/daophot/psf/showpsf.key new file mode 100644 index 00000000..7551477e --- /dev/null +++ b/noao/digiphot/daophot/psf/showpsf.key @@ -0,0 +1,24 @@ + Interactive Graphics Keystroke Commands + +? Print help +p Print the photometry for this star +t Print the plot parameters and data minimum and maximum +a Accept star and proceed +d Reject star and select another with image cursor +m Plot the default mesh plot for this star +n Increase vertical angle by 15 degrees (mesh plot only) +s Decrease vertical angle by 15 degrees (mesh plot only) +w Decrease horizontal angle by 15 degrees (mesh plot only) +e Increase horizontal angle by 15 degrees (mesh plot only) +c Plot the default contour plot for this star +r Plot the radial profile for this star + + + Colon Graphics Commands + +:m [val] [val] Set the mesh plot vertical and horizontal viewing angles +:v [val] Set the mesh plot vertical viewing angle +:h [val] Set the mesh plot horizontal viewing angle +:c [val] [val] Set the contour plot floor and ceiling levels +:l [value] Set the contour plot floor level +:u [value] Set the contour plot ceiling level diff --git a/noao/digiphot/daophot/psf/t_psf.x b/noao/digiphot/daophot/psf/t_psf.x new file mode 100644 index 00000000..b7444b3e --- /dev/null +++ b/noao/digiphot/daophot/psf/t_psf.x @@ -0,0 +1,509 @@ +include <fset.h> +include <imhdr.h> +include <gset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +# T_PSF -- Generate a point spread function from one or more stars in the +# image frame. + +procedure t_psf () + +pointer image # the input image +pointer photfile # input aperture photometry file +pointer pstarfile # input psf star file +pointer psfimage # output psf image +pointer groupfile # output psf group file +pointer opstfile # output psf star file +pointer graphics # pointer to graphics device name +pointer plotfile # pointer to plotfile name +int cache # cache the input image pixels +pointer display # pointer to display device name +bool matchbyid # match psf stars by id or position +bool interactive # the mode of task operation +bool showplots # display plots of the psf stars +pointer plottype # type of psf plot +bool mkstars # mark deleted and accepted psf stars + +pointer sp, im, apd, psfim, dao, mgd, gd, id +pointer outfname, curfile, str +int imlist, limlist, alist, lalist, clist, lclist, pimlist, lpimlist +int olist, lolist, oclist, loclist, up, verify, update, wcs +int root, min_lenuserarea, pltype, pfd, pst, psfgr, opst +int req_size, old_size, buf_size, memstat +bool ap_text, pst_text + +pointer immap(), tbtopn(), gopen() +int fnldir(), strlen(), strncmp(), btoi(), envfind(), ctoi(), clgwrd() +int strdic(), open(), access(), fstati(), dp_stati(), dp_pstati() +int imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb(), fntgfnb() +int sizeof(), dp_memstat() +bool streq(), clgetb(), itob(), dp_updatepsf() +errchk gopen + +begin + # Set the standard output to flush on newline. + if (fstati (STDOUT, F_REDIR) == NO) + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get some working memory. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (photfile, SZ_FNAME, TY_CHAR) + call salloc (pstarfile, SZ_FNAME, TY_CHAR) + call salloc (psfimage, SZ_FNAME, TY_CHAR) + call salloc (groupfile, SZ_FNAME, TY_CHAR) + call salloc (opstfile, SZ_FNAME, TY_CHAR) + call salloc (plottype, SZ_FNAME, TY_CHAR) + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (plotfile, SZ_FNAME, TY_CHAR) + call salloc (curfile, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the various task parameters. + call clgstr ("image", Memc[image], SZ_FNAME) + call clgstr ("photfile", Memc[photfile], SZ_FNAME) + call clgstr ("pstfile", Memc[pstarfile], SZ_FNAME) + call clgstr ("psfimage", Memc[psfimage], SZ_FNAME) + call clgstr ("opstfile", Memc[opstfile], SZ_FNAME) + call clgstr ("groupfile", Memc[groupfile], SZ_FNAME) + cache = btoi (clgetb ("cache")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + + # Get the lists. + imlist = imtopen (Memc[image]) + limlist = imtlen (imlist) + alist = fntopnb (Memc[photfile], NO) + lalist = fntlenb (alist) + clist = fntopnb (Memc[pstarfile], NO) + lclist = fntlenb (clist) + pimlist = imtopen (Memc[psfimage]) + lpimlist = imtlen (pimlist) + olist = fntopnb (Memc[groupfile], NO) + lolist = fntlenb (olist) + oclist = fntopnb (Memc[opstfile], NO) + loclist = fntlenb (oclist) + + # Test that the lengths of the photometry file, psf image, and + # output file lists are the same as the length of the input image + # list. + + # Compare the image and photometry file list lengths. + if ((limlist != lalist) && (strncmp (Memc[photfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + call sfree (sp) + call error (0, + "Incompatable image and photometry file list lengths") + } + + # Compare the image and psf star list lengths. + if ((lclist != 0) && (limlist != lclist) && (strncmp (Memc[pstarfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + call sfree (sp) + call error (0, + "Incompatable image and psf star file list lengths") + } + + # Compare the image and psf image list lengths. + if ((limlist != lpimlist) && (strncmp (Memc[psfimage], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + call sfree (sp) + call error (0, + "Incompatable image and psf file list lengths") + } + + # Compare the image and groupfile list lengths. + if ((limlist != lolist) && (strncmp (Memc[groupfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + call sfree (sp) + call error (0, + "Incompatable image and group file list lengths") + } + + # Compare the image and output pstfile list lengths. + if ((limlist != loclist) && (strncmp (Memc[opstfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + call sfree (sp) + call error (0, + "Incompatable image and output psf file list lengths") + } + + # Initialize DAOPHOT main structure, get pset parameters. + call dp_gppars (dao) + + # Verify the critical parameters and update if appropriate. + if (verify == YES) { + call dp_pconfirm (dao) + if (update == YES) + call dp_pppars (dao) + } + + # Get the wcs information. + wcs = clgwrd ("wcsin", Memc[str], SZ_FNAME, WCSINSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the input coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSIN, wcs) + wcs = clgwrd ("wcsout", Memc[str], SZ_FNAME, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSOUT, wcs) + + # Initialize the photometry structure. + call dp_apselsetup (dao) + + # Initialize the PSF structure. + call dp_fitsetup (dao) + + # Intialize the PSF fitting structure. + call dp_psfsetup (dao) + + # Matching algorithm for stars in the psf star list. + matchbyid = clgetb ("matchbyid") + + # Is the task interactive or not? + call clgstr ("icommands.p_filename", Memc[curfile], SZ_FNAME) + if (Memc[curfile] == EOS) + interactive = clgetb ("interactive") + else + interactive = false + + # Get the graphics, display and plot file devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + call clgstr ("plottype", Memc[plottype], SZ_FNAME) + call clgstr ("plotfile", Memc[plotfile], SZ_FNAME) + + # Open graphics and display devices if appropriate. + if (interactive) { + call dp_seti (dao, VERBOSE, YES) + showplots = clgetb ("showplots") + if (Memc[graphics] == EOS) + gd = NULL + else { + iferr { + gd = gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH) + } then { + call eprintf ( + "Warning: Error opening graphics device. \n") + gd = NULL + } + } + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) + id = gd + else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + if (id != NULL) + mkstars = clgetb ("mkstars") + else + mkstars = false + } else { + gd = NULL + id = NULL + call dp_seti (dao, VERBOSE, btoi (clgetb ("verbose"))) + showplots = false + mkstars = false + } + + # Open the plot file. + if (Memc[plotfile] == EOS) + pfd = NULL + else + pfd = open (Memc[plotfile], APPEND, BINARY_FILE) + if (pfd != NULL) + mgd = gopen (Memc[graphics], NEW_FILE, pfd) + else + mgd = NULL + + # Set the default plot type. + pltype = strdic (Memc[plottype], Memc[plottype], SZ_FNAME, PSF_PLOTS) + call dp_pseti (dao, PLOTTYPE, pltype) + + # Loop over the list of images. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open input image + im = immap (Memc[image], READ_ONLY, 0) + call dp_imkeys (dao, im) + call dp_sets (dao, INIMAGE, Memc[image]) + + # Set up the display coordinate system. + if ((id != NULL) && (id != gd)) + call dp_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = dp_memstat (cache, req_size, old_size) + if (memstat == YES) + call dp_pcache (im, INDEFI, buf_size) + + # Open the input photometry list and store the descriptor. + # PSF can read either an APPHOT PHOT file or an ST TABLE + # file. + + if (fntgfnb (alist, Memc[photfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[photfile], SZ_FNAME) + root = fnldir (Memc[photfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[photfile+root], DEF_LENDEFNAME) == + 0 || root == strlen (Memc[photfile])) + call dp_inname (Memc[image], Memc[outfname], "mag", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[photfile], Memc[outfname], SZ_FNAME) + ap_text = itob (access (Memc[outfname], 0, TEXT_FILE)) + if (ap_text) + apd = open (Memc[outfname], READ_ONLY, TEXT_FILE) + else + apd = tbtopn (Memc[outfname], READ_ONLY, 0) + call dp_wgetapert (dao, im, apd, dp_stati (dao, MAXNSTAR), ap_text) + call dp_sets (dao, INPHOTFILE, Memc[outfname]) + + # Open the input photometry list and store the descriptor. + # PSF can read either an APPHOT PHOT file or an ST TABLE + # file. + + if (lclist == 0) { + pst = NULL + Memc[outfname] = EOS + } else { + if (fntgfnb (clist, Memc[pstarfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[pstarfile], SZ_FNAME) + root = fnldir (Memc[pstarfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[pstarfile+root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[pstarfile])) + call dp_inname (Memc[image], Memc[outfname], "pst", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[pstarfile], Memc[outfname], SZ_FNAME) + pst_text = itob (access (Memc[outfname], 0, TEXT_FILE)) + if (pst_text) + pst = open (Memc[outfname], READ_ONLY, TEXT_FILE) + else + pst = tbtopn (Memc[outfname], READ_ONLY, 0) + } + call dp_sets (dao, COORDS, Memc[outfname]) + + # Open output image containing PSF, output file containing list + # of PSF stars actually used, and the file for PSF neighbors + # If the output is "default", dir$default or a directory + # specification then the extension "psf" is added to the PSF image + # and "psg" to the neighbors file. A suitable version number is + # added to the output name as well. Check that there is enough + # space in the image header user area to hold many PSF stars. + + if (imtgetim (pimlist, Memc[psfimage], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[psfimage], SZ_FNAME) + root = fnldir (Memc[psfimage], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[psfimage + root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[psfimage])) { + call dp_oimname (Memc[image], Memc[outfname], "psf", + Memc[outfname], SZ_FNAME) + } else + call strcpy (Memc[psfimage], Memc[outfname], SZ_FNAME) + if (envfind ("min_lenuserarea", Memc[str], SZ_FNAME) > 0) { + up = 1 + if (ctoi (Memc[str], up, min_lenuserarea) <= 0) + min_lenuserarea = MIN_LENUSERAREA + else + min_lenuserarea = max (MIN_LENUSERAREA, min_lenuserarea) + } else + min_lenuserarea = MIN_LENUSERAREA + call dp_pseti (dao, LENUSERAREA, min_lenuserarea) + psfim = immap (Memc[outfname], NEW_IMAGE, min_lenuserarea) + call dp_sets (dao, PSFIMAGE, Memc[outfname]) + + if (fntgfnb (olist, Memc[groupfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[groupfile], SZ_FNAME) + root = fnldir (Memc[groupfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[groupfile + root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[groupfile])) { + call dp_outname (Memc[image], Memc[outfname], "psg", + Memc[outfname], SZ_FNAME) + } else + call strcpy (Memc[groupfile], Memc[outfname], SZ_FNAME) + if (dp_stati (dao, TEXT) == YES) + psfgr = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + psfgr = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTPHOTFILE, Memc[outfname]) + + + if (fntgfnb (oclist, Memc[opstfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[opstfile], SZ_FNAME) + root = fnldir (Memc[opstfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[opstfile+root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[opstfile])) { + call dp_outname (Memc[image], Memc[outfname], "pst", + Memc[outfname], SZ_FNAME) + } else + call strcpy (Memc[opstfile], Memc[outfname], SZ_FNAME) + if (dp_stati (dao, TEXT) == YES) + opst = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + opst = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTREJFILE, Memc[outfname]) + + # Print banner. + if (DP_VERBOSE(dao) == YES) { + call printf ("\nComputing PSF for image: %s\n") + call pargstr (Memc[image]) + call dp_stats (dao, INPHOTFILE, Memc[str], SZ_FNAME) + call printf ("%d stars read from %s\n\n") + call pargi (dp_stati (dao, APNUM)) + call pargstr (Memc[str]) + } + + # Read in the PSF star list. + call dp_pseti (dao, PNUM, 0) + if (pst != NULL) { + call dp_rpstars (dao, im, pst, pst_text, gd, mgd, id, mkstars, + matchbyid, showplots) + if (DP_VERBOSE(dao) == YES) { + call dp_stats (dao, COORDS, Memc[str], SZ_FNAME) + call printf ("\n%d PSF stars read from %s\n\n") + call pargi (dp_pstati (dao, PNUM)) + call pargstr (Memc[str]) + } + } + + # Make the PSF. + if (Memc[curfile] != EOS) + call dp_mkpsf (dao, im, psfim, opst, psfgr, gd, mgd, id, false, + false, false) + else if (interactive) + call dp_mkpsf (dao, im, psfim, opst, psfgr, gd, mgd, id, + mkstars, true, true) + else if (! dp_updatepsf (dao, im, psfim, opst, psfgr, false, false, + false)) + call dp_rmpsf (dao, psfim, opst, psfgr) + + # Close the input image. + call imunmap (im) + + # Close the input photometry file. + if (apd != NULL) { + if (ap_text) + call close (apd) + else + call tbtclo (apd) + } + + # Close the input psf star file. + if (pst != NULL) { + if (pst_text) + call close (pst) + else + call tbtclo (pst) + } + + # Close PSF image. + if (psfim != NULL) + call imunmap (psfim) + + # Close the output PSF star file. + if (opst != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (opst) + else + call tbtclo (opst) + } + + # Close the group file. + if (psfgr != NULL) { + if (dp_stati (dao, TEXT) == YES) + call close (psfgr) + else + call tbtclo (psfgr) + } + + # Uncache memory. + call fixmem (old_size) + } + + # Close up the graphics and display streams. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + # Close the metacode plot files. + if (mgd != NULL) + call gclose (mgd) + if (pfd != NULL) + call close (pfd) + + # Close the image / file lists. + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (clist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (oclist) + + # Free the photometry structure. + call dp_apclose (dao) + + # Free the PSF fitting structure. + call dp_psfclose (dao) + + # Free the PSF structure. + call dp_fitclose (dao) + + # Close up the daophot structures. + call dp_free (dao) + + call sfree (sp) +end diff --git a/noao/digiphot/daophot/psf/t_pstselect.x b/noao/digiphot/daophot/psf/t_pstselect.x new file mode 100644 index 00000000..497ab374 --- /dev/null +++ b/noao/digiphot/daophot/psf/t_pstselect.x @@ -0,0 +1,329 @@ +include <fset.h> +include <imhdr.h> +include <gset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/psfdef.h" + +# T_PSTSELECT -- Select good candidate PSF stars from a DAOPHOT photometry +# file. + +procedure t_pstselect () + +pointer image # input image +pointer photfile # input rough photometry +pointer pstfile # output PSTFILE table +int maxnpsf # maximimum number of psf stars +pointer plotfile # pointer to the plot metacode file +bool interactive # interactive mode +pointer plottype # default plot type +int cache # cache the input image pixels +int verify # verify the critical parameters +int update # update the critical parameters +pointer graphics # the graphics device +pointer display # the display device +bool mkstars # mark deleted and accepted psf stars + +pointer sp, pfd, dao, outfname, curfile, str, im, gd, id, mgd +int imlist, limlist, alist, lalist, olist, lolist, root, apd, pmgd, pltype +int wcs, req_size, old_size, buf_size, memstat() +bool ap_text + +pointer immap(), gopen() +int tbtopn(), open(), fnldir(), strlen(), strncmp(), fstati(), btoi() +int access(), fntopnb(), fntlenb(), clgeti(), imtopen(), imtlen() +int fntgfnb(), imtgetim(), strdic(), dp_stati(), clgwrd(), sizeof() +int dp_memstat() +bool clgetb(), itob(), streq() + +errchk gopen() + +begin + # Set the standard output to flush on newline. + if (fstati (STDOUT, F_REDIR) == NO) + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get some working memory. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (photfile, SZ_FNAME, TY_CHAR) + call salloc (pstfile, SZ_FNAME, TY_CHAR) + call salloc (plotfile, SZ_FNAME, TY_CHAR) + call salloc (plottype, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (curfile, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the various task parameters. + call clgstr ("image", Memc[image], SZ_FNAME) + call clgstr ("photfile", Memc[photfile], SZ_FNAME) + call clgstr ("pstfile", Memc[pstfile], SZ_FNAME) + maxnpsf = clgeti ("maxnpsf") + cache = btoi (clgetb ("cache")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + + # Get the lists. + imlist = imtopen (Memc[image]) + limlist = imtlen (imlist) + alist = fntopnb (Memc[photfile], NO) + lalist = fntlenb (alist) + olist = fntopnb (Memc[pstfile], NO) + lolist = fntlenb (olist) + + # Test that the lengths of the photometry file and psf star file + # lists are the same as the input image list. + + if ((limlist != lalist) && (strncmp (Memc[photfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (olist) + call sfree (sp) + call error (0, + "Incompatable image and photometry file list lengths\n") + } + + if ((limlist != lolist) && (strncmp (Memc[pstfile], DEF_DEFNAME, + DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (olist) + call sfree (sp) + call error (0, + "Incompatable image and photometry file list lengths\n") + } + + # Initialize the DAOPHOT structure, and get the pset parameters. + call dp_gppars (dao) + + # Confirm the parameters. + if (verify == YES) { + call dp_ptconfirm (dao) + if (update == YES) + call dp_pppars (dao) + } + + # Get the wcs information. + wcs = clgwrd ("wcsin", Memc[str], SZ_LINE, WCSINSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the input coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSIN, wcs) + wcs = clgwrd ("wcsout", Memc[str], SZ_LINE, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSOUT, wcs) + + # Initialize the photometry structure. + call dp_apsetup (dao) + + # Initialize the PSF fitting structure. + call dp_psfsetup (dao) + + # Initialize the PSF structure. + call dp_fitsetup (dao) + + # Is the task interactive or not? + call clgstr ("icommands.p_filename", Memc[curfile], SZ_FNAME) + if (Memc[curfile] == EOS) + interactive = clgetb ("interactive") + else + interactive = false + + # Get the graphics display and plot file devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + call clgstr ("plottype", Memc[plottype], SZ_FNAME) + call clgstr ("plotfile", Memc[plotfile], SZ_FNAME) + + # Open graphics and display devices if appropriate. + if (interactive) { + call dp_seti (dao, VERBOSE, YES) + if (Memc[graphics] == EOS) + gd = NULL + else { + iferr { + gd = gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH) + } then { + call eprintf ("Warning: Error opening graphics device\n") + gd = NULL + } + } + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) { + id = gd + } else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device\n") + id = NULL + } + } + if (id != NULL) + mkstars = clgetb ("mkstars") + else + mkstars = false + } else { + gd = NULL + id = NULL + call dp_seti (dao, VERBOSE, btoi (clgetb ("verbose"))) + mkstars = false + } + + # Open the plot file. + if (Memc[plotfile] == EOS) + pmgd = NULL + else + pmgd = open (Memc[plotfile], APPEND, BINARY_FILE) + if (pmgd != NULL) + mgd = gopen (Memc[graphics], NEW_FILE, pmgd) + else + mgd = NULL + + # Set the default plot type. + pltype = strdic (Memc[plottype], Memc[plottype], SZ_FNAME, PSF_PLOTS) + call dp_pseti (dao, PLOTTYPE, pltype) + + # Loop over the list of input files + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the input image. + im = immap (Memc[image], READ_ONLY, 0) + call dp_imkeys (dao, im) + call dp_sets (dao, INIMAGE, Memc[image]) + + # Set up the display coordinate system. + if ((id != NULL) && (id != gd)) + call dp_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = dp_memstat (cache, req_size, old_size) + if (memstat == YES) + call dp_pcache (im, INDEFI, buf_size) + + # Open the input photometry table and read in the photometry. + if (fntgfnb (alist, Memc[photfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[photfile], SZ_FNAME) + root = fnldir (Memc[photfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[photfile+root], DEF_LENDEFNAME) == + 0 || root == strlen (Memc[photfile])) + call dp_inname (Memc[image], Memc[outfname], "mag", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[photfile], Memc[outfname], SZ_FNAME) + ap_text = itob (access (Memc[outfname], 0, TEXT_FILE)) + if (ap_text) + apd = open (Memc[outfname], READ_ONLY, TEXT_FILE) + else + apd = tbtopn (Memc[outfname], READ_ONLY, 0) + call dp_wgetapert (dao, im, apd, DP_MAXNSTAR(dao), ap_text) + call dp_sets (dao, INPHOTFILE, Memc[outfname]) + + # Open the output PSTSELECT file. If the output is "default", + # dir$default or a directory specification then the extension .pst + # is added to the image name and a suitable version number is + # appended to the output name. + + if (fntgfnb (olist, Memc[pstfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[pstfile], SZ_FNAME) + root = fnldir (Memc[pstfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[pstfile + root], DEF_LENDEFNAME) == + 0 || root == strlen (Memc[pstfile])) + call dp_outname (Memc[image], Memc[outfname], "pst", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[pstfile], Memc[outfname], SZ_FNAME) + if (ap_text) + pfd = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + pfd = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTPHOTFILE, Memc[outfname]) + + if (DP_VERBOSE(dao) == YES) { + call printf ("\nSelecting PSF stars for image %s\n") + call pargstr (Memc[image]) + call dp_stats (dao, INPHOTFILE, Memc[outfname], SZ_FNAME) + call printf ("\t%d stars read from file %s\n\n") + call pargi (dp_stati (dao, APNUM)) + call pargstr (Memc[outfname]) + } + + # Now select the PSF stars. + if (Memc[curfile] != EOS) + call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, NULL, + mgd, NULL, false, false, true) + else if (interactive) + call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, gd, mgd, + id, mkstars, true, false) + else + call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, NULL, + mgd, NULL, false, false, false) + + # Close the input image. + call imunmap (im) + + # Close the photometry file. + if (ap_text) + call close (apd) + else + call tbtclo (apd) + + # Close the output table. + if (ap_text) + call close (pfd) + else + call tbtclo (pfd) + + # Uncache memory. + call fixmem (old_size) + + } + + # Close up the graphics and display streams. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + # Close the metacode plot files. + if (mgd != NULL) + call gclose (mgd) + if (pmgd != NULL) + call close (pmgd) + + # Close the image/file lists. + call imtclose (imlist) + call fntclsb (alist) + call fntclsb (olist) + + # Close the PSF structure. + call dp_fitclose (dao) + + # Close the PSF fitting structure. + call dp_psfclose (dao) + + # Free the photometry structure. + call dp_apclose (dao) + + # Free the daophot structure. + call dp_free (dao) + + call sfree(sp) +end |