1 files changed, 62 insertions, 0 deletions
diff --git a/math/bevington/pchisq.f b/math/bevington/pchisq.f
new file mode 100644
index 00000000..0dbe6d9b
--- /dev/null
+++ b/math/bevington/pchisq.f
@@ -0,0 +1,62 @@
+c function pchisq.f
+c
+c source
+c   Bevington, pages 192-193.
+c
+c purpose
+c   evaluate probability for exceeding chi square
+c
+c usage
+c   result = pchisq (chisqr, nfree)
+c
+c description of parameters
+c   chisqr - comparison value of reduced chi square
+c   nfree  - number of degrees of freedom
+c
+c subroutines and function subprograms required
+c   gamma (x)
+c      calculates gamma function for integers and half-integers
+c
+c comments
+c   calculation is approximate for nfree odd and
+c      chi square greater than 50
+c
+	function pchisq (chisqr,nfree)
+	double precision z,term,sum
+11	if (nfree) 12,12,14
+12	pchisq=0.
+	goto 60
+14	free=nfree
+	z=chisqr*free/2.
+	neven=2*(nfree/2)
+	if (nfree-neven) 21,21,41
+c
+c number of degrees of freedom is even
+c
+21	imax=nfree/2
+	term=1.
+	sum=0.
+31	do 34 i=1,imax
+	fi=i
+	sum=sum+term
+34	term=term*z/fi
+35	pchisq=sum*dexp(z)
+	goto 60
+c
+c number of degrees of freedom is odd
+c
+41	if (z-25) 44,44,42
+42	z=chisqr*(free-1.)/2.
+	goto 21
+44	pwr=free/2.
+	term=1.
+	sum=term/pwr
+51	do 56 i=1,1000
+	fi=i
+	term=-term*z/fi
+	sum=sum+term/(pwr+fi)
+55	if (dabs(term/sum)-0.00001) 57,57,56
+56	continue
+57	pchisq=1.-(z**pwr)*sum/gamma(pwr)
+60	return
+	end