Repatch (from linux) of OSX IRAF

31 files changed, 2710 insertions, 0 deletions

diff --git a/noao/astutil/pdm/README b/noao/astutil/pdm/README
new file mode 100644
index 00000000..b6102f37
--- /dev/null
+++ b/noao/astutil/pdm/README
@@ -0,0 +1,9 @@
+ASTUTIL.PDM -- This directory contains the sources for the PDM program.
+
+PDM finds periods in light curve data by a method called
+Phase Dispersion Minimization.
+
+Reference:
+
+Stellingwerf, R. F., 1978, "Period Determination by Phase Dispersion
+    Minimization", The Astrophysical Journal, 224, pp. 953-960.
diff --git a/noao/astutil/pdm/TODO b/noao/astutil/pdm/TODO
new file mode 100644
index 00000000..4e668a13
--- /dev/null
+++ b/noao/astutil/pdm/TODO
@@ -0,0 +1,4 @@
+The following are things which would be nice to add.
+
+- Compute RMS in theta plots
+- Text output of data points for all plots
diff --git a/noao/astutil/pdm/mkpkg b/noao/astutil/pdm/mkpkg
new file mode 100644
index 00000000..aeb97b97
--- /dev/null
+++ b/noao/astutil/pdm/mkpkg
@@ -0,0 +1,40 @@
+# PDM -- Make the Phase Dispersion Minimization package.
+
+update:
+	$checkout libpkg.a ../
+	$update   libpkg.a
+	$checkin  libpkg.a ../
+	;
+
+
+libpkg.a:
+	pdmalltheta.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmampep.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmautorang.x	pdm.h <ctype.h> <error.h> <mach.h> <pkg/rg.h>
+	pdmbatch.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmclose.x	pdm.h <mach.h>
+	pdmcolon.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmcompare.x	pdm.h <mach.h>
+	pdmcursor.x	<fset.h> <gset.h> <math/curfit.h> pdm.h <ctype.h>\
+			<error.h> <mach.h>
+	pdmdelete.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmdplot.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmfindmin.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmfitphase.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmgdata.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmminmaxp.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmopen.x	pdm.h <mach.h>
+	pdmphase.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmpplot.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmranperm.x	<gset.h> pdm.h <ctype.h> <error.h> <mach.h> <pkg/rg.h>
+	pdmshow.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmsignif.x	<gset.h> pdm.h <ctype.h> <error.h> <mach.h> <pkg/rg.h>
+	pdmsort.x	pdm.h <mach.h>
+	pdmstats.x	pdm.h <ctype.h> <error.h> <mach.h>
+	pdmtheta.x	pdm.h <ctype.h> <error.h> <mach.h> <pkg/rg.h>
+	pdmthetaran.x	pdm.h <ctype.h> <error.h> <mach.h> <pkg/rg.h>
+	pdmtplot.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	pdmundelete.x	pdm.h <ctype.h> <error.h> <gset.h> <mach.h>
+	t_pdm.x         pdm.h <ctype.h> <error.h> <mach.h>
+	;
+
diff --git a/noao/astutil/pdm/pdm.h b/noao/astutil/pdm/pdm.h
new file mode 100644
index 00000000..7ad38d1e
--- /dev/null
+++ b/noao/astutil/pdm/pdm.h
@@ -0,0 +1,77 @@
+# The PDM data structure and other definitions
+
+define	PDM_LENSTRUCT	49		# Length of PDM structure
+
+# Double precision.
+define	PDM_PMIN	Memd[P2D($1)]		# Minimum period to search
+define	PDM_PMAX	Memd[P2D($1+2)]		# Maximum period to search
+define	PDM_FMIN	Memd[P2D($1+4)]		# Minimum frequency to search
+define	PDM_FMAX	Memd[P2D($1+6)]		# Maximum frequency to search
+define	PDM_MINR	Memd[P2D($1+8)]		# Period to remember (min)
+define	PDM_NSIGMA	Memd[P2D($1+10)]	# num std dev. for range break
+define	PDM_SUMSQ	Memd[P2D($1+12)]	# Sum of squares of the data
+define	PDM_DVAR	Memd[P2D($1+14)]	# Variance (s ** 2) of the data
+define	PDM_AMPL	Memd[P2D($1+16)]	# Amplitude of light curve
+define	PDM_EPOCH	Memd[P2D($1+18)]	# Epoch of first maxima in data
+
+# Pointers
+define	PDM_ICD		Memi[$1+20]	# ICFIT pointer for data fit
+define	PDM_ICP		Memi[$1+21]	# ICFIT pointer for phasecurve fit
+define	PDM_CVD		Memi[$1+22]	# CURFIT pointer for data fit
+define	PDM_CVP		Memi[$1+23]	# CURFIT pointer for phasecurve fit
+define	PDM_GP		Memi[$1+24]	# PDM graphics GIO pointer
+define	PDM_LFD		Memi[$1+25]	# Log file descriptor
+define	PDM_PFD		Memi[$1+26]	# Plot file descriptor
+define	PDM_GT		Memi[$1+27]	# PDM gtools pointer
+define	PDM_XP		Memi[$1+28]	# Pointer to data ordinates
+define	PDM_ODYP	Memi[$1+29]	# Pointer to original data abscissas
+define	PDM_DYP		Memi[$1+30]	# Pointer to working data abscissas
+define	PDM_INUSEP	Memi[$1+31]	# Pointer to in-use array
+define	PDM_ERRP	Memi[$1+32]	# Pointer to error bar array
+define	PDM_XTHP	Memi[$1+33]	# Pointer to theta plot ordinates
+define	PDM_YTHP	Memi[$1+34]	# Pointer to theta plot abscissas
+define	PDM_XPHP	Memi[$1+35]	# Pointer to phasecurve plot ordinates
+define	PDM_YPHP	Memi[$1+36]	# Pointer to phasecurve plot abscissas
+define	PDM_PHERRP	Memi[$1+37]	# Pointer to phasecurve plot errors
+define	PDM_SORTP	Memi[$1+38]	# Pointer to array defining sort
+define	PDM_SAMPLEP	Memi[$1+39]	# Pointer to sample (range) string
+define	PDM_RG		Memi[$1+40]	# Pointer to range (sample) structure
+
+# Constants
+define	PDM_NPT		Memi[$1+41]	# Number of data points
+define	PDM_NTHPT	Memi[$1+42]	# Number of theta points
+define	PDM_NRANGE	Memi[$1+43]	# Number of ranges.
+
+# Other
+define	PDM_RESID	Memi[$1+44]	# Using residuals? flag
+define	PDM_RANGE	Memi[$1+45]	# Using ranges? flag
+define	PDM_DEBUG	Memb[$1+46]	# Debug? flag
+define	PDM_PLUSPOINT	Memi[$1+47]	# Threshold data number to use plus
+define	PDM_EB		Memi[$1+48]	# Use error bars? flag
+
+# Macro definitions
+define	PDM_X		Memd[PDM_XP($1)+$2-1]		# data ordinates
+define	PDM_ODY		Memd[PDM_ODYP($1)+$2-1]		# orig data abscissas
+define	PDM_DY		Memd[PDM_DYP($1)+$2-1]		# working data abscissas
+define	PDM_ERR		Memr[PDM_ERRP($1)+$2-1]		# error bars
+define	PDM_INUSE	Memi[PDM_INUSEP($1)+$2-1]	# in-use array
+define	PDM_XTH		Memd[PDM_XTHP($1)+$2-1]		# theta plot ordinates
+define	PDM_YTH		Memd[PDM_YTHP($1)+$2-1]		# theta plot abscissas
+define	PDM_XPH		Memd[PDM_XPHP($1)+$2-1]		# phase plot ordinates
+define	PDM_YPH		Memd[PDM_YPHP($1)+$2-1]		# phase plot abscissas
+define	PDM_PHERR	Memr[PDM_PHERRP($1)+$2-1]	# phase plot errors
+define	PDM_SORT	Memi[PDM_SORTP($1)+$2-1]	# array defining sort
+define	PDM_SAMPLE	Memc[PDM_SAMPLEP($1)]		# sample (range) string
+
+
+# Plot types. (ptype)
+define	DATAPLOT	0		# data plot
+define	THETAPPLOT	1		# theta period plot
+define	THETAFPLOT	2		# theta frequency plot
+define	PHASEPLOT	3		# phase plot
+
+define	MAX_RANGES	20		# maximum number of range segments
+define	PDM_SZ_TITLE	(4*SZ_LINE)	# size of pdm plot title buffer
+define	BIN10		100		# numpts for bin 5/10 split
+					# otherwise, plot pluses
+define	HELP	"noao$lib/scr/pdm.key"	# where help file is
diff --git a/noao/astutil/pdm/pdmalltheta.x b/noao/astutil/pdm/pdmalltheta.x
new file mode 100644
index 00000000..b7a48889
--- /dev/null
+++ b/noao/astutil/pdm/pdmalltheta.x
@@ -0,0 +1,104 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_ALLTHETA -- Calculate the theta statistic for the period window.
+
+procedure pdm_alltheta (pdmp, porf)
+
+pointer	pdmp				# structure pointer
+int	porf				# period or frequency flag
+
+int	i
+double	factor, period, pdm_theta(), mmean, mmin, mmax, scale
+pointer rg, rg_xrangesd()
+errchk	calloc, realloc, rg_xrangesd
+
+begin
+	# Check to see that maxp > minp.
+	if (PDM_PMAX(pdmp) <= PDM_PMIN(pdmp))
+	    call error (0,"maxp smaller or equal minp in alltheta")
+
+	# Calculate the mean of the Y data, remove it.
+	mmean = 0.0d+0
+	do i = 1, PDM_NPT(pdmp) {
+	    mmean = mmean + PDM_DY(pdmp,i)
+	}
+	mmean = mmean/PDM_NPT(pdmp)
+	do i = 1, PDM_NPT(pdmp) {
+	    PDM_DY(pdmp,i) = PDM_DY(pdmp,i) - mmean
+	}
+
+	# Find the min and max of the data, scale the data 0 to 1.
+	mmin = PDM_DY(pdmp,1)
+	mmax = PDM_DY(pdmp,1)
+	do i = 1, PDM_NPT(pdmp) {
+	    if (PDM_DY(pdmp,i) > mmax)
+		mmax = PDM_DY(pdmp,i)
+	    if (PDM_DY(pdmp,i) < mmin)
+		mmin = PDM_DY(pdmp,i)
+	}
+	scale = mmax - mmin
+	do i = 1, PDM_NPT(pdmp) {
+	    PDM_DY(pdmp,i) = (PDM_DY(pdmp,i) - mmin)/scale
+	}
+
+	# Bring the statistics up to date.
+	iferr (call pdm_statistics (pdmp))
+	    call error (0, "alltheta: error in statistics")
+
+	# Allocate space for the ordinate and abscissa vectors for theta
+	# in the pdm data structure.
+
+	if (PDM_XTHP(pdmp) == NULL) {
+	    call calloc (PDM_XTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE)
+	    call calloc (PDM_YTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE)
+	} else {
+	    call realloc (PDM_XTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE)
+	    call realloc (PDM_YTHP(pdmp), PDM_NTHPT(pdmp), TY_DOUBLE)
+	}
+
+	# Calculate constant for incrementing the period.  Give equally
+	# spaced frequencies.
+
+	if (PDM_PMIN(pdmp) <= EPSILOND)
+	    if (PDM_NTHPT(pdmp) >= 1)
+	        PDM_PMIN(pdmp) = PDM_PMAX(PDMP)/PDM_NTHPT(pdmp)
+	    else
+		call error (1, "alltheta: num thpts < 1")
+
+	if (PDM_PMIN(pdmp) >= EPSILONR && PDM_NTHPT(pdmp) != 1)
+	    factor = (PDM_PMAX(pdmp)/PDM_PMIN(pdmp))**
+		(1.0/(PDM_NTHPT(pdmp)-1))
+
+	# Calculate the ranges information from the sample string.
+	rg = rg_xrangesd (PDM_SAMPLE(pdmp), PDM_X(pdmp,1), PDM_NPT(pdmp))
+
+	# Call pdm_theta for each period and store the thetas and periods
+	# in the pdm data structure, calculate frequencies (1/p) as we go.
+
+	period = PDM_PMIN(pdmp)
+	do i = 1, PDM_NTHPT(pdmp) {
+	    PDM_YTH(pdmp,i) = pdm_theta (pdmp, rg, period)
+	    if (porf == THETAPPLOT)
+	        PDM_XTH(pdmp,i) = period
+	    else {
+	        if (period > EPSILOND)
+	            PDM_XTH(pdmp,i) = 1.0d+0/period
+		else
+		    call error (0, "alltheta: period very close to zero")
+	    }
+	    period = period * factor
+	}
+
+	# Remove the scaling.
+	do i = 1, PDM_NPT(pdmp) {
+	    PDM_DY(pdmp,i) = (PDM_DY(pdmp,i) * scale) + mmin
+	}
+
+	# Put the mean back in.
+	do i = 1, PDM_NPT(pdmp) {
+	    PDM_DY(pdmp,i) = PDM_DY(pdmp,i) + mmean
+	}
+end
diff --git a/noao/astutil/pdm/pdmampep.x b/noao/astutil/pdm/pdmampep.x
new file mode 100644
index 00000000..08c2c4e1
--- /dev/null
+++ b/noao/astutil/pdm/pdmampep.x
@@ -0,0 +1,38 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_AMPEP -- Calculate the amplitude and epoch for this data.
+
+procedure pdm_ampep (pdmp, period)
+
+pointer	pdmp			# PDM structure pointer
+double	period			# period for which to calculate
+
+int	i, isave
+double	npt, ymin, ymax
+errchk	pdm_phase
+
+begin
+	npt = PDM_NPT(pdmp)
+
+	# Find the maximum and minimum values in the data.
+	# The difference is the amplitude.
+
+	ymax = -MAX_DOUBLE
+	ymin = MAX_DOUBLE
+	do i = 1, npt {
+	    if (PDM_INUSE(pdmp,i) == 0)
+		next
+	    if (PDM_DY(pdmp,i) < ymin)
+		ymin = PDM_DY(pdmp,i)
+	    if (PDM_DY(pdmp,i) > ymax) {
+		ymax = PDM_DY(pdmp,i)
+		isave = i
+	    }
+	}
+
+	PDM_AMPL(pdmp) = ymax - ymin
+	PDM_EPOCH(pdmp) = PDM_X(pdmp, isave)
+end
diff --git a/noao/astutil/pdm/pdmautorang.x b/noao/astutil/pdm/pdmautorang.x
new file mode 100644
index 00000000..ba27d628
--- /dev/null
+++ b/noao/astutil/pdm/pdmautorang.x
@@ -0,0 +1,101 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <pkg/rg.h>
+include "pdm.h"
+
+# PDM_AUTORANG -- Calculate the ranges division of the data automatically.
+
+int procedure pdm_autorang (pdmp)
+
+pointer pdmp			# PDM structure pointer
+
+int	npt, i, nrng
+double	sumsq, sum, var, sd, maxdif, meansep, rbegin, rend
+int	rngstrt
+pointer	sep, command, sp
+
+begin
+	call smark (sp)
+	call salloc (command, SZ_LINE, TY_CHAR)
+
+	npt = PDM_NPT(pdmp)
+
+	# Calculate the mean and standard deviation of the x-axis
+	# separation of the data points.  (time intervals)
+	# Allocate an array of separations and fill it.
+
+	call salloc (sep, npt-1, TY_DOUBLE)
+	do i = 1, npt-1 {
+	    Memd[sep+i-1] = PDM_X(pdmp,i+1) - PDM_X(pdmp,i)
+	}
+
+	sumsq = 0.0
+	sum = 0.0
+	for (i=1; i<=(npt-1); i=i+1) {
+	    sumsq = sumsq + Memd[sep+i-1]**2	# Sum of squares.
+	    sum = sum + Memd[sep+i-1]
+	}
+	if (npt != 1) {
+	    var = (sumsq - sum**2/npt)/(npt - 1)	# Variance.
+	    sd = var**.5				# Standard Deviation.
+	}
+
+	# Mean separation, maximum time diff.
+	if (npt != 1) {
+	    meansep = (PDM_X(pdmp,npt) - PDM_X(pdmp,1)) / double(npt-1)
+	    maxdif = meansep + sd * PDM_NSIGMA(pdmp)
+	}
+
+	# Look through the separations and if we find one that is more
+	# than nsigma away from the mean on the plus side, divide the
+	# data at this point into another range.
+	
+	nrng = 0
+	rngstrt = 1
+	PDM_SAMPLE(pdmp) = EOS
+	do i = 1, npt - 1 {
+	    if (Memd[sep+i-1] > maxdif) {
+		nrng = nrng + 1
+		if (nrng > MAX_RANGES) {
+		    call sfree (sp)
+		    call error (0,"Max num ranges exceeded in autorange")
+		    break
+		}
+
+		# End of last range = x(i)
+		# If (i+1 != npts) beginning of next range = x(i+1)
+		# Remember where the next range starts.
+
+		rbegin = PDM_X(pdmp,rngstrt)
+		rend = PDM_X(pdmp,i)
+		if ((i+1) < npt)
+		    rngstrt = i+1
+
+		# Sprintf range info at end of string.
+		call sprintf (Memc[command], SZ_LINE, " %g:%g")
+		    call pargd (rbegin)
+		    call pargd (rend)
+		call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE)
+	    }
+	}
+
+	# Finish up last range if needed.
+	If (rngstrt <= npt && rngstrt > 1) {
+	    rbegin = PDM_X(pdmp,rngstrt)
+	    rend = PDM_X(pdmp,i)
+
+	    # Sprintf range info at end of string.
+	    call sprintf (Memc[command], SZ_LINE, " %g:%g")
+		call pargd (rbegin)
+		call pargd (rend)
+	    call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE)
+	}
+
+	# If no ranges found, set the sample string to '*', all data.
+	if (nrng == 0)
+	    call sprintf (PDM_SAMPLE(pdmp), SZ_LINE, "*")
+
+	call sfree (sp)
+	return (nrng)
+end
diff --git a/noao/astutil/pdm/pdmbatch.x b/noao/astutil/pdm/pdmbatch.x
new file mode 100644
index 00000000..6828b783
--- /dev/null
+++ b/noao/astutil/pdm/pdmbatch.x
@@ -0,0 +1,49 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_BATCH -- Batch mode calculation.
+
+procedure pdm_batch (pdmp, file, infile, flip)
+
+pointer	pdmp				# PDM structure pointer
+char	file[SZ_FNAME]			# file to store information
+char	infile[SZ_FNAME]		# input data file name
+bool	flip				# flip y-axis scale
+
+int	fd
+double	pdm_signif(), signif
+bool	verbose
+errchk	pdm_alltheta, pdm_signif, pdm_ampep, pdm_phase
+
+begin
+	# Plot the data.
+	call pdm_dplot (pdmp, infile, flip)
+
+	# Call pdm_alltheta to get the theta array; plot to metafile.
+	call pdm_alltheta (pdmp, THETAPPLOT)
+	call pdm_tplot (pdmp, THETAPPLOT, infile)
+
+	# Call pdm_signif to calculate the significance of the theta minimum.
+	signif = pdm_signif (pdmp, PDM_MINR(pdmp))
+
+	# Write this significance out to the file.
+	fd = PDM_LFD(pdmp)
+	call fprintf (fd, "significance at this minimum = %g\n")
+	    call pargd (signif)
+	call close (fd)
+
+	# Call pdm_amplitudeepoch to get the amplitude and epoch.
+	call pdm_ampep (pdmp, PDM_MINR(pdmp))
+
+	# Call pdm_phase to get the phase curve.
+	# Plot to metafile.
+	
+	call pdm_phase(pdmp, PDM_MINR(pdmp), PDM_EPOCH(pdmp))
+	call pdm_pplot (pdmp, PDM_MINR(pdmp), infile, flip)
+
+	# Call pdm_show to print the parameter information to a file.
+	verbose = TRUE
+	call pdm_show (pdmp, file, verbose)
+end
diff --git a/noao/astutil/pdm/pdmclose.x b/noao/astutil/pdm/pdmclose.x
new file mode 100644
index 00000000..b7e4541e
--- /dev/null
+++ b/noao/astutil/pdm/pdmclose.x
@@ -0,0 +1,63 @@
+include <mach.h>
+include "pdm.h"
+
+# PDM_CLOSE -- Close a PDM data structure.
+
+procedure pdm_close (pdmp, interactive)
+
+pointer	pdmp			# PDM structure pointer
+bool	interactive		# interactive flag
+
+begin
+	# Close icfit pointers, curfit pointers, and graphics pointers,
+	# and the ranges pointer.
+
+	if (PDM_ICD(pdmp) != NULL)
+	    call ic_closed (PDM_ICD(pdmp))
+	if (PDM_ICP(pdmp) != NULL)
+	    call ic_closed (PDM_ICP(pdmp))
+	if (PDM_CVD(pdmp) != NULL)
+	    call dcvfree (PDM_CVD(pdmp))
+	if (PDM_CVP(pdmp) != NULL)
+	    call dcvfree (PDM_CVP(pdmp))
+	if (PDM_GT(pdmp) != NULL)
+	    call gt_free (PDM_GT(pdmp))
+	if (PDM_GP(pdmp) != NULL)
+	    call gclose (PDM_GP(pdmp))
+	if (PDM_RG(pdmp) != NULL)
+	    call rg_free (PDM_RG(pdmp))
+	if (PDM_LFD(pdmp) != NULL)
+	    call close (PDM_LFD(pdmp))
+	if (!interactive)
+	    if (PDM_PFD(pdmp) != NULL)
+		call close (PDM_PFD(pdmp))
+
+	# Free the data vectors.
+	if (PDM_XP(pdmp) != NULL)
+	    call mfree (PDM_XP(pdmp), TY_DOUBLE)
+	if (PDM_ODYP(pdmp) != NULL)
+	    call mfree (PDM_ODYP(pdmp), TY_DOUBLE)
+	if (PDM_DYP(pdmp) != NULL)
+	    call mfree (PDM_DYP(pdmp), TY_DOUBLE)
+	if (PDM_ERRP(pdmp) != NULL)
+	    call mfree (PDM_ERRP(pdmp), TY_REAL)
+	if (PDM_INUSEP(pdmp) != NULL)
+	    call mfree (PDM_INUSEP(pdmp), TY_INT)
+	if (PDM_XTHP(pdmp) != NULL)
+	    call mfree (PDM_XTHP(pdmp), TY_DOUBLE)
+	if (PDM_YTHP(pdmp) != NULL)
+	    call mfree (PDM_YTHP(pdmp), TY_DOUBLE)
+	if (PDM_XPHP(pdmp) != NULL)
+	    call mfree (PDM_XPHP(pdmp), TY_DOUBLE)
+	if (PDM_YPHP(pdmp) != NULL)
+	    call mfree (PDM_YPHP(pdmp), TY_DOUBLE)
+	if (PDM_PHERRP(pdmp) != NULL)
+	    call mfree (PDM_PHERRP(pdmp), TY_REAL)
+	if (PDM_SORTP(pdmp) != NULL)
+	    call mfree (PDM_SORTP(pdmp), TY_INT)
+	if (PDM_SAMPLEP(pdmp) != NULL)
+	    call mfree (PDM_SAMPLEP(pdmp), TY_CHAR)
+
+	# Free the pdm data structure.
+	call mfree (pdmp, TY_STRUCT)
+end
diff --git a/noao/astutil/pdm/pdmcolon.x b/noao/astutil/pdm/pdmcolon.x
new file mode 100644
index 00000000..6d438515
--- /dev/null
+++ b/noao/astutil/pdm/pdmcolon.x
@@ -0,0 +1,292 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+define PDM_KEYWORDS "|show|vshow|minp|maxp|minf|maxf|ntheta|sample|signif\
+		     |ampep|phase|unreject|alldata|origdata|"
+
+define	SHOW		1	# Show parameter settings
+define	VSHOW		2	# Show verbose information
+define	PMIN		3	# Set min search period
+define	PMAX		4	# Set max search period
+define	FMIN		5	# Set min search frequency
+define	FMAX		6	# Set max search frequency
+define	NTHETA		7	# Set number of points for theta
+define	SAMPLE		8	# Set/show the sample ranges
+define	SIGNIF		9	# Find theta significance
+define	AMPEP		10	# Amplitude and Epoch
+define	PHASE		11	# Graph phase curve
+define	UNREJECT	12	# Unreject all the rejected data points
+define	ALLDATA		13	# Reset range to entire dataset
+define	ORIGDATA	14	# Reset data to origional dataset
+
+define	SLOWTHRESH	500	# Threshold above which theta calc gets slow
+
+# PDM_COLON -- Decode colon commands.
+
+procedure pdm_colon (pdmp, cmdstr, ptype, infile, period, flip)
+
+pointer	pdmp				# PDM structure pointer
+char	cmdstr[ARB]			# Command string
+int	ptype				# plot type
+char	infile[SZ_FNAME]		# input file name
+double	period				# current working period
+bool	flip				# flip the y-axis scale
+
+int	nscan(), strdic()
+int	itemp, i
+double	temp, p1, signif, pdm_signif()
+bool	verbose
+pointer	cmd, sp
+errchk	pdm_signif, pdm_ampep, pdm_phase
+
+string	keywords PDM_KEYWORDS
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+
+	# Unpack keyword from string, look up command in dictionary.
+	# Switch on command.  Call the appropriate subroutines.
+
+	switch (strdic (Memc[cmd], Memc[cmd], SZ_FNAME, keywords)) {
+	case SHOW:
+	    # Show parameter settings.
+	    call gargstr (Memc[cmd], SZ_LINE)
+	    verbose = false
+	    if (Memc[cmd] == EOS) {
+		call gdeactivate (PDM_GP(pdmp), AW_CLEAR)
+		iferr (call pdm_show (pdmp, "STDOUT", verbose))
+		    call erract (EA_WARN)
+		call greactivate (PDM_GP(pdmp), AW_PAUSE)
+	    } else {
+		iferr (call pdm_show (pdmp, Memc[cmd], verbose))
+		    call erract (EA_WARN)
+	    }
+	case VSHOW:
+	    # Show verbose information.
+	    call gargstr (Memc[cmd], SZ_LINE)
+	    verbose = true
+	    if (Memc[cmd] == EOS) {
+		call gdeactivate (PDM_GP(pdmp), AW_CLEAR)
+		iferr (call pdm_show (pdmp, "STDOUT", verbose))
+		    call erract (EA_WARN)
+		call greactivate (PDM_GP(pdmp), AW_PAUSE)
+	    } else {
+		iferr (call pdm_show (pdmp, Memc[cmd], verbose))
+		    call erract (EA_WARN)
+	    }
+	case PMIN:
+	    # List or set minimum period.
+	    call gargd (temp)
+	    if (nscan() == 2) {
+		# Set the period minimum in structure.
+		PDM_PMIN(pdmp) = temp
+		if (temp >= EPSILOND)
+		    PDM_FMAX(pdmp) = 1.0d+0/temp
+		# Save this minp out to the parameter.
+		call clputd ("minp", temp)
+	    } else {
+		# Print out the period minimum from the structure.
+		call printf ("Current minimum period = %g\n")
+		    call pargd (PDM_PMIN(pdmp))
+		call flush (STDOUT)
+	    }
+	case PMAX:
+	    # List or set maximum period.
+	    call gargd (temp)
+	    if (nscan() == 2) {
+		# Set the period maximum in structure.
+		PDM_PMAX(pdmp) = temp
+		if (temp >= EPSILOND)
+		    PDM_FMIN(pdmp) = 1.0d+0/temp
+		# Save this minp out to the parameter.
+		call clputd ("maxp", temp)
+	    } else {
+		# Print out the period maximum from the structure.
+		call printf ("Current maximum period = %g\n")
+		    call pargd (PDM_PMAX(pdmp))
+		call flush (STDOUT)
+	    }
+	case FMIN:
+	    # List or set minimum frequency.
+	    call gargd (temp)
+	    if (nscan() == 2) {
+		# Set the frequency minimum in structure.
+		PDM_FMIN(pdmp) = temp
+		if (temp >= EPSILOND)
+		    PDM_PMAX(pdmp) = 1.0d+0/temp
+		# Save this minp out to the parameter.
+		if (temp >= EPSILOND)
+		    call clputd ("maxp", 1.0d+0/temp)
+	    } else {
+		# Print out the frequency minimum from the structure.
+		call printf ("Current minimum frequency = %g\n")
+		    call pargd (PDM_FMIN(pdmp))
+		call flush (STDOUT)
+	    }
+	case FMAX:
+	    # List or set maximum frequency.
+	    call gargd (temp)
+	    if (nscan() == 2) {
+		# Set the frequency maximum in structure.
+		PDM_FMAX(pdmp) = temp
+		if (temp >= EPSILOND)
+		    PDM_PMIN(pdmp) = 1.0d+0/temp
+		# Save this minp out to the parameter.
+		if (temp >= EPSILOND)
+		    call clputd ("minp", 1.0d+0/temp)
+	    } else {
+		# Print out the frequency maximum from the structure.
+		call printf ("Current maximum frequency = %g\n")
+		    call pargd (PDM_FMAX(pdmp))
+		call flush (STDOUT)
+	    }
+	case NTHETA:
+	    # Set/show number of theta points
+	    call gargi (itemp)
+	    if (nscan() == 2) {
+		# Set ntheta in structure.
+		PDM_NTHPT(pdmp) = itemp
+		if (itemp > SLOWTHRESH)
+		    # Give message saying that with this ntheta, the
+		    # theta calculation will take quite a while.
+		    call printf ("Large ntheta => long calculation time \007\n")
+	    } else {
+		# Print out the value of ntheta from the structure.
+		call printf ("Number of theta points = %g\n")
+		    call pargi (PDM_NTHPT(pdmp))
+		call flush (STDOUT)
+	    }
+	case SAMPLE:
+	    # List or set the sample points.
+	    call gargstr (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+	        call printf ("sample = %s\n")
+		    call pargstr (PDM_SAMPLE(pdmp))
+	    } else {
+		if (ptype == DATAPLOT) {
+		    call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp),
+			PDM_X(pdmp,1), PDM_NPT(pdmp), 0)
+		}
+		call strcpy (Memc[cmd], PDM_SAMPLE(pdmp), SZ_LINE)
+		if (ptype == DATAPLOT) {
+		    call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp),
+			PDM_X(pdmp,1), PDM_NPT(pdmp), 1)
+		}
+	    }
+	case SIGNIF:
+	    # Calculate the significance of theta at period.
+	    p1 = 0.0
+	    call gargd (temp)
+	    if (nscan() == 2)		# User entered a period.
+		p1 = temp
+	    else {
+		# Use remembered period.
+		if (PDM_MINR(pdmp) >= EPSILOND)
+		    p1 = PDM_MINR(pdmp)
+		else {
+		    call printf ("No remembered minimum period. \007\n")
+		    call flush (STDOUT)
+		}
+	    }
+	    # Calculate significance at cursor x position (per).
+	    if (p1 >= EPSILOND) {
+	        signif = pdm_signif (pdmp, p1)
+	        # Print at bottom of screen.
+	        call printf ("Significance at period %g = %g\n")
+		    call pargd (p1)
+		    call pargd (signif)
+	        call flush (STDOUT)
+	    }
+	case AMPEP:
+	    # Calculate the amplitude and epoch for the data.
+	    p1 = 0.0
+	    call gargd (temp)
+	    if (nscan() == 2)		# User entered a period.
+		p1 = temp
+	    else {
+		# Use remembered period.
+		if (PDM_MINR(pdmp) >= EPSILOND)
+		    p1 = PDM_MINR(pdmp)
+		else {
+		    call printf ("No remembered minimum period. \007\n")
+		    call flush (STDOUT)
+		}
+	    }
+	    # Calculate ampl & epoch at p1.
+	    if (p1 >= EPSILOND) {
+	        call pdm_ampep (pdmp, p1)
+	        # Print at bottom of screen.
+	        call printf("amplitude of data at period %g = %g, epoch = %g\n")
+		    call pargd (p1)
+		    call pargd (PDM_AMPL(pdmp))
+		    call pargd (PDM_EPOCH(pdmp))
+	        call flush (STDOUT)
+	    }
+	case PHASE:
+	    # Phase curve plot.
+	    call gargd (temp)
+	    if (nscan() == 2) {
+		# Calculate the phase curve, then make the plot.
+		call pdm_ampep (pdmp, temp)
+		call pdm_phase(pdmp, temp, PDM_EPOCH(pdmp))
+		call pdm_pplot (pdmp, temp, infile, flip)
+		period = temp
+	    } else {
+		# Use remembered period.
+		if (PDM_MINR(pdmp) >= EPSILOND) {
+		    call pdm_ampep (pdmp, PDM_MINR(pdmp))
+		    call pdm_phase(pdmp, PDM_MINR(pdmp), PDM_EPOCH(pdmp))
+		    call pdm_pplot (pdmp, PDM_MINR(pdmp), infile, flip)
+		} else {
+		    call printf ("No remembered minimum period. \007\n")
+		    call flush (STDOUT)
+		}
+	    }
+	    ptype = PHASEPLOT
+	case UNREJECT:
+	    # Copy original data vector into working data vector and
+	    # set inuse array to all ones.
+
+	    do i = 1, PDM_NPT(pdmp) {
+	        PDM_INUSE(pdmp,i) = 1
+	        PDM_DY(pdmp,i) = PDM_ODY(pdmp,i)
+	    }
+	    PDM_RESID(pdmp) = NO
+
+	    # Check type of plot and replot.
+	    if (ptype == DATAPLOT)
+		call pdm_dplot (pdmp, infile, flip)
+	    if (ptype == PHASEPLOT)
+		call pdm_pplot (pdmp, period, infile, flip)
+	case ALLDATA:
+	    # Initialize the sample string and erase from the graph.
+	    if (ptype == DATAPLOT) {
+		call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp),
+		    PDM_X(pdmp,1), PDM_NPT(pdmp), 0)
+	    }
+	    call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE)
+	case ORIGDATA:
+	    # Copy the original data vector into the working data vector.
+	    do i = 1, PDM_NPT(pdmp)
+		PDM_DY(pdmp,i) = PDM_ODY(pdmp,i)
+	    PDM_RESID(pdmp) = NO
+
+	    # Replot data.
+	    call pdm_dplot (pdmp, infile, flip)
+	    ptype = DATAPLOT
+	default:
+	    # Error, unknown colon command; ring bell.
+	    call printf ("\007")
+	    call printf ("v/show minp/f maxp/f ntheta sample ")
+	    call printf ("signif ampep phase unreject all/origdata\n")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/astutil/pdm/pdmcompare.x b/noao/astutil/pdm/pdmcompare.x
new file mode 100644
index 00000000..b2d60978
--- /dev/null
+++ b/noao/astutil/pdm/pdmcompare.x
@@ -0,0 +1,25 @@
+include <mach.h>
+include "pdm.h"
+
+# Compare procedure for qsort.
+
+int procedure pdm_compare (item1, item2)
+
+int	item1		# index of first phase
+int	item2		# index of second phase
+
+pointer	comarray
+common /sortcom/ comarray
+double	p1, p2
+
+begin
+	p1 = Memd[comarray+item1-1]
+	p2 = Memd[comarray+item2-1]
+
+	if (p1 > p2)
+	    return (1)
+	if (p1 == p2)
+	    return (0)
+	if (p1 < p2)
+	    return (-1)
+end
diff --git a/noao/astutil/pdm/pdmcursor.x b/noao/astutil/pdm/pdmcursor.x
new file mode 100644
index 00000000..067ce966
--- /dev/null
+++ b/noao/astutil/pdm/pdmcursor.x
@@ -0,0 +1,383 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include	<gset.h>
+include	<fset.h>
+include	<math/curfit.h>
+include "pdm.h"
+
+define	PROMPT		"pdm options"
+
+# PDM_CURSOR -- Get the next command from the user in a graphics cursor loop.
+# Perform the requested function.
+
+procedure pdm_cursor (pdmp, ptype, infile, flip)
+
+pointer	pdmp			# pointer to PDM data structure
+int	ptype			# type of plot on the screen
+char	infile[SZ_FNAME]	# input file name
+bool	flip			# flip the y-axis scale
+
+real	xx, yy
+double	x, y			# cursor coordinates
+double	xmax, xmin
+double	period
+int	wcs			# wcs to which coordinates belong
+int	key			# keystroke value of cursor event
+int	stridxs(), gqverify(), ier
+int	clgcur(), i
+int	pdm_delete(), index
+int	pdm_undelete(), pdm_findmin()
+double	rx1, rx2
+double	dcveval()
+double	pdm_signif(), signif
+pointer	weights			# pointer to temporary weights array for icfit
+pointer	command			# string value, if any
+pointer	sp, sptemp
+errchk	pdm_colon, icg_fit, pdm_fitphase, pdm_alltheta
+errchk	pdm_phase, pdm_signif, pdm_ampep, pdm_findmin
+
+begin
+	call smark (sp)
+	call salloc (command, SZ_LINE, TY_CHAR)
+
+	while (clgcur ("cursor", xx, yy, wcs, key, Memc[command],
+	    SZ_LINE) != EOF) {
+
+	    x = double(xx)
+	    y = double(yy)
+
+	    # Switch on command, take appropriate action.
+	    switch (key) {
+	    case '?':
+	        # List options.
+		call gpagefile (PDM_GP(pdmp), HELP, PROMPT)
+	    case ':':
+		# Colon command.
+		call pdm_colon (pdmp, Memc[command], ptype, infile,
+		    period, flip)
+	    case 'h':
+		# Graph the data.
+		call pdm_dplot (pdmp, infile, flip)
+		ptype = DATAPLOT
+	    case 'f':
+		# Call icfit on the data.
+		if (ptype == DATAPLOT) {
+		    # Set the min/max ordinate values for icfit.
+		    call alimd (PDM_X(pdmp,1), PDM_NPT(pdmp), xmin, xmax)
+		    call ic_putr (PDM_ICD(pdmp), "xmin", real(xmin))
+		    call ic_putr (PDM_ICD(pdmp), "xmax", real(xmax))
+
+		    # Allocate a temporary weights array and fill it with the
+		    # in-use array values.
+
+		    call smark (sptemp)
+		    call salloc (weights, PDM_NPT(pdmp), TY_DOUBLE)
+		    do i = 1, PDM_NPT(pdmp)
+			Memd[weights+i-1] = double(PDM_INUSE(pdmp,i))
+
+		    # Call icfit.
+		    if (PDM_NPT(pdmp) >= 2)
+	    		call icg_fitd (PDM_ICD(pdmp), PDM_GP(pdmp), "cursor",
+			    PDM_GT(pdmp), PDM_CVD(pdmp), PDM_X(pdmp,1),
+			    PDM_DY(pdmp,1), Memd[weights], PDM_NPT(pdmp))
+
+		    # Recover the weights array back into the in-use array.
+		    do i = 1, PDM_NPT(pdmp)
+			 PDM_INUSE(pdmp,i) = int(Memd[weights+i-1])
+
+		    call sfree (sptemp)
+
+		    # Replot.
+		    call pdm_dplot (pdmp, infile, flip)
+		    call sfree (sptemp)
+		} else if (ptype == PHASEPLOT) {
+		    # Call icfit on the phases.
+		    call pdm_fitphase (pdmp)
+
+		    # Replot.
+		    call pdm_pplot (pdmp, period, infile, flip)
+		} else
+		    call printf ("Can't fit a THETA plot \007\n")
+	    case 'i':
+		# Theta vs. frequency plot.
+		if (PDM_PMIN(pdmp) < EPSILOND && PDM_PMAX(pdmp) < EPSILOND)
+		    call pdm_minmaxp(pdmp)
+
+		# Calculate theta.
+		call pdm_alltheta (pdmp, THETAFPLOT)
+
+		# Graph theta vs frequency.
+		call pdm_tplot (pdmp, THETAFPLOT, infile)
+		ptype = THETAFPLOT
+	    case 'k':
+		# Theta vs. period plot.
+		if (PDM_PMIN(pdmp) < EPSILOND && PDM_PMAX(pdmp) < EPSILOND)
+		    call pdm_minmaxp(pdmp)
+
+		# Calculate theta.
+		call pdm_alltheta (pdmp, THETAPPLOT)
+
+		# Graph theta vs frequency.
+		call pdm_tplot (pdmp, THETAPPLOT, infile)
+		ptype = THETAPPLOT
+	    case 'p':
+		# Phase curve plot.
+		if (ptype == THETAPPLOT) {
+		    # Find the epoch, calculate the phases w.r.t. this epoch.
+		    call pdm_ampep (pdmp, x)
+		    call pdm_phase(pdmp, x, PDM_EPOCH(pdmp))
+		    call pdm_pplot (pdmp, x, infile, flip)
+		    ptype = PHASEPLOT
+		    period = x
+		} else if (ptype == THETAFPLOT) {
+		    # Find the epoch, calculate the phases w.r.t. this epoch.
+		    call pdm_ampep (pdmp, 1.0d+0/x)
+		    call pdm_phase(pdmp, 1.0d+0/x, PDM_EPOCH(pdmp))
+		    call pdm_pplot (pdmp, 1.0d+0/x, infile, flip)
+		    ptype = PHASEPLOT
+		    period = 1.0d+0/x
+		} else
+		    call printf ("Wrong type of plot on screen for p key\007\n")
+	    case 'd':
+		# Delete the point nearest the cursor.
+		index = pdm_delete (pdmp, x, y, ptype)
+	    case 'u':
+		# Undelete the point nearest the cursor.
+		index = pdm_undelete (pdmp, x, y, ptype)
+	    case 'j':
+		# Subtract the fit from the data and use residuals.
+		if (PDM_CVD(pdmp) == NULL)
+		    call printf ("Fit has not been done. \007\n")
+		else if (PDM_RESID(pdmp) == YES)
+		    call printf ("Already using residuals. \007\n")
+		else {
+		    # For each point, calculate the fit function and subtract
+		    # it from the data.
+
+		    do i = 1, PDM_NPT(pdmp) {
+			PDM_DY(pdmp,i) = PDM_DY(pdmp,i) -
+			    dcveval (PDM_CVD(pdmp), PDM_X(pdmp,i))
+		    }
+		    PDM_RESID(pdmp) = YES
+		    if (ptype == DATAPLOT)
+		        call pdm_dplot (pdmp, infile, flip)
+		}
+	    case 's':
+		# Set sample regions with the cursor.
+		if (ptype == DATAPLOT) {
+		    if (stridxs ("*", PDM_SAMPLE(pdmp)) > 0)
+		        PDM_SAMPLE(pdmp) = EOS
+
+		    rx1 = x
+		    call printf ("again:\n")
+		    if (clgcur ("cursor", xx, yy, wcs, key, Memc[command],
+			SZ_LINE) == EOF)
+		        break
+		    rx2 = double(xx)
+
+		    call sprintf (Memc[command], SZ_LINE, " %g:%g")
+		        call pargd (rx1)
+		        call pargd (rx2)
+
+		    call strcat (Memc[command], PDM_SAMPLE(pdmp), SZ_LINE)
+		    call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp),
+			PDM_X(pdmp,1), PDM_NPT(pdmp), 1)
+		    call printf (" \n")
+		    call gflush (PDM_GP(pdmp))
+		} else
+		    call printf ("Wrong type of plot for s key \007\n")
+	    case 't':
+		# Initialize the sample string and erase from the graph.
+		if (ptype == DATAPLOT) {
+		    call rg_gxmarkd (PDM_GP(pdmp), PDM_SAMPLE(pdmp),
+			PDM_X(pdmp,1), PDM_NPT(pdmp), 0)
+		}
+		call gflush (PDM_GP(pdmp))
+		call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE)
+		call gflush (PDM_GP(pdmp))
+	    case 'g':
+		# Significance of theta at cursor x position.
+		if (ptype == THETAPPLOT) {
+	    	    # Calculate significance at cursor x position (per).
+		        signif = pdm_signif (pdmp, x)
+	    	    # Print at bottom of screen.
+		    call printf ("Significance at cursor = %g\n")
+			call pargd (signif)
+		} else if (ptype == THETAFPLOT) {
+	    	    # Calculate significance at cursor x position (per).
+		        signif = pdm_signif (pdmp, 1.0d+0/x)
+		    # Print at bottom of screen.
+		    call printf ("Significance at cursor = %g\n")
+			call pargd (signif)
+		} else if (ptype == PHASEPLOT) {
+		    # Calculate significance at current period/frequency
+		        signif = pdm_signif (pdmp, period)
+		    # Print at bottom of screen
+		    call printf ("Significance of this period = %g\n")
+			call pargd (signif)
+		} else {  
+		    # Data plot.
+		    call printf ("Wrong type of plot for g key \007\n")
+		}
+	    case 'a':
+		# Amplitude and epoch at cursor x position.
+		if (ptype == THETAPPLOT) {
+		    # Calculate ampl & epoch at cursor x position.
+		    call pdm_ampep (pdmp, x)
+		    # Print at bottom of screen.
+		    call printf ("amplitude of data = %g, epoch = %g\n")
+			call pargd (PDM_AMPL(pdmp))
+			call pargd (PDM_EPOCH(pdmp))
+		} else if (ptype == THETAFPLOT) {
+		    # Calculate ampl & epoch at cursor x position.
+		    call pdm_ampep (pdmp, 1.0d+0/x)
+		    # Print at bottom of screen.
+		    call printf ("amplitude of data = %g, epoch = %g\n")
+			call pargd (PDM_AMPL(pdmp))
+			call pargd (PDM_EPOCH(pdmp))
+		} else if (ptype == PHASEPLOT) {
+		    # Calculate ampl & epoch at current period/frequency
+		    call pdm_ampep (pdmp, period)
+		    # Print at bottom of screen.
+		    call printf ("amplitude of data = %g, epoch = %g\n")
+			call pargd (PDM_AMPL(pdmp))
+			call pargd (PDM_EPOCH(pdmp))
+		} else {
+		    # Data plot.
+		    call printf ("Wrong type of plot for a key \007\n")
+		}
+	    case ',':
+		# Set minp or minf to cursor x position.
+		if (ptype == THETAPPLOT) {
+		    PDM_PMIN(pdmp) = x
+		    PDM_FMIN(pdmp) = 1.0d+0/x
+		    # Print at bottom of screen.
+		    call printf ("minp now %g\n")
+			call pargd (PDM_PMIN(pdmp))
+		} else if (ptype == THETAFPLOT) {
+		    PDM_FMAX(pdmp) = x
+		    PDM_PMAX(pdmp) = 1.0d+0/x
+		    # Print at bottom of screen.
+		    call printf ("minf now %g\n")
+			call pargd (PDM_FMAX(pdmp))
+		} else {
+		    # Data plot or phase plot.
+		    call printf ("Wrong type of plot for , key \007\n")
+		}
+	    case '.':
+		# Set maxp or maxf to cursor x position.
+		if (ptype == THETAPPLOT) {
+		    PDM_PMAX(pdmp) = x
+		    PDM_FMAX(pdmp) = 1.0d+0/x
+		    # Print at bottom of screen.
+		    call printf ("maxp now %g\n")
+			call pargd (PDM_PMAX(pdmp))
+		} else if (ptype == THETAFPLOT) {
+		    PDM_FMIN(pdmp) = x
+		    PDM_PMIN(pdmp) = 1.0d+0/x
+		    # Print at bottom of screen.
+		    call printf ("maxf now %g\n")
+			call pargd (PDM_FMIN(pdmp))
+		} else {
+		    # Data plot or phase plot.
+		    call printf ("Wrong type of plot for . key \007\n")
+		}
+	    case 'm':
+		# Mark range and find minimum in this range.
+		if (ptype == THETAFPLOT || ptype == THETAPPLOT) {
+		    rx1 = x
+		    call printf ("again:\n")
+		    if (clgcur ("cursor", xx, yy, wcs, key, Memc[command],
+			SZ_LINE) == EOF)
+		        break
+		    rx2 = double(xx)
+		    index = pdm_findmin(pdmp, ptype, rx1, rx2, 1,
+			PDM_NTHPT(pdmp))
+	    	    PDM_MINR(pdmp) = PDM_XTH(pdmp,index)
+		    call printf ("period at minimum = %g, frequency = %g\n")
+			call pargd (PDM_XTH(pdmp,index))
+			call pargd (1.0d+0/PDM_XTH(pdmp,index))
+
+		} else
+		    call printf ("Wrong type of plot for m key. \007\n")
+	    case 'r':
+		# Check type of plot and replot.
+		if (ptype == DATAPLOT)
+		    call pdm_dplot (pdmp, infile, flip)
+		if (ptype == THETAPPLOT)
+		    call pdm_tplot (pdmp, THETAPPLOT, infile)
+		if (ptype == THETAFPLOT)
+		    call pdm_tplot (pdmp, THETAFPLOT, infile)
+		if (ptype == PHASEPLOT)
+		    call pdm_pplot (pdmp, period, infile, flip)
+	    case 'e':
+		# Toggle error bars.
+		if (PDM_EB(pdmp) == NO)
+		    PDM_EB(pdmp) = YES
+		else
+		    PDM_EB(pdmp) = NO
+
+		# Check type of plot and replot.
+		if (ptype == PHASEPLOT)
+		    call pdm_pplot (pdmp, period, infile, flip)
+	    case 'x':
+		# Remove a trend from the data by fitting a straight line to
+		# the data and removing it.
+
+		# Set the min/max ordinate values for icfit.
+		call alimd (PDM_X(pdmp,1), PDM_NPT(pdmp), xmin, xmax)
+		call dcvinit (PDM_CVD(pdmp), SPLINE1, 1, xmin, xmax)
+
+		# Allocate a weights array and fill it.
+		call smark (sptemp)
+		call salloc (weights, PDM_NPT(pdmp), TY_DOUBLE)
+		do i = 1, PDM_NPT(pdmp)
+		    Memd[weights+i-1] = PDM_INUSE(pdmp,i)
+
+		call dcvfit (PDM_CVD(pdmp), PDM_X(pdmp,1), PDM_DY(pdmp,1),
+		    Memd[weights], PDM_NPT(pdmp), WTS_USER, ier)
+		if (ier != 0) {
+		    call eprintf ("error in dcvfit\n")
+		    call erract (EA_WARN)
+		}
+
+		# Subtract the fit from the data.
+		if (PDM_RESID(pdmp) == YES) {
+		    call printf ("Already using residuals. \007\n")
+		    next
+		} else {
+		    # For each point, calculate the fit function and subtract
+		    # it from the data.
+
+		    do i = 1, PDM_NPT(pdmp) {
+			PDM_DY(pdmp,i) = PDM_DY(pdmp,i) -
+			    dcveval (PDM_CVD(pdmp), PDM_X(pdmp,i))
+		    }
+		    PDM_RESID(pdmp) = YES
+		}
+
+		call dcvfree (PDM_CVD(pdmp))
+		call sfree (sptemp)
+
+		# Replot.
+		call pdm_dplot (pdmp, infile, flip)
+		ptype = DATAPLOT
+	    case 'z':
+		# Flip the y-axis scale.
+		flip = !flip
+	    case 'q':
+		# Quit. Exit PDM.
+		if (gqverify() == YES)
+		    break
+	    default:
+		# Error: unknown command: ring bell.
+		call printf ("\007\n")
+		call printf ("? for help or (h,f,i,k,p")
+		call printf (",d,u,j,s,t,g,a,m,r,x,q,:)\n")
+	    }
+	}
+
+	call flush (STDOUT)
+	call sfree (sp)
+end
diff --git a/noao/astutil/pdm/pdmdelete.x b/noao/astutil/pdm/pdmdelete.x
new file mode 100644
index 00000000..5052a2e7
--- /dev/null
+++ b/noao/astutil/pdm/pdmdelete.x
@@ -0,0 +1,103 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+define	MSIZE		2.0		# Mark size
+
+# PDM_DELETE -- Delete the point from the plot and set it's inuse array
+# entry to zero.
+
+int procedure pdm_delete (pdmp, cx, cy, ptype)
+
+pointer	pdmp			# pointer to PDM data structure
+double	cx, cy			# device cursor coordinates
+int	ptype			# plot type
+
+pointer	gp
+real	x, y
+int	npts, i, j, index
+real	x0, y0, r2, r2min
+
+begin
+	gp = PDM_GP(pdmp)
+	npts = PDM_NPT(pdmp)
+
+	if (ptype == DATAPLOT) {
+	    # Transform world cursor coordinates to NDC.
+	    call gctran (gp, real(cx), real(cy), x, y, 1, 0)
+
+	    # Search for nearest point in-use.
+	    j = 0
+	    r2min = MAX_REAL
+	    do i = 1, npts {
+	        if (PDM_INUSE(pdmp,i) == 0)
+		    next
+
+	        call gctran (gp, real(PDM_X(pdmp,i)), real(PDM_DY(pdmp,i)),
+			     x0, y0, 1, 0)
+		
+	        r2 = (x0 - x) ** 2 + (y0 - y) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+
+	    # Mark the deleted point with a cross and set the weight to zero.
+	    if (j != 0) {
+	        call gscur (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)))
+	        call gmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)),
+		    GM_CROSS, MSIZE, MSIZE)
+	        PDM_INUSE(pdmp,j) = 0
+		call gflush(gp)
+	    }
+
+	    return (j)
+
+	} else if (ptype == PHASEPLOT) {
+	    # Transform world cursor coordinates to NDC.
+	    call gctran (gp, real(cx), real(cy), x, y, 1, 0)
+
+	    # Search for nearest point in-use.
+	    j = 0
+	    r2min = MAX_REAL
+	    do i = 1, npts {
+		index = PDM_SORT(pdmp,i)
+	        if (PDM_INUSE(pdmp,index) == 0)
+		    next
+
+	        call gctran (gp, real(PDM_XPH(pdmp,i)), real(PDM_YPH(pdmp,i)),
+			     x0, y0, 1, 0)
+		
+	        r2 = (x0 - x) ** 2 + (y0 - y) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+
+	    # Mark the deleted point with a cross and set the weight to zero.
+	    # We mark two points since all of the points are displayed twice
+	    # on the phase plot.
+	    
+	    if (j != 0) {
+	        call gscur (gp, real(PDM_XPH(pdmp,j))+1.0,
+				real(PDM_YPH(pdmp,j)))
+	        call gmark (gp, real(PDM_XPH(pdmp,j))+1.0,
+			        real(PDM_YPH(pdmp,j)), GM_CROSS, MSIZE, MSIZE)
+	        call gscur (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)))
+	        call gmark (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)),
+		    GM_CROSS, MSIZE, MSIZE)
+		
+		# Calculate which point this corresponds to.
+		index = PDM_SORT(pdmp,j)
+	        PDM_INUSE(pdmp,index) = 0
+		call gflush (gp)
+	    }
+
+	    return (index)
+	} else
+	    return (0)
+end
diff --git a/noao/astutil/pdm/pdmdplot.x b/noao/astutil/pdm/pdmdplot.x
new file mode 100644
index 00000000..c0b53d19
--- /dev/null
+++ b/noao/astutil/pdm/pdmdplot.x
@@ -0,0 +1,101 @@
+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_DPLOT -- Plot the data on the screen.
+
+procedure pdm_dplot (pdmp, filename, flip)
+
+pointer	pdmp				# pointer to PDM data structure
+char	filename[SZ_FNAME]		# name of the input data file
+bool	flip				# flip the y-axis scale
+
+int	npt, i
+real	x1, x2, y1, y2, scldif, sclspc
+pointer	gp, xtemp, ytemp
+pointer	title
+pointer	system_id, sp
+
+begin
+	call smark (sp)
+	call salloc (system_id, SZ_LINE, TY_CHAR)
+	call salloc (title, 4*SZ_LINE, TY_CHAR)
+	call salloc (xtemp, PDM_NPT(pdmp), TY_REAL)
+	call salloc (ytemp, PDM_NPT(pdmp), TY_REAL)
+
+	npt = PDM_NPT(pdmp)
+	gp = PDM_GP(pdmp)
+	call gclear (gp)
+
+	# Scale the wcs.
+	do i = 1, PDM_NPT(pdmp) {
+	    Memr[xtemp+i-1] = real(PDM_X(pdmp,i))
+	    Memr[ytemp+i-1] = real(PDM_DY(pdmp,i))
+	}
+	call gascale (gp, Memr[xtemp], npt, 1)
+	call gascale (gp, Memr[ytemp], npt, 2)
+
+	# Get the X and Y boundaries.
+	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)
+	call sprintf (Memc[title], 4*SZ_LINE,
+	    "%s\nFile = %s\n%s\nnumpts = %d")
+	    call pargstr (Memc[system_id])
+	    call pargstr (filename)
+	    if (PDM_RESID(pdmp) == YES)
+	        call pargstr ("Data with fit removed")
+	    else
+	        call pargstr ("Data")
+	    call pargi (npt)
+
+	# Draw the axes.
+	call glabax (gp, Memc[title], "obs time", "magnitude")
+
+	# Make the plot.
+	if (npt <= PDM_PLUSPOINT(pdmp)) {
+	    call gpmark (gp, Memr[xtemp], Memr[ytemp], npt,
+			GM_PLUS, MSIZE, MSIZE)
+	} else {
+	    call gpmark (gp, Memr[xtemp], Memr[ytemp], npt,
+			GM_POINT, 1.0, 1.0)
+	}
+
+	# Call the routine to mark the ranges if they are in effect.
+	call rg_gxmarkd (gp, PDM_SAMPLE(pdmp), PDM_X(pdmp,1), npt, 1)
+
+	# Draw an x over any deleted points.
+	do i = 1, npt {
+	    if (PDM_INUSE(pdmp,i) == 0) {
+	        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
diff --git a/noao/astutil/pdm/pdmfindmin.x b/noao/astutil/pdm/pdmfindmin.x
new file mode 100644
index 00000000..4053bf43
--- /dev/null
+++ b/noao/astutil/pdm/pdmfindmin.x
@@ -0,0 +1,57 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_FINDMIN -- Find the minimum value of the abscissa.
+
+int procedure pdm_findmin (pdmp, ptype, startint, endint, is, ie)
+
+pointer	pdmp			# pointer to PDM data structure
+int	ptype			# plot type
+double	startint, endint	# start and end ordinates
+int	is, ie			# start and end indexs
+
+double	miny, dy, dx
+int	i, isave, npt
+pointer	xpt, ypt
+
+begin
+	# Dereference npt.
+	npt = PDM_NPT(pdmp)
+
+	# Dereference the appropriate abcissa and ordinate.
+	switch (ptype) {
+	case DATAPLOT:
+	    xpt = PDM_XP(pdmp)
+	    ypt = PDM_DYP(pdmp)
+	case THETAPPLOT:
+	    xpt = PDM_XTHP(pdmp)
+	    ypt = PDM_YTHP(pdmp)
+	case THETAFPLOT:
+	    ypt = PDM_YTHP(pdmp)
+	    xpt = PDM_XTHP(pdmp)
+	case PHASEPLOT:
+	    ypt = PDM_YPHP(pdmp)
+	    xpt = PDM_XPHP(pdmp)
+	}
+
+	# Search the abscissas between startint and endint
+	# for the minimum value.
+
+	isave = 1
+	miny = MAX_DOUBLE
+	do i = is, ie {
+	    dx = Memd[xpt+i-1]
+	    dy = Memd[ypt+i-1]
+	    if (dx > startint && dx < endint) {
+		if (dy < miny) {
+	    	    miny = dy
+	    	    isave = i
+		}
+	    }
+	}
+
+	# Return the corresponding index value.
+	return (isave)
+end
diff --git a/noao/astutil/pdm/pdmfitphase.x b/noao/astutil/pdm/pdmfitphase.x
new file mode 100644
index 00000000..6d9aa901
--- /dev/null
+++ b/noao/astutil/pdm/pdmfitphase.x
@@ -0,0 +1,43 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+# PDM_FITPHASE -- Call ICFIT on the Phase Curve.
+
+procedure pdm_fitphase (pdmp)
+
+pointer	pdmp			# pointer to PDM data structure
+
+double	xmin, xmax
+pointer	weights
+int	i, npt
+errchk	calloc, icg_fit
+
+begin
+	# Dereference some pointers.
+	npt = PDM_NPT(pdmp)
+
+	# Calloc a phase in-use array.
+	call alimd (PDM_XPH(pdmp,1), npt, xmin, xmax)
+	call ic_putr (PDM_ICP(pdmp), "xmin", real(xmin))
+	call ic_putr (PDM_ICP(pdmp), "xmax", real(xmax))
+	call calloc (weights, npt, TY_DOUBLE)
+
+	# Permute the data in-use array and save it in the phase weights array.
+	do i = 1, npt
+	    Memd[weights+i-1] = double(PDM_INUSE(pdmp,PDM_SORT(pdmp,i)))
+
+	# Call icfit on the phase curve (pass the phase x, y, wts)
+	if (npt >= 2)
+	    call icg_fitd (PDM_ICP(pdmp), PDM_GP(pdmp), "cursor",
+		PDM_GT(pdmp), PDM_CVP(pdmp), PDM_XPH(pdmp,1),
+		PDM_YPH(pdmp,1), Memd[weights], npt)
+
+	# Update the data in-use array appropriately.
+	do i = 1, npt
+	    PDM_INUSE(pdmp,PDM_SORT(pdmp,i)) = int(Memd[weights+i-1])
+
+	call mfree (weights, TY_DOUBLE)
+end
diff --git a/noao/astutil/pdm/pdmgdata.x b/noao/astutil/pdm/pdmgdata.x
new file mode 100644
index 00000000..79e8d49f
--- /dev/null
+++ b/noao/astutil/pdm/pdmgdata.x
@@ -0,0 +1,136 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+define	SZ_BUF	100
+
+# PDM_GDATA -- Get Data from the input files.
+
+int procedure pdm_gdata (pdmp, infile)
+
+pointer	pdmp			# pointer to PDM data structure
+char	infile[SZ_LINE]		# input data file name
+
+int	fntopnb(), list, clgfil()
+int	n, ncols, lineno, buflen
+int	open(), getline(), nscan()
+int	fd
+pointer	nextfile
+pointer	lbuf, ip, sp
+errchk	realloc, fntopnb, open
+
+begin
+	# Get a line buffer.
+	call smark (sp)
+	call salloc (lbuf, SZ_LINE, TY_CHAR)
+	call salloc (nextfile, SZ_LINE, TY_CHAR)
+
+	# Open the input file as a list of files.
+	list = fntopnb (infile, 0)
+
+	# Initialize some variables.
+	n = 0
+	ncols = 0
+	lineno = 0
+
+	# For each input file in the list, read the data.
+	while (clgfil (list, Memc[nextfile], SZ_FNAME) != EOF) {
+
+	    # Open this input file.
+	    fd = open (Memc[nextfile], READ_ONLY, TEXT_FILE)
+
+	    # Read in the data from this file.
+	    while (getline (fd, Memc[lbuf]) != EOF) {
+	        # Skip white space and blank lines.
+	        lineno = lineno + 1
+	        for (ip = lbuf; IS_WHITE(Memc[ip]); ip = ip + 1)
+		    ;
+	        if (Memc[ip] == '\n' || Memc[ip] == EOS)
+		    next
+
+	        if (n == 0) {
+		    buflen = SZ_BUF
+		    iferr {
+		        call calloc (PDM_XP(pdmp), buflen, TY_DOUBLE)
+		        call calloc (PDM_DYP(pdmp), buflen, TY_DOUBLE)
+		        call calloc (PDM_ODYP(pdmp), buflen, TY_DOUBLE)
+		        call calloc (PDM_INUSEP(pdmp), buflen, TY_INT)
+		        call calloc (PDM_ERRP(pdmp), buflen, TY_REAL)
+		    } then
+		        call erract (EA_FATAL)
+	        } else if (n + 1 > buflen) {
+		    buflen = buflen + SZ_BUF
+		    call realloc (PDM_XP(pdmp), buflen, TY_DOUBLE)
+		    call realloc (PDM_DYP(pdmp), buflen, TY_DOUBLE)
+		    call realloc (PDM_ODYP(pdmp), buflen, TY_DOUBLE)
+		    call realloc (PDM_INUSEP(pdmp), buflen, TY_INT)
+		    call realloc (PDM_ERRP(pdmp), buflen, TY_REAL)
+	        }
+
+		# Read data from the file, put it in the data structure.
+	        call sscan (Memc[ip])
+	        call gargd (PDM_X(pdmp,n+1))
+	        call gargd (PDM_ODY(pdmp,n+1))
+	        call gargr (PDM_ERR(pdmp,n+1))
+	        PDM_INUSE(pdmp,n+1) = 1
+	        PDM_DY(pdmp,n+1) = PDM_ODY(pdmp,n+1)
+
+	        # If this is line one, then determine the number of columns.
+	        if (ncols == 0 && nscan() > 0)
+		    ncols = nscan()
+
+	        # Check this line against the number of columns and do the
+	        # appropriate thing.
+
+	        switch (nscan()) {
+	        case 0:
+		    call printf ("no args; %s, line %d: %s\n")
+		        call pargstr (Memc[nextfile])
+		        call pargi (lineno)
+		        call pargstr (Memc[lbuf])
+		    next
+	        case 1:
+		    if (ncols >= 2) {
+		        call eprintf ("only one arg; %s, line %d: %s\n")
+		            call pargstr (Memc[nextfile])
+		            call pargi (lineno)
+		            call pargstr (Memc[lbuf])
+		        next
+		    } else {
+		        PDM_ODY(pdmp,n+1) = PDM_X(pdmp,n+1)
+		        PDM_DY(pdmp,n+1) = PDM_X(pdmp,n+1)
+		        PDM_X(pdmp,n+1) = n + 1.0d+0
+			PDM_ERR(pdmp,n+1) = 0.0
+		    }
+		case 2:
+		    if (ncols == 3) {
+			call eprintf ("only two args; %s, line %d: %s\n")
+			    call pargstr (Memc[nextfile])
+			    call pargi (lineno)
+			    call pargstr (Memc[lbuf])
+			next
+		    } else {
+		        PDM_ODY(pdmp,n+1) = PDM_ODY(pdmp,n+1)
+		        PDM_DY(pdmp,n+1) = PDM_ODY(pdmp,n+1)
+		        PDM_X(pdmp,n+1) = PDM_X(pdmp,n+1)
+			PDM_ERR(pdmp,n+1) = 0.0
+		    }
+
+	        }
+
+	        n = n + 1
+	    }
+	    call close (fd)
+	}
+
+	call realloc (PDM_XP(pdmp), n, TY_DOUBLE)
+	call realloc (PDM_DYP(pdmp), n, TY_DOUBLE)
+	call realloc (PDM_ODYP(pdmp), n, TY_DOUBLE)
+	call realloc (PDM_INUSEP(pdmp), n, TY_INT)
+	call realloc (PDM_ERRP(pdmp), n, TY_REAL)
+
+	call fntclsb (list)
+	call sfree (sp)
+	return (n)
+end
diff --git a/noao/astutil/pdm/pdmminmaxp.x b/noao/astutil/pdm/pdmminmaxp.x
new file mode 100644
index 00000000..01cdcd2c
--- /dev/null
+++ b/noao/astutil/pdm/pdmminmaxp.x
@@ -0,0 +1,43 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_MINMAXP -- Calculate the minimum and maximum periods automatically.
+
+procedure pdm_minmaxp (pdmp)
+
+pointer pdmp			# pointer to PDM data structure
+
+int	npt, i
+double	minp, maxp
+pointer	sep, sp
+
+begin
+	call smark (sp)
+	npt = PDM_NPT(pdmp)
+
+	# Allocate an array of separations and fill it.  Find the minimum
+	# separation as we go.
+
+	call salloc (sep, npt-1, TY_DOUBLE)
+	maxp = PDM_X(pdmp,npt) - PDM_X(pdmp,1)
+	minp = maxp
+	do i = 1, npt-1 {
+	    Memd[sep+i-1] = PDM_X(pdmp,i+1) - PDM_X(pdmp,i)
+	    if (Memd[sep+i-1] < minp)
+		minp = Memd[sep+i-1]
+	}
+
+	# Set minp equal to twice this minimum (Nyquist criterion).  Set fmax.
+	PDM_PMIN(pdmp) = 2.0d+0 * minp
+	if (PDM_PMIN(pdmp) != 0.0d+0)
+	    PDM_FMAX(pdmp) = 1.0d+0/PDM_PMIN(pdmp)
+
+	# Set maxp equal to 4 times maxp.  Set fmin.
+	PDM_PMAX(pdmp) = 4.0d+0 * maxp
+	if (PDM_PMAX(pdmp) != 0.0d+0)
+	    PDM_FMIN(pdmp) = 1.0d+0/PDM_PMAX(pdmp)
+
+	call sfree (sp)
+end
diff --git a/noao/astutil/pdm/pdmopen.x b/noao/astutil/pdm/pdmopen.x
new file mode 100644
index 00000000..aa58ad5f
--- /dev/null
+++ b/noao/astutil/pdm/pdmopen.x
@@ -0,0 +1,47 @@
+include <mach.h>
+include "pdm.h"
+
+# PDM_OPEN -- Open a new PDM data structure.
+
+int procedure pdm_open (device, batchfile, metafile, interactive)
+
+char	device[SZ_FNAME]		# graphics device
+char	batchfile[SZ_FNAME]		# file to store batch information
+char	metafile[SZ_FNAME]		# file to store plots
+bool	interactive			# interactive flag
+
+pointer	pdmp
+pointer	gt_init(), gopen()
+int	open()
+errchk	ic_open, gopen, open, calloc, malloc
+
+begin
+	# Allocate a pdm data structure.
+	call calloc (pdmp, PDM_LENSTRUCT, TY_STRUCT)
+
+	# Set up icfit structure pointers for data and phase curve fits.
+	call ic_open (PDM_ICD(pdmp)) 
+	call ic_open (PDM_ICP(pdmp)) 
+
+ 	# Set up gtools and the gio structure pointers.
+	PDM_GT(pdmp) = gt_init ()
+	if (interactive) {
+	    PDM_GP(pdmp) = gopen (device, NEW_FILE, STDGRAPH)
+	    PDM_LFD(pdmp) = open ("STDOUT", APPEND, TEXT_FILE)
+	} else {
+	    PDM_LFD(pdmp) = open (batchfile, APPEND, TEXT_FILE)
+	    PDM_PFD(pdmp) = open (metafile, APPEND, BINARY_FILE)
+	    PDM_GP(pdmp) = gopen ("stdvdm", NEW_FILE, PDM_PFD(pdmp))
+	}
+
+	# Allocate space for the sample string and put a '*' in it.
+	call malloc (PDM_SAMPLEP(pdmp), SZ_LINE, TY_CHAR)
+	call strcpy ("*", PDM_SAMPLE(pdmp), SZ_LINE)
+
+	# Booleans.
+	PDM_RESID(pdmp) = NO
+	PDM_RANGE(pdmp) = YES
+	PDM_EB(pdmp) = NO
+
+	return (pdmp)
+end
diff --git a/noao/astutil/pdm/pdmphase.x b/noao/astutil/pdm/pdmphase.x
new file mode 100644
index 00000000..76b3e0d2
--- /dev/null
+++ b/noao/astutil/pdm/pdmphase.x
@@ -0,0 +1,72 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+# PDM_PHASE -- Calculate Phase Curve for period (compressed light curve).
+
+procedure pdm_phase (pdmp, period, epoch)
+
+pointer	pdmp			# pointer to PDM data structure
+double	period			# period to calculate the phase for
+double	epoch			# epoch of this data
+
+int	j, offset, temp
+double	p
+pointer	npt, phaseint, sp
+errchk	calloc, realloc
+
+begin
+	call smark (sp)
+	npt = PDM_NPT(pdmp)
+	call salloc (phaseint, npt, TY_DOUBLE)
+
+	
+	# Allocate space for the output phase data (ordinate and abscissa)
+	# in the pdm data structure.
+
+	if (PDM_XPHP(pdmp) == NULL) {
+	    call calloc (PDM_XPHP(pdmp), 2*npt, TY_DOUBLE)
+	    call calloc (PDM_YPHP(pdmp), 2*npt, TY_DOUBLE)
+	    call calloc (PDM_PHERRP(pdmp), 2*npt, TY_REAL)
+	} else {
+	    call realloc (PDM_XPHP(pdmp), 2*npt, TY_DOUBLE)
+	    call realloc (PDM_YPHP(pdmp), 2*npt, TY_DOUBLE)
+	    call realloc (PDM_PHERRP(pdmp), 2*npt, TY_REAL)
+	}
+
+	# Set up the sort array and a temporary array for the phases.
+	if (PDM_SORTP(pdmp) == NULL)
+	    call calloc (PDM_SORTP(pdmp), npt, TY_INT)
+	else
+	    call realloc (PDM_SORTP(pdmp), npt, TY_INT)
+
+	# Calculate the phases for all the points.
+	for (j=1; j<=npt; j=j+1) {
+	    PDM_SORT(pdmp,j) = j
+	    if (period > EPSILONR) {
+		temp = (int(epoch/period)+1)*period
+	        p = (PDM_X(pdmp,j) - epoch + temp)/period
+	    }
+	    Memd[phaseint+j-1] = double(p - int(p))
+	}
+
+	# Sort the phase array into ascending order and permute
+	# the index array (sort).
+
+	call pdm_sort (phaseint, PDM_SORTP(pdmp), npt)
+
+	# Store the data in the pdm data structure.
+	do j = 1, npt {
+	    offset = PDM_SORT(pdmp,j)
+	    PDM_YPH(pdmp,j) = PDM_DY(pdmp,offset)
+	    PDM_XPH(pdmp,j) = Memd[phaseint+offset-1]
+	    PDM_PHERR(pdmp,j) = PDM_ERR(pdmp,offset)
+	    PDM_YPH(pdmp,j+npt) = PDM_DY(pdmp,offset)
+	    PDM_XPH(pdmp,j+npt) = Memd[phaseint+offset-1] + 1.0
+	    PDM_PHERR(pdmp,j+npt) = PDM_ERR(pdmp,offset)
+	}
+
+	call sfree (sp)
+end
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
diff --git a/noao/astutil/pdm/pdmranperm.x b/noao/astutil/pdm/pdmranperm.x
new file mode 100644
index 00000000..1d80264e
--- /dev/null
+++ b/noao/astutil/pdm/pdmranperm.x
@@ -0,0 +1,56 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include	<gset.h>
+include <pkg/rg.h>
+include "pdm.h"
+
+# PDM_RANPERM -- Make a random permutation of the data vector.
+# This is the algorithm:
+# Starting at the beginning of the input array,
+# Do this numdatapoints times {
+#   Find the next random position.
+#   Call urand for a random number between one and numdatapoints(hash).
+#   Move forward in the input array this number of places (mod numpts).
+#   While (marker array for this position has a zero) {
+#       Move to the next position (linear rehash).
+#   }
+#   Put the input array value found at this position into the
+#   output array, set the marker array value corresponding
+#   to this position to zero.
+# }
+
+
+procedure pdm_ranperm (inarray, inuse, outarray, outinuse, numpts, seed)
+
+int	numpts			# number of points in the data
+double	inarray[numpts]		# data to be permuted
+double	outarray[numpts]	# output permuted data to this array
+int	inuse[numpts]		# the PDM in-use array
+int	outinuse[numpts]	# output scrambled in-use array
+long	seed			# a seed for the random number generator
+
+int	count, pos
+real	urand()
+pointer	p, sp			# array to keep track of which have been used
+
+begin
+	# Allocate the output array and the marker array.
+	# Fill the marker array with ones (amovki)
+
+	call smark (sp)
+	call salloc (p, numpts, TY_INT)
+	call amovki (1, Memi[p], numpts)	# Set this array to all ones.
+
+	pos = 0
+	do count = 1, numpts {
+	    pos = mod(pos+int(urand(seed)*numpts)+1, numpts) # Hash.
+	    while (Memi[p+pos] == 0)			     # Linear rehash.
+		pos = mod(pos+1, numpts)
+	    outarray[count] = inarray[pos+1]
+	    outinuse[count] = inuse[pos+1]
+	    Memi[p+pos] = 0
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/astutil/pdm/pdmshow.x b/noao/astutil/pdm/pdmshow.x
new file mode 100644
index 00000000..6f8af37a
--- /dev/null
+++ b/noao/astutil/pdm/pdmshow.x
@@ -0,0 +1,56 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_SHOW -- Print information to file.
+
+procedure pdm_show (pdmp, file, verbose)
+
+pointer	pdmp			# pointer to PDM data structure
+char	file[ARB]		# file to put the show information in
+bool	verbose			# verbose output flag
+
+int	fd, open(), i
+errchk	open()
+
+begin
+	# Open the output file.
+	fd = open (file, APPEND, TEXT_FILE)
+
+	# Print information from the data structure.
+	call fprintf (fd,
+		    "minimum period searched = %12.12g, maximum = %g12.12\n")
+	    call pargd (PDM_PMIN(pdmp))
+	    call pargd (PDM_PMAX(pdmp))
+	call fprintf (fd,
+		    "period = %12.12g, amplitude = %12.12g, epoch = %12.12g\n")
+	    call pargd (PDM_MINR(pdmp))
+	    call pargd (PDM_AMPL(pdmp))
+	    call pargd (PDM_EPOCH(pdmp))
+
+	if (verbose) {
+	    # Print the working data set out as x,y,in-use triplets.
+	    call fprintf (fd, "The working data vector is as follows: \n")
+	    do i = 1, PDM_NPT(pdmp) {
+		call fprintf ( fd, "index = %d, x = %12.12g, y = %12.12g\n")
+		    call pargi (i)
+		    call pargd (PDM_X(pdmp,i))
+		    call pargd (PDM_DY(pdmp,i))
+	    }
+
+	    if (PDM_XPHP(pdmp) != NULL) {
+	        # Print the phasecurve out as x,y pairs
+	        call fprintf (fd, "\nThe phase curve vector is as follows: \n")
+	        do i = 1, PDM_NPT(pdmp) {
+		    call fprintf ( fd, "index = %d, x = %12.12g, y = %12.12g\n")
+		        call pargi (i)
+		        call pargd (PDM_XPH(pdmp,i))
+		        call pargd (PDM_YPH(pdmp,i))
+	        }
+	    }
+	}
+
+	# Close the output file.
+	call close (fd)
+end
diff --git a/noao/astutil/pdm/pdmsignif.x b/noao/astutil/pdm/pdmsignif.x
new file mode 100644
index 00000000..13a17acc
--- /dev/null
+++ b/noao/astutil/pdm/pdmsignif.x
@@ -0,0 +1,61 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include	<gset.h>
+include <pkg/rg.h>
+include "pdm.h"
+
+define	NUMTRIES	100
+
+# PDM_SIGNIF -- Calculate the significance of the theta statistic for
+# a certain period.  Use the "Method of Randomization".
+
+double procedure pdm_signif (pdmp, period)
+
+pointer	pdmp			# pointer to PDM data structure
+double	period			# period at which to find significance
+
+int	lesscount, i, npt
+double	otheta, theta, pdm_theta(), pdm_thetaran()
+long	seed
+pointer rg, pt, inuse, oinuse, rg_xrangesd(), sp
+errchk	pdm_statistics, pdm_theta(), pdm_ranperm, pdm_thetaran()
+
+begin
+        # Do NUMTRIES random permutations on the data, calculate the theta
+        # statistic on the scrambled data for this period.  Return the
+        # fraction of these permutaions which yield thetas less than
+        # the unmixed data.
+
+	call smark (sp)
+	npt = PDM_NPT(pdmp)
+
+	# Make sure the statistics are up to date.
+	call pdm_statistics (pdmp)
+
+	# Allocate a temporary array for the scrambled data, allocate a 
+	# temporary copy of the inuse array and copy the real inuse array
+	# into it.
+
+	call salloc (pt, npt, TY_DOUBLE)
+	call salloc (inuse, npt, TY_INT)
+	call salloc (oinuse, npt, TY_INT)
+	call amovi (PDM_INUSE(pdmp,1), Memi[inuse], npt)
+	lesscount = 0
+
+	# Calculate the ranges information from the sample string.
+	rg = rg_xrangesd (PDM_SAMPLE(pdmp), PDM_X(pdmp,1), npt)
+	otheta = pdm_theta (pdmp, rg, period)
+	seed = 1.0
+
+	do i = 1, NUMTRIES {
+	    call pdm_ranperm (PDM_DY(pdmp,1), Memi[inuse], Memd[pt],
+		Memi[oinuse], npt, seed)
+	    theta = pdm_thetaran (pdmp, pt, oinuse, rg, period)
+	    if (theta < otheta)
+		lesscount = lesscount + 1
+	}
+
+	call sfree (sp)
+	return (1.0d+0 - (double(lesscount)/double(NUMTRIES)))
+end
diff --git a/noao/astutil/pdm/pdmsort.x b/noao/astutil/pdm/pdmsort.x
new file mode 100644
index 00000000..93c59aec
--- /dev/null
+++ b/noao/astutil/pdm/pdmsort.x
@@ -0,0 +1,20 @@
+include <mach.h>
+include "pdm.h"
+
+# PDM_SORT -- Sort the phases into ascending order.
+
+procedure pdm_sort (array, sort, numpts)
+
+pointer	array			# array to sort
+pointer	sort			# array to contain the sort indexes
+int	numpts			# number of points to sort
+
+pointer	comarray
+common /sortcom/ comarray
+extern	pdm_compare()
+
+begin
+	comarray = array
+	if (numpts > 1)
+	    call qsort (Memi[sort], numpts, pdm_compare)
+end
diff --git a/noao/astutil/pdm/pdmstats.x b/noao/astutil/pdm/pdmstats.x
new file mode 100644
index 00000000..aeabacaa
--- /dev/null
+++ b/noao/astutil/pdm/pdmstats.x
@@ -0,0 +1,37 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include "pdm.h"
+
+# PDM_STATISTICS -- Calculate the sum of squares and the variance of the data.
+
+procedure pdm_statistics (pdmp)
+
+pointer	pdmp			# pointer to PDM data structure
+
+int	npt, i, j
+double	var, sumx2, sum
+
+begin
+	npt = PDM_NPT(pdmp)
+
+	# Calculate the sum of squares and the variance of the data.
+	sumx2 = 0.0
+	sum = 0.0
+	j = 0
+
+	do i = 1, npt {
+	    if (PDM_INUSE(pdmp,i) == 1) {
+	        sumx2 = sumx2 + PDM_DY(pdmp,i)**2	    # Sum of squares.
+	        sum = sum + PDM_DY(pdmp,i)
+		j = j + 1
+	    }
+	}
+
+	if (j != 1)
+	    var = (sumx2 - sum**2/double(j))/double(j - 1)    # Variance.
+
+	# Put these two values in the data structure.
+	PDM_SUMSQ(pdmp) = sumx2
+	PDM_DVAR(pdmp) = var
+end
diff --git a/noao/astutil/pdm/pdmtheta.x b/noao/astutil/pdm/pdmtheta.x
new file mode 100644
index 00000000..b2b295bc
--- /dev/null
+++ b/noao/astutil/pdm/pdmtheta.x
@@ -0,0 +1,120 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <pkg/rg.h>
+include "pdm.h"
+
+# PDM_THETA -- Calculate the theta statistic for this period.
+# Theta is a measure of the dispersion of data phases about a mean
+# light curve.
+
+double procedure pdm_theta (pdmp, rg, period)
+
+pointer	pdmp			# pointer to PDM data structure
+pointer	rg			# segments structure pointer
+double	period			# period at which to find theta
+
+int	i, j, k, l, m
+double	s2
+int	ndof, segst, segend
+bool	bin10
+double	theta
+pointer	sumbin, sum2bin, numbin, sp
+errchk	binem
+
+begin
+	# Allocate bin storage.
+	call smark (sp)
+	call salloc (sumbin, 10, TY_DOUBLE)
+	call salloc (sum2bin, 10, TY_DOUBLE)
+	call salloc (numbin, 10, TY_INT)
+
+	ndof = 0
+	s2 = 0
+
+	# Do loop on the segments.
+	do i = 1, RG_NRGS(rg) {
+
+	    # Calculate segst, segend, bin10.
+	    segst = min(RG_X2(rg,i),RG_X1(rg,i))
+	    segend = max(RG_X2(rg,i),RG_X1(rg,i))
+	    bin10 = ((segend - segst) >= BIN10)
+
+	    # Calculate the number of points in each bin and the sum of 
+	    # the bins.
+
+	    call binem (period, bin10, PDM_XP(pdmp), PDM_DYP(pdmp),
+		segst, segend, PDM_INUSEP(pdmp), sumbin, sum2bin, numbin)
+
+	    # Calculate sigma**2 for this period.
+	    for (j=0; j<=9; j=j+1) {
+	        k = numbin+j
+	        l = sumbin+j
+	        m = sum2bin+j
+	        if (Memi[k] > 1) {
+		    ndof = ndof + Memi[k] - 1
+		    s2 = s2 +  (Memd[m] - Memd[l]**2 / Memi[k])
+	        }
+	    }
+	}
+
+	# Calculate theta.
+	theta = (s2 / double(ndof)) / PDM_DVAR(pdmp)
+
+	call sfree (sp)
+
+	return (theta)
+
+end
+
+
+# BINEM -- Put the data points into the appropriate bins.
+
+procedure binem(incper, bin10, x, y, segst, segend, inuse, sumbin, sum2bin,
+	numbin)
+
+double	incper			# period increment
+bool	bin10			# use 5 or 10 bins flag
+pointer	x			# ordinates
+pointer	y			# abcissas
+pointer	inuse			# PDM in-use array
+int	segst, segend		# segment start and segment end
+pointer	sumbin, sum2bin, numbin	# pointers to bins of sum, sum2, and number
+
+int	bin1, bin2, j, k, l, m
+double	p, phase, p0
+
+begin
+	do j = 1, 10 {
+	    Memi[numbin+j-1] = 0
+	    Memd[sumbin+j-1] = 0.0
+	    Memd[sum2bin+j-1] = 0.0
+	}
+
+	#p0 = Memd[x]
+	call alimd (Memd[x+segst-1], segend-segst+1, p0, p)
+	do j = segst, segend {
+	    if (Memi[inuse+j-1] == 0)
+		next
+	    p = (Memd[x+j-1] - p0)/incper
+	    phase = double(p - int(p))
+	    if (bin10) {
+		bin1 = mod(int(10.*phase+0.5), 10)
+	    } else {
+		bin1 = 2 * int(5. * phase) + 1
+		bin2 = 2 * (mod(int(5. * phase + 0.5), 5))
+		k = numbin+bin2
+		l = sumbin+bin2
+		m = sum2bin+bin2
+		Memi[k] = Memi[k] + 1
+		Memd[l] = Memd[l] + Memd[y+j-1]
+		Memd[m] = Memd[m] + Memd[y+j-1] ** 2
+	    }
+	    k = numbin+bin1
+	    l = sumbin+bin1
+	    m = sum2bin+bin1
+	    Memi[k] = Memi[k] + 1
+	    Memd[l] = Memd[l] + Memd[y+j-1]
+	    Memd[m] = Memd[m] + Memd[y+j-1] ** 2
+	}
+end
diff --git a/noao/astutil/pdm/pdmthetaran.x b/noao/astutil/pdm/pdmthetaran.x
new file mode 100644
index 00000000..28e63b34
--- /dev/null
+++ b/noao/astutil/pdm/pdmthetaran.x
@@ -0,0 +1,118 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <pkg/rg.h>
+include "pdm.h"
+
+# PDMTHETARAN -- This program is a copy of pdmtheta but can be used on
+# scrambled data.
+
+double procedure pdm_thetaran (pdmp, y, inuse, rg, period)
+
+pointer	pdmp			# pointer to PDM data structure
+pointer	y			# pointer to abcissas
+pointer	inuse			# pointer to PDM in-use array 
+pointer	rg			# pointer to ranges structure
+double	period			# period to calculate theta for
+
+int	i, j, k, l, m
+double	s2
+int	ndof, segst, segend
+bool	bin10
+double	theta
+pointer	sumbin, sum2bin, numbin, sp
+errchk	binemran
+
+begin
+	# Allocate bin storage.
+	call smark (sp)
+	call salloc (sumbin, 10, TY_DOUBLE)
+	call salloc (sum2bin, 10, TY_DOUBLE)
+	call salloc (numbin, 10, TY_INT)
+
+	s2 = 0
+	ndof = 0
+
+	# Do loop on the segments.
+	do i = 1, RG_NRGS(rg) {
+
+	    # Calculate segst, segend, bin10.
+	    segst = min(RG_X2(rg,i),RG_X1(rg,i))
+	    segend = max(RG_X2(rg,i),RG_X1(rg,i))
+	    bin10 = ((segend - segst) >= BIN10)
+
+	    # Calculate the number of points in each bin and the sum of 
+	    # the bins.
+
+	    call binemran (period, bin10, PDM_XP(pdmp), y, segst,
+		segend, inuse, sumbin, sum2bin, numbin)
+
+	    # Calculate sigma**2 for this period.
+	    for (j=0; j<=9; j=j+1) {
+	        k = numbin+j
+	        l = sumbin+j
+	        m = sum2bin+j
+	        if (Memi[k] > 1) {
+		    ndof = ndof + Memi[k] - 1
+		    s2 = s2 +  (Memd[m] - Memd[l]**2 / Memi[k])
+	        }
+	    }
+	}
+
+	# Calculate theta.
+	theta = (s2 / double(ndof)) / PDM_DVAR(pdmp)
+
+	call sfree (sp)
+	return (theta)
+end
+
+
+# BINEMRAN -- Put the data points into the appropriate bins (scrambled data).
+
+procedure binemran (incper, bin10, x, y, segst, segend, inuse, sumbin, sum2bin,
+	numbin)
+
+double	incper
+bool	bin10
+pointer	x
+pointer	y
+pointer	inuse
+int	segst, segend
+pointer	sumbin, sum2bin, numbin
+
+int	bin1, bin2, j, k, l, m
+double	p, phase, p0
+
+begin
+	do j = 1, 10 {
+	    Memi[numbin+j-1] = 0
+	    Memd[sumbin+j-1] = 0.0
+	    Memd[sum2bin+j-1] = 0.0
+	}
+
+	p0 = Memd[x]
+	do j = segst, segend {
+	    if (Memi[inuse+j-1] == 0)
+		next
+	    p = (Memd[x+j-1] - p0)/incper
+	    phase = double(p - int(p))
+	    if (bin10) {
+		bin1 = mod(int(10.*phase+0.5d+0), 10)
+	    } else {
+		bin1 = 2 * int(5.0d+0 * phase) + 1
+		bin2 = 2 * (mod(int(5.0d+0 * phase + 0.5d+0), 5))
+		k = numbin+bin2
+		l = sumbin+bin2
+		m = sum2bin+bin2
+		Memi[k] = Memi[k] + 1
+		Memd[l] = Memd[l] + Memd[y+j-1]
+		Memd[m] = Memd[m] + Memd[y+j-1] ** 2
+	    }
+	    k = numbin+bin1
+	    l = sumbin+bin1
+	    m = sum2bin+bin1
+	    Memi[k] = Memi[k] + 1
+	    Memd[l] = Memd[l] + Memd[y+j-1]
+	    Memd[m] = Memd[m] + Memd[y+j-1] ** 2
+	}
+end
diff --git a/noao/astutil/pdm/pdmtplot.x b/noao/astutil/pdm/pdmtplot.x
new file mode 100644
index 00000000..e8c3ce8a
--- /dev/null
+++ b/noao/astutil/pdm/pdmtplot.x
@@ -0,0 +1,139 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+define	EDGESCL	5.0			# Percent extra space around plot data
+
+# PDM_TPLOT -- Plot the data on the screen.
+
+procedure pdm_tplot (pdmp, porf, filename)
+
+pointer	pdmp			# pointer to PDM data structure
+int	porf			# period or frequency flag
+char	filename[SZ_FNAME]	# name of the input data file
+
+int	nthpt
+real	x1, x2, y1, y2, scldif, sclspc
+pointer	gp, xtemp, ytemp
+pointer	title
+char	system_id[SZ_LINE]
+int	indx, pdm_findmin()
+errchk	malloc, pdm_findmin()
+
+begin
+	# Dereference some structure stuff.
+	nthpt = PDM_NTHPT(pdmp)
+	gp = PDM_GP(pdmp)
+
+	call malloc (title, PDM_SZ_TITLE, TY_CHAR)
+	call malloc (xtemp, nthpt, TY_REAL)
+	call malloc (ytemp, nthpt, TY_REAL)
+	call gclear (gp)
+
+	do indx = 1, nthpt {
+	    Memr[xtemp+indx-1] = real(PDM_XTH(pdmp,indx))
+	    Memr[ytemp+indx-1] = real(PDM_YTH(pdmp,indx))
+	}
+
+	if (porf == THETAPPLOT) {
+	    # Scale the wcs.
+	    call gascale (gp, Memr[xtemp], nthpt, 1)
+	    call gascale (gp, Memr[ytemp], nthpt, 2)
+
+	    # Get the X and Y boundary values.
+	    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
+
+	    call gswind (gp, x1, x2, y1, y2)
+
+	    # Find the minimum value and save it in the remembered minimum.
+	    indx = pdm_findmin (pdmp, THETAPPLOT, PDM_PMIN(pdmp),
+		PDM_PMAX(pdmp), 1, nthpt)
+	    PDM_MINR(pdmp) = PDM_XTH(pdmp,indx)
+
+	    # Multiline title, save in an array and sprintf to it.
+	    # Get the system identification.
+
+	    call sysid (system_id, SZ_LINE)
+	    call sprintf (Memc[title], PDM_SZ_TITLE,
+		"%s\nFile = %s, minimum = %12.12g\n%s\nnumpts = %d")
+		call pargstr (system_id)
+	        call pargstr (filename)
+		call pargd (PDM_MINR(pdmp))
+	        call pargstr ("Theta vs Period")
+		call pargi (nthpt)
+
+	    # Draw the axes.
+	    call glabax (gp, Memc[title], "period", "theta")
+
+	    # Make the plot.
+	    call gpline (gp, Memr[xtemp], Memr[ytemp], nthpt)
+
+	    # Put the cursor at the minimum.
+	    call gscur (gp, Memr[xtemp+indx-1], Memr[ytemp+indx-1])
+	} else {
+	    # Scale the wcs.
+	    call gascale (gp, Memr[xtemp], nthpt, 1)
+	    call gascale (gp, Memr[ytemp], nthpt, 2)
+
+	    # Get the X and Y boundary values.
+	    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
+
+	    call gswind (gp, x1, x2, y1, y2)
+
+	    # Find the minimum value and save it in the remembered minimum.
+	    indx = pdm_findmin (pdmp, THETAFPLOT, PDM_FMIN(pdmp),
+		PDM_FMAX(pdmp), 1, nthpt)
+	    if (PDM_XTH(pdmp,indx) > EPSILOND)
+	        PDM_MINR(pdmp) = 1./PDM_XTH(pdmp,indx)
+
+	    # Multiline title, save in an array and sprintf to it.
+	    # Get the system identification.
+
+	    call sysid (system_id, SZ_LINE)
+	    call sprintf (Memc[title], PDM_SZ_TITLE,
+	        "%s\nFile = %s, minimum = %12.12g\n%s\nnumpts = %d")
+		call pargstr (system_id)
+	        call pargstr (filename)
+		if (PDM_MINR(pdmp) > EPSILOND)
+		    call pargd (1.0d+0/PDM_MINR(pdmp))
+		else
+		    call pargd (0.0d+0)
+	        call pargstr ("Theta vs Frequency")
+		call pargi (nthpt)
+
+	    # Draw the axes.
+	    call glabax (gp, Memc[title], "frequency", "theta")
+
+	    # Make the plot.
+	    call gpline (gp, Memr[xtemp], Memr[ytemp], nthpt)
+
+	    # Put the cursor at the minimum.
+	    call gscur (gp, Memr[xtemp+indx-1], Memr[ytemp+indx-1])
+        }
+
+	call mfree (title, TY_CHAR)
+	call mfree (xtemp, TY_REAL)
+	call mfree (ytemp, TY_REAL)
+end
diff --git a/noao/astutil/pdm/pdmundelete.x b/noao/astutil/pdm/pdmundelete.x
new file mode 100644
index 00000000..2f5f0e31
--- /dev/null
+++ b/noao/astutil/pdm/pdmundelete.x
@@ -0,0 +1,124 @@
+include <mach.h>
+include <ctype.h>
+include <error.h>
+include <gset.h>
+include "pdm.h"
+
+define	MSIZE		2.0		# Mark size
+
+# PDM_UNDELETE -- Undelete the nearest deleted point from the plot
+# and set it's inuse array entry to one (in-use).
+
+int procedure pdm_undelete (pdmp, cx, cy, ptype)
+
+pointer	pdmp			# pointer to PDM data structure
+double	cx, cy			# device coordinates of point to undelete
+int	ptype			# plot type flag
+
+pointer	gp
+real	x, y
+int	npt, i, j, index
+real	x0, y0, r2, r2min
+
+begin
+	gp = PDM_GP(pdmp)
+	npt = PDM_NPT(pdmp)
+
+	if (ptype == DATAPLOT) {
+	    # Transform world cursor coordinates to NDC.
+	    call gctran (gp, real(cx), real(cy), x, y, 1, 0)
+
+	    # Search for nearest point not in-use.
+	    j = 0
+	    r2min = MAX_REAL
+	    do i = 1, npt {
+	        if (PDM_INUSE(pdmp,i) == 1)
+		    next
+
+	        call gctran (gp, real(PDM_X(pdmp,i)), real(PDM_DY(pdmp,i)),
+			     x0, y0, 1, 0)
+		
+	        r2 = (x0 - x) ** 2 + (y0 - y) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+
+	    # Unmark the deleted point and reset the weight.
+	    if (j != 0) {
+	        call gscur (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)))
+	        call gseti (gp, G_PMLTYPE, GL_CLEAR)
+	        call gmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)),
+			    GM_CROSS, MSIZE, MSIZE)
+	        call gseti (gp, G_PMLTYPE, GL_SOLID)
+		if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp))
+	            call gpmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)),
+		        1, GM_PLUS, MSIZE, MSIZE)
+		else
+	            call gpmark (gp, real(PDM_X(pdmp,j)), real(PDM_DY(pdmp,j)),
+		        1, GM_POINT, 1.0, 1.0)
+	        PDM_INUSE(pdmp,j) = 1
+		call gflush (gp)
+	    }
+
+	    return (j)
+	} else if (ptype == PHASEPLOT) {
+	    # Transform world cursor coordinates to NDC.
+	    call gctran (gp, real(cx), real(cy), x, y, 1, 0)
+
+	    # Search for nearest point not in-use.
+	    j = 0
+	    r2min = MAX_REAL
+	    do i = 1, npt {
+		index = PDM_SORT(pdmp,i)
+	        if (PDM_INUSE(pdmp,index) == 1)
+		    next
+
+	        call gctran (gp, real(PDM_XPH(pdmp,i)), real(PDM_YPH(pdmp,i)),
+			     x0, y0, 1, 0)
+		
+	        r2 = (x0 - x) ** 2 + (y0 - y) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    j = i
+	        }
+	    }
+
+	    # Unmark the deleted point and reset the weight.
+	    if (j != 0) {
+	        call gscur (gp, real(PDM_XPH(pdmp,j))+1.0,
+		    real(PDM_YPH(pdmp,j)))
+	        call gseti (gp, G_PMLTYPE, GL_CLEAR)
+	        call gmark (gp, real(PDM_XPH(pdmp,j))+1.0,
+		    real(PDM_YPH(pdmp,j)), GM_CROSS, MSIZE, MSIZE)
+	        call gseti (gp, G_PMLTYPE, GL_SOLID)
+		if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp))
+	            call gpmark (gp, real(PDM_XPH(pdmp,j))+1.0,
+		    real(PDM_YPH(pdmp,j)), 1, GM_PLUS, MSIZE, MSIZE)
+		else
+	            call gpmark (gp, real(PDM_XPH(pdmp,j))+1.0,
+		    real(PDM_YPH(pdmp,j)), 1, GM_POINT, 1.0, 1.0)
+
+	        call gscur (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)))
+	        call gseti (gp, G_PMLTYPE, GL_CLEAR)
+	        call gmark (gp, real(PDM_XPH(pdmp,j)), real(PDM_YPH(pdmp,j)),
+			    GM_CROSS, MSIZE, MSIZE)
+	        call gseti (gp, G_PMLTYPE, GL_SOLID)
+		if (PDM_NPT(pdmp) <= PDM_PLUSPOINT(pdmp))
+	            call gpmark (gp, real(PDM_XPH(pdmp,j)),
+		    real(PDM_YPH(pdmp,j)), 1, GM_PLUS, MSIZE, MSIZE)
+		else
+	            call gpmark (gp, real(PDM_XPH(pdmp,j)),
+		    real(PDM_YPH(pdmp,j)), 1, GM_POINT, 1.0, 1.0)
+		
+		# Calculate which point this corresponds to.
+		index = PDM_SORT(pdmp,j)
+	        PDM_INUSE(pdmp,index) = 1
+		call gflush (gp)
+	    }
+
+	    return (index)
+	} else
+	    return (0)
+end
diff --git a/noao/astutil/pdm/t_pdm.x b/noao/astutil/pdm/t_pdm.x
new file mode 100644
index 00000000..7498804a
--- /dev/null
+++ b/noao/astutil/pdm/t_pdm.x
@@ -0,0 +1,72 @@
+
+include <ctype.h>
+include <error.h>
+include <mach.h>
+include "pdm.h"
+
+# PDM -- Phase Dispersion Minimization.  Find periodicities in light
+# curve data. Root program.
+
+procedure t_pdm()
+
+char	infile[SZ_FNAME]		# input data file
+char	metafile[SZ_FNAME]		# batch metafile
+char	batchfile[SZ_FNAME]		# batch log file
+char	device[SZ_FNAME]		# plot device
+bool	interactive			# interactive or batch flag
+bool	autoranges			# autoranges flag
+
+int	ptype				# plot type
+pointer	pdmp, pdm_open()		# structure stuff
+bool	flip, clgetb()
+real	clgetr()
+int	clgeti(), pdm_gdata(), pdm_autorang()
+errchk	pdm_open, pdm_autorang, pdm_batch, pdm_cursor
+
+begin
+	# Get some of the CL parameters and open the pdm data structure.
+	call clgstr ("infiles", infile, SZ_FNAME)
+	call clgstr ("metafile", metafile, SZ_FNAME)
+	call clgstr ("batchfile", batchfile, SZ_FNAME)
+	call clgstr ("device", device, SZ_FNAME)
+	interactive = clgetb ("interactive")
+	flip = clgetb ("flip")
+	pdmp = pdm_open (device, batchfile, metafile, interactive)
+
+	# Clio to get more parameters.
+	PDM_PMIN(pdmp) = double(clgetr ("minp"))
+	call clputr ("minp", real(PDM_PMIN(pdmp)))
+	if (PDM_PMIN(pdmp) > EPSILOND)
+	    PDM_FMAX(pdmp) = 1.0d+0/PDM_PMIN(pdmp)
+	else
+	    PDM_FMAX(pdmp) = 0.0d+0
+	PDM_PMAX(pdmp) = double(clgetr ("maxp"))
+	call clputr ("maxp", real(PDM_PMAX(pdmp)))
+	if (PDM_PMAX(pdmp) > EPSILOND)
+	    PDM_FMIN(pdmp) = 1.0d+0/PDM_PMAX(pdmp)
+	else
+	    PDM_FMIN(pdmp) = 0.0d+0
+	PDM_NTHPT(pdmp) = clgeti ("ntheta")
+	autoranges = clgetb ("autoranges")
+	PDM_NSIGMA(pdmp) = double(clgetr ("nsigma"))
+	PDM_PLUSPOINT(pdmp) = clgeti ("pluspoint")
+
+	# Read in the data.
+	PDM_NPT(pdmp) = pdm_gdata (pdmp, infile)
+
+	# If the autoranges flag is set, call the autorange subroutine.
+	if (autoranges)
+	    PDM_NRANGE(pdmp) = pdm_autorang (pdmp)
+
+	if (!interactive)
+	    call pdm_batch (pdmp, batchfile, infile, flip)
+	else {
+	    # Plot the data on the screen and call the cursor loop. 
+	    call pdm_dplot (pdmp, infile, flip)
+	    ptype = DATAPLOT
+	    call pdm_cursor(pdmp, ptype, infile, flip)
+	}
+
+	# Close the pdm data structure.
+	call pdm_close (pdmp, interactive)
+end