#run energy_spin_vs_time to create all_Energy_spin_time.dat, all_Energy_vs_time_spin.dat and phase_vs_time.dat in all config _## directories
array subdir[5]
array tit[5]

npoints=5
tit[3]='{/Symbol D}B/B=0'
tit[1]='{/Symbol D}B/B=-0.007'
tit[2]='{/Symbol D}B/B=-0.0035'
tit[4]='{/Symbol D}B/B=+0.0035'
tit[5]='{/Symbol D}B/B=+0.007'

#subdir[2] = 'config_078'; subdir[1]='config_079';subdir[3]='config_080'
#subdir[2] = 'config_072'; subdir[1]='config_074';subdir[3]='config_075';config_slope = -13.1
#subdir[2] = 'config_073'; subdir[1]='config_076';subdir[3]='config_077';config_slope = 0.79

#subdir[2] = 'config_072'; subdir[1]='config_074a';subdir[3]='config_075a';config_slope = -13.1
subdir[3] = 'config_072'; subdir[1]='config_074a';subdir[2] = 'config_074b';subdir[4]='config_075b';subdir[5]='config_075a';config_slope = -13.1
subdir[3] = 'config_081'; subdir[1]='config_082a'; subdir[2]='null';subdir[4]='null';subdir[5]='config_083a';config_slope = -1.72


set term pdf fontscale 0.5 size 5in,3in
#set term x11
outfile = subdir[3].'/phase_vs_energy.pdf'
set output outfile

#set timestamp offset graph 1.18,0.7 rotate font 'Verdana,6'
set timestamp font 'Verdana,6'
set label 1 'plotting_scripts/phase_vs_energy.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
!pwd > directory
!sed -n -e 's/^/name="/' -e 's/$/"/p' directory > name
load 'name'
print 'directory='.name
set label at graph 0.0,1.02 name font 'Verdana,6' noenhance

#plot energy of stored particles in each of subdirectories and compute average
#measured phase vs momentum
run1_slope = -1.; run1_slope_err = 0.16
latest_slope = -1.35; latest_slope_err = 0.14


array average_e[npoints]
array a[npoints]
array b[npoints]
array omega[npoints]

a[1] = 2; a[2]=0; a[3]=2; a[4]=0; a[5]=2


do for [i=1:npoints]{
if(a[i]!=0){
set ylabel 'dp/p'
set xlabel 'muon'
set xrange [0:]
set key left Left
plot subdir[i].'/all_Energy_spin_time.dat' u 3 t tit[i]                                    # energy of stored beam
stats subdir[i].'/all_Energy_spin_time.dat' u 3                                      # energy of stored beam
average_e[i]= STATS_mean

#pause -1

b[i]=1.4
f(x)  = a[i] + b[i]*x
set title subdir[i] noenhance

fit f(x)  subdir[i].'/phase_vs_time.dat' u ($1*1.e6):2 via a[i], b[i]
set ylabel '{/Symbol f} [rad]'
set xlabel 'time [us]'
plot subdir[i].'/phase_vs_time.dat' u ($1*1.e6):2 not, f(x) t tit[i]
#pause -1
omega[i]=b[i]
print subdir[i],' ',a[i],' rad','    ', b[i], 'X 10^6  [rad/s]'
}
}

f(x) = aa+bb*x

set autoscale
set key right Left
set title subdir[3] noenhance
set xlabel '{/Symbol D}B/B'
set ylabel '{/Symbol Df}(t=0) [rad]'
set yrange [-0.04:0.02]
set ytics -0.04,0.005
fit f(x) a u (($1-3.)*0.0035):($2 != 0 ? ($2-a[3]):1/0) via aa,bb
g(x) = run1_slope *x
h(x) = latest_slope*x
k(x) = config_slope*x
plot a u (($1-3.)*0.0035):($2 != 0 ? ($2-a[3]):1/0) pt 4  t 'BMAD "measured": d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", bb).' [rad]',f(x) lt 1 not
plot a u (($1-3.)*0.0035):($2 != 0 ? ($2-a[3]):1/0) pt 4  t 'BMAD "measured": d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", bb).' [rad]',f(x) lt 1 not,\
g(x) t 'LaBounty Run1: d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", run1_slope).' [rad]' ,\
h(x) t 'LaBounty Now:  d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", latest_slope).' [rad]' 

set autoscale y
set yrange [-0.05:0.05]
set ytics -0.1,0.02
plot a u (($1-3.)*0.0035):($2 != 0 ? ($2-a[3]):1/0) pt 4  t 'BMAD "measured": d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", bb).' [rad]',f(x) not,\
g(x) t 'LaBounty Run1: d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", run1_slope).' [rad]' ,\
h(x) t 'LaBounty Now:  d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", latest_slope).' [rad]' ,\
k(x) lt -1 dt 2 t 'Stored beam: d{/Symbol f}/d{/Symbol d}} ='.sprintf("%3.1f", config_slope).' [rad]' 

set autoscale
set key left Left
set title subdir[3] noenhance
set ylabel '{/Symbol Dw/w}'
set yrange [-0.01:0.01]
set ytics -0.01,0.005
set xlabel '{/Symbol D}B/B'
fit f(x) b u (($1-3)*0.0035):($2 !=0 ? (($2-b[3])/b[3]):1/0) via aa,bb
plot b u (($1-3)*0.0035):($2 !=0 ? (($2-b[3])/b[3]):1/0) pt 4 not, f(x) t '(d{/Symbol w}/d{/Symbol w})/d{/Symbol d} ='.sprintf("%3.1f", bb)

set autoscale
set title subdir[3] noenhance
set ylabel '<dp/p>'
set yrange [-0.01:0.01]
set ytics -0.01,0.005
set xlabel '{/Symbol D}B/B'
fit f(x) average_e u (($1-3)*0.0035):($2 != 0 ? $2:1/0) via aa,bb
plot average_e u (($1-3)*0.0035):($2 != 0 ? $2:1/0) pt 4 not, f(x)  t 'd{/Symbol d}/d({/Symbol D}B/B) ='.sprintf("%3.1f", bb)
do for [i=1:5]{print subdir[i],a[i],' rad','    ', b[i], 'X 10^6  [rad/s]'}
print outfile