working on file output a bit

examples/KernelV0.cpp

@@ -22,7 +22,10 @@ using namespace Lemma;
 
 std::shared_ptr<PolygonalWireAntenna> CircularLoop ( int nd, Real radius, Real Offsetx, Real Offsety ) ;
 
-int main() {
+int main(int argc, char** argv) {
+
+    Real offset = atof(argv[1]);
+        std::cout << offset << std::endl;
 
 	auto earth = LayeredEarthEM::NewSP();
 		earth->SetNumberOfLayers(3);
@@ -34,7 +37,7 @@ int main() {
 
     // Transmitter loops
     auto Tx1 = CircularLoop(21, 15, 100, 100);
-    auto Tx2 = CircularLoop(21, 15, 100, 124.8);
+    auto Tx2 = CircularLoop(21, 15, 100, 100 + offset); // 100, 115, 124.8, 130
     //auto Tx1 = CircularLoop(60, 15, 0, 0); // was 60
 
     auto Kern = KernelV0::NewSP();
@@ -44,7 +47,7 @@ int main() {
 
         Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
         Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
-        Kern->SetTolerance( 1e-9 );
+        Kern->SetTolerance( 1e-12 );
 
         Kern->SetPulseDuration(0.020);
         VectorXr I(36);
@@ -66,24 +69,25 @@ int main() {
             interfaces(ilay) = interfaces(ilay-1) + thick;
             thick *= 1.1;
         }
-        // TODO log spacing results in a strange transposition of the matrix? Difficult to understand
-        //10**np.linspace(np.log10(10),np.log10(19),10)
-//         VectorXr interfaces2 = VectorXr::LinSpaced(31, std::log10(2), std::log10(150)); // 30 log spaced
-//         for (int i=0; i<interfaces2.size(); ++i) {
-//             interfaces(i) = std::pow(10, interfaces2(i));
-//         }
-//         std::cout << interfaces << std::endl;
-
         Kern->SetDepthLayerInterfaces( interfaces ); // nlay, low, high
 
     // We could, I suppose, take the earth model in here? For non-linear that
     // may be more natural to work with?
-    //std::vector<std::string> tx = {std::string("Coil 1"), std::string("Coil 2") };
-    std::vector<std::string> tx = {std::string("Coil 1")};
+    std::vector<std::string> tx = {std::string("Coil 1"), std::string("Coil 2") };
+    //std::vector<std::string> tx = {std::string("Coil 1")};
     std::vector<std::string> rx = {std::string("Coil 1")};
-    Kern->CalculateK0( tx, rx, true );
-
-    std::ofstream out = std::ofstream("k.yaml");
+    Kern->CalculateK0( tx, rx, false );
+
+    std::ofstream dout = std::ofstream(std::string("k-")+ std::string(argv[1])+ std::string(".dat"));
+        dout << interfaces.transpose() << std::endl;
+        dout << I.transpose() << std::endl;
+        dout << "#real\n";
+        dout << Kern->GetKernel().real() << std::endl;
+        dout << "#imag\n";
+        dout << Kern->GetKernel().imag() << std::endl;
+        dout.close();
+
+    std::ofstream out = std::ofstream(std::string("k-")+std::string(argv[1])+std::string(".yaml"));
     out << *Kern;
     out.close();
 }

examples/plotkernel.py

@@ -2,6 +2,11 @@ import numpy as np
 import matplotlib.pyplot as plt
 import sys
 from pylab import meshgrid
+from matplotlib.colors import LightSource
+from matplotlib.ticker import ScalarFormatter
+from matplotlib.ticker import MaxNLocator
+from matplotlib.ticker import AutoMinorLocator
+from matplotlib.ticker import LogLocator
 
 kf = open(sys.argv[1])
 ifaces = np.array( kf.readline().split(), dtype=np.float ) 
@@ -9,41 +14,39 @@ q = np.array( kf.readline().split(), dtype=np.float )
 q = np.append( q, (q[-1]+q[-2]) ) # for pcolor mesh
 Y,X = meshgrid( ifaces, q )   
 
-K = np.loadtxt(sys.argv[1], comments="#", skiprows=2)
+K = 1e9*np.loadtxt(sys.argv[1], comments="#", skiprows=2)
 nx, ny = np.shape(K)
 
 fig = plt.figure( )
-ax1 = fig.add_axes( [.10,.10,.35,.80] )
-ax2 = fig.add_axes( [.5,.10,.35,.80], sharex=ax1, sharey=ax1 )
-axcb = fig.add_axes( [.9,.10,.05,.80] )
+ax1 =  fig.add_axes( [.100,.125,.335,.775] )
+ax2 =  fig.add_axes( [.465,.125,.335,.775] , sharex=ax1, sharey=ax1 )
+axcb = fig.add_axes( [.835,.125,.040,.775] )
+
+ccmap = "seismic" # RdBu
 
 # Real plot
-ax1.pcolormesh(X, Y, K[0:nx//2].T , cmap="RdBu", vmin=-np.max(np.abs(K)), vmax=np.max(np.abs(K)) )
+ax1.pcolormesh(X, Y, K[0:nx//2].T, cmap=ccmap, vmin=-np.max(np.abs(K)), vmax=np.max(np.abs(K)))
 ax1.set_ylim( ifaces[-1], ifaces[0] )
 ax1.set_xlim( q[-1], q[0] )
 ax1.set_xscale("log", nonposx='clip')
-#plt.colorbar()
+ax1.minorticks_off()
+ax1.set_title(r"$\mathrm{Re} \left( \mathcal{K}_N^{1D} \right)$")
+
+ax1.xaxis.set_major_locator(LogLocator(base = 2.0))
+ax1.get_xaxis().set_major_formatter(ScalarFormatter())
 
 # imaginary 
-im = ax2.pcolormesh(X, Y, K[nx//2:].T , cmap="RdBu", vmin=-np.max(np.abs(K)), vmax=np.max(np.abs(K)) )
+im = ax2.pcolormesh(X, Y, K[nx//2:].T, cmap=ccmap, vmin=-np.max(np.abs(K)), vmax=np.max(np.abs(K)))
 plt.setp(ax2.get_yticklabels(), visible=False)
-plt.colorbar(im, axcb)
+cb = plt.colorbar(im, axcb)
+cb.set_label(r"$\overline{V}^{\left(0\right)}_N$  (nV)")
 
 ax1.set_ylabel("Depth (m)")
-ax1.set_xlabel("Pulse moment (A$\cdot$s)")
-ax2.set_xlabel("Pulse moment (A$\cdot$s)")
-
-plt.show()
-exit()
-
+ax1.set_xlabel(r"Pulse moment (A$\cdot$s)")
+ax2.set_xlabel(r"Pulse moment (A$\cdot$s)")
+ax2.set_title(r"$\mathrm{Im} \left( \mathcal{K}_N^{1D} \right)$")
 
-plt.matshow(K[0:nx//2])
-plt.colorbar()
+plt.suptitle("exp 0")
 
-plt.matshow(K[nx//2:])
-plt.colorbar()
-
-KA = np.abs(K[0:nx//2] + 1j*K[nx//2:])
-plt.matshow(1e9*KA, cmap='viridis')
-plt.colorbar()
 plt.show()
+exit()

include/KernelV0.h

@@ -127,6 +127,13 @@ namespace Lemma {
         }		// -----  end of method KernelV0::get_SigmaModel  -----
 
         /**
+         * @return the kernel matrix
+         */
+        inline MatrixXcr GetKernel ( ) {
+            return Kern;
+        }
+
+        /**
          * @param[in] value the 1D-EM model used for calculations
          */
         inline void SetLayeredEarthEM ( std::shared_ptr< LayeredEarthEM > value ) {

src/KernelV0.cpp

@@ -154,21 +154,12 @@ namespace Lemma {
         std::cout << "Calculating K0 kernel\n";
         Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
         for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
-            std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "<< Interfaces(ilay+1) << std::endl; //<< " q " << iq << std::endl;
+            std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "<< Interfaces(ilay+1) << std::endl;
             Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
             Origin(2) = Interfaces(ilay);
             IntegrateOnOctreeGrid( vtkOutput );
         }
-        std::cout << "\rFinished KERNEL\n";
-
-        std::ofstream out = std::ofstream("k.dat");
-        out << Interfaces.transpose() << std::endl;
-        out << PulseI.transpose() << std::endl;
-        out << "#real\n";
-        out << Kern.real() << std::endl;
-        out << "#imag\n";
-        out << Kern.imag() << std::endl;
-        out.close();
+        std::cout << "\nFinished KERNEL\n";
     }
 
     //--------------------------------------------------------------------------------------