mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge remote-tracking branch 'walshjon/iso-lab' into iso-lab
This commit is contained in:
Will Boyd
commit
b6090e26c6
4 changed files with 328 additions and 25 deletions
|
|
@ -7,7 +7,7 @@ module input_xml
|
|||
use error, only: fatal_error, warning
|
||||
use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice
|
||||
use global
|
||||
use list_header, only: ListChar, ListReal
|
||||
use list_header, only: ListChar, ListLog, ListReal
|
||||
use mesh_header, only: StructuredMesh
|
||||
use output, only: write_message
|
||||
use plot_header
|
||||
|
|
@ -1581,6 +1581,7 @@ contains
|
|||
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
|
||||
type(ListChar) :: list_names ! temporary list of nuclide names
|
||||
type(ListReal) :: list_density ! temporary list of nuclide densities
|
||||
type(ListLog) :: list_iso_lab ! temporary list of isotropic lab scatterers
|
||||
type(Material), pointer :: mat => null()
|
||||
type(Node), pointer :: doc => null()
|
||||
type(Node), pointer :: node_mat => null()
|
||||
|
|
@ -1735,6 +1736,20 @@ contains
|
|||
end if
|
||||
end if
|
||||
|
||||
! Check enforced isotropic lab scattering
|
||||
if (check_for_node(node_nuc, "lab")) then
|
||||
call get_node_value(node_nuc, "lab", temp_str)
|
||||
if (trim(adjustl(to_lower(temp_str))) == "true") then
|
||||
call list_iso_lab % append(.true.)
|
||||
else if (trim(adjustl(to_lower(temp_str))) == "false") then
|
||||
call list_iso_lab % append(.false.)
|
||||
else
|
||||
call fatal_error("Isotropic lab scattering must be true or false")
|
||||
end if
|
||||
else
|
||||
call list_iso_lab % append(.false.)
|
||||
end if
|
||||
|
||||
! store full name
|
||||
call get_node_value(node_nuc, "name", temp_str)
|
||||
if (check_for_node(node_nuc, "xs")) &
|
||||
|
|
@ -1826,6 +1841,7 @@ contains
|
|||
allocate(mat % names(n))
|
||||
allocate(mat % nuclide(n))
|
||||
allocate(mat % atom_density(n))
|
||||
allocate(mat % p0(n))
|
||||
|
||||
ALL_NUCLIDES: do j = 1, mat % n_nuclides
|
||||
! Check that this nuclide is listed in the cross_sections.xml file
|
||||
|
|
@ -1862,6 +1878,10 @@ contains
|
|||
! Copy name and atom/weight percent
|
||||
mat % names(j) = name
|
||||
mat % atom_density(j) = list_density % get_item(j)
|
||||
|
||||
! Copy isotropic lab scattering flag
|
||||
mat % p0(j) = list_iso_lab % get_item(j)
|
||||
|
||||
end do ALL_NUCLIDES
|
||||
|
||||
! Check to make sure either all atom percents or all weight percents are
|
||||
|
|
@ -1878,6 +1898,7 @@ contains
|
|||
! Clear lists
|
||||
call list_names % clear()
|
||||
call list_density % clear()
|
||||
call list_iso_lab % clear()
|
||||
|
||||
! =======================================================================
|
||||
! READ AND PARSE <sab> TAG FOR S(a,b) DATA
|
||||
|
|
|
|||
|
|
@ -34,6 +34,13 @@ module list_header
|
|||
type(ListElemChar), pointer :: prev => null()
|
||||
end type ListElemChar
|
||||
|
||||
type :: ListElemLog
|
||||
logical :: data
|
||||
type(ListElemLog), pointer :: next => null()
|
||||
type(ListElemLog), pointer :: prev => null()
|
||||
end type ListElemLog
|
||||
|
||||
|
||||
!===============================================================================
|
||||
! LIST* types contain the linked list with convenience methods. We originally
|
||||
! considered using unlimited polymorphism to provide a single type, but compiler
|
||||
|
|
@ -107,6 +114,28 @@ module list_header
|
|||
procedure :: size => list_size_char ! Size of list
|
||||
end type ListChar
|
||||
|
||||
type, public :: ListLog
|
||||
private
|
||||
integer :: count = 0 ! Number of elements in list
|
||||
|
||||
! Used in get_item for fast sequential lookups
|
||||
integer :: last_index = huge(0)
|
||||
type(ListElemLog), pointer :: last_elem => null()
|
||||
|
||||
! Pointers to beginning and end of list
|
||||
type(ListElemLog), public, pointer :: head => null()
|
||||
type(ListElemLog), public, pointer :: tail => null()
|
||||
contains
|
||||
procedure :: append => list_append_log ! Add item to end of list
|
||||
procedure :: clear => list_clear_log ! Remove all items
|
||||
procedure :: contains => list_contains_log ! Does list contain?
|
||||
procedure :: get_item => list_get_item_log ! Get i-th item in list
|
||||
procedure :: index => list_index_log ! Determine index of given item
|
||||
procedure :: insert => list_insert_log ! Insert item in i-th position
|
||||
procedure :: remove => list_remove_log ! Remove specified item
|
||||
procedure :: size => list_size_log ! Size of list
|
||||
end type ListLog
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -189,6 +218,31 @@ contains
|
|||
|
||||
end subroutine list_append_char
|
||||
|
||||
subroutine list_append_log(this, data)
|
||||
class(ListLog) :: this
|
||||
logical :: data
|
||||
|
||||
type(ListElemLog), pointer :: elem
|
||||
|
||||
! Create element and set dat
|
||||
allocate(elem)
|
||||
elem % data = data
|
||||
|
||||
if (.not. associated(this % head)) then
|
||||
! If list is empty, set head and tail to new element
|
||||
this % head => elem
|
||||
this % tail => elem
|
||||
else
|
||||
! Otherwise append element at end of list
|
||||
this % tail % next => elem
|
||||
elem % prev => this % tail
|
||||
this % tail => this % tail % next
|
||||
end if
|
||||
|
||||
this % count = this % count + 1
|
||||
|
||||
end subroutine list_append_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_CLEAR removes all elements from the list
|
||||
!===============================================================================
|
||||
|
|
@ -271,6 +325,32 @@ contains
|
|||
|
||||
end subroutine list_clear_char
|
||||
|
||||
subroutine list_clear_log(this)
|
||||
class(ListLog) :: this
|
||||
|
||||
type(ListElemLog), pointer :: current => null()
|
||||
type(ListElemLog), pointer :: next => null()
|
||||
|
||||
if (this % count > 0) then
|
||||
current => this % head
|
||||
do while (associated(current))
|
||||
! Set pointer to next element
|
||||
next => current % next
|
||||
|
||||
! Deallocate memory for current element
|
||||
deallocate(current)
|
||||
|
||||
! Move to next element
|
||||
current => next
|
||||
end do
|
||||
|
||||
nullify(this % head)
|
||||
nullify(this % tail)
|
||||
this % count = 0
|
||||
end if
|
||||
|
||||
end subroutine list_clear_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_CONTAINS determines whether the list contains a specified item. Since it
|
||||
! relies on the index method, it is O(n).
|
||||
|
|
@ -303,6 +383,15 @@ contains
|
|||
|
||||
end function list_contains_char
|
||||
|
||||
function list_contains_log(this, data) result(in_list)
|
||||
class(ListLog) :: this
|
||||
logical :: data
|
||||
logical :: in_list
|
||||
|
||||
in_list = (this % index(data) > 0)
|
||||
|
||||
end function list_contains_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_GET_ITEM returns the item in the list at position 'i_list'. If the index
|
||||
! is out of bounds, an error code is returned.
|
||||
|
|
@ -407,6 +496,39 @@ contains
|
|||
|
||||
end function list_get_item_char
|
||||
|
||||
function list_get_item_log(this, i_list) result(data)
|
||||
class(ListLog) :: this
|
||||
integer :: i_list
|
||||
logical :: data
|
||||
|
||||
integer :: last_index
|
||||
|
||||
if (i_list < 1 .or. i_list > this % count) then
|
||||
! Check for index out of bounds
|
||||
data = .false.
|
||||
elseif (i_list == 1) then
|
||||
data = this % head % data
|
||||
this % last_index = 1
|
||||
this % last_elem => this % head
|
||||
elseif (i_list == this % count) then
|
||||
data = this % tail % data
|
||||
this % last_index = this % count
|
||||
this % last_elem => this % tail
|
||||
else
|
||||
if (i_list < this % last_index) then
|
||||
this % last_index = 1
|
||||
this % last_elem => this % head
|
||||
end if
|
||||
|
||||
do last_index = this % last_index + 1, i_list
|
||||
this % last_elem => this % last_elem % next
|
||||
this % last_index = last_index
|
||||
end do
|
||||
data = this % last_elem % data
|
||||
end if
|
||||
|
||||
end function list_get_item_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_INDEX determines the first index in the list that contains 'data'. If
|
||||
! 'data' is not present in the list, the return value is -1.
|
||||
|
|
@ -475,6 +597,27 @@ contains
|
|||
|
||||
end function list_index_char
|
||||
|
||||
function list_index_log(this, data) result(i_list)
|
||||
|
||||
class(ListLog) :: this
|
||||
logical :: data
|
||||
integer :: i_list
|
||||
|
||||
type(ListElemLog), pointer :: elem
|
||||
|
||||
i_list = 0
|
||||
elem => this % head
|
||||
do while (associated(elem))
|
||||
i_list = i_list + 1
|
||||
if (data .eqv. elem % data) exit
|
||||
elem => elem % next
|
||||
end do
|
||||
|
||||
! Check if we reached the end of the list
|
||||
if (.not. associated(elem)) i_list = -1
|
||||
|
||||
end function list_index_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_INSERT inserts 'data' at index 'i_list' within the list. If 'i_list'
|
||||
! exceeds the size of the list, the data is appends at the end of the list.
|
||||
|
|
@ -646,6 +789,62 @@ contains
|
|||
|
||||
end subroutine list_insert_char
|
||||
|
||||
subroutine list_insert_log(this, i_list, data)
|
||||
|
||||
class(ListLog) :: this
|
||||
integer :: i_list
|
||||
logical :: data
|
||||
|
||||
integer :: i
|
||||
type(ListElemLog), pointer :: elem => null()
|
||||
type(ListElemLog), pointer :: new_elem => null()
|
||||
|
||||
if (i_list > this % count) then
|
||||
! Check whether specified index is greater than number of elements -- if
|
||||
! so, just append it to the end of the list
|
||||
call this % append(data)
|
||||
|
||||
else if (i_list == 1) then
|
||||
! Check for new head element
|
||||
allocate(new_elem)
|
||||
new_elem % data = data
|
||||
new_elem % next => this % head
|
||||
this % head => new_elem
|
||||
this % count = this % count + 1
|
||||
|
||||
else
|
||||
! Default case with new element somewhere in middle of list
|
||||
if (i_list >= this % last_index) then
|
||||
i = this % last_index
|
||||
elem => this % last_elem
|
||||
else
|
||||
i = 0
|
||||
elem => this % head
|
||||
end if
|
||||
do while (associated(elem))
|
||||
i = i + 1
|
||||
if (i == i_list - 1) then
|
||||
! Allocate new element
|
||||
allocate(new_elem)
|
||||
new_elem % data = data
|
||||
|
||||
! Put it before the i-th element
|
||||
new_elem % prev => elem % prev
|
||||
new_elem % next => elem
|
||||
new_elem % prev % next => new_elem
|
||||
new_elem % next % prev => new_elem
|
||||
this % count = this % count + 1
|
||||
this % last_index = i_list
|
||||
this % last_elem => new_elem
|
||||
exit
|
||||
end if
|
||||
i = i + 1
|
||||
elem => elem % next
|
||||
end do
|
||||
end if
|
||||
|
||||
end subroutine list_insert_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_REMOVE removes the first item in the list that contains 'data'. If 'data'
|
||||
! is not in the list, no action is taken.
|
||||
|
|
@ -768,6 +967,45 @@ contains
|
|||
|
||||
end subroutine list_remove_char
|
||||
|
||||
subroutine list_remove_log(this, data)
|
||||
|
||||
class(ListLog) :: this
|
||||
logical :: data
|
||||
|
||||
type(ListElemLog), pointer :: elem => null()
|
||||
|
||||
elem => this % head
|
||||
do while (associated(elem))
|
||||
! Check for matching data
|
||||
if (elem % data .eqv. data) then
|
||||
|
||||
! Determine whether the current element is the head, tail, or a middle
|
||||
! element
|
||||
if (associated(elem, this % head)) then
|
||||
this % head => elem % next
|
||||
if (associated(elem, this % tail)) nullify(this % tail)
|
||||
if (associated(this % head)) nullify(this % head % prev)
|
||||
deallocate(elem)
|
||||
else if (associated(elem, this % tail)) then
|
||||
this % tail => elem % prev
|
||||
deallocate(this % tail % next)
|
||||
else
|
||||
elem % prev % next => elem % next
|
||||
elem % next % prev => elem % prev
|
||||
deallocate(elem)
|
||||
end if
|
||||
|
||||
! Decrease count and exit
|
||||
this % count = this % count - 1
|
||||
exit
|
||||
end if
|
||||
|
||||
! Advance pointers
|
||||
elem => elem % next
|
||||
end do
|
||||
|
||||
end subroutine list_remove_log
|
||||
|
||||
!===============================================================================
|
||||
! LIST_SIZE returns the number of elements in the list
|
||||
!===============================================================================
|
||||
|
|
@ -799,4 +1037,13 @@ contains
|
|||
|
||||
end function list_size_char
|
||||
|
||||
function list_size_log(this) result(size)
|
||||
|
||||
class(ListLog) :: this
|
||||
integer :: size
|
||||
|
||||
size = this % count
|
||||
|
||||
end function list_size_log
|
||||
|
||||
end module list_header
|
||||
|
|
|
|||
|
|
@ -32,6 +32,9 @@ module material_header
|
|||
! Does this material contain fissionable nuclides?
|
||||
logical :: fissionable = .false.
|
||||
|
||||
! enforce isotropic scattering in lab
|
||||
logical, allocatable :: p0(:)
|
||||
|
||||
end type Material
|
||||
|
||||
end module material_header
|
||||
|
|
|
|||
|
|
@ -73,11 +73,12 @@ contains
|
|||
type(Particle), intent(inout) :: p
|
||||
|
||||
integer :: i_nuclide ! index in nuclides array
|
||||
integer :: i_nuc_mat ! index in material's nuclides array
|
||||
integer :: i_reaction ! index in nuc % reactions array
|
||||
type(Nuclide), pointer, save :: nuc => null()
|
||||
!$omp threadprivate(nuc)
|
||||
|
||||
i_nuclide = sample_nuclide(p, 'total ')
|
||||
call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat)
|
||||
|
||||
! Get pointer to table
|
||||
nuc => nuclides(i_nuclide)
|
||||
|
|
@ -107,7 +108,7 @@ contains
|
|||
|
||||
! Sample a scattering reaction and determine the secondary energy of the
|
||||
! exiting neutron
|
||||
call scatter(p, i_nuclide)
|
||||
call scatter(p, i_nuclide, i_nuc_mat)
|
||||
|
||||
! Play russian roulette if survival biasing is turned on
|
||||
|
||||
|
|
@ -122,13 +123,13 @@ contains
|
|||
! SAMPLE_NUCLIDE
|
||||
!===============================================================================
|
||||
|
||||
function sample_nuclide(p, base) result(i_nuclide)
|
||||
subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat)
|
||||
|
||||
type(Particle), intent(in) :: p
|
||||
character(7), intent(in) :: base ! which reaction to sample based on
|
||||
integer :: i_nuclide
|
||||
integer, intent(out) :: i_nuclide
|
||||
integer, intent(out) :: i_nuc_mat
|
||||
|
||||
integer :: i
|
||||
real(8) :: prob
|
||||
real(8) :: cutoff
|
||||
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
|
||||
|
|
@ -149,20 +150,20 @@ contains
|
|||
cutoff = prn() * material_xs % fission
|
||||
end select
|
||||
|
||||
i = 0
|
||||
i_nuc_mat = 0
|
||||
prob = ZERO
|
||||
do while (prob < cutoff)
|
||||
i = i + 1
|
||||
i_nuc_mat = i_nuc_mat + 1
|
||||
|
||||
! Check to make sure that a nuclide was sampled
|
||||
if (i > mat % n_nuclides) then
|
||||
if (i_nuc_mat > mat % n_nuclides) then
|
||||
call write_particle_restart(p)
|
||||
call fatal_error("Did not sample any nuclide during collision.")
|
||||
end if
|
||||
|
||||
! Find atom density
|
||||
i_nuclide = mat % nuclide(i)
|
||||
atom_density = mat % atom_density(i)
|
||||
i_nuclide = mat % nuclide(i_nuc_mat)
|
||||
atom_density = mat % atom_density(i_nuc_mat)
|
||||
|
||||
! Determine microscopic cross section
|
||||
select case (base)
|
||||
|
|
@ -179,7 +180,7 @@ contains
|
|||
prob = prob + sigma
|
||||
end do
|
||||
|
||||
end function sample_nuclide
|
||||
end subroutine sample_nuclide
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_FISSION
|
||||
|
|
@ -303,10 +304,11 @@ contains
|
|||
! SCATTER
|
||||
!===============================================================================
|
||||
|
||||
subroutine scatter(p, i_nuclide)
|
||||
subroutine scatter(p, i_nuclide, i_nuc_mat)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
integer, intent(in) :: i_nuclide
|
||||
integer, intent(in) :: i_nuc_mat
|
||||
|
||||
integer :: i
|
||||
integer :: i_grid
|
||||
|
|
@ -334,6 +336,10 @@ contains
|
|||
! ELASTIC SCATTERING
|
||||
|
||||
if (micro_xs(i_nuclide) % index_sab /= NONE) then
|
||||
if (materials(p % material) % p0(i_nuc_mat)) &
|
||||
& call fatal_error("thermal scattering law data and isotropic lab&
|
||||
& scattering specified for the same nuclide")
|
||||
|
||||
! S(a,b) scattering
|
||||
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, &
|
||||
p % E, p % coord0 % uvw, p % mu)
|
||||
|
|
@ -344,7 +350,8 @@ contains
|
|||
|
||||
! Perform collision physics for elastic scattering
|
||||
call elastic_scatter(i_nuclide, rxn, &
|
||||
p % E, p % coord0 % uvw, p % mu, p % wgt)
|
||||
p % E, p % coord0 % uvw, p % mu, p % wgt, &
|
||||
& materials(p % material) % p0(i_nuc_mat))
|
||||
end if
|
||||
|
||||
p % event_MT = ELASTIC
|
||||
|
|
@ -387,7 +394,7 @@ contains
|
|||
|
||||
! Perform collision physics for inelastic scattering
|
||||
call inelastic_scatter(nuc, rxn, p % E, p % coord0 % uvw, &
|
||||
p % mu, p % wgt)
|
||||
p % mu, p % wgt, materials(p % material) % p0(i_nuc_mat))
|
||||
p % event_MT = rxn % MT
|
||||
|
||||
end if
|
||||
|
|
@ -402,21 +409,24 @@ contains
|
|||
! target.
|
||||
!===============================================================================
|
||||
|
||||
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt)
|
||||
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt, iso_lab)
|
||||
|
||||
integer, intent(in) :: i_nuclide
|
||||
type(Reaction), pointer :: rxn
|
||||
real(8), intent(inout) :: E
|
||||
real(8), intent(inout) :: uvw(3)
|
||||
real(8), intent(out) :: mu_lab
|
||||
real(8), intent(inout) :: wgt
|
||||
real(8), intent(out) :: mu_lab ! cosine of polar angle in lab system
|
||||
logical, intent(in) :: iso_lab
|
||||
|
||||
real(8) :: awr ! atomic weight ratio of target
|
||||
real(8) :: mu_cm ! cosine of polar angle in center-of-mass
|
||||
real(8) :: phi ! azimuthal angle
|
||||
real(8) :: vel ! magnitude of velocity
|
||||
real(8) :: v_n(3) ! velocity of neutron
|
||||
real(8) :: v_cm(3) ! velocity of center-of-mass
|
||||
real(8) :: v_t(3) ! velocity of target nucleus
|
||||
real(8) :: uvw_in(3) ! incoming direction
|
||||
real(8) :: uvw_cm(3) ! directional cosines in center-of-mass
|
||||
type(Nuclide), pointer, save :: nuc => null()
|
||||
!$omp threadprivate(nuc)
|
||||
|
|
@ -430,6 +440,9 @@ contains
|
|||
! Neutron velocity in LAB
|
||||
v_n = vel * uvw
|
||||
|
||||
! incoming direction
|
||||
uvw_in(:) = uvw(:)
|
||||
|
||||
! Sample velocity of target nucleus
|
||||
if (.not. micro_xs(i_nuclide) % use_ptable) then
|
||||
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
|
||||
|
|
@ -465,11 +478,17 @@ contains
|
|||
vel = sqrt(E)
|
||||
|
||||
! compute cosine of scattering angle in LAB frame by taking dot product of
|
||||
! neutron's pre- and post-collision angle
|
||||
mu_lab = dot_product(uvw, v_n) / vel
|
||||
! neutron's pre- and post-collision unit vectors
|
||||
if (iso_lab) then
|
||||
uvw(1) = TWO * prn() - ONE
|
||||
phi = TWO * PI * prn()
|
||||
uvw(2) = cos(phi) * sqrt(ONE - uvw(1)*uvw(1))
|
||||
uvw(3) = sin(phi) * sqrt(ONE - uvw(1)*uvw(1))
|
||||
else
|
||||
uvw = v_n / vel
|
||||
end if
|
||||
|
||||
! Set energy and direction of particle in LAB frame
|
||||
uvw = v_n / vel
|
||||
mu_lab = dot_product(uvw_in, uvw)
|
||||
|
||||
end subroutine elastic_scatter
|
||||
|
||||
|
|
@ -1277,7 +1296,7 @@ contains
|
|||
! than fission), i.e. level scattering, (n,np), (n,na), etc.
|
||||
!===============================================================================
|
||||
|
||||
subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt)
|
||||
subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt, iso_lab)
|
||||
|
||||
type(Nuclide), pointer :: nuc
|
||||
type(Reaction), pointer :: rxn
|
||||
|
|
@ -1285,6 +1304,7 @@ contains
|
|||
real(8), intent(inout) :: uvw(3) ! directional cosines
|
||||
real(8), intent(out) :: mu ! cosine of scattering angle in lab
|
||||
real(8), intent(inout) :: wgt ! particle weight
|
||||
logical, intent(in) :: iso_lab
|
||||
|
||||
integer :: law ! secondary energy distribution law
|
||||
real(8) :: A ! atomic weight ratio of nuclide
|
||||
|
|
@ -1292,9 +1312,12 @@ contains
|
|||
real(8) :: E_cm ! outgoing energy in center-of-mass
|
||||
real(8) :: Q ! Q-value of reaction
|
||||
real(8) :: yield ! neutron yield
|
||||
real(8) :: uvw_in(3) ! incoming direction
|
||||
real(8) :: phi ! azimuthal angle
|
||||
|
||||
! copy energy of neutron
|
||||
! copy energy, direction of neutron
|
||||
E_in = E
|
||||
uvw_in(:) = uvw(:)
|
||||
|
||||
! determine A and Q
|
||||
A = nuc % awr
|
||||
|
|
@ -1328,8 +1351,17 @@ contains
|
|||
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
|
||||
end if
|
||||
|
||||
! change direction of particle
|
||||
uvw = rotate_angle(uvw, mu)
|
||||
! compute cosine of scattering angle in LAB frame by taking dot product of
|
||||
! neutron's pre- and post-collision unit vectors
|
||||
if (iso_lab) then
|
||||
uvw(1) = TWO * prn() - ONE
|
||||
phi = TWO * PI * prn()
|
||||
uvw(2) = cos(phi) * sqrt(ONE - uvw(1)*uvw(1))
|
||||
uvw(3) = sin(phi) * sqrt(ONE - uvw(1)*uvw(1))
|
||||
mu = dot_product(uvw_in, uvw)
|
||||
else
|
||||
uvw = rotate_angle(uvw_in, mu)
|
||||
end if
|
||||
|
||||
! change weight of particle based on yield
|
||||
if (rxn % multiplicity_with_E) then
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue