diff options
author | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
---|---|---|
committer | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
commit | 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch) | |
tree | 4464880c571602d54f6ae114729bf62a89518057 /math/bevington/interp.f | |
download | iraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz |
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'math/bevington/interp.f')
-rw-r--r-- | math/bevington/interp.f | 85 |
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 |