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..34b08d6 100644
--- a/inres1/aegis_settings.json
+++ b/inres1/aegis_settings.json
@@ -1,40 +1,41 @@
 {
-  "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": "fine_inres1.h5m",
+        "step_size": 0.01,
+        "max_steps": 10000,
+        "launch_pos": "fixed",
+  
+        "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": "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
+    }
+   }
+  
diff --git a/src/aegis_lib/CoordTransform.cpp b/src/aegis_lib/CoordTransform.cpp
index effa031..c5690a2 100644
--- a/src/aegis_lib/CoordTransform.cpp
+++ b/src/aegis_lib/CoordTransform.cpp
@@ -9,177 +9,85 @@
 #include "SimpleLogger.h"
 
 std::vector<double>
-CoordTransform::cart_to_polar(std::vector<double> inputVector)
+CoordTransform::cart_to_polar(const 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
-  phi = atan2(-y, x);              // calculate phi
-
-  outputVector[0] = r;
-  outputVector[1] = z;
-  outputVector[2] = phi;
+  double x = inputVector[0];
+  double y = inputVector[1];
+  double z = inputVector[2];
+  double r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
+  double phi = atan2(-y, x);              // calculate phi
+  std::vector<double> outputVector = {r, z, phi};
   return outputVector;
 }
 
 std::vector<double>
-CoordTransform::polar_to_cart(std::vector<double> inputVector)
+CoordTransform::polar_to_cart(const 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
-
-  r = inputVector[0];
-  z = inputVector[1];
-  phi = inputVector[2];
-
-  x = r * cos(phi);  // calculate x
-  y = -r * sin(phi); // calculate y
-
-  outputVector[0] = x;
-  outputVector[1] = y;
-  outputVector[2] = z;
-
+  double r = inputVector[0];
+  double z = inputVector[1];
+  double phi = inputVector[2];
+  double x = r * cos(phi);  // calculate x
+  double y = -r * sin(phi); // calculate y
+  std::vector<double> outputVector = {x, y, z};
   return outputVector;
 }
 
 std::vector<double>
-CoordTransform::cart_to_polar(double e0, double e1, double e2)
+CoordTransform::cart_to_polar(const double & e0, const double & e1, const double & e2)
 {
-  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 = e0;
-  y = e1;
-  z = e2;
-
-  r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
-  phi = atan2(-y, x);              // calculate phi
-
-  outputVector[0] = r;
-  outputVector[1] = z;
-  outputVector[2] = phi;
+  double x = e0;
+  double y = e1;
+  double z = e2;
+  double r = sqrt(pow(x, 2) + pow(y, 2)); // calculate
+  double phi = atan2(-y, x);              // calculate phi
+  std::vector<double> outputVector = {r, z, phi};
   return outputVector;
 }
 
 std::vector<double>
-CoordTransform::polar_to_cart(double e0, double e1, double e2)
+CoordTransform::polar_to_cart(const double & e0, const double & e1, const double & e2)
 {
-  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
-
-  r = e0;
-  z = e1;
-  phi = e2;
-
-  x = r * cos(phi); // calculate x
-
-  y = -r * sin(phi); // calculate y
-
-  outputVector[0] = x;
-  outputVector[1] = y;
-  outputVector[2] = z;
-
+  double r = e0;
+  double z = e1;
+  double phi = e2;
+  double x = r * cos(phi);  // calculate x
+  double y = -r * sin(phi); // calculate y
+  std::vector<double> outputVector = {x, y, z};
   return outputVector;
 }
 
 std::vector<double>
-CoordTransform::polar_to_flux(std::vector<double> inputVector,
+CoordTransform::polar_to_flux(const std::vector<double> & inputVector,
                               const std::shared_ptr<EquilData> & equilibrium)
 {
-  std::vector<double> outputVector(3);
-  double r;     // local polar r
-  double z;     // local polar z
-  double phi;   // local polar phi
-  double psi;   // local psi
-  double theta; // local theta
-
-  r = inputVector[0];
-  z = inputVector[1];
-  phi = inputVector[2];
-
-  psi = alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
-
-  theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
+  double r = inputVector[0];
+  double z = inputVector[1];
+  double phi = inputVector[2];
+  double psi =
+      alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
+  double theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
   if (theta < -M_PI_2)
   {
     theta = 2 * M_PI + theta;
   }
-
-  outputVector[0] = -psi;
-  outputVector[1] = theta;
-  outputVector[2] = phi;
-
+  std::vector<double> outputVector = {-psi, theta, phi};
   return outputVector;
 }
 
-// TODO
-// std::vector<double>
-// CoordTransform::flux_to_polar(std::vector<double> inputVector,
-//                               const std::shared_ptr<EquilData> & equilibrium)
-// {
-//   std::vector<double> outputVector(3);
-//   double r;     // local polar r
-//   double z;     // local polar z
-//   double phi;   // local polar phi
-//   double psi;   // local psi
-//   double theta; // local theta
-
-//   if (direction == "backwards") // backwards transform (flux -> polar coords)
-//   {
-//     psi = inputVector[0];
-//     theta = inputVector[1];
-//     phi = inputVector[2];
-//   }
-
-//   return outputVector;
-// }
-
 std::vector<double>
-CoordTransform::polar_to_flux(double e0, double e1, double e2,
+CoordTransform::polar_to_flux(const double & e0, const double & e1, const double & e2,
                               const std::shared_ptr<EquilData> & equilibrium)
 {
-  std::vector<double> outputVector(3);
-  double r;     // local polar r
-  double z;     // local polar z
-  double phi;   // local polar phi
-  double psi;   // local psi
-  double theta; // local theta
-
-  r = e0;
-  z = e1;
-  phi = e2;
-
-  psi = alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
-
-  theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
+  double r = e0;
+  double z = e1;
+  double phi = e2;
+  double psi =
+      alglib::spline2dcalc(equilibrium->psiSpline, r, z); // spline interpolation of psi(R,Z)
+  double theta = atan2(z - equilibrium->zcen, r - equilibrium->rcen);
   if (theta < -M_PI_2)
   {
     theta = 2 * M_PI + theta;
   }
-
-  outputVector[0] = -psi;
-  outputVector[1] = theta;
-  outputVector[2] = phi;
-
+  std::vector<double> outputVector = {-psi, theta, phi};
   return outputVector;
 }
diff --git a/src/aegis_lib/CoordTransform.h b/src/aegis_lib/CoordTransform.h
index c43f8b1..18746a7 100644
--- a/src/aegis_lib/CoordTransform.h
+++ b/src/aegis_lib/CoordTransform.h
@@ -6,22 +6,44 @@
 #include <cmath>
 #include <string>
 #include "EquilData.h"
+#include "alglib/interpolation.h"
+
 
 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(double e0, double e1, double e2);
+std::vector<double> cart_to_polar(const std::vector<double> &xyz);
+
+// polar (R,Z) coords from (X,Y,Z)
+inline std::vector<double> cart_to_polar_xy(const std::vector<double> &xyz)
+{
+  double r = sqrt(pow(xyz[0], 2) + pow(xyz[1], 2));
+  std::vector<double> rz = {r, xyz[2]};
+  return rz;
+};
+
+// get polar angle (PHI) from (X,Y)
+inline double cart_to_polar_phi(const std::vector<double> &xyz)
+{
+  double phi = atan2(-xyz[1], xyz[0]);
+  return phi;
+};
+
+
+std::vector<double> cart_to_polar(const double &e0, const double &e1, const double &e2);
 
-std::vector<double> polar_to_cart(std::vector<double> inputVector);
-std::vector<double> polar_to_cart(double e0, double e1, double e2);
+std::vector<double> polar_to_cart(const std::vector<double> &inputVector);
+std::vector<double> polar_to_cart(const double &e0, const double &e1, const double &e2);
 
 // Polar toroidal to flux coordinates defined by psi spline
-std::vector<double> polar_to_flux(std::vector<double> inputVector,
+std::vector<double> polar_to_flux(const 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);
+
+
+std::vector<double> polar_to_flux(const double &e0, const double &e1, const 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..5291b6b 100644
--- a/src/aegis_lib/EquilData.cpp
+++ b/src/aegis_lib/EquilData.cpp
@@ -746,7 +746,7 @@ EquilData::r_extrema()
         work1[j] = 2;
         break;
       }
-      zpsi = alglib::spline2dcalc(psiSpline, re, ze);
+      zpsi = get_psi(re, ze);
     }
   }
 
@@ -803,7 +803,7 @@ EquilData::rz_splines()
 
 // Caculate B field vector (in toroidal polars) at given position
 std::vector<double>
-EquilData::b_field(std::vector<double> position, std::string startingFrom)
+EquilData::b_field(const std::vector<double> & position, std::string startingFrom)
 {
   std::vector<double> bVector(3); // B in toroidal polars
   double zr;                      // local R from position vector supplied
@@ -824,8 +824,7 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
   // otherwise transform from cartesian -> polar before calculating B
   else
   {
-    std::vector<double> polarPosition;
-    polarPosition = CoordTransform::cart_to_polar(position);
+    auto polarPosition = CoordTransform::cart_to_polar_xy(position);
     zr = polarPosition[0];
     zz = polarPosition[1];
   }
@@ -851,35 +850,14 @@ EquilData::b_field(std::vector<double> position, std::string startingFrom)
 }
 
 std::vector<double>
-EquilData::b_field_cart(std::vector<double> polarBVector, double phi, int normalise)
+EquilData::b_field_cart(const std::vector<double> & polarBVector, const double & phi)
 {
-  std::vector<double> cartBVector(3);
-  double zbx; // cartesian Bx component
-  double zby; // cartesian By component
-  double zbz; // Bz component
-  double zbr; // polar Br component
-  double zbt; // polar toroidal Bt component
-
-  zbr = polarBVector[0];
-  zbz = polarBVector[1];
-  zbt = polarBVector[2];
-
-  zbx = zbr * cos(phi) - zbt * sin(phi);
-  zby = -(zbr * sin(phi) + zbt * cos(phi));
-
-  cartBVector[0] = zbx;
-  cartBVector[1] = zby;
-  cartBVector[2] = zbz;
-
-  if (normalise == 1)
-  {
-    double norm;
-    norm = sqrt(pow(cartBVector[0], 2) + pow(cartBVector[1], 2) + pow(cartBVector[2], 2));
-    cartBVector[0] = cartBVector[0] / norm;
-    cartBVector[1] = cartBVector[1] / norm;
-    cartBVector[2] = cartBVector[2] / norm;
-  }
-
+  double zbr = polarBVector[0];
+  double zbz = polarBVector[1];
+  double zbt = polarBVector[2];
+  double zbx = zbr * cos(phi) - zbt * sin(phi);
+  double zby = -(zbr * sin(phi) + zbt * cos(phi));
+  std::vector<double> cartBVector = {zbx, zby, zbz};
   return cartBVector;
 }
 
@@ -958,7 +936,7 @@ EquilData::write_bfield(int phiSamples)
                                         "Fixup eqdsk required");
           }
           polarB = b_field(polarPos, "polar");               // calculate B(R,Z,phi)
-          cartB = b_field_cart(polarB, polarPos[2], 0);      // transform to B(x,y,z)
+          cartB = b_field_cart(polarB, polarPos[2]);         // transform to B(x,y,z)
           cartPos = CoordTransform::polar_to_cart(polarPos); // transform position to cartesian
 
           // write out magnetic field data for plotting
@@ -1074,7 +1052,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)
     {
@@ -1249,7 +1227,7 @@ EquilData::boundary_rb()
 }
 
 double
-EquilData::omp_power_dep(double psi, double bn, std::string formula)
+EquilData::omp_power_dep(const double & psi, const double & bn, std::string formula)
 {
   double heatFlux;
   double exponential, fpfac;
@@ -1305,13 +1283,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)
@@ -1407,9 +1384,9 @@ EquilData::get_midplane_params()
 
 // check if in b_field from polar coords (R,Z)
 bool
-EquilData::check_if_in_bfield(std::vector<double> xyzPos)
+EquilData::check_if_in_bfield(const 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 +1395,9 @@ EquilData::check_if_in_bfield(std::vector<double> xyzPos)
   }
   return true;
 }
+
+double
+EquilData::get_psi(const double & r, const 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..bf92129 100644
--- a/src/aegis_lib/EquilData.h
+++ b/src/aegis_lib/EquilData.h
@@ -91,6 +91,9 @@ class EquilData : public AegisBase
   // Return eqdsk struct
   eqdskData get_eqdsk_struct();
 
+  // return psi value at given (R,Z) coords
+  double get_psi(const double &r, const double &z);
+
   // Read eqdsk file
   void read_eqdsk(std::string filename);
 
@@ -114,10 +117,10 @@ class EquilData : public AegisBase
 
   // Caculate B field vector (in toroidal polars) at given position
   // set string to "polar" if position vector already in polars
-  std::vector<double> b_field(std::vector<double> position, std::string startingFrom);
+  std::vector<double> b_field(const std::vector<double> &position, std::string startingFrom);
   
   // Convert B Field vectors to cartesian given polar form and value of angle 
-  std::vector<double> b_field_cart(std::vector<double> polarBVector, double phi, int normalise); 
+  std::vector<double> b_field_cart(const std::vector<double> &polarBVector,const double &phi); 
 
   // Write out positions and associated BField vectors in cartesian and/or polar toroidal
   void write_bfield(int phiSamples = 12);
@@ -127,7 +130,7 @@ class EquilData : public AegisBase
   // Determine Rm and Bpm (R and Bpol at omp)
   void boundary_rb();
 
-  double omp_power_dep(double psi, double bn, std::string formula);
+  double omp_power_dep(const double &psi, const double &bn, std::string formula);
 
   void psi_limiter(std::vector<std::vector<double>> vertices);
 
@@ -135,7 +138,7 @@ class EquilData : public AegisBase
 
   void psiref_override();
 
-  bool check_if_in_bfield(std::vector<double> xyzPos);
+  bool check_if_in_bfield(const std::vector<double> &xyzPos);
 
   std::array<double, 3> get_midplane_params(); // return rInnerMidplane, rOuterMidplane, zMidplane
 
diff --git a/src/aegis_lib/Particle.cpp b/src/aegis_lib/Particle.cpp
index d02de88..fa52239 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;
 }
 
@@ -95,8 +91,6 @@ ParticleBase::get_psi(const std::shared_ptr<EquilData> & equilibrium)
 void
 ParticleBase::set_dir(const std::shared_ptr<EquilData> & equilibrium)
 {
-  double norm = 0;              // normalisation constant for magnetic field
-  std::vector<double> normB(3); // normalised magnetic field
   std::vector<double> polarPos; // polar position
   std::vector<double> magn(3);
 
@@ -105,7 +99,7 @@ ParticleBase::set_dir(const std::shared_ptr<EquilData> & equilibrium)
     case coordinateSystem::CARTESIAN:
       polarPos = CoordTransform::cart_to_polar(pos);
       magn = equilibrium->b_field(pos, "cart"); // magnetic field
-      magn = equilibrium->b_field_cart(magn, polarPos[2], 0);
+      magn = equilibrium->b_field_cart(magn, polarPos[2]);
       break;
 
     case coordinateSystem::POLAR:
@@ -116,8 +110,10 @@ ParticleBase::set_dir(const std::shared_ptr<EquilData> & equilibrium)
       break;
   }
 
-  norm = sqrt(pow(magn[0], 2) + pow(magn[1], 2) + pow(magn[2], 2));
+  // normalise particle direction vector
+  double norm = sqrt(pow(magn[0], 2) + pow(magn[1], 2) + pow(magn[2], 2));
   Bfield = magn;
+  std::vector<double> normB(3);
   normB[0] = magn[0] / norm;
   normB[1] = magn[1] / norm;
   normB[2] = magn[2] / norm;
@@ -150,7 +146,7 @@ ParticleBase::align_dir_to_surf(double Bn)
 
 // update position vector of particle with set distance travelled
 void
-ParticleBase::update_vectors(double distanceTravelled)
+ParticleBase::update_vectors(const double & distanceTravelled)
 {
   // This will update the position
   std::vector<double> newPt(3);
@@ -172,7 +168,7 @@ ParticleBase::update_vectors(double distanceTravelled)
 // update position vector of particle with set distance travelled and update
 // the particle direction at new position
 void
-ParticleBase::update_vectors(double distanceTravelled,
+ParticleBase::update_vectors(const double & distanceTravelled,
                              const std::shared_ptr<EquilData> & equilibrium)
 {
   // This will update the position and direction vector of the particle
diff --git a/src/aegis_lib/Particle.h b/src/aegis_lib/Particle.h
index 24b5345..49c405a 100644
--- a/src/aegis_lib/Particle.h
+++ b/src/aegis_lib/Particle.h
@@ -73,8 +73,8 @@ class ParticleBase : public AegisBase
   // return STL vector of unit direction
   std::vector<double> get_dir(); 
   void align_dir_to_surf(double Bn); // align particle dir to surface normal
-  void update_vectors(double distanceTravelled); // update position  
-  void update_vectors(double distanceTravelled, const std::shared_ptr<EquilData> &EquData); // overload to update dir as well
+  void update_vectors(const double &distanceTravelled); // update position  
+  void update_vectors(const double &distanceTravelled, const std::shared_ptr<EquilData> &EquData); // overload to update dir as well
   void check_if_midplane_crossed(const std::array<double, 3> &midplaneParameters); // check if particle has reached inner or outer midplane and set value of atMidplane
   void set_intersection_threshold(double distanceThreshold);
   bool check_if_threshold_crossed();
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..24d9151 100644
--- a/src/aegis_lib/Source.cpp
+++ b/src/aegis_lib/Source.cpp
@@ -101,11 +101,10 @@ 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");
 
-  bField = equilibrium->b_field_cart(bField, polarPos[2], 0);
+  bField = equilibrium->b_field_cart(bField, polarPos[2]);
   double product = 0;
   for (int i = 0; i < 3; i++)
   {
@@ -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
 }