diff options
Diffstat (limited to 'noao/digiphot/apphot/phot')
31 files changed, 4359 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/phot/apbphot.x b/noao/digiphot/apphot/phot/apbphot.x new file mode 100644 index 00000000..ea6ee325 --- /dev/null +++ b/noao/digiphot/apphot/phot/apbphot.x @@ -0,0 +1,115 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +# APBPHOT -- Procedure to compute the magnitudes for a list of objects +# interactively. + +procedure apbphot (ap, im, cl, sd, out, id, ld, gd, mgd, gid, interactive) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +int sd # sky file descriptor +int out # output file descriptor +int id, ld # sequence and list numbers +pointer gd # pointer to stdgraph stream +pointer mgd # pointer to plot metacode file +pointer gid # pointer to image display stream +int interactive # interactive pr batch mode + +int stdin, ild, radius, cier, sier, pier +pointer sp, str +real wx, wy +int fscan(), nscan(), apfitsky(), apfitcenter(), apmag(), strncmp() +int apstati() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call fstats (cl, F_FILENAME, Memc[str], SZ_FNAME) + + # Initialize + ild = ld + radius = 4 * (apstatr (ap, ANNULUS) + apstatr (ap, DANNULUS) + 1.0) * + apstatr (ap, SCALE) + call ap_imbuf (ap, radius, YES) + + # Print query. + if (strncmp ("STDIN", Memc[str], 5) == 0) { + stdin = YES + call printf ("Type object x and y coordinates (^D or ^Z to end): ") + call flush (STDOUT) + } else + stdin = NO + + # Loop over the coordinate file. + while (fscan (cl) != EOF) { + + # Get and store the coordinates. + call gargr (wx) + call gargr (wy) + if (nscan () != 2) { + if (stdin == YES) { + call printf ( + "Type object x and y coordinates (^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: + ; + } + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Center the coordinates, fit the sky and compute magnitudes. + cier = apfitcenter (ap, im, wx, wy) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), sd, gd) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, YCENTER), + apstati (ap, POSITIVE), apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + + # Print the results if interactive mode. + if (interactive == YES) { + call ap_qpmag (ap, cier, sier, pier) + if (gid != NULL) + call apmark (ap, gid, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + } + + # Write the results. + if (id == 1) + call ap_param (ap, out, "phot") + call ap_pmag (ap, out, id, ild, cier, sier, pier) + + # Make plots if mgd is enabled. + call ap_pplot (ap, im, id, mgd, YES) + + # Prepare 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 object x and y coordinates (^Z or ^D to end): ") + call flush (STDOUT) + } + } + + call ap_imbuf (ap, 0, YES) + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apcomags.x b/noao/digiphot/apphot/phot/apcomags.x new file mode 100644 index 00000000..c6d15a16 --- /dev/null +++ b/noao/digiphot/apphot/phot/apcomags.x @@ -0,0 +1,118 @@ +include <mach.h> +include "../lib/noise.h" + +# APCOPMAGS -- Procedure to compute the magnitudes from the aperture sums, +# areas and sky values. + +procedure apcopmags (sums, areas, mags, magerrs, naperts, sky, sigma, nsky, + zmag, noise, padu) + +double sums[ARB] # aperture sums +double areas[ARB] # aperture areas +real mags[ARB] # output magnitudes +real magerrs[ARB] # errors in the magnitudes +int naperts # number of apertures +real sky # sky value +real sigma # sigma of the sky values +int nsky # number of sky pixels +real zmag # sky value, sigma and zero point magnitudes +int noise # noise model +real padu # photons per adu + +int i +real err1, err2, err3, err + +begin + # Compute the magnitudes and errors + do i = 1, naperts { + mags[i] = sums[i] - areas[i] * sky + if (mags[i] <= 0.0) + mags[i] = INDEFR + else { + if (IS_INDEFR(sigma)) + err1 = 0.0 + else + err1 = areas[i] * sigma ** 2 + switch (noise) { + case AP_NCONSTANT: + err2 = 0.0 + case AP_NPOISSON: + err2 = mags[i] / padu + default: + err2 = 0.0 + } + if (nsky <= 0) + err3 = 0.0 + else if (IS_INDEFR(sigma)) + err3 = 0.0 + else + err3 = sigma ** 2 * areas[i] ** 2 / nsky + err = err1 + err2 + err3 + if (err <= 0.0) + magerrs[i] = 0.0 + else { + magerrs[i] = 1.0857 * sqrt (err) / mags[i] + } + mags[i] = zmag - 2.5 * log10 (mags[i]) + } + } +end + + +# APCONMAGS -- Procedure to compute the magnitudes from the aperture sums, +# areas and sky values. + +procedure apconmags (sums, areas, mags, magerrs, naperts, sky, sigma, nsky, + zmag, noise, padu, readnoise) + +double sums[ARB] # aperture sums +double areas[ARB] # aperture areas +real mags[ARB] # output magnitudes +real magerrs[ARB] # errors in the magnitudes +int naperts # number of apertures +real sky # sky value +real readnoise # readout noise in electrons +real sigma # sigma of the sky values +int nsky # number of sky pixels +real zmag # sky value, sigma and zero point magnitudes +int noise # noise model +real padu # photons per adu + +int i +real err1, err2, err3, err + +begin + # Compute the magnitudes and errors + do i = 1, naperts { + mags[i] = areas[i] * sky - sums[i] + if (mags[i] <= 0.0) + mags[i] = INDEFR + else { + if (IS_INDEFR(readnoise)) + err1 = 0.0 + else + err1 = areas[i] * (readnoise / padu) ** 2 + switch (noise) { + case AP_NCONSTANT: + err2 = 0.0 + case AP_NPOISSON: + err2 = abs (sums[i]) / padu + default: + err2 = 0.0 + } + if (nsky <= 0) + err3 = 0.0 + else if (IS_INDEFR(sigma)) + err3 = 0.0 + else + err3 = sigma ** 2 * areas[i] ** 2 / nsky + err = err1 + err2 + err3 + if (err <= 0.0) + magerrs[i] = 0.0 + else { + magerrs[i] = 1.0857 * sqrt (err) / abs (sums[i]) + } + mags[i] = zmag - 2.5 * log10 (mags[i]) + } + } +end diff --git a/noao/digiphot/apphot/phot/apgppars.x b/noao/digiphot/apphot/phot/apgppars.x new file mode 100644 index 00000000..09961678 --- /dev/null +++ b/noao/digiphot/apphot/phot/apgppars.x @@ -0,0 +1,35 @@ +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# AP_GPPARS -- Procedure to fetch the phot task parameters. + +procedure ap_gppars (ap) + +pointer ap # pointer to apphot structure + +bool clgetb() +int btoi() + +begin + # Open the apphot strucuture. + call appinit (ap, AP_CENTROID1D, 2.5, AP_MODE, 10.0, 10.0, 3.0, 1, + AP_PWCONSTANT, 2.0, AP_NPOISSON) + + # Get the data dependent parameters. + call ap_gdapars (ap) + + # Get the centering algorithm parameters. + call ap_gcepars (ap) + + # Get the sky fitting parameters. + call ap_gsapars (ap) + + # Get the photometry parameters. + call ap_gphpars (ap) + + # Get the plotting parameters. + call apseti (ap, RADPLOTS, btoi (clgetb ("radplots"))) +end diff --git a/noao/digiphot/apphot/phot/apgqppars.x b/noao/digiphot/apphot/phot/apgqppars.x new file mode 100644 index 00000000..dea31951 --- /dev/null +++ b/noao/digiphot/apphot/phot/apgqppars.x @@ -0,0 +1,75 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# AP_GQPPARS -- Procedure to fetch the phot task parameters. + +procedure ap_gqppars (ap) + +pointer ap # pointer to apphot structure + +int naperts +pointer mp, aperts, str, apstr +real cbox, annulus, dannulus +bool clgetb() +int ap_getaperts(), btoi() +real clgetr() + +begin + call smark (mp) + call salloc (aperts, MAX_NAPERTS, TY_REAL) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (apstr, SZ_LINE, TY_CHAR) + + # Get the center, sky fitting and photometry apertures. + cbox = clgetr ("cbox") / 2.0 + annulus = clgetr ("annulus") + dannulus = clgetr ("dannulus") + call clgstr ("apertures", Memc[apstr], SZ_LINE) + naperts = ap_getaperts (Memc[apstr], Memr[aperts], MAX_NAPERTS) + + # Open the apphot structure. + if (naperts <= 0.0) + call appinit (ap, AP_CENTROID1D, cbox, AP_CENTROID, annulus, + dannulus, 0.0, 1, AP_PWCONSTANT, 2.5, AP_NPOISSON) + else + call appinit (ap, AP_CENTROID1D, cbox, AP_CENTROID, annulus, + dannulus, Memr[aperts], naperts, AP_PWCONSTANT, 2.5, + AP_NPOISSON) + + # Set remaining parameters. + call apseti (ap, SMOOTH, YES) + + if (naperts > 0) + call apsets (ap, APSTRING, Memc[apstr]) + call apsetr (ap, ZMAG, clgetr ("zmag")) + + call clgstr ("exposure", Memc[str], SZ_FNAME) + call apsets (ap, EXPOSURE, Memc[str]) + + call clgstr ("airmass", Memc[str], SZ_FNAME) + call apsets (ap, AIRMASS, Memc[str]) + call apsetr (ap, XAIRMASS, INDEFR) + + call clgstr ("filter", Memc[str], SZ_FNAME) + call apsets (ap, FILTER, Memc[str]) + call apsets (ap, FILTERID, "INDEF") + + call clgstr ("obstime", Memc[str], SZ_FNAME) + call apsets (ap, OBSTIME, Memc[str]) + call apsets (ap, OTIME, "INDEF") + + call apsetr (ap, EPADU, clgetr ("epadu")) + + # Print the display parameters. + call apseti (ap, MKCENTER, btoi (true)) + call apseti (ap, MKSKY, btoi (true)) + call apseti (ap, MKAPERT, btoi (true)) + call apseti (ap, RADPLOTS, btoi (clgetb ("radplots"))) + + # Close the pset files. + call sfree (mp) +end diff --git a/noao/digiphot/apphot/phot/apmag.x b/noao/digiphot/apphot/phot/apmag.x new file mode 100644 index 00000000..e1490be2 --- /dev/null +++ b/noao/digiphot/apphot/phot/apmag.x @@ -0,0 +1,118 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noisedef.h" +include "../lib/photdef.h" +include "../lib/phot.h" + +# APMAG -- Procedure to compute the magnitudes inside a set of apertures for +# a single of object. + +int procedure apmag (ap, im, wx, wy, positive, skyval, skysig, nsky) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +real wx, wy # object coordinates +int positive # emission or absorption features +real skyval # sky value +real skysig # sky sigma +int nsky # number of sky pixels + +int c1, c2, l1, l2, ier, nap +pointer sp, nse, phot, temp +real datamin, datamax, zmag +int apmagbuf() + +begin + # Initalize. + phot = AP_PPHOT(ap) + nse = AP_NOISE(ap) + AP_PXCUR(phot) = wx + AP_PYCUR(phot) = wy + if (IS_INDEFR(wx) || IS_INDEFR(wy)) { + AP_OPXCUR(phot) = INDEFR + AP_OPYCUR(phot) = INDEFR + } else { + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_OPXCUR(phot), AP_OPYCUR(phot), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_OPXCUR(phot), AP_OPYCUR(phot), 1) + default: + AP_OPXCUR(phot) = wx + AP_OPYCUR(phot) = wy + } + } + + call amovkd (0.0d0, Memd[AP_SUMS(phot)], AP_NAPERTS(phot)) + call amovkd (0.0d0, Memd[AP_AREA(phot)], AP_NAPERTS(phot)) + call amovkr (INDEFR, Memr[AP_MAGS(phot)], AP_NAPERTS(phot)) + call amovkr (INDEFR, Memr[AP_MAGERRS(phot)], AP_NAPERTS(phot)) + + # Make sure the center is defined. + if (IS_INDEFR(wx) || IS_INDEFR(wy)) + return (AP_APERT_NOAPERT) + + # Fetch the aperture pixels. + ier = apmagbuf (ap, im, wx, wy, c1, c2, l1, l2) + if (ier == AP_APERT_NOAPERT) + return (AP_APERT_NOAPERT) + + call smark (sp) + call salloc (temp, AP_NAPERTS(phot), TY_REAL) + + # Do photometry for all the apertures. + call amulkr (Memr[AP_APERTS(phot)], AP_SCALE(ap), Memr[temp], + AP_NAPERTS(phot)] + if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap))) { + call apmeasure (im, wx, wy, c1, c2, l1, l2, Memr[temp], + Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot)) + AP_NMINAP(phot) = AP_NMAXAP(phot) + 1 + } else { + 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) + call apbmeasure (im, wx, wy, c1, c2, l1, l2, datamin, datamax, + Memr[temp], Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], + AP_NMAXAP(phot), AP_NMINAP(phot)) + } + + # Make sure that the sky value has been defined. + if (IS_INDEFR(skyval)) + ier = AP_APERT_NOSKYMODE + else { + + # Check for bad pixels. + if ((ier == AP_OK) && (AP_NMINAP(phot) <= AP_NMAXAP(phot))) + ier = AP_APERT_BADDATA + + nap = min (AP_NMINAP(phot) - 1, AP_NMAXAP(phot)) + + # Compute the magnitudes and errors. + if (positive == YES) + call apcopmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], + Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], + nap, skyval, skysig, nsky, AP_ZMAG(phot), + 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, skyval, skysig, nsky, AP_ZMAG(phot), + AP_NOISEFUNCTION(nse), AP_EPADU(nse), AP_READNOISE(nse)) + + # Compute correction for itime. + zmag = 2.5 * log10 (AP_ITIME(ap)) + call aaddkr (Memr[AP_MAGS(phot)], zmag, Memr[AP_MAGS(phot)], nap) + } + + call sfree (sp) + if (ier != AP_OK) + return (ier) + else + return (AP_OK) +end diff --git a/noao/digiphot/apphot/phot/apmagbuf.x b/noao/digiphot/apphot/phot/apmagbuf.x new file mode 100644 index 00000000..15034561 --- /dev/null +++ b/noao/digiphot/apphot/phot/apmagbuf.x @@ -0,0 +1,89 @@ +include <imhdr.h> +include "../lib/apphotdef.h" +include "../lib/photdef.h" +include "../lib/phot.h" + +# APMAGBUF -- Procedure to determine the mapping of the aperture list +# into the input image. + +int procedure apmagbuf (ap, im, wx, wy, c1, c2, l1, l2) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # center coordinates +int c1, c2 # column limits +int l1, l2 # line limits + +int i +pointer phot +real rbuf +int ap_photpix() + +begin + # Check for 0 radius aperture. + phot = AP_PPHOT(ap) + if (Memr[AP_APERTS(phot)] <= 0.0) + return (AP_APERT_NOAPERT) + + # Compute the maximum aperture size + AP_APIX(phot) = NULL + for (i = AP_NAPERTS(phot); AP_APIX(phot) == NULL && i >= 1; i = i - 1) { + rbuf = 2. * Memr[AP_APERTS(phot)+i-1] * AP_SCALE(ap) + AP_APIX(phot) = ap_photpix (im, wx, wy, rbuf, c1, c2, l1, l2) + AP_AXC(phot) = wx - c1 + 1 + AP_AYC(phot) = wy - l1 + 1 + AP_ANX(phot) = c2 - c1 + 1 + AP_ANY(phot) = l2 - l1 + 1 + AP_NMAXAP(phot) = i + } + + # Return the appropriate error code. + if (AP_APIX(phot) == NULL) { + return (AP_APERT_NOAPERT) + } else if (AP_NMAXAP(phot) < AP_NAPERTS(phot)) { + return (AP_APERT_OUTOFBOUNDS) + } else { + return (AP_OK) + } +end + + +# AP_PHOTPIX -- Procedure to determine the line and column limits of the +# required subraster. + +int procedure ap_photpix (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 + +begin + # Check for 0 radius aperture. + half_papert = papert / 2. + if (half_papert <= 0) + return (0) + 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 < 0.5) || (xc2 > (real (ncols) + 0.5)) || + (xl1 < 0.5) || (xl2 > (real (nlines) + 0.5))) + return (0) + + # 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 + 0.5)) + l1 = max (1.0, min (real (nlines), xl1)) + l2 = min (real (nlines), max (1.0, xl2 + 0.5)) + return ((c2 - c1 + 1) * (l2 - l1 + 1)) +end diff --git a/noao/digiphot/apphot/phot/apmeasure.x b/noao/digiphot/apphot/phot/apmeasure.x new file mode 100644 index 00000000..6d64db5e --- /dev/null +++ b/noao/digiphot/apphot/phot/apmeasure.x @@ -0,0 +1,116 @@ +# APMEASURE -- Procedure to measure the fluxes and effective areas of a set of +# apertures. + +procedure apmeasure (im, wx, wy, c1, c2, l1, l2, aperts, sums, areas, naperts) + +pointer im # pointer to image +real wx, wy # center of subraster +int c1, c2 # column limits +int l1, l2 # line limits +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, nx, yindex +double fctn +pointer buf +real xc, yc, apmaxsq, dy2, r2, r +pointer imgs2r() + +begin + # Initialize. + nx = c2 - c1 + 1 + xc = wx - c1 + 1 + yc = wy - l1 + 1 + apmaxsq = (aperts[naperts] + 0.5) ** 2 + call aclrd (sums, naperts) + call aclrd (areas, naperts) + + # Loop over the pixels. + do j = l1, l2 { + buf = imgs2r (im, c1, c2, j, j) + if (buf == NULL) + return + yindex = j - l1 + 1 + dy2 = (yindex - yc) ** 2 + do i = 1, nx { + r2 = (i - xc) ** 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 * Memr[buf+i-1] + areas[k] = areas[k] + fctn + } + } + } +end + + +# APBMEASURE -- Procedure to measure the fluxes and effective areas of a set of +# apertures. + +procedure apbmeasure (im, wx, wy, c1, c2, l1, l2, datamin, datamax, aperts, + sums, areas, naperts, minapert) + +pointer im # pointer to image +real wx, wy # center of subraster +int c1, c2 # column limits +int l1, l2 # line limits +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 apertures + +int i, j, k, nx, yindex, kindex +double fctn +pointer buf +real xc, yc, apmaxsq, dy2, r2, r, pixval +pointer imgs2r() + +begin + # Initialize. + nx = c2 - c1 + 1 + xc = wx - c1 + 1 + yc = wy - l1 + 1 + apmaxsq = (aperts[naperts] + 0.5) ** 2 + call aclrd (sums, naperts) + call aclrd (areas, naperts) + minapert = naperts + 1 + + # Loop over the pixels. + do j = l1, l2 { + buf = imgs2r (im, c1, c2, j, j) + if (buf == NULL) + return + yindex = j - l1 + 1 + dy2 = (yindex - yc) ** 2 + do i = 1, nx { + r2 = (i - xc) ** 2 + dy2 + if (r2 > apmaxsq) + next + r = sqrt (r2) - 0.5 + pixval = Memr[buf+i-1] + kindex = naperts + 1 + 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/phot/appconfirm.x b/noao/digiphot/apphot/phot/appconfirm.x new file mode 100644 index 00000000..53d24034 --- /dev/null +++ b/noao/digiphot/apphot/phot/appconfirm.x @@ -0,0 +1,98 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# AP_PCONFIRM -- Procedure to confirm the critical phot parameters. + +procedure ap_pconfirm (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 +int apstati() +real apstatr(), ap_vcapert(), ap_vannulus(), ap_vdannulus(), ap_vsigma() +real ap_vfwhmpsf(), ap_vdatamin(), ap_vdatamax + +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) + + # Confirm the sky annulus. + if (apstati (ap, SKYFUNCTION) != AP_CONSTANT && + apstati (ap, SKYFUNCTION) != AP_SKYFILE) { + annulus = ap_vannulus (ap) + dannulus = ap_vdannulus (ap) + } else { + annulus = apstatr (ap, ANNULUS) + dannulus = apstatr (ap, DANNULUS) + } + + # Confirm the aperture radii parameter. + call ap_vaperts (ap, Memc[aperts], SZ_LINE) + + # Confirm the datamin and datamax parameters. + skysigma = ap_vsigma (ap) + 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 aperture") + 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 sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apperrors.x b/noao/digiphot/apphot/phot/apperrors.x new file mode 100644 index 00000000..ff67e3ee --- /dev/null +++ b/noao/digiphot/apphot/phot/apperrors.x @@ -0,0 +1,48 @@ +include "../lib/phot.h" + +# AP_PERRORS -- Print the photometry task errors. + +procedure ap_perrors (ap, cier, sier, pier) + +pointer ap # apphot structure +int cier # centering error code +int sier # sky fitting error code +int pier # photometry error code + +begin + # Print the centering error message. + call ap_cerrors (ap, cier) + + # Print the sky fitting error message. + call ap_serrors (ap, sier) + + # Print the photometry error message. + call ap_merrors (ap, pier) + +end + + +# AP_MERRORS -- Print the photometry errors. + +procedure ap_merrors (ap, pier) + +pointer ap # pointer to the apphot structure (unused) +int pier # photometry error + +begin + # Print the photometry error message. + switch (pier) { + case AP_APERT_NOAPERT: + call printf ("Photometry apertures are outside of the image.\n") + case AP_APERT_OUTOFBOUNDS: + call printf ("Photometry apertures are partially outside the image.\n") + case AP_APERT_NOSKYMODE: + call printf ("The sky value is undefined.\n") + case AP_APERT_NEGMAG: + call printf ("The total flux inside the aperture is negative.\n") + case AP_APERT_BADDATA: + call printf ("Bad data in the aperture(s).\n") + default: + call printf ("") + } +end diff --git a/noao/digiphot/apphot/phot/appfree.x b/noao/digiphot/apphot/phot/appfree.x new file mode 100644 index 00000000..6d588f8a --- /dev/null +++ b/noao/digiphot/apphot/phot/appfree.x @@ -0,0 +1,66 @@ +include "../lib/apphotdef.h" +include "../lib/photdef.h" + +# APFREE -- Free the apphot structure. + +procedure appfree (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_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_PHOTCLS -- Procedure to close up the photometry structure and arrays. + +procedure ap_photcls (ap) + +pointer ap # pointer to apphot structure + +pointer phot + +begin + if (AP_PPHOT(ap) == NULL) + return + phot = AP_PPHOT(ap) + if (AP_APERTS(phot) != NULL) + call mfree (AP_APERTS(phot), TY_REAL) + if (AP_MAGS(phot) != NULL) + call mfree (AP_MAGS(phot), TY_REAL) + if (AP_MAGERRS(phot) != NULL) + call mfree (AP_MAGERRS(phot), TY_REAL) + if (AP_SUMS(phot) != NULL) + call mfree (AP_SUMS(phot), TY_DOUBLE) + if (AP_AREA(phot) != NULL) + call mfree (AP_AREA(phot), TY_DOUBLE) + + #if (AP_APIX(phot) != NULL) + #call mfree (AP_APIX(phot), TY_REAL) + #if (AP_XAPIX(phot) != NULL) + #call mfree (AP_XAPIX(phot), TY_REAL) + #if (AP_YAPIX(phot) != NULL) + #call mfree (AP_YAPIX(phot), TY_REAL) + + call mfree (AP_PPHOT(ap), TY_STRUCT) +end diff --git a/noao/digiphot/apphot/phot/apphot.x b/noao/digiphot/apphot/phot/apphot.x new file mode 100644 index 00000000..1f8b942d --- /dev/null +++ b/noao/digiphot/apphot/phot/apphot.x @@ -0,0 +1,511 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/phot.h" +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +define HELPFILE "apphot$phot/phot.key" + +# APPHOT -- Procedure to compute magnitudes for a list of objects + +int procedure apphot (ap, im, cl, sd, gd, mgd, id, out, stid, interactive, + cache) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +int sd # the sky file descriptor +pointer gd # pointer to graphcis descriptor +pointer mgd # pointer to the metacode file +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 newimage, newskybuf, newsky, newcenterbuf, newcenter, newmagbuf, newmag +int colonkey, prev_num, req_num, ip, cier, sier, pier, oid, req_size +int old_size, buf_size, memstat, wcs, key, ltid, newlist + +real apstatr() +int clgcur(), apfitsky(), aprefitsky(), apfitcenter(), aprefitcenter() +int apmag(), apremag(), apgscur(), ctoi(), apstati(), apgqverify() +int apgtverify(), apnew(), ap_avsky(), ap_memstat(), sizeof() +bool fp_equalr() + +define endswitch_ 99 + +begin + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize cursor command. + key = ' ' + Memc[cmd] = EOS + + # Initialize the fitting parameters. + newimage = NO + newcenterbuf = YES + newcenter = YES + newskybuf = YES + newsky = YES + newmagbuf = YES + newmag = YES + cier = AP_OK + sier = AP_OK + pier = AP_OK + + # Intialize 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 current cursor coordinates. + 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 moved. + if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) { + newcenterbuf = YES + newcenter = YES + newskybuf = YES + newsky = YES + newmagbuf = YES + newmag = YES + } + + # Store previous cursor coordinates. + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + # Loop over the colon keystroke commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("phot", ap, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + # Print out error messages. + case 'e': + if (interactive == YES) + call ap_perrors (ap, cier, sier, pier) + + # Print out the help page(s). + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Plot a centered stellar radial profile. + case 'd': + if (interactive == YES) { + call ap_qrad (ap, im, wx, wy, gd) + newmagbuf = YES; newmag = YES + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + } + + # Rewind the list. + case 'r': + if (cl != NULL) { + call seek (cl, BOFL) + ltid = 0 + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Get, measure the 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 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 next object. + newlist = YES + if (key == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + goto endswitch_ + } + + # Measure next object. + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), sd, gd) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + if (stid == 1) + call ap_param (ap, out, "phot") + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Process the remainder of the list. + case 'l': + if (cl != NULL) { + ltid = ltid + 1 + oid = stid + call apbphot (ap, im, cl, sd, out, stid, ltid, gd, mgd, id, + 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 apphot 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 a phot parameter. + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call apphotcolon (ap, im, cl, out, stid, ltid, + Memc[cmd], newimage, newcenterbuf, newcenter, + newskybuf, newsky, newmagbuf, newmag) + goto endswitch_ + } + + # No coordinate list. + 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 the 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 the next object. + newlist = YES + if (colonkey == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + goto endswitch_ + } + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), sd, gd) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), + apstati (ap, MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + if (stid == 1) + call ap_param (ap, out, "phot") + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Show/set a phot parameter. + default: + call apphotcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newcenterbuf, newcenter, newskybuf, newsky, + newmagbuf, newmag) + } + + # Reestablish the image display viewport if necessary. + 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 current parameters in the pset files. + case 'w': + call ap_ppars (ap) + + # Setup phot parameters interactively. + case 'i': + if (interactive == YES) { + call ap_radsetup (ap, im, wx, wy, gd, out, stid) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + + # Verify the critical PHOT parameters. + case 'v': + call ap_pconfirm (ap, out, stid) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + + # Fit the center around the cursor position. + case 'c': + if (newcenterbuf == YES) + cier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + cier = aprefitcenter (ap, im, cier) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), NO, NO) + 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 cursor position. + case 't': + if (newskybuf == YES || ! fp_equalr (wx, + apstatr (ap, SXCUR)) || ! fp_equalr (wy, apstatr (ap, + SYCUR))) + sier = apfitsky (ap, im, wx, wy, sd, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + if (id != NULL) { + call apmark (ap, id, NO, apstati (ap, MKSKY), NO) + 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, sd, id, gd, interactive) + if (interactive == YES) + call ap_qaspsky (ap, sier) + newskybuf = NO; newsky = NO + + # Fit the sky around the current center position. + 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), sd, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + if (id != NULL) { + call apmark (ap, id, NO, apstati (ap, MKSKY), NO) + 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 star 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) + if (newmagbuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, PXCUR)) || ! fp_equalr (apstatr (ap, + PYCUR), apstatr (ap, YCENTER))) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + else + pier = apremag (ap, im, apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, NO, NO, apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + newcenterbuf = NO; newcenter = NO + newmagbuf = NO; newmag = NO + + if (key == 'o') { + if (stid == 1) + call ap_param (ap, out, "phot") + if (newlist == YES) + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + else + call ap_pmag (ap, out, stid, 0, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + } + + # Compute the center, sky, and 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), sd, gd) + else if (newsky == YES) + sier = aprefitsky (ap, im, gd) + if (newmagbuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, PXCUR)) || ! fp_equalr (apstatr (ap, YCENTER), + apstatr (ap, PYCUR))) + pier = apmag (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), apstati (ap, POSITIVE), + apstatr (ap, SKY_MODE), apstatr (ap, SKY_SIGMA), + apstati (ap, NSKY)) + else + pier = apremag (ap, im, apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + if (key == ' ') { + if (stid == 1) + call ap_param (ap, out, "phot") + if (newlist == YES) + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + else + call ap_pmag (ap, out, stid, 0, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + } + + default: + call printf ("Unknown or ambiguous keystroke command\n") + } + +endswitch_ + # Setup for the next object. + key = ' ' + Memc[cmd] = EOS + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apphotcolon.x b/noao/digiphot/apphot/phot/apphotcolon.x new file mode 100644 index 00000000..825459b2 --- /dev/null +++ b/noao/digiphot/apphot/phot/apphotcolon.x @@ -0,0 +1,209 @@ +include <gset.h> +include "../lib/apphot.h" +include "../lib/fitsky.h" +include "../lib/center.h" +include "../lib/phot.h" +include "../lib/display.h" +include "../lib/noise.h" + +# APPHOTCOLON -- Process phot colon commands. + +procedure apphotcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + newcenterbuf, newcenter, newskybuf, newsky, newmagbuf, newmag) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the 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 center buffer ? new center fit ? +int newskybuf, newsky # new sky buffer ? new sky fit ? +int newmagbuf, newmag # new aperture buffer ? new fit ? + +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 command. + 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, PCMDS) != 0) + call apmagcolon (ap, out, stid, cmdstr, newmagbuf, newmag) + 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, newmagbuf, newmag) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0) + call apnscolon (ap, im, out, stid, cmdstr, newcenterbuf, + newcenter, newskybuf, newsky, newmagbuf, newmag) + else + call ap_himcolon (ap, cmdstr) + + call sfree (sp) +end + + +# APMAGCOLON -- Procedure to display and edit the photometry parameters. + +procedure apmagcolon (ap, out, stid, cmdstr, newmagbuf, newmag) + +pointer ap # pointer to apphot structure +int out # output file descriptor +int stid # output number +char cmdstr[ARB] # command string +int newmagbuf # new aperture buffers +int newmag # compute new magnitudes + +bool bval +int ncmd +pointer sp, cmd +real rval +bool itob() +int btoi(), 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 command + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PCMDS) + switch (ncmd) { + case PCMD_APERTURES: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, APERTS, Memc[cmd], SZ_LINE) + call printf ("%s = %s %s\n") + call pargstr (KY_APERTS) + call pargstr (Memc[cmd]) + call pargstr (UN_PSCALEUNIT) + } else { + call apsets (ap, APERTS, Memc[cmd]) + if (stid > 1) + call ap_sparam (out, KY_APERTS, Memc[cmd], UN_PSCALEUNIT, + "list of aperture radii") + newmag = YES + newmagbuf = YES + } + case PCMD_ZMAG: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_ZMAG) + call pargr (apstatr (ap, ZMAG)) + } else { + call apsetr (ap, ZMAG, rval) + if (stid > 1) + call ap_rparam (out, KY_ZMAG, rval, UN_PZMAG, + "zero point of magnitude scale") + newmag = YES + } + case PCMD_MKAPERT: + call gargb (bval) + if (nscan () == 1) { + call printf ("%s = %b\n") + call pargstr (KY_MKAPERT) + call pargb (itob (apstati (ap, MKAPERT))) + } else { + call apseti (ap, MKAPERT, btoi (bval)) + } + default: + call printf ("Unknown command\7\n") + } + + call sfree (sp) +end + + +# AP_HIMCOLON -- Procedure to process commands which alter the centering, sky +# fitting and photometry buffers. + +procedure ap_himcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr[ARB] # command string + +bool bval +int ncmd +pointer sp, cmd +bool itob() +int strdic(), nscan(), btoi(), apstati() + +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, PSHOWARGS) + switch (ncmd) { + case PCMD_CENTER: + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + case PCMD_SKY: + call printf ("\n") + call ap_spshow (ap) + call printf ("\n") + case PCMD_PHOT: + call printf ("\n") + call ap_mpshow (ap) + call printf ("\n") + case PCMD_DATA: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_pshow (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 ambigous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/appinit.x b/noao/digiphot/apphot/phot/appinit.x new file mode 100644 index 00000000..b9f558c2 --- /dev/null +++ b/noao/digiphot/apphot/phot/appinit.x @@ -0,0 +1,100 @@ +include "../lib/apphotdef.h" +include "../lib/photdef.h" +include "../lib/phot.h" + +# APPINIT - Procedure to initialize apphot structure. + +procedure appinit (ap, cfunction, cbox, sfunction, annulus, dannulus, + aperts, napert, weight, 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] # apertures +int napert # number of apertures +int weight # weight for photometry +real fwhmpsf # FWHM of the PSF +int noise # noise model + +begin + # Set the image parameters. + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set up the apphot package defaults. + call ap_defsetup (ap, fwhmpsf) + + # Set up noise structure. + call ap_noisesetup (ap, noise) + + # Set up centering structure. + call ap_ctrsetup (ap, cfunction, cbox) + + # Set up sky fitting structure. + call ap_skysetup (ap, sfunction, annulus, dannulus) + + # Set up photometry structure. + call ap_photsetup (ap, aperts, napert, weight) + + # Set the display options. + call ap_dispsetup (ap) + + # Unused structures are set to null. + AP_PPSF(ap) = NULL + AP_POLY(ap) = NULL + AP_RPROF(ap) = NULL +end + + +# AP_PHOTSETUP -- Procedure to set up the photometry parameters. + +procedure ap_photsetup (ap, aperts, napert, weight) + +pointer ap # pointer to apphot structure +real aperts[ARB] # array of apertures +int napert # number of apertures +int weight # weighting function for photometry + +pointer phot + +begin + # phot structure + call malloc (AP_PPHOT(ap), LEN_PHOTSTRUCT, TY_STRUCT) + phot = AP_PPHOT(ap) + + # Set the default values forthe photometry parameters. + AP_PXCUR(phot) = INDEFR + AP_PYCUR(phot) = INDEFR + AP_NAPERTS(phot) = napert + AP_ZMAG(phot) = DEF_ZMAG + AP_PWEIGHTS(phot) = weight + AP_APSTRING(phot) = EOS + switch (weight) { + case AP_PWCONSTANT: + call strcpy ("constant", AP_PWSTRING(phot), SZ_FNAME) + case AP_PWCONE: + call strcpy ("cone", AP_PWSTRING(phot), SZ_FNAME) + case AP_PWGAUSS: + call strcpy ("gauss", AP_PWSTRING(phot), SZ_FNAME) + default: + call strcpy ("constant", AP_PWSTRING(phot), SZ_FNAME) + } + + # Initialize buffers. + AP_LENABUF(phot) = 0 + AP_NAPIX(phot) = 0 + AP_APIX(phot) = NULL + AP_XAPIX(phot) = NULL + AP_YAPIX(phot) = NULL + + # Allocate the buffers to hold the answers. + call malloc (AP_APERTS(phot), napert, TY_REAL) + call malloc (AP_MAGS(phot), napert, TY_REAL) + call malloc (AP_MAGERRS(phot), napert, TY_REAL) + call malloc (AP_AREA(phot), napert, TY_DOUBLE) + call malloc (AP_SUMS(phot), napert, TY_DOUBLE) + call amovr (aperts, Memr[AP_APERTS(phot)], napert) + call asrtr (Memr[AP_APERTS(phot)], Memr[AP_APERTS(phot)], napert) +end diff --git a/noao/digiphot/apphot/phot/appmag.x b/noao/digiphot/apphot/phot/appmag.x new file mode 100644 index 00000000..9bc7e43c --- /dev/null +++ b/noao/digiphot/apphot/phot/appmag.x @@ -0,0 +1,120 @@ +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/photdef.h" +include "../lib/phot.h" + +# AP_PMAG -- Procedure to write the results of the phot task to the output +# file. + +procedure ap_pmag (ap, fd, id, lid, cier, sier, pier) + +pointer ap # pointer to apphot structure +int fd # output text file +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 i, naperts +int apstati() +real apstatr() + +begin + # Initialize. + if (fd == NULL) + return + + # Write out the object id parameters. + call ap_wid (ap, fd, apstatr (ap, OXINIT), apstatr(ap, OYINIT), + id, lid, '\\') + + # Write out the centering results. + call ap_wcres (ap, fd, cier, '\\') + + # Write out the sky fitting results. + call ap_wsres (ap, fd, sier, '\\') + + # Write out the photometry results. + naperts = apstati (ap, NAPERTS) + if (naperts == 0) + call ap_wmres (ap, fd, 0, pier, " ") + else { + do i = 1, naperts { + if (naperts == 1) + call ap_wmres (ap, fd, i, pier, " ") + else if (i == naperts) + call ap_wmres (ap, fd, i, pier, "* ") + else + call ap_wmres (ap, fd, i, pier, "*\\") + } + } +end + + +# AP_QPMAG -- Procedure to print a quick summary of the phot output on the +# standard output. + +procedure ap_qpmag (ap, cier, sier, pier) + +pointer ap # pointer to apphot structure +int cier # centering error +int sier # sky fitting error +int pier # photometry error + +int i +pointer sp, imname, phot +real apstatr() + +begin + # Initialize. + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + phot = AP_PPHOT(ap) + + # Print the center and sky value. + call apstats (ap, IMROOT, Memc[imname], SZ_FNAME) + call printf ("%s %8.2f %8.2f %8g ") + call pargstr (Memc[imname]) + call pargr (apstatr (ap, OPXCUR)) + call pargr (apstatr (ap, OPYCUR)) + call pargr (apstatr (ap, SKY_MODE)) + + # Print out the magnitudes and errors. + do i = 1, AP_NAPERTS(phot) { + if (i == AP_NAPERTS(phot)) + call printf ("%7.3f ") + else + call printf ("%7.3f ") + call pargr (Memr[AP_MAGS(phot)+i-1]) + } + + # Print out the error codes. + if (cier != AP_OK || sier != AP_OK || pier != AP_OK) { + call printf ("err\n") + } else { + call printf ("ok\n") + } + + call sfree (sp) +end + + +# AP_MAGHDR -- Procedure to print the banner for the phot task on the +# standard output. + +procedure ap_maghdr (ap, fd) + +pointer ap # pointer to apphot strucuture +int fd # output file descriptor + +begin + if (fd == NULL) + return + call ap_idhdr (ap, fd) + call ap_chdr (ap, fd) + call ap_shdr (ap, fd) + call ap_mhdr (ap, fd) +end diff --git a/noao/digiphot/apphot/phot/apppars.x b/noao/digiphot/apphot/phot/apppars.x new file mode 100644 index 00000000..68e08aeb --- /dev/null +++ b/noao/digiphot/apphot/phot/apppars.x @@ -0,0 +1,26 @@ +include "../lib/display.h" + +# AP_PPARS -- Procedure to write out the phot task parameters. + +procedure ap_ppars (ap) + +pointer ap # pointer to apphot structure + +bool itob() +int apstati() + +begin + # Write the data dependent parameters. + call ap_dapars (ap) + + # Write the centering parameters. + call ap_cepars (ap) + + # Write the sky fitting paameters. + call ap_sapars (ap) + + # Write the photometry parameters. + call ap_phpars (ap) + + call clputb ("radplots", itob (apstati (ap, RADPLOTS))) +end diff --git a/noao/digiphot/apphot/phot/appplot.x b/noao/digiphot/apphot/phot/appplot.x new file mode 100644 index 00000000..5172533d --- /dev/null +++ b/noao/digiphot/apphot/phot/appplot.x @@ -0,0 +1,273 @@ +include <pkg/gtools.h> +include <gset.h> +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# AP_PPLOT -- Procedure to compute radial profile plots for the centering +# routine. + +procedure ap_pplot (ap, im, sid, gd, makeplot) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the iraf image +int sid # id number of the star +pointer gd # graphics stream +int makeplot # make a plot ? + +int apert, nx, ny +pointer buf, sp, str, r, gt +real xcenter, ycenter, xc, yc, rmin, rmax, imin, imax +real u1, u2, v1, v2, x1, x2, y1, y2 +int ap_ctrpix() +pointer ap_gtinit() +real apstatr() + +begin + # Initialize + if (gd == NULL || makeplot == NO) + return + + # Check for defined center and get the pixels. + xcenter = apstatr (ap, XCENTER) + ycenter = apstatr (ap, YCENTER) + if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter)) + return + + # Fetch the pixels. + apert = 2 * int (apstatr (ap, SCALE) * (apstatr (ap, ANNULUS) + + apstatr (ap, DANNULUS))) + 1 + buf = ap_ctrpix (im, xcenter, ycenter, apert, xc, yc, nx, ny) + if (buf == NULL) + return + + # Reactivate the work station. + call greactivate (gd, 0) + + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (r, nx * ny, TY_REAL) + + # Compute the radii and the plot limits. + call ap_ijtor2 (Memr[r], nx, ny, xc, yc) + call alimr (Memr[r], nx * ny, rmin, rmax) + call alimr (Memr[buf], nx * ny, imin, imax) + + # Store the viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Initialize the plot. + call apstats (ap, IMROOT, Memc[str], SZ_LINE) + call sprintf (Memc[str], SZ_LINE, "%s Star %d") + call pargstr (Memc[str]) + call pargi (sid) + gt = ap_gtinit (Memc[str], apstatr (ap, OXINIT), apstatr (ap, OYINIT)) + call gclear (gd) + + # Label and annotate the plot. + call ap_ppset (gd, gt, ap, rmin, rmax, imin, imax) + call ap_ppannotate (gd, ap, rmin, rmax, imin, imax) + + # Plot the coordinates. + call ap_plotrad (gd, gt, Memr[r], Memr[buf], nx * ny, "plus") + + # Store the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call ap_gtfree (gt) + call gdeactivate (gd, 0) + call sfree (sp) +end + + +# AP_PPSET -- Procedure to set up the parameters for the phot radial profile +# plot. + +procedure ap_ppset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # the graphics stream +pointer gt # the gtools pointer +pointer ap # the apphot pointer +real xmin, xmax # the minimum and maximum radial distance +real ymin, ymax # the minimum and maximum of the y axes + +int fd, naperts +pointer sp, str, title, temp +real aspect, scale, vx1, vx2, vy1, vy2 +int stropen(), apstati() +real apstatr(), gstatr() + +begin + call smark (sp) + call salloc (str, 5 * SZ_LINE, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + naperts = apstati (ap, NAPERTS) + call salloc (temp, naperts, TY_REAL) + + fd = stropen (Memc[str], 5 * SZ_LINE, WRITE_ONLY) + call sysid (Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + + # Encode the center parameter 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 fitting parameter string + call fprintf (fd, + "Sky: value=%0.2f sigma=%0.2f skew=%0.2f nsky=%d nrej=%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 apertures and magnitudes. + call ap_arrayr (ap, APERTS, Memr[temp]) + call amulkr (Memr[temp], apstatr (ap, SCALE), Memr[temp], naperts) + call fprintf (fd, "Photometry: maxapert=") + call fprintf (fd, "%0.2f mag=") + call pargr (Memr[temp+naperts-1]) + call ap_arrayr (ap, MAGS, Memr[temp]) + call fprintf (fd, "%0.3f merr=") + call pargr (Memr[temp+naperts-1]) + call ap_arrayr (ap, MAGERRS, Memr[temp]) + call fprintf (fd, "%0.3f\n") + call pargr (Memr[temp+naperts-1]) + + # Encode the title. + call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + + call strclose (fd) + + aspect = gstatr (gd, G_ASPECT) + call gsetr (gd, G_ASPECT, 0.70) + scale = apstatr (ap, SCALE) + + # Draw three axes. + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, ymin, ymax) + call glabax (gd, Memc[str], "", "Intensity") + + # Draw the bottom x axis. + 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)", "") + + # Restore the default draw axis parameters. + call gseti (gd, G_YDRAWAXES, 3) + call gseti (gd, G_XDRAWAXES, 3) + call gsetr (gd, G_ASPECT, aspect) + + # Set the mark type. + call gt_sets (gt, GTTYPE, "mark") + + call sfree (sp) +end + + +# AP_PPANNOTATE -- Procedure to annotate the radial plot in phot. + +procedure ap_ppannotate (gd, ap, xmin, xmax, ymin, ymax) + +pointer gd # graphics stream +pointer ap # apphot structure +real xmin, xmax # minimum and maximum of the x axis +real ymin, ymax # minimum and maximum of the y axis + +int i, naperts +pointer sp, str, temp +real annulus, dannulus, skyval, skysigma +int apstati() +real apstatr () + +begin + naperts = apstati (ap, NAPERTS) + call smark (sp) + call salloc (str, SZ_LINE + 1, TY_CHAR) + call salloc (temp, naperts, TY_REAL) + + # Define some temporary variables + annulus = apstatr (ap, SCALE) * apstatr (ap, ANNULUS) + dannulus = annulus + apstatr (ap, SCALE) * apstatr (ap, DANNULUS) + + # Mark the inner sky annulus. + 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, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the outer sky annulus. + 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, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the sky value. + call gseti (gd, G_PLTYPE, GL_SOLID) + skyval = apstatr (ap, SKY_MODE) + if (skyval >= ymin && skyval <= ymax) { + call gamove (gd, xmin, skyval) + call gadraw (gd, xmax, skyval) + } + + # Mark the upper sky sigma. + call gseti (gd, G_PLTYPE, GL_DASHED) + if (! IS_INDEFR(apstatr (ap, SKY_SIGMA))) + skysigma = skyval + apstatr (ap, SHIREJECT) * apstatr (ap, + SKY_SIGMA) + else + skysigma = INDEFR + if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) { + call gamove (gd, xmin, skysigma) + call gadraw (gd, xmax, skysigma) + } + + # Mark the lower sky sigma + if (! IS_INDEFR(apstatr (ap, SKY_SIGMA))) + skysigma = skyval - apstatr (ap, SLOREJECT) * apstatr (ap, + SKY_SIGMA) + else + skysigma = INDEFR + if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) { + call gamove (gd, xmin, skysigma) + call gadraw (gd, xmax, skysigma) + } + + # Mark the appertures. + call gseti (gd, G_PLTYPE, GL_SOLID) + call ap_arrayr (ap, APERTS, Memr[temp]) + call amulkr (Memr[temp], apstatr (ap, SCALE), Memr[temp], naperts) + do i = 1, naperts { + call gamove (gd, Memr[temp+i-1], ymin) + call gadraw (gd, Memr[temp+i-1], ymax) + call sprintf (Memc[str], SZ_LINE, "apert[%d] = %0.2f") + call pargi (i) + call pargr (Memr[temp+i-1]) + call gtext (gd, Memr[temp+i-1], ymax, Memc[str], + "q=h;u=180;v=t;p=r") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/appshow.x b/noao/digiphot/apphot/phot/appshow.x new file mode 100644 index 00000000..5e52b42d --- /dev/null +++ b/noao/digiphot/apphot/phot/appshow.x @@ -0,0 +1,83 @@ +include "../lib/display.h" +include "../lib/phot.h" + +# AP_PSHOW -- Procedure to print the photometry parameters. + +procedure ap_pshow (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 (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) +end + + +# AP_MSHOW -- Procedure to print the photometry parameters on the standard +# output. + +procedure ap_mshow (ap) + +pointer ap # pointer to apphot structure + +bool itob() +int apstati() +begin + call ap_nshow (ap) + call printf ("\n") + call ap_mpshow (ap) + call printf (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) +end + + +# AP_MPSHOW -- Procedure to print the photometry parameters on the standard +# output. + +procedure ap_mpshow (ap) + +pointer ap # pointer to apphot structure + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + # Write out the image and cursor position. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Write out the phot parameters. + call printf ("Photometry Parameters\n") + call apstats (ap, PWSTRING, Memc[str], SZ_LINE) + call printf (" %s = %s %s\n") + call pargstr (KY_PWSTRING) + call pargstr (Memc[str]) + call pargstr (UN_PMODEL) + call apstats (ap, APERTS, Memc[str], SZ_LINE) + call printf (" %s = %s %s\n") + call pargstr (KY_APERTS) + call pargstr (Memc[str]) + call pargstr (UN_PSCALEUNIT) + call printf (" %s = %g\n") + call pargstr (KY_ZMAG) + call pargr (apstatr (ap, ZMAG)) + call printf (" %s = %b\n") + call pargstr (KY_MKAPERT) + call pargb (itob (apstati (ap, MKAPERT))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apqcolon.x b/noao/digiphot/apphot/phot/apqcolon.x new file mode 100644 index 00000000..c3050190 --- /dev/null +++ b/noao/digiphot/apphot/phot/apqcolon.x @@ -0,0 +1,331 @@ +include <error.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# APQCOLON -- Procedure to display and edit the quick photometry parameters. + +procedure apqcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, newcbuf, + newcenter, newsbuf, newsky, newmagbuf, newmag) + +pointer ap # pointer to apphot structure +pointer im # pointer to the iraf image +int cl # coordinate file descriptor +int out # output file descriptor +int stid # output file sequence number +int ltid # coordinate file sequence number +char cmdstr[ARB] # command string +int newimage # new image ? +int newcbuf # new centering buffers ? +int newcenter # compute new center ? +int newsbuf # new sky fitting buffers ? +int newsky # compute new sky ? +int newmagbuf # new aperture buffers ? +int newmag # compute new magnitudes ? + +bool bval +int ip, ncmd +pointer sp, cmd, str +real rval + +bool streq(), itob() +int btoi(), strdic(), nscan(), apstati(), ctowrd(), open() +pointer immap() +real apstatr() +errchk immap, open + +begin + # Allocate temporary space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, 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, QCMDS) + switch (ncmd) { + case QCMD_SHOW: + call printf ("\n") + call ap_qshow (ap) + call printf ("\n") + + case QCMD_CBOXWIDTH: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_CAPERT) + call pargr (2.0 * apstatr (ap, CAPERT)) + call pargstr ("pixels") + } else { + call apsetr (ap, CAPERT, rval / 2.0) + if (stid > 1) + call ap_rparam (out, KY_CAPERT, rval, UN_CSCALEUNIT, + "width of the centering box") + newcbuf = YES + newcenter = YES + } + + case QCMD_ANNULUS: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_ANNULUS) + call pargr (apstatr (ap, ANNULUS)) + call pargstr ("pixels") + } else { + call apsetr (ap, ANNULUS, rval) + if (stid > 1) + call ap_rparam (out, KY_ANNULUS, rval, UN_SSCALEUNIT, + "inner radius of the sky annulus") + newsbuf = YES + newsky = YES + } + + case QCMD_DANNULUS: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_DANNULUS) + call pargr (apstatr (ap, DANNULUS)) + call pargstr ("pixels") + } else { + call apsetr (ap, DANNULUS, rval) + if (stid > 1) + call ap_rparam (out, KY_DANNULUS, rval, UN_SSCALEUNIT, + "width of the sky annulus") + newsbuf = YES + newsky = YES + } + + case QCMD_APERTURES: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, APERTS, Memc[cmd], SZ_LINE) + call printf ("%s = %s %s\n") + call pargstr (KY_APERTS) + call pargstr (Memc[cmd]) + call pargstr ("pixels") + } else { + call apsets (ap, APERTS, Memc[cmd]) + if (stid > 1) + call ap_sparam (out, KY_APERTS, Memc[cmd], UN_PSCALEUNIT, + "list of aperture radii") + newmag = YES + newmagbuf = YES + } + + case QCMD_ZMAG: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_ZMAG) + call pargr (apstatr (ap, ZMAG)) + } else { + call apsetr (ap, ZMAG, rval) + if (stid > 1) + call ap_rparam (out, KY_ZMAG, rval, UN_PZMAG, + "zero point of magnitude scale") + newmag = YES + } + + case QCMD_EPADU: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_EPADU) + call pargr (apstatr (ap, EPADU)) + } else { + call apsetr (ap, EPADU, rval) + if (stid > 1) + call ap_rparam (out, KY_EPADU, rval, UN_NEPADU, "gain") + newmag = YES + } + + case QCMD_EXPOSURE: + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, EXPOSURE, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_EXPOSURE) + call pargstr (Memc[str]) + } else { + ip = 1 + if (ctowrd (Memc[cmd], ip, Memc[str], SZ_FNAME) <= 0) + Memc[str] = EOS + call apsets (ap, EXPOSURE, Memc[str]) + if (im != NULL) + call ap_itime (im, ap) + if (stid > 1) + call ap_sparam (out, KY_EXPOSURE, Memc[str], + UN_AKEYWORD, "exposure time keyword") + } + + case QCMD_AIRMASS: + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, AIRMASS, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_AIRMASS) + call pargstr (Memc[str]) + } else { + ip = 1 + if (ctowrd (Memc[cmd], ip, Memc[str], SZ_FNAME) <= 0) + Memc[str] = EOS + call apsets (ap, AIRMASS, Memc[str]) + if (im != NULL) + call ap_airmass (im, ap) + if (stid > 1) + call ap_sparam (out, KY_AIRMASS, Memc[str], + UN_AKEYWORD, "airmass keyword") + } + + case QCMD_FILTER: + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, FILTER, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_FILTER) + call pargstr (Memc[str]) + } else { + ip = 1 + if (ctowrd (Memc[cmd], ip, Memc[str], SZ_FNAME) <= 0) + Memc[str] = EOS + call apsets (ap, FILTER, Memc[str]) + if (im != NULL) + call ap_filter (im, ap) + if (stid > 1) + call ap_sparam (out, KY_FILTER, Memc[str], + UN_AKEYWORD, "filter keyword") + } + + case QCMD_OBSTIME: + call gargstr (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, OBSTIME, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_OBSTIME) + call pargstr (Memc[str]) + } else { + ip = 1 + if (ctowrd (Memc[cmd], ip, Memc[str], SZ_FNAME) <= 0) + Memc[str] = EOS + call apsets (ap, OBSTIME, Memc[str]) + if (im != NULL) + call ap_otime (im, ap) + if (stid > 1) + call ap_sparam (out, KY_OBSTIME, Memc[str], + UN_AKEYWORD, "filter keyword") + } + + case QCMD_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)) + } + + case QCMD_IMAGE: + call gargwrd (Memc[cmd], SZ_LINE) + call apstats (ap, IMNAME, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s = %s\n") + call pargstr (KY_IMNAME) + call pargstr (Memc[str]) + } else { + if (im != NULL) { + call imunmap (im) + im = NULL + } + iferr { + im = immap (Memc[cmd], READ_ONLY, 0) + } then { + call erract (EA_WARN) + call printf ("Reopening image %s.\n") + call pargstr (Memc[str]) + im = immap (Memc[str], READ_ONLY, 0) + } else { + call apimkeys (ap, im, Memc[cmd]) + newimage = YES + newcbuf = YES; newcenter = YES + newsbuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + } + + case QCMD_COORDS: + call gargwrd (Memc[cmd], SZ_LINE) + call apstats (ap, CLNAME, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s: %s\n") + call pargstr (KY_CLNAME) + call pargstr (Memc[str]) + } else { + if (cl != NULL) { + call close( cl) + cl = NULL + } + iferr { + cl = open (Memc[cmd], READ_ONLY, TEXT_FILE) + } then { + cl = NULL + call erract (EA_WARN) + call apsets (ap, CLNAME, "") + call apsets (ap, CLROOT, "") + call printf ("Coordinate file is undefined.\n") + } else { + call apsets (ap, CLNAME, Memc[cmd]) + call apfroot (Memc[cmd], Memc[str], SZ_FNAME) + call apsets (ap, CLROOT, Memc[str]) + ltid = 0 + } + } + + case QCMD_OUTPUT: + call gargwrd (Memc[cmd], SZ_LINE) + call apstats (ap, OUTNAME, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s = %s\n") + call pargstr (KY_OUTNAME) + call pargstr (Memc[str]) + } else { + if (out != NULL) { + call close (out) + out = NULL + if (stid <= 1) + call delete (Memc[str]) + } + iferr { + out = open (Memc[cmd], READ_ONLY, TEXT_FILE) + } then { + call erract (EA_WARN) + call printf ("Reopening output file: %s\n") + call pargstr (Memc[str]) + if (Memc[str] != EOS) + out = open (Memc[str], APPEND, TEXT_FILE) + else + out = NULL + } else { + call apsets (ap, OUTNAME, Memc[cmd]) + stid = 1 + } + } + + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apqphot.x b/noao/digiphot/apphot/phot/apqphot.x new file mode 100644 index 00000000..b3c80b73 --- /dev/null +++ b/noao/digiphot/apphot/phot/apqphot.x @@ -0,0 +1,494 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/phot.h" +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +define HELPFILE "apphot$phot/qphot.key" + +# APQPHOT -- Procedure to compute quick magnitudes for a list of objects + +int procedure apqphot (ap, im, cl, gd, mgd, id, out, stid, interactive, cache) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # the coordinate list descriptor +pointer gd # pointer to graphcis descriptor +pointer mgd # pointer to the metacode file +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 newskybuf, newsky, newcenterbuf, newcenter, newmagbuf, newmag +int newimage, newlist, ip, wcs, key, colonkey, cier, sier, pier +int ltid, oid, prev_num, req_num, req_size, old_size, buf_size, memstat + +real apstatr() +int clgcur(), apfitsky(), aprefitsky(), apfitcenter(), aprefitcenter() +int apmag(), apremag(), apstati(), apgqverify(), apnew(), ctoi() +int apgtverify(), apgscur(), ap_avsky(), ap_memstat(), sizeof() +bool fp_equalr() + +define endswitch_ 99 + +begin + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize cursor command. + key = ' ' + Memc[cmd] = EOS + + # Initialize the fitting parameters. + newimage = NO + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + cier = AP_OK; sier = AP_OK; pier = 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 current cursor coordinates. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + + # Store previous cursor coordinates. + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + # Loop over the colon keystroke commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("qphot", ap, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + # Print out error messages. + case 'e': + if (interactive == YES) + call ap_perrors (ap, cier, sier, pier) + + # Print out the help page(s). + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (gd != NULL) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Process apphot 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 a qphot parameter. + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call apqcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newcenterbuf, newcenter, newskybuf, + newsky, newmagbuf, newmag) + goto endswitch_ + } + + # No coordinate list. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Get the next object. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + + # Fetch the 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 the next object. + newlist = YES + if (colonkey == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + goto endswitch_ + } + + # Measure object. + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), NULL, gd) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr(ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), + apstati (ap, MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + if (stid == 1) + call ap_param (ap, out, "qphot") + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + default: + call apqcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newcenterbuf, newcenter, newskybuf, newsky, + newmagbuf, newmag) + } + + # Reestablish the correct viewport. + 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 current parameters in the pset files. + case 'w': + call ap_qppars (ap) + + # Plot a centered stellar radial profile. + case 'd': + if (interactive == YES) { + call ap_qrad (ap, im, wx, wy, gd) + newmagbuf = YES; newmag = YES + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + } + + # Rewind the list. + case 'r': + if (cl != NULL) { + call seek (cl, BOF) + ltid = 0 + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Process the remainder of the list. + case 'l': + if (cl != NULL) { + ltid = ltid + 1 + oid = stid + call apbphot (ap, im, cl, NULL, out, stid, ltid, gd, mgd, + id, 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") + + # Get, measure next object in the 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 next object. + newlist = YES + if (key == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + goto endswitch_ + } + + # Measure next object. + cier = apfitcenter (ap, im, xlist, ylist) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), NULL, gd) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + if (stid == 1) + call ap_param (ap, out, "qphot") + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Setup phot parameters interactively. + case 'i': + if (interactive == YES) { + call ap_qradsetup (ap, im, wx, wy, gd, out, stid) + newmagbuf = YES; newmag = YES + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + } + + # Fit the center around the cursor position. + case 'c': + if (newcenterbuf == YES) + cier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + cier = aprefitcenter (ap, im, cier) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), NO, NO) + 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 cursor position. + 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) + if (id != NULL) { + call apmark (ap, id, NO, apstati (ap, MKSKY), NO) + 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 the current center position. + 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) + if (id != NULL) { + call apmark (ap, id, NO, apstati (ap, MKSKY), NO) + 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 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) + if (newmagbuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, PXCUR)) || ! fp_equalr (apstatr (ap, + PYCUR), apstatr (ap, YCENTER))) + pier = apmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + else + pier = apremag (ap, im, apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + if (id != NULL) { + call apmark (ap, id, NO, NO, apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + newmagbuf = NO; newmag = NO + + if (key == 'o') { + if (stid == 1) + call ap_param (ap, out, "qphot") + if (newlist == YES) + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + else + call ap_pmag (ap, out, stid, 0, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + } + + # Compute the center, sky, and 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 (newmagbuf == YES || ! fp_equalr (apstatr (ap, XCENTER), + apstatr (ap, PXCUR)) || ! fp_equalr (apstatr (ap, YCENTER), + apstatr (ap, PYCUR))) + pier = apmag (ap, im, apstatr (ap, XCENTER), + apstatr (ap, YCENTER), apstati (ap, POSITIVE), + apstatr (ap, SKY_MODE), apstatr (ap, SKY_SIGMA), + apstati (ap, NSKY)) + else + pier = apremag (ap, im, apstati (ap, POSITIVE), apstatr (ap, + SKY_MODE), apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), apstati (ap, + MKSKY), apstati (ap, MKAPERT)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_pplot (ap, im, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qpmag (ap, cier, sier, pier) + + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + if (key == ' ') { + if (stid == 1) + call ap_param (ap, out, "qphot") + if (newlist == YES) + call ap_pmag (ap, out, stid, ltid, cier, sier, pier) + else + call ap_pmag (ap, out, stid, 0, cier, sier, pier) + call ap_pplot (ap, im, stid, mgd, YES) + stid = stid + 1 + } + + default: + call printf ("Unknown or ambiguous keystroke command\n") + } + +endswitch_ + # Setup for the next object. + key = ' ' + Memc[cmd] = EOS + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apqppars.x b/noao/digiphot/apphot/phot/apqppars.x new file mode 100644 index 00000000..cb87383b --- /dev/null +++ b/noao/digiphot/apphot/phot/apqppars.x @@ -0,0 +1,42 @@ +include "../lib/apphot.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" +include "../lib/noise.h" +include "../lib/display.h" + +# AP_QPPARS -- Procedure to write out qthe phot task parameters. + +procedure ap_qppars (ap) + +pointer ap # pointer to apphot structure + +pointer mp, str +bool itob() +int apstati() +real apstatr() + +begin + call smark (mp) + call salloc (str, SZ_LINE, TY_CHAR) + + call clputr ("cbox", 2.0 * apstatr (ap, CAPERT)) + call clputr ("annulus", apstatr (ap, ANNULUS)) + call clputr ("dannulus", apstatr (ap, DANNULUS)) + call apstats (ap, APERTS, Memc[str], SZ_LINE) + call clpstr ("apertures", Memc[str]) + + call apstats (ap, EXPOSURE, Memc[str], SZ_LINE) + call clpstr ("exposure", Memc[str]) + call apstats (ap, AIRMASS, Memc[str], SZ_LINE) + call clpstr ("airmass", Memc[str]) + call apstats (ap, FILTER, Memc[str], SZ_LINE) + call clpstr ("filter", Memc[str]) + call apstats (ap, OBSTIME, Memc[str], SZ_LINE) + call clpstr ("obstime", Memc[str]) + call clputr ("epadu", apstatr (ap, EPADU)) + call clputr ("zmag", apstatr (ap, ZMAG)) + call clputb ("radplots", itob (apstati (ap, RADPLOTS))) + + call sfree (mp) +end diff --git a/noao/digiphot/apphot/phot/apqradsetup.x b/noao/digiphot/apphot/phot/apqradsetup.x new file mode 100644 index 00000000..233cd8a7 --- /dev/null +++ b/noao/digiphot/apphot/phot/apqradsetup.x @@ -0,0 +1,105 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/fitsky.h" +include "../lib/center.h" + +define HELPFILE "apphot$phot/iqphot.key" + +# AP_QRADSETUP -- Procedure to set up phot interactively using a radial profile +# plot of a bright star. + +procedure ap_qradsetup (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, wcs, key +pointer sp, str, cmd +real xcenter, ycenter, rmin, rmax, imin, imax, xc, yc, rval +real u1, u2, v1, v2, x1, x2, y1, y2 + +int apfitcenter(), apfitsky(), ap_wmag(), apstati(), clgcur(), ap_showplot() +real apstatr(), ap_ccapert(), ap_cannulus(), ap_cdannulus() + +begin + # Check for open graphics stream + if (gd == NULL) + return + call greactivate (gd, 0) + + # Store the viewport and window coordinates. + 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 (cmd, SZ_LINE, TY_CHAR) + call salloc (str, 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 'c': + rval = ap_ccapert (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 'r': + call ap_caper (ap, gd, out, stid, Memc[str], rmin, rmax, + imin, imax) + case 'v': + 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) + call ap_caper (ap, gd, out, stid, Memc[str], 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 the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call gdeactivate (gd, 0) + call sfree (sp) + + # Print the answer. + cier = apfitcenter (ap, im, xcenter, ycenter) + if (! IS_INDEFR (apstatr (ap, XCENTER)) && + ! IS_INDEFR (apstatr (ap, YCENTER))) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), NULL, gd) + if (! IS_INDEFR (apstatr (ap, SKY_MODE))) + pier = ap_wmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + call ap_pplot (ap, im, 0, gd, apstati (ap, RADPLOTS)) + call ap_qpmag (ap, cier, sier, pier) +end diff --git a/noao/digiphot/apphot/phot/apqshow.x b/noao/digiphot/apphot/phot/apqshow.x new file mode 100644 index 00000000..d8a4b6db --- /dev/null +++ b/noao/digiphot/apphot/phot/apqshow.x @@ -0,0 +1,86 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/phot.h" + +# AP_QSHOW -- Procedure to display the current data parameters. + +procedure ap_qshow (ap) + +pointer ap # pointer to the apphot strucuture + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Set the object charactersitics. + call printf ("Quick Phot Parameters\n") + call printf (" %s: %s (%.2f,%.2f) %s: %s\n") + call apstats (ap, IMNAME, Memc[str], SZ_FNAME) + call pargstr (KY_IMNAME) + call pargstr (Memc[str]) + call pargr (apstatr (ap, CWX)) + call pargr (apstatr (ap, CWY)) + call apstats (ap, OUTNAME, Memc[str], SZ_FNAME) + call pargstr (KY_OUTNAME) + call pargstr (Memc[str]) + + call printf (" %s: %s\n") + call apstats (ap, CLNAME, Memc[str], SZ_FNAME) + call pargstr (KY_CLNAME) + call pargstr (Memc[str]) + + call printf (" %s = %g %s\n") + call pargstr (KY_CAPERT) + call pargr (2.0* apstatr (ap, CAPERT)) + call pargstr ("pixels") + + call printf (" %s = %g %s %s = %g %s\n") + call pargstr (KY_ANNULUS) + call pargr (apstatr (ap, ANNULUS)) + call pargstr ("pixels") + call pargstr (KY_DANNULUS) + call pargr (apstatr (ap, DANNULUS)) + call pargstr ("pixels") + + call apstats (ap, APERTS, Memc[str], SZ_FNAME) + call printf (" %s = %s %s %s = %g %s\n") + call pargstr (KY_APERTS) + call pargstr (Memc[str]) + call pargstr ("pixels") + call pargstr (KY_ZMAG) + call pargr (apstatr (ap, ZMAG)) + call pargstr (UN_PZMAG) + + call printf (" %s = %g %s\n") + call pargstr (KY_EPADU) + call pargr (apstatr (ap, EPADU)) + call pargstr (UN_NEPADU) + call printf (" %s = %s %s = %s\n") + call apstats (ap, EXPOSURE, Memc[str], SZ_FNAME) + call pargstr (KY_EXPOSURE) + call pargstr (Memc[str]) + call apstats (ap, AIRMASS, Memc[str], SZ_FNAME) + call pargstr (KY_AIRMASS) + call pargstr (Memc[str]) + call printf (" %s = %s %s = %s\n") + call apstats (ap, FILTER, Memc[str], SZ_FNAME) + call pargstr (KY_FILTER) + call pargstr (Memc[str]) + call apstats (ap, OBSTIME, Memc[str], SZ_FNAME) + call pargstr (KY_OBSTIME) + call pargstr (Memc[str]) + + call printf (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/apradsetup.x b/noao/digiphot/apphot/phot/apradsetup.x new file mode 100644 index 00000000..b7f47f41 --- /dev/null +++ b/noao/digiphot/apphot/phot/apradsetup.x @@ -0,0 +1,126 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +define HELPFILE "apphot$phot/iphot.key" + +# AP_RADSETUP -- Procedure to set up phot interactively using a radial profile +# plot of a bright star. + +procedure ap_radsetup (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, key, wcs +pointer sp, cmd, str +real xcenter, ycenter, xc, yc, rmin, rmax, imin, imax +real u1, u2, v1, v2, x1, x2, y1, y2, rval + +int apfitcenter(), apfitsky(), ap_wmag(), apstati(), clgcur(), ap_showplot() +real apstatr(), ap_cfwhmpsf(), ap_ccapert(), ap_cannulus(), ap_cdannulus() +real ap_csigma(), ap_crgrow(), ap_crclean(), ap_crclip() +real ap_cdatamin(), ap_cdatamax() + +begin + # Check for open graphics stream + if (gd == NULL) + return + call greactivate (gd, 0) + call gclear (gd) + + # Store the viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Make the plot. + if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax, + imin, imax) == ERR) { + call gdeactivate (gd, 0) + return + } + + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, 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 '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_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + #rval = ap_ccthresh (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) + call ap_caper (ap, gd, out, stid, Memc[str], 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 the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call gdeactivate (gd, 0) + call sfree (sp) + + # Compute and print the answer. + cier = apfitcenter (ap, im, xcenter, ycenter) + if (! IS_INDEFR (apstatr (ap, XCENTER)) && + ! IS_INDEFR (apstatr (ap, YCENTER))) + sier = apfitsky (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), NULL, gd) + if (! IS_INDEFR (apstatr (ap, SKY_MODE))) + pier = ap_wmag (ap, im, apstatr (ap, XCENTER), apstatr (ap, + YCENTER), apstati (ap, POSITIVE), apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + call ap_pplot (ap, im, 0, gd, apstati (ap, RADPLOTS)) + call ap_qpmag (ap, cier, sier, pier) +end diff --git a/noao/digiphot/apphot/phot/apremag.x b/noao/digiphot/apphot/phot/apremag.x new file mode 100644 index 00000000..9cf705d4 --- /dev/null +++ b/noao/digiphot/apphot/phot/apremag.x @@ -0,0 +1,77 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noisedef.h" +include "../lib/photdef.h" +include "../lib/phot.h" + +# APREMAG -- Procedure to recompute the magnitudes inside a set of apertures +# given that the sums and effective areas have already been computed. + +int procedure apremag (ap, im, positive, skyval, skysig, nsky) + +pointer ap # pointer to the apphot structure +pointer im # the input image descriptor +int positive # emission and absorption features +real skyval # sky value +real skysig # sigma of sky +int nsky # number of sky pixels + +int nap +pointer nse, phot +real zmag + +begin + # Initalize. + phot = AP_PPHOT(ap) + nse = AP_NOISE(ap) + call amovkr (INDEFR, Memr[AP_MAGS(phot)], AP_NAPERTS(phot)) + call amovkr (INDEFR, Memr[AP_MAGERRS(phot)], AP_NAPERTS(phot)) + if (IS_INDEFR(AP_PXCUR(phot)) || IS_INDEFR(AP_PYCUR(phot))) { + AP_OPXCUR(phot) = AP_PXCUR(phot) + AP_OPYCUR(phot) = AP_PYCUR(phot) + } else { + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, AP_PXCUR(phot), AP_PYCUR(phot), + AP_OPXCUR(phot), AP_OPYCUR(phot), 1) + case WCS_TV: + call ap_ltov (im, AP_PXCUR(phot), AP_PYCUR(phot), + AP_OPXCUR(phot), AP_OPYCUR(phot), 1) + default: + AP_OPXCUR(phot) = AP_PXCUR(phot) + AP_OPYCUR(phot) = AP_PYCUR(phot) + } + } + + # Check for errors. + if (IS_INDEFR(AP_PXCUR(phot)) || IS_INDEFR(AP_PYCUR(phot))) + return (AP_APERT_NOAPERT) + if (IS_INDEFR(skyval)) + return (AP_APERT_NOSKYMODE) + + nap = min (AP_NMINAP(phot) - 1, AP_NMAXAP(phot)) + + # Compute the magnitudes and errors. + if (positive == YES) + call apcopmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], + Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], nap, + skyval, skysig, nsky, AP_ZMAG(phot), 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, + skyval, skysig, nsky, AP_ZMAG(phot), AP_NOISEFUNCTION(nse), + AP_EPADU(nse), AP_READNOISE(nse)) + + # Correct for itime. + zmag = 2.5 * log10 (AP_ITIME(ap)) + call aaddkr (Memr[AP_MAGS(phot)], zmag, Memr[AP_MAGS(phot)], nap) + + if (AP_NMAXAP(phot) < AP_NAPERTS(phot)) + return (AP_APERT_OUTOFBOUNDS) + else if (AP_NMINAP(phot) <= AP_NMAXAP(phot)) + return (AP_APERT_BADDATA) + else + return (AP_OK) +end diff --git a/noao/digiphot/apphot/phot/iphot.key b/noao/digiphot/apphot/phot/iphot.key new file mode 100644 index 00000000..679f2d50 --- /dev/null +++ b/noao/digiphot/apphot/phot/iphot.key @@ -0,0 +1,18 @@ + Interactive Phot Setup Menu + + v Mark and verify the critical parameters (f,s,c,a,d,r) + + f Mark and verify the full-width half-maximum of psf + 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 aperture radii diff --git a/noao/digiphot/apphot/phot/iqphot.key b/noao/digiphot/apphot/phot/iqphot.key new file mode 100644 index 00000000..426bd523 --- /dev/null +++ b/noao/digiphot/apphot/phot/iqphot.key @@ -0,0 +1,8 @@ + Interactive Qphot Setup Menu + + v Mark and verify the critical parameters (c,a,d,r) + + c Mark and verify the centering box width + a Mark and verify the inner radius of the sky annulus + d Mark and verify the width of the sky annulus + r Mark and verify the aperture radii diff --git a/noao/digiphot/apphot/phot/mkpkg b/noao/digiphot/apphot/phot/mkpkg new file mode 100644 index 00000000..8ca684f8 --- /dev/null +++ b/noao/digiphot/apphot/phot/mkpkg @@ -0,0 +1,81 @@ +# PHOT,QPHOT tasks + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + # general photometry and PHOT task routines + + apbphot.x <fset.h> ../lib/apphot.h \ + ../lib/center.h ../lib/display.h \ + ../lib/fitsky.h + apcomags.x <mach.h> ../lib/noise.h + apgppars.x ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h ../lib/display.h \ + ../lib/noise.h + apmag.x <mach.h> ../lib/apphotdef.h \ + ../lib/noisedef.h ../lib/photdef.h \ + ../lib/apphot.h ../lib/phot.h + apmagbuf.x <imhdr.h> ../lib/apphotdef.h \ + ../lib/photdef.h ../lib/phot.h + apmeasure.x + apppars.x ../lib/display.h + appconfirm.x ../lib/apphot.h ../lib/noise.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h + apperrors.x ../lib/phot.h + appfree.x ../lib/apphotdef.h ../lib/photdef.h + apphot.x <ctype.h> <gset.h> \ + ../lib/apphot.h ../lib/display.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h <imhdr.h> + apphotcolon.x ../lib/noise.h ../lib/apphot.h \ + ../lib/display.h ../lib/fitsky.h \ + ../lib/center.h ../lib/phot.h \ + <gset.h> + appinit.x ../lib/apphotdef.h ../lib/photdef.h \ + ../lib/phot.h + appplot.x <gset.h> <pkg/gtools.h> \ + ../lib/apphot.h ../lib/noise.h \ + ../lib/center.h ../lib/fitsky.h \ + ../lib/phot.h + appmag.x ../lib/apphot.h ../lib/apphotdef.h \ + ../lib/fitsky.h ../lib/center.h \ + ../lib/photdef.h ../lib/phot.h + appshow.x ../lib/display.h ../lib/phot.h + apradsetup.x ../lib/fitsky.h ../lib/center.h \ + ../lib/display.h ../lib/apphot.h + apremag.x <mach.h> ../lib/apphotdef.h \ + ../lib/noisedef.h ../lib/photdef.h \ + ../lib/apphot.h ../lib/phot.h + t_phot.x <fset.h> <gset.h> \ + <lexnum.h> ../lib/apphot.h \ + ../lib/fitsky.h <imhdr.h> + + # QPHOT TASK SPECIFIC ROUTINES + + apgqppars.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/noise.h" "../lib/center.h" \ + "../lib/fitsky.h" "../lib/phot.h" + apqppars.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/noise.h" "../lib/center.h" \ + "../lib/fitsky.h" "../lib/phot.h" + apqphot.x <ctype.h> <gset.h> \ + "../lib/apphot.h" "../lib/display.h" \ + "../lib/center.h" "../lib/fitsky.h" \ + "../lib/phot.h" <imhdr.h> + apqradsetup.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/center.h" "../lib/fitsky.h" + apqcolon.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/noise.h" "../lib/center.h" \ + "../lib/fitsky.h" "../lib/phot.h" \ + <error.h> + apqshow.x "../lib/apphot.h" "../lib/display.h" \ + "../lib/noise.h" "../lib/center.h" \ + "../lib/fitsky.h" "../lib/phot.h" + t_qphot.x <fset.h> <gset.h> \ + <lexnum.h> ../lib/apphot.h \ + ../lib/fitsky.h <imhdr.h> + ; diff --git a/noao/digiphot/apphot/phot/phot.key b/noao/digiphot/apphot/phot/phot.key new file mode 100644 index 00000000..83e577fd --- /dev/null +++ b/noao/digiphot/apphot/phot/phot.key @@ -0,0 +1,109 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Save the current parameters +d Plot radial profile of current star +i Interactively set parameters using current star +c Fit center for current star +t Fit sky around cursor +a Average sky values fit around several cursor positions +s Fit sky around current centered star +p Do photometry for current star, using current sky +o Do photometry for current star, using current sky, output results +f Do photometry for current star +spbar Do photometry for current star, output results +m Move to next star in coordinate list +n Do photometry for next star in coordinate list, output results +l Do photometry for remaining stars in coordinate list, output results +e Print error messages +r Rewind coordinate list +q Exit task + + +Photometry parameters are listed or set with the following commands. + + Colon commands + +:show [data/center/sky/phot] List the parameters +:m [n] Move to next [nth] star in coordinate list +:n [n] Do photometry for next [nth] star 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 [value] Image scale (units per pixel) +:fwhmpsf [value] Full width half maximum of PSF (scale units) +:emission [y/n] Emission feature (y), absorption (n) +:sigma [value] Standard deviation of sky (counts) +:datamin [value] Minimum good data value (counts) +:datamax [value] Maximum good data 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) + +# Observations parameters + +:exposure [string] 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 (+/- sky sigma) +:binsize [value] Resolution of sky histogram (sky sigma) +:smooth [y/n] Lucy smooth the sky histogram +:sloclip [value] Low-side clipping factor in percent +:shiclip [value] High-side clipping factor in percent +:smaxiter [value] Maximum number of iterations +: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 aperture radii (scale units) +:zmag [value] Zero point of magnitude scale + +# Plotting and marking 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 radial profile of object diff --git a/noao/digiphot/apphot/phot/qphot.key b/noao/digiphot/apphot/phot/qphot.key new file mode 100644 index 00000000..ff4dd5b7 --- /dev/null +++ b/noao/digiphot/apphot/phot/qphot.key @@ -0,0 +1,48 @@ + Interactive Photometry Commands + +? Print help +: Colon commands +w Save 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 for current centered star +p Do photometry for current star, using current sky +o Do photometry for current star, using current sky, output results +f Do photometry for current star +spbar Do photometry for current star, output results +e Print error messages +m Move to next star in coordinate list +n Do photometry for next star in coordinate list, output results +l Do photometry for remaining stars in coordinate list, output results +r Rewind the coordinate list +q Exit task + + + Colon Commands + +:show List the parameters +:m [n] Move to next [nth] star in coordinate list +:n [n] Do photometry for next [nth] star in coordinate list, output results + + Colon Parameter Editing Commands + +:image [string] Image name +:output [string] Output file name +:coords [string] Coords file name + +:cbox [value] Width of the centering box (pixels) +:annulus [value] Inner radius of sky annulus (pixels) +:dannulus [value] Width of sky annulus (pixels) +:apertures [string] List of aperture radii (pixels) +:zmag [value] Zero point of magnitude scale (magnitudes) +:epadu [value] Gain (electrons per adu) + +:exposure [string] 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 + +:radplot [y/n] Plot radial profile of object diff --git a/noao/digiphot/apphot/phot/t_phot.x b/noao/digiphot/apphot/phot/t_phot.x new file mode 100644 index 00000000..9d35c994 --- /dev/null +++ b/noao/digiphot/apphot/phot/t_phot.x @@ -0,0 +1,341 @@ +include <fset.h> +include <gset.h> +include <lexnum.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/fitsky.h" + +# T_PHOT -- Procedure to measure magnitudes inside a set of apertures for a list +# of stars in a list of images. + +procedure t_phot () + +pointer image # pointer name of the image +pointer output # pointer output file name +pointer coords # pointer to name of coords file +pointer skyfile # pointer to name of file with sky values +pointer plotfile # file of plot metacode +pointer graphics # graphics display device +pointer display # display device +int interactive # mode of use +int cache # cache the input image pixels in memory +int verify # verify critical parameters in batch mode +int update # update the critical parameters +int verbose # type messages on the terminal + +pointer sp, cname, outfname, str, ap, im, gd, mgd, id +int limlist, lclist, lolist, lslist, sid, lid, sd, out, cl, root, stat, pfd +int imlist, clist, olist, slist, memstat, old_size, wcs, req_size +int buf_size + +pointer immap(), gopen() +int imtlen(), imtgetim(), clplen(), clgfil(), btoi(), apstati(), strncmp() +int fnldir(), strlen(), apphot(), imtopenp(), clpopnu(), open(), clgwrd() +int ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate temporary 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 (skyfile, 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_LINE, TY_CHAR) + + # Set the standard output to flush on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get the task parameters. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + clist = clpopnu ("coords") + lclist = clplen (clist) + olist = clpopnu ("output") + lolist = clplen (olist) + + # Check that 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")) + verbose = btoi (clgetb ("verbose")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + + # Get the graphics and display devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + + # Open the graphics and display devices. + 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 + + # Intialize the phot structure. + call ap_gppars (ap) + + # Confirm the algorithm parameters. + if (verify == YES && interactive == NO) { + call ap_pconfirm (ap, NULL, 1) + if (update == YES) + call ap_ppars (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 file name for the sky values. + if (apstati (ap, SKYFUNCTION) == AP_SKYFILE) { + slist = clpopnu ("skyfile") + lslist = clplen (slist) + if (limlist < 1 || (lslist > 1 && lslist != limlist)) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call clpcls (slist) + call error (0, "Imcompatible image and sky file list lengths") + } + } else + sd = NULL + + # Begin looping over the image list. + sid = 1 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the image and store image parameters. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (ap, im, Memc[image]) + + # Set the image display viewport. + 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_LINE) + call apsets (ap, CLROOT, Memc[str]) + + # Open the skys file. + if (lslist <= 0) { + sd = NULL + call strcpy ("", Memc[skyfile], SZ_FNAME) + } else if (clgfil (slist, Memc[skyfile], SZ_FNAME) != EOF) + sd = open (Memc[skyfile], READ_ONLY, TEXT_FILE) + else + call seek (sd, BOF) + #call apsets (ap, SKYNAME, Memc[skyfile]) + + # Open the output text file, if output is "default", dir$default or + # a directory specification then the extension "mag" is added to the + # image name and a suitable version number is appended to the output + # name. If output is the null string 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], "mag", + 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]) + + # Do aperture photometry. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = apphot (ap, im, cl, sd, NULL, mgd, NULL, out, sid, + NO, cache) + else if (cl != NULL) { + lid = 1 + call apbphot (ap, im, cl, sd, out, sid, lid, gd, mgd, id, + verbose) + stat = NO + } else + stat = NO + } else + stat = apphot (ap, im, cl, sd, gd, mgd, id, out, sid, YES, + cache) + + # Cleanup. + call imunmap (im) + if (cl != NULL) { + if (clplen(clist) > 1) + call close (cl) + } + if (sd != NULL) { + if (lslist > 1) + call close (sd) + } + 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 + } + + + # Close the coordinate, sky and output files. + if (cl != NULL && lclist == 1) + call close (cl) + if (sd != NULL && lslist == 1) + call close (sd) + 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 the 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 apphot data structures. + call appfree (ap) + + # Close the coord, sky and image lists. + call imtclose (imlist) + call clpcls (clist) + if (sd != NULL) + call clpcls (slist) + call clpcls (olist) + + # Free working space. + call sfree (sp) +end diff --git a/noao/digiphot/apphot/phot/t_qphot.x b/noao/digiphot/apphot/phot/t_qphot.x new file mode 100644 index 00000000..91174e78 --- /dev/null +++ b/noao/digiphot/apphot/phot/t_qphot.x @@ -0,0 +1,293 @@ +include <fset.h> +include <gset.h> +include <lexnum.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/fitsky.h" + +# T_QPHOT -- Procedure to measure magnitudes inside a set of apertures for a +# list of stars in a list of images. + +procedure t_qphot () + +pointer image # pointer name of the image +pointer output # pointer output file name +pointer coords # pointer to the coordinate file +pointer plotfile # file of plot metacode +pointer graphics # graphics display device +pointer display # display device +int cache # cache input image pixels in memory +int verbose # verbose mode + +pointer sp, outfname, cname, ap, im, gd, mgd, id, str +int limlist, lclist, lolist, cl, sid, lid, out, root, stat, pfd, interactive +int imlist, olist, clist, memstat, wcs, req_size, old_size, buf_size + +pointer immap(), gopen() +int imtlen(), imtgetim(), clplen(), clgfil(), btoi(), strncmp() +int fnldir(), strlen(), apqphot(), imtopenp(), clpopnu(), open() +int clgwrd(), ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate temporary space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (output, 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 the standard output to flush on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get the task parameters. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + clist = clpopnu ("coords") + lclist = clplen (clist) + olist = clpopnu ("output") + lolist = clplen (olist) + + # Check that 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 (olist) + call error (0, "Imcompatible image and output list lengths") + } + + # Intialize the phot structure. + 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")) + verbose = btoi (clgetb ("verbose")) + + # Get the parameters. + call ap_gqppars (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 and store image parameters. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (ap, im, Memc[image]) + + # Set the image display viewport. + 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 the output text file, if output is "default", dir$default or + # a directory specification then the extension "mag" is added to the + # image name and a suitable version number is appended to the output + # name. If output is the null string 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], "mag", + 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]) + + # Do aperture photometry. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = apqphot (ap, im, cl, NULL, mgd, NULL, out, sid, NO, + cache) + else if (cl != NULL) { + lid = 1 + call apbphot (ap, im, cl, NULL, out, sid, lid, gd, mgd, id, + verbose) + stat = NO + } else + stat = NO + } else + stat = apqphot (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 + } + + + # Close the coordinate, sky and output files. + if (cl != NULL && 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]) + } + } + + # Close up the 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 apphot data structure. + call appfree (ap) + + # Close the coord, sky and image lists. + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + + # Free working space. + call sfree (sp) +end |