!*********************************************************************
!
!    Subroutine:  Grid
!
!    Purpose:    Initialize parameters, smolyak points, integration 
!                weights, etc.
!
!    Author:   Ben Malin
!
!    Date:   6/04
!    Revised: 9/06/07
!*********************************************************************

subroutine Grid
  
  ! inputs:
  !
  ! output:
  ! 
  ! remarks:
  !
  
  use Params
  
  implicit none
  
 !Variable Declaration
  integer :: h,i,j,i1,j1 
  
 !Initialize run-time variable
  policy_runtime = 0.0
  
 !Construct weights used in weighting of low-dimensional cheb polynomials for smolyak interpolation
  if (d == 2) then
     coc1 = 1
     coc2 = 1
  else 
     coc1 = d-1
     coc2 = (d-1)*(d-2)/2
  endif
  
 !Specify lower and upper bounds for technological shock, e^(z)
  zmin = -0.8 * sigma / (1.0 - rho)    !Note: this may truncate innovations.
  zmax = 0.8 * sigma / (1.0 - rho)
  
 !Compute deterministic steady state (no adjustment costs), 
 !and specify lower and upper bounds for capital stock.
  kss = 1.0	! Note: kss = ((1.0/betah-1.0)/(A_tfp*alpha))**(1.0/(alpha-1.0)) = 1.0
 
  do j = 1,N
     kmin(j) = 0.8*kss 
     kmax(j)  = 1.2*kss 
  end do
!  do i = 1,N
!     kmin(i) = 0.45*kss - (i-1.0)/(N-1.0)*0.1*kss	![kmin, kmax] centered around [0.4, 1.6] * kss 
!     kmax(i) = 1.55*kss + (i-1.0)/(N-1.0)*0.1*kss
!  end do
!kmin(1) = 0.44*kss
!kmax(1) = 1.56*kss
!kmin(2) = 0.40
!kmax(2) = 1.60
!kmin(N) = 0.36*kss
!kmax(N) = 1.64*kss
!print*, kmin
!print*, kmax
!pause

 !Construct the grid for knorm, znorm -- Used in DecRule subroutine
  do i = 1,nk
     do j = 1,N
        knorm(i,j) = kmin(j) + ((i-1.0)/(nk-1.0))*(kmax(j)-kmin(j))
     end do
  end do

  do i = 1,nz
     do j = 1,N
        znorm(i,j) = zmin + ((i-1.0)/(nz-1.0))*(zmax-zmin)
     end do
  end do
  
 !Defines grid of length (2**(q-d)+1) (on [-1,1],[zmin,zmax] and [kmin,kmax])
 !and defines Chebyshev polynomial at those points.
  do i = 1, 2**(q-d)+1
     x(i) = -cos (real( i - 1 ,prec)/real( 2**(q-d) ,prec)*pie)	!Sets grid of Chebyshev extrema on [-1,1]
     
     do h = 1,N
        k(i,h) = (x(i) + 1.0)*(kmax(h) - kmin(h))/(2.0) + kmin(h)  !Chebyshev extrema on [kmin,kmax] x N
        z(i,h) = (x(i) + 1.0)*(zmax - zmin)/(2.0) + zmin  !Chebyshev extrema on [zmin,zmax] x N
        
        T3(i,1,h) = 1.0		!Calculates Chebyshev polynomial at x(i)
        T3(i,1,h+N) = 1.0
        T3(i,2,h) = x(i)
        T3(i,2,h+N) = x(i)
        do j = 3, (2**(q-d)+1)
           T3(i,j,h) = 2.0 * x(i) * T3(i,j-1,h) - T3(i,j-2,h)
           T3(i,j,h+N) = 2.0 * x(i) * T3(i,j-1,h+N) - T3(i,j-2,h+N)
        end do
        
        if (i == 1 .or. i == 3 .or. i == 5) then	!Using q = d+2
           i1 = (i-1)/2 + 1
           T2(i1,1,h) = 1.0
           T2(i1,1,h+N) = 1.0
           T2(i1,2,h) = x(i)
           T2(i1,2,h+N) = x(i)
           T2(i1,3,h) = 2.0 * x(i) * T2(i1,2,h) - T2(i1,1,h)
           T2(i1,3,h+N) = 2.0 * x(i) * T2(i1,2,h+N) - T2(i1,1,h+N)
        endif
     end do
  end do
  
  !**************************************************************************
  ! Defines Smolyak points on ([kmin,kmax] x N) Union ([zmin,zmax] x N) grid
  !**************************************************************************
  !The origin
  do h = 1,N
     point0(h) = k(3,h)		!Note: 3 = 2**(q-d-1)+1
     point0(N+h) = z(3,h)
  enddo
  
  !Extrema (excluding origin) of 4th degree polynomial in ith direction, holding other directions constant
  do i = 1, d
     do j = 1,2**(q-d)
        do h = 1,d
           if (i <= N) then
              if (i == h .and. j <= 2**(q-d-1)) then
                 points3(i,j,h) = k(j,h)
              elseif (i == h .and. j > 2**(q-d-1)) then
                 points3(i,j,h) = k(j+1,h)
              else
                 points3(i,j,h) = point0(h)
              endif
           else
              if (i == h .and. j <= 2**(q-d-1)) then
                 points3(i,j,h) = z(j,h-N)
              elseif (i == h .and. j > 2**(q-d-1)) then
                 points3(i,j,h) = z(j+1,h-N)
              else 
                 points3(i,j,h) = point0(h)
              endif
           endif
        enddo
     enddo
  enddo

  !Points used in 2-dimensional tensor product (excluding those defined above)
  do i = 1, d-1
     do i1 = i+1,d
        do j = 1,2**(q-d-1)
           do j1 = 1,2**(q-d-1)
              do h = 1,d
                 if (i <= N .and. i1 <= N) then
                    if (h == i .and. j == 1) then	!Here, I'm using j == 1 since q = d+2
                       points2(i,i1,j,j1,h) = k(1,h)	!Also, q = d+2 allows me to use k(1,h)
                    elseif (h == i .and. j == 2) then
                       points2(i,i1,j,j1,h) = k(5,h)	!k(5,h), etc.
                    elseif (h == i1 .and. j1 == 1) then
                       points2(i,i1,j,j1,h) = k(1,h)
                    elseif (h == i1 .and. j1 == 2) then
                       points2(i,i1,j,j1,h) = k(5,h)
                    elseif (h <= N) then
                       points2(i,i1,j,j1,h) = k(3,h)
                    else
                       points2(i,i1,j,j1,h) = z(3,h-N)
                    endif
                 elseif (i <= N .and. i1 > N) then
                    if (h == i .and. j == 1) then	!Here, I'm using j == 1 since q = d+2
                       points2(i,i1,j,j1,h) = k(1,h)	!Also, q = d+2 allows me to use k(1,h)
                    elseif (h == i .and. j == 2) then
                       points2(i,i1,j,j1,h) = k(5,h)	!k(5,h), etc.
                    elseif (h == i1 .and. j1 == 1) then
                       points2(i,i1,j,j1,h) = z(1,h-N)
                    elseif (h == i1 .and. j1 == 2) then
                       points2(i,i1,j,j1,h) = z(5,h-N)
                    elseif (h <= N) then
                       points2(i,i1,j,j1,h) = k(3,h)
                    else 
                       points2(i,i1,j,j1,h) = z(3,h-N)
                    endif
                 else
                    if (h == i .and. j == 1) then	!Here, I'm using j == 1 since q = d+2
                       points2(i,i1,j,j1,h) = z(1,h-N)	!Also, q = d+2 allows me to use k(1,h)
                    elseif (h == i .and. j == 2) then
                       points2(i,i1,j,j1,h) = z(5,h-N)	!k(5,h), etc.
                    elseif (h == i1 .and. j1 == 1) then
                       points2(i,i1,j,j1,h) = z(1,h-N)
                    elseif (h == i1 .and. j1 == 2) then
                       points2(i,i1,j,j1,h) = z(5,h-N)
                    elseif (h <= N) then
                       points2(i,i1,j,j1,h) = k(3,h)
                    else 
                       points2(i,i1,j,j1,h) = z(3,h-N)
                    endif
                 endif
              enddo
           enddo
        enddo
     enddo
  enddo
  
  !*******************************************************
  !Constructs grid of points at which integrand will be 
  !evaluated.  Also assigns values to weights as specified
  !by formula (7.5.11).
  !*******************************************************
 
  !Weights used in Integration
  V_mono = pie**(0.5*(N+1))
  A_mono = 2.0 * V_mono / (real(N+1,prec) + 2.0)
  B_mono = (4.0 - real(N+1,prec)) * V_mono / (2.0 * (real(N+1,prec) + 2.0)*(real(N+1,prec)+2.0))
  C_mono = V_mono / (real(N+1,prec) + 2.0) / (real(N+1,prec) + 2.0)
  
  !Construct points used in integration
  e_mono = 0.0
  do i = 1,N+1
     e_mono(i,i) = 1.0
  enddo
  
  r_mono = (1.0 + 0.5 * real(N+1,prec))**(.5)
  s_mono = (0.5 + 0.25 * real(N+1,prec))**(.5)
  
  !The Origin
  point0_mono = 0.0
  
  !"2nd term" points
  do i = 1,N+1
     points1_mono(i,:) = r_mono * e_mono(i,:)
  enddo
  
  !"3rd term" points
  do i = 1, N
     do j = i+1, N+1
        points2_mono(i,j,1,:) = s_mono * e_mono(i,:) + s_mono * e_mono(j,:)
        points2_mono(i,j,2,:) = s_mono * e_mono(i,:) - s_mono * e_mono(j,:)
     enddo
  enddo

end subroutine Grid
