diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h
index 6cf34e0c..d4790df2 100644
--- a/Grid/qcd/utils/WilsonLoops.h
+++ b/Grid/qcd/utils/WilsonLoops.h
@@ -6,10 +6,12 @@
 
     Copyright (C) 2015
 
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
-Author: paboyle <paboyle@ph.ed.ac.uk>
+    Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+    Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+    Author: neo <cossu@post.kek.jp>
+    Author: paboyle <paboyle@ph.ed.ac.uk>
+    Author: James Harrison <J.Harrison@soton.ac.uk>
+    Author: Antonin Portelli <antonin.portelli@me.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -645,6 +647,184 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
       }
     }
   }
+
+  //////////////////////////////////////////////////
+  // Wilson loop of size (R1, R2), oriented in mu,nu plane
+  //////////////////////////////////////////////////
+  static void wilsonLoop(GaugeMat &wl, const std::vector<GaugeMat> &U,
+                           const int Rmu, const int Rnu,
+                           const int mu, const int nu) {
+    wl = U[nu];
+
+    for(int i = 0; i < Rnu-1; i++){
+      wl = Gimpl::CovShiftForward(U[nu], nu, wl);
+    }
+
+    for(int i = 0; i < Rmu; i++){
+      wl = Gimpl::CovShiftForward(U[mu], mu, wl);
+    }
+
+    for(int i = 0; i < Rnu; i++){
+      wl = Gimpl::CovShiftBackward(U[nu], nu, wl);
+    }
+
+    for(int i = 0; i < Rmu; i++){
+      wl = Gimpl::CovShiftBackward(U[mu], mu, wl);
+    }
+  }
+  //////////////////////////////////////////////////
+  // trace of Wilson Loop oriented in mu,nu plane
+  //////////////////////////////////////////////////
+  static void traceWilsonLoop(LatticeComplex &wl,
+                                const std::vector<GaugeMat> &U,
+                                const int Rmu, const int Rnu,
+                                const int mu, const int nu) {
+    GaugeMat sp(U[0]._grid);
+    wilsonLoop(sp, U, Rmu, Rnu, mu, nu);
+    wl = trace(sp);
+  }
+  //////////////////////////////////////////////////
+  // sum over all planes of Wilson loop
+  //////////////////////////////////////////////////
+  static void siteWilsonLoop(LatticeComplex &Wl,
+                            const std::vector<GaugeMat> &U,
+                            const int R1, const int R2) {
+    LatticeComplex siteWl(U[0]._grid);
+    Wl = zero;
+    for (int mu = 1; mu < U[0]._grid->_ndimension; mu++) {
+      for (int nu = 0; nu < mu; nu++) {
+        traceWilsonLoop(siteWl, U, R1, R2, mu, nu);
+        Wl = Wl + siteWl;
+        traceWilsonLoop(siteWl, U, R2, R1, mu, nu);
+        Wl = Wl + siteWl;
+      }
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum over planes of Wilson loop with length R1
+  // in the time direction
+  //////////////////////////////////////////////////
+  static void siteTimelikeWilsonLoop(LatticeComplex &Wl,
+                            const std::vector<GaugeMat> &U,
+                            const int R1, const int R2) {
+    LatticeComplex siteWl(U[0]._grid);
+
+    int ndim = U[0]._grid->_ndimension;
+
+    Wl = zero;
+    for (int nu = 0; nu < ndim - 1; nu++) {
+      traceWilsonLoop(siteWl, U, R1, R2, ndim-1, nu);
+      Wl = Wl + siteWl;
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum Wilson loop over all planes orthogonal to the time direction
+  //////////////////////////////////////////////////
+  static void siteSpatialWilsonLoop(LatticeComplex &Wl,
+                            const std::vector<GaugeMat> &U,
+                            const int R1, const int R2) {
+    LatticeComplex siteWl(U[0]._grid);
+
+    Wl = zero;
+    for (int mu = 1; mu < U[0]._grid->_ndimension - 1; mu++) {
+      for (int nu = 0; nu < mu; nu++) {
+        traceWilsonLoop(siteWl, U, R1, R2, mu, nu);
+        Wl = Wl + siteWl;
+        traceWilsonLoop(siteWl, U, R2, R1, mu, nu);
+        Wl = Wl + siteWl;
+      }
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of Wilson loop
+  //////////////////////////////////////////////////
+  static Real sumWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    std::vector<GaugeMat> U(4, Umu._grid);
+
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex Wl(Umu._grid);
+
+    siteWilsonLoop(Wl, U, R1, R2);
+
+    TComplex Tp = sum(Wl);
+    Complex p = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of timelike Wilson loop
+  //////////////////////////////////////////////////
+  static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    std::vector<GaugeMat> U(4, Umu._grid);
+
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex Wl(Umu._grid);
+
+    siteTimelikeWilsonLoop(Wl, U, R1, R2);
+
+    TComplex Tp = sum(Wl);
+    Complex p = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of spatial Wilson loop
+  //////////////////////////////////////////////////
+  static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    std::vector<GaugeMat> U(4, Umu._grid);
+
+    for (int mu = 0; mu < Umu._grid->_ndimension; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex Wl(Umu._grid);
+
+    siteSpatialWilsonLoop(Wl, U, R1, R2);
+
+    TComplex Tp = sum(Wl);
+    Complex p = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of Wilson loop
+  //////////////////////////////////////////////////
+  static Real avgWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
+    Real sumWl = sumWilsonLoop(Umu, R1, R2);
+    Real vol = Umu._grid->gSites();
+    Real faces = 1.0 * ndim * (ndim - 1);
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of timelike Wilson loop
+  //////////////////////////////////////////////////
+  static Real avgTimelikeWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
+    Real sumWl = sumTimelikeWilsonLoop(Umu, R1, R2);
+    Real vol = Umu._grid->gSites();
+    Real faces = 1.0 * (ndim - 1);
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of spatial Wilson loop
+  //////////////////////////////////////////////////
+  static Real avgSpatialWilsonLoop(const GaugeLorentz &Umu,
+                            const int R1, const int R2) {
+    int ndim = Umu._grid->_ndimension;
+    Real sumWl = sumSpatialWilsonLoop(Umu, R1, R2);
+    Real vol = Umu._grid->gSites();
+    Real faces = 1.0 * (ndim - 1) * (ndim - 2);
+    return sumWl / vol / faces / Nc; // Nc dependent... FIXME
+  }
 };
 
 typedef WilsonLoops<PeriodicGimplR> ColourWilsonLoops;