module evalTijk
! Last updated 9/2/2007
!
! This module includes 5 subroutines used to evaluate a
! model of the three-point third order velocity correlation 
! tensor for high Reynolds number isotropic turbulence
! as described in the Physics of Fluids paper entitled
!"An inertial range model for the three-point third-order velocity correlation"
! by Chang & Moser (2007). The model is expressed as the
! sum of a large number of terms which are rational in the magnitudes
! of the separation vectors. The coefficients and rational powers
! of these terms are read in from a data file, which is provided with
! these routines (basisdata.txt). See the document accompanying this 
! module (Tijk_define.pdf) for details of the data representation.
!
! This set of subroutines is intended to provide an example of the evaluation
! of the 3-point 3rd-order correlation model. There is certainly room for
! improvement in these routines. Two obvious improvements are:
! 1) When evaluating multiple components of the tensor at the same r_ and s_
!    it would be more efficient to evaluate them all at once.
! 2) It would be useful to implement an asymptotic form for the model that is
!    valid when one of the separation vectors is small, so that under such
!    conditions, evaluations are not polluted by roundoff error.
!
! The subroutines provided in this module (in the order in which 
! they would normally be called) are:
!
! read_basis(datafile)  : Reads the data describing the 5 basis tensors for
!                         the 3-point 3rd-order correlation. "datafile" is
!                         a string containing the name of the data file to read
!
! print_basis           : Prints the basis characterization data to stdout
!                         in the same format as the input data file, though 
!                         without comments. This is a debugging diagnostic.
!
! build_composite(weightsin,normalize)
!                       : Constructs a representation for a linear combination
!                         of the 5 basis tensors specified by weightsin
!                         an array of dimension 5 of double precision reals.
!                         This specifies a tensor in the 5-dimensional space
!                         of correlation tensors determined by the model.
!                         Normalize is a logical, if true, the weights will be
!                         normalized to sum to one, resulting in a correlation
!                         consistent with a dissipation of one. If not
!                         normalized, the dissipation related to the
!                         correlation is the sum of the coefficients in
!                         weightsin. If an individual basis is to be evaluated
!                         then weightsin should include all zeros, except for
!                         a one for the basis tensor to be evaluated.
!
! print_composite       : prints the characterization of the composite,
!                         similar to print_basis. This is intended as a
!                         debugging diagnostic
!
!
! eval_composite(i,j,k,r_,s_,checkprecision)
!                       : is a function returning the value of
!                         T_{ijk}(r_,s_)=<v_i(x_) v_j(x_+r_) v_k(x_+s_)>
!                         for the composite tensor defined by build_composite.
!                         i,j,k are integers, while r_ and s_ are 
!                         three-dimensional double precision real vectors.
!                         checkprecision is a logical. If true, eval_composite
!                         will return a NaN if the magnitude the
!                         smallest separation vector
!                         is more than a factor of 100 smaller than
!                         the largest separation vector. This is useful because
!                         for such large aspect ratios, the evaluation
!                         of the model is subject to pollution by round-off
!                         errors. When this aspect ratio is less than 100,
!                         the roundoff error will be less than a part in 10^5.

  implicit none
  integer :: basis(64,3,14,5) !stores triple (coeff,rpow,spow) for 
                              !up to 64 rational terms in r and s, for 14 
                              !tensor components of 5 basis functions
  integer :: nterms(14,5)     !stores the number of rational terms for 14 
                              !tensor components of 5 basis functions
  integer :: basisnorms(5)    !stores norms for 5 basis functions
  double precision :: compositecoeffs(64,14) !stores coeff for up to 64 
                              !rational terms in r and s, for 14 tensor 
                              !components of the composite function
  integer :: compositepowers(64,2,14) !stores (rpow,spow) for up to 64 
                              !rational terms in r and s, for 14 tensor 
                              !components of the composite function
  integer :: compositeterms(14) !stores the number of rational terms for 14 
                              !tensor components of the composite function
  double precision :: w_sum

contains

  subroutine read_comments
  !internal subroutine used by read_basis
  !ignores all lines that start with "!", carriage return lines, 
  !and lines with nothing but spaces
    implicit none
    character(500) :: str
    do
      str = ''
      read(unit=7,fmt='(a)') str
      if (scan(str,'!') == verify(str,' ') .or. len(trim(str)) == 0) then
!        write (*,*) "ignoring: "//str
      else
        backspace(unit=7)
        exit
      end if
    end do
  end subroutine read_comments
  
  subroutine read_basis(datafile)
  !reads in basis functions from datafile, such as "basisdata.txt"
  !comment lines (starting with "!") and carriage returns may be 
  !added to the datafile but order/format of the basis functions 
  !and order/format/spacing of tensor components must remain fixed
    implicit none
    integer :: i,j,k,ios,baseterms
    character(*), intent(in) :: datafile

    open(unit=7,file=datafile,iostat=ios,status="old",form="formatted",action="read")
    if (ios/=0) then
      print *,"unable to open file "
    else
      do k=1,5 !loop over basis functions
        call read_comments
        read (unit=7,fmt='(i8)') basisnorms(k)
        baseterms = 0
        do j=1,14 !loop over tensor components
          call read_comments
          read (unit=7,fmt='(i8)') nterms(j,k)
          baseterms = baseterms + nterms(j,k)
          call read_comments
          do i=1,nterms(j,k) !loop over rational terms
            read (unit=7,fmt='(i8)',advance="no") basis(i,1,j,k) !coeff
          end do
          read (unit=7,fmt=*)
          call read_comments
          do i=1,nterms(j,k) !loop over rational terms
            read (unit=7,fmt='(i8)',advance="no") basis(i,2,j,k) !rpow
          end do
          read (unit=7,fmt=*)
          call read_comments
          do i=1,nterms(j,k) !loop over rational terms
            read (unit=7,fmt='(i8)',advance="no") basis(i,3,j,k) !spow
          end do
          read (unit=7,fmt=*)
        end do
        write (*,*) "basis",k," total terms: ",baseterms
      end do
      close(unit=7)
    end if
  end subroutine read_basis


  subroutine print_basis
  !dumps to screen all basis data
  !this should be the same as datafile but uncommented
  !This is a debugging diagnostic
    implicit none
    integer :: i,j,k

    do k=1,5
      write (unit=*,fmt=*)
      write (unit=*,fmt=*)
      write (unit=*,fmt='(i8)') basisnorms(k)
      do j=1,14
        write (unit=*,fmt=*)
        write (unit=*,fmt='(i8)') nterms(j,k)
        do i=1,nterms(j,k)
          write (unit=*,fmt='(i8)',advance="no") basis(i,1,j,k)
        end do
        write (unit=*,fmt=*)
        do i=1,nterms(j,k)
          write (unit=*,fmt='(i8)',advance="no") basis(i,2,j,k)
        end do
        write (unit=*,fmt=*)
        do i=1,nterms(j,k)
          write (unit=*,fmt='(i8)',advance="no") basis(i,3,j,k)
        end do
        write (unit=*,fmt=*)
      end do
    end do
  end subroutine print_basis


  subroutine build_composite(weightsin,normalize)
  !This routine builds the description of the composite
  !specified by the linear combination of the bases with
  !coefficients in weightsin.
  !input: weights for composite function, where 
  !compositefunction = SUM(weights(n)*basisfunction(n))
  !output: compositecoeffs and compositepowers, coefficients 
  !and powers for rational functions in composite function

    implicit none
    double precision, intent(in) :: weightsin(5)
    logical, intent(in) :: normalize
    double precision :: weights(5)
    double precision :: composite(14,-7:10,-7:10) !intermediate array of 
                   !coefficents for 14 tensor components, 18 powers of r, 
                   !and 18 powers of s
    integer :: i,j,k,p,q

    weights = weightsin
    w_sum = weights(1)+weights(2)+weights(3)+weights(4)+weights(5)
    
    !if necessary, normalize weights so that epsilon=1
    if (w_sum /= 1. .and. w_sum /= 0. .and. normalize) then
      weights = weights/w_sum
      print *,"weights were normalized: ",weights
      w_sum = 1
    end if

    composite = 0.
    do k=1,5 !loop over basis functions
      if (weights(k) /= 0.) then
        do j=1,14 !loop over tensor components
          do i=1,nterms(j,k) !loop over rational terms
            !sum up coefficients in like terms (those with same tensor component, rpow, and spow) over 5 basis functions
            composite(j,basis(i,2,j,k),basis(i,3,j,k)) = composite(j,basis(i,2,j,k),basis(i,3,j,k)) + basis(i,1,j,k)*weights(k)/basisnorms(k)
          end do
        end do
      end if
    end do
    
    !build compositecoeffs and compositepowers from data in composite
    do j=1,14 !loop over tensor components
      i = 0
      do p=-7,10 !loop over rpow
        do q=-7,10 !loop over spow
          if (composite(j,p,q) /= 0) then
            i = i+1
            if (i > 64) then
              print *,"i > 64: compositecoeffs,compositepowers matrices too small"
            else
              compositecoeffs(i,j) = composite(j,p,q) !coeff
              compositepowers(i,1,j) = p !rpow
              compositepowers(i,2,j) = q !spow
            end if
          end if
        end do
      end do
      compositeterms(j) = i
    end do
   
  end subroutine build_composite


  subroutine print_composite
  !dump composite description (uncommented)
  !this if for diagnostics
    implicit none
    integer :: i,j

    do j=1,14
      write (unit=*,fmt=*)
      write (unit=*,fmt='(i8)') compositeterms(j)
      do i=1,compositeterms(j)
        write (unit=*,fmt='(i8)',advance="no") compositecoeffs(i,j)
      end do
      write (unit=*,fmt=*)
      do i=1,compositeterms(j)
        write (unit=*,fmt='(i8)',advance="no") compositepowers(i,1,j)
      end do
      write (unit=*,fmt=*)
      do i=1,compositeterms(j)
        write (unit=*,fmt='(i8)',advance="no") compositepowers(i,2,j)
      end do
      write (unit=*,fmt=*)
    end do
  end subroutine print_composite


  double precision function eval_composite(iin,jin,kin,rin_,sin_,checkprecision)
  !input: i,j,k components for three velcities
  !input: r_,s_ separation vectors
  !output: numerical value of composite function at (r_,s_): 
  !T_{ijk}(r_,s_)=<v_i(x_) v_j(x_+r_) v_k(x_+s_)>

    implicit none
    integer, intent(in) :: iin,jin,kin
    double precision, intent(in) :: rin_(3),sin_(3)
    logical, intent(in) :: checkprecision
    integer :: i,j,k !components of three vectors
    double precision :: r_(3),s_(3),t_(3) !separation vectors
    double precision :: r,s,t !lengths of separation vectors
    double precision :: rpow(-7:10), spow(-7:10) !r and s to 18 possible powers: -7,-6,...10
    integer :: DELTA_ij,DELTA_jk,DELTA_ik !kronecker delta evaluated for input i,j,k
    double precision :: sum1,sum2,minrst,maxrst,comp,tempreal
    integer :: m,n,p,tempint

    i=iin
    j=jin
    k=kin
    r_=rin_
    s_=sin_

    !calculate lengths r, s, t, and vector t_
    r = sqrt(r_(1)*r_(1) + r_(2)*r_(2) + r_(3)*r_(3))
    s = sqrt(s_(1)*s_(1) + s_(2)*s_(2) + s_(3)*s_(3))
    t_ = r_-s_
    t = sqrt(t_(1)*t_(1) + t_(2)*t_(2) + t_(3)*t_(3))

    if (i==j) then
      DELTA_ij = 1
    else
      DELTA_ij = 0
    end if
    if (j==k) then
      DELTA_jk = 1
    else
      DELTA_jk = 0
    end if
    if (i==k) then
      DELTA_ik = 1
    else
      DELTA_ik = 0
    end if

    if (r==0.) then !evaluate 2pt correlations (much simpler)
      eval_composite = w_sum/15*(s_(k)*DELTA_ij-1.5*(s_(i)*DELTA_jk+s_(j)*DELTA_ik))
    else if (s==0.) then
      eval_composite = w_sum/15*(r_(j)*DELTA_ik-1.5*(r_(k)*DELTA_ij+r_(i)*DELTA_jk))
    else if (t==0.) then
      eval_composite = -w_sum/15*(r_(i)*DELTA_jk-1.5*(r_(j)*DELTA_ik+r_(k)*DELTA_ij))
    else
     !evaluate full 3pt correlation

     minrst = min(r,s,t)
     maxrst = max(r,s,t)
     if (minrst/maxrst < 0.01 .and. checkprecision) then
       !print *,"minrst/maxrst < 0.01"
       eval_composite = 0./0
     else
      if (minrst/=r .and. minrst/=s) then
        !t is the smallest separation
        !we found that evaluation as t->0 is less accurate than as r->0 or s->0
        !therefore use symmetry T_{ijk}(r,s) = T_{jik}(-r,-t)
        !swap subscripts and variables
        tempint=i
        i=j
        j=tempint
        r_ = -r_
        s_ = -t_
        !note we don't set t_ again since it is not used below
        tempreal=s
        s=t
        t=tempreal
        if (i==j) then
          DELTA_ij = 1
        else
          DELTA_ij = 0
        end if
        if (j==k) then
          DELTA_jk = 1
        else
          DELTA_jk = 0
        end if
        if (i==k) then
          DELTA_ik = 1
        else
          DELTA_ik = 0
        end if
      end if

      !normalize variables by t (T_{ijk} will be un-normalized at the end)
      r_ = r_/t
      s_ = s_/t
      r = r/t
      s = s/t

      !pre-calculate r to different powers
      rpow(0) = 1.
      do p = 1,10
        rpow(p) = rpow(p-1)*r
      end do
      do p = 1,7
        rpow(-p) = 1./rpow(p)
      end do

      !pre-calculate s to different powers
      spow(0) = 1.
      do p = 1,10
        spow(p) = spow(p-1)*s
      end do
      do p = 1,7
        spow(-p) = 1./spow(p)
      end do


      sum2 = 0. !accumulates value of composite function

      do n=1,14 !loop over tensor components

        !determine which tensor component we are dealing with
        select case (n)
        case (1)
          comp = r_(i)*r_(j)*r_(k)
        case (2)
          comp = r_(i)*r_(j)*s_(k)
        case (3)
          comp = r_(i)*s_(j)*r_(k)
        case (4)
          comp = s_(i)*r_(j)*r_(k)
        case (5)
          comp = s_(i)*s_(j)*s_(k)
        case (6)
          comp = s_(i)*s_(j)*r_(k)
        case (7)
          comp = s_(i)*r_(j)*s_(k)
        case (8)
          comp = r_(i)*s_(j)*s_(k)
        case (9)
          comp = r_(i)*DELTA_jk
        case (10)
          comp = r_(j)*DELTA_ik
        case (11)
          comp = r_(k)*DELTA_ij
        case (12)
          comp = s_(i)*DELTA_jk
        case (13)
          comp = s_(j)*DELTA_ik
        case (14)
          comp = s_(k)*DELTA_ij
        end select
      
        if (comp /= 0.) then
          sum1 = 0. !accumulates value of rational function
          do m=1,compositeterms(n) !sum over rational terms
            sum1 = sum1 + compositecoeffs(m,n)*rpow(compositepowers(m,1,n))*spow(compositepowers(m,2,n))
          end do
          sum2 = sum2 + comp*sum1 !sum up tensor components times their respective rational functions
        end if
      end do
      !un-normalized by t
      eval_composite = sum2*t
     end if
    end if

  end function eval_composite

end module evalTijk