diff --git a/CMakeLists.txt b/CMakeLists.txt
index d5af8f9..6a38aab 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -28,7 +28,7 @@ include_directories(externals/json/include)
 #include(${CMAKE_SOURCE_DIR}/cmake/SetBuildType.cmake)
 
 # Default to cpmpile with debugging symbols
-set(CMAKE_BUILD_TYPE DEBUG)
+set(CMAKE_BUILD_TYPE RELEASE)
 # set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
 # set (CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
 # set(BUILD_TYPE "PROFILING")
@@ -92,7 +92,7 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/lib)
 ############################################################
 # Set the build type environment variable
 ############################################################
-if(BUILD_TYPE STREQUAL "TEST")
+if(CMAKE_BUILD_TYPE STREQUAL "TEST")
   set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;--coverage;")
   set(TEST_LIBRARIES "gcov")
   message(STATUS "BUILD_TYPE SET TO TEST")
@@ -102,18 +102,18 @@ elseif(CMAKE_BUILD_TYPE STREQUAL "DEBUG")
   set(TEST_LIBRARIES "")
   message(STATUS "BUILD_TYPE SET TO DEBUG")
 
-elseif(BUILD_TYPE STREQUAL "RELEASE")
-  set(BUILD_TYPE_COMPILE_FLAGS "-O2")
+elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE")
+  set(BUILD_TYPE_COMPILE_FLAGS "-g;-O2")
   set(TEST_LIBRARIES "")
   message(STATUS "BUILD_TYPE SET TO RELEASE")
 
-elseif(BUILD_TYPE STREQUAL "PROFILING")
+elseif(CMAKE_BUILD_TYPE STREQUAL "PROFILING")
   set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;-pg")
   set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg")
   set (CMAKE_EXE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -pg")
   message(STATUS "BUILD_TYPE SET TO PROFILING")
 
-elseif(BUILD_TYPE STREQUAL "ADDRESS-SANITIZE")
+elseif(CMAKE_BUILD_TYPE STREQUAL "ADDRESS-SANITIZE")
   set(BUILD_TYPE_COMPILE_FLAGS "-g;-O0;")
   set(TEST_LIBRARIES "")  
   set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
diff --git a/inres1/aegis_settings.json b/inres1/aegis_settings.json
index 52cd7a5..ac6b718 100644
--- a/inres1/aegis_settings.json
+++ b/inres1/aegis_settings.json
@@ -1,40 +1,42 @@
 {
-  "description": "inres1 test case to compare against SMARDDA",
-  "aegis_params":{
-      "DAGMC": "inres1_shad.h5m",
-      "step_size": 0.01,
-      "max_steps": 10000,
-      "launch_pos": "mc",
-
-      "monte_carlo_params":{
-            "number_of_particles_per_facet":5,
-            "write_particle_launch_positions":false
-      },
-
-      "target_surfs": [1,2,3,4,5,6],
-      "force_no_deposition": false,
-      "coordinate_system": "flux",
-      "execution_type": "dynamic",
-
-      "dynamic_batch_params":{
-          "batch_size":16,
-          "worker_profiling": false,
-	  "debug":false
-      } 
-  },
-  "equil_params":{
-      "eqdsk": "EQ3.eqdsk",
-      "power_sol": 7.5e+06,
-      "lambda_q": 0.012,
-      "r_outrbdry": 8.07385,
-      "cenopt": 2,
-      "psiref": -2.1628794,
-      "rmove": -0.006,
-      "draw_equil_rz": false,
-      "draw_equil_xyz": false,
-      "print_debug_info": false
-  },
-  "vtk_params":{
-      "draw_particle_tracks": false
-  }
- }
+    "description": "inres1 test case to compare against SMARDDA",
+    "aegis_params":{
+        "DAGMC": "inres1_shad.h5m",
+        "step_size": 0.01,
+        "max_steps": 10000,
+        "launch_pos": "mc",
+  
+        "monte_carlo_params":{
+              "number_of_particles_per_facet":1,
+              "write_particle_launch_positions":false
+        },
+  
+        "target_surfs": [1,2,3,4,5,6],
+        "force_no_deposition": false,
+        "coordinate_system": "cart",
+        "execution_type": "dynamic",
+  
+        "dynamic_batch_params":{
+            "batch_size":16,
+            "worker_profiling": false,
+        "debug":false
+        } 
+    },
+    "equil_params":{
+        "eqdsk": "EQ3.eqdsk",
+        "power_sol": 7.5e+06,
+        "lambda_q": 0.012,
+        "r_outrbdry": 8.07385,
+        "cenopt": 2,
+        "psiref": -2.1628794,
+        "rmove": -0.006,
+        "spline_type":"cubic",
+        "draw_equil_rz": false,
+        "draw_equil_xyz": false,
+        "print_debug_info": false
+    },
+    "vtk_params":{
+        "draw_particle_tracks": false
+    }
+   }
+  
diff --git a/src/aegis_lib/CoordTransform.cpp b/src/aegis_lib/CoordTransform.cpp
index effa031..fec2cc1 100644
--- a/src/aegis_lib/CoordTransform.cpp
+++ b/src/aegis_lib/CoordTransform.cpp
@@ -31,6 +31,28 @@ CoordTransform::cart_to_polar(std::vector<double> inputVector)
   return outputVector;
 }
 
+// cart_to_polar without the calculation of phi
+std::vector<double>
+CoordTransform::cart_to_polar_xy(std::vector<double> inputVector)
+{
+  std::vector<double> outputVector(3);
+  double r;   // polar r
+  double phi; // polar phi
+  double x;   // cart x
+  double y;   // cart y
+  double z;   // cart z
+
+  x = inputVector[0];
+  y = inputVector[1];
+  z = inputVector[2];
+
+  r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
+
+  outputVector[0] = r;
+  outputVector[1] = z;
+  return outputVector;
+}
+
 std::vector<double>
 CoordTransform::polar_to_cart(std::vector<double> inputVector)
 {
@@ -182,4 +204,4 @@ CoordTransform::polar_to_flux(double e0, double e1, double e2,
   outputVector[2] = phi;
 
   return outputVector;
-}
+}
\ No newline at end of file
diff --git a/src/aegis_lib/CoordTransform.h b/src/aegis_lib/CoordTransform.h
index c43f8b1..d358350 100644
--- a/src/aegis_lib/CoordTransform.h
+++ b/src/aegis_lib/CoordTransform.h
@@ -6,6 +6,8 @@
 #include <cmath>
 #include <string>
 #include "EquilData.h"
+#include "alglib/interpolation.h"
+
 
 namespace CoordTransform
 {
@@ -13,6 +15,7 @@ namespace CoordTransform
 
 // Cartesian to polar toroidal (direction = "backwards" for polar->cartesian)
 std::vector<double> cart_to_polar(std::vector<double> inputVector);
+std::vector<double> cart_to_polar_xy(std::vector<double> inputVector);
 std::vector<double> cart_to_polar(double e0, double e1, double e2);
 
 std::vector<double> polar_to_cart(std::vector<double> inputVector);
@@ -21,7 +24,10 @@ std::vector<double> polar_to_cart(double e0, double e1, double e2);
 // Polar toroidal to flux coordinates defined by psi spline
 std::vector<double> polar_to_flux(std::vector<double> inputVector,
                                   const std::shared_ptr<EquilData>& equilibrium);
+
+
 std::vector<double> polar_to_flux(double e0, double e1, double e2, const std::shared_ptr<EquilData>& equilibrium);
+
 };
 
 #endif
\ No newline at end of file
diff --git a/src/aegis_lib/EquilData.cpp b/src/aegis_lib/EquilData.cpp
index 7d4e18e..054bb02 100644
--- a/src/aegis_lib/EquilData.cpp
+++ b/src/aegis_lib/EquilData.cpp
@@ -51,6 +51,9 @@ EquilData::read_params(const std::shared_ptr<JsonHandler> & configFile)
     drawEquRZ = equilParams.get_optional<bool>("draw_equil_rz").value_or(drawEquRZ);
     drawEquXYZ = equilParams.get_optional<bool>("draw_equil_xyz").value_or(drawEquXYZ);
     debug = equilParams.get_optional<bool>("debug").value_or(debug);
+
+    auto splineTypeStr = equilParams.get_optional<std::string>("spline_type").value();
+    _splineType = (splineTypeStr == "linear") ? SplineType::LINEAR : SplineType::CUBIC;
   }
   else
   {
@@ -438,7 +441,7 @@ EquilData::init_interp_splines()
 
   // set_rsig();
 
-  // loop over Z creating spline knots
+  // loop over R creating spline knots
   for (int i = 1; i < nw; i++)
   {
     rPts[i] = rPts[i - 1] + dr;
@@ -481,24 +484,62 @@ EquilData::init_interp_splines()
   psi_1dgrid.setcontent(nw, psi1dpts.data());
   f_grid.setcontent(nw, eqdsk.fpol.data());
 
-  //  for (int i=0 ; i<nh; i++)
-  //  {
-  //    psi1d_out << psi_1dgrid[i] << std::endl;
-  //  }
-
   // Construct the spline interpolant for flux function psi(R,Z)
-  alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline);
+  switch (_splineType)
+  {
+    case SplineType::CUBIC:
+      alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline);
+      alglib::spline1dbuildcubic(psi_1dgrid, f_grid, fSpline);
+
+      log_string(LogLevel::WARNING,
+                 "Splines initialised as cubic splines (cubic splines are default, use "
+                 "spline_type:'linear' to switch to linear splines)");
+      break;
+
+    case SplineType::LINEAR:
+      alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, psi_grid, 1, psiSpline); // psi spline
+      alglib::spline1dbuildlinear(psi_1dgrid, f_grid, fSpline); // I(psi) aka f(psi)
 
-  // Construct the spline interpolant for toroidal component flux function
-  // I(psi) aka f(psi)
-  alglib::spline1dbuildcubic(psi_1dgrid, f_grid, fSpline);
+      log_string(LogLevel::WARNING, "Splines intialised as linear splines");
+      break;
+  }
+
+  int count2 = 0;
+  std::vector<double> dpsidrPts(nw * nh);
+  std::vector<double> dpsidzPts(nw * nh);
+  double psi, dpsidr, dpsidz, dpsi2drdz;
+  for (int j = 0; j < nh; j++) // flatten Psi(R,Z) array into Psi(index)
+  {
+    for (int i = 0; i < nw; i++)
+    {
+      alglib::spline2ddiff(psiSpline, rPts[i], zPts[j], psi, dpsidr, dpsidz, dpsi2drdz);
+      dpsidrPts[count2] = dpsidr;
+      dpsidzPts[count2] = dpsidz;
+      count2 += 1;
+    }
+  }
+
+  alglib::real_1d_array dpsidr_grid;
+  alglib::real_1d_array dpsidz_grid;
+  dpsidr_grid.setcontent(nw * nh, dpsidrPts.data());
+  dpsidz_grid.setcontent(nw * nh, dpsidzPts.data());
+
+  switch (_splineType)
+  {
+    case SplineType::CUBIC:
+      alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, dpsidr_grid, 1, dpsidrSpline);
+      alglib::spline2dbuildbicubicv(r_grid, nw, z_grid, nh, dpsidz_grid, 1, dpsidzSpline);
+      break;
+
+    case SplineType::LINEAR:
+      alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, dpsidr_grid, 1, dpsidrSpline);
+      alglib::spline2dbuildbilinearv(r_grid, nw, z_grid, nh, dpsidz_grid, 1, dpsidzSpline);
+      break;
+  }
 
   // Alglib::spline2ddiff function can return a value of spline at (R,Z) as well
   // as derivatives. I.e no need to have a separate spline for each derivative
   // dPsidR and dPsidZ
-
-  // create 1d spline for f(psi) aka eqdsk.fpol
-  // alglib::spline1dbuildlinear(f_grid,)
 }
 
 // Write out psi(R,Z) data for gnuplotting
@@ -746,7 +787,7 @@ EquilData::r_extrema()
         work1[j] = 2;
         break;
       }
-      zpsi = alglib::spline2dcalc(psiSpline, re, ze);
+      zpsi = get_psi(re, ze);
     }
   }
 
@@ -825,7 +866,7 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
   else
   {
     std::vector<double> polarPosition;
-    polarPosition = CoordTransform::cart_to_polar(position);
+    polarPosition = CoordTransform::cart_to_polar_xy(position);
     zr = polarPosition[0];
     zz = polarPosition[1];
   }
@@ -836,7 +877,10 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
   }
 
   // evaluate psi and psi derivs at given (R,Z) coords
-  alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
+  zpsi = alglib::spline2dcalc(psiSpline, zr, zz);
+  zdpdr = alglib::spline2dcalc(dpsidrSpline, zr, zz);
+  zdpdz = alglib::spline2dcalc(dpsidzSpline, zr, zz);
+  // alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
 
   // evaluate I aka f at psi (I aka f is the flux function)
   zf = alglib::spline1dcalc(fSpline, zpsi);
@@ -1074,7 +1118,7 @@ EquilData::boundary_rb()
       // extremal Z reached
       break;
     }
-    zpsi = alglib::spline2dcalc(psiSpline, re, ze);
+    zpsi = get_psi(re, ze);
 
     if (rsig > 0) // psiaxis < psiqbdry (psi increasing outwards)
     {
@@ -1115,6 +1159,7 @@ EquilData::boundary_rb()
 
     // check monotone
     alglib::spline2ddiff(psiSpline, re, ze, zpsi, zdpdr, zdpdz, zddpdz);
+
     zdpdsr = zdpdr * zcostheta + zdpdz * zsintheta;
     if (idplset == 1)
     {
@@ -1305,13 +1350,12 @@ EquilData::psi_limiter(std::vector<std::vector<double>> vertices)
   for (auto & i : vertices)
   {
     std::vector<double> cartPos = i;
-    std::vector<double> polarPos = CoordTransform::cart_to_polar(i);
+    std::vector<double> polarPos = CoordTransform::cart_to_polar_xy(i);
 
     zr = polarPos[0];
     zz = polarPos[1];
 
-    zpsi = alglib::spline2dcalc(this->psiSpline, zr,
-                                zz); // spline interpolation of psi(R,Z)
+    zpsi = get_psi(zr, zz);
 
     ztheta = atan2(zz - zcen, zr - rcen);
     if (ztheta < -M_PI_2)
@@ -1338,7 +1382,10 @@ EquilData::psi_limiter(std::vector<std::vector<double>> vertices)
 
   double zdpdr;
   double zdpdz;
-  alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
+  // alglib::spline2ddiff(psiSpline, zr, zz, zpsi, zdpdr, zdpdz, null);
+  zpsi = alglib::spline2dcalc(psiSpline, zr, zz);
+  zdpdr = alglib::spline2dcalc(dpsidrSpline, zr, zz);
+  zdpdz = alglib::spline2dcalc(dpsidzSpline, zr, zz);
 
   double zbpbdry, zbtotbdry;
   zbpbdry = (1 / zr) * sqrt(std::max(0.0, (pow(zdpdr, 2) + pow(zdpdz, 2))));
@@ -1409,7 +1456,7 @@ EquilData::get_midplane_params()
 bool
 EquilData::check_if_in_bfield(std::vector<double> xyzPos)
 {
-  std::vector<double> polarPos = CoordTransform::cart_to_polar(xyzPos);
+  std::vector<double> polarPos = CoordTransform::cart_to_polar_xy(xyzPos);
   double r = polarPos[0];
   double z = polarPos[1];
   if (r < rmin || r > rmax || z < zmin || z > zmax)
@@ -1418,3 +1465,9 @@ EquilData::check_if_in_bfield(std::vector<double> xyzPos)
   }
   return true;
 }
+
+double
+EquilData::get_psi(double r, double z)
+{
+  return -(alglib::spline2dcalc(psiSpline, r, z)); // spline returns sign flipped psi
+}
\ No newline at end of file
diff --git a/src/aegis_lib/EquilData.h b/src/aegis_lib/EquilData.h
index 7e6565e..6ec578e 100644
--- a/src/aegis_lib/EquilData.h
+++ b/src/aegis_lib/EquilData.h
@@ -10,6 +10,12 @@
 #include "SimpleLogger.h"
 #include "AegisBase.h"
 
+enum class SplineType
+{
+  CUBIC,
+  LINEAR
+};
+
 
 struct eqdskData
 {
@@ -91,6 +97,9 @@ class EquilData : public AegisBase
   // Return eqdsk struct
   eqdskData get_eqdsk_struct();
 
+  // return psi value at given (R,Z) coords
+  double get_psi(double r, double z);
+
   // Read eqdsk file
   void read_eqdsk(std::string filename);
 
@@ -202,7 +211,12 @@ class EquilData : public AegisBase
   alglib::spline2dinterpolant psiSpline; // 2d spline interpolant for Psi(R,Z)
   alglib::spline1dinterpolant fSpline; // 1d spline interpolant for f(psi) or I(psi) toroidal component
 
+  alglib::spline2dinterpolant dpsidrSpline; // 2d spline interpolant for dpsidr(R,Z)
+  alglib::spline2dinterpolant dpsidzSpline; // 2d spline interpolant for dpsidz(R,Z)
+  
+
   private:
+  SplineType _splineType = SplineType::CUBIC; 
 
   bool debug = false;
   bool drawEquRZ = false;
diff --git a/src/aegis_lib/Particle.cpp b/src/aegis_lib/Particle.cpp
index d02de88..9d12f49 100644
--- a/src/aegis_lib/Particle.cpp
+++ b/src/aegis_lib/Particle.cpp
@@ -73,21 +73,17 @@ ParticleBase::get_xyz_pos()
 double
 ParticleBase::get_psi(const std::shared_ptr<EquilData> & equilibrium)
 {
-  std::vector<double> fluxPos;
   std::vector<double> polarPos;
   double psi;
   switch (coordSystem)
   {
     case coordinateSystem::CARTESIAN:
-      polarPos = CoordTransform::cart_to_polar(pos);
-      fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium);
-      psi = fluxPos[0];
+      polarPos = CoordTransform::cart_to_polar_xy(pos);
+      psi = equilibrium->get_psi(polarPos[0], polarPos[1]);
       break;
     case coordinateSystem::POLAR:
-      fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium);
-      psi = fluxPos[0];
+      psi = equilibrium->get_psi(pos[0], pos[1]);
   }
-
   return psi;
 }
 
diff --git a/src/aegis_lib/ParticleSimulation.cpp b/src/aegis_lib/ParticleSimulation.cpp
index 4481bf7..eb24300 100644
--- a/src/aegis_lib/ParticleSimulation.cpp
+++ b/src/aegis_lib/ParticleSimulation.cpp
@@ -64,7 +64,6 @@ ParticleSimulation::Execute_dynamic_mpi()
 
   MPI_Status mpiStatus;
 
-
   const int noMoreWork = -1;
 
   std::vector<double> particleHeatfluxes(num_particles_launched());
@@ -442,8 +441,7 @@ ParticleSimulation::init_geometry()
   }
   dagmcMeshReadTime = MPI_Wtime() - dagmcMeshReadStart;
 
-
-  // setup B Field data 
+  // setup B Field data
   double prepSurfacesStart = MPI_Wtime();
 
   std::vector<double> vertexCoordinates;
@@ -483,7 +481,7 @@ ParticleSimulation::init_geometry()
   targetFacets = select_target_surface();
   integrator = std::make_unique<SurfaceIntegrator>(targetFacets);
   prepSurfacesTime = MPI_Wtime() - prepSurfacesStart;
-  
+
   double setupArrayOfParticlesTimeStart = MPI_Wtime();
   setup_sources();
   setupArrayOfParticlesTime = MPI_Wtime() - setupArrayOfParticlesTimeStart;
@@ -819,7 +817,6 @@ ParticleSimulation::Execute_serial()
   vtkInterface->init();
   moab::ErrorCode rval;
 
-
   int start = 0;
   int end = num_particles_launched();
 
@@ -870,7 +867,6 @@ ParticleSimulation::Execute_serial()
   vtkInterface->write_multiBlockData("particle_tracks.vtm");
 
   print_particle_stats(particleStats);
-
 }
 
 void
@@ -882,7 +878,6 @@ ParticleSimulation::Execute_mpi()
   vtkInterface->init();
   MPI_Status mpiStatus;
 
-
   int totalFacets = target_num_facets();
   int nFacetsPerProc = totalFacets / nprocs;
   int remainder = totalFacets % nprocs;
@@ -955,7 +950,6 @@ ParticleSimulation::Execute_mpi()
     attach_mesh_attribute("Heatflux", targetFacets, rootQvalues);
     write_out_mesh(meshWriteOptions::BOTH, targetFacets);
     mainParticleLoopTime = MPI_Wtime() - mainParticleLoopStart;
-
   }
 
   mpi_particle_stats();
@@ -1171,9 +1165,12 @@ ParticleSimulation::get_profiling_times()
 {
   std::map<std::string, double> profilingTimes;
   profilingTimes.insert(std::make_pair("DAGMC Mesh read runtime = ", dagmcMeshReadTime));
-  profilingTimes.insert(std::make_pair("Preparing surfaces for particles runtime = ", prepSurfacesTime));  
-  profilingTimes.insert(std::make_pair("Pool of particles generation runtime = ", setupArrayOfParticlesTime));
-  profilingTimes.insert(std::make_pair("Main particle tracking loop runtime = ", mainParticleLoopTime));
+  profilingTimes.insert(
+      std::make_pair("Preparing surfaces for particles runtime = ", prepSurfacesTime));
+  profilingTimes.insert(
+      std::make_pair("Pool of particles generation runtime = ", setupArrayOfParticlesTime));
+  profilingTimes.insert(
+      std::make_pair("Main particle tracking loop runtime = ", mainParticleLoopTime));
   profilingTimes.insert(std::make_pair("Mesh Write out runtime = ", aegisMeshWriteTime));
 
   return profilingTimes;
diff --git a/src/aegis_lib/Source.cpp b/src/aegis_lib/Source.cpp
index 944ca8e..e0a10b7 100644
--- a/src/aegis_lib/Source.cpp
+++ b/src/aegis_lib/Source.cpp
@@ -101,7 +101,6 @@ Sources::set_heatflux_params(const std::shared_ptr<EquilData> & equilibrium,
   std::vector<double> bFieldXYZ;
 
   polarPos = CoordTransform::cart_to_polar(launchPos);
-  fluxPos = CoordTransform::polar_to_flux(polarPos, equilibrium);
 
   bField = equilibrium->b_field(polarPos, "forwards");
 
@@ -114,7 +113,7 @@ Sources::set_heatflux_params(const std::shared_ptr<EquilData> & equilibrium,
 
   Bn = product; // store B.n
 
-  psi = fluxPos[0]; // store psi at particle start
+  psi = equilibrium->get_psi(polarPos[0], polarPos[1]); // store psi at particle start
   psid = psi + equilibrium->psibdry;
   _totalHeatflux = equilibrium->omp_power_dep(psid, Bn, "exp"); // store Q at particle start
 }