Initial commit

1 files changed, 306 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/polyphot/apyget.x b/noao/digiphot/apphot/polyphot/apyget.x
new file mode 100644
index 00000000..a990b6f8
--- /dev/null
+++ b/noao/digiphot/apphot/polyphot/apyget.x
@@ -0,0 +1,306 @@
+include <fset.h>
+include "../lib/apphot.h"
+include "../lib/fitsky.h"
+include "../lib/polyphot.h"
+
+# AP_YGET -- Procedure to fetch the coordinates of the vertices of a
+# polygon from a text file.
+
+int procedure ap_yget (py, im, fd, delim, x, y, max_nvertices)
+
+pointer	py		# polyphot structure
+pointer	im		# the input image descriptor
+int	fd		# polygon file descriptor
+int	delim		# delimiter character
+real	x[ARB]		# x coords of vertices
+real	y[ARB]		# y coords of vertices
+int	max_nvertices	# maximum number of vertices
+
+real	xc, yc, r2max, r2, xtemp, ytemp
+pointer	sp, str
+int	i, nvertices, nreal, stat
+char	marker
+real	asumr()
+int	fscan(), nscan(), strncmp(), apstati()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Print prompts if input file is STDIN.
+	call fstats (fd, F_FILENAME, Memc[str], SZ_LINE)
+	if (strncmp ("STDIN", Memc[str], 5) == 0) {
+	    call printf ("Type x and y coordinates of vertex\n")
+	    call printf ("\t1 vertex per line\n")
+	    call printf ("\t; to end polygon\n")
+	    call printf ("\t^Z to end list\n")
+	    call flush (STDOUT)
+	}
+
+	# Get the polygon.
+	nvertices = 0
+	stat = fscan (fd)
+	while (stat != EOF) {
+
+	    # Read the vertices from the file
+	    call gargr (xtemp)
+	    call gargr (ytemp)
+	    nreal = nscan ()
+	    if (nreal != 2) {
+		call reset_scan ()
+		call gargc (marker)
+		if (int (marker) == delim)
+		    break
+		nreal = nscan ()
+	    }
+
+	    # Store the vertices.
+	    if (nreal >= 2) {
+	        nvertices = nvertices + 1
+	        if (nvertices <= max_nvertices) {
+		    x[nvertices] = xtemp
+		    y[nvertices] = ytemp
+	        }
+	    }
+
+	    stat = fscan (fd)
+	}
+
+	if (nvertices == 0 && stat == EOF) {
+
+	    nvertices = EOF
+	    call apsetr (py, PYXMEAN, INDEFR)
+	    call apsetr (py, PYYMEAN, INDEFR)
+	    call apsetr (py, OPYXMEAN, INDEFR)
+	    call apsetr (py, OPYYMEAN, INDEFR)
+	    call apsetr (py, PYCX, INDEFR)
+	    call apsetr (py, PYCY, INDEFR)
+	    call apsetr (py, OPYCX, INDEFR)
+	    call apsetr (py, OPYCY, INDEFR)
+	    call apsetr (py, PYMINRAD, INDEFR)
+	    call apseti (py, PYNVER, 0)
+
+	} else if (nvertices <= 2) {
+
+	    call apsetr (py, PYXMEAN, INDEFR)
+	    call apsetr (py, PYYMEAN, INDEFR)
+	    call apsetr (py, OPYXMEAN, INDEFR)
+	    call apsetr (py, OPYYMEAN, INDEFR)
+	    call apsetr (py, PYCX, INDEFR)
+	    call apsetr (py, PYCY, INDEFR)
+	    call apsetr (py, OPYCX, INDEFR)
+	    call apsetr (py, OPYCY, INDEFR)
+	    call apsetr (py, PYMINRAD, INDEFR)
+	    call apseti (py, PYNVER, 0)
+	    nvertices = 0
+
+	} else {
+
+	    # Add in the last vertex.
+	    x[nvertices+1] = x[1]
+	    y[nvertices+1] = y[1]
+
+            # Transform the input coordinates.
+            switch (apstati(py,WCSIN)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_itol (py, x, y, x, y, nvertices + 1)
+            case WCS_TV:
+                call ap_vtol (im, x, y, x, y, nvertices + 1)
+            default:
+                ;
+            }
+
+	    # Compute the mean polygon coordinates.
+	    xc = asumr (x, nvertices) / nvertices
+	    yc = asumr (y, nvertices) / nvertices
+	    call apsetr (py, PYXMEAN, xc)
+	    call apsetr (py, PYYMEAN, yc)
+	    call apsetr (py, PYCX, xc)
+	    call apsetr (py, PYCY, yc)
+
+	    # Set the minimum size of the sky annulus.
+	    r2max = 0.0
+	    do i = 1, nvertices {
+		r2 = (x[i] - xc) ** 2 + (y[i] - yc) ** 2
+		if (r2 > r2max)
+		    r2max = r2
+	    }
+
+            switch (apstati(py,WCSOUT)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_ltoo (py, xc, yc, xc, yc, 1)
+            case WCS_TV:
+                call ap_ltov (im, xc, yc, xc, yc, 1)
+            default:
+    		;
+	    }
+	    call apsetr (py, OPYXMEAN, xc)
+	    call apsetr (py, OPYYMEAN, yc)
+	    call apsetr (py, OPYCX, xc)
+	    call apsetr (py, OPYCY, yc)
+
+	    call apsetr (py, PYMINRAD, (sqrt (r2max) + 1.1))
+	    call apseti (py, PYNVER, nvertices)
+	}
+
+	call sfree (sp)
+	return (nvertices)
+end
+
+
+# AP_YMKPOLY -- Mark the coordinates of a polygon on the display device.
+
+int procedure ap_ymkpoly (py, im, id, x, y, max_nvertices, idflush)
+
+pointer	py		# polyphot structure
+pointer	im		# the input image descriptor
+pointer	id		# display pointer
+real	x[ARB]		# x coords of vertices
+real	y[ARB]		# y coords of vertices
+int	max_nvertices	# maximum number of vertices
+int	idflush		# flush the imd interface ?
+
+int	i, nvertices, stat, wcs, key
+pointer	sp, cmd
+real	xtemp, ytemp, xc, yc, r2max, r2
+real	apstatr(), asumr()
+int	clgcur(), apstati()
+errchk	gscur
+
+begin
+	# Reopen the device.
+	if (id != NULL)
+	    call greactivate (id, 0)
+
+	# Initialize.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Type prompt.
+	call printf (
+	    "Mark polygon vertex [space=mark,q=quit].\n")
+	stat = clgcur ("icommands", xtemp, ytemp, wcs, key, Memc[cmd], SZ_LINE)
+	call ap_vtol (im, xtemp, ytemp, xtemp, ytemp, 1)
+
+
+	# Fetch the polygon and draw it on the display.
+	nvertices = 0
+	while (stat != EOF) {
+
+	    if (key == 'q')
+		break
+
+	    # Decode and draw vertices.
+	    nvertices = nvertices + 1
+	    if (nvertices <= max_nvertices) {
+		x[nvertices] = xtemp
+		y[nvertices] = ytemp
+		if (id != NULL) {
+		    if (nvertices == 1)
+		        call gamove (id, x[1], y[1])
+		    else {
+		        call gadraw (id, x[nvertices], y[nvertices])
+			if (idflush == YES)
+			    call gframe (id)
+			else
+			    call gflush (id)
+		    }
+		}
+	    } else
+		break
+
+	    # Type prompt.
+	    call printf (
+	        "Mark polygon vertex [space=mark,q=quit].\n")
+	    stat = clgcur ("icommands", xtemp, ytemp, wcs, key, Memc[cmd],
+	        SZ_LINE)
+	    call ap_vtol (im, xtemp, ytemp, xtemp, ytemp, 1)
+	}
+	call printf ("\n")
+
+	call sfree (sp)
+
+	# Return EOF or the number of vertices in the polygon.
+	if (stat == EOF) {
+
+	    call apsetr (py, PYXMEAN, INDEFR)
+	    call apsetr (py, PYYMEAN, INDEFR)
+	    call apsetr (py, OPYXMEAN, INDEFR)
+	    call apsetr (py, OPYYMEAN, INDEFR)
+	    call apsetr (py, PYCX, INDEFR)
+	    call apsetr (py, PYCY, INDEFR)
+	    call apsetr (py, OPYCX, INDEFR)
+	    call apsetr (py, OPYCY, INDEFR)
+	    call apsetr (py, PYMINRAD, INDEFR)
+	    call apseti (py, PYNVER, 0)
+	    nvertices = EOF
+
+	} else if (nvertices  <= 2) {
+
+	    call apsetr (py, PYXMEAN, INDEFR)
+	    call apsetr (py, PYYMEAN, INDEFR)
+	    call apsetr (py, OPYXMEAN, INDEFR)
+	    call apsetr (py, OPYYMEAN, INDEFR)
+	    call apsetr (py, PYCX, INDEFR)
+	    call apsetr (py, PYCY, INDEFR)
+	    call apsetr (py, OPYCX, INDEFR)
+	    call apsetr (py, OPYCY, INDEFR)
+	    call apsetr (py, PYMINRAD, INDEFR)
+	    call apseti (py, PYNVER, 0)
+	    nvertices = 0
+
+	} else {
+
+	    # Add the last vertex and close the polygon.
+	    x[nvertices+1] = x[1]
+	    y[nvertices+1] = y[1]
+	    if (id != NULL) {
+	        call gadraw (id, x[1], y[1])
+		if (idflush == YES)
+		    call gframe (id)
+		else
+	            call gflush (id)
+	    }
+
+	    # Compute and save the mean polygon coords.
+	    xc = asumr (x, nvertices) / nvertices
+	    yc = asumr (y, nvertices) / nvertices
+	    call apsetr (py, PYXMEAN, xc)
+	    call apsetr (py, PYYMEAN, yc)
+	    call apsetr (py, PYCX, xc)
+	    call apsetr (py, PYCY, yc)
+	    if (id != NULL)
+	        call gscur (id, apstatr (py, PYCX), apstatr (py, PYCY))
+
+	    # Compute the mean sky annulus.
+	    r2max = 0.0
+	    do i = 1, nvertices {
+		r2 = (x[i] - xc) ** 2 + (y[i] - yc) ** 2
+		if (r2 > r2max)
+		    r2max = r2
+	    }
+
+	    # Compute output coordinate centers.
+            switch (apstati(py,WCSOUT)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_ltoo (py, xc, yc, xc, yc, 1)
+            case WCS_TV:
+                call ap_ltov (im, xc, yc, xc, yc, 1)
+            default:
+                ;
+            }
+            call apsetr (py, OPYXMEAN, xc)
+            call apsetr (py, OPYYMEAN, yc)
+            call apsetr (py, OPYCX, xc)
+            call apsetr (py, OPYCY, yc)
+
+	    call apseti (py, PYNVER, nvertices)
+	    call apsetr (py, PYMINRAD, (sqrt (r2max) + 1.1))
+
+	}
+
+	if (id != NULL)
+	    call gdeactivate (id, 0)
+
+	return (nvertices)
+end