* Nom du fichier: coulomb.kumac
*
* KUMAC pour simuler la diffusion multiple de Coulomb
*                                      MP 2/95
*                             Histogrammes pour la presentation
his/cre/2d 500 'Distribution y vs. Theta ' 20 -1. 1. 20 -1. 1.
his/cre/proy 500
his/cre/prox 500
his/cre/slix 500 10
his/cre/sliy 500 10
*                             Generation de l'echantillon
ne = 10000
sigma t  = 2.            
sigma om = 2.*0.100/0.021
sigma st = sqrt(2.*t)/om
sigma sy = t*st/sqrt(3.)
sigma  y = rndm(array([ne]))*10.*sy-5.*sy
sigma th = rndm(array([ne]))*10.*st-5.*st
sigma pp = array([ne])
sigma pp = exp(-om**2*(th**2/t - 3.*th*y/t**2 + 3.*y**2/t**3))
sigma pm = rndm(array([ne]))*vmax(pp)

* create a vector >0 if pm > pp and <0 else

sigma po = pm - pp
ve/print po(1:10)

* order po

sigma j = array([ne],1#[ne])
ve/print j(1:10)
sigma k = order(j,po)
ve/print k(1:10)

* count the number of hits <0

sigma l = order(po,po)

ns = 0
repeat
       ns = [ns] + 1
until  (l([ns])>0)
ns = [ns] - 1

ve/print l(1:[ns])

* fill yo and to

vec/cre yo([ns])
vec/cre to([ns])

do i=1,[ns]
    m = k([i])     
    vec/copy y([m])  yo([i])
    vec/copy th([m]) to([i])
    call hfill(500,y([m]),th([m]),1.)
enddo    




*                             Sauvegarde de l'echantillon
vec/write yo(1:[ns]) yotst.txt
vec/write to(1:[ns]) totst.txt
*                             Sauvegarde du graphique sur fichier eps
for/file 20 'coulomb1.eps'
meta     20 -113
*              Options de presentation graphique
opt stat
*                             Presentation bidimensionelle
zone 1 1
his/plot 500                  ***** DEMO OPTIONS *****
atit 'y' 'theta'
for/close 20                                                                   
wait
*                             Distributions marginales
for/file 20 'coulomb2.eps'
meta     20 -113
zone 1 2
his/plot 500.prox
atit 'y' 'evenements'
his/plot 500.proy             
atit 'theta' 'evenements'
for/close 20
wait       
*                             Distributions conditionelles
for/file 20 'coulomb3.eps'
meta     20 -113
zone 2 5
his/plot 500.slix.1
his/plot 500.slix.2
his/plot 500.slix.3
his/plot 500.slix.4
his/plot 500.slix.5
his/plot 500.slix.6
his/plot 500.slix.7
his/plot 500.slix.8
his/plot 500.slix.9
his/plot 500.slix.10
for/close 20
wait
*                             Distributions conditionelles
for/file 20 'coulomb4.eps'
meta     20 -113
his/plot 500.sliy.1
his/plot 500.sliy.2
his/plot 500.sliy.3
his/plot 500.sliy.4
his/plot 500.sliy.5
his/plot 500.sliy.6
his/plot 500.sliy.7
his/plot 500.sliy.8
his/plot 500.sliy.9
his/plot 500.sliy.10
for/close 20
wait