program main
! Last updated 9/2/2007
!
! prints out S (longitudinal 2-pt third-order correlation) and Sf 
! (longitudinal 2-pt third-order correlation of filtered velocity, with a
! tensor product top hat filter) as functions of separation for epsilon=1
! user is asked to enter the separation, and the quantities
! are printed to stdout
!
  use evalTijk
  implicit none
  double precision :: coeff(5),r_(3),s_(3)
  integer :: x1,x2,x3,y1,y2,y3,z1,z2,z3
  double precision :: S,Sf,sep,g5pts(-2:2),g5wts(-2:2),weight

   coeff = (/0.884,-2.692,-6.099,-5.853,14.760/) ! coefficients for model
  call read_basis("basisdata.txt")
!  call print_basis
  call build_composite(coeff,.true.)
!  call print_composite

  write (*,*) "enter separation:"
  read (*,*) sep
  r_ = (/0.,0.,0./)
  s_ = (/sep,0.,0./)
  S = eval_composite(1,1,1,r_,s_,.true.)

!integrate Tijk using 5-point Gauss quadrature to find Sf
  Sf = 0.
  g5pts = (/-.906180,-.538469,0,.538469,.906180/)
  g5wts = (/.236927,.478629,.568889,.478629,.236927/)
  g5pts = g5pts/2 !adjust points for filter width=1
  g5wts = g5wts/2 !adjust weights for filter width=1
  do x1 = -2,2
   do x2 = -2,2
    do x3 = -2,2
     do y1 = -2,2
      do y2 = -2,2
       do y3 = -2,2
        do z1 = -2,2
         do z2 = -2,2
          do z3 = -2,2
           r_(1) = g5pts(y1)-g5pts(x1)+sep
           r_(2) = g5pts(y2)-g5pts(x2)
           r_(3) = g5pts(y3)-g5pts(x3)
           s_(1) = g5pts(z1)-g5pts(x1)
           s_(2) = g5pts(z2)-g5pts(x2)
           s_(3) = g5pts(z3)-g5pts(x3)
           weight = g5wts(x1)*g5wts(x2)*g5wts(x3)*g5wts(y1)*g5wts(y2)*g5wts(y3)*g5wts(z1)*g5wts(z2)*g5wts(z3)
           Sf = Sf + eval_composite(1,1,1,r_,s_,.true.)*weight
          end do
         end do
        end do
       end do
      end do
     end do
    end do
   end do
  end do
  Sf = Sf

  print *, sep, S, Sf, (S-Sf)

end program