Fluid autocorrelation

module correlator
  implicit none

  integer, parameter :: rk = 8

  complex, allocatable :: vc(:)
  complex, allocatable :: sc(:)

  real, allocatable :: work(:)
  real, allocatable :: wsave(:)

  integer :: m, ier, lenwrk, lensav
contains
  function autocfft(v) ! autocorrelation via fft
    real(rk), allocatable :: autocfft(:)
    real(rk), allocatable :: v(:)

    m = size(v)

    call cfr(vc, v) ! create complex

    lenwrk = 2*m
    lensav = 2*m + int(log(real(m, kind = 4)) / log(2.0E+00)) + 4

    allocate(work(lenwrk))
    allocate(wsave(lensav))

    call cfft1i(m, wsave, lensav, ier) ! init fft

    ! fft of v
    call cfft1f(m, 1, vc, 1*(m-1) + 1, wsave, lensav, work, lenwrk, ier)

    sc = vc

    call conjugate(vc) ! complex conjugate of fft

    sc = sc * vc ! S = F x F* (Wiener-Khinchin)

    ! inverse fft
    call cfft1b(m, 1, sc, 1*(m-1) + 1, wsave, lensav, work, lenwrk, ier)

    call rfc(autocfft, sc) ! back to real

    deallocate(work)
    deallocate(wsave)
    deallocate(vc)
    deallocate(sc)
  end function autocfft

  function autoc(v) ! regular autocorrelation
    real(rk), allocatable :: autoc(:)
    real(rk), allocatable :: v(:)

    integer :: n, i

    n = size(v)

    allocate(autoc(n))

    do i = 1, n, 1
       autoc(i) = rs(v, real(i, kind = rk))
    end do
  end function autoc

  real(rk) function rs(v, dt) ! R(s)
    real(rk), allocatable :: v(:)
    real(rk) :: dt

    integer :: n, t
    real(rk) :: sd, s

    sd = stdv(v)
    n = size(v)
    s = 0
    do t = 1, n - dt, 1
       s = s + v(t)*v(t+dt)
    end do

    rs = s * (1.0_8/((real(n) - dt)*sd**2))
  end function rs

  function eulertime(v) ! time scale
    real(rk) :: eulertime
    real(rk), allocatable :: v(:)
    real(rk), allocatable :: r(:)

    integer :: i, n

    n = size(v)

    r = autocfft(v)

    eulertime = 0
    do i = 1, n - 2, 1
       if (isnan(r(i)) .eqv. .false.) then
          eulertime = eulertime + r(i)
       endif
    end do
  end function eulertime

  function eulerlength(v) ! length scale
    real(rk) :: eulerlength
    real(rk), allocatable :: v(:)

    eulerlength = eulertime(v) * mean(v)
  end function eulerlength

  subroutine mka(a, min, max, diff) ! make array
    real(rk), allocatable :: a(:)
    real(rk) :: max, min, diff

    integer :: n, i

    n = (max - min)/diff

    allocate(a(n))

    a(1) = min
    do i = 2, n, 1
       a(i) = a(i - 1) + diff
    end do
  end subroutine mka

  subroutine cfr(c, r) ! complex from real array
    real(rk), allocatable :: r(:)
    complex, allocatable :: c(:)

    integer :: i, n

    n = size(r)

    allocate(c(n))

    do i = 1, n, 1
       c(i) = complex(r(i), 0)
    end do
  end subroutine cfr

  subroutine rfc(r, c) ! real from complex array
    real(rk), allocatable :: r(:)
    complex, allocatable :: c(:)

    integer :: i, n

    n = size(c)

    allocate(r(n))

    do i = 1, n, 1
       r(i) = real(c(i))
    end do
  end subroutine rfc

  subroutine conjugate(c) ! conjugate of complex array
    complex, allocatable :: c(:)

    integer :: i, n

    n = size(c)

    do i = 1, n, 1
       c(i) = complex(real(c(i)), -aimag(c(i)))
    end do
  end subroutine conjugate

  real function mean(a) ! arithmetic mean of array
    real(rk), allocatable :: a(:)

    integer :: i, n
    real(rk) :: s

    n = size(a)
    s = 0.0
    do i = 1, n, 1
       s = s + a(i)
    end do

    mean = s / real(n, kind = rk)
  end function mean

  real function stdv(a) ! standard deviation of array
    real(rk), allocatable :: a(:)

    integer :: i, n, m
    real(rk) :: s

    m = mean(a)

    n = size(a)
    s = 0.0
    do i = 1, n, 1
       s = s + (a(i) - m)**2
    end do

    s = s / real(n, kind = rk)

    stdv = sqrt(s)
  end function stdv

  function parser(p) ! create array from file
    character(32) :: p
    real(rk), allocatable :: parser(:)

    integer :: i, n, io
    real(rk), allocatable :: a(:)
    real(rk) :: s

    open(8, file = p, action = 'read')

    allocate(a(0))

    n = 2
    i = 1
    do while (i < n)
       read(8, *, iostat = io) s
       if (io == 0) then
          n = n + 1
          call append(a, s)
       end if
       i = i + 1
    end do

    allocate(parser(n))

    close(8)

    parser = a
  end function parser

  real(rk) function c2r(c) ! complex to real
    character(32) :: c
    read(c, *) c2r
  end function c2r

  subroutine append(a, v) ! append real to array
    real(rk), allocatable :: a(:)
    real(rk) :: v

    real(rk), allocatable :: b(:)
    integer :: n

    n = size(a)

    allocate(b(n))

    b = a
    deallocate(a)
    allocate(a(n + 1))

    a(1:n) = b
    a(n + 1) = v
  end subroutine append
end module correlator

program main
  use correlator
  use gnufor2 ! gnuplot

  implicit none

  character(32) :: s
  integer :: n

  real(rk), allocatable :: v(:) ! velocity (m/s)
  real(rk), allocatable :: t(:) ! time (s)
  real(rk), allocatable :: r(:) ! result

  call get_command_argument(1, s)

  print *, "Reading data..."

  v = parser(s) ! read data from file in command line arg

  n = size(v)

  call mka(t, 0.0_8, real(n, kind = 8), 1.0_8)

  print *, "Calculating autocorrelation function..."

  r = autocfft(v)

  print *, "Eulerian integral time scale: ", eulertime(v)
  print *, "Eulerian integral length scale: ", eulerlength(v)

  call plot(t, v) ! plot velocity
  call plot(t, r) ! plot autocorrelation
end program main