!*****************************************************************************
!
!    Subroutine: MonIntegrate
!
!    Purpose: Calculates an integral using monomial formula (7.5.11) in Judd
!              (pg. 275).
!
!    Author: Ben Malin
!
!    Date: 6/04
!    Revised: 6/07
!
!*****************************************************************************

subroutine MonIntegrate(nc, IntResult, check)

  ! inputs:
  !   nc -- Index of country.
  !   check -- Denotes whether MCIntegrate is called in policy function iteration
  !             or in simulation.
  !
  ! output:
  !   IntResult -- The value of the integral. 
  !
  ! remarks:
  !   The subroutines constructs the integral on the RHS of the Euler Equation
  !   for country nc.
  !
  
  use Params
  
  implicit none

  !Declaration of Variables
  integer,intent(in) :: nc
  real(prec),intent(out) :: IntResult 
  integer,intent(in) :: check
  integer :: h,i,j
  
  !Store values of integrand
  real(prec) :: Integrand0				! Integrand at the origin
  real(prec),dimension(N+1,2) :: Integrand1		! Integrand at 2nd term points
  real(prec),dimension(N, N+1, 4) :: Integrand2	! Integrand at 3rd term points
  
  !Returns Integrand on RHS of Euler Equation at specified points
  real(prec), external :: ExpectFOC
  real(prec), external :: ExpectFOC_Simul
  
  if (check == 1) then
     !Calculate values of integrand at specified points
     !The origin
     Integrand0 = ExpectFOC(N, point0_mono, nc)
     
     !2nd term points
     do i = 1, N+1
        Integrand1(i,1) = ExpectFOC(N, points1_mono(i,:),nc)
        Integrand1(i,2) = ExpectFOC(N, -1.0 * points1_mono(i,:),nc)
     enddo
     
     !3rd term points
     do i = 1, N
        do j = i+1, N+1
           Integrand2(i,j,1) = ExpectFOC(N, points2_mono(i,j,1,:),nc)
           Integrand2(i,j,2) = ExpectFOC(N, points2_mono(i,j,2,:),nc)
           Integrand2(i,j,3) = ExpectFOC(N, -1.0 * points2_mono(i,j,1,:),nc)
           Integrand2(i,j,4) = ExpectFOC(N, -1.0 * points2_mono(i,j,2,:),nc)
        enddo
     enddo
     
  elseif (check == 3) then
     !Calculate values of integrand at specified points
     !The origin
     Integrand0 = ExpectFOC_Simul(N, point0_mono, nc)
     
     !2nd term points
     do i = 1, N+1
        Integrand1(i,1) = ExpectFOC_Simul(N, points1_mono(i,:),nc)
        Integrand1(i,2) = ExpectFOC_Simul(N, -1.0 * points1_mono(i,:),nc)
     enddo
     
     !3rd term points
     do i = 1, N 
        do j = i+1, N+1
           Integrand2(i,j,1) = ExpectFOC_Simul(N, points2_mono(i,j,1,:),nc)
           Integrand2(i,j,2) = ExpectFOC_Simul(N, points2_mono(i,j,2,:),nc)
           Integrand2(i,j,3) = ExpectFOC_Simul(N, -1.0 * points2_mono(i,j,1,:),nc)
           Integrand2(i,j,4) = ExpectFOC_Simul(N, -1.0 * points2_mono(i,j,2,:),nc)
        enddo
     enddo
  endif
  
  !Weight and Sum values of integrand  -- See formula (7.5.11) on pg. 275 of Judd 
  IntResult = 0.0
  
  !The Origin
  IntResult = A_mono * Integrand0 + IntResult
  
  !"2nd term" points
  do i = 1, N+1
     IntResult = B_mono * (Integrand1(i,1) + Integrand1(i,2)) + IntResult
  enddo
  
  !"3rd term" points
  do i = 1, N
     do j = i + 1, N+1
        do h = 1,4
           IntResult = C_mono * Integrand2(i,j,h) + IntResult
        enddo
     enddo
  enddo
  
end subroutine MonIntegrate
