#!/usr/bin/gnuplot
infile1 = 'histo.dat'
infile2 = '2017-04-20/bi_tx_2017_04_20'
infile3 = '2017-04-20/blue_tx_2017_04_20'
infile4 = '2017-04-20/red_tx_2017_04_20'
infile5 = 'histo_bi.dat'

blue_defined = '#0078FF'

set terminal postscript eps color enhanced solid font "Helvetica,20" size 20cm, 9cm
set output "hist_spec.eps"

###### FITTING ########################################
lorentz_b(x) = 1 - tx_b / ( 4*( (x0_b-x) / (gamma_b) )**2 + 1 )
tx_b = 0.35
lorentz_r(x) = 1 - tx_r / ( 4*( (x0_r-x) / (gamma_r) )**2 + 1 )
tx_r = 0.35
lorentz_sum(x) = 1 - tx_sum / ( 4*( (x0_sum-x) / (gamma_sum) )**2 + 1 )
tx_sum = 0.35

lorentz_bs(x) = 1 - tx_bs / ( 4*( (x0_bs-x) / (gamma_bs) )**2 + 1 )
tx_bs = 0.15
lorentz_rs(x) = 1 - tx_rs / ( 4*( (x0_rs-x) / (gamma_rs) )**2 + 1 )
tx_rs = 0.15

lorentz_dis_b(x) = (1 + (4*lambda_b**2 - 4*lambda_b*cos(phi_b))*gamma_b**2/ ( 4*(x0_b-x)**2 + (gamma_b)**2 ) - ( 8*lambda_b*sin(phi_b)*(x-x0_b)*gamma_b ) / ( 4*(x0_b-x)**2 + gamma_b**2 ) ) 
gamma_b = 6.9
x0_b = 40
phi_b = 0.01
lambda_b = 0.047

lorentz_dis_r(x) = (1 + (4*lambda_r**2 - 4*lambda_r*cos(phi_r))*gamma_r**2/ ( 4*(x0_r-x)**2 + (gamma_r)**2 ) - ( 8*lambda_r*sin(phi_r)*(x-x0_r)*gamma_r ) / ( 4*(x0_r-x)**2 + gamma_r**2 ) ) 
gamma_r = 6.9
x0_r = 40
phi_r = 0.01
lambda_r = 0.047

lorentz_dis_bs(x) = (1 + (4*lambda_bs**2 - 4*lambda_bs*cos(phi_bs))*gamma_bs**2/ ( 4*(x0_bs-x)**2 + (gamma_bs)**2 ) - ( 8*lambda_bs*sin(phi_bs)*(x-x0_bs)*gamma_bs ) / ( 4*(x0_bs-x)**2 + gamma_bs**2 ) ) 
gamma_bs = 6.9
x0_bs = 40
phi_bs = 0.01
lambda_bs = 0.047

lorentz_dis_rs(x) = (1 + (4*lambda_rs**2 - 4*lambda_rs*cos(phi_rs))*gamma_rs**2/ ( 4*(x0_rs-x)**2 + (gamma_rs)**2 ) - ( 8*lambda_rs*sin(phi_rs)*(x-x0_rs)*gamma_rs ) / ( 4*(x0_rs-x)**2 + gamma_rs**2 ) ) 
gamma_rs = 6.9
x0_rs = 40
phi_rs = 0.01
lambda_rs = 0.047

lorentz_dis_sum(x) = (1 + (4*lambda_sum**2 - 4*lambda_sum*cos(phi_sum))*gamma_sum**2/ ( 4*(x0_sum-x)**2 + (gamma_sum)**2 ) - ( 8*lambda_sum*sin(phi_sum)*(x-x0_sum)*gamma_sum ) / ( 4*(x0_sum-x)**2 + gamma_sum**2 ) ) 
gamma_sum = 6.9
x0_sum = 40
phi_sum = 0.01
lambda_sum = 0.047

#fit lorentz_dis_sum(x) infile2 u (($1-190000)/1000*2):11:12 via x0_sum, gamma_sum, lambda_sum, phi_sum
#fit lorentz_dis_b(x) infile2 u (($1-190000)/1000*2):7:8 via x0_b, gamma_b, lambda_b, phi_b
#fit lorentz_dis_r(x) infile2 u (($1-190000)/1000*2):9:10 via x0_r, gamma_r, lambda_r, phi_r
#fit lorentz_dis_rs(x) infile4 u (($1-190000)/1000*2):4:5 via x0_rs, gamma_rs, lambda_rs, phi_rs
#fit lorentz_dis_bs(x) infile3 u (($1-190000)/1000*2):2:3 via x0_bs, gamma_bs, lambda_bs, phi_bs

fit lorentz_sum(x) infile2 u (($1-190000)/1000*2):11:12 via x0_sum, gamma_sum, tx_sum
fit lorentz_b(x) infile2 u (($1-190000)/1000*2):7:8 via x0_b, gamma_b, tx_b
fit lorentz_r(x) infile2 u (($1-190000)/1000*2):9:10 via x0_r, gamma_r, tx_r
fit lorentz_rs(x) infile4 u (($1-190000)/1000*2):4:5 via x0_rs, gamma_rs, tx_rs
fit lorentz_bs(x) infile3 u (($1-190000)/1000*2):2:3 via x0_bs, gamma_bs, tx_bs

#################################
#######
####GENERAL SETTINGS####
set border linewidth 1.5
set samples 5000
set multiplot 
set pointsize 1.5
unset bar
unset key
set style fill solid 0.3 noborder

###########################
#####     ONESIDE SPEC       #####
set xrange [23:58]
set yrange [58:102]
set ytics 10 nomirror 
set xtics 10 nomirror
set xlabel 'Probe detuning (MHz)'
set ylabel 'Transmission (%)' offset 1.5,0
set title 'One-sided illumination'

set lmargin at screen 0.08
set rmargin at screen 0.48

set label 1 sprintf("(a)") at graph -0.1,1.1 font "Helvetica, 24"

plot lorentz_rs(x)*100 lc rgb 'red' lw 2 not, \
	lorentz_bs(x)*100 lc rgb blue_defined lw 2 not, \
	infile4 u (($1-190000)/1000*2):($4*100):($5*100) w yerr ps 1.3 pt 5 lw 3 lc rgb 'red' not, \
	infile3 u (($1-190000)/1000*2):($2*100):($3*100) w yerr ps 1.8 pt 13 lw 3 lc rgb blue_defined not
	
#BI-SPEC

set title '4Pi illumination'
set lmargin at screen 0.58
set rmargin at screen 0.98

set label 1 sprintf("(b)")

plot lorentz_r(x)*100 lc rgb 'red' lw 2 not, \
lorentz_b(x)*100 lc rgb blue_defined lw 2 not, \
lorentz_sum(x)*100 lt 1 lc rgb 'black' lw 2 not, \
infile2 u (($1-190000)/1000*2):($9*100):($10*100) w yerr ps 1.3 pt 5 lc rgb 'red' lw 3 not, \
infile2 u (($1-190000)/1000*2):($7*100):($8*100) w yerr ps 1.8 pt 13 lc rgb blue_defined lw 3 not, \
infile2 u (($1-190000)/1000*2):($11*100):($12*100) w yerr pt 7 lc rgb 'black' lw 3 not


set border 15
###########################
#####      HISTO      #####
set lmargin at screen 0.35
set rmargin at screen 0.45
set xtics nomirror
set ytics nomirror
unset ylabel
unset title
unset label

#RED HISTO
set xrange [0:100]
set yrange [0:0.18]
#set ytics 0.1 offset 0.5
unset ytics
set xtics 50 nomirror out offset 0,0.2
set ylabel 'Occurences' offset 0,1
set xlabel 'Photons' offset 0,0.7

set tmargin at screen 0.35
set bmargin at screen 0.25
plot infile1  u 1:($5) lc -1 with boxes , '' u 1:($4) lc rgb 'red' lw 3 with steps

#BLUE HISTO
unset xlabel
unset xtics
unset ylabel
set tmargin at screen 0.45
set bmargin at screen 0.35
plot infile1  u 1:($3) lc -1 with boxes , '' u 1:($2) lc rgb blue_defined lw 3 with steps

#BI HISTO
set lmargin at screen 0.85
set rmargin at screen 0.95
set xrange [0:200]
set yrange [0:0.34]
set xtics 100 nomirror out offset 0,0.2
set xlabel 'Photons' offset 0,0.7
set ylabel 'Occurences' offset 0,2

set tmargin at screen 0.35
set bmargin at screen 0.25
plot infile5  u 1:($7) lc -1 with boxes , '' u 1:($6) lc rgb 'black' lw 3 with steps

#red
unset xlabel
unset xtics
unset ylabel
set tmargin at screen 0.45
set bmargin at screen 0.35
plot infile5  u 1:($5) lc -1 with boxes , '' u 1:($4) lc rgb 'red' lw 3 with steps

#blue
set tmargin at screen 0.55
set bmargin at screen 0.45
plot infile5  u 1:($3) lc -1 with boxes , '' u 1:($2) lc rgb blue_defined lw 3 with steps

unset multiplot   
unset output