Initial commit

1 files changed, 329 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/psf/t_pstselect.x b/noao/digiphot/daophot/psf/t_pstselect.x
new file mode 100644
index 00000000..497ab374
--- /dev/null
+++ b/noao/digiphot/daophot/psf/t_pstselect.x
@@ -0,0 +1,329 @@
+include	<fset.h>
+include <imhdr.h>
+include <gset.h>
+include "../lib/daophotdef.h"
+include "../lib/apseldef.h"
+include "../lib/psfdef.h"
+
+# T_PSTSELECT -- Select good candidate PSF stars from a DAOPHOT photometry
+# file.
+
+procedure t_pstselect ()
+
+pointer	image			# input image
+pointer	photfile		# input rough photometry
+pointer	pstfile			# output PSTFILE table
+int	maxnpsf			# maximimum number of psf stars
+pointer	plotfile		# pointer to the plot metacode file
+bool	interactive		# interactive mode
+pointer	plottype		# default plot type
+int	cache			# cache the input image pixels
+int	verify			# verify the critical parameters
+int	update			# update the critical parameters
+pointer	graphics		# the graphics device
+pointer	display			# the display device
+bool	mkstars			# mark deleted and accepted psf stars
+
+pointer	sp, pfd, dao, outfname, curfile, str, im, gd, id, mgd
+int	imlist, limlist, alist, lalist, olist, lolist, root, apd, pmgd, pltype
+int	wcs, req_size, old_size, buf_size, memstat()
+bool	ap_text
+
+pointer	immap(), gopen()
+int	tbtopn(), open(), fnldir(), strlen(), strncmp(), fstati(), btoi()
+int	access(), fntopnb(), fntlenb(), clgeti(), imtopen(), imtlen()
+int	fntgfnb(), imtgetim(), strdic(), dp_stati(), clgwrd(), sizeof()
+int	dp_memstat()
+bool	clgetb(), itob(), streq()
+
+errchk	gopen()
+
+begin
+	# Set the standard output to flush on newline.
+	if (fstati (STDOUT, F_REDIR) == NO)
+	    call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Get some working memory.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (photfile, SZ_FNAME, TY_CHAR)
+	call salloc (pstfile, SZ_FNAME, TY_CHAR)
+	call salloc (plotfile, SZ_FNAME, TY_CHAR)
+	call salloc (plottype, 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 (curfile, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the various task parameters.
+	call clgstr ("image", Memc[image], SZ_FNAME)
+	call clgstr ("photfile", Memc[photfile], SZ_FNAME)
+	call clgstr ("pstfile", Memc[pstfile], SZ_FNAME)
+	maxnpsf = clgeti ("maxnpsf")
+	cache = btoi (clgetb ("cache"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+
+	# Get the lists.
+	imlist = imtopen (Memc[image])
+	limlist = imtlen (imlist)
+	alist = fntopnb (Memc[photfile], NO)
+	lalist = fntlenb (alist)
+	olist = fntopnb (Memc[pstfile], NO)
+	lolist =  fntlenb (olist)
+
+	# Test that the lengths of the photometry file and psf star file
+	# lists are the same as the input image list.
+
+	if ((limlist != lalist) && (strncmp (Memc[photfile], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call fntclsb (olist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and photometry file list lengths\n")
+	}
+
+	if ((limlist != lolist) && (strncmp (Memc[pstfile], DEF_DEFNAME,
+	    DEF_LENDEFNAME) != 0)) {
+	    call imtclose (imlist)
+	    call fntclsb (alist)
+	    call fntclsb (olist)
+	    call sfree (sp)
+	    call error (0,
+	        "Incompatable image and photometry file list lengths\n")
+	}
+
+	# Initialize the DAOPHOT structure, and get the pset parameters.
+	call dp_gppars (dao)	
+
+	# Confirm the parameters.
+	if (verify == YES) {
+	    call dp_ptconfirm (dao)
+	    if (update == YES)
+		call dp_pppars (dao)
+	}
+
+        # 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 dp_seti (dao, 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 dp_seti (dao, WCSOUT, wcs)
+
+	# Initialize the photometry structure.
+	call dp_apsetup (dao)
+
+	# Initialize the PSF fitting structure.
+	call dp_psfsetup (dao)
+
+	# Initialize the PSF structure.
+	call dp_fitsetup (dao)
+
+	# Is the task interactive or not?
+	call clgstr ("icommands.p_filename", Memc[curfile], SZ_FNAME)
+	if (Memc[curfile] == EOS)
+	    interactive = clgetb ("interactive")
+	else
+	    interactive = false
+
+	# Get the graphics display and plot file devices.
+	call clgstr ("graphics", Memc[graphics], SZ_FNAME)
+	call clgstr ("display", Memc[display], SZ_FNAME)
+	call clgstr ("plottype", Memc[plottype], SZ_FNAME)
+	call clgstr ("plotfile", Memc[plotfile], SZ_FNAME)
+
+	# Open graphics and display devices if appropriate.
+	if (interactive) {
+	    call dp_seti (dao, VERBOSE, 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
+		}
+	    }
+	    if (id != NULL)
+	        mkstars = clgetb ("mkstars")
+	    else
+		mkstars = false
+	} else {
+	    gd = NULL
+	    id = NULL
+	    call dp_seti (dao, VERBOSE, btoi (clgetb ("verbose")))
+	    mkstars = false
+	}
+
+	# Open the plot file.
+	if (Memc[plotfile] == EOS)
+	    pmgd = NULL
+	else
+	    pmgd = open (Memc[plotfile], APPEND, BINARY_FILE)
+	if (pmgd != NULL)
+	    mgd = gopen (Memc[graphics], NEW_FILE, pmgd)
+	else
+	    mgd = NULL
+
+	# Set the default plot type.
+	pltype = strdic (Memc[plottype], Memc[plottype], SZ_FNAME, PSF_PLOTS)
+	call dp_pseti (dao, PLOTTYPE, pltype)
+
+	# Loop over the list of input files
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open the input image.
+	    im = immap (Memc[image], READ_ONLY, 0)
+	    call dp_imkeys (dao, im)
+	    call dp_sets (dao, INIMAGE, Memc[image])
+
+            # Set up the display coordinate system.
+            if ((id != NULL) && (id != gd))
+                call dp_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 = dp_memstat (cache, req_size, old_size)
+            if (memstat == YES)
+                call dp_pcache (im, INDEFI, buf_size)
+
+	    # Open the input photometry table and read in the photometry.
+	    if (fntgfnb (alist, Memc[photfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[photfile], SZ_FNAME)
+	    root = fnldir (Memc[photfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[photfile+root], DEF_LENDEFNAME) ==
+	        0 || root == strlen (Memc[photfile]))
+	        call dp_inname (Memc[image], Memc[outfname], "mag",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[photfile], Memc[outfname], SZ_FNAME)
+	    ap_text =  itob (access (Memc[outfname], 0, TEXT_FILE))
+	    if (ap_text)
+	        apd = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+	    else
+	        apd = tbtopn (Memc[outfname], READ_ONLY, 0)
+	    call dp_wgetapert (dao, im, apd, DP_MAXNSTAR(dao), ap_text)
+	    call dp_sets (dao, INPHOTFILE, Memc[outfname])
+
+	    # Open the output PSTSELECT file. If the output is "default",
+	    # dir$default or a directory specification then the extension .pst
+	    # is added to the image name and a suitable version number is
+	    # appended to the output name.
+
+	    if (fntgfnb (olist, Memc[pstfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[pstfile], SZ_FNAME)
+	    root = fnldir (Memc[pstfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[pstfile + root], DEF_LENDEFNAME) ==
+	        0 || root == strlen (Memc[pstfile]))
+	        call dp_outname (Memc[image], Memc[outfname], "pst",
+		    Memc[outfname], SZ_FNAME)
+	    else
+	        call strcpy (Memc[pstfile], Memc[outfname], SZ_FNAME)
+	    if (ap_text)
+	        pfd = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        pfd = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTPHOTFILE, Memc[outfname])
+
+	    if (DP_VERBOSE(dao) == YES) {
+		call printf ("\nSelecting PSF stars for image %s\n")
+		    call pargstr (Memc[image])
+	        call dp_stats (dao, INPHOTFILE, Memc[outfname], SZ_FNAME)
+		call printf ("\t%d stars read from file %s\n\n")
+		    call pargi (dp_stati (dao, APNUM))
+		    call pargstr (Memc[outfname])
+	    }
+
+	    # Now select the PSF stars.
+	    if (Memc[curfile] != EOS)
+	        call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, NULL,
+		    mgd, NULL, false, false, true)
+	    else if (interactive)
+	        call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, gd, mgd,
+		    id, mkstars, true, false)
+	    else
+	        call dp_gpstars (dao, im, apd, pfd, ap_text, maxnpsf, NULL,
+		    mgd, NULL, false, false, false)
+
+	    # Close the input image.
+	    call imunmap (im)
+
+	    # Close the photometry file. 
+	    if (ap_text)
+	        call close (apd)
+	    else
+	        call tbtclo (apd)
+
+	    # Close the output table.
+	    if (ap_text)
+	        call close (pfd)
+	    else
+	        call tbtclo (pfd)
+
+            # Uncache memory.
+            call fixmem (old_size)
+
+	}
+
+	# Close up the graphics and display streams.
+	if (id == gd && id != NULL)
+	    call gclose (id)
+	else {
+	    if (gd != NULL)
+		call gclose (gd)
+	    if (id != NULL)
+		call gclose (id)
+	}
+
+	# Close the metacode plot files.
+	if (mgd != NULL)
+	    call gclose (mgd)
+	if (pmgd != NULL)
+	    call close (pmgd)
+
+	# Close the image/file lists.
+	call imtclose (imlist)
+	call fntclsb (alist)
+	call fntclsb (olist)
+
+	# Close the PSF structure.
+	call dp_fitclose (dao)
+
+	# Close the PSF fitting structure.
+	call dp_psfclose (dao)
+
+	# Free the photometry structure.
+	call dp_apclose (dao)
+	
+	# Free the daophot structure.
+	call dp_free (dao)
+
+	call sfree(sp)
+end