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]{

i=1
#E[1]=310.
E[1]=3100.
gamma = E[i]/mpi
beta = sqrt(1. - 1/gamma**2)
ppi = gamma*beta*mpi


sleng(x) = 1+gamma*gamma_r*(beta+beta_r*cos(x))
polarization3(x)= (cos(x)*(1+gamma*gamma_r)+gamma*gamma_r*beta*beta_r)/sleng(x)
polarization1(x) = sin(x)*(gamma+gamma_r)/sleng(x)

thetapol(x)=atan(polarization1(x)/polarization3(x))

psquare(x) = sqrt((gamma(beta+beta_r*cos(x)))**2 + (beta_r*sin(x))**2)

velocity3(x) = gamma*(beta+beta_r*cos(x))/psquare(x)
velocity1(x) = beta_r*sin(x)/psquare(x)

thetao(x) = atan(velocity1(x)/velocity3(x))

sindif(x) = sin(thetaexpol(x)-thetao(x))


gammap(x) = gamma*gamma_r*(1+beta*beta_r*cos(x))
betap(x) = 1-1/gammap(x)**2

sinphi(x)=beta*gamma/gammap(x)/betap(x)*sin(x)

tantheta2(x) = tan(x/2)*(sqrt(gamma+1)*(gamma_r+1)-sqrt(gamma-1)*gamma_r*beta_r)/(sqrt(gamma+1)*(gamma_r+1)+sqrt(gamma-1)*gamma_r*beta_r)
theta2(x) = atan(tantheta2(x))
thetap(x)=2*theta2(x)


sigmadotp(x) = gamma*gamma_r*(beta+beta_r*cos(x))/(1+gamma*gamma_r*(1+beta*beta_r*cos(x)))
acossdp(x) = acos(sigmadotp(x))

set xrange [-3*pi/4:3*pi/4]