set term pdf fontscale 0.5 size 6in,2in
set timestamp offset graph 1.111,0.7 rotate font ',6'
set label 2 'example10/fixed_pitch/plot_longitudinal.gnu' at graph 1.01,0.02 rotate left font ',6' noenhance

file_long_0percent = 'tbt_long_0percent.'
file_long_1percent = 'tbt_long_1percent.'
file_long_10percent = 'tbt_long_10percent.'

amu=0.00116592
c_light=2.99792458e2
R=7.112
gamma0= 29.303464
beta= (1.-1/gamma0**2)**.5
omega = beta*c_light/R

omega0=omega*gamma0*(amu+1/gamma0)
omegaa = omega *gamma0*amu
blz = 0.01 #Bl/Bz
bomega(blz)=omegaa*blz*(1/gamma0*(1+amu))/amu
omegap(blz) = ((omega-omega0)**2+bomega(blz)**2)**.5
beta = (1-1./gamma0**2)**.5
beta_x(x)=beta*cos(omega*x)
beta_y(x)=-beta*sin(omega*x)
beta_z=0

sx(x,blz)= ((cos(omegap(blz) * x/2))**2-((omega-omega0)**2-bomega(blz)**2)/(omegap(blz)**2) *(sin(omegap(blz) * x/2))**2)*cos(omega *x)+(2/omegap(blz)*cos(omegap(blz)*x/2)*sin(omegap(blz)*x/2)*(omega-omega0))*sin(omega*x)
sx1(x,blz)= ((cos(omegap(blz) * x/2))**2-((omega-omega0)**2-bomega(blz)**2)/(omegap(blz)**2) *(sin(omegap(blz) * x/2))**2)*cos(omega *x)
sx2(x,blz)= (2/omegap(blz)*cos(omegap(blz)*x/2)*sin(omegap(blz)*x/2)*(omega-omega0))*sin(omega*x)
sy(x,blz)= -((cos(omegap(blz) * x/2))**2-((omega-omega0)**2-bomega(blz)**2)/(omegap(blz)**2) *(sin(omegap(blz) * x/2))**2)*sin(omega *x)+(2/omegap(blz)*cos(omegap(blz)*x/2)*sin(omegap(blz)*x/2)*(omega-omega0))*cos(omega*x)
sz(x,blz) =bomega(blz)/omegap(blz)*(sin(omegap(blz)*x))/2

sdotb(x,blz) = -(sx(x,blz)*beta_x(x)+ sy(x,blz)*beta_y(x)+ sz(x,blz)*beta_z)
set yrange [-1.2:1.2]
set xrange [0:50]

 #plot 'tbt_long_1percent.' u (1.e6*$14):22 w l, '' u (1.e6*$14):13 w l, 'tbt_long_1percent.' u (1.e6*$14):22 w l, '' u (1.e6*$14):13 w l, 'tbt_long_0percent.' u (1.e6*$14):22 w l, '' u (1.e6*$14):13 w l
# plot 'tbt_long_1percent.'  u (1.e6*$14):13 w l, 'tbt_long_1percent.'  u (1.e6*$14):13 w l, 'tbt_long_0percent.'  u (1.e6*$14):13 w l

# plot 'tbt_long_1percent.'  u (1.e6*$14):($13/($19**2+$20**2+$21**2)**.5) w l, 'tbt_long_1percent.'  u (1.e6*$14):13 w l, 'tbt_long_0percent.'  u (1.e6*$14):13 w l
set xlabel 'time [us]'
set ylabel 'Amplitude'
set samples 100000
set label 1 at graph 0.02, 0.95 'simulation'  rotate by 0 font ',10'
#set term qt 1
outfile = 'plot_longitudinal1.pdf'
set output outfile
plot 'tbt_long_1percent.'  u (1.e6*$14):($13/($19**2+$20**2+$21**2)**.5) w l t 'B_l=B_z/100', 'tbt_long_0percent.'  u (1.e6*$14):13 w l t 'B_l=0'
print 'write '.outfile
set label 1 at graph 0.02, 0.95 'analytic' rotate by 0 font ',10'
#set term qt 2;
outfile = 'plot_longitudinal2.pdf'
set output outfile
phi = 0.
plot sdotb(x+phi,0.01) t 'B_l=B_z/100', sdotb(x+phi,0.) t 'B_l=0'
print 'write '.outfile
unset label 1
#set term qt 3;
outfile = 'plot_longitudinal3.pdf'
set output outfile

plot 'tbt_long_1percent.'  u (1.e6*$14):(int($14*1.e8/10)*10 == int($14*1.e8) ? ($13/($19**2+$20**2+$21**2)**.5):1/0) w p lc 1 pt 6 ps 0.25 t 'B_l=B_z/100 - simulation', sdotb(x,0.01) lc 2 t 'B_l=B_z/100 analytic' 
# 'tbt_long_0percent.'  u (1.e6*$14):13 w l lc 1 t 'B_l=0 simulation', -sdotb(x,0.) lc 2 dt 2  t 'B_l=0 analytic'


print 'write '.outfile