From 56a95fb4aa719d70080f301a13a9f7ef7a728b6e Mon Sep 17 00:00:00 2001
From: Brown <2mx@ornl.gov>
Date: Thu, 27 Feb 2020 17:02:41 -0500
Subject: [PATCH] Implement Yield class
Added classes to store yield variables for the ssm module
---
.gitignore | 1 +
docs/tex/sammy_main.tex | 1 +
docs/tex/sammyrefs.bib | 10 +-
docs/tex/scattering-theory.tex | 127 ++++++-
...ions_common.f90 => Convolution_common.f90} | 9 +-
sammy/src/blk/Observable_common.f90 | 12 +
sammy/src/convolution/CMakeLists.txt | 9 +-
sammy/src/convolution/CapYieldCorrections.cpp | 15 +-
sammy/src/convolution/YieldNormalization.cpp | 142 +++++++
sammy/src/convolution/YieldNormalization.h | 60 +++
.../src/convolution/cmake/Dependencies.cmake | 2 +-
.../cix/YieldNormalization.cpp2f.xml | 15 +
.../cpp/YieldNormalizationInterface.cpp | 26 ++
.../cpp/YieldNormalizationInterface.h | 22 ++
.../fortran/YieldNormalization_I.f90 | 32 ++
.../fortran/YieldNormalization_M.f90 | 43 +++
sammy/src/convolution/tests/CMakeLists.txt | 1 +
.../tests/YieldNormalizationTest.cpp | 353 ++++++++++++++++++
sammy/src/endf/CovarianceData.h | 2 +-
sammy/src/inp/minp04.F | 46 +--
sammy/src/inp/minp05.F | 7 +-
sammy/src/mso/mmso5.f | 40 +-
sammy/src/par/mpar0.f | 20 +-
sammy/src/ref/mcon.F | 2 +
sammy/src/salmon/CMakeLists.txt | 11 +
sammy/src/salmon/ExperimentalParameters.cpp | 2 +
sammy/src/salmon/ExperimentalParameters.h | 18 +
sammy/src/salmon/MultipleScattering.cpp | 2 +
sammy/src/salmon/MultipleScattering.h | 36 +-
sammy/src/salmon/Yield.cpp | 32 ++
sammy/src/salmon/Yield.h | 287 ++++++++++++++
.../cix/ExperimentalParameters.cpp2f.xml | 4 +
.../cix/MultipleScattering.cpp2f.xml | 8 +-
.../src/salmon/interface/cix/Yield.cpp2f.xml | 70 ++++
.../cpp/ExperimentalParametersInterface.cpp | 12 +-
.../cpp/ExperimentalParametersInterface.h | 4 +-
.../cpp/MultipleScatteringInterface.cpp | 12 +-
.../cpp/MultipleScatteringInterface.h | 4 +-
.../salmon/interface/cpp/YieldInterface.cpp | 163 ++++++++
.../src/salmon/interface/cpp/YieldInterface.h | 48 +++
.../fortran/ExperimentalParameters_I.f90 | 13 +-
.../fortran/ExperimentalParameters_M.f90 | 16 +-
.../fortran/MultipleScattering_I.f90 | 13 +-
.../fortran/MultipleScattering_M.f90 | 16 +-
.../src/salmon/interface/fortran/Yield_I.f90 | 178 +++++++++
.../src/salmon/interface/fortran/Yield_M.f90 | 226 +++++++++++
sammy/src/salmon/tests/CMakeLists.txt | 1 +
.../tests/ExperimentalParametersTest.cpp | 31 +-
sammy/src/salmon/tests/YieldTest.cpp | 58 +++
sammy/src/sam/msam.F | 2 +
sammy/src/sammy/CMakeLists.txt | 3 +-
sammy/src/ssm/mssm04.f90 | 1 +
sammy/src/ssm/mssm19.f90 | 59 ++-
sammy/src/ssm/mssm20.f90 | 100 +++--
sammy/src/ssm/mssm21.f90 | 130 ++++---
55 files changed, 2332 insertions(+), 225 deletions(-)
rename sammy/src/blk/{CapYCorrections_common.f90 => Convolution_common.f90} (69%)
create mode 100644 sammy/src/blk/Observable_common.f90
create mode 100644 sammy/src/convolution/YieldNormalization.cpp
create mode 100644 sammy/src/convolution/YieldNormalization.h
create mode 100644 sammy/src/convolution/interface/cix/YieldNormalization.cpp2f.xml
create mode 100644 sammy/src/convolution/interface/cpp/YieldNormalizationInterface.cpp
create mode 100644 sammy/src/convolution/interface/cpp/YieldNormalizationInterface.h
create mode 100644 sammy/src/convolution/interface/fortran/YieldNormalization_I.f90
create mode 100644 sammy/src/convolution/interface/fortran/YieldNormalization_M.f90
create mode 100644 sammy/src/convolution/tests/YieldNormalizationTest.cpp
create mode 100644 sammy/src/salmon/Yield.cpp
create mode 100644 sammy/src/salmon/Yield.h
create mode 100644 sammy/src/salmon/interface/cix/Yield.cpp2f.xml
create mode 100644 sammy/src/salmon/interface/cpp/YieldInterface.cpp
create mode 100644 sammy/src/salmon/interface/cpp/YieldInterface.h
create mode 100644 sammy/src/salmon/interface/fortran/Yield_I.f90
create mode 100644 sammy/src/salmon/interface/fortran/Yield_M.f90
create mode 100644 sammy/src/salmon/tests/YieldTest.cpp
diff --git a/.gitignore b/.gitignore
index 90400f6c7..2f77a0a9a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,4 +16,5 @@ sammy/TriBITS
.DS_Store
sammy/[bB][uU][iI][lL][dD]*
/install
+[Ss][Cc][Aa][Ll][Ee]
.vscode
diff --git a/docs/tex/sammy_main.tex b/docs/tex/sammy_main.tex
index da8f82f45..8a7b1c831 100644
--- a/docs/tex/sammy_main.tex
+++ b/docs/tex/sammy_main.tex
@@ -50,6 +50,7 @@
\usepackage{lipsum}
\usepackage[hidelinks]{hyperref} % ,backref
\usepackage[all]{nowidow}
+\usepackage{threeparttable}
%\DeclareMathAlphabet{\pazocal}{OMS}{zplm}{m}{n}
\newcommand{\Ma}{\mathcal{M}}
\newcommand{\Mb}{\pazocal{M}}
diff --git a/docs/tex/sammyrefs.bib b/docs/tex/sammyrefs.bib
index aca4cfa1e..fdd0f044f 100644
--- a/docs/tex/sammyrefs.bib
+++ b/docs/tex/sammyrefs.bib
@@ -1161,7 +1161,6 @@ abstract = "The response of the RPI 16-segment NaI(Tl) multiplicity detector sys
number = {CONF-880546},
}
-<<<<<<< HEAD
%----------- SAMMY specific references ---------------
@conference{larson_multscat_2002,
@@ -1191,9 +1190,12 @@ abstract = "The response of the RPI 16-segment NaI(Tl) multiplicity detector sys
publisher = {Los Alamos Scientific Laboratory},
location = {Los Alamos, New Mexico, USA},
}
-=======
-
->>>>>>> f90b3e53... Added tex files for the manual multiple scattering documentation
+@article{lynn_1968,
+ title={The theory of neutron resonance reactions},
+ author={Lynn, John Eric},
+ year={1968},
+ publisher={Clarendon Press, Oxford}
+}
@Comment{jabref-meta: databaseType:bibtex;}
diff --git a/docs/tex/scattering-theory.tex b/docs/tex/scattering-theory.tex
index d1d27cc67..8bc8d065e 100644
--- a/docs/tex/scattering-theory.tex
+++ b/docs/tex/scattering-theory.tex
@@ -26,26 +26,16 @@ R-matrix theory is expressed in terms of channels, where a channel is defined as
\begin{figure}
\centering
- \includegraphics[width=0.5\textwidth]{figures/rmatrix_inc_channel.pdf} \\[\smallskipamount]
- \includegraphics[width=0.5\textwidth]{figures/rmatrix_interior_region.pdf}\hfill
- \includegraphics[width=0.5\textwidth]{figures/rmatrix_exit_channel.pdf}
+ \includegraphics[width=0.4\textwidth]{figures/rmatrix_inc_channel.pdf} \\[\smallskipamount]
+ \includegraphics[width=0.4\textwidth]{figures/rmatrix_interior_region.pdf}\hfill
+ \includegraphics[width=0.4\textwidth]{figures/rmatrix_exit_channel.pdf}
\caption{Schematic of entrance and exit channels as used in scattering theory. For the interior region (with separation distance $r < a$), no assumptions are made about the nature of the interaction. In the figure, $m$, $i$, and $z$ refer to the mass, spin, and charge of the incident particle while $M$, $I$ and $Z$ refer to the target particle. Orbital angular momentum is denoted by $l$ and velocity by $v$. Primes are used for post-collision quantities.}
\label{scattering-theory_rmatrix_channel_diagram}
\end{figure}
-In Section II.A\ref{}, general equations of scattering theory are presented and their derivations
-discussed. The fundamental R-matrix equations are presented. Section II.A.1\ref{} gives a detailed
-derivation of the equations for a simple case. Section II.A.2\ref{} shows the relationship between the
-R-matrix and the A-matrix, which is another common representation of scattering theory.
+In Section \ref{sec:equations-for-scattering-theory}, general equations of scattering theory are presented and their derivations discussed. The fundamental R-matrix equations are presented. Section II.A.1\ref{} gives a detailed derivation of the equations for a simple case. Section II.A.2\ref{} shows the relationship between the R-matrix and the A-matrix, which is another common representation of scattering theory.
-The approximations to R-matrix theory available in the SAMMY code are detailed in
-Section II.B\ref{}. The recommended choice for most applications is the Reich-Moore approximation,
-described in Section II.B.1\ref{}. For some applications, the Reich-Moore approximation is
-inadequate; for those cases, a method for using SAMMY's Reich-Moore approximation to
-mimic the full (exact) R-matrix is presented Section II.B.2\ref{}. Two historically useful but now
-obsolete approximations are single-level and multilevel Breit Wigner (SLBW and MLBW),
-discussed in Section II.B.3\ref{}. Provisions for including non-compound (direct) effects are
-discussed in Section II.B.4\ref{}.
+The approximations to R-matrix theory available in the SAMMY code are detailed in Section II.B\ref{}. The recommended choice for most applications is the Reich-Moore approximation, described in Section II.B.1\ref{}. For some applications, the Reich-Moore approximation is inadequate; for those cases, a method for using SAMMY's Reich-Moore approximation to mimic the full (exact) R-matrix is presented Section II.B.2\ref{}. Two historically useful but now obsolete approximations are single-level and multilevel Breit Wigner (SLBW and MLBW), discussed in Section II.B.3\ref{}. Provisions for including non-compound (direct) effects are discussed in Section II.B.4\ref{}.
In Section II.C\ref{}, details are given for the SAMMY nomenclature and other conventions,
for transformations to the center-of-momentum system, and for the calculation of penetrability,
@@ -75,6 +65,113 @@ vector sum of the spins of the two particles of the pair: $\vec{s} = \vec{i} + \
of $l$ and $s$: $\vec{J} = \vec{l} + \vec{s}$.
\end{itemize}
+Only $J$ and its associated parity $\pi$ are conserved for any given interaction. The other quantum numbers may differ from channel to channel, as long as the sum rules for spin and parity are obeyed. Within this document and within the SAMMY code, the set of all channels with the same $J$ and $\pi$ are called a ``spin group.''
+
+In all formulae given below, spin quantum numbers (e.g., $J$ ) are implicitly assumed to include the associated parity. Quantized vector sum rules are implicitly assumed to be obeyed. Readers unfamiliar with these sum rules are referred to Section II.C.1.a\ref{} for a mini-tutorial on the subject.
+
+Let the angle-integrated cross sections from entrance channel $c$ to exit channel $c'$ with total angular momentum $J$ be represented by $\sigma_{cc'}$. This cross section is given in terms of the scattering matrix $U_{cc'}$ as
+
+\begin{equation}\label{eq:sigma-ccprime}
+\sigma_{cc'} = \frac{\pi}{k_\alpha^2}g_{J\alpha}\left|e^{2iw_c}\delta_{cc'}-U_{cc'}\right|^2\delta_{JJ'} \:,
+\end{equation}
+
+\noindent
+where $k_\alpha$ is the wave number (and $K_\alpha = \hbar k_\alpha =$ center-of-mass momentum) associated with incident particle pair $\alpha$, $g_{J\alpha}$ is the spin statistical factor, and $w_c$ is the Coulomb phase-shift difference. Note that $w_c$ is zero for non-Coulomb channels. (Details for the charged-particle case are presented in Section II.C.4.\ref{}) The spin statistical factor $g_{J\alpha}$ is given by
+
+\begin{equation}\label{eq:spin-stat-factor}
+g_{J\alpha} = \frac{2J+1}{(2i+1)(2I+1)} \:,
+\end{equation}
+
+\noindent
+and center-of-mass momentum $K_\alpha$ by
+
+\begin{equation}\label{eq:center-of-mass-momentum}
+K_\alpha^2 = (\hbar k_\alpha)^2 = \frac{2mM^2}{(m+M)^2}E \:.
+\end{equation}
+
+\noindent
+Here $E$ is the \textbf{laboratory} kinetic energy of the incident (moving) particle. A derivation of this value for $K_\alpha$ is given in Section II.C.2\ref{}.
+
+The scattering matrix $U$ can be written in terms of matrix $W$ as
+
+\begin{equation}\label{eq:scat-matrix}
+U_{cc'} = \Omega_cW_{cc'}\Omega_{c'} \:,
+\end{equation}
+
+\noindent
+where $\Omega$ is given by
+
+\begin{equation}\label{eq:omega}
+\Omega_c = e^{i(w_c-\phi_c)} \:.
+\end{equation}
+
+\noindent
+Here again, $w_c$ is zero for non-Coulomb channels, and the potential scattering phase shifts for non-Coulomb interactions $\phi_c$ are defined in many references (e.g., \cite{lane_thomas_1958}) and shown in Table \ref{}. The matrix $W$ in Eq. \ref{eq:scat-matrix} is related to the R-matrix (in matrix notation with indices suppressed) via
+
+\begin{equation}\label{eq:W-matrix}
+W = P^{1/2}\left(I-RL\right)^{-1}\left(I-RL^*\right)P^{-1/2} \:.
+\end{equation}
+
+\noindent
+The quantity $I$ in this equation represents the identity matrix, and superscript $*$ indicates a complex conjugate. The form of the R-matrix given in Section IIA.1\ref{} in general Section II.B\ref{} for the versions used in SAMMY. The quantity $L$ in Eq. \ref{eq:W-matrix} is given by
+
+\begin{equation}\label{eq:L-matrix}
+L = (S-B) + iP \:,
+\end{equation}
+
+\noindent
+with $P$ being the penetration factor (penetrability) $S$ the shift factor, and $B$ the arbitrary boundary constant at the channel radius $a_c$. $P$ and $S$ are functions of energy $E$, and also depend on the orbital angular momentum $l$ and the channel radius $a_c$. Formulae for $P$ and $S$ are found in many references (see, for example Eq. (2.9) in \cite{lynn_1968}).
+
+For non-Coulomb interactions, the penetrability and shift factor have the form
+
+\begin{equation}\label{eq:pen-shift-func-of-rho}
+ P\rightarrow P_l(\rho) \qquad \text{and} \qquad S \rightarrow S_l(\rho) \:,
+\end{equation}
+
+\noindent
+where $\rho$ is related to the center-of-mass momentum which in turn is related to the laboratory energy of the incident particle $(E)$. For arbitrary channel $c$ with a particle pair $\alpha$, orbital angular momentum $l$, and channel radius $a_c$, $\rho$ has the form
+
+\begin{equation}\label{eq:de-Broglie-radius}
+ \rho = k_\alpha a_c = \frac{1}{\hbar} \sqrt{\frac{2m_\alpha M_\alpha}{m_\alpha+M_\alpha} \frac{M}{m+M}} \sqrt{(E-\Xi_\alpha)}\: a_c \:,
+\end{equation}
+
+as shown in Section II.C.2\ref{}. Here $\Xi_\alpha$ is the energy threshold for the particle pair $\alpha$, $m_\alpha$ and $M_\alpha$ are the masses of the two particles of particle pair $\alpha$, and $m$ and $M$ are the masses of the incident particle and target nuclide, respectively.
+
+Appropriate formulae for $P$, $S$, and $\phi$ in the non-Coulomb case are shown in Table IIA.1. For two charged particles, formulae for the penetrabilities are given in Section II.C.4\ref{}.
+
+The energy dependence of fission and capture widths is negligible over the energy range of these calculations. Therefore, a penetrability of unity may be used.
+
+% multiline cell: \begin{tabular}{@{}c@{}} line1 \\ line2 \end{tabular}
+\begin{threeparttable}
+\centering
+\caption{Hard-sphere penetrability (penetration factor) $P$, level shift factor $S$, and potential-scattering phase shift $\phi$ for orbital angular momentum $l$, wave number $k$, and channel radius $a_c$, with $\rho=ka_c$} \label{tab:penetrabilities}
+\begin{tabular}{ l c c c }
+ \hline\hline
+ $\mathbf{l}$ & $\mathbf{P_l}$ & $\mathbf{S_l}$ & $\mathbf{\phi_l}$ \cr
+ \hline\hline
+ 0 & $\rho$ & 0 & $\rho$ \cr
+ 1 & $\rho^3/(1+\rho^2)$ & $-1/(1+\rho^2)$ & $\rho-tan^{-1}\rho$ \cr
+ 2 & $\rho^5/(9+3\rho^2+\rho^4)$ & $-(18+3\rho^2)/(9+3\rho^2+\rho^4)$ & $\rho-tan^{-1}\left[3\rho/(3-\rho^2)\right]$ \cr
+ 3 & \begin{tabular}{@{}c@{}}$\rho^7/$ \\ $(225+45\rho^2+6\rho^4+\rho^6)$\end{tabular} & \begin{tabular}{@{}c@{}} $-(675+90\rho^2+6\rho^4)/$ \\ $(225+45\rho^2+6\rho^4+\rho^6)$ \end{tabular} & $\rho-tan^{-1}\left[\rho(15-\rho^2)/(15-6\rho^2)\right]$ \cr
+ \hline
+ $l$ & $\frac{\rho^2P_{l-1}}{(1-S_{l-1})^2+P_{l-1}^2}$ & $\frac{\rho^2(l-S_{l-1})}{(1-S_{l-1})^2+P_{l-1}^2}-l$ & $\phi_{l-1} - tan^{-1}\left(P_{l-1}/(l-S_{l-1})\right)$\tnote{\textdagger} \cr
+
+ \hline\hline
+
+\end{tabular}
+\begin{tablenotes}
+ \item[\textdagger] \footnotesize{The iterative formula for $\phi_l$ could also be defined by $B_l = (B_{l-1}+X_l)/(1-B_{l-1}X_l)$ where $B_l=tan(\rho-\phi_l)$ and $X_l=P_{l-1}/(l-S_l)$}
+\end{tablenotes}
+\end{threeparttable}
+
+
+
+
+
+
+
+
+
diff --git a/sammy/src/blk/CapYCorrections_common.f90 b/sammy/src/blk/Convolution_common.f90
similarity index 69%
rename from sammy/src/blk/CapYCorrections_common.f90
rename to sammy/src/blk/Convolution_common.f90
index e4427788e..43c937565 100644
--- a/sammy/src/blk/CapYCorrections_common.f90
+++ b/sammy/src/blk/Convolution_common.f90
@@ -3,10 +3,13 @@
! * used to pass the variables contained in C++ classes until
! * the code has been fully modernized.
! **********************************************************/
-module CapYCorrections_common_m
- use CapYieldCorrections_M
+module Convolution_common_m
+ use CapYieldCorrections_M
+ use YieldNormalization_M
IMPLICIT NONE
type(CapYieldCorrections)::capYieldCor
-end module CapYCorrections_common_m
+ type(YieldNormalization)::yNorm
+
+end module Convolution_common_m
diff --git a/sammy/src/blk/Observable_common.f90 b/sammy/src/blk/Observable_common.f90
new file mode 100644
index 000000000..7fd5eb21f
--- /dev/null
+++ b/sammy/src/blk/Observable_common.f90
@@ -0,0 +1,12 @@
+! /**********************************************************
+! * This common block style module is a temporary solution
+! * used to pass the variables contained in C++ classes until
+! * the code has been fully modernized.
+! **********************************************************/
+module Observable_common_m
+ use Yield_M
+ IMPLICIT NONE
+
+ type(Yield)::yld
+
+end module Observable_common_m
diff --git a/sammy/src/convolution/CMakeLists.txt b/sammy/src/convolution/CMakeLists.txt
index 81c874b58..f576ac814 100644
--- a/sammy/src/convolution/CMakeLists.txt
+++ b/sammy/src/convolution/CMakeLists.txt
@@ -11,8 +11,10 @@ INCLUDE_DIRECTORIES(${CMAKE_CURRENT_BINARY_DIR})
# Set headers
SET(HEADERS
CapYieldCorrections.h
-
interface/cpp/CapYieldCorrectionsInterface.h
+
+ YieldNormalization.h
+ interface/cpp/YieldNormalizationInterface.h
)
APPEND_SET(CONVOLUTION_SOURCES
@@ -21,6 +23,11 @@ APPEND_SET(CONVOLUTION_SOURCES
interface/cpp/CapYieldCorrectionsInterface.cpp
interface/fortran/CapYieldCorrections_I.f90
interface/fortran/CapYieldCorrections_M.f90
+
+ YieldNormalization.cpp
+ interface/cpp/YieldNormalizationInterface.cpp
+ interface/fortran/YieldNormalization_I.f90
+ interface/fortran/YieldNormalization_M.f90
)
diff --git a/sammy/src/convolution/CapYieldCorrections.cpp b/sammy/src/convolution/CapYieldCorrections.cpp
index ace8fdc6c..8c8f736ac 100644
--- a/sammy/src/convolution/CapYieldCorrections.cpp
+++ b/sammy/src/convolution/CapYieldCorrections.cpp
@@ -86,6 +86,9 @@ namespace sammy{
if( N_sigma <= 0.0 ){
std::ostringstream err;
err << "(sample thickness [at/b]) * (tot. x.s. [b]) is less than zero" << std::endl;
+ err << "--------------------------------" << std::endl;
+ err << "N * sigma_tot = " << N_sigma << std::endl;
+ err << "--------------------------------" << std::endl;
throw std::runtime_error(err.str());
}
@@ -94,10 +97,11 @@ namespace sammy{
std::ostringstream err;
err << "Under or overflow may occur in exponent. Check sample thickness" << std::endl;
err << " and total cross section values." << std::endl;
+ err << "--------------------------------" << std::endl;
+ err << "N * sigma_tot = " << N_sigma << std::endl;
+ err << "--------------------------------" << std::endl;
throw std::runtime_error(err.str());
}
-
-
}
void CapYieldCorrections::check_attenFactor(const double &attenFactor)
@@ -106,7 +110,9 @@ namespace sammy{
if( attenFactor <= 0.0 ){
std::ostringstream err;
err << "Gamma attenuation factor C1 (from documentation) is less than zero" << std::endl;
+ err << "--------------------------------" << std::endl;
err << "C1 = " << attenFactor << std::endl;
+ err << "--------------------------------" << std::endl;
throw std::runtime_error(err.str());
}
@@ -115,7 +121,9 @@ namespace sammy{
std::ostringstream err;
err << "Gamma attenuation factor C1 (from documentation) is too large, and" << std::endl;
err << " may cause under/overflow issues." << std::endl;
+ err << "--------------------------------" << std::endl;
err << "C1 = " << attenFactor << std::endl;
+ err << "--------------------------------" << std::endl;
throw std::runtime_error(err.str());
}
}
@@ -126,9 +134,10 @@ namespace sammy{
std::ostringstream err;
err << "The gamma attenuation correction to the yield is less than or" << std::endl;
err << "equal to zero or greater than 1.0" << std::endl;
+ err << "--------------------------------" << std::endl;
err << "k_a = " << gammaCor << std::endl;
+ err << "--------------------------------" << std::endl;
throw std::runtime_error(err.str());
}
}
-
}
\ No newline at end of file
diff --git a/sammy/src/convolution/YieldNormalization.cpp b/sammy/src/convolution/YieldNormalization.cpp
new file mode 100644
index 000000000..1c0ed3720
--- /dev/null
+++ b/sammy/src/convolution/YieldNormalization.cpp
@@ -0,0 +1,142 @@
+
+#include "YieldNormalization.h"
+
+#include <iostream>
+namespace sammy{
+
+ YieldNormalization::YieldNormalization(){}
+
+ void YieldNormalization::norm(sammy::Yield & yld, sammy::SampleDimensions & sd,
+ sammy::MultipleScattering &multScat, double & sigTot)
+ {
+ // How many scatters considered before neutron absorption?
+ // (We change the behavior of normTypes == 1 and 2 for different scatters)
+ int NS = multScat.getNumScatters();
+
+ if( yld.normTypeIsSelfShielded() )
+ {
+ selfShielded(yld,sd);
+ }
+ else if( yld.normTypeIsDivideByN() )
+ {
+ divideByN(yld, sd, NS);
+ }
+ else if( yld.normTypeIsTimesSigTot() )
+ {
+ timesSigmaTot(yld, sd, sigTot, NS);
+ }
+ else
+ {
+ std::ostringstream err;
+ err << "------------------------------------------" << std::endl;
+ err << "The yield normalization type has not been set, check input file." << std::endl;
+ err << "------------------------------------------" << std::endl;
+ throw std::runtime_error(err.str());
+ }
+ }
+
+ void YieldNormalization::selfShielded(sammy::Yield &yld, sammy::SampleDimensions &sd)
+ {
+ yld.setY0( yld.getY0() * sd.getArealDensity() );
+ }
+
+ void YieldNormalization::divideByN(sammy::Yield &yld, sammy::SampleDimensions &sd, int numScatters)
+ {
+ if( numScatters == 0 )
+ {
+ // do nothing
+ }
+
+ if( numScatters == 1 or numScatters == 2 )
+ {
+ // Yield
+ yld.setY1( yld.getY1() / sd.getArealDensity() );
+
+ // if numParam == 0, we don't need to calculate
+ for( int i=0; i<yld.getNumTotalParam(); i++ )
+ {
+ // Derivative of yield
+ yld.setDerivY1(i, yld.getDerivY1(i) / sd.getArealDensity() );
+ }
+ }
+
+ // Really if the number of scatters is >= 2
+ if( numScatters == 2 )
+ {
+ // Yield
+ yld.setY2( yld.getY2() / sd.getArealDensity() );
+
+ // if numParam == 0, we don't need to calculate
+ for( int i=0; i<yld.getNumTotalParam(); i++ )
+ {
+ // Derivative of yield
+ yld.setDerivY2(i, yld.getDerivY2(i) / sd.getArealDensity() );
+ }
+ }
+ }
+
+ void YieldNormalization::timesSigmaTot(sammy::Yield &yld, sammy::SampleDimensions &sd,
+ double & sigTot, int numScatters)
+ {
+ // for all cases
+ yld.setY0( yld.getY0() * sigTot * sd.getArealDensity() );
+
+ if( numScatters == 0 )
+ {
+ // nothing extra
+ }
+
+ if( numScatters == 1 or numScatters == 2 )
+ {
+ yld.setY1( yld.getY1() * sigTot );
+
+ // if numParam == 0, we don't need to calculate
+ for( int i=0; i<yld.getNumTotalParam(); i++ )
+ {
+ yld.setDerivY1(i, yld.getDerivY1(i) * sigTot );
+ }
+ }
+
+ if( numScatters == 2 )
+ {
+ yld.setY2( yld.getY2() * sigTot );
+
+ // if numParam == 0, we don't need to calculate
+ for( int i=0; i<yld.getNumTotalParam(); i++ )
+ {
+ yld.setDerivY2(i, yld.getDerivY2(i) * sigTot );
+ }
+ }
+ }
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sammy/src/convolution/YieldNormalization.h b/sammy/src/convolution/YieldNormalization.h
new file mode 100644
index 000000000..b592d3f83
--- /dev/null
+++ b/sammy/src/convolution/YieldNormalization.h
@@ -0,0 +1,60 @@
+#ifndef YIELDNORMALIZATION_H
+#define YIELDNORMALIZATION_H
+
+#include "../salmon/Yield.h"
+#include "../salmon/ExperimentalParameters.h"
+#include "../salmon/MultipleScattering.h"
+
+namespace sammy{
+
+ /** @brief Normalization of the observed yield
+
+ This class applies normalization to yields (capture, fission, scattering...)
+ @date February 2020
+ */
+ class YieldNormalization
+ {
+ public:
+ YieldNormalization();
+ YieldNormalization(const YieldNormalization& orig);
+ virtual ~YieldNormalization(){};
+
+ /**
+ * Master normalization function
+ */
+ void norm(sammy::Yield & yld, sammy::SampleDimensions & sd,
+ sammy::MultipleScattering &multScat, double & sigTot );
+
+ /**
+ * Normalize as self-shielded yield \f$ Y \f$
+ * \f[ = \frac{\sigma_x}{\sigma_t} \left(1-e^{-n\sigma_t}\right) \f]
+ */
+ void selfShielded(sammy::Yield & yld, sammy::SampleDimensions &sd);
+
+ /**
+ * Normalize as self-shielded yield divided by sample thickness \f$ Y/n \f$ which, in the
+ * thin sample approximation, is \f$ \approx \sigma_x \f$
+ * \f[ = \frac{\sigma_\gamma}{\sigma_t n} \left(1-e^{-n\sigma_t}\right) \f]
+ */
+ void divideByN(sammy::Yield & yld, sammy::SampleDimensions &sd, int numScatters);
+
+ /**
+ * Normalize as self-shielded yield \f$ Y \f$ multiplied by total cross section \f$ \sigma_t \f$
+ * \f[ = \sigma_x \left(1-e^{-n\sigma_t}\right) \f]
+ */
+ void timesSigmaTot(sammy::Yield & yld, sammy::SampleDimensions &sd,
+ double & sigTot, int numScatters);
+
+ /**
+ * Check that the number of varied parameters is greater than zero if
+ * Bayes' equations are to be solved
+ * @param yld A sammy::Yield object
+ */
+ //void checkNumParam(sammy::Yield & yld);
+
+ };
+
+}
+
+
+#endif // YIELDNORMALIZATION_H
\ No newline at end of file
diff --git a/sammy/src/convolution/cmake/Dependencies.cmake b/sammy/src/convolution/cmake/Dependencies.cmake
index 13244a48a..1cbd56532 100644
--- a/sammy/src/convolution/cmake/Dependencies.cmake
+++ b/sammy/src/convolution/cmake/Dependencies.cmake
@@ -1,4 +1,4 @@
-SET(LIB_REQUIRED_DEP_PACKAGES)
+SET(LIB_REQUIRED_DEP_PACKAGES salmon)
SET(LIB_OPTIONAL_DEP_PACKAGES)
SET(TEST_REQUIRED_DEP_PACKAGES)
SET(TEST_OPTIONAL_DEP_PACKAGES)
diff --git a/sammy/src/convolution/interface/cix/YieldNormalization.cpp2f.xml b/sammy/src/convolution/interface/cix/YieldNormalization.cpp2f.xml
new file mode 100644
index 000000000..c3fc0a03e
--- /dev/null
+++ b/sammy/src/convolution/interface/cix/YieldNormalization.cpp2f.xml
@@ -0,0 +1,15 @@
+<generate name="YieldNormalization">
+ <include_relative name="../../YieldNormalization.h"/>
+ <using_namespace name="sammy"/>
+
+ <class name="YieldNormalization">
+
+ <method name="norm">
+ <param name="yld" type="Yield"/>
+ <param name="sampDim" type="SampleDimensions"/>
+ <param name="multscat" type="MultipleScattering"/>
+ <param name="sigTot" type="double"/>
+ </method>
+
+ </class>
+</generate>
\ No newline at end of file
diff --git a/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.cpp b/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.cpp
new file mode 100644
index 000000000..908f0efc9
--- /dev/null
+++ b/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.cpp
@@ -0,0 +1,26 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX)
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Mon Mar 02 10:18:42 EST 2020
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#include <string.h>
+#include "YieldNormalizationInterface.h"
+using namespace sammy;
+void YieldNormalization_norm(void * YieldNormalization_ptr,Yield * yld,SampleDimensions * sampDim,MultipleScattering * multscat,double * sigTot)
+{
+ ((YieldNormalization*)YieldNormalization_ptr)->norm(*yld,*sampDim,*multscat,*sigTot);
+}
+
+void* YieldNormalization_initialize()
+{
+ return new YieldNormalization();
+}
+
+void YieldNormalization_destroy(void * YieldNormalization_ptr)
+{
+ delete (YieldNormalization*)YieldNormalization_ptr;
+}
+
diff --git a/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.h b/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.h
new file mode 100644
index 000000000..36db444e9
--- /dev/null
+++ b/sammy/src/convolution/interface/cpp/YieldNormalizationInterface.h
@@ -0,0 +1,22 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX)
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Mon Mar 02 10:18:42 EST 2020
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#ifndef YIELDNORMALIZATIONINTERFACE_H
+#define YIELDNORMALIZATIONINTERFACE_H
+#include "../../YieldNormalization.h"
+using namespace sammy;
+#ifdef __cplusplus
+extern "C" {
+#endif
+void YieldNormalization_norm(void * YieldNormalization_ptr,Yield * yld,SampleDimensions * sampDim,MultipleScattering * multscat,double * sigTot);
+void* YieldNormalization_initialize();
+void YieldNormalization_destroy(void * YieldNormalization_ptr);
+#ifdef __cplusplus
+}
+#endif
+#endif /* YIELDNORMALIZATIONINTERFACE_H */
diff --git a/sammy/src/convolution/interface/fortran/YieldNormalization_I.f90 b/sammy/src/convolution/interface/fortran/YieldNormalization_I.f90
new file mode 100644
index 000000000..4610ddb66
--- /dev/null
+++ b/sammy/src/convolution/interface/fortran/YieldNormalization_I.f90
@@ -0,0 +1,32 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX)
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Mon Mar 02 10:18:42 EST 2020
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module YieldNormalization_I
+use, intrinsic :: ISO_C_BINDING
+interface
+subroutine f_YieldNormalization_norm(YieldNormalization_ptr, yld,sampDim,multscat,sigTot ) BIND(C,name="YieldNormalization_norm")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: YieldNormalization_ptr;
+ type(C_PTR), value :: yld;
+ type(C_PTR), value :: sampDim;
+ type(C_PTR), value :: multscat;
+ real(C_DOUBLE) :: sigTot;
+end subroutine
+type(C_PTR) function f_YieldNormalization_initialize( )BIND(C,name="YieldNormalization_initialize")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR) :: YieldNormalization_ptr;
+end function
+subroutine f_YieldNormalization_destroy(this) BIND(C,name="YieldNormalization_destroy")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: this;
+end subroutine
+end interface
+end module YieldNormalization_I
diff --git a/sammy/src/convolution/interface/fortran/YieldNormalization_M.f90 b/sammy/src/convolution/interface/fortran/YieldNormalization_M.f90
new file mode 100644
index 000000000..0473ffa31
--- /dev/null
+++ b/sammy/src/convolution/interface/fortran/YieldNormalization_M.f90
@@ -0,0 +1,43 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX)
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Mon Mar 02 10:18:42 EST 2020
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module YieldNormalization_M
+use, intrinsic :: ISO_C_BINDING
+use YieldNormalization_I
+use MultipleScattering_M
+use ExperimentalParameters_M
+use Yield_M
+type YieldNormalization
+ type(C_PTR) :: instance_ptr=C_NULL_PTR
+ contains
+ procedure, pass(this) :: norm => YieldNormalization_norm
+ procedure, pass(this) :: initialize => YieldNormalization_initialize
+ procedure, pass(this) :: destroy => YieldNormalization_destroy
+end type YieldNormalization
+contains
+subroutine YieldNormalization_norm(this, yld, sampDim, multscat, sigTot)
+ implicit none
+ class(YieldNormalization)::this
+ class(Yield)::yld
+ class(SampleDimensions)::sampDim
+ class(MultipleScattering)::multscat
+ real(C_DOUBLE)::sigTot
+ call f_YieldNormalization_norm(this%instance_ptr, yld%instance_ptr,sampDim%instance_ptr,multscat%instance_ptr,sigTot)
+end subroutine
+subroutine YieldNormalization_initialize(this)
+ implicit none
+ class(YieldNormalization) :: this
+ this%instance_ptr = f_YieldNormalization_initialize()
+end subroutine
+subroutine YieldNormalization_destroy(this)
+ implicit none
+ class(YieldNormalization) :: this
+ call f_YieldNormalization_destroy(this%instance_ptr)
+ this%instance_ptr = C_NULL_PTR
+end subroutine
+end module YieldNormalization_M
diff --git a/sammy/src/convolution/tests/CMakeLists.txt b/sammy/src/convolution/tests/CMakeLists.txt
index 50399bf91..1fb53c922 100644
--- a/sammy/src/convolution/tests/CMakeLists.txt
+++ b/sammy/src/convolution/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
include(NemesisTest)
ADD_NEMESIS_TEST(CapYieldCorrectionsTest.cpp TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(YieldNormalizationTest.cpp TIMEOUT 60000 NP 1 )
diff --git a/sammy/src/convolution/tests/YieldNormalizationTest.cpp b/sammy/src/convolution/tests/YieldNormalizationTest.cpp
new file mode 100644
index 000000000..6dd8ee1a9
--- /dev/null
+++ b/sammy/src/convolution/tests/YieldNormalizationTest.cpp
@@ -0,0 +1,353 @@
+#include "Nemesis/gtest/Gtest_Functions.hh"
+#include "Nemesis/gtest/nemesis_gtest.hh"
+#include "../YieldNormalization.h"
+#include <exception>
+#include <stdexcept>
+
+using namespace sammy;
+
+/**************************************************************************************************
+* TODO: need to include tests for the Observable class members of yield from Observable<Yield>
+**************************************************************************************************/
+
+/**************************************************************************************************
+* Test norm function: multiply by tot x.s., scattering = 0
+**************************************************************************************************/
+TEST( YieldNormalization,multTot_scat0 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeTimesSigTot();
+ ms.setNumScatters(0);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setDerivY0(0,-0.01);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * sigTot * 0.01;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+
+}
+
+/**************************************************************************************************
+* Test norm function: multiply by tot x.s., scattering = 1
+**************************************************************************************************/
+TEST( YieldNormalization,multTot_scat1 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeTimesSigTot();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(1);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * sigTot * 0.01;
+ double gold_Y1 = 0.2 * sigTot;
+ double gold_DerivY1 = -0.02 * sigTot;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: multiply by tot x.s., scattering = 2
+**************************************************************************************************/
+TEST( YieldNormalization,multTot_scat2 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeTimesSigTot();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(2);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ yld.setY2(0.3);
+ yld.setDerivY2(0,-0.03);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * sigTot * 0.01;
+ double gold_Y1 = 0.2 * sigTot;
+ double gold_Y2 = 0.3 * sigTot;
+
+ double gold_DerivY1 = -0.02 * sigTot;
+ double gold_DerivY2 = -0.03 * sigTot;
+
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_Y2, yld.getY2());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+ ASSERT_EQ(gold_DerivY2,yld.getDerivY2(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: divide by thick, scattering = 0
+**************************************************************************************************/
+TEST( YieldNormalization,divideByN_scat0 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeDivideByN();
+ ms.setNumScatters(0);
+
+ yld.setY0(0.1);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+
+}
+
+/**************************************************************************************************
+* Test norm function: divide by thick, scattering = 1
+**************************************************************************************************/
+TEST( YieldNormalization,divideByN_scat1 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeDivideByN();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(1);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1;
+ double gold_Y1 = 0.2 / 0.01;
+ double gold_DerivY1 = -0.02 / 0.01;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: divide by thick, scattering = 2
+**************************************************************************************************/
+TEST( YieldNormalization,divideByN_scat2 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeDivideByN();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(2);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ yld.setY2(0.3);
+ yld.setDerivY2(0,-0.03);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1;
+ double gold_Y1 = 0.2 / 0.01;
+ double gold_DerivY1 = -0.02 / 0.01;
+ double gold_Y2 = 0.3 / 0.01;
+ double gold_DerivY2 = -0.03 / 0.01;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+ ASSERT_EQ(gold_Y2, yld.getY2());
+ ASSERT_EQ(gold_DerivY2,yld.getDerivY2(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: self-shielded norm, scattering = 0
+**************************************************************************************************/
+TEST( YieldNormalization,selfShielded_scat0 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeSelfShielded();
+ ms.setNumScatters(0);
+
+ yld.setY0(0.1);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * 0.01;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+
+}
+
+/**************************************************************************************************
+* Test norm function: self-shielded norm, scattering = 1
+**************************************************************************************************/
+TEST( YieldNormalization,selfShielded_scat1 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeSelfShielded();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(1);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * 0.01;
+ double gold_Y1 = 0.2;
+ double gold_DerivY1 = -0.02;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: self-shielded norm, scattering = 2
+**************************************************************************************************/
+TEST( YieldNormalization,selfShielded_scat2 ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ YieldNormalization yNorm;
+
+ // These change test-to-test
+ yld.setNormTypeSelfShielded();
+ yld.setCalcDerivs(true);
+ ms.setNumScatters(2);
+ yld.setNumTotalParam(1);
+
+ yld.setY0(0.1);
+ yld.setY1(0.2);
+ yld.setDerivY1(0,-0.02);
+ yld.setY2(0.3);
+ yld.setDerivY2(0,-0.03);
+ double sigTot = 3.0;
+ sd.setArealDensity(0.01);
+
+ yNorm.norm(yld, sd, ms, sigTot);
+
+ double gold_Y0 = 0.1 * 0.01;
+ double gold_Y1 = 0.2;
+ double gold_DerivY1 = -0.02;
+ double gold_Y2 = 0.3;
+ double gold_DerivY2 = -0.03;
+
+ ASSERT_EQ(gold_Y0, yld.getY0());
+ ASSERT_EQ(gold_Y1, yld.getY1());
+ ASSERT_EQ(gold_DerivY1,yld.getDerivY1(0));
+ ASSERT_EQ(gold_Y2, yld.getY2());
+ ASSERT_EQ(gold_DerivY2,yld.getDerivY2(0));
+
+}
+
+/**************************************************************************************************
+* Test norm function: normalization not set, throw the user an error
+**************************************************************************************************/
+TEST( YieldNormalization,normNotSet ){
+
+ Yield yld;
+ SampleDimensions sd;
+ MultipleScattering ms;
+
+ double sigTot = 3.0;
+
+ YieldNormalization yNorm;
+
+ // Expect to throw error
+ int gold_result = 0;
+ int result;
+
+ try{
+ yNorm.norm(yld, sd, ms, sigTot);
+ result = 1;
+ }
+ catch(...){
+ result = 0;
+ }
+
+
+ ASSERT_EQ(gold_result,result);
+}
diff --git a/sammy/src/endf/CovarianceData.h b/sammy/src/endf/CovarianceData.h
index 9c1031ea5..602ae91a6 100644
--- a/sammy/src/endf/CovarianceData.h
+++ b/sammy/src/endf/CovarianceData.h
@@ -173,7 +173,7 @@ namespace sammy{
/**
- * Get the number of total parameters
+ * Get the number of total parameters: The pup'd and non-pup'd
*/
int getNumTotalParam() const{
return (int)paramIndex.size();
diff --git a/sammy/src/inp/minp04.F b/sammy/src/inp/minp04.F
index e58f420c3..2e68ef8e7 100644
--- a/sammy/src/inp/minp04.F
+++ b/sammy/src/inp/minp04.F
@@ -14,6 +14,7 @@ C
use lbro_common_m
use misccc_common_m
use ntyp_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
C
@@ -81,25 +82,26 @@ C
C
IF (Kdebug.EQ.0) Debug = .False.
IF (Kdebug.EQ.1) Debug = .True.
- IF (Kyield.EQ.-1 .AND. (Kssmsc.GT.0 .OR. Kssmsc.EQ.-2)) THEN
- WRITE (6,10000)
- WRITE (21,10000)
-10000 FORMAT (/, ' ### CAUTION -- Normalization method has not been
- * specified in INPut file. ###', /,
- * ' ### See page 56r of the SAMMY users manual for a des
- *cription ###', /,
- * ' ### of the various options. Please add the appropria
- *te command ###', /,
- * ' ### to your INPut file and resubmit. Possibilities x
- *include ###', /,
- * ' ### NORMALIZE AS CROSS SECTION or CROSS SECTION x
- * ###', /,
- * ' ### NORMALIZE AS YIELD RATHER or YIELD x
- * ###', /,
- * ' ### NORMALIZE AS (1-E)SIGMA x
- * ###')
- STOP
- END IF
+C IF (yld%getNormType().EQ.-1 .AND.
+C * (Kssmsc.GT.0 .OR. Kssmsc.EQ.-2)) THEN
+C WRITE (6,10000)
+C WRITE (21,10000)
+C 10000 FORMAT (/, ' ### CAUTION -- Normalization method has not been
+C * specified in INPut file. ###', /,
+C * ' ### See page 56r of the SAMMY users manual for a des
+C *cription ###', /,
+C * ' ### of the various options. Please add the appropria
+C *te command ###', /,
+C * ' ### to your INPut file and resubmit. Possibilities x
+C *include ###', /,
+C * ' ### NORMALIZE AS CROSS SECTION or CROSS SECTION x
+C * ###', /,
+C * ' ### NORMALIZE AS YIELD RATHER or YIELD x
+C * ###', /,
+C * ' ### NORMALIZE AS (1-E)SIGMA x
+C * ###')
+C STOP
+C END IF
C
IF (Mretro.EQ.1) THEN
IF (Nretro.EQ.1) THEN
@@ -189,6 +191,7 @@ C
use broad_common_m
use namfil_common_m
use lbro_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
C
#ifdef COMPILE_WITH_SAMINT
@@ -533,11 +536,12 @@ C
C *** Do not make plot file of multiple-scattering pieces (-268,269)
Kwssms = 0
C
-C *** Normalize as cross section rather than yield
+C *** Normalize as cross section rather than yield
+C *** Kyield=yld%getNormType()
C *** for Kyield=1: (270,271) cross section
C *** for Kyield=0: (272,273) yield
C *** for Kyield=2: (274 ) 1-e
- Kyield = -1
+ !call yld%setNormType(-1)
C
C *** Do not prepare ssm comparison file (-275)
C *** (i.e., do not print multiple-scattering corrections)
diff --git a/sammy/src/inp/minp05.F b/sammy/src/inp/minp05.F
index 2c6bd23e0..6292686e5 100644
--- a/sammy/src/inp/minp05.F
+++ b/sammy/src/inp/minp05.F
@@ -15,6 +15,7 @@ C
use namfil_common_m
use lbro_common_m
use misccc_common_m
+ use Observable_common_m
use mdf5_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
@@ -416,11 +417,11 @@ C
ELSE IF (Number.EQ.268 .OR. Number.EQ.269) THEN
Kwssms = 1
ELSE IF (Number.EQ.270 .OR. Number.EQ.271) THEN
- Kyield = 1
+ call yld%setNormTypeDivideByN()
ELSE IF (Number.EQ.272 .OR. Number.EQ.273) THEN
- Kyield = 0
+ call yld%setNormTypeSelfShielded()
ELSE IF (Number.EQ.274) THEN
- Kyield = 2
+ call yld%setNormTypeTimesSigTot()
ELSE IF (Number.EQ.275) THEN
Kssmpr = 1
ELSE IF (Number.EQ.276 .OR. Number.eq.277) THEN
diff --git a/sammy/src/mso/mmso5.f b/sammy/src/mso/mmso5.f
index 792fd5a5f..23fff9a61 100644
--- a/sammy/src/mso/mmso5.f
+++ b/sammy/src/mso/mmso5.f
@@ -487,13 +487,14 @@ C
use fixedi_m
use ifwrit_m
use fixedr_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION Capsig(*), Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*),
* Dy2(*), Sigxxx(*), Dasigx(*), Dbsigx(*)
DATA Zero /0.0d0/, Half /0.5d0/, One /1.0d0/
C
C *** Default normalization is such that Y0 = (1-e)Capture /(n Total)
-C *** Which corresponds to Kyield=1
+C *** Which corresponds to yld%normTypeIsDivideByN()==true
C
Cap = Capsig(Nn)
IF (Kssmsc.NE.-2) THEN
@@ -522,15 +523,15 @@ C
Sigxxx(1) = Y0
IF (Kssmsc.NE.-2) Sigxxx(1) = Sigxxx(1) + Y1
IF (Kssdbl.EQ. 1) Sigxxx(1) = Sigxxx(1) + Y2
- IF (Kyield.EQ. 0) Sigxxx(1) = Sigxxx(1)*Dthick
- IF (Kyield.EQ. 2) Sigxxx(1) = Sigxxx(1)*Total*Dthick
+ IF (yld%normTypeIsSelfShielded() ) Sigxxx(1) = Sigxxx(1)*Dthick
+ IF (yld%normTypeIsTimesSigTot()) Sigxxx(1) =Sigxxx(1)*Total*Dthick
C
IF (Kssmpr.EQ.1) THEN
- IF (Kyield.EQ.0) THEN
+ IF ( yld%normTypeIsSelfShielded() ) THEN
Y0 = Y0 * Dthick
Y1 = Y1 * Dthick
Y2 = Y2 * Dthick
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF ( yld%normTypeIsTimesSigTot() ) THEN
Y0 = Y0 * Dthick * Total
Y1 = Y1 * Dthick * Total
Y2 = Y2 * Dthick * Total
@@ -551,9 +552,9 @@ C
IF (Kvthck.GT.0) THEN
C *** Add deriv of Y0 wrt Thickness; modify deriv of Y1 & Y2 wrt Thick
Kv = Kvthck - Nvadif - Ndasig
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Dbsigx(Kv) = Exp1*Cap + Dbsigx(Kv)
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Sn = Dthick*Total
IF (Sn.LT.0.3d0) THEN
Abc = Abcexp (-Sn, A, B, C, D, Ijklmn)
@@ -564,14 +565,14 @@ C *** Add deriv of Y0 wrt Thickness; modify deriv of Y1 & Y2 wrt Thick
Dbsigx(Kv) = F + Dbsigx(Kv)
IF (Kssmsc.NE.-2) Dbsigx(Kv) = Dbsigx(Kv) - Y1/Dthick
IF (Kssdbl.EQ. 1) Dbsigx(Kv) = Dbsigx(Kv) - Y2/Dthick
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Dbsigx(Kv)= Exp1*Cap*Total + Dbsigx(Kv)
END IF
END IF
C
IF (Sensin.NE.Zero .AND. Ksensc.GT.0) THEN
- IF (Kyield.EQ.0) Asensn = Asensn * Dthick
- IF (Kyield.EQ.2) Asensn = Asensn * Dthick * Total
+ IF (yld%normTypeIsSelfShielded()) Asensn = Asensn * Dthick
+ IF (yld%normTypeIsTimesSigTot()) Asensn = Asensn *Dthick *Total
C *** Add deriv of Y0 wrt neutron sensitivity multiplier
Kv = Ksensc - Nvadif - Ndasig
Dbsigx(Kv) = Asensn
@@ -586,6 +587,7 @@ C
* Dy2, Y0, Y1, Y2, Dthick, Total, Exp1, Dddyy0, Nn, Nd, Na)
use fixedi_m
use ifwrit_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*), Dy2(*),ds(*)
DATA One /1.0d0/
@@ -594,40 +596,40 @@ C
Iipar = Ipar + Na
C
C *** First, the derivative of Y0 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) +
* (Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total)
* + Dtotsi(Iipar,Nn)*Exp1*Cap/Total ) * Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) +
* Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total)
* + Dtotsi(Iipar,Nn)*Exp1*Cap/Total
IF (Iipar.EQ.1) Dddyy0 = Dddyy0 +
* Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total)
* + Dtotsi(Iipar,Nn)*Exp1*Cap/Total
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + (One-Exp1)*Dcapsi(Iipar,Nn)
* + Dtotsi(Iipar,Nn)*Exp1*Cap*Dthick
END IF
C
IF (Kssmsc.NE.-2) THEN
C *** Next, the derivative of Y1 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)*Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)*Total*Dthick +
* Y1*Dtotsi(Iipar,Nn)*Dthick
END IF
C
IF (Kssdbl.EQ.1) THEN
C *** Finally, the derivative of Y2 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar)*Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar)*Total*Dthick
* + Y2*Dtotsi(Iipar,Nn)*Dthick
END IF
diff --git a/sammy/src/par/mpar0.f b/sammy/src/par/mpar0.f
index afaa8d445..3650a3b7e 100644
--- a/sammy/src/par/mpar0.f
+++ b/sammy/src/par/mpar0.f
@@ -427,7 +427,9 @@ C
use fixedr_m
use ExpPars_common_m
use EndfData_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ character(len=80)::norm_string
DATA Zero /0.0d0/
real(kind=8) Ybeam
@@ -448,17 +450,15 @@ C
10700 FORMAT ('sigma.dat')
WRITE (12,10800) Lllmax-1, Numiso
10800 FORMAT (2I3)
- IF (Kyield.EQ.0) THEN
- WRITE (12,10900)
-10900 FORMAT ('Yield = Normalization type')
- ELSE IF (Kyield.EQ.1) THEN
- WRITE (12,11000)
-11000 FORMAT ('Cross section = Normalization type')
- ELSE IF (Kyield.EQ.2) THEN
- WRITE (12,11100)
-11100 FORMAT ('Other (1-e)sigma = Normalization type')
+ call yld%getNormString(norm_string)
+ IF ( yld%normTypeIsSelfShielded() ) THEN
+ WRITE (12,*) norm_string
+ ELSE IF ( yld%normTypeIsDivideByN() ) THEN
+ WRITE (12,*) norm_string
+ ELSE IF ( yld%normTypeIsTimesSigTot() ) THEN
+ WRITE (12,*) norm_string
ELSE
- STOP '[STOP in mpar0.f, Kyield is unknown]'
+ STOP '[STOP in mpar0.f, Kyield/( normaliz. type) is unknown]'
END IF
IF (Numiso.GT.0) THEN
WRITE (12,11200) (resParData%getMassForIsotope(I),I=1,Numiso)
diff --git a/sammy/src/ref/mcon.F b/sammy/src/ref/mcon.F
index 270828393..7a9852c99 100644
--- a/sammy/src/ref/mcon.F
+++ b/sammy/src/ref/mcon.F
@@ -19,6 +19,7 @@ C These are here to replace the common block
use zzzzz_common_m
use EndfData_common_m
use Samdat_0_M
+ use Observable_common_m
C this is not a common block, this contains functions
use mold4_m
IMPLICIT DOUBLE PRECISION (A-h,o-z)
@@ -31,6 +32,7 @@ CC ****************************************************
C
real(kind=8)::A(-1*(SAMMY_ARRAY_SIZE):(SAMMY_ARRAY_SIZE))
+ call yld%initialize()
call resParData%initialize()
call covData%initialize()
call covData%makeNewCovariance()
diff --git a/sammy/src/salmon/CMakeLists.txt b/sammy/src/salmon/CMakeLists.txt
index f80412a95..c26baf928 100644
--- a/sammy/src/salmon/CMakeLists.txt
+++ b/sammy/src/salmon/CMakeLists.txt
@@ -1,5 +1,8 @@
INCLUDE(SammyPackageSetup)
+# --------------------------------------
+# SAmmy Legacy Memory Obsolete Now
+# --------------------------------------
TRIBITS_PACKAGE(salmon)
@@ -14,6 +17,9 @@ SET(HEADERS ExperimentalParameters.h
MultipleScattering.h
interface/cpp/MultipleScatteringInterface.h
+
+ Yield.h
+ interface/cpp/YieldInterface.h
)
APPEND_SET(SALMON_SOURCES
@@ -27,6 +33,11 @@ APPEND_SET(SALMON_SOURCES
interface/cpp/MultipleScatteringInterface.cpp
interface/fortran/MultipleScattering_M.f90
interface/fortran/MultipleScattering_I.f90
+
+ Yield.cpp
+ interface/cpp/YieldInterface.cpp
+ interface/fortran/Yield_M.f90
+ interface/fortran/Yield_I.f90
)
diff --git a/sammy/src/salmon/ExperimentalParameters.cpp b/sammy/src/salmon/ExperimentalParameters.cpp
index b1f6a370c..5c0ed733e 100644
--- a/sammy/src/salmon/ExperimentalParameters.cpp
+++ b/sammy/src/salmon/ExperimentalParameters.cpp
@@ -8,6 +8,7 @@ namespace sammy{
radius = 0.0;
height = 0.0;
thickness = 0.0;
+ arealDensity = 0.0;
beamRadius = 0.0;
}
@@ -15,6 +16,7 @@ namespace sammy{
radius = orig.radius;
height = orig.height;
thickness = orig.thickness;
+ arealDensity = orig.arealDensity;
beamRadius = orig.beamRadius;
}
diff --git a/sammy/src/salmon/ExperimentalParameters.h b/sammy/src/salmon/ExperimentalParameters.h
index 59a60b524..17eed0c47 100644
--- a/sammy/src/salmon/ExperimentalParameters.h
+++ b/sammy/src/salmon/ExperimentalParameters.h
@@ -35,6 +35,13 @@ namespace sammy{
void setThickness(double t){
thickness = t;
}
+ /**
+ * Set the areal density
+ * @param ad the areal density of the sample [at/b]
+ */
+ void setArealDensity(double ad){
+ arealDensity = ad;
+ }
/**
* Set the beam radius
* @param br the beam radius [cm]
@@ -72,6 +79,12 @@ namespace sammy{
double getThickness() const{
return thickness;
}
+ /**
+ * Get the thickness of the sample [at/b]
+ */
+ double getArealDensity() const{
+ return arealDensity;
+ }
/**
* Get the beam radius [cm]
*/
@@ -106,6 +119,11 @@ namespace sammy{
******************************/
double thickness; /**< thickness [cm] */
+ /******************************
+ * The thickness of the sample in atom/b
+ ******************************/
+ double arealDensity; /**< thickness [atom/b] */
+
/******************************
* The radius of the beam striking the sample in cm
******************************/
diff --git a/sammy/src/salmon/MultipleScattering.cpp b/sammy/src/salmon/MultipleScattering.cpp
index 0257a6664..c83322264 100644
--- a/sammy/src/salmon/MultipleScattering.cpp
+++ b/sammy/src/salmon/MultipleScattering.cpp
@@ -5,6 +5,7 @@
namespace sammy{
MultipleScattering::MultipleScattering(){
+ numScatters = 0;
numTheta = 0;
numZHatIntegration = 0;
numBeamAreaIntegration = 0;
@@ -18,6 +19,7 @@ namespace sammy{
}
MultipleScattering::MultipleScattering(const MultipleScattering &orig){
+ numScatters = orig.numScatters;
numTheta = orig.numTheta;
numZHatIntegration = orig.numZHatIntegration;
numBeamAreaIntegration = orig.numBeamAreaIntegration;
diff --git a/sammy/src/salmon/MultipleScattering.h b/sammy/src/salmon/MultipleScattering.h
index 7d6d18ce5..36fd5b950 100644
--- a/sammy/src/salmon/MultipleScattering.h
+++ b/sammy/src/salmon/MultipleScattering.h
@@ -14,54 +14,59 @@ namespace sammy{
* Setters
******************************/
+ /**
+ * Set number of scatters considered before neutron absorption
+ * @param numScat the number of scatters considered before neutron absorption
+ */
+ void setNumScatters(int numScat) { numScatters = numScat;}
/**
* Set number of integration points over theta
- * @ param Ntheta number of integration points over theta
+ * @param Ntheta number of integration points over theta
*/
void setNumTheta(int Ntheta) { numTheta = Ntheta;}
/**
* Set number of integration points in the z-direction
- * @ param Ngausz number of integration points in the z-direction (4, 8, or 16 Default = 16)
+ * @param Ngausz number of integration points in the z-direction (4, 8, or 16 Default = 16)
*/
void setNumZHatIntegration(int Ngausz) { numZHatIntegration = Ngausz;}
/**
* Set number of integration points over cross section of beam
- * @ param Ngaus number of integration points over cross section of beam (4, 8, or 16 Default = 16)
+ * @param Ngaus number of integration points over cross section of beam (4, 8, or 16 Default = 16)
*/
void setNumBeamAreaIntegration(int Ngaus) { numBeamAreaIntegration = Ngaus;}
/**
* Set number of theta integration points near the end points
- * @ param Mtheta number of integration points theta integration near the end points
+ * @param Mtheta number of integration points theta integration near the end points
*/
void setNumThetaEndPoints(int Mtheta) { numThetaEndPoints = Mtheta;}
/**
* Set number of integration points for interpolation of sigma
- * @ param Nxtpt number of integration points for interpolation of sigma
+ * @param Nxtpt number of integration points for interpolation of sigma
*/
void setNumInterpSigma(int Nxtpt) { numInterpSigma = Nxtpt;}
/**
* Set number of integration points for interpolation over sigma prime (scattered)
- * @ param Mxtpt number of integration points for interpolation over sigma prime (scattered)
+ * @param Mxtpt number of integration points for interpolation over sigma prime (scattered)
*/
void setNumInterpSigmaPrime(int Mxtpt) { numInterpSigmaPrime = Mxtpt;}
/**
* Set minimum value of log of total cross section
- * @ param Xtmin minimum value of log of total cross section
+ * @param Xtmin minimum value of log of total cross section
*/
void setLogSigmaTotMin(double Xtmin) { logSigmaTotMin = Xtmin;}
/**
* Set maximum value of log of total cross section
- * @ param Xtmax maximum value of log of total cross section
+ * @param Xtmax maximum value of log of total cross section
*/
void setLogSigmaTotMax(double Xtmax) { logSigmaTotMax = Xtmax;}
/**
* Set number of integration points in theta integration near cos(theta) = 1
- * @ param Jtheta number of integration points in theta integration near cos(theta) = 1
+ * @param Jtheta number of integration points in theta integration near cos(theta) = 1
*/
void setNumThetaNearOne(int Jtheta) { numThetaNearOne = Jtheta;}
/**
* Set number of integration points in theta integration near cos(theta) = 0
- * @ param Ktheta number of integration points in theta integration near cos(theta) = 0
+ * @param Ktheta number of integration points in theta integration near cos(theta) = 0
*/
void setNumThetaNearZero(int Ktheta) { numThetaNearZero = Ktheta;}
@@ -69,6 +74,12 @@ namespace sammy{
Getters
******************************/
+ /**
+ * Get number of scatters considered before neutron absorption
+ */
+ int getNumScatters() const{
+ return numScatters;
+ }
/**
* Get number of integration points over theta
*/
@@ -134,6 +145,11 @@ namespace sammy{
private:
+ /**
+ * How many scatters to consider before neutron absorption?
+ */
+ int numScatters;
+
/******************************
* The number of points to integrate over angle theta (Default = 33)
******************************/
diff --git a/sammy/src/salmon/Yield.cpp b/sammy/src/salmon/Yield.cpp
new file mode 100644
index 000000000..88e4588db
--- /dev/null
+++ b/sammy/src/salmon/Yield.cpp
@@ -0,0 +1,32 @@
+#include "Yield.h"
+#include <iostream>
+
+namespace sammy{
+
+ Yield::Yield()
+ {
+ normString = "";
+ normType = normalization::NOTSET;
+ calcDerivs = false;
+ numParam = 0;
+ Y0 = 0.0;
+ Y1 = 0.0;
+ Y2 = 0.0;
+ totY = 0.0;
+ }
+
+ Yield::Yield(const Yield &orig)
+ {
+ normType = orig.normType;
+ calcDerivs = orig.calcDerivs;
+ numParam = orig.numParam;
+ Y0 = orig.Y0;
+ Y1 = orig.Y1;
+ Y2 = orig.Y2;
+ totY = orig.totY;
+ derivY0.insert( derivY0.begin(), orig.derivY0.begin(), orig.derivY0.end() );
+ derivY1.insert( derivY1.begin(), orig.derivY1.begin(), orig.derivY1.end() );
+ derivY2.insert( derivY2.begin(), orig.derivY2.begin(), orig.derivY2.end() );
+ derivTotY.insert( derivTotY.begin(), orig.derivTotY.begin(), orig.derivTotY.end() );
+ }
+}
\ No newline at end of file
diff --git a/sammy/src/salmon/Yield.h b/sammy/src/salmon/Yield.h
new file mode 100644
index 000000000..8eca77bfd
--- /dev/null
+++ b/sammy/src/salmon/Yield.h
@@ -0,0 +1,287 @@
+#ifndef YIELD_H
+#define YIELD_H
+
+#include <vector>
+#include <string>
+#include <memory>
+
+namespace sammy{
+
+ /** @brief Storage for the observed quantity of yield.
+
+ Store the values specific to yields (capture, fission, scattering...)
+ @date February 2020
+ */
+ class Yield{
+
+ public:
+
+
+ /******************************
+ * The type of normalization applied to the yield
+ ******************************/
+ enum normalization {SELFSHIELDED, DIVIDEBYN, TIMESSIGTOT, NOTSET} normType;
+
+ Yield();
+ Yield(const Yield& orig);
+ virtual ~Yield(){};
+
+ /******************************
+ * Setters
+ ******************************/
+
+ /**
+ * Set normalization type for yield to: self-shielded
+ */
+ void setNormTypeSelfShielded() {
+ normType = normalization::SELFSHIELDED;
+ normString = "Yield = Normalization type";
+ }
+
+ /**
+ * Set normalization type for yield to: divide by sample thickness \f$ N \f$
+ */
+ void setNormTypeDivideByN() {
+ normType = normalization::DIVIDEBYN;
+ normString = "Cross section = Normalization type";
+ }
+
+ /**
+ * Set normalization type for yield to: multiply by \f$ \sigma_t \f$
+ */
+ void setNormTypeTimesSigTot() {
+ normType = normalization::TIMESSIGTOT;
+ normString = "Other (1-e)sigma = Normalization type";
+ }
+
+ /**
+ * Set the boolean for whether to solve Bayes (calc. derivs.)
+ * @param sb boolean for whether we are solving Bayes eqtns.
+ */
+ void setCalcDerivs(bool cd) { calcDerivs = cd;}
+
+ /**
+ * Set the number of varied parameters (this should be temporary...we should take
+ * the info from a CovarianceData object or resonance parameter object)
+ */
+ void setNumTotalParam(int np){
+ numParam = np;
+ }
+
+ /**
+ * Set the primary yield
+ * @param y0 the primary yield
+ */
+ void setY0(double y0) { Y0 = y0;}
+
+ /**
+ * Set the secondary yield
+ * @param y1 the secondary yield
+ */
+ void setY1(double y1) { Y1 = y1;}
+
+ /**
+ * Set the tertiary yield
+ * @param y2 the tertiary yield
+ */
+ void setY2(double y2) { Y2 = y2;}
+
+ /**
+ * Set the total yield
+ * @param yt the total yield
+ */
+ void setTotY(double yt) { totY = yt; }
+
+ /**
+ * Set the derivative of primary yield w.r.t. parameter at indice iPar
+ * @param iPar index for the parameter of which deriv. was taken w.r.t.
+ * @param dy0 the primary yield
+ */
+ void setDerivY0(int iPar, double dy0){
+ if ((int)derivY0.size() < (iPar+1)){
+ derivY0.resize(iPar+1);
+ }
+ derivY0[iPar] = dy0;
+ }
+
+ /**
+ * Set the derivative of secondary yield w.r.t. parameter at indice iPar
+ * @param iPar index for the parameter of which deriv. was taken w.r.t.
+ * @param dy1 the secondary yield
+ */
+ void setDerivY1(int iPar, double dy1){
+ if ((int)derivY1.size() < (iPar+1)){
+ derivY1.resize(iPar+1);
+ }
+ derivY1[iPar] = dy1;
+ }
+
+ /**
+ * Set the derivative of tertiary yield w.r.t. parameter at indice iPar
+ * @param iPar index for the parameter of which deriv. was taken w.r.t.
+ * @param dy2 the tertiary yield
+ */
+ void setDerivY2(int iPar, double dy2){
+ if ((int)derivY2.size() < (iPar+1)){
+ derivY2.resize(iPar+1);
+ }
+ derivY2[iPar] = dy2;
+ }
+
+ /**
+ * Set the derivative of total yield w.r.t. parameter at indice iPar
+ * @param iPar index for the parameter of which deriv. was taken w.r.t.
+ * @param dyt the total yield
+ */
+ void setDerivTotY(int iPar, double dyt){
+ if ((int)derivTotY.size() < (iPar+1)){
+ derivTotY.resize(iPar+1);
+ }
+ derivTotY[iPar] = dyt;
+ }
+
+ /******************************
+ Getters
+ ******************************/
+
+ /**
+ * Get the normalization string for printing
+ */
+ std::string getNormString() const{ return normString; }
+
+ /**
+ * Get boolean whether normalization type for yield is self-shielded
+ */
+ bool normTypeIsSelfShielded() const{
+ return normType == normalization::SELFSHIELDED;
+ }
+
+ /**
+ * Get boolean whether normalization type for yield is divide by sample thickness \f$ N \f$
+ */
+ bool normTypeIsDivideByN() const{
+ return normType == normalization::DIVIDEBYN;
+ }
+
+ /**
+ * Get boolean whether normalization type for yield is to multiply by \f$ \sigma_t \f$
+ */
+ bool normTypeIsTimesSigTot() const{
+ return normType == normalization::TIMESSIGTOT;
+ }
+
+ /**
+ * Get boolean whether normalization type for yield has not been set
+ */
+ bool normTypeIsNotSet() const{
+ return normType == normalization::NOTSET;
+ }
+
+ /**
+ * Get decision on whether we solve Bayes (calc. derivs.)
+ */
+ bool getCalcDerivs() const{
+ return calcDerivs;
+ }
+
+ /**
+ * Get the total number of varied parameters (this should be temporary...we want to take
+ * the info from a CovarianceData object or resonance parameter object)
+ */
+ int getNumTotalParam() const{
+ return numParam;
+ }
+
+ /**
+ * Get the primary yield
+ */
+ double getY0() const{
+ return Y0;
+ }
+
+ /**
+ * Get the secondary yield
+ */
+ double getY1() const{
+ return Y1;
+ }
+
+ /**
+ * Get the tertiary yield
+ */
+ double getY2() const{
+ return Y2;
+ }
+
+ /**
+ * Get the total yield
+ */
+ double getTotY() const{
+ return totY;
+ }
+
+ /**
+ * Get the derivative of primary yield w.r.t. param. iPar
+ * @param iPar the indice of the parameter of which der. was taken w.r.t.
+ */
+ double getDerivY0(int iPar) const{
+ return derivY0[iPar];
+ }
+
+ /**
+ * Get the derivative of secondary yield w.r.t. param. iPar
+ * @param iPar the indice of the parameter of which der. was taken w.r.t.
+ */
+ double getDerivY1(int iPar) const{
+ return derivY1[iPar];
+ }
+
+ /**
+ * Get the derivative of tertiary yield w.r.t. param. iPar
+ * @param iPar the indice of the parameter of which der. was taken w.r.t.
+ */
+ double getDerivY2(int iPar) const{
+ return derivY2[iPar];
+ }
+
+ /**
+ * Get the derivative oftotal yield
+ */
+ double getDerivTotY(int iPar) const{
+ return derivTotY[iPar];
+ }
+
+ private:
+
+ /******************************
+ * String to print to files, stdout, etc.
+ *****************************/
+
+ std::string normString;
+
+ /******************************
+ * Should we calculate derivatives?
+ *****************************/
+
+ bool calcDerivs;
+
+ /******************************
+ * How many varied parameters are there?
+ *****************************/
+ int numParam;
+
+ /******************************
+ * The primary, secondary, tertiary, and total (sum) yields
+ ******************************/
+ double Y0, Y1, Y2, totY;
+
+ /******************************
+ * The derivatives of the primary, secondary, tertiary, and total (sum) yields
+ ******************************/
+ std::vector<double> derivY0, derivY1, derivY2, derivTotY;
+ };
+
+}
+
+
+#endif // YIELD_H
diff --git a/sammy/src/salmon/interface/cix/ExperimentalParameters.cpp2f.xml b/sammy/src/salmon/interface/cix/ExperimentalParameters.cpp2f.xml
index 4169973e6..40d773d09 100644
--- a/sammy/src/salmon/interface/cix/ExperimentalParameters.cpp2f.xml
+++ b/sammy/src/salmon/interface/cix/ExperimentalParameters.cpp2f.xml
@@ -13,6 +13,9 @@
<method name="setThickness">
<param name="t" type="double"/>
</method>
+ <method name="setArealDensity">
+ <param name="t" type="double"/>
+ </method>
<method name="setBeamRadius">
<param name="br" type="double"/>
</method>
@@ -23,6 +26,7 @@
<method name="getRadius" return_type="double"/>
<method name="getHeight" return_type="double"/>
<method name="getThickness" return_type="double"/>
+ <method name="getArealDensity" return_type="double"/>
<method name="getBeamRadius" return_type="double"/>
<method name="getBeamHeight" return_type="double"/>
diff --git a/sammy/src/salmon/interface/cix/MultipleScattering.cpp2f.xml b/sammy/src/salmon/interface/cix/MultipleScattering.cpp2f.xml
index ae7d1a398..6262f21a1 100644
--- a/sammy/src/salmon/interface/cix/MultipleScattering.cpp2f.xml
+++ b/sammy/src/salmon/interface/cix/MultipleScattering.cpp2f.xml
@@ -1,9 +1,12 @@
<generate name="MultipleScattering">
- <include_relative name="../MultipleScattering.h"/>
+ <include_relative name="../../MultipleScattering.h"/>
<using_namespace name="sammy"/>
<class name="MultipleScattering">
+ <method name="setNumScatters">
+ <param name="ns" type="int"/>
+ </method>
<method name="setNumTheta">
<param name="nt" type="int"/>
</method>
@@ -35,6 +38,7 @@
<param name="ntnearzero" type="int"/>
</method>
+ <method name="getNumScatters" return_type="int"/>
<method name="getNumTheta" return_type="int"/>
<method name="getNumZHatIntegration" return_type="int"/>
<method name="getNumBeamAreaIntegration" return_type="int"/>
@@ -47,4 +51,4 @@
<method name="getNumThetaNearZero" return_type="int"/>
</class>
-</generate>
\ No newline at end of file
+</generate>
diff --git a/sammy/src/salmon/interface/cix/Yield.cpp2f.xml b/sammy/src/salmon/interface/cix/Yield.cpp2f.xml
new file mode 100644
index 000000000..1fe843295
--- /dev/null
+++ b/sammy/src/salmon/interface/cix/Yield.cpp2f.xml
@@ -0,0 +1,70 @@
+<generate name="Yield">
+ <include_relative name="../../Yield.h"/>
+ <using_namespace name="sammy"/>
+
+ <class name="Yield">
+
+ <method name="setNormTypeSelfShielded"/>
+ <method name="setNormTypeDivideByN"/>
+ <method name="setNormTypeTimesSigTot"/>
+ <method name="setCalcDerivs">
+ <param name="sb" type="bool"/>
+ </method>
+ <method name="setNumTotalParam">
+ <param name="np" type="int"/>
+ </method>
+ <method name="setY0">
+ <param name="y0" type="double"/>
+ </method>
+ <method name="setY1">
+ <param name="y1" type="double"/>
+ </method>
+ <method name="setY2">
+ <param name="y2" type="double"/>
+ </method>
+ <method name="setTotY">
+ <param name="toty" type="double"/>
+ </method>
+ <method name="setDerivY0">
+ <param name="iPar" type="int" offset="-1"/>
+ <param name="dy0" type="double"/>
+ </method>
+ <method name="setDerivY1">
+ <param name="iPar" type="int" offset="-1"/>
+ <param name="dy1" type="double"/>
+ </method>
+ <method name="setDerivY2">
+ <param name="iPar" type="int" offset="-1"/>
+ <param name="dy2" type="double"/>
+ </method>
+ <method name="setDerivTotY">
+ <param name="iPar" type="int" offset="-1"/>
+ <param name="dy2" type="double"/>
+ </method>
+
+ <method name="getNormString" return_type="string"/>
+ <method name="normTypeIsSelfShielded" return_type="bool"/>
+ <method name="normTypeIsDivideByN" return_type="bool"/>
+ <method name="normTypeIsTimesSigTot" return_type="bool"/>
+ <method name="getCalcDerivs" return_type="bool"/>
+ <method name="getNumTotalParam" return_type="int"/>
+ <method name="getY0" return_type="double"/>
+ <method name="getY1" return_type="double"/>
+ <method name="getY2" return_type="double"/>
+ <method name="getTotY" return_type="double"/>
+
+ <method name="getDerivY0" return_type="double">
+ <param name="iPar" type="int" offset="-1"/>
+ </method>
+ <method name="getDerivY1" return_type="double">
+ <param name="iPar" type="int" offset="-1"/>
+ </method>
+ <method name="getDerivY2" return_type="double">
+ <param name="iPar" type="int" offset="-1"/>
+ </method>
+ <method name="getDerivTotY" return_type="double">
+ <param name="iPar" type="int" offset="-1"/>
+ </method>
+
+ </class>
+</generate>
\ No newline at end of file
diff --git a/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.cpp b/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.cpp
index c3eacbacf..96e23ad4c 100644
--- a/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.cpp
+++ b/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.cpp
@@ -2,7 +2,7 @@
* This file has been dynamically generated by Class Interface Xml (CIX)
* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
* If changes need to occur, modify the appropriate CIX xml file
-* Date Generated: Fri Dec 20 13:08:25 EST 2019
+* Date Generated: Tue Mar 03 09:03:32 EST 2020
* If any issues are experiences with this generated file that cannot be fixed
* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
*/
@@ -24,6 +24,11 @@ void SampleDimensions_setThickness(void * SampleDimensions_ptr,double * t)
((SampleDimensions*)SampleDimensions_ptr)->setThickness(*t);
}
+void SampleDimensions_setArealDensity(void * SampleDimensions_ptr,double * t)
+{
+ ((SampleDimensions*)SampleDimensions_ptr)->setArealDensity(*t);
+}
+
void SampleDimensions_setBeamRadius(void * SampleDimensions_ptr,double * br)
{
((SampleDimensions*)SampleDimensions_ptr)->setBeamRadius(*br);
@@ -49,6 +54,11 @@ double SampleDimensions_getThickness(void * SampleDimensions_ptr)
return ((SampleDimensions*)SampleDimensions_ptr)->getThickness();
}
+double SampleDimensions_getArealDensity(void * SampleDimensions_ptr)
+{
+ return ((SampleDimensions*)SampleDimensions_ptr)->getArealDensity();
+}
+
double SampleDimensions_getBeamRadius(void * SampleDimensions_ptr)
{
return ((SampleDimensions*)SampleDimensions_ptr)->getBeamRadius();
diff --git a/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.h b/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.h
index a33085c0e..43d5407dc 100644
--- a/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.h
+++ b/sammy/src/salmon/interface/cpp/ExperimentalParametersInterface.h
@@ -2,7 +2,7 @@
* This file has been dynamically generated by Class Interface Xml (CIX)
* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
* If changes need to occur, modify the appropriate CIX xml file
-* Date Generated: Fri Dec 20 13:08:25 EST 2019
+* Date Generated: Tue Mar 03 09:03:32 EST 2020
* If any issues are experiences with this generated file that cannot be fixed
* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
*/
@@ -16,11 +16,13 @@ extern "C" {
void SampleDimensions_setRadius(void * SampleDimensions_ptr,double * r);
void SampleDimensions_setHeight(void * SampleDimensions_ptr,double * h);
void SampleDimensions_setThickness(void * SampleDimensions_ptr,double * t);
+void SampleDimensions_setArealDensity(void * SampleDimensions_ptr,double * t);
void SampleDimensions_setBeamRadius(void * SampleDimensions_ptr,double * br);
void SampleDimensions_setBeamHeight(void * SampleDimensions_ptr,double * br);
double SampleDimensions_getRadius(void * SampleDimensions_ptr);
double SampleDimensions_getHeight(void * SampleDimensions_ptr);
double SampleDimensions_getThickness(void * SampleDimensions_ptr);
+double SampleDimensions_getArealDensity(void * SampleDimensions_ptr);
double SampleDimensions_getBeamRadius(void * SampleDimensions_ptr);
double SampleDimensions_getBeamHeight(void * SampleDimensions_ptr);
void* SampleDimensions_initialize();
diff --git a/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.cpp b/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.cpp
index 1e88a180d..65c95b66f 100644
--- a/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.cpp
+++ b/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.cpp
@@ -2,13 +2,18 @@
* This file has been dynamically generated by Class Interface Xml (CIX)
* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
* If changes need to occur, modify the appropriate CIX xml file
-* Date Generated: Mon Dec 09 09:48:43 EST 2019
+* Date Generated: Mon Mar 02 10:33:41 EST 2020
* If any issues are experiences with this generated file that cannot be fixed
* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
*/
#include <string.h>
#include "MultipleScatteringInterface.h"
using namespace sammy;
+void MultipleScattering_setNumScatters(void * MultipleScattering_ptr,int * ns)
+{
+ ((MultipleScattering*)MultipleScattering_ptr)->setNumScatters(*ns);
+}
+
void MultipleScattering_setNumTheta(void * MultipleScattering_ptr,int * nt)
{
((MultipleScattering*)MultipleScattering_ptr)->setNumTheta(*nt);
@@ -59,6 +64,11 @@ void MultipleScattering_setNumThetaNearZero(void * MultipleScattering_ptr,int *
((MultipleScattering*)MultipleScattering_ptr)->setNumThetaNearZero(*ntnearzero);
}
+int MultipleScattering_getNumScatters(void * MultipleScattering_ptr)
+{
+ return ((MultipleScattering*)MultipleScattering_ptr)->getNumScatters();
+}
+
int MultipleScattering_getNumTheta(void * MultipleScattering_ptr)
{
return ((MultipleScattering*)MultipleScattering_ptr)->getNumTheta();
diff --git a/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.h b/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.h
index 269bd916d..778c51f06 100644
--- a/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.h
+++ b/sammy/src/salmon/interface/cpp/MultipleScatteringInterface.h
@@ -2,7 +2,7 @@
* This file has been dynamically generated by Class Interface Xml (CIX)
* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
* If changes need to occur, modify the appropriate CIX xml file
-* Date Generated: Mon Dec 09 09:48:43 EST 2019
+* Date Generated: Mon Mar 02 10:33:41 EST 2020
* If any issues are experiences with this generated file that cannot be fixed
* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
*/
@@ -13,6 +13,7 @@ using namespace sammy;
#ifdef __cplusplus
extern "C" {
#endif
+void MultipleScattering_setNumScatters(void * MultipleScattering_ptr,int * ns);
void MultipleScattering_setNumTheta(void * MultipleScattering_ptr,int * nt);
void MultipleScattering_setNumZHatIntegration(void * MultipleScattering_ptr,int * nz);
void MultipleScattering_setNumBeamAreaIntegration(void * MultipleScattering_ptr,int * nba);
@@ -23,6 +24,7 @@ void MultipleScattering_setLogSigmaTotMin(void * MultipleScattering_ptr,double *
void MultipleScattering_setLogSigmaTotMax(void * MultipleScattering_ptr,double * nsigmax);
void MultipleScattering_setNumThetaNearOne(void * MultipleScattering_ptr,int * ntnearone);
void MultipleScattering_setNumThetaNearZero(void * MultipleScattering_ptr,int * ntnearzero);
+int MultipleScattering_getNumScatters(void * MultipleScattering_ptr);
int MultipleScattering_getNumTheta(void * MultipleScattering_ptr);
int MultipleScattering_getNumZHatIntegration(void * MultipleScattering_ptr);
int MultipleScattering_getNumBeamAreaIntegration(void * MultipleScattering_ptr);
diff --git a/sammy/src/salmon/interface/cpp/YieldInterface.cpp b/sammy/src/salmon/interface/cpp/YieldInterface.cpp
new file mode 100644
index 000000000..e118efb1b
--- /dev/null
+++ b/sammy/src/salmon/interface/cpp/YieldInterface.cpp
@@ -0,0 +1,163 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX)
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Wed Apr 08 13:29:36 EDT 2020
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#include <string.h>
+#include "YieldInterface.h"
+using namespace sammy;
+void Yield_setNormTypeSelfShielded(void * Yield_ptr)
+{
+ ((Yield*)Yield_ptr)->setNormTypeSelfShielded();
+}
+
+void Yield_setNormTypeDivideByN(void * Yield_ptr)
+{
+ ((Yield*)Yield_ptr)->setNormTypeDivideByN();
+}
+
+void Yield_setNormTypeTimesSigTot(void * Yield_ptr)
+{
+ ((Yield*)Yield_ptr)->setNormTypeTimesSigTot();
+}
+
+void Yield_setCalcDerivs(void * Yield_ptr,bool * sb)
+{
+ ((Yield*)Yield_ptr)->setCalcDerivs(*sb);
+}
+
+void Yield_setNumTotalParam(void * Yield_ptr,int * np)
+{
+ ((Yield*)Yield_ptr)->setNumTotalParam(*np);
+}
+
+void Yield_setY0(void * Yield_ptr,double * y0)
+{
+ ((Yield*)Yield_ptr)->setY0(*y0);
+}
+
+void Yield_setY1(void * Yield_ptr,double * y1)
+{
+ ((Yield*)Yield_ptr)->setY1(*y1);
+}
+
+void Yield_setY2(void * Yield_ptr,double * y2)
+{
+ ((Yield*)Yield_ptr)->setY2(*y2);
+}
+
+void Yield_setTotY(void * Yield_ptr,double * toty)
+{
+ ((Yield*)Yield_ptr)->setTotY(*toty);
+}
+
+void Yield_setDerivY0(void * Yield_ptr,int * iPar,double * dy0)
+{
+ ((Yield*)Yield_ptr)->setDerivY0(*iPar,*dy0);
+}
+
+void Yield_setDerivY1(void * Yield_ptr,int * iPar,double * dy1)
+{
+ ((Yield*)Yield_ptr)->setDerivY1(*iPar,*dy1);
+}
+
+void Yield_setDerivY2(void * Yield_ptr,int * iPar,double * dy2)
+{
+ ((Yield*)Yield_ptr)->setDerivY2(*iPar,*dy2);
+}
+
+void Yield_setDerivTotY(void * Yield_ptr,int * iPar,double * dy2)
+{
+ ((Yield*)Yield_ptr)->setDerivTotY(*iPar,*dy2);
+}
+
+void Yield_getNormString(void * Yield_ptr, char * string2Populate, int * lengthOfString2Populate)
+{
+ const std::string& cix_string = ((Yield*)Yield_ptr)->getNormString();
+ int t_length = *lengthOfString2Populate < static_cast<int>(cix_string.length())
+ ? *lengthOfString2Populate : cix_string.length();
+ strncpy(string2Populate, cix_string.c_str(), t_length);
+ for(int cix_i = t_length; cix_i < *lengthOfString2Populate; cix_i++){
+ string2Populate[cix_i] = ' ';
+ }
+
+}
+
+bool Yield_normTypeIsSelfShielded(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->normTypeIsSelfShielded();
+}
+
+bool Yield_normTypeIsDivideByN(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->normTypeIsDivideByN();
+}
+
+bool Yield_normTypeIsTimesSigTot(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->normTypeIsTimesSigTot();
+}
+
+bool Yield_getCalcDerivs(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getCalcDerivs();
+}
+
+int Yield_getNumTotalParam(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getNumTotalParam();
+}
+
+double Yield_getY0(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getY0();
+}
+
+double Yield_getY1(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getY1();
+}
+
+double Yield_getY2(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getY2();
+}
+
+double Yield_getTotY(void * Yield_ptr)
+{
+ return ((Yield*)Yield_ptr)->getTotY();
+}
+
+double Yield_getDerivY0(void * Yield_ptr,int * iPar)
+{
+ return ((Yield*)Yield_ptr)->getDerivY0(*iPar);
+}
+
+double Yield_getDerivY1(void * Yield_ptr,int * iPar)
+{
+ return ((Yield*)Yield_ptr)->getDerivY1(*iPar);
+}
+
+double Yield_getDerivY2(void * Yield_ptr,int * iPar)
+{
+ return ((Yield*)Yield_ptr)->getDerivY2(*iPar);
+}
+
+double Yield_getDerivTotY(void * Yield_ptr,int * iPar)
+{
+ return ((Yield*)Yield_ptr)->getDerivTotY(*iPar);
+}
+
+void* Yield_initialize()
+{
+ return new Yield();
+}
+
+void Yield_destroy(void * Yield_ptr)
+{
+ delete (Yield*)Yield_ptr;
+}
+
diff --git a/sammy/src/salmon/interface/cpp/YieldInterface.h b/sammy/src/salmon/interface/cpp/YieldInterface.h
new file mode 100644
index 000000000..64065a15a
--- /dev/null
+++ b/sammy/src/salmon/interface/cpp/YieldInterface.h
@@ -0,0 +1,48 @@
+/*!
+* This file has been dynamically generated by Class Interface Xml (CIX)
+* DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+* If changes need to occur, modify the appropriate CIX xml file
+* Date Generated: Wed Apr 08 13:29:36 EDT 2020
+* If any issues are experiences with this generated file that cannot be fixed
+* with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+*/
+#ifndef YIELDINTERFACE_H
+#define YIELDINTERFACE_H
+#include "../../Yield.h"
+using namespace sammy;
+#ifdef __cplusplus
+extern "C" {
+#endif
+void Yield_setNormTypeSelfShielded(void * Yield_ptr);
+void Yield_setNormTypeDivideByN(void * Yield_ptr);
+void Yield_setNormTypeTimesSigTot(void * Yield_ptr);
+void Yield_setCalcDerivs(void * Yield_ptr,bool * sb);
+void Yield_setNumTotalParam(void * Yield_ptr,int * np);
+void Yield_setY0(void * Yield_ptr,double * y0);
+void Yield_setY1(void * Yield_ptr,double * y1);
+void Yield_setY2(void * Yield_ptr,double * y2);
+void Yield_setTotY(void * Yield_ptr,double * toty);
+void Yield_setDerivY0(void * Yield_ptr,int * iPar,double * dy0);
+void Yield_setDerivY1(void * Yield_ptr,int * iPar,double * dy1);
+void Yield_setDerivY2(void * Yield_ptr,int * iPar,double * dy2);
+void Yield_setDerivTotY(void * Yield_ptr,int * iPar,double * dy2);
+void Yield_getNormString(void * Yield_ptr, char * string2Populate, int * lengthOfString2Populate);
+bool Yield_normTypeIsSelfShielded(void * Yield_ptr);
+bool Yield_normTypeIsDivideByN(void * Yield_ptr);
+bool Yield_normTypeIsTimesSigTot(void * Yield_ptr);
+bool Yield_getCalcDerivs(void * Yield_ptr);
+int Yield_getNumTotalParam(void * Yield_ptr);
+double Yield_getY0(void * Yield_ptr);
+double Yield_getY1(void * Yield_ptr);
+double Yield_getY2(void * Yield_ptr);
+double Yield_getTotY(void * Yield_ptr);
+double Yield_getDerivY0(void * Yield_ptr,int * iPar);
+double Yield_getDerivY1(void * Yield_ptr,int * iPar);
+double Yield_getDerivY2(void * Yield_ptr,int * iPar);
+double Yield_getDerivTotY(void * Yield_ptr,int * iPar);
+void* Yield_initialize();
+void Yield_destroy(void * Yield_ptr);
+#ifdef __cplusplus
+}
+#endif
+#endif /* YIELDINTERFACE_H */
diff --git a/sammy/src/salmon/interface/fortran/ExperimentalParameters_I.f90 b/sammy/src/salmon/interface/fortran/ExperimentalParameters_I.f90
index 060172f10..f51a57315 100644
--- a/sammy/src/salmon/interface/fortran/ExperimentalParameters_I.f90
+++ b/sammy/src/salmon/interface/fortran/ExperimentalParameters_I.f90
@@ -2,7 +2,7 @@
!! This file has been dynamically generated by Class Interface Xml (CIX)
!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
!! If changes need to occur, modify the appropriate CIX xml file
-!! Date Generated: Fri Dec 20 13:08:25 EST 2019
+!! Date Generated: Tue Mar 03 09:03:32 EST 2020
!! If any issues are experiences with this generated file that cannot be fixed
!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
!!/
@@ -27,6 +27,12 @@ subroutine f_SampleDimensions_setThickness(SampleDimensions_ptr, t ) BIND(C,name
type(C_PTR), value :: SampleDimensions_ptr;
real(C_DOUBLE) :: t;
end subroutine
+subroutine f_SampleDimensions_setArealDensity(SampleDimensions_ptr, t ) BIND(C,name="SampleDimensions_setArealDensity")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: SampleDimensions_ptr;
+ real(C_DOUBLE) :: t;
+end subroutine
subroutine f_SampleDimensions_setBeamRadius(SampleDimensions_ptr, br ) BIND(C,name="SampleDimensions_setBeamRadius")
use,intrinsic :: ISO_C_BINDING
implicit none
@@ -54,6 +60,11 @@ real(C_DOUBLE) function f_SampleDimensions_getThickness(SampleDimensions_ptr ) B
implicit none
type(C_PTR), value :: SampleDimensions_ptr;
end function
+real(C_DOUBLE) function f_SampleDimensions_getArealDensity(SampleDimensions_ptr ) BIND(C,name="SampleDimensions_getArealDensity")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: SampleDimensions_ptr;
+end function
real(C_DOUBLE) function f_SampleDimensions_getBeamRadius(SampleDimensions_ptr ) BIND(C,name="SampleDimensions_getBeamRadius")
use,intrinsic :: ISO_C_BINDING
implicit none
diff --git a/sammy/src/salmon/interface/fortran/ExperimentalParameters_M.f90 b/sammy/src/salmon/interface/fortran/ExperimentalParameters_M.f90
index b24a121f3..b608ed5de 100644
--- a/sammy/src/salmon/interface/fortran/ExperimentalParameters_M.f90
+++ b/sammy/src/salmon/interface/fortran/ExperimentalParameters_M.f90
@@ -2,7 +2,7 @@
!! This file has been dynamically generated by Class Interface Xml (CIX)
!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
!! If changes need to occur, modify the appropriate CIX xml file
-!! Date Generated: Fri Dec 20 13:08:25 EST 2019
+!! Date Generated: Tue Mar 03 09:03:32 EST 2020
!! If any issues are experiences with this generated file that cannot be fixed
!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
!!/
@@ -15,11 +15,13 @@ type SampleDimensions
procedure, pass(this) :: setRadius => SampleDimensions_setRadius
procedure, pass(this) :: setHeight => SampleDimensions_setHeight
procedure, pass(this) :: setThickness => SampleDimensions_setThickness
+ procedure, pass(this) :: setArealDensity => SampleDimensions_setArealDensity
procedure, pass(this) :: setBeamRadius => SampleDimensions_setBeamRadius
procedure, pass(this) :: setBeamHeight => SampleDimensions_setBeamHeight
procedure, pass(this) :: getRadius => SampleDimensions_getRadius
procedure, pass(this) :: getHeight => SampleDimensions_getHeight
procedure, pass(this) :: getThickness => SampleDimensions_getThickness
+ procedure, pass(this) :: getArealDensity => SampleDimensions_getArealDensity
procedure, pass(this) :: getBeamRadius => SampleDimensions_getBeamRadius
procedure, pass(this) :: getBeamHeight => SampleDimensions_getBeamHeight
procedure, pass(this) :: initialize => SampleDimensions_initialize
@@ -54,6 +56,12 @@ subroutine SampleDimensions_setThickness(this, t)
real(C_DOUBLE)::t
call f_SampleDimensions_setThickness(this%instance_ptr, t)
end subroutine
+subroutine SampleDimensions_setArealDensity(this, t)
+ implicit none
+ class(SampleDimensions)::this
+ real(C_DOUBLE)::t
+ call f_SampleDimensions_setArealDensity(this%instance_ptr, t)
+end subroutine
subroutine SampleDimensions_setBeamRadius(this, br)
implicit none
class(SampleDimensions)::this
@@ -84,6 +92,12 @@ function SampleDimensions_getThickness(this) result(result2Return)
real(C_DOUBLE):: result2Return
result2Return=f_SampleDimensions_getThickness(this%instance_ptr)
end function
+function SampleDimensions_getArealDensity(this) result(result2Return)
+ implicit none
+ class(SampleDimensions)::this
+ real(C_DOUBLE):: result2Return
+ result2Return=f_SampleDimensions_getArealDensity(this%instance_ptr)
+end function
function SampleDimensions_getBeamRadius(this) result(result2Return)
implicit none
class(SampleDimensions)::this
diff --git a/sammy/src/salmon/interface/fortran/MultipleScattering_I.f90 b/sammy/src/salmon/interface/fortran/MultipleScattering_I.f90
index f25f04dc2..5422e5029 100644
--- a/sammy/src/salmon/interface/fortran/MultipleScattering_I.f90
+++ b/sammy/src/salmon/interface/fortran/MultipleScattering_I.f90
@@ -2,13 +2,19 @@
!! This file has been dynamically generated by Class Interface Xml (CIX)
!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
!! If changes need to occur, modify the appropriate CIX xml file
-!! Date Generated: Mon Dec 09 09:48:43 EST 2019
+!! Date Generated: Mon Mar 02 10:33:41 EST 2020
!! If any issues are experiences with this generated file that cannot be fixed
!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
!!/
module MultipleScattering_I
use, intrinsic :: ISO_C_BINDING
interface
+subroutine f_MultipleScattering_setNumScatters(MultipleScattering_ptr, ns ) BIND(C,name="MultipleScattering_setNumScatters")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: MultipleScattering_ptr;
+ integer(C_INT) :: ns;
+end subroutine
subroutine f_MultipleScattering_setNumTheta(MultipleScattering_ptr, nt ) BIND(C,name="MultipleScattering_setNumTheta")
use,intrinsic :: ISO_C_BINDING
implicit none
@@ -69,6 +75,11 @@ subroutine f_MultipleScattering_setNumThetaNearZero(MultipleScattering_ptr, ntne
type(C_PTR), value :: MultipleScattering_ptr;
integer(C_INT) :: ntnearzero;
end subroutine
+integer(C_INT) function f_MultipleScattering_getNumScatters(MultipleScattering_ptr ) BIND(C,name="MultipleScattering_getNumScatters")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: MultipleScattering_ptr;
+end function
integer(C_INT) function f_MultipleScattering_getNumTheta(MultipleScattering_ptr ) BIND(C,name="MultipleScattering_getNumTheta")
use,intrinsic :: ISO_C_BINDING
implicit none
diff --git a/sammy/src/salmon/interface/fortran/MultipleScattering_M.f90 b/sammy/src/salmon/interface/fortran/MultipleScattering_M.f90
index aa10aee1d..d6ed8dd22 100644
--- a/sammy/src/salmon/interface/fortran/MultipleScattering_M.f90
+++ b/sammy/src/salmon/interface/fortran/MultipleScattering_M.f90
@@ -2,7 +2,7 @@
!! This file has been dynamically generated by Class Interface Xml (CIX)
!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
!! If changes need to occur, modify the appropriate CIX xml file
-!! Date Generated: Mon Dec 09 09:48:43 EST 2019
+!! Date Generated: Mon Mar 02 10:33:41 EST 2020
!! If any issues are experiences with this generated file that cannot be fixed
!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
!!/
@@ -12,6 +12,7 @@ use MultipleScattering_I
type MultipleScattering
type(C_PTR) :: instance_ptr=C_NULL_PTR
contains
+ procedure, pass(this) :: setNumScatters => MultipleScattering_setNumScatters
procedure, pass(this) :: setNumTheta => MultipleScattering_setNumTheta
procedure, pass(this) :: setNumZHatIntegration => MultipleScattering_setNumZHatIntegration
procedure, pass(this) :: setNumBeamAreaIntegration => MultipleScattering_setNumBeamAreaIntegration
@@ -22,6 +23,7 @@ type MultipleScattering
procedure, pass(this) :: setLogSigmaTotMax => MultipleScattering_setLogSigmaTotMax
procedure, pass(this) :: setNumThetaNearOne => MultipleScattering_setNumThetaNearOne
procedure, pass(this) :: setNumThetaNearZero => MultipleScattering_setNumThetaNearZero
+ procedure, pass(this) :: getNumScatters => MultipleScattering_getNumScatters
procedure, pass(this) :: getNumTheta => MultipleScattering_getNumTheta
procedure, pass(this) :: getNumZHatIntegration => MultipleScattering_getNumZHatIntegration
procedure, pass(this) :: getNumBeamAreaIntegration => MultipleScattering_getNumBeamAreaIntegration
@@ -36,6 +38,12 @@ type MultipleScattering
procedure, pass(this) :: destroy => MultipleScattering_destroy
end type MultipleScattering
contains
+subroutine MultipleScattering_setNumScatters(this, ns)
+ implicit none
+ class(MultipleScattering)::this
+ integer(C_INT)::ns
+ call f_MultipleScattering_setNumScatters(this%instance_ptr, ns)
+end subroutine
subroutine MultipleScattering_setNumTheta(this, nt)
implicit none
class(MultipleScattering)::this
@@ -96,6 +104,12 @@ subroutine MultipleScattering_setNumThetaNearZero(this, ntnearzero)
integer(C_INT)::ntnearzero
call f_MultipleScattering_setNumThetaNearZero(this%instance_ptr, ntnearzero)
end subroutine
+function MultipleScattering_getNumScatters(this) result(result2Return)
+ implicit none
+ class(MultipleScattering)::this
+ integer(C_INT):: result2Return
+ result2Return=f_MultipleScattering_getNumScatters(this%instance_ptr)
+end function
function MultipleScattering_getNumTheta(this) result(result2Return)
implicit none
class(MultipleScattering)::this
diff --git a/sammy/src/salmon/interface/fortran/Yield_I.f90 b/sammy/src/salmon/interface/fortran/Yield_I.f90
new file mode 100644
index 000000000..695d2115f
--- /dev/null
+++ b/sammy/src/salmon/interface/fortran/Yield_I.f90
@@ -0,0 +1,178 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX)
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Wed Apr 08 13:29:36 EDT 2020
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module Yield_I
+use, intrinsic :: ISO_C_BINDING
+interface
+subroutine f_Yield_setNormTypeSelfShielded(Yield_ptr ) BIND(C,name="Yield_setNormTypeSelfShielded")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end subroutine
+subroutine f_Yield_setNormTypeDivideByN(Yield_ptr ) BIND(C,name="Yield_setNormTypeDivideByN")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end subroutine
+subroutine f_Yield_setNormTypeTimesSigTot(Yield_ptr ) BIND(C,name="Yield_setNormTypeTimesSigTot")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end subroutine
+subroutine f_Yield_setCalcDerivs(Yield_ptr, sb ) BIND(C,name="Yield_setCalcDerivs")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ logical(C_BOOL) :: sb;
+end subroutine
+subroutine f_Yield_setNumTotalParam(Yield_ptr, np ) BIND(C,name="Yield_setNumTotalParam")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: np;
+end subroutine
+subroutine f_Yield_setY0(Yield_ptr, y0 ) BIND(C,name="Yield_setY0")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ real(C_DOUBLE) :: y0;
+end subroutine
+subroutine f_Yield_setY1(Yield_ptr, y1 ) BIND(C,name="Yield_setY1")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ real(C_DOUBLE) :: y1;
+end subroutine
+subroutine f_Yield_setY2(Yield_ptr, y2 ) BIND(C,name="Yield_setY2")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ real(C_DOUBLE) :: y2;
+end subroutine
+subroutine f_Yield_setTotY(Yield_ptr, toty ) BIND(C,name="Yield_setTotY")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ real(C_DOUBLE) :: toty;
+end subroutine
+subroutine f_Yield_setDerivY0(Yield_ptr, iPar,dy0 ) BIND(C,name="Yield_setDerivY0")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+ real(C_DOUBLE) :: dy0;
+end subroutine
+subroutine f_Yield_setDerivY1(Yield_ptr, iPar,dy1 ) BIND(C,name="Yield_setDerivY1")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+ real(C_DOUBLE) :: dy1;
+end subroutine
+subroutine f_Yield_setDerivY2(Yield_ptr, iPar,dy2 ) BIND(C,name="Yield_setDerivY2")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+ real(C_DOUBLE) :: dy2;
+end subroutine
+subroutine f_Yield_setDerivTotY(Yield_ptr, iPar,dy2 ) BIND(C,name="Yield_setDerivTotY")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+ real(C_DOUBLE) :: dy2;
+end subroutine
+subroutine f_Yield_getNormString(Yield_ptr, string2Populate, lengthOfString2Populate) BIND(C,name="Yield_getNormString")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ character(C_CHAR) :: string2Populate;
+ integer(C_INT) :: lengthOfString2Populate;
+end subroutine
+logical(C_BOOL) function f_Yield_normTypeIsSelfShielded(Yield_ptr ) BIND(C,name="Yield_normTypeIsSelfShielded")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+logical(C_BOOL) function f_Yield_normTypeIsDivideByN(Yield_ptr ) BIND(C,name="Yield_normTypeIsDivideByN")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+logical(C_BOOL) function f_Yield_normTypeIsTimesSigTot(Yield_ptr ) BIND(C,name="Yield_normTypeIsTimesSigTot")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+logical(C_BOOL) function f_Yield_getCalcDerivs(Yield_ptr ) BIND(C,name="Yield_getCalcDerivs")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+integer(C_INT) function f_Yield_getNumTotalParam(Yield_ptr ) BIND(C,name="Yield_getNumTotalParam")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+real(C_DOUBLE) function f_Yield_getY0(Yield_ptr ) BIND(C,name="Yield_getY0")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+real(C_DOUBLE) function f_Yield_getY1(Yield_ptr ) BIND(C,name="Yield_getY1")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+real(C_DOUBLE) function f_Yield_getY2(Yield_ptr ) BIND(C,name="Yield_getY2")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+real(C_DOUBLE) function f_Yield_getTotY(Yield_ptr ) BIND(C,name="Yield_getTotY")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+end function
+real(C_DOUBLE) function f_Yield_getDerivY0(Yield_ptr, iPar ) BIND(C,name="Yield_getDerivY0")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+end function
+real(C_DOUBLE) function f_Yield_getDerivY1(Yield_ptr, iPar ) BIND(C,name="Yield_getDerivY1")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+end function
+real(C_DOUBLE) function f_Yield_getDerivY2(Yield_ptr, iPar ) BIND(C,name="Yield_getDerivY2")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+end function
+real(C_DOUBLE) function f_Yield_getDerivTotY(Yield_ptr, iPar ) BIND(C,name="Yield_getDerivTotY")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: Yield_ptr;
+ integer(C_INT) :: iPar;
+end function
+type(C_PTR) function f_Yield_initialize( )BIND(C,name="Yield_initialize")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR) :: Yield_ptr;
+end function
+subroutine f_Yield_destroy(this) BIND(C,name="Yield_destroy")
+ use,intrinsic :: ISO_C_BINDING
+ implicit none
+ type(C_PTR), value :: this;
+end subroutine
+end interface
+end module Yield_I
diff --git a/sammy/src/salmon/interface/fortran/Yield_M.f90 b/sammy/src/salmon/interface/fortran/Yield_M.f90
new file mode 100644
index 000000000..30eded3bb
--- /dev/null
+++ b/sammy/src/salmon/interface/fortran/Yield_M.f90
@@ -0,0 +1,226 @@
+!>
+!! This file has been dynamically generated by Class Interface Xml (CIX)
+!! DO NOT MODIFY THIS FILE -- CHANGES WILL BE OVERWRITTEN UPON REGENERATION
+!! If changes need to occur, modify the appropriate CIX xml file
+!! Date Generated: Wed Apr 08 13:29:36 EDT 2020
+!! If any issues are experiences with this generated file that cannot be fixed
+!! with adjustment of the CIX xml file, please contact Robert A. Lefebvre, raq@ornl.gov
+!!/
+module Yield_M
+use, intrinsic :: ISO_C_BINDING
+use Yield_I
+type Yield
+ type(C_PTR) :: instance_ptr=C_NULL_PTR
+ contains
+ procedure, pass(this) :: setNormTypeSelfShielded => Yield_setNormTypeSelfShielded
+ procedure, pass(this) :: setNormTypeDivideByN => Yield_setNormTypeDivideByN
+ procedure, pass(this) :: setNormTypeTimesSigTot => Yield_setNormTypeTimesSigTot
+ procedure, pass(this) :: setCalcDerivs => Yield_setCalcDerivs
+ procedure, pass(this) :: setNumTotalParam => Yield_setNumTotalParam
+ procedure, pass(this) :: setY0 => Yield_setY0
+ procedure, pass(this) :: setY1 => Yield_setY1
+ procedure, pass(this) :: setY2 => Yield_setY2
+ procedure, pass(this) :: setTotY => Yield_setTotY
+ procedure, pass(this) :: setDerivY0 => Yield_setDerivY0
+ procedure, pass(this) :: setDerivY1 => Yield_setDerivY1
+ procedure, pass(this) :: setDerivY2 => Yield_setDerivY2
+ procedure, pass(this) :: setDerivTotY => Yield_setDerivTotY
+ procedure, pass(this) :: getNormString => Yield_getNormString
+ procedure, pass(this) :: normTypeIsSelfShielded => Yield_normTypeIsSelfShielded
+ procedure, pass(this) :: normTypeIsDivideByN => Yield_normTypeIsDivideByN
+ procedure, pass(this) :: normTypeIsTimesSigTot => Yield_normTypeIsTimesSigTot
+ procedure, pass(this) :: getCalcDerivs => Yield_getCalcDerivs
+ procedure, pass(this) :: getNumTotalParam => Yield_getNumTotalParam
+ procedure, pass(this) :: getY0 => Yield_getY0
+ procedure, pass(this) :: getY1 => Yield_getY1
+ procedure, pass(this) :: getY2 => Yield_getY2
+ procedure, pass(this) :: getTotY => Yield_getTotY
+ procedure, pass(this) :: getDerivY0 => Yield_getDerivY0
+ procedure, pass(this) :: getDerivY1 => Yield_getDerivY1
+ procedure, pass(this) :: getDerivY2 => Yield_getDerivY2
+ procedure, pass(this) :: getDerivTotY => Yield_getDerivTotY
+ procedure, pass(this) :: initialize => Yield_initialize
+ procedure, pass(this) :: destroy => Yield_destroy
+end type Yield
+contains
+subroutine Yield_setNormTypeSelfShielded(this)
+ implicit none
+ class(Yield)::this
+ call f_Yield_setNormTypeSelfShielded(this%instance_ptr)
+end subroutine
+subroutine Yield_setNormTypeDivideByN(this)
+ implicit none
+ class(Yield)::this
+ call f_Yield_setNormTypeDivideByN(this%instance_ptr)
+end subroutine
+subroutine Yield_setNormTypeTimesSigTot(this)
+ implicit none
+ class(Yield)::this
+ call f_Yield_setNormTypeTimesSigTot(this%instance_ptr)
+end subroutine
+subroutine Yield_setCalcDerivs(this, sb)
+ implicit none
+ class(Yield)::this
+ logical(C_BOOL)::sb
+ call f_Yield_setCalcDerivs(this%instance_ptr, sb)
+end subroutine
+subroutine Yield_setNumTotalParam(this, np)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::np
+ call f_Yield_setNumTotalParam(this%instance_ptr, np)
+end subroutine
+subroutine Yield_setY0(this, y0)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE)::y0
+ call f_Yield_setY0(this%instance_ptr, y0)
+end subroutine
+subroutine Yield_setY1(this, y1)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE)::y1
+ call f_Yield_setY1(this%instance_ptr, y1)
+end subroutine
+subroutine Yield_setY2(this, y2)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE)::y2
+ call f_Yield_setY2(this%instance_ptr, y2)
+end subroutine
+subroutine Yield_setTotY(this, toty)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE)::toty
+ call f_Yield_setTotY(this%instance_ptr, toty)
+end subroutine
+subroutine Yield_setDerivY0(this, iPar, dy0)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE)::dy0
+ call f_Yield_setDerivY0(this%instance_ptr, iPar-1,dy0)
+end subroutine
+subroutine Yield_setDerivY1(this, iPar, dy1)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE)::dy1
+ call f_Yield_setDerivY1(this%instance_ptr, iPar-1,dy1)
+end subroutine
+subroutine Yield_setDerivY2(this, iPar, dy2)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE)::dy2
+ call f_Yield_setDerivY2(this%instance_ptr, iPar-1,dy2)
+end subroutine
+subroutine Yield_setDerivTotY(this, iPar, dy2)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE)::dy2
+ call f_Yield_setDerivTotY(this%instance_ptr, iPar-1,dy2)
+end subroutine
+subroutine Yield_getNormString(this, string2Populate)
+ implicit none
+ class(Yield)::this
+ character(len=*) :: string2Populate
+ integer(C_INT) :: lengthOfString2Populate
+ lengthOfString2Populate = len(string2Populate)
+ call f_Yield_getNormString(this%instance_ptr, string2Populate, lengthOfString2Populate)
+end subroutine
+function Yield_normTypeIsSelfShielded(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ logical(C_BOOL):: result2Return
+ result2Return=f_Yield_normTypeIsSelfShielded(this%instance_ptr)
+end function
+function Yield_normTypeIsDivideByN(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ logical(C_BOOL):: result2Return
+ result2Return=f_Yield_normTypeIsDivideByN(this%instance_ptr)
+end function
+function Yield_normTypeIsTimesSigTot(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ logical(C_BOOL):: result2Return
+ result2Return=f_Yield_normTypeIsTimesSigTot(this%instance_ptr)
+end function
+function Yield_getCalcDerivs(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ logical(C_BOOL):: result2Return
+ result2Return=f_Yield_getCalcDerivs(this%instance_ptr)
+end function
+function Yield_getNumTotalParam(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ integer(C_INT):: result2Return
+ result2Return=f_Yield_getNumTotalParam(this%instance_ptr)
+end function
+function Yield_getY0(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getY0(this%instance_ptr)
+end function
+function Yield_getY1(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getY1(this%instance_ptr)
+end function
+function Yield_getY2(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getY2(this%instance_ptr)
+end function
+function Yield_getTotY(this) result(result2Return)
+ implicit none
+ class(Yield)::this
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getTotY(this%instance_ptr)
+end function
+function Yield_getDerivY0(this, iPar) result(result2Return)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getDerivY0(this%instance_ptr, iPar-1)
+end function
+function Yield_getDerivY1(this, iPar) result(result2Return)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getDerivY1(this%instance_ptr, iPar-1)
+end function
+function Yield_getDerivY2(this, iPar) result(result2Return)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getDerivY2(this%instance_ptr, iPar-1)
+end function
+function Yield_getDerivTotY(this, iPar) result(result2Return)
+ implicit none
+ class(Yield)::this
+ integer(C_INT)::iPar
+ real(C_DOUBLE):: result2Return
+ result2Return=f_Yield_getDerivTotY(this%instance_ptr, iPar-1)
+end function
+subroutine Yield_initialize(this)
+ implicit none
+ class(Yield) :: this
+ this%instance_ptr = f_Yield_initialize()
+end subroutine
+subroutine Yield_destroy(this)
+ implicit none
+ class(Yield) :: this
+ call f_Yield_destroy(this%instance_ptr)
+ this%instance_ptr = C_NULL_PTR
+end subroutine
+end module Yield_M
diff --git a/sammy/src/salmon/tests/CMakeLists.txt b/sammy/src/salmon/tests/CMakeLists.txt
index 7f3bf5f57..9bbab0a95 100644
--- a/sammy/src/salmon/tests/CMakeLists.txt
+++ b/sammy/src/salmon/tests/CMakeLists.txt
@@ -2,3 +2,4 @@ include(NemesisTest)
ADD_NEMESIS_TEST(ExperimentalParametersTest.cpp TIMEOUT 60000 NP 1 )
ADD_NEMESIS_TEST(MultipleScatteringTest.cpp TIMEOUT 60000 NP 1 )
+ADD_NEMESIS_TEST(YieldTest.cpp TIMEOUT 60000 NP 1 )
diff --git a/sammy/src/salmon/tests/ExperimentalParametersTest.cpp b/sammy/src/salmon/tests/ExperimentalParametersTest.cpp
index 2d6d9c192..96f48c0b2 100644
--- a/sammy/src/salmon/tests/ExperimentalParametersTest.cpp
+++ b/sammy/src/salmon/tests/ExperimentalParametersTest.cpp
@@ -11,20 +11,23 @@ TEST(SampleDimensions,setterGetter){
sammy::SampleDimensions result;
- result.setRadius( 1.0);
- result.setHeight( 2.0);
- result.setThickness( 3.0);
- result.setBeamRadius(4.0);
-
- double goldRadius = 1.0;
- double goldHeight = 2.0;
- double goldThickness = 3.0;
- double goldBeamRadius = 4.0;
-
- ASSERT_EQ(goldRadius,result.getRadius());
- ASSERT_EQ(goldHeight,result.getHeight());
- ASSERT_EQ(goldThickness,result.getThickness());
- ASSERT_EQ(goldBeamRadius,result.getBeamRadius());
+ result.setRadius( 1.0);
+ result.setHeight( 2.0);
+ result.setThickness( 3.0);
+ result.setArealDensity(4.0);
+ result.setBeamRadius( 5.0);
+
+ double goldRadius = 1.0;
+ double goldHeight = 2.0;
+ double goldThickness = 3.0;
+ double goldArealDensity = 4.0;
+ double goldBeamRadius = 5.0;
+
+ ASSERT_EQ(goldRadius, result.getRadius());
+ ASSERT_EQ(goldHeight, result.getHeight());
+ ASSERT_EQ(goldThickness, result.getThickness());
+ ASSERT_EQ(goldArealDensity,result.getArealDensity());
+ ASSERT_EQ(goldBeamRadius, result.getBeamRadius());
}
/**********************************************************
diff --git a/sammy/src/salmon/tests/YieldTest.cpp b/sammy/src/salmon/tests/YieldTest.cpp
new file mode 100644
index 000000000..3ce0bb696
--- /dev/null
+++ b/sammy/src/salmon/tests/YieldTest.cpp
@@ -0,0 +1,58 @@
+#include "Nemesis/gtest/Gtest_Functions.hh"
+#include "Nemesis/gtest/nemesis_gtest.hh"
+#include "../Yield.h"
+
+
+
+/**************************************************************************************************
+* Test the getting and setting functions
+**************************************************************************************************/
+TEST(Yield,setterGetter){
+
+ sammy::Yield yield;
+
+ int gold_def_normType_NOTSET = true;
+ int gold_normType_tst = true; // <-- Normtype is now set to normalization::TIMESSIGTOT
+ double gold_Y0 = 0.1;
+ double gold_Y1 = 0.2;
+ double gold_Y2 = 0.3;
+ double gold_derivY0 = 0.1;
+ double gold_derivY1 = 0.2;
+ double gold_derivY2 = 0.3;
+ double gold_totY = gold_Y0 + gold_Y1 + gold_Y2;
+ double gold_derivTotY = -0.1;
+ bool gold_def_calcDerivs = false;
+ bool gold_calcDerivs = true;
+ int gold_numParam = 20;
+
+
+ yield.setY0( gold_Y0);
+ yield.setY1( gold_Y1);
+ yield.setY2( gold_Y2);
+ yield.setTotY( gold_totY);
+ yield.setDerivY0(0, gold_derivY0);
+ yield.setDerivY1(0, gold_derivY1);
+ yield.setDerivY2(0, gold_derivY2);
+ yield.setDerivTotY(0, gold_derivTotY);
+ yield.setNumTotalParam( gold_numParam);
+
+
+ ASSERT_EQ(gold_def_normType_NOTSET, yield.normTypeIsNotSet()); // <-- Default normalization type
+ ASSERT_EQ(gold_def_calcDerivs, yield.getCalcDerivs());
+
+ yield.setNormTypeTimesSigTot();
+ yield.setCalcDerivs(gold_calcDerivs);
+
+ ASSERT_EQ(gold_normType_tst, yield.normTypeIsTimesSigTot());
+ ASSERT_EQ(gold_Y0, yield.getY0());
+ ASSERT_EQ(gold_Y1, yield.getY1());
+ ASSERT_EQ(gold_Y2, yield.getY2());
+ ASSERT_EQ(gold_derivY0, yield.getDerivY0(0));
+ ASSERT_EQ(gold_derivY1, yield.getDerivY1(0));
+ ASSERT_EQ(gold_derivY2, yield.getDerivY2(0));
+ ASSERT_EQ(gold_totY, yield.getTotY());
+ ASSERT_EQ(gold_derivTotY, yield.getDerivTotY(0));
+ ASSERT_EQ(gold_numParam, yield.getNumTotalParam());
+ ASSERT_EQ(gold_calcDerivs, yield.getCalcDerivs());
+
+}
diff --git a/sammy/src/sam/msam.F b/sammy/src/sam/msam.F
index 602134bd4..cce141037 100755
--- a/sammy/src/sam/msam.F
+++ b/sammy/src/sam/msam.F
@@ -16,6 +16,7 @@ C
use EndfData_common_m
use MultScatPars_common_m
use ExpPars_common_m
+ use Observable_common_m
use ssm_m
use Samdat_0_M
use Samodf_0_m
@@ -43,6 +44,7 @@ C
call multScat%initialize()
call sampleDim%initialize()
call covData%initialize()
+ call yld%initialize()
call covData%makeNewUCovariance()
call covData%makeNewCovariance()
call covData%makeNewInverseUCovariance()
diff --git a/sammy/src/sammy/CMakeLists.txt b/sammy/src/sammy/CMakeLists.txt
index 8eae11683..bc1bfa2dd 100644
--- a/sammy/src/sammy/CMakeLists.txt
+++ b/sammy/src/sammy/CMakeLists.txt
@@ -572,7 +572,8 @@ APPEND_SET(SAMMY_SOURCES
../blk/AllocateFunctions.f90
../blk/ExpPars_common.f90
../blk/MultScatPars_common.f90
- ../blk/CapYCorrections_common.f90
+ ../blk/Convolution_common.f90
+ ../blk/Observable_common.f90
)
TRIBITS_ADD_LIBRARY(
diff --git a/sammy/src/ssm/mssm04.f90 b/sammy/src/ssm/mssm04.f90
index 12a4ea9c1..df05f554a 100644
--- a/sammy/src/ssm/mssm04.f90
+++ b/sammy/src/ssm/mssm04.f90
@@ -40,6 +40,7 @@ module ssm_4_m
Itimes = 1
Kountr = 0
!
+! JMB: The a and b here are Y's at energy a and b read from a file
IF (Mcy2.EQ.1) CALL Gety2zzz (Emy2_a, Emy2_b, Y2tab_a, Y2tab_b)
!
CALL Findpr (multScat%getNumTheta(), Knthet)
diff --git a/sammy/src/ssm/mssm19.f90 b/sammy/src/ssm/mssm19.f90
index 7449e786a..f885014f8 100644
--- a/sammy/src/ssm/mssm19.f90
+++ b/sammy/src/ssm/mssm19.f90
@@ -11,32 +11,47 @@ module ssm_19_m
use ifwrit_m
use fixedr_m
use logic_ssm_common_m
- use CapYCorrections_common_m
+ use Convolution_common_m
+ use Observable_common_m
+ use MultScatPars_common_m
+ use ExpPars_common_m
+ use Yield_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ logical(C_BOOL) calc_derivs
DIMENSION Capsig(*), Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Sigxxx(*), Dasigx(*), Dbsigx(*)
DATA Zero /0.0d0/, One /1.0d0/
call capYieldCor%initialize()
+ call yNorm%initialize()
!
! *** Default normalization is such that Y0 = (1-e)Capture /(n Total)
-! *** which corresponds to Kyield=1 ("cross section", not "yield")
+! *** which corresponds to yld%normTypeIsDivideByN()==true ("cross section", not "yield")
!
IF (Noyzer.NE.0) Y0 = Zero
+
+ ! -----------------------------------
+ ! Put things into C++ memory here
+ ! -----------------------------------
+ call yld%setY0(Y0)
+ calc_derivs = Ksolve.NE.2
+ call yld%setCalcDerivs(calc_derivs)
+ call multScat%setNumScatters(0)
+ call yld%setNumTotalParam(Nnpar)
+ call sampleDim%setArealDensity(Dthick)
!
Cap = Capsig(Nn)
- IF (Kyield.EQ.0) THEN
- Y0 = Y0*Dthick
- ELSE IF (Kyield.EQ.1) THEN
- ELSE IF (Kyield.EQ.2) THEN
- Y0 = Y0*Total*Dthick
- END IF
+
+ ! -----------------------------------
+ ! Now normalize the yield
+ ! -----------------------------------
+ call yNorm%norm(yld,sampleDim,multScat,Total)
!
- Sigxxx(1) = Y0
+ call yld%setTotY(yld%getY0())
+ Sigxxx(1) = yld%getTotY()
+ Y0 = yld%getY0()
!
! If the gamma attenuation factor C1 is not zero
IF (Gam_att_C1.NE.0.0) THEN
- dC1 = dGam_att_C1
- Td = Total*Dthick
! Sigxxx(1) is the place we store the output to ODF/LST files
! -- Here we are modifying Sigxxx(1) with the gamma attenuation
! -- correction
@@ -61,13 +76,13 @@ module ssm_19_m
IF (Kvthck.GT.0) THEN
! *** Add deriv of Y0 wrt Thickness
Kv = Kvthck - Nvadif - Ndasig
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
IF (Noyzer.EQ.0) THEN
Dbsigx(Kv) = Exp1*Cap + Dbsigx(Kv)
! ELSE
! disabling Y0 portion of calculation
END IF
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Sn = Dthick*Total
IF (Noyzer.EQ.0) THEN
IF (Sn.LT.0.3d0) THEN
@@ -81,18 +96,21 @@ module ssm_19_m
F = Zero
END IF
Dbsigx(Kv) = F + Dbsigx(Kv)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
IF (Noyzer.EQ.0) Dbsigx(Kv)= Exp1*Cap*Total + Dbsigx(Kv)
END IF
END IF
!
IF (Sensin.NE.Zero .AND. Ksensc.GT.0) THEN
- IF (Kyield.EQ.0) Asensn = Asensn * Dthick
- IF (Kyield.EQ.2) Asensn = Asensn * Dthick * Total
+ IF (yld%normTypeIsSelfShielded()) Asensn = Asensn * Dthick
+ IF (yld%normTypeIsTimesSigTot()) Asensn = Asensn * Dthick * Total
! *** Add deriv of Y0 wrt neutron sensitivity multiplier
Kv = Ksensc - Nvadif - Ndasig
Dbsigx(Kv) = Asensn
END IF
+
+ call capYieldCor%destroy()
+
RETURN
END
!
@@ -104,7 +122,8 @@ module ssm_19_m
use fixedi_m
use ifwrit_m
use logic_ssm_common_m
- use CapYCorrections_common_m
+ use Convolution_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Ds(*)
DATA One /1.0d0/
@@ -120,18 +139,18 @@ module ssm_19_m
IF (Noyzer.EQ.0) THEN
! *** First, the derivative of Y0 wrt u(Iipar)
! *** (Nn=Energy) JMB
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total * Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
IF (Iipar.EQ.1) Dddyy0 = Dddyy0 + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + (One-Exp1)*Dcapsi(Iipar,Nn) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap*Dthick
END IF
diff --git a/sammy/src/ssm/mssm20.f90 b/sammy/src/ssm/mssm20.f90
index 8ce787e29..9c3c3ce70 100644
--- a/sammy/src/ssm/mssm20.f90
+++ b/sammy/src/ssm/mssm20.f90
@@ -149,54 +149,72 @@ module ssm_20_m
use ifwrit_m
use fixedr_m
use logic_ssm_common_m
- use CapYCorrections_common_m
+ use Convolution_common_m
+ use Observable_common_m
+ use MultScatPars_common_m
+ use ExpPars_common_m
+ use Yield_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ logical(C_BOOL) calc_derivs
DIMENSION Capsig(*), Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*), &
Sigxxx(*), Dasigx(*), Dbsigx(*)
DATA Zero /0.0d0/, One /1.0d0/
call capYieldCor%initialize()
+ call yNorm%initialize()
!
! *** Default normalization is such that Y0 = (1-e)Capture /(n Total)
-! *** which corresponds to Kyield=1 ("cross section", not "yield")
+! *** which corresponds to yld%normTypeIsDivideByN()==true ("cross section", not "yield")
!
IF (Noyzer.NE.0) Y0 = Zero
+
+ ! -----------------------------------
+ ! Put things into C++ memory here
+ ! -----------------------------------
+ call yld%setY0(Y0)
+ call yld%setY1(Y1)
+ if( Ksolve.ne.2 ) then
+ do j=1,Nnpar
+ call yld%setDerivY1(J,Dy1(J))
+ end do
+ end if
+ calc_derivs = Ksolve.NE.2
+ call yld%setCalcDerivs(calc_derivs)
+ call multScat%setNumScatters(1)
+ call yld%setNumTotalParam(Nnpar)
+ call sampleDim%setArealDensity(Dthick)
+
+ ! -----------------------------------
+ ! Now normalize the yield
+ ! -----------------------------------
+ call yNorm%norm(yld,sampleDim,multScat,Total)
!
Cap = Capsig(Nn)
- IF (Kyield.EQ.0) THEN
- Y0 = Y0*Dthick
- ELSE IF (Kyield.EQ.1) THEN
- Y1 = Y1/Dthick
- IF (Ksolve.NE.2) THEN
- DO J=1,Nnpar
- Dy1(J) = Dy1(J)/Dthick
- END DO
- END IF
- ELSE IF (Kyield.EQ.2) THEN
- Y0 = Y0*Total*Dthick
- Y1 = Y1*Total
- IF (Ksolve.NE.2) THEN
- DO J=1,Nnpar
- Dy1(J) = Dy1(J)*Total
- END DO
- END IF
- END IF
-!
-!
- Sigxxx(1) = Y0
- Sigxxx(1) = Sigxxx(1) + Y1
+
+ ! -----------------------------------
+ ! Return the values to the old places
+ ! -----------------------------------
+ call yld%setTotY(yld%getY0() + yld%getY1())
+ Sigxxx(1) = yld%getTotY()
+ Y0 = yld%getY0()
+ Y1 = yld%getY1()
+ if( yld%getCalcDerivs() ) then
+ do j=1,yld%getNumTotalParam()
+ Dy1(j) = yld%getDerivY1(j)
+ end do
+ end if
!
! If the gamma attenuation factor C1 is not zero
IF (Gam_att_C1.NE.0.0) THEN
- dC1 = dGam_att_C1
- Td = Total*Dthick
! Sigxxx(1) is the place we store the output to ODF/LST files
! -- Here we are modifying Sigxxx(1) with the gamma attenuation
! -- correction
call capYieldCor%gammaAttenCorr( Dthick,Total,Gam_att_C1,Sigxxx(1) )
END IF
!
-!
+ ! -----------------------------------
+ ! Write yields to file
+ ! -----------------------------------
IF (Kssmpr.EQ.1) THEN
WRITE (60,12000) Em, Y0, Y1, Sigxxx(1)
12000 FORMAT (1X, f15.8, 1P5E15.6)
@@ -213,13 +231,13 @@ module ssm_20_m
IF (Kvthck.GT.0) THEN
! *** Add deriv of Y0 wrt Thickness; modify deriv of Y1 & Y2 wrt Thick
Kv = Kvthck - Nvadif - Ndasig
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
IF (Noyzer.EQ.0) THEN
Dbsigx(Kv) = Exp1*Cap + Dbsigx(Kv)
! ELSE
! disabling Y0 portion of calculation
END IF
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Sn = Dthick*Total
IF (Noyzer.EQ.0) THEN
IF (Sn.LT.0.3d0) THEN
@@ -233,18 +251,21 @@ module ssm_20_m
F = Zero
END IF
Dbsigx(Kv) = F + Dbsigx(Kv) - Y1/Dthick
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
IF (Noyzer.EQ.0) Dbsigx(Kv)= Exp1*Cap*Total + Dbsigx(Kv)
END IF
END IF
!
IF (Sensin.NE.Zero .AND. Ksensc.GT.0) THEN
- IF (Kyield.EQ.0) Asensn = Asensn * Dthick
- IF (Kyield.EQ.2) Asensn = Asensn * Dthick * Total
+ IF (yld%normTypeIsSelfShielded()) Asensn = Asensn * Dthick
+ IF (yld%normTypeIsTimesSigTot()) Asensn = Asensn * Dthick * Total
! *** Add deriv of Y0 wrt neutron sensitivity multiplier
Kv = Ksensc - Nvadif - Ndasig
Dbsigx(Kv) = Asensn
END IF
+
+ call capYieldCor%destroy()
+
RETURN
END
!
@@ -256,7 +277,8 @@ module ssm_20_m
use fixedi_m
use ifwrit_m
use logic_ssm_common_m
- use CapYCorrections_common_m
+ use Convolution_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*), Ds(*)
DATA One /1.0d0/
@@ -272,29 +294,29 @@ module ssm_20_m
IF (Noyzer.EQ.0) THEN
! *** First, the derivative of Y0 wrt u(Iipar)
! *** (Nn=Energy) JMB
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total * Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
IF (Iipar.EQ.1) Dddyy0 = Dddyy0 + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + (One-Exp1)*Dcapsi(Iipar,Nn) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap*Dthick
END IF
END IF
!
! *** Next, the derivative of Y1 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar) + &
Y1*Dtotsi(Iipar,Nn)/Total
END IF
diff --git a/sammy/src/ssm/mssm21.f90 b/sammy/src/ssm/mssm21.f90
index 49047ddeb..0bcd285b6 100644
--- a/sammy/src/ssm/mssm21.f90
+++ b/sammy/src/ssm/mssm21.f90
@@ -49,9 +49,16 @@ module ssm_21_m
use fixedi_m
use ifwrit_m
use fixedr_m
- use CapYCorrections_common_m
+ use Convolution_common_m
use logic_ssm_common_m
+ use Observable_common_m
+ use MultScatPars_common_m
+ use ExpPars_common_m
+ use EndfData_common_m
+ use Yield_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
+ logical(C_BOOL) calc_derivs
+ double precision, dimension(:), allocatable :: Dabsigx
DIMENSION Capsig(*), Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*),&
Dy2(*), Sigxxx(*), Dasigx(*), Dbsigx(*)
DATA Zero /0.0d0/, One /1.0d0/
@@ -59,54 +66,72 @@ module ssm_21_m
!type(CapYieldCorrections)::capYieldCor
call capYieldCor%initialize()
+ call yNorm%initialize()
+
+ ! ------- Nnpar, Ndasig, and Ndbsig are set to zero when Ksolve = 2 -----------
+
!
! *** Default normalization is such that Y0 = (1-e)Capture /(n Total)
-! *** which corresponds to Kyield=1 ("cross section", not "yield")
+! *** which corresponds to yld%normTypeIsDivideByN()==true ("cross section", not "yield")
!
IF (Noyzer.NE.0) Y0 = Zero
+
+ ! -----------------------------------
+ ! Put things into C++ memory here
+ ! -----------------------------------
+ call yld%setY0(Y0)
+ call yld%setY1(Y1)
+ call yld%setY2(Y2)
+ if( Ksolve.ne.2 ) then
+ do j=1,Nnpar
+ call yld%setDerivY1(J,Dy1(J))
+ call yld%setDerivY2(J,Dy2(J))
+ end do
+ end if
+ calc_derivs = Ksolve.NE.2
+ call yld%setCalcDerivs(calc_derivs)
+ call multScat%setNumScatters(2)
+ call yld%setNumTotalParam(Nnpar)
+ call sampleDim%setArealDensity(Dthick)
+
+ ! -----------------------------------
+ ! Now normalize the yield
+ ! -----------------------------------
+ call yNorm%norm(yld,sampleDim,multScat,Total)
!
Cap = Capsig(Nn)
- IF (Kyield.EQ.0) THEN
- Y0 = Y0*Dthick
- ELSE IF (Kyield.EQ.1) THEN
- Y1 = Y1/Dthick
- Y2 = Y2/Dthick
- IF (Ksolve.NE.2) THEN
- DO J=1,Nnpar
- Dy1(J) = Dy1(J)/Dthick
- END DO
- DO J=1,Nnpar
- Dy2(J) = Dy2(J)/Dthick
- END DO
- END IF
- ELSE IF (Kyield.EQ.2) THEN
- Y0 = Y0*Total*Dthick
- Y1 = Y1*Total
- Y2 = Y2*Total
- IF (Ksolve.NE.2) THEN
- DO J=1,Nnpar
- Dy1(J) = Dy1(J)*Total
- END DO
- DO J=1,Nnpar
- Dy2(J) = Dy2(J)*Total
- END DO
- END IF
- END IF
- IF (Mcy2.NE.0) Y2 = Y2tab
-!
- Sigxxx(1) = Y0 + Y1 + Y2
+
+ IF (Mcy2.NE.0) then
+ Y2 = Y2tab
+ call yld%setY2(Y2tab)
+ end if
+!
+ ! -----------------------------------
+ ! Return the values to the old places
+ ! -----------------------------------
+ call yld%setTotY(yld%getY0() + yld%getY1() + yld%getY2())
+ Sigxxx(1) = yld%getTotY()
+ Y0 = yld%getY0()
+ Y1 = yld%getY1()
+ Y2 = yld%getY2()
+ if( yld%getCalcDerivs() ) then
+ do j=1,yld%getNumTotalParam()
+ Dy1(j) = yld%getDerivY1(j)
+ Dy2(j) = yld%getDerivY2(j)
+ end do
+ end if
!
! If the gamma attenuation factor C1 is not zero
IF (Gam_att_C1.NE.0.0) THEN
- dC1 = dGam_att_C1
- Td = Total*Dthick
! Sigxxx(1) is the place we store the output to ODF/LST files
! -- Here we are modifying Sigxxx(1) with the gamma attenuation
! -- correction
call capYieldCor%gammaAttenCorr( Dthick,Total,Gam_att_C1,Sigxxx(1) )
END IF
!
-!
+ ! -----------------------------------
+ ! Write yields to file
+ ! -----------------------------------
IF (Kssmpr.EQ.1) THEN
WRITE (60,12000) Em, Y0, Y1, Y2, Sigxxx(1)
12000 FORMAT (1X, f15.8, 1P5E15.6)
@@ -115,7 +140,7 @@ module ssm_21_m
IF (Ksolve.EQ.2) RETURN
! ######################################
! It seems that Ndasig is a subset of varied parameters and that
-! Ndbsig is the total number of varied parameters minus Ndasig
+! Ndbsig is another subset of varied parameters minus Ndasig
! ######################################
IF (Ndasig.GT.0) CALL Ynrm_2_Derivative (Cap, Dtotsi, Dcapsi, &
Dasigx, Dy1, Dy2, Y0, Y1, Y2, Dthick, Total, Exp1, Dddyy0, &
@@ -127,13 +152,13 @@ module ssm_21_m
IF (Kvthck.GT.0) THEN
! *** Add deriv of Y0 wrt Thickness; modify deriv of Y1 & Y2 wrt Thick
Kv = Kvthck - Nvadif - Ndasig
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
IF (Noyzer.EQ.0) THEN
Dbsigx(Kv) = Exp1*Cap + Dbsigx(Kv)
! ELSE
! disabling Y0 portion of calculation
END IF
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Sn = Dthick*Total
IF (Noyzer.EQ.0) THEN
IF (Sn.LT.0.3d0) THEN
@@ -148,18 +173,21 @@ module ssm_21_m
END IF
Dbsigx(Kv) = F + Dbsigx(Kv) - Y1/Dthick
IF (Double_Plus) Dbsigx(Kv) = Dbsigx(Kv)-Y2/Dthick
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
IF (Noyzer.EQ.0) Dbsigx(Kv)= Exp1*Cap*Total + Dbsigx(Kv)
END IF
END IF
!
IF (Sensin.NE.Zero .AND. Ksensc.GT.0) THEN
- IF (Kyield.EQ.0) Asensn = Asensn * Dthick
- IF (Kyield.EQ.2) Asensn = Asensn * Dthick * Total
+ IF (yld%normTypeIsSelfShielded()) Asensn = Asensn * Dthick
+ IF (yld%normTypeIsTimesSigTot()) Asensn = Asensn * Dthick * Total
! *** Add deriv of Y0 wrt neutron sensitivity multiplier
Kv = Ksensc - Nvadif - Ndasig
Dbsigx(Kv) = Asensn
END IF
+
+ call capYieldCor%destroy()
+
RETURN
END
!
@@ -171,10 +199,12 @@ module ssm_21_m
use fixedi_m
use ifwrit_m
use logic_ssm_common_m
- use CapYCorrections_common_m
+ use Convolution_common_m
+ use Observable_common_m
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION Dtotsi(Nnparx,*), Dcapsi(Nnparx,*), Dy1(*), Dy2(*),Ds(*)
DATA One /1.0d0/
+
!
! Dthick = N [at/b]
! Cap = sigma_gamma
@@ -187,39 +217,39 @@ module ssm_21_m
IF (Noyzer.EQ.0) THEN
! *** First, the derivative of Y0 wrt u(Iipar)
! *** (Nn=Energy) JMB
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total * Dthick
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
IF (Iipar.EQ.1) Dddyy0 = Dddyy0 + &
Y0* (Dcapsi(Iipar,Nn)/Cap-Dtotsi(Iipar,Nn)/Total) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap/Total
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + (One-Exp1)*Dcapsi(Iipar,Nn) &
+ Dtotsi(Iipar,Nn)*Exp1*Cap*Dthick
END IF
END IF
!
! *** Next, the derivative of Y1 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + Dy1(Iipar) + &
Y1*Dtotsi(Iipar,Nn)/Total
END IF
!
! *** Finally, the derivative of Y2 wrt u(Iipar)
- IF (Kyield.EQ.0) THEN
+ IF (yld%normTypeIsSelfShielded()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar)
- ELSE IF (Kyield.EQ.1) THEN
+ ELSE IF (yld%normTypeIsDivideByN()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar)
- ELSE IF (Kyield.EQ.2) THEN
+ ELSE IF (yld%normTypeIsTimesSigTot()) THEN
Ds(Ipar) = Ds(Ipar) + Dy2(Iipar) &
+ Y2*Dtotsi(Iipar,Nn)/Total
END IF
--
GitLab