1 |        subroutine tirb(p,numbande)
   2 | c.......................................................................
   3 | c     anciennement ntir
   4 | c     distribution d'un nombre de paquets a tirer sur les bandes
   5 | c     en fonction d'une probabilite
   6 | c
   7 | c     tirage du numero numbande conformement a la proba discrete p
   8 | c
   9 | c
  10 | c.......................................................................
  11 | c
  12 | c  in :  * probabilites (une fonction == vecteur)    ---> p(nband_mx)
  13 | c
  14 | c  out : * numero de la bande tiree bande            ---> numbande
  15 | c
  16 | c.......................................................................
  17 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  18 |       IMPLICIT NONE
  19 | c.......................................................................
  20 |       include 'cecile.inc'
  21 | c     cecile.inc sert a la declaration de ngaz_mx
  22 |       INCLUDE 'propradia.inc'
  23 | c.......................................................................
  24 |       DOUBLE PRECISION p(nbande_mx)
  25 | c.......................................................................
  26 |       DOUBLE PRECISION sp(nbande_mx)
  27 | c.......................................................................
  28 |       DOUBLE PRECISION tmp1
  29 |       INTEGER itmp1,itmp2,itmp3,ibandeloc, numbande
  30 | c.......................................................................
  31 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  32 | c construction de la fonction de repartition discrete de p ==> sp
  33 | 
  34 |       sp(1)=p(1)
  35 |       DO ibandeloc=2,nbande
  36 |         sp(ibandeloc)=sp(ibandeloc-1)+p(ibandeloc)
  37 |       ENDDO
  38 | 
  39 | c tirage avec sp pour trouver numbande
  40 |       call rand_uniforme(tmp1)
  41 |       itmp1=0
  42 |       itmp2=nbande+1
  43 | 10000 CONTINUE
  44 |       IF (itmp2-itmp1.GT.1) THEN
  45 |           itmp3=(itmp2+itmp1)/2
  46 |           IF (tmp1.GT.sp(itmp3)) THEN
  47 |             itmp1=itmp3
  48 |           ELSE
  49 |             itmp2=itmp3
  50 |           ENDIF
  51 |       GOTO 10000
  52 |       ENDIF
  53 |       numbande=itmp1+1
  54 | 
  55 | c modif du poids
  56 | c aucune precaution si p(numbande est nul)
  57 |       poids = poids / (nbande*p(numbande))
  58 | c      if (in.eq. 110) then
  59 | c         print *, 'poids, nbande, p(numbande)'
  60 | c         print *, poids, nbande, p(numbande)
  61 | c         read *
  62 | c         endif
  63 |      
  64 | c.......................................................................
  65 |       RETURN
  66 |       END