diff options
Diffstat (limited to 'noao/digiphot/apphot/polyphot')
29 files changed, 3947 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/polyphot/apgypars.x b/noao/digiphot/apphot/polyphot/apgypars.x new file mode 100644 index 00000000..960f9cbc --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apgypars.x @@ -0,0 +1,29 @@ +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +# AP_GYPARS -- Procedure to fetch the polyphot task parameters. + +procedure ap_gypars (ap) + +pointer ap # pointer to apphot fitting structure + +begin + # Open the apphot strucuture. + call apyinit (ap, AP_CENTROID1D, 2.5, AP_MODE, 10.0, 10.0, 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_gpopars (ap) +end diff --git a/noao/digiphot/apphot/polyphot/apmkpylist.x b/noao/digiphot/apphot/polyphot/apmkpylist.x new file mode 100644 index 00000000..6b9d1ea0 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apmkpylist.x @@ -0,0 +1,279 @@ +include <ctype.h> +include <gset.h> +include "../lib/apphot.h" +include "../lib/polyphot.h" + +define HELPFILE "apphot$polyphot/polymark.key" + +# AP_MKPYLIST -- Procedure to make polygon and accompanying coordinate list. + +int procedure ap_mkpylist (im, py, pl, cl, id, gd, pid, cid) + +pointer im # pointer to IRAF image +pointer py # pointer to the POLYPHOT structure +int pl # starlist file descriptor +int cl # coordinate file list +pointer id # pointer to image display stream +pointer gd # pointer to graphics display stream +int pid # polygon id sequence number +int cid # coordinate list sequence number + +int key, nvertices, wcs, ptid, ltid, prev_num, req_num +int ip, colonkey, firstpoly, newpoly, delim +pointer sp, cmd, x, y, xshift, yshift +real wx, wy, xmean, ymean + +int clgcur(), apgqverify(), apgtverify(), ap_ymkpoly(), ctoi() +int ap_ynextobj() +real apstatr() +data delim /';'/ + +define endswitch_ 99 + +begin + # Allocate temporary space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (x, MAX_NVERTICES + 1, TY_REAL) + call salloc (y, MAX_NVERTICES + 1, TY_REAL) + call salloc (xshift, MAX_NVERTICES + 1, TY_REAL) + call salloc (yshift, MAX_NVERTICES + 1, TY_REAL) + + # Initialize the cursor read. + key = ' ' + Memc[cmd] = EOS + + # Initialize the sequencing. + firstpoly = YES + newpoly = NO + ptid = 0 + ltid = 0 + xmean = INDEFR + ymean = INDEFR + + # Loop over the polygon file. + nvertices = 0 + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Set the current cursor coordinates. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (py, CWX, wx) + call apsetr (py, CWY, wy) + + # Loop over the colon commands. + switch (key) { + + # Quit. + case 'q': + if (apgqverify ("polymark", NULL, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + + # Plot a centered stellar radial profile. + case 'd': + call ap_qrad (py, im, wx, wy, gd) + + # Print the help page. + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Colon escape commands. + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + colonkey = Memc[cmd+ip-1] + + switch (colonkey) { + case 'm': + + # Decode a polymark colon command. + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call printf ( + "Unknown or ambigous keystroke command\n") + goto endswitch_ + } + + # The polygon file is undefined. + if (pl == NULL) { + call printf ("The polygon list is undefined\n") + goto endswitch_ + } + + # Read the next polygon. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, + Memr[xshift], Memr[yshift], MAX_NVERTICES, prev_num, + req_num, ltid, ptid) + if (nvertices == EOF) { + call printf ( + "End of polygon list, use r key to rewind\n") + goto endswitch_ + } + + # The polygon is undefined. + if (nvertices < 3) { + call printf ("The polygon has fewer than 3 vertices\n") + goto endswitch_ + } + + # Mark the polygon. + if (id != NULL) { + call appymark (py, id, Memr[xshift], Memr[yshift], + nvertices + 1, NO, NO, YES) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + default: + call printf ("Unknown or ambigous keystroke command\n") + } + + # Draw the next polygon in the list. + case 'm': + + # No polygon file. + if (pl == NULL) { + call printf ("The polygon list is undefined\n") + goto endswitch_ + } + + # Read the next polygon. + prev_num = ltid + req_num = ltid + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, + Memr[xshift], Memr[yshift], MAX_NVERTICES, prev_num, + req_num, ltid, ptid) + + # The polygon is undefined. + if (nvertices == EOF) { + call printf ("End of polygon list, use r key to rewind\n") + goto endswitch_ + } + + # The polygon is ill-defined. + if (nvertices < 3) { + call printf ("The polygon has fewer than 3 vertices\n") + goto endswitch_ + } + + # Mark the polygon. + if (id != NULL) { + call appymark (py, id, Memr[xshift], Memr[yshift], + nvertices + 1, NO, NO, YES) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # Rewind the polygon and coordinate lists. + case 'r': + if (pl != NULL) { + call seek (pl, BOF) + if (cl != NULL) + call seek (cl, BOF) + ptid = 0 + ltid = 0 + } else + call printf ("The polygon list is undefined\n") + + # Draw the remaining polygons on the display. + case 'l': + if (pl == NULL) { + call printf ("The polygon list is undefined\n") + } else if (id != NULL) { + call ap_ydraw (py, im, cl, pl, ltid, ptid, id) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + + # Define the polygon interactively. + case 'g': + if (gd == id) + nvertices = ap_ymkpoly (py, im, id, Memr[x], Memr[y], + MAX_NVERTICES, NO) + else + nvertices = ap_ymkpoly (py, im, id, Memr[x], Memr[y], + MAX_NVERTICES, YES) + xmean = apstatr (py, PYXMEAN) + ymean = apstatr (py, PYYMEAN) + if (nvertices == EOF) { + newpoly = NO + call printf ("The polygon is undefined\n") + } else if (nvertices <= 2) { + newpoly = NO + call printf ( + "The polygon has fewer then 3 vertices\n") + } else { + newpoly = YES + if (id != NULL) { + if (gd == id) + call gflush (id) + else + call gframe (id) + } + } + + # Mark the current polygon on the display. + case 'f': + if (id != NULL) { + if (! IS_INDEFR(xmean) && ! IS_INDEFR(ymean)) { + call aaddkr (Memr[x], wx - xmean, Memr[xshift], + nvertices + 1) + call aaddkr (Memr[y], wy - ymean, Memr[yshift], + nvertices + 1) + call appymark (py, id, Memr[xshift], Memr[yshift], + nvertices + 1, NO, NO, YES) + if (gd == id) + call gflush (id) + else + call gframe (id) + } else + call printf ("The polygon is undefined\n") + } + + + # Mark the current polygon on the display and write to file. + case ' ': + if (! IS_INDEFR(xmean) && ! IS_INDEFR(ymean)) { + call ap_ywrite (py, im, cl, pl, Memr[x], Memr[y], nvertices, + cid, pid, firstpoly, newpoly) + if (id != NULL) { + call aaddkr (Memr[x], wx - xmean, Memr[xshift], + nvertices + 1) + call aaddkr (Memr[y], wy - ymean, Memr[yshift], + nvertices + 1) + call appymark (py, id, Memr[xshift], Memr[yshift], + nvertices + 1, NO, NO, YES) + if (gd == id) + call gflush (id) + else + call gframe (id) + } + } else + call printf ("The polygon is undefined\n") + + default: + call printf ("Unknown or ambigous keystroke command\n") + } + +endswitch_ + # Reset the keystroke and command defaults. + call apsetr (py, WX, apstatr (py, CWX)) + call apsetr (py, WY, apstatr (py, CWY)) + key = ' ' + Memc[cmd] = EOS + } +end diff --git a/noao/digiphot/apphot/polyphot/appyerrors.x b/noao/digiphot/apphot/polyphot/appyerrors.x new file mode 100644 index 00000000..96a7bc93 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/appyerrors.x @@ -0,0 +1,35 @@ +include "../lib/polyphot.h" + +# AP_PYERRORS -- Procedure to print out polyphot error messages when the +# task is run in interactive mode. + +procedure ap_pyerrors (ap, cier, sier, pier) + +pointer ap # apphot structure +int cier # centering error code +int sier # sky fitting error code +int pier # polyphot error code + +begin + # Print centering errors. + call ap_cerrors (ap, cier) + + # Print sky fitting errors. + call ap_serrors (ap, sier) + + # Print the polyphot errors. + switch (pier) { + case PY_NOPOLYGON: + call printf ("The polygon is undefined or too few vertices.\n") + case PY_OUTOFBOUNDS: + call printf ("The polygon is partially outside the image.\n") + case PY_NOPIX: + call printf ("The effective polygon area is 0.0.\n") + case PY_NOSKYMODE: + call printf ("The sky value is undefined.\n") + case PY_BADDATA: + call printf ("Bad pixels inside the polygon.\n") + default: + call printf ("") + } +end diff --git a/noao/digiphot/apphot/polyphot/apybphot.x b/noao/digiphot/apphot/polyphot/apybphot.x new file mode 100644 index 00000000..8a4783e8 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apybphot.x @@ -0,0 +1,94 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +# AP_YBPHOT -- Measure the flux inside a list of polygons. + +procedure ap_ybphot (py, im, cl, pl, out, id, ld, pd, gid, interactive) + +pointer py # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # coordinates file descriptor +int pl # vertices list file descriptor +int out # output file descriptor +int id # output file sequence number +int ld # coordinate list number +int pd # polygon list number +pointer gid # pointer to image display stream +int interactive # interactive or batch mode + +real apstatr() +pointer sp, x, y, xout, yout +int req_num, prev_num, cier, sier, pier, nvertices, delim +int ap_ynextobj(), ap_ycenter(), apfitsky(), ap_yfit(), apstati() +data delim /';'/ + +begin + # Allocate temporary space for arrays. + call smark (sp) + call salloc (x, MAX_NVERTICES + 1, TY_REAL) + call salloc (y, MAX_NVERTICES + 1, TY_REAL) + call salloc (xout, MAX_NVERTICES + 1, TY_REAL) + call salloc (yout, MAX_NVERTICES + 1, TY_REAL) + + # Initialize + if (pl != NULL) + call seek (pl, BOF) + if (cl != NULL) + call seek (cl, BOF) + + # Get the first polygon. + pd = 0 + prev_num = 0 + req_num = ld + 1 + nvertices = ap_ynextobj (py, im, gid, pl, cl, delim, Memr[x], Memr[y], + MAX_NVERTICES, prev_num, req_num, ld, pd) + + while (nvertices != EOF) { + + # Fit the center, sky and measure the polygon. + cier = ap_ycenter (py, im, apstatr (py, PYCX), apstatr (py, PYCY), + Memr[x], Memr[y], nvertices + 1) + sier = apfitsky (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), NULL, NULL) + pier = ap_yfit (py, im, Memr[x], Memr[y], + nvertices + 1, apstatr (py, SKY_MODE), apstatr (py, + SKY_SIGMA), apstati (py, NSKY)) + + # Write the output to the standard output. + if (interactive == YES) { + call ap_qyprint (py, cier, sier, pier) + if (gid != NULL) + call appymark (py, gid, Memr[x], Memr[y], nvertices + 1, + YES, apstati (py, MKSKY), apstati (py, MKPOLYGON)) + } + + # Write the output to a file. + if (id == 1) + call ap_param (py, out, "polyphot") + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, Memr[x], Memr[y], Memr[xout], Memr[yout], + nvertices + 1) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[xout], Memr[yout], + nvertices + 1) + default: + call amovr (Memr[x], Memr[xout], nvertices + 1) + call amovr (Memr[y], Memr[yout], nvertices + 1) + } + call ap_yprint (py, out, Memr[xout], Memr[yout], nvertices, id, ld, + pd, cier, sier, pier) + id = id + 1 + + # Setup for next polygon. + prev_num = ld + req_num = ld + 1 + nvertices = ap_ynextobj (py, im, gid, pl, cl, delim, Memr[x], + Memr[y], MAX_NVERTICES, prev_num, req_num, ld, pd) + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/polyphot/apycenter.x b/noao/digiphot/apphot/polyphot/apycenter.x new file mode 100644 index 00000000..faa2d884 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apycenter.x @@ -0,0 +1,59 @@ +include "../lib/center.h" +include "../lib/polyphot.h" + +# AP_YCENTER -- Center the polygon using the centering package routines. + +int procedure ap_ycenter (py, im, wx, wy, x, y, nver) + +pointer py # pointer to polyphot structure +pointer im # pointer to the IRAF image +real wx # initial x center +real wy # initial y center +real x[ARB] # x coordinates of the polygon +real y[ARB] # y coordinates of the polygon +int nver # number of vertices + +int cier +int apfitcenter() +real apstatr() + +begin + if (IS_INDEFR(apstatr (py, PYCX)) || IS_INDEFR(apstatr (py, PYCY))) + cier = apfitcenter (py, im, INDEFR, INDEFR) + else + cier = apfitcenter (py, im, wx, wy) + + if (! IS_INDEFR (apstatr (py, XCENTER)) && ! IS_INDEFR (apstatr (py, + YCENTER)) && ! IS_INDEFR (apstatr (py, PYCX)) && + ! IS_INDEFR (apstatr (py, PYCY))) + call ap_yshift (py, im, x, y, nver, apstatr (py, XCENTER), + apstatr (py, YCENTER)) + + return (cier) +end + + +# AP_YRECENTER -- Recenter the polygon using the existing buffer of pixels +# and the centering package routines. + +int procedure ap_yrecenter (py, im, x, y, nver, cier) + +pointer py # pointer to polyphot structure +pointer im # the input image descriptor +real x[ARB] # x coordinates of the polygon +real y[ARB] # y coordinates of the polygon +int nver # number of vertices +int cier # original centering error + +int aprefitcenter() +real apstatr() + +begin + cier = aprefitcenter (py, im, cier) + if (! IS_INDEFR (apstatr (py, XCENTER)) && ! IS_INDEFR (apstatr (py, + YCENTER)) && ! IS_INDEFR (apstatr (py, PYCX)) && ! + IS_INDEFR (apstatr (py, PYCY))) + call ap_yshift (py, im, x, y, nver, apstatr (py, XCENTER), + apstatr (py, YCENTER)) + return (cier) +end diff --git a/noao/digiphot/apphot/polyphot/apycolon.x b/noao/digiphot/apphot/polyphot/apycolon.x new file mode 100644 index 00000000..9b70199b --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apycolon.x @@ -0,0 +1,226 @@ +include <error.h> +include <gset.h> +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/display.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +# AP_YCOLON -- Process polyphot task colon commands + +procedure ap_ycolon (ap, im, pl, cl, out, stid, ptid, ltid, cmdstr, newimage, + newcenterbuf, newcenter, newskybuf, newsky, newmagbuf, newmag) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the iraf image +int pl # polygon file descriptor +int cl # coord file descriptor +int out # output file descriptor +int stid # output file number +int ptid # polygon file sequence number +int ltid # coord file 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 + # Fetch the command. + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + # Set the command to the appropriate routine. + 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, 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 ap_nscolon (ap, im, out, stid, cmdstr, newcenterbuf, + newcenter, newskybuf, newsky, newmagbuf, newmag) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, PYCMDS) != 0) + call ap_yycolon (ap, pl, out, stid, ptid, ltid, cmdstr, newmagbuf, + newmag) + else + call ap_yimcolon (ap, cmdstr) + + call sfree (sp) +end + + +# AP_YIMCOLON -- Procedure to process remaining polyphot commands which + +procedure ap_yimcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr[ARB] # command string + +int ncmd +pointer sp, cmd +int strdic() + +begin + # Fetch 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, MISC1) + switch (ncmd) { + case 1: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PYSHOWARGS) + switch (ncmd) { + case 1: + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + case 2: + call printf ("\n") + call ap_spshow (ap) + call printf ("\n") + case 3: + call printf ("\n") + call ap_ypshow (ap) + call printf ("\n") + case 4: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_yshow (ap) + call printf ("\n") + } + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end + + +# AP_YYCOLON -- Procedure to process polyphot commands + +procedure ap_yycolon (ap, pl, out, stid, ptid, ltid, cmdstr, newmagbuf, newmag) + +pointer ap # pointer to apphot structure +int pl # polygon file descriptor +int out # output file descriptor +int stid # output file sequence number +int ptid # polygon file sequence number +int ltid # coords file sequence number +char cmdstr[ARB] # command string +int newmagbuf # new aperture buffers ? +int newmag # compute new magnitudes ? + +bool bval +int ncmd +pointer sp, cmd, str +real rval +string cmds PYCMDS + +bool itob(), streq() +int apstati(), btoi(), strdic(), nscan(), open() +real apstatr() +errchk open, close + +begin + # Allocate the space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the colon command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, cmds) + switch (ncmd) { + + case PLCMD_ZMAG: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_PYZMAG) + call pargr (apstatr (ap, PYZMAG)) + } else { + call apsetr (ap, PYZMAG, rval) + if (stid > 1) + call ap_rparam (out, KY_PYZMAG, rval, UN_PYZMAG, + "zero point of magnitude scale") + newmag = YES + } + + case PLCMD_MKPOLYGON: + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_MKPOLYGON) + call pargb (itob (apstati (ap, MKPOLYGON))) + } else { + call apseti (ap, MKPOLYGON, btoi (bval)) + } + + case PLCMD_POLYGONS: + call gargwrd (Memc[cmd], SZ_LINE) + call apstats (ap, PYNAME, Memc[str], SZ_FNAME) + if (Memc[cmd] == EOS || streq (Memc[cmd], Memc[str])) { + call printf ("%s = %s\n") + call pargstr (KY_PYNAME) + call pargstr (Memc[str]) + } else { + if (pl != NULL) { + call close (pl) + pl = NULL + } + iferr { + pl = open (Memc[cmd], READ_ONLY, TEXT_FILE) + } then { + pl = NULL + call erract (EA_WARN) + call apsets (ap, PYNAME, "") + call apsets (ap, PYROOT, "") + ptid = 0 + ltid = 0 + } else { + call apsets (ap, PYNAME, Memc[cmd]) + call apfroot (Memc[cmd], Memc[str], SZ_FNAME) + call apsets (ap, PYROOT, Memc[str]) + ptid = 0 + ltid = 0 + } + } + + default: + call printf ("Unknown or ambiguous colon command\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/polyphot/apyconfirm.x b/noao/digiphot/apphot/polyphot/apyconfirm.x new file mode 100644 index 00000000..fb10ad32 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyconfirm.x @@ -0,0 +1,99 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" +include "../lib/fitsky.h" + +# AP_YCONFIRM -- Procedure to confirm the critical polyphot parameters. + +procedure ap_yconfirm (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 +real fwhmpsf, capert, annulus, dannulus, skysigma +real datamin, datamax +int apstati() +real apstatr(), ap_vfwhmpsf(), ap_vcapert(), ap_vsigma() +real ap_vannulus(), ap_vdannulus(), ap_vdatamin(), ap_vdatamax() + +begin + call smark (sp) + call salloc (cstr, SZ_FNAME, TY_CHAR) + call salloc (sstr, SZ_FNAME, TY_CHAR) + + call printf ("\n") + + # Confirm the centering algorithm. + call ap_vcstring (ap, Memc[cstr], SZ_FNAME) + + if (apstati (ap, CENTERFUNCTION) != AP_NONE) { + + # Confirm the fwhmpsf. + if (apstati (ap, CENTERFUNCTION) != AP_CENTROID1D) + fwhmpsf = ap_vfwhmpsf (ap) + else + fwhmpsf = apstatr (ap, FWHMPSF) + + # Confirm the centering box. + capert = 2.0 * ap_vcapert (ap) + + } else { + fwhmpsf = apstatr (ap, FWHMPSF) + capert = 2.0 * apstatr (ap, CAPERT) + } + + # Confirm the sky fitting algorithm. + call ap_vsstring (ap, Memc[sstr], SZ_FNAME) + + if (apstati (ap, SKYFUNCTION) != AP_CONSTANT && + apstati (ap, SKYFUNCTION) != AP_SKYFILE) { + + # Confirm the sky annulus parameter. + annulus = ap_vannulus (ap) + + # Confirm the width of the sky annulus. + dannulus = ap_vdannulus (ap) + + } else { + annulus = apstatr (ap, ANNULUS) + dannulus = apstatr (ap, DANNULUS) + } + + # Confirm the sky sigma parameter. + if (apstati (ap, SKYFUNCTION) != AP_SKYFILE) + skysigma = ap_vsigma (ap) + else + skysigma = apstatr (ap, SKYSIGMA) + + # Confirm the datamin and datamax parameters. + datamin = ap_vdatamin (ap) + datamax = ap_vdatamax (ap) + + call printf ("\n") + + # Update the database file. + if (out != NULL && stid > 1) { + call ap_sparam (out, KY_CSTRING, Memc[cstr], UN_CALGORITHM, + "centering algorithm") + call ap_rparam (out, KY_FWHMPSF, fwhmpsf, UN_ASCALEUNIT, + "full width half maximum of the psf") + call ap_rparam (out, KY_CAPERT, capert, UN_CSCALEUNIT, + "width of the centering box") + 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_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/polyphot/apycoords.x b/noao/digiphot/apphot/polyphot/apycoords.x new file mode 100644 index 00000000..0298b0f8 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apycoords.x @@ -0,0 +1,43 @@ +include "../lib/polyphot.h" + +# AP_YCOORDS -- Procedure to fetch the next center coordinates for the polygon. + +int procedure ap_ycoords (cl, delim, xshift, yshift, stdin) + +int cl # coordinates file descriptor +int delim # delimiter character +real xshift # new x coordinate +real yshift # new y coordinate +int stdin # is cl STDIN? + +char marker +int stat +int fscan(), nscan() + +begin + if (stdin == YES) { + call printf ("Type object x and y coordinates (^D or^Z to end): ") + call flush (STDOUT) + } + + stat = fscan (cl) + while (stat != EOF) { + + call gargr (xshift) + call gargr (yshift) + + if (nscan() != 2) { + call reset_scan() + call gargc (marker) + if (int (marker) == delim) + return (NEXT_POLYGON) + #else + #return (NEXT_OBJECT) + } else + return (THIS_OBJECT) + + stat = fscan (cl) + } + + return (stat) +end diff --git a/noao/digiphot/apphot/polyphot/apydraw.x b/noao/digiphot/apphot/polyphot/apydraw.x new file mode 100644 index 00000000..aadcbed0 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apydraw.x @@ -0,0 +1,55 @@ +include "../lib/polyphot.h" + +# APYDRAW -- Procedure to draw polygons interactively on the display. + +procedure ap_ydraw (py, im, cl, pl, ld, pd, id) + +pointer py # pointer to apphot structure +pointer im # theinput image descriptor +int cl # coordinates file descriptor +int pl # vertices list file descriptor +int ld # coordinate list number +int pd # polygon list number +pointer id # pointer to image display stream + +int req_num, prev_num, nvertices, delim +pointer sp, x, y +int ap_ynextobj() +data delim /';'/ + +begin + # Allocate temporary space for arrays. + call smark (sp) + call salloc (x, MAX_NVERTICES + 1, TY_REAL) + call salloc (y, MAX_NVERTICES + 1, TY_REAL) + + # Rewind the polygon and coordinate files. + if (pl != NULL) + call seek (pl, BOF) + if (cl != NULL) + call seek (cl, BOF) + + # Initialize. + pd = 0 + prev_num = 0 + req_num = ld + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, Memr[x], Memr[y], + MAX_NVERTICES, prev_num, req_num, ld, pd) + + # Loop over the coordinate and polygon file. + while (nvertices != EOF) { + + # Draw the polygon if it has more than three vertices. + if (id != NULL && nvertices >= 3) + call appymark (py, id, Memr[x], Memr[y], nvertices + 1, + NO, NO, YES) + + # Setup for next polygon. + prev_num = ld + req_num = ld + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, Memr[x], + Memr[y], MAX_NVERTICES, prev_num, req_num, ld, pd) + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/polyphot/apyfit.x b/noao/digiphot/apphot/polyphot/apyfit.x new file mode 100644 index 00000000..abb76a6d --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyfit.x @@ -0,0 +1,578 @@ +include <imhdr.h> +include <mach.h> +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/polyphot.h" + +# AP_YFIT -- Procedure to compute the magnitude of an object inside a polygonal +# aperture. + +int procedure ap_yfit (py, im, xver, yver, nver, skyval, skysig, nsky) + +pointer py # pointer to polyphot strucuture +pointer im # pointer to IRAF image +real xver[ARB] # x vertices coords +real yver[ARB] # y vertices coords +int nver # number of vertices +real skyval # sky value +real skysig # sigma of sky pixels +int nsky # number of sky pixels + +double flux, area +int noise, badpix, ier +real datamin, datamax, mag, magerr, zmag, padu, itime, readnoise +int apstati(), ap_yyfit(), ap_byyfit() +real apstatr() + +begin + # Initialize. + call apsetd (py, PYFLUX, 0.0d0) + call apsetd (py, PYNPIX, 0.0d0) + call apsetr (py, PYMAG, INDEFR) + call apsetr (py, PYMAGERR, INDEFR) + + # Compute the flux inside the polygon. + if (IS_INDEFR(apstatr (py, DATAMIN)) && IS_INDEFR(apstatr(py, + DATAMAX))) { + ier = ap_yyfit (im, xver, yver, nver, flux, area) + badpix = NO + } else { + if (IS_INDEFR(apstatr (py, DATAMIN))) + datamin = -MAX_REAL + else + datamin = apstatr (py, DATAMIN) + if (IS_INDEFR(apstatr (py, DATAMAX))) + datamax = MAX_REAL + else + datamax = apstatr (py, DATAMAX) + ier = ap_byyfit (im, xver, yver, nver, datamin, datamax, flux, + area, badpix) + } + + if (ier == PY_NOPOLYGON) + return (PY_NOPOLYGON) + else if (ier == PY_NOPIX) + return (PY_NOPIX) + + # Store the results. + call apsetd (py, PYFLUX, flux) + call apsetd (py, PYNPIX, area) + call apseti (py, PYBADPIX, badpix) + + if (IS_INDEFR(skyval)) + return (PY_NOSKYMODE) + + # Get the photometry parameters. + zmag = apstatr (py, PYZMAG) + itime = apstatr (py, ITIME) + noise = apstati (py, NOISEFUNCTION) + padu = apstatr (py, EPADU) + readnoise = apstatr (py, READNOISE) + + # Compute the magnitude and error. + if (badpix == NO) { + if (apstati (py, POSITIVE) == YES) + call apcopmags (flux, area, mag, magerr, 1, skyval, + skysig, nsky, zmag, noise, padu) + else + call apconmags (flux, area, mag, magerr, 1, skyval, + skysig, nsky, zmag, noise, padu, readnoise) + mag = mag + 2.5 * log10 (itime) + + call apsetr (py, PYMAG, mag) + call apsetr (py, PYMAGERR, magerr) + } + + return (ier) +end + + +# AP_YYFIT -- Measure the total flux inside a polygon. + +int procedure ap_yyfit (im, xver, yver, nver, flux, area) + +pointer im # pointer to IRAF image +real xver[ARB] # x coordinates of the vertices +real yver[ARB] # y coordinates of the vertices: +int nver # number of vertices +double flux # flux interior to the polygon +double area # approximate area of polygon + +double fluxx, areax, fctnx, fctny +real xmin, xmax, ymin, ymax, x1, x2, lx, ld +pointer sp, work1, work2, xintr, buf +int i, j, k, linemin, linemax, colmin, colmax, nintr, ier +int ap_yclip() +pointer imgl2r() + +begin + # Check that polygon has at least 3 vertices plus the closing vertex. + if (nver < 4) { + flux = INDEFD + area = 0.0d0 + return (PY_NOPOLYGON) + } + + # Allocate working space. + call smark (sp) + call salloc (work1, nver, TY_REAL) + call salloc (work2, nver, TY_REAL) + call salloc (xintr, nver, TY_REAL) + + # Find the minimum and maximum x and y values of the polygon + # and detemine whether the polygon is partially off the image. + + call alimr (xver, nver, xmin, xmax) + call alimr (yver, nver, ymin, ymax) + if (xmin < 0.5 || xmax > (IM_LEN(im,1) + 0.5) || ymin < 0.5 || ymax > + (IM_LEN(im,2) + 0.5)) + ier = PY_OUTOFBOUNDS + else + ier = PY_OK + + # Find the minimum and maximum image line numbers. + ymin = max (0.5, min (real (IM_LEN (im,2) + 0.5), ymin)) + ymax = min (real (IM_LEN(im,2) + 0.5), max (0.5, ymax)) + linemin = min (int (ymin + 0.5), int (IM_LEN(im,2))) + linemax = min (int (ymax + 0.5), int (IM_LEN(im,2))) + + # Set up the line segment parameters and initialize fluxes and areas. + x1 = 0.5 + x2 = IM_LEN(im,1) + 0.5 + lx = x2 - x1 + flux = 0.0d0 + area = 0.0d0 + + # Loop over the range of lines of interest. + do i = linemin, linemax { + + # Read in image line. + buf = imgl2r (im, i) + if (buf == EOF) + next + + # Find all the x intersection points of image line and polygon. + if (ymin > i) + ld = min (i + 1, linemax) + else if (ymax < i) + ld = max (i - 1, linemin) + else + ld = i + nintr = ap_yclip (xver, yver, Memr[work1], Memr[work2], + Memr[xintr], nver, lx, ld) + if (nintr <= 0) + next + fctny = min (i + 0.5, ymax) - max (i - 0.5, ymin) + + # Sort the x intersection points + call asrtr (Memr[xintr], Memr[xintr], nintr) + + # Integrate the flux in each line segment. + fluxx = 0.0d0 + areax = 0.0d0 + do j = 1, nintr, 2 { + + # Compute the line segment limits. + xmin = min (real (IM_LEN(im,1) + 0.5), max (0.5, + Memr[xintr+j-1])) + xmax = min (real (IM_LEN(im,1) + 0.5), max (0.5, + Memr[xintr+j])) + colmin = min (int (xmin + 0.5), int (IM_LEN(im,1))) + colmax = min (int (xmax + 0.5), int (IM_LEN(im,1))) + + # Sum the contribution from a particular line segment. + do k = colmin, colmax { + fctnx = min (k + 0.5, xmax) - max (k - 0.5, xmin) + fluxx = fluxx + fctnx * Memr[buf+k-1] + areax = areax + fctnx + } + } + + # Add the line sum to the total. + area = area + areax * fctny + flux = flux + fluxx * fctny + } + + call sfree (sp) + + # Return the appropriate error code. + if (area <= 0.0d0) + return (PY_NOPIX) + else if (ier != PY_OK) + return (ier) + else + return (PY_OK) +end + + +# AP_BYYFIT -- Measure the total flux inside a polygon while searching for +# bad pixels at the same time. + +int procedure ap_byyfit (im, xver, yver, nver, datamin, datamax, flux, area, + badpix) + +pointer im # pointer to IRAF image +real xver[ARB] # x coordinates of the vertices +real yver[ARB] # y coordinates of the vertices: +int nver # number of vertices +real datamin # minimum good data value +real datamax # maximum good data value +double flux # flux interior to the polygon +double area # approximate area of polygon +int badpix # are there bad pixels + +int i, j, k, linemin, linemax, colmin, colmax, nintr, ier +pointer sp, work1, work2, xintr, buf +real xmin, xmax, ymin, ymax, x1, x2, lx, ld +double fluxx, areax, fctnx, fctny +int ap_yclip() +pointer imgl2r() + +begin + # Check that polygon has at least 3 vertices plus a closing vertex. + if (nver < 4) { + flux = INDEFD + area = 0.0d0 + badpix = NO + return (PY_NOPOLYGON) + } + + # Allocate working space. + call smark (sp) + call salloc (work1, nver, TY_REAL) + call salloc (work2, nver, TY_REAL) + call salloc (xintr, nver, TY_REAL) + + # Find minimum and maximum y values of the polygon vertices and + # compute the minimum and maximum image line limits. + + call alimr (xver, nver, xmin, xmax) + call alimr (yver, nver, ymin, ymax) + if (xmin < 0.5 || xmax > (IM_LEN(im,1) + 0.5) || ymin < 0.5 || ymax > + (IM_LEN(im,2) + 0.5)) + ier = PY_OUTOFBOUNDS + else + ier = PY_OK + + # Find the min and max image line numbers. + ymin = max (0.5, min (real (IM_LEN (im,2) + 0.5), ymin)) + ymax = min (real (IM_LEN(im,2) + 0.5), max (0.5, ymax)) + linemin = max (1, min (int (ymin + 0.5), int (IM_LEN(im,2)))) + linemax = max (1, min (int (ymax + 0.5), int (IM_LEN(im,2)))) + + # Set up line segment parameters and initialize fluxes. + x1 = 0.5 + x2 = IM_LEN(im,1) + 0.5 + lx = x2 - x1 + flux = 0.0d0 + area = 0.0d0 + + # Loop over the range of lines of interest. + badpix = NO + do i = linemin, linemax { + + # Read in the image line. + buf = imgl2r (im, i) + if (buf == EOF) + next + + # Find all the intersection points. + if (ymin > i) + ld = min (i + 1, linemax) + else if (ymax < i) + ld = max (i - 1, linemin) + else + ld = i + nintr = ap_yclip (xver, yver, Memr[work1], Memr[work2], + Memr[xintr], nver, lx, ld) + if (nintr <= 0) + next + fctny = min (i + 0.5, ymax) - max (i - 0.5, ymin) + + # Sort the x intersection points + call asrtr (Memr[xintr], Memr[xintr], nintr) + + # Integrate the flux in the line segment + fluxx = 0.0d0 + areax = 0.0d0 + do j = 1, nintr, 2 { + + # Compute the line segment limits. + xmin = min (real (IM_LEN(im,1) + 0.5), max (0.5, + Memr[xintr+j-1])) + xmax = min (real (IM_LEN(im,1) + 0.5), max (0.5, Memr[xintr+j])) + colmin = min (int (xmin + 0.5), int (IM_LEN(im,1))) + colmax = min (int (xmax + 0.5), int (IM_LEN(im,1))) + + # Sum the contribution from a particular line segment. + do k = colmin, colmax { + fctnx = min (k + 0.5, xmax) - max (k - 0.5, xmin) + fluxx = fluxx + fctnx * Memr[buf+k-1] + areax = areax + fctnx + if (Memr[buf+k-1] < datamin || Memr[buf+k-1] > datamax) + badpix = YES + } + } + + # Add the line sum to the total. + area = area + areax * fctny + flux = flux + fluxx * fctny + } + + call sfree (sp) + + if (area <= 0.0d0) + return (PY_NOPIX) + if (badpix == YES) + return (PY_BADDATA) + else if (ier != PY_OK) + return (ier) + else + return (PY_OK) +end + + +# AP_YCLIP -- Compute the intersection of an image line with a polygon defined +# by a list of vertices. The output is a list of ranges stored in the array +# xranges. Two work additional work arrays xintr and slope are required for +# the computation. + +int procedure ap_yclip (xver, yver, xintr, slope, xranges, nver, lx, ld) + +real xver[ARB] # x vertex coords +real yver[ARB] # y vertex coords +real xintr[ARB] # work array of x intersection points +real slope[ARB] # work array of y slopes at intersection points +real xranges[ARB] # x line segments +int nver # number of vertices +real lx, ld # equation of image line + +bool collinear +int i, j, nintr, nplus, nzero, nneg, imin, imax, nadd +real u1, u2, u1u2, dx, dy, dd, xa, wa + +begin + # Compute the intersection points of the image line and the polygon. + collinear = false + nplus = 0 + nzero = 0 + nneg = 0 + nintr = 0 + #u1 = - lx * yver[1] + ld + u1 = lx * (- yver[1] + ld) + do i = 2, nver { + + #u2 = - lx * yver[i] + ld + u2 = lx * (- yver[i] + ld) + u1u2 = u1 * u2 + + # Does the polygon side intersect the image line ? + if (u1u2 <= 0.0) { + + + # Compute the x intersection coordinate if the point of + # intersection is not a vertex. + + if ((u1 != 0.0) && (u2 != 0.0)) { + + dy = yver[i-1] - yver[i] + dx = xver[i-1] - xver[i] + dd = xver[i-1] * yver[i] - yver[i-1] * xver[i] + #xa = (dx * ld - lx * dd) + xa = lx * (dx * ld - dd) + wa = dy * lx + nintr = nintr + 1 + xranges[nintr] = xa / wa + slope[nintr] = -dy + if (slope[nintr] < 0.0) + nneg = nneg + 1 + else if (slope[nintr] > 0.0) + nplus = nplus + 1 + else + nzero = nzero + 1 + collinear = false + + # For each collinear line segment add two intersection + # points. Remove interior collinear intersection points. + + } else if (u1 == 0.0 && u2 == 0.0) { + + if (! collinear) { + nintr = nintr + 1 + xranges[nintr] = xver[i-1] + if (i == 2) + slope[nintr] = yver[1] - yver[nver-1] + else + slope[nintr] = yver[i-1] - yver[i-2] + if (slope[nintr] < 0.0) + nneg = nneg + 1 + else if (slope[nintr] > 0.0) + nplus = nplus + 1 + else + nzero = nzero + 1 + nintr = nintr + 1 + xranges[nintr] = xver[i] + slope[nintr] = 0.0 + nzero = nzero + 1 + } else { + xranges[nintr] = xver[i] + slope[nintr] = 0.0 + nzero = nzero + 1 + } + collinear = true + + # If the intersection point is a vertex add it to the + # list if it is not collinear with the next point. Add + # another point to the list if the vertex is at the + # apex of an acute angle. + + } else if (u1 != 0.0) { + + if (i == nver) { + dx = (xver[2] - xver[nver]) + dy = (yver[2] - yver[nver]) + dd = dy * (yver[nver-1] - yver[nver]) + } else { + dx = (xver[i+1] - xver[i]) + dy = (yver[i+1] - yver[i]) + dd = dy * (yver[i-1] - yver[i]) + } + + # Test whether the point is collinear with the point + # ahead. If it is not include the intersection point. + + if (dy != 0.0) { + nintr = nintr + 1 + xranges[nintr] = xver[i] + slope[nintr] = yver[i] - yver[i-1] + if (slope[nintr] < 0.0) + nneg = nneg + 1 + else if (slope[nintr] > 0.0) + nplus = nplus + 1 + else + nzero = nzero + 1 + } + + # If the intersection point is an isolated vertex add + # another point to the list. + + if (dd > 0.0) { + nintr = nintr + 1 + xranges[nintr] = xver[i] + slope[nintr] = dy + if (slope[nintr] < 0.0) + nneg = nneg + 1 + else if (slope[nintr] > 0.0) + nplus = nplus + 1 + else + nzero = nzero + 1 + } + + collinear = false + + } else + collinear = false + } else + collinear = false + + u1 = u2 + } + + # Join up any split collinear line segments. + if (collinear && (slope[1] == 0.0)) { + xranges[1] = xranges[nintr-1] + slope[1] = slope[nintr-1] + nintr = nintr - 2 + nzero = nzero - 2 + } + + # Return the number of intersection points if there are no interior + # collinear line segments. + if (nzero == 0 || nplus == 0 || nneg == 0) + return (nintr) + + # Find the minimum and maximum intersection points. + call ap_alimr (xranges, nintr, u1, u2, imin, imax) + + # Check for vertices at the ends of the ranges. + + u1 = xranges[min(imin,imax)] - xranges[1] + u2 = xranges[nintr] - xranges[max(imin,imax)] + + # Vertices were traversed in order of increasing x. + if ((u1 >= 0.0 && u2 > 0.0) || (u1 > 0.0 && u2 >= 0.0) || + (u1 == u2 && imax > imin)) { + do i = imax + 1, nintr { + if (xranges[i] != xranges[i-1]) + break + imax = i + } + do i = imin - 1, 1, -1 { + if (xranges[i] != xranges[i+1]) + break + imin = i + } + } + + # Vertices were traversed in order of decreasing x. + if ((u1 <= 0.0 && u2 < 0.0) || (u1 < 0.0 && u2 <= 0.0) || + (u1 == u2 && imax < imin)) { + do i = imin + 1, nintr { + if (xranges[i] != xranges[i-1]) + break + imin = i + } + do i = imax - 1, 1, -1 { + if (xranges[i] != xranges[i+1]) + break + imax = i + } + } + + # Reorder the x ranges and slopes if necessary. + if ((imax < imin) && ! (imin == nintr && imax == 1)) { + call amovr (xranges, xintr, nintr) + do i = 1, imax + xranges[nintr-imax+i] = xintr[i] + do i = imin, nintr + xranges[i-imax] = xintr[i] + call amovr (slope, xintr, nintr) + do i = 1, imax + slope[nintr-imax+i] = xintr[i] + do i = imin, nintr + slope[i-imax] = xintr[i] + } else if ((imin < imax) && ! (imin == 1 && imax == nintr)) { + call amovr (xranges, xintr, nintr) + do i = 1, imin + xranges[nintr-imin+i] = xintr[i] + do i = imax, nintr + xranges[i-imin] = xintr[i] + call amovr (slope, xintr, nintr) + do i = 1, imin + slope[nintr-imin+i] = xintr[i] + do i = imax, nintr + slope[i-imin] = xintr[i] + } + + # Add any extra intersection points that are required to deal with + # the collinear line segments. + + nadd = 0 + for (i = 1; i <= nintr-2; ) { + if (slope[i] * slope[i+2] > 0.0) { + i = i + 2 + } else { + nadd = nadd + 1 + xranges[nintr+nadd] = xranges[i+1] + for (j = i + 3; j <= nintr; j = j + 1) { + if (slope[i] * slope[j] > 0) + break + nadd = nadd + 1 + xranges[nintr+nadd] = xranges[j-1] + } + i = j + } + } + + return (nintr + nadd) +end diff --git a/noao/digiphot/apphot/polyphot/apyfree.x b/noao/digiphot/apphot/polyphot/apyfree.x new file mode 100644 index 00000000..0bede579 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyfree.x @@ -0,0 +1,41 @@ +include "../lib/apphotdef.h" +include "../lib/polyphotdef.h" + +# AP_YFREE -- Procedure to free the polyphot structure. + +procedure ap_yfree (ap) + +pointer ap # pointer to the apphot structure + +begin + if (ap == NULL) + return + if (AP_NOISE(ap) != NULL) + call ap_noisecls (ap) + if (AP_PDISPLAY(ap) != NULL) + call ap_dispcls (ap) + if (AP_POLY(ap) != NULL) + call ap_ycls (ap) + if (AP_PSKY(ap) != NULL) + call ap_skycls (ap) + if (AP_PCENTER(ap) != NULL) + call ap_ctrcls (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_YCLS -- Procedure to close up the polyphot structure and arrays. + +procedure ap_ycls (ap) + +pointer ap # pointer to the apphot structure + +begin + if (AP_POLY(ap) == NULL) + return + call mfree (AP_POLY(ap), TY_STRUCT) +end diff --git a/noao/digiphot/apphot/polyphot/apyget.x b/noao/digiphot/apphot/polyphot/apyget.x new file mode 100644 index 00000000..a990b6f8 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyget.x @@ -0,0 +1,306 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +# AP_YGET -- Procedure to fetch the coordinates of the vertices of a +# polygon from a text file. + +int procedure ap_yget (py, im, fd, delim, x, y, max_nvertices) + +pointer py # polyphot structure +pointer im # the input image descriptor +int fd # polygon file descriptor +int delim # delimiter character +real x[ARB] # x coords of vertices +real y[ARB] # y coords of vertices +int max_nvertices # maximum number of vertices + +real xc, yc, r2max, r2, xtemp, ytemp +pointer sp, str +int i, nvertices, nreal, stat +char marker +real asumr() +int fscan(), nscan(), strncmp(), apstati() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Print prompts if input file is STDIN. + call fstats (fd, F_FILENAME, Memc[str], SZ_LINE) + if (strncmp ("STDIN", Memc[str], 5) == 0) { + call printf ("Type x and y coordinates of vertex\n") + call printf ("\t1 vertex per line\n") + call printf ("\t; to end polygon\n") + call printf ("\t^Z to end list\n") + call flush (STDOUT) + } + + # Get the polygon. + nvertices = 0 + stat = fscan (fd) + while (stat != EOF) { + + # Read the vertices from the file + call gargr (xtemp) + call gargr (ytemp) + nreal = nscan () + if (nreal != 2) { + call reset_scan () + call gargc (marker) + if (int (marker) == delim) + break + nreal = nscan () + } + + # Store the vertices. + if (nreal >= 2) { + nvertices = nvertices + 1 + if (nvertices <= max_nvertices) { + x[nvertices] = xtemp + y[nvertices] = ytemp + } + } + + stat = fscan (fd) + } + + if (nvertices == 0 && stat == EOF) { + + nvertices = EOF + call apsetr (py, PYXMEAN, INDEFR) + call apsetr (py, PYYMEAN, INDEFR) + call apsetr (py, OPYXMEAN, INDEFR) + call apsetr (py, OPYYMEAN, INDEFR) + call apsetr (py, PYCX, INDEFR) + call apsetr (py, PYCY, INDEFR) + call apsetr (py, OPYCX, INDEFR) + call apsetr (py, OPYCY, INDEFR) + call apsetr (py, PYMINRAD, INDEFR) + call apseti (py, PYNVER, 0) + + } else if (nvertices <= 2) { + + call apsetr (py, PYXMEAN, INDEFR) + call apsetr (py, PYYMEAN, INDEFR) + call apsetr (py, OPYXMEAN, INDEFR) + call apsetr (py, OPYYMEAN, INDEFR) + call apsetr (py, PYCX, INDEFR) + call apsetr (py, PYCY, INDEFR) + call apsetr (py, OPYCX, INDEFR) + call apsetr (py, OPYCY, INDEFR) + call apsetr (py, PYMINRAD, INDEFR) + call apseti (py, PYNVER, 0) + nvertices = 0 + + } else { + + # Add in the last vertex. + x[nvertices+1] = x[1] + y[nvertices+1] = y[1] + + # Transform the input coordinates. + switch (apstati(py,WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_itol (py, x, y, x, y, nvertices + 1) + case WCS_TV: + call ap_vtol (im, x, y, x, y, nvertices + 1) + default: + ; + } + + # Compute the mean polygon coordinates. + xc = asumr (x, nvertices) / nvertices + yc = asumr (y, nvertices) / nvertices + call apsetr (py, PYXMEAN, xc) + call apsetr (py, PYYMEAN, yc) + call apsetr (py, PYCX, xc) + call apsetr (py, PYCY, yc) + + # Set the minimum size of the sky annulus. + r2max = 0.0 + do i = 1, nvertices { + r2 = (x[i] - xc) ** 2 + (y[i] - yc) ** 2 + if (r2 > r2max) + r2max = r2 + } + + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, xc, yc, xc, yc, 1) + case WCS_TV: + call ap_ltov (im, xc, yc, xc, yc, 1) + default: + ; + } + call apsetr (py, OPYXMEAN, xc) + call apsetr (py, OPYYMEAN, yc) + call apsetr (py, OPYCX, xc) + call apsetr (py, OPYCY, yc) + + call apsetr (py, PYMINRAD, (sqrt (r2max) + 1.1)) + call apseti (py, PYNVER, nvertices) + } + + call sfree (sp) + return (nvertices) +end + + +# AP_YMKPOLY -- Mark the coordinates of a polygon on the display device. + +int procedure ap_ymkpoly (py, im, id, x, y, max_nvertices, idflush) + +pointer py # polyphot structure +pointer im # the input image descriptor +pointer id # display pointer +real x[ARB] # x coords of vertices +real y[ARB] # y coords of vertices +int max_nvertices # maximum number of vertices +int idflush # flush the imd interface ? + +int i, nvertices, stat, wcs, key +pointer sp, cmd +real xtemp, ytemp, xc, yc, r2max, r2 +real apstatr(), asumr() +int clgcur(), apstati() +errchk gscur + +begin + # Reopen the device. + if (id != NULL) + call greactivate (id, 0) + + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Type prompt. + call printf ( + "Mark polygon vertex [space=mark,q=quit].\n") + stat = clgcur ("icommands", xtemp, ytemp, wcs, key, Memc[cmd], SZ_LINE) + call ap_vtol (im, xtemp, ytemp, xtemp, ytemp, 1) + + + # Fetch the polygon and draw it on the display. + nvertices = 0 + while (stat != EOF) { + + if (key == 'q') + break + + # Decode and draw vertices. + nvertices = nvertices + 1 + if (nvertices <= max_nvertices) { + x[nvertices] = xtemp + y[nvertices] = ytemp + if (id != NULL) { + if (nvertices == 1) + call gamove (id, x[1], y[1]) + else { + call gadraw (id, x[nvertices], y[nvertices]) + if (idflush == YES) + call gframe (id) + else + call gflush (id) + } + } + } else + break + + # Type prompt. + call printf ( + "Mark polygon vertex [space=mark,q=quit].\n") + stat = clgcur ("icommands", xtemp, ytemp, wcs, key, Memc[cmd], + SZ_LINE) + call ap_vtol (im, xtemp, ytemp, xtemp, ytemp, 1) + } + call printf ("\n") + + call sfree (sp) + + # Return EOF or the number of vertices in the polygon. + if (stat == EOF) { + + call apsetr (py, PYXMEAN, INDEFR) + call apsetr (py, PYYMEAN, INDEFR) + call apsetr (py, OPYXMEAN, INDEFR) + call apsetr (py, OPYYMEAN, INDEFR) + call apsetr (py, PYCX, INDEFR) + call apsetr (py, PYCY, INDEFR) + call apsetr (py, OPYCX, INDEFR) + call apsetr (py, OPYCY, INDEFR) + call apsetr (py, PYMINRAD, INDEFR) + call apseti (py, PYNVER, 0) + nvertices = EOF + + } else if (nvertices <= 2) { + + call apsetr (py, PYXMEAN, INDEFR) + call apsetr (py, PYYMEAN, INDEFR) + call apsetr (py, OPYXMEAN, INDEFR) + call apsetr (py, OPYYMEAN, INDEFR) + call apsetr (py, PYCX, INDEFR) + call apsetr (py, PYCY, INDEFR) + call apsetr (py, OPYCX, INDEFR) + call apsetr (py, OPYCY, INDEFR) + call apsetr (py, PYMINRAD, INDEFR) + call apseti (py, PYNVER, 0) + nvertices = 0 + + } else { + + # Add the last vertex and close the polygon. + x[nvertices+1] = x[1] + y[nvertices+1] = y[1] + if (id != NULL) { + call gadraw (id, x[1], y[1]) + if (idflush == YES) + call gframe (id) + else + call gflush (id) + } + + # Compute and save the mean polygon coords. + xc = asumr (x, nvertices) / nvertices + yc = asumr (y, nvertices) / nvertices + call apsetr (py, PYXMEAN, xc) + call apsetr (py, PYYMEAN, yc) + call apsetr (py, PYCX, xc) + call apsetr (py, PYCY, yc) + if (id != NULL) + call gscur (id, apstatr (py, PYCX), apstatr (py, PYCY)) + + # Compute the mean sky annulus. + r2max = 0.0 + do i = 1, nvertices { + r2 = (x[i] - xc) ** 2 + (y[i] - yc) ** 2 + if (r2 > r2max) + r2max = r2 + } + + # Compute output coordinate centers. + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, xc, yc, xc, yc, 1) + case WCS_TV: + call ap_ltov (im, xc, yc, xc, yc, 1) + default: + ; + } + call apsetr (py, OPYXMEAN, xc) + call apsetr (py, OPYYMEAN, yc) + call apsetr (py, OPYCX, xc) + call apsetr (py, OPYCY, yc) + + call apseti (py, PYNVER, nvertices) + call apsetr (py, PYMINRAD, (sqrt (r2max) + 1.1)) + + } + + if (id != NULL) + call gdeactivate (id, 0) + + return (nvertices) +end diff --git a/noao/digiphot/apphot/polyphot/apyinit.x b/noao/digiphot/apphot/polyphot/apyinit.x new file mode 100644 index 00000000..6d9c2fe7 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyinit.x @@ -0,0 +1,64 @@ +include "../lib/apphotdef.h" +include "../lib/polyphotdef.h" + +# AP_YINIT - Initialize the polyphot structure. + +procedure ap_yinit (ap, cfunction, cbox, sfunction, annulus, dannulus, fwhmpsf, + noise) + +pointer ap # pointer to the apphot structure +int cfunction # centering algorithm +real cbox # centering box half-width +int sfunction # sky fitting algorithm +real annulus # inner radius of sky annulus +real dannulus # width of sky annulus +real fwhmpsf # fwhmpsf +int noise # Noise function + +begin + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set the global apphot package parameters. + call ap_defsetup (ap, fwhmpsf) + + # Setup noise parameters. + call ap_noisesetup (ap, noise) + + # Set display options. + call ap_dispsetup (ap) + + # Setup the centering parameters. + call ap_ctrsetup (ap, cfunction, cbox) + + # Set up the sky fitting parameters. + call ap_skysetup (ap, sfunction, annulus, dannulus) + + # Setup the polyphot parameters. + call ap_ysetup (ap) + + # Set unused structure pointers to null. + AP_PPSF(ap) = NULL + AP_PPHOT(ap) = NULL + AP_RPROF(ap) = NULL +end + + +# AP_YSETUP -- Procedure to set parameters for polygonal aperture photometry. + +procedure ap_ysetup (ap) + +pointer ap # pointer to apphot strucuture + +pointer ply + +begin + call malloc (AP_POLY(ap), LEN_PYSTRUCT, TY_STRUCT) + ply = AP_POLY(ap) + AP_PYCX(ply) = INDEFR + AP_PYCY(ply) = INDEFR + AP_PYX(ply) = INDEFR + AP_PYY(ply) = INDEFR + AP_PYZMAG(ply) = DEF_PYZMAG + AP_PYNAME(ply) = EOS + AP_PYROOT(ply) = EOS +end diff --git a/noao/digiphot/apphot/polyphot/apymkfree.x b/noao/digiphot/apphot/polyphot/apymkfree.x new file mode 100644 index 00000000..424b94e3 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apymkfree.x @@ -0,0 +1,23 @@ +include "../lib/apphotdef.h" + +# AP_YMKFREE -- Procedure to free the polyphot structure. + +procedure ap_ymkfree (ap) + +pointer ap # pointer to the apphot structure + +begin + if (ap == NULL) + return + if (AP_NOISE(ap) != NULL) + call ap_noisecls (ap) + if (AP_PDISPLAY(ap) != NULL) + call ap_dispcls (ap) + if (AP_POLY(ap) != NULL) + call ap_ycls (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 diff --git a/noao/digiphot/apphot/polyphot/apymkinit.x b/noao/digiphot/apphot/polyphot/apymkinit.x new file mode 100644 index 00000000..215d3e9c --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apymkinit.x @@ -0,0 +1,31 @@ +include "../lib/apphotdef.h" +include "../lib/noise.h" + +# AP_YMKINIT - Procedure to initialize the polymark structure. + +procedure ap_ymkinit (ap) + +pointer ap # pointer to the apphot structure + +begin + call calloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set the main structure parameters. + call ap_defsetup (ap, 2.5) + + # Set the noise options. + call ap_noisesetup (ap, AP_NPOISSON) + + # Set display options. + call ap_dispsetup (ap) + + # Setup the polyphot parameters. + call ap_ysetup (ap) + + # Set unused structure pointers to null. + AP_PCENTER(ap) = NULL + AP_PSKY(ap) = NULL + AP_PPHOT(ap) = NULL + AP_RPROF(ap) = NULL + AP_PPSF(ap) = NULL +end diff --git a/noao/digiphot/apphot/polyphot/apynextobj.x b/noao/digiphot/apphot/polyphot/apynextobj.x new file mode 100644 index 00000000..292d84bd --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apynextobj.x @@ -0,0 +1,190 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/polyphot.h" + +# AP_YNEXTOBJ -- Read the next polygon from the list polygon and/or coordinate +# list. + +int procedure ap_ynextobj (py, im, id, pl, cl, delim, x, y, maxnver, prev_num, + req_num, ld, pd) + +pointer py # pointer to the apphot structre +pointer im # the input image descriptor +pointer id # pointer to image display stream +int pl # polygons file descriptor +int cl # coordinates file descriptor +int delim # delimiter character +real x[ARB] # x coordinates of the polygon vertices +real y[ARB] # y coordinates of the polygon vertices +int maxnver # maximum number of vertices +int prev_num # previous object +int req_num # requested object +int ld # current object +int pd # current polygon + +real xshift, yshift +pointer sp, fname +int stdin, nskip, ncount, nver, stat +real apstatr() +int strncmp(), ap_yget(), ap_ycoords(), apstati() +errchk greactivate, gdeactivate, gscur + +begin + # The coordinate list is undefined. + if (cl == NULL) { + + # Get the input file name. + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call fstats (pl, F_FILENAME, Memc[fname], SZ_FNAME) + + # Compute the number of polygons that must be skipped. + if (strncmp ("STDIN", Memc[fname], 5) == 0) { + stdin = YES + nskip = 1 + } else { + stdin = NO + if (req_num <= prev_num) { + call seek (pl, BOF) + nskip = req_num + } else + nskip = req_num - prev_num + } + + # Initialize the polygon search. + ncount = 0 + pd = prev_num + + # Find the correct polygon. + repeat { + + call apsetr (py, PYX, apstatr (py, PYCX)) + call apsetr (py, PYY, apstatr (py, PYCY)) + nver = ap_yget (py, im, pl, delim, x, y, maxnver) + + if (nver == EOF) { + ncount = EOF + } else if (nver > 0) { + ncount = ncount + 1 + pd = pd + 1 + } + + } until (ncount == EOF || ncount == nskip) + + # Return the polygon number. + if (req_num <= prev_num) + pd = ncount + ld = pd + + call sfree (sp) + + if (ncount == EOF) { + return (EOF) + } else if (id != NULL) { + iferr { + call greactivate (id, 0) + call gscur (id, apstatr (py, PYCX), apstatr (py, PYCY)) + call gdeactivate (id, 0) + } then + ; + return (nver) + } else { + return (nver) + } + + # Both the polygon and coordinate file are defined. + } else { + + # Get the input file name. + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call fstats (cl, F_FILENAME, Memc[fname], SZ_FNAME) + + # Compute the number of objects that must be skipped. + if (strncmp ("STDIN", Memc[fname], 5) == 0) { + stdin = YES + nskip = 1 + } else { + stdin = NO + if (req_num <= prev_num) { + call seek (cl, BOF) + call seek (pl, BOF) + pd = 0 + nskip = req_num + } else + nskip = req_num - prev_num + } + + # Initialize the polygon id. + ncount = 0 + ld = prev_num + + # Find the correct object and shift the coordinates. + repeat { + + # Read the coordinates and the polygon. + stat = ap_ycoords (cl, delim, xshift, yshift, stdin) + if (stat == EOF) { + ncount = EOF + } else if (stat == NEXT_POLYGON || pd == 0) { + call apsetr (py, PYX, apstatr (py, PYCX)) + call apsetr (py, PYY, apstatr (py, PYCY)) + nver = ap_yget (py, im, pl, delim, x, y, maxnver) + if (nver == EOF) + ncount = EOF + else if (nver > 0) + pd = pd + 1 + } + + # Shift the polygon coordinates. + if (stat == THIS_OBJECT && ncount != EOF && nver > 0) { + switch (apstati(py,WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_itol (py, xshift, yshift, xshift, yshift, 1) + case WCS_TV: + call ap_vtol (im, xshift, yshift, xshift, yshift, 1) + default: + ; + } + call aaddkr (x, (xshift - apstatr (py, PYCX)), x, nver + 1) + call aaddkr (y, (yshift - apstatr (py, PYCY)), y, nver + 1) + call apsetr (py, PYCX, xshift) + call apsetr (py, PYCY, yshift) + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, xshift, yshift, xshift, yshift, 1) + case WCS_TV: + call ap_ltov (im, xshift, yshift, xshift, yshift, 1) + default: + ; + } + call apsetr (py, OPYCX, xshift) + call apsetr (py, OPYCY, yshift) + ncount = ncount + 1 + ld = ld + 1 + } + + } until (ncount == EOF || ncount == nskip) + + # Get the correct id. + if (req_num <= prev_num) + ld = ncount + + call sfree (sp) + + if (ncount == EOF) { + return (EOF) + } else if (id != NULL) { + iferr { + call greactivate (id, 0) + call gscur (id, apstatr (py, PYCX), apstatr (py, PYCY)) + call gdeactivate (id, 0) + } then + ; + return (nver) + } else { + return (nver) + } + } + +end diff --git a/noao/digiphot/apphot/polyphot/apypars.x b/noao/digiphot/apphot/polyphot/apypars.x new file mode 100644 index 00000000..4f988323 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apypars.x @@ -0,0 +1,20 @@ +# AP_PYPARS -- Procedure to write the current polyphot parameters to the +# output file. + +procedure ap_pypars (ap) + +pointer ap # pointer to apphot structure + +begin + # Write the data dependent parameters. + call ap_dapars (ap) + + # Write the centering parameters. + call ap_cepars (ap) + + # Write the sky fitting parameters. + call ap_sapars (ap) + + # Set the photometry parameters. + call ap_popars (ap) +end diff --git a/noao/digiphot/apphot/polyphot/apyphot.x b/noao/digiphot/apphot/polyphot/apyphot.x new file mode 100644 index 00000000..cd7bac0a --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyphot.x @@ -0,0 +1,577 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +define HELPFILE "apphot$polyphot/polyphot.key" + +# AP_YPHOT -- Compute the flux inside polygonal apertures. + +int procedure ap_yphot (py, im, cl, pl, gd, id, out, stid, interactive, cache) + +pointer py # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # coordinate file descriptor +int pl # starlist file descriptor +pointer gd # pointer to graphcis stream +pointer id # pointer to image display stream +int out # output file descriptor +int stid # output file sequence number +int interactive # interactive mode +int cache # cache the input image pixels + +real wx, wy +pointer sp, x, y, xout, yout, cmd +int nvertices, cier, sier, pier, wcs, key, ip, colonkey +int prev_num, req_num, ptid, ltid, delim, newlist, newimage +int newcenterbuf, newcenter, newskybuf, newsky, newmagbuf, newmag +int req_size, old_size, buf_size, memstat + +real apstatr() +int ap_ymkpoly(), ap_yfit(), clgcur(), apfitsky(), aprefitsky() +int apstati(), ctoi(), ap_ynextobj(), ap_ycenter(), ap_yrecenter() +int apgqverify(), apgtverify(), ap_yradsetup(), apnew(), ap_avsky() +int ap_memstat(), sizeof() +bool fp_equalr() +data delim /';'/ + +define endswitch_ 99 + +begin + # Allocate temporary space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (x, MAX_NVERTICES + 1, TY_REAL) + call salloc (y, MAX_NVERTICES + 1, TY_REAL) + call salloc (xout, MAX_NVERTICES + 1, TY_REAL) + call salloc (yout, MAX_NVERTICES + 1, TY_REAL) + + # Initialize the cursor read. + key = ' ' + Memc[cmd] = EOS + + # Initialize the fit. + 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 + ptid = 0 + ltid = 0 + + # Loop over the polygon file. + nvertices = 0 + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Store the current coordinates. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (py, CWX, wx) + call apsetr (py, CWY, wy) + + # Has the polygon moved ? + if (apnew (py, wx, wy, apstatr (py, PYCX), apstatr (py, PYCY), + newlist) == YES) { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + + # Loop over the colon commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("polyphot", py, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + # Print the help page. + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Plot a centered stellar radial profile. + case 'd': + if (interactive == YES) + call ap_qrad (py, im, wx, wy, gd) + + # Setup the parameters using a radial profile plot + case 'i': + if (interactive == YES) { + nvertices = ap_yradsetup (py, im, id, gd, out, stid, + Memr[x], Memr[y], MAX_NVERTICES) + newlist = NO + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + + # Verify the critical parameters. + case 'v': + call ap_yconfirm (py, out, stid) + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + + # Define a polygon interactively. + case 'g': + if (interactive == YES) { + if (id == gd) + nvertices = ap_ymkpoly (py, im, id, Memr[x], Memr[y], + MAX_NVERTICES, NO) + else + nvertices = ap_ymkpoly (py, im, id, Memr[x], Memr[y], + MAX_NVERTICES, YES) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + newlist = NO + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + } + + # Print error messages. + case 'e': + if (interactive == YES) + call ap_pyerrors (py, cier, sier, pier) + + # Rewind the polygon and coordinate lists. + case 'r': + if (pl != NULL) { + call seek (pl, BOF) + if (cl != NULL) + call seek (cl, BOF) + ptid = 0 + ltid = 0 + } else if (interactive == YES) + call printf ("No polygon list\n") + + # Show/set polyphot parameters. + 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 polyphot parameters + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call ap_ycolon (py, im, pl, cl, out, stid, ptid, ltid, + Memc[cmd], newimage, newcenterbuf, newcenter, + newskybuf, newsky, newmagbuf, newmag) + goto endswitch_ + } + + # Measure the next object in the list. + if (pl != NULL) { + if (interactive == YES) + call printf ("No polygon list\n") + goto endswitch_ + } + + # Get the next polygon. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, Memr[x], + Memr[y], MAX_NVERTICES, prev_num, req_num, ltid, ptid) + if (nvertices == EOF) { + if (interactive == YES) + call printf ( + "End of polygon list, use r key to rewind\n") + goto endswitch_ + } + + # Move to the next object. + newlist = YES + if (colonkey == 'm') { + newcenterbuf = YES; newcenter = YES + newskybuf = YES; newsky = YES + newmagbuf = YES; newmag = YES + goto endswitch_ + } + + # Fit the center, sky and measure polygon. + cier = ap_ycenter (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), Memr[x], Memr[y], nvertices + 1) + sier = apfitsky (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), NULL, gd) + pier = ap_yfit (py, im, Memr[x], Memr[y], nvertices + 1, + apstatr (py, SKY_MODE), apstatr (py, SKY_SIGMA), + apstati (py, NSKY)) + if (interactive == YES) + call ap_qyprint (py, cier, sier, pier) + + # Draw the polygon. + if (id != NULL) { + call appymark (py, id, Memr[x], Memr[y], nvertices + 1, + apstati (py, MKCENTER), apstati (py, MKSKY), + apstati (py, MKPOLYGON)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + # Write the results to output. + if (stid == 1) + call ap_param (py, out, "polyphot") + + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + default: + call amovr (Memr[x], Memr[xout], nvertices + 1) + call amovr (Memr[y], Memr[yout], nvertices + 1) + } + if (newlist == YES) + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, ltid, ptid, cier, sier, pier) + else + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, 0, ptid, cier, sier, pier) + + # Set up for the next object. + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Show/set polyphot parameters. + default: + call ap_ycolon (py, im, pl, cl, out, stid, ptid, ltid, + Memc[cmd], newimage, newcenterbuf, newcenter, + newskybuf, newsky, newmagbuf, newmag) + } + + # Reestablish the image display viewport and window. + if (newimage == YES) { + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[cmd], im, 4) + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + } + newimage = NO + + # Get measure next object in the list. + case 'm', 'n': + if (pl == NULL) { + if (interactive == YES) + call printf ("No polygon list\n") + goto endswitch_ + } + + # Get the next polygon. + prev_num = ltid + req_num = ltid + 1 + nvertices = ap_ynextobj (py, im, id, pl, cl, delim, Memr[x], + Memr[y], MAX_NVERTICES, prev_num, req_num, ltid, ptid) + if (nvertices == EOF) { + if (interactive == YES) + call printf ( + "End of polygon list, use r key to rewind\n") + goto endswitch_ + } + + # 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 = ap_ycenter (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), Memr[x], Memr[y], nvertices + 1) + sier = apfitsky (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), NULL, gd) + pier = ap_yfit (py, im, Memr[x], Memr[y], nvertices + 1, + apstatr (py, SKY_MODE), apstatr (py, SKY_SIGMA), + apstati (py, NSKY)) + if (interactive == YES) + call ap_qyprint (py, cier, sier, pier) + + # Draw the polygon. + if (id != NULL) { + call appymark (py, id, Memr[x], Memr[y], nvertices + 1, + apstati (py, MKCENTER), apstati (py, MKSKY), + apstati (py, MKPOLYGON)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + # Write the results to output. + if (stid == 1) + call ap_param (py, out, "polyphot") + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, Memr[x], Memr[y], Memr[xout], Memr[yout], + nvertices + 1) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[xout], Memr[yout], + nvertices + 1) + default: + call amovr (Memr[x], Memr[xout], nvertices + 1) + call amovr (Memr[y], Memr[yout], nvertices + 1) + } + if (newlist == YES) + call ap_yprint (py, out, Memr[xout], Memr[yout], nvertices, + stid, ltid, ptid, cier, sier, pier) + else + call ap_yprint (py, out, Memr[xout], Memr[yout], nvertices, + stid, 0, ptid, cier, sier, pier) + + # Set up for the next object. + stid = stid + 1 + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Process the remainder of the list. + case 'l': + if (pl != NULL) { + call ap_ybphot (py, im, cl, pl, out, stid, ltid, ptid, id, + YES) + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + } else if (interactive == YES) + call printf ("No polygon list\n") + + # Save current parameters in the pset file. + case 'w': + call ap_pypars (py) + + # Center the current polygon around the cursor poition. + case 'c': + if (newcenterbuf == YES) + cier = ap_ycenter (py, im, wx, wy, Memr[x], Memr[y], + nvertices + 1) + else if (newcenter == YES) + cier = ap_yrecenter (py, im, Memr[x], Memr[y], + nvertices + 1, cier) + if (interactive == YES) + call ap_qcenter (py, cier) + if (id != NULL) { + call apmark (py, id, apstati (py, MKCENTER), NO, NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + newcenterbuf = NO; newcenter = NO + + # Fit the sky around the cursor position. + case 't': + if (newskybuf == YES || ! fp_equalr (wx, + apstatr (py, SXCUR)) || ! fp_equalr (wy, apstatr (py, + SYCUR))) + sier = apfitsky (py, im, wx, wy, NULL, gd) + else if (newsky == YES) + sier = aprefitsky (py, im, gd) + if (interactive == YES) + call ap_qspsky (py, sier) + if (id != NULL) { + call apmark (py, id, NO, apstati (py, MKSKY), NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + newskybuf = NO; newsky = NO + + # Compute the average of several sky measurements around + # different cursor postions. + case 'a': + sier = ap_avsky (py, im, stid, NULL, id, gd, interactive) + if (interactive == YES) + call ap_qaspsky (py, sier) + newskybuf = NO; newsky = NO + + # Fit sky in an annulus around the current polygon. + case 's': + if (newskybuf == YES || ! fp_equalr (apstatr (py, PYCX), + apstatr (py, SXCUR)) || ! fp_equalr (apstatr (py, PYCY), + apstatr (py, SYCUR))) + sier = apfitsky (py, im, apstatr (py, PYCX), apstatr (py, + PYCY), NULL, gd) + else if (newsky == YES) + sier = aprefitsky (py, im, gd) + if (interactive == YES) + call ap_qspsky (py, sier) + if (id != NULL) { + call apmark (py, id, NO, apstati (py, MKSKY), NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + newskybuf = NO; newsky = NO + + # Compute the magnitudes of current polygon using the current sky. + case 'p', 'o': + if (newcenterbuf == YES) + cier = ap_ycenter (py, im, wx, wy, Memr[x], Memr[y], + nvertices + 1) + else if (newcenter == YES) + cier = ap_yrecenter (py, im, Memr[x], Memr[y], + nvertices + 1, cier) + pier = ap_yfit (py, im, Memr[x], Memr[y], nvertices + 1, + apstatr (py, SKY_MODE), apstatr (py, SKY_SIGMA), + apstati (py, NSKY)) + if (interactive == YES) + call ap_qyprint (py, cier, sier, pier) + if (id != NULL) { + call appymark (py, id, Memr[x], Memr[y], nvertices + 1, + NO, NO, apstati (py, MKPOLYGON)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + newcenterbuf = NO; newcenter = NO + newmagbuf = NO; newmag = NO + + # Write the results. + if (key == 'o') { + if (stid == 1) + call ap_param (py, out, "polyphot") + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + default: + call amovr (Memr[x], Memr[xout], nvertices + 1) + call amovr (Memr[y], Memr[yout], nvertices + 1) + } + if (newlist == YES) + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, ltid, ptid, cier, sier, pier) + else + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, 0, ptid, cier, sier, pier) + stid = stid + 1 + } + + # Center, compute the sky and the magnitudes and save the results. + case 'h', 'j', 'f', ' ': + + # Compute the centers. + if (key == 'f' || key == ' ') { + if (newcenterbuf == YES) + cier = ap_ycenter (py, im, wx, wy, Memr[x], Memr[y], + nvertices + 1) + else if (newcenter == YES) + cier = ap_yrecenter (py, im, Memr[x], Memr[y], + nvertices + 1, cier) + } else { + if (newcenterbuf == YES) + cier = ap_ycenter (py, im, apstatr (py, PYCX), + apstatr (py, PYCY), Memr[x], Memr[y], nvertices + 1) + else if (newcenter == YES) + cier = ap_yrecenter (py, im, Memr[x], Memr[y], + nvertices + 1, cier) + } + + # Compute the sky values and the magnitudes. + if (newskybuf == YES || ! fp_equalr (apstatr (py, + PYCX), apstatr (py, SXCUR)) || ! fp_equalr (apstatr (py, + PYCY), apstatr (py, SYCUR))) sier = apfitsky (py, im, + apstatr (py, PYCX), apstatr (py, PYCY), NULL, gd) + else if (newsky == YES) + sier = aprefitsky (py, im, gd) + pier = ap_yfit (py, im, Memr[x], Memr[y], nvertices + 1, + apstatr (py, SKY_MODE), apstatr (py, SKY_SIGMA), + apstati (py, NSKY)) + + if (interactive == YES) + call ap_qyprint (py, cier, sier, pier) + if (id != NULL) { + call appymark (py, id, Memr[x], Memr[y], nvertices + 1, + apstati (py, MKCENTER), apstati (py, MKSKY), + apstati (py, MKPOLYGON)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + newcenterbuf = NO; newcenter = NO + newskybuf = NO; newsky = NO + newmagbuf = NO; newmag = NO + + # Write the results to an output file. + if (key == ' ' || key == 'j') { + if (stid == 1) + call ap_param (py, out, "polyphot") + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + case WCS_TV: + call ap_ltov (im, Memr[x], Memr[y], Memr[xout], + Memr[yout], nvertices + 1) + default: + call amovr (Memr[x], Memr[xout], nvertices + 1) + call amovr (Memr[y], Memr[yout], nvertices + 1) + } + if (newlist == YES) + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, ltid, ptid, cier, sier, pier) + else + call ap_yprint (py, out, Memr[xout], Memr[yout], + nvertices, stid, 0, ptid, cier, sier, pier) + stid = stid + 1 + } + + default: + call printf ("Unknown or ambigous keystroke command\n") + } + +endswitch_ + # Reset the keystroke and command defaults. + key = ' ' + Memc[cmd] = EOS + call apsetr (py, WX, apstatr (py, CWX)) + call apsetr (py, WY, apstatr (py, CWY)) + } + +end diff --git a/noao/digiphot/apphot/polyphot/apyprint.x b/noao/digiphot/apphot/polyphot/apyprint.x new file mode 100644 index 00000000..dda24c3f --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyprint.x @@ -0,0 +1,94 @@ +include "../lib/apphot.h" +include "../lib/center.h" +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +# AP_YPRINT -- Write the results of the polyphot task to the output file. + +procedure ap_yprint (py, fd, xver, yver, nver, id, lid, pid, cier, sier, pier) + +pointer py # pointer to apphot structure +int fd # output file descriptor +real xver[ARB] # coords of x vertices +real yver[ARB] # coords of y vertices +int nver # number of vertices +int id # id number of str +int lid # list id number of star +int pid # polygon number +int cier # centering error +int sier # sky fitting error +int pier # photometric error + +real apstatr() + +begin + if (fd == NULL) + return + + # Write object id parameters. + call ap_wid (py, fd, apstatr (py, OXINIT), apstatr (py, OYINIT), id, + lid, '\\') + + # Write the centering parameters + call ap_wcres (py, fd, cier, '\\') + + # Write the fitsky parameters. + call ap_wsres (py, fd, sier, '\\') + + # Write the polyphot parameters + call ap_wlres (py, fd, xver, yver, nver, pid, pier) +end + + +# AP_QYPRINT -- Print a quick summary of the polyphot results on the standard +# output. + +procedure ap_qyprint (py, cier, sier, pier) + +pointer py # pointer to apphot structure +int cier # centering error +int sier # sky fitting error +int pier # photometry error + +pointer sp, imname +real apstatr() + +begin + # Allocate space. + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + # Print polyphot magnitudes. + #call apstats (py, IMNAME, Memc[imname], SZ_FNAME) + call apstats (py, IMROOT, Memc[imname], SZ_FNAME) + call printf ("%s %8.2f %8.2f %8g ") + call pargstr (Memc[imname]) + call pargr (apstatr (py, OPYCX)) + call pargr (apstatr (py, OPYCY)) + call pargr (apstatr (py, SKY_MODE)) + call printf ("%7.3f %s\n") + call pargr (apstatr (py, PYMAG)) + if (cier != AP_OK || sier != AP_OK || pier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + + call sfree (sp) +end + + +# AP_YHDR -- Write the polyphot header banner to the output file. + +procedure ap_yhdr (ap, fd) + +pointer ap # apphot structure +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_plhdr (ap, fd) +end diff --git a/noao/digiphot/apphot/polyphot/apyradsetup.x b/noao/digiphot/apphot/polyphot/apyradsetup.x new file mode 100644 index 00000000..a966110b --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyradsetup.x @@ -0,0 +1,135 @@ +include "../lib/fitsky.h" +include "../lib/polyphot.h" + +define HELPFILE "apphot$polyphot/ipolyphot.key" + +# AP_YRADSETUP - Set up polyphot interactively using a radial profile +# plot of a bright star. + +int procedure ap_yradsetup (ap, im, id, gd, out, stid, x, y, max_nvertices) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +pointer id # pointer to the image display +pointer gd # pointer to graphics stream +int out # output file descriptor +int stid # output file sequence number +real x[ARB] # array of x vertices +real y[ARB] # array of y vertices +int max_nvertices # maximum number of vertices + +int nvertices, cier, sier, pier, key, wcs +pointer sp, str, cmd +real rmin, rmax, imin, imax, u1, u2, v1, v2, x1, x2, y1, y2 +real xcenter, ycenter, xc, yc, rval + +int ap_ycenter(), clgcur(), ap_showplot() +int apfitsky(), ap_yfit(), apstati(), ap_ymkpoly() +real apstatr(), ap_cfwhmpsf(), ap_ccapert(), ap_cannulus(), ap_csigma() +real ap_cdannulus(), ap_cdatamin(), ap_cdatamax() +real ap_crgrow(), ap_crclean(), ap_crclip() + +begin + # Mark the polygon interactively. + if (id == gd) + nvertices = ap_ymkpoly (ap, im, id, x, y, max_nvertices, NO) + else + nvertices = ap_ymkpoly (ap, im, id, x, y, max_nvertices, YES) + if (id != NULL) { + if (gd == id) + call gflush (id) + else + call gframe (id) + } + if (nvertices <= 0) + return (nvertices) + + # Store the viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Check for open display device and graphics stream. + if (gd == NULL) + return (0) + call greactivate (gd, 0) + + # Show the plot. + if (ap_showplot (ap, im, apstatr (ap, PYCX), apstatr (ap, PYCY), gd, + xcenter, ycenter, rmin, rmax, imin, imax) == ERR) { + call gdeactivate (gd, 0) + return (nvertices) + } + + # Allocate memory. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + while (clgcur ("gcommands", xc, yc, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + switch (key) { + + case 'q': + break + case '?': + call gpagefile (gd, HELPFILE, "") + case 'f': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'h': + #rval = ap_ccthresh (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'c': + rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax) + case 's': + rval = ap_csigma (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 '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 '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 '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_cannulus (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_cdannulus (ap, gd, out, stid, apstatr (ap, ANNULUS), + 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) + + # Clean up memory space. + call gdeactivate (gd, 0) + call sfree (sp) + + # Print the answer. + cier = ap_ycenter (ap, im, xcenter, ycenter, x, y, nvertices + 1) + sier = apfitsky (ap, im, apstatr (ap, PYCX), apstatr (ap, PYCY), + NULL, gd) + pier = ap_yfit (ap, im, x, y, nvertices + 1, apstatr (ap, SKY_MODE), + apstatr (ap, SKY_SIGMA), apstati (ap, NSKY)) + call ap_qyprint (ap, cier, sier, pier) + + return (nvertices) +end diff --git a/noao/digiphot/apphot/polyphot/apyshift.x b/noao/digiphot/apphot/polyphot/apyshift.x new file mode 100644 index 00000000..4fd56cbe --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apyshift.x @@ -0,0 +1,41 @@ +include "../lib/apphot.h" +include "../lib/polyphot.h" + +# AP_YSHIFT -- Shift a polygon to a new center position. + +procedure ap_yshift (py, im, x, y, nver, newx, newy) + +pointer py # pointer to the polyphot structure +pointer im # the input image descriptor +real x[ARB] # x coordinates of the vertices +real y[ARB] # y coordinates of the vertices +int nver # number of vertices +real newx, newy # new x and y coordinates of center + +real xshift, yshift +real apstatr() +int apstati() + +begin + call apsetr (py, PYX, apstatr (py, PYCX)) + call apsetr (py, PYY, apstatr (py, PYCY)) + + xshift = newx - apstatr (py, PYCX) + yshift = newy - apstatr (py, PYCY) + call aaddkr (x, xshift, x, nver) + call aaddkr (y, yshift, y, nver) + call apsetr (py, PYCX, newx) + call apsetr (py, PYCY, newy) + + switch (apstati(py,WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, newx, newy, xshift, yshift, 1) + case WCS_TV: + call ap_ltov (im, newx, newy, xshift, yshift, 1) + default: + xshift = newx + yshift = newy + } + call apsetr (py, OPYCX, xshift) + call apsetr (py, OPYCY, yshift) +end diff --git a/noao/digiphot/apphot/polyphot/apywrite.x b/noao/digiphot/apphot/polyphot/apywrite.x new file mode 100644 index 00000000..b246aac3 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/apywrite.x @@ -0,0 +1,92 @@ +include "../lib/apphot.h" +include "../lib/polyphot.h" + +# AP_YWRITE -- Procedure to write a polygon to a file. + +procedure ap_ywrite (py, im, cl, pl, x, y, nvertices, cid, pid, firstpoly, + newpoly) + +pointer py # pointer to the apphot structure +pointer im # the input image descriptor +int cl # pointer to the coordinate list file +int pl # pointer to the polygon list file +real x[ARB] # x coordinates of the vertices +real y[ARB] # y coordinates of the vertices +int nvertices # number of vertices +int cid # coordinate list index +int pid # polygon list index +int firstpoly # first polygon measured +int newpoly # new polygon + +real xtemp, ytemp +int i +real apstatr() +int apstati() + +begin + # Make sure the output files are at EOF. + if (pl == NULL) + return + call seek (pl, EOF) + if (cl != NULL) + call seek (cl, EOF) + + if (newpoly == YES) { + + # Terminate the coord list that belongs with the first polygon. + if (firstpoly == NO && cl != NULL) + call fprintf (cl, ";\n") + + # Write out the coordinates. + do i = 1, nvertices { + switch (apstati(py, WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, x[i], y[i], xtemp, ytemp, 1) + case WCS_TV: + call ap_ltov (im, x[i], y[i], xtemp, ytemp, 1) + default: + xtemp = x[i] + ytemp = y[i] + } + call fprintf (pl, "%g %g\n") + call pargr (xtemp) + call pargr (ytemp) + } + if (nvertices > 0) + call fprintf (pl, ";\n") + + pid = pid + 1 + + # Reset polygon parameters. + newpoly = NO + if (firstpoly == YES) + firstpoly = NO + } + + # Write out the central coordinates of the polygon. + if (firstpoly == NO && cl != NULL) { + + switch (apstati(py, WCSOUT)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (py, apstatr(py,CWX), apstatr(py,CWY), xtemp, + ytemp, 1) + case WCS_TV: + call ap_ltov (im, apstatr(py,CWX), apstatr(py,CWY), xtemp, + ytemp, 1) + default: + xtemp = apstatr (py, CWX) + ytemp = apstatr (py, CWY) + } + call fprintf (cl, "%g %g\n") + call pargr (xtemp) + call pargr (ytemp) + + cid = cid + 1 + } + + # Flush the output files. + if (pl != NULL) + call flush (pl) + if (cl != NULL) + call flush (cl) +end diff --git a/noao/digiphot/apphot/polyphot/ipolyphot.key b/noao/digiphot/apphot/polyphot/ipolyphot.key new file mode 100644 index 00000000..f8ef878c --- /dev/null +++ b/noao/digiphot/apphot/polyphot/ipolyphot.key @@ -0,0 +1,16 @@ + Interactive Photometry Setup Menu + + v Mark and verify the critical parameters (f,c,s,a,d) + + f Mark and verify the psf full-width half-maximum + s Mark and verify the standard deviation of the background + l Mark and verify the minimum good data value + u Mark and verify the maximum good data value + + c Mark and verify the centering box width + n Mark and verify the cleaning radius + p Mark and verify the clipping radius + + a Mark and verify the inner radius of the sky annulus + d Mark and verify the width of the sky annulus + g Mark and verify the region growing radius diff --git a/noao/digiphot/apphot/polyphot/mkpkg b/noao/digiphot/apphot/polyphot/mkpkg new file mode 100644 index 00000000..b98b6803 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/mkpkg @@ -0,0 +1,60 @@ +# POLYPHOT and POLYMARK routines + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + # POLYPHOT routines + + appyerrors.x ../lib/polyphot.h + apgypars.x ../lib/display.h ../lib/fitsky.h \ + ../lib/noise.h ../lib/polyphot.h \ + ../lib/center.h + apycenter.x ../lib/center.h ../lib/polyphot.h + apyconfirm.x ../lib/apphot.h ../lib/noise.h \ + ../lib/center.h ../lib/fitsky.h + apycoords.x ../lib/polyphot.h + apyfree.x ../lib/apphotdef.h ../lib/polyphotdef.h + apyinit.x ../lib/apphotdef.h ../lib/polyphotdef.h + apynextobj.x <fset.h> ../lib/apphot.h \ + ../lib/polyphot.h + apypars.x + apyradsetup.x ../lib/fitsky.h ../lib/polyphot.h + apyshift.x ../lib/apphot.h ../lib/polyphot.h + apycolon.x ../lib/polyphot.h ../lib/center.h \ + ../lib/apphot.h ../lib/fitsky.h \ + ../lib/noise.h ../lib/display.h \ + <gset.h> <error.h> + apyfit.x ../lib/polyphot.h ../lib/noise.h \ + ../lib/apphot.h <imhdr.h> \ + <mach.h> + apyget.x <fset.h> ../lib/apphot.h \ + ../lib/fitsky.h ../lib/polyphot.h + polyshow.x ../lib/display.h ../lib/polyphot.h + apybphot.x <fset.h> ../lib/display.h \ + ../lib/polyphot.h ../lib/fitsky.h \ + ../lib/apphot.h + apyphot.x <ctype.h> <gset.h> \ + ../lib/display.h ../lib/polyphot.h \ + ../lib/fitsky.h ../lib/apphot.h \ + <imhdr.h> + apyprint.x ../lib/fitsky.h ../lib/center.h \ + ../lib/apphot.h ../lib/polyphot.h + t_polyphot.x <gset.h> <fset.h> \ + ../lib/apphot.h ../lib/polyphot.h \ + <imhdr.h> + + # POLYMARK routines + + t_polymark.x <gset.h> <fset.h> \ + "../lib/apphot.h" "../lib/polyphot.h" \ + <imhdr.h> + apmkpylist.x <gset.h> <ctype.h> \ + "../lib/apphot.h" "../lib/polyphot.h" + apydraw.x "../lib/polyphot.h" + apymkinit.x "../lib/apphotdef.h" "../lib/noise.h" + apymkfree.x "../lib/apphotdef.h" + apywrite.x "../lib/apphot.h" "../lib/polyphot.h" + ; diff --git a/noao/digiphot/apphot/polyphot/polymark.key b/noao/digiphot/apphot/polyphot/polymark.key new file mode 100644 index 00000000..9e1901a1 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/polymark.key @@ -0,0 +1,16 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +d Plot radial profile of star near cursor +g Define the current polygonal aperture +f Draw the current polygon on the display +spbar Draw the current polygon on the display, output the polygon +r Rewind the polygon list +m Draw the next polygon in the polygon list on the display +l Draw all the remaining polygons in the list on the display +q Exit + + Colon commands + +:m [n] Draw the next [nth] polygon in the polygon list on the display diff --git a/noao/digiphot/apphot/polyphot/polyphot.key b/noao/digiphot/apphot/polyphot/polyphot.key new file mode 100644 index 00000000..e2c92da2 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/polyphot.key @@ -0,0 +1,110 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Store the current parameters +d Plot radial profile of current object +i Define current polygon, graphically set parameters using current object +g Define current polygon +c Fit center for current object +t Fit sky around cursor +a Average sky values fit around several cursor positions +s Fit sky around current object +h Do photometry for current polygon +j Do photometry for current polygon, output results +p Do photometry for current object using current sky +o Do photometry for current object using current sky, output results +f Do photometry for current object +spbar Do photometry for current object, output results +m Move to next object in coordinate list +n Do photometry for next object in coordinate list, output results +l Do photometry for remaining objects in coordinate list, output results +r Rewind the the coordinate list +e Print error messages +q Exit task + + + Colon Commands + +:show [data/center/sky/phot] List the parameters +:m [n] Move to next [nth] object in coordinate list +:n [n] Do photometry for next [nth] object in coordinate list, output results + + + Colon Parameter Editing Commands + +# Image and file name parameters + +:image [string] Image name +:polygon [string] Polygon file +:coords [string] Coordinate file +:output [string] Results file + +# 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 pixel value (counts) +:datamax [value] Maximum good pixel value (counts) + +# Noise parameters + +:noise [string] Noise model (constant|poisson) +:gain [string] Gain image header keyword +:ccdread [string] Readout noise image header keyword +:epadu [value] Gain (electrons per count) +:epadu [value] Readout noise (electrons) + +# Observing 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 centering box (scale units) +:cthreshold [value] Centering intensity threshold (sigma) +:cmaxiter [value] Maximum number of iterations +:maxshift [value] Maximum center shift (scale units) +:minsnratio [value] Minimum S/N ratio for centering +:clean [y/n] Clean subraster before centering +:rclean [value] Cleaning radius (scale units) +:rclip [value] Clipping radius (scale units) +:kclean [value] Clean K-sigma rejection limit (sigma) + +# Sky fitting algorithm parameters + +:salgorithm [string] Sky fitting algorithm +:skyvalue [value] User supplied sky value (counts) +:annulus [value] Inner radius of sky annulus (scale units) +:dannulus [value] Width of sky annulus (scale units) +:khist [value] Sky histogram extent (+/- sigma) +:binsize [value] Resolution of sky histogram (sigma) +:sloclip [value] Low-side clipping factor in percent +:shiclip [value] High-side clipping factor in percent +:smooth [y/n] Lucy smooth the sky histogram +: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 + +:zmag [value] Zero point of magnitude scale + +# Plotting and marking parameters + +:mkcenter [y/n] Mark computed centers on the display +:mksky [y/n] Mark the sky annuli on the display +:mkpolygon [y/n] Mark the polygon on the display diff --git a/noao/digiphot/apphot/polyphot/polyshow.x b/noao/digiphot/apphot/polyphot/polyshow.x new file mode 100644 index 00000000..9e3da070 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/polyshow.x @@ -0,0 +1,65 @@ +include "../lib/display.h" +include "../lib/polyphot.h" + +# AP_YSHOW -- Print the all the polyphot parameters on the standard output. + +procedure ap_yshow (ap) + +pointer ap # pointer to the apphot strucuture + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + call ap_spshow (ap) + call printf ("\n") + call ap_yyshow (ap) +end + + +# AP_YPSHOW -- Print the noise and polypars parameters on the standard output. + +procedure ap_ypshow (ap) + +pointer ap # pointer to apphot structure + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_yyshow (ap) +end + + +# AP_YYSHOW -- Print the polypars parameters on the standard output. + +procedure ap_yyshow (ap) + +pointer ap # pointer to apphot structure + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call printf ("Photometry Parameters\n") + + call apstats (ap, PYNAME, Memc[str], SZ_LINE) + call printf (" %s = %s\n") + call pargstr (KY_PYNAME) + call pargstr (Memc[str]) + + call printf (" %s = %g\n") + call pargstr (KY_PYZMAG) + call pargr (apstatr (ap, PYZMAG)) + + call printf (" %s = %b\n") + call pargstr (KY_MKPOLYGON) + call pargb (itob (apstati (ap, MKPOLYGON))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/polyphot/t_polymark.x b/noao/digiphot/apphot/polyphot/t_polymark.x new file mode 100644 index 00000000..4eea5207 --- /dev/null +++ b/noao/digiphot/apphot/polyphot/t_polymark.x @@ -0,0 +1,240 @@ +include <gset.h> +include <fset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/polyphot.h" + +define MAX_NVERTICES 100 + +# T_POLYMARK -- Create a polygons list file interactively. + +procedure t_polymark() + +pointer image # pointer to name of image +pointer polygons # pointer to the polygons file +pointer coords # pointer to the coords file +int cache # cache the input image pixels in memory +pointer display # pointer to display device name +pointer graphics # pointer to the graphics device + +pointer sp, cfname, pfname, im, py, id, gd, str +int limlist, lplist, lclist, stat, pl, cl, root, pid, cid, wcs +int imlist, plist, clist, newpy, newcoo, memstat, req_size, old_size +int buf_size + +pointer gopen(), immap() +int imtlen(), clplen(), imtgetim(), clgfil(), strncmp(), strlen() +int fnldir(), ap_mkpylist(), imtopenp(), clpopnu(), open(), clgwrd() +int access(), btoi(), ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen() + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (polygons, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (cfname, SZ_FNAME, TY_CHAR) + call salloc (pfname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Set STDOUT. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get input and output file lists. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + plist = clpopnu ("polygons") + lplist = clplen (plist) + clist = clpopnu ("coords") + lclist = clplen (clist) + + # Check that image input and coords file list lengths match. + if (lclist > 0 && lclist != limlist) { + call imtclose (imlist) + call clpcls (plist) + call clpcls (clist) + call error (0, "Imcompatible image and coord file list lengths") + } + + # Check that image input and polygons file list lengths match. + if (lplist != limlist) { + call imtclose (imlist) + call clpcls (plist) + call clpcls (clist) + call error (0, "Imcompatible image and polygons file list lengths") + } + + cache = btoi (clgetb ("cache")) + + # Open the graphics device. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + 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 + } + } + + # Open the display device. + call clgstr ("display", Memc[display], SZ_FNAME) + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[display], Memc[graphics])) { + id = gd + } else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + + # Open the polymark structure. + call ap_ymkinit (py) + + # 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 (py, 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 (py, WCSOUT, wcs) + + # Mark polygons. + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open image. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (py, im, Memc[image]) + + # Establish the image display viewport and window. + 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) + + # Set the polygon file name. + if (lplist == 0) + pl = NULL + else { + stat = clgfil (plist, Memc[polygons], SZ_FNAME) + root = fnldir (Memc[polygons], Memc[pfname], SZ_FNAME) + if (strncmp ("default", Memc[polygons+root], 7) == 0 || root == + strlen (Memc[polygons])) { + call apoutname (Memc[image], Memc[pfname], "ver", + Memc[pfname], SZ_FNAME) + lplist = limlist + } else if (stat != EOF) + call strcpy (Memc[polygons], Memc[pfname], SZ_FNAME) + if (access (Memc[pfname], READ_WRITE, TEXT_FILE) == YES) { + newpy = NO + pl = open (Memc[pfname], READ_WRITE, TEXT_FILE) + } else { + newpy = YES + pl = open (Memc[pfname], NEW_FILE, TEXT_FILE) + call close (pl) + pl = open (Memc[pfname], READ_WRITE, TEXT_FILE) + } + } + call apsets (py, PYNAME, Memc[pfname]) + + # Set the coord file name. + if (lclist == 0) + cl = NULL + else { + stat = clgfil (clist, Memc[coords], SZ_FNAME) + root = fnldir (Memc[coords], Memc[cfname], SZ_FNAME) + if (strncmp ("default", Memc[coords+root], 7) == 0 || root == + strlen (Memc[coords])) { + call apoutname (Memc[image], Memc[cfname], "coo", + Memc[cfname], SZ_FNAME) + lclist = limlist + } else if (stat != EOF) + call strcpy (Memc[coords], Memc[cfname], SZ_FNAME) + if (access (Memc[cfname], READ_WRITE, TEXT_FILE) == YES) { + newcoo = NO + cl = open (Memc[cfname], READ_WRITE, TEXT_FILE) + } else { + newcoo = YES + cl = open (Memc[cfname], NEW_FILE, TEXT_FILE) + call close (cl) + cl = open (Memc[cfname], READ_WRITE, TEXT_FILE) + } + } + call apsets (py, CLNAME, Memc[cfname]) + + # Mark polygons on each image. + pid = 1 + cid = 1 + stat = ap_mkpylist (im, py, pl, cl, id, gd, pid, cid) + + # Unmap the input image. + call imunmap (im) + + # Close the polygon file. + if (pl != NULL) { + call close (pl) + if (newpy == YES && pid <= 1) + call delete (Memc[pfname]) + } + + # Close the coordinate file. + if (cl != NULL) { + if (cid > 1) { + call seek (cl, EOF) + call fprintf (cl, ";\n") + } + call close (cl) + if (newcoo == YES && cid <= 1) + call delete (Memc[cfname]) + } + + # Uncache memory. + call fixmem (old_size) + + if (stat == YES) + break + } + + # Close the graphics and image display devices. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + # Free the polymark structure. + call ap_ymkfree (py) + + # Close image, polygon, and coordinate lists. + call imtclose (imlist) + call clpcls (plist) + call clpcls (clist) + call sfree (sp) +end diff --git a/noao/digiphot/apphot/polyphot/t_polyphot.x b/noao/digiphot/apphot/polyphot/t_polyphot.x new file mode 100644 index 00000000..5242d02e --- /dev/null +++ b/noao/digiphot/apphot/polyphot/t_polyphot.x @@ -0,0 +1,329 @@ +include <gset.h> +include <fset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/polyphot.h" + +# T_POLYPHOT -- Measure the total magnitudes inside a list of polygons. + +procedure t_polyphot() + +pointer image # pointer to name of image +pointer output # pointer to the results file +pointer coords # pointer to file of coordinates +pointer polygon # pointer to file containing polygon +pointer graphics # pointer to graphics device name +pointer display # pointer to display device name +int interactive # mode of use +int cache # cache the input image pixels +int verify # verify critical parameters +int update # update the critical parameters +int verbose # print messages + +pointer sp, outfname, cname, im, py, id, gd, str +int limlist, lplist, lolist, lclist, sid, lid, pid, pl, cl, out, root, stat +int imlist, plist, olist, clist, memstat, wcs, req_size, old_size, buf_size + +pointer immap(), gopen() +int imtlen(), imtgetim(), clplen(), clgfil(), btoi(), strncmp() +int fnldir(), strlen(), ap_yphot(), open(), imtopenp(), clpopnu() +int clgwrd(), ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (polygon, 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 STDOUT. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get input and output file lists. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + plist = clpopnu ("polygons") + lplist = clplen (plist) + olist = clpopnu ("output") + lolist = clplen (olist) + clist = clpopnu ("coords") + lclist = clplen (clist) + + # Check that image and polygon list lengths match. + if (limlist < 1 || (lplist > 1 && lplist != limlist)) { + call imtclose (imlist) + call clpcls (plist) + call clpcls (olist) + call clpcls (clist) + call error (0, "Imcompatible image and polygon list lengths") + } + + # Check that image and coordinates list lengths match. + if (limlist < 1 || (lclist > 1 && lclist != limlist)) { + call imtclose (imlist) + call clpcls (plist) + call clpcls (olist) + call clpcls (clist) + call error (0, "Imcompatible image and coordinate list lengths") + } + + # Check that image input and output list lengths match. + if (lolist > 1 && lolist != limlist) { + call imtclose (imlist) + call clpcls (plist) + call clpcls (olist) + call clpcls (clist) + 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 + interactive = btoi (clgetb ("interactive")) + cache = btoi (clgetb ("cache")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + verbose = btoi (clgetb ("verbose")) + + # Get polygon fitting parameters. + call ap_gypars (py) + + # Confirm the algorithm parameters. + if (verify == YES && interactive == NO) { + call ap_yconfirm (py, NULL, 1) + if (update == YES) + call ap_pypars (py) + } + + # 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 (py, 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 (py, 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 + } + + # Measure flux in a polygon. + sid = 1 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open image. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (py, im, Memc[image]) + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + + # Open the polygons file. + if (lplist <= 0) { + pl = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (plist, Memc[polygon], SZ_FNAME) + root = fnldir (Memc[polygon], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[polygon+root], 7) == 0 || root == + strlen (Memc[polygon])) { + call ap_inname (Memc[image], Memc[outfname], "ver", + Memc[outfname], SZ_FNAME) + lplist = limlist + pl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else if (stat != EOF) { + call strcpy (Memc[polygon], Memc[outfname], SZ_FNAME) + pl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else { + call apstats (py, PYNAME, Memc[outfname], SZ_FNAME) + call seek (pl, BOF) + } + } + call apsets (py, PYNAME, Memc[outfname]) + call apfroot (Memc[outfname], Memc[str], SZ_LINE) + call apsets (py, PYROOT, Memc[str]) + + # Open the coordinates file. + 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 (py, CLNAME, Memc[outfname], SZ_FNAME) + call seek (cl, BOF) + } + } + call apsets (py, CLNAME, Memc[outfname]) + call apfroot (Memc[outfname], Memc[str], SZ_LINE) + call apsets (py, CLROOT, Memc[str]) + + # Set output file name. + 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], "ply", + 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 (py, OUTNAME, Memc[outfname], SZ_FNAME) + } + call apsets (py, OUTNAME, Memc[outfname]) + + # Do the photometry. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = ap_yphot (py, im, cl, pl, NULL, NULL, out, sid, NO, + cache) + else if (pl != NULL) { + lid = 0 + pid = 0 + call ap_ybphot (py, im, cl, pl, out, sid, lid, pid, NULL, + verbose) + stat = NO + } else + stat = NO + } else + stat = ap_yphot (py, im, cl, pl, gd, id, out, sid, YES, cache) + + # Unmap the input image. + call imunmap (im) + + # Close the polygon file. + if (pl != NULL) { + if (lplist > 1) + call close (pl) + } + + # Close the coordinate file. + if (cl != NULL) { + if (lclist > 1) + call close (cl) + } + + # Close the output file. + if (out != NULL && lolist != 1) { + call close (out) + if (sid <= 1) { + call apstats (py, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + sid = 1 + } + + # Uncache memory. + call fixmem (old_size) + + if (stat == YES) + break + } + + # Close plot files. + if (id == gd && id != NULL) + call gclose (id) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + + # Close the single coords and polygon file. + if (pl != NULL && lplist == 1) + call close (pl) + if (cl != NULL && lclist == 1) + call close (cl) + + # Close the singled output file. + if (out != NULL && lolist == 1) { + call close (out) + if (sid <= 1) { + call apstats (py, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + } + + # Close up the files. + call ap_yfree (py) + + # Close image, coord and shift lists. + call imtclose (imlist) + call clpcls (plist) + call clpcls (clist) + call clpcls (olist) + + call sfree (sp) +end |