Beginning work on adding lagged convolution to Key FHT classes.

Modules/FDEM1D/examples/Hantenna.cpp

@@ -202,7 +202,9 @@ const char *buildString = __DATE__ ", " __TIME__;
 	for (int iy=0; iy<ny; ++iy) {
 	for (int ix=0; ix<nx; ++ix) {
         hreal << receivers->GetLocation(i).transpose() << "\t";
- 		hreal << receivers->GetHfield(0, i).transpose() << "\n";
+ 		//hreal << receivers->GetHfield(0, i).transpose() << "\n"; ( complex, notation )
+ 		hreal << receivers->GetHfield(0, i).transpose().real() << "\t";
+ 		hreal << receivers->GetHfield(0, i).transpose().imag() << "\n";
         ++i;
     }
     }

Modules/FDEM1D/include/FHTKey201.h

@@ -20,6 +20,7 @@
 #define  FHTKEY201_INC
 
 #include "HankelTransform.h"
+#include "CubicSplineInterpolator.h"
 
 namespace Lemma {
 
@@ -27,8 +28,8 @@ namespace Lemma {
       \ingroup FDEM1D
       \brief   Impliments the fast Hankel transform as outlines by Key 2011
       \details This filter uses 51, 101, or 201 filter points and supports both lagged and related
-                kernel arguments. This algorithm is a port of code carrying the following copyright and
-                restriction:
+                kernel arguments. This algorithm is a port of code carrying the following copyright
+                and restriction:
                 %------------------------------------------------------------------%
                 % Copyright (c) 2012 by the Society of Exploration Geophysicists.  %
                 % For more information, go to http://software.seg.org/2012/0003 .  %
@@ -83,10 +84,16 @@ namespace Lemma {
 
         void ComputeRelated(const Real& rho, std::shared_ptr<KernelEM1DManager> KernelManager);
 
+        void ComputeLaggedRelated(const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager);
+
         // ====================  ACCESS        =======================
 
         // ====================  INQUIRY       =======================
 
+        Real GetABSER();
+
+        void SetLaggedArg(const Real& rho);
+
         /** Returns the name of the underlying class, similiar to Python's type */
         virtual std::string GetName() const;
 
@@ -101,6 +108,12 @@ namespace Lemma {
         // Shared Filter Weights
 		static const Eigen::Matrix<Real, 201, 3>  WT201;
 
+        /// Spines for lagged convolutions (real part)
+        std::vector <std::shared_ptr<CubicSplineInterpolator> > splineVecReal;
+
+        /// Spines for lagged convolutions (imaginary part)
+        std::vector < std::shared_ptr<CubicSplineInterpolator> > splineVecImag;
+
         /// Holds answer, dimensions are NumConv, and NumberRelated.
         Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> Zans;
 

Modules/FDEM1D/src/EMEarth1D.cpp

@@ -231,7 +231,7 @@ namespace Lemma {
         if (Antenna->GetName() == std::string("PolygonalWireAntenna") || Antenna->GetName() == std::string("TEMTransmitter") ) {
             icalc += 1;
             // Check to see if they are all on a plane? If so we can do this fast
-            if (Antenna->IsHorizontallyPlanar() && HankelType == ANDERSON801) {
+            if (Antenna->IsHorizontallyPlanar() && 1==2 && ( HankelType == ANDERSON801 || HankelType== FHTKEY201  )) {
                 #ifdef HAVE_BOOST_PROGRESS
                 if (progressbar) {
                     disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
@@ -243,7 +243,12 @@ namespace Lemma {
                     #pragma omp parallel
                     {
                     #endif
-                    auto Hankel = FHTAnderson801::NewSP();
+                    //if (HankelType == ANDERSON801) {
+                        auto Hankel = FHTAnderson801::NewSP();
+                    //}
+                    //else if(HankelType == FHTKEY201) {
+                    //    auto Hankel = FHTKey201::NewSP();
+                    //}
                     #ifdef LEMMAUSEOMP
                     #pragma omp for schedule(static, 1)
                     #endif
@@ -797,9 +802,7 @@ namespace Lemma {
         //std::cout << Hankel->GetAnswer() << std::endl;
         //std::cout << Hankel->GetArg() << std::endl;
 
-
         // Sort the dipoles by rho
-
         for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
         //for (int idip=0; idip<1; ++idip) {
             auto tDipole = antenna->GetDipoleSource(idip);

Modules/FDEM1D/src/FHTKey201.cpp

@@ -371,4 +371,65 @@ namespace Lemma {
         return ;
     }		// -----  end of method FHTKey201::ComputeRelated  -----
 
+    //--------------------------------------------------------------------------------------
+    //       Class:  FHTKey201
+    //      Method:  ComputeLaggedRelated
+    //--------------------------------------------------------------------------------------
+    void FHTKey201::ComputeLaggedRelated ( const Real& rho, const int& nlag, std::shared_ptr<KernelEM1DManager> KernelManager ) {
+
+        int nrel = (int)(KernelManager->GetSTLVector().size());
+        Eigen::Matrix<Complex, 201, Eigen::Dynamic > Zwork;
+        Zans= Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(nlag, nrel);
+        Zwork.resize(201, nrel);
+        VectorXr lambda = WT201.col(0).array()/rho;
+        int NumFun = 0;
+        int idx = 0;
+
+        // Get Kernel values
+        for (int ir=0; ir<lambda.size(); ++ir) {
+            // irelated loop
+            ++NumFun;
+            KernelManager->ComputeReflectionCoeffs(lambda(ir), idx, rho);
+            for (int ir2=0; ir2<nrel; ++ir2) {
+                // Zwork* needed due to sign convention of filter weights
+ 			    Zwork(ir, ir2) = std::conj(KernelManager->GetSTLVector()[ir2]->RelBesselArg(lambda(ir)));
+            }
+
+        }
+
+        // We diverge slightly from Key here, each kernel is evaluated seperately, whereby instead
+        // they prefer to sum them. The reason is that all those terms have been removed from the kernels
+        // in the interests of making them as generic and reusable as possible. This approach requires slightly
+        // more multiplies, but the same number of kernel evaluations, which is the expensive part.
+        // Inner product and scale
+        for (int ir2=0; ir2<nrel; ++ir2) {
+            Zans(0, ir2) = Zwork.col(ir2).dot(WT201.col(KernelManager->GetSTLVector()[ir2]->GetBesselOrder() + 1))/rho;
+        }
+
+        return ;
+    }		// -----  end of method FHTKey201::ComputeLaggedRelated  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  FHTKey201
+    //      Method:  GetABSER
+    //--------------------------------------------------------------------------------------
+    Real FHTKey201::GetABSER (  ) {
+        return WT201(1,0)/WT201(0,0);
+    }		// -----  end of method FHTKey201::GetABSER  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  FHTKey201
+    //      Method:  SetLaggedArg
+    //--------------------------------------------------------------------------------------
+    void FHTKey201::SetLaggedArg ( const Real& rho  ) {
+        for (int i=0; i<Zans.cols(); ++ i) {
+            Zans(0, i) = Complex( splineVecReal[i]->Interpolate(rho),
+                                  splineVecImag[i]->Interpolate(rho) );
+        }
+        return ;
+    }		// -----  end of method FHTKey201::SetLaggedArg  -----
+
+
 }		// -----  end of Lemma  name  -----