!     Last change:  DOS  29 Jul 2000   12:05 pm
!                                       *** copyright 2000 ***
! *** filename orthplz.f95              *** John F. Monahan **
!                                       **********************
      PROGRAM ORTHPLZ
! DEMO OF ORTHOGONAL POLYNOMIAL ROUTINES
!  JUST LAGUERRE AND HERMITE HERE
      IMPLICIT NONE
      INTEGER, PARAMETER :: NMAX = 6
      INTEGER, PARAMETER :: NPI = 1024  ! INTEGRATION POINTS
      INTEGER, PARAMETER :: NPE = 25    ! EVALUATION
      REAL, DIMENSION(0:NMAX) :: C,V
      REAL, DIMENSION(2,0:NPE) :: FHAT
      REAL X,FX,W,FNPI
      REAL ORTPE,WEIGHT
      INTEGER KIND,I,J
!
  20  FORMAT(' TEST OF TWO TYPES OF ORTHOGONAL FUNCTIONS'//&
     & ' KIND=5  LAGUERRE POLYNOMIALS'/                    &
     & ' KIND=6  HERMITE POLYNOMIALS'//                    &
     & '    FOURIER COEFFICIENTS')
  21  FORMAT(/' KIND=',I2,4(I4,F12.8)/6X,4(I4,F12.8)/4X,4(I4,F12.8))
  22  FORMAT(/4X,'X',10X,'F(X)',6X,'KIND',I4)
  23  FORMAT(F8.3,F14.8,2X,F14.8)
!
!
      OPEN( UNIT=6, FILE='orthplz.out' )
!
      WRITE(6,20)
      FNPI = REAL(NPI)
      DO KIND = 5,6,1
      C = 0.                  ! INITIALIZE TO ZERO
!                   INTEGRATE TO GET FOURIER COEFS
!                    USING MIDPOINT RULE -- NOTE WEIGHT ON KIND=3
      DO I = 0,NPI
      IF( KIND .EQ. 5 ) THEN     ! LAGUERRE
      X = 32.*REAL(I)/FNPI
      FX = 3.*X*X - 2.*X -4. ! SIMPLE POLYNOMIAL FUNCTION
                     ELSE     ! HERMITE
      X = 4.*REAL(2*I-NPI)/FNPI
      FX = 3.*X*X - 2.*X -4. ! SIMPLE POLYNOMIAL FUNCTION
                     END IF   ! ( KIND = 5 )
      CALL ORTPV(X,KIND,NMAX,V)
      W = WEIGHT(X,KIND)
!                   SIMPSON'S RULE WEIGHTS
      IF( (I .NE. 0).AND.(I .NE. NPI) ) W = W*2.*REAL(1+MOD(I,2))
      DO J = 0,NMAX
      C(J) = C(J) + W*V(J)*FX
      END DO  ! LOOP ON J
      END DO  ! LOOP ON I
!                   FINISH INTEGRATION
!
!                   GET NORMALIZATION CONSTANTS
      CALL ORTPN(KIND,NMAX,V)
!                   NORMALIZE & APPLY INTEGRATION STEP SIZE
      IF( KIND .EQ. 5 ) THEN     ! LAGUERRE
      W = 32./FNPI
                     ELSE     ! HERMITE
      W = 4.*2./FNPI
                     END IF   ! ( KIND = 5 )
      W = W / 3.              ! SIMPSON WEIGHTS
      DO J = 0,NMAX
      C(J) = C(J)*W/V(J)
      END DO
!                   HAVE FOURIER COEFFS
      WRITE(6,21) KIND,(J,C(J),J=0,NMAX)
!                   CHECK IT OUT ON THESE POINTS
      DO I = 0,NPE
      IF( KIND .EQ. 5 ) THEN     ! LAGUERRE
      X = 32.*REAL(I)/REAL(NPE)
                     ELSE     ! HERMITE
      X = 4.*REAL(2*I-NPE)/REAL(NPE)
                     END IF   ! ( KIND = 5 )
      FHAT(KIND-4,I) = ORTPE(X,KIND,NMAX,C)
      END DO  ! LOOP ON I
!
      END DO  ! LOOP ON KIND
!                   LOOK AT LAGUERRE
      KIND = 5
      WRITE(6,22) KIND
      DO I = 0,NPE
      X = 32.*REAL(I)/REAL(NPE)
      FX = 3.*X*X - 2.*X -4.  ! SIMPLE POLYNOMIAL FUNCTION
      WRITE(6,23) X,FX,FHAT(KIND-4,I)
      END DO  ! LOOP ON I
!                   LOOK AT HERMITE
      KIND = 6
      WRITE(6,22) KIND
      DO I = 0,NPE
      X = 4.*REAL(2*I-NPE)/REAL(NPE)
      FX = 3.*X*X - 2.*X -4.  ! SIMPLE POLYNOMIAL FUNCTION
      WRITE(6,23) X,FX,FHAT(KIND-4,I)
      END DO  ! LOOP ON I
      STOP
      END PROGRAM ORTHPLZ
      REAL FUNCTION ORTPE(X,KIND,NMAX,C)
!   USES CLENSHAW'S METHOD TO EVALUATE A FUNCTION GIVEN AS A FINITE
!   EXPANSION IN ORTHOGONAL POLYNOMIALS WITH COEFFICIENTS STORED
!   IN C(0),...,C(NMAX)
!
!  KIND=1   LEGENDRE                         W(X)=1
!  KIND=2   CHEBYSHEV, FIRST                 W(X)=1/SQRT(1-X**2)
!  KIND=3   CHEBYSHEV, SECOND                W(X)=SQRT(1-X**2)
!  KIND=4   STANDARD FOURIER SERIES          W(X)=1.
!  KIND=5   LAGUERRE  (STANDARD, ALPHA=0)    W(X)=EXP(-X)
!  KIND=6   HERMITE                          W(X)=EXP(-X**2)
!  J F MONAHAN    AUGUST 1981   DEPT OF STAT, NCSU, RALEIGH, NC 27650
!    RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      INTEGER, INTENT(IN) :: KIND,NMAX
      REAL, DIMENSION(0:NMAX), INTENT(IN) :: C
      REAL, PARAMETER :: PI = 3.141593       ! 3.14159265358979
      REAL YNEW,YOLD,YOLDR
      INTEGER L,LL,K
!
      YOLD = 0.
      YOLDR = 0.
      ORTPE = 0.
      SELECT CASE (KIND)
!  LEGENDRE
      CASE (1)
      DO L = 0,NMAX
      K = NMAX - L
      YNEW = REAL(2*K+1)*X*YOLD/REAL(K+1) -               &
     &          REAL(K+1)*YOLDR/REAL(K+2) + C(K)
      YOLDR = YOLD
      YOLD =  YNEW
      END DO  ! LOOP ON L
      ORTPE = YNEW
!  CHEBYSHEV, FIRST
      CASE (2)
      DO L = 1,NMAX
      K = NMAX + 1 - L
      YNEW = 2.*X*YOLD - YOLDR + C(K)
      YOLDR = YOLD
      YOLD =  YNEW
      END DO  ! LOOP ON L
      ORTPE = X*YOLD - YOLDR + C(0)
!  CHEBYSHEV, SECOND
      CASE (3)
      DO L = 0,NMAX
      K = NMAX - L
      YNEW = 2.*X*YOLD - YOLDR + C(K)
      YOLDR = YOLD
      YOLD =  YNEW
      END DO  ! LOOP ON L
      ORTPE = YNEW
!  STANDARD FOURIER SERIES
      CASE (4)
      ORTPE = 0.
      DO LL = 1,NMAX
      L = NMAX + 1 - LL
      K = (L+1)/2
      IF( L .EQ. 2*K ) THEN   ! EVEN
      ORTPE = ORTPE + C(L)*SIN(K*PI*X)
                       ELSE
      ORTPE = ORTPE + C(L)*COS(K*PI*X)
                       END IF ! ( L .EQ. 2*K )
      END DO  ! LOOP ON LL
      ORTPE = ORTPE + C(0)
!  LAGUERRE
      CASE (5)
      DO L = 0,NMAX
      K = NMAX - L
      YNEW = (REAL(2*K+1)-X)*YOLD/REAL(K+1) -             &
     &            REAL(K+1)*YOLDR/REAL(K+2) + C(K)
      YOLDR = YOLD
      YOLD =  YNEW
      END DO  ! LOOP ON L
      ORTPE = YNEW
!  HERMITE
      CASE (6)
      DO L = 0,NMAX
      K = NMAX - L
      YNEW = 2.*X*YOLD - 2.*REAL(K+1)*YOLDR + C(K)
      YOLDR = YOLD
      YOLD =  YNEW
      END DO  ! LOOP ON L
      ORTPE = YNEW
      CASE DEFAULT
      END SELECT
      RETURN                                                            
      END FUNCTION ORTPE
      SUBROUTINE ORTPV(X,KIND,NMAX,V)                                    
!  COMPUTES VALUES OF ORTHOGONAL POLYS FROM 0 THRU NMAX
!    IN V(0),,.V(NMAX)
!
!  KIND=1   LEGENDRE                         W(X)=1
!  KIND=2   CHEBYSHEV, FIRST                 W(X)=1/SQRT(1-X**2)
!  KIND=3   CHEBYSHEV, SECOND                W(X)=SQRT(1-X**2)
!  KIND=4   STANDARD FOURIER SERIES          W(X)=1.
!  KIND=5   LAGUERRE  (STANDARD, ALPHA=0)    W(X)=EXP(-X)
!  KIND=6   HERMITE                          W(X)=EXP(-X**2)
!  J F MONAHAN    AUGUST 1981   DEPT OF STAT, NCSU, RALEIGH, NC 27650
!    RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      INTEGER, INTENT(IN) :: KIND, NMAX
      REAL, DIMENSION(0:NMAX), INTENT(OUT) :: V
      REAL, PARAMETER :: PI = 3.141593       ! 3.14159265358979
      INTEGER N,K
      SELECT CASE (KIND)
!  LEGENDRE
      CASE(1)
      V(0) = 1.
      V(1) = X
      DO N = 2,NMAX
      V(N) = ( REAL(2*N-1)*X*V(N-1) - REAL(N-1)*V(N-2) )/REAL(N)
      END DO  ! LOOP ON N
!  CHEBYSHEV, FIRST
      CASE (2)
      V(0) = 1.
      V(1) = X
      DO N = 2,NMAX
      V(N) = 2.*X*V(N-1) - V(N-2)
      END DO  ! LOOP ON N
!  CHEBYSHEV, SECOND
      CASE (3)
      V(0) = 1.
      V(1) = 2.*X
      DO N = 2,NMAX
      V(N) = 2.*X*V(N-1) - V(N-2)
      END DO  ! LOOP ON N
!  STANDARD FOURIER SERIES
      CASE (4)
      V(0) = 1.
      DO N = 1,NMAX
      K = (N+1)/2
      IF( N .EQ. 2*K ) THEN   ! EVEN
      V(N) = SIN(K*PI*X)
                       ELSE
      V(N) = COS(K*PI*X)
                       END IF ! ( N .EQ. 2*K )
      END DO  ! LOOP ON N
!  LAGUERRE
      CASE (5)
      V(0) = 1.
      V(1) = 1.-X
      DO N = 2,NMAX
      V(N) = ( (REAL(2*N-1)-X)*V(N-1) - REAL(N-1)*V(N-2) ) / REAL(N)
      END DO  ! LOOP ON N
!  HERMITE
      CASE (6)
      V(0) = 1.
      V(1) = 2.*X
      DO N = 2,NMAX
      V(N) = 2.*( X*V(N-1) - REAL(N-1)*V(N-2) )
      END DO  ! LOOP ON N
      CASE DEFAULT
      END SELECT
      RETURN
      END SUBROUTINE ORTPV
      REAL FUNCTION WEIGHT(X,KIND)
!  COMPUTES WEIGHT FUNCTION FOR ORTHOGONAL POLYS
!  KIND=1   LEGENDRE                         W(X)=1
!  KIND=2   CHEBYSHEV, FIRST                 W(X)=1/SQRT(1-X**2)
!  KIND=3   CHEBYSHEV, SECOND                W(X)=SQRT(1-X**2)
!  KIND=4   STANDARD FOURIER SERIES          W(X)=1.
!  KIND=5   LAGUERRE  (STANDARD, ALPHA=0)    W(X)=EXP(-X)
!  KIND=6   HERMITE                          W(X)=EXP(-X**2)
!  J F MONAHAN    AUGUST 1981   DEPT OF STAT, NCSU, RALEIGH, NC 27650
!    RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      INTEGER, INTENT(IN) :: KIND
      SELECT CASE (KIND)
!  LEGENDRE
      CASE (1)
      WEIGHT = 1.
!  CHEBYSHEV, FIRST
      CASE (2)
      WEIGHT = 1./ SQRT(1.-X*X)
!  CHEBYSHEV, SECOND
      CASE (3)
      WEIGHT = SQRT(1.-X*X)
!  STANDARD FOURIER SERIES
      CASE (4)
      WEIGHT = 1.
!  LAGUERRE
      CASE (5)
      IF( X .GT. 85. ) THEN
      WEIGHT = 0.             ! UNDERFLOW
                       ELSE
      WEIGHT = EXP(-X)        ! OK
                       END IF ! ( X .GT. 85. )
!  HERMITE
      CASE (6)
      IF( X .GT. 9. ) THEN
      WEIGHT = 0.             ! UNDERFLOW
                      ELSE
      WEIGHT = EXP(-X*X)      ! OK
                      END IF  ! ( X .GT. 9. )
      CASE DEFAULT
      END SELECT
      RETURN                                                            
      END FUNCTION WEIGHT
      SUBROUTINE ORTPN(KIND,NMAX,V)                                      
!  COMPUTES NORMALIZATION CONSTANTS FOR ORTHOGONAL POLYS FROM
!    0 THRU NMAX AND STORES THEM IN V(0),...,V(NMAX)
!  KIND=1   LEGENDRE                         W(X)=1
!  KIND=2   CHEBYSHEV, FIRST                 W(X)=1/SQRT(1-X**2)
!  KIND=3   CHEBYSHEV, SECOND                W(X)=SQRT(1-X**2)
!  KIND=4   STANDARD FOURIER SERIES          W(X)=1.
!  KIND=5   LAGUERRE  (STANDARD, ALPHA=0)    W(X)=EXP(-X)
!  KIND=6   HERMITE                          W(X)=EXP(-X**2)
!  J F MONAHAN    AUGUST 1981   DEPT OF STAT, NCSU, RALEIGH, NC 27650
!     RECODED FEBRUARY 2000, APRIL 2000 FOR FORTRAN 95
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: KIND
      INTEGER, INTENT(IN) :: NMAX
      REAL, DIMENSION(0:NMAX), INTENT(OUT) :: V
      REAL, PARAMETER :: PI = 3.141593       ! 3.14159265358979
      REAL, PARAMETER :: PIO2 = 1.570796     ! 1.5707963267949
      REAL, PARAMETER :: SQRTPI = 1.772454   ! 1.77245385090552
      INTEGER N
      SELECT CASE (KIND)
!  LEGENDRE
      CASE (1)
      V(0) = 2.
      DO N = 1,NMAX
      V(N) = 2./REAL(2*N+1)
      END DO  ! LOOP ON N
!  CHEBYSHEV, FIRST
      CASE (2)
      V(0) = PI
      DO N = 1,NMAX
      V(N) = PIO2
      END DO  ! LOOP ON N
!  CHEBYSHEV, SECOND
      CASE (3)
      V(0) = PIO2
      DO N = 1,NMAX
      V(N) = PIO2
      END DO  ! LOOP ON N
!  FOURIER
      CASE (4)
      V(0) = 2.
      DO N = 1,NMAX
      V(N) = 1.
      END DO  ! LOOP ON N
!  LAGUERRE
      CASE (5)
      V(0) = 1.
      DO N = 1,NMAX
      V(N) = 1.
      END DO  ! LOOP ON N
!  HERMITE                                                              
      CASE (6)
      V(0) = SQRTPI
      DO N = 1,NMAX
      V(N) = 2.*N*V(N-1)
      END DO  ! LOOP ON N
      CASE DEFAULT
      END SELECT
      RETURN
      END SUBROUTINE ORTPN
! *** end of filename orthplz.f95       *********************