diff options
Diffstat (limited to 'noao/digiphot/apphot/fitpsf')
22 files changed, 2475 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/fitpsf/apbfitpsf.x b/noao/digiphot/apphot/fitpsf/apbfitpsf.x new file mode 100644 index 00000000..9854ab9f --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apbfitpsf.x @@ -0,0 +1,95 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/display.h" + +# APBFITPSF -- Procedure to fit a functional form to the PSF for a list of +# objects. + +procedure apbfitpsf (ap, im, cl, out, gid, id, ld, interactive) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +int out # output file descriptor +pointer gid # image display stream +int id # initial id +int ld # list number +int interactive # interactive mode + +int ier, ild, stdin +pointer sp, str +real wx, wy +int fscan, nscan(), strncmp(), apsffit(), apstati() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call fstats (cl, F_FILENAME, Memc[str], SZ_FNAME) + + # Initialize. + ild = ld + + # Print query if input is STDIN. + if (strncmp ("STDIN", Memc[str], 5) == 0) { + stdin = YES + call printf ("Type object x and y coordinates (^Z or ^D to end): ") + call flush (STDOUT) + } else + stdin = NO + + # Loop over the interesting objects. + while (fscan(cl) != EOF) { + + # Decode the centers. + call gargr (wx) + call gargr (wy) + if (nscan () != 2) { + if (stdin == YES) { + call printf ( + "Type object x and y coordinates (^Z or ^D to end): ") + call flush (STDOUT) + } + next + } + + # Transform the input coordinates. + switch (apstati(ap,WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_itol (ap, wx, wy, wx, wy, 1) + case WCS_TV: + call ap_vtol (im, wx, wy, wx, wy, 1) + default: + ; + } + + # Store the current cursor coordinates. + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Fit the psf. + ier = apsffit (ap, im, wx, wy) + + # Output the results. + if (interactive == YES) { + call ap_qppsf (ap, ier) + if (gid != NULL) + call appfmark (ap, gid, apstati (ap, MKPSFBOX)) + } + if (id == 1) + call ap_param (ap, out, "fitpsf") + call ap_ppsf (ap, out, id, ild, ier) + + # Prepare for the next object. + id = id + 1 + ild = ild + 1 + call apsetr (ap, WX, wx) + call apsetr (ap, WY, wy) + if (stdin == YES) { + call printf ( + "Type object x and y coordinates (^D or ^Z to end): ") + call flush (STDOUT) + } + + } + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/apfbuf.x b/noao/digiphot/apphot/fitpsf/apfbuf.x new file mode 100644 index 00000000..428c30cb --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apfbuf.x @@ -0,0 +1,88 @@ +include <imhdr.h> +include "../lib/apphotdef.h" +include "../lib/fitpsfdef.h" +include "../lib/fitpsf.h" + +# APFBUF -- Procedure to fetch the subraster of pixels to be fitted given the +# pointer to the IRAF image, the coordinates of the center and the apphot +# structure. + +int procedure apfbuf (ap, im, wx, wy) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # center coordinates + +int ippix +pointer psf +real ppix +pointer ap_psfpix() + +begin + # Check for 0 sized aperture. + psf = AP_PPSF(ap) + if (AP_PSFAPERT(psf) <= 0.0) + return (AP_NOPSFAREA) + + # Get the PSF pixels. + ppix = max (1.0, AP_PSFAPERT(psf) * AP_SCALE(ap)) + ippix = 2 * int (ppix) + 1 + AP_PSFPIX(psf) = ap_psfpix (im, wx, wy, ippix, AP_PXC(psf), AP_PYC(psf), + AP_PNX(psf), AP_PNY(psf)) + if (AP_PSFPIX(psf) == NULL) + return (AP_NOPSFAREA) + else if (AP_PNX(psf) < ippix || AP_PNY(psf) < ippix) + return (AP_PSF_OUTOFBOUNDS) + else + return (AP_OK) +end + + +# AP_PSFPIX -- Procedure to fetch the pixels to be used for centering. + +pointer procedure ap_psfpix (im, wx, wy, papert, xc, yc, nx, ny) + +pointer im # pointer to IRAF image +real wx, wy # center of subraster to be extracted +int papert # width of subraster to be extracted +real xc, yc # center of extracted subraster +int nx, ny # dimensions of extracted subraster + +int ncols, nlines, c1, c2, l1, l2, half_papert +pointer buf +real xc1, xc2, xl1, xl2 +pointer imgs2r() + +begin + # Check for nonsensical input. + half_papert = (papert - 1) / 2 + if (half_papert <= 0) + return (NULL) + + # Test for out of bounds pixels + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + xc1 = wx - half_papert + xc2 = wx + half_papert + xl1 = wy - half_papert + xl2 = wy + half_papert + if (xc1 > real (ncols) || xc2 < 1.0 || xl1 > real (nlines) || xl2 < 1.0) + return (NULL) + + # Get column and line limits, dimensions and center of subraster. + c1 = max (1.0, min (real (ncols), xc1)) + 0.5 + c2 = min (real (ncols), max (1.0, xc2)) + 0.5 + l1 = max (1.0, min (real (nlines), xl1)) + 0.5 + l2 = min (real (nlines), max (1.0, xl2)) + 0.5 + nx = c2 - c1 + 1 + ny = l2 - l1 + 1 + xc = wx - c1 + 1 + yc = wy - l1 + 1 + + # Get pixels and return. + buf = imgs2r (im, c1, c2, l1, l2) + if (buf == EOF) + return (NULL) + else + return (buf) +end diff --git a/noao/digiphot/apphot/fitpsf/apfitpsf.x b/noao/digiphot/apphot/fitpsf/apfitpsf.x new file mode 100644 index 00000000..94957b6a --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apfitpsf.x @@ -0,0 +1,344 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/fitpsf.h" + +define HELPFILE "apphot$fitpsf/fitpsf.key" + +# APFITPSF -- Procedure to fit a functional form to the PSF for a list of +# objects in interactive mode. + +int procedure apfitpsf (ap, im, cl, gd, id, out, stid, interactive, cache) + +pointer ap # pointer to apphot structure +pointer im # pointer to IRAF image +int cl # starlist file descriptor +pointer gd # graphics pointer +pointer id # display pointer +int out # output file descriptor +int stid # output file sequence number +int interactive # interactive mode +int cache # cache the input image pixels + +real wx, wy, xlist, ylist +pointer sp, cmd +int wcs, key, newimage, newbuf, newfit, newlist, ltid, ier +int ip, oid, colonkey, prev_num, req_num, buf_size, req_size, old_size +int memstat + +real apstatr() +int clgcur(), apgscur(), apsffit(), apsfrefit(), apstati() +int apgqverify(), apgtverify(), ctoi(), apnew(), ap_memstat(), sizeof() + +define endswitch_ 99 + +begin + # Initialize. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize the cursor commands. + key = ' ' + Memc[cmd] = EOS + + # Initialize the fitting commands. + newimage = NO + newbuf = YES + newfit = YES + ier = AP_OK + + # Initialize the sequencing. + newlist = NO + ltid = 0 + + # Loop over the cursor commands. + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Store the cursor coordinates + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Test to see if the cursor moved. + if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) { + newbuf = YES + newfit = YES + } + + # Loop over the cursor commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("fitpsf", ap, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + + # Print error messages. + case 'e': + if (interactive == YES) + call ap_pferrors (ap, ier) + + # Print the help pages. + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]") + + # Rewind the list. + case 'r': + if (cl != NULL) { + call seek (cl, BOFL) + ltid = 0 + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Process the fitpsf colon commands. + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + colonkey = Memc[cmd+ip-1] + + switch (colonkey) { + case 'm', 'n': + + # Show/set fitpsf parameters. + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call apsfcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newbuf, newfit) + goto endswitch_ + } + + # No coordinate list. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Get next object from the list. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + + # Fetch the next object from the list. + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, + ltid) == EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + } + + # Convert the coordinates. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to the next object. + newlist = YES + if (colonkey == 'm') { + newbuf = YES + newfit = YES + goto endswitch_ + } + + # Measure the next object. + ier = apsffit (ap, im, xlist, ylist) + if (id != NULL) { + call appfmark (ap, id, apstati (ap, MKPSFBOX)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (interactive == YES) + call ap_qppsf (ap, ier) + if (stid == 1) + call ap_param (ap, out, "fitpsf") + call ap_ppsf (ap, out, stid, ltid, ier) + stid = stid + 1 + newbuf = NO + newfit = NO + + default: + call apsfcolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newbuf, newfit) + } + + # Reestablish the image viewport and window. + if (newimage == YES) { + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[cmd], im, 4) + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + } + + newimage = NO + + # Plot a centered stellar radial profile. + case 'd': + if (interactive == YES) { + call ap_qrad (ap, im, wx, wy, gd) + newbuf = YES + newfit = YES + } + + # Verify the parameters interactively. + case 'v': + call ap_pfconfirm (ap, out, stid) + newbuf = YES + newfit = YES + + # Interactively set up fitpsf parameters. + case 'i': + if (interactive == YES) { + call ap_pfradsetup (ap, im, wx, wy, gd, out, stid) + newbuf = YES + newfit = YES + } + + # Save the parameters. + case 'w': + call ap_ppfpars (ap) + + # Fit the PSF and save the results. + case 'f', ' ': + if (newbuf == YES) + ier = apsffit (ap, im, wx, wy) + else if (newfit == YES) + ier = apsfrefit (ap, im) + if (id != NULL) { + call appfmark (ap, id, apstati (ap, MKPSFBOX)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (interactive == YES) + call ap_qppsf (ap, ier) + newbuf = NO + newfit = NO + + if (key == ' ') { + if (stid == 1) + call ap_param (ap, out, "fitpsf") + if (newlist == YES) + call ap_ppsf (ap, out, stid, ltid, ier) + else + call ap_ppsf (ap, out, stid, 0, ier) + stid = stid + 1 + } + + + # Get, measure the next star in the coordinate list. + case 'm', 'n': + + # No coordinate file. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + prev_num = ltid + req_num = ltid + 1 + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, + ltid) == EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + } + + # Convert coordinates if necessary. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to next object. + newlist = YES + if (key == 'm') { + newbuf = YES + newfit = YES + goto endswitch_ + } + + # Measure next object. + ier = apsffit (ap, im, xlist, ylist) + if (id != NULL) { + call appfmark (ap, id, apstati (ap, MKPSFBOX)) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + if (interactive == YES) + call ap_qppsf (ap, ier) + if (stid == 1) + call ap_param (ap, out, "fitpsf") + call ap_ppsf (ap, out, stid, ltid, ier) + stid = stid + 1 + newbuf = NO + newfit = NO + + # Process the remainder of the coordinate list. + case 'l': + if (cl != NULL) { + oid = stid + ltid = ltid + 1 + call apbfitpsf (ap, im, cl, out, id, stid, ltid, YES) + ltid = ltid + stid - oid + 1 + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + } else if (interactive == YES) + call printf ("No coordinate list\n") + + default: + # do nothing + call printf ("Unknown keystroke command\n") + } + +endswitch_ + # Setup for the next object by setting the default keystroke + # command and storing the old cursor coordinates in the + # centering structure. + + key = ' ' + Memc[cmd] = EOS + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/apgpfpars.x b/noao/digiphot/apphot/fitpsf/apgpfpars.x new file mode 100644 index 00000000..ed4b93ab --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apgpfpars.x @@ -0,0 +1,40 @@ +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/fitpsf.h" + +# AP_GPFPARS -- Procedure to fetch the fitpsf parameters. + +procedure ap_gpfpars (ap) + +pointer ap # pointer to the apphot structure + +int function +pointer sp, str +bool clgetb() +int clgeti(), btoi(), clgwrd() +real clgetr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Open the apphot structure. + call apsfinit (ap, AP_RADGAUSS, 3.5, 2.0, AP_NPOISSON) + call apsetr (ap, PSFAPERT, clgetr ("box") / 2.0) + + # Get the data dependent parameters. + call ap_gdapars (ap) + + # Get the rest of the FITPSF fitting parameters. + function = clgwrd ("function", Memc[str], SZ_LINE, PSFFUNCS) + call apsets (ap, PSFSTRING, Memc[str]) + call apseti (ap, PSFUNCTION, function) + call apseti (ap, PMAXITER, clgeti ("maxiter")) + call apseti (ap, PNREJECT, clgeti ("nreject")) + call apsetr (ap, PK2, clgetr ("kreject")) + + # Get the plotting parameters. + call apseti (ap, MKPSFBOX, btoi (clgetb ("mkbox"))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/appfconfirm.x b/noao/digiphot/apphot/fitpsf/appfconfirm.x new file mode 100644 index 00000000..30e111b2 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/appfconfirm.x @@ -0,0 +1,58 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/fitpsf.h" + +# AP_PFCONFIRM -- Procedure to confirm the critical fitpsf parameters. + +procedure ap_pfconfirm (ap, out, stid) + +pointer ap # pointer to the apphot structure +int out # the output file descriptor +int stid # the output file sequence number + +pointer sp, str +real fwhmpsf, psfapert, datamin, datamax +int apstati() +real apstatr(), ap_vfwhmpsf(), ap_vpsfapert() +real ap_vdatamin(), ap_vdatamax() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + call printf ("\n") + + # Confirm the fitting function. + call ap_vpfstring (ap, Memc[str], SZ_FNAME) + + # Confirm the fwhmpsf. + if (apstati (ap, PSFUNCTION) != AP_MOMENTS) + fwhmpsf = ap_vfwhmpsf (ap) + else + fwhmpsf = apstatr (ap, FWHMPSF) + + # Confirm the fitting box. + psfapert = 2.0 * ap_vpsfapert (ap) + + # Confirm the good data minimum and maximum values. + datamin = ap_vdatamin (ap) + datamax = ap_vdatamax (ap) + + call printf ("\n") + + # Update the database file. + if (out != NULL && stid > 1) { + call ap_sparam (out, KY_PSFSTRING, Memc[str], UN_PSFMODEL, + "psf fitting function") + call ap_rparam (out, KY_FWHMPSF, fwhmpsf, UN_ASCALEUNIT, + "full width half maximum of the psf") + call ap_rparam (out, KY_PSFAPERT, psfapert, UN_PSFSCALEUNIT, + "width of fitting box") + 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/fitpsf/appferrors.x b/noao/digiphot/apphot/fitpsf/appferrors.x new file mode 100644 index 00000000..20ed8820 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/appferrors.x @@ -0,0 +1,26 @@ +include "../lib/fitpsf.h" + +# AP_PFERRORS -- Print detailed fitpsf error messages on the standard output. + +procedure ap_pferrors (ap, ier) + +pointer ap # pointer to apphot structure (unused) +int ier # integer error code + +begin + switch (ier) { + case AP_NOPSFAREA: + call printf ( + "The psf fitting aperture is outside the image.\n") + case AP_PSF_OUTOFBOUNDS: + call printf ( + "The psf fitting aperture is partially outside the image.\n") + case AP_NPSF_TOO_SMALL: + call printf ( + "The number of data points is too few for psf fitting.\n") + case AP_PSF_SINGULAR: + call printf ("The psf fitting solution is singular.\n") + case AP_PSF_NOCONVERGE: + call printf ("The psf fitting solution did not converge.\n") + } +end diff --git a/noao/digiphot/apphot/fitpsf/appfradsetup.x b/noao/digiphot/apphot/fitpsf/appfradsetup.x new file mode 100644 index 00000000..aaae9be6 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/appfradsetup.x @@ -0,0 +1,86 @@ +define HELPFILE "apphot$fitpsf/ifitpsf.key" + +# AP_PFRADSETUP -- Procedure to set up fitpsf interactively + +procedure ap_pfradsetup (ap, im, wx, wy, gd, out, stid) + +pointer ap # pointer to apphot structure +pointer im # pointer to to the IRAF image +real wx, wy # cursor coordinates +pointer gd # pointer to graphics stream +int out # output file descriptor +int stid # output file sequence number + +int ier, wcs, key +pointer sp, cmd +real xcenter, ycenter, rmin, rmax, imin, imax, xc, yc, rval +real u1, u2, v1, v2, x1, x2, y1, y2 + +int apsffit(), clgcur(), ap_showplot() +real ap_cfwhmpsf(), ap_cpapert(), ap_cdatamin(), ap_cdatamax() + +begin + if (gd == NULL) + return + call greactivate (gd, 0) + call gclear (gd) + + # Save old viewport and window coordinates + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # plot the radial profile. + if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax, + imin, imax) == ERR) { + call gdeactivate (gd, 0) + return + } + + # Allocate temporary space. + call smark (sp) + 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 '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 'b': + rval = ap_cpapert (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'v': + rval = ap_cpapert (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + default: + call printf ("Unknown or ambiguous keystroke command\007\n") + } + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + } + + # Interactive setup is complete. + call printf ( + "Interactive setup is complete. Type w to save parameters.\n") + + # Restore old viewport and window coords. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + call gdeactivate (gd, 0) + call sfree (sp) + + # Fit the object. + ier = apsffit (ap, im, xcenter, ycenter) + call ap_qppsf (ap, ier) +end diff --git a/noao/digiphot/apphot/fitpsf/apppfpars.x b/noao/digiphot/apphot/fitpsf/apppfpars.x new file mode 100644 index 00000000..58f40bcf --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apppfpars.x @@ -0,0 +1,34 @@ +include "../lib/fitpsf.h" +include "../lib/display.h" + +# AP_PPFPARS -- Procedure to write the fitpsf parameters to a parameter file. + +procedure ap_ppfpars (ap) + +pointer ap # pointer to apphot structure + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + # Initialize and open psets. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Store the psf fitting parameters. + call clputr ("box", 2.0 * apstatr (ap, PSFAPERT)) + call apstats (ap, PSFSTRING, Memc[str], SZ_FNAME) + call clpstr ("function", Memc[str]) + call clputi ("maxiter", apstati (ap, PMAXITER)) + call clputr ("kreject", apstatr (ap, PK2)) + call clputi ("nreject", apstati (ap, PNREJECT)) + call clputb ("mkbox", itob (apstati (ap, MKPSFBOX))) + + # Store the data dependent parameters. + call ap_dapars (ap) + + # Closeup. + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/apppsf.x b/noao/digiphot/apphot/fitpsf/apppsf.x new file mode 100644 index 00000000..c3604c87 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apppsf.x @@ -0,0 +1,123 @@ +include "../lib/apphotdef.h" +include "../lib/fitpsfdef.h" +include "../lib/apphot.h" +include "../lib/fitpsf.h" + +# AP_PPSF -- Procedure to write the results of fitpsf to the output file. + +procedure ap_ppsf (ap, fd, id, lid, ier) + +pointer ap # pointer to apphot structure +int fd # output file descriptor +int id # sequence number of star +int lid # list id of star +int ier # comment string + +real apstatr() + +begin + # Initialize. + if (fd == NULL) + return + + # Print the stars id. + call ap_wid (ap, fd, apstatr (ap, OPFXCUR), apstatr (ap, OPFYCUR), id, + lid, '\\') + + # Print the parameters. + call ap_wfres (ap, fd, ier) +end + + +# AP_QPPSF -- Procedure to print the results of fitpsf on the standard output +# in short form. + +procedure ap_qppsf (ap, ier) + +pointer ap # pointer to apphot structure +int ier # comment string + +pointer sp, imname, psf + +begin + # Initialize. + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + psf = AP_PPSF(ap) + + # Print the parameters on the standard output. + call apstats (ap, IMROOT, Memc[imname], SZ_FNAME) + switch (AP_PSFUNCTION(psf)) { + case AP_RADGAUSS: + call printf ( + "%s %8.2f %8.2f %6.3f %8g %8g %s\n") + call pargstr (Memc[imname]) + call pargr (Memr[AP_PPARS(psf)+1]) + call pargr (Memr[AP_PPARS(psf)+2]) + call pargr (Memr[AP_PPARS(psf)+3]) + call pargr (Memr[AP_PPARS(psf)]) + call pargr (Memr[AP_PPARS(psf)+4]) + if (ier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + case AP_ELLGAUSS: + call printf ("%s %8.2f %8.2f %6.3f %6.3f ") + call pargstr (Memc[imname]) + call pargr (Memr[AP_PPARS(psf)+1]) + call pargr (Memr[AP_PPARS(psf)+2]) + call pargr (Memr[AP_PPARS(psf)+3]) + call pargr (Memr[AP_PPARS(psf)+4]) + call printf ("%6.1f %8g %8g %s\n") + call pargr (Memr[AP_PPARS(psf)+5]) + call pargr (Memr[AP_PPARS(psf)]) + call pargr (Memr[AP_PPARS(psf)+6]) + if (ier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + case AP_MOMENTS: + call printf ("%s %8.2f %8.2f %6.3f %6.3f ") + call pargstr (Memc[imname]) + call pargr (Memr[AP_PPARS(psf)+1]) + call pargr (Memr[AP_PPARS(psf)+2]) + call pargr (Memr[AP_PPARS(psf)+3]) + call pargr (Memr[AP_PPARS(psf)+4]) + call printf ("%6.1f %8g %8g %s\n") + call pargr (Memr[AP_PPARS(psf)+5]) + call pargr (Memr[AP_PPARS(psf)]) + call pargr (Memr[AP_PPARS(psf)+6]) + if (ier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + } + + call sfree (sp) +end + + +# AP_PFHDR -- Procedure to write the banner for the fitpsf output file. + +procedure ap_pfhdr (ap, fd) + +pointer ap # pointer to the apphot strucuture +int fd # output file descriptor + +int apstati() + +begin + if (fd == NULL) + return + + switch (apstati (ap, PSFUNCTION)) { + case AP_RADGAUSS: + call radhdr (ap, fd) + case AP_ELLGAUSS: + call elhdr (ap, fd) + case AP_MOMENTS: + call momhdr (ap, fd) + default: + ; + } +end diff --git a/noao/digiphot/apphot/fitpsf/appsfshow.x b/noao/digiphot/apphot/fitpsf/appsfshow.x new file mode 100644 index 00000000..3317bd22 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/appsfshow.x @@ -0,0 +1,59 @@ +include "../lib/display.h" +include "../lib/fitpsf.h" + +# AP_PSFSHOW -- Procedure to print the fitpsf parameters on the standard output. + +procedure ap_psfshow (ap) + +pointer ap # pointer to the apphot structure + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_pfshow (ap) +end + + +# AP_PFSHOW -- Procedure to print out the fitting parameters on the standard +# output. + +procedure ap_pfshow (ap) + +pointer ap # pointer to apphot structure + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + # Print the PSF fitting parameters. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + call printf ("Psf Modelling Parameters\n") + call apstats (ap, PSFSTRING, Memc[str], SZ_FNAME) + call printf (" %s = %s\n") + call pargstr (KY_PSFSTRING) + call pargstr (Memc[str]) + + call printf (" %s = %g %s %s = %d\n") + call pargstr (KY_PSFAPERT) + call pargr (2.0 * apstatr (ap, PSFAPERT)) + call pargstr (UN_PSFSCALEUNIT) + call pargstr (KY_PMAXITER) + call pargi (apstati (ap, PMAXITER)) + + call printf (" %s = %g %s %s = %d\n") + call pargstr (KY_PK2) + call pargr (apstatr (ap, PK2)) + call pargstr (UN_PSFSIGMA) + call pargstr (KY_PNREJECT) + call pargi (apstati (ap, PNREJECT)) + + call printf (" %s = %b\n") + call pargstr (KY_MKPSFBOX) + call pargb (itob (apstati (ap, MKPSFBOX))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfcolon.x b/noao/digiphot/apphot/fitpsf/apsfcolon.x new file mode 100644 index 00000000..cb4d4521 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfcolon.x @@ -0,0 +1,238 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/fitpsf.h" + +# APSFCOLON -- Process the fitpsf task colon commands. + +procedure apsfcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, newbuf, + newfit) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the iraf image +int cl # coordinate file descriptor +int out # output file descriptor +int stid # output file sequence number +int ltid # coord file sequence number +char cmdstr # command string +int newimage # new image ? +int newbuf # new psf buffer ? +int newfit # new psf fit ? + +int junk +pointer sp, incmd, outcmd +int strdic() + +begin + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, PSFCMDS) != 0) + call ap_fitcolon (ap, out, stid, cmdstr, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0) + call ap_apcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + junk, junk, junk, junk, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0) + call apnscolon (ap, im, out, stid, cmdstr, junk, junk, + junk, junk, newbuf, newfit) + else + call ap_pfimcolon (ap, cmdstr) + + call sfree (sp) +end + + +# AP_FITCOLON -- Procedure to show/set the fitpsf parameters. + +procedure ap_fitcolon (ap, out, stid, cmdstr, newbuf, newfit) + +pointer ap # pointer to the apphot structure +int out # output file descriptor +int stid # output file sequence number +char cmdstr # command string +int newbuf # new psf buffer +int newfit # new psf fit + +bool bval +int ncmd, ival, stat +pointer sp, cmd, str +real rval + +bool itob() +int nscan(), strdic(), apstati(), btoi() +real apstatr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PSFCMDS) + switch (ncmd) { + case PFCMD_FUNCTION: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + call apstats (ap, PSFSTRING, Memc[str], SZ_FNAME) + call printf ("%s = %s\n") + call pargstr (KY_PSFSTRING) + call pargstr (Memc[str]) + } else { + stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PSFFUNCS) + if (stat > 0) { + call apseti (ap, PSFUNCTION, stat) + call apsets (ap, PSFSTRING, Memc[cmd]) + switch (stat) { + case AP_RADGAUSS: + call apseti (ap, NPARS, 5) + case AP_ELLGAUSS: + call apseti (ap, NPARS, 7) + case AP_MOMENTS: + call apseti (ap, NPARS, 7) + } + if (stid > 1) + call ap_sparam (out, "FUNCTION", Memc[cmd], "model", + "fitting function") + newfit = YES + } + } + case PFCMD_BOX: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_PSFAPERT) + call pargr (2.0 * apstatr (ap, PSFAPERT)) + call pargstr (UN_PSFSCALEUNIT) + } else { + call apsetr (ap, PSFAPERT, rval / 2.0) + if (stid > 1) + call ap_rparam (out, KY_PSFAPERT, rval, UN_PSFSCALEUNIT, + "fitting box width") + newbuf = YES + newfit = YES + } + case PFCMD_KREJECT: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_PK2) + call pargr (apstatr (ap, PK2)) + call pargstr (UN_PSFSIGMA) + } else { + call apsetr (ap, PK2, rval) + if (stid > 1) + call ap_rparam (out, KY_PK2, rval, UN_PSFSIGMA, + "k-sigma rejection criterion") + newfit = YES + } + case PFCMD_MAXITER: + call gargi (ival) + if (nscan () == 1) { + call printf ("%s = %d\n") + call pargstr (KY_PMAXITER) + call pargi (apstati (ap, PMAXITER)) + } else { + call apseti (ap, PMAXITER, ival) + if (stid > 1) + call ap_iparam (out, KY_PMAXITER, ival, UN_PSFNUMBER, + "maximum number of iterations") + newfit = YES + } + case PFCMD_NREJECT: + call gargi (ival) + if (nscan () == 1) { + call printf ("%s = %d\n") + call pargstr (KY_PNREJECT) + call pargi (apstati (ap, PNREJECT)) + } else { + call apseti (ap, PNREJECT, ival) + if (stid > 1) + call ap_iparam (out, KY_PNREJECT, ival, UN_PSFNUMBER, + "maximum number of rejection cycles") + newfit = YES + } + case PFCMD_MKBOX: + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (KY_MKPSFBOX) + call pargb (itob (apstati (ap, MKPSFBOX))) + } else + call apseti (ap, MKPSFBOX, btoi (bval)) + default: + # do nothing gracefully + call printf ("Unknown or ambiguous colon command\n") + } + + call sfree (sp) +end + + +# AP_PFIMCOLON -- Procedure to process fitpsf commands that are not fitpsf +# parameters. + +procedure ap_pfimcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr # command string + +int ncmd +pointer sp, cmd +int strdic() + +begin + # Get the command. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, MISC1) + switch (ncmd) { + case ACMD_SHOW: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, PFSHOWARGS) + switch (ncmd) { + case PFCMD_DATA: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + case PFCMD_FIT: + call printf ("\n") + call ap_pfshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_psfshow (ap) + call printf ("\n") + } + default: + call printf ("Unknown or ambigous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfelgauss.x b/noao/digiphot/apphot/fitpsf/apsfelgauss.x new file mode 100644 index 00000000..df1a9ef4 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfelgauss.x @@ -0,0 +1,185 @@ +include <math.h> +include <math/nlfit.h> +include "../lib/noise.h" +include "../lib/fitpsf.h" + +define NPARAMETERS 7 +define TOL 0.001 + +# APSFELGAUSS -- Procedure to fit an elliptical Gaussian function to the +# stellar data. + +int procedure apsfelgauss (ctrpix, nx, ny, emission, fwhmpsf, datamin, + datamax, noise, gain, sigma, maxiter, k2, nreject, par, perr, npar) + +real ctrpix[nx,ny] # object to be centered +int nx, ny # dimensions of subarray +int emission # emission or absorption object +real fwhmpsf # full width half max of the psf +real datamin # minimum good data value +real datamax # maximum good data value +int noise # noise model +real gain # the gain parameter +real sigma # constant term to noise +int maxiter # maximum number of iterations +real k2 # k-sigma rejection criterion +int nreject # maximum number of rejection cycles +real par[ARB] # parameters +real perr[ARB] # errors in parameters +int npar # number of parameters + +extern elgauss, delgauss +int i, j, npts, fier, imin,imax +pointer sp, x, w, list, zfit, nl, ptr +real sumw, dummy, chisqr, locut, hicut, ptemp +int locpr(), apreject() +real asumr(), apwssqr() + +begin + # Initialize. + npts = nx * ny + if (npts < NPARAMETERS) + return (AP_NPSF_TOO_SMALL) + + # Allocate working space. + call smark (sp) + call salloc (x, 2 * npts, TY_REAL) + call salloc (w, npts, TY_REAL) + call salloc (zfit, npts, TY_REAL) + call salloc (list, NPARAMETERS, TY_INT) + + # Define the active parameters. + do i = 1, NPARAMETERS + Memi[list+i-1] = i + + # Set up the varaibles array. + ptr = x + do j = 1, ny { + do i = 1, nx { + Memr[ptr] = i + Memr[ptr+1] = j + ptr = ptr + 2 + } + } + + # Set up the weight array. + switch (noise) { + case AP_NCONSTANT: + call amovkr (1.0, Memr[w], npts) + case AP_NPOISSON: + call amaxkr (ctrpix, 0.0, Memr[w], npts) + if (gain > 0.0) + call adivkr (Memr[w], gain, Memr[w], npts) + if (! IS_INDEFR(sigma)) + call aaddkr (Memr[w], sigma ** 2, Memr[w], npts) + call apreciprocal (Memr[w], Memr[w], npts, 1.0) + default: + call amovkr (1.0, Memr[w], npts) + } + + # Make an initial guess at the fitting parameters. + if (emission == YES) + call ap_wlimr (ctrpix, Memr[w], nx * ny, datamin, datamax, + par[7], par[1], imin, imax) + else + call ap_wlimr (ctrpix, Memr[w], nx * ny, datamin, datamax, + par[7], par[1], imax, imin) + par[1] = par[1] - par[7] + if (mod (imax, nx) == 0) + imin = imax / nx + else + imin = imax / nx + 1 + par[3] = imin + imin = imax - (imin - 1) * nx + par[2] = imin + par[4] = (fwhmpsf ** 2 / 4.0) + par[5] = (fwhmpsf ** 2 / 4.0) + par[6] = 0.0 + + # Get the centers and errors. + call nlinitr (nl, locpr (elgauss), locpr (delgauss), par, perr, + NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, fier) + + # Perform the rejection cycle. + if (nreject > 0 && k2 > 0.0) { + do i = 1, nreject { + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call asubr (ctrpix, Memr[zfit], Memr[zfit], npts) + chisqr = apwssqr (Memr[zfit], Memr[w], npts) + sumw = asumr (Memr[w], npts) + if (sumw <= 0.0) + break + else if (chisqr <= 0.0) + break + else + chisqr = sqrt (chisqr / sumw) + locut = - k2 * chisqr + hicut = k2 * chisqr + if (apreject (Memr[zfit], Memr[w], npts, locut, hicut) == 0) + break + call nlpgetr (nl, par, npar) + call nlfreer (nl) + call nlinitr (nl, locpr (elgauss), locpr (delgauss), par, + perr, NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, + fier) + } + } + + # Fetch the parameters. + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call nlpgetr (nl, par, npar) + par[4] = sqrt (abs(par[4])) + par[5] = sqrt (abs(par[5])) + + # Fetch the errors. + call nlerrorsr (nl, ctrpix, Memr[zfit], Memr[w], npts, dummy, + chisqr, perr) + perr[4] = sqrt (perr[4]) + perr[5] = sqrt (perr[5]) + + # Compute the mean errors. + dummy = 0.0 + do i = 1, npts { + if (Memr[w+i-1] > 0.0) + dummy = dummy + 1.0 + } + dummy = sqrt (dummy) + if (dummy > 0.0) + call adivkr (perr, dummy, perr, npar) + + # Transform the parameters. + par[6] = mod (RADTODEG(par[6]), 360.0) + if (par[6] < 0.0) + par[6] = 360.0 + par[6] + if (par[6] > 90.0 && par[6] <= 270.0) + par[6] = par[6] - 180.0 + else if (par[6] > 270.0 && par[6] <= 360.0) + par[6] = par[6] - 360.0 + if (par[5] > par[4]) { + if (par[6] > 0.0) + par[6] = par[6] - 90.0 + else if (par[6] < 0.0) + par[6] = par[6] + 90.0 + ptemp = par[4] + par[4] = par[5] + par[5] = ptemp + } + perr[6] = mod (RADTODEG(perr[6]), 360.0) + + call nlfreer (nl) + + call sfree (sp) + + # Return the appropriate error code. + if (fier == NO_DEG_FREEDOM) { + return (AP_NPSF_TOO_SMALL) + } else if (fier == SINGULAR) { + return (AP_PSF_SINGULAR) + } else if (fier == NOT_DONE) { + return (AP_PSF_NOCONVERGE) + } else { + return (AP_OK) + } +end diff --git a/noao/digiphot/apphot/fitpsf/apsffit.x b/noao/digiphot/apphot/fitpsf/apsffit.x new file mode 100644 index 00000000..2e6d6dbc --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsffit.x @@ -0,0 +1,136 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/fitpsfdef.h" +include "../lib/noisedef.h" +include "../lib/fitpsf.h" + +# APSFFIT -- Procedure to fit an analytic function to the PSF. + +int procedure apsffit (ap, im, wx, wy) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +real wx, wy # object coordinates + +int ier, fier +pointer psf, nse +real datamin, datamax, dmin, dmax, threshold +int apfbuf(), apsfradgauss(), apsfelgauss(), apsfmoments() + +begin + # Initialize. + psf = AP_PPSF(ap) + nse = AP_NOISE(ap) + AP_PFXCUR(psf) = wx + AP_PFYCUR(psf) = wy + if (IS_INDEFR(wx) || IS_INDEFR(wy)) { + AP_OPFXCUR(psf) = INDEFR + AP_OPFYCUR(psf) = INDEFR + } else { + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_OPFXCUR(psf), AP_OPFYCUR(psf), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_OPFXCUR(psf), AP_OPFYCUR(psf), 1) + default: + AP_OPFXCUR(psf) = wx + AP_OPFYCUR(psf) = wy + } + } + call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_MAXNPARS(psf)) + call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_MAXNPARS(psf)) + + # Fetch the buffer of pixels. + ier = apfbuf (ap, im, wx, wy) + if (ier == AP_NOPSFAREA) + return (AP_NOPSFAREA) + + # Compute the min and max of the data subraster. + if (IS_INDEFR(AP_DATAMIN(ap))) + datamin = -MAX_REAL + else + datamin = AP_DATAMIN(ap) + if (IS_INDEFR(AP_DATAMAX(ap))) + datamax = MAX_REAL + else + datamax = AP_DATAMAX(ap) + + switch (AP_PSFUNCTION(psf)) { + + case AP_RADGAUSS: + + fier = apsfradgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf), + AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin, + datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse), + AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf), + AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf) + + case AP_ELLGAUSS: + + fier = apsfelgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf), + AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin, + datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse), + AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf), + AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf) + + case AP_MOMENTS: + + call alimr (Memr[AP_PSFPIX(psf)], AP_PNX(psf) * AP_PNY(psf), + dmin, dmax) + dmin = max (dmin, datamin) + dmax = min (dmax, datamax) + threshold = 0.0 + + if (AP_POSITIVE(ap) == YES) + fier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf), + AP_PNY(psf), dmin + threshold, dmax, + AP_POSITIVE(ap), Memr[AP_PPARS(psf)], Memr[AP_PPERRS(psf)], + AP_PSFNPARS(psf)) + else + fier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf), + AP_PNY(psf), dmax - threshold, dmin, + AP_POSITIVE(ap), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf) + + default: + + # do nothing gracefully + + } + + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], + Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1) + case WCS_TV: + call ap_ltov (im, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], + Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1) + default: + ; + } + + # Return the appropriate error code. + if (fier == AP_OK) { + if (ier == AP_PSF_OUTOFBOUNDS) + return (AP_PSF_OUTOFBOUNDS) + else + return (AP_OK) + } else if (fier == AP_NPSF_TOO_SMALL) { + call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_PSFNPARS(psf)) + call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + return (AP_NPSF_TOO_SMALL) + } else + return (fier) +end diff --git a/noao/digiphot/apphot/fitpsf/apsffree.x b/noao/digiphot/apphot/fitpsf/apsffree.x new file mode 100644 index 00000000..6bac4607 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsffree.x @@ -0,0 +1,52 @@ +include "../lib/apphotdef.h" +include "../lib/fitpsfdef.h" + +# APSFFREE -- Procedure to free the point spread fitting structure. + +procedure apsffree (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_PPSF(ap) != NULL) + call ap_psfcls (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_PSFCLS -- Procedure to close up the point spread fitting function +# and arrays + +procedure ap_psfcls (ap) + +pointer ap # pointer to apphot structure + +pointer psf + +begin + if (AP_PPSF(ap) == NULL) + return + psf = AP_PPSF(ap) + #if (AP_PSFPIX(psf) != NULL) + #call mfree (AP_PSFPIX(psf), TY_REAL) + if (AP_PSFXPIX(psf) != NULL) + call mfree (AP_PSFXPIX(psf), TY_REAL) + if (AP_PSFYPIX(psf) != NULL) + call mfree (AP_PSFYPIX(psf), TY_REAL) + if (AP_PPARS(psf) != NULL) + call mfree (AP_PPARS(psf), TY_REAL) + if (AP_PPERRS(psf) != NULL) + call mfree (AP_PPERRS(psf), TY_REAL) + if (AP_PPSF(ap) != NULL) + call mfree (AP_PPSF(ap), TY_STRUCT) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfinit.x b/noao/digiphot/apphot/fitpsf/apsfinit.x new file mode 100644 index 00000000..19e5a42d --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfinit.x @@ -0,0 +1,94 @@ +include "../lib/apphotdef.h" +include "../lib/fitpsfdef.h" +include "../lib/fitpsf.h" + +# APSFINIT - Procedure to initialize the point spread modelling structure. + +procedure apsfinit (ap, function, rbox, fwhmpsf, noise) + +pointer ap # pointer to the apphot structure +int function # fitting function +real rbox # fitting radius +real fwhmpsf # full width half max of psf +int noise # noise model + +begin + # Initialize the image parameters. + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set up the global apphot package parameters. + call ap_defsetup (ap, fwhmpsf) + + # Setup noise model. + call ap_noisesetup (ap, noise) + + # Set up the point spread fitting function. + call ap_psfsetup (ap, function, rbox) + + # Set display options. + call ap_dispsetup (ap) + + # Set remaining unused structure pointers to NULL. + AP_PCENTER(ap) = NULL + AP_PSKY(ap) = NULL + AP_PPHOT(ap) = NULL + AP_POLY(ap) = NULL + AP_RPROF(ap) = NULL +end + + +# AP_PSFSETUP -- Procedure to define the PSF fitting parameters. + +procedure ap_psfsetup (ap, function, rbox) + +pointer ap # pointer to apphot structure +int function # fitting function +real rbox # fitting aperture + +pointer psf + +begin + call malloc (AP_PPSF(ap), LEN_PSFSTRUCT, TY_STRUCT) + psf = AP_PPSF(ap) + + # Set PSF fitting function. + AP_PSFUNCTION(psf) = function + switch (function) { + case AP_RADGAUSS: + call strcpy ("radgauss", AP_PSFSTRING(psf), SZ_FNAME) + case AP_ELLGAUSS: + call strcpy ("elgauss", AP_PSFSTRING(psf), SZ_FNAME) + case AP_MOMENTS: + call strcpy ("moments", AP_PSFSTRING(psf), SZ_FNAME) + default: + call strcpy ("radgauss", AP_PSFSTRING(psf), SZ_FNAME) + } + AP_PFXCUR(psf) = INDEFR + AP_PFYCUR(psf) = INDEFR + switch (function) { + case AP_RADGAUSS: + AP_PSFNPARS(psf) = 5 + case AP_ELLGAUSS: + AP_PSFNPARS(psf) = 7 + case AP_MOMENTS: + AP_PSFNPARS(psf) = 7 + } + + # Set remaining PSF parameters. + AP_PSFAPERT(psf) = rbox + AP_MAXNPARS(psf) = DEF_MAXNPARS + AP_PK2(psf) = DEF_PK2 + AP_PMAXITER(psf) = DEF_PMAXITER + AP_PNREJECT(psf) = DEF_PNREJECT + + # Initialize buffers. + AP_LENPSFBUF(psf) = 0 + AP_NPSFPIX(psf) = 0 + AP_PSFPIX(psf) = NULL + AP_PSFXPIX(psf) = NULL + AP_PSFYPIX(psf) = NULL + + # Allocate space for computed parameters. + call calloc (AP_PPARS(psf), AP_MAXNPARS(psf), TY_REAL) + call calloc (AP_PPERRS(psf), AP_MAXNPARS(psf), TY_REAL) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfmoments.x b/noao/digiphot/apphot/fitpsf/apsfmoments.x new file mode 100644 index 00000000..ddae00f2 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfmoments.x @@ -0,0 +1,106 @@ +include <math.h> +include "../lib/fitpsf.h" + +define NPARAMETERS 7 + +# APSFMOMENTS -- Procedure to compute the 0, 1st and second moments of an +# image and estimate the x,y center, the ellipticity and the position angle. + +int procedure apsfmoments (ctrpix, nx, ny, lthreshold, uthreshold, positive, + par, perr, npar) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +real lthreshold # lower threshold for moment computation +real uthreshold # upper threshold for moment computation +int positive # emission feature +real par[ARB] # parameters +real perr[ARB] # errors in parameters +int npar # number of parameters + +int i, j +real temp, sumi, sumxi, sumyi, sumx2i, sumy2i, sumxyi, r2, varx, vary, varxy +bool fp_equalr() + +begin + # Initialize the sums. + sumi = 0.0 + sumxi = 0.0 + sumyi = 0.0 + sumxyi = 0.0 + sumx2i = 0.0 + sumy2i = 0.0 + + # Accumulate the moments. + if (positive == YES) { + do j = 1, ny { + do i = 1, nx { + if (ctrpix[i,j] > uthreshold) + next + temp = ctrpix[i,j] - lthreshold + if (temp <= 0.0) + next + sumi = sumi + temp + sumxi = sumxi + i * temp + sumyi = sumyi + j * temp + sumxyi = sumxyi + i * j * temp + sumx2i = sumx2i + i * i * temp + sumy2i = sumy2i + j * j * temp + } + } + } else { + do j = 1, ny { + do i = 1, nx { + if (ctrpix[i,j] < uthreshold) + next + temp = lthreshold - ctrpix[i,j] + if (temp <= 0.0) + next + sumi = sumi + temp + sumxi = sumxi + i * temp + sumyi = sumyi + j * temp + sumxyi = sumxyi + i * j * temp + sumx2i = sumx2i + i * i * temp + sumy2i = sumy2i + j * j * temp + } + } + } + + # Compute the parameters. + if (fp_equalr (sumi, 0.0)) { + par[1] = 0.0 + par[2] = (1 + nx) / 2.0 + par[3] = (1 + ny) / 2.0 + par[4] = 0.0 + par[5] = 0.0 + par[6] = 0.0 + par[7] = lthreshold + } else { + par[1] = sumi + par[2] = sumxi / sumi + par[3] = sumyi / sumi + varx = max (0.0, sumx2i / sumi - par[2] ** 2) + vary = max (0.0, sumy2i / sumi - par[3] ** 2) + r2 = varx + vary + if (r2 <= 0.0) { + par[4] = 0.0 + par[5] = 0.0 + par[6] = 0.0 + } else { + par[4] = sqrt (r2) + varxy = sumxyi / sumi - par[2] * par[3] + par[5] = sqrt ((varx - vary) ** 2 + 4.0 * varxy ** 2) / r2 + par[6] = (0.5 * r2 * (1.0 + par[5]) - vary) + if (fp_equalr (par[6], 0.0)) + par[6] = 90.0 + else + par[6] = RADTODEG (atan (varxy / par[6])) + } + par[7] = lthreshold + } + + # Compute the errors. + npar = NPARAMETERS + call amovkr (0.0, perr, NPARAMETERS) + return (AP_OK) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfradgauss.x b/noao/digiphot/apphot/fitpsf/apsfradgauss.x new file mode 100644 index 00000000..0b4a93c6 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfradgauss.x @@ -0,0 +1,183 @@ +include <math/nlfit.h> +include "../lib/noise.h" +include "../lib/fitpsf.h" + +define NPARAMETERS 5 +define TOL 0.001 + +# APSFRADGAUSS -- Fit a radial Gaussian function to the data. + +int procedure apsfradgauss (ctrpix, nx, ny, emission, fwhmpsf, datamin, + datamax, noise, gain, sigma, maxiter, k2, nreject, par, perr, npar) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +int emission # emission or absorption version +real fwhmpsf # full width half max of the psf +real datamin # minimum good data value +real datamax # maximum good data value +int noise # noise model to be used +real gain # the gain in the data +real sigma # sigma of constant noise term +int maxiter # maximum number of iterations +real k2 # k-sigma rejection criterion +int nreject # maximum number of rejection cycles +real par[ARB] # parameters +real perr[ARB] # errors in parameters +int npar # number of parameters + +extern gaussr, dgaussr +int i, j, npts, list, imin, imax, fier +pointer sp, x, w, zfit, nl, ptr +real sumw, dummy, chisqr, locut, hicut +int locpr(), apreject() +real asumr(), apwssqr() + +begin + # Initialize. + npts = nx * ny + if (npts < NPARAMETERS) + return (AP_NPSF_TOO_SMALL) + + call smark (sp) + call salloc (x, 2 * npts, TY_REAL) + call salloc (w, npts, TY_REAL) + call salloc (zfit, npts, TY_REAL) + call salloc (list, NPARAMETERS, TY_INT) + + # Define the active parameters. + do i = 1, NPARAMETERS + Memi[list+i-1] = i + + # Set variables array. + ptr = x + do j = 1, ny { + do i = 1, nx { + Memr[ptr] = i + Memr[ptr+1] = j + ptr = ptr + 2 + } + } + + # Define the weight array. + switch (noise) { + case AP_NCONSTANT: + call amovkr (1.0, Memr[w], npts) + case AP_NPOISSON: + call amaxkr (ctrpix, 0.0, Memr[w], npts) + if (gain > 0.0) + call adivkr (Memr[w], gain, Memr[w], npts) + if (! IS_INDEFR(sigma)) + call aaddkr (Memr[w], sigma ** 2, Memr[w], npts) + call apreciprocal (Memr[w], Memr[w], npts, 1.0) + default: + call amovkr (1.0, Memr[w], npts) + } + + # Make an initial guess at the parameters. + if (emission == YES) + call ap_wlimr (ctrpix, Memr[w], npts, datamin, datamax, + par[5], par[1], imin, imax) + else + call ap_wlimr (ctrpix, Memr[w], npts, datamin, datamax, + par[1], par[5], imax, imin) + par[1] = par[1] - par[5] + if (mod (imax, nx) == 0) + imin = imax / nx + else + imin = imax / nx + 1 + par[3] = imin + imin = imax - (imin - 1) * nx + par[2] = imin + par[4] = (fwhmpsf ** 2 / 4.0) + + # Fit the function and the errors. + call nlinitr (nl, locpr (gaussr), locpr (dgaussr), par, perr, + NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, fier) + + # Perform the rejection cycle. + if ((nreject > 0) && (k2 > 0.0)) { + do i = 1, nreject { + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call asubr (ctrpix, Memr[zfit], Memr[zfit], npts) + chisqr = apwssqr (Memr[zfit], Memr[w], npts) + sumw = asumr (Memr[w], npts) + if (sumw <= 0.0) + break + else if (chisqr <= 0.0) + break + else + chisqr = sqrt (chisqr / sumw) + locut = - k2 * chisqr + hicut = k2 * chisqr + if (apreject (Memr[zfit], Memr[w], npts, locut, hicut) == 0) + break + call nlpgetr (nl, par, npar) + call nlfreer (nl) + call nlinitr (nl, locpr (gaussr), locpr (dgaussr), par, perr, + NPARAMETERS, Memi[list], NPARAMETERS, TOL, maxiter) + call nlfitr (nl, Memr[x], ctrpix, Memr[w], npts, 2, WTS_USER, + fier) + } + } + + # Get the parameters. + call nlpgetr (nl, par, npar) + par[4] = sqrt (abs(par[4])) + + # Get the errors. + call nlvectorr (nl, Memr[x], Memr[zfit], npts, 2) + call nlerrorsr (nl, ctrpix, Memr[zfit], Memr[w], npts, dummy, + chisqr, perr) + perr[4] = sqrt (abs(perr[4])) + + # Compute the mean errors in the parameters. + dummy = 0.0 + do i = 1, npts { + if (Memr[w+i-1] > 0.0) + dummy = dummy + 1.0 + } + dummy = sqrt (dummy) + if (dummy > 0.0) + call adivkr (perr, dummy, perr, npar) + + call nlfreer (nl) + + call sfree (sp) + + # Return the appropriate error code. + if (fier == NO_DEG_FREEDOM) { + return (AP_NPSF_TOO_SMALL) + } else if (fier == SINGULAR) { + return (AP_PSF_SINGULAR) + } else if (fier == NOT_DONE) { + return (AP_PSF_NOCONVERGE) + } else { + return (AP_OK) + } +end + + +# APREJECT -- Reject points outside of the specified intensity limits by +# setting their weights to zero. + +int procedure apreject (pix, w, npts, locut, hicut) + +real pix[ARB] # data +real w[ARB] # weights +int npts # number of data points +real locut, hicut # data limits + +int i, nreject + +begin + nreject = 0 + do i = 1, npts { + if ((pix[i] < locut || pix[i] > hicut) && w[i] > 0.0) { + w[i] = 0.0 + nreject = nreject + 1 + } + } + return (nreject) +end diff --git a/noao/digiphot/apphot/fitpsf/apsfrefit.x b/noao/digiphot/apphot/fitpsf/apsfrefit.x new file mode 100644 index 00000000..a1b950fa --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/apsfrefit.x @@ -0,0 +1,113 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/fitpsfdef.h" +include "../lib/fitpsf.h" +include "../lib/noisedef.h" + +# APSFREFIT -- Procedure to fit a function to the data already in the +# data buffers. + +int procedure apsfrefit (ap, im) + +pointer ap # pointer to the apphot structure +pointer im # the input image descriptor + +int ier +pointer psf, nse +int apsfradgauss(), apsfelgauss(), apsfmoments() +real datamin, datamax, dmin, dmax, threshold + +begin + psf = AP_PPSF(ap) + nse = AP_NOISE(ap) + call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_MAXNPARS(psf)) + call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_MAXNPARS(psf)) + + # Compute the minimum and maximum good data value. + if (IS_INDEFR(AP_DATAMIN(ap))) + datamin = -MAX_REAL + else + datamin = AP_DATAMIN(ap) + if (IS_INDEFR(AP_DATAMAX(ap))) + datamax = MAX_REAL + else + datamax = AP_DATAMAX(ap) + + switch (AP_PSFUNCTION(psf)) { + + case AP_RADGAUSS: + + ier = apsfradgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf), + AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin, + datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse), + AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf), + AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) - + AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) - + AP_PYC(psf) + + case AP_ELLGAUSS: + + ier = apsfelgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf), + AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin, + datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse), + AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf), + AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) - + AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) - + AP_PYC(psf) + + case AP_MOMENTS: + + call alimr (Memr[AP_PSFPIX(psf)], AP_PNX(psf) * AP_PNY(psf), + dmin, dmax) + dmin = max (dmin, datamin) + dmax = min (dmax, datamax) + threshold = 0.0 + + if (AP_POSITIVE(ap) == YES) + ier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf), + AP_PNY(psf), dmin + threshold, dmax, + AP_POSITIVE(ap), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + else + ier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf), + AP_PNY(psf), dmax + threshold, dmin, + AP_POSITIVE(ap), Memr[AP_PPARS(psf)], + Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + + Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + AP_PFXCUR(psf) - + AP_PXC(psf) + Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + AP_PFYCUR(psf) - + AP_PYC(psf) + + default: + + # do nothing gracefully + } + + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], + Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1) + case WCS_TV: + call ap_ltov (im, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], + Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1) + default: + ; + } + + if (ier == AP_NPSF_TOO_SMALL) { + call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_PSFNPARS(psf)) + call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf)) + } + + return (ier) +end diff --git a/noao/digiphot/apphot/fitpsf/fitpsf.key b/noao/digiphot/apphot/fitpsf/fitpsf.key new file mode 100644 index 00000000..f7ab725e --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/fitpsf.key @@ -0,0 +1,73 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Save the current parameters +d Plot radial profile of current star +i Interactively set parameters using current star +f Fit current star +spbar Fit current star, output results +m Move to next star in coordinate list +n Fit next star in coordinate list, output results +l Fit remaining stars in coordinate list, output results +e Print error messages +r Rewind the coordinate list +q Exit task + + + + Colon Commands + +:show [data/fit] List the parameters +:m [n] Move to next [nth] star in coordinate list +:n [n] Fit next [nth] star in coordinate list, output results + + + Colon Parameter Editing Commands + +# Image and file name parameters + +:image [string] Image name +:coords [string] Coordinate file name +:output [string] Output file name + +# Data dependent parameters + +:scale [value] Image scale (units per pixel) +:fwhmpsf [value] Scale factor (scale units) +:emission [y/n] Emission feature (y), absorption (n) +:sigma [value] Standard deviation of sky (counts) +:datamin [value] Minimum good data value (counts) +:datamax [value] Maximum good data value (counts) + +# Noise description parameters + +:noise [string] Noise model (constant|poisson) +:gain [string] Gain image header keyword +:ccdread [string] Readout noise image header keyword +:epadu [value] Gain (electrons per adu) +:readnoise [value] Readnoise (electrons) + +# Observation 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] Exposure time (time units) +:xairmass [value] Airmass value (number) +:ifilter [string] Filter id string +:otime [string] Time of observation (time units) + +# Fitting parameters + +:function [string] PSF model (radgauss|elgauss|moments) +:box [value] Width of the fitting box (scale units) +:maxiter [value] Maximum number of iterations +:nreject [value] Maximum number of rejection cycles +:kreject [value] Rejection limit (sigma) + +# Plotting and marking functions + +:mkbox [y/n] Mark the fitting box on the display diff --git a/noao/digiphot/apphot/fitpsf/ifitpsf.key b/noao/digiphot/apphot/fitpsf/ifitpsf.key new file mode 100644 index 00000000..9bbb4ba9 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/ifitpsf.key @@ -0,0 +1,10 @@ + Interactive Fitpsf Setup Menu + + v Mark and verify the critical fitpsf parameters (f,s,b) + + f Mark and verify the full-width half-maximum of the psf + s Mark and verify the standard deviation of the background + l Mark and verify the minimum good data value + u Mark and verify the maximum good data value + + b Mark and verify the half-width of the fitting box diff --git a/noao/digiphot/apphot/fitpsf/mkpkg b/noao/digiphot/apphot/fitpsf/mkpkg new file mode 100644 index 00000000..8e8efb12 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/mkpkg @@ -0,0 +1,44 @@ +# FITPSF Task Routines + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + apbfitpsf.x <fset.h> ../lib/apphot.h \ + ../lib/display.h + apfbuf.x <imhdr.h> ../lib/apphotdef.h \ + ../lib/fitpsfdef.h ../lib/fitpsf.h + apfitpsf.x <ctype.h> <gset.h> \ + ../lib/apphot.h ../lib/display.h \ + ../lib/fitpsf.h <imhdr.h> + apppfpars.x ../lib/display.h ../lib/fitpsf.h + apgpfpars.x ../lib/display.h ../lib/fitpsf.h \ + ../lib/noise.h + appfconfirm.x ../lib/apphot.h ../lib/noise.h \ + ../lib/fitpsf.h + appferrors.x ../lib/fitpsf.h + apppsf.x ../lib/apphotdef.h ../lib/fitpsfdef.h \ + ../lib/fitpsf.h ../lib/apphot.h + appsfshow.x ../lib/display.h ../lib/fitpsf.h + appfradsetup.x + apsfcolon.x ../lib/display.h ../lib/fitpsf.h \ + ../lib/apphot.h ../lib/noise.h + apsfelgauss.x <math.h> ../lib/fitpsf.h \ + <math/nlfit.h> ../lib/noise.h + apsffree.x ../lib/apphotdef.h ../lib/fitpsfdef.h + apsfinit.x ../lib/apphotdef.h ../lib/fitpsfdef.h \ + ../lib/fitpsf.h + apsffit.x ../lib/apphotdef.h ../lib/fitpsfdef.h \ + ../lib/noisedef.h ../lib/fitpsf.h \ + ../lib/apphot.h <mach.h> + apsfmoments.x <math.h> ../lib/fitpsf.h + apsfradgauss.x ../lib/fitpsf.h <math/nlfit.h> \ + ../lib/noise.h + apsfrefit.x ../lib/apphotdef.h ../lib/fitpsfdef.h \ + ../lib/fitpsf.h ../lib/noisedef.h \ + ../lib/apphot.h <mach.h> + t_fitpsf.x <fset.h> <gset.h> \ + ../lib/apphot.h <imhdr.h> + ; diff --git a/noao/digiphot/apphot/fitpsf/t_fitpsf.x b/noao/digiphot/apphot/fitpsf/t_fitpsf.x new file mode 100644 index 00000000..28345a09 --- /dev/null +++ b/noao/digiphot/apphot/fitpsf/t_fitpsf.x @@ -0,0 +1,288 @@ +include <gset.h> +include <fset.h> +include <imhdr.h> +include "../lib/apphot.h" + +# T_FITPSF -- Procedure to fit an analytic function to the PSF for a list +# of objects in a list of images. + +procedure t_fitpsf () + +pointer image # pointer to the name of the image +pointer output # pointer to the output file name +pointer coords # pointer to the coordinate file +pointer graphics # pointer to the graphics display device +pointer display # pointer to the display device +int interactive # mode of use +int cache # cache the input image pixels +int verify # verify critical parameters +int update # update the critical parameter +int verbose # verbose mode + +pointer sp, outfname, ap, im, gd, id, cname, str +int cl, out, limlist, lclist, lolist, lid, sid, root, stat, memstat +int imlist, clist, olist, wcs, req_size, buf_size, old_size + +pointer gopen(), immap() +int imtlen(), imtgetim(), clplen(), btoi(), clgfil(), fnldir() +int open(), strncmp(), strlen(), apfitpsf(), imtopenp(), clpopnu() +int clgwrd(), ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate workin 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 (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 standard output to flush on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get input image list. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + + # Get input coordinate lists. + clist = clpopnu ("coords") + lclist = clplen (clist) + + # Get output file list and check for zero length list. + olist = clpopnu ("output") + lolist = clplen (olist) + + # Check that image and coordinate list lengths match. + if (limlist < 1 || (lclist > 1 && lclist != limlist)) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatable image and coordinate list lengths") + } + + # Check that image and output list lengths match. + if (lolist > 1 && lolist != limlist) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatable image and output list lengths") + } + + call clgstr ("icommands.p_filename", Memc[cname], SZ_FNAME) + if (Memc[cname] != EOS) + interactive = NO + #else if (lclist == 0) + #interactive = YES + else + interactive = btoi (clgetb ("interactive")) + cache = btoi (clgetb("cache")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + verbose = btoi (clgetb ("verbose")) + + # Get the parameters. + call ap_gpfpars (ap) + if (verify == YES && interactive == NO) { + call ap_pfconfirm (ap, NULL, 1) + if (update == YES) + call ap_ppfpars (ap) + } + + # Get the wcs information. + wcs = clgwrd ("wcsin", Memc[str], SZ_LINE, WCSINSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the input coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call apseti (ap, WCSIN, wcs) + wcs = clgwrd ("wcsout", Memc[str], SZ_LINE, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call apseti (ap, WCSOUT, wcs) + + # Get the graphics and display devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + + # Open the graphics and display devices. + if (interactive == YES) { + if (Memc[graphics] == EOS) + gd = NULL + else { + iferr { + gd = gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH) + } then { + call eprintf ( + "Warning: Error opening graphics device.\n") + gd = NULL + } + } + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) + id = gd + else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + } else { + id = NULL + gd = NULL + } + + # Begin looping over the image list. + sid = 1 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open the image. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (ap, im, Memc[image]) + + # Set the image display viewport. + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + + # Open coordinate file, where coords is assumed to be a simple text + # file in which the x and y positions are in columns 1 and 2 + # respectively and all remaining fields are ignored. + + if (lclist <= 0) { + cl = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (clist, Memc[coords], SZ_FNAME) + root = fnldir (Memc[coords], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[coords+root], 7) == 0 || root == + strlen (Memc[coords])) { + call ap_inname (Memc[image], Memc[outfname], "coo", + Memc[outfname], SZ_FNAME) + lclist = limlist + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else if (stat != EOF) { + call strcpy (Memc[coords], Memc[outfname], SZ_FNAME) + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else { + call apstats (ap, CLNAME, Memc[outfname], SZ_FNAME) + call seek (cl, BOF) + } + } + call apsets (ap, CLNAME, Memc[outfname]) + call apfroot (Memc[outfname], Memc[str], SZ_LINE) + call apsets (ap, CLROOT, Memc[str]) + + # Open output text file, if output is "default", dir$default or + # a directory specification then the extension "psf" is added on + # to the image name and a suitable version number is appended to + # the output name. If the output string is null then no output + # file is written. + + 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], "psf", + Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + lolist = limlist + } else if (stat != EOF) { + call strcpy (Memc[output], Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + } else + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + } + call apsets (ap, OUTNAME, Memc[outfname]) + + # Fit the PSF. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = apfitpsf (ap, im, cl, NULL, NULL, out, sid, NO, + cache) + else if (cl != NULL) { + lid = 1 + call apbfitpsf (ap, im, cl, out, id, sid, lid, verbose) + stat = NO + } else + stat = NO + } else + stat = apfitpsf (ap, im, cl, gd, id, out, sid, YES, cache) + + # Cleanup. + call imunmap (im) + if (cl != NULL) { + if (lclist > 1) + call close (cl) + } + if (out != NULL && lolist != 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + sid = 1 + } + + # Uncache memory. + call fixmem (old_size) + + if (stat == YES) + break + } + + # Close the plot files. + if (id == gd && id != NULL) + call gclose (id) + else { + if (id != NULL) + call gclose (id) + if (gd != NULL) + call gclose (gd) + } + + # If only one coordinate file for a list of images close file. + if (cl != NULL && lclist == 1) + call close (cl) + + # If only one output file for a list of images close file. + if (out != NULL && lolist == 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + } + + # Close up the PSF fitting structure. + call apsffree (ap) + + # Close up the lists. + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + + call sfree (sp) +end |