anchor.diff:          t=1.2*sq*(1.+y**2)*dexp(oldg-gammln(em+1.) 
        +-gammln(en-em+1.)+em*plog+(en-em)*pclog) 
   C       Reject. This happens about 1.5 times per deviate, on average.
  -        call random_number(rn)
  +!         call random_number(rn)
  +          rn = rand()
           if (rn.gt.t) goto 2 
           x=em 
           endif 
  
positive.code: C       require up to 25 calls to ran1.
          x=0. 
          do 11 j=1,n 
            call random_number(rn)
            if(rn.lt.p) x=x+1. 
     11   enddo 
        else if (am.lt.1.) then 
  
negative.code:         g=dexp(-am) 
          t=1. 
          do 12 j=0,n
  !         call random_number(rn) 
            rn = rand()
          t=t*rn
          if (t.lt.g) goto 1 
     12   enddo  
  