From c902f022795157b28ef01cffcc499753e2212c28 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 22 Jun 2017 23:36:35 -0400
Subject: added correct support for negative external fields

---
 lib/wolff.h       | 2 ++
 lib/wolff_tools.c | 8 ++++++--
 src/wolff.c       | 4 ++--
 3 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/lib/wolff.h b/lib/wolff.h
index 4270f33..428ce85 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -31,6 +31,8 @@ typedef struct {
   int32_t dHb;
 } cluster_t;
 
+int32_t sign(double x);
+
 cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
                         gsl_rng *r);
 
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index b3758cc..7886bab 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -2,6 +2,10 @@
 #include "queue.h"
 #include "wolff.h"
 
+int32_t sign(double x) {
+  return x > 0 ? 1 : -1;
+}
+
 graph_t *graph_add_ext(const graph_t *g) {
   graph_t *h = (graph_t *)calloc(1, sizeof(graph_t));
 
@@ -131,8 +135,8 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
 
   cluster_t *c = flip_cluster(s->g, ps, s->spins, r);
 
-  s->M += -2 * c->dHb;
-  s->H += 2 * (c->dJb + H * c->dHb);
+  s->M += - sign(H) * 2 * c->dHb;
+  s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
 
   uint32_t n_flips = c->nv;
 
diff --git a/src/wolff.c b/src/wolff.c
index a8266ff..1176fde 100644
--- a/src/wolff.c
+++ b/src/wolff.c
@@ -84,8 +84,8 @@ int main(int argc, char *argv[]) {
 
   s->g = graph_add_ext(h);
   s->spins = (bool *)calloc(h->nv + 1, sizeof(bool));
-  s->M = h->nv;
-  s->H = -(1.0 * h->ne) - H * h->nv;
+  s->M = sign(H) * h->nv;
+  s->H = -(1.0 * h->ne) - sign (H) * H * h->nv;
 
   double *bond_probs = (double *)malloc(2 * sizeof(double));
   bond_probs[0] = 1 - exp(-2 / T);
-- 
cgit v1.2.3-70-g09d2