first bug in scattering fixed

code/include/physics.h

@@ -31,11 +31,17 @@ class Physics
     /** @brief GetScatteringXS gives back vector of vectors of scattering cross sections for materials defined by density and energies in vector
      * energy
      * @param energies is vector with energies
-     * @param density is vector with patient densities (at different spatial cells)
-     * @param Omega are scattering angles
+     * @param angle are scattering angles (mu)
      */
     VectorVector GetScatteringXS( Vector energies, Vector angle );
 
+    /** @brief GetScatteringXS gives back vector of Matrix of scattering cross sections for materials defined by density and energies in vector
+     * energy
+     * @param energies is vector with energies
+     * @param angle is matrix of scattering angles (<Omega_i,Omega_j>)
+     */
+    std::vector<Matrix> GetScatteringXS( const Vector& energies, const Matrix& angle );
+
     /**
      * @brief GetTotalXS gives back vector of vectors of total cross sections for materials defined by density and energies in vector energy
      * @param energies is vector with energies

code/include/problems/electronrt.h

@@ -22,7 +22,7 @@ class ElectronRT : public ProblemBase
 
     virtual VectorVector GetScatteringXS( const Vector& energies );
     virtual VectorVector GetTotalXS( const Vector& energies );
-    virtual VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles );
+    virtual std::vector<Matrix> GetScatteringXSE( const Vector& energies, const Matrix& angles );
     virtual Vector GetTotalXSE( const Vector& energies );
     virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
     virtual VectorVector SetupIC();

code/include/problems/linesource.h

@@ -37,7 +37,7 @@ class LineSource_SN_Pseudo1D_Physics : public LineSource_SN_Pseudo1D
   public:
     LineSource_SN_Pseudo1D_Physics( Config* settings, Mesh* mesh );
 
-    VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles ) override;
+    std::vector<Matrix> GetScatteringXSE( const Vector& energies, const Matrix& angles ) override;
     Vector GetTotalXSE( const Vector& energies ) override;
 };
 

code/include/problems/problembase.h

@@ -49,8 +49,8 @@ class ProblemBase
      *        in vector energy
      * @param energy is the energy the cross section is queried for
      */
-    virtual VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles ) {
-        return VectorVector( energies.size(), Vector( angles.size() ) );
+    virtual std::vector<Matrix> GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
+        return std::vector<Matrix>( energies.size(), Matrix( angles.rows(), angles.columns() ) );
     }
 
     /**

code/include/solvers/csdsnsolver.h

@@ -15,8 +15,8 @@ class CSDSNSolver : public SNSolver
     Vector _angle;    /*! @brief: angles for SN */
     Vector _density;  /*! @brief: patient density for each grid cell */
 
-    VectorVector _sigmaSE; /*!  @brief scattering cross section for all energies*/
-    Vector _sigmaTE;       /*!  @brief total cross section for all energies*/
+    std::vector<Matrix> _sigmaSE; /*!  @brief scattering cross section for all energies*/
+    Vector _sigmaTE;              /*!  @brief total cross section for all energies*/
 
   public:
     /**

code/src/physics.cpp

@@ -75,6 +75,28 @@ void Physics::LoadDatabase( std::string fileName_H, std::string fileName_O, std:
     _stpowH2O = ReadStoppingPowers( fileName_stppower );
 }
 
+std::vector<Matrix> Physics::GetScatteringXS( const Vector& energies, const Matrix& angle ) {
+    std::vector<Matrix> out( energies.size(), Matrix( angle.rows(), angle.columns() ) );
+    Vector angleVec( angle.columns() * angle.rows() );
+    // store Matrix with mu values in vector format to call GetScatteringXS
+    for( unsigned i = 0; i < angle.rows(); ++i ) {
+        for( unsigned j = 0; j < angle.columns(); ++j ) {
+            angleVec[i * angle.columns() + j] = angle( i, j );
+        }
+    }
+    VectorVector outVec = GetScatteringXS( energies, angleVec );
+
+    // rearrange output to matrix format
+    for( unsigned n = 0; n < energies.size(); ++n ) {
+        for( unsigned i = 0; i < angle.rows(); ++i ) {
+            for( unsigned j = 0; j < angle.columns(); ++j ) {
+                out[n]( i, j ) = outVec[n][i * angle.columns() + j];
+            }
+        }
+    }
+    return out;
+}
+
 VectorVector Physics::GetScatteringXS( Vector energies, Vector angle ) {
     std::vector<std::vector<double>> tmpH, tmpO;          // vectorvector which stores data at fixed energies
     std::vector<std::vector<double>> xsHGrid, xsOGrid;    // matrix which stores tensorized data for given angular grid, original energy grid

code/src/problems/electronrt.cpp

@@ -18,7 +18,7 @@ VectorVector ElectronRT::GetTotalXS( const Vector& energies ) {
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
-VectorVector ElectronRT::GetScatteringXSE( const Vector& energies, const Vector& angles ) {
+std::vector<Matrix> ElectronRT::GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
     // @TODO
     return _physics->GetScatteringXS( energies, angles );
 }

code/src/problems/linesource.cpp

@@ -50,7 +50,7 @@ LineSource_SN_Pseudo1D_Physics::LineSource_SN_Pseudo1D_Physics( Config* settings
     _physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" );
 }
 
-VectorVector LineSource_SN_Pseudo1D_Physics::GetScatteringXSE( const Vector& energies, const Vector& angles ) {
+std::vector<Matrix> LineSource_SN_Pseudo1D_Physics::GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
     return _physics->GetScatteringXS( energies, angles );
 }
 

code/src/solvers/csdsnsolver.cpp

@@ -26,11 +26,18 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
         _energies[n] = energyMin + ( energyMax - energyMin ) / ( _nEnergies - 1 ) * n;
     }
     // write mu grid
-    for( unsigned k = 0; k < _settings->GetNQuadPoints(); ++k ) {
-        _angle[k] = _quadPoints[k][2];
+    Matrix muMatrix( _settings->GetNQuadPoints(), _settings->GetNQuadPoints() );
+    for( unsigned l = 0; l < _settings->GetNQuadPoints(); ++l ) {
+        for( unsigned k = 0; k < _settings->GetNQuadPoints(); ++k ) {
+            double inner = 0.0;
+            for( unsigned j = 0; j < 3; ++j ) {
+                inner += _quadPoints[l][j] * _quadPoints[k][j];    // compute mu at Omega_l*Omega_k
+            }
+            muMatrix( l, k ) = inner;
+        }
     }
 
-    _sigmaSE = _problem->GetScatteringXSE( _energies, _angle );
+    _sigmaSE = _problem->GetScatteringXSE( _energies, muMatrix );
     _sigmaTE = _problem->GetTotalXSE( _energies );
     _s       = _problem->GetStoppingPower( _energies );
     _Q       = _problem->GetExternalSource( _energies );
@@ -105,7 +112,7 @@ void CSDSNSolver::Solve() {
                 psiNew[j][i] = _sol[j][i] - ( _dE / _areas[j] ) * psiNew[j][i] - _dE * _sigmaTE[n] * _sol[j][i];
             }
             // compute scattering effects (_scatteringKernel is simply multiplication with quad weights)
-            psiNew[j] += _dE * ( blaze::trans( _sigmaSE[n] ) * _scatteringKernel * _sol[j] );    // multiply scattering matrix with psi
+            psiNew[j] += _dE * ( _sigmaSE[n] * _scatteringKernel * _sol[j] );    // multiply scattering matrix with psi
             // TODO: Check if _sigmaS^T*psi is correct
 
             // TODO: figure out a more elegant way