PROGRAM LEASTSQ C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 5.1 (Least Squares Line). C Section 5.1, Least-Squares Line, Page 264 C PARAMETER(MaxN=200) INTEGER N REAL A,B,X,Y CHARACTER ANS*1 DIMENSION X(1:MaxN),Y(1:MaxN) 10 CALL GETDATA(X,Y,N) CALL REGRESS(X,Y,N,A,B) CALL RESULTS(X,Y,N,A,B) WRITE(9,*)' ' WRITE(9,*)'WANT TO FIND ANOTHER LEAST SQUARES LINE ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END SUBROUTINE REGRESS(X,Y,N,A,B) PARAMETER(MaxN=200) INTEGER K,N REAL A,B,SumX,SumXY,X,Xmean,Y,Ymean DIMENSION X(1:MaxN),Y(1:MaxN) Xmean=0 DO K=1,N Xmean=Xmean+X(K) ENDDO Xmean=Xmean/N Ymean=0 DO K=1,N Ymean=Ymean+Y(K) ENDDO Ymean=Ymean/N SumX=0 DO K=1,N SumX=SumX+(X(K)-Xmean)*(X(K)-Xmean) ENDDO SumXY=0 DO K=1,N SumXY=SumXY+(X(K)-Xmean)*(Y(K)-Ymean) ENDDO A=SumXY/SumX B=Ymean-A*Xmean RETURN END SUBROUTINE XREGRESS(X,Y,N,A,B) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=200) INTEGER K,N REAL A,B,SumX,SumXY,X,Xmean,Y,Ymean DIMENSION X(1:MaxN),Y(1:MaxN) Xmean=0 DO 10 K=1,N Xmean=Xmean+X(K) 10 CONTINUE Xmean=Xmean/N Ymean=0 DO 20 K=1,N Ymean=Ymean+Y(K) 20 CONTINUE Ymean=Ymean/N SumX=0 DO 30 K=1,N SumX=SumX+(X(K)-Xmean)*(X(K)-Xmean) 30 CONTINUE SumXY=0 DO 40 K=1,N SumXY=SumXY+(X(K)-Xmean)*(Y(K)-Ymean) 40 CONTINUE A=SumXY/SumX B=Ymean-A*Xmean RETURN END SUBROUTINE GETDATA(X,Y,N) PARAMETER(MaxN=200) INTEGER I,K,N REAL X,Y DIMENSION X(1:MaxN),Y(1:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE LEAST SQUARES LINE Y = AX + B IS CONSTRUCTED.' WRITE(9,*)' ' WRITE(9,*)'IT FITS THE POINTS (X ,Y ), (X ,Y ),..., (X ,Y )' WRITE(9,*)' 1 1 2 2 N N ' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER OF POINTS N = ' READ(9,*) N IF (N.LT.2) THEN N=2 ENDIF WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',N,' POINTS' WRITE(9,*)' ' DO 20 K=1,N WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) WRITE(9,*)'Y(',K,') = ' READ(9,*) Y(K) WRITE(9,*)' ' 20 CONTINUE RETURN END SUBROUTINE RESULTS(X,Y,N,A,B) PARAMETER(MaxN=200) INTEGER I,K,N REAL A,B,Err,X,Y,Z DIMENSION X(1:MaxN),Y(1:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'THE LEAST SQUARES LINE IS Y = A*X + B, WHERE:' WRITE(9,*)' ' WRITE(9,*)'A = ',A,' B = ',B WRITE(9,*)' ' WRITE(9,*)' K X Y A*X + B Error' WRITE(9,*)' K K K' WRITE(9,*)' ---------------------------------------------------- +----' DO 20 K=1,N Z=A*X(K)+B Err=Y(K)-Z WRITE(9,999) K,X(K),Y(K),Z,Err 20 CONTINUE 999 FORMAT(I5,2F12.5,2F14.7) RETURN END PROGRAM LPOLYNOM C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 5.2 (Least Squares Polynomial). C Section 5.2, Curve Fitting, Page 278 C PARAMETER(MaxN=12,MaxM=100) INTEGER M,N REAL A,B,C,P,X,Y DIMENSION A(1:MaxN,1:MaxN),B(1:MaxN),C(1:MaxN),X(1:MaxM),Y(1:MaxM) CHARACTER ANS*1 EXTERNAL P 10 CALL GETDATA(X,Y,M,N) CALL FILLMAT(X,Y,M,A,B,N) CALL SOLVELI(A,B,N,C) 20 CALL RESULTS(C,X,Y,M,N) CALL EVALUATE(C,N) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE P(T) AGAIN ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO FIND ANOTHER LEAST SQUARES POLYNOMIAL ? ' READ(9,'(A)') ANS IF(ANS .EQ. 'Y' .OR. ANS .EQ. 'y') GOTO 10 STOP END SUBROUTINE FILLMAT(X,Y,M,A,B,N) PARAMETER(MaxN=12,MaxM=100) INTEGER Col,J,K,M,N,R REAL A,B,Pow,Prod,X,XK,Y,YK DIMENSION A(1:MaxN,1:MaxN),B(1:MaxN),X(1:MaxM),Y(1:MaxM) DIMENSION Pow(0:MaxN) DO R=1,(N+1) B(R)=0 ENDDO DO K=1,M YK=Y(K) XK=X(K) Prod=1 DO R=1,(N+1) B(R)=B(R)+YK*Prod Prod=Prod*XK ENDDO ENDDO DO J=1,(2*N) Pow(J)=0 ENDDO Pow(0)=M DO K=1,M XK=X(K) Prod=X(K) DO J=1,(2*N) Pow(J)=Pow(J)+Prod Prod=Prod*XK ENDDO ENDDO DO R=1,(N+1) DO Col=1,(N+1) A(R,Col)=Pow(R+Col-2) ENDDO ENDDO RETURN END SUBROUTINE SOLVELI(A,B,N,C) PARAMETER(MaxN=12) INTEGER Col,J,K,N,P,Row,T REAL A,B,C,Sum DIMENSION A(1:MaxN,1:MaxN),B(1:MaxN),C(1:MaxN) DIMENSION Row(1:MaxN),Z(1:MaxN) DO J=1,(N+1) Row(J)=J ENDDO DO P=1,N DO K=(P+1),(N+1) IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T=Row(P) Row(P)=Row(K) Row(K)=T ENDIF ENDDO DO K=(P+1),(N+1) A(Row(K),P)=A(Row(K),P)/A(Row(P),P) DO Col=(P+1),(N+1) A(Row(K),Col)=A(Row(K),Col)-A(Row(K),P)*A(Row(P),Col) ENDDO ENDDO ENDDO Z(1)=B(Row(1)) DO K=2,(N+1) Sum=0 DO Col=1,(K-1) Sum=Sum+A(Row(K),Col)*Z(Col) ENDDO Z(K)=B(Row(K))-Sum ENDDO C(N+1)=Z(N+1)/A(Row(N+1),N+1) DO K=N,1,-1 Sum=0 DO Col=(K+1),(N+1) Sum=Sum+A(Row(K),Col)*C(Col) ENDDO C(K)=(Z(K)-Sum)/A(Row(K),K) ENDDO RETURN END REAL FUNCTION P(C,N,T) PARAMETER(MaxN=12) INTEGER K,N REAL C,Sum,T DIMENSION C(1:MaxN) Sum=C(N+1) DO K=N,1,-1 Sum=C(K)+Sum*T ENDDO P=Sum RETURN END SUBROUTINE XFILLMAT(X,Y,M,A,B,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=12,MaxM=100) INTEGER Col,J,K,M,N,R REAL A,B,Pow,Prod,X,XK,Y,YK DIMENSION A(1:MaxN,1:MaxN),B(1:MaxN),X(1:MaxM),Y(1:MaxM) DIMENSION Pow(0:MaxN) DO 10 R=1,(N+1) B(R)=0 10 CONTINUE DO 30 K=1,M YK=Y(K) XK=X(K) Prod=1 DO 20 R=1,(N+1) B(R)=B(R)+YK*Prod Prod=Prod*XK 20 CONTINUE 30 CONTINUE DO 40 J=1,(2*N) Pow(J)=0 40 CONTINUE Pow(0)=M DO 60 K=1,M XK=X(K) Prod=X(K) DO 50 J=1,(2*N) Pow(J)=Pow(J)+Prod Prod=Prod*XK 50 CONTINUE 60 CONTINUE DO 80 R=1,(N+1) DO 70 Col=1,(N+1) A(R,Col)=Pow(R+Col-2) 70 CONTINUE 80 CONTINUE RETURN END SUBROUTINE XSOLVELI(A,B,N,C) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=12) INTEGER Col,J,K,N,P,Row,T REAL A,B,C,Sum DIMENSION A(1:MaxN,1:MaxN),B(1:MaxN),C(1:MaxN) DIMENSION Row(1:MaxN),Z(1:MaxN) DO 10 J=1,(N+1) Row(J)=J 10 CONTINUE DO 50 P=1,N DO 20 K=(P+1),(N+1) IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN T=Row(P) Row(P)=Row(K) Row(K)=T ENDIF 20 CONTINUE DO 40 K=(P+1),(N+1) A(Row(K),P)=A(Row(K),P)/A(Row(P),P) DO 30 Col=(P+1),(N+1) A(Row(K),Col)=A(Row(K),Col)-A(Row(K),P)*A(Row(P),Col) 30 CONTINUE 40 CONTINUE 50 CONTINUE Z(1)=B(Row(1)) DO 70 K=2,(N+1) Sum=0 DO 60 Col=1,(K-1) Sum=Sum+A(Row(K),Col)*Z(Col) 60 CONTINUE Z(K)=B(Row(K))-Sum 70 CONTINUE C(N+1)=Z(N+1)/A(Row(N+1),N+1) DO 90 K=N,1,-1 Sum=0 DO 80 Col=(K+1),(N+1) Sum=Sum+A(Row(K),Col)*C(Col) 80 CONTINUE C(K)=(Z(K)-Sum)/A(Row(K),K) 90 CONTINUE RETURN END REAL FUNCTION XP(C,N,T) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=12) INTEGER K,N REAL C,Sum,T DIMENSION C(1:MaxN) Sum=C(N+1) DO 10 K=N,1,-1 Sum=C(K)+Sum*T 10 CONTINUE P=Sum RETURN END SUBROUTINE PRINTPOL(C,N) PARAMETER(MaxN=12) INTEGER I,K,N REAL C DIMENSION C(1:MaxN) IF (N.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'P(X) = C X + 2' WRITE(9,*)' 2 1' ELSEIF (N.EQ.2) THEN WRITE(9,*)' 2' WRITE(9,*)'P(X) = C X + C X + C' WRITE(9,*)' 3 2 1' ELSEIF (N.EQ.3) THEN WRITE(9,*)' 3 2' WRITE(9,*)'P(X) = C X + C X + C X + C' WRITE(9,*)' 4 3 2 1' ELSE WRITE(9,*)' ',N,' 2' WRITE(9,*)'P(X) = A C +...+ C X + C X + C' WRITE(9,*) N,' 3 2 1' ENDIF WRITE(9,*)' ' DO 10 K=(N+1),1,-2 IF (K.GE.2) THEN WRITE(9,*)'C(',K,') =',C(K),' C(',(K-1),') =',C(K-1) ELSE WRITE(9,*)'C(',K,') =',C(K) ENDIF 10 CONTINUE RETURN END SUBROUTINE GETDATA(X,Y,M,N) PARAMETER(MaxM=100) INTEGER I,K,M,N REAL X,Y DIMENSION X(1:MaxM),Y(1:MaxM) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE LEAST SQUARES POLYNOMIAL OF DEGREE N IS CONSTRUCTED +' WRITE(9,*)' ' WRITE(9,*)'WHICH FITS THE M POINTS (X ,Y ),(X ,Y ),...,(X ,Y ).' WRITE(9,*)' 1 1 2 2 M M' WRITE(9,*)' ' WRITE(9,*)'ENTER DEGREE N OF THE LEAST SQUARES POLYNOMIAL.' WRITE(9,*)' ' WRITE(9,*)'N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER OF POINTS M = ' READ(9,*) M IF (M.LT.(N+1)) THEN WRITE(9,*)' ' WRITE(9,*)'YOU MUST HAVE AT LEAST ',N+1,' POINTS.' M=N+1 ENDIF WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',M,' POINTS' WRITE(9,*)' ' WRITE(9,*)'(X ,Y ),(X ,Y ),...,(X ,Y )' WRITE(9,*)' 1 1 2 2 M M' WRITE(9,*)' ' DO 20 K=1,M WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) WRITE(9,*)'Y(',K,') = ' READ(9,*) Y(K) WRITE(9,*)' ' 20 CONTINUE RETURN END SUBROUTINE RESULTS(C,X,Y,M,N) PARAMETER(MaxN=12) INTEGER I,K,M,N REAL C,Err,P,X,Y,Z DIMENSION C(1:MaxN),X(1:MaxM),Y(1:MaxM) EXTERNAL P DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'THE LEAST SQUARES POLYNOMIAL IS:' WRITE(9,*)' ' CALL PRINTPOL(C,N) WRITE(9,*)' ' WRITE(9,*)' K X Y P(X ) Error' WRITE(9,*)' K K K' WRITE(9,*)' ----------------------------------------------------- +----' DO 20 K=1,M Z=P(C,N,X(K)) Err=Y(K)-Z WRITE(9,999) K,X(K),Y(K),Z,Err 20 CONTINUE 999 FORMAT(I5,2F12.5,2F14.7) RETURN END SUBROUTINE EVALUATE(C,N) PARAMETER(MaxN=12) INTEGER N REAL C,P,T DIMENSION C(1:MaxN) EXTERNAL P WRITE(9,*)' ' WRITE(9,*)'NOW EVALUATE P(T) ENTER A VALUE T = ' READ(9,*) T WRITE(9,*)'P(',T,' ) = ',P(C,N,T) RETURN END PROGRAM CURVEFIT C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 5.3 (Non-linear Curve Fitting). C Section 5.2, Curve Fitting, Page 280 C PARAMETER(MaxN=100) INTEGER Ctype,N REAL A,B,C,D,F,L,X,Y,X1,Y1 CHARACTER ANS*1 DIMENSION X(1:MaxN),X1(1:MaxN),Y(1:MaxN),Y1(1:MaxN) EXTERNAL F 10 CALL GETDATA(X,Y,X1,Y1,N) 20 CALL CURTYPE(Y1,N,Ctype,L) CALL CHANGEV(X1,Y1,X,Y,N,Ctype,L) CALL REGRESS(X,Y,N,A,B) CALL CONSTAN(A,B,C,D,Ctype) CALL RESULTS(X1,Y1,N,A,B,C,D,Ctype,L) WRITE(9,*)' ' WRITE(9,*)'WANT TO FIT ANOTHER CURVE TO THIS DATA ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO FIT A CURVE FOR SOME NEW DATA ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 STOP END REAL FUNCTION F(X,A,B,C,D,Ctype,L) INTEGER Ctype REAL A,B,C,D,L,X IF (Ctype.EQ.1) THEN F=A*X+B ELSEIF (Ctype.EQ.2) THEN F=A/X+B ELSEIF (Ctype.EQ.3) THEN F=D/(X+C) ELSEIF (Ctype.EQ.4) THEN F=1/(A*X+B) ELSEIF (Ctype.EQ.5) THEN F=X/(A+B*X) ELSEIF (Ctype.EQ.6) THEN F=A*ALOG(X)+B ELSEIF (Ctype.EQ.7) THEN F=C*EXP(A*X) ELSEIF (Ctype.EQ.8) THEN F=C*EXP(A*ALOG(X)) ELSEIF (Ctype.EQ.9) THEN F=1/((A*X+B)*(A*X+B)) ELSEIF (Ctype.EQ.10) THEN F=C*X*EXP(-D*X) ELSEIF (Ctype.EQ.11) THEN F=L/(1.0+C*EXP(A*X)) ENDIF RETURN END SUBROUTINE CHANGEV(X1,Y1,X,Y,N,Ctype,L) PARAMETER(MaxN=100) INTEGER Ctype,K,N REAL L,X,Y,X1,Y1 DIMENSION X(1:MaxN),X1(1:MaxN),Y(1:MaxN),Y1(1:MaxN) IF (Ctype.EQ.1) THEN DO K=1,N X(K)=X1(K) Y(K)=Y1(K) ENDDO ELSEIF (Ctype.EQ.2) THEN DO K=1,N X(K)=1.0/X1(K) Y(K)=Y1(K) ENDDO ELSEIF (Ctype.EQ.3) THEN DO K=1,N X(K)=X1(K)*Y1(K) Y(K)=Y1(K) ENDDO ELSEIF (Ctype.EQ.4) THEN DO K=1,N X(K)=X1(K) Y(K)=1.0/Y1(K) ENDDO ELSEIF (Ctype.EQ.5) THEN DO K=1,N X(K)=1.0/X1(K) Y(K)=1.0/Y1(K) ENDDO ELSEIF (Ctype.EQ.6) THEN DO K=1,N X(K)=ALOG(X1(K)) Y(K)=Y1(K) ENDDO ELSEIF (Ctype.EQ.7) THEN DO K=1,N X(K)=X1(K) Y(K)=ALOG(Y1(K)) ENDDO ELSEIF (Ctype.EQ.8) THEN DO K=1,N X(K)=ALOG(X1(K)) Y(K)=ALOG(Y1(K)) ENDDO ELSEIF (Ctype.EQ.9) THEN DO K=1,N X(K)=X1(K) Y(K)=1.0/SQRT(Y1(K)) ENDDO ELSEIF (Ctype.EQ.10) THEN DO K=1,N X(K)=X1(K) Y(K)=ALOG(Y1(K)/X1(K)) ENDDO ELSEIF (Ctype.EQ.11) THEN DO K=1,N X(K)=X1(K) Y(K)=ALOG(L/Y1(K)-1.0) ENDDO ENDIF RETURN END SUBROUTINE REGRESS(X,Y,N,A,B) PARAMETER(MaxN=100) INTEGER K,N REAL A,B,SumX,SumXY,X,Xmean,Y,Ymean DIMENSION X(1:MaxN),Y(1:MaxN) Xmean=0.0 DO K=1,N Xmean=Xmean+X(K) ENDDO Xmean=Xmean/N Ymean=0.0 DO K=1,N Ymean=Ymean+Y(K) ENDDO Ymean=Ymean/N SumX=0.0 DO K=1,N SumX=SumX+(X(K)-Xmean)*(X(K)-Xmean) ENDDO SumXY=0.0 DO K=1,N SumXY=SumXY+(X(K)-Xmean)*(Y(K)-Ymean) ENDDO A=SumXY/SumX B=Ymean-A*Xmean RETURN END SUBROUTINE CONSTAN(A,B,C,D,Ctype) INTEGER Ctype REAL A,B,C,D IF (Ctype.EQ.3) THEN C=-1.0/A D=-B/A ELSEIF (Ctype.EQ.7) THEN C=EXP(B) ELSEIF (Ctype.EQ.8) THEN C=EXP(B) ELSEIF (Ctype.EQ.10) THEN C=EXP(B) D=-A ELSEIF (Ctype.EQ.11) THEN C=EXP(B) ENDIF RETURN END REAL FUNCTION XF(X,A,B,C,D,Ctype,L) INTEGER Ctype REAL A,B,C,D,L,X IF (Ctype.EQ.1) THEN F=A*X+B ELSEIF (Ctype.EQ.2) THEN F=A/X+B ELSEIF (Ctype.EQ.3) THEN F=D/(X+C) ELSEIF (Ctype.EQ.4) THEN F=1/(A*X+B) ELSEIF (Ctype.EQ.5) THEN F=X/(A+B*X) ELSEIF (Ctype.EQ.6) THEN F=A*ALOG(X)+B ELSEIF (Ctype.EQ.7) THEN F=C*EXP(A*X) ELSEIF (Ctype.EQ.8) THEN F=C*EXP(A*ALOG(X)) ELSEIF (Ctype.EQ.9) THEN F=1/((A*X+B)*(A*X+B)) ELSEIF (Ctype.EQ.10) THEN F=C*X*EXP(-D*X) ELSEIF (Ctype.EQ.11) THEN F=L/(1.0+C*EXP(A*X)) ENDIF RETURN END SUBROUTINE XCHANGEV(X1,Y1,X,Y,N,Ctype,L) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=100) INTEGER Ctype,K,N REAL L,X,Y,X1,Y1 DIMENSION X(1:MaxN),X1(1:MaxN),Y(1:MaxN),Y1(1:MaxN) IF (Ctype.EQ.1) THEN DO 10 K=1,N X(K)=X1(K) Y(K)=Y1(K) 10 CONTINUE ELSEIF (Ctype.EQ.2) THEN DO 20 K=1,N X(K)=1.0/X1(K) Y(K)=Y1(K) 20 CONTINUE ELSEIF (Ctype.EQ.3) THEN DO 30 K=1,N X(K)=X1(K)*Y1(K) Y(K)=Y1(K) 30 CONTINUE ELSEIF (Ctype.EQ.4) THEN DO 40 K=1,N X(K)=X1(K) Y(K)=1.0/Y1(K) 40 CONTINUE ELSEIF (Ctype.EQ.5) THEN DO 50 K=1,N X(K)=1.0/X1(K) Y(K)=1.0/Y1(K) 50 CONTINUE ELSEIF (Ctype.EQ.6) THEN DO 60 K=1,N X(K)=ALOG(X1(K)) Y(K)=Y1(K) 60 CONTINUE ELSEIF (Ctype.EQ.7) THEN DO 70 K=1,N X(K)=X1(K) Y(K)=ALOG(Y1(K)) 70 CONTINUE ELSEIF (Ctype.EQ.8) THEN DO 80 K=1,N X(K)=ALOG(X1(K)) Y(K)=ALOG(Y1(K)) 80 CONTINUE ELSEIF (Ctype.EQ.9) THEN DO 90 K=1,N X(K)=X1(K) Y(K)=1.0/SQRT(Y1(K)) 90 CONTINUE ELSEIF (Ctype.EQ.10) THEN DO 100 K=1,N X(K)=X1(K) Y(K)=ALOG(Y1(K)/X1(K)) 100 CONTINUE ELSEIF (Ctype.EQ.11) THEN DO 110 K=1,N X(K)=X1(K) Y(K)=ALOG(L/Y1(K)-1.0) 110 CONTINUE ENDIF RETURN END SUBROUTINE XREGRESS(X,Y,N,A,B) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=100) INTEGER K,N REAL A,B,SumX,SumXY,X,Xmean,Y,Ymean DIMENSION X(1:MaxN),Y(1:MaxN) Xmean=0.0 DO 10 K=1,N Xmean=Xmean+X(K) 10 CONTINUE Xmean=Xmean/N Ymean=0.0 DO 20 K=1,N Ymean=Ymean+Y(K) 20 CONTINUE Ymean=Ymean/N SumX=0.0 DO 30 K=1,N SumX=SumX+(X(K)-Xmean)*(X(K)-Xmean) 30 CONTINUE SumXY=0.0 DO 40 K=1,N SumXY=SumXY+(X(K)-Xmean)*(Y(K)-Ymean) 40 CONTINUE A=SumXY/SumX B=Ymean-A*Xmean RETURN END SUBROUTINE XCONSTAN(A,B,C,D,Ctype) INTEGER Ctype REAL A,B,C,D IF (Ctype.EQ.3) THEN C=-1.0/A D=-B/A ELSEIF (Ctype.EQ.7) THEN C=EXP(B) ELSEIF (Ctype.EQ.8) THEN C=EXP(B) ELSEIF (Ctype.EQ.10) THEN C=EXP(B) D=-A ELSEIF (Ctype.EQ.11) THEN C=EXP(B) ENDIF RETURN END SUBROUTINE GETDATA(X,Y,X1,Y1,N) PARAMETER(MaxN=100) INTEGER I,K,N REAL X,Y,X1,Y1 DIMENSION X(1:MaxN),X1(1:MaxN),Y(1:MaxN),Y1(1:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' THE METHOD OF "DATA LINEARIZATION" IS USED' WRITE(9,*)' ' WRITE(9,*)'TO FIT ONE OF THE CURVES Y=A*X+B, Y=A/X+B' WRITE(9,*)' ' WRITE(9,*)'Y=D/(X+C), Y=1/(A*X+B), Y=X/(A+B*X), Y=A*Ln(X)+B' WRITE(9,*)' ' WRITE(9,*)'Y=C*Exp(A*X), Y=C*X^A, Y=(A*X+B)^-2' WRITE(9,*)' ' WRITE(9,*)'Y=C*X*Exp(-D*X), Y=L/(1+C*Exp(A*X))' WRITE(9,*)' ' WRITE(9,*)'TO THE N DATA POINTS (X ,Y ), (X ,Y ),..., (X ,Y ).' WRITE(9,*)' 1 1 2 2 N N ' WRITE(9,*)' ' WRITE(9,*)'ENTER NUMBER OF POINTS N = ' READ(9,*) N IF (N.LT.2) THEN N=2 ENDIF WRITE(9,*)' ' DO 20 K=1,N WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) X1(K)=X(K) WRITE(9,*)'Y(',K,') = ' READ(9,*) Y(K) Y1(K)=Y(K) WRITE(9,*)' ' 20 CONTINUE RETURN END SUBROUTINE CURTYPE(Y,N,Ctype,L) PARAMETER(MaxN=100) INTEGER Ctype,K,N REAL L,LM,Y DIMENSION Y(1:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' WHICH CURVE DO YOU WANT TO FIT ?' WRITE(9,*)' ' WRITE(9,*)' 1: Y = A*X + B' WRITE(9,*)' ' WRITE(9,*)' 2: Y = A/X + B' WRITE(9,*)' ' WRITE(9,*)' 3: Y = D/(X + C)' WRITE(9,*)' ' WRITE(9,*)' 4: Y = 1/(A*X + B)' WRITE(9,*)' ' WRITE(9,*)' 5: Y = X/(A + B*X)' WRITE(9,*)' ' WRITE(9,*)' 6: Y = A*Ln(X) + B' WRITE(9,*)' ' WRITE(9,*)' 7: Y = C*Exp(A*X)' WRITE(9,*)' ' WRITE(9,*)' 8: Y = C*X^A' WRITE(9,*)' ' WRITE(9,*)' 9: Y = (A*X + B)^-2' WRITE(9,*)' ' WRITE(9,*)'10: Y = C*X*EXP(-D*X)' WRITE(9,*)' ' WRITE(9,*)'11: Y = L/(1+C*EXP(A*X))' WRITE(9,*)' ' WRITE(9,*)' CHOOSE TYPE OF CURVE <1-11> ' READ(9,*) Ctype IF (Ctype.EQ.11) THEN DO 20 I=1,20 WRITE(9,*)' ' 20 CONTINUE WRITE(9,*)' ' WRITE(9,*)' YOU CHOSE THE CURVE Y = L/(1+C*EXP(A*X))' WRITE(9,*)' ' WRITE(9,*)'THE LIMITING VALUE L = LIM Y(X) MUST BE GIVEN.' WRITE(9,*)' X ->INF' WRITE(9,*)' ' WRITE(9,*)'ENTER THE VALUE L = ' READ(9,*) L LM=ABS(Y(1)) DO 30 K=2,N IF (LM.LT.ABS(Y(K))) THEN LM=ABS(Y(K)) ENDIF 30 CONTINUE IF (L.LE.LM) THEN L=LM*1.0001 ENDIF ENDIF RETURN END SUBROUTINE RESULTS(X,Y,N,A,B,C,D,Ctype,L) PARAMETER(MaxN=100) INTEGER Ctype,I,K,N REAL A,B,C,D,Err,F,L,X,Y,Z DIMENSION X(1:MaxN),Y(1:MaxN) EXTERNAL F DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE IF (Ctype.EQ.1) THEN WRITE(9,*)' F(X) = A*X +B' WRITE(9,*)' ' WRITE(9,*)' A = ',A WRITE(9,*)' ' WRITE(9,*)' B = ',B ELSEIF (Ctype.EQ.2) THEN WRITE(9,*)' F(X) = A/X + B' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' B =',B ELSEIF (Ctype.EQ.3) THEN WRITE(9,*)' F(X) = D/(X + C)' WRITE(9,*)' ' WRITE(9,*)' C =',C WRITE(9,*)' ' WRITE(9,*)' D =',D ELSEIF (Ctype.EQ.4) THEN WRITE(9,*)' F(X) = 1/(A*X + B)' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' B =',B ELSEIF (Ctype.EQ.5) THEN WRITE(9,*)' F(X) = X/(A + B*X)' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' B =',B ELSEIF (Ctype.EQ.6) THEN WRITE(9,*)' F(X) = A*Ln(X) + B' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' B =',B ELSEIF (Ctype.EQ.7) THEN WRITE(9,*)' F(X) = C*Exp(A*X)' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' C =',C ELSEIF (Ctype.EQ.8) THEN WRITE(9,*)' F(X) = C*X^A' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' C =',C ELSEIF (Ctype.EQ.9) THEN WRITE(9,*)' F(X) = (A*X + B)^-2' WRITE(9,*)' ' WRITE(9,*)' A =',A WRITE(9,*)' ' WRITE(9,*)' B =',B ELSEIF (Ctype.EQ.10) THEN WRITE(9,*)' F(X) = C*X*EXP(-D*X)' WRITE(9,*)' ' WRITE(9,*)' C =',C WRITE(9,*)' ' WRITE(9,*)' D =',D ELSEIF (Ctype.EQ.11) THEN WRITE(9,*)' F(X) = L/(1+C*EXP(A*X))' WRITE(9,*)' ' WRITE(9,*)' L =',L WRITE(9,*)' ' WRITE(9,*)' C =',C WRITE(9,*)' ' WRITE(9,*)' A =',A ENDIF WRITE(9,*)' ' WRITE(9,*)' K X Y F(X ) Error' WRITE(9,*)' K K K' WRITE(9,*)' ----------------------------------------------------- +----' DO 20 K=1,N Z=F(X(K),A,B,C,D,Ctype,L) Err=Y(K)-Z WRITE(9,999) K,X(K),Y(K),Z,Err 20 CONTINUE 999 FORMAT(I5,2F12.5,2F14.7) RETURN END PROGRAM CUBICSPL C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 5.4 (Cubic Splines). C Section 5.3, Interpolation by Spline Functions, Page 297 C PARAMETER(MaxN=50) INTEGER N,Stype REAL A,B,C,D,DX0,DXN,DDX0,DDXN,H,M,S,T,V,X,Y DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN),D(0:MaxN),H(0:MaxN) DIMENSION M(0:MaxN),V(0:MaxN),X(0:MaxN),Y(0:MaxN),S(0:MaxN,0:3) CHARACTER ANS*1 EXTERNAL CS 10 CALL GETDATA(X,Y,N) 20 CALL DIFFERS(X,Y,N,A,B,C,D,H,V) CALL SPLTYPE(Stype) CALL MODIFYS(A,B,C,D,H,M,V,X,N,Stype,DX0,DXN,DDX0,DDXN) CALL TRIDIAG(A,B,C,V,M,N) CALL ENDCOEF(D,H,M,N,Stype,DX0,DXN,DDX0,DDXN) CALL COMCOEF(D,H,M,S,Y,N) CALL RESULTS(S,X,Y,N,Stype,DX0,DXN,DDX0,DDXN) 30 CALL EVALUAT(CS,S,X,Y,N,Stype,DX0,DXN,DDX0,DDXN) WRITE(9,*)' ' WRITE(9,*)'WANT TO EVALUATE THIS SPLINE AT ANOTHER POINT ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 30 WRITE(9,*)' ' WRITE(9,*)'WANT TO FIT ANOTHER CUBIC SPLINE TO THIS DATA SET ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 20 WRITE(9,*)' ' WRITE(9,*)'WANT TO FIND A CUBIC SPLINE FOR SOME NEW DATA SET ? ' READ(9,'(A)') ANS IF (ANS.EQ.'Y' .OR. ANS.EQ.'y') GOTO 10 STOP END SUBROUTINE DIFFERS(X,Y,N,A,B,C,D,H,V) PARAMETER(MaxN=50) INTEGER K,N REAL A,B,C,D,H,V,X,Y DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN),D(0:MaxN) DIMENSION H(0:MaxN),V(0:MaxN),X(0:MaxN),Y(0:MaxN) H(0)= X(1)-X(0) D(0)=(Y(1)-Y(0))/H(0) DO K=1,(N-1) H(K)= X(K+1)-X(K) D(K)=(Y(K+1)-Y(K))/H(K) A(K)=H(K) B(K)=2*(H(K-1)+H(K)) C(K)=H(K) V(K)=6*(D(K)-D(K-1)) ENDDO RETURN END SUBROUTINE MODIFYS(A,B,C,D,H,M,V,X,N,Stype,DX0,DXN,DDX0,DDXN) PARAMETER(MaxN=50) INTEGER N,Stype REAL A,B,C,D,DX0,DXN,DDX0,DDXN,H,M,V DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN),D(0:MaxN) DIMENSION H(0:MaxN),M(0:MaxN),V(0:MaxN),X(0:MaxN) IF (Stype.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'ENTER THE FIRST DERIVATIVES.' WRITE(9,*)' ' WRITE(9,*)'S`(',X(0),' ) = ' READ(9,*) DX0 WRITE(9,*)' ' WRITE(9,*)'S`(',X(N),' ) = ' READ(9,*) DXN WRITE(9,*)' ' B(1)=B(1)-H(0)/2 V(1)=V(1)-3*(D(0)-DX0) B(N-1)=B(N-1)-H(N-1)/2 V(N-1)=V(N-1)-3*(DXN-D(N-1)) ELSEIF (Stype.EQ.2) THEN M(0)=0 M(N)=0 ELSEIF (Stype.EQ.3) THEN B(1)=B(1)+H(0)+H(0)*H(0)/H(1) C(1)=C(1)-H(0)*H(0)/H(1) B(N-1)=B(N-1)+H(N-1)+H(N-1)*H(N-1)/H(N-2) A(N-2)=A(N-2)-H(N-1)*H(N-1)/H(N-2) ELSEIF (Stype.EQ.4) THEN B(1)=B(1)+H(0) B(N-1)=B(N-1)+H(N-1) ELSEIF (Stype.EQ.5) THEN WRITE(9,*)' ' WRITE(9,*)'ENTER THE SECOND DERIVATIVES.' WRITE(9,*)' ' WRITE(9,*)'S``(',X(0),' ) = ' READ(9,*) DDX0 WRITE(9,*)' ' WRITE(9,*)'S``(',X(N),' ) = ' READ(9,*) DDXN WRITE(9,*)' ' V(1)=V(1)-H(0)*DDX0 V(N-1)=V(N-1)-H(N-1)*DDXN ENDIF RETURN END SUBROUTINE TRIDIAG(A,B,C,V,M,N) PARAMETER(MaxN=50) INTEGER K,N REAL A,B,C,M,T,V DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN) DIMENSION M(0:MaxN),V(0:MaxN) DO K=2,(N-1) T=A(K-1)/B(K-1) B(K)=B(K)-T*C(K-1) V(K)=V(K)-T*V(K-1) ENDDO M(N-1)=V(N-1)/B(N-1) DO K=(N-2),1,-1 M(K)=(V(K)-C(K)*M(K+1))/B(K) ENDDO RETURN END SUBROUTINE ENDCOEF(D,H,M,N,Stype,DX0,DXN,DDX0,DDXN) PARAMETER(MaxN=50) INTEGER N,Stype REAL D,DX0,DXN,DDX0,DDXN,H,M DIMENSION D(0:MaxN),H(0:MaxN),M(0:MaxN) IF (Stype.EQ.1) THEN M(0)=3*(D(0)-DX0)/H(0)-M(1)/2 M(N)=3*(DXN-D(N-1))/H(N-1)-M(N-1)/2 ELSEIF (Stype.EQ.2) THEN M(0)=0 M(N)=0 ELSEIF (Stype.EQ.3) THEN M(0)=M(1)-H(0)*(M(2)-M(1))/H(1) M(N)=M(N-1)+H(N-1)*(M(N-1)-M(N-2))/H(N-2) ELSEIF (Stype.EQ.4) THEN M(0)=M(1) M(N)=M(N-1) ELSEIF (Stype.EQ.5) THEN M(0)=DDX0 M(N)=DDXN ENDIF RETURN END SUBROUTINE COMCOEF(D,H,M,S,Y,N) PARAMETER(MaxN=50) INTEGER I,K,N REAL D,H,M,S DIMENSION D(0:MaxN),H(0:MaxN),M(0:MaxN),Y(0:MaxN),S(0:MaxN,0:3) DO K=0,(N-1) S(K,0)=Y(K) S(K,1)=D(K)-H(K)*(2*M(K)+M(K+1))/6 S(K,2)=M(K)/2 S(K,3)=(M(K+1)-M(K))/(6*H(K)) ENDDO RETURN END REAL FUNCTION CS(S,X,N,T) PARAMETER(MaxN=50) INTEGER J,K,N REAL S,T,W,X DIMENSION X(0:MaxN),S(0:MaxN,0:3) J=0 K=N-1 WHILE (J.LE.N) J=J+1 IF (T.LE.X(J)) THEN K=J-1 J=N+1 ENDIF REPEAT IF (T.LE.X(0)) K=0 W=T-X(K) CS=((S(K,3)*W + S(K,2))*W + S(K,1))*W + S(K,0) RETURN END SUBROUTINE XDIFFERS(X,Y,N,A,B,C,D,H,V) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER K,N REAL A,B,C,D,H,V,X,Y DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN),D(0:MaxN) DIMENSION H(0:MaxN),V(0:MaxN),X(0:MaxN),Y(0:MaxN) H(0)= X(1)-X(0) D(0)=(Y(1)-Y(0))/H(0) DO 10 K=1,(N-1) H(K)= X(K+1)-X(K) D(K)=(Y(K+1)-Y(K))/H(K) A(K)=H(K) B(K)=2*(H(K-1)+H(K)) C(K)=H(K) V(K)=6*(D(K)-D(K-1)) 10 CONTINUE RETURN END SUBROUTINE XMODIFYS(A,B,C,D,H,M,V,X,N,Stype,DX0,DXN,DDX0,DDXN) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER N,Stype REAL A,B,C,D,DX0,DXN,DDX0,DDXN,H,M,V DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN),D(0:MaxN) DIMENSION H(0:MaxN),M(0:MaxN),V(0:MaxN),X(0:MaxN) IF (Stype.EQ.1) THEN WRITE(9,*)' ' WRITE(9,*)'ENTER THE FIRST DERIVATIVES.' WRITE(9,*)' ' WRITE(9,*)'S`(',X(0),' ) = ' READ(9,*) DX0 WRITE(9,*)' ' WRITE(9,*)'S`(',X(N),' ) = ' READ(9,*) DXN WRITE(9,*)' ' B(1)=B(1)-H(0)/2 V(1)=V(1)-3*(D(0)-DX0) B(N-1)=B(N-1)-H(N-1)/2 V(N-1)=V(N-1)-3*(DXN-D(N-1)) ELSEIF (Stype.EQ.2) THEN M(0)=0 M(N)=0 ELSEIF (Stype.EQ.3) THEN B(1)=B(1)+H(0)+H(0)*H(0)/H(1) C(1)=C(1)-H(0)*H(0)/H(1) B(N-1)=B(N-1)+H(N-1)+H(N-1)*H(N-1)/H(N-2) A(N-2)=A(N-2)-H(N-1)*H(N-1)/H(N-2) ELSEIF (Stype.EQ.4) THEN B(1)=B(1)+H(0) B(N-1)=B(N-1)+H(N-1) ELSEIF (Stype.EQ.5) THEN WRITE(9,*)' ' WRITE(9,*)'ENTER THE SECOND DERIVATIVES.' WRITE(9,*)' ' WRITE(9,*)'S``(',X(0),' ) = ' READ(9,*) DDX0 WRITE(9,*)' ' WRITE(9,*)'S``(',X(N),' ) = ' READ(9,*) DDXN WRITE(9,*)' ' V(1)=V(1)-H(0)*DDX0 V(N-1)=V(N-1)-H(N-1)*DDXN ENDIF RETURN END SUBROUTINE XTRIDIAG(A,B,C,V,M,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER K,N REAL A,B,C,M,T,V DIMENSION A(0:MaxN),B(0:MaxN),C(0:MaxN) DIMENSION M(0:MaxN),V(0:MaxN) DO 10 K=2,(N-1) T=A(K-1)/B(K-1) B(K)=B(K)-T*C(K-1) V(K)=V(K)-T*V(K-1) 10 CONTINUE M(N-1)=V(N-1)/B(N-1) DO 20 K=(N-2),1,-1 M(K)=(V(K)-C(K)*M(K+1))/B(K) 20 CONTINUE RETURN END SUBROUTINE XENDCOEF(D,H,M,N,Stype,DX0,DXN,DDX0,DDXN) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER N,Stype REAL D,DX0,DXN,DDX0,DDXN,H,M DIMENSION D(0:MaxN),H(0:MaxN),M(0:MaxN) IF (Stype.EQ.1) THEN M(0)=3*(D(0)-DX0)/H(0)-M(1)/2 M(N)=3*(DXN-D(N-1))/H(N-1)-M(N-1)/2 ELSEIF (Stype.EQ.2) THEN M(0)=0 M(N)=0 ELSEIF (Stype.EQ.3) THEN M(0)=M(1)-H(0)*(M(2)-M(1))/H(1) M(N)=M(N-1)+H(N-1)*(M(N-1)-M(N-2))/H(N-2) ELSEIF (Stype.EQ.4) THEN M(0)=M(1) M(N)=M(N-1) ELSEIF (Stype.EQ.5) THEN M(0)=DDX0 M(N)=DDXN ENDIF RETURN END SUBROUTINE XCOMCOEF(D,H,M,S,Y,N) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER I,K,N REAL D,H,M,S DIMENSION D(0:MaxN),H(0:MaxN),M(0:MaxN),Y(0:MaxN),S(0:MaxN,0:3) DO 10 K=0,(N-1) S(K,0)=Y(K) S(K,1)=D(K)-H(K)*(2*M(K)+M(K+1))/6 S(K,2)=M(K)/2 S(K,3)=(M(K+1)-M(K))/(6*H(K)) 10 CONTINUE RETURN END REAL FUNCTION XCS(S,X,N,T) C This subroutine uses labeled DO loop(s). PARAMETER(MaxN=50) INTEGER J,K,N REAL S,T,W,X DIMENSION X(0:MaxN),S(0:MaxN,0:3) K=N-1 DO 10 J=1,N IF (T.LE.X(J)) THEN K=J-1 GOTO 20 ENDIF 10 CONTINUE 20 CONTINUE IF (T.LE.X(0)) THEN K=0 ENDIF W=T-X(K) CS=((S(K,3)*W + S(K,2))*W + S(K,1))*W + S(K,0) RETURN END SUBROUTINE GETDATA(X,Y,N) PARAMETER(MaxN=50) INTEGER I,K,N REAL X,Y DIMENSION X(0:MaxN),Y(0:MaxN) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' ' WRITE(9,*)'A CUBIC SPLINE IS CONSTRUCTED WHICH PASSES THROUGH' WRITE(9,*)' ' WRITE(9,*)'THE N+1 POINTS (X ,Y ), (X ,Y ),..., (X ,Y )' WRITE(9,*)' 0 0 1 1 N N ' WRITE(9,*)' ' WRITE(9,*)'THE SOLUTION IS A SET OF N PIECEWISE CUBIC FUNCTIONS' WRITE(9,*)' ' WRITE(9,*)' 3 2' WRITE(9,*)'S (X) = S (X-X ) + S (X-X ) + S (X-X ) + S' WRITE(9,*)' K K,3 K K,2 K K,1 K K,0' WRITE(9,*)' ' WRITE(9,*)' WHERE X IS IN [X ,X ]' WRITE(9,*)' K K+1 ' WRITE(9,*)' ' WRITE(9,*)' AND K = 0,1,...,N-1.' WRITE(9,*)' ' WRITE(9,*)'ENTER THE NUMBER N = ' READ(9,*) N WRITE(9,*)' ' WRITE(9,*)'NOW ENTER THE ',(N+1),' POINTS:' WRITE(9,*)' ' DO 20 K=0,N WRITE(9,*)'X(',K,') = ' READ(9,*) X(K) WRITE(9,*)'Y(',K,') = ' READ(9,*) Y(K) WRITE(9,*)' ' 20 CONTINUE RETURN END SUBROUTINE SPLTYPE(Stype) INTEGER I,Stype DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)'WHICH SPLINE TYPE DO YOU WANT? <1-5> ' WRITE(9,*)' ' WRITE(9,*)'1 - CLAMPED SPLINE' WRITE(9,*)' ' WRITE(9,*)'2 - NATURAL SPLINE' WRITE(9,*)' ' WRITE(9,*)'3 - EXTRAPOLATE S``(X) NEAR THE ENDPOINTS' WRITE(9,*)' ' WRITE(9,*)'4 - PARABOLIC NEAR THE ENDPOINTS' WRITE(9,*)' ' WRITE(9,*)'5 - SPECIFY S``(X(0)) AND S``(X(N))' WRITE(9,*)' ' WRITE(9,*)' CHOOSE TYPE <1-5> ' READ(9,*) Stype RETURN END SUBROUTINE RESULTS(S,X,Y,N,Stype,DX0,DXN,DDX0,DDXN) PARAMETER(MaxN=50) INTEGER I,K,N,Stype REAL DX0,DXN,DDX0,DDXN,S,X,Y DIMENSION X(0:MaxN),Y(0:MaxN),S(0:MaxN,0:3) DO 10 I=1,18 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' THE DATA POINTS ARE:' WRITE(9,*)' ' IF (Stype.EQ.1) THEN WRITE(9,*)' K X(K) Y(K) S`(X(K))' ELSEIF (Stype.EQ.5) THEN WRITE(9,*)' K X(K) Y(K) S``(X(K))' ELSE WRITE(9,*)' K X(K) Y(K)' ENDIF DO 20 K=0,N IF (Stype.EQ.1 .AND. K.EQ.0) THEN WRITE(9,999) K,X(K),Y(K),DX0 ELSEIF (Stype.EQ.1 .AND. K.EQ.N) THEN WRITE(9,999) K,X(K),Y(K),DXN ELSEIF (Stype.EQ.5 .AND. K.EQ.0) THEN WRITE(9,999) K,X(K),Y(K),DDX0 ELSEIF (Stype.EQ.5 .AND. K.EQ.N) THEN WRITE(9,999) K,X(K),Y(K),DDXN ELSE WRITE(9,999) K,X(K),Y(K) ENDIF 20 CONTINUE WRITE(9,*)' ' IF (Stype.EQ.1) THEN WRITE(9,*)'YOU CHOSE THE CLAMPED SPLINE WITH S`(X(0)) AND S`(X(N +)) SPECIFIED.' ELSEIF (Stype.EQ.2) THEN WRITE(9,*)'YOU CHOSE THE NATURAL SPLINE.' ELSEIF (Stype.EQ.3) THEN WRITE(9,*)'YOU CHOSE THE SPLINE WITH S``(X) EXTRAPOLATED NEAR TH +E ENDPOINTS.' ELSEIF (Stype.EQ.4) THEN WRITE(9,*)'YOU CHOSE THE SPLINE WHICH IS PARABOLIC NEAR THE ENDP +OINTS.' ELSEIF (Stype.EQ.5) THEN WRITE(9,*)'YOU CHOSE THE SPLINE WITH S``(X(0)) AND S``(X(N)) SPE +CIFIED.' ENDIF WRITE(9,*)' ' WRITE(9,*)' 3 2' WRITE(9,*)'S (X) = S(K,3)*(X-X(K)) + S(K,2)*(X-X(K)) + S(K,1)*(X +-X(K)) + S(K,0)' WRITE(9,*)' K' WRITE(9,*)' ' WRITE(9,*)' K X(K) S(K,3) S(K,2) +S(K,1) S(K,0)' DO 30 K=0,(N-1) WRITE(9,999) K,X(K),S(K,3),S(K,2),S(K,1),S(K,0) 30 CONTINUE WRITE(9,*)' ' 999 FORMAT(1X,I2,5F15.7) RETURN END SUBROUTINE EVALUAT(CS,S,X,Y,N,Stype,DX0,DXN,DDX0,DDXN) PARAMETER(MaxN=50) INTEGER N,Stype REAL DX0,DXN,DDX0,DDXN,S,T,X,Y,Z DIMENSION X(0:MaxN),Y(0:MaxN),S(0:MaxN,0:3) EXTERNAL CS WRITE(9,*)' ' WRITE(9,*)'ENTER A VALUE T = ' READ(9,*) T Z=CS(S,X,N,T) WRITE(9,*)'THE SPLINE`S VALUE IS S(',T,' ) =',Z RETURN END PROGRAM TRIGPOLY C NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994 C To accompany the text: C NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 C Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A. C This free software is complements of the author. C C Algorithm 5.5 (Trigonometric Polynomials). C Section 5.4, Fourier Series and Trigonometric Polynomials, Page 311 C PARAMETER(MaxN=361,MaxM=181,Pi=3.14159265358979323846,FunMax=7) INTEGER Ftype,FunType,Itype,M,M0,N,N0 REAL A,B,L,X,Y DIMENSION A(0:361),B(0:361),X(0:361),Y(0:361) EXTERNAL F,FX,TP CALL GETDATA(F,X,Y,M,M0,N,N0,A0,B0,L,Ftype,Itype,FunType,FunMax) CALL COEFFIC(X,Y,A,B,M,M0,N,N0,L) CALL RESULTS(TP,A,B,X,Y,XK,M,M0,N,N0,L,Itype,Ftype,FunType) STOP END FUNCTION FX(X,FunType) PARAMETER(Pi=3.14159265358979323846) INTEGER FunType REAL X IF (FunType.EQ.1) THEN IF (X.LT.-Pi) FX=X+2*Pi IF ((-Pi.LE.X).AND.(X.LT.0)) FX=0 IF ((0.LE.X).AND.(X.LT.Pi)) FX=X IF (Pi.LE.X) FX=0 ELSEIF (FunType.EQ.2) THEN IF (X.LT.-Pi) FX=1 IF ((-Pi.LE.X).AND.(X.LT.0)) FX=-1 IF ((0.LE.X).AND.(X.LT.Pi)) FX=1 IF (Pi.LE.X) FX=-1 ELSEIF (FunType.EQ.3) THEN IF (X.LT.-Pi) FX=-3*Pi/2 - X IF ((-Pi.LE.X).AND.(X.LT.0)) FX=Pi/2 + X IF ((0.LE.X).AND.(X.LT.Pi)) FX=Pi/2 - X IF (Pi.LE.X) FX=-3*Pi/2 + X ELSEIF (FunType.EQ.4) THEN IF (X.LT.-Pi) FX=-1 IF ((-Pi.LE.X).AND.(X.LT.-Pi/2)) FX=-1 IF ((-Pi/2.LE.X).AND.(X.LT.Pi/2)) FX=1 IF ((Pi/2.LE.X).AND.(X.LT.Pi)) FX=-1 IF (Pi.LE.X) FX=-1 ELSEIF (FunType.EQ.5) THEN IF (X.LT.-Pi) FX=Pi + X IF ((-Pi.LE.X).AND.(X.LT.-Pi/2)) FX=-Pi - X IF ((-Pi/2.LE.X).AND.(X.LT.Pi/2)) FX=X IF ((Pi/2.LE.X).AND.(X.LT.Pi)) FX=Pi - X IF (Pi.LE.X) FX=-Pi + X ELSEIF (FunType.EQ.6) THEN IF (X.LT.-Pi) FX=(X+2*Pi)*(X+2*Pi)/4 IF ((-Pi.LE.X).AND.(X.LT.Pi)) FX=X*X/4 IF (Pi.LE.X) FX=(X-2*Pi)*(X-2*Pi)/4 ELSEIF (FunType.EQ.7) THEN IF (X.LT.-Pi) FX=X/2+Pi IF ((-Pi.LE.X).AND.(X.LT.Pi)) FX=X/2 IF (Pi.LE.X) FX=X/2-Pi ENDIF RETURN END SUBROUTINE PRINTFUN(FunType) INTEGER FunType IF (FunType.EQ.1) THEN WRITE(9,*)'<1> F(X) = X + 2*Pi for X < -Pi ' WRITE(9,*)' F(X) = 0 for -Pi <= X < 0 ' WRITE(9,*)' F(X) = X for 0 <= X < Pi' WRITE(9,*)' F(X) = 0 for Pi <= X ' ELSEIF (FunType.EQ.2) THEN WRITE(9,*)'<2> F(X) = 1 for X <= -Pi ' WRITE(9,*)' F(X) = -1 for -Pi <= X < 0 ' WRITE(9,*)' F(X) = 1 for 0 <= X < Pi ' WRITE(9,*)' F(X) = -1 for Pi <= X ' ELSEIF (FunType.EQ.3) THEN WRITE(9,*)'<3> F(X) = -3*Pi/2 - X for X < -Pi' WRITE(9,*)' F(X) = Pi/2 + X for -Pi <= X < 0' WRITE(9,*)' F(X) = Pi/2 - X for 0 <= X < Pi' WRITE(9,*)' F(X) = -3Pi/2 + X for Pi <= X ' ELSEIF (FunType.EQ.4) THEN WRITE(9,*)'<4> F(X) = -1 for X < -Pi ' WRITE(9,*)' F(X) = -1 for -Pi <= X < -Pi/2 ' WRITE(9,*)' F(X) = 1 for -Pi/2 <= X < Pi/2 ' WRITE(9,*)' F(X) = -1 for Pi/2 <= X < Pi ' WRITE(9,*)' F(X) = -1 for Pi <= X ' ELSEIF (FunType.EQ.5) THEN WRITE(9,*)'<5> F(X) = Pi + X for X < -Pi ' WRITE(9,*)' F(X) = -Pi - X for -Pi <= X < -Pi/2' WRITE(9,*)' F(X) = X for -Pi/2 <= X < Pi/2' WRITE(9,*)' F(X) = Pi - X for Pi/2 <= X < Pi ' WRITE(9,*)' F(X) = -Pi + X for Pi <= X ' ELSEIF (FunType.EQ.6) THEN WRITE(9,*)'<6> F(X) = (X+2Pi)^2/4 for X < -Pi' WRITE(9,*)' F(X) = X*X/4 for -Pi <= X < Pi' WRITE(9,*)' F(X) = (X-2Pi)^2/4 for Pi <= X ' ELSEIF (FunType.EQ.7) THEN WRITE(9,*)'<7> F(X) = X/2 + Pi for X < -Pi ' WRITE(9,*)' F(X) = X/2 for -Pi <= X < Pi ' WRITE(9,*)' F(X) = X/2 - Pi for Pi <= X ' ENDIF RETURN END FUNCTION F(X,FunType) REAL X,XL,XR EXTERNAL FX XL=X - 0.000001 XR=X + 0.000001 F=(FX(XL,FunType)+FX(XR,FunType))/2 RETURN END SUBROUTINE COEFFIC(X,Y,A,B,M,M0,N,N0,L) PARAMETER(Pi=3.14159265358979323846) INTEGER J,K,M,M0,N,N0 REAL A,B,L,T,X,Y DIMENSION A(0:361),B(0:361),X(0:361),Y(0:361) DO K=0,M0 A(K)=0 B(K)=0 ENDDO DO K=1,N A(0)=A(0)+Y(K) DO J=1,M0 T=J*Pi*X(K)/L A(J)=A(J)+Y(K)*COS(T) B(J)=B(J)+Y(K)*SIN(T) ENDDO ENDDO DO J=0,M0 A(J)=2*A(J)/N B(J)=2*B(J)/N ENDDO RETURN END FUNCTION TP(A,B,M,M0,L,X) PARAMETER(Pi=3.14159265358979323846) INTEGER J,M,M0 REAL A,B,L,P,X,Z DIMENSION A(0:361),B(0:361) P=A(0)/2 DO J=1,M Z=J*Pi*X/L P=P+A(J)*COS(Z)+B(J)*SIN(Z) ENDDO IF (M.LT.M0) THEN Z=M0*Pi*X/L P=P+A(M0)*COS(Z)/2 ENDIF TP=P RETURN END SUBROUTINE XCOEFFIC(X,Y,A,B,M,M0,N,N0,L) C This subroutine uses labeled DO loop(s). PARAMETER(Pi=3.14159265358979323846) INTEGER J,K,M,M0,N,N0 REAL A,B,L,T,X,Y DIMENSION A(0:361),B(0:361),X(0:361),Y(0:361) DO 10 K=0,M0 A(K)=0 B(K)=0 10 CONTINUE DO 30 K=1,N A(0)=A(0)+Y(K) DO 30 J=1,M0 T=J*Pi*X(K)/L A(J)=A(J)+Y(K)*COS(T) B(J)=B(J)+Y(K)*SIN(T) 20 CONTINUE 30 CONTINUE DO 40 J=0,M0 A(J)=2*A(J)/N B(J)=2*B(J)/N 40 CONTINUE RETURN END FUNCTION XTP(A,B,M,M0,L,X) C This subroutine uses labeled DO loop(s). PARAMETER(Pi=3.14159265358979323846) INTEGER J,M,M0 REAL A,B,L,P,X,Z DIMENSION A(0:361),B(0:361) P=A(0)/2 DO 10 J=1,M Z=J*Pi*X/L P=P+A(J)*COS(Z)+B(J)*SIN(Z) 10 CONTINUE IF (M.LT.M0) THEN Z=M0*Pi*X/L P=P+A(M0)*COS(Z)/2 ENDIF TP=P RETURN END SUBROUTINE PRTPOLY(M,M0,Itype) INTEGER Itype,K,M,M0,U,V IF ((Itype.EQ.1).OR.(Itype.EQ.2)) THEN IF (M.EQ.0) THEN WRITE(9,*)'P(X) = A /2' WRITE(9,*)' 0' ELSEIF (M.EQ.1) THEN WRITE(9,*)'P(X) = A /2 + A COS(X) + B SIN(X)' WRITE(9,*)' 0 1 1' ELSEIF (M.EQ.2) THEN WRITE(9,*)'P(X) = A /2 + A COS( X ) + B SIN( X ) + A COS( +2X ) + B SIN( 2X )' WRITE(9,*)' 0 1 1 2 + 2' ELSE WRITE(9,*)'P(X) = A /2 + A COS( X ) + B SIN( X ) +...+ A C +OS(',M,'X ) + B SIN(',M,'X )' WRITE(9,*)' 0 1 1 ',M +,' ',M ENDIF IF (M.LT.M0) THEN WRITE(9,*)' ' WRITE(9,*)' + A /2 COS(',M0,'X )' WRITE(9,*)' ',M0 ENDIF ENDIF IF (Itype.GT.2) THEN IF (M.EQ.0) THEN WRITE(9,*)'P(X) = A /2' WRITE(9,*)' 0' ELSEIF (M.EQ.1) THEN WRITE(9,*)' Pi*X Pi*X' WRITE(9,*)'P(X) = A /2 + A COS(----) + B SIN(----)' WRITE(9,*)' 0 1 L 1 L' ELSEIF (M.EQ.2) THEN WRITE(9,*)' Pi*X Pi*X + 2Pi*X 2Pi*X' WRITE(9,*)'P(X) = A /2 + A COS(----) + B SIN(----) + A COS +(-----) + B SIN(-----)' WRITE(9,*)' 0 1 L 1 L 2 + L 2 L' ELSE WRITE(9,*)' Pi*X Pi*X + ',M,'Pi*X ',M,'Pi*X' WRITE(9,*)'P(X) = A /2 + A COS(----) + B SIN(----) +...+ A + COS(------) + B SIN(------)' WRITE(9,*)' 0 1 L 1 L ' +,M,' L ',M,' L' ENDIF IF (M.LT.M0) THEN WRITE(9,*)' ' WRITE(9,*)' ',M0,'Pi*X ' WRITE(9,*)' + A /2 COS(------)' WRITE(9,*)' ',M0,' L' ENDIF ENDIF RETURN END SUBROUTINE INPUTFUN(FunType,FunMax) INTEGER FunType,FunMax,I,K REAL Resp DO 10 I=1,20 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' Choose your function:' WRITE(9,*)' ' FunMax=7 DO 20 K=1,FunMax WRITE(9,*)' ' CALL PRINTFUN(K) 20 CONTINUE WRITE(9,*)' ' WRITE(9,*)' SELECT < 1 - ',FunMax,' > ' READ(9,*) Resp FunType=INT(Resp) IF ((FunType.LT.1).OR.(FunType.GT.FunMax)) FunType=1 RETURN END SUBROUTINE GETDATA(F,X,Y,M,M0,N,N0,A0,B0,L,Ftype,Itype, +FunType,FunMax) PARAMETER(Pi=3.14159265358979323846) INTEGER Ftype,FunMax,FunType,Itype,K,M,M0,MAX,N,N0 REAL A0,B0,H,L,X,Y DIMENSION X(0:361),Y(0:361) EXTERNAL F,FX DO 10 I=1,20 WRITE(9,*)' ' 10 CONTINUE WRITE(9,*)' The trigonometric polynomial of degree M is constr +ucted' WRITE(9,*)' ' WRITE(9,*)'P(X) = A /2 + A COS(X) + B SIN(X) +...+ A COS(MX) + B S +IN(MX)' WRITE(9,*)' 0 1 1 M M' WRITE(9,*)' ' WRITE(9,*)'It is the `least squares fit` for the N data points (X +,Y ).' WRITE(9,*)' K + K' WRITE(9,*)'The abscissas are equally spaced: X = A + (B-A)/N for +K=0,..,N' WRITE(9,*)' K' WRITE(9,*)'and the ordinates are periodic, i.e. Y = Y ' WRITE(9,*)' N 0' WRITE(9,*)' ' WRITE(9,*)' Choose the interval type:' WRITE(9,*)' ' WRITE(9,*)' <1> [-Pi,Pi] <2> [0,2Pi] <3> [0,Pi]' WRITE(9,*)' ' WRITE(9,*)' <4> [-L,L] <5> [0,L] <6> [A,B]' WRITE(9,*)' ' WRITE(9,*)' Select <1-6> ' READ(9,*) Itype IF (Itype.EQ.1) THEN A0=-Pi B0=Pi L=Pi ELSEIF (Itype.EQ.2) THEN A0=0 B0=2*Pi L=Pi ELSEIF (Itype.EQ.3) THEN A0=0 B0=Pi L=Pi/2 ELSEIF (Itype.EQ.4) THEN WRITE(9,*)' ' WRITE(9,*)' The interval is [-L,L], Enter L = ' READ(9,*) L A0=-L B0=L ELSEIF (Itype.EQ.5) THEN WRITE(9,*)' ' WRITE(9,*)' The interval is [0,L], Enter L = ' READ(9,*) L A0=0 B0=L L=L/2 ELSEIF (Itype.EQ.6) THEN WRITE(9,*)' ' WRITE(9,*)' The interval is [A,B].' WRITE(9,*)' ' WRITE(9,*)' Enter the left endpoint A = ' READ(9,*) A0 WRITE(9,*)' ' WRITE(9,*)' Enter the right endpoint B = ' READ(9,*) B0 L=(B0-A0)/2 ENDIF WRITE(9,*)' ' WRITE(9,*)' There are N `periodic` data points, enter N = ' READ(9,*) N N0=N WRITE(9,*)' Degree of the trigonometric polynomial is M = ' READ(9,*) M M0=M WRITE(9,*)' ' MAX=INT((N-1)/2) IF (M.GT.MAX) M=MAX IF (INT(N/2).NE.0) THEN M0=M ELSE IF (N.LT.(1+2*M0)) THEN M0=M+1 ELSE M0=M ENDIF ENDIF H=(B0-A0)/N DO 20 K=0,(N-1) X(K)=A0+K*H 20 CONTINUE X(N)=B0 DO 30 I=1,20 WRITE(9,*)' ' 30 CONTINUE IF (M.EQ.M0) THEN WRITE(9,*)'The trigonometric polynomial of degree ',M,' is con +structed.' ENDIF IF (M.LT.M0) THEN WRITE(9,*)'The trigonometric polynomial of degree ',M0,' is co +nstructed.' ENDIF WRITE(9,*)' ' CALL PRTPOLY(M,M0,Itype) WRITE(9,*)' ' WRITE(9,*)'Over the interval [ ',A0,' , ',B0,' ]' WRITE(9,*)' ' IF (Itype.EQ.1 .OR. Itype.EQ.2) THEN WRITE(9,*)'of width 2Pi.' ELSEIF (Itype.EQ.1 .OR. Itype.EQ.2) THEN WRITE(9,*)'of width Pi.' ELSEIF (Itype.EQ.4 .OR. Itype.EQ.5 .OR. Itype.EQ.5) THEN WRITE(9,*)'where L =',L ENDIF WRITE(9,*)' ' WRITE(9,*)' Do you want to enter data points or use Y = F(X ) + ?' WRITE(9,*)' K K' WRITE(9,*)' ' WRITE(9,*)'<1> Enter data points.' WRITE(9,*)' ' WRITE(9,*)'<2> Use Y = F(X ).' WRITE(9,*)' K K' WRITE(9,*)' ' WRITE(9,*)' Select <1 - 2> ' READ(9,*) Ftype WRITE(9,*)' ' IF (Ftype.EQ.2) THEN FunMax=7 CALL INPUTFUN(FunType,FunMax) DO 40 K=0,(N-1) Y(K)=F(X(K),FunType) 40 CONTINUE Y(N)=Y(0) ELSE WRITE(9,*)' ' WRITE(9,*)'ENTER THE ',N,' ORDINATES Y , Y ,...,Y . ' WRITE(9,*)' 0 1 ',N-1 WRITE(9,*)' ' DO 50 K=0,(N-1) WRITE(9,*)'X[',K,' ] = ',X(K),' Y[',K,' ] = ' READ(9,*) Y(K) 50 CONTINUE Y(N)=Y(0) ENDIF DO 60 I=1,20 WRITE(9,*)' ' 60 CONTINUE WRITE(9,*)' Let me think for a while.' WRITE(9,*)' ' WRITE(9,*)'It takes a considerable amount of computing effort to f +ind:' WRITE(9,*)' ' CALL PRTPOLY(M,M0,Itype) IF (Ftype.EQ.2) THEN WRITE(9,*)' ' WRITE(9,*)'For the function:' WRITE(9,*)' ' WRITE(9,*)' ' CALL PRINTFUN(FunType) ENDIF RETURN END SUBROUTINE PRTCOEF(A,B,M,M0,Itype) INTEGER Itype,K,M,M0 REAL A,B DIMENSION A(0:361),B(0:361) WRITE(9,*)' ' IF (Itype.GT.2) THEN WRITE(9,*)' L =',L ELSE WRITE(9,*)' ' ENDIF DO 10 K=0,M WRITE(9,*)' A(',K,') =',A(K) 10 CONTINUE IF (M.LT.M0) WRITE(9,*)' A(',M0,') =',A(M0) DO 20 K=1,M WRITE(9,*)' B(',K,') =',B(K) 20 CONTINUE RETURN END SUBROUTINE RESULTS(TP,A,B,X,Y,XK,M,M0,N,N0,L,Itype,Ftype,FunType) INTEGER Ftype,FunType,Itype,K,M,M0,N,N0 REAL A,B,Err,L,PK,X,XK,Y,YK DIMENSION A(0:361),B(0:361),X(0:361),Y(0:361) EXTERNAL TP DO 10 I=1,20 WRITE(9,*)' ' 10 CONTINUE IF (Ftype.EQ.2) THEN WRITE(9,*)' ' CALL PRINTFUN(FunType) ENDIF WRITE(9,*)' ' CALL PRTPOLY(M,M0,Itype) CALL PRTCOEF(A,B,M,M0,Itype) WRITE(9,*)' ' WRITE(9,*)' K X(K) Y(K) P(X(K)) + Error' DO 20 K=0,N XK=X(K) YK=Y(K) PK=TP(A,B,M,M0,L,XK) Err=YK-PK WRITE(9,1000) K,XK,YK,PK,Err 20 CONTINUE PAUSE RETURN 1000 FORMAT(I4,1X,F15.7,1X,F15.7,1X,F15.7,1X,F15.7) END