#!/usr/bin/gnuplot
# Plan: 2x2 plot, top two plots show diff of tx curves, bottom plots show rx curves

set border linewidth 1.5
set border 3
set samples 5000

### 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_tx_decaying_diff_theory='tx_diff_model_binning_decaying_2ns'
infile_tx_rising_diff_theory='tx_diff_model_binning_rising_2ns'
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_rx_decaying_theory='decaying_exp_reflection_model_binning_5ns'
infile_rx_rising_theory='rising_exp_reflection_model_binning_5ns'
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'

### useful values ###
decaying_data_integral=0.00368197041753 #obtained from ipython notebook
rising_data_integral=0.00370053312007
gamma0=38124285
eta=3.68e-3*2*2
rx_det_eff=0.0078
timebin=5e-9
timebin_trans=2e-9
rx_scalefactor=gamma0*eta*rx_det_eff*timebin
rx_decaying_offset=3.11e-7
rx_decaying_offset_err=4.67e-9
rx_rising_offset=3.2e-7
rx_rising_offset_err=4.95e-9
##########################
## theory :
# parameters:
gamma0 = 1/26.2e-9
gammap = 1/13.5e-9 
lambda = 0.03
#----
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)


    #integrated versions
p_decay_integral(x) = x<0 ? 0 : 4*lambda*gamma0*gammap / (gammap-gamma0)**2 * ( (1-exp(-gamma0*x))/gamma0 + (1-exp(-gammap*x))/gammap + 4*( exp(-(gamma0+gammap)*x/2) -1 )/(gamma0+gammap) )

p_rise_integral(x) = x<0 ? 4*lambda*gamma0*exp(gammap*x)/(gamma0+gammap)**2 : 4*lambda/(gamma0+gammap)**2 * (gamma0 + gammap - gammap*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)

    #integrated versions, gonna be real ugly...
diff_decay_integral(x) = x<0 ? 0: -4*gamma0**2*lambda**2*gammap/(gamma0-gammap)**2 * ( (1-exp(-gamma0*x))/gamma0 + (1-exp(-gammap*x))/gammap + 4*(exp(-(gamma0+gammap)*x/2) -1)/(gamma0+gammap) ) + 4*gammap*gamma0*lambda/sqrt((gamma0-gammap)**2) * ( (exp(-gammap*x) -1)/gammap + 2*(1-exp(-(gamma0+gammap)*x/2))/(gamma0+gammap) )

diff_rise_integral(x) = x<0 ? 4*gamma0*lambda*(gamma0+gammap-(gamma0*lambda))/(gamma0+gammap)**2 * exp(gammap*x) : 4*gamma0*lambda/(gamma0+gammap)**2 * exp(-gamma0*x) * ( -exp(gamma0*x)*(gamma0+gammap)*(lambda-1) + lambda*gammap) 

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

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

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

set multiplot
set nobar
unset key

##### tx diff plots ####
set origin 0,0.48
set size 0.56,0.52
set lmargin 9

unset xlabel
# set ylabel "P(scattered) (x 10^{6})"
set ylabel "{/Symbol-Oblique d} (x 10^{6})"
set xr [-40:200]
set xtics 40 nomirror 
# format " "
set yr [-0.3:4]
set ytics 2 nomirror 

plot infile_tx_decaying u ($2-881.5):(($9-$5)/decaying_data_integral/timebin_trans*1e-6):(sqrt($10**2+$6**2)/decaying_data_integral/timebin_trans*1e-6) w yerrorbars pt 6 lc rgb "red" lw 2,\
(diff_decay_integral(x*1e-9 + timebin_trans/2)-diff_decay_integral(x*1e-9 - timebin_trans/2))/timebin_trans*1e-6 w l lt 1 lc rgb "blue" lw 2.5
#(diff_decay(x*1e-9)*1e-6) w l lt 1 lc rgb "blue" lw 2.5
#infile_tx_decaying_diff_theory u 1:($2*1000) w l lt 1 lc rgb "blue" lw 2.5

#######
set origin 0.45,0.48
set size 0.56,0.52


unset xlabel
unset ylabel
set xr [-120:120]
set xtics 40 
# nomirror format " "
#set yr [-0.7:4]
set ytics 2     nomirror

plot infile_tx_rising u ($2-878.5):(($9-$5)/rising_data_integral/timebin_trans*1e-6):(sqrt($10**2+$6**2)/rising_data_integral/timebin_trans*1e-6) w yerrorbars pt 6 lc rgb "red" lw 2,\
(diff_rise_integral(x*1e-9 + timebin_trans/2)-diff_rise_integral(x*1e-9 - timebin_trans/2))/timebin_trans*1e-6 w l lt 1 lc rgb "blue" lw 2.5
#(diff_rise(x*1e-9)*1e-6) w l lt 1 lc rgb "blue" lw 2.5
#infile_tx_rising_diff_theory u 1:($2*1000) w l lt 1 lc rgb "blue" lw 2.5


#### rx plots ####
set origin 0,0
set size 0.56,0.52
set xlabel "Time from heralding event, {/Times-Italic t_i} (ns)" offset graph 0.6,-0.05
#set ylabel "Coincidence \nprobability (x 10^{-7})" offset graph -0.05,0
set ylabel "{/Times-Italic P_e} (%)" offset graph -0.05,0
set bmargin 3

set xr [-40:200]
set xtics 40 format "%g"
set ytics 1 format "%g" 
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 - (881.5+12)):(($5-rx_decaying_offset)/rx_scalefactor*100):(($6+rx_decaying_offset_err)/rx_scalefactor*100) w yerrorbars pt 12 lc rgb "dark-green" lw 2,\
p_decay(x*1e-9)*100 w lines lt 1 lc rgb "blue" lw 2.5
#(p_decay_integral(x*1e-9 + timebin/2)-p_decay_integral(x*1e-9 - timebin/2))/timebin*100 w lines lt 1 lc rgb "blue" lw 2.5
#infile_rx_decaying_theory u ($1):((($2*1e-7)-rx_decaying_offset)/rx_scalefactor*100) w lines lt 1 lc rgb "web-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.56,0.52
unset xlabel
unset ylabel

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

plot infile_ODE_rising u 1:($2*100):($3*100) w yerrorbars pt 6 lc rgb "red" lw 2,\
infile_rx_rising u ($2 - (878.5+12)):(($5-rx_rising_offset)/rx_scalefactor*100):(($6+rx_rising_offset_err)/rx_scalefactor*100) w yerrorbars pt 12 lc rgb "dark-green" lw 2,\
p_rise(x*1e-9)*100 w lines lt 1 lc rgb "blue" lw 2.5
#(p_rise_integral(x*1e-9 + timebin/2)-p_rise_integral(x*1e-9 - timebin/2))/timebin*100 w lines lt 1 lc rgb "blue" lw 2.5
#infile_rx_rising_theory u ($1):((($2*1e-7)-rx_rising_offset)/rx_scalefactor*100) w lines lt 1 lc rgb "web-blue" lw 2.5


####### INSETS #######
set origin 0.25,0.7
set size 0.28,0.26
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.69,0.7
set size 0.28,0.26
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