#!/usr/bin/gnuplot

set border linewidth 1.5
set encoding utf8


## input file ##
infile_1='ODE_13.dat'
infile_2='ODE_9.dat'
infile_3='ODE_6.dat'
infile_4='ODE_5.dat'
infile_5='ODE_4.dat'


##########################
## theory :
# parameters:
gamma0 = 1/26.2e-9
gammap_1 = 1/ 13.3e-9
gammap_2 = 1/ 10.0e-9
gammap_3 = 1/ 6.97e-9
gammap_4 = 1/ 5.50e-9
gammap_5 = 1/ 4.31e-9

lambda = 0.033
#----

# solutions of diff eq for excited state population P_e(t)
p_decay_1(x) = x<0 ? 0 : 4*lambda*gamma0*gammap_1 / (gammap_1-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap_1/2*x) )**2 
p_decay_2(x) = x<0 ? 0 : 4*lambda*gamma0*gammap_2 / (gammap_2-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap_2/2*x) )**2 
p_decay_3(x) = x<0 ? 0 : 4*lambda*gamma0*gammap_3 / (gammap_3-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap_3/2*x) )**2 
p_decay_4(x) = x<0 ? 0 : 4*lambda*gamma0*gammap_4 / (gammap_4-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap_4/2*x) )**2 
p_decay_5(x) = x<0 ? 0 : 4*lambda*gamma0*gammap_5 / (gammap_5-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap_5/2*x) )**2 
#------------------


set terminal postscript eps color enhanced solid font "Helvetica,22" size 12.6cm,16cm
# set terminal pdf color enhanced solid font "Helvetica,22" size 25.2cm,32cm

set output "pe.eps"

set border linewidth 1.5
set border 3

set multiplot

set origin 0.0,0.5
set size 1,0.45
set label 1 "a)" at graph -0.10,1.15 left
# set nobar
unset key

set key at graph 1.01,1 spacing 1.1 samplen 0


set xtics nomirror
set ytics nomirror
set ytics 1
set ytics add ("0" 0,"1" 1, "2" 2)
set xr [-6:150]
set yr [-0.1:2]
set xtics 20
set xlabel "Time from heralding event, {/Times-Italic t_i} (ns)"
set ylabel "excited state population (\%)"  #(x10^{-2})"
unset bars

# linestyles to programmatically change line styles
load "dark2.pal"
col_seq = (8, 6, 3, 2, 1)
set style fill transparent solid 0.3 noborder

th_offset = 3e-9#6e-9
pt_size = 1.1

# plot infile_1 u 1:(($2-$3)*100):(($2+$3)*100) with filledcurves  ls 8 notitle,\
#      infile_1 u ($1):($2*100) every 4 with points  ls 8 pt 7 ps pt_size t "1.97(1) {/Symbol G}_0",\
#      infile_2 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 6 notitle,\
#      infile_2 u ($1):($2*100) every 4 with points  ls 6 pt 7 ps pt_size t "2.62(1) {/Symbol G}_0",\
#      infile_3 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 3 notitle,\
#      infile_3 u ($1):($2*100) every 4 with points  ls 3 pt 7 ps pt_size t "3.76(2) {/Symbol G}_0",\
#      infile_4 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 2 notitle,\
#      infile_4 u ($1):($2*100) every 4 with points  ls 2 pt 7 ps pt_size t "4.76(3) {/Symbol G}_0",\
#      infile_5 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 1 notitle,\
#      infile_5 u ($1):($2*100) every 4 with points  ls 1 pt 7 ps pt_size t "6.08(5) {/Symbol G}_0" ,\
#      p_decay_1(x*1e-9+th_offset)*100 ls 8 lw 5 t '',\
#      p_decay_2(x*1e-9+th_offset)*100 ls 6 lw 5 t '',\
#      p_decay_3(x*1e-9+th_offset)*100 ls 3 lw 5 t '',\
#      p_decay_4(x*1e-9+th_offset)*100 ls 2 lw 5 t '',\
#      p_decay_5(x*1e-9+th_offset)*100 ls 1 lw 5 t ''

# plot infile_1 u 1:(($2-$3)*100):(($2+$3)*100) with filledcurves  ls 8 notitle,\
#      infile_2 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 6 notitle,\
#      infile_3 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 3 notitle,\
#      infile_4 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 2 notitle,\
#      infile_5 u 1:(($2-$3)*100):(($2+$3)*100)  with filledcurves ls 1 notitle,\
#      p_decay_1(x*1e-9+th_offset)*100 ls 8 lw 5 t '',\
#      p_decay_2(x*1e-9+th_offset)*100 ls 6 lw 5 t '',\
#      p_decay_3(x*1e-9+th_offset)*100 ls 3 lw 5 t '',\
#      p_decay_4(x*1e-9+th_offset)*100 ls 2 lw 5 t '',\
#      p_decay_5(x*1e-9+th_offset)*100 ls 1 lw 5 t ''

ev_step = 3

plot \
     infile_1 u 1:($2*100):($3*100) every ev_step with yerrorbars  ls 8  lw 3 pt 7 ps pt_size notitle,\
     infile_2 u 1:($2*100):($3*100) every ev_step with yerrorbars  ls 6  lw 3 pt 7 ps pt_size notitle,\
     infile_3 u 1:($2*100):($3*100) every ev_step with yerrorbars  ls 3  lw 3 pt 7 ps pt_size notitle,\
     infile_4 u 1:($2*100):($3*100) every ev_step with yerrorbars  ls 2  lw 3 pt 7 ps pt_size notitle,\
     infile_5 u 1:($2*100):($3*100) every ev_step with yerrorbars  ls 1  lw 3  pt 7 ps pt_size notitle,\
     infile_1 u 1:($2*100):($3*100) every ev_step with points  ls 8 pt 7 ps pt_size t "1.97(1) {/Symbol G}_0",\
     infile_2 u 1:($2*100):($3*100) every ev_step with points  ls 6 pt 7 ps pt_size t "2.62(1) {/Symbol G}_0",\
     infile_3 u 1:($2*100):($3*100) every ev_step with points  ls 3 pt 7 ps pt_size t "3.76(2) {/Symbol G}_0",\
     infile_4 u 1:($2*100):($3*100) every ev_step with points  ls 2 pt 7 ps pt_size t "4.76(3) {/Symbol G}_0",\
     infile_5 u 1:($2*100):($3*100) every ev_step with points  ls 1 pt 7 ps pt_size t "6.08(5) {/Symbol G}_0" ,\
     p_decay_1(x*1e-9+th_offset)*100 ls 8 lw 5 t '',\
     p_decay_2(x*1e-9+th_offset)*100 ls 6 lw 5 t '',\
     p_decay_3(x*1e-9+th_offset)*100 ls 3 lw 5 t '',\
     p_decay_4(x*1e-9+th_offset)*100 ls 2 lw 5 t '',\
     p_decay_5(x*1e-9+th_offset)*100 ls 1 lw 5 t ''

#########################
set origin 0.0,0.0
set size 1,0.45
set label 1 "b)" at graph -0.1,1.15 left
# theory:
lambda = 0.033
tau_0 = 26.23
pemax(x) = 4*lambda*(tau_0/x)**((tau_0+x)/(x-tau_0))
#-------------------------

infile_extinction_data_decaying='../tx_safwm_full_analysis/tx_safwm_summary.dat'

set ytics 1 nomirror format "%g"
set ytics 1
set ytics add ("0" 0,"1" 1, "2" 2)
set yr [-0.1:2]

set xtics autofreq nomirror
set xr [0:6.5]
#set xr [4:14]
set xlabel "Relative bandwidth {/Symbol G}_p / {/Symbol G}_0"
set ylabel "{/Times-Italic P_{e,max}} (\%)"  #(x10^{-2})"
set border 3
set key right top
unset grid
unset bars
#
set key invert
unset key

plot pemax(tau_0/x)*100 lw 3 lc rgb "blue" t "Eq.6 with {/Symbol L}=3.3%", \
    infile_extinction_data_decaying u (26.23/$1):(100*$5):(($2*26.23/($1*$1))):(100*$6)  w xyerrorbars pt 7 ps 2 lc rgb "red" lw 3  t "experiment"

unset multiplot
unset output

# pause -1