diff options
Diffstat (limited to 'noao/digiphot/daophot/peak')
| -rw-r--r-- | noao/digiphot/daophot/peak/dpmempk.x | 72 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/dppeakphot.x | 277 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/dppkconfirm.x | 25 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/dppkfit.x | 411 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/dppkwrite.x | 365 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/dprrphot.x | 98 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/mkpkg | 22 | ||||
| -rw-r--r-- | noao/digiphot/daophot/peak/t_peak.x | 299 |
8 files changed, 1569 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/peak/dpmempk.x b/noao/digiphot/daophot/peak/dpmempk.x new file mode 100644 index 00000000..97c7ff48 --- /dev/null +++ b/noao/digiphot/daophot/peak/dpmempk.x @@ -0,0 +1,72 @@ +include "../lib/daophotdef.h" +include "../lib/peakdef.h" + + +# DP_PKSETUP -- Initialize the PEAK fitting structure. + +procedure dp_pksetup (dao) + +pointer dao # pointer to the daophot structure + +pointer peak + +begin + call malloc (DP_PEAK(dao), LEN_PKSTRUCT, TY_STRUCT) + peak = DP_PEAK(dao) + + DP_PKNTERM(peak) = 0 + DP_PKCLAMP(peak) = NULL + DP_PKNORMAL(peak) = NULL + DP_PKRESID(peak) = NULL + DP_PKDERIV(peak) = NULL + DP_PKRESULT(peak) = NULL + DP_PKOLDRESULT(peak) = NULL +end + + +# DP_MEMPK -- Allocate memory for the PEAK fitting arrays. + +procedure dp_mempk (dao, nterm) + +pointer dao # pointer to the daophot structure +int nterm # the number of terms to be fit + +pointer peak + +begin + peak = DP_PEAK(dao) + + call malloc (DP_PKCLAMP(peak), nterm, TY_REAL) + call malloc (DP_PKNORMAL(peak), nterm * nterm, TY_REAL) + call malloc (DP_PKRESID(peak), nterm * nterm, TY_REAL) + call malloc (DP_PKDERIV(peak), nterm * nterm, TY_REAL) + call malloc (DP_PKRESULT(peak), nterm * nterm, TY_REAL) + call malloc (DP_PKOLDRESULT(peak), nterm * nterm, TY_REAL) + DP_PKNTERM(peak) = nterm +end + + +# DP_PKCLOSE -- Free the PEAK fitting structure. + +procedure dp_pkclose (dao) + +pointer dao # pointer to the daophot structure + +pointer peak + +begin + peak = DP_PEAK(dao) + if (DP_PKCLAMP(peak) != NULL) + call mfree (DP_PKCLAMP(peak), TY_REAL) + if (DP_PKNORMAL(peak) != NULL) + call mfree (DP_PKNORMAL(peak), TY_REAL) + if (DP_PKRESID(peak) != NULL) + call mfree (DP_PKRESID(peak), TY_REAL) + if (DP_PKDERIV(peak) != NULL) + call mfree (DP_PKDERIV(peak), TY_REAL) + if (DP_PKRESULT(peak) != NULL) + call mfree (DP_PKRESULT(peak), TY_REAL) + if (DP_PKOLDRESULT(peak) != NULL) + call mfree (DP_PKOLDRESULT(peak), TY_REAL) + call mfree (peak, TY_STRUCT) +end diff --git a/noao/digiphot/daophot/peak/dppeakphot.x b/noao/digiphot/daophot/peak/dppeakphot.x new file mode 100644 index 00000000..182b10d8 --- /dev/null +++ b/noao/digiphot/daophot/peak/dppeakphot.x @@ -0,0 +1,277 @@ +include <imhdr.h> +include <tbset.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/peakdef.h" + + +# DP_PEAKPHOT -- Fit the PSF to a single star. + +procedure dp_peakphot (dao, im, tp, tpout, tprej, ap_text) + +pointer dao # pointer to the daophot structure +pointer im # the input image descriptor +int tp # the input photometry file descriptor +int tpout # the ouput photometry file descriptor +int tprej # the rejections file descriptor +bool ap_text # which style of photometry + +real rel_bright, xold, yold, x, y, dx, dy, mag, sky, errmag +real chi, sharp, radius, itx, ity, otx, oty +pointer psffit, key, sp, subim, colpoint, indices, fields, perror +int id, in_nrow, instar, lowx, lowy, nxpix, nypix, niter, out_nrow +int rout_nrow, nterm, ier, plen + +int tbpsta(), dp_rrphot(), dp_pkfit(), dp_gpkerr() +pointer dp_gsubrast() + +begin + # Get the daophot pointers. + psffit = DP_PSFFIT (dao) + + # Store the original fitting radius. + radius = DP_FITRAD(dao) + + # Check that the fitting radius is less than the psf radius. + DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao)) + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) + + # Allocate working space. + call smark (sp) + call salloc (indices, NAPPAR, TY_INT) + call salloc (fields, SZ_LINE, TY_CHAR) + call salloc (colpoint, PK_NOUTCOL, TY_INT) + call salloc (perror, SZ_FNAME, TY_CHAR) + + # Initialze the output table. + if (DP_TEXT(dao) == YES) { + call dp_xnewpeak (dao, tpout) + if (tprej != NULL) + call dp_xnewpeak (dao, tprej) + } else { + call dp_tnewpeak (dao, tpout, Memi[colpoint]) + if (tprej != NULL) + call dp_tnewpeak (dao, tprej, Memi[colpoint]) + } + + # Intialize the input table. + if (ap_text) { + call pt_kyinit (key) + Memi[indices] = DP_PAPID + Memi[indices+1] = DP_PAPXCEN + Memi[indices+2] = DP_PAPYCEN + Memi[indices+3] = DP_PAPMAG1 + Memi[indices+4] = DP_PAPSKY + call dp_gappsf (Memi[indices], Memc[fields], NAPRESULT) + in_nrow = 0 + } else { + call dp_tpkinit (tp, Memi[indices]) + in_nrow = tbpsta (tp, TBL_NROWS) + } + + # Initialize the photometry file reading code. + instar = 0 + + # Initialize the fitting code. + if (DP_RECENTER(dao) == YES) + nterm = 3 + else + nterm = 1 + if (DP_FITSKY(dao) == YES) + nterm = nterm + 1 + call dp_mempk (dao, nterm) + + out_nrow = 0 + rout_nrow = 0 + repeat { + + # Read in the photometry for a single star. + if (dp_rrphot (tp, key, Memc[fields], Memi[indices], id, itx, + ity, sky, mag, instar, in_nrow) == EOF) + break + + # Convert to and from logical coordinates. + call dp_win (dao, im, itx, ity, x, y, 1) + call dp_wout (dao, im, x, y, otx, oty, 1) + + if (DP_VERBOSE(dao) == YES) { + call printf ( + "Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky: %8.2f\n") + call pargi (id) + call pargr (otx) + call pargr (oty) + call pargr (mag) + call pargr (sky) + } + + # Check that the center is defined. + if (IS_INDEFR(x) || IS_INDEFR(y)) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "\tWarning: X and/or Y for star %d are undefined\n") + call pargi (id) + } + next + } + + # Read in the subraster. + subim = dp_gsubrast (im, x, y, DP_FITRAD(dao), lowx, lowy, + nxpix, nypix) + if (subim == NULL) { + if (DP_VERBOSE(dao) == YES) { + call printf ( + "\tWarning: Cannot read in image data for star %d\n") + call pargi (id) + } + next + } + + # Save the old x and y values for use with the variable psf + # option. + xold = x + yold = y + call dp_wpsf (dao, im, xold, yold, xold, yold, 1) + + # Compute the relative centers and the relative brightness and + # fit the star. + if (IS_INDEFR(sky)) { + + ier = PKERR_INDEFSKY + + } else { + + x = x - lowx + 1.0 + y = y - lowy + 1.0 + dx = (xold - 1.0) / DP_PSFX(psffit) - 1.0 + dy = (yold - 1.0) / DP_PSFY(psffit) - 1.0 + if (IS_INDEFR(mag)) + mag = DP_PSFMAG (psffit) + DELTA_MAG + rel_bright = DAO_RELBRIGHT (psffit, mag) + ier = dp_pkfit (dao, Memr[subim], nxpix, nypix, DP_FITRAD(dao), + x, y, dx, dy, rel_bright, sky, errmag, chi, sharp, niter) + x = x + lowx - 1.0 + y = y + lowy - 1.0 + } + + call dp_wout (dao, im, x, y, otx, oty, 1) + + if (ier != PKERR_OK) { + + # Set fitting parameters to INDEF. + mag = INDEFR + niter = 0 + errmag = INDEFR + chi = INDEFR + sharp = INDEFR + + if (DP_VERBOSE(dao) == YES) { + switch (ier) { + case PKERR_INDEFSKY: + call printf ( + "\tWarning: The sky value for star %d is undefined\n") + call pargi (id) + case PKERR_NOPIX: + call printf ( + "\tWarning: Too few pixels to fit star %d\n") + call pargi (id) + case PKERR_SINGULAR: + call printf ( + "\tWarning: Singular matrix computed for star %d\n") + call pargi (id) + case PKERR_FAINT: + call printf ("\tWarning: Star %d too faint\n") + call pargi (id) + default: + call printf ( + "\tWarning: Unknown error detected for star %d\n") + call pargi (id) + } + } + + } else { + + # Compute the results. + mag = DP_PSFMAG (psffit) - 2.5 * log10 (rel_bright) + errmag = 1.085736 * errmag / rel_bright + if (errmag >= 2.0) + errmag = INDEFR + + if (DP_VERBOSE (dao) == YES) { + call printf ( + "\tFIT: Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky =%8.2f\n") + call pargi (id) + call pargr (otx) + call pargr (oty) + call pargr (mag) + call pargr (sky) + } + } + + # Get the error code. + plen = dp_gpkerr (ier, Memc[perror], SZ_FNAME) + + # Now write the results to the output photometry or rejections + # file. + if (DP_TEXT(dao) == YES) { + if ((tprej != NULL) && (ier != PKERR_OK)) + call dp_xpkwrite (tprej, id, otx, oty, mag, errmag, sky, + niter, chi, sharp, ier, Memc[perror], plen) + else + call dp_xpkwrite (tpout, id, otx, oty, mag, errmag, sky, + niter, chi, sharp, ier, Memc[perror], plen) + } else { + if ((tprej != NULL) && (ier != PKERR_OK)) { + rout_nrow = rout_nrow + 1 + call dp_tpkwrite (tprej, Memi[colpoint], id, otx, oty, mag, + errmag, sky, niter, chi, sharp, ier, + Memc[perror], plen, rout_nrow) + } else { + out_nrow = out_nrow + 1 + call dp_tpkwrite (tpout, Memi[colpoint], id, otx, oty, mag, + errmag, sky, niter, chi, sharp, ier, Memc[perror], + plen, out_nrow) + } + } + } + + # Free the text file descriptor. + if (ap_text) + call pt_kyfree (key) + + # Restore the original fitting radius. + DP_FITRAD(dao) = radius + DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao) + + call sfree (sp) +end + + +# DP_GPKERR -- Set the PEAK task error code string. + +int procedure dp_gpkerr (ier, perror, maxch) + +int ier # the integer error code +char perror # the output error code string +int maxch # the maximum size of the error code string + +int plen +int gstrcpy() + +begin + switch (ier) { + case PKERR_OK: + plen = gstrcpy ("No_error", perror, maxch) + case PKERR_INDEFSKY: + plen = gstrcpy ("Bad_sky", perror, maxch) + case PKERR_NOPIX: + plen = gstrcpy ("Npix_too_few", perror, maxch) + case PKERR_SINGULAR: + plen = gstrcpy ("Singular", perror, maxch) + case PKERR_FAINT: + plen = gstrcpy ("Too_faint", perror, maxch) + default: + plen = gstrcpy ("No_error", perror, maxch) + } + + return (plen) +end diff --git a/noao/digiphot/daophot/peak/dppkconfirm.x b/noao/digiphot/daophot/peak/dppkconfirm.x new file mode 100644 index 00000000..c7dc6333 --- /dev/null +++ b/noao/digiphot/daophot/peak/dppkconfirm.x @@ -0,0 +1,25 @@ +# DP_PKCONFIRM -- Procedure to confirm the critical peak parameters. + +procedure dp_pkconfirm (dao) + +pointer dao # pointer to the daophot structure + +begin + call printf ("\n") + + # Confirm recentering and sky fitting. + call dp_vrecenter (dao) + call dp_vfitsky (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/peak/dppkfit.x b/noao/digiphot/daophot/peak/dppkfit.x new file mode 100644 index 00000000..7818100c --- /dev/null +++ b/noao/digiphot/daophot/peak/dppkfit.x @@ -0,0 +1,411 @@ +include <mach.h> +include "../lib/daophotdef.h" +include "../lib/peakdef.h" + +# DP_PKFIT -- Do the actual fit. + +int procedure dp_pkfit (dao, subim, nx, ny, radius, x, y, dx, dy, rel_bright, + sky, errmag, chi, sharp, iter) + +pointer dao # pointer to the DAOPHOT Structure +real subim[nx,ny] # pointer to the image subraster +int nx, ny # size of the image subraster +real radius # the fitting radius +real x, y # initial estimate of the stellar position +real dx, dy # distance of star from the psf position +real rel_bright # initial estimate of stellar brightness +real sky # the sky value associated with this star +real errmag # error estimate for this star +real chi # estimated goodness of fit parameter +real sharp # broadness of the profile compared to the PSF +int iter # number of iterations needed for a fit. + +bool clip, redo +int i, flag, npix +pointer psffit, peak +real ronoise, numer, denom, chiold, sum_weight, noise, wcrit +int dp_fitbuild() + +begin + # Get the pointer to the PSF structure. + psffit = DP_PSFFIT (dao) + peak = DP_PEAK(dao) + + # Initialize the parameters which control the fit. + + chiold = 1.0 + sharp = INDEFR + clip = false + ronoise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2 + + call amovkr (1.0, Memr[DP_PKCLAMP(peak)], DP_PKNTERM(peak)) + call amovkr (0.0, Memr[DP_PKOLDRESULT(peak)], DP_PKNTERM(peak)) + + # Iterate until a solution is found. + + for (iter = 1; iter <= DP_MAXITER(dao); iter = iter + 1) { + + # Initialize the matrices and vectors required by the fit. + + chi = 0.0 + numer = 0.0 + denom = 0.0 + sum_weight = 0.0 + call aclrr (Memr[DP_PKRESID(peak)], DP_PKNTERM(peak)) + call aclrr (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak) * + DP_PKNTERM(peak)) + + # Set up the critical error limit. + if (iter >= WCRIT_NMAX) + wcrit = WCRIT_MAX + else if (iter >= WCRIT_NMED) + wcrit = WCRIT_MED + else if (iter >= WCRIT_NMIN) + wcrit = WCRIT_MIN + else + wcrit = MAX_REAL + + + # Build up the vector of residuals and the normal matrix. + + npix = dp_fitbuild (dao, subim, nx, ny, radius, x, y, dx, dy, + rel_bright, sky, chiold, chi, clip, iter, + Memr[DP_PKCLAMP(peak)], Memr[DP_PKNORMAL(peak)], + Memr[DP_PKRESID(peak)], Memr[DP_PKDERIV(peak)], + DP_PKNTERM[(peak)], numer, denom, sum_weight) + + # Return if the iteration was unsuccessful. + + if (npix < MIN_NPIX) + return (PKERR_NOPIX) + + # Invert the matrix. Return if inversion was unsuccessful or + # if any of the diagonal elements are less or equal to 0.0. + + call invers (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak), + DP_PKNTERM(peak), flag) + if (flag != 0) + return (PKERR_SINGULAR) + do i = 1, DP_PKNTERM(peak) { + if (Memr[DP_PKNORMAL(peak)+(i-1)*DP_PKNTERM(peak)+i-1] <= 0.0) + return (PKERR_SINGULAR) + } + # Solve the matrix. + + call mvmul (Memr[DP_PKNORMAL(peak)], DP_PKNTERM(peak), + DP_PKNTERM(peak), Memr[DP_PKRESID(peak)], + Memr[DP_PKRESULT(peak)]) + + # In the beginning the brightness of the star will not be + # permitted to change by more than two magnitudes per + # iteration (that is to say if the estimate is getting + # brighter, it may not get brighter by more than 525% + # per iteration, and if it is getting fainter, it may not + # get fainter by more than 84% per iteration). The x and y + # coordinates of the centroid will be allowed to change by + # no more than one-half pixel per iteration. Any time + # that a parameter correction changes sign, the maximum + # permissible change in that parameter will be reduced + # by a factor of 2. + + # Perform the sign check. + + do i = 1, DP_PKNTERM(peak) { + if ((Memr[DP_PKOLDRESULT(peak)+i-1] * + Memr[DP_PKRESULT(peak)+i-1]) < 0.0) + Memr[DP_PKCLAMP(peak)+i-1] = 0.5 * + Memr[DP_PKCLAMP(peak)+i-1] + else + Memr[DP_PKCLAMP(peak)+i-1] = min (1.0, 1.1 * + Memr[DP_PKCLAMP(peak)+i-1]) + Memr[DP_PKOLDRESULT(peak)+i-1] = Memr[DP_PKRESULT(peak)+i-1] + } + + # Compute the new x, y, sky, and magnitude. + rel_bright = rel_bright + Memr[DP_PKRESULT(peak)] / + (1.0 + max (Memr[DP_PKRESULT(peak)] / + (MAX_DELTA_BRIGHTER * rel_bright), -Memr[DP_PKRESULT(peak)] / + (MAX_DELTA_FAINTER * rel_bright)) / Memr[DP_PKCLAMP(peak)]) + + # Return if the star becomes too faint) + if (rel_bright < MIN_REL_BRIGHT) + return (PKERR_FAINT) + + if (DP_RECENTER(dao) == YES) { + x = x + Memr[DP_PKRESULT(peak)+1] / (1.0 + + abs (Memr[DP_PKRESULT(peak)+1]) / (MAX_DELTA_PIX * + Memr[DP_PKCLAMP(peak)+1])) + y = y + Memr[DP_PKRESULT(peak)+2] / (1.0 + + abs (Memr[DP_PKRESULT(peak)+2]) / (MAX_DELTA_PIX * + Memr[DP_PKCLAMP(peak)+2])) + } + + if (DP_FITSKY(dao) == YES) { + noise = sqrt ((abs (sky / DP_PHOTADU(dao)) + ronoise)) + sky = sky + max (-3.0 * noise, min (Memr[DP_PKRESULT(peak)+ + DP_PKNTERM(peak)-1], 3.0 * noise)) + } + + # Force at least one iteration before checking for convergence. + if (iter <= 1) + next + + # Start on the assumption the fit has converged. + + redo = false + + # Check for convergence. If the most recent computed correction + # to the brightness is larger than 0.1% of the brightness or + # 0.05 * sigma of the brightness whichever is larger, convergence + # has not been achieved. + + errmag = chiold * sqrt (Memr[DP_PKNORMAL(peak)]) + if (clip) { + if (abs (Memr[DP_PKRESULT(peak)]) > max ((MAX_NEW_ERRMAG * + errmag), (MAX_NEW_RELBRIGHT1 * rel_bright))) { + redo = true + } else { + if (DP_RECENTER(dao) == YES) { + if (max (abs (Memr[DP_PKRESULT(peak)+1]), + abs (Memr[DP_PKRESULT(peak)+2])) > MAX_PIXERR1) + redo = true + } + if (DP_FITSKY(dao) == YES) { + if (abs (Memr[DP_PKRESULT(peak)+DP_PKNTERM(peak)-1]) > + 1.0e-4 * sky) + redo = true + } + } + } else { + if (abs (Memr[DP_PKRESULT(peak)]) > max (errmag, + (MAX_NEW_RELBRIGHT2 * rel_bright))) { + redo = true + } else { + if (DP_RECENTER(dao) == YES) { + if (max (abs (Memr[DP_PKRESULT(peak)+1]), + abs (Memr[DP_PKRESULT(peak)+2])) > MAX_PIXERR2) + redo = true + } + if (DP_FITSKY(dao) == YES) { + if (abs (Memr[DP_PKRESULT(peak)+DP_PKNTERM(peak)-1]) > + 1.0e-4 * sky) + redo = true + } + } + } + if (redo) + next + + if ((iter < DP_MAXITER(dao)) && (! clip)) { + if (DP_CLIPEXP(dao) > 0) + clip = true + call aclrr (Memr[DP_PKOLDRESULT(peak)], DP_PKNTERM(peak)) + call amaxkr (Memr[DP_PKCLAMP(peak)], MAX_CLAMPFACTOR, + Memr[DP_PKCLAMP(peak)], DP_PKNTERM(peak)) + } else { + sharp = 1.4427 * Memr[DP_PSFPARS(psffit)] * + Memr[DP_PSFPARS(psffit)+1] * numer / (DP_PSFHEIGHT(psffit) * + rel_bright * denom) + if (sharp <= MIN_SHARP || sharp >= MAX_SHARP) + sharp = INDEFR + break + } + } + + if (iter > DP_MAXITER(dao)) + iter = iter - 1 + if ((errmag / rel_bright) >= wcrit) + return (PKERR_FAINT) + else + return (PKERR_OK) +end + + +# DP_FITBUILD -- Build the normal and vector of residuals for the fit. + +int procedure dp_fitbuild (dao, subim, nx, ny, radius, x, y, xfrom_psf, + yfrom_psf, rel_bright, sky, chiold, chi, clip, iter, clamp, normal, + resid, deriv, nterm, numer, denom, sum_weight) + +pointer dao # pointer to the DAOPHOT Structure +real subim[nx,ny] # subimage containing star +int nx, ny # size of the image subraster +real radius # the fitting radius +real x, y # initial estimate of the position +real xfrom_psf, yfrom_psf # distance from the psf star +real rel_bright # initial estimate of stellar brightness +real sky # the sky value associated with this star +real chiold # previous estimate of gof +real chi # estimated goodness of fit parameter +bool clip # clip the weights ? +int iter # current iteration number +real clamp[ARB] # clamp array +real normal[nterm,ARB] # normal matrix +real resid[ARB] # residual matrix +real deriv[ARB] # derivative matrix +int nterm # the number of terms to be fit +real numer, denom # used in sharpness calculation +real sum_weight # sum of the weights + +int i, j, ix, iy, lowx, lowy, highx, highy, npix +pointer psffit +real fitradsq, pererr, peakerr, datamin, datamax, read_noise +real dx, dy, dxsq, dysq, radsq, dvdx, dvdy, d_pixval +real pred_pixval, sigma, sigmasq, relerr, weight +real rhosq, pixval, dfdsig + +real dp_usepsf() + +begin + # Get the pointer to the PSF structure. + psffit = DP_PSFFIT (dao) + + # Set up some constants. + fitradsq = radius * radius + read_noise = (DP_READNOISE(dao) / DP_PHOTADU(dao)) ** 2 + pererr = 0.01 * DP_FLATERR(dao) + peakerr = 0.01 * DP_PROFERR(dao) / + (Memr[DP_PSFPARS(psffit)] * Memr[DP_PSFPARS(psffit)+1]) + 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) + + # Define the size of the subraster to be used in the fit. + lowx = max (1, int (x - radius)) + lowy = max (1, int (y - radius)) + highx = min (nx, int (x + radius) + 1) + highy = min (ny, int (y + radius) + 1) + + npix = 0 + do iy = lowy, highy { + + dy = real (iy) - y + dysq = dy * dy + + do ix = lowx, highx { + + # Is this pixel in the good data range. + pixval = subim[ix,iy] + if (pixval < datamin || pixval > datamax) + next + + # Is this pixel inside the fitting radius? + dx = real(ix) - x + dxsq = dx * dx + radsq = (dxsq + dysq) / fitradsq + if ((1.0 - radsq) <= PEAK_EPSILONR) + next + + # We have a good pixel within the fitting radius. + deriv[1] = dp_usepsf (DP_PSFUNCTION(psffit), dx, dy, + DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)], + Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), + DP_NVLTABLE(psffit), DP_NFEXTABLE(psffit), xfrom_psf, + yfrom_psf, dvdx, dvdy) + if (((rel_bright * deriv[1] + sky) > datamax) && (iter >= 3)) + next + if (DP_RECENTER(dao) == YES) { + deriv[2] = rel_bright * dvdx + deriv[3] = rel_bright * dvdy + } + if (DP_FITSKY(dao) == YES) + deriv[nterm] = 1.0 + + # Get the residual from the PSF fit and the pixel + # intensity as predicted by the fit + d_pixval = (pixval - sky) - rel_bright * deriv[1] + pred_pixval = max (0.0, pixval - d_pixval) + + # Calculate the anticipated error in the intensity of + # in this pixel including READOUT noise, photon statistics + # and the error of interpolating within the PSF. + + sigmasq = pred_pixval / DP_PHOTADU (dao) + read_noise + + (pererr * pred_pixval) ** 2 + (peakerr * + (pred_pixval - sky)) ** 2 + if (sigmasq <= 0.0) + next + sigma = sqrt (sigmasq) + relerr = abs (d_pixval / sigma) + + # Compute the radial wweighting function. + weight = 5.0 / (5.0 + radsq / (1.0 - radsq)) + + # Now add this pixel into the quantities which go to make + # up the sharpness index. + rhosq = dxsq / Memr[DP_PSFPARS(psffit)] ** 2 + dysq / + Memr[DP_PSFPARS(psffit)+1] ** 2 + + # Include in the sharpness calculation only those pixels + # which are within NCORE_SIGMASQ core-sigmasq of the + # centroid. This saves time and floating underflows + # bu excluding pixels which contribute less than one + # part in a million to the fit. + + if (rhosq <= NCORE_SIGMASQ) { + rhosq = 0.6931472 * rhosq + dfdsig = exp (-rhosq) * (rhosq - 1.0) + #pred_pixval = max (0.0, pixval - sky) + sky + numer = numer + dfdsig * d_pixval / sigmasq + denom = denom + (dfdsig ** 2) / sigmasq + } + + + # Derive the weight of this pixel. First of all the weight + # depends on the distance of the pixel from the centroid of + # the star --- it is determined from a function which is very + # nearly unity for radii much smaller than the fitting radius, + # and which goes to zero for radii very near the fitting + # radius. Then reject any pixels with 10 sigma residuals + # after the first iteration. + + chi = chi + weight * relerr + sum_weight = sum_weight + weight + + # Now the weight is scaled to the inverse square of the + # expected mean error. + + weight = weight / sigmasq + + # Reduce the weight of a bad pixel. A pixel having a residual + # of 2.5 sigma gets reduced to half weight; a pixel having + # a residual of of 5.0 sigma gets weight 1/257. + + if (clip) + weight = weight / (1.0 + (relerr / (DP_CLIPRANGE(dao) * + chiold)) ** DP_CLIPEXP(dao)) + + # Now add this pixel into the vector of residuals + # and the normal matrix. + do i = 1, nterm { + resid[i] = resid[i] + d_pixval * deriv[i] * weight + do j = 1, nterm { + normal[i,j] = normal[i,j] + deriv[i] * deriv[j] * + weight + } + } + + npix = npix + 1 + } + } + + # Compute the goodness of fit index CHI. CHI is pulled toward its + # expected value of unity before being stored in CHIOLD to keep the + # statistics of a small number of pixels from dominating the error + # analysis. + + if (sum_weight > nterm) { + chi = CHI_NORM * chi / sqrt (sum_weight * (sum_weight - nterm)) + chiold = ((sum_weight - MIN_SUMWT) * chi + MIN_SUMWT) / sum_weight + } else { + chi = 1.0 + chiold = 1.0 + } + + return (npix) +end diff --git a/noao/digiphot/daophot/peak/dppkwrite.x b/noao/digiphot/daophot/peak/dppkwrite.x new file mode 100644 index 00000000..97f770bf --- /dev/null +++ b/noao/digiphot/daophot/peak/dppkwrite.x @@ -0,0 +1,365 @@ +include <tbset.h> +include <time.h> +include "../lib/daophotdef.h" +include "../lib/apseldef.h" +include "../lib/peakdef.h" + +# DP_TNEWPEAK -- Create a new output PEAK ST table. + +procedure dp_tnewpeak (dao, tp, colpoint) + +pointer dao # pointer to the daophot strucuture +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, PK_NOUTCOL * (SZ_COLNAME + 1), TY_CHAR) + call salloc (colunits, PK_NOUTCOL * (SZ_COLUNITS + 1), TY_CHAR) + call salloc (colformat, PK_NOUTCOL * (SZ_COLFMT + 1), TY_CHAR) + call salloc (col_dtype, PK_NOUTCOL, TY_INT) + call salloc (col_len, PK_NOUTCOL, 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 (MAGERR, Memc[colnames+4*SZ_COLNAME+4], SZ_COLNAME) + call strcpy (SKY, Memc[colnames+5*SZ_COLNAME+5], SZ_COLNAME) + call strcpy (NITER, Memc[colnames+6*SZ_COLNAME+6], SZ_COLNAME) + call strcpy (SHARP, Memc[colnames+7*SZ_COLNAME+7], SZ_COLNAME) + call strcpy (CHI, Memc[colnames+8*SZ_COLNAME+8], SZ_COLNAME) + call strcpy (PIER, Memc[colnames+9*SZ_COLNAME+9], SZ_COLNAME) + call strcpy (PERROR, Memc[colnames+10*SZ_COLNAME+10], 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 ("%14.3f", Memc[colformat+4*SZ_COLFMT+4], SZ_COLFMT) + call strcpy ("%15.7g", Memc[colformat+5*SZ_COLFMT+5], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+6*SZ_COLFMT+6], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+7*SZ_COLFMT+7], SZ_COLFMT) + call strcpy ("%12.3f", Memc[colformat+8*SZ_COLFMT+8], SZ_COLFMT) + call strcpy ("%6d", Memc[colformat+9*SZ_COLFMT+9], SZ_COLFMT) + call strcpy ("%13s", Memc[colformat+10*SZ_COLFMT+10], 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 ("MAGNITUDES", Memc[colunits+4*SZ_COLUNITS+4], SZ_COLUNITS) + call strcpy ("ADU", Memc[colunits+5*SZ_COLUNITS+5], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+6*SZ_COLUNITS+6], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+7*SZ_COLUNITS+7], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+8*SZ_COLUNITS+8], SZ_COLUNITS) + call strcpy ("NUMBER", Memc[colunits+9*SZ_COLUNITS+9], SZ_COLUNITS) + call strcpy ("PERRORS", Memc[colunits+10*SZ_COLUNITS+10], 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 + Memi[col_dtype+5] = TY_REAL + Memi[col_dtype+6] = TY_INT + Memi[col_dtype+7] = TY_REAL + Memi[col_dtype+8] = TY_REAL + Memi[col_dtype+9] = TY_INT + Memi[col_dtype+10] = -13 + + # Define the column lengths. + do i = 1, PK_NOUTCOL + 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], PK_NOUTCOL) + call tbtcre (tp) + + # Write out some header parameters. + call dp_tpeakpars (dao, tp) + + call sfree (sp) +end + + +define PK_NAME1STR "#N%4tID%10tXCENTER%20tYCENTER%30tMAG%42tMERR%56tMSKY%71t\ +NITER%80t\\\n" +define PK_UNIT1STR "#U%4t##%10tpixels%20tpixels%30tmagnitudes%42t\ +magnitudes%56tcounts%71t##%80t\\\n" +define PK_FORMAT1STR "#F%4t%%-9d%10t%%-10.3f%20t%%-10.3f%30t%%-12.3f%42t\ +%%-14.3f%56t%%-15.7g%71t%%-6d%80t \n" + +define PK_NAME2STR "#N%12tSHARPNESS%24tCHI%36tPIER%42tPERROR%80t\\\n" +define PK_UNIT2STR "#U%12t##%24t##%36t##%42tperrors%80t\\\n" +define PK_FORMAT2STR "#F%12t%%-23.3f%24t%%-12.3f%36t%%-6d%42t%%-13s%80t \n" + + +# DP_XNEWPEAK -- Create a new PEAK output text file. + +procedure dp_xnewpeak (dao, tp) + +pointer dao # pointer to the daophot structure +int tp # the output file descriptor + +begin + # Write out some header parameters. + call dp_xpeakpars (dao, tp) + + # Print out the banner file. + call fprintf (tp, "#\n") + call fprintf (tp, PK_NAME1STR) + call fprintf (tp, PK_UNIT1STR) + call fprintf (tp, PK_FORMAT1STR) + call fprintf (tp, "#\n") + call fprintf (tp, PK_NAME2STR) + call fprintf (tp, PK_UNIT2STR) + call fprintf (tp, PK_FORMAT2STR) + call fprintf (tp, "#\n") +end + + +# DP_TPEAKPARS -- Add parameters to the header of the PEAK ST table. + +procedure dp_tpeakpars (dao, tp) + +pointer dao # pointer to the daophot structure +pointer tp # pointer to the output table + +pointer sp, psffit, outstr, date, time +bool itob() +int envfind() + +begin + # Define some daophot pointers. + psffit = DP_PSFFIT(dao) + + # 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) + + # Write the id. + 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]) + call tbhadt (tp, "PACKAGE", "daophot") + call tbhadt (tp, "TASK", "peak") + + # Write out the files names. + 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]) + call dp_imroot (DP_PSFIMAGE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PSFIMAGE", Memc[outstr]) + call dp_froot (DP_OUTPHOTFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "PEAKFILE", Memc[outstr]) + if (DP_OUTREJFILE(dao) == EOS) + call tbhadt (tp, "REJFILE", "\"\"") + else { + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE) + call tbhadt (tp, "REJFILE", Memc[outstr]) + } + + # Define the data characteristics. + call tbhadr (tp, "SCALE", DP_SCALE(dao)) + call tbhadr (tp, "DATAMIN", DP_MINGDATA(dao)) + call tbhadr (tp, "DATAMAX", DP_MAXGDATA(dao)) + call tbhadr (tp, "GAIN", DP_PHOTADU(dao)) + call tbhadr (tp, "READNOISE", DP_READNOISE(dao)) + + # Define the observing parameters. + call tbhadt (tp, "OTIME", DP_OTIME(dao)) + call tbhadr (tp, "XAIRMASS", DP_XAIRMASS(dao)) + call tbhadt (tp, "IFILTER", DP_IFILTER(dao)) + + # Define the fitting parameters. + call tbhadb (tp, "RECENTER", itob (DP_RECENTER(dao))) + call tbhadb (tp, "FITSKY", itob (DP_FITSKY(dao))) + call tbhadr (tp, "PSFRAD", DP_SPSFRAD(dao)) + call tbhadr (tp, "FITRAD", DP_SFITRAD(dao)) + call tbhadr (tp, "PSFMAG", DP_PSFMAG(psffit)) + call tbhadi (tp, "MAXITER", DP_MAXITER(dao)) + call tbhadr (tp, "FLATERROR", DP_FLATERR(dao)) + call tbhadr (tp, "PROFERROR", DP_PROFERR(dao)) + call tbhadi (tp, "CLIPEXP", DP_CLIPEXP(dao)) + call tbhadr (tp, "CLIPRANGE", DP_CLIPRANGE(dao)) + #call tbhadi (tp, "MAXNSTAR", DP_MAXNSTAR(dao)) + + call sfree(sp) +end + + +# DP_XPEAKPARS -- Add parameters to the header of the output PEAK text file. + +procedure dp_xpeakpars (dao, tp) + +pointer dao # pointer to the daophot structure +int tp # the output file descriptor + +pointer sp, psffit, outstr, date, time +bool itob() +int envfind() + +begin + # Define some daophot pointers. + psffit = DP_PSFFIT(dao) + + # 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) + + # Write the id. + 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", "") + call dp_sparam (tp, "PACKAGE", "daophot", "name", "") + call dp_sparam (tp, "TASK", "peak", "name", "") + + # Write out the files names. + 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", "") + 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, "PEAKFILE", Memc[outstr], "filename", "") + + if (DP_OUTREJFILE(dao) == EOS) + call dp_sparam (tp, "REJFILE", "\"\"", "filename", "") + else { + call dp_froot (DP_OUTREJFILE(dao), Memc[outstr], SZ_LINE) + call dp_sparam (tp, "REJFILE", Memc[outstr], "filename", "") + } + + # Define the data characteristics. + call dp_rparam (tp, "SCALE", DP_SCALE(dao), "units/pix", "") + call dp_rparam (tp, "DATAMIN", DP_MINGDATA(dao), "counts", "") + call dp_rparam (tp, "DATAMAX", DP_MAXGDATA(dao), "counts", "") + call dp_rparam (tp, "GAIN", DP_PHOTADU(dao), "e-/adu", "") + call dp_rparam (tp, "READNOISE", DP_READNOISE(dao), "electrons", "") + + # Define the observing parameters. + call dp_sparam (tp, "OTIME", DP_OTIME(dao), "timeunit", "") + call dp_rparam (tp, "XAIRMASS", DP_XAIRMASS(dao), "number", "") + call dp_sparam (tp, "IFILTER", DP_IFILTER(dao), "filter", "") + + # Define the fitting parameters. + call dp_bparam (tp, "RECENTER", itob (DP_RECENTER(dao)), "switch", "") + call dp_bparam (tp, "FITSKY", itob (DP_FITSKY(dao)), "switch", "") + call dp_rparam (tp, "PSFRAD", DP_SPSFRAD(dao), "scaleunit", "") + call dp_rparam (tp, "FITRAD", DP_SFITRAD(dao), "scaleunit", "") + call dp_rparam (tp, "PSFMAG", DP_PSFMAG(psffit), "magnitudes", "") + call dp_iparam (tp, "MAXITER", DP_MAXITER(dao), "number", "") + call dp_rparam (tp, "FLATERROR", DP_FLATERR(dao), "percentage", "") + call dp_rparam (tp, "PROFERROR", DP_PROFERR(dao), "percentage", "") + call dp_iparam (tp, "CLIPEXP", DP_CLIPEXP(dao), "number", "") + call dp_rparam (tp, "CLIPRANGE", DP_CLIPRANGE(dao), "sigma", "") + #call dp_iparam (tp, "MAXNSTAR", DP_MAXNSTAR(dao), "number", "") + + call sfree(sp) +end + + +# DP_TPKWRITE -- Write the output photometry record to a PEAK ST table. + +procedure dp_tpkwrite (tpout, colpoint, id, x, y, mag, errmag, sky, niter, chi, + sharp, pier, perror, plen, star) + +pointer tpout # pointer to the output table +int colpoint[ARB] # array of column pointers +int id # id of the star +real x, y # position of the star +real mag # magnitude of the star +real errmag # error magnitude of the star +real sky # value of sky +int niter # number of iterations +real chi # chi squared value +real sharp # sharpness characteristic +int pier # the error code +char perror[ARB] # the error string +int plen # the length of the error code string +int star # row number + +begin + call tbrpti (tpout, colpoint[1], id, 1, star) + call tbrptr (tpout, colpoint[2], x, 1, star) + call tbrptr (tpout, colpoint[3], y, 1, star) + call tbrptr (tpout, colpoint[4], mag, 1, star) + call tbrptr (tpout, colpoint[5], errmag, 1, star) + call tbrptr (tpout, colpoint[6], sky, 1, star) + call tbrpti (tpout, colpoint[7], niter, 1, star) + call tbrptr (tpout, colpoint[8], sharp, 1, star) + call tbrptr (tpout, colpoint[9], chi, 1, star) + call tbrpti (tpout, colpoint[10], pier, 1, star) + call tbrptt (tpout, colpoint[11], perror, plen, 1, star) +end + + +define PK_DATA1STR "%-9d%10t%-10.3f%20t%-10.3f%30t%-12.3f%42t%-14.3f%56t\ +%-15.7g%71t%-6d%80t\\\n" +define PK_DATA2STR "%12t%-12.3f%24t%-12.3f%36t%-6d%42t%-13.13s%80t \n" + +# DP_XPKWRITE -- Write the output photometry record to a PEAK text file. + +procedure dp_xpkwrite (tpout, id, x, y, mag, errmag, sky, niter, chi, sharp, + pier, perror, plen) + +int tpout # the output file descriptor +int id # id of the star +real x, y # position of the star +real mag # magnitude of the star +real errmag # error magnitude of the star +real sky # value of sky +int niter # number of iterations +real chi # chi squared value +real sharp # sharpness characteristic +int pier # the error code +char perror[ARB] # the error string +int plen # the length of the error code string + +begin + call fprintf (tpout, PK_DATA1STR) + call pargi (id) + call pargr (x) + call pargr (y) + call pargr (mag) + call pargr (errmag) + call pargr (sky) + call pargi (niter) + call fprintf (tpout, PK_DATA2STR) + call pargr (sharp) + call pargr (chi) + call pargi (pier) + call pargstr (perror) +end diff --git a/noao/digiphot/daophot/peak/dprrphot.x b/noao/digiphot/daophot/peak/dprrphot.x new file mode 100644 index 00000000..36aed88f --- /dev/null +++ b/noao/digiphot/daophot/peak/dprrphot.x @@ -0,0 +1,98 @@ +include "../lib/apseldef.h" + +# DP_TPKINIT -- Procedure to initialize for reading the "standard" fields from +# a photometry table. The "standard" fields being ID, X, Y, MAG, ERR, and SKY + +procedure dp_tpkinit (tp, colpoint) + +pointer tp # the table descriptor +pointer colpoint[ARB] # the column descriptor + +begin + # Get the results one by one + # First the ID + call tbcfnd (tp, ID, colpoint[1], 1) + if (colpoint[1] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (ID) + } + + # Then the position + call tbcfnd (tp, XCENTER, colpoint[2], 1) + if (colpoint[2] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (XCENTER) + } + + call tbcfnd (tp, YCENTER, colpoint[3], 1) + if (colpoint[3] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (YCENTER) + } + + # Now the Magnitude + call tbcfnd (tp, MAG, colpoint[4], 1) + if (colpoint[4] == NULL) # No column + call tbcfnd (tp, APMAG, colpoint[4], 1) + if (colpoint[4] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (APMAG) + } + + # The sky + call tbcfnd (tp, SKY, colpoint[5], 1) + if (colpoint[5] == NULL) + call tbcfnd (tp, SKY, colpoint[5], 1) + if (colpoint[5] == NULL) { + call eprintf ("Column %s not found\n") + call pargstr (SKY) + } +end + + +# DP_RRPHOT -- Fetch the photometry for a single star from either a +# table or a text photometry file. + +int procedure dp_rrphot (tp, key, fields, indices, id, x, y, sky, mag, + instar, nrow) + +int tp # the input file descriptor +pointer key # pointer to text apphot structure +char fields[ARB] # character fields +int indices[ARB] # columns pointers +int id # star id +real x # x center value +real y # y center value +real sky # sky value +real mag # magnitude +int instar # current record +int nrow # maximum number of rows for ST table + +bool nullflag +int nrec +int dp_apsel() + +begin + # If nrow is 0 the file file is a text file otherwise it is a table. + + if (nrow == 0) { + + nrec = dp_apsel (key, tp, fields, indices, id, x, y, sky, mag) + if (nrec != EOF) + instar = instar + 1 + + } else if ((instar + 1) <= nrow) { + + instar = instar + 1 + call tbrgti (tp, indices[1], id, nullflag, 1, instar) + call tbrgtr (tp, indices[2], x, nullflag, 1, instar) + call tbrgtr (tp, indices[3], y, nullflag, 1, instar) + call tbrgtr (tp, indices[4], mag, nullflag, 1, instar) + call tbrgtr (tp, indices[5], sky, nullflag, 1, instar) + nrec = instar + + } else + nrec = EOF + + return (nrec) +end diff --git a/noao/digiphot/daophot/peak/mkpkg b/noao/digiphot/daophot/peak/mkpkg new file mode 100644 index 00000000..f3570728 --- /dev/null +++ b/noao/digiphot/daophot/peak/mkpkg @@ -0,0 +1,22 @@ +# PEAK task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + dpmempk.x ../lib/daophotdef.h ../lib/peakdef.h + dppeakphot.x <imhdr.h> <tbset.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/peakdef.h + dppkconfirm.x + dppkfit.x <mach.h> ../lib/daophotdef.h \ + ../lib/peakdef.h + dppkwrite.x <tbset.h> <time.h> \ + ../lib/daophotdef.h ../lib/apseldef.h \ + ../lib/peakdef.h + dprrphot.x ../lib/apseldef.h + t_peak.x <fset.h> <imhdr.h> \ + ../lib/daophotdef.h + ; diff --git a/noao/digiphot/daophot/peak/t_peak.x b/noao/digiphot/daophot/peak/t_peak.x new file mode 100644 index 00000000..2022fe16 --- /dev/null +++ b/noao/digiphot/daophot/peak/t_peak.x @@ -0,0 +1,299 @@ +include <fset.h> +include <imhdr.h> +include "../lib/daophotdef.h" + +# T_PEAK -- Procedure to fit the PSF to single stars. + +procedure t_peak () + +pointer image # name of the image +pointer photfile # input photometry file +pointer psfimage # the PSF image +pointer peakfile # output PEAK photometry file +pointer rejfile # output PEAK rejections file + +pointer sp, im, psfim, outfname, dao, str +int apd, root, verbose, verify, cache, update, pkfd, rejfd +int imlist, limlist, alist, lalist, pimlist, lpimlist, olist, lolist +int rlist, lrlist, wcs, req_size, old_size, buf_size, memstat +bool ap_text + +pointer immap(), tbtopn() +int open(), fnldir(), strlen(), strncmp(), fstati(), btoi() +int access(), imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb() +int fntgfnb(), clgwrd(), sizeof(), dp_memstat() +bool clgetb(), itob() + +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 (psfimage, SZ_FNAME, TY_CHAR) + call salloc (peakfile, SZ_FNAME, TY_CHAR) + call salloc (rejfile, SZ_FNAME, TY_CHAR) + call salloc (outfname, 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 ("psfimage", Memc[psfimage], SZ_FNAME) + call clgstr ("peakfile", Memc[peakfile], SZ_FNAME) + call clgstr ("rejfile", Memc[rejfile], SZ_FNAME) + verbose = btoi (clgetb ("verbose")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + cache = btoi (clgetb ("cache")) + + # Get the lists. + imlist = imtopen (Memc[image]) + limlist = imtlen (imlist) + alist = fntopnb (Memc[photfile], NO) + lalist = fntlenb (alist) + pimlist = imtopen (Memc[psfimage]) + lpimlist = imtlen (pimlist) + olist = fntopnb (Memc[peakfile], NO) + lolist = fntlenb (olist) + rlist = fntopnb (Memc[rejfile], NO) + lrlist = fntlenb (rlist) + + # 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. + + if ((limlist != lalist) && (strncmp (Memc[photfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and input photometry file list lengths") + } + + if ((limlist != lpimlist) && (strncmp (Memc[psfimage], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and psf image list lengths") + } + + if ((limlist != lolist) && (strncmp (Memc[peakfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and output photometry file list lengths") + } + + if ((lrlist != 0) && (limlist != lolist) && (strncmp (Memc[rejfile], + DEF_DEFNAME, DEF_LENDEFNAME) != 0)) { + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + call sfree (sp) + call error (0, + "Incompatable image and rejections file list lengths") + } + + # Initialize the DAOPHOT structure, and get the pset parameters. + call dp_gppars (dao) + call dp_seti (dao, VERBOSE, verbose) + + # Verify the standard algorithm parameters. + if (verify == YES) { + call dp_pkconfirm (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) + wcs = clgwrd ("wcspsf", Memc[str], SZ_FNAME, WCSPSFSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the psf coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call dp_seti (dao, WCSPSF, wcs) + + # Initialize the PSF structure. + call dp_fitsetup (dao) + + # Initialize the PEAK structure. + call dp_pksetup (dao) + + # Loop over the list of images. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the image. + im = immap (Memc[image], READ_ONLY, 0) + call dp_imkeys (dao, im) + call dp_sets (dao, INIMAGE, Memc[image]) + + # 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 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_sets (dao, INPHOTFILE, Memc[outfname]) + + # Read in the PSF function. + 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_iimname (Memc[image], Memc[outfname], "psf", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[psfimage], Memc[outfname], SZ_FNAME) + psfim = immap (Memc[outfname], READ_ONLY, 0) + call dp_readpsf (dao, psfim) + call dp_sets (dao, PSFIMAGE, Memc[outfname]) + + # Open the output photometry file. If the output is "default", + # dir$default or a directory specification then the extension .pk + # is added to the image name and a suitable version number is + # appended to the output name. + + if (fntgfnb (olist, Memc[peakfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[peakfile], SZ_FNAME) + root = fnldir (Memc[peakfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[peakfile + root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[peakfile])) + call dp_outname (Memc[image], Memc[outfname], "pk", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[peakfile], Memc[outfname], SZ_FNAME) + if (DP_TEXT(dao) == YES) + pkfd = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + pkfd = tbtopn (Memc[outfname], NEW_FILE, 0) + call dp_sets (dao, OUTPHOTFILE, Memc[outfname]) + + # Open the output rejections file if any are defined. If the + # output is "default", dir$default or a directory specification + # then the extension .prj is added to the image name and a + # suitable version number is appended to the output file name. + + if (lrlist <= 0) { + rejfd = NULL + Memc[outfname] = EOS + } else { + if (fntgfnb (rlist, Memc[rejfile], SZ_FNAME) == EOF) + call strcpy (DEF_DEFNAME, Memc[rejfile], SZ_FNAME) + root = fnldir (Memc[rejfile], Memc[outfname], SZ_FNAME) + if (strncmp (DEF_DEFNAME, Memc[rejfile+root], + DEF_LENDEFNAME) == 0 || root == strlen (Memc[rejfile])) + call dp_outname (Memc[image], Memc[outfname], "prj", + Memc[outfname], SZ_FNAME) + else + call strcpy (Memc[rejfile], Memc[outfname], SZ_FNAME) + if (DP_TEXT(dao) == YES) + rejfd = open (Memc[outfname], NEW_FILE, TEXT_FILE) + else + rejfd = tbtopn (Memc[outfname], NEW_FILE, 0) + } + call dp_sets (dao, OUTREJFILE, Memc[outfname]) + + # Now go and do the PSF fitting. + call dp_peakphot (dao, im, apd, pkfd, rejfd, ap_text) + + # Close the input image. + call imunmap (im) + + # Close the input photometry file. + if (ap_text) + call close (apd) + else + call tbtclo (apd) + + # Close the PSF image. + call imunmap (psfim) + + # Close the output photometry file. + if (DP_TEXT(dao) == YES) + call close (pkfd) + else + call tbtclo (pkfd) + + # Close the output rejections file. + if (rejfd != NULL) { + if (DP_TEXT(dao) == YES) + call close (rejfd) + else + call tbtclo (rejfd) + } + + # Uncache memory. + call fixmem (old_size) + + } + + # Close the image/file lists. + call imtclose (imlist) + call fntclsb (alist) + call imtclose (pimlist) + call fntclsb (olist) + call fntclsb (rlist) + + # Free the PEAK structure. + call dp_pkclose (dao) + + # Free the PSF structure. + call dp_fitclose (dao) + + # Free the daophot structure. + call dp_free (dao) + + call sfree(sp) +end |