!**********************************************************************************
!
!	FUNCTION:  ChebValue
!
!	PURPOSE: Given the state variables, returns the value of the 
!                 policy function.
!			 
!	AUTHOR:	 Ben Malin 
!
!	DATE:    6/2004
!       REVISED: 6/2007
!
!**********************************************************************************

function ChebValue(Ttemp,tempCoeffs1,tempCoeffs2,tempCoeffs3,ChebValueOrigin)
  
  ! inputs:
  !    Ttemp -- Chebyshev polynomials of state variables
  !    tempCoeffs1 -- coefficients of policy function
  !    tempCoeffs2 -- coefficients of policy function
  !    tempCoeffs3 -- coefficients of policy function
  !    ChebValueOrigin -- coefficient of policy function
  !
  ! output:
  ! 
  ! remarks:
  !    The relevant information concerning the state is included in Ttemp.
  !    The formula for constructing the value is given in Krueger/Kubler (2003)
  !     pg. 10.
  !

  use Params
  
  implicit none
  
  !Variable Declaration
  real(prec),dimension(2**(q-d)+1,d),intent(in) :: Ttemp
  real(prec),dimension(d,2**(q-d-1)+1),intent(in) :: tempCoeffs1		
  real(prec),dimension(d-1,d,2**(q-d-1)+1,2**(q-d-1)+1),intent(in) :: tempCoeffs2		
  real(prec),dimension(d,2**(q-d)+1),intent(in) :: tempCoeffs3		
  real(prec),intent(in) :: ChebValueOrigin

  real(prec) :: ChebValue
  real(prec),dimension(d) :: ChebValue1, ChebValue3
  real(prec),dimension(d-1,d) :: ChebValue2
  integer :: i,j,i1,j1
  
  !Construct 4th order polynomials in each of d dimensions
  ChebValue3 = 0.0
  do i = 1,d
     do j = 1,5		!Note that 5 = 2**(q-d) + 1, since (q-d) = 2
        ChebValue3(i) = tempCoeffs3(i,j)*Ttemp(j,i) + ChebValue3(i) 
     end do
  end do
  
  !Construct 2nd order polynomials in each of d dimensions
  ChebValue1 = 0.0
  do i = 1,d
     do j = 1,3		!Note that 3 = 2**(q-d-1) + 1, since (q-d) = 2
        ChebValue1(i) = tempCoeffs1(i,j)*Ttemp(j,i) + ChebValue1(i) 
     end do
  end do
  
  !Construct 2-dimensional tensor products
  ChebValue2 = 0.0
  do i = 1,d-1
     do i1 = i+1,d
        do j = 1,3		!Note that 3 = 2**(q-d) + 1, since (q-d) = 2
           do j1 = 1,3		!Note that 3 = 2**(q-d) + 1, since (q-d) = 2
              ChebValue2(i,i1) = tempCoeffs2(i,i1,j,j1)*Ttemp(j,i)*Ttemp(j1,i1) + ChebValue2(i,i1)
           enddo
        enddo
     enddo
  enddo
  
  !Constuct ChebValue ... Using formula on pg. 10 of Kubler/Krueger
  ChebValue = 0.0
  do i = 1,d
     ChebValue = ChebValue3(i) + ChebValue
  enddo
  
  do i = 1,d
     ChebValue = (-1.0)*coc1*ChebValue1(i) + ChebValue
  enddo
  
  do i = 1,d-1
     do i1 = i+1,d
        ChebValue = ChebValue2(i,i1) + ChebValue
     enddo
  enddo
  
  if (d >= 3) then
     ChebValue = coc2*ChebValueOrigin + ChebValue
  endif

end function ChebValue
