Fixes elliptic field calculation bug

examples/KernelV0-2.cpp

         Kern->PushCoil( "Coil 2", Rx1 );
         Kern->SetLayeredEarthEM( earth );
 
-        Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
+        Kern->SetIntegrationSize( (Vector3r() << 200,200,100).finished() );
         Kern->SetIntegrationOrigin( (Vector3r() << -100, -100, .5).finished() );
-        Real tol(1e-9);
+        Real tol(1e-11); // 13
         Kern->SetTolerance( tol ); // 1e-12
 
 //         Kern->AlignWithAkvoDataset( YAML::LoadFile(argv[2]) );
         Kern->SetPulseCurrent( I ); // nbins, low, high
 
         //VectorXr interfaces = VectorXr::LinSpaced( 41, .5, 45.5 ); // nlay, low, high
-        VectorXr interfaces = VectorXr::LinSpaced( 61, .5, 45.5 ); // nlay, low, high
+        //VectorXr interfaces = VectorXr::LinSpaced( 61, .5, 45.5 ); // nlay, low, high
+        VectorXr interfaces = VectorXr::LinSpaced( 2, .5, 45.5 ); // nlay, low, high
         Real thick = .5;
         for (int ilay=1; ilay<interfaces.size(); ++ilay) {
             interfaces(ilay) = interfaces(ilay-1) + thick;
     // 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")};
-    Kern->CalculateK0( tx, rx, false ); // 3rd argument is vtk output
+    Kern->CalculateK0( tx, rx, true ); // 3rd argument is vtk output
 
     std::ofstream dout = std::ofstream(std::string("Rx-")+std::string(argv[3])+std::string(".dat"));
     dout << "# Transmitters: ";

include/KernelV0.h

 namespace Lemma {
 
     struct EllipticB {
-        Vector3r  bhat;
-        Vector3r  bhatp;
         Real      alpha;
         Real      beta;
-        Complex   eizt;
         Real      zeta;
+        Real      err;
+        Complex   eizt;
+//        Complex   BperpdotB;
+        Vector3r  bhat;
+        Vector3r  bhatp;
+//        Vector3cr Bperp;
     };
 
     template <typename T> int sgn(T val) {
 
         int                                       ilay;
         int                                       nleaves;
-        int                                       minLevel=4;
-        int                                       maxLevel=10;
+        int                                       minLevel=0;
+        int                                       maxLevel=12;
 
         Real                                      VOLSUM;
         Real                                      tol=1e-11;
 
         // Physical constants and conversion factors
         static constexpr Real GAMMA = 2.67518e8;                  // MKS units
-        static constexpr Real INVSQRT2 = 0.70710678118654746;
+        static constexpr Real INVSQRT2 = 0.70710678118654746;     // 1/sqrt(2)
         static constexpr Real HBAR = 1.05457148e-34;              // m2 kg / s
         static constexpr Real NH2O = 6.692e28;                    // [m^3]
         static constexpr Real KB = 1.3805e-23;                    // m^2 kg s-2 K-1

src/KernelV0.cpp

     //       Class:  KernelV0
     //      Method:  f
     //--------------------------------------------------------------------------------------
+#if 1
     VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
 
         // Compute the elliptic fields
 
         // Compute Mn0
         Vector3r Mn0 = ComputeMn0(1.0, B0);
+        //std::cout << "Mn0\t" << Mn0.transpose() << std::endl;
         Real Mn0Abs = Mn0.norm();
 
         // Compute phase delay
         VectorXcr F = VectorXcr::Zero( PulseI.size() );
         for (int iq=0; iq<PulseI.size(); ++iq) {
             // Compute the tipping angle
-            Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta)); // why std::abs
+            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;
+            //TODO TEST FOR ASYMETRY
+            //Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta));
+            //F(iq) = volume * Complex(EBT.Bperp.real().norm(), EBT.Bperp.imag().norm()); //Complex(sintheta, EBT.Bperp.norm() );
+            //F(iq) = volume * Complex(EBT.alpha, EBT.beta);
+            //F(iq) = volume * EBT.err;
+            //F(iq) = volume * sintheta;
         }
+
+        return F;
+    }
+#endif
+#if 0
+    VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
+        VectorXcr F = VectorXcr::Ones( PulseI.size() );
+        F.array() *= volume * Complex(Ht.norm(), Hr.norm()); //*Ht.dot(Hr);
         return F;
     }
+#endif
 
 //     //--------------------------------------------------------------------------------------
 //     //       Class:  KernelV0
     //      Method:  ComputeV0Cell
     //--------------------------------------------------------------------------------------
     EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
+        // This all follows Weichman et. al., 2000. There are some numerical stability isseus
+        // below. Reformulating may be welcome, may be in B field calculation too.
         EllipticB ElipB = EllipticB();
-        Vector3cr Bperp = B.array() - B0hat.dot(B)*B0hat.array();
-        Real BperpNorm = Bperp.norm();
+        Vector3cr Bperp = B - B0hat.dot(B)*B0hat; // complex - real??
+        //ElipB.BperpdotB = Bperp.dot(B0hat);       // TODO remove
+        Real BperpNorm  = Bperp.norm();
         Complex Bp2 = Bperp.transpose() * Bperp;
         VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
         ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
         ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
-        ElipB.beta = sgn(std::real(iB0.dot(Bperp.cross(Bperp.conjugate())))) *
-                (INVSQRT2)*std::sqrt(std::abs(BperpNorm*BperpNorm-std::abs(Bp2)));
+        ElipB.beta = std::copysign(1, std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
+                     (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2)));
         ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
         ElipB.bhatp = B0hat.cross(ElipB.bhat);
         ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
+        /* use as an error check decomposed field - computed actual */
+        Vector3cr Bperp2 = ElipB.eizt * (ElipB.alpha * ElipB.bhat
+                       + (Complex(0,1) * ElipB.beta * ElipB.bhatp) );
+        ElipB.err = (Bperp-Bperp2).norm();
+        if (ElipB.err > 1e-12) {
+            std::cout << "Elip error\n";
+            std::cout << "Bperp \t" << Bperp.transpose() << std::endl;
+            std::cout << "Bperp2\t" << Bperp2.transpose() << std::endl;
+            std::cout << "err   \t" << ElipB.err << std::endl;
+        }
+        //std::cout << "B0\t" << B0hat.transpose() << std::endl;
         return ElipB;
     }
 
         Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
         Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
         for ( auto EMCalc : EMEarths ) {
+
+
+
             EMCalc.second->GetFieldPoints()->ClearFields();
             EMCalc.second->CalculateWireAntennaFields();
             switch (EMCalc.second->GetTxRxMode()) {