source: PSPA/parmelaPSPA/trunk/distriphase.f @ 496

        subroutine distriphase(xmax,sigma,npoints)
c       for phase in input10 ran=2 and 3 (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*
        sqtwopi=sqrt(twopi)
        fac=xmax/sigma
        if(fac.lt.0.5)then
          do 10, i=1,npoints
            distp(i)=-xmax+(i-1)*2*xmax/npoints
            distp(i)=-distp(i)
 10       continue
        else
        ech=xmax/1.4142/sigma
        const=npoints/sqtwopi/erf(ech)/sigma
        sigsq=2*sigma**2
        nbin=int(min(npoints/5,100))
        if (mod(nbin,2).ne.0.)then
          nbin=nbin+1
        endif
        nhbin=nbin/2
        pas=xmax/nhbin
        k=0
        ntot=0
        do 20,i=1,nhbin
          n(i)=int(const*exp(-((nhbin-i+0.5)*pas)**2/sigsq)*pas+0.5)
          ntot=ntot+n(i)
  20    continue
        if(ntot.lt.int(npoints/2)+1)then
          nextra=int(npoints/2)+1-ntot
          do 25,l=nhbin,nhbin-nextra+1,-1
             n(l)=n(l)+1
  25      continue
        endif
        nj = 0   ! ajout glm 29/10/2012
        do 30,i=1,nhbin
          do 35,j=k+1,k+n(i)
            if(j.gt.npoints/2+1)then 
              nj=j
              go to 40
            endif
            distp(j)=(n(i)-j+k)*pas/n(i)+(nhbin-i)*pas
  35      continue
          k=k+n(i)
  30      continue
  40      continue
          if(nj.eq.0) nj=j
          do 50, i=nj,npoints
             distp(i)=-distp(2*nj-i-2)
  50      continue
        endif
        return
        end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*