c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      REAL*4 FUNCTION GASDEV(ISEED)
c
c Standard Normal
c
      SAVE
      INTEGER ISEED
C      REAL*4 RAN2
      DATA ISET/0/
      IF (ISET.EQ.0) THEN
1       V1=2.*RAN2(ISEED)-1.
        V2=2.*RAN2(ISEED)-1.
        R=V1**2+V2**2
        IF(R.GE.1.)GO TO 1
        FAC=SQRT(-2.*LOG(R)/R)
        GSET=V1*FAC
        GASDEV=V2*FAC
        ISET=1
      ELSE
        GASDEV=GSET
        ISET=0
      ENDIF
      RETURN
      END

c#######################################################################
C#######################################################################
C
c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      real*4 function rgam(a,b,iseed)
c     generates from the Gamma(a,b) distribution
      integer iseed
      real*4 a,b
      real*4 gamdev
      rgam=gamdev(a,iseed)/b
      return
      end
c +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
        real*4 function gamdev(a,ix)
c
c       Generates a random gamma variable with shape a>0 and scale=1
c       ix : random seed
c       requires uniform random generator, ran2
c
        integer ix
        real*4 aa,ea,u0,u1,c1,c2,c3,c4,c5,w,ran2

        real*8 a
        aa=sngl(a)
        if ( aa .eq. 1.0 ) then
                gamdev = -alog( ran2( ix ) )
                return
        endif
        if ( aa .lt. 1.0 ) then
                ea = 2.7182818 / ( aa + 2.7182818 )
 11             u0 = ran2( ix )
                u1 = ran2( ix )
                if ( u0 .le. ea ) then   
                        gamdev = ( u0 / ea ) ** (1.0/aa)
                        if ( u1 .gt. exp( -gamdev ) ) then
                                goto 11
                        else
                                return
                        endif
                else
                        gamdev = -alog( (1.0-u0)/(ea*aa) )
                        if( u1 .gt. gamdev**(aa-1.0) ) then   
                                go to 11
                        else
                                return
                        endif
                endif
        else
                c1 = aa - 1.0
                c2 = ( aa- 1.0/(6.0*aa) ) / c1
                c3 = 2.0 / c1
                c4 = c3 + 2.0
                c5 = 1.0 / sqrt(aa)
 12             u1 = ran2( ix ) 
                u2 = ran2( ix )
                u1 = u2 + c5 * ( 1.0-1.86 * u1 )
                if ( abs( u1- 0.5 ) .ge. 0.5 ) goto 12
                w = c2 * u2 / u1
          if ( ( c3*alog( u1 ) - alog( w ) + w) .ge. 1.0 ) then 
                goto 12
          else      
                gamdev = c1 * w
          endif
        endif  
        return
        end
C###################################################################################
C##########################################################################
C
c
      real*8 function resnorm(a,b,iseed)
c samples standard normal truncated to (a,b) with a<b
c for large |a|, uses MW's use of Mill's + Taylor series 
      integer iseed
      real*8 a,b,y,aa,bb,sg
      real*4 ran2, invnor, pnorm, x, z
c      
      sg=a*b
      IF (sg.LE.0) THEN 
         x=pnorm(a)
         y=dble(invnor(x+ran2(iseed)*(pnorm(b)-x)))
      ENDIF   

c      
      sg=a*b
      IF (sg.LE.0) THEN 
         x=pnorm(a)
         y=dble(invnor(x+ran2(iseed)*(pnorm(b)-x)))
      ENDIF   
c
      IF (sg.GT.0) THEN
         aa=a
         bb=b
         if (a.LT.0) then 
            aa=dabs(a)
            bb=dabs(b)
         endif     
         IF (aa.LT.6.0) THEN 
            if (bb.LT.6.0) then 
               z=pnorm(bb)
            else 
               z=1.0
            endif   
            x=pnorm(aa)
            y=dble(invnor(x+ran2(iseed)*(z-x)))
         ELSE 
            y=3.D0/(dsqrt(12.D0+aa*aa)+5.D0*aa)+
     +           1.D0/(dsqrt( 4.D0+aa*aa)+aa)
            y=aa+dble(ran2(iseed))*y     
            if (y.GT.bb) then 
               y=bb
            endif
         ENDIF
         if (a.LT.0) then      
            y=-y
         endif
      ENDIF  
c
      if (y.LT.a) then
         y=a
      else if (y.GT.b) then
         y=b
      endif
      resnorm=y     
c
      return
      end
c
C##############################################################################
C#############################################################################
C
      real function invnor(pr)
c
c Evaluates inverse normal CDF 
c
      real pr
      real p,eta,term,f1,f2,f3,f4,f5,f6
      data f1/.010328/
      data f2/.802853/
      data f3/2.515517/
      data f4/.001308/
      data f5/.189269/
      data f6/1.432788/
      if (pr.le.0..or.pr.ge.1.) then 
c         write(*,*) 'probability out of range (0,1)'
         continue
      else   
         p = pr
         if (p .gt. 0.5) then 
            p = 1.0 - pr
         endif   
         eta=sqrt(-2.*alog(p))
         term=((f1*eta+f2)*eta+f3)/(((f4*eta+f5)*eta+f6)*eta+1.0)
         if (pr.le..5) then 
            p = ( term - eta )
         else
            p = ( eta - term )
         endif   
      endif   
      invnor=p
      return
      end
c
C#######################################################################
C
        real*4 function pnorm(dx)
c
c       cdf for the standard Gaussian distribution
c       Using Kennedy & Gentle sec 5.3 pg 93
c
       
        parameter(zero = 0.00, one = 1.00, half = 0.5e0)
        parameter(two = 2.0e0, ohalf = 1.5e0, three = 3.0e0)
        parameter(ten = 10.0e0)
        parameter(four = 4.0e0, fty2 = 52.0e0, thr45 = 345.0e0)
        parameter(ety5 = 85.0e0, two55 = 255.0e0)
        parameter(rt2pi= 2.506628275e0)
        real*8 dx
        real*4 x
        x=sngl(dx)       
        dax=dabs(dx)
        phi=exp(-x*x/two)/rt2pi
        if ((x .gt.-ohalf) .and. (x .lt. three)) then
          is=int(ten + ten*abs(x))
          w=x
          do 10 ir=is,1,-1
            w=w*x*x/(two*ir+one) + x
10        continue
          pnorm=phi*w + half
        else if (x .lt. zero) then
          s=four - fty2/x + thr45/(x*x)
          is=iabs(int(s))
          w=x*x
          do 20 ir=is,1,-1
            r=dble(ir)
            r2=two*r
            w=x*x + four*r - one - r2*(r2+one)/w
20        continue
          pnorm=phi*(one/w - one)/x
        else
          s=-ten + ety5/x + two55/(x*x)
          is=int(abs(s))
          w=x*x
          do 30 ir=is,1,-1
            r=dble(ir)
            r2=two*r
            w=x*x + four*r - one - r2*(r2+one)/w
30        continue
          pnorm=one + phi*(one/w - one)/x
        endif
c
        return
        end
C
C############################################################################
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      REAL*4 FUNCTION RAN2(ISEED)
c
c U(0,1) generator --- ran2 from 2nd Edn of Numerical Recipes
c
      integer iseed,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
      real*4 am,eps,rnmx
      parameter(im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1,
     &  ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,
     &  ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.-eps)
      integer idum2,j,k,iv(ntab),iy
c
      SAVE IV,IY,IDUM2
c
      data idum2/123456789/,iv/ntab*0/,iy/0/
      if (iseed.le.0) then
         iseed=max(-iseed,1)
         idum2=iseed
         do j=ntab+8,1,-1
            k=iseed/iq1
            iseed=ia1*(iseed-k*iq1)-k*ir1
            if (iseed.lt.0) iseed=iseed+im1
            if (j.le.ntab) iv(j)=iseed
         enddo
         iy=iv(1)
      endif
      k=iseed/iq1
      iseed=ia1*(iseed-k*iq1)-k*ir1
      if (iseed.lt.0) iseed=iseed+im1
      k=idum2/iq2
      idum2=ia2*(idum2-k*iq2)-k*ir2
      if (idum2.lt.0) idum2=idum2+im2
      j=1+iy/ndiv
      iy=iv(j)-idum2
      iv(j)=iseed
      if (iy.lt.1) iy=iy+imm1
      ran2=min(am*iy,rnmx)

      RETURN
      END
c
c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
