#!/bin/bash

if ! which gnuplot > /dev/null 2>&1
then
    echo 'gnuplot not found - skipping graph creation' >&2
    exit 1
fi

gnuplot<<EOF


TIME    = system('foamListTimes -latestTime -case ..')
FILE    = "../postProcessing/graph/".TIME."/line.csv"

EFILE(p) = 'experiment/F250/'.p.'.dat'

# Reported mean water level by Fabre (1987)
ymeas = real(system("foamDictionary ../caseParameterDict -entry hWater -value"))

# Simulated interface region (0.01 < alpha < 0.99) (green area on plots)
fsmin = real(system(sprintf("tail -n+2 %s | awk '\$13>0.01{print \$1; exit}'", FILE)))
fsmax = real(system(sprintf("tail -n+2 %s | awk '\$13>0.99{print \$1; exit}'", FILE)))

set linetype 1 lc rgb "black" lw 1.5 dt 1
set linetype 2 lc rgb "red"   lw 1   dt 1
set linetype 3 lc rgb "blue"  lw 2.5 dt 2
set linetype 4 lc rgb "green" lw 1.5 dt 2

set style fill solid 0.1 noborder

set term pngcairo enhanced
set grid


#-------------------------Streamwise velocity profile--------------------------
set output 'U_x.png'

set key top left
set ylabel "y [m]"
set xlabel "U_x [m/s]"

set x2label "{/Symbol a} [-]"
set x2range [0:1]
set x2tics 0, 0.1, 1

plot fsmax                                                       axis x2y1 w filledcurves y=fsmin ls 4                  notitle                 , \
     EFILE('gasVel')   u (column("UMean.air_x")):(column("y"))             w p                    ls 2  pt 6            t "Fabre (1987): air"   , \
     EFILE('liqVel_x') u (column("UMean.water_x")):(column("y"))           w p                    ls 2  pt 7            t "water"               , \
     ymeas                                                       axis x2y1 w l                    ls 2  lw 1.5 dt 4     t "water level"         , \
     1/0                                                                                          ls -1 lc rgb("white") t " "                   , \
     FILE              u (column("U.airC_x")):(column("y"))                w l                    ls 1                  t "Simulated: airC"     , \
     FILE              u (column("U.water_x")):(column("y"))               w l                    ls 3                  t "water"               , \
     1/0                                                         axis x2y1 w filledcurves y=0     ls 4                  t "interface region"    , \
     FILE              u (column("alpha.airC")):(column("y"))    axis x2y1 w l                    ls 4  dt 4            t "{/Symbol a}_{airC}"


#--------------------------Vertical velocity profile---------------------------
set output 'U_y.png'

set key top right
set xlabel "U_y [m/s]"
# set xtics  1e-3

plot fsmax                                                       axis x2y1 w filledcurves y=fsmin ls 4                  notitle                 , \
     EFILE('gasVel')   u (column("UMean.air_y")  ):(column("y"))           w p                    ls 2  pt 6            t "Fabre (1987): air"   , \
     ymeas                                                       axis x2y1 w l                    ls 2  lw 1.5 dt 4     t "water level"         , \
     FILE              u (column("U.airC_y")     ):(column("y"))           w l                    ls 1                  t "Simulated: airC"     , \
     FILE              u (column("U.water_y")    ):(column("y"))           w l                    ls 3                  t "water"               , \
     1/0                                                         axis x2y1 w filledcurves y=0     ls 4                  t "interface region"    , \
     FILE              u (column("alpha.airC")   ):(column("y")) axis x2y1 w l                    ls 4  dt 4            t "{/Symbol a}_{airC}"


#----------------------Turbulence kinetic energy profile-----------------------
set output 'k.png'

set key top left
set xlabel "k [m^2/s^2]"
set logscale x
set format x '10^{%T}'
set xrange[1e-7:10]

set xtics 10
set xtics nomirror
set x2range [0:1]
set x2tics 0,0.1,1 nomirror

kmix = "(column('alpha.airC')*column('k.airC') + column('alpha.water')*column('k.water'))"

plot fsmax                                                  axis x2y1  w filledcurves y=fsmin ls 4                  notitle                 , \
     EFILE('gasKe') u (column("kMean.air")  ):(column("y"))            w p                    ls 2  pt 6            t "Fabre (1987): air"   , \
     EFILE('liqKe') u (column("kMean.water")):(column("y"))            w p                    ls 2  pt 7            t "water"               , \
     ymeas                                                  axis x2y1  w l                    ls 2  lw 1.5 dt 4     t "water level"         , \
     1/0                                                                                      ls -1 lc rgb("white") t " "                   , \
     FILE           u @kmix                  :(column("y"))            w l                    ls 1                  t "Simulated: mixture"  , \
     1/0                                                    axis x2y1  w filledcurves y=0     ls 4                  t "interface region"    , \
     FILE           u (column("alpha.airC") ):(column("y")) axis x2y1  w l                    ls 4  dt 4            t "{/Symbol a}_{airC}"

EOF

#------------------------------------------------------------------------------
