*     Collection of random number programs ... by V. Blobel

      FUNCTION RNUM()
*
*             A simple portable random number generator
*     There are only 714025 different values (spacing  =  1.4 E-6),  but
*     the period is effectively infinite!
*
*     U = RNUM()     random U(0,1) (uniform) distributed number
*     G = RNGA()     random N(0,1) (Gaussian with mean 0 and variance 1)
*
      PARAMETER (M=714025,IA=1366,IC=150889)
      INTEGER IR(97)
      LOGICAL START
      REAL RN(2)
      SAVE START,K,SN,CS,AL,RM,IY,IR,JDUM 
      DATA START/.TRUE./,K/1/
*     ...
*     random U(0,1) number
      IF(START) THEN
*        init
         START=.FALSE.
         RM=1.0/FLOAT(M)
         JDUM=-1
         JDUM=MOD(IC-JDUM,M)
         DO J=1,97
          JDUM=MOD(IA*JDUM+IC,M)
          IR(J)=JDUM
         END DO
         JDUM=MOD(IA*JDUM+IC,M)
         IY=JDUM
      END IF
*     next U(0,1) random number from array ...
      J=1+(97*IY)/M
      IY=IR(J)
      RNUM=FLOAT(IY)*RM
*     ... and store new random number in array
      JDUM=MOD(IA*JDUM+IC,M)
      IR(J)=JDUM
      GOTO 100
*     ...
      ENTRY RNGA()
*     random N(0,1) standard gaussian number
      IF(START) THEN
*        init
         START=.FALSE.
         RM=1.0/FLOAT(M)
         JDUM=-1
         JDUM=MOD(IC-JDUM,M)
         DO J=1,97
          JDUM=MOD(IA*JDUM+IC,M)
          IR(J)=JDUM
         END DO
         JDUM=MOD(IA*JDUM+IC,M)
         IY=JDUM
      END IF
*
      IF(K.LE.1) THEN
*        one U(0,1) random number
         J=1+(97*IY)/M
         IY=IR(J)
         RNDL=FLOAT(IY)*RM
         JDUM=MOD(IA*JDUM+IC,M)
         IR(J)=JDUM
*        two U(-1,+1) random numbers
   10    DO K=1,2
          J=1+(97*IY)/M
          IY=IR(J)
          RN(K)=FLOAT(2*IY)*RM-1.0
          JDUM=MOD(IA*JDUM+IC,M)
          IR(J)=JDUM
         END DO
         RADSQ=RN(1)*RN(1)+RN(2)*RN(2)
*        test point inside sircle
         IF(RADSQ.GT.1.0) GOTO 10
*        sine and cosine for random phi
         SN=RN(1)/SQRT(RADSQ)
         CS=RN(2)/SQRT(RADSQ)
*        transform to gaussians
         AL=SQRT(-2.0*ALOG(RNDL))
         RNGA=SN*AL
         K=2
      ELSE
         RNGA=CS*AL
         K=1
      END IF
  100 RETURN
      END

      FUNCTION NRPOIS(XMU)
*     random number of Poisson distribution
      REAL XLAST,EXPXMU
      SAVE XLAST,EXPXMU 
      DATA XLAST/0.0/,EXPXMU/1.0/
*     ...     
      IF(XMU.LE.0.0) THEN
         NRPOIS=0.0
         RETURN
      END IF
      IF(XMU.LT.12.0) THEN
*        small mean value
         IF(XMU.NE.XLAST) THEN
            EXPXMU=EXP(XMU)
            XLAST=XMU
         END IF
         K=0
         A=1.0
 10      A=A*RNUM()
         IF(A.GE.EXPXMU) THEN
            K=K+1
            GOTO 10
         END IF 
      ELSE
*        large value, normal approximation 
         K=MAX(0.0,XMU+RNGA()*SQRT(XMU)+0.5)  
      END IF
      NRPOIS=K
      END