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

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

snux(x)=sin(x); snuz(x)=cos(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) = (Emu_r+pmu_r/beta*cos(x))/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 sperp(x) t '{/Symbol S}_T (dlr)', scpperp(x) t '{/Symbol S}_T (cp)', xx(x) t 'x=p^{/Symbol m}_{||}/p_{/Symbol p}'

#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)
}

print outfile