Ability to sample monodirectional, monoenergetic, maxwell, watt sources.

src/input_xml.F90

@@ -170,7 +170,7 @@ contains
        end select
 
        ! Read parameters for spatial distribution
-       if (type /= 'file') then
+       if (type /= 'file' .and. associated(source_ % space % parameters)) then
           n = size(source_ % space % parameters)
           if (n < coeffs_reqd) then
              message = "Not enough parameters specified for spatial " &
@@ -198,9 +198,10 @@ contains
        select case (type)
        case ('isotropic')
           external_source % type_angle = SRC_ANGLE_ISOTROPIC
-       case ('monoenergtic')
+          coeffs_reqd = 0
+       case ('monodirectional')
           external_source % type_angle = SRC_ANGLE_MONO
-          coeffs_reqd = 1
+          coeffs_reqd = 3
        case ('tabular')
           external_source % type_angle = SRC_ANGLE_TABULAR
        case default
@@ -208,6 +209,23 @@ contains
                // trim(source_ % angle % type)
           call fatal_error()
        end select
+
+       ! Read parameters for angle distribution
+       if (associated(source_ % angle % parameters)) then
+          n = size(source_ % angle % parameters)
+          if (n < coeffs_reqd) then
+             message = "Not enough parameters specified for angle " &
+                  // "distribution of external source."
+             call fatal_error()
+          elseif (n > coeffs_reqd) then
+             message = "Too many parameters specified for angle " &
+                  // "distribution of external source."
+             call fatal_error()
+          else
+             allocate(external_source % params_angle(n))
+             external_source % params_angle = source_ % angle % parameters
+          end if
+       end if
     else
        ! Set default angular distribution isotropic
        external_source % type_angle  = SRC_ANGLE_ISOTROPIC
@@ -221,11 +239,13 @@ contains
        select case (type)
        case ('monoenergetic')
           external_source % type_energy = SRC_ENERGY_MONO
+          coeffs_reqd = 1
        case ('maxwell')
           external_source % type_energy = SRC_ENERGY_MAXWELL
           coeffs_reqd = 1
        case ('watt')
           external_source % type_energy = SRC_ENERGY_WATT
+          coeffs_reqd = 2
        case ('tabular')
           external_source % type_energy = SRC_ENERGY_TABULAR
        case default
@@ -233,6 +253,23 @@ contains
                // trim(source_ % energy % type)
           call fatal_error()
        end select
+
+       ! Read parameters for energy distribution
+       if (associated(source_ % energy % parameters)) then
+          n = size(source_ % energy % parameters)
+          if (n < coeffs_reqd) then
+             message = "Not enough parameters specified for energy " &
+                  // "distribution of external source."
+             call fatal_error()
+          elseif (n > coeffs_reqd) then
+             message = "Too many parameters specified for energy " &
+                  // "distribution of external source."
+             call fatal_error()
+          else
+             allocate(external_source % params_energy(n))
+             external_source % params_energy = source_ % energy % parameters
+          end if
+       end if
     else
        ! Set default energy distribution to Watt fission spectrum
        external_source % type_energy = SRC_ENERGY_WATT

src/source.F90

@@ -7,7 +7,7 @@ module source
   use global
   use output,          only: write_message
   use particle_header, only: deallocate_coord
-  use physics,         only: watt_spectrum
+  use physics,         only: maxwell_spectrum, watt_spectrum
   use random_lcg,      only: prn, set_particle_seed
   use string,          only: to_str
 
@@ -108,8 +108,8 @@ contains
     real(8) :: mu         ! cosine of polar angle
     real(8) :: p_min(3)   ! minimum coordinates of source
     real(8) :: p_max(3)   ! maximum coordinates of source
-    real(8) :: Watt_a     ! Watt spectrum parameter 'a'
-    real(8) :: Watt_b     ! Watt spectrum parameter 'b'
+    real(8) :: a          ! Arbitrary parameter 'a'
+    real(8) :: b          ! Arbitrary parameter 'b'
 
     ! Sample position
     select case (external_source % type_space)
@@ -125,29 +125,48 @@ contains
     ! Sample angle
     select case (external_source % type_angle)
     case (SRC_ANGLE_ISOTROPIC)
+       ! Sample isotropic distribution
        phi = TWO*PI*prn()
        mu = TWO*prn() - ONE
        site % uvw(1) = mu
        site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
        site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
+
+    case (SRC_ANGLE_MONO)
+       site % uvw = external_source % params_angle
+
     case default
        message = "No angle distribution specified for external source!"
        call fatal_error()
     end select
           
-    ! Sample energy
+    ! Sample energy distribution
     select case (external_source % type_energy)
-    case (SRC_ENERGY_WATT)
-       Watt_a = external_source % params_energy(1)
-       Watt_b = external_source % params_energy(2)
+    case (SRC_ENERGY_MONO)
+       ! Monoenergtic source
+       site % E = external_source % params_energy(1)
 
+    case (SRC_ENERGY_MAXWELL)
+       a = external_source % params_energy(1)
        do
-          ! Sample Watt fission spectrum
-          site % E = watt_spectrum(Watt_a, Watt_b)
+          ! Sample Maxwellian fission spectrum
+          site % E = maxwell_spectrum(a)
 
           ! resample if energy is >= 20 MeV
           if (site % E < 20) exit
        end do
+
+    case (SRC_ENERGY_WATT)
+       a = external_source % params_energy(1)
+       b = external_source % params_energy(2)
+       do
+          ! Sample Watt fission spectrum
+          site % E = watt_spectrum(a, b)
+
+          ! resample if energy is >= 20 MeV
+          if (site % E < 20) exit
+       end do
+
     case default
        message = "No energy distribution specified for external source!"
        call fatal_error()

src/xml-fortran/templates/settings_t.xml

@@ -11,7 +11,7 @@
   </typedef>
 
   <typedef name="distribution_xml">
-    <component name="type" type="word" length="10" />
+    <component name="type" type="word" length="16" />
     <component name="length" type="integer" />
     <component name="interpolation" type="word" length="10" />
     <component name="parameters" type="double-array" />