diff options
Diffstat (limited to 'noao/digiphot/apphot/radprof')
21 files changed, 3004 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/radprof/apbradprof.x b/noao/digiphot/apphot/radprof/apbradprof.x new file mode 100644 index 00000000..08e5547b --- /dev/null +++ b/noao/digiphot/apphot/radprof/apbradprof.x @@ -0,0 +1,107 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +# AP_BRADPROF -- Procedure to fit the radial profiles of a list of objects +# in a list of images. + +procedure ap_bradprof (ap, im, cl, gid, gd, mgd, out, id, ld, interactive) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +pointer gid # display device stream +pointer gd # pointer to graphics stream +pointer mgd # pointer to metacode file +int out # output file descriptor +int id, ld # sequence and list numbers +int interactive # interactive pr batch mode + +int stdin, ild, cier, sier, pier, rier +pointer sp, str +real wx, wy +int fscan(), nscan(), apfitsky(), apfitcenter(), ap_frprof(), apstati() +int strncmp() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME + 1, TY_CHAR) + call fstats (cl, F_FILENAME, Memc[str], SZ_FNAME) + + # Initialize. + ild = ld + + # Print query if coords is equal to STDIN. + if (strncmp ("STDIN", Memc[str], 5) == 0) { + stdin = YES + call printf ( + "Type x and y coordinates of object (^D or ^Z to end): ") + call flush (STDOUT) + } else + stdin = NO + + # Loop over the coordinate file. + while (fscan (cl) != EOF) { + + # Get the coordinates. + call gargr (wx) + call gargr (wy) + if (nscan () != 2) { + if (stdin == YES) { + call printf ( + "Type x and y coordinates of object (^D or ^Z to end): ") + call flush (STDOUT) + } + next + } + + # Transform the input coordinates. + switch (apstati(ap,WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_itol (ap, wx, wy, wx, wy, 1) + case WCS_TV: + call ap_vtol (im, wx, wy, wx, wy, 1) + default: + ; + } + + # Set the current cursor position. + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Fit the center, sky and radial profile. + cier = apfitcenter (ap, im, wx, wy) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), NULL, gd) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), pier) + + # Output the results. + if (interactive == YES) { + call ap_qprprof (ap, cier, sier, pier, rier) + if (id != NULL) + call aprmark (ap, gid, apstati (ap, MKCENTER), + apstati (ap, MKSKY), apstati (ap, MKAPERT)) + } + if (id == 1) + call ap_param (ap, out, "radprof") + call ap_prprof (ap, out, id, ild, cier, sier, pier, rier) + call ap_rpplot (ap, id, mgd, YES) + + # Setup for the next object. + id = id + 1 + ild = ild + 1 + call apsetr (ap, WX, wx) + call apsetr (ap, WY, wy) + if (stdin == YES) { + call printf ( + "Type x and y coordinates of object (^D or ^Z to end): ") + call flush (STDOUT) + } + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/apfrprof.x b/noao/digiphot/apphot/radprof/apfrprof.x new file mode 100644 index 00000000..093e90bd --- /dev/null +++ b/noao/digiphot/apphot/radprof/apfrprof.x @@ -0,0 +1,417 @@ +include <mach.h> +include <gset.h> +include <math.h> + +include <math/curfit.h> +include <math/iminterp.h> +include "../lib/apphotdef.h" +include "../lib/noisedef.h" +include "../lib/fitskydef.h" +include "../lib/photdef.h" +include "../lib/radprofdef.h" +include "../lib/apphot.h" +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_FRPROF -- Compute the radial profile of an object. + +int procedure ap_frprof (ap, im, wx, wy, pier) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +real wx, wy # object coordinates +int pier # photometry error + +int i, ier, fier, nxpts, nypts, nrpts, order +pointer sky, rprof, cv, asi +real datamin, datamax, step, rmin, rmax, inorm, tinorm +int ap_rpbuf(), ap_rmag(), ap_rpmeasure(), ap_rpiter() +real asigrl(), cveval(), ap_rphalf() + +errchk cvinit(), cvfit(), cvvector(), cveval(), cvfree() +errchk asinit(), asifit(), asigrl(), asifree() +errchk ap_rpmeasure(), ap_rpiter() + +begin + # Set up some apphot pointers. + sky = AP_PSKY(ap) + rprof = AP_RPROF(ap) + + # Initialize. + AP_RPXCUR(rprof) = wx + AP_RPYCUR(rprof) = wy + call ap_rpindef (ap) + if (IS_INDEFR(wx) || IS_INDEFR(wy)) { + AP_ORPXCUR(rprof) = INDEFR + AP_ORPYCUR(rprof) = INDEFR + } else { + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_ORPXCUR(rprof), + AP_ORPYCUR(rprof), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_ORPXCUR(rprof), + AP_ORPYCUR(rprof), 1) + default: + AP_ORPXCUR(rprof) = wx + AP_ORPYCUR(rprof) = wy + } + } + + # Get the pixels and check for error conditions. + if (IS_INDEFR(wx) || IS_INDEFR(wy)) { + pier = AP_APERT_NOAPERT + return (AP_RP_NOPROFILE) + } else if (IS_INDEFR(AP_SKY_MODE(sky))) { + pier = AP_APERT_NOSKYMODE + return (AP_RP_NOSKYMODE) + } else if (ap_rpbuf (ap, im, wx, wy) == AP_RP_NOPROFILE) { + pier = AP_APERT_NOAPERT + return (AP_RP_NOPROFILE) + } + + # Do the photometry. + pier = ap_rmag (ap) + + # Initialize some common variables. + nxpts = AP_RPNX(rprof) + nypts = AP_RPNY(rprof) + nrpts = AP_RPNPTS(rprof) + + # Initialize the radial profile curve fitting. + step = AP_SCALE(ap) * AP_RPSTEP(rprof) + order = max (1, min (AP_RPORDER(rprof), nxpts * nypts - 3)) + rmin = 0.0 + rmax = (nrpts - 1) * step + if (IS_INDEFR(AP_DATAMIN(ap))) + datamin = -MAX_REAL + else + datamin = AP_DATAMIN(ap) - AP_SKY_MODE(sky) + if (IS_INDEFR(AP_DATAMAX(ap))) + datamax = MAX_REAL + else + datamax = AP_DATAMAX(ap) - AP_SKY_MODE(sky) + + # Fit the curve. + call cvinit (cv, SPLINE3, order, rmin, rmax) + call asubkr (Memr[AP_RPIX(rprof)], AP_SKY_MODE(sky), + Memr[AP_RPIX(rprof)], nxpts * nypts) + fier = ap_rpmeasure (cv, Memr[AP_RPIX(rprof)], nxpts, nypts, + AP_RPXC(rprof), AP_RPYC(rprof), datamin, datamax, rmax, + AP_RPNDATA(rprof), AP_RPNBAD(rprof)) + + # Perform the rejection cycle. + if (fier != NO_DEG_FREEDOM && AP_RPNREJECT(rprof) > 0 && + AP_RPKSIGMA(rprof) > 0.0) + AP_RPNDATAREJ(rprof) = ap_rpiter (cv, Memr[AP_RPIX(rprof)], nxpts, + nypts, AP_RPXC(rprof), AP_RPYC(rprof), rmax, datamin, datamax, + AP_RPNREJECT(rprof), AP_RPKSIGMA(rprof), fier) + else + AP_RPNDATAREJ(rprof) = 0 + + AP_RPNDATA(rprof) = AP_RPNDATA(rprof) - AP_RPNDATAREJ(rprof) + AP_RPNDATAREJ(rprof) = AP_RPNDATAREJ(rprof) + AP_RPNBAD(rprof) + + # Evaluate the fit. + if (fier != NO_DEG_FREEDOM) { + + # Evaluate the profile. + do i = 1, nrpts + Memr[AP_RPDIST(rprof)+i-1] = (i - 1) * step + call cvvector (cv, Memr[AP_RPDIST(rprof)], + Memr[AP_INTENSITY(rprof)], nrpts) + + # Evaluate the integral. + call asiinit (asi, II_SPLINE3) + call amulr (Memr[AP_RPDIST(rprof)], Memr[AP_INTENSITY(rprof)], + Memr[AP_TINTENSITY(rprof)], nrpts) + call asifit (asi, Memr[AP_TINTENSITY(rprof)], nrpts) + Memr[AP_TINTENSITY(rprof)] = 0.0 + do i = 2, nrpts + Memr[AP_TINTENSITY(rprof)+i-1] = Memr[AP_TINTENSITY(rprof)+ + i-2] + asigrl (asi, real (i - 1), real (i)) + call amulkr (Memr[AP_TINTENSITY(rprof)], real (TWOPI) * step, + Memr[AP_TINTENSITY(rprof)], nrpts) + call asifree (asi) + + # Normalize the radial profile. + inorm = cveval (cv, 0.0) + if (inorm != 0.0) + call adivkr (Memr[AP_INTENSITY(rprof)], inorm, + Memr[AP_INTENSITY(rprof)], nrpts) + call apsetr (ap, INORM, inorm) + + # Normalize the total intensity. + tinorm = Memr[AP_TINTENSITY(rprof)+AP_RPNPTS(rprof)-1] + if (tinorm != 0.0) + call adivkr (Memr[AP_TINTENSITY(rprof)], tinorm, + Memr[AP_TINTENSITY(rprof)], nrpts) + call apsetr (ap, TNORM, tinorm) + + # Compute the FWHMPSF. + call apsetr (ap, RPFWHM, 2.0 * ap_rphalf (Memr[AP_RPDIST(rprof)], + Memr[AP_INTENSITY(rprof)], nrpts)) + } + + # Set the error code and return. + call cvfree (cv) + if (fier == NO_DEG_FREEDOM) + ier = AP_RP_NPTS_TOO_SMALL + else if (fier == SINGULAR) + ier = AP_RP_SINGULAR + + # Free space. + return (ier) +end + + +# AP_RMAG -- Compute the magnitudes for the radial profile + +int procedure ap_rmag (ap) + +pointer ap # pointer to the apphot structure + +int pier, nap +pointer sp, nse, sky, phot, rprof, aperts +real datamin, datamax, zmag + +begin + # Initialize some apphot pointers. + nse = AP_NOISE(ap) + sky = AP_PSKY(ap) + phot = AP_PPHOT(ap) + rprof = AP_RPROF(ap) + + # Allocate working space. + call smark (sp) + call salloc (aperts, AP_NAPERTS(phot), TY_REAL) + + # Check for out of bounds apertures. + call ap_maxap (ap, pier) + + # Define the good data minimum and maximum for photometry. + if (IS_INDEFR(AP_DATAMIN(ap))) + datamin = -MAX_REAL + else + datamin = AP_DATAMIN(ap) + if (IS_INDEFR(AP_DATAMAX(ap))) + datamax = MAX_REAL + else + datamax = AP_DATAMAX(ap) + + # Compute the sums. + call ap_arrayr (ap, APERTS, Memr[aperts]) + call amulkr (Memr[aperts], AP_SCALE(ap), Memr[aperts], AP_NAPERTS(phot)) + if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap))) { + call ap_rmmeasure (Memr[AP_RPIX(rprof)], AP_RPNX(rprof), + AP_RPNY(rprof), AP_RPXC(rprof), AP_RPYC(rprof), Memr[aperts], + Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot)) + AP_NMINAP(phot) = AP_NMAXAP(phot) + 1 + } else + call ap_brmmeasure (Memr[AP_RPIX(rprof)], AP_RPNX(rprof), + AP_RPNY(rprof), AP_RPXC(rprof), AP_RPYC(rprof), datamin, + datamax, Memr[aperts], Memd[AP_SUMS(phot)], + Memd[AP_AREA(phot)], AP_NMAXAP(phot), AP_NMINAP(phot)) + + # Check for bad pixels. + if ((pier == AP_OK) && (AP_NMINAP(phot) <= AP_NMAXAP(phot))) + pier = AP_APERT_BADDATA + nap = min (AP_NMINAP(phot) - 1, AP_NMAXAP(phot)) + + # Compute the magnitudes. + zmag = AP_ZMAG(phot) + 2.5 * log10 (AP_ITIME(ap)) + if (AP_POSITIVE(ap) == YES) + call apcopmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], + Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], nap, + AP_SKY_MODE(sky), AP_SKY_SIG(sky), AP_NSKY(sky), zmag, + AP_NOISEFUNCTION(nse), AP_EPADU(nse)) + else + call apconmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], + Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], nap, + AP_SKY_MODE(sky), AP_SKY_SIG(sky), AP_NSKY(sky), zmag, + AP_NOISEFUNCTION(nse), AP_EPADU(nse), AP_READNOISE(nse)) + + call sfree (sp) + + return (pier) +end + + +# AP_RPMEASURE -- Procedure to measure the flux and effective area in a set of +# apertures and accumulate the fit to the radial profile. + +int procedure ap_rpmeasure (cv, pixels, nx, ny, wx, wy, datamin, datamax, + rmax, ndata, nbad) + +pointer cv # pointer to curfit structure +real pixels[nx,ARB] # subraster pixel values +int nx, ny # dimensions of the subraster +real wx, wy # center of subraster +real datamin # minimum good data value +real datamax # maximum good data value +real rmax # the maximum radius +int ndata # the number of ok data points +int nbad # the number of bad data points + +int i, j, ier +real weight, dy2, r2, rcv + +begin + # Initialize. + ndata = 0 + nbad = 0 + call cvzero (cv) + + # Loop over the pixels. + do j = 1, ny { + dy2 = (j - wy) ** 2 + do i = 1, nx { + r2 = (i - wx) ** 2 + dy2 + rcv = sqrt (r2) + if (rcv > rmax) + next + if (pixels[i,j] < datamin || pixels[i,j] > datamax) { + nbad = nbad + 1 + } else { + call cvaccum (cv, rcv, pixels[i,j], weight, WTS_UNIFORM) + ndata = ndata + 1 + } + } + } + + call cvsolve (cv, ier) + return (ier) +end + + +# AP_RPITER -- Procedure to reject pixels from the fit. + +int procedure ap_rpiter (cv, pixels, nx, ny, wx, wy, rmax, datamin, datamax, + niter, ksigma, fier) + +pointer cv # pointer to the curfit structure +real pixels[nx,ARB] # pixel values +int nx, ny # dimensions of image subraster +real wx, wy # x and y coordinates of the center +real rmax # maximum radius value +real datamin # minimum good data value +real datamax # maximum good data value +int niter # maximum number of rejection cycles +real ksigma # ksigma rejection limit +int fier # fitting error code + +int i, j, k, npts, ntreject, nreject +pointer sp, rtemp, w, ptr +real chisqr, diff, locut, hicut +real cveval() +errchk cveval, cvrject, cvsolve + +begin + # Allocate the necessary space. + call smark (sp) + call salloc (rtemp, nx, TY_REAL) + call salloc (w, nx * ny, TY_REAL) + call amovkr (1.0, Memr[w], nx * ny) + + # Set the weights of out of range and bad data points to 0.0. + ptr = w + do j = 1, ny { + call ap_ijtor (Memr[rtemp], nx, j, wx, wy) + do k = 1, nx { + if (Memr[rtemp+k-1] > rmax || pixels[k,j] < datamin || + pixels[k,j] > datamax) + Memr[ptr+k-1] = 0.0 + } + ptr = ptr + nx + } + + ntreject = 0 + do i = 1, niter { + + # Compute the chisqr. + chisqr = 0.0 + npts = 0 + ptr = w + do j = 1, ny { + call ap_ijtor (Memr[rtemp], nx, j, wx, wy) + do k = 1, nx { + if (Memr[ptr+k-1] <= 0.0) + next + chisqr = chisqr + (pixels[k,j] - cveval (cv, + Memr[rtemp+k-1])) ** 2 + npts = npts + 1 + } + ptr = ptr + nx + } + + # Compute the new limits. + if (npts > 1) + chisqr = sqrt (chisqr / (npts - 1)) + else + chisqr = 0.0 + locut = - ksigma * chisqr + hicut = ksigma * chisqr + + # Reject pixels from the fit. + nreject = 0 + ptr = w + do j = 1, ny { + call ap_ijtor (Memr[rtemp], nx, j, wx, wy) + do k = 1, nx { + if (Memr[ptr+k-1] <= 0.0) + next + diff = pixels[k,j] - cveval (cv, Memr[rtemp+k-1]) + if (diff >= locut && diff <= hicut) + next + call cvrject (cv, Memr[rtemp+k-1], pixels[k,j], 1.0) + nreject = nreject + 1 + Memr[ptr+k-1] = 0.0 + } + ptr = ptr + nx + } + + if (nreject == 0) + break + ntreject = ntreject + nreject + + # Recompute the fit. + call cvsolve (cv, fier) + if (fier == NO_DEG_FREEDOM) + break + } + + + call sfree (sp) + + return (ntreject) +end + + +# AP_RPHALF -- Compute the FWHM of the PSF. + +real procedure ap_rphalf (radius, intensity, npts) + +real radius[ARB] # radius in pixels +real intensity[ARB] # profile intensity +int npts # number of points + +int i +real halfp + +begin + # Seach for the appropriate interval. + do i = 1, npts + if (intensity[i] < 0.5) + break + + # Compute the full width half maximum. + if (i == 1) + halfp = radius[1] + else if (i == npts && intensity[npts] >= 0.5) + halfp = radius[npts] + else + halfp = (radius[i] * (0.5 - intensity[i-1]) + radius[i-1] * + (intensity[i] - 0.5)) / (intensity[i] - intensity[i-1]) + + return (halfp) +end diff --git a/noao/digiphot/apphot/radprof/apgrpars.x b/noao/digiphot/apphot/radprof/apgrpars.x new file mode 100644 index 00000000..c7442121 --- /dev/null +++ b/noao/digiphot/apphot/radprof/apgrpars.x @@ -0,0 +1,46 @@ +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_GRPPARS -- Procedure to fetch the radprof parameters. + +procedure ap_grppars (ap) + +pointer ap # pointer to apphot structure + +bool clgetb() +int clgeti(), btoi() +real clgetr() + +begin + # Open the apphot strucuture. + call ap_rpinit (ap, AP_CENTROID1D, 2.5, AP_MODE, 10.0, 10.0, + 3.0, 1, 8.0, 0.5, 2.0, AP_NPOISSON) + + # Get the radial profile parameters. + call apsetr (ap, RPRADIUS, clgetr ("radius")) + call apsetr (ap, RPSTEP, clgetr ("step")) + + # Get the data dependent parameters. + call ap_gdapars (ap) + + # Get the centering algorithm parameters. + call ap_gcepars (ap) + + # Get the sky fitting algorithm parameters. + call ap_gsapars (ap) + + # Get the photometry parameters. + call ap_gphpars (ap) + + # Set remainder of the radprof parameters. + call apsetr (ap, RPKSIGMA, clgetr ("kreject")) + call apseti (ap, RPNREJECT, clgeti ("nreject")) + call apseti (ap, RPORDER, clgeti ("order")) + + # Set the plotting parameters. + call apseti (ap, RADPLOTS, btoi (clgetb ("radplots"))) +end diff --git a/noao/digiphot/apphot/radprof/approfsetup.x b/noao/digiphot/apphot/radprof/approfsetup.x new file mode 100644 index 00000000..826bdd21 --- /dev/null +++ b/noao/digiphot/apphot/radprof/approfsetup.x @@ -0,0 +1,123 @@ +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +define HELPFILE "apphot$radprof/iradprof.key" + +# AP_PROFSETUP -- Procedure to set up radprof interactively using a radial +# profile plot. + +procedure ap_profsetup (ap, im, wx, wy, gd, out, stid) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # cursor coordinates +pointer gd # pointer to graphics stream +int out # output file descriptor +int stid # output file sequence number + +int cier, sier, pier, rier, wcs, key +pointer sp, str, cmd +real xc, yc, xcenter, ycenter, rmin, rmax, imin, imax, rval +real u1, u2, v1, v2, x1, x2, y1, y2 + +int apfitcenter(), clgcur(), apfitsky(), ap_frprof(), apstati() +int ap_showplot() +real apstatr(), ap_cfwhmpsf(), ap_ccapert(), ap_cannulus(), ap_cdannulus() +real ap_crprof(), ap_crpstep(), ap_cdatamin() +real ap_cdatamax(), ap_crgrow(), ap_crclip(), ap_crclean(), ap_csigma() + +begin + if (gd == NULL) + return + call greactivate (gd, 0) + + # Save old viewport and window + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Plot the radial profile. + if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax, + imin, imax) == ERR) { + call gdeactivate (gd, 0) + return + } + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + while (clgcur ("gcommands", xc, yc, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + switch (key) { + + case 'q': + break + case '?': + call gpagefile (gd, HELPFILE, "") + case 'f': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + case 's': + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'l': + rval = ap_cdatamin (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'u': + rval = ap_cdatamax (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'c': + rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'n': + rval = ap_crclean (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'p': + rval = ap_crclip (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'a': + rval = ap_cannulus (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'd': + rval = ap_cdannulus (ap, gd, out, stid, apstatr (ap, ANNULUS), + rmin, rmax, imin, imax) + case 'g': + rval = ap_crgrow (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'r': + call ap_caper (ap, gd, out, stid, Memc[str], rmin, rmax, + imin, imax) + case 'w': + rval = ap_crprof (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'x': + rval = ap_crpstep (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'v': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_cannulus (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_cdannulus (ap, gd, out, stid, apstatr (ap, ANNULUS), + rmin, rmax, imin, imax) + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + call ap_caper (ap, gd, out, stid, Memc[str], rmin, rmax, + imin, imax) + rval = ap_crprof (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_crpstep (ap, gd, out, stid, rmin, rmax, imin, imax) + default: + call printf ("Unknown or ambiguous keystroke command\007\n") + } + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + } + call printf ( + "Interactive setup is complete. Type w to save parameters.\n") + + # Restore old view port and window + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call gdeactivate (gd, 0) + call sfree (sp) + + # Compute the answer. + cier = apfitcenter (ap, im, xcenter, ycenter) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, YCENTER), + NULL, gd) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), apstatr (ap, YCENTER), + pier) + call ap_rpplot (ap, 0, gd, apstati (ap, RADPLOTS)) + call ap_qprprof (ap, cier, sier, pier, rier) +end diff --git a/noao/digiphot/apphot/radprof/apprprof.x b/noao/digiphot/apphot/radprof/apprprof.x new file mode 100644 index 00000000..1d7461ee --- /dev/null +++ b/noao/digiphot/apphot/radprof/apprprof.x @@ -0,0 +1,109 @@ +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/photdef.h" +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_PRPROF -- Procedure to write the results of radprof to the output file. + +procedure ap_prprof (ap, fd, id, lid, cier, sier, pier, rier) + +pointer ap # pointer to apphot structure +int fd # output text file descriptor +int id # id number of str +int lid # list id of star +int cier # centering error +int sier # sky fitting error +int pier # photometric error +int rier # radial profile error + +int i, naperts +int apstati() +real apstatr() + +begin + if (fd == NULL) + return + + # Print the id parameters. + call ap_wid (ap, fd, apstatr (ap, OXINIT), apstatr (ap, OYINIT), + id, lid, '\\') + + # Print the center parameters. + call ap_wcres (ap, fd, cier, '\\') + + # Print the sky values. + call ap_wsres (ap, fd, sier, '\\') + + # Print photometry parameters. + naperts = apstati (ap, NAPERTS) + do i = 1, naperts { + if (naperts == 1) + call ap_wmres (ap, fd, i, pier, " \\") + else + call ap_wmres (ap, fd, i, pier, "*\\") + } + + # Print the radprof parameters. + call ap_wrres (ap, fd, rier) +end + + +# AP_RPHDR -- Procedure to write the radprof banner header to the output file. + +procedure ap_rphdr (ap, out) + +pointer ap # apphot descriptor +int out # output file descriptor + +begin + if (out == NULL) + return + + # Print out the keywords. + call ap_idhdr (ap, out) + call ap_chdr (ap, out) + call ap_shdr (ap, out) + call ap_mhdr (ap, out) + call ap_rhdr (ap, out) +end + + +# AP_QPRPROF -- Procedure to print a short version of the radprof results +# on the standard output. + +procedure ap_qprprof (ap, cier, sier, pier, rier) + +pointer ap # pointer to apphot structure +int cier # centering error +int sier # sky fitting error +int pier # phot error +int rier # radprof error + +pointer sp, imname, phot +real apstatr() + +begin + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + phot = AP_PPHOT(ap) + + # Print quick summary of radprof results on the standard output. + call apstats (ap, IMROOT, Memc[imname], SZ_FNAME) + call printf ("%s %8.2f %8.2f %8g %5.2f ") + call pargstr (Memc[imname]) + call pargr (apstatr (ap, ORPXCUR)) + call pargr (apstatr (ap, ORPYCUR)) + call pargr (apstatr (ap, SKY_MODE)) + call pargr (apstatr (ap, RPFWHM) / apstatr (ap, SCALE)) + call printf ("%7.3f %s\n") + call pargr (Memr[AP_MAGS(phot)+AP_NAPERTS(phot)-1]) + if (cier != AP_OK || sier != AP_OK || pier != AP_OK || rier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/apradprof.x b/noao/digiphot/apphot/radprof/apradprof.x new file mode 100644 index 00000000..10448af2 --- /dev/null +++ b/noao/digiphot/apphot/radprof/apradprof.x @@ -0,0 +1,500 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/radprof.h" + +define HELPFILE "apphot$radprof/radprof.key" + +# AP_RADPROF -- Procedure to determine radial profiles for a list of objects +# in a list of images. + +int procedure ap_radprof (ap, im, cl, gd, mgd, id, out, stid, interactive, + cache) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +pointer gd # pointer to graphcis descriptor +pointer mgd # pointer to plot metacode stream +pointer id # pointer to image display stream +int out # output file descriptor +int stid # output file sequence number +int interactive # interactive mode +int cache # cache the input image pixels + +real wx, wy, xlist, ylist +pointer sp, cmd +int wcs, key, oid, ltid, newlist, colonkey, req_size, buf_size, old_size +int newimage, newskybuf, newsky, newcenterbuf, newcenter, newbuf, newfit +int ip, prev_num, req_num, cier, sier, pier, rier, memstat + +real apstatr() +int clgcur(), apfitsky(), aprefitsky(), apfitcenter(), aprefitcenter() +int apstati(), apgscur(), ap_frprof(), ctoi(), apgqverify(), apgtverify() +int apnew(), ap_avsky(), sizeof(), ap_memstat() +bool fp_equalr() + +define endswitch_ 99 + +begin + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize the cursor command. + key = ' ' + Memc[cmd] = EOS + + # Initialize the fit. + newimage = NO + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + cier = AP_OK; sier = AP_OK; pier = AP_OK; rier = AP_OK + + # Initialize the sequencing. + newlist = NO + ltid = 0 + + # Loop over the coordinate file. + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Store the cursor coords. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Test to see if the cursor has moved. + if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + } + + # Switch on the keystroke commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("radprof", ap, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + # Print the errors. + case 'e': + if (interactive == YES) + call ap_rferrors (ap, cier, sier, pier, rier) + + # Print the help page. + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Rewind the list. + case 'r': + if (cl != NULL) { + call seek (cl, BOF) + ltid = 0 + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Move, measure next object in the coordinate list. + case 'm', 'n': + + # No coordinate file. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Need to rewind the coordinate file. + prev_num = ltid + req_num = ltid + 1 + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, + ltid) == EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + } + + # Convert coordinates if necessary. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to the next object. + newlist = YES + if (key == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + goto endswitch_ + } + + # Measure next object. + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), NULL, gd) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), pier) + call aprmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_rpplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qprprof (ap, cier, sier, pier, rier) + if (stid == 1) + call ap_param (ap, out, "radprof") + call ap_prprof (ap, out, stid, ltid, cier, sier, pier, rier) + call ap_rpplot (ap, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newbuf = NO; newfit = NO + + # Process the remainder of the list. + case 'l': + if (cl != NULL) { + oid = stid + ltid = ltid + 1 + call ap_bradprof (ap, im, cl, id, gd, mgd, out, stid, ltid, + YES) + ltid = ltid + stid - oid + 1 + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Process radprof colon commands. + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + colonkey = Memc[cmd+ip-1] + switch (colonkey) { + case 'm', 'n': + + # Show/set radprof parameters. + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call ap_rpcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newcenterbuf, newcenter, newskybuf, + newsky, newbuf, newfit) + goto endswitch_ + } + + # Process the next object. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Get next object from the list. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + + # Fetch next object from the list. + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, + ltid) == EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + + } + + # Convert the coordinates. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to next object. + newlist = YES + if (colonkey == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + } + + # Measure the next object. + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), NULL, gd) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), pier) + call aprmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_rpplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qprprof (ap, cier, sier, pier, rier) + + if (stid == 1) + call ap_param (ap, out, "radprof") + call ap_prprof (ap, out, stid, ltid, cier, sier, pier, rier) + call ap_rpplot (ap, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newbuf = NO; newfit = NO + + default: + call ap_rpcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newcenterbuf, newcenter, newskybuf, newsky, + newbuf, newfit) + } + + # Reestablish the image viewport and window. + if (newimage == YES) { + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[cmd], im, 4) + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + } + + newimage = NO + + # Save the parameters. + case 'w': + call ap_rpars (ap) + + # Plot a simple centered radial profile. + case 'd': + if (interactive == YES) { + call ap_qrad (ap, im, wx, wy, gd) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + } + + + # Setup the radial profile fitting parameters interactively. + case 'i': + if (interactive == YES) { + call ap_profsetup (ap, im, wx, wy, gd, out, stid) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + } + + # Verify the critical radprof parameters. + case 'v': + call ap_rconfirm (ap, out, stid) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newbuf = YES; newfit = YES + + # Fit the center around the current cursor value. + case 'c': + if (newcenterbuf == YES) + cier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + cier = aprefitcenter (ap, im, cier) + call aprmark (ap, id, apstati (ap, MKCENTER), NO, NO) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qcenter (ap, cier) + newcenterbuf = NO; newcenter = NO + + # Fit the sky around the current cursor value. + case 't': + if (newskybuf == YES || ! fp_equalr (wx, + apstatr (ap, SXCUR)) || ! fp_equalr (wy, apstatr (ap, + SYCUR))) + sier = apfitsky (ap, im, wx, wy, NULL, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + call aprmark (ap, id, NO, apstati (ap, MKSKY), NO) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qspsky (ap, sier) + newskybuf = NO; newsky = NO + + # Compute the average of several sky measurements around + # different cursor postions. + case 'a': + sier = ap_avsky (ap, im, stid, NULL, id, gd, interactive) + if (interactive == YES) + call ap_qaspsky (ap, sier) + newskybuf = NO; newsky = NO + + # Fit the sky around derived center value. + case 's': + if (newskybuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, SXCUR)) || ! fp_equalr (apstatr (ap, + SYCUR), apstatr (ap, YCENTER))) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), NULL, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + call aprmark (ap, id, NO, apstati (ap, MKSKY), NO) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qspsky (ap, sier) + newskybuf = NO; newsky = NO + + # Compute magnitudes around the current cursor position using + # the current sky. + case 'p', 'o': + if (newcenterbuf == YES) + cier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + cier = aprefitcenter (ap, im, cier) + call aprmark (ap, id, apstati (ap, MKCENTER), NO, + apstati (ap, MKAPERT)) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (newfit == YES || newbuf == YES || ! fp_equalr (apstatr (ap, + XCENTER), apstatr (ap, RPXCUR)) || + ! fp_equalr (apstatr (ap, RPYCUR), apstatr (ap, YCENTER))) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), pier) + if (interactive == YES) + call ap_qprprof (ap, cier, sier, pier, rier) + newcenterbuf = NO; newcenter = NO + newbuf = NO; newfit = NO + + if (key == 'o') { + if (stid == 1) + call ap_param (ap, out, "radprof") + if (newlist == YES) + call ap_prprof (ap, out, stid, ltid, cier, sier, pier, + rier) + else + call ap_prprof (ap, out, stid, 0, cier, sier, pier, + rier) + call ap_rpplot (ap, stid, mgd, YES) + stid = stid + 1 + } + + # Center, fit the sky, and compute magnitudes. + # Compute the centers, fit the sky, compute the magnitudes + # and save the results. + case 'f', ' ': + if (newcenterbuf == YES) + cier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + cier = aprefitcenter (ap, im, cier) + if (newskybuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, SXCUR)) || ! fp_equalr (apstatr (ap, YCENTER), + apstatr (ap, SYCUR))) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), NULL, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + + if (newfit == YES || newbuf == YES || ! fp_equalr (apstatr (ap, + XCENTER), apstatr (ap, RPXCUR)) || ! fp_equalr (apstatr (ap, + YCENTER), apstatr (ap, RPYCUR))) + rier = ap_frprof (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), pier) + call aprmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_rpplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qprprof (ap, cier, sier, pier, rier) + + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newbuf = NO; newfit = NO + + if (key == ' ') { + if (stid == 1) + call ap_param (ap, out, "radprof") + if (newlist == YES) + call ap_prprof (ap, out, stid, ltid, cier, sier, pier, + rier) + else + call ap_prprof (ap, out, stid, 0, cier, sier, pier, + rier) + call ap_rpplot (ap, stid, mgd, YES) + stid = stid + 1 + } + + default: + # do nothing + call printf ("Print unknown or ambiguous colon command\n") + } + +endswitch_ + # Prepare for the next object. + key = ' ' + Memc[cmd] = EOS + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/aprconfirm.x b/noao/digiphot/apphot/radprof/aprconfirm.x new file mode 100644 index 00000000..3c775e79 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprconfirm.x @@ -0,0 +1,118 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_RCONFIRM -- Procedure to confirm the critical phot parameters. + +procedure ap_rconfirm (ap, out, stid) + +pointer ap # pointer to the apphot structure +int out # output file descriptor +int stid # output file sequence number + +pointer sp, cstr, sstr, aperts +real fwhmpsf, capert, annulus, dannulus, skysigma +real datamin, datamax, radius, step +int apstati() +real apstatr(), ap_vfwhmpsf(), ap_vcapert() +real ap_vannulus(), ap_vdannulus(), ap_vsigma(), ap_vstep() +real ap_vdatamin(), ap_vdatamax(), ap_vrpradius() + +begin + call smark (sp) + call salloc (cstr, SZ_FNAME, TY_CHAR) + call salloc (sstr, SZ_FNAME, TY_CHAR) + call salloc (aperts, SZ_LINE, TY_CHAR) + + call printf ("\n") + + # Confirm the centering algorithm. + call ap_vcstring (ap, Memc[cstr], SZ_FNAME) + + if (apstati (ap, CENTERFUNCTION) != AP_NONE) { + + # Confirm the fwhmpsf. + if (apstati (ap, CENTERFUNCTION) != AP_CENTROID1D) + fwhmpsf = ap_vfwhmpsf (ap) + else + fwhmpsf = apstatr (ap, FWHMPSF) + + # Confirm the centering box. + capert = 2.0 * ap_vcapert (ap) + + } else { + fwhmpsf = apstatr (ap, FWHMPSF) + capert = 2.0 * apstatr (ap, CAPERT) + } + + # Confirm the sky fitting algorithm. + call ap_vsstring (ap, Memc[sstr], SZ_FNAME) + + if (apstati (ap, SKYFUNCTION) != AP_CONSTANT && + apstati (ap, SKYFUNCTION) != AP_SKYFILE) { + + # Confirm the sky annulus parameter. + annulus = ap_vannulus (ap) + + # Confirm the width of the sky annulus. + dannulus = ap_vdannulus (ap) + + } else { + annulus = apstatr (ap, ANNULUS) + dannulus = apstatr (ap, DANNULUS) + } + + # Confirm the sky sigma parameter. + if (apstati (ap, SKYFUNCTION) != AP_SKYFILE) + skysigma = ap_vsigma (ap) + else + skysigma = apstatr (ap, SKYSIGMA) + + # Confirm the aperture radii parameter. + call ap_vaperts (ap, Memc[aperts], SZ_LINE) + + # Confirm the radius of profile + radius = ap_vrpradius (ap) + + # Confirm the step size of profile + step = ap_vstep (ap) + + # Confirm the minimum and maximum good data values. + datamin = ap_vdatamin (ap) + datamax = ap_vdatamax (ap) + + call printf ("\n") + + # Update the database file. + if (out != NULL && stid > 1) { + call ap_sparam (out, KY_CSTRING, Memc[cstr], UN_CALGORITHM, + "centering algorithm") + call ap_rparam (out, KY_FWHMPSF, fwhmpsf, UN_ASCALEUNIT, + "full width half maximum of the psf") + call ap_rparam (out, KY_CAPERT, capert, UN_CSCALEUNIT, + "centering box width") + call ap_sparam (out, KY_SSTRING, Memc[sstr], UN_SALGORITHM, + "sky fitting algorithm") + call ap_rparam (out, KY_ANNULUS, annulus, UN_SSCALEUNIT, + "inner radius of the sky annulus") + call ap_rparam (out, KY_DANNULUS, dannulus, UN_SSCALEUNIT, + "width of the sky annulus") + call ap_rparam (out, KY_SKYSIGMA, skysigma, UN_NCOUNTS, + "standard deviation of 1 sky pixel") + call ap_sparam (out, KY_APERTS, Memc[aperts], UN_PSCALEUNIT, + "list of apertures") + call ap_rparam (out, KY_DATAMIN, datamin, UN_ACOUNTS, + "minimum good data value") + call ap_rparam (out, KY_DATAMAX, datamax, UN_ACOUNTS, + "maximum good data value") + call ap_rparam (out, KY_RPRADIUS, radius, UN_RSCALEUNIT, + "fitting radius") + call ap_rparam (out, KY_RPSTEP, step, UN_RSCALEUNIT, + "step size in radius") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/aprferrors.x b/noao/digiphot/apphot/radprof/aprferrors.x new file mode 100644 index 00000000..f72586f3 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprferrors.x @@ -0,0 +1,40 @@ +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_RFERRORS -- Procedure to print a short form of the output of radprof on +# the standard output. + +procedure ap_rferrors (ap, cier, sier, pier, rier) + +pointer ap # pointer to apphot structure +int cier # centering error +int sier # sky fitting error +int pier # photmetry error +int rier # photometric error + +begin + # Print the centering errors. + call ap_cerrors (ap, cier) + + # Print the sky fitting errors. + call ap_serrors (ap, sier) + + # Print the photometry errors. + call ap_merrors (ap, pier) + + # Print the radial profile fitting errors. + switch (rier) { + case AP_RP_NOPROFILE: + call printf ("The profile fitting region is outside the image.\n") + case AP_RP_OUTOFBOUNDS: + call printf ( + "The profile fitting region is partially outside the image.\n") + case AP_RP_NPTS_TOO_SMALL: + call printf ( + "There are too few points in the profile fitting region.\n") + case AP_RP_SINGULAR: + call printf ("The profile fit is singular.\n") + default: + call printf ("") + } +end diff --git a/noao/digiphot/apphot/radprof/aprmmeasure.x b/noao/digiphot/apphot/radprof/aprmmeasure.x new file mode 100644 index 00000000..fe455008 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprmmeasure.x @@ -0,0 +1,96 @@ +include <math/curfit.h> + +# AP_RMMEASURE -- Procedure to compute the sums and areas inside the apertures. + +procedure ap_rmmeasure (pixels, nx, ny, wx, wy, aperts, sums, areas, naperts) + +real pixels[nx,ARB] # subraster pixel values +int nx, ny # dimensions of the subraster +real wx, wy # center of subraster +real aperts[ARB] # array of apertures +double sums[ARB] # array of sums +double areas[ARB] # aperture areas +int naperts # number of apertures + +int i, j, k +double fctn +real apmaxsq, dy2, r2, r + +begin + # Initialize. + apmaxsq = (aperts[naperts] + 0.5) ** 2 + call aclrd (sums, naperts) + call aclrd (areas, naperts) + + # Loop over the pixels. + do j = 1, ny { + dy2 = (j - wy) ** 2 + do i = 1, nx { + r2 = (i - wx) ** 2 + dy2 + if (r2 > apmaxsq) + next + r = sqrt (r2) - 0.5 + do k = 1, naperts { + if (r > aperts[k]) + next + fctn = max (0.0, min (1.0, aperts[k] - r)) + sums[k] = sums[k] + fctn * pixels[i,j] + areas[k] = areas[k] + fctn + } + } + } +end + + +# AP_BRMMEASURE -- Procedure to compute the sums and areas inside the apertures. + +procedure ap_brmmeasure (pixels, nx, ny, wx, wy, datamin, datamax, aperts, + sums, areas, naperts, minapert) + +real pixels[nx,ARB] # subraster pixel values +int nx, ny # dimensions of the subraster +real wx, wy # center of subraster +real datamin # minimum good data value +real datamax # maximum good data value +real aperts[ARB] # array of apertures +double sums[ARB] # array of sums +double areas[ARB] # aperture areas +int naperts # number of apertures +int minapert # minimum number of apertures + +int i, j, k, kindex +double fctn +real apmaxsq, dy2, r2, r, pixval + +begin + # Initialize. + apmaxsq = (aperts[naperts] + 0.5) ** 2 + call aclrd (sums, naperts) + call aclrd (areas, naperts) + minapert = naperts + 1 + + # Loop over the pixels. + do j = 1, ny { + dy2 = (j - wy) ** 2 + do i = 1, nx { + r2 = (i - wx) ** 2 + dy2 + if (r2 > apmaxsq) + next + r = sqrt (r2) - 0.5 + kindex = naperts + 1 + pixval = pixels[i,j] + do k = 1, naperts { + if (r > aperts[k]) + next + kindex = min (k, kindex) + fctn = max (0.0, min (1.0, aperts[k] - r)) + sums[k] = sums[k] + fctn * pixval + areas[k] = areas[k] + fctn + } + if (kindex < minapert) { + if (pixval < datamin || pixval > datamax) + minapert = kindex + } + } + } +end diff --git a/noao/digiphot/apphot/radprof/aprpars.x b/noao/digiphot/apphot/radprof/aprpars.x new file mode 100644 index 00000000..fb4ed401 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpars.x @@ -0,0 +1,37 @@ +include "../lib/display.h" +include "../lib/radprof.h" + +# AP_RPARS -- Procedure to write the radprof parameters to the parameter +# file. + +procedure ap_rpars (ap) + +pointer ap # pointer to apphot structure + +bool itob() +int apstati() +real apstatr() + +begin + # Write the data dependent parameters. + call ap_dapars (ap) + + # Write the centering parameters. + call ap_cepars (ap) + + # Write the sky fitting parameters. + call ap_sapars (ap) + + # Write out the photometry parameters. + call ap_phpars (ap) + + # Set the radphot parameters + call clputr ("radius", apstatr (ap, RPRADIUS)) + call clputr ("step", apstatr (ap, RPSTEP)) + call clputi ("order", apstati (ap, RPORDER)) + call clputr ("kreject", apstatr (ap, RPKSIGMA)) + call clputi ("nreject", apstati (ap, RPNREJECT)) + + # Set radial profile plots + call clputb ("radplots", itob (apstati (ap, RADPLOTS))) +end diff --git a/noao/digiphot/apphot/radprof/aprpbuf.x b/noao/digiphot/apphot/radprof/aprpbuf.x new file mode 100644 index 00000000..3a2f73ed --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpbuf.x @@ -0,0 +1,82 @@ +include <imhdr.h> +include "../lib/apphotdef.h" +include "../lib/radprofdef.h" +include "../lib/radprof.h" + +# AP_RPBUF -- Procedure to determine the mapping of the of the radial +# profile size into the apertures. + +int procedure ap_rpbuf (ap, im, wx, wy) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # center coordinates + +int c1, c2, l1, l2 +pointer rprof +real rbuf +pointer ap_rppix() + +begin + # Check for 0 radius aperture. + rprof = AP_RPROF(ap) + if (AP_RPRADIUS(rprof) <= 0.0) + return (AP_RP_NOPROFILE) + + # Compute the maximum aperture size + rbuf = 2. * AP_RPRADIUS(rprof) * AP_SCALE(ap) + 1. + AP_RPIX(rprof) = ap_rppix (im, wx, wy, rbuf, c1, c2, l1, l2) + AP_RPXC(rprof) = wx - c1 + 1 + AP_RPYC(rprof) = wy - l1 + 1 + AP_RPNX(rprof) = c2 - c1 + 1 + AP_RPNY(rprof) = l2 - l1 + 1 + + # Return the appropriate error code. + if (AP_RPIX(rprof) == NULL) { + return (AP_RP_NOPROFILE) + } else if (AP_RPNX(rprof) < rbuf || AP_RPNY(rprof) < rbuf) { + return (AP_RP_OUTOFBOUNDS) + } else { + return (AP_OK) + } +end + + +# AP_RPPIX -- Procedure to read in the aperture pixels + +pointer procedure ap_rppix (im, wx, wy, papert, c1, c2, l1, l2) + +pointer im # pointer to IRAF image +real wx, wy # center of centering subraster annulus +real papert # centering radius +int c1, c2 # column limits +int l1, l2 # line limits + +int ncols, nlines +real half_papert, xc1, xc2, xl1, xl2 +pointer imgs2r() + +begin + # Check for 0 radius aperture. + half_papert = papert / 2. + if (half_papert <= 0.) + return (NULL) + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + + # Test for an out of bounds aperture. + xc1 = wx - half_papert + xc2 = wx + half_papert + xl1 = wy - half_papert + xl2 = wy + half_papert + if (xc1 > real (ncols) || xc2 < 1.0 || xl1 > real (nlines) || xl2 < 1.0) + return (NULL) + + # Get the column and line limits, dimensions and center of the subraster + # to be extracted. + c1 = max (1.0, min (real (ncols), xc1)) + c2 = min (real (ncols), max (1.0, xc2)) + l1 = max (1.0, min (real (nlines), xl1)) + l2 = min (real (nlines), max (1.0, xl2)) + return (imgs2r (im, c1, c2, l1, l2)) +end diff --git a/noao/digiphot/apphot/radprof/aprpcolon.x b/noao/digiphot/apphot/radprof/aprpcolon.x new file mode 100644 index 00000000..a5bf0eab --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpcolon.x @@ -0,0 +1,241 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" +include "../lib/radprof.h" +include "../lib/display.h" + +# AP_RPCOLON -- Show/set radprof parameters. + +procedure ap_rpcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + newcenterbuf, newcenter, newskybuf, newsky, newbuf, newfit) + +pointer ap # pointer to the apphot structure +pointer im # pointer to iraf image +int cl # coord file descriptor +int out # output file descriptor +int stid # output file sequence number +int ltid # coord list sequence number +char cmdstr[ARB] # command string +int newimage # new image +int newcenterbuf, newcenter # new sky fit +int newskybuf, newsky # new sky buffer +int newbuf, newfit # new aperture + +pointer sp, incmd, outcmd +int strdic() + +begin + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + # Get the commands. + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, CCMDS) != 0) + call apccolon (ap, out, stid, cmdstr, newcenterbuf, newcenter) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, SCMDS) != 0) + call apscolon (ap, out, stid, cmdstr, newskybuf, newsky) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, RPCMDS) != 0) + call ap_profcolon (ap, out, stid, cmdstr, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, PCMDS) != 0) + call apmagcolon (ap, out, stid, cmdstr, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0) + call ap_apcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + newcenterbuf, newcenter, newskybuf, newsky, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0) + call ap_nscolon (ap, im, out, stid, cmdstr, newcenterbuf, + newcenter, newskybuf, newsky, newbuf, newfit) + else + call ap_rpimcolon (ap, cmdstr) + + call sfree (sp) +end + + +# AP_PROFCOLON -- Procedure to display and modify radprof parameters. + +procedure ap_profcolon (ap, out, stid, cmdstr, newbuf, newfit) + +pointer ap # pointer to apphot structure +int out # output file descriptor +int stid # output file number +char cmdstr[ARB] # command string +int newbuf # new aperture buffers +int newfit # compute new magnitudes + +int ival, ncmd +pointer sp, cmd +real rval +int strdic(), nscan(), apstati() +real apstatr() + +begin + # Get the command. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the colon command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, RPCMDS) + switch (ncmd) { + case RCMD_RADIUS: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RPRADIUS) + call pargr (apstatr (ap, RPRADIUS)) + call pargstr (UN_RSCALEUNIT) + } else { + call apsetr (ap, RPRADIUS, rval) + if (stid > 1) + call ap_rparam (out, KY_RPRADIUS, rval, UN_RSCALEUNIT, + "fitting radius") + newbuf = YES; newfit = YES + } + case RCMD_STEPSIZE: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RPSTEP) + call pargr (apstatr (ap, RPSTEP)) + call pargstr (UN_RSCALEUNIT) + } else { + call apsetr (ap, RPSTEP, rval) + if (stid > 1) + call ap_rparam (out, KY_RPSTEP, rval, UN_RSCALEUNIT, + "step size in radius") + newfit = YES + } + case RCMD_ORDER: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_RPORDER) + call pargi (apstati (ap, RPORDER)) + } else { + call apseti (ap, RPORDER, ival) + if (stid > 1) + call ap_iparam (out, KY_RPORDER, ival, UN_RNUMBER, + "maximum number of rejection cycels") + newfit = YES + } + case RCMD_KREJECT: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RPKSIGMA) + call pargr (apstatr (ap, RPKSIGMA)) + call pargstr (UN_RSIGMA) + } else { + call apsetr (ap, RPKSIGMA, rval) + if (stid > 1) + call ap_rparam (out, KY_RPKSIGMA, rval, UN_RSIGMA, + "k-sigma rejection criteron") + newfit = YES + } + case RCMD_NREJECT: + call gargi (ival) + if (nscan() == 1) { + call printf ("%s = %d\n") + call pargstr (KY_RPNREJECT) + call pargi (apstati (ap, RPNREJECT)) + } else { + call apseti (ap, RPNREJECT, ival) + if (stid > 1) + call ap_iparam (out, KY_RPNREJECT, ival, UN_RNUMBER, + "maximum number of rejection cycles") + newfit = YES + } + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + call sfree (sp) +end + + +# AP_RPIMCOLON -- Show/set quantities which are not radprof parameters. + +procedure ap_rpimcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr[ARB] # command string + +bool bval +int ncmd +pointer sp, cmd +bool itob() +int apstati(), strdic(), nscan(), btoi() + +begin + call smark (sp) + call salloc (cmd, 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, MISC) + switch (ncmd) { + case ACMD_SHOW: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, RPSHOWARGS) + switch (ncmd) { + case RCMD_CENTER: + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + case RCMD_SKY: + call printf ("\n") + call ap_spshow (ap) + call printf ("\n") + case RCMD_PHOT: + call printf ("\n") + call ap_mpshow (ap) + call printf ("\n") + case RCMD_FIT: + call printf ("\n") + call ap_rppshow (ap) + call printf ("\n") + case RCMD_DATA: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_rprofshow (ap) + call printf ("\n") + } + case ACMD_RADPLOTS: + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) + } else + call apseti (ap, RADPLOTS, btoi (bval)) + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/aprpfree.x b/noao/digiphot/apphot/radprof/aprpfree.x new file mode 100644 index 00000000..2f16e03a --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpfree.x @@ -0,0 +1,60 @@ +include "../lib/apphotdef.h" +include "../lib/radprofdef.h" + +# AP_RPFREE -- Procedure to free the radial profile fitting structure. + +procedure ap_rpfree (ap) + +pointer ap # pointer to the apphot structure + +begin + if (ap == NULL) + return + if (AP_NOISE(ap) != NULL) + call ap_noisecls (ap) + if (AP_PCENTER(ap) != NULL) + call ap_ctrcls (ap) + if (AP_PDISPLAY(ap) != NULL) + call ap_dispcls (ap) + if (AP_POLY(ap) != NULL) + call ap_ycls (ap) + if (AP_PPHOT(ap) != NULL) + call ap_photcls (ap) + if (AP_PPSF(ap) != NULL) + call ap_psfcls (ap) + if (AP_RPROF(ap) != NULL) + call ap_rpcls (ap) + if (AP_PSKY(ap) != NULL) + call ap_skycls (ap) + if (AP_IMBUF(ap) != NULL) + call mfree (AP_IMBUF(ap), TY_REAL) + if (AP_MW(ap) != NULL) + call mw_close (AP_MW(ap)) + call mfree (ap, TY_STRUCT) +end + + +# AP_RPCLS -- Procedure to closee up the radial profile fitting arrays. + +procedure ap_rpcls (ap) + +pointer ap # pointer to apphot structure + +pointer rprof + +begin + rprof = AP_RPROF(ap) + if (rprof == NULL) + return + #if (AP_RPIX(rprof) != NULL) + #call mfree (AP_RPIX(rprof), TY_REAL) + if (AP_RPDIST(rprof) != NULL) + call mfree (AP_RPDIST(rprof), TY_REAL) + if (AP_INTENSITY(rprof) != NULL) + call mfree (AP_INTENSITY(rprof), TY_REAL) + if (AP_DINTENSITY(rprof) != NULL) + call mfree (AP_DINTENSITY(rprof), TY_REAL) + if (AP_TINTENSITY(rprof) != NULL) + call mfree (AP_TINTENSITY(rprof), TY_REAL) + call mfree (rprof, TY_STRUCT) +end diff --git a/noao/digiphot/apphot/radprof/aprpindef.x b/noao/digiphot/apphot/radprof/aprpindef.x new file mode 100644 index 00000000..4d92d615 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpindef.x @@ -0,0 +1,69 @@ +include "../lib/apphotdef.h" +include "../lib/radprofdef.h" +include "../lib/photdef.h" +include "../lib/phot.h" +include "../lib/radprof.h" + +# AP_RPINDEF -- Routine to return INDEF valued photometry and radial profile +# buffers. + +procedure ap_rpindef (ap) + +pointer ap # pointer to the apphot structure + +pointer phot, rprof + +begin + phot = AP_PPHOT(ap) + rprof = AP_RPROF(ap) + + AP_RPFWHM(rprof) = INDEFR + AP_INORM(rprof) = INDEFR + AP_TINORM(rprof) = INDEFR + + call amovkr (INDEFR, Memr[AP_INTENSITY(rprof)], AP_RPNPTS(rprof)) + call amovkr (INDEFR, Memr[AP_TINTENSITY(rprof)], AP_RPNPTS(rprof)) + call amovkd (0.0d0, Memd[AP_AREA(phot)], AP_NAPERTS(phot)) + call amovkd (0.0d0, Memd[AP_SUMS(phot)], AP_NAPERTS(phot)) + call amovkr (INDEFR, Memr[AP_MAGS(phot)], AP_NAPERTS(phot)) + call amovkr (INDEFR, Memr[AP_MAGERRS(phot)], AP_NAPERTS(phot)) +end + + +# AP_MAXAP -- Procedure to setup the maximum number of apertures for phot. + +procedure ap_maxap (ap, pier) + +pointer ap # pointer to the apphot structure +int pier # photometric error + +int i +pointer phot, rprof +real dxc1, dxc2, dyc1, dyc2, rdist, rapert + +begin + phot = AP_PPHOT(ap) + rprof = AP_RPROF(ap) + + dxc1 = AP_RPXC(rprof) - 0.5 + dxc2 = AP_RPNX(rprof) - AP_RPXC(rprof) + 0.5 + dyc1 = AP_RPYC(rprof) - 0.5 + dyc2 = AP_RPNY(rprof) - AP_RPYC(rprof) + 0.5 + + # Compute the maximum aperture. + AP_NMAXAP(phot) = 0 + rdist = min (abs (dxc1), abs (dxc2), abs (dyc1), abs (dyc2)) + do i = 1, AP_NAPERTS(phot) { + rapert = AP_SCALE(ap) * Memr[AP_APERTS(phot)+i-1] + if (rapert <= rdist) { + AP_NMAXAP(phot) = i + } else { + break + } + } + + if (AP_NMAXAP(phot) < AP_NAPERTS(phot)) + pier = AP_APERT_OUTOFBOUNDS + else + pier = AP_OK +end diff --git a/noao/digiphot/apphot/radprof/aprpinit.x b/noao/digiphot/apphot/radprof/aprpinit.x new file mode 100644 index 00000000..49bea62c --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpinit.x @@ -0,0 +1,77 @@ +include "../lib/apphotdef.h" +include "../lib/radprofdef.h" +include "../lib/phot.h" + +# AP_RPINIT - Procedure to initialize the radial profile fitting structure. + +procedure ap_rpinit (ap, cfunction, cbox, sfunction, annulus, dannulus, + aperts, napert, radius, step, fwhmpsf, noise) + +pointer ap # pointer to the apphot structure +int cfunction # centering algorithm +real cbox # half width of the centering box +int sfunction # sky fitting algorithm +real annulus # radius of sky annulus +real dannulus # width of sky annulus +real aperts[ARB] # array of apertures +int napert # number of apertures +real radius # radius of fitting region +real step # step size of output +real fwhmpsf # FWHM of the PSF +int noise # Noise model + +begin + # Set the image dependent parameters. + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set up the global apphot package defaults. + call ap_defsetup (ap, fwhmpsf) + + # Set up the noise model parameters. + call ap_noisesetup (ap, noise) + + # Set up the centering algorithm parameters. + call ap_ctrsetup (ap, cfunction, cbox) + + # Set up the sky fitting parameters. + call ap_skysetup (ap, sfunction, annulus, dannulus) + + # Set up the photometry parameters. + call ap_photsetup (ap, aperts, napert, AP_PWCONSTANT) + + # Set up the radial profile fitting parameters. + call ap_rpsetup (ap, radius, step) + + # Set up the display options. + call ap_dispsetup (ap) + + # Set psf fitting and polyphot structures to null. + AP_PPSF(ap) = NULL + AP_POLY(ap) = NULL +end + + +# AP_RPSETUP -- Procedure to set up the radial profle fitting parameters. + +procedure ap_rpsetup (ap, radius, step) + +pointer ap # pointer to apphot structure +real radius # radius of psf to be fit +real step # step size + +pointer rprof + +begin + call malloc (AP_RPROF(ap), LEN_RPSTRUCT, TY_STRUCT) + rprof = AP_RPROF(ap) + AP_RPXCUR(rprof) = INDEFR + AP_RPYCUR(rprof) = INDEFR + AP_RPRADIUS(rprof) = radius + AP_RPSTEP(rprof) = step + AP_RPIX(rprof) = NULL + AP_RPNPTS(rprof) = int (AP_RPRADIUS(rprof) / AP_RPSTEP(rprof)) + 1 + call malloc (AP_RPDIST(rprof), AP_RPNPTS(rprof), TY_REAL) + call malloc (AP_INTENSITY(rprof), AP_RPNPTS(rprof), TY_REAL) + call malloc (AP_DINTENSITY(rprof), AP_RPNPTS(rprof), TY_REAL) + call malloc (AP_TINTENSITY(rprof), AP_RPNPTS(rprof), TY_REAL) +end diff --git a/noao/digiphot/apphot/radprof/aprpplot.x b/noao/digiphot/apphot/radprof/aprpplot.x new file mode 100644 index 00000000..dee441f1 --- /dev/null +++ b/noao/digiphot/apphot/radprof/aprpplot.x @@ -0,0 +1,307 @@ +include <gset.h> +include <pkg/gtools.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/photdef.h" +include "../lib/phot.h" +include "../lib/radprofdef.h" +include "../lib/radprof.h" + +define IMIN -0.1 # Minimum intensity value for plot +define IMAX 1.1 # Maximum intensity value for plot + +# AP_RPPLOT -- Procedure to plot the radial profile. + +procedure ap_rpplot (ap, sid, gd, makeplots) + +pointer ap # pointer to the apphot structure +int sid # output file id number (not used) +pointer gd # pointer to the plot stream +int makeplots # make plots on the screen ? + +int nxpts, nypts, nrpts +pointer rprof, gt +real rmin, rmax, inorm, x1, x2, y1, y2, u1, u2, v1, v2 +int apstati() +pointer ap_gtinit() +real apstatr() + +begin + # Return if no graphics stream. + if (gd == NULL || makeplots == NO) + return + + # Return if the center or sky is undefined. + if (IS_INDEFR(apstatr (ap, XCENTER)) || IS_INDEFR(apstatr (ap, + YCENTER)) || IS_INDEFR (apstatr (ap, SKY_MODE))) + return + + # Return if there are no pixels. + rprof = AP_RPROF(ap) + if (AP_RPIX(rprof) == NULL) + return + + # Set up some useful constants. + rmin = 0.0 + rmax = apstatr (ap, SCALE) * apstatr (ap, RPRADIUS) + nxpts = AP_RPNX(rprof) + nypts = AP_RPNY(rprof) + nrpts = apstati (ap, RPNPTS) + inorm = apstatr (ap, INORM) + + # Reopen the work station. + call greactivate (gd, 0) + + # Save old viewport and coordinates. + call ggwind (gd, x1, x2, y1, y2) + call ggview (gd, u1, u2, v1, v2) + + # Set up the labels and annotate the plot. + gt = ap_gtinit (AP_IMROOT(ap), apstatr (ap, OXINIT), apstatr (ap, + OYINIT)) + call gclear (gd) + call ap_rpset (gd, gt, ap, rmin, rmax, IMIN, IMAX) + call ap_rpannotate (gd, ap, rmin, rmax, IMIN, IMAX) + + # Plot the intensity and total intensity. + call ap_plothist (gd, gt, Memr[AP_RPDIST(rprof)], + Memr[AP_INTENSITY(rprof)], nrpts, "line", GL_SOLID) + call ap_plothist (gd, gt, Memr[AP_RPDIST(rprof)], + Memr[AP_TINTENSITY(rprof)], nrpts, "line", GL_DOTDASH) + + # Plot the points. + call gswind (gd, rmin, rmax, (IMIN * inorm), (IMAX * inorm)) + call rp_ptsplot (gd, gt, Memr[AP_RPIX(rprof)], nxpts, nypts, + AP_RPXC(rprof), AP_RPYC(rprof)) + + # Restore old viewport and world coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call ap_gtfree (gt) + call gdeactivate (gd, 0) +end + + +# AP_RPSET -- Procedure to set up the parameters for the radial profile +# plot. + +procedure ap_rpset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # graphics stream +pointer gt # gtools pointer +pointer ap # apphot pointer +real xmin, xmax # minimum and maximum radial distance +real ymin, ymax # minimum and maximum of the y axis + +int fd, naperts +pointer sp, str, tstr, temp +real scale, aspect, vx1, vx2, vy1, vy2 +int stropen(), apstati() +real apstatr(), gstatr() + +begin + call smark (sp) + call salloc (str, 6 * SZ_LINE, TY_CHAR) + call salloc (tstr, SZ_LINE, TY_CHAR) + naperts = apstati (ap, NAPERTS) + call salloc (temp, naperts, TY_REAL) + + # Open the title string. + fd = stropen (Memc[str], 6 * SZ_LINE, WRITE_ONLY) + + # Encode the sysid. + call sysid (Memc[tstr], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[tstr]) + + # Encode the center string + call fprintf (fd, + "Center: xc=%0.2f yc=%0.2f xerr=%0.2f yerr=%0.2f\n") + call pargr (apstatr (ap, OXCENTER)) + call pargr (apstatr (ap, OYCENTER)) + call pargr (apstatr (ap, XERR)) + call pargr (apstatr (ap, YERR)) + + # Encode the sky string + call fprintf (fd, + "Sky: value=%0.2f sigma=%0.2f skew=%0.2f nsky=%d nreject=%d\n") + call pargr (apstatr (ap, SKY_MODE)) + call pargr (apstatr (ap, SKY_SIGMA)) + call pargr (apstatr (ap, SKY_SKEW)) + call pargi (apstati (ap, NSKY)) + call pargi (apstati (ap, NSKY_REJECT)) + + # Encode the value of the magnitude at the maximum aperture. + call fprintf (fd, "Photometry: fwhmpsf=%0.3f maxapert=") + call pargr (apstatr (ap, RPFWHM)) + call ap_arrayr (ap, APERTS, Memr[temp]) + call amulkr (Memr[temp], apstatr (ap, SCALE), Memr[temp], naperts) + call fprintf (fd, "%0.2f mags=") + call pargr (Memr[temp+naperts-1]) + call ap_arrayr (ap, MAGS, Memr[temp]) + call fprintf (fd, "%0.3f\n") + call pargr (Memr[temp+naperts-1]) + + # Encode the parameter string. + call fprintf (fd, + "Fit: spline3 order=%d krej=%0.1f sigma np=%d nprej=%d\n") + call pargi (apstati (ap, RPORDER)) + call pargr (apstatr (ap, RPKSIGMA)) + call pargi (apstati (ap, RPNDATA)) + call pargi (apstati (ap, RPNDATAREJ)) + + # Encode the title. + call gt_gets (gt, GTTITLE, Memc[tstr], SZ_LINE) + call fprintf (fd, "%s\n\n") + call pargstr (Memc[tstr]) + + call strclose (fd) + + aspect = gstatr (gd, G_ASPECT) + scale = apstatr (ap, SCALE) + call gsetr (gd, G_ASPECT, 0.70) + call gseti (gd, G_WCS, 2) + + # Draw and label the axes. + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, ymin, ymax) + call glabax (gd, Memc[str], "", "Intensity") + call gseti (gd, G_YDRAWAXES, 0) + call gseti (gd, G_XDRAWAXES, 1) + call ggview (gd, vx1, vx2, vy1, vy2) + call gsview (gd, vx1, vx2, vy1, vy2) + call gswind (gd, xmin, xmax, ymin, ymax) + call glabax (gd, "", + "Radial Distance (lower-pixels, upper-scale units)", "") + + # Reset the axes parameters. + call gseti (gd, G_YDRAWAXES, 3) + call gseti (gd, G_XDRAWAXES, 3) + call gsetr (gd, G_ASPECT, aspect) + + # Set the plot type. + call gt_sets (gt, GTTYPE, "line") + + call sfree (sp) +end + + +# AP_RPANNOTATE -- Procedure to annotate the radial plot in radprof. + +procedure ap_rpannotate (gd, ap, xmin, xmax, ymin, ymax) + +pointer gd # graphics stream +pointer ap # apphot structure +real xmin, xmax # min and max of x axis +real ymin, ymax # min and max of y axis + +int naperts +pointer sp, str, temp +real scale, rpfwhm, annulus, dannulus, inorm, tinorm +int apstati() +real apstatr() + +begin + # Allocate working space. + naperts = apstati (ap, NAPERTS) + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (temp, naperts, TY_REAL) + call gseti (gd, G_PLTYPE, GL_DASHED) + + # Draw the zero level line. + call gamove (gd, xmin, 0.0) + call gadraw (gd, xmax, 0.0) + + # Draw the half power point. + call gamove (gd, xmin, 0.5) + call gadraw (gd, xmax, 0.5) + + # Draw the unit normalization value. + call gamove (gd, xmin, 1.0) + call gadraw (gd, xmax, 1.0) + + # Plot the full width half maximum of the radial profile. + scale = apstatr (ap, SCALE) + rpfwhm = apstatr (ap, RPFWHM) / 2.0 + if (rpfwhm >= xmin && rpfwhm <= xmax) { + call gamove (gd, rpfwhm, ymin) + call gadraw (gd, rpfwhm, ymax) + call sprintf (Memc[str], SZ_LINE, "hwhm = %0.2f") + call pargr (rpfwhm) + call gtext (gd, rpfwhm, 0.0, Memc[str], "q=h;u=180;p=r") + } + + # Mark the sky annuli. + annulus = scale * apstatr (ap, ANNULUS) + dannulus = scale * (apstatr (ap, ANNULUS) + apstatr (ap, DANNULUS)) + if (annulus >= xmin && annulus <= xmax) { + call gamove (gd, annulus, ymin) + call gadraw (gd, annulus, ymax) + call sprintf (Memc[str], SZ_LINE, "inner sky radius = %0.2f") + call pargr (annulus) + call gtext (gd, annulus, 0.0, Memc[str], "q=h;u=180;p=r") + } + if (dannulus >= xmin && dannulus <= xmax) { + call gamove (gd, dannulus, ymin) + call gadraw (gd, dannulus, ymax) + call sprintf (Memc[str], SZ_LINE, "outer sky radius = %0.2f") + call pargr (dannulus) + call gtext (gd, dannulus, 0.0, Memc[str], "q=h;u=180;p=r") + } + + # Plot the aperture value. + call ap_arrayr (ap, APERTS, Memr[temp]) + call amulkr (Memr[temp], scale, Memr[temp], naperts) + call gseti (gd, G_PLTYPE, GL_SOLID) + if (Memr[temp+naperts-1] >= xmin && Memr[temp+naperts-1] <= xmax) { + call gamove (gd, Memr[temp+naperts-1], ymin) + call gadraw (gd, Memr[temp+naperts-1], ymax) + call sprintf (Memc[str], SZ_LINE, "maxapert = %0.2f") + call pargr (Memr[temp+naperts-1]) + call gtext (gd, Memr[temp+naperts-1], 0.0, Memc[str], + "q=h;u=180;p=r") + } + call gseti (gd, G_PLTYPE, GL_DASHED) + + # Plot the inorm value. + inorm = apstatr (ap, INORM) + call sprintf (Memc[str], SZ_LINE, "inorm = %0.2f") + call pargr (inorm) + call gtext (gd, 0.0, 1.0, Memc[str], "q=h") + + # Plot the tinorm value. + tinorm = apstatr (ap, TNORM) + call sprintf (Memc[str], SZ_LINE, "tinorm = %0.2f") + call pargr (tinorm) + call gtext (gd, xmax, 1.0, Memc[str], "q=h;h=r") + + call sfree (sp) +end + + +# RP_PTSPLOT -- Plot the radial profile plots. + +procedure rp_ptsplot (gd, gt, pixels, nx, ny, wx, wy) + +pointer gd # pointer to the graphics stream +pointer gt # pointer to the gtools structure +real pixels[nx,ARB] # subraster of pixel values +int nx, ny # dimensions of the pixel subraster +real wx, wy # x and y coordinates of the center + +int i +pointer sp, rtemp + +begin + call smark (sp) + call salloc (rtemp, nx, TY_REAL) + do i = 1, ny { + call ap_ijtor (Memr[rtemp], nx, i, wx, wy) + call ap_plotrad (gd, gt, Memr[rtemp], pixels[1,i], nx, "plus") + } + call sfree (sp) +end diff --git a/noao/digiphot/apphot/radprof/iradprof.key b/noao/digiphot/apphot/radprof/iradprof.key new file mode 100644 index 00000000..1c5a8974 --- /dev/null +++ b/noao/digiphot/apphot/radprof/iradprof.key @@ -0,0 +1,20 @@ + Interactive Radprof Setup Menu + + v Mark and verify the critical parameters (f,c,s,a,d,r,w,x) + + f Mark and verify the psf full-width half-maximum + s Mark and verify the standard deviation of the background + l Mark and verify the minimum good data value + u Mark and verify the maximum good data value + + c Mark and verify the centering box width + n Mark and verify the cleaning radius + p Mark and verify the clipping radius + + a Mark and verify the inner radius of the sky annulus + d Mark and verify the width of the sky annulus + g Mark and verify the region growing radius + + r Mark and verify the photometry aperture radii + w Mark and verify the radius of the radial profile + x Mark and verify the step size of radial profile diff --git a/noao/digiphot/apphot/radprof/mkpkg b/noao/digiphot/apphot/radprof/mkpkg new file mode 100644 index 00000000..8ef786bc --- /dev/null +++ b/noao/digiphot/apphot/radprof/mkpkg @@ -0,0 +1,58 @@ +RADPROF Fitting Task Routines + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + apbradprof.x <fset.h> ../lib/apphot.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/display.h + apgrpars.x ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h ../lib/display.h \ + ../lib/noise.h ../lib/radprof.h + approfsetup.x ../lib/center.h ../lib/fitsky.h \ + ../lib/display.h + apprprof.x ../lib/apphotdef.h ../lib/apphot.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/photdef.h ../lib/phot.h \ + ../lib/radprof.h + apradprof.x <ctype.h> <gset.h> \ + ../lib/apphot.h ../lib/center.h \ + ../lib/fitsky.h ../lib/radprof.h \ + ../lib/display.h <imhdr.h> + aprconfirm.x ../lib/apphot.h ../lib/noise.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h ../lib/radprof.h + aprferrors.x ../lib/phot.h ../lib/radprof.h + aprmmeasure.x <math/curfit.h> + aprpars.x ../lib/radprof.h ../lib/display.h + aprpbuf.x <imhdr.h> ../lib/apphotdef.h \ + ../lib/radprofdef.h ../lib/radprof.h + aprpcolon.x ../lib/apphot.h ../lib/fitsky.h \ + ../lib/center.h ../lib/phot.h \ + ../lib/radprof.h ../lib/noise.h \ + ../lib/display.h + aprpfree.x ../lib/apphotdef.h ../lib/radprofdef.h + aprpindef.x ../lib/apphotdef.h ../lib/radprofdef.h \ + ../lib/radprof.h ../lib/photdef.h \ + ../lib/phot.h + aprpinit.x ../lib/apphotdef.h ../lib/radprofdef.h \ + ../lib/phot.h + aprpplot.x <gset.h> <pkg/gtools.h> \ + ../lib/apphotdef.h ../lib/photdef.h \ + ../lib/radprof.h ../lib/phot.h \ + ../lib/fitsky.h ../lib/radprofdef.h \ + ../lib/center.h ../lib/apphot.h + apfrprof.x <math/curfit.h> <math/iminterp.h> \ + <gset.h> ../lib/fitskydef.h \ + ../lib/apphotdef.h ../lib/radprofdef.h \ + ../lib/phot.h ../lib/photdef.h \ + ../lib/radprof.h <math.h> \ + <mach.h> "../lib/noisedef.h" \ + ../lib/apphot.h + rprofshow.x ../lib/display.h ../lib/radprof.h + t_radprof.x <fset.h> <gset.h> \ + ../lib/apphot.h <imhdr.h> + ; diff --git a/noao/digiphot/apphot/radprof/radprof.key b/noao/digiphot/apphot/radprof/radprof.key new file mode 100644 index 00000000..4204743f --- /dev/null +++ b/noao/digiphot/apphot/radprof/radprof.key @@ -0,0 +1,116 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Store the current parameters +d Plot radial profile of current star +i Interactively set parameters using current star +c Fit center of current star +t Fit sky around cursor +a Average sky values fit around several cursor positions +s Fit sky around the current star +p Fit star using current sky +o Fit star using current sky, output results +f Fit current star +spbar Fit current star, output results +m Move to next star in coordinate list +n Fit next star in coordinate list, output results +l Fit remaining stars in coordinate list, output results +r Rewind the coordinate list +e Print error messages +q Exit task + + + Colon Commands + +:show [data/center/sky/fit] List the parameters +:m [n] Move to next [nth] object in coordinate list +:n [n] Fit next [nth] object in coordinate list, output results + + + Colon Parameter Editing Commands + +# Image and file name parameters + +:image [string] Image name +:coords [string] Coordinate file name +:output [string] Output file name + +# Data dependent parameters + +:scale [vlaue] Image scale (units per pixel) +:fwhmpsf [value] Full-width half-maximum of psf (scale units) +:emission [y/n] Emission features (y), absorption (n) +:sigma [value] Standard deviation of sky (counts) +:datamin [value] Minimum good pixel value (counts) +:datamax [value] Maximum good pixel value (counts) + +# Noise parameters + +:noise [string] Noise model (constant|poisson) +:gain [string] Gain image header keyword +:ccdread [string] Readout noise image header keyword +:epadu [value] Gain (electrons per adu) +:readnoise [value] Readout noise (electrons) + +# Observing parameters + +:exposure [value] Exposure time image header keyword +:airmass [string] Airmass image header keyword +:filter [string] Filter image header keyword +:obstime [string] Time of observation image header keyword +:itime [value] Integration time (time units) +:xairmass [value] Airmass value (number) +:ifilter [string] Filter id string +:otime [string] Time of observation (time units) + +# Centering algorithm parameters + +:calgorithm [string] Centering algorithm +:cbox [value] Width of the centering box (scale units) +:cthreshold [value] Centering intensity threshold (sigma) +:cmaxiter [value] Maximum number of iterations +:maxshift [value] Maximum center shift (scale units) +:minsnratio [value] Minimum S/N ratio for centering +:clean [y/n] Clean subraster before centering +:rclean [value] Cleaning radius (scale units) +:rclip [value] Clipping radius (scale units) +:kclean [value] Clean K-sigma rejection limit (sigma) + +# Sky fitting algorithm parameters + +:salgorithm [string] Sky fitting algorithm +:skyvalue [value] User supplied sky value (counts) +:annulus [value] Inner radius of sky annulus (scale units) +:dannulus [value] Width of sky annulus (scale units) +:khist [value] Sky histogram extent (+/- sigma) +:binsize [value] Resolution of sky histogram (sigma) +:sloclip [value] Low-side clipping factor in percent +:shiclip [value] High-side clipping factor in percent +:smaxiter [value] Maximum number of iterations +:smooth [y/n] Lucy smooth the sky histogram +:snreject [value] Maximum number of rejection cycles +:sloreject [value] Low-side pixel rejection limits (sky sigma) +:shireject [value] High-side pixel rejection limits (sky sigma) +:rgrow [value] Region growing radius (scale units) + +# Photometry parameters + +:apertures [string] List of apertures (scale units) +:zmag [value] Zero point of magnitude scale + +# Profile fitting parameters + +:radius [value] Maximum profile radius (scale units) +:step [value] Step size for computed profile (scale units) +:order [value] Number of spline pieces in fit +:kreject [value] K-sigma rejection for fit (fit sigma) +:nreject [value] Maximum number of rejection cycles + +# Marking and plotting parameters + +:mkcenter [y/n] Mark computed centers on display +:mksky [y/n] Mark the sky annuli on the display +:mkapert [y/n] Mark apertures on the display +:radplot [y/n] Plot the radial profile diff --git a/noao/digiphot/apphot/radprof/rprofshow.x b/noao/digiphot/apphot/radprof/rprofshow.x new file mode 100644 index 00000000..c4264ede --- /dev/null +++ b/noao/digiphot/apphot/radprof/rprofshow.x @@ -0,0 +1,75 @@ +include "../lib/display.h" +include "../lib/radprof.h" + +# AP_RPROFSHOW -- Procedure to display all the radprof parameters. + +procedure ap_rprofshow (ap) + +pointer ap # pointer to the apphot strucuture + +bool itob() +int apstati() + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + call ap_spshow (ap) + call printf ("\n") + call ap_mpshow (ap) + call printf ("\n") + call ap_rppshow (ap) + call printf (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) +end + + +# AP_RPSHOW -- Procedure to display radprof parameters. + +procedure ap_rpshow (ap) + +pointer ap # pointer to apphot structure + +bool itob() +int apstati() + +begin + call ap_nshow (ap) + call ap_rppshow (ap) + call printf (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) +end + + +# AP_RPPSHOW -- Procedure to display radprof parameters. + +procedure ap_rppshow (ap) + +pointer ap # pointer to apphot structure + +int apstati() +real apstatr() + +begin + # Print the radial profile characteristics. + call printf ("Radial Profile Fitting Parameters\n") + call printf (" %s = %g %s %s = %g %s\n") + call pargstr (KY_RPRADIUS) + call pargr (apstatr (ap, RPRADIUS)) + call pargstr (UN_RSCALEUNIT) + call pargstr (KY_RPSTEP) + call pargr (apstatr (ap, RPSTEP)) + call pargstr (UN_RSCALEUNIT) + + call printf (" %s = %d %s = %g %s %s = %d\n") + call pargstr (KY_RPORDER) + call pargi (apstati (ap, RPORDER)) + call pargstr (KY_RPKSIGMA) + call pargr (apstatr (ap, RPKSIGMA)) + call pargstr (UN_RSIGMA) + call pargstr (KY_RPNREJECT) + call pargi (apstati (ap, RPNREJECT)) +end diff --git a/noao/digiphot/apphot/radprof/t_radprof.x b/noao/digiphot/apphot/radprof/t_radprof.x new file mode 100644 index 00000000..fda4c459 --- /dev/null +++ b/noao/digiphot/apphot/radprof/t_radprof.x @@ -0,0 +1,306 @@ +include <fset.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" + +# T_RADPROF -- Procedure to compute the radial profiles of a list of +# objects. + +procedure t_radprof () + +pointer image # pointer to the name of the image +pointer output # pointer to the output file name +pointer coords # pointer to the coordinate file +pointer plotfile # pointer to the metacode plot file +pointer graphics # pointer to the graphics device +pointer display # pointer to the display device +int interactive # interactive mode +int cache # cache the input image +int verify # verify the critical parameters +int update # update the critical parameters +int verbose # verbose mode + +pointer sp, cname, outfname, ap, im, id, gd, mgd, str +int limlist, lclist, lolist, sid, lid, cl, out, pfd, root, stat +int imlist, clist, olist, wcs, memstat, req_size, old_size, buf_size + +pointer immap(), gopen() +int imtlen(), imtgetim(), clplen(), clgfil(), btoi(), fnldir(), strncmp() +int strlen(), ap_radprof(), imtopenp(), clpopnu(), open(), clgwrd() +int ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (plotfile, 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 (cname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Set standard output to flush on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get input and output file names. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + clist = clpopnu ("coords") + lclist = clplen (clist) + olist = clpopnu ("output") + lolist = clplen (olist) + + # Check that the image and coordinate list lengths match. + if (limlist < 1 || (lclist > 1 && lclist != limlist)) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatible image and coordinate list lengths") + } + + # Check that image and output list lengths match. + if (lolist > 1 && lolist != limlist) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatible image and output list lengths") + } + + call clgstr ("icommands.p_filename", Memc[cname], SZ_FNAME) + if (Memc[cname] != EOS) + interactive = NO + #else if (lclist == 0) + #interactive = YES + else + interactive = btoi (clgetb ("interactive")) + cache = btoi (clgetb ("cache")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + verbose = btoi (clgetb ("verbose")) + + # Get the radial profile fitting parameters. + call ap_grppars (ap) + + # Verify the radial profile fitting parameters. + if (verify == YES && interactive == NO) { + call ap_rconfirm (ap, NULL, 1) + if (update == YES) + call ap_rpars (ap) + } + + # 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 apseti (ap, 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 apseti (ap, WCSOUT, wcs) + + # Get the graphics and display devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + + # Open the plot files. + if (interactive == 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 + } + } + } else { + gd = NULL + id = NULL + } + + # Open the plot metacode file. + call clgstr ("plotfile", Memc[plotfile], SZ_FNAME) + 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 + + # Begin looping over the image list. + sid = 1 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the image. + im = immap (Memc[image], READ_ONLY, 0) + + # Set the image display viewport. + call apimkeys (ap, im, Memc[image]) + if ((id != NULL) && (id != gd)) + call ap_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 = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + + # Open the coordinate file, where coords is assumed to be a simple + # text file in which the x and y positions are in columns 1 and 2 + # respectively and all remaining fields are ignored. + + if (lclist <= 0) { + cl = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (clist, Memc[coords], SZ_FNAME) + root = fnldir (Memc[coords], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[coords+root], 7) == 0 || root == + strlen (Memc[coords])) { + call ap_inname (Memc[image], Memc[outfname], "coo", + Memc[outfname], SZ_FNAME) + lclist = limlist + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else if (stat != EOF) { + call strcpy (Memc[coords], Memc[outfname], SZ_FNAME) + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else { + call apstats (ap, CLNAME, Memc[outfname], SZ_FNAME) + call seek (cl, BOF) + } + } + call apsets (ap, CLNAME, Memc[outfname]) + call apfroot (Memc[outfname], Memc[str], SZ_FNAME) + call apsets (ap, CLROOT, Memc[str]) + + # Open output text file, if output is "default", dir$default, + # or a directory specification then the extension "prf" is added + # to the image name and a suitable version number is appended to + # the output name. If the output name is null then no output + # file is created. + + if (lolist == 0) { + out = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (olist, Memc[output], SZ_FNAME) + root = fnldir (Memc[output], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[output+root], 7) == 0 || root == + strlen (Memc[output])) { + call apoutname (Memc[image], Memc[outfname], "prf", + Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + lolist = limlist + } else if (stat != EOF) { + call strcpy (Memc[output], Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + } else + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + } + call apsets (ap, OUTNAME, Memc[outfname]) + + # Fit the radial profiles. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = ap_radprof (ap, im, cl, NULL, mgd, NULL, out, + sid, NO, cache) + else if (cl != NULL) { + lid = 1 + call ap_bradprof (ap, im, cl, id, gd, mgd, out, sid, lid, + verbose) + stat = NO + } else + stat = NO + } else + stat = ap_radprof (ap, im, cl, gd, mgd, id, out, sid, YES, + cache) + + # Cleanup. + call imunmap (im) + if (cl != NULL) { + if (lclist > 1) + call close (cl) + } + if (out != NULL && lolist != 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + sid = 1 + } + + # Uncache memory. + call fixmem (old_size) + + if (stat == YES) + break + + } + + # If only one coordinate file for a list of images close list. + if (cl != NULL && lclist == 1) + call close (cl) + + # If only one output file for a list of images close list. + if (out != NULL && lolist == 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + } + + # Close up plot files. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + if (mgd != NULL) + call gclose (mgd) + if (pfd != NULL) + call close (pfd) + + # Free the radial profile fitting structure. + call ap_rpfree (ap) + + # Close up the lists. + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + + call sfree (sp) +end |