1 files changed, 118 insertions, 0 deletions
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