Support for discrete univariate probability distribution.

src/distribution_univariate.F90

@@ -30,12 +30,14 @@ module distribution_univariate
 ! Derived classes of Distribution
 !===============================================================================
 
-  ! delta function at a single point
-  type, extends(Distribution) :: Delta
-    real(8) :: x0
+  ! Discrete distribution
+  type, extends(Distribution) :: Discrete
+    real(8), allocatable :: x(:)
+    real(8), allocatable :: p(:)
   contains
-    procedure :: sample => delta_sample
-  end type Delta
+    procedure :: sample => discrete_sample
+    procedure :: initialize => discrete_initialize
+  end type Discrete
 
   ! Uniform distribution over the interval [a,b]
   type, extends(Distribution) :: Uniform
@@ -73,12 +75,50 @@ module distribution_univariate
 
 contains
 
-  function delta_sample(this) result(x)
-    class(Delta), intent(in) :: this
+  function discrete_sample(this) result(x)
+    class(Discrete), intent(in) :: this
     real(8) :: x
 
-    x = this%x0
-  end function delta_sample
+    integer :: n  ! size of distribution
+    real(8) :: c  ! cumulative frequency
+    real(8) :: xi ! sampled CDF value
+
+    n = size(this%x)
+    if (n > 1) then
+      xi = prn()
+      c = ZERO
+      do i = 1, size(this%x)
+        c = c + this%p(i)
+        if (xi < c) exit
+      end do
+      x = this%x(i)
+    else
+      x = this%x(1)
+    end if
+  end function discrete_sample
+
+  subroutine discrete_initialize(this, x, p)
+    class(Discrete), intent(inout) :: this
+    real(8), intent(in) :: x(:)
+    real(8), intent(in) :: p(:)
+
+    integer :: n
+
+    ! Check length of x, p arrays
+    if (size(x) /= size(p)) then
+      call fatal_error('Tabulated probabilities not of same length as &
+           &independent variable.')
+    end if
+
+    ! Copy probability density function
+    n = size(x)
+    allocate(this%x(n), this%p(n))
+    this%x(:) = x(:)
+    this%p(:) = p(:)
+
+    ! Normalize density function
+    this%p(:) = this%p(:)/sum(this%p)
+  end subroutine
 
   function uniform_sample(this) result(x)
     class(Uniform), intent(in) :: this

src/input_xml.F90

@@ -529,7 +529,7 @@ contains
              call get_node_value(node_dist, "type", type)
         select case (to_lower(type))
         case ('monoenergetic')
-          allocate(Delta :: external_source%energy)
+          allocate(Discrete :: external_source%energy)
           if (n /= 1) then
             call fatal_error('Monoenergetic energy distribution must have one &
                  &parameter specified.')
@@ -559,8 +559,10 @@ contains
 
         ! Read parameters for energy distribution
         select type(energy => external_source%energy)
-        type is (Delta)
-          call get_node_value(node_dist, "parameters", energy%x0)
+        type is (Discrete)
+          allocate(energy%x(1), energy%p(1))
+          call get_node_value(node_dist, "parameters", energy%x(1))
+          energy%p(1) = ONE
 
         type is (Maxwell)
           call get_node_value(node_dist, "parameters", energy%theta)

src/source.F90

@@ -2,7 +2,7 @@ module source
 
   use bank_header,      only: Bank
   use constants
-  use distribution_univariate, only: Delta
+  use distribution_univariate, only: Discrete
   use error,            only: fatal_error
   use geometry,         only: find_cell
   use geometry_header,  only: BASE_UNIVERSE
@@ -189,8 +189,8 @@ contains
 
     ! Check for monoenergetic source above maximum neutron energy
     select type (energy => external_source%energy)
-    type is (Delta)
-      if (energy%x0 >= energy_max_neutron) then
+    type is (Discrete)
+      if (any(energy%x >= energy_max_neutron)) then
         call fatal_error("Source energy above range of energies of at least &
              &one cross section table")
       end if