#!/usr/bin/gnuplot

## 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,8cm
set output "pe_01.eps"

set border linewidth 1.5
set border 3
set encoding utf8

set key at graph 1.01,1 spacing 1.1 samplen 0

set xr [-6:150]
set yr [-0.1:2]
set xtics 20 nomirror
set ytics 1 nomirror
set ytics add ("0" 0,"1" 1, "2" 2)

set xlabel "Time from heralding event {/Times-Italic t_i} (ns)"
set ylabel "Excited state population {/Times-Italic P_e}(\%)"  #(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

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.96(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.61(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.09(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 output "pe_02.eps"
# 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 add ("0" 0,"1" 1, "2" 2)
set yr [-0.1:2]

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

unset bars

plot pemax(tau_0/x)*100 lw 3 lc rgb "blue" notitle, \
    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  notitle