function erf(z)
       implicit none
       real*8 erf,z,gammp
c returns the error function erf(x) with fractional error everywhere less
c than 1.2E-07. From Numerical Recipes routine erfcc.f
c      erf=gammp(.5d0,z**2)
       real*8 t
c for z>0 only
       t=1.d0/(1.d0+0.5d0*z)
       erf=1.d0-t*dexp(-z*z-1.26551223d0+t*(1.00002368d0+t*(.37409196d0+
     &  t*(.09678418d0+t*(-.18628806d0+t*(.27886807d0+t*(-1.13520398d0+
     &  t*(1.48851587d0+t*(-.82215223d0+t*.17087277d0)))))))))
       return
       end

       function erfc(z)
       implicit none
       real*8 erfc,z,gammp
c returns the complementary error function erf(x) with fractional 
c error everywhere less
c than 1.2E-07. From Numerical Recipes routine erfcc.f
c      erfc=gammp(.5d0,z**2)
       real*8 t
c for z>0 only
       t=1.d0/(1.d0+0.5d0*z)
       erfc=t*dexp(-z*z-1.26551223d0+t*(1.00002368d0+t*(.37409196d0+
     &  t*(.09678418d0+t*(-.18628806d0+t*(.27886807d0+t*(-1.13520398d0+
     &  t*(1.48851587d0+t*(-.82215223d0+t*.17087277d0)))))))))
       return
       end
       

      FUNCTION BESSI0(X)
	implicit none
	REAL*8 x,Y,P1,P2,P3,P4,P5,P6,P7,
     *	  Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,bessi0,ax
        DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0,1.2067492D
     *	  0,
     *	  0.2659732D0,0.360768D-1,0.45813D-2/
	DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1,
     *	  0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1,
     *	  0.2635537D-1,-0.1647633D-1,0.392377D-2/
	IF (ABS(X).LT.3.75) THEN
	Y=(X/3.75)**2
	BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
	BESSI0=BESSI0*DEXP(-X)
	ELSE
	AX=ABS(X)
	Y=3.75/AX
	BESSI0=(1.d0/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
     *	    +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
	ENDIF
      RETURN
      END

		FUNCTION BESSI1(X)
	 implicit none
	 REAL*8 x,Y,P1,P2,P3,P4,P5,P6,P7,
     *	  Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,bessi1,ax
	 DATA P1,P2,P3,P4,P5,P6,P7/0.5D0,0.87890594D0,0.51498869D0,
     *	  0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/
	 DATA Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,-0.3988024D-1,
     *	  -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1,
     *	  -0.2895312D-1,0.1787654D-1,-0.420059D-2/
	 IF (ABS(X).LT.3.75) THEN
	 Y=(X/3.75)**2
	 BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
	 BESSI1=BESSI1*DEXP(-X)
        ELSE
	 AX=ABS(X)
	 Y=3.75/AX
	 BESSI1=(1.d0/SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+
     *	    Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
        ENDIF
      RETURN
      END

	function gaussn()
c	real *8 r1, r2, dlog, dcos, dsqrt
	real *8 r1, r2
	real*8 gaussn
10	continue
        call rmar(r1)
        call rmar(r2)
	 if(abs(r1).lt.0.0000000000001d0) goto 10
	gaussn= dsqrt(-2.0d0*dlog(r1))*dcos(6.2831853d0*r2)
c	print *,r1,r2,gaussn
	 if(dabs(gaussn).gt.100.d0) stop
	return
	end


      subroutine rmaset(iunit)
c     initializing routine for rmar, must be called before 
c     beginning to generate pseudorandom numbers by means of the
c     subroutine rmar
c     ranges : 0<=ij<=31328 and 0<=kl<=30081.
c     implicit real*8 (a-h,o-z)
      implicit none
      real*8 u,c,cd,cm,s,t
      integer i,j,k,m,iunit,n,kl,ij,ii,jj

      common/raset1/u(97),c,cd,cm,i,j
      data ij,kl/1802,9373/
c
      if(iunit.eq.0)                            then
      i=mod(ij/177, 177)+2
      j=mod(ij,     177)+2
      k=mod(kl/169, 178)+1
      m=mod(kl,     169)
      print '(a,2i7,4i4)',
     & ' marsaglia initialized (seed=0): ',ij,kl,i,j,k,m
c
      do ii=1,97
      s=0.0d00
      t=0.5d00
       do jj=1,24
       n=mod(mod(i*j,179)*k, 179)
       i=j
       j=k
       k=n
       m=mod(53*m+1, 169)
       if(mod(m*n,64).ge.32) s=s+t
       t=0.5d00*t
       enddo
      u(ii)=s
      enddo
c
      c =  362436.d00/16777216.d00
      cd= 7654321.d00/16777216.d00
      cm=16777213.d00/16777216.d00
c
      i=97
      j=33
c
						   else
      
      print*,'marsaglia initialized, seed=',iunit
      open(41,file='rmaseed.dat',
     & status='old',form='unformatted')
      rewind 41
      read(41) u,c,cd,cm,i,j
      close(41)
						   end if
c
      return
      end


      subroutine rmar(xr)          
c     pseudo random number generator 
c     proposed by marsaglia, zaman and tsang
c     implicit real*8 (a-h,o-z)
      implicit none
      real*8 xr,u,c,cd,cm,uni
      integer i,j
      common/raset1/u(97),c,cd,cm,i,j
c
      uni=u(i)-u(j)
      if(uni.lt. 0.0e00) uni=uni+1.0e00
      u(i)=uni           
      i=i-1              
      if(i.eq.0) i=97    
      j=j-1              
      if(j.eq.0) j=97   
      c=c-cd            
      if(c  .lt. 0.0e00)   c=c+cm
      uni=uni-c
      if(uni.lt. 0) uni=uni+1
      xr=uni
c
      return
      end
	
	   

      subroutine  rmaget(iunit)
c     implicit real*8 (a-h,o-z)
      implicit none
      real*8 u,c,cd,cm
      integer iunit,i,j
      common/raset1/u(97),c,cd,cm,i,j
      open(41,file='rmaseed.dat',
     & status='unknown',form='unformatted')
      write(41) u,c,cd,cm,i,j
      return
      end