#!/usr/bin/gnuplot
# Plan: 1x2 plot, two plots show diff of tx curves


### input files ###
infile_tx_decaying='non_reversed_with_alternating_09_26_onwards_combined_histo_tx_2ns'
infile_tx_rising='reversed_with_alternating_09_26_onwards_combined_histo_tx_2ns'
infile_decaying_pulse='decaying_exp_pulse_fewerpoints'
infile_rising_pulse='rising_exp_pulse_fewerpoints'


### useful values ###
# the following values for the heralding efficiencies and accidental background are obtained from the 
# two files: non_reversed_with_alternating_09_26_onwards_combined_histo_tx_2ns and reversed_with_alternating_09_26_onwards_combined_histo_tx_2ns via python script: safwm_full_analysis.py.
eta_f = 0.00370 # accidental corrected, python script safwm_full_analysis.py
bg_decaying_atom = 0.000005147 # accidental events per time bin, python script safwm_full_analysis.py
bg_decaying_no_atom = 0.000005279 # accidental events per time bin, python script safwm_full_analysis.py
error_bg_decaying_atom = 0.0000000159
error_bg_decaying_no_atom = 0.0000000140
bg_rising_atom = 0.000005219 # accidental events per time bin, python script safwm_full_analysis.py
bg_rising_no_atom = 0.000005366 # accidental events per time bin, python script safwm_full_analysis.py
error_bg_rising_atom = 0.0000000164
error_bg_rising_no_atom = 0.0000000191

timebin=5e-9
timebin_trans=2e-9
offset_tx = 879
##########################
## theory :
# parameters:
gamma0 = 1/26.2e-9
gammap = 1/13.3e-9 
lambda = 0.033
#----
exp_rise(x)  =  x<0 ? exp(gammap/2*x) : 0
exp_decay(x)  =  x<0 ? 0: exp(-gammap/2*x) 

# solutions of diff eq for excited state population P_e(t)
p_decay(x) = x<0 ? 0 : 4*lambda*gamma0*gammap / (gammap-gamma0)**2 * ( exp(-gamma0/2*x) - exp(-gammap/2*x) )**2 
p_rise(x) =  x<0 ? 4*lambda*gamma0*gammap / (gammap+gamma0)**2 *exp(gammap*x) : 4*lambda*gamma0*gammap / (gammap+gamma0)**2 *exp(-gamma0*x)

# difference in transmission
diff_decay(x) = 2*sqrt( lambda*gamma0*gammap*p_decay(x))*exp_decay(x) - lambda*gamma0*p_decay(x)
diff_rise(x) = 2*sqrt( lambda*gamma0*gammap*p_rise(x))*exp_rise(x) - lambda*gamma0*p_rise(x)


###################

###################

set terminal postscript eps color enhanced solid "Helvetica,18" size 12.6cm, 8cm
set output "tx_diff.eps"

set border linewidth 1.5
set border 3
set samples 5000
set multiplot
set nobar
unset key

##### tx diff plots ####

set origin 0,0
set size 0.54,1
#set tmargin 1.5
set lmargin 9
set bmargin 4
#set rmargin 0
set xtics nomirror

set label 1 "a" at graph 0.1,0.96 font "Helvetica-Bold,22"


unset xlabel
set xlabel "Time from heralding event {/Times-Italic t_i} (ns)" offset graph 0.55,0
set ylabel "{/Symbol-Oblique d} ( x 10^{6} ) s^{-1}"
set xr [-40:200]
set xtics 40 nomirror 
set yr [-0.3:4]
set ytics 2 nomirror 

plot infile_tx_decaying u ($2-offset_tx ):(( ($9-bg_decaying_atom) - ($5-bg_decaying_no_atom) )/eta_f/timebin_trans*1e-6):(sqrt($10**2 + $6**2 + error_bg_decaying_atom**2 +error_bg_decaying_no_atom**2  )/eta_f/timebin_trans*1e-6) w yerrorbars pt 6 lc rgb "red" lw 2,\
(diff_decay(x*1e-9)*1e-6) w l lt 1 lc rgb "blue" lw 2


#######
set origin 0.45,0
set size 0.54,1
unset key
set border 3
unset xlabel
unset ylabel
set xr [-120:120]
set xtics 40 
# nomirror format " "
#set yr [-0.7:4]
set ytics 2     nomirror

set label 1 "b" at graph 0.11,0.96 font "Helvetica-Bold,22"


plot infile_tx_rising u ($2-offset_tx ):( ( ($9-bg_rising_atom) - ($5-bg_rising_no_atom) ) /eta_f/timebin_trans*1e-6):(sqrt($10**2+$6**2 + error_bg_rising_atom**2 + error_bg_rising_no_atom**2  )/eta_f/timebin_trans*1e-6) w yerrorbars pt 6 lc rgb "red" lw 2,\
(diff_rise(x*1e-9)*1e-6) w l lt 1 lc rgb "blue" lw 2


####### INSETS #######
set origin 0.35,0.75
set size 0.14,0.18
set lmargin 0
set bmargin 0
unset key
unset ylabel
unset xlabel
set xr [-0.3:1.5]
set yr [0:1]
unset xtics
#set xtics (0) offset -0.1,0.4 font "Helvetica,18"
unset border
unset ytics
set label 1 "{/Symbol-Oblique x}" at -0.8,0.5
set label 2 "t" at 1.6,0 font "Times-Italic,22"
set label 3 "0" at -0.15,-0.25 font "Helvetica,18"

#set style arrow 1 head filled size screen 0.005,40,40
#set arrow 1 from -0.85,0.4 to -0.85,0.1 as 1

set style arrow 1 lw 1.2 head filled size screen 0.01,15
set arrow 1 from -0.3,0 to 1.5,0 as 1 
set arrow 2 from -0.3,0 to -0.3,1.2 as 1

plot infile_decaying_pulse u ($1*1e7):($2*1e-4) w l lc rgb "blue" lw 2



set origin 0.85,0.75
set size 0.14,0.18
unset ylabel
unset xlabel
set xr [-1.3:0.3]
set yr [0:1]
unset xtics
unset ytics
#set xtics (0) offset -0.1,0.4 font "Helvetica,18"
set border 2 lw 1
unset ytics
set label 1 "{/Symbol-Oblique x}" at -1.8,0.5
set label 2 "t" at 0.65,0 font "Times-Italic,22"
set style arrow 1 lw 1.2 head filled size screen 0.01,15
set label 3 "0" at -0.15,-0.25 font "Helvetica,18"

set arrow 1 from -1.3,0 to 0.55,0 as 1 
set arrow 2 from -1.3,0 to -1.3,1.2 as 1
plot infile_rising_pulse u ($1*1e7):($2*1e-4) w l lc rgb "blue" lw 2

unset multiplot
unset output

# pause -1