#!/usr/bin/gnuplot
# plots the tx curves, with and without atom, for 7ns

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 output "pe.eps"

set border linewidth 1.5
set border 3

set multiplot
set nobar

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.05,1 spacing 1.1 samplen 0
set key

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
plot infile_1 u ($1):($2*100):($3*100)  with yerrorbars pt 7 ps 0.9 lc rgb "green" t "1.97(1) {/Symbol G}_0", infile_2 u ($1):($2*100):($3*100)  with yerrorbars pt 7 ps 0.9 lc rgb "orange" t "2.62(1) {/Symbol G}_0", infile_3 u ($1):($2*100):($3*100)  with yerrorbars pt 7 ps 0.9 lc rgb "blue" t "3.76(2) {/Symbol G}_0", infile_4 u ($1):($2*100):($3*100)  with yerrorbars pt 7 ps 0.9 lc rgb "red" t "4.76(3) {/Symbol G}_0",  infile_5 u ($1):($2*100):($3*100)  with yerrorbars pt 7 ps 0.9 lc rgb "black" t "6.08(5) {/Symbol G}_0", p_decay_1(x*1e-9+3e-9)*100 lw 3 lc rgb "green" t '', p_decay_2(x*1e-9+3e-9)*100 lw 3 lc rgb "orange" t '', p_decay_3(x*1e-9+3e-9)*100 lw 3 lc rgb "blue" t '', p_decay_4(x*1e-9+3e-9)*100 lw 3 lc rgb "red" t '', p_decay_5(x*1e-9+3e-9)*100  lw 3 lc rgb "black" 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