* 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 = 100000
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)
ns = 0
vec/cre yo(5000)
vec/cre to(5000)
do i=1,[ne]
*  sigma print th([i]),y([i]),pp([i]),pm([i])
   if (pp([i])>pm([i])) then
      ns = [ns] + 1
      vec/copy y([i])  yo([ns])
      vec/copy th([i]) to([ns]) 
      call hfill(500,y([i]),th([i]),1.)
   endif
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