program bivar_prb
!
!***********************************************************************
!
!! BIVAR_PRB runs tests on the BIVAR bivariate scattered data routines.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BIVAR_PRB'
  write ( *, '(a)' ) '  Tests for BIVAR, which interpolates scattered '
  write ( *, '(a)' ) '  (X,Y) data.'

  call test01
  call test02

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'BIVAR_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!***********************************************************************
!
!! TEST01 tests IDBVIP, interpolation at scattered output points.
!
  implicit none
!
  integer, parameter :: nd = 20
  integer, parameter :: ni = 5
!
  integer, parameter :: niwk = 31 * nd + ni
  integer, parameter :: nwk = 8 * nd
!
  real error
  integer i
  integer iwk(niwk)
  integer md
  real temp
  real wk(nwk)
  real xd(nd)
  real xi(ni)
  real yd(nd)
  real yi(ni)
  real zd(nd)
  real zi(ni)
!
!  X coordinates of the data points.
!
  data xd / &
     0.632000E+00, 0.295109E+00, 0.578769E+00, 0.140772E+00, 0.985549E+00, &
     0.209717E+00, 0.893022E+00, 0.928088E+00, 0.766803E+00, 0.819351E+00, &
     0.959191E+00, 0.516312E+00, 0.052801E+00, 0.626174E+00, 0.275320E+00, &
     0.659645E+00, 0.396438E+00, 0.204824E+00, 0.906341E+00, 0.880356E+00/
!
!  Y coordinates of the data points.
!
  data yd / &
     0.918675E+00, 0.821677E+00, 0.962992E+00, 0.477774E+00, 0.056257E+00, &
     0.654656E+00, 0.377358E+00, 0.263626E+00, 0.706323E+00, 0.584658E+00, & 
     0.154619E+00, 0.986449E+00, 0.197552E+00, 0.924616E+00, 0.788099E+00, &
     0.886828E+00, 0.945136E+00, 0.643341E+00, 0.335303E+00, 0.416050E+00/
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  IDBVIP accepts scattered Z(X,Y) data, and interpolates '
  write ( *, '(a)' ) '  values of Z at new points.'
  write ( *, '(a)' ) ' '
!
!  Set up the Z data.
!
  do i = 1, nd
    zd(i) = exp ( xd(i) ) * sin ( yd(i) )
  end do
!
!  Set the points at which we want interpolated values.
!
  do i = 1, ni
    xi(i) = real(i-1) / real(ni-1)
    yi(i) = 1.0E+00 - xi(i)
  end do
!
!  Call IDBVIP to estimate the value of Z at the points (XI,YI).
!
  md = 1
 
  call idbvip ( md, nd, xd, yd, zd, ni, xi, yi, zi, iwk, wk )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   X     Y    Z    Exact Z    Error'
  write ( *, '(a)' ) ' '

  do i = 1, ni
    temp = exp ( xi(i) ) * sin ( yi(i) )
    error = zi(i) - temp
    write ( *, '(5g14.6)' ) xi(i), yi(i), zi(i), temp, error
  end do
 
  return
end
subroutine test02
!
!***********************************************************************
!
!! TEST02 tests IDSFFT, interpolation on a regular grid.
!
  implicit none
!
  integer, parameter :: nd = 20
  integer, parameter :: nx = 5
  integer, parameter :: ny = 4
!
  integer, parameter :: nz = nx
  integer, parameter :: niwk = 31*nd+nx*ny
  integer, parameter :: nwk = 6*nd
!
  real error
  integer i
  integer iwk(niwk)
  integer j
  integer md
  real temp
  real wk(nwk)
  real xd(nd)
  real xi(nx)
  real yd(nd)
  real yi(ny)
  real zd(nd)
  real zi(nz,ny)
!
!  X coordinates of the data points.
!
  data xd / &
     0.632000E+00, 0.295109E+00, 0.578769E+00, 0.140772E+00, 0.985549E+00, &
     0.209717E+00, 0.893022E+00, 0.928088E+00, 0.766803E+00, 0.819351E+00, &
     0.959191E+00, 0.516312E+00, 0.052801E+00, 0.626174E+00, 0.275320E+00, &
     0.659645E+00, 0.396438E+00, 0.204824E+00, 0.906341E+00, 0.880356E+00/
!
!  Y coordinates of the data points.
!
  data yd / &
      0.918675E+00, 0.821677E+00, 0.962992E+00, 0.477774E+00, 0.056257E+00, &
      0.654656E+00, 0.377358E+00, 0.263626E+00, 0.706323E+00, 0.584658E+00, &
      0.154619E+00, 0.986449E+00, 0.197552E+00, 0.924616E+00, 0.788099E+00, &
      0.886828E+00, 0.945136E+00, 0.643341E+00, 0.335303E+00, 0.416050E+00/
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  IDSFFT accepts scattered Z(X,Y) data, and interpolates'
  write ( *, '(a)' ) '  values of Z on a regular grid of points.'
  write ( *, '(a)' ) ' '
!
!  Set up the Z data
!
  do i = 1, nd
    zd(i) = exp ( xd(i) ) * sin ( yd(i) )
  end do
!
!  Set up the points at which we want interpolated values.
!
  do i = 1, nx
    xi(i) = real(i-1) / real(nx-1)
  end do
 
  do j = 1, ny
    yi(j) = real(j-1) / real(ny-1)
  end do
!
!  Call IDSFFT to estimate the value of Z.
!
  md = 1
 
  call idsfft ( md, nd, xd, yd, zd, nx, ny, nz, xi, yi, zi, iwk, wk )
 
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '   X     Y    Z    Exact Z    Error'
  write ( *, '(a)' ) ' '

  do i = 1, nx
    do j = 1, ny
      temp = exp ( xi(i) ) * sin ( yi(j) )
      error = zi(i,j) - temp
      write ( *, '(5g14.6)' ) xi(i), yi(j), zi(i,j), temp, error
    end do
    write ( *, * ) ' '
  end do
 
  return
end