Working on kernel matrix calculation

examples/KernelV0.cpp

 		earth->SetNumberOfLayers(3);
 		earth->SetLayerConductivity( (VectorXcr(3) << Complex(0.,0), Complex(1./50.,0), Complex(1./100.)).finished() );
 		earth->SetLayerThickness( (VectorXr(1) << 10).finished() );
+        // Set mag field info
+        // From NOAA, Laramie WY, June 9 2016, aligned with mag. north
+        earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
 
     // Transmitter loops
-    auto Tx1 = CircularLoop(60, 15, 100, 100);
-    auto Tx2 = CircularLoop(60, 15, 100, 120);
+    auto Tx1 = CircularLoop(65, 15, 100, 100);
+    auto Tx2 = CircularLoop(65, 15, 100, 120);
     //auto Tx1 = CircularLoop(60, 15, 0, 0); // was 60
 
     auto Kern = KernelV0::NewSP();
         Kern->SetLayeredEarthEM( earth );
         // std::cout << *Kern << std::endl;
 
+        // Kern->SetPulseDuration();
+        // Kern->SetPulseCurrent();
+        // Kern->SetPulseMoments();
+        // Kern->SetDepthLayers();
+        // Kern->SetIntegrationOrigin();
+        // Kern->SetIntegrationSize();
+
+
     // 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::vector<std::string> rx = {std::string("Coil 2")};
+    std::vector<std::string> rx = {std::string("Coil 1")};
     Kern->CalculateK0( tx, rx , true ); //, false );
     //Kern->CalculateK0( "Coil 1", "Coil 1" );
 
 }
 
 std::shared_ptr<Lemma::PolygonalWireAntenna> CircularLoop ( int nd, Real Radius, Real Offsetx, Real Offsety ) {
+
     auto Tx1 = Lemma::PolygonalWireAntenna::NewSP();
          Tx1->SetNumberOfPoints(nd);
 
     VectorXr range = VectorXr::LinSpaced(nd, 0, 2*PI);
     int ii;
-    for (ii=0; ii<nd-1; ++ii) {
+    for (ii=0; ii<nd; ++ii) {
         Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*std::cos(range(ii)), Offsety+Radius*std::sin(range(ii)),  -1e-3));
     }
-    Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety,  -1e-3));
+    //Tx1->SetPoint(ii, Vector3r(Offsetx+Radius*1, Offsety,  -1e-3));
 
     Tx1->SetCurrent(1.);
     Tx1->SetNumberOfTurns(1);

include/KernelV0.h

         void CalculateK0 (const std::vector< std::string >& tx, const std::vector< std::string >& rx,
                 bool vtkOutput=false );
 
+        /**
+         *  Sets the temperature, which has implications in calculation of \f$ M_N^{(0)}\f$. Units in
+         *  Kelvin.
+         */
+        void SetTemperature(const Real& tempK) {
+            Temperature = tempK;
+        }
+
         // ====================  INQUIRY       =======================
         /**
          *  Returns the name of the underlying class, similiar to Python's type
 
         Real                                      VOLSUM;
         Real                                      tol=1e-3;
-        //Real                                        Temperature;
+        Real                                      Temperature=283.;
 
         Complex                                   SUM;
 
 
         #ifdef LEMMAUSEVTK
         std::map< int, Complex  >                 LeafDict;
+        std::map< int, Real     >                 LeafDictErr;
         #endif
 
         // Physical constants and conversion factors

src/KernelV0.cpp

                     EMEarths[rx]->SetTxRxMode(RX);
             }
         }
-        IntegrateOnOctreeGrid( 1e-7, vtkOutput );
-
+        IntegrateOnOctreeGrid( 1e-13, vtkOutput );
     }
 
     //--------------------------------------------------------------------------------------
 
         this->tol = tolerance;
         //Vector3r                Size;
-            Size << 200,200,20;
+            Size << 200,200,2;
         //Vector3r                Origin;
             Origin << 0,0,5.0;
         Vector3r                cpos;  // centre position
                 ki->SetNumberOfComponents(1);
                 ki->SetName("Im($K_0$)");
                 ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkDoubleArray* km = vtkDoubleArray::New();
+                km->SetNumberOfComponents(1);
+                km->SetName("mod($K_0$)");
+                km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            vtkIntArray* kid = vtkIntArray::New();
+                kid->SetNumberOfComponents(1);
+                kid->SetName("ID");
+                kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
+            //vtkDoubleArray* kerr = vtkDoubleArray::New();
+            //    kerr->SetNumberOfComponents(1);
+            //    kerr->SetName("Error");
+
             for (auto leaf : LeafDict) {
                 kr->InsertTuple1( leaf.first, std::real(leaf.second) );
                 ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
+                kid->InsertTuple1( leaf.first, leaf.first );
             }
+
+            //for (auto leaf : LeafDictErr) {
+            //    kerr->InsertTuple1( leaf.first, std::imag(leaf.second) );
+            //}
+
             oct->GetLeafData()->AddArray(kr);
             oct->GetLeafData()->AddArray(ki);
+            oct->GetLeafData()->AddArray(km);
+            oct->GetLeafData()->AddArray(kid);
+            //oct->GetLeafData()->AddArray(kerr);
 
             auto write = vtkXMLHyperOctreeWriter::New();
                 //write.SetDataModeToAscii()
                 write->Write();
                 write->Delete();
 
+            //kerr->Delete();
+            kid->Delete();
             kr->Delete();
             ki->Delete();
+            km->Delete();
             curse->Delete();
             oct->Delete();
         #else
     //      Method:  f
     //--------------------------------------------------------------------------------------
     Complex KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
-        return Complex(volume*Ht.dot(Hr));
-        //return ComputeV0Cell(MU0*Ht, MU0*Hr, volume, 1.0);
+        //return Complex(volume*Ht.dot(Hr));
+        return ComputeV0Cell(MU0*Ht, MU0*Hr, volume, 1.0);
     }
 
     //--------------------------------------------------------------------------------------
                 const Vector3cr& Br, const Real& vol, const Real& phi) {
 
         // Compute the elliptic fields
-        Vector3r B0hat = {1,0,0};
-        Vector3r B0 = 53000 * B0hat; // nT
+        Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
+        Vector3r B0 = SigmaModel->GetMagneticField();
+
+        // Elliptic representation
         EllipticB EBT = EllipticFieldRep(Bt, B0hat);
         EllipticB EBR = EllipticFieldRep(Br, B0hat);
 
         Real Mn0Abs = Mn0.norm();
 
         Real Taup = 0.020; // s
-        Real Ip = 10;      // A
+        Real Ip   = 10;      // A
 
         // Compute the tipping angle
         Real sintheta = std::sin(0.5*GAMMA*Ip*Taup*std::abs(EBT.alpha-EBT.beta));
         return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
     }
 
+    //--------------------------------------------------------------------------------------
+    //       Class:  KernelV0
+    //      Method:  ComputeV0Cell
+    //--------------------------------------------------------------------------------------
     Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
-        Real Temperature = 283; // in K
         Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
         return chi_n*Porosity*B0;
     }