diff --git a/CMakeLists.txt b/CMakeLists.txt
index 48765a9f..07949a9c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.5)
 
-project(Windows-CalcEngine VERSION 1.0.32 LANGUAGES CXX)
+project(Windows-CalcEngine VERSION 1.0.33 LANGUAGES CXX)
 
 set(target_name ${PROJECT_NAME})
 
diff --git a/src/Common/src/MatrixSeries.cpp b/src/Common/src/MatrixSeries.cpp
index 695f8fe9..1cca1525 100644
--- a/src/Common/src/MatrixSeries.cpp
+++ b/src/Common/src/MatrixSeries.cpp
@@ -1,16 +1,20 @@
 #include <cassert>
 #include <stdexcept>
 
+#include <thread>
+
 #include "MatrixSeries.hpp"
 #include "SquareMatrix.hpp"
 #include "Series.hpp"
 #include "IntegratorStrategy.hpp"
+#include "Utility.hpp"
 
 
 namespace FenestrationCommon
 {
     CMatrixSeries::CMatrixSeries(const size_t t_Size1, const size_t t_Size2, size_t seriesSize) :
-        m_Size1(t_Size1), m_Size2(t_Size2)
+        m_Size1(t_Size1),
+        m_Size2(t_Size2)
     {
         m_Matrix = std::vector<std::vector<CSeries>>(m_Size1);
         for(size_t i = 0; i < m_Size1; ++i)
@@ -106,6 +110,7 @@ namespace FenestrationCommon
 
     void CMatrixSeries::mMult(const CSeries & t_Series)
     {
+        // Parallelization here did not show any improvements. Program was performing even slower.
         for(size_t i = 0; i < m_Matrix.size(); ++i)
         {
             for(size_t j = 0; j < m_Matrix[i].size(); ++j)
@@ -118,11 +123,11 @@ namespace FenestrationCommon
 
     void CMatrixSeries::mMult(const std::vector<CSeries> & t_Series)
     {
+        // Parallelization here did not show any improvements. Program was performing even slower.
         for(size_t i = 0; i < m_Matrix.size(); ++i)
         {
             for(size_t j = 0; j < m_Matrix[i].size(); ++j)
             {
-                // assert( t_Series[ i ]->size() == ( *m_Matrix[ i ][ j ] ).size() );
                 m_Matrix[i][j] = m_Matrix[i][j] * t_Series[i];
             }
         }
@@ -139,21 +144,57 @@ namespace FenestrationCommon
     {
         for(size_t i = 0; i < m_Matrix.size(); ++i)
         {
-            for(size_t j = 0; j < m_Matrix[i].size(); ++j)
+            auto numberOfThreads{1u};
+#if MULTITHREADING
+            numberOfThreads = std::thread::hardware_concurrency();
+#endif
+            const auto chunks{
+              FenestrationCommon::chunkIt(0u, m_Matrix[i].size() - 1u, numberOfThreads)};
+
+            std::vector<std::thread> workers;
+            for(const auto & chunk : chunks)
+            {
+                workers.emplace_back([&]() {
+                    for(size_t j = chunk.start; j < chunk.end; ++j)
+                    {
+                        m_Matrix[i][j] = m_Matrix[i][j].integrate(
+                          t_Integration, normalizationCoefficient, integrationPoints);
+                    }
+                });
+            }
+
+            for(auto & worker : workers)
             {
-                m_Matrix[i][j] = m_Matrix[i][j].integrate(
-                  t_Integration, normalizationCoefficient, integrationPoints);
+                worker.join();
             }
         }
     }
 
     void CMatrixSeries::interpolate(const std::vector<double> & t_Wavelengths)
     {
-        for(size_t i = 0; i < m_Matrix.size(); ++i)
+        for(size_t i = 0u; i < m_Matrix.size(); ++i)
         {
-            for(size_t j = 0; j < m_Matrix[i].size(); ++j)
+            auto numberOfThreads{1u};
+#if MULTITHREADING
+            numberOfThreads = std::thread::hardware_concurrency();
+#endif
+            const auto chunks{
+              FenestrationCommon::chunkIt(0u, m_Matrix[i].size() - 1u, numberOfThreads)};
+
+            std::vector<std::thread> workers;
+            for(const auto & chunk : chunks)
+            {
+                workers.emplace_back([&]() {
+                    for(size_t j = chunk.start; j < chunk.end; ++j)
+                    {
+                        m_Matrix[i][j] = m_Matrix[i][j].interpolate(t_Wavelengths);
+                    }
+                });
+            }
+
+            for(auto & worker : workers)
             {
-                m_Matrix[i][j] = m_Matrix[i][j].interpolate(t_Wavelengths);
+                worker.join();
             }
         }
     }
diff --git a/src/MultiLayerOptics/src/EquivalentBSDFLayer.cpp b/src/MultiLayerOptics/src/EquivalentBSDFLayer.cpp
index 8e272de9..189d2364 100644
--- a/src/MultiLayerOptics/src/EquivalentBSDFLayer.cpp
+++ b/src/MultiLayerOptics/src/EquivalentBSDFLayer.cpp
@@ -102,15 +102,6 @@ namespace MultiLayerOptics
         return m_Layer.size();
     }
 
-    void CEquivalentBSDFLayer::setMatrixLayerWavelengths(const std::vector<double> & wavelenghts)
-    {
-        m_CombinedLayerWavelengths = wavelenghts;
-        for(const auto & layer : m_Layer)
-        {
-            layer->setBandWavelengths(wavelenghts);
-        }
-    }
-
     void CEquivalentBSDFLayer::calculate()
     {
         const size_t matrixSize = m_Lambda.size();
@@ -134,12 +125,6 @@ namespace MultiLayerOptics
 
     void CEquivalentBSDFLayer::calculateWavelengthByWavelengthProperties()
     {
-        std::vector<size_t> wavelengthIndexes;
-        for(size_t i = 0; i < m_CombinedLayerWavelengths.size(); ++i)
-        {
-            wavelengthIndexes.push_back(i);
-        }
-
         std::mutex absorptanceMutex;
         std::mutex jscMutex;
         std::mutex totMutex;
diff --git a/src/MultiLayerOptics/src/EquivalentBSDFLayer.hpp b/src/MultiLayerOptics/src/EquivalentBSDFLayer.hpp
index 948b58d8..aa107293 100644
--- a/src/MultiLayerOptics/src/EquivalentBSDFLayer.hpp
+++ b/src/MultiLayerOptics/src/EquivalentBSDFLayer.hpp
@@ -32,43 +32,25 @@ namespace MultiLayerOptics
         [[nodiscard]] double getMaxLambda() const;
 
         // Absorptance wavelength by wavelength matrices
-        FenestrationCommon::CMatrixSeries getTotalA(const FenestrationCommon::Side t_Side);
+        FenestrationCommon::CMatrixSeries getTotalA(FenestrationCommon::Side t_Side);
 
         // Photovoltaic current (scaled to income irradiance equal to one)
         FenestrationCommon::CMatrixSeries getTotalJSC(FenestrationCommon::Side t_Side);
 
         // Transmittance and reflectance wavelength by wavelength matrices
-        FenestrationCommon::CMatrixSeries
-          getTotal(const FenestrationCommon::Side t_Side,
-                   const FenestrationCommon::PropertySimple t_Property);
+        FenestrationCommon::CMatrixSeries getTotal(FenestrationCommon::Side t_Side,
+                                                   FenestrationCommon::PropertySimple t_Property);
 
         void setSolarRadiation(FenestrationCommon::CSeries & t_SolarRadiation);
 
         [[nodiscard]] std::vector<std::shared_ptr<SingleLayerOptics::CBSDFLayer>> & getLayers();
         [[nodiscard]] size_t numberOfLayers() const;
 
-        void setMatrixLayerWavelengths(const std::vector<double> & wavelenghts);
-
-    private:
-        struct wavelenghtData
-        {
-            wavelenghtData(double wl, CEquivalentBSDFLayerSingleBand & layerWl) :
-                wavelength(wl),
-                layer(layerWl)
-            {}
-
-            double wavelength;
-            CEquivalentBSDFLayerSingleBand & layer;
-            std::map<std::pair<FenestrationCommon::Side, size_t>, std::vector<double>> totA;
-            std::map<std::pair<FenestrationCommon::Side, size_t>, std::vector<double>> totJSC;
-            std::map<std::pair<FenestrationCommon::Side, FenestrationCommon::PropertySimple>,
-                     FenestrationCommon::SquareMatrix>
-              tot;
-        };
-
         void calculate();
 
-        CEquivalentBSDFLayerSingleBand getEquivalentLayerAtWavelength(size_t wavelengthIndex) const;
+    private:
+        [[nodiscard]] CEquivalentBSDFLayerSingleBand
+          getEquivalentLayerAtWavelength(size_t wavelengthIndex) const;
 
         static std::vector<double> unionOfLayerWavelengths(
           const std::vector<std::shared_ptr<SingleLayerOptics::CBSDFLayer>> & t_Layer);
@@ -97,7 +79,6 @@ namespace MultiLayerOptics
         bool m_Calculated;
 
         void calculateWavelengthByWavelengthProperties();
-        
     };
 
 }   // namespace MultiLayerOptics
diff --git a/src/MultiLayerOptics/src/MultiPaneBSDF.cpp b/src/MultiLayerOptics/src/MultiPaneBSDF.cpp
index 88236c5a..29b1cc2f 100644
--- a/src/MultiLayerOptics/src/MultiPaneBSDF.cpp
+++ b/src/MultiLayerOptics/src/MultiPaneBSDF.cpp
@@ -1,11 +1,9 @@
 #include <cmath>
 #include <numeric>
-#include <algorithm>
 #include <cassert>
 #include <utility>
 
 #include "MultiPaneBSDF.hpp"
-#include "EquivalentBSDFLayer.hpp"
 #include "EquivalentBSDFLayerSingleBand.hpp"
 #include "WCESingleLayerOptics.hpp"
 #include "WCECommon.hpp"
@@ -112,6 +110,7 @@ namespace MultiLayerOptics
                 m_IncomingSolar.push_back(iTotalSolar.sum(minLambda, maxLambda));
             }
 
+
             // Produce local results matrices for each side and property
             std::map<std::pair<Side, PropertySimple>, SquareMatrix> aResults;