Initial commit

1 files changed, 509 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/psf/t_psf.x b/noao/digiphot/daophot/psf/t_psf.x
new file mode 100644
index 00000000..b7444b3e
--- /dev/null
+++ b/noao/digiphot/daophot/psf/t_psf.x
@@ -0,0 +1,509 @@
+include	<fset.h>
+include <imhdr.h>
+include <gset.h>
+include	"../lib/daophotdef.h"
+include	"../lib/apseldef.h"
+include	"../lib/psfdef.h"
+
+# T_PSF  -- Generate a point spread function from one or more stars in the
+# image frame.
+
+procedure t_psf ()
+
+pointer	image				# the input image
+pointer	photfile			# input aperture photometry file
+pointer	pstarfile			# input psf star file
+pointer	psfimage			# output psf image
+pointer	groupfile			# output psf group file
+pointer	opstfile			# output psf star file
+pointer	graphics			# pointer to graphics device name
+pointer	plotfile			# pointer to plotfile name
+int	cache				# cache the input image pixels
+pointer	display				# pointer to display device name
+bool	matchbyid			# match psf stars by id or position
+bool	interactive			# the mode of task operation
+bool	showplots			# display plots of the psf stars
+pointer	plottype			# type of psf plot
+bool	mkstars				# mark deleted and accepted psf stars
+
+pointer	sp, im, apd, psfim, dao, mgd, gd, id
+pointer	outfname, curfile, str
+int	imlist, limlist, alist, lalist, clist, lclist, pimlist, lpimlist
+int	olist, lolist, oclist, loclist, up, verify, update, wcs
+int	root, min_lenuserarea, pltype, pfd, pst, psfgr, opst
+int	req_size, old_size, buf_size, memstat
+bool	ap_text, pst_text
+
+pointer	immap(), tbtopn(), gopen()
+int	fnldir(), strlen(), strncmp(), btoi(), envfind(), ctoi(), clgwrd()
+int	strdic(), open(), access(), fstati(), dp_stati(), dp_pstati()
+int	imtopen(), imtlen(), imtgetim(), fntopnb(), fntlenb(), fntgfnb()
+int	sizeof(), dp_memstat()
+bool	streq(), clgetb(), itob(), dp_updatepsf()
+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 (pstarfile, SZ_FNAME, TY_CHAR)
+	call salloc (psfimage, SZ_FNAME, TY_CHAR)
+	call salloc (groupfile, SZ_FNAME, TY_CHAR)
+	call salloc (opstfile, SZ_FNAME, TY_CHAR)
+	call salloc (plottype, SZ_FNAME, TY_CHAR)
+	call salloc (outfname, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+	call salloc (display, SZ_FNAME, TY_CHAR)
+	call salloc (plotfile, 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[pstarfile], SZ_FNAME)
+	call clgstr ("psfimage", Memc[psfimage], SZ_FNAME)
+	call clgstr ("opstfile", Memc[opstfile], SZ_FNAME)
+	call clgstr ("groupfile", Memc[groupfile], SZ_FNAME)
+	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)
+	clist = fntopnb (Memc[pstarfile], NO)
+	lclist = fntlenb (clist)
+	pimlist = imtopen (Memc[psfimage])
+	lpimlist = imtlen (pimlist)
+	olist = fntopnb (Memc[groupfile], NO)
+	lolist =  fntlenb (olist)
+	oclist = fntopnb (Memc[opstfile], NO)
+	loclist = fntlenb (oclist)
+
+	# Test that the lengths of the photometry file, psf image, and
+        # output file lists are the same as the length of the input image
+        # list.
+
+	# Compare the image and photometry file list lengths.
+        if ((limlist != lalist) && (strncmp (Memc[photfile],
+            DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+            call imtclose (imlist)
+            call fntclsb (alist)
+            call fntclsb (clist)
+            call imtclose (pimlist)
+            call fntclsb (olist)
+            call fntclsb (oclist)
+            call sfree (sp)
+            call error (0,
+                "Incompatable image and photometry file list lengths")
+        }
+
+	# Compare the image and psf star list lengths.
+        if ((lclist != 0) && (limlist != lclist) && (strncmp (Memc[pstarfile],
+            DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+            call imtclose (imlist)
+            call fntclsb (alist)
+            call fntclsb (clist)
+            call imtclose (pimlist)
+            call fntclsb (olist)
+            call fntclsb (oclist)
+            call sfree (sp)
+            call error (0,
+                "Incompatable image and psf star file list lengths")
+        }
+
+	# Compare the image and psf image list lengths.
+        if ((limlist != lpimlist) && (strncmp (Memc[psfimage],
+            DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+            call imtclose (imlist)
+            call fntclsb (alist)
+            call fntclsb (clist)
+            call imtclose (pimlist)
+            call fntclsb (olist)
+            call fntclsb (oclist)
+            call sfree (sp)
+            call error (0,
+                "Incompatable image and psf file list lengths")
+        }
+
+	# Compare the image and groupfile list lengths.
+        if ((limlist != lolist) && (strncmp (Memc[groupfile],
+            DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+            call imtclose (imlist)
+            call fntclsb (alist)
+            call fntclsb (clist)
+            call imtclose (pimlist)
+            call fntclsb (olist)
+            call fntclsb (oclist)
+            call sfree (sp)
+            call error (0,
+                "Incompatable image and group file list lengths")
+        }
+
+	# Compare the image and output pstfile list lengths.
+        if ((limlist != loclist) && (strncmp (Memc[opstfile],
+            DEF_DEFNAME, DEF_LENDEFNAME) != 0)) {
+            call imtclose (imlist)
+            call fntclsb (alist)
+            call fntclsb (clist)
+            call imtclose (pimlist)
+            call fntclsb (olist)
+            call fntclsb (oclist)
+            call sfree (sp)
+            call error (0,
+                "Incompatable image and output psf file list lengths")
+        }
+
+	# Initialize DAOPHOT main structure, get pset parameters. 
+	call dp_gppars (dao)	
+
+	# Verify the critical parameters and update if appropriate.
+	if (verify == YES) {
+	    call dp_pconfirm (dao)
+	    if (update == YES)
+	        call dp_pppars (dao)
+	}
+
+        # Get the wcs information.
+        wcs = clgwrd ("wcsin", Memc[str], SZ_FNAME, 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_FNAME, 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_apselsetup (dao)
+
+	# Initialize the PSF structure.
+	call dp_fitsetup (dao)
+
+	# Intialize the PSF fitting structure.
+	call dp_psfsetup (dao)
+
+	# Matching algorithm for stars in the psf star list.
+	matchbyid = clgetb ("matchbyid")
+
+	# 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)
+	    showplots = clgetb ("showplots")
+	    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")))
+	    showplots = false
+	    mkstars = false
+	}
+
+	# Open the plot file.
+	if (Memc[plotfile] == EOS)
+	    pfd = NULL
+	else
+	    pfd = open (Memc[plotfile], APPEND, BINARY_FILE)
+	if (pfd != NULL)
+	    mgd = gopen (Memc[graphics], NEW_FILE, pfd)
+	else
+	    mgd = NULL
+
+	# 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 images.
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open 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 list and store the descriptor.
+	    # PSF can read either an APPHOT PHOT file or an ST TABLE 
+	    # file.
+
+	    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_stati (dao, MAXNSTAR), ap_text)
+	    call dp_sets (dao, INPHOTFILE, Memc[outfname])
+
+	    # Open the input photometry list and store the descriptor.
+	    # PSF can read either an APPHOT PHOT file or an ST TABLE 
+	    # file.
+
+	    if (lclist == 0) {
+		pst = NULL
+		Memc[outfname] = EOS
+	    } else {
+	        if (fntgfnb (clist, Memc[pstarfile], SZ_FNAME) == EOF)
+		    call strcpy (DEF_DEFNAME, Memc[pstarfile], SZ_FNAME)
+	        root = fnldir (Memc[pstarfile], Memc[outfname], SZ_FNAME)
+	        if (strncmp (DEF_DEFNAME, Memc[pstarfile+root],
+		    DEF_LENDEFNAME) == 0 || root == strlen (Memc[pstarfile])) 
+	            call dp_inname (Memc[image], Memc[outfname], "pst",
+		        Memc[outfname], SZ_FNAME)
+	        else
+	            call strcpy (Memc[pstarfile], Memc[outfname], SZ_FNAME)
+	        pst_text = itob (access (Memc[outfname], 0, TEXT_FILE))
+	        if (pst_text)
+	            pst = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+	        else
+	            pst = tbtopn (Memc[outfname], READ_ONLY, 0)
+	    }
+	    call dp_sets (dao, COORDS, Memc[outfname])
+
+	    # Open output image containing PSF, output file containing list
+	    # of PSF stars actually used, and the file for PSF neighbors
+	    # If the output is "default", dir$default or a directory
+	    # specification then the extension "psf" is added to the PSF image
+	    # and "psg" to the neighbors file. A suitable version number is
+	    # added to the output name as well. Check that there is enough
+	    # space in the image header user area to hold many PSF stars.
+
+	    if (imtgetim (pimlist, Memc[psfimage], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[psfimage], SZ_FNAME)
+	    root = fnldir (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[psfimage + root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[psfimage])) {
+	        call dp_oimname (Memc[image], Memc[outfname], "psf",
+		    Memc[outfname], SZ_FNAME)
+	    } else
+	        call strcpy (Memc[psfimage], Memc[outfname], SZ_FNAME)
+	    if (envfind ("min_lenuserarea", Memc[str], SZ_FNAME) > 0) {
+	        up = 1
+	        if (ctoi (Memc[str], up, min_lenuserarea) <= 0)
+		    min_lenuserarea = MIN_LENUSERAREA
+	        else
+		    min_lenuserarea = max (MIN_LENUSERAREA, min_lenuserarea)
+            } else
+	        min_lenuserarea = MIN_LENUSERAREA
+	    call dp_pseti (dao, LENUSERAREA, min_lenuserarea)
+	    psfim = immap (Memc[outfname], NEW_IMAGE, min_lenuserarea)
+	    call dp_sets (dao, PSFIMAGE, Memc[outfname])
+
+	    if (fntgfnb (olist, Memc[groupfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[groupfile], SZ_FNAME)
+	    root = fnldir (Memc[groupfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[groupfile + root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[groupfile])) {
+	        call dp_outname (Memc[image], Memc[outfname], "psg",
+		    Memc[outfname], SZ_FNAME)
+	    } else
+	        call strcpy (Memc[groupfile], Memc[outfname], SZ_FNAME)
+	    if (dp_stati (dao, TEXT) == YES)
+	        psfgr = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        psfgr = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTPHOTFILE, Memc[outfname])
+
+
+	    if (fntgfnb (oclist, Memc[opstfile], SZ_FNAME) == EOF)
+		call strcpy (DEF_DEFNAME, Memc[opstfile], SZ_FNAME)
+	    root = fnldir (Memc[opstfile], Memc[outfname], SZ_FNAME)
+	    if (strncmp (DEF_DEFNAME, Memc[opstfile+root],
+	        DEF_LENDEFNAME) == 0 || root == strlen (Memc[opstfile])) {
+	        call dp_outname (Memc[image], Memc[outfname], "pst",
+		    Memc[outfname], SZ_FNAME)
+	    } else
+	        call strcpy (Memc[opstfile], Memc[outfname], SZ_FNAME)
+	    if (dp_stati (dao, TEXT) == YES)
+	        opst = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+	    else
+	        opst = tbtopn (Memc[outfname], NEW_FILE, 0)
+	    call dp_sets (dao, OUTREJFILE, Memc[outfname])
+
+	    # Print banner.
+	    if (DP_VERBOSE(dao) == YES) {
+		call printf ("\nComputing PSF for image: %s\n")
+		    call pargstr (Memc[image])
+		call dp_stats (dao, INPHOTFILE, Memc[str], SZ_FNAME)
+		call printf ("%d stars read from %s\n\n")
+		    call pargi (dp_stati (dao, APNUM))
+		    call pargstr (Memc[str])
+	    }
+
+	    # Read in the PSF star list.
+	    call dp_pseti (dao, PNUM, 0)
+	    if (pst != NULL) {
+	        call dp_rpstars (dao, im, pst, pst_text, gd, mgd, id, mkstars,
+		    matchbyid, showplots)
+	        if (DP_VERBOSE(dao) == YES) {
+		    call dp_stats (dao, COORDS, Memc[str], SZ_FNAME)
+		    call printf ("\n%d PSF stars read from %s\n\n")
+		        call pargi (dp_pstati (dao, PNUM))
+		        call pargstr (Memc[str])
+	        }
+	    }
+
+	    # Make the PSF.
+	    if (Memc[curfile] != EOS)
+	        call dp_mkpsf (dao, im, psfim, opst, psfgr, gd, mgd, id, false,
+	            false, false)
+	    else if (interactive)
+	        call dp_mkpsf (dao, im, psfim, opst, psfgr, gd, mgd, id,
+		    mkstars, true, true)
+	    else if (! dp_updatepsf (dao, im, psfim, opst, psfgr, false, false,
+	        false))
+		call dp_rmpsf (dao, psfim, opst, psfgr)
+
+	    # Close the input image.
+	    call imunmap (im)
+
+	    # Close the input photometry file.
+	    if (apd != NULL) {
+	        if (ap_text)
+	    	    call close (apd)
+	        else
+		    call tbtclo (apd)
+	    }
+
+	    # Close the input psf star file.
+	    if (pst != NULL) {
+	        if (pst_text)
+	    	    call close (pst)
+	        else
+		    call tbtclo (pst)
+	    }
+
+	    # Close PSF image.
+	    if (psfim != NULL)
+	        call imunmap (psfim)
+
+	    # Close the output PSF star file.
+	    if (opst != NULL) {
+	        if (dp_stati (dao, TEXT) == YES)
+		    call close (opst)
+	        else
+	            call tbtclo (opst)
+	    }
+
+	    # Close the group file.
+	    if (psfgr != NULL) {
+	        if (dp_stati (dao, TEXT) == YES)
+		    call close (psfgr)
+	        else
+	            call tbtclo (psfgr)
+	    }
+
+            # 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 (pfd != NULL)
+	    call close (pfd)
+
+	# Close the image / file lists.
+	call imtclose (imlist)
+	call fntclsb (alist)
+	call fntclsb (clist)
+	call imtclose (pimlist)
+	call fntclsb (olist)
+	call fntclsb (oclist)
+
+	# Free the photometry structure.
+	call dp_apclose (dao)
+
+	# Free the PSF fitting structure.
+	call dp_psfclose (dao)
+
+	# Free the PSF structure.
+	call dp_fitclose (dao)
+
+	# Close up the daophot structures.
+	call dp_free (dao)
+
+	call sfree (sp)
+end