now using const passing of object reference instead of pointers

1 files changed, 34 insertions, 31 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index f948586..f3271b0 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -12,45 +12,45 @@
 #include "graph.h"
 
 template <class R_t, class X_t>
-void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
+void flip_cluster(state_t<R_t, X_t>& s, v_t v0, const R_t& r, gsl_rng *rand) {
   v_t nv = 0;
 
   std::stack<v_t> stack;
   stack.push(v0);
 
-  std::vector<bool> marks(state->g.nv, false);
+  std::vector<bool> marks(s.g.nv, false);
 
   while (!stack.empty()) {
     v_t v = stack.top();
     stack.pop();
 
-    if (!marks[v]) {
+    if (!marks[v]) { // don't reprocess anyone we've already visited!
       X_t si_old, si_new;
       R_t R_old, R_new;
 
-      R_old = state->R;
+      R_old = s.R;
       marks[v] = true;
 
-      bool v_is_ghost = (v == state->nv);
+      bool v_is_ghost = (v == s.nv); // ghost site has the last index
 
       if (v_is_ghost) {
-        R_new = r.act(R_old);
+        R_new = r.act(R_old); // if we are, then we're moving the transformation
       } else {
-        si_old = state->spins[v];
-        si_new = r.act(si_old);
+        si_old = s.spins[v];
+        si_new = r.act(si_old); // otherwise, we're moving the spin at our site
       }
 
-      for (const v_t &vn : state->g.v_adj[v]) {
+      for (const v_t &vn : s.g.v_adj[v]) {
         X_t sj;
-        bool vn_is_ghost = (vn == state->nv);
+        bool vn_is_ghost = (vn == s.nv); // any of our neighbors could be the ghost
 
         if (!vn_is_ghost) {
-          sj = state->spins[vn];
+          sj = s.spins[vn];
         }
 
         double prob;
 
-        if (v_is_ghost || vn_is_ghost) {
+        if (v_is_ghost || vn_is_ghost) { // if this is a ghost-involved bond...
           X_t rs_old, rs_new;
           v_t non_ghost;
           if (vn_is_ghost) {
@@ -63,51 +63,54 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
             non_ghost = vn;
           }
 
-          double dE = state->H(rs_old) - state->H(rs_new);
+          double dE = s.H(rs_old) - s.H(rs_new);
+
 #ifdef FINITE_STATES
           prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
 #else
-          prob = 1.0 - exp(-dE / state->T);
+          prob = 1.0 - exp(-dE / s.T);
 #endif
 
-          state->M -= rs_old;
-          state->M += rs_new;
+          s.M -= rs_old;
+          s.M += rs_new;
 
-          state->E += dE;
+          s.E += dE;
 
-          for (D_t i = 0; i < state->D; i++) {
-            L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L;
+          for (D_t i = 0; i < s.D; i++) {
+            L_t x = (non_ghost / (v_t)pow(s.L, s.D - i - 1)) % s.L;
 
-            state->ReF[i] -= rs_old * state->precomputed_cos[x];
-            state->ReF[i] += rs_new * state->precomputed_cos[x];
+            s.ReF[i] -= rs_old * s.precomputed_cos[x];
+            s.ReF[i] += rs_new * s.precomputed_cos[x];
 
-            state->ImF[i] -= rs_old * state->precomputed_sin[x];
-            state->ImF[i] += rs_new * state->precomputed_sin[x];
+            s.ImF[i] -= rs_old * s.precomputed_sin[x];
+            s.ImF[i] += rs_new * s.precomputed_sin[x];
           }
-        } else {
-          double dE = state->J(si_old, sj) - state->J(si_new, sj);
+        } else { // otherwise, we're at a perfectly normal bond!
+          double dE = s.J(si_old, sj) - s.J(si_new, sj);
+
 #ifdef FINITE_STATES
           prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
 #else
-          prob = 1.0 - exp(-dE / state->T);
+          prob = 1.0 - exp(-dE / s.T);
 #endif
-          state->E += dE;
+
+          s.E += dE;
         }
 
-        if (gsl_rng_uniform(rand) < prob) { // and with probability...
+        if (gsl_rng_uniform(rand) < prob) {
           stack.push(vn); // push the neighboring vertex to the stack
         }
       }
 
       if (v_is_ghost) {
-        state->R = R_new;
+        s.R = R_new;
       } else {
-        state->spins[v] = si_new;
+        s.spins[v] = si_new;
         nv++;
       }
     }
   }
 
-  state->last_cluster_size = nv;
+  s.last_cluster_size = nv;
 }