anchor.diff:         ! sample outgoing energy based on evaporation spectrum probability
          ! density function
          do
  -          xi1 = rang()
  -          xi2 = rang()
  -          E_out = -T * log(xi1*xi2)
  +          r1 = rang()
  +          r2 = rang()
  +          E_out = -T * log(r1*r2)
             if (E_out <= E_in - U) exit
          end do
          
  
positive.code:     ! sample outgoing energy
      if (law == 44 .or. law == 61) then
         call sample_energy(rxn%edist, E_in, E, mu)
      else
         call sample_energy(rxn%edist, E_in, E)
      end if
  
negative.code:        do i = 1, n/2
            x = x * rang()
         end do
         x = -2./n * log(x)
  
      case (1)
         ! Odd number of degrees of freedom can be sampled via rule C64. We can
  