* Nom du fichier: helicity.kumac
*
* KUMAC pour simuler une distrubution angulaire
*                                      MP 5/01
zone 1 2
opt stat
set stat 11
opt htit
his/cre/2d 500 'Distribution Phi vs. Theta ' 20 -3.14 3.14 20 0. 3.14
his/cre/1d 501 'Distribution g1' 50 -1. 1. 
his/cre/1d 502 'Distribution g2' 50 -1. 1.
his/cre/1d 503 'Distribution g2' 50 -1. 1.
 
* Generation de l'echantillon
ne = 100000
sigma r00 = 0.50           
sigma r11 = 0.15
sigma r10 = 0.20
fun/fun2 10 3./4./3.14*((1-.5)/2.+(3.*.5-1.)*cos(x)**2/2.-.15*sin(x)**2*cos(2.*y)-sqrt(2.)*.2*sin(2.*x)*cos(y)) 20 -3.14 3.14 20 0. 3.14
* 
*               generer echantillon
* 
sigma th = array([ne])
sigma th = (th-2.*rndm(th))*3.14
sigma ph = rndm(array([ne]))*3.14
sigma pp = (1.-r00)/2.*array([ne])+(3.*r00-1.)*cos(th)**2/2.-r11*sin(th)**2*cos(2.*ph)-sqrt(2.)*r10*sin(2.*th)*cos(ph)
sigma pp = pp*3./4./3.14
sigma pm = rndm(array([ne]))*vmax(pp)
ns = 0
vec/cre tho([ne])
vec/cre pho([ne])
do i=1,[ne]
   if (pp([i])>pm([i])) then
      ns = [ns] + 1
      vec/copy th([i]) tho([ns])
      vec/copy ph([i]) pho([ns]) 
      call hfill(500,th([i]),ph([i]),1.)
   endif
enddo
h/pl 500 lego
pic/pri helicity_1.eps
wait
* 
* projecteurs helicite
* 
vec/cre g1([ne])
vec/cre g2([ne])
vec/cre g3([ne])
sigma g1 = cos(tho)**2
sigma g2 = sin(tho)**2*cos(2.*pho)
sigma g3 = sin(2.*tho)*cos(pho)
vec/hfi g1(1:[ns]) 501
vec/hfi g2(1:[ns]) 502
vec/hfi g3(1:[ns]) 503
zone 1 3
set stat 1111
h/pl 501
h/pl 502
h/pl 503
pic/pri helicity_2.eps
wait