anchor.diff:      case (66)
          ! =======================================================================
          ! N-BODY PHASE SPACE DISTRIBUTION
  +
  +       ! read number of bodies in phase space and total mass ratio
  +       n_bodies = edist % data(1)
  +       Ap       = edist % data(2)
  +
  +       ! determine E_max parameter
  +       E_max = (Ap - ONE)/Ap * (A/(A+ONE)*E_in + Q)
  +
  +       ! x is essentially a Maxwellian distribution
  +       x = maxwell_spectrum(ONE)
  +
  +       select case (n_bodies)
  +       case (3)
  +          y = maxwell_spectrum(ONE)
  +       case (4)
  +          r1 = rang()
  +          r2 = rang()
  +          r3 = rang()
  +          y = -log(r1*r2*r3)
  +       case (5)
  +          r1 = rang()
  +          r2 = rang()
  +          r3 = rang()
  +          r4 = rang()
  +          r5 = rang()
  +          r6 = rang()
  +          y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/2.*r6)**2
  +       end select
  +
  +       ! now determine v and E_out
  +       v = x/(x+y)
  +       E_out = E_max * v
   
       case (67)
          ! =======================================================================
  
positive.code: ! SAMPLE_ENERGY
  !===============================================================================
  
    subroutine sample_energy(edist, E_in, E_out, mu_out)
  
      type(AceDistEnergy), intent(inout):: edist
      real(8), intent(in)               :: E_in
      real(8), intent(out)              :: E_out
      real(8), intent(inout), optional  :: mu_out
  
      integer :: i           ! index on incoming energy grid
      integer :: k           ! sampled index on outgoing grid
  
negative.code:        end if
  
         ! Sample between the ith and (i+1)th bin
         r2 = rang()
         if (r > r2) then
            l = i + 1
         else
            l = i
  