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(*),Y(*)
	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