Added kiha library check, things seem to be working, need more verification and benching

CMakeLists.txt

 # Lemma Configuration
 #####################
 
+#####################
+# Look for Ki Ha Lee 
+#####################
+FIND_LIBRARY( KIHA_EM1D kihaem1d ) # PATHS ${/home/tirons/local/lib}  )
+if ( KIHA_EM1D )
+	add_compile_options(-DKIHALEE_EM1D)
+endif()
+
 ####################
 # Add the c++11 flag 
 # TODO add compiler specific instructions

Modules/FDEM1D/CMakeLists.txt

 # Linking
 target_link_libraries(fdem1d "lemmacore")
 
+if ( KIHA_EM1D )
+	target_link_libraries(fdem1d "kihaem1d")
+	target_link_libraries(fdem1d "gfortran")
+endif()
+
 # Testing
 if (LEMMA_ENABLE_TESTING)
 	add_subdirectory(testing)

Modules/FDEM1D/examples/CMakeLists.txt

 add_executable( PolygonalWireAntenna PolygonalWireAntenna.cpp  )
 target_link_libraries(  PolygonalWireAntenna  "lemmacore" "fdem1d")
 
+add_executable( EMDipEarth1D EMDipEarth1D.cpp  )
+target_link_libraries(  EMDipEarth1D  "lemmacore" "fdem1d")
+
 INSTALL_TARGETS( "/share/FEM1D/"
 	LayeredEarthEM 
 	FieldPoints
 	PolygonalWireAntenna
+	EMDipEarth1D
 )

Modules/FDEM1D/examples/emdipearth1d.cpp → Modules/FDEM1D/examples/EMDipEarth1D.cpp

 //  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 //
 // ===========================================================================
-#include "emearth1d.h"
-#include "dipolesource.h"
-#include "layeredearthem.h"
-#include "receiverpoints.h"
+#include <FDEM1D>
 
 using namespace Lemma;
 
 int main() {
 
 	// Test with a single dipole
-	DipoleSource *dipole = DipoleSource::New();
+	auto dipole = DipoleSource::NewSP();
 		dipole->SetType(GROUNDEDELECTRICDIPOLE);
 		dipole->SetPolarisation(XPOLARISATION);
 		dipole->SetNumberOfFrequencies(1);
 	VectorXr  thick(1);
 		thick << 10;//, 10, 10;
 
-	LayeredEarthEM *earth = LayeredEarthEM::New();
+	auto earth = LayeredEarthEM::NewSP();
         earth->SetNumberOfLayers(2);
 		earth->SetLayerConductivity(sigma);
 		//earth->SetLayerThickness(thick);
 
 	// Receivers
-	ReceiverPoints *receivers = ReceiverPoints::New();
+	auto receivers = FieldPoints::NewSP();
 		Vector3r loc;
 		//Real ox = -5.1456;
 		//Real oy =  2.2350;
 		Real dz = 2.6;
 		int  nz = 1;
 
-		receivers->SetNumberOfReceivers(nz);
+		receivers->SetNumberOfPoints(nz);
 		int ir = 0;
 		for (int iz=0; iz<nz; ++iz) {
 			loc << ox, oy, depth;
 		}
 	//receivers->SetLocation(1, ox, oy, depth);
 
-	EMEarth1D *EmEarth = EMEarth1D::New();
+    auto EmEarth = EMEarth1D::NewSP();
         //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
         EmEarth->SetFieldsToCalculate(BOTH); // Fortran needs this
 		EmEarth->AttachDipoleSource(dipole);
 		EmEarth->AttachLayeredEarthEM(earth);
-		EmEarth->AttachReceiverPoints(receivers);
+		EmEarth->AttachFieldPoints(receivers);
 
         //dipole->SetType(ELECTRICDIPOLE);
 //         receivers->SetNumberOfReceivers(1);
 	receivers->ClearFields();
 
     // swap tx rx posigion
+    /*
     receivers->SetLocation(0, dipole->GetLocation());
     dipole->SetLocation(loc);
 	EmEarth->MakeCalc3();
 	std::cout << receivers->GetEfield(0,0) << std::endl;
-
 	receivers->ClearFields();
+    */
+
  #ifdef KIHALEE_EM1D
-// 	//std::cout << "\nFORTRAN\n";
+    std::cout << "\nFORTRAN\n";
  	EmEarth->MakeCalc();
 // 	std::cout << receivers->GetHfield(0,0) << std::endl;
  	std::cout << receivers->GetEfield(0,0) << std::endl;
  #endif
 
-    dipole->Delete();
-    EmEarth->Delete();
-    receivers->Delete();
-    earth->Delete();
-
-// 	try {
-//
-// 	}
-//
-// 	// Catch exceptions
-// 	catch (NullEarth) {
-// 			std::cout << "Earth model not set dummy\n";
-// 	}
-//
-// 	catch (NullReceivers) {
-// 			std::cout << "Receivers not set dummy\n";
-// 	}
-
-	//EmEarth->TestBess();
-
 }

Modules/FDEM1D/include/EMEarth1D.h

 
             /// Calculates the field(s) due to an ungrounded dipole source
             /// Calls FORTRAN library em1d (em1dnew.for)
-#ifdef COMPILE_FORTRAN
+#ifdef KIHALEE_EM1D
             void MakeCalc();
 #endif
 

Modules/FDEM1D/include/FDEM1D

 
 #include "DipoleSource.h"
 
+#include "EMEarth1D.h"
+
 // internal
 //#include "KernelEM1DManager.h"
 //#include "KernelEM1DSpec.h"

Modules/FDEM1D/src/EMEarth1D.cpp

+/* This file is part of Lemma, a geophysical modelling and inversion API */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     12/02/2009
+ **/
+
+#include "EMEarth1D.h"
+#include "FieldPoints.h"
+#include "WireAntenna.h"
+#include "PolygonalWireAntenna.h"
+
+#ifdef LEMMAUSEOMP
+#include "omp.h"
+#endif
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const EMEarth1D &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc ---
+        return stream;
+    }
+
+#ifdef KIHALEE_EM1D
+    // Wrapper function for Fortran subroutine Em1D bi kihand
+    // Returns E or H fields (SLOW)
+    extern "C" { void em1dcall_(int &itype,   // source
+                            int     &ipol,    // source
+                            int     &nlay,    // Earth
+                            int     &nfreq,   // source
+                            int     &nfield,  // Calculator
+                            int     &nres,    // Receivers
+                            int     &jtype,   // N/A
+                            int     &jgamma,  // Controller
+                            double  &acc,     // Controller
+                            double  *dep,     // Earth
+                            std::complex<double> *sig,     // Earth
+                            double  *susl,    // Earth
+                            double  *sush,    // Earth
+                            double  *sustau,  // Earth
+                            double  *susalp,  // Earth
+                            double  *eprl,    // Earth
+                            double  *eprh,    // Earth
+                            double  *eprtau,  // Earth
+                            double  *epralp,  // Earth
+                            double  &finit,   // N/A
+                            double  &flimit,  // N/A
+                            double  &dlimit,  // N/A
+                            double  &lfinc,   // N/A
+                            double  &tx,      // Source
+                            double  &ty,      // Source
+                            double  &tz,      // Source
+                            double  *rxx,     // Receivers
+                            double  *rxy,     // Receivers
+                            double  *rxz,     // Receivers
+                            std::complex<double> *ex,      // Receivers
+                            std::complex<double> *ey,      //    |
+                            std::complex<double> *ez,      //    |
+                            std::complex<double> *hx,      //    |
+                            std::complex<double> *hy,      //    V
+                            std::complex<double> *hz );    //   ___
+    }
+#endif
+
+    // ====================  LIFECYCLE     ===================================
+
+    // TODO init large arrays here.
+    EMEarth1D::EMEarth1D( const ctor_key& ) : LemmaObject( ),
+            Dipole(nullptr), Earth(nullptr), Receivers(nullptr), Antenna(nullptr),
+            FieldsToCalculate(BOTH), HankelType(ANDERSON801), icalcinner(0), icalc(0)
+        //#ifdef HAVEBOOSTPROGRESS
+        //    , disp(0)
+        //#endif
+        {
+    }
+
+    EMEarth1D::~EMEarth1D() {
+    }
+
+    std::shared_ptr<EMEarth1D> EMEarth1D::NewSP() {
+        return std::make_shared<EMEarth1D>(ctor_key());
+    }
+
+    YAML::Node EMEarth1D::Serialize() const {
+        YAML::Node node = LemmaObject::Serialize();
+        node["FieldsToCalculate"] = enum2String(FieldsToCalculate);
+        node["HankelType"] = enum2String(HankelType);
+        //if (Dipole != nullptr) node["Dipole"] = Dipole->Serialize();
+        if (Earth != nullptr)     node["Earth"] = Earth->Serialize();
+        //if (Receivers != nullptr) node["Receivers"] = Receivers->Serialize(); Can be huge?
+        if (Antenna != nullptr)   node["Antenna"] = Antenna->Serialize();
+        node.SetTag( this->GetName() );
+        return node;
+    }
+
+    // ====================  ACCESS        ===================================
+    void EMEarth1D::AttachDipoleSource( std::shared_ptr<DipoleSource> dipoleptr) {
+        Dipole = dipoleptr;
+    }
+
+    void EMEarth1D::AttachLayeredEarthEM( std::shared_ptr<LayeredEarthEM> earthptr) {
+        Earth = earthptr;
+    }
+
+    void EMEarth1D::AttachFieldPoints( std::shared_ptr<FieldPoints> recptr) {
+
+        Receivers = recptr;
+        if (Receivers == nullptr) {
+            std::cout << "nullptr Receivers in emearth1d.cpp " << std::endl;
+            return;
+        }
+
+        // This has an implicid need to first set a source before receivers, users
+        // will not expect this. Fix
+        if (Dipole != nullptr) {
+            switch (FieldsToCalculate) {
+                case E:
+                    Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
+                    break;
+                case H:
+                    Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
+                    break;
+                case BOTH:
+                    Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
+                    Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
+                    break;
+            }
+        } else if (Antenna != nullptr) {
+            switch (FieldsToCalculate) {
+                case E:
+                    Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
+                    break;
+                case H:
+                    Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
+                    break;
+                case BOTH:
+                    Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
+                    Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
+                    break;
+            }
+        }
+    }
+
+    void EMEarth1D::AttachWireAntenna(std::shared_ptr<WireAntenna> antennae) {
+        this->Antenna = antennae;
+    }
+
+    void EMEarth1D::SetFieldsToCalculate(const FIELDCALCULATIONS &calc) {
+        FieldsToCalculate = calc;
+    }
+
+    void EMEarth1D::SetHankelTransformMethod( const HANKELTRANSFORMTYPE &type) {
+        HankelType = type;
+    }
+
+    void EMEarth1D::Query() {
+        std::cout << "EmEarth1D::Query()" << std::endl;
+
+        std::cout << "Dipole " << Dipole;
+        if (Dipole) std::cout << *Dipole << std::endl;
+
+        std::cout << "Earth " << Earth;
+        if (Earth) std::cout << *Earth << std::endl;
+
+        std::cout << "Receivers " << Earth;
+        if (Earth) std::cout << *Receivers << std::endl;
+
+        std::cout << "Antenna " << Earth;
+        if (Antenna) std::cout << *Antenna << std::endl;
+
+        std::cout << "icalc " << icalc << std::endl;
+
+        std::cout << "icalcinner " << icalcinner << std::endl;
+    }
+
+    // ====================  OPERATIONS    ===================================
+
+    void EMEarth1D::CalculateWireAntennaFields(bool progressbar) {
+
+        #ifdef HAVEBOOSTPROGRESS
+        boost::progress_display *disp;
+        #endif
+
+        if (Earth == nullptr) {
+            throw NullEarth();
+        }
+        if (Receivers == nullptr) {
+            throw NullReceivers();
+        }
+        if (Antenna == nullptr) {
+            throw NullAntenna();
+        }
+        if (Dipole != nullptr) {
+            throw DipoleSourceSpecifiedForWireAntennaCalc();
+        }
+
+        Receivers->ClearFields();
+
+        // Check to make sure Receivers are set up for all calculations
+        switch(FieldsToCalculate) {
+            case E:
+                if (Receivers->NumberOfBinsE != Antenna->GetNumberOfFrequencies())
+                    Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
+                break;
+            case H:
+                if (Receivers->NumberOfBinsH != Antenna->GetNumberOfFrequencies())
+                    Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
+                break;
+            case BOTH:
+                if (Receivers->NumberOfBinsH != Antenna->GetNumberOfFrequencies())
+                    Receivers->SetNumberOfBinsH(Antenna->GetNumberOfFrequencies());
+                if (Receivers->NumberOfBinsE != Antenna->GetNumberOfFrequencies())
+                    Receivers->SetNumberOfBinsE(Antenna->GetNumberOfFrequencies());
+                break;
+        }
+
+        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
+            /* TODO FIX THIS ISSUES */
+            if (Antenna->IsHorizontallyPlanar() && HankelType == ANDERSON801) {
+                //std::cout << "Lag baby lag" << std::endl;
+                for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies();++ifreq) {
+                    //std::cout << "Num Recs" <<  Receivers->GetNumberOfPoints() << std::endl;
+                    Real wavef = 2.*PI* Antenna->GetFrequency(ifreq);
+                    #ifdef LEMMAUSEOMP
+                    #pragma omp parallel
+                    {
+                    #endif
+                    auto Hankel = FHTAnderson801::NewSP();
+                    #ifdef LEMMAUSEOMP
+                    #pragma omp for schedule(static, 1)
+                    #endif
+                    for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
+                    //for (int irec=0; irec<2; ++irec) { // TODO FIXME BELO
+                        auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
+                        SolveLaggedTxRxPair(irec, Hankel.get(), wavef, ifreq, AntCopy.get());
+                        //exit(0);
+                    }
+                    //Receivers->ClearFields(); // FIXME DEBUG TODO
+                    #ifdef LEMMAUSEOMP
+                    }
+                    #endif
+                }
+            } else
+            if (Receivers->GetNumberOfPoints() > Antenna->GetNumberOfFrequencies()) {
+
+                //std::cout << "freq parallel #1" << std::endl;
+                //** Progress display bar for long calculations */
+                #ifdef HAVEBOOSTPROGRESS
+                if (progressbar) {
+                    disp = new boost::progress_display( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
+                }
+                #endif
+
+                // parallelise across receivers
+                #ifdef LEMMAUSEOMP
+                #pragma omp parallel
+                #endif
+                { // OpenMP Parallel Block
+                    // Since these antennas change we need a local copy for each
+                    // thread.
+                    auto AntCopy = static_cast<PolygonalWireAntenna*>(Antenna.get())->ClonePA();
+
+                    std::shared_ptr<HankelTransform> Hankel;
+                    switch (HankelType) {
+                        case ANDERSON801:
+                            Hankel = FHTAnderson801::NewSP();
+                            break;
+                        case CHAVE:
+                            Hankel = GQChave::NewSP();
+                            break;
+                        case FHTKEY201:
+                            Hankel = FHTKey201::NewSP();
+                            break;
+                        case FHTKEY101:
+                            Hankel = FHTKey101::NewSP();
+                            break;
+                        case FHTKEY51:
+                            Hankel = FHTKey51::NewSP();
+                            break;
+                        case QWEKEY:
+                            Hankel = QWEKey::NewSP();
+                            break;
+                        default:
+                            std::cerr << "Hankel transform cannot be created\n";
+                            exit(EXIT_FAILURE);
+                    }
+
+                    //for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
+                    #ifdef LEMMAUSEOMP
+                    #pragma omp for schedule(static, 1) //nowait
+                    #endif
+                    for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
+                        if (!Receivers->GetMask(irec)) {
+                            AntCopy->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
+                            for (int idip=0; idip<AntCopy->GetNumberOfDipoles(); ++idip) {
+                                auto tDipole = AntCopy->GetDipoleSource(idip);
+                                //#ifdef LEMMAUSEOMP
+                                //#pragma omp for schedule(static, 1)
+                                //#endif
+                                for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies();
+                                        ++ifreq) {
+                                    // Propogation constant in free space
+                                    Real wavef   = tDipole->GetAngularFrequency(ifreq) *
+                                          std::sqrt(MU0*EPSILON0);
+                                    SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
+                                } // freq loop
+                            } // dipole loop
+                        } // mask
+                        //std::cout << "Normal Path\n";
+                        //std::cout << Receivers->GetHfield(0, irec) << std::endl;
+                        //if (irec == 1) exit(0);
+                        #ifdef HAVEBOOSTPROGRESS
+                        if (progressbar) ++(*disp);
+                        #endif
+                    } // receiver loop
+                } // OMP_PARALLEL BLOCK
+            } else if (Antenna->GetNumberOfFrequencies() > 8) {
+                // parallel across frequencies
+                //std::cout << "freq parallel #2" << std::endl;
+                for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
+                    if (!Receivers->GetMask(irec)) {
+                        static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
+                        #ifdef LEMMAUSEOMP
+                        #pragma omp parallel
+                        #endif
+                        { // OpenMP Parallel Block
+
+                            std::shared_ptr<HankelTransform> Hankel;
+                            switch (HankelType) {
+                                case ANDERSON801:
+                                    Hankel = FHTAnderson801::NewSP();
+                                    break;
+                                case CHAVE:
+                                    Hankel = GQChave::NewSP();
+                                    break;
+                                case FHTKEY201:
+                                    Hankel = FHTKey201::NewSP();
+                                    break;
+                                case FHTKEY101:
+                                    Hankel = FHTKey101::NewSP();
+                                    break;
+                                case FHTKEY51:
+                                    Hankel = FHTKey51::NewSP();
+                                    break;
+                                case QWEKEY:
+                                    Hankel = QWEKey::NewSP();
+                                    break;
+                                default:
+                                    std::cerr << "Hankel transform cannot be created\n";
+                                    exit(EXIT_FAILURE);
+                            }
+                            #ifdef LEMMAUSEOMP
+                            #pragma omp for schedule(static, 1)
+                            #endif
+                            for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
+                                for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
+                                    auto tDipole = Antenna->GetDipoleSource(idip);
+                                    // Propogation constant in free space
+                                    Real wavef   = tDipole->GetAngularFrequency(ifreq) *
+                                          std::sqrt(MU0*EPSILON0);
+                                    SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
+                                } // dipole loop
+                            } // frequency loop
+                        } // OMP_PARALLEL BLOCK
+                    } // mask loop
+                    #ifdef HAVEBOOSTPROGRESS
+                    //if (Receivers->GetNumberOfPoints() > 100) {
+                    //    ++ disp;
+                    //}
+                    #endif
+                } // receiver loop
+                //std::cout << "End freq parallel " << std::endl;
+            } // Frequency Parallel
+              else {
+                //std::cout << "parallel across #3 " << std::endl;
+                for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
+                    if (!Receivers->GetMask(irec)) {
+
+                        static_cast<PolygonalWireAntenna*>(Antenna.get())->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
+//                         std::cout << "Not Masked " << std::endl;
+//                         std::cout << "n Freqs " << Antenna->GetNumberOfFrequencies() << std::endl;
+//                         std::cout << "n Dipoles " << Antenna->GetNumberOfDipoles() << std::endl;
+//                         if ( !Antenna->GetNumberOfDipoles() ) {
+//                             std::cout << "NO DIPOLES!!!!!!!!!!!!!!!!!!!!!!!!!!\n";
+// //                            std::cout << "rec location " << Receivers->GetLocation(irec) << std::endl;
+// //                        }
+
+                        #ifdef LEMMAUSEOMP
+                        #pragma omp parallel
+                        #endif
+                        { // OpenMP Parallel Block
+                            std::shared_ptr<HankelTransform> Hankel;
+                            switch (HankelType) {
+                                case ANDERSON801:
+                                    Hankel = FHTAnderson801::NewSP();
+                                    break;
+                                case CHAVE:
+                                    Hankel = GQChave::NewSP();
+                                    break;
+                                case FHTKEY201:
+                                    Hankel = FHTKey201::NewSP();
+                                    break;
+                                case FHTKEY101:
+                                    Hankel = FHTKey101::NewSP();
+                                    break;
+                                case FHTKEY51:
+                                    Hankel = FHTKey51::NewSP();
+                                    break;
+                                case QWEKEY:
+                                    Hankel = QWEKey::NewSP();
+                                    break;
+                                default:
+                                    std::cerr << "Hankel transform cannot be created\n";
+                                    exit(EXIT_FAILURE);
+                            }
+                            for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
+                                #ifdef LEMMAUSEOMP
+                                #pragma omp for schedule(static, 1)
+                                #endif
+                                for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
+                                    //#pragma omp critical
+                                    //{
+                                    //cout << "idip=" << idip << "\tthread num=" << omp_get_thread_num() << '\n';
+                                    //}
+                                    auto tDipole = Antenna->GetDipoleSource(idip);
+                                    // Propogation constant in free space
+                                    Real wavef   = tDipole->GetAngularFrequency(ifreq) *
+                                          std::sqrt(MU0*EPSILON0);
+                                    SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
+                                } // dipole loop
+                            } // frequency loop
+                        } // OMP_PARALLEL BLOCK
+                    } // mask loop
+                    #ifdef HAVEBOOSTPROGRESS
+                    //if (Receivers->GetNumberOfPoints() > 100) {
+                    //    ++ disp;
+                    //}
+                    #endif
+                } // receiver loop
+           } // Polygonal parallel logic
+        } else {
+             std::cerr << "Lemma with WireAntenna class is currently broken"
+                  << " fix or use PolygonalWireAntenna\n" << std::endl;
+             exit(EXIT_FAILURE);
+            // TODO, getting wrong answer, curiously worKernel->GetKs() with MakeCalc, maybe
+            // a threading issue, use SolveSingleTxRxPair maype instead of call
+            // to MakeCalc3? !!!
+            for (int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
+                this->Dipole = Antenna->GetDipoleSource(idip);
+                MakeCalc3();
+                //++disp;
+            }
+            this->Dipole = nullptr;
+        }
+
+        #ifdef HAVEBOOSTPROGRESS
+        if (progressbar) {
+            delete disp;
+        }
+        #endif
+    }
+
+    #ifdef KIHALEE_EM1D
+    void EMEarth1D::MakeCalc() {
+
+        int itype;       // 1 = elec, 2 = mag
+        switch (this->Dipole->GetDipoleSourceType()) {
+            case (GROUNDEDELECTRICDIPOLE) :
+                itype = 1;
+                break;
+            case (MAGNETICDIPOLE) :
+                itype = 2;
+                break;
+            case (UNGROUNDEDELECTRICDIPOLE) :
+                std::cerr << "Fortran routine cannot calculate ungrounded"
+                             "electric dipole\n";
+            default:
+                throw NonValidDipoleType();
+        }
+
+        int ipol ;
+        Vector3r Pol = this->Dipole->GetPolarisation();
+        if (std::abs(Pol[0]-1) < 1e-5) {
+            ipol = 1;
+        } else if (std::abs(Pol[1]-1) < 1e-5) {
+            ipol = 2;
+        } else if (std::abs(Pol[2]-1) < 1e-5) {
+            ipol = 3;
+        } else {
+            std::cerr << "Fortran routine cannot calculate arbitrary "
+                             "dipole polarisation, set to x, y, or z\n";
+        }
+
+        int nlay   =  Earth->GetNumberOfNonAirLayers();
+
+        if (nlay > MAXLAYERS) {
+            std::cerr << "FORTRAN CODE CAN ONLY HANDLE " << MAXLAYERS
+                      << " LAYERS\n";
+            throw  EarthModelWithMoreThanMaxLayers();
+        }
+
+        int nfreq  = 1;     // number of freqs
+
+        int nfield;     // field output 1 = elec, 2 = mag, 3 = both
+        switch (FieldsToCalculate) {
+            case E:
+                nfield = 1;
+                break;
+            case H:
+                nfield = 2;
+                break;
+            case BOTH:
+                nfield = 3;
+                break;
+            default:
+                throw 7;
+        }
+
+        int nres   = Receivers->GetNumberOfPoints();
+        int jtype  = 3;     // form ouf output,
+                            // 1 = horizontal,
+                            // 2 = down hole,
+                            // 3 = freq sounding
+                            // 4 = down hole logging
+
+        int jgamma = 0;     // Units 0 = MKS (H->A/m and E->V/m)
+                            //       1 = h->Gammas   E->V/m
+
+        double acc = 0.;    // Tolerance
+
+        // TODO, fix FORTRAN calls so these arrays can be nlay long, not
+        // MAXLAYERS.
+
+        // Model Parameters
+        double *dep    = new double[MAXLAYERS];
+        dep[0] = 0.;          // We always say air starts at 0
+        for (int ilay=1; ilay<Earth->GetNumberOfLayers(); ++ilay) {
+            dep[ilay] = dep[ilay-1] + Earth->GetLayerThickness(ilay);
+            //std::cout << "Depth " << dep[ilay] << std::endl;
+        }
+
+        std::complex<double> *sig = new std::complex<double> [MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            sig[ilay-1] = (std::complex<double>)(Earth->GetLayerConductivity(ilay));
+        }
+
+        // TODO, pass these into Fortran call, and return Cole-Cole model
+        // parameters. Right now this does nothing
+        //std::complex<double> *sus   = new std::complex<double>[MAXLAYERS];
+        //std::complex<double> *epr   = new std::complex<double>[MAXLAYERS];
+
+        // Cole-Cole model stuff
+        double *susl   = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            susl[ilay-1] = Earth->GetLayerLowFreqSusceptibility(ilay);
+        }
+
+        double *sush   = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            sush[ilay-1] = Earth->GetLayerHighFreqSusceptibility(ilay);
+        }
+
+        double *sustau = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            sustau[ilay-1] = Earth->GetLayerTauSusceptibility(ilay);
+        }
+
+        double *susalp = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            susalp[ilay-1] = Earth->GetLayerBreathSusceptibility(ilay);
+        }
+
+        double *eprl   = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            eprl[ilay-1] = Earth->GetLayerLowFreqPermitivity(ilay);
+        }
+
+        double *eprh   = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            eprh[ilay-1] = Earth->GetLayerHighFreqPermitivity(ilay);
+        }
+
+        double *eprtau = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            eprtau[ilay-1] = Earth->GetLayerTauPermitivity(ilay);
+        }
+
+        double *epralp = new double[MAXLAYERS];
+        for (int ilay=1; ilay<=nlay; ++ilay) {
+            epralp[ilay-1] = Earth->GetLayerBreathPermitivity(ilay);
+        }
+
+        // Freq stuff
+        double finit  = Dipole->GetFrequency(0); //(1000);   // Starting freq
+        double flimit = Dipole->GetFrequency(0); //(1000);   // max freq
+        double dlimit = Dipole->GetFrequency(0); //(1000);   // difusion limit
+        double lfinc(1);       // no. freq per decade
+
+        // tx location jtype != 4
+        double txx = Dipole->GetLocation(0); // (0.);
+        double txy = Dipole->GetLocation(1); // (0.);
+        double txz = Dipole->GetLocation(2); // (0.);
+
+        // rx position
+        // TODO, fix Fortran program to not waste this memory
+        // maybe
+        const int MAXREC = 15;
+        double *rxx = new double [MAXREC];
+        double *rxy = new double [MAXREC];
+        double *rxz = new double [MAXREC];
+
+        std::complex<double> *ex = new std::complex<double>[MAXREC];
+        std::complex<double> *ey = new std::complex<double>[MAXREC];
+        std::complex<double> *ez = new std::complex<double>[MAXREC];
+
+        std::complex<double> *hx = new std::complex<double>[MAXREC];
+        std::complex<double> *hy = new std::complex<double>[MAXREC];
+        std::complex<double> *hz = new std::complex<double>[MAXREC];
+
+        int nres2 = MAXREC;
+        int ii=0;
+
+        for (ii=0; ii<nres-MAXREC; ii+=MAXREC) {
+
+            for (int ir=0; ir<MAXREC; ++ir) {
+                //Vector3r pos = Receivers->GetLocation(ii+ir);
+                rxx[ir] = Receivers->GetLocation(ii+ir)[0];
+                rxy[ir] = Receivers->GetLocation(ii+ir)[1];
+                rxz[ir] = Receivers->GetLocation(ii+ir)[2];
+            }
+
+            em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
+                  jgamma, acc, dep, sig, susl, sush, sustau, susalp,
+                  eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
+                  lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
+                  hx, hy, hz);
+
+            // Scale By Moment
+            for (int ir=0; ir<MAXREC; ++ir) {
+
+                ex[ir] *= Dipole->GetMoment();
+                ey[ir] *= Dipole->GetMoment();
+                ez[ir] *= Dipole->GetMoment();
+
+                hx[ir] *= Dipole->GetMoment();
+                hy[ir] *= Dipole->GetMoment();
+                hz[ir] *= Dipole->GetMoment();
+
+                // Append values instead of setting them
+                this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
+                                                        (Complex)(ey[ir]),
+                                                        (Complex)(ez[ir]) );
+                this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
+                                                        (Complex)(hy[ir]),
+                                                        (Complex)(hz[ir]) );
+            }
+        }
+
+        //ii += MAXREC;
+        nres2 = 0;
+        // Perform last positions
+        for (int ir=0; ir<nres-ii; ++ir) {
+            rxx[ir] = Receivers->GetLocation(ii+ir)[0];
+            rxy[ir] = Receivers->GetLocation(ii+ir)[1];
+            rxz[ir] = Receivers->GetLocation(ii+ir)[2];
+            ++nres2;
+        }
+
+        em1dcall_(itype, ipol, nlay, nfreq, nfield, nres2, jtype,
+              jgamma, acc, dep, sig, susl, sush, sustau, susalp,
+              eprl, eprh, eprtau, epralp, finit, flimit, dlimit,
+              lfinc, txx, txy, txz, rxx, rxy, rxz, ex, ey, ez,
+              hx, hy, hz);
+
+        // Scale By Moment
+        for (int ir=0; ir<nres-ii; ++ir) {
+
+            ex[ir] *= Dipole->GetMoment();
+            ey[ir] *= Dipole->GetMoment();
+            ez[ir] *= Dipole->GetMoment();
+
+            hx[ir] *= Dipole->GetMoment();
+            hy[ir] *= Dipole->GetMoment();
+            hz[ir] *= Dipole->GetMoment();
+
+            // Append values instead of setting them
+            this->Receivers->AppendEfield(0, ii+ir, (Complex)(ex[ir]),
+                                                    (Complex)(ey[ir]),
+                                                    (Complex)(ez[ir]) );
+            this->Receivers->AppendHfield(0, ii+ir, (Complex)(hx[ir]),
+                                                    (Complex)(hy[ir]),
+                                                    (Complex)(hz[ir]) );
+
+        }
+
+        delete [] sig;
+        delete [] dep;
+
+        //delete [] sus;
+        //delete [] epr;
+
+        delete [] susl;
+        delete [] sush;
+        delete [] susalp;
+        delete [] sustau;
+
+        delete [] eprl;
+        delete [] eprh;
+        delete [] epralp;
+        delete [] eprtau;
+
+        delete [] rxx;
+        delete [] rxy;
+        delete [] rxz;
+
+        delete [] ex;
+        delete [] ey;
+        delete [] ez;
+
+        delete [] hx;
+        delete [] hy;
+        delete [] hz;
+
+    }
+#endif
+
+
+    void EMEarth1D::SolveSingleTxRxPair (const int &irec, HankelTransform *Hankel, const Real &wavef, const int &ifreq,
+                   DipoleSource *tDipole) {
+        ++icalcinner;
+
+        Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
+
+        tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
+        Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
+        tDipole->UpdateFields( ifreq,  Hankel, wavef );
+    }
+
+    void EMEarth1D::SolveLaggedTxRxPair(const int &irec, FHTAnderson801* Hankel,
+                    const Real &wavef, const int &ifreq, PolygonalWireAntenna* antenna) {
+
+        antenna->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
+
+        // Determine the min and max arguments
+        Real rhomin = 1e9;
+        Real rhomax = 1e-9;
+        for (int idip=0; idip<antenna->GetNumberOfDipoles(); ++idip) {
+            auto tDipole = antenna->GetDipoleSource(idip);
+            Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
+            rhomin = std::min(rhomin, rho);
+            rhomax = std::max(rhomax, rho);
+        }
+        //std::cout << "rhomin\t" << rhomin << "\trhomax" << rhomax << std::endl;
+
+        // Determine number of lagged convolutions to do
+        // TODO, can Hankel2 adjust the lagg spacing safely?
+        int nlag = 1; // We need an extra for some reason for stability
+        Real lrho ( 1.01* rhomax );
+        while ( lrho > rhomin ) {
+            nlag += 1;
+            lrho *= Hankel->GetABSER();
+        }
+
+        //int nlag = rhomin
+        auto tDipole = antenna->GetDipoleSource(0);
+        tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
+
+        // Instead we should pass the antenna into this so that Hankel hass all the rho arguments...
+        Hankel->ComputeLaggedRelated( 1.01* rhomax, nlag, tDipole->GetKernelManager() );
+
+        //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);
+            tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
+            // Pass Hankel2 a message here so it knows which one to return in Zgauss!
+            Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
+            //std::cout << " in Lagged " <<  rho << "\t" << rhomin << "\t" << rhomax << std::endl;
+            Hankel->SetLaggedArg( rho );
+            //std::cout << "out Lagged" << std::endl;
+            tDipole->UpdateFields( ifreq,  Hankel, wavef );
+        }
+        //std::cout << "Spline\n";
+        //std::cout << Receivers->GetHfield(0, irec) << std::endl;
+    }
+
+    //////////////////////////////////////////////////////////
+    // Thread safe OO Reimplimentation of KiHand's
+    // EM1DNEW.for programme
+    void EMEarth1D::MakeCalc3() {
+
+        if ( Dipole == nullptr ) throw NullDipoleSource();
+
+        if (Earth == nullptr) throw NullEarth();
+
+        if (Receivers == nullptr) throw NullReceivers();
+
+        #ifdef LEMMAUSEOMP
+        #pragma omp parallel
+        #endif
+        { // OpenMP Parallel Block
+
+            #ifdef LEMMAUSEOMP
+            int tid = omp_get_thread_num();
+            int nthreads = omp_get_num_threads();
+            #else
+            int tid=0;
+            int nthreads=1;
+            #endif
+
+            auto tDipole = Dipole->Clone();
+
+            std::shared_ptr<HankelTransform> Hankel;
+            switch (HankelType) {
+                case ANDERSON801:
+                    Hankel = FHTAnderson801::NewSP();
+                    break;
+                case CHAVE:
+                    Hankel = GQChave::NewSP();
+                    break;
+                case FHTKEY201:
+                    Hankel = FHTKey201::NewSP();
+                    break;
+                case FHTKEY101:
+                    Hankel = FHTKey101::NewSP();
+                    break;
+                case FHTKEY51:
+                    Hankel = FHTKey51::NewSP();
+                    break;
+                case QWEKEY:
+                    Hankel = QWEKey::NewSP();
+                    break;
+                default:
+                    std::cerr << "Hankel transform cannot be created\n";
+                    exit(EXIT_FAILURE);
+            }
+
+            if ( tDipole->GetNumberOfFrequencies() < Receivers->GetNumberOfPoints() ) {
+                for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies(); ++ifreq) {
+                    // Propogation constant in free space being input to Hankel
+                    Real wavef   = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
+                    for (int irec=tid; irec<Receivers->GetNumberOfPoints(); irec+=nthreads) {
+                        SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
+                    }
+                }
+            } else {
+                for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
+                    for (int ifreq=tid; ifreq<tDipole->GetNumberOfFrequencies(); ifreq+=nthreads) {
+                        // Propogation constant in free space being input to Hankel
+                        Real wavef   = tDipole->GetAngularFrequency(ifreq) * std::sqrt(MU0*EPSILON0);
+                        SolveSingleTxRxPair(irec, Hankel.get(), wavef, ifreq, tDipole.get());
+                    }
+                }
+            }
+        } // OpenMP Parallel Block
+    }
+
+    NullReceivers::NullReceivers() :
+        runtime_error("nullptr RECEIVERS") {}
+
+    NullAntenna::NullAntenna() :
+        runtime_error("nullptr ANTENNA") {}
+
+    NullInstrument::NullInstrument(LemmaObject* ptr) :
+        runtime_error("nullptr INSTRUMENT") {
+        std::cout << "Thrown by instance of "
+                  << ptr->GetName() << std::endl;
+    }
+
+    DipoleSourceSpecifiedForWireAntennaCalc::
+        DipoleSourceSpecifiedForWireAntennaCalc() :
+        runtime_error("DIPOLE SOURCE SPECIFIED FOR WIRE ANTENNA CALC"){}
+
+} // end of Lemma Namespace

vim/c.vim

 
 " Lemma types
 highlight leType ctermfg=Yellow guifg=Yellow
-syn keyword leType HankelTransform FHTAnderson801 FHTKey201 FHTKey101 FHTKey51 GQChave QWEKey KernelEM1D KernelEM1DSpec KernelEM1DBase KernelEM1DManager DipoleSource EarthModel LayeredEarth LayeredEarthEM TEMSurvey TEMSurveyLine TEMSurveyLineRecord TEMInductiveReceiver WireAntenna PolygonalWireAntenna TEMTransmitter TEMReceiver Instrument InstrumentTem LemmaObject FieldPoints DCIPElectrode TEMSurveyData TEMSurveyLineData TEMSurveyLineRecordData  ASCIIParser CubicSplineInterpolator RectilinearGrid GridReader RectilinearGridReader RectilinearGridVTKExporter Filter WindowFilter DEMParticle DEMGrain Data DataReader 
+syn keyword leType HankelTransform FHTAnderson801 FHTKey201 FHTKey101 FHTKey51 GQChave QWEKey KernelEM1D KernelEM1DSpec KernelEM1DBase KernelEM1DManager DipoleSource EarthModel LayeredEarth LayeredEarthEM TEMSurvey TEMSurveyLine TEMSurveyLineRecord TEMInductiveReceiver WireAntenna PolygonalWireAntenna TEMTransmitter TEMReceiver Instrument InstrumentTem LemmaObject FieldPoints DCIPElectrode TEMSurveyData TEMSurveyLineData TEMSurveyLineRecordData  ASCIIParser CubicSplineInterpolator RectilinearGrid GridReader RectilinearGridReader RectilinearGridVTKExporter Filter WindowFilter DEMParticle DEMGrain Data DataReader EMEarth1D 
 
 " Deprecated Lemma Types
 highlight dleType ctermfg=Blue guifg=Blue