#!/usr/bin/gnuplot

set grid
set terminal pngcairo size 800,600

# analytical solution
mu2 = 1.8479e-05
mu1 = 8.904e-04
K = 0.105
Ly = 0.02
h = 0.5*Ly

set parametric
set trange [0:Ly]
set autoscale yfix
y(t) = t - h
Ux2(t) =  K/mu2*( 0.5*y(t)*y(t) - h*y(t) )
Ux1(t) = -K/mu1*( 0.5*y(t)*y(t) + h*y(t) )
Ux(t) = y(t) < 0 ? Ux1(t) : Ux2(t)

PATH = '../postProcessing/graph/'
LIST     = system('ls -1 '.PATH.' | sort -g')
NUM      = words(LIST)
TIMES(n) = word(LIST,n)
FILES(n) = PATH.TIMES(n).'/line.csv'
FILE = FILES(NUM)

REFR = './postProcessing/graph/25/line.csv'

#---Alpha---
set output 'alpha.png'
set title 't = '.TIMES(NUM).' s'
set xlabel '{/Symbol a}'
set yrange [0:Ly]
set ylabel 'y'
plot FILE u (column("alpha.water")):1 w lp t "water"           , \
     REFR u (column("alpha.water")):1 w lp t "reference"

#---Velocity---
set output 'velocity-x.png'
set title 't = '.TIMES(NUM).' s'
set xlabel 'U_x'
set yrange [0:Ly]
plot FILE u (column("U.water_x")):1 w lp t "water"             , \
     FILE u (column("U.air_x")  ):1 w lp t "air"               , \
     REFR u (column("U.water_x")):1 w lp t "reference water"   , \
     REFR u (column("U.air_x")  ):1 w lp t "reference air"     , \
     Ux(t), y(t)+h w l lw 2 lc 0 dt 2 title "analytical"

set output 'velocity-x_evolution.png'
set title 'evolution'
set cbrange [0:TIMES(NUM)]
set cblabel 't'
set xlabel 'U_x'
plot for [i=1:NUM] FILES(i) u (column("U.water_x")):1 w lp lc palette cb TIMES(i) notitle, \
     for [i=1:NUM] FILES(i) u (column("U.air_x")  ):1 w lp lc palette cb TIMES(i) notitle, \
     Ux(t), y(t)+h w l lw 2 lc 0 dt 2 title "analytical"

#---Drag-Force---
set output 'dragForce.png'
set title 't = '.TIMES(NUM).' s'
set xlabel 'drag_x'
set yrange [0:Ly]
plot FILE u (column("dragForce.water_x")):1 w lp t "water"     , \
     REFR u (column("dragForce.water_x")):1 w lp t "reference"

#---Dimless-Interface-Thickness---
set output 'yPlusI.png'
set xlabel 'y^{+I}'
set xrange [0:]
plot FILE u (column("yPlusI.air_water")):1 w lp t "interface"   , \
     REFR u (column("yPlusI.air_water")):1 w lp t "reference"
set autoscale x

print 'Finished'

