diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp
index c83c478..df88b23 100644
--- a/src/poly/solver/private/polySolver.cpp
+++ b/src/poly/solver/private/polySolver.cpp
@@ -70,7 +70,7 @@ namespace laneEmden {
 }
 
 
-PolySolver::PolySolver(double n, double order) {
+PolySolver::PolySolver(const double n, const double order) {
 
     // --- Check the polytropic index ---
     if (n > 4.99 || n < 0.0) {
@@ -81,14 +81,15 @@ PolySolver::PolySolver(double n, double order) {
     m_polytropicIndex = n;
     m_feOrder = order;
 
-    ResourceManager& rm = ResourceManager::getInstance();
-    const Resource& resource = rm.getResource("mesh:polySphere");
-    const auto &meshIO = std::get<std::unique_ptr<MeshIO>>(resource);
+    // Load and rescale the mesh
+    const ResourceManager& rm = ResourceManager::getInstance();
+    const Resource& genericResource = rm.getResource("mesh:polySphere");
+    const auto &meshIO = std::get<std::unique_ptr<MeshIO>>(genericResource);
     meshIO->LinearRescale(polycoeff::x1(m_polytropicIndex));
     m_mesh = std::make_unique<mfem::Mesh>(meshIO->GetMesh());
 
     // Use feOrder - 1 for the RT space to satisfy Brezzi-Babuska condition
-    // for the H1 and RT spaces
+    // for the H1 and RT [H(div)] spaces
     m_fecH1 = std::make_unique<mfem::H1_FECollection>(m_feOrder, m_mesh->SpaceDimension());
     m_fecRT = std::make_unique<mfem::RT_FECollection>(m_feOrder - 1, m_mesh->SpaceDimension());
 
@@ -120,7 +121,7 @@ void PolySolver::assembleBlockSystem() {
     blockOffsets[0] = 0;
     blockOffsets[1] = feSpaces[0]->GetVSize();
     blockOffsets[2] = feSpaces[1]->GetVSize();
-    blockOffsets.PartialSum();
+    blockOffsets.PartialSum(); // Cumulative sum to get the offsets
 
     // Add integrators to block form one by one
     // We add integrators corresponding to each term in the weak form
@@ -159,9 +160,10 @@ void PolySolver::assembleBlockSystem() {
     Dform->Finalize();
 
     // --- Assemble the NonlinearForm (f) ---
+    // Note that the nonlinear form is built here but its essential true dofs (boundary conditions) are
+    // not set until later (when solve is called). They are applied through the operator rather than directly.
     auto fform = std::make_unique<mfem::NonlinearForm>(m_feTheta.get());
     fform->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex));
-    // TODO: Add essential boundary conditions to the nonlinear form
 
     // -- Build the BlockOperator --
     m_polytropOperator = std::make_unique<PolytropeOperator>(
@@ -189,7 +191,7 @@ void PolySolver::solve() const {
     mfem::Vector zero_rhs(block_offsets.Last());
     zero_rhs = 0.0;
 
-    solverBundle sb = setupNewtonSolver();
+    const solverBundle sb = setupNewtonSolver();
 
     // EMB 2025: Calling Mult gets the gradient of the operator for the NewtonSolver
     // This then is set as the operator for the solver for NewtonSolver
@@ -261,11 +263,11 @@ void PolySolver::setInitialGuess() const {
     // --- Set the initial guess for the solution ---
     mfem::FunctionCoefficient thetaInitGuess (
         [this](const mfem::Vector &x) {
-            double r = x.Norml2();
-            double radius = Probe::getMeshRadius(*m_mesh);
-            double u = 1/radius;
+            const double r = x.Norml2();
+            const double radius = Probe::getMeshRadius(*m_mesh);
+            const double u = 1/radius;
 
-            return -std::pow((u*r), 2)+1.0;
+            return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not
         }
     );
     m_theta->ProjectCoefficient(thetaInitGuess);
diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h
index ca1d02a..d496684 100644
--- a/src/poly/solver/public/polySolver.h
+++ b/src/poly/solver/public/polySolver.h
@@ -46,7 +46,7 @@ struct solverBundle {
 
 class PolySolver {
 public: // Public methods
-    PolySolver(double n, double order);
+    PolySolver(const double n, const double order);
     ~PolySolver();
 
     void solve() const;