# input filename. Could be 'domega_ddelta_vs_time.dat' or 'combined_domega_ddelta_vs_time.dat'
# run energy_vs_time to compute average and rms momentum. Reads 'all_Energy_vs_time_0.dat'
set term pdf fontscale 0.75 size 6in,3in
#set term x11
outfile = 'dphi_ddelta_vs_time.pdf'
set output outfile

set timestamp offset graph 1.18,0.7 rotate font 'Verdana,6'
set label 4 'plotting_scripts/dphi_ddelta_vs_time.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 1 at graph 0.0, 1.02 name font 'Verdana,6' noenhance
set xlabel 'time [{/Symbol m}s]'
set ylabel '{/Symbol D}p/p [%]'
plot 'all_Energy_vs_time_0.dat' u ($3*1.e9):($2*100) not
stats 'all_Energy_vs_time_0.dat' u ($3*1.e9):($2*100)

#pause -1


 tau=2.197; gamma_0=29.3; t_rev=0.1492
stdev_energy = STATS_stddev_y/100.
average_energy=STATS_mean_y/100.

 dgamma_dt = stdev_energy**2/(1+average_energy)**2/tau
 omega_a = 2*pi/t_rev * gamma_0/(gamma_0**2-1)
print ' dgamma_dt =', dgamma_dt
set ylabel 'd{/Symbol f}/d{/Symbol d} [rad]'
set xlabel 'time [{/Symbol m}s]'
set xrange [5:22]

stats filename u ($1*1.e6):2
average_dphi_ddelta = STATS_mean_y

set yrange [-30:0]#config_067 [0:20]#config_065 set yrange [-20:15]#config_049   [-20:0]#config_050  [-10:20]#config_048[-10:10]#config_051 [-20:0]#config_045 
set xrange [0:22]
set label 2 at graph 0.63,0.95 '<{d{/Symbol f}/d{/Symbol d}>='.sprintf("%3.1f",average_dphi_ddelta).'[rad]'
plot filename u ($1*1.e6):2:3 w e not
set xrange [:10]
plot filename u ($1*1.e6):2:3 w e not

set xrange [18:20]
plot filename u ($1*1.e6):2:3 w e not


#pause -1

set ylabel '{/Symbol D w}/{/Symbol w} [ppb]' 

set xrange [5:22]
set yrange [-300:100]#config_066 set yrange [-100:300]#config_065 set yrange [-300:100]#config_050 [-200:100]#config_049    [-100:200]#config_048 [-150:150]#config_051 [-300:100]#config_045 [-150:150]#config_047 
set xrange [18:20]
stats filename u ($1*1.e6):($2*dgamma_dt/gamma_0/omega_a*1.e9)
average_domega_omega = STATS_mean_y
set autoscale x
set xrange [5:22]

unset label 2
set label 3 at graph 0.63,0.95 '<{/Symbol Dw}_a/{/Symbol w}_a>='.sprintf("%3.1f",average_domega_omega).'[ppb]'
plot filename u ($1*1.e6):($2*dgamma_dt/gamma_0/omega_a*1.e9):($3*dgamma_dt/gamma_0/omega_a*1.e9) w e not

set xrange [0:10]
plot filename u ($1*1.e6):($2*dgamma_dt/gamma_0/omega_a*1.e9):($3*dgamma_dt/gamma_0/omega_a*1.e9) w e not
set xrange [18:20]
plot filename u ($1*1.e6):($2*dgamma_dt/gamma_0/omega_a*1.e9):($3*dgamma_dt/gamma_0/omega_a*1.e9) w e not

load 'fort.401'
reset
#set yrange [2:6]#config_048 #[-30:15]#config_050 
set xrange [15:20]
set xlabel 'time [us]'
set ylabel 'd{/Symbol f}/d{/Symbol d} [rad]'
set pointsize 0.2
plot f01 u ($1*1.e6):2 not, f02 u ($1*1.e6):2 not, f03 u ($1*1.e6):2 not, f04 u ($1*1.e6):2 not, f05 u ($1*1.e6):2 not, f06 u ($1*1.e6):2 not, f07 u ($1*1.e6):2 not, f08 u ($1*1.e6):2 not #, f09 u ($1*1.e6):2 not, f10 u ($1*1.e6):2 not\
f11 u 1:2, f12 u 1:2,f13 u 1:2, f14 u 1:2,f15 u 1:2, f16 u 1:2,

print ' t_rev[us] = ', t_rev
print ' dgamma_dt =', dgamma_dt
print ' omega_a = ', omega_a*1.e6
print ' dgamma_dt/gamma_0/omega_a[1/rad] = ', dgamma_dt/gamma_0/omega_a*1.e9,' [ppb/rad]' 



print outfile