KiHa testing tweaks

Modules/FDEM1D/examples/EMDipEarth1D.cpp

@@ -45,7 +45,7 @@ int main() {
 		dipole->SetFrequency(0, 4400.1000);
 		//dipole->SetPhase(0);
 		//dipole->SetLocation( (VectorXr(3) << 49, 49, -1e-4).finished() );
-		dipole->SetLocation( 49, 49, -1e-4  );
+		dipole->SetLocation( 0, 0, -1e-4  );
 
 
 	// Define model
@@ -74,6 +74,7 @@ int main() {
 		int ir = 0;
 		for (int iz=0; iz<nz; ++iz) {
 			loc << ox, oy, depth;
+            std::cout << "Receiver location " << loc.transpose() << std::endl;
 			receivers->SetLocation(ir, loc);
 			depth += dz;
 			++ ir;
@@ -113,7 +114,7 @@ int main() {
 
  #ifdef KIHALEE_EM1D
 	receivers->ClearFields();
-    std::cout << "\nFORTRAN\n";
+    std::cout << "\nFORTRAN KiHa\n";
  	EmEarth->MakeCalc();
     auto fc = receivers->GetEfield(0,0);
 // 	std::cout << receivers->GetHfield(0,0) << std::endl;

Modules/FDEM1D/src/DipoleSource.cpp

@@ -770,6 +770,7 @@ namespace Lemma {
     }
 
     void DipoleSource::UpdateFields( const int& ifreq, HankelTransform* Hankel, const Real& wavef) {
+
         Vector3r Pol = Phat;
 
         switch (Type) {
@@ -813,7 +814,7 @@ namespace Lemma {
                                      0. );
                     } // Fields to calculate Z polarity Electric dipole
                 }
-                if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // y dipole
+                if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y dipole
                     switch(FieldsToCalculate) {
                         case E:
                             f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
@@ -823,9 +824,9 @@ namespace Lemma {
                             f(4) = Hankel->Zgauss(4, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4])) / KernelManager->GetRAWKernel(ik[4])->GetYm();
                             if (std::abs(Pol[1]) > 0) {
                                 this->Receivers->AppendEfield(ifreq, irec,
-                                    Pol[1]*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho))*Moment,
-                                    Pol[1]*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho))*Moment,
-                                    Pol[1]*QPI*sp*f(4)*Moment);
+                                    Pol[1]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
+                                    Pol[1]*Moment*QPI*((sps*f(0)+c2p*f(1)/rho)-(cps*f(2)-c2p*f(3)/rho)),
+                                    Pol[1]*Moment*QPI*sp*f(4));
                             }
                             if (std::abs(Pol[0]) > 0) {
                                 this->Receivers->AppendEfield(ifreq, irec,

Modules/FDEM1D/src/EMEarth1D.cpp

@@ -755,9 +755,12 @@ namespace Lemma {
 
     void EMEarth1D::SolveSingleTxRxPair (const int &irec, HankelTransform *Hankel, const Real &wavef, const int &ifreq,
                    DipoleSource *tDipole) {
-        //std::cout << "SolveSingleTxRxPair" << std::endl;
         ++icalcinner;
+
+        // The PGI compilers fail on the below line, and others like it.
         Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
+        //Real rho = ( ((Receivers->GetLocation(irec) - tDipole->GetLocation()).head(2)).eval() ).norm();
+
         tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
         Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
         tDipole->UpdateFields( ifreq,  Hankel, wavef );

Modules/FDEM1D/src/FHTAnderson801.cpp

@@ -998,7 +998,8 @@ namespace Lemma {
         // in release.
         #ifndef NDEBUG
 		if (rho<=0) {
-			throw std::runtime_error("In Hankel 2 Argument rho < 0.");
+            //std::cout << "rho= " << rho << std::endl;
+			throw std::runtime_error("In Hankel 2 Argument rho <= 0; rho=" + to_string(rho) );
 		}
 
 		if (this->NumConv<1) {

Modules/FDEM1D/src/GQChave.cpp

@@ -325,10 +325,9 @@ namespace Lemma{
         return CName;
     }		// -----  end of method GQChave::GetName  -----
 
-    Complex GQChave::
-            Zgauss(const int &ikk, const EMMODE &mode,
-                    const int &itype, const Real &rho, const Real &wavef,
-                    KernelEM1DBase* Kernel){
+    Complex GQChave::Zgauss(const int &ikk, const EMMODE &mode,
+            const int &itype, const Real &rho, const Real &wavef,
+            KernelEM1DBase* Kernel){
 
         // TI, TODO, change calls to Zgauss to reflect this, go and fix so we
         // dont subract 1 from this everywhere