From 9f755e0379e192826a593d758a53124bd28b781a Mon Sep 17 00:00:00 2001
From: James Harrison <jch1g10@soton.ac.uk>
Date: Mon, 27 Mar 2017 16:49:18 +0100
Subject: [PATCH] Add functions momD1 and momD2 to ScalarVP

---
 extras/Hadrons/Modules/MScalar/ScalarVP.cc  | 68 ++++++++++++++++++++-
 extras/Hadrons/Modules/MScalar/ScalarVP.hpp |  5 ++
 2 files changed, 72 insertions(+), 1 deletion(-)

diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.cc b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
index dff636dd..423fb1a2 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.cc
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
@@ -15,7 +15,7 @@ TScalarVP::TScalarVP(const std::string name)
 // dependencies/products ///////////////////////////////////////////////////////
 std::vector<std::string> TScalarVP::getInput(void)
 {
-    prop0Name_ = par().scalarProp + "_0";
+	prop0Name_ = par().scalarProp + "_0";
 	propD1Name_ = par().scalarProp + "_D1";
     propD1D1Name_ = par().scalarProp + "_D1D1";
     propD2Name_ = par().scalarProp + "_D2";
@@ -36,11 +36,77 @@ std::vector<std::string> TScalarVP::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TScalarVP::setup(void)
 {
+	phaseName_.clear();
+	for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        phaseName_.push_back("_shiftphase_" + std::to_string(mu));
+    }
 
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TScalarVP::execute(void)
 {
+	for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
+    }
 
 }
+
+void TScalarVP::momD1(ScalarField &s, EmField &A, FFT &fft)
+{
+    ScalarField buf(env().getGrid()), result(env().getGrid()),
+                Amu(env().getGrid());
+    Complex     ci(0.0,1.0);
+
+    result = zero;
+
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        Amu = peekLorentz(A, mu);
+        buf = (*phase_[mu])*s;
+        fft.FFT_all_dim(buf, buf, FFT::backward);
+        buf = Amu*buf;
+        fft.FFT_all_dim(buf, buf, FFT::forward);
+        result = result - ci*buf;
+    }
+    fft.FFT_all_dim(s, s, FFT::backward);
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        Amu = peekLorentz(A, mu);
+        buf = Amu*s;
+        fft.FFT_all_dim(buf, buf, FFT::forward);
+        result = result + ci*adj(*phase_[mu])*buf;
+    }
+
+    s = result;
+}
+
+void TScalarVP::momD2(ScalarField &s, EmField &Asquared, FFT &fft)
+{
+    ScalarField buf(env().getGrid()), result(env().getGrid()),
+                A2mu(env().getGrid());
+
+    result = zero;
+    
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        A2mu = peekLorentz(Asquared, mu);
+        buf = (*phase_[mu])*s;
+        fft.FFT_all_dim(buf, buf, FFT::backward);
+        buf = A2mu*buf;
+        fft.FFT_all_dim(buf, buf, FFT::forward);
+        result = result + .5*buf;
+    }
+    fft.FFT_all_dim(s, s, FFT::backward);
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        A2mu = peekLorentz(Asquared, mu);        
+        buf = A2mu*s;
+        fft.FFT_all_dim(buf, buf, FFT::forward);
+        result = result + .5*adj(*phase_[mu])*buf;
+    }
+
+    s = result;
+}
diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
index 3c3be434..92a4f246 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
@@ -41,8 +41,13 @@ public:
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    void momD1(ScalarField &s, EmField &A, FFT &fft);
+    void momD2(ScalarField &s, EmField &Asquared, FFT &fft);
 private:
     std::string  prop0Name_, propD1Name_, propD1D1Name_, propD2Name_;
+    std::vector<std::string>   phaseName_;
+    std::vector<ScalarField *> phase_;
 };
 
 MODULE_REGISTER_NS(ScalarVP, TScalarVP, MScalar);