amu=0.00116592
gamma0= 29.3
betaz_0=0.6
betaz=betaz_0
gammaz=1./sqrt(1-betaz**2)

gammap=gamma0/gammaz
kappa = (amu+1/(1+gamma0))

beta=(1-1/gamma0**2)**.5
clight=2.99792458e2
R=7.112
omega = clight/R  # with betaz=0.5
omega0(betaz)=omega*gamma0*(amu+1/gamma0-amu*gamma0/(1+gamma0)*betaz**2)
omegaa = omega *gamma0*amu
bomega(betaz)=omegaa*gamma0/(gamma0+1)*betaz*beta
omegap(betaz) = ((omega-omega0(betaz))**2+bomega(betaz)**2)**.5
omegav = 0.34*omega

beta_x(x)=beta*cos(omega*x)
beta_y(x)=beta*sin(omega*x)

beta_xye(x,gamma0) = (1-1/gamma0**2)**.5*cos(betaz_0*sin(omegav*x))
beta_ze(x) = sin(betaz_0*sin(omegav*x))


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

sxy_ve(x,gamma0) = cos(kappa*gamma0*(1-1/gamma0**2)**.5*betaz_0*sin(omegav*x))
sz_ve(x,gamma0) = sin(kappa*gamma0*(1-1/gamma0**2)**.5*betaz_0*sin(omegav*x))

sdotb(x,betaz) = sx(x,betaz)*beta_x(x) + sy(x,betaz)*beta_y(x) + sz(x,betaz)*beta_z
bz(x) = betaz_0*sin(omegav*x)
bpi(x) = omegaa*(amu+1/(gamma0+1))*omegav/omega * bz(x) * beta
omegap_me(betaz) = ((omega-omega0(betaz))**2+bomega(betaz)**2 + bpi(betaz)**2)**.5

sx_v(x)= ((cos(omegap_me(bz(x)) * x/2))**2-((omega-omega0(bz(x)))**2-(bomega(bz(x))**2+bpi(x))**2))/(omegap_me(bz(x))**2) *(sin(omegap_me(bz(x))*x/2))**2)*cos(omega *x)-(2/omegap_me(bz(x))*cos(omega*x/2)*sin(omega*x/2)*(omega-omega0(bz(x))))*sin(omega*x)
sy_v(x)= ((cos(omegap_me(bz(x)) * x/2))**2-((omega-omega0(bz(x)))**2-(bomega(bz(x))**2+bpi(x)**2))/(omegap_me(bz(x))**2) *(sin(omegap_me(bz(x))*x/2))**2)*sin(omega *x)+(2/omegap_me(bz(x))*cos(omega*x/2)*sin(omega*x/2)*(omega-omega0(bz(x))))*cos(omega*x)
sz_v(x) = 2*(omega-omega0(bz(x)))*bomega(bz(x))/omegap_me(bz(x))**2*(sin(omegap_me(bz(x))*x/2))**2
sdotb_v(x) = (sx_v(x)*beta_x(x) + sy_v(x)*beta_y(x))*beta_xye(x,gamma0) + sz_v(x)*bz(x)*beta_ze(x)
lengths(x) =((beta*cos(bz(x)))**2 + (sin(bz(x)))**2)**.5

sdotbv_v(x,gamma0) = (sxy_ve(x,gamma0)*beta_xye(x,gamma0) + sz_ve(x,gamma0)*beta_ze(x))/lengths(x)