mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Added cubic spline interpolation/integration; other fixes
This commit is contained in:
amandalund
parent
a2b5cfb2d7
commit
62985f6198
5 changed files with 225 additions and 56 deletions
|
|
@ -11,7 +11,7 @@ import pandas as pd
|
|||
from openmc.mixin import EqualityMixin
|
||||
import openmc.checkvalue as cv
|
||||
from . import HDF5_VERSION
|
||||
from .data import ATOMIC_SYMBOL
|
||||
from .data import ATOMIC_SYMBOL, EV_PER_MEV
|
||||
from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
|
||||
from .function import Tabulated1D
|
||||
|
||||
|
|
@ -297,18 +297,18 @@ class IncidentPhoton(EqualityMixin):
|
|||
:math:`p_z` for each subshell). Note that subshell occupancies may not
|
||||
match the atomic relaxation data.
|
||||
stopping_powers : dict
|
||||
Dictionary of stopping power data with keys 'energy', 'density' (mass
|
||||
density in g/cm:sup:`3`), 'I' (mean excitation energy), 's_collision'
|
||||
(collision stopping power in MeV cm:sup:`2`/g), 's_radiative'
|
||||
(radiative stopping power in MeV cm:sup:`2`/g), and 'density_effect'
|
||||
(density effect parameter).
|
||||
Dictionary of stopping power data with keys 'energy' (in eV), 'density'
|
||||
(mass density in g/cm:sup:`3`), 'I' (mean excitation energy),
|
||||
's_collision' (collision stopping power in eV cm:sup:`2`/g),
|
||||
's_radiative' (radiative stopping power in eV cm:sup:`2`/g), and
|
||||
'density_effect' (density effect parameter).
|
||||
bremsstrahlung : dict
|
||||
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
|
||||
(incident electron kinetic energy values in MeV), 'photon_energy'
|
||||
(incident electron kinetic energy values in eV), 'photon_energy'
|
||||
(ratio of the energy of the emitted photon to the incident electron
|
||||
kinetic energy), and 'dcs' (cross sectin values). The cross sections
|
||||
are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`, where
|
||||
:math:`E_k` is the energy of the emitted photon.
|
||||
kinetic energy), and 'dcs' (cross sectin values in mb). The cross
|
||||
sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`,
|
||||
where :math:`E_k` is the energy of the emitted photon.
|
||||
reactions : collections.OrderedDict
|
||||
Contains the cross sections for each photon reaction. The keys are MT
|
||||
values and the values are instances of :class:`PhotonReaction`.
|
||||
|
|
@ -429,7 +429,8 @@ class IncidentPhoton(EqualityMixin):
|
|||
if not _STOPPING_POWERS:
|
||||
filename = os.path.join(os.path.dirname(__file__), 'stopping_powers.h5')
|
||||
with h5py.File(filename, 'r') as f:
|
||||
_STOPPING_POWERS['energy'] = f['energy'].value
|
||||
# Units are in MeV; convert to eV
|
||||
_STOPPING_POWERS['energy'] = f['energy'].value*EV_PER_MEV
|
||||
for i in range(1, 99):
|
||||
group = f['{:03}'.format(i)]
|
||||
_STOPPING_POWERS[i] = {'density': group.attrs['density'],
|
||||
|
|
@ -438,6 +439,10 @@ class IncidentPhoton(EqualityMixin):
|
|||
's_radiative': group['s_radiative'].value,
|
||||
'density_effect': group['density_effect'].value}
|
||||
|
||||
# Units are in MeV cm^2/g; convert to eV cm^2/g
|
||||
_STOPPING_POWERS[i]['s_collision'] *= EV_PER_MEV
|
||||
_STOPPING_POWERS[i]['s_radiative'] *= EV_PER_MEV
|
||||
|
||||
# Add stopping power data
|
||||
if Z < 99:
|
||||
data.stopping_powers['energy'] = _STOPPING_POWERS['energy']
|
||||
|
|
@ -448,8 +453,8 @@ class IncidentPhoton(EqualityMixin):
|
|||
filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT')
|
||||
brem = open(filename, 'r').read().split()
|
||||
|
||||
# Incident electron kinetic energy grid
|
||||
_BREMSSTRAHLUNG['electron_energy'] = np.logspace(-3, 3, 200)
|
||||
# Incident electron kinetic energy grid in eV
|
||||
_BREMSSTRAHLUNG['electron_energy'] = np.logspace(3, 9, 200)
|
||||
log_energy = np.log(_BREMSSTRAHLUNG['electron_energy'])
|
||||
|
||||
# Get number of tabulated electron and photon energy values
|
||||
|
|
@ -460,8 +465,8 @@ class IncidentPhoton(EqualityMixin):
|
|||
p = 39
|
||||
|
||||
# Get log of incident electron kinetic energy values, used for cubic
|
||||
# spline interpolation in log energy
|
||||
logx = np.log(np.fromiter(brem[p:p+n], float, n))
|
||||
# spline interpolation in log energy. Units are in MeV, so convert to eV.
|
||||
logx = np.log(np.fromiter(brem[p:p+n], float, n)*EV_PER_MEV)
|
||||
p += n
|
||||
|
||||
# Get reduced photon energy values
|
||||
|
|
|
|||
145
src/math.F90
145
src/math.F90
|
|
@ -2,6 +2,7 @@ module math
|
|||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use constants
|
||||
use random_lcg, only: prn
|
||||
|
||||
|
|
@ -827,4 +828,148 @@ contains
|
|||
end do
|
||||
end subroutine broaden_wmp_polynomials
|
||||
|
||||
!===============================================================================
|
||||
! SPLINE
|
||||
!===============================================================================
|
||||
|
||||
subroutine spline(x, y, z, n)
|
||||
|
||||
integer, intent(in) :: n
|
||||
real(8), intent(in) :: x(n)
|
||||
real(8), intent(in) :: y(n)
|
||||
real(8), intent(out) :: z(n)
|
||||
|
||||
integer :: i
|
||||
real(8) :: a, b, c, d
|
||||
real(8), allocatable :: c_new(:)
|
||||
|
||||
allocate(c_new(n-1))
|
||||
|
||||
! Set natural boundary conditions
|
||||
c_new(1) = ZERO
|
||||
z(1) = ZERO
|
||||
z(n) = ZERO
|
||||
|
||||
! Solve using tridiagonal matrix algorithm; first do forward sweep
|
||||
do i = 2, n - 1
|
||||
a = x(i) - x(i-1)
|
||||
c = x(i+1) - x(i)
|
||||
b = TWO * (a + c)
|
||||
d = 6.0_8 * ((y(i+1) - y(i)) / c - (y(i) - y(i-1)) / a)
|
||||
|
||||
c_new(i) = c / (b - a * c_new(i-1))
|
||||
z(i) = (d - a * z(i-1)) / (b - a * c_new(i-1))
|
||||
end do
|
||||
|
||||
! Back substitution
|
||||
do i = n - 1, 1, -1
|
||||
z(i) = z(i) - c_new(i) * z(i+1)
|
||||
end do
|
||||
|
||||
deallocate(c_new)
|
||||
|
||||
end subroutine spline
|
||||
|
||||
!===============================================================================
|
||||
! SPLINE_INTERPOLATE
|
||||
!===============================================================================
|
||||
|
||||
function spline_interpolate(x, y, z, n, xint) result(yint)
|
||||
|
||||
integer, intent(in) :: n
|
||||
real(8), intent(in) :: x(n)
|
||||
real(8), intent(in) :: y(n)
|
||||
real(8), intent(in) :: z(n)
|
||||
real(8), intent(in) :: xint
|
||||
real(8) :: yint
|
||||
|
||||
integer :: i
|
||||
real(8) :: h, r
|
||||
real(8) :: b, c, d
|
||||
|
||||
! Find the lower bounding index in x of xint
|
||||
if (xint < x(1)) then
|
||||
i = 1
|
||||
else if (xint > x(n)) then
|
||||
i = n - 1
|
||||
else
|
||||
i = binary_search(x, n, xint)
|
||||
end if
|
||||
|
||||
h = x(i+1) - x(i)
|
||||
r = xint - x(i)
|
||||
|
||||
! Compute the coefficients
|
||||
b = (y(i+1) - y(i)) / h - (h / 6.0_8) * (z(i+1) + TWO * z(i))
|
||||
c = HALF * z(i)
|
||||
d = (z(i+1) - z(i)) / (h * 6.0_8)
|
||||
|
||||
yint = y(i) + b * r + c * r**2 + d * r**3
|
||||
|
||||
end function spline_interpolate
|
||||
|
||||
!===============================================================================
|
||||
! SPLINE_INTEGRATE
|
||||
!===============================================================================
|
||||
|
||||
function spline_integrate(x, y, z, n, xa, xb) result(s)
|
||||
|
||||
integer, intent(in) :: n
|
||||
real(8), intent(in) :: x(n)
|
||||
real(8), intent(in) :: y(n)
|
||||
real(8), intent(in) :: z(n)
|
||||
real(8), intent(in) :: xa ! Lower limit of integration
|
||||
real(8), intent(in) :: xb ! Upper limit of integration
|
||||
real(8) :: s
|
||||
|
||||
integer :: i
|
||||
integer :: ia, ib
|
||||
real(8) :: h, r
|
||||
real(8) :: a, b, c, d
|
||||
|
||||
! Find the lower bounding index in x of the lower limit of integration.
|
||||
if (xa < x(1)) then
|
||||
ia = 1
|
||||
else if (xa > x(n)) then
|
||||
ia = n - 1
|
||||
else
|
||||
ia = binary_search(x, n, xa)
|
||||
end if
|
||||
|
||||
! Find the lower bounding index in x of the upper limit of integration.
|
||||
if (xb < x(1)) then
|
||||
ib = 1
|
||||
else if (xb > x(n)) then
|
||||
ib = n - 1
|
||||
else
|
||||
ib = binary_search(x, n, xb)
|
||||
end if
|
||||
|
||||
! Evaluate the integral
|
||||
s = ZERO
|
||||
do i = ia, ib
|
||||
h = x(i+1) - x(i)
|
||||
|
||||
! Compute the coefficients
|
||||
b = (y(i+1) - y(i)) / h - (h / 6.0_8) * (z(i+1) + TWO * z(i))
|
||||
c = HALF * z(i)
|
||||
d = (z(i+1) - z(i)) / (h * 6.0_8)
|
||||
|
||||
! Subtract the integral from x(ia) to xa
|
||||
if (i == ia) then
|
||||
r = xa - x(ia)
|
||||
s = s - (y(i) * r + b * r**2 * HALF + c * r**3 / THREE + d * r**4 / FOUR)
|
||||
end if
|
||||
|
||||
! Integrate from x(ib) to xb in final interval
|
||||
if (i == ib) then
|
||||
h = xb - x(ib)
|
||||
end if
|
||||
|
||||
! Accumulate the integral
|
||||
s = s + y(i) * h + b * h**2 * HALF + c * h**3 / THREE + d * h**4 / FOUR
|
||||
end do
|
||||
|
||||
end function spline_integrate
|
||||
|
||||
end module math
|
||||
|
|
|
|||
|
|
@ -2,16 +2,19 @@ module photon_header
|
|||
|
||||
use hdf5, only: HID_T, HSIZE_T, SIZE_T
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use constants, only: ZERO, HALF, SUBSHELLS
|
||||
use dict_header, only: DictIntInt, DictCharInt
|
||||
use endf_header, only: Tabulated1D
|
||||
use algorithm, only: binary_search
|
||||
use constants, only: ZERO, HALF, SUBSHELLS
|
||||
use dict_header, only: DictIntInt, DictCharInt
|
||||
use endf_header, only: Tabulated1D
|
||||
use hdf5_interface
|
||||
use math, only: spline, spline_integrate
|
||||
use material_header, only: Material, materials
|
||||
use nuclide_header, only: nuclides
|
||||
use settings
|
||||
|
||||
real(8), allocatable :: compton_profile_pz(:)
|
||||
real(8), allocatable :: ttb_energy_electron(:) ! incident electron energy grid
|
||||
real(8), allocatable :: ttb_energy_photon(:) ! reduced photon energy grid
|
||||
real(8) :: ttb_energy_cutoff
|
||||
|
||||
type ElectronSubshell
|
||||
integer :: index_subshell ! index in SUBSHELLS
|
||||
|
|
@ -72,7 +75,7 @@ module photon_header
|
|||
type Bremsstrahlung
|
||||
integer :: i_material ! Index in materials array
|
||||
|
||||
real(8), allocatable :: yield(:,:) ! Photon number yield
|
||||
real(8), allocatable :: yield(:) ! Photon number yield
|
||||
real(8), allocatable :: dcs(:,:) ! Scaled bremsstrahlung DCS
|
||||
real(8), allocatable :: cdf(:,:) ! Bremsstrahlung energy CDF
|
||||
|
||||
|
|
@ -80,12 +83,12 @@ module photon_header
|
|||
procedure :: init => bremsstrahlung_init
|
||||
end type Bremsstrahlung
|
||||
|
||||
type(PhotonInteraction), allocatable :: elements(:) ! Photon cross sections
|
||||
type(PhotonInteraction), allocatable, target :: elements(:) ! Photon cross sections
|
||||
integer :: n_elements ! Number of photon cross section tables
|
||||
|
||||
type(DictCharInt) :: element_dict
|
||||
|
||||
type(Bremsstrahlung), allocatable :: ttb(:) ! Bremsstrahlung cross sections
|
||||
type(Bremsstrahlung), allocatable, target :: ttb(:) ! Bremsstrahlung cross sections
|
||||
|
||||
!===============================================================================
|
||||
! ELEMENTMICROXS contains cached microscopic photon cross sections for a
|
||||
|
|
@ -302,9 +305,6 @@ contains
|
|||
this % pair_production_electron
|
||||
|
||||
if (electron_treatment == ELECTRON_TTB) then
|
||||
! TODO: read in
|
||||
ttb_energy_cutoff = 1.e-3
|
||||
|
||||
! Read bremsstrahlung scaled DCS
|
||||
rgroup = open_group(group_id, 'bremsstrahlung')
|
||||
dset_id = open_dataset(rgroup, 'dcs')
|
||||
|
|
@ -385,9 +385,14 @@ contains
|
|||
real(8) :: x_l, x_r, x_c
|
||||
real(8) :: Z_eq_sq
|
||||
real(8) :: beta
|
||||
real(8) :: atom_fraction
|
||||
real(8), allocatable :: mass_fraction(:)
|
||||
real(8), allocatable :: stopping_power(:)
|
||||
real(8), allocatable :: mfp_inv(:)
|
||||
real(8), allocatable :: x(:)
|
||||
real(8), allocatable :: y(:)
|
||||
real(8), allocatable :: z(:)
|
||||
type(DictIntInt) :: nuc_dict
|
||||
type(Material), pointer :: mat
|
||||
type(PhotonInteraction), pointer :: elm
|
||||
|
||||
|
|
@ -398,24 +403,29 @@ contains
|
|||
! Allocate and initialize arrays
|
||||
n_k = size(ttb_energy_photon)
|
||||
n_e = size(ttb_energy_electron)
|
||||
allocate(mass_fraction(size(mat % element)))
|
||||
allocate(mass_fraction(mat % n_nuclides))
|
||||
allocate(stopping_power(n_e))
|
||||
allocate(mfp_inv(n_e))
|
||||
allocate(this % yield(n_e))
|
||||
allocate(this % dcs(n_k, n_e))
|
||||
allocate(this % cdf(n_k, n_e))
|
||||
allocate(x(n_e))
|
||||
allocate(y(n_e))
|
||||
allocate(z(n_e))
|
||||
stopping_power(:) = ZERO
|
||||
mfp_inv(:) = ZERO
|
||||
this % dcs(:,:) = ZERO
|
||||
this % cdf(:,:) = ZERO
|
||||
|
||||
! TODO
|
||||
! Calculate the "equivalent" atomic number Zeq and get the mass fraction of
|
||||
! Calculate the "equivalent" atomic number Zeq and the mass fraction of
|
||||
! each element
|
||||
Z_eq_sq = 0
|
||||
Z_eq_sq = ZERO
|
||||
do i = 1, mat % n_nuclides
|
||||
! Zeq**2 = (atomic fraction of x)*Zx**2 + (atomic fraction of y)*Zy**2 + ...
|
||||
atom_fraction = mat % atom_density(i) / mat % density
|
||||
mass_fraction(i) = atom_fraction * nuclides(mat % nuclide(i)) % awr
|
||||
Z_eq_sq = Z_eq_sq + atom_fraction * nuclides(mat % nuclide(i)) % Z**2
|
||||
end do
|
||||
mass_fraction = mass_fraction / sum(mass_fraction)
|
||||
|
||||
! Calculate the molecular DCS and the molecular total stopping power using
|
||||
! Bragg's additivity rule. Note: the collision stopping power cannot be
|
||||
|
|
@ -426,13 +436,17 @@ contains
|
|||
! form at normal temperature and pressure (at which the NIST stopping
|
||||
! powers are given). It will be used to approximate the collision stopping
|
||||
! powers for now, but should be fixed in the future.
|
||||
do i = 1, size(mat % element)
|
||||
do i = 1, mat % n_nuclides
|
||||
! Get pointer to current element
|
||||
elm => mat % element(i)
|
||||
elm => elements(mat % element(i))
|
||||
|
||||
! TODO: n_atoms
|
||||
! Get atomic fraction
|
||||
atom_fraction = mat % atom_density(i) / mat % density
|
||||
|
||||
! TODO: for molecular DCS, atom_fraction should actually be the number of
|
||||
! atoms in the molecule.
|
||||
! Accumulate material DCS
|
||||
this % dcs = this % dcs + n_atoms(i) * elm % Z**2 / Z_eq_sq * elm % dcs
|
||||
this % dcs = this % dcs + atom_fraction * elm % Z**2 / Z_eq_sq * elm % dcs
|
||||
|
||||
! Accumulate material total stopping power
|
||||
stopping_power = stopping_power + mass_fraction(i) * elm % density * &
|
||||
|
|
@ -442,13 +456,13 @@ contains
|
|||
! Calculate inverse bremsstrahlung mean free path
|
||||
do i = 1, n_e
|
||||
e = exp(ttb_energy_electron(i))
|
||||
if (e <= ttb_energy_cutoff) cycle
|
||||
if (e <= energy_cutoff(PHOTON)) cycle
|
||||
|
||||
! Ratio of the velocity of the charged particle to the speed of light
|
||||
beta = sqrt(e*(e + TWO*MASS_ELECTRON/1.e6)) / (e + MASS_ELECTRON/1.e6)
|
||||
beta = sqrt(e*(e + TWO*MASS_ELECTRON/1.e6_8)) / (e + MASS_ELECTRON/1.e6_8)
|
||||
|
||||
! Integration lower bound
|
||||
k_c = ttb_energy_cutoff / e
|
||||
k_c = energy_cutoff(PHOTON) / e
|
||||
|
||||
! Find the upper bounding index of the reduced photon cutoff energy
|
||||
i_k = binary_search(ttb_energy_photon, n_k, k_c) + 1
|
||||
|
|
@ -470,18 +484,24 @@ contains
|
|||
c = c + HALF * (log(ttb_energy_photon(j+1)) - &
|
||||
log(ttb_energy_photon(j))) * (this % dcs(j,i) + this % dcs(j+1,i))
|
||||
this % cdf(j+1,i) = c
|
||||
end do
|
||||
|
||||
! TODO: density in atom/cm^3
|
||||
mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-27_8
|
||||
! Calculate the inverse bremsstrahlung mean free path
|
||||
mfp_inv(i) = c * mat % density * Z_eq_sq / beta**2 * 1.0e-3_8
|
||||
end do
|
||||
|
||||
! TODO:
|
||||
! Calculate photon number yield
|
||||
! cs = cubic_spline(exp(ttb_energy_electron), mfp_inv / stopping_power)
|
||||
! do i = 1, n_e
|
||||
! yield(i) = cs % integrate(ttb_energy_cutoff, exp(ttb_energy_electron(i)))
|
||||
x = exp(ttb_energy_electron)
|
||||
y = mfp_inv / stopping_power
|
||||
call spline(x, y, z, n_e)
|
||||
do i = 1, n_e
|
||||
this % yield(i) = spline_integrate(x, y, z, n_e, energy_cutoff(PHOTON), x(i))
|
||||
end do
|
||||
|
||||
! Use logarithm of number yield since it is log-log interpolated
|
||||
! yield = log(yield)
|
||||
this % yield = log(this % yield)
|
||||
|
||||
deallocate(mass_fraction, stopping_power, mfp_inv, x, y, z)
|
||||
|
||||
end subroutine bremsstrahlung_init
|
||||
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@ module photon_physics
|
|||
use particle_header, only: Particle
|
||||
use photon_header, only: PhotonInteraction, Bremsstrahlung, &
|
||||
compton_profile_pz, ttb_energy_electron, &
|
||||
ttb_energy_photon, ttb_energy_cutoff
|
||||
ttb_energy_photon, ttb
|
||||
use random_lcg, only: prn
|
||||
|
||||
contains
|
||||
|
|
@ -382,7 +382,7 @@ contains
|
|||
integer :: n
|
||||
integer :: n_e, n_k
|
||||
real(8) :: c_max
|
||||
real(8) :: pi
|
||||
real(8) :: w
|
||||
real(8) :: r
|
||||
real(8) :: e, e_l, e_r
|
||||
real(8) :: y, y_l, y_r
|
||||
|
|
@ -413,17 +413,17 @@ contains
|
|||
! Sample number of secondary bremsstrahlung photons
|
||||
n = floor(y + prn())
|
||||
|
||||
! Calculate the interpolation weight pi_j of the bremsstrahlung energy PDF
|
||||
! Calculate the interpolation weight w_j of the bremsstrahlung energy PDF
|
||||
! interpolated in log energy, which can be interpreted as the probability
|
||||
! of index j
|
||||
pi = (e_r - e) / (e_r - e_l)
|
||||
w = (e_r - e) / (e_r - e_l)
|
||||
|
||||
! Sample the energies of the emitted photons
|
||||
do i = 1, n
|
||||
|
||||
! Sample index of the tabulated PDF in the energy grid, j or j+1
|
||||
i_e = j
|
||||
if (prn() > pi)
|
||||
if (prn() > w) then
|
||||
i_e = i_e + 1
|
||||
end if
|
||||
|
||||
|
|
@ -451,7 +451,7 @@ contains
|
|||
x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l)
|
||||
|
||||
! Determine whether to deliver k
|
||||
if (prn()*max(x_l, x_r) < x) exit
|
||||
if (prn() * max(x_l, x_r) < x) exit
|
||||
end do
|
||||
|
||||
! Create secondary photon
|
||||
|
|
|
|||
|
|
@ -15,7 +15,8 @@ module physics
|
|||
use particle_restart_write, only: write_particle_restart
|
||||
use photon_header
|
||||
use photon_physics, only: rayleigh_scatter, compton_scatter, &
|
||||
atomic_relaxation
|
||||
atomic_relaxation, &
|
||||
thick_target_bremsstrahlung
|
||||
use physics_common
|
||||
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
|
||||
use reaction_header, only: Reaction
|
||||
|
|
@ -304,9 +305,7 @@ contains
|
|||
! TODO: create reaction types
|
||||
|
||||
if (electron_treatment == ELECTRON_TTB) then
|
||||
! TODO: implement thick-target bremsstrahlung model
|
||||
call fatal_error("Thick-target bremsstrahlung treatment of electrons &
|
||||
&is not yet implemented.")
|
||||
call thick_target_bremsstrahlung(p)
|
||||
end if
|
||||
|
||||
p % E = ZERO
|
||||
|
|
@ -333,7 +332,7 @@ contains
|
|||
|
||||
if (electron_treatment == ELECTRON_TTB) then
|
||||
! TODO: implement thick-target bremsstrahlung model
|
||||
call fatal_error("Thick-target bremsstrahlung treatment of electrons &
|
||||
call fatal_error("Thick-target bremsstrahlung treatment of positrons &
|
||||
&is not yet implemented.")
|
||||
end if
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue