Initial commit

1 files changed, 120 insertions, 0 deletions

diff --git a/noao/astutil/pdm/pdmpplot.x b/noao/astutil/pdm/pdmpplot.x
new file mode 100644
index 00000000..78d0ee04
--- /dev/null
+++ b/noao/astutil/pdm/pdmpplot.x
@@ -0,0 +1,120 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+define	MSIZE		2.0		# Mark size
+define	EDGESCL	5.0			# Percent extra space around plot data
+
+# PDM_PPLOT -- Plot the Phase curve.
+
+procedure pdm_pplot (pdmp, period, filename, flip)
+
+pointer	pdmp			# pointer to PDM data structure
+double	period			# period for which to plot the phase curve
+char	filename[SZ_FNAME]	# name of the input data file
+bool	flip			# flip the y-axis scale
+
+int	npt, i, index
+double	frequency
+real	x1, x2, y1, y2, scldif, sclspc
+pointer	gp, xtemp, ytemp, etemp
+pointer	title, system_id, sp
+
+begin
+	npt = PDM_NPT(pdmp)
+	gp = PDM_GP(pdmp)
+	call gclear (PDM_GP(pdmp))
+
+	call smark (sp)
+	call salloc (system_id, SZ_LINE, TY_CHAR)
+	call salloc (title, PDM_SZ_TITLE, TY_CHAR)
+	call salloc (xtemp, 2*npt, TY_REAL)
+	call salloc (ytemp, 2*npt, TY_REAL)
+	call salloc (etemp, 2*npt, TY_REAL)
+
+	do i = 1, 2*npt {
+	    Memr[xtemp+i-1] = PDM_XPH(pdmp,i)
+	    Memr[ytemp+i-1] = PDM_YPH(pdmp,i)
+	    Memr[etemp+i-1] = PDM_PHERR(pdmp,i)
+	}
+
+	# Scale the wcs.
+	call gascale (gp, Memr[xtemp], 2*npt, 1)
+	call gascale (gp, Memr[ytemp], 2*npt, 2)
+
+	# Get the boundaries in X and Y.
+	call ggwind (gp, x1, x2, y1, y2)
+
+	# Add boundry space.
+	scldif = x2 - x1
+	sclspc = scldif * (EDGESCL / 100.)
+	x1 = x1 - sclspc
+	x2 = x2 + sclspc
+	scldif = y2 - y1
+	sclspc = scldif * (EDGESCL / 100.)
+	y1 = y1 - sclspc
+	y2 = y2 + sclspc
+
+	# Flip the y-axis scale if flip = TRUE.
+	if (flip)
+	    call gswind (gp, x1, x2, y2, y1)
+	else
+	    call gswind (gp, x1, x2, y1, y2)
+
+	# Multiline title, save in an array and sprintf to it.
+	# Get the system identification.
+
+	call sysid (Memc[system_id], SZ_LINE)
+
+	# Calculate the frequency.
+	if (period != 0.0)
+	    frequency = 1.0d+0/period
+
+	call sprintf (Memc[title], PDM_SZ_TITLE,
+"%s\nPhase curve at period %12.12g, frequency %12.12g\nFile = %s, numpts = %d")
+	    call pargstr (Memc[system_id])
+	    call pargd (period)
+	    call pargd (frequency)
+	    call pargstr (filename)
+	    call pargi (npt)
+
+	# Draw the axes.
+	call glabax (gp, Memc[title], "phase", "magnitude")
+
+	# Make the plot. If error bars are turned on, draw them.
+	if (PDM_EB(pdmp) == YES && Memr[etemp] > EPSILONR) {
+	    call gpmark (gp, Memr[xtemp], Memr[ytemp],
+			2*npt, GM_CIRCLE, 1.0, 1.0)
+	    call gpmark (gp, Memr[xtemp], Memr[ytemp],
+			2*npt, GM_CIRCLE+GM_FILL, 1.0, 1.0)
+	    do i = 1, 2*npt
+		call gpmark (gp, Memr[xtemp+i-1], Memr[ytemp+i-1],
+			     1, GM_VEBAR, MSIZE, -(2.0*Memr[etemp+i-1]))
+	} else {
+
+	    if (npt <= PDM_PLUSPOINT(pdmp)) {
+		call gpmark (gp, Memr[xtemp], Memr[ytemp],
+			    2*npt, GM_PLUS, MSIZE, MSIZE)
+	    } else {
+		call gpmark (gp, Memr[xtemp], Memr[ytemp],
+			    2*npt, GM_POINT, 1.0, 1.0)
+	    }
+	}
+
+	# Draw an x over any deleted points.
+	do i = 1, npt {
+	    index = PDM_SORT(pdmp,i)
+	    if (PDM_INUSE(pdmp,index) == 0) {
+	        call gscur (gp, Memr[xtemp+i-1]+1, Memr[ytemp+i-1])
+	        call gmark (gp, Memr[xtemp+i-1]+1, Memr[ytemp+i-1], GM_CROSS,
+		    MSIZE, MSIZE)
+	        call gscur (gp, Memr[xtemp+i-1], Memr[ytemp+i-1])
+	        call gmark (gp, Memr[xtemp+i-1], Memr[ytemp+i-1], GM_CROSS,
+		    MSIZE, MSIZE)
+	    }
+	}
+
+	call sfree (sp)
+end