C*********************************************************************** C * C NATURAL CUBIC SPLINE * C * C*********************************************************************** C C C C C TO CONSTRUCT THE CUBIC SPLINE INTERPOLANT S FOR C THE FUNCTION F, DEFINED AT THE NUMBERS C X(0) < X(1) < ... < X(N), SATISFYING C S''(X(0)) = S''(X(N)) = 0: C C INPUT: N; X(O),X(1), ...,X(N); EITHER GENERATE C A(I) = F(X(I)) FOR I = 0,1, ...,N OR INPUT C A(I) FOR I = 0,1, ...,N. C C OUTPUT: A(J),B(J),C(J),D(J) FOR J = 0,1, ...,N-1. C C NOTE: S(X) = A(J) + B(J)*(X-X(J)) + C(J)*(X-X(J))**2 + C D(J)*(X-X(J))**3 FOR X(J) < X < X(J+1) C DIMENSION X(10),A(10),B(10),C(10),D(10),H(10),XA(10),XL(10), 1XU(10),XZ(10) INTEGER N, I, J, M CHARACTER NAME*14,NAME1*14,AA*1 INTEGER INP,OUP,FLAG LOGICAL OK F(XX)=EXP(2.0*XX) OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL') OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL') WRITE(6,*) 'This is the natural cubic spline interpolation.' OK = .FALSE. 10 IF ( .NOT. OK) THEN WRITE(6,*) 'Choice of input method: ' WRITE(6,*) '1. Input entry by entry from keyboard ' WRITE(6,*) '2. Input data from a text file ' WRITE(6,*) '3. Generate data using a function F with nodes' WRITE(6,*) 'entered from keyboard ' WRITE(6,*) '4. Generate data using a function F with nodes' WRITE(6,*) 'from a text file ' WRITE(6,*) 'Choose 1, 2, 3, or 4 please ' WRITE(6,*) ' ' READ(5,*) FLAG IF( ( FLAG .GE. 1 ) .AND. ( FLAG .LE. 4 )) OK = .TRUE. GOTO 10 ENDIF IF (FLAG .EQ. 1) THEN OK = .FALSE. 20 IF (.NOT. OK ) THEN WRITE(6,*) 'Input number N ' WRITE(6,*) ' ' READ(5,*) N IF (N .GT. 0 ) THEN OK = .TRUE. N=N+1 DO 30 I = 1, N J=I-1 WRITE(6,*) 'Input X(',J,') and F(X(',J,')) ' WRITE(6,*) 'separated by space - include ' WRITE(6,*) 'decimal point ' WRITE(6,*) ' ' READ(5,*) X(I), A(I) 30 CONTINUE ELSE WRITE(6,*) 'Number must be a positive integer ' ENDIF GOTO 20 ENDIF ENDIF IF (FLAG .EQ. 2) THEN WRITE(6,*) 'Has a text file been created with data in two ' WRITE(6,*) 'columns? ' WRITE(6,*) 'Enter Y or N - letter within quotes ' WRITE(6,*) ' ' READ(5,*) AA IF (( AA .EQ. 'Y' ) .OR.( AA .EQ. 'y' )) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext contained in quotes' WRITE(6,*) 'as example: ''A:DATA.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME INP = 4 OPEN(UNIT=INP,FILE=NAME,ACCESS='SEQUENTIAL') OK = .FALSE. 40 IF (.NOT. OK) THEN WRITE(6,*) 'Input number N ' WRITE(6,*) ' ' READ(5,*) N IF ( N .GT. 0) THEN OK = .TRUE. N=N+1 DO 50 I = 1, N READ(4,*) X(I) , A(I) 50 CONTINUE CLOSE(UNIT=4) ELSE WRITE(6,*) 'Number must be a positive integer ' ENDIF GOTO 40 ENDIF ELSE WRITE(6,*) 'Please create the input file in two column ' WRITE(6,*) 'form with the ' WRITE(6,*) 'X values and F(X) values in the ' WRITE(6,*) 'corresponding columns ' WRITE(6,*) 'The program will end so the input file can ' WRITE(6,*) 'be created. ' OK = .FALSE. ENDIF ENDIF IF (FLAG .EQ. 3) THEN WRITE(6,*) 'Has the function F been created in the program?' WRITE(6,*) 'Enter Y or N - letter within quotes' WRITE(6,*) ' ' READ(5,*) AA IF (( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN OK = .FALSE. 60 IF (.NOT. OK) THEN WRITE(6,*) 'Input number N ' WRITE(6,*) ' ' READ(5,*) N IF ( N .GT. 0 ) THEN OK = .TRUE. N=N+1 DO 70 I = 1, N WRITE(6,*) 'Input X(',I,') - include ' WRITE(6,*) 'decimal point ' WRITE(6,*) ' ' READ(5,*) X(I) A(I) = F( X(I) ) 70 CONTINUE ELSE WRITE(6,*) 'Number must be a positive integer ' ENDIF GOTO 60 ENDIF ELSE WRITE(6,*) 'The program will end so that the function F ' WRITE(6,*) 'can be created ' OK = .FALSE. ENDIF ENDIF IF ( FLAG .EQ. 4) THEN WRITE(6,*) 'Has the text file with X-values been created and ' WRITE(6,*) 'has the function F been created in the program?' WRITE(6,*) 'Enter Y or N - letter within quotes' WRITE(6,*) ' ' READ(5,*) AA IF (( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y')) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext contained within quotes' WRITE(6,*) 'as example: ''A:DATA.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME INP = 4 OPEN(UNIT=4,FILE=NAME,ACCESS='SEQUENTIAL') OK = .FALSE. 80 IF (.NOT. OK ) THEN WRITE(6,*) 'Input number N ' WRITE(6,*) ' ' READ(5,*) N IF ( N .GT. 0 ) THEN OK = .TRUE. N=N+1 DO 90 I = 1, N READ(4,*) X(I) A(I) = F(X(I)) 90 CONTINUE CLOSE(UNIT=4) ELSE WRITE(6,*) 'Number must be a positive integer ' ENDIF GOTO 80 ENDIF ELSE WRITE(6,*) 'Please create the input file with one entry ' WRITE(6,*) 'per row and create the function F. ' WRITE(6,*) 'The program will end so the input file and ' WRITE(6,*) 'F can be created. ' OK = .FALSE. ENDIF ENDIF IF(.NOT.OK) GOTO 400 WRITE(6,*) 'Select output destinations: ' WRITE(6,*) '1. Screen ' WRITE(6,*) '2. Text file ' WRITE(6,*) 'Enter 1 or 2 ' WRITE(6,*) ' ' READ(5,*) FLAG IF ( FLAG .EQ. 2 ) THEN WRITE(6,*) 'Input the file name in the form - ' WRITE(6,*) 'drive:name.ext' WRITE(6,*) 'with the name contained within quotes' WRITE(6,*) 'as example: ''A:OUTPUT.DTA'' ' WRITE(6,*) ' ' READ(5,*) NAME1 OUP = 3 OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW') ELSE OUP = 6 ENDIF WRITE(OUP,*) 'NATURAL CUBIC SPLIINE INTERPOLATION ' WRITE(OUP,*) 'The numbers X(0), ..., X(N) are: ' WRITE(OUP,1) (X(I),I=1,N) 1 FORMAT(5(2X,E15.8)) C STEP 1 M=N-1 DO 11 I=1,M 11 H(I) = X(I+1)-X(I) C A(N)=F(X(N)), IF AVAILABLE C STEP 2 DO 21 I=2,M C USE XA INSTEAD OF ALPHA 21 XA(I) = 3*(A(I+1)*H(I-1)-A(I)*(X(I+1)-X(I-1))+A(I-1)*H(I))/ *(H(I)*H(I-1)) C STEP 3 C STEPS 3,4,5 AND PART OF 6 SOLVE THE TRIDIAGONAL SYSTEM USING C ALGORITHM 6.7. C C USE XL, XU, XZ IN PLACE OF L, MU, Z RESP. XL(1) = 1.0 XU(1) = 0.0 XZ(1) = 0.0 C STEP 4 DO 31 I=2,M XL(I) = 2*(X(I+1)-X(I-1))-H(I-1)*XU(I-1) XU(I) = H(I)/XL(I) 31 XZ(I) = (XA(I)-H(I-1)*XZ(I-1))/XL(I) C STEP 5 XL(N) = 1.0 XZ(N) = 0.0 C(N) = XZ(N) C STEP 6 DO 41 I=1,M J = N-I C(J) = XZ(J)-XU(J)*C(J+1) B(J) = (A(J+1)-A(J))/H(J)-H(J)*(C(J+1)+2*C(J))/3 41 D(J) = (C(J+1)-C(J))/(3*H(J)) C STEP 7 C OUTPUT WRITE(OUP,3) M WRITE(OUP,4)(A(I),B(I),C(I),D(I),I=1,M) 400 CLOSE(UNIT=5) CLOSE(UNIT=OUP) IF(OUP.NE.6) CLOSE(UNIT=6) STOP 3 FORMAT(1X,'SPLINE CONSISTS OF',3X,I3,3X,'PIECES WITH COEF. * A(I),B(I),C(I),D(I)') 4 FORMAT((1X,4(E15.8,3X))) END