Opalstack -> GitHub sync test

README.md

@@ -1,4 +1,4 @@
 ==About
-Merlin is a NMR modelling module for Lemma. 
+Merlin is a surface NMR modelling module for Lemma. 
 
 =Capabilities 

include/KernelV0.h

@@ -216,6 +216,20 @@ namespace Lemma {
         }
 
         /**
+         * @param[in] dfin is the transmitter offset from Larmor in Hz. Calculated as Txf - Lf
+         */
+        inline void SetDf( const Real& dfin ) {
+            df = dfin;
+        }
+
+        /**
+         * @return the pulse offset from the Larmor frequency
+         */
+        inline Real GetDf ( ) {
+            return df;
+        }
+
+        /**
          * @param[in] value the 1D-EM model used for calculations
          */
         inline void SetLayeredEarthEM ( std::shared_ptr< LayeredEarthEM > value ) {
@@ -376,6 +390,7 @@ namespace Lemma {
         Real                                      Temperature=283.;
         Real                                      Taup = .020;  // Sec
         Real                                      Larmor;
+        Real                                      df=0;
 
         Vector3r                                  Size;
         Vector3r                                  Origin;

python/pyMerlin.cpp

@@ -56,6 +56,7 @@ PYBIND11_MODULE(Merlin, m) {
             .def("SetPulseDuration", &Lemma::KernelV0::SetPulseDuration, "Sets the duration of the pulse")
             .def("SetDepthLayerInterfaces", &Lemma::KernelV0::SetDepthLayerInterfaces, "Sets the layer depth interfaces")
             .def("SetHankelTransformType", &Lemma::KernelV0::SetHankelTransformType, "Sets the Hankel transform type")
+            .def("SetDf", &Lemma::KernelV0::SetDf, "Sets the transmitter offset frequncy in Hz, (Txf - Lf)")
 
             // accessors
             .def("GetName", &Lemma::KernelV0::GetName, "Returns the name of the class")
@@ -65,6 +66,7 @@ PYBIND11_MODULE(Merlin, m) {
             .def("GetInterfaces", &Lemma::KernelV0::GetInterfaces, "Returns the layer interfaces")
             .def("GetPulseCurrent", &Lemma::KernelV0::GetPulseCurrent, "Returns the pulse current")
             .def("GetPulseDuration", &Lemma::KernelV0::GetPulseDuration, "Returns the length of the pulse moment")
+            .def("GetDf", &Lemma::KernelV0::GetDf, "Returns the transmitter offset in Hz")
 
             // operations
             .def("CalculateK0", &Lemma::KernelV0::CalculateK0, "Calculates an intial amplitude kernel")

src/KernelV0.cpp

@@ -218,6 +218,10 @@ namespace Lemma {
                 EMEarths[tx]->SetHankelTransformMethod(HankelType);
                 EMEarths[tx]->SetTxRxMode(TX);
                 TxRx[tx]->SetCurrent(1.);
+
+                // calculate df TODO, fix for multiple tx frequencies
+                df = TxRx[tx]->GetFrequency(0) - Larmor/(2.*PI);
+                std::cout << "df=" << df << std::endl;
         }
         for (auto rx : Rx) {
             if (EMEarths.count(rx)) {
@@ -482,9 +486,35 @@ namespace Lemma {
         // Calcuate vector of all responses
         VectorXcr F = VectorXcr::Zero( PulseI.size() );
         for (int iq=0; iq<PulseI.size(); ++iq) {
-            // Compute the tipping angle
+
+            /////////////////////////////////////////////////////////////////////////////
+            // Compute the tipping angle for on-resonance
+            // Weichman formulation
+            /*
             Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta));
             F(iq) = -volume*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
+            */
+
+            /////////////////////////////////////////////////////////////////////////////
+            // compute tipping ange for off-resonance, from MRSMatlab (Mueller, et. al)
+            // but based on seperate works by Grombacher and Walbrecher.
+            //theta    = atan2(0.5*gamma*pm_vec(n)/taup*(Bcomps.alpha - Bcomps.beta),(2*pi*df));
+            //flip_eff = sqrt((0.5*gamma*pm_vec(n)*(Bcomps.alpha - Bcomps.beta)).^2 + ...
+            //            (2*pi*df*taup).^2 );
+            //m     = sin(flip_eff) .* sin(theta) + ...
+            //    1i*(-1)*sin(theta).*cos(theta) .* (cos(flip_eff) - 1);
+            //kern = gamma * earth.erdt^2 * 3.29e-3 * Px .* Bcomps.e_zeta.^2 .* ...
+            //            (Bcomps.alpha + Bcomps.beta) .* m;
+            //K(n,:) = sum(sum(kern.*dh*dz));
+            // TODO, benchmark calls to pow below...modern compilers should optimize this
+            // Real df = 25;  // Hz? df is a class data member
+            Real theta = std::atan2( 0.5*GAMMA*PulseI(iq)*(EBT.alpha-EBT.beta), 2*PI*df );
+            Real flip_eff = std::sqrt( std::pow(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta),2) +
+                                       std::pow(2*PI*df*Taup,2));
+            Complex m = std::sin(flip_eff)*std::sin(theta) +
+                Complex(0,-1)*std::sin(theta)*std::cos(theta) * (std::cos(flip_eff) - 1.);
+
+            F(iq) = -volume*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*m*PhaseTerm;
         }
         return F;
     }