Initial commit

1 files changed, 85 insertions, 0 deletions

diff --git a/math/bevington/interp.f b/math/bevington/interp.f
new file mode 100644
index 00000000..75ca1b23
--- /dev/null
+++ b/math/bevington/interp.f
@@ -0,0 +1,85 @@
+C SUBROUTINE INTERP.F 
+C
+C SOURCE
+C   BEVINGTON, PAGES 266-267.
+C
+C PURPOSE
+C   INTERPOLATE BETWEEN DATA POINTS TO EVALUATE A FUNCTION
+C
+C USAGE 
+C   CALL INTERP (X, Y, NPTS, NTERMS, XIN, YOUT) 
+C
+C DESCRIPTION OF PARAMETERS
+C   X	   - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
+C   Y	   - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
+C   NPTS   - NUMBER OF PAIRS OF DATA POINTS
+C   NTERMS - NUMBER OF TERMS IN FITTING POLYNOMIAL
+C   XIN    - INPUT VALUE OF X
+C   YOUT   - INTERPOLATED VALUE OF Y
+C
+C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
+C   NONE
+C
+C COMMENTS
+C   DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
+C   VALUE OF NTERMS MAY BE MODIFIED BY THE PROGRAM
+C
+	SUBROUTINE INTERP (X,Y,NPTS,NTERMS,XIN,YOUT)
+	DOUBLE PRECISION DELTAX,DELTA,A,PROD,SUM
+	DIMENSION X(1),Y(1)
+	DIMENSION DELTA(10),A(10)
+C
+C SEARCH FOR APPROPRIATE VALUE OF X(1)
+C
+11	DO 19 I=1,NPTS
+	IF (XIN-X(I)) 13,17,19
+13	I1=I-NTERMS/2
+	IF (I1) 15,15,21
+15	I1=1
+	GOTO 21 
+17	YOUT=Y(I)
+18	GOTO 61 
+19	CONTINUE
+	I1=NPTS-NTERMS+1
+21	I2=I1+NTERMS-1
+	IF (NPTS-I2) 23,31,31
+23	I2=NPTS 
+	I1=I2-NTERMS+1
+25	IF (I1) 26,26,31
+26	I1=1
+27	NTERMS=I2-I1+1
+C
+C EVALUATE DEVIATIONS DELTA
+C
+31	DENOM=X(I1+1)-X(I1)
+	DELTAX=(XIN-X(I1))/DENOM
+	DO 35 I=1,NTERMS
+	IX=I1+I-1
+35	DELTA(I)=(X(IX)-X(I1))/DENOM
+C
+C ACCUMULATE COEFFICIENTS A
+C
+40	A(1)=Y(I1)
+41	DO 50 K=2,NTERMS
+	PROD=1. 
+	SUM=0.
+	IMAX=K-1
+	IXMAX=I1+IMAX
+	DO 49 I=1,IMAX
+	J=K-I
+	PROD=PROD*(DELTA(K)-DELTA(J))
+49	SUM=SUM-A(J)/PROD
+50	A(K)=SUM+Y(IXMAX)/PROD
+C
+C ACCUMULATE SUM OF EXPANSION
+C
+51	SUM=A(1)
+	DO 57 J=2,NTERMS
+	PROD=1. 
+	IMAX=J-1
+	DO 56 I=1,IMAX
+56	PROD=PROD*(DELTAX-DELTA(I))
+57	SUM=SUM+A(J)*PROD
+60	YOUT=SUM
+61	RETURN
+	END