Initial commit

1 files changed, 72 insertions, 0 deletions

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