#set term pdf fontscale 0.75 size 6in,4in enhanced
#outfile='sigma_perp_vs_theta.pdf'
#set output outfile 
set timestamp font 'Verdana,6' #offset graph 1.13,0.7 rotate font 'Verdana,6'
set label 4 'sigma_perp_vs_theta.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



mpi= 139.57
mmu= 105.658
b=mmu/mpi
pmu_r = (mpi**2-mmu**2)/2./mpi
Emu_r = sqrt(pmu_r**2+mmu**2)
gamma_r = Emu_r/mmu
beta_r = pmu_r/Emu_r
set samples 10000
array E[3]
E[1] = 310.; E[2]=3100.; E[3]=31000.

do for [i=1:3]{
gamma = E[i]/mpi
beta = sqrt(1. - 1/gamma**2)
ppi = gamma*beta*mpi

coshxp = sqrt(gamma+1)
sinhxp = sqrt(gamma-1)
coshx = sqrt(gamma_r+1)
sinhx = sqrt(gamma_r-1)

spin1(x) = (coshxp*coshx + sinhxp*sinhx)*cos(x/2)
spin2(x) = (coshxp*coshx - sinhxp*sinhx)*sin(x/2)
spin3(x) = (-sinhxp*coshx - coshxp*sinhx)*cos(x/2)
spin4(x) = (sinhxp*coshx - coshxp*sinhx)*sin(x/2)

tan_thetap_s(x) = spin2(x)/spin1(x)
tan_thetap_o(x) = -spin4(x)/spin3(x)

tan_thetap_direct(x) = sin(x)/(gamma*(beta+beta_r*cos(x)))

thetap_s(x) = atan(tan_thetap_s(x))*2
thetap_o(x) = atan(tan_thetap_o(x))*2

dtheta(x) = thetap_s(x) - thetap_o(x)
sin_spin_theta(x) = sin(dtheta(x))

snux(x)=sin(x); snuz(x)=cos(x)

gammap(x) = gamma_r*gamma*(1-beta*beta_r*cos(x))

beta_rx(x) = beta_r*sin(x)
beta_rz(x) = beta_r*cos(x)
gammap(x) = gamma_r*gamma*(1+beta*beta_rz(x))

g2_g1 = gamma_r**2/(1+gamma_r)

betaxp(x)=gamma_r*beta_r*sin(x)/gammap(x)
betazp(x)=gamma_r*gamma*(beta+beta_r*cos(x))/gammap(x)

betap(x)=sqrt(betaxp(x)**2+betazp(x)**2)

nbetaxp(x) = betaxp(x)/betap(x)
nbetazp(x) = betazp(x)/betap(x)

sx(x) = cos(x)*nbetaxp(x)*nbetazp(x)*(gammap(x)-1) + 0.5*sin(x)*((gammap(x)+1)+(nbetaxp(x)**2-nbetazp(x)**2)*(gammap(x)-1))
sz(x) = 0.5*cos(x)*((gammap(x)+1)+(nbetazp(x)**2-nbetaxp(x)**2)*(gammap(x)-1))+sin(x)*nbetaxp(x)*nbetazp(x)*(gammap(x)-1)

slength(x)=sqrt(sx(x)**2+sz(x)**2)
sin_alpha(x) =(sx(x)*betazp(x)-sz(x)*betaxp(x))/slength(x)/betap(x)

svecx(x) = (1+(gammap(x)-1)*nbetaxp(x)**2)*sin(x) + (gammap(x)-1)*nbetaxp(x)*nbetazp(x)*cos(x)
svecz(x)= (gammap(x)-1)*nbetaxp(x)*nbetazp(x)*sin(x) + (1+(gammap(x)-1)*nbetazp(x)**2)*cos(x)

tan_alpha_w(x) = sin(x)/gamma_r/(cos(x)+beta_r/beta)
sin_alpha_w(x) = tan_alpha(x)/(1+tan_alpha(x)**2)**.5
 
length_svec(x) = sqrt(svecx(x)**2+svecz(x)**2)
sin_alpha_vec(x)= (svecx(x)*betazp(x)-svecz(x)*betaxp(x))/length_svec(x)/betap(x)




#pmux(x)=pmu_r*sin(x); pmuz(x)=gamma*(beta*Emu_r+pmu_r*cos(x))
#sperp(x) = (snux(x)*pmuz(x)-snuz(x)*pmux(x))/(pmux(x)**2+pmuz(x)**2)**.5


xx(x) = gammap(x)*betazp(x)*mmu/(beta*gamma*mpi)
scpperp(x) = 2*b/xx(x)/(1-b**2)*sqrt((1-xx(x))*(xx(x)-b**2))*x/abs(x)

print gamma_r, beta_r, beta,gamma
#sdotp(x)=gamma_r**2*(-beta_r*sin(x)*sin(x)+gamma**2*(-beta_r*beta**2+(beta_r**2*beta+beta)*cos(x) - beta_r*cos(x)**2))
#sp(x)=gamma_r**2*(sin(x)**2+gamma**2*(-beta_r*beta+cos(x))**2)**.5 * (beta_r**2*sin(x)**2+gamma**2*(beta-beta_r*cos(x))**2)**.5
#costheta(x)=sdotp(x)/sp(x)


#sdotp(x)=gamma_r**2*(-beta_r*sin(x)*sin(x)+gamma*(beta*cos(x) - beta_r*cos(x)**2))
#sp(x)=gamma_r**2*(sin(x)**2+cos(x)**2)**.5 * (beta_r**2*sin(x)**2+gamma**2*(beta-beta_r*cos(x))**2)**.5

set key bottom r L
set label 7 at graph 0.05,0.95 'E_{/Symbol p}='.sprintf("%3.1f", E[i]).' MeV'
set label 5 at graph 0.05,0.87 '{/Symbol g}='.sprintf("%3.1f",gamma)
set label 6 at graph 0.05,0.8 '{/Symbol b}='.sprintf("%4.3f",beta)
set xrange [-pi:pi]
set yrange [-1:1.2]
set grid
set xlabel '{/Symbol q} [rad]'
set xtics ("-pi" -pi, "-pi/2" -pi/2, "0" 0, "pi/2" pi/2, "pi" pi)
set ylabel '{/Symbol S}_T'
plot sin_alpha(x) t '{/Symbol S}_T (dlr)', scpperp(x) t '{/Symbol S}_T (cp)', xx(x) t 'x=p^{/Symbol m}_{||}/p_{/Symbol p}',sin_alpha_w(x) w lp t 'werbrouck', sin(thetap_s(x)-x) lt -1 t 'spinor - orbit'
pause -1

set xrange [-pi/8:pi/8]
set xtics -pi/8,pi/40
set autoscale y
plot sin_alpha(x) t '{/Symbol S}_T (dlr)',sin_alpha_vec(x) t '4-vector'
pause -1
#set xrange [-pi:0]
#set autoscale y
#plot sdotp(x)
#plot sp(x)
#plot costheta(x)
#plot (1-costheta(x)**2)**.5*x/abs(x)

reset
set xrange [-pi:pi]
plot sx(x) t 'spinor', svecx(x) t 'vector'
pause -1
}

#print outfile