PROGRAM IntegrateTest
USE Numeric_Kinds_Module
USE Precision_Module
USE FileUnits_Module
IMPLICIT NONE
INTEGER (KIND=IW_Kind)            :: i, n, Kase, m
REAL (KIND=WP_Kind) a, b, delta, quad, exact, res(1)
REAL (KIND=WP_Kind) , ALLOCATABLE :: xpt(:), wht(:), fun(:), junk(:)

m=2
exact=SQRT(Pi)/EXP(m**2/Four)

OPEN(Unit=Out_Unit,FILE='IntegrateResults.txt')
CALL NumericInfo

WRITE(*,*)'Input Number of Quadrature Points'
READ(*,*) n
ALLOCATE(xpt(n),wht(n),fun(n),junk(n))

b= Six
a=-Six
delta=(b-a)/(n-1)


!Trapezoidal Test Case
CALL Trapezoidal(n, xpt, wht, a, b)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Trapezoidal Quadrature'
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO
junk=wht*COS(m*xpt)*EXP(-xpt*xpt)
junk(1:3)=wht(1:3)*COS(m*xpt(1:3))*EXP(xpt(1:3)*xpt(1:3)) +1.d+12
res=SUM(xpt*xpt)

quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)
   !quad=quad+wht(i)*(one+two*xpt(i)+Three*xpt(i)**2)
   write(88,*)xpt(i),junk(i)
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Trap_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Trap_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


!Simpson Test Case
IF(2*(n/2).ne.n)THEN
   CALL Simpson(n, xpt, wht, a, b)
   WRITE(Out_Unit,*)
   WRITE(Out_Unit,*)
   WRITE(Out_Unit,*) 'Simpson Quadrature'
   DO i=1,n
      WRITE(Out_Unit,*)xpt(i), wht(i)
   ENDDO
ENDIF
quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)
   !quad=quad+wht(i)*(one+two*xpt(i)+Three*xpt(i)**2)
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Simp_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Simp_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
   

! Gauss Legendre Test Case
CALL Glegen(n, xpt, wht, a, b)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Gauss-Legendre Quadrature'
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO
quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Gleg_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'Gleg_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


! Regular Gauss-Hermite
Kase=0
CALL Hermit(n, xpt, wht, Kase)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Gauss-Hermite Quadrature, Kase=',Kase
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO
quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'RGHQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'RGHQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


! Weighted Gauss-Hermite
Kase=1
CALL Hermit(n, xpt, wht, Kase)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Gauss-Hermite Quadrature, Kase=',Kase
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO
quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'WGHQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
!WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'WGHQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


! Regular Gauss-Lagurre
Kase=0
CALL Laguer(n, xpt, wht, Kase)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Gauss-Lagurre Quadrature, Kase=',Kase
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO
quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)*EXP(xpt(i))*Two
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'RGLQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'RGLQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


! Weighted Gauss-Lagurre
Kase=1
CALL Laguer(n, xpt, wht, Kase)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*)
WRITE(Out_Unit,*) 'Gauss-Lagurre Quadrature, Kase=',Kase
DO i=1,n
   WRITE(Out_Unit,*)xpt(i), wht(i)
ENDDO

quad=Zero
DO i=1,n
   quad=quad+wht(i)*COS(m*xpt(i))*EXP(-xpt(i)**2)*Two
ENDDO
WRITE(Out_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'WGLQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad
!WRITE(Msg_Unit,'(A,e15.7,A,e15.7,A,e15.7)')'WGLQ_quad=',quad,'  Exact=', exact, '   Diff=', exact-quad


DEALLOCATE(xpt,wht)


STOP 'Integration Test Complete'
END PROGRAM IntegrateTest