anchor.diff:        DOUBLE PRECISION X, Y, S, ONE, TWO
         DATA ONE /1.0/, TWO /2.0/
   C
  -    1 call random_number(x)
  -      call random_number(y)
  +!     1 call random_number(x)
  +    1 X = rand()
  +!       call random_number(y)
  +      Y = rand()
         X = TWO * X - ONE
         Y = TWO * Y - ONE
         S = X * X + Y * Y
  
positive.code:         sq=sqrt(2.*am*pc) 
  C       The following code should by now seem familiar: rejection 
  C       method with a Lorentzian comparison function.
          call random_number(rn)
      2   y=tan(PI*rn) 
          em=sq*y+am 
  C       Reject.
  
negative.code:         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)
            rn = rand()
          if (rn.gt.t) goto 2 
          x=em 
          endif 
  