file = 'tbt_-0.1'
file = 'tbt_0.001'
file = 'tbt_0.05'

amu=0.00116592
file = 'tbt_0.5'
gamma_0= 32.758453
gammaz=1./sqrt(1-betaz**2)
gammap=gamma_0/gammaz
omegac = 37.6847  # with betaz=0.5
betap = sqrt(1-1./gammap**2)
a=betap/gammaz
phi=0
f(x) = -a*sin(omegac*x)
time_step= 1.e6
fit f(x) file u ($14*time_step):($19) via a, omegac

beta_0=sqrt(1.-1/gamma_0**2)
omegas=  38.9755352065633 #omegac*(amu*gamma_0+1)
omegaa=omegas-omegac
betaz=0.5
sx(x) = cos(omegas*x) *(1-(betaz*betap*(sin(omegaa*x)))**2)**.5
sy(x) = sin(omegas*x) *(1-(betaz*betap*(sin(omegaa*x)))**2)**.5
sz(x) = ((betaz*betap)*(sin(omegaa*x) )) #/(1+(betaz*betap*(sin(omegaa*x)))**2)**.5
betax(x) = -betap/gammaz*sin(omegac*x)
betay(x) = betap/gammaz*cos(omegac*x)
betaz(x)= betaz
dsx_dt(x) = -omegas*sin(omegas*x)*(1-(betaz*betap*(sin(omegaa*x)))**2)**.5 + cos(omegas*x)*(betaz*betap)**2*omegaa*cos(omegaa*x)*sin(omegaa*x)/(1-(betaz*betap*(sin(omegaa*x)))**2)**.5
dsx_dt_rename1(x) = omegas * sy(x)
dsx_dt_rename2(x) = -sx(x)*(betaz*betap)**2*omegaa*cos(omegaa*x)/(1-(betaz*betap*(sin(omegaa*x)))**2)
d=1250
rhsx1(x)= -d*((amu+1/gamma_0)*sy(x))
rhsx2(x)= -d*(- amu*gamma_0/(gamma_0+1)*betaz*(sy(x)*betaz(x)-sz(x)*betay(x)))
reset
set xrange [0:5]

set samples 10000
set xlabel 'time [us]'
# set term qt 0; plot file u ($14*time_step):($10/100) w l,''  u ($14*time_step):11 w l, ''  u ($14*time_step):($12/100) w l lc -1
set title  'beta cartesian'
set term qt 0; plot file u ($14*time_step):($19) w l ,''  u ($14*time_step):($20) w l, ''  u ($14*time_step):($21) w l lc -1  , f(x)
#set term qt 1; plot file u ($14*time_step):($4) w l,''  u ($14*time_step):6 w l, '' u ($14*time_step):7 w l
set title 'spin cartesian'
  set term qt 1; plot file u ($14*time_step):($22) w l ,''  u ($14*time_step):($23) w l, ''  u ($14*time_step):($24) w l lc -1 , ((betaz/gammaz)*(sin(omegaa*x) ))
#set term qt 2; plot file u ($14*time_step):($13) w l
set title 's X beta'
  set term qt 2; plot file u ($14*time_step):($25) w l,''  u ($14*time_step):($26) w l, ''  u ($14*time_step):($27) w l lc -1 
set title 's X B'
  set term qt 3; plot file u ($14*time_step):($28) w l,''  u ($14*time_step):($29) w l, ''  u ($14*time_step):($30) w l lc -1 

set title 's X B analytic'
set term qt 4; plot sy(x), -sx(x)

#set term qt 5; plot cos(omegas*x)*(1-((betaz*gammaz*betap)**2*cos(omegaa*x)))
#set title 'sx/N'
#set term qt 6; plot cos(omegas*x)* (1-(0.63*betaz*gammaz*betap*(cos(omegaa*x) -1))**2)**.5,sin(omegas*x) *(1-(0.63*betaz*gammaz*betap*(cos(omegaa*x) -1))**2)**.5
set title 'spin analytic'; set term qt 6; plot sx(x) , sy(x), sz(x)
set title 'sz '
set term qt 7;  plot file  u ($14*time_step):($23) w l, ((betaz/gammaz)*(sin(omegaa*x) ))
set title 's X beta analytic'
phi=pi/2
a=.9
c=1.05
set term qt 8; plot sy(x)*betaz(x)-sz(x)*betay(x) t 'x', sz(x)*betax(x)-sx(x)*betaz(x) t 'y',sx(x)*betay(x)-sy(x)*betax(x) t 'z', file u ($14*time_step):(-$26) w l , a*sin(omegaa*x+phi) lt -1, c*(sx(x)*betay(x)-sy(x)*betax(x)) 
set title 'betax analytic'
set term qt 9; plot betax(x), betay(x)
set title 'd sx/dt'; set term qt 10; plot dsx_dt(x) , rhsx1(x)+rhsx2(x), 100*rhsx2(x), dsx_dt_rename1(x),  100*dsx_dt_rename2(x)
set title 'd sx1/dt, d sx2/dt'; set term qt 11; plot   100*rhsx2(x),   100*dsx_dt_rename2(x)
# set term qt 3; plot file u ($14*time_step):($12) w l
 #set term qt 3; plot file u ($14*time_step):($7*$12) w l
 #set term qt 0; plot file u ($14*time_step):($10/100) w l,''  u ($14*time_step):11 w l, ''  u ($14*time_step):($12/100) w l lc -1
#set xrange [0:10]
#set term qt 4; plot file u ($14*time_step):($10*sin($18)+$12*cos($18)) w l
#set term qt 5; plot file u($14*time_step):($17*cos($18)) w l