#!/usr/bin/gnuplot
# Plan: 1x2 plot,  bottom plots show P_e curves

set border linewidth 1.5
set border 3
set samples 5000

### input files ###
infile_rx_decaying='non_reversed_with_alternating_09_26_onwards_combined_histo_rx_5ns'
infile_rx_rising='reversed_with_alternating_09_26_onwards_combined_histo_rx_5ns_4'
infile_decaying_pulse='decaying_exp_pulse_fewerpoints'
infile_rising_pulse='rising_exp_pulse_fewerpoints'
infile_ODE_rising='ODE_rising.dat'
infile_ODE_decaying='ODE_decaying.dat'

##########################
## theory :
# parameters:
gamma0 = 1/26.2e-9
gammap = 1/13.3e-9 
lambda = 0.033

### useful values ###
# values for efficiency and errors are calculated in python script: safwm_full_analysis.py
eta=3.7e-3/(0.46*0.519) #overall eff for two shapes, coupling of tx mode into fiber = 46% (fiber coupling + optical path losses), efficiency of tx APD = 51.9%
rx_det_eff=0.0078*0.56/0.6186 # 0.0078 was from saturation curve measurement done with APD Tauhuay [eff 61.86%]. Expt APD eff is 56%
timebin=5e-9
timebin_trans=2e-9
rx_scalefactor=gamma0*eta*rx_det_eff*timebin
rx_decaying_offset=0.000000310 #
rx_decaying_offset_err= 0.0000000048
rx_rising_offset= 0.000000325
rx_rising_offset_err= 0.0000000051
offset_tx = 879
#----
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)

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

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

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

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


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"


set xlabel "Time from heralding event {/Times-Italic t_i} (ns)" offset graph 0.55,0
#set ylabel "Coincidence \nprobability (x 10^{-7})" offset graph -0.05,0
set ylabel "{/Times-Italic P_{/Times-Roman e} (%)}" offset graph -0.05,0
set bmargin 3

set xr [-40:200]
set xtics 40 format "%g"
set ytics 1 format "%g" nomirror
set yr [-0.6:3]
plot infile_ODE_decaying u 1:($2*100):($3*100) w yerrorbars pt 6 lc rgb "red" lw 2,\
infile_rx_decaying u ($2 - (offset_tx+12)):(($5-rx_decaying_offset)/rx_scalefactor*100):(sqrt($6**2+rx_decaying_offset_err**2)/rx_scalefactor*100) w yerrorbars pt 13 ps 1.1 lc rgb "dark-green" lw 2,\
p_decay(x*1e-9)*100 w lines lt 1 lc rgb "blue" lw 2.5

##NOTE: the offset of rx curves wrt tx curves is set to 12ns here instead of 11ns from the g2 measurements, mainly because 1) there is an uncertainty of +/-1ns in the g2 measurement, and 2) the bin offsets had already been chosen, and it just looks a lot nicer if the exp rising photon has its theory curve peaked where the centre of a point is... I can of course redo it all and choose a different bin offset etc, but I think this is fair...

#######
set origin 0.45,0
set size 0.54,1
unset xlabel
unset ylabel

set xr [-120:120]
set xtics 40 format "%g"
set yr [-0.6:3]
set ytics 1

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


plot infile_ODE_rising u 1:($2*100):($3*100) w yerrorbars pt 6 lc rgb "red" lw 2,\
infile_rx_rising u ($2 - (offset_tx+12)):(($5-rx_rising_offset)/rx_scalefactor*100):(sqrt($6**2+rx_rising_offset_err**2)/rx_scalefactor*100) w yerrorbars pt 13 ps 1.1 lc rgb "dark-green" lw 2,\
p_rise(x*1e-9)*100 w lines lt 1 lc rgb "blue" lw 2.5


####### INSETS #######
set origin 0.35,0.75
set size 0.14,0.18
set lmargin 0
set bmargin 0
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