       PROGRAM POISSONdr

       INTEGER idum,n,j
       REAL xm,p,xx

C       open(5,file='pered.dat')

       write(*,*)'n:'
       read(*,*)n
       write(*,*)'xm(valos):'
       read(*,*)xm

       do j=-1,-n,-1
	idum=j
       p=poidev(xm,idum)
       write(*,*)p
       enddo
C       close(5)
       stop
       end



       FUNCTION poidev(xm,idum)
       REAL PI,poidev
       PARAMETER (PI=3.141592654)
       REAL alxm,em,g,oldm,sq,t,y,gammln,ran1
       SAVE alxm,g,oldm,sq
       DATA oldm /-1/
       if (xm.lt.12)then
	   if (xm.ne.oldm)then
	      oldm=xm
	      g=exp(-xm)
	   endif
	   em=-1
	   t=1
 2         em=em+1
	   t=t*ran1(idum)
	   if (t.gt.g)goto2
       else
	   if (xm.ne.oldm)then
	     oldm=xm
	     sq=sqrt(2*xm)
	     alxm=log(xm)
	     g=xm*alxm-gammln(xm+1)
	   endif
 1         y=tan(PI*ran1(idum))
	   em=sq*y+xm
	   if(em.lt.0)goto1
	   em=int(em)
	   t=0.9*(1+y**2)*exp(em*alxm-gammln(em+1)-g)
	   if (ran1(idum).gt.t)goto1
       endif
       poidev=em
       return

       END



       FUNCTION gammln(xx)
       REAL gammln
       INTEGER j
       DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
       SAVE cof,stp
       DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     * 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     * -.5395239384953d-5,2.5066282746310005d0/
       x=xx
       y=x
       tmp=x+5.5d0
       tmp=(x+0.5d0)*log(tmp)-tmp
       ser=1.000000000190015d0
       do j=1,6
	y=y+1.d0
	ser=ser+cof(j)/y
       enddo
       gammln=tmp+log(stp*ser/x)
       return
       END



      FUNCTION ran1(idum)

      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     * NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv/NTAB*0/,iy/0/
      if (idum.le.0.or.iy.eq.0)then
	  idum=max(-idum,1)
	  do j=NTAB+8,1,-1
	    k=idum/IQ
	    idum=IA*(idum-k*IQ)-IR*k
	    if (idum.lt.0)idum=idum+IM
	    if (j.le.NTAB)iv(j)=idum
	    enddo
	    iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0)idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
c      return
      END










