Initial commit

1 files changed, 100 insertions, 0 deletions

diff --git a/math/bevington/polfit.f b/math/bevington/polfit.f
new file mode 100644
index 00000000..cfbdb178
--- /dev/null
+++ b/math/bevington/polfit.f
@@ -0,0 +1,100 @@
+C SUBROUTINE POLFIT.F 
+C
+C SOURCE
+C   BEVINGTON, PAGES 140-142.
+C
+C PURPOSE
+C   MAKE A LEAST-SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
+C      Y = A(1) + A(2)*X + A(3)*X**2 + A(3)*X**3 + . . .
+C
+C USAGE 
+C   CALL POLFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, CHISQR)
+C
+C DESCRIPTION OF PARAMETERS
+C   X	   - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
+C   Y	   - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
+C   SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS
+C   NPTS   - NUMBER OF PAIRS OF DATA POINTS
+C   NTERMS - NUMBER OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
+C   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
+C	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
+C	      0 (NO WEIGHTING) WEIGHT(I) = 1.
+C	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
+C   A	   - ARRAY OF COEFFICIENTS OF POLYNOMIAL
+C   CHISQR - REDUCED CHI SQUARE FOR FIT 
+C
+C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
+C   DETERM (ARRAY, NORDER)
+C      EVALUATES THE DETERMINANT OF A SYMETRICAL
+C      TWO-DIMENSIONAL MATRIX OF ORDER NORDER
+C
+C COMMENTS
+C   DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
+C
+	SUBROUTINE POLFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)
+	DOUBLE PRECISION SUMX,SUMY,XTERM,YTERM,ARRAY,CHISQ
+	DIMENSION X(1),Y(1),SIGMAY(1),A(1)
+	DIMENSION SUMX(19),SUMY(10),ARRAY(10,10)
+C
+C ACCUMULATE WEIGHTED SUMS
+C
+11	NMAX=2*NTERMS-1 
+	DO 13 N=1,NMAX
+13	SUMX(N)=0.
+	DO 15 J=1,NTERMS
+15	SUMY(J)=0.
+	CHISQ=0.
+21	DO 50 I=1,NPTS
+	XI=X(I) 
+	YI=Y(I) 
+31	IF (MODE) 32,37,39
+32	IF (YI) 35,37,33
+33	WEIGHT=1./YI
+	GOTO 41 
+35	WEIGHT=1./(-YI) 
+	GOTO 41 
+37	WEIGHT=1.
+	GOTO 41 
+39	WEIGHT=1./SIGMAY(I)**2
+41	XTERM=WEIGHT
+	DO 44 N=1,NMAX
+	SUMX(N)=SUMX(N)+XTERM
+44	XTERM=XTERM*XI
+45	YTERM=WEIGHT*YI 
+	DO 48 N=1,NTERMS
+	SUMY(N)=SUMY(N)+YTERM
+48	YTERM=YTERM*XI
+49	CHISQ=CHISQ+WEIGHT*YI**2
+50	CONTINUE
+C
+C CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS 
+C
+51	DO 54 J=1,NTERMS
+	DO 54 K=1,NTERMS
+	N=J+K-1 
+54	ARRAY(J,K)=SUMX(N)
+	DELTA=DETERM (ARRAY,NTERMS)
+	IF (DELTA) 61,57,61
+57	CHISQR=0.
+	DO 59 J=1,NTERMS
+59	A(J)=0. 
+	GOTO 80 
+61	DO 70 L=1,NTERMS
+62	DO 66 J=1,NTERMS
+	DO 65 K=1,NTERMS
+	N=J+K-1 
+65	ARRAY(J,K)=SUMX(N)
+66	ARRAY(J,L)=SUMY(J)
+70	A(L)=DETERM (ARRAY,NTERMS) /DELTA
+C
+C CALCULATE CHI SQUARE
+C
+71	DO 75 J=1,NTERMS
+	CHISQ=CHISQ-2.*A(J)*SUMY(J)
+	DO 75 K=1,NTERMS
+	N=J+K-1 
+75	CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
+76	FREE=NPTS-NTERMS
+77	CHISQR=CHISQ/FREE
+80	RETURN
+	END