!Module to do lagrange interpolation

module wei_lagrange
use wei_process
implicit none


contains 
subroutine lagrange_init(index_map,frac,zloc)
integer,intent (out) :: index_map(:,:)
real,intent (out) :: frac(:,:)
real,intent (in) :: zloc(:)
real, allocatable :: dist(:),xpos(:)
integer :: iloc,n1in,icoef,ncoef,i1
real    :: loc,mymax,mymin
integer :: myloc





allocate(xpos(size(index_map,1)),dist(size(zloc)))
n1in=size(zloc);
ncoef=size(frac,1)

do iloc=1,pro%az%n
  loc=pro%az%o+pro%az%d*(iloc-1)
  dist=abs(loc-zloc)
  mymax=maxval(dist)
  if(loc > zloc(size(zloc))) then
    frac(:,iloc)=0.
    index_map(:,iloc)=size(zloc)
    if(loc - zloc(size(zloc)) < pro%az%d) frac(1,iloc)=0.
  else
    do icoef=1,ncoef
      mymin=mymax+3.
      do i1=1,n1in 
        if(dist(i1)< mymin)  then
          mymin=dist(i1)
          myloc=i1
        end if
      end do
      index_map(icoef,iloc)=myloc
      xpos(icoef)=zloc(myloc)
      dist(myloc)=mymax+1
    end do
    if(.not. find_lagrange_coefs(loc,xpos,frac(:,iloc))) &
      call erexit("trouble finding lagrange coefs")
  end if
end do
deallocate(dist,xpos)

end subroutine lagrange_init



logical function find_lagrange_coefs(x_out,x_in,wts)
real :: x_out
real :: x_in(:)
real :: wts(:)
integer n1,i,j

if(size(x_in)/=size(wts)) then;
  write(0,*) "input coef size and wt size not same"
  find_lagrange_coefs=.false.
  return
end if
  find_lagrange_coefs=.true.
  n1=size(wts)
  do i=1,n1
    wts(i)=1
    do j=1,n1
      if(i/=j) wts(i)=wts(i)*(x_out-x_in(j))/(x_in(i)-x_in(j))
    end do
  end do
end function



subroutine lagrange_interp(adj,iw,k,iz,image_buf,image_out,image_temp,tag)
logical :: adj
integer :: iw,k,iz,ij
integer :: i,ib,stat
integer, pointer :: status(:)
integer, pointer :: iloc(:)
real, pointer :: frac(:)
complex :: image_buf(:,:,:,:,:),image_out(:,:,:,:),image_temp(:,:,:,:)
logical :: found
character(len=*) :: tag

allocate(frac(size(pro%z_samp(iw)%iloc,1)))
allocate(iloc(size(pro%z_samp(iw)%iloc,1)))

iloc=pro%z_samp(iw)%iloc(:,iz)
frac=pro%z_samp(iw)%frac(:,iz)
status=>pro%z_samp(iw)%status

do i=1,size(iloc)
  if(abs(frac(i)) > .001) then
    if(adj) then
      !find the image component
      ij=iloc(i)-(k-1)*pro%nzs
      if(ij >0 .and. ij <= size(image_buf,5)) then !slice is in image_buf
        image_temp=image_buf(:,:,:,:,ij)
      else  if(status(iloc(i))>0) then !slice is in variable memory buffer
        image_temp=pro%v_i_buf(:,:,:,:,status(iloc(i)))
      else if(status(iloc(i)) <0) then !slice is on disk
        stat=sseek_block(tag,(-status(iloc(i))-1),size(image_temp)*8,0)
        call sreed(tag,image_temp,size(image_temp)*8)
      else
       call seperr("error error status set to 0 ")
      end if
      image_out=image_out+image_temp*frac(i)
    else 
      call seperr("unimplemented modeling variable sampling option")
   end if
  end if
end do

deallocate(iloc,frac)
end subroutine lagrange_interp






end module wei_lagrange
