!*****************************************************************************
!
!    Subroutine: Interpolate
!
!    Purpose: Takes function values at the smolyak points as input
!             and (along with subroutine ChebyshevRegression) constucts the
!             coefficients of the polynomials used in the smolyak
!             algorithm
!
!    Author: Ben Malin
!
!    Date: 6/04
!    Revised: 5/31/07
!
!*****************************************************************************

subroutine Interpolate(data0,data2,data3,tempCoeffs1,tempCoeffs2,tempCoeffs3)
  
  ! inputs:
  !   data0 -- Function value at "origin" 
  !   data2 -- Function values at points only used in tensor product
  !   data3 -- Function values at 4th-order polynomial points, excluding origin.
  !
  ! output:
  !   tempCoeffs1 -- Coefficients of 2nd-order polynomial in each dimension
  !   tempCoeffs2 -- Coefficients of tensor-products
  !   tempCoeffs3 -- Coefficients of 4th-order polynomial in each dimension
  !
  ! remarks:
  !  Output will be stored in CoeffsCap0 (etc.) arrays.
  !
  
  use Params
  use Regression
  
  implicit none
 
  !Variable Declaration
  real(prec),intent(in) :: data0
  real(prec),dimension(d-1,d,2**(q-d-1),2**(q-d-1)),intent(in) :: data2
  real(prec),dimension(d,2**(q-d)),intent(in) :: data3
  real(prec),dimension(d,2**(q-d-1)+1),intent(out) :: tempCoeffs1
  real(prec),dimension(d-1,d,2**(q-d-1)+1,2**(q-d-1)+1),intent(out) :: tempCoeffs2
  real(prec),dimension(d,2**(q-d)+1),intent(out) :: tempCoeffs3
  
  real(prec),dimension(d,2**(q-d)+1) :: alldata3
  real(prec),dimension(d,2**(q-d-1)+1) :: alldata1
  real(prec),dimension(d-1,d,9) :: alldata2	!9 == (2**(q-d-1)+1) * (2**(q-d-1)+1) 
  
  real(prec),dimension(d-1,d,9) :: tempCoeffs2vec	!9 == (2**(q-d-1)+1) * (2**(q-d-1)+1)
  
  integer :: i,i1,j,j1
  
  do i = 1,d
     !Function values for fourth order polynomials			
     do j = 1,5		!5 = 2**(q-d)+1 
        if (j < 3) then
           alldata3(i,j) = data3(i,j)
        elseif (j == 3) then
           alldata3(i,j) = data0
        else
           alldata3(i,j) = data3(i,j-1)
        endif
     end do
     
     !Function values for second order, one-dimensional polynomials
     alldata1(i,1) = data3(i,1)
     alldata1(i,2) = data0
     alldata1(i,3) = data3(i,4)
  end do
  
  do i = 1,d-1
     do i1 = i+1,d
        !Function values for 2-dimensional tensor products
        alldata2(i,i1,1) = data2(i,i1,1,1)
        alldata2(i,i1,2) = data3(i,1)
        alldata2(i,i1,3) = data2(i,i1,1,2)
        alldata2(i,i1,4) = data3(i1,1)
        alldata2(i,i1,5) = data0
        alldata2(i,i1,6) = data3(i1,4)
        alldata2(i,i1,7) = data2(i,i1,2,1)
        alldata2(i,i1,8) = data3(i,4)
        alldata2(i,i1,9) = data2(i,i1,2,2)
     end do
  end do
  
  !Calculates Coefficients of Chebyshev Polynomials
  do i = 1,d
     call ChebyshevRegression(1,5,T3,alldata3(i,:),tempCoeffs3(i,:))
     call ChebyshevRegression(1,3,T2,alldata1(i,:),tempCoeffs1(i,:))
  end do
  
  do i = 1,d-1
     do i1 = i+1,d
        call ChebyshevRegression(2,3,T2,alldata2(i,i1,:),tempCoeffs2vec(i,i1,:))
     end do
  end do
  
  do i = 1,d-1
     do i1 = i+1,d
        do j = 1,3
           do j1 = 1,3
              tempCoeffs2(i,i1,j,j1) = tempCoeffs2vec(i,i1,3*(j-1)+j1)
           enddo
        enddo
     enddo
  enddo
  
end subroutine Interpolate
