#!/bin/bash

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

gnuplot<<EOF

    set terminal pngcairo enhanced color
    set output '../numberDensity.png'
    set decimalsign '.'

    set grid
    set key box
    set key font ",11"

    set format xy '%g'
    set xrange [0:4]
    set xlabel 'd [mm]'
    set ylabel 'n [mm^-{1} m^{-3}]'

    set key top left

    plot '../postProcessing/numberDensity/0/numberDensity.xy' u (\$1*1000):2 w l ls 1 dt 2 lc rgb 'black'  t 'Initial condition',\
    '../postProcessing/numberDensity/10/numberDensity.xy' u (\$1*1000):2 w l lc rgb 'black' t 'Result'

EOF

gnuplot<<EOF

    set term pngcairo size 550,500
    set terminal pngcairo enhanced color
    set output '../CO2.png'

    set grid
    set key box vertical width 0.5 height 1 maxcols 1 spacing 1.5
    set key font ",11"

    set xrange [0:10]
    set xlabel 't [s]'
    set ylabel 'Mass per unit volume'

    set key right center

    plot  '../postProcessing/probes/0/multiply(alpha.gas3,CO2.gas3,rho.gas3)' w l t '{CO_2}.gas3',\
    '../postProcessing/probes/0/multiply(alpha.liquid,CO2.liquid,rho.liquid)' w l t '{CO_2}.liquid',\
    '../postProcessing/probes/0/add(multiply(alpha.gas3,CO2.gas3,rho.gas3),multiply(alpha.liquid,CO2.liquid,rho.liquid))' w l t '{CO_2} total mass'

EOF

gnuplot<<EOF

    set term pngcairo size 550,500
    set terminal pngcairo enhanced color

    set grid
    set key box
    set key font ",11"

    set xrange [0:10]
    set xlabel 't [s]'

## Equations for analytical solution
# Initial values

    Y_G0 = 0.9        #initial gas CO2 mass fraction
    a_G0 = 1e-3       #initial gas fraction
    Y_L0 = 0          #initial liquid CO2 mass fraction
    d_B0 = 3.0        #initial bubble diameter in mm
    S = 1.0           #time scale in sec (T^-1)
    R = 1.0 / 1000    #ratio of gas and liquid densities

# Function for Y_G(t)

    Y_G(x) = 1 - ((1 - Y_G0) / (Y_G0 * exp(-S * x) + (1 - Y_G0)))

# Function for alpha_G(t)

    a_G1(x) = a_G0 * (1 - Y_G0 * (1 - exp(-S * x)))
    a_G(x) = a_G1(x) * 1e3   #Scale a_G by 1e3

# Function for Y_L(t)

    Y_L1(x) = (1 - (1 - Y_L0) * ((1 - a_G0) / ((1 - a_G0 + a_G0 * Y_G0) - a_G0 * Y_G0 * exp(-S * x)))**R)
    Y_L(x) = Y_L1(x) * 1e6   #Scale Y_L by 1e6

# Function for d_B(t)

    d_B(x) = ((a_G1(x) * (d_B0**(3.0))) / a_G0)**(1.0 / 3.0)

    set output '../alphaG.png'

    set ylabel '{/Symbol a}_G [-] {/Symbol \264} 10^{-3}'
    set yrange [0:]

    plot a_G(x) w l dt "-.." lw 4 lc rgb 'black' t 'Analy.',\
    '../postProcessing/probes/0/alpha.gas3' u 1:(\$2/1e-3) w l lw 2 lc rgb 'red' t 'Sim.'

    set output '../db.png'

    set ylabel 'd_B [mm]'
    set yrange [1:4]
    set format y "%g"
    set ytics 1, 0.6, 4

    plot d_B(x) w l dt "-.." lw 4 lc rgb 'black' t 'Analy.',\
    '../postProcessing/probes/0/d.gas3' u 1:(\$2*1000) w l lw 2 lc rgb 'blue' t 'Sim.'

    set output '../Y.png'

    set ylabel '{Y_G}^{A} [-]'
    set y2tics  nomirror
    set y2label '{Y_L}^{A} [-] {/Symbol \264} 10^{-6}'
    set yrange [0:1]
    set y2range [0:1]

    set key center right

    plot Y_G(x) w l dt "-.." lw 4 lc rgb 'black' t 'Analy.',\
    '../postProcessing/probes/0/CO2.gas3' u 1:2 w l lw 2 lc rgb 'red' t 'Sim.',\
    Y_L(x) w l dt "-.." lw 4 lc rgb 'black' t 'Analy.' axis x1y2,\
    '../postProcessing/probes/0/CO2.liquid' u 1:(\$2*1000000) w l lw 2 lc rgb 'blue' t 'Sim.' axis x1y2

EOF

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