Started with post-processing of physics data

code/CMakeLists.txt

@@ -117,7 +117,8 @@ include_directories( ${ParMETIS_PATH}/include )
 include_directories( ${ParMETIS_PATH}/metis/include )
 add_subdirectory( ${ParMETIS_PATH} )
 
-find_package( VTK REQUIRED COMPONENTS vtkIOGeometry vtkFiltersCore )
+find_package( VTK REQUIRED COMPONENTS vtkIOGeometry vtkFiltersCore HINTS /home/pia/vtk-inst/include/vtk-8.2)
+include(${VTK_USE_FILE})
 
 include_directories( ${CMAKE_SOURCE_DIR}/ext/cpptoml/include )
 include_directories( ${CMAKE_SOURCE_DIR}/ext/spdlog/include )

code/include/physics.h

@@ -5,34 +5,41 @@
 #include "settings/config.h"
 #include "settings/typedef.h"
 
+#include "math.h"
+#include <fstream>
+
+
 class Physics
 {
 
   private:
-    Matrix _xsH2O;
-    Matrix _totalxsH2O;
+    VectorVector _xsH2O;
+    VectorVector _xsTotalH2O;
+    VectorVector _xsTransportH2O;
+	VectorVector _stpowH2O;
     /**
-     * @brief LoadXSH2O loads (total) scattering cross sections for water from file and saves them in _xsH2O and _totalxsH2O
-     * @param fileName1 is name of cross section file
-     * @param fileName2 is name of total cross section file
+     * @brief LoadDatabase loads required physics data from file in ENDL format (default IAEA EEDL database)
+     * @param fileName is name of cross section file
      */
-    void LoadXSH2O( std::string fileName1, std::string fileName2 );
 
   public:
-    /**
+ 
+    void LoadDatabase (std::string fileName);
+    
+     /**
      * @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
      */
-    VectorVector GetScatteringXS( std::vector<double> energies, std::vector<double> density, std::vector<double> Omegas );
+    VectorVector GetScatteringXS (Vector energies,Vector density, Vector Omegas);
 
     /**
      * @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
      * @param density is vector with patient densities (at different spatial cells)
      */
-    VectorVector GetTotalXS( std::vector<double> energies, std::vector<double> density );
+    VectorVector GetTotalXS (Vector energies,Vector density);
 
     /**
      * @brief GetStoppingPower gives back vector of vectors of stopping powers for materials defined by density and energies in vector energy
@@ -40,20 +47,27 @@ class Physics
      * @param density is vector with patient densities (at different spatial cells)
      * @param sH2O is vector of stopping powers in water
      */
-    VectorVector GetStoppingPower( std::vector<double> energies, std::vector<double> sH2O );
+    VectorVector GetStoppingPower (Vector energies,Vector density,Vector sH2O);
+
+    /**
+     * @brief GetTransportXS gives back vector of vectors of stopping powers 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)
+     */
+    VectorVector GetTransportXS (Vector energies,Vector density);
 
     /**
      * @brief Physics constructor
      * @param settings stores all needed user information
      */
-    Physics( Config* settings );
+    Physics ();
 
     /**
      * @brief Create constructor
      * @param settings stores all needed information
      * @return pointer to Physics
      */
-    static Physics* Create( Config* settings );
+    static Physics* Create ();
 };
 
 #endif

code/src/physics.cpp

 #include "physics.h"
+using namespace std;
+using blaze::CompressedVector;
+
 
 Physics::Physics( Config* settings ) {
     // @TODO set up physical model from settings?
@@ -7,19 +10,67 @@ Physics::Physics( Config* settings ) {
     LoadXSH2O( "test1", "test2" );
 }
 
+Physics::Physics (){
+    //LoadXSH2O("test1","test2");
+    //LoadTransportXSH2O("test3","test4");
+
+}
+
 void Physics::LoadXSH2O( std::string fileName1, std::string fileName2 ) {
     // @TODO
 }
 
-VectorVector Physics::GetScatteringXS( std::vector<double> energies, std::vector<double> density, std::vector<double> Omegas ) {
-    // @TODO
+void Physics::LoadDatabase( std::string fileName){
+    VectorVector transport_XS_H;
+    VectorVector transport_XS_O;
+	
+	VectorVector scattering_XS_H;
+    VectorVector scattering_XS_O;
+	
+	VectorVector total_scat_XS_H;
+    VectorVector total_scat_XS_O;
+	
+	VectorVector stopp_pow_H;
+    VectorVector stopp_pow_O;
+	
+    //Read data for H
+
+
+
+    //Read data for O
+	
+	//Combine values for H and O according to mixture ratio in water 
+    Physics::_xsTransportH2O = 0.11189400*transport_XS_H + 0.88810600*transport_XS_O;
+	Physics::_xsH2O = 0.11189400*scattering_XS_H + 0.88810600*scattering_XS_O;
+	Physics::_xsTotalH2O = 0.11189400*total_scat_XS_H + 0.88810600*total_scat_XS_O;
+	Physics::_stpowH2O = 0.11189400*stopp_pow_XS_H + 0.88810600*stopp_pow_XS_O;
+}
+
+VectorVector Physics::GetScatteringXS (Vector energies,Vector density, Vector Theta){
+
     VectorVector scattering_XS;
     return scattering_XS;
 }
 
-VectorVector Physics::GetTotalXS( std::vector<double> energies, std::vector<double> density ) {
-    // @TODO
+
+VectorVector Physics::GetTotalXS (Vector energies,Vector density) {
     VectorVector total_XS;
+	
+	if (energies < 1.000000000E-05) {
+	
+	}
+	else if (energies > 1.000000000E+05) {
+		
+	}
+	else {
+	//Find upper and lower bound of energy value in table	
+	CompressedVector<int>::Iterator upper( Physics::_xsTotalH2O(0).upperBound(energies));
+	CompressedVector<int>::Iterator lower( Physics::_xsTotalH2O(0).lowerBound(energies));
+	
+	//Linear interpolation between upper and lower bound
+	
+	}
+	
     return total_XS;
 }
 
@@ -29,4 +80,11 @@ VectorVector Physics::GetStoppingPower( std::vector<double> energies, std::vecto
     return stopping_power;
 }
 
-Physics* Physics::Create( Config* settings ) { return new Physics( settings ); }
+
+VectorVector Physics::GetTransportXS (Vector energies,Vector density){
+    // @TODO
+    VectorVector transport_XS;
+    return transport_XS;
+}
+
+Physics* Physics::Create() { return new Physics(); }