      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