source: PSPA/parmelaPSPA/trunk/distrir.f @ 13

        subroutine distrir(xmax,sigma,npoints)
c       for x-y in input10 ran=2  (LAL)
c-----------------------------------------------------------------------
c
        include 'param_sz.h'
        include 'constcom.h'
c
        common/dist/distp(imaa),xr(imaa),yr(imaa)
        dimension n(imaa/5)
c-----------------------------------------------------------------------
c*
        fac=xmax/sigma
        xmaxsq=xmax**2
        sigsq=sigma**2
        const1=2*npoints/xmaxsq
        const2=npoints/(sigsq*(1.-exp(-xmaxsq/(2*sigsq))))
        nbin=int(min(npoints/5,100))
        pas=xmax/nbin
        ntot=0
c-------
        do 20 i=1,nbin
          if(fac.lt.0.5) then
c  distribution uniforme
              n(i)=int(const1*(i-0.5)*pas*pas+0.5)
              ntot=ntot+n(i)
              if(ntot.gt.npoints) then
                 n(i)=n(i)-ntot+npoints
                 ntot=npoints
                 go to 21
              endif
          else
c  distribution gaussienne
              n(i)=int(const2*exp(-((i-0.5)*pas)**2/(2*sigsq))
     1             *(i-0.5)*pas*pas+0.5)
          ntot=ntot+n(i)
              if(ntot.gt.npoints) then
                 n(i)=n(i)-ntot+npoints
                 ntot=npoints
                 go to 21
              endif
          endif
  20    continue
c-------
        if(ntot.lt.npoints)then
          nextra=npoints-ntot
          do 25 l=1,nextra
             n(l)=n(l)+1
  25      continue
        endif
  21    continue
        k=0
        do 30 i=1,nbin
          if(n(i).eq.0) go to 30
          theta=twopi/n(i)
          theta0=i*twopi/nbin
          do 35,j=k+1,k+n(i)
              xr(j)=(i-0.5)*pas*cos((j-1)*theta+theta0)
              yr(j)=(i-0.5)*pas*sin((j-1)*theta+theta0)
  35      continue
          k=k+n(i)
  30      continue
        return
        end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*