reworked interpolation + added unit test; icru77 xs are now validated; icru77...

reworked interpolation + added unit test; icru77 xs are now validated; icru77 xs added to the csd solvers needs to be checked

code/include/cubic2dspline.h

-#ifndef CUBIC2DSPLINE_H
-#define CUBIC2DSPLINE_H
-
-#include "common/typedef.h"
-
-class Cubic2DSpline
-{
-  private:
-    Vector _x;
-    Vector _y;
-    Matrix _data;
-
-    /*
-     * calculates a cubic 1D interpolation at point x
-     * @param[in]: double param[4] - contains the 4 points needed for the purpose of cubic interpolation
-     * @param[in]: double x - the desired interpolation point
-     * @param[out]: double - interpolation result
-     */
-    inline double interpolate1D( double param[4], double x );
-
-    /*
-     * finds index of the closest distance to provided value in a given vector
-     * @param[in]: double value - value to search for
-     * @param[in]: Vector v - sorted vector holding the data
-     * @param[out]: unsigned - resulting index of closest value in v
-     */
-    unsigned indexOfClosestValue( double value, const Vector& v );
-
-    Cubic2DSpline() = delete;
-
-  public:
-    /*
-     * calculates a cubic 2D interpolation at point (paramX, paramY)
-     * @param[in]: double x - x coordinate of the desired interpolation point
-     * @param[in]: double y - y coordinate of the desired interpolation point
-     * @param[out]: double - interpolation result
-     */
-    double operator()( double x, double y );
-
-    Cubic2DSpline( const Vector& x, const Vector& y, const Matrix& data );
-    ~Cubic2DSpline();
-};
-
-#endif

code/include/icru.h

@@ -6,6 +6,7 @@
 #include <sstream>
 #include <vector>
 
+#include "common/globalconstants.h"
 #include "common/typedef.h"
 #include "interpolation.h"
 #include "toolboxes/errormessages.h"
@@ -140,7 +141,7 @@ class ICRU
     }
     ~ICRU(){};
 
-    void GetAngularScatteringXS( Matrix& elastic, Matrix& inelastic, Matrix& total );
+    void GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS );
 };
 
 #endif    // ICRU_H

code/include/interpolation.h

@@ -6,29 +6,29 @@
 class Interpolation
 {
   public:
-    enum BOUNDARY { first_deriv, second_deriv };
     enum TYPE { linear, loglinear, cubic };
 
   private:
+    unsigned _dim;
     Vector _x, _y;
-    Vector _a, _b, _c;
-    double _b0, _c0;
-    BOUNDARY _left, _right;
+    Matrix _data;
     TYPE _type;
-    double _left_value, _right_value;
-    bool _force_linear_extrapolation;
 
     Interpolation() = delete;
 
+    void Setup();
+    unsigned IndexOfClosestValue( double value, const Vector& v ) const;
+    inline double EvalCubic1DSpline( double param[4], double x ) const;
+
   public:
     Interpolation( const std::vector<double>& x, const std::vector<double>& y, TYPE type = linear );
     Interpolation( const Vector& x, const Vector& y, TYPE type = linear );
+    Interpolation( const Vector& x, const Vector& y, const Matrix& data, TYPE type = cubic );
 
-    void set_boundary( BOUNDARY left, double left_value, BOUNDARY right, double right_value, bool force_linear_extrapolation = false );
-    void set_points( const std::vector<double>& x, const std::vector<double>& y, TYPE type );
-    void set_points( const Vector& x, const Vector& y, TYPE type );
     double operator()( double x ) const;
     Vector operator()( Vector v ) const;
+    std::vector<double> operator()( std::vector<double> v ) const;
+    double operator()( double x, double y ) const;
 };
 
 #endif

code/include/physics.h

@@ -7,7 +7,6 @@
 
 #include "common/globalconstants.h"
 #include "common/typedef.h"
-#include "cubic2dspline.h"
 #include "interpolation.h"
 #include "toolboxes/errormessages.h"
 

code/include/solvers/csdsnsolver.h

 #ifndef CSDSNSOLVER_H
 #define CSDSNSOLVER_H
 
+#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;

code/include/solvers/csdsnsolvernotrafo.h

 #ifndef CSDSNSOLVERNOTRAFO_H
 #define CSDSNSOLVERNOTRAFO_H
 
+#include "icru.h"
 #include "solvers/snsolver.h"
 
 class Physics;
@@ -35,4 +36,4 @@ class CSDSNSolverNoTrafo : public SNSolver
     virtual void Save( int currEnergy ) const;
 };
 
-#endif // CSDSNSOLVERNOTRAFO_H
+#endif    // CSDSNSOLVERNOTRAFO_H

code/input/example_csd.cfg

@@ -8,11 +8,11 @@ CONTINUOUS_SLOWING_DOWN = YES
 HYDROGEN_FILE = ENDL_H.txt
 OXYGEN_FILE = ENDL_O.txt
 KERNEL = ISOTROPIC_1D
-CFL_NUMBER = 0.1
+CFL_NUMBER = 0.025
 TIME_FINAL = 1.0
 CLEAN_FLUX_MATRICES = NO
 
 BC_DIRICHLET = ( void )
 BC_NEUMANN = ( adiabatic )
 QUAD_TYPE = GAUSS_LEGENDRE_1D
-QUAD_ORDER = 21
+QUAD_ORDER = 101

code/src/cubic2dspline.cpp

-#include "cubic2dspline.h"
-
-Cubic2DSpline::Cubic2DSpline( const Vector& x, const Vector& y, const Matrix& data ) : _x( x ), _y( y ), _data( data ) {}
-Cubic2DSpline::~Cubic2DSpline() {}
-
-inline double Cubic2DSpline::interpolate1D( double param[4], double x ) {
-    return ( param[1] + 0.5 * x *
-                            ( param[2] - param[0] +
-                              x * ( 2.0 * param[0] - 5.0 * param[1] + 4.0 * param[2] - param[3] +
-                                    x * ( 3.0 * ( param[1] - param[2] ) + param[3] - param[0] ) ) ) );
-}
-
-double Cubic2DSpline::operator()( double x, double y ) {
-    unsigned xId = indexOfClosestValue( x, _x );
-    unsigned yId = indexOfClosestValue( y, _y );
-
-    // store all 16 interpolation points needed
-    double points[4][4];
-    for( int i = -1; i < 3; ++i ) {
-        unsigned idx_y;
-        idx_y = yId + i < 0 ? 0 : yId + i;
-        idx_y = yId + i > _data.rows() - 1 ? _data.rows() - 1 : yId + i;
-        for( int j = -1; j < 3; ++j ) {
-            unsigned idx_x;
-            idx_x = xId + j < 0 ? 0 : xId + j;
-            idx_x = xId + j > _data.columns() - 1 ? _data.columns() - 1 : xId + j;
-
-            points[i + 1][j + 1] = _data( idx_x, idx_y );
-        }
-    }
-
-    // rescale data to [0,1]
-    double t = ( x - _x[xId] ) / ( _x[xId + 1] - _x[xId] );
-    double u = ( y - _y[yId] ) / ( _y[yId + 1] - _y[yId] );
-
-    // first interpolate in x-direction and store the results on which the final interpolation in y will be done
-    double interpolationBuffer[4];
-    interpolationBuffer[0] = interpolate1D( points[0], t );
-    interpolationBuffer[1] = interpolate1D( points[1], t );
-    interpolationBuffer[2] = interpolate1D( points[2], t );
-    interpolationBuffer[3] = interpolate1D( points[3], t );
-    return interpolate1D( interpolationBuffer, u );
-}
-
-unsigned Cubic2DSpline::indexOfClosestValue( double value, const Vector& v ) {
-    auto i = std::min_element( begin( v ), end( v ), [=]( double x, double y ) { return abs( x - value ) < abs( y - value ); } );
-    return std::distance( begin( v ), i );
-}

code/src/icru.cpp

@@ -447,18 +447,18 @@ double ICRU::DODCSC( double RMU, double EL ) {
     return res;
 }
 
-void ICRU::GetAngularScatteringXS( Matrix& elastic, Matrix& inelastic, Matrix& total ) {
-    elastic.resize( _QMU.size(), _E.size() );
-    inelastic.resize( _QMU.size(), _E.size() );
-    total.resize( _QMU.size(), _E.size() );
+void ICRU::GetAngularScatteringXS( Matrix& angularXS, Vector& integratedXS ) {
+    angularXS.resize( _QMU.size(), _E.size() );
+    integratedXS.resize( _E.size() );
     for( unsigned i = 0; i < _E.size(); ++i ) {
         std::vector<double> dxse, dxsi, dxs;
         angdcs( 3u, _E[i], dxse, dxsi, dxs );
-        Interpolation interpDXSE( _XMU, dxse );
-        Interpolation interpDXSI( _XMU, dxsi );
         Interpolation interpDXS( _XMU, dxs );
-        blaze::column( elastic, i )   = interpDXSE( _QMU );
-        blaze::column( inelastic, i ) = interpDXSI( _QMU );
-        blaze::column( total, i )     = interpDXS( _QMU );
+        blaze::column( angularXS, i ) = interpDXS( _QMU );
+        for( unsigned j = 0; j < _XMU.size() - 1; ++j ) {
+            integratedXS[i] += 0.5 * ( _XMU[j + 1] - _XMU[j] ) * ( dxs[j + 1] + dxs[j] );
+        }
     }
+    angularXS *= H2OMolecularDensity;
+    integratedXS *= H2OMolecularDensity;
 }

code/src/interpolation.cpp

@@ -2,136 +2,108 @@
 #include "toolboxes/errormessages.h"
 #include <blaze/math/lapack/posv.h>
 
-Interpolation::Interpolation( const std::vector<double>& x, const std::vector<double>& y, TYPE type )
-    : _left( first_deriv ), _right( first_deriv ), _type( type ), _left_value( 0.0 ), _right_value( 0.0 ), _force_linear_extrapolation( false ) {
-    this->set_points( x, y, type );
+Interpolation::Interpolation( const std::vector<double>& x, const std::vector<double>& y, TYPE type ) : _dim( 1u ), _type( type ) {
+    _x = Vector( x.size(), x.data() );
+    _y = Vector( y.size(), y.data() );
+    Setup();
 }
 
-Interpolation::Interpolation( const Vector& x, const Vector& y, TYPE type )
-    : _left( first_deriv ), _right( first_deriv ), _type( type ), _left_value( 0.0 ), _right_value( 0.0 ), _force_linear_extrapolation( false ) {
-    this->set_points( x, y, type );
-}
+Interpolation::Interpolation( const Vector& x, const Vector& y, TYPE type ) : _dim( 1u ), _x( x ), _y( y ), _type( type ) { Setup(); }
 
-void Interpolation::set_boundary( BOUNDARY left, double left_value, BOUNDARY right, double right_value, bool force_linear_extrapolation ) {
-    _left                       = left;
-    _right                      = right;
-    _left_value                 = left_value;
-    _right_value                = right_value;
-    _force_linear_extrapolation = force_linear_extrapolation;
+Interpolation::Interpolation( const Vector& x, const Vector& y, const Matrix& data, TYPE type )
+    : _dim( 2u ), _x( x ), _y( y ), _data( data ), _type( type ) {
+    Setup();
 }
 
-void Interpolation::set_points( const std::vector<double>& x, const std::vector<double>& y, TYPE type ) {
-    Vector xn( x.size(), 0.0 );
-    Vector yn( y.size(), 0.0 );
-    for( unsigned i = 0; i < x.size(); ++i ) {
-        xn[i] = x[i];
-        yn[i] = y[i];
+void Interpolation::Setup() {
+    if( _dim == 1u ) {
+        if( _x.size() != _y.size() ) ErrorMessages::Error( "Vectors are of unequal length!", CURRENT_FUNCTION );
+        if( _type == loglinear ) _y = blaze::log( _y );
+
+        for( unsigned i = 0; i < _x.size() - 1u; i++ ) {
+            if( !( _x[i] < _x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
+        }
+    }
+    else if( _dim == 2u ) {
+        if( _x.size() != _data.rows() ) ErrorMessages::Error( "x and data are of unequal length!", CURRENT_FUNCTION );
+        if( _y.size() != _data.columns() ) ErrorMessages::Error( "y and data are of unequal length!", CURRENT_FUNCTION );
+        for( unsigned i = 0; i < _x.size() - 1u; i++ ) {
+            if( !( _x[i] < _x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
+        }
+        for( unsigned i = 0; i < _y.size() - 1u; i++ ) {
+            if( !( _y[i] < _y[i + 1] ) ) ErrorMessages::Error( "y is not sorted ascendingly!", CURRENT_FUNCTION );
+        }
     }
-    this->set_points( xn, yn, type );
 }
 
-void Interpolation::set_points( const Vector& x, const Vector& y, TYPE type ) {
-    if( x.size() != y.size() ) ErrorMessages::Error( "Vectors are of unequal length!", CURRENT_FUNCTION );
-    _x = x;
-    if( type == loglinear )
-        _y = log( y );
-    else
-        _y = y;
-    int n = x.size();
-    for( int i = 0; i < n - 1; i++ ) {
-        if( !( _x[i] < _x[i + 1] ) ) ErrorMessages::Error( "x is not sorted ascendingly!", CURRENT_FUNCTION );
+double Interpolation::operator()( double x ) const {
+    if( _dim != 1u ) ErrorMessages::Error( "Invalid data dimension for operator(x)!", CURRENT_FUNCTION );
+
+    if( x < _x[0] || x > _x[_x.size() - 1u] ) {
+        ErrorMessages::Error( "Extrapolation is not supported!", CURRENT_FUNCTION );
     }
 
-    if( type == cubic ) {
-        Matrix A( n, n, 0.0 );    // TODO: should be a sparse matrix!
-        Vector rhs( n, 0.0 );
-        for( int i = 1; i < n - 1; i++ ) {
-            A( i, i - 1 ) = 1.0 / 3.0 * ( x[i] - x[i - 1] );
-            A( i, i )     = 2.0 / 3.0 * ( x[i + 1] - x[i - 1] );
-            A( i, i + 1 ) = 1.0 / 3.0 * ( x[i + 1] - x[i] );
-            rhs[i]        = ( y[i + 1] - y[i] ) / ( x[i + 1] - x[i] ) - ( y[i] - y[i - 1] ) / ( x[i] - x[i - 1] );
-        }
-        if( _left == Interpolation::second_deriv ) {
-            A( 0, 0 ) = 2.0;
-            A( 0, 1 ) = 0.0;
-            rhs[0]    = _left_value;
-        }
-        else if( _left == Interpolation::first_deriv ) {
-            A( 0, 0 ) = 2.0 * ( x[1] - x[0] );
-            A( 0, 1 ) = 1.0 * ( x[1] - x[0] );
-            rhs[0]    = 3.0 * ( ( y[1] - y[0] ) / ( x[1] - x[0] ) - _left_value );
-        }
-        else {
-            ErrorMessages::Error( "Invalid bd_type", CURRENT_FUNCTION );
-        }
-        if( _right == Interpolation::second_deriv ) {
-            A( n - 1, n - 1 ) = 2.0;
-            A( n - 1, n - 2 ) = 0.0;
-            rhs[n - 1]        = _right_value;
-        }
-        else if( _right == Interpolation::first_deriv ) {
-            A( n - 1, n - 1 ) = 2.0 * ( x[n - 1] - x[n - 2] );
-            A( n - 1, n - 2 ) = 1.0 * ( x[n - 1] - x[n - 2] );
-            rhs[n - 1]        = 3.0 * ( _right_value - ( y[n - 1] - y[n - 2] ) / ( x[n - 1] - x[n - 2] ) );
-        }
-        else {
-            ErrorMessages::Error( "Invalid boundary type!", CURRENT_FUNCTION );
-        }
+    Vector::ConstIterator it = std::lower_bound( _x.begin(), _x.end(), x );
 
-        _b = rhs;
-        blaze::posv( A, _b, 'U' );
+    unsigned idx = static_cast<unsigned>( std::max( int( it - _x.begin() ) - 1, 0 ) );
 
-        _a.resize( n );
-        _c.resize( n );
-        for( int i = 0; i < n - 1; i++ ) {
-            _a[i] = 1.0 / 3.0 * ( _b[i + 1] - _b[i] ) / ( x[i + 1] - x[i] );
-            _c[i] = ( y[i + 1] - y[i] ) / ( x[i + 1] - x[i] ) - 1.0 / 3.0 * ( 2.0 * _b[i] + _b[i + 1] ) * ( x[i + 1] - x[i] );
-        }
+    if( _type == linear || _type == loglinear ) {
+        double res = _y[idx] + ( _y[idx + 1] - _y[idx] ) / ( _x[idx + 1] - _x[idx] ) * ( x - _x[idx] );
+
+        if( _type == loglinear )
+            return std::exp( res );
+        else
+            return res;
     }
-    else if( type == linear || type == loglinear ) {
-        _a.resize( n );
-        _b.resize( n );
-        _c.resize( n );
-        for( int i = 0; i < n - 1; i++ ) {
-            _a[i] = 0.0;
-            _b[i] = 0.0;
-            _c[i] = ( _y[i + 1] - _y[i] ) / ( _x[i + 1] - _x[i] );
-        }
+    else if( _type == cubic ) {
+        double param[4] = { _y[idx - 1], _y[idx], _y[idx + 1], _y[idx + 2] };
+        double t        = ( x - _x[idx] ) / ( _x[idx + 1] - _x[idx] );
+        return EvalCubic1DSpline( param, t );
     }
     else {
-        ErrorMessages::Error( "Invalid interpolation type!", CURRENT_FUNCTION );
+        ErrorMessages::Error( "Invalid type!", CURRENT_FUNCTION );
+        return -1.0;
     }
+}
 
-    _b0 = ( _force_linear_extrapolation == false ) ? _b[0] : 0.0;
-    _c0 = _c[0];
+double Interpolation::operator()( double x, double y ) const {
+    if( _dim != 2u ) ErrorMessages::Error( "Invalid data dimension for operator(x,y)!", CURRENT_FUNCTION );
+    if( _type == cubic ) {
 
-    double h  = x[n - 1] - x[n - 2];
-    _a[n - 1] = 0.0;
-    _c[n - 1] = 3.0 * _a[n - 2] * h * h + 2.0 * _b[n - 2] * h + _c[n - 2];
-    if( _force_linear_extrapolation == true ) _b[n - 1] = 0.0;
-}
+        unsigned xId = IndexOfClosestValue( x, _x );
+        unsigned yId = IndexOfClosestValue( y, _y );
 
-double Interpolation::operator()( double x ) const {
-    size_t n = _x.size();
-    Vector::ConstIterator it;
-    it      = std::lower_bound( _x.begin(), _x.end(), x );
-    int idx = std::max( int( it - _x.begin() ) - 1, 0 );
-
-    double h = x - _x[idx];
-    double interpol;
-    if( x < _x[0] ) {
-        interpol = ( _b0 * h + _c0 ) * h + _y[0];
-    }
-    else if( x > _x[n - 1] ) {
-        interpol = ( _b[n - 1] * h + _c[n - 1] ) * h + _y[n - 1];
+        // store all 16 interpolation points needed
+        double points[4][4];
+        for( int i = -1; i < 3; ++i ) {
+            unsigned idx_y;
+            idx_y = yId + i < 0 ? 0 : yId + i;
+            idx_y = yId + i > _data.rows() - 1 ? _data.rows() - 1 : yId + i;
+            for( int j = -1; j < 3; ++j ) {
+                unsigned idx_x;
+                idx_x = xId + j < 0 ? 0 : xId + j;
+                idx_x = xId + j > _data.columns() - 1 ? _data.columns() - 1 : xId + j;
+
+                points[i + 1][j + 1] = _data( idx_x, idx_y );
+            }
+        }
+
+        // rescale data to [0,1]
+        double t = ( x - _x[xId] ) / ( _x[xId + 1] - _x[xId] );
+        double u = ( y - _y[yId] ) / ( _y[yId + 1] - _y[yId] );
+
+        // first interpolate in x-direction and store the results on which the final interpolation in y will be done
+        double buffer[4];
+        buffer[0] = EvalCubic1DSpline( points[0], t );
+        buffer[1] = EvalCubic1DSpline( points[1], t );
+        buffer[2] = EvalCubic1DSpline( points[2], t );
+        buffer[3] = EvalCubic1DSpline( points[3], t );
+        return EvalCubic1DSpline( buffer, u );
     }
     else {
-        interpol = ( ( _a[idx] * h + _b[idx] ) * h + _c[idx] ) * h + _y[idx];
+        ErrorMessages::Error( "Unsupported interpolation type!", CURRENT_FUNCTION );
+        return -1.0;
     }
-    if( _type == loglinear )
-        return std::exp( interpol );
-    else
-        return interpol;
 }
 
 Vector Interpolation::operator()( Vector v ) const {
@@ -141,3 +113,23 @@ Vector Interpolation::operator()( Vector v ) const {
     }
     return res;
 }
+
+std::vector<double> Interpolation::operator()( std::vector<double> v ) const {
+    std::vector<double> res( v.size() );
+    for( unsigned i = 0; i < v.size(); ++i ) {
+        res[i] = this->operator()( v[i] );
+    }
+    return res;
+}
+
+inline double Interpolation::EvalCubic1DSpline( double param[4], double x ) const {
+    return ( param[1] + 0.5 * x *
+                            ( param[2] - param[0] +
+                              x * ( 2.0 * param[0] - 5.0 * param[1] + 4.0 * param[2] - param[3] +
+                                    x * ( 3.0 * ( param[1] - param[2] ) + param[3] - param[0] ) ) ) );
+}
+
+unsigned Interpolation::IndexOfClosestValue( double value, const Vector& v ) const {
+    auto i = std::min_element( begin( v ), end( v ), [=]( double x, double y ) { return abs( x - value ) < abs( y - value ); } );
+    return std::distance( begin( v ), i );
+}

code/src/solvers/csdsnsolver.cpp

@@ -25,6 +25,13 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
     for( unsigned n = 0; n < _nEnergies; ++n ) {
         _energies[n] = energyMin + ( energyMax - energyMin ) / ( _nEnergies - 1 ) * n;
     }
+
+    std::vector<double> buffer;
+    for( auto q : _quadPoints )
+        if( q[0] > 0.0 ) buffer.push_back( q[0] );
+    std::sort( buffer.begin(), buffer.end() );
+    Vector posMu( buffer.size(), buffer.data() );
+
     // write mu grid
     Matrix muMatrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints() );
     for( unsigned l = 0; l < _settings->GetNQuadPoints(); ++l ) {
@@ -37,9 +44,30 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
         }
     }
 
-    _sigmaSE = _problem->GetScatteringXSE( _energies, muMatrix );
-    _sigmaTE = M_PI * 1e19 * _problem->GetTotalXSE( _energies );    // 5e19 *
+    _sigmaSE = std::vector<Matrix>( _energies.size(), Matrix( muMatrix.rows(), muMatrix.columns(), 0.0 ) );
+    Vector angleVec( muMatrix.columns() * muMatrix.rows() );
+    // store Matrix with mu values in vector format to call GetScatteringXS
+    for( unsigned i = 0; i < muMatrix.rows(); ++i ) {
+        for( unsigned j = 0; j < muMatrix.columns(); ++j ) {
+            angleVec[i * muMatrix.columns() + j] = std::fabs( muMatrix( i, j ) );    // icru xs go from 0 to 1 due to symmetry
+        }
+    }
+
+    ICRU database( angleVec, _energies );
+    Matrix total;
+    database.GetAngularScatteringXS( total, _sigmaTE );
+
+    // rearrange output to matrix format
+    for( unsigned n = 0; n < _energies.size(); ++n ) {
+        for( unsigned i = 0; i < muMatrix.rows(); ++i ) {
+            for( unsigned j = 0; j < muMatrix.columns(); ++j ) {
+                _sigmaSE[n]( i, j ) = total( i * muMatrix.columns() + j, n );
+            }
+        }
+    }
+
     // compute scaling s.t. scattering kernel integrates to one for chosen quadrature
+    /*
     for( unsigned n = 0; n < _nEnergies; ++n ) {
         for( unsigned p = 0; p < _nq; ++p ) {
             double cp = 0.0;
@@ -51,6 +79,7 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
             }
         }
     }
+    */
 
     _s = _problem->GetStoppingPower( _energies );
     _Q = _problem->GetExternalSource( _energies );
@@ -62,6 +91,7 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
         }
         _scatteringKernel( p, p ) = _weights[p];
     }
+
     /*
         // test 1
         for( unsigned n = 0; n < _nEnergies; ++n ) {
@@ -191,7 +221,6 @@ void CSDSNSolver::Solve() {
         // Save( n );
         dFlux   = blaze::l2Norm( fluxNew - fluxOld );
         fluxOld = fluxNew;
-        std::cout << _energies[n] << " " << dFlux << std::endl;
         if( rank == 0 ) log->info( "{:03.8f}  {:01.5e}  {:01.5e}", _energies[n], _dE / densityMin, dFlux );
         if( std::isinf( dFlux ) || std::isnan( dFlux ) ) break;
     }

code/src/solvers/csdsnsolvernotrafo.cpp

@@ -37,8 +37,29 @@ CSDSNSolverNoTrafo::CSDSNSolverNoTrafo( Config* settings ) : SNSolver( settings
         }
     }
 
-    _sigmaSE = _problem->GetScatteringXSE( _energies, muMatrix );
-    _sigmaTE = _problem->GetTotalXSE( _energies );    // M_PI * 1e19 *
+    _sigmaSE = std::vector<Matrix>( _energies.size(), Matrix( muMatrix.rows(), muMatrix.columns(), 0.0 ) );
+    Vector angleVec( muMatrix.columns() * muMatrix.rows() );
+    // store Matrix with mu values in vector format to call GetScatteringXS
+    for( unsigned i = 0; i < muMatrix.rows(); ++i ) {
+        for( unsigned j = 0; j < muMatrix.columns(); ++j ) {
+            angleVec[i * muMatrix.columns() + j] = std::fabs( muMatrix( i, j ) );    // icru xs go from 0 to 1 due to symmetry
+        }
+    }
+
+    ICRU database( angleVec, _energies );
+    Matrix total;
+    database.GetAngularScatteringXS( total, _sigmaTE );
+
+    // rearrange output to matrix format
+    for( unsigned n = 0; n < _energies.size(); ++n ) {
+        for( unsigned i = 0; i < muMatrix.rows(); ++i ) {
+            for( unsigned j = 0; j < muMatrix.columns(); ++j ) {
+                _sigmaSE[n]( i, j ) = total( i * muMatrix.columns() + j, n );
+            }
+        }
+    }
+
+    /*
     // compute scaling s.t. scattering kernel integrates to one for chosen quadrature
     for( unsigned n = 0; n < _nEnergies; ++n ) {
         for( unsigned p = 0; p < _nq; ++p ) {
@@ -51,6 +72,7 @@ CSDSNSolverNoTrafo::CSDSNSolverNoTrafo( Config* settings ) : SNSolver( settings
             }
         }
     }
+    */
 
     _s = _problem->GetStoppingPower( _energies );
     _Q = _problem->GetExternalSource( _energies );
@@ -84,12 +106,12 @@ void CSDSNSolverNoTrafo::Solve() {