Merge branch 'master' of https://lemma.codes/Lemma

CMake/FindUmfpack.cmake

@@ -47,7 +47,7 @@ if(UMFPACK_LIBRARIES)
 endif(UMFPACK_LIBRARIES)
 
 include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(UMFPACK DEFAULT_MSG
+find_package_handle_standard_args(Umfpack DEFAULT_MSG
                                   UMFPACK_INCLUDES UMFPACK_LIBRARIES)
 
 mark_as_advanced(UMFPACK_INCLUDES UMFPACK_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY CHOLMOD_LIBRARY SUITESPARSE_LIBRARY)

CMake/SuperBuild.cmake

@@ -8,7 +8,7 @@ else()
         ExternalProject_Add(EIGEN
 	    GIT_REPOSITORY "https://gitlab.com/libeigen/eigen.git"
         UPDATE_COMMAND "" 
-	    GIT_TAG "3.3.7" #"default"
+	    GIT_TAG "3.3.7" 
    	    PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/eigen
    	    CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX}
 		#CONFIGURE_COMMAND ""
@@ -24,12 +24,12 @@ else()
 
     ExternalProject_Add(YAML_CPP
         GIT_REPOSITORY  "https://github.com/jbeder/yaml-cpp.git" 
-        GIT_TAG  "yaml-cpp-0.6.2"  # "master" 
+        GIT_TAG  "yaml-cpp-0.6.3"  # "master" 
         UPDATE_COMMAND ""
         PATCH_COMMAND ""
         PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/yaml-cpp
         CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} 
-                   -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS}
+                   -DYAML_BUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} 
                    -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} 
     	           -DYAML_CPP_BUILD_TESTS=OFF
                    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}
@@ -41,13 +41,13 @@ else()
     )
 endif()
 
-if ( LEMMA_VTK8_SUPPORT )
-    if ( NOT VTK_FOUND OR  VTK_VERSION VERSION_GREATER "8.2.0" )
+if ( LEMMA_VTK9_SUPPORT )
+    if ( NOT VTK_FOUND OR  VTK_VERSION VERSION_LESS "9.0.0" )
         message( STATUS "VTK > 8.20.0 was found! Version found: " ${VTK_VERSION}, ${VTK_USE_FILE} )
         message( STATUS "External build of VTK 8 has been added, this may take some time to build." )
-        ExternalProject_Add(VTK8
+        ExternalProject_Add(VTK9
         GIT_REPOSITORY "https://gitlab.kitware.com/vtk/vtk.git"
-        GIT_TAG  "v8.2.0"
+        GIT_TAG  "v9.0.1"
         PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/vtk8
         CMAKE_ARGS
             -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX}
@@ -77,9 +77,11 @@ if (LEMMA_PYTHON3_BINDINGS)
         message( STATUS "pybind11 was found" )
     else()
         message( STATUS "pybind11 was NOT found, please build or remove LEMMA_PYTHON3_BINDINGS" )
-	    ExternalProject_Add(pybind11
+        #find_package(Python COMPONENTS Interpreter Development REQUIRED)
+	#find_package(PythonLibs 3.0 REQUIRED)
+        ExternalProject_Add(pybind11
 		    GIT_REPOSITORY "https://github.com/pybind/pybind11.git"
-		    GIT_TAG "v2.3.0" #"master"
+		    GIT_TAG "v2.5.0" # "master" #"v2.4.3" #"master"
 		    UPDATE_COMMAND ""
 		    PATCH_COMMAND ""
     	    PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/pybind11

CMakeLists.txt

@@ -13,8 +13,8 @@ SET_PROPERTY(GLOBAL PROPERTY TARGET_SUPPORTS_SHARED_LIBS TRUE)
 ####################################################################################################
 # Lemma versioning, set Major, minor, and patch here                                               #
 set(LEMMA_VERSION_MAJOR "0")                                                                       #
-set(LEMMA_VERSION_MINOR "3")                                                                       #
-set(LEMMA_VERSION_PATCH "3")                                                                       #
+set(LEMMA_VERSION_MINOR "4")                                                                       #
+set(LEMMA_VERSION_PATCH "0")                                                                       #
 set(LEMMA_VERSION "\"${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}.${LEMMA_VERSION_PATCH}\"")      #
 set(LEMMA_VERSION_NOQUOTES "${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}.${LEMMA_VERSION_PATCH}") #
 ####################################################################################################
@@ -47,6 +47,7 @@ option ( LEMMA_BUILD_EXAMPLES       "Compile example Lemma applications" OFF )
 option ( LEMMA_USE_OPENMP           "Use OpenMP in Lemma" OFF )
 option ( LEMMA_BUILD_DOCUMENTATION  "Build Doxygen man pages" OFF )
 option ( LEMMA_VTK8_SUPPORT         "VTK 8.x library for visualisation and grids" OFF )
+option ( LEMMA_VTK9_SUPPORT         "VTK 9.x library for visualisation and grids" OFF )
 option ( LEMMA_PYTHON3_BINDINGS     "Compile Python 3 bindings" OFF )
 
 # We end up using this for all builds, TODO remove this variable but follow same path
@@ -99,11 +100,32 @@ if (CMAKE_CROSSCOMPILING)
         QUIET
         ) 
     endif()
+    if (LEMMA_VTK9_SUPPORT)                                                            # Visualisation 
+        find_package (VTK  9.0.1      
+        COMPONENTS 
+        CommonCore 
+        RenderingCore 
+        FiltersCore 
+        FiltersSources 
+        CommonDataModel 
+        FiltersHyperTree 
+        IOXML 
+        IOImage 
+        IOLegacy 
+        IOGeometry 
+        InteractionStyle 
+        RenderingAnnotation 
+        FiltersHybrid 
+        FiltersModeling 
+        RenderingVolumeOpenGL2
+        QUIET
+        ) 
+    endif()
 else()
     find_package (Eigen3   3.3 QUIET )  # Matrix/Vector & Math
     find_package (yaml-cpp 0.6 QUIET )  # Serialisation of classes 
-    if (LEMMA_VTK8_SUPPORT)
-        find_package (VTK  8.2.0 COMPONENTS
+    if (LEMMA_VTK9_SUPPORT)
+        find_package (VTK  9.0.1 COMPONENTS
         vtkCommonCore 
         vtkRenderingCore 
         vtkFiltersCore 
@@ -125,14 +147,19 @@ endif()
 
 INCLUDE_DIRECTORIES(${YAML_CPP_INCLUDE_DIR})
 
-if (VTK_VERSION VERSION_GREATER "8.2.0")
-    message( STATUS "${VTK_VERSION} is compatible: ${VTK_VERSION_COMPATIBLE} exact? ${VTK_VERSION_EXACT}" )
-    set (VTK_FOUND False) 
-endif()
+#if (VTK_VERSION VERSION_GREATER "8.2.0")
+#    message( STATUS "${VTK_VERSION} is compatible: ${VTK_VERSION_COMPATIBLE} exact? ${VTK_VERSION_EXACT}" )
+#    set (VTK_FOUND False) 
+#endif()
 
 if (VTK_FOUND)
     set(volumeRenderer volumerenderer.cxx)
 	add_definitions(-DLEMMAUSEVTK) 
+	if (LEMMA_VTK8_SUPPORT)
+    add_definitions(-DLEMMA_VTK8_SUPPORT)
+    else (LEMMA_VTK9_SUPPORT) 
+    add_definitions(-DLEMMA_VTK9_SUPPORT)
+    endif()
 endif()
 
 if (LEMMA_BUILD_DOCUMENTATION)
@@ -160,7 +187,7 @@ message ( STATUS "VTK Found? " ${VTK_FOUND} )
 if ( NOT Eigen3_FOUND OR 
      NOT yaml-cpp_FOUND OR 
      (LEMMA_PYTHON3_BINDINGS AND NOT pybind11_FOUND) OR 
-     (LEMMA_VTK8_SUPPORT AND NOT VTK_FOUND) OR
+     (LEMMA_VTK9_SUPPORT AND NOT VTK_FOUND) OR
      (LEMMA_ENABLE_TESTING AND NOT CxxTest_FOUND) )
     message ( STATUS "Missing hard dependencies have been found, these will be downloaded any compiled." )
     message ( STATUS "This necessitates a two step build." )
@@ -247,13 +274,16 @@ if (LEMMA_USE_BOOST)
 endif()
 
 if (LEMMA_PYTHON3_BINDINGS)
-    find_package(PythonLibs 3.0 REQUIRED)
+    #find_package(PythonLibs 3.0 REQUIRED)
+    #find_package(pythoninterp REQUIRED)
     find_package(pybind11 REQUIRED)
     #    find_package(Boost 1.64.0   REQUIRED COMPONENTS system python3 numpy3)
     # Is this OK for portability of non-python libraries
     #include_directories(${PYTHON_INCLUDE_DIRS})
     #include_directories(${Boost_INCLUDE_DIRS})
     install ( FILES python/setup.py DESTINATION ${CMAKE_INSTALL_PREFIX}/pyLemma/ ) 
+    install ( FILES python/long.md DESTINATION ${CMAKE_INSTALL_PREFIX}/pyLemma/ ) 
+    install ( FILES python/publish.sh DESTINATION ${CMAKE_INSTALL_PREFIX}/pyLemma/ ) 
     install ( DIRECTORY python/pyLemma DESTINATION ${CMAKE_INSTALL_PREFIX}/pyLemma/ ) 
 endif()
 
@@ -353,9 +383,10 @@ endif()
 ########################################################################
 # add a target to generate API documentation with Doxyfile.in 
 # ALL make documentation build by default if possible
-find_package(Doxygen)
+if (LEMMA_BUILD_DOCUMENTATION)
+    find_package(Doxygen)
 	if(DOXYGEN_FOUND)
-	if (LEMMA_BUILD_DOCUMENTATION)
+    message( STATUS "LEMMA_BUILD_DOCUMENTATION must be positive" )
 # Custom header and footer option, enable in Doxygen 
 #	file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/Documentation/header.html
 #        DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/Documentation/header.html
@@ -370,7 +401,7 @@ find_package(Doxygen)
 			WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
 			COMMENT "Generating API documentation with Doxygen" VERBATIM
 		)
-	endif (LEMMA_BUILD_DOCUMENTATION)
 	endif(DOXYGEN_FOUND)
+endif (LEMMA_BUILD_DOCUMENTATION)
 
 # vim: set tabstop=4 shiftwidth=4 expandtab: 

Modules/FDEM1D/CMakeLists.txt

@@ -10,6 +10,11 @@ set_target_properties(fdem1d PROPERTIES
     CXX_STANDARD_REQUIRED ON
 )
 
+if ( LEMMA_VTK9_SUPPORT ) 
+	target_link_libraries(fdem1d ${visibility}VTK::RenderingCore VTK::CommonDataModel)
+	vtk_module_autoinit(TARGETS fdem1d MODULES VTK::RenderingCore)
+endif()
+
 # Linking
 target_link_libraries(fdem1d "lemmacore")
 target_link_libraries(fdem1d ${YAML_CPP_LIBRARIES}) 

Modules/FDEM1D/include/KernelEM1DSpec.h

@@ -1338,7 +1338,7 @@ namespace Lemma {
         static bool called = false;
         if (!called) {
             std::cout << "WARNING\n";
-            std::cout << "Unspecialised PotentialBelowSourceLayer <" << Mode << " " << Ikernel << " "
+            std::cout << "Unspecialised RelPotentialBelowSourceLayer <" << Mode << " " << Ikernel << " "
                   << Isource << " " << Irecv << ">...this function will be slow\n\n";
             called = true;
         }

Modules/FDEM1D/include/PolygonalWireAntenna.h

@@ -118,6 +118,12 @@ namespace Lemma {
             /// Maximum dipole moment allowed
             Real maxDipoleMoment;
 
+            /// Adds grounding point to non-closed wire loop
+            void AddGroundingPoint(
+                            const Vector3r &cp,
+                            const Vector3r &dir,
+                            std::vector< std::shared_ptr<DipoleSource> > &Dipoles );
+
             /// appends
             void PushXYZDipoles( const Vector3r &step, const Vector3r &cp,
                             const Vector3r &dir,

Modules/FDEM1D/python/pyFDEM1D.cpp

@@ -249,14 +249,14 @@ PYBIND11_MODULE(FDEM1D, m) {
         // modifiers
         .def("AttachWireAntenna", &Lemma::EMEarth1D::AttachWireAntenna,
             "Sets the wire antenna to use for calculations")
-        .def("AttachDipoleSOurce", &Lemma::EMEarth1D::AttachDipoleSource,
+        .def("AttachDipoleSource", &Lemma::EMEarth1D::AttachDipoleSource,
             "Sets a DipoleSource to use for calculations")
         .def("AttachFieldPoints", &Lemma::EMEarth1D::AttachFieldPoints,
             "Sets the FieldPoints to use for calculations")
         .def("AttachLayeredEarthEM", &Lemma::EMEarth1D::AttachLayeredEarthEM,
             "Sets the LayeredEarthEM to use for calculations")
 
-        .def("SetFieldToCalculate", &Lemma::EMEarth1D::SetFieldsToCalculate,
+        .def("SetFieldsToCalculate", &Lemma::EMEarth1D::SetFieldsToCalculate,
             "Sets which fields to calculate")
         .def("SetHankelTransformMethod", &Lemma::EMEarth1D::SetHankelTransformMethod,
             "Sets which Hankel transform to use")

Modules/FDEM1D/src/DipoleSource.cpp

@@ -183,6 +183,9 @@ namespace Lemma {
             case (UNGROUNDEDELECTRICDIPOLE):
                 this->Type = stype;
                 break;
+            case (GROUNDINGPOINT):
+                this->Type = stype;
+                break;
             case (MAGNETICDIPOLE):
                 this->Type = stype;
                 break;
@@ -485,6 +488,172 @@ namespace Lemma {
                     }
                 break;
 
+            case (GROUNDINGPOINT):
+
+                if (std::abs(Pol[2]) > 0) { // z dipole
+
+                    switch(FieldsToCalculate) {
+
+                        case E:
+                            if (lays == 0 && layr == 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INAIR>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INGROUND>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INAIR>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
+                            } else {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INGROUND>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
+                            }
+                            break;
+
+                        case H:
+                            if (lays == 0 && layr == 0) {
+                                ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
+                            } else {
+                                ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
+                            }
+                            break;
+
+
+                        case BOTH:
+                            if ( lays == 0 && layr == 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INAIR>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INAIR>( );
+                                ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INAIR, INGROUND>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INAIR, INGROUND>( );
+                                ik[12] = KernelManager->AddKernel<TM, 12, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INAIR>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INAIR>( );
+                                ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INAIR>( );
+                            } else {
+                                ik[10] = KernelManager->AddKernel<TM, 10, INGROUND, INGROUND>( );
+                                ik[11] = KernelManager->AddKernel<TM, 11, INGROUND, INGROUND>( );
+                                ik[12] = KernelManager->AddKernel<TM, 12, INGROUND, INGROUND>( );
+                            }
+                        }
+                }
+                if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y grounded HED dipole
+
+                    switch(FieldsToCalculate) {
+
+                        case E:
+                            if ( lays == 0 && layr == 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INAIR>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INAIR>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INAIR>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INGROUND>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INGROUND>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INGROUND>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INAIR>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INAIR>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INAIR>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
+                            } else {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INGROUND>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INGROUND>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INGROUND>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
+                            }
+                            break;
+
+                        case H:
+                            if (lays == 0 && layr == 0) {
+                                ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
+                            } else {
+                                ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
+                            }
+                            break;
+
+                        case BOTH:
+                            if (lays == 0 && layr == 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INAIR>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INAIR>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INAIR>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INAIR>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INAIR>( );
+                                ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INAIR>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INAIR>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INAIR>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INAIR>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INAIR>( );
+                            } else if (lays == 0 && layr > 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INAIR, INGROUND>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INAIR, INGROUND>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INAIR, INGROUND>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INAIR, INGROUND>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INAIR, INGROUND>( );
+                                ik[5] = KernelManager->AddKernel<TM, 5, INAIR, INGROUND>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INAIR, INGROUND>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INAIR, INGROUND>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INAIR, INGROUND>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INAIR, INGROUND>( );
+                            } else if (lays > 0 && layr == 0) {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INAIR>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INAIR>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INAIR>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INAIR>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INAIR>( );
+                                ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INAIR>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INAIR>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INAIR>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INAIR>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INAIR>( );
+                            } else {
+                                ik[0] = KernelManager->AddKernel<TM, 0, INGROUND, INGROUND>( );
+                                ik[1] = KernelManager->AddKernel<TM, 1, INGROUND, INGROUND>( );
+                                ik[4] = KernelManager->AddKernel<TM, 4, INGROUND, INGROUND>( );
+                                //ik[2] = KernelManager->AddKernel<TE, 2, INGROUND, INGROUND>( );
+                                //ik[3] = KernelManager->AddKernel<TE, 3, INGROUND, INGROUND>( );
+                                ik[5] = KernelManager->AddKernel<TM, 5, INGROUND, INGROUND>( );
+                                ik[6] = KernelManager->AddKernel<TM, 6, INGROUND, INGROUND>( );
+                                //ik[7] = KernelManager->AddKernel<TE, 7, INGROUND, INGROUND>( );
+                                //ik[8] = KernelManager->AddKernel<TE, 8, INGROUND, INGROUND>( );
+                                //ik[9] = KernelManager->AddKernel<TE, 9, INGROUND, INGROUND>( );
+                            }
+                            break;
+                        }
+                    }
+                break;
+
             case (UNGROUNDEDELECTRICDIPOLE):
 
                 if (std::abs(Pol[2]) > 0) { // z dipole
@@ -786,8 +955,6 @@ namespace Lemma {
                 exit(EXIT_FAILURE);
 
         }
-
-
     }
 
     void DipoleSource::UpdateFields( const int& ifreq, HankelTransform* Hankel, const Real& wavef) {
@@ -874,6 +1041,7 @@ namespace Lemma {
                                     Pol[0]*Moment*QPI*sp*f(9) );
                             }
                             break;
+
                         case BOTH:
                             f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
                             f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
@@ -913,6 +1081,128 @@ namespace Lemma {
                 }
                 break; // GROUNDEDELECTRICDIPOLE
 
+            case (GROUNDINGPOINT):
+                if (std::abs(Pol[2]) > 0) { // z dipole
+                    switch(FieldsToCalculate) {
+                        case E:
+                            f(10) = Hankel->Zgauss(10, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10])) / KernelManager->GetRAWKernel(ik[10])->GetYm();
+                            f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
+                            this->Receivers->AppendEfield(ifreq, irec,
+                                -Pol[2]*QPI*cp*f(10)*Moment,
+                                -Pol[2]*QPI*sp*f(10)*Moment,
+                                 Pol[2]*QPI*f(11)*Moment);
+                            break;
+
+                        case H:
+                            f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
+                            this->Receivers->AppendHfield(ifreq, irec,
+                                -Pol[2]*QPI*sp*f(12)*Moment,
+                                 Pol[2]*QPI*cp*f(12)*Moment,
+                                 0. );
+                            break;
+
+                        case BOTH:
+                            f(10) = Hankel->Zgauss(10, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[10])) / KernelManager->GetRAWKernel(ik[10])->GetYm();
+                            f(11) = Hankel->Zgauss(11, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[11])) / KernelManager->GetRAWKernel(ik[11])->GetYm();
+                            this->Receivers->AppendEfield(ifreq, irec,
+                                    -Pol[2]*QPI*cp*f(10)*Moment,
+                                    -Pol[2]*QPI*sp*f(10)*Moment,
+                                     Pol[2]*QPI*f(11)*Moment   );
+
+                            f(12) = Hankel->Zgauss(12, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[12]));
+                            this->Receivers->AppendHfield(ifreq, irec,
+                                    -Pol[2]*QPI*sp*f(12)*Moment,
+                                     Pol[2]*QPI*cp*f(12)*Moment,
+                                     0. );
+                    } // Fields to calculate Z polarity Electric dipole
+                }
+                if (std::abs(Pol[1]) > 0 || std::abs(Pol[0]) > 0) { // x or y dipole
+                    switch(FieldsToCalculate) {
+                        case E:
+                            f(2) = 0;//Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
+                            f(3) = 0;//Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
+                            f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
+                            f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
+                            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,
+                                    0,0,0);
+                                    //Pol[1]*QPI*Moment*scp*(f(0)+f(2)),
+                                    //Pol[1]*QPI*Moment*((sps*f(0)+c2p*f(1)/rho)),
+                                    //Pol[1]*QPI*Moment*(f(0)+f(2)),
+                                    //Pol[1]*QPI*Moment*((f(0)+f(1)/rho)),
+                                    //Pol[1]*QPI*sp*f(4)*Moment);
+                                    // std dipole
+                                    //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);
+                            }
+                            if (std::abs(Pol[0]) > 0) {
+                                this->Receivers->AppendEfield(ifreq, irec,
+                                    Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)),          //-(sps*f(2)+c2p*f(3)/rho)),
+                                    Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)),   //+(f(2)-(Real)(2.)*f(3)/rho)),
+                                    Pol[0]*Moment*QPI*cp*f(4) );
+                            }
+                            break;
+                        case H:
+                            f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
+                            f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
+                            f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
+                            f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
+                            f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
+                            if (std::abs(Pol[1]) > 0) {
+                                this->Receivers->AppendHfield(ifreq, irec,
+                                        Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
+                                        Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
+                                        -Pol[1]*QPI*cp*f(9)*Moment );
+                            }
+                            if (std::abs(Pol[0]) > 0) {
+                                this->Receivers->AppendHfield(ifreq, irec,
+                                    Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
+                                    Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
+                                    Pol[0]*Moment*QPI*sp*f(9) );
+                            }
+                            break;
+
+                        case BOTH:
+                            f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
+                            f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
+                            f(4) = Hankel->Zgauss(4, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[4])) / KernelManager->GetRAWKernel(ik[4])->GetYm();
+                            f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
+                            f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
+                            f(5) = Hankel->Zgauss(5, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[5]));
+                            f(6) = Hankel->Zgauss(6, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[6]));
+                            f(7) = Hankel->Zgauss(7, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[7]))*KernelManager->GetRAWKernel(ik[7])->GetZs()/KernelManager->GetRAWKernel(ik[7])->GetZm();
+                            f(8) = Hankel->Zgauss(8, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[8]))*KernelManager->GetRAWKernel(ik[8])->GetZs()/KernelManager->GetRAWKernel(ik[8])->GetZm();
+                            f(9) = Hankel->Zgauss(9, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[9]))*KernelManager->GetRAWKernel(ik[9])->GetZs()/KernelManager->GetRAWKernel(ik[9])->GetZm();
+
+                            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);
+
+                                this->Receivers->AppendHfield(ifreq, irec,
+                                        Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
+                                        Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
+                                       -Pol[1]*QPI*cp*f(9)*Moment );
+                            }
+                            if (std::abs(Pol[0]) > 0) {
+                                this->Receivers->AppendEfield(ifreq, irec,
+                                    Pol[0]*Moment*QPI*((cps*f(0)-c2p*f(1)/rho)-(sps*f(2)+c2p*f(3)/rho)),
+                                    Pol[0]*Moment*QPI*scp*((f(0)-(Real)(2.)*f(1)/rho)+(f(2)-(Real)(2.)*f(3)/rho)),
+                                    Pol[0]*Moment*QPI*cp*f(4) );
+
+                                this->Receivers->AppendHfield(ifreq, irec,
+                                    Pol[0]*Moment*QPI*scp*(f(5)-(Real)(2.)*f(6)/rho+f(7)-(Real)(2.)*f(8)/rho),
+                                    Pol[0]*Moment*QPI*(-cps*f(5)+c2p*f(6)/rho+sps*f(7)+c2p*f(8)/rho),
+                                    Pol[0]*Moment*QPI*sp*f(9) );
+                            }
+                            break;
+                    }
+                }
+                break; // GROUNDINGPOINT
+
 
             case UNGROUNDEDELECTRICDIPOLE:
 
@@ -955,11 +1245,16 @@ namespace Lemma {
 
                     switch(FieldsToCalculate) {
                         case E:
-                            f(0) = 0;
-                            f(1) = 0;
+                            //f(0) = 0;
+                            //f(1) = 0;
+                            //f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
+                            //f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
+                            //f(4) = 0;
                             f(2) = Hankel->Zgauss(2, TE, 0, rho, wavef, KernelManager->GetRAWKernel(ik[2])) * KernelManager->GetRAWKernel(ik[2])->GetZs();
                             f(3) = Hankel->Zgauss(3, TE, 1, rho, wavef, KernelManager->GetRAWKernel(ik[3])) * KernelManager->GetRAWKernel(ik[3])->GetZs();
-                            f(4) = 0;
+                            f(0) = Hankel->Zgauss(0, TM, 0, rho, wavef, KernelManager->GetRAWKernel(ik[0])) / KernelManager->GetRAWKernel(ik[0])->GetYm();
+                            f(1) = Hankel->Zgauss(1, TM, 1, rho, wavef, KernelManager->GetRAWKernel(ik[1])) / KernelManager->GetRAWKernel(ik[1])->GetYm();
+                            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,
@@ -984,7 +1279,7 @@ namespace Lemma {
                                 this->Receivers->AppendHfield(ifreq, irec,
                                         Pol[1]*QPI*(sps*f(5)+c2p*f(6)/rho-cps*f(7)+c2p*f(8)/rho)*Moment,
                                         Pol[1]*QPI*scp*(-f(5)+(Real)(2.)*f(6)/rho-f(7)+(Real)(2.)*f(8)/rho)*Moment,
-                                        -Pol[1]*QPI*cp*f(9)*Moment );
+                                       -Pol[1]*QPI*cp*f(9)*Moment );
                             }
                             if (std::abs(Pol[0]) > 0) {
                                 this->Receivers->AppendHfield(ifreq, irec,
@@ -1036,7 +1331,6 @@ namespace Lemma {
                 break; // UNGROUNDEDELECTRICDIPOLE
 
             case MAGNETICDIPOLE:
-
                 //Hankel->ComputeRelated(rho, KernelManager);
                 if (std::abs(Pol[2]) > 0) { // z dipole
                     switch(FieldsToCalculate) {

Modules/FDEM1D/src/EMEarth1D.cpp

@@ -72,8 +72,7 @@ namespace Lemma {
     // TODO init large arrays here.
     EMEarth1D::EMEarth1D( const ctor_key& key ) : LemmaObject( key ),
             Dipole(nullptr), Earth(nullptr), Receivers(nullptr), Antenna(nullptr),
-            FieldsToCalculate(BOTH), HankelType(ANDERSON801), icalcinner(0), icalc(0)
-        {
+            FieldsToCalculate(BOTH), HankelType(ANDERSON801), icalcinner(0), icalc(0) {
     }
 
     EMEarth1D::~EMEarth1D() {
@@ -107,21 +106,36 @@ namespace Lemma {
     // ====================  ACCESS        ===================================
     void EMEarth1D::AttachDipoleSource( std::shared_ptr<DipoleSource> dipoleptr) {
         Dipole = dipoleptr;
-    }
-
-    void EMEarth1D::AttachLayeredEarthEM( std::shared_ptr<LayeredEarthEM> earthptr) {
-        Earth = earthptr;
+        if (Receivers != nullptr) {
+            // Check to make sure Receivers are set up for all calculations
+            switch(FieldsToCalculate) {
+                case E:
+                    if (Receivers->NumberOfBinsE != Dipole->GetNumberOfFrequencies())
+                        Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
+                    break;
+                case H:
+                    if (Receivers->NumberOfBinsH != Dipole->GetNumberOfFrequencies())
+                        Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
+                    break;
+                case BOTH:
+                    if (Receivers->NumberOfBinsH != Dipole->GetNumberOfFrequencies())
+                        Receivers->SetNumberOfBinsH(Dipole->GetNumberOfFrequencies());
+                    if (Receivers->NumberOfBinsE != Dipole->GetNumberOfFrequencies())
+                        Receivers->SetNumberOfBinsE(Dipole->GetNumberOfFrequencies());
+                    break;
+            }
+        }
     }
 
     void EMEarth1D::AttachFieldPoints( std::shared_ptr<FieldPoints> recptr) {
 
         Receivers = recptr;
         if (Receivers == nullptr) {
-            std::cout << "nullptr Receivers in emearth1d.cpp " << std::endl;
+            std::cerr << "nullptr Receivers in emearth1d.cpp " << std::endl;
             return;
         }
 
-        // This has an implicid need to first set a source before receivers, users
+        // This has an implicit need to first set a source before receivers, users
         // will not expect this. Fix
         if (Dipole != nullptr) {
             switch (FieldsToCalculate) {
@@ -152,6 +166,10 @@ namespace Lemma {
         }
     }
 
+    void EMEarth1D::AttachLayeredEarthEM( std::shared_ptr<LayeredEarthEM> earthptr) {
+        Earth = earthptr;
+    }
+
     void EMEarth1D::AttachWireAntenna(std::shared_ptr<WireAntenna> antennae) {
         this->Antenna = antennae;
     }
@@ -227,6 +245,7 @@ namespace Lemma {
             if ( Antenna->IsHorizontallyPlanar() && ( HankelType == ANDERSON801 || HankelType == FHTKEY201  || HankelType==FHTKEY101 ||
                                                       HankelType == FHTKEY51    || HankelType == FHTKONG61  || HankelType == FHTKONG121 ||
                                                       HankelType == FHTKONG241  || HankelType == IRONS )) {
+
                 std::unique_ptr<ProgressBar> mdisp;
                 if (progressbar) {
                     mdisp = std::make_unique< ProgressBar >( Receivers->GetNumberOfPoints()*Antenna->GetNumberOfFrequencies() );
@@ -255,7 +274,6 @@ namespace Lemma {
                     #endif
                 }
 
-
             } else if (Receivers->GetNumberOfPoints() > Antenna->GetNumberOfFrequencies()) {
 
                 //** Progress display bar for long calculations */
@@ -280,7 +298,7 @@ namespace Lemma {
                     for (int irec=0; irec<Receivers->GetNumberOfPoints(); ++irec) {
                         if (!Receivers->GetMask(irec)) {
                             AntCopy->ApproximateWithElectricDipoles(Receivers->GetLocation(irec));
-                            for (unsigned int idip=0; idip<AntCopy->GetNumberOfDipoles(); ++idip) {
+                            for (int idip=0; idip < static_cast<int>(AntCopy->GetNumberOfDipoles()); ++idip) {
                                 auto tDipole = AntCopy->GetDipoleSource(idip);
                                 //#ifdef LEMMAUSEOMP
                                 //#pragma omp for schedule(static, 1)
@@ -318,7 +336,7 @@ namespace Lemma {
                             #pragma omp for schedule(static, 1)
                             #endif
                             for (int ifreq=0; ifreq<Antenna->GetNumberOfFrequencies(); ++ifreq) {
-                                for (unsigned int idip=0; idip<Antenna->GetNumberOfDipoles(); ++idip) {
+                                for (int idip=0; idip< static_cast<int>(Antenna->GetNumberOfDipoles()); ++idip) {
                                     auto tDipole = Antenna->GetDipoleSource(idip);
                                     // Propogation constant in free space
                                     Real wavef   = tDipole->GetAngularFrequency(ifreq) *
@@ -354,7 +372,7 @@ namespace Lemma {
                                 #ifdef LEMMAUSEOMP
                                 #pragma omp for schedule(static, 1)
                                 #endif
-                                for (int idip=0; idip < (int)Antenna->GetNumberOfDipoles(); ++idip) {
+                                for (int idip=0; idip<static_cast<int>(Antenna->GetNumberOfDipoles()); ++idip) {
                                     //#pragma omp critical
                                     //{
                                     //cout << "idip=" << idip << "\tthread num=" << omp_get_thread_num() << '\n';

Modules/FDEM1D/src/FHTAnderson801.cpp

@@ -940,7 +940,8 @@ namespace Lemma {
         icount = 0;
         Manager = KernelManager;
         this->kernelVec = KernelManager->GetSTLVector();
-        this->SetNumConv(nlag);
+        this->SetNumConv(nlag+1);
+
 #ifdef LEMMA_SINGLE_PRECISION
  		Compute(rho, 1, 1e-8);
 #else
@@ -958,7 +959,6 @@ namespace Lemma {
             SplineI->SetKnots( Arg, Zans.col(ii).imag() );
             splineVecImag.push_back(SplineI);
         }
-
     }
 
     void FHTAnderson801::SetLaggedArg(const Real& rho) {
@@ -998,7 +998,6 @@ namespace Lemma {
         // in release.
         #ifndef NDEBUG
 		if (rho<=0) {
-            //std::cout << "rho= " << rho << std::endl;
 			throw std::runtime_error("In Hankel 2 Argument rho <= 0; rho=" + to_string(rho) );
 		}
 

Modules/FDEM1D/src/PolygonalWireAntenna.cpp

@@ -117,25 +117,43 @@ namespace Lemma {
 	// ====================  OPERATIONS    =======================
 
 	void PolygonalWireAntenna::ApproximateWithElectricDipoles(const Vector3r &rp) {
+
         // Only resplit if necessary. Save a few cycles if repeated
         if ( (rRepeat-rp).norm() > 1e-16 ) {
-		    Dipoles.clear();
+
+            Dipoles.clear();
 
             ///////////////////
 		    // dipole array, this has impoved performance over directly pushing
 		    std::vector< std::shared_ptr<DipoleSource> >       xDipoles;
 
-		    // loop over all segments
+            // check to see if loop is closed
+            /*
+            int lastPoint = Points.cols()-1;
+            if ( (Points.col(0) - Points.col(lastPoint)).norm() > 1e-3 ) {
+                AddGroundingPoint( Points.col(0), Points.col(1)-Points.col(0), xDipoles );
+            }
+            */
+
+		    // loop over all segments, TODO fix for closed loops
 		    int iseg;
-            for (iseg=0; iseg<NumberOfPoints-1; ++iseg) {
+            for (iseg=0; iseg<NumberOfPoints-1; ++iseg) { // closed loop
+            //for (iseg=1; iseg<NumberOfPoints-2; ++iseg) { // grounded wire
 			    InterpolateLineSegment(Points.col(iseg), Points.col(iseg+1), rp, xDipoles);
 		    }
 
-            // Check to see if the loop is closed, if not, assume its grounded on ends,
-            if ( (Points.col(0)-Points.col(iseg)).norm() > 1e-3 ) {
-                xDipoles[0]->SetType(GROUNDEDELECTRICDIPOLE);
-                xDipoles.back()->SetType(GROUNDEDELECTRICDIPOLE);
+            // check to see if loop is closed
+            /*
+            if ( (Points.col(lastPoint)-Points.col( lastPoint-1 )).norm() > 1e-3 ) {
+                AddGroundingPoint( Points.col(lastPoint), Points.col(lastPoint)-Points.col(lastPoint-1), xDipoles );
             }
+            */
+
+            // Check to see if the loop is closed, if not, assume its grounded on ends,
+            //if ( (Points.col(0)-Points.col(iseg)).norm() > 1e-3 ) {
+            //    xDipoles[0]->SetType(GROUNDEDELECTRICDIPOLE);
+            //    xDipoles.back()->SetType(GROUNDEDELECTRICDIPOLE);
+            //}
 
             Dipoles = std::move(xDipoles);
             rRepeat = rp;
@@ -147,6 +165,23 @@ namespace Lemma {
         }
 	}
 
+    void PolygonalWireAntenna::AddGroundingPoint( const Vector3r &cp, const Vector3r& dir,
+                                std::vector< std::shared_ptr<DipoleSource> > &xDipoles) {
+		Real scale = (Real)(NumberOfTurns)*Current;
+        auto tx = DipoleSource::NewSP();
+		    tx->SetLocation(cp);
+		    tx->SetType(GROUNDINGPOINT);
+		    //tx->SetType(GROUNDEDELECTRICDIPOLE);
+		    //tx->SetType(MAGNETICDIPOLE);
+			tx->SetPolarisation(dir.array()); //dir.norm());
+			tx->SetFrequencies(Freqs);
+			tx->SetMoment(scale);
+			xDipoles.push_back(tx);
+        //std::cout << "cp " << cp.transpose() << std::endl;
+        //std::cout << "pol " << tx->GetPolarisation().transpose() << std::endl;
+        //std::cout << "moment " << tx->GetMoment() << std::endl;
+    }
+
 	Vector3r PolygonalWireAntenna::ClosestPointOnLine(const Vector3r &p1,
 					const Vector3r &p2, const Vector3r &tp) {
 

Modules/LemmaCore/CMakeLists.txt

@@ -33,6 +33,12 @@ if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT OR LEMMA_VTK8_SUPPORT )
 #	target_link_libraries(lemmacore "matplot")
 endif()
 
+
+if ( LEMMA_VTK9_SUPPORT ) 
+	target_link_libraries(lemmacore ${visibility}VTK::CommonCore VTK::IOXML VTK::FiltersHyperTree)
+	vtk_module_autoinit(TARGETS lemmacore MODULES VTK::CommonCore VTK::IOXML VTK::FiltersHyperTree)
+endif()
+
 # find_package(yaml-cpp) does not seem to properly define library names...
 # a better solution than this is welcome 
 

Modules/LemmaCore/include/RectilinearGridVTKExporter.h

@@ -26,7 +26,7 @@
 #include "LemmaObject.h"
 #include "RectilinearGrid.h"
 
-#include <vtkSmartPointer.h>
+#include "vtkSmartPointer.h"
 #include "vtkXMLRectilinearGridWriter.h"
 #include "vtkRectilinearGrid.h"
 #include "vtkDoubleArray.h"

Modules/LemmaCore/include/lemma.h

@@ -301,9 +301,10 @@
          @param NOSOURCETYPE is default.
          @param GROUNDEDELECTRICDIPOLE is an grounded electric dipole
          @param UNGROUNDEDELECTRICDIPOLE is an ungrounded electric dipole
+         @param GROUNDINGPOINT is a point of grounding in a bipole (or similiar) transmitter
          @param MAGNETICDIPOLE is a magnetic dipole
         */
-        enum DIPOLESOURCETYPE {NOSOURCETYPE, GROUNDEDELECTRICDIPOLE, UNGROUNDEDELECTRICDIPOLE, MAGNETICDIPOLE};
+        enum DIPOLESOURCETYPE {NOSOURCETYPE, GROUNDEDELECTRICDIPOLE, UNGROUNDEDELECTRICDIPOLE, GROUNDINGPOINT, MAGNETICDIPOLE};
 
         /// Only three polarizations are supported. They may be summed to
         /// approximate others

Modules/LemmaCore/src/CubicSplineInterpolator.cpp

@@ -173,17 +173,21 @@ namespace Lemma {
     //      Method:  Interpolate
     //--------------------------------------------------------------------------------------
     Real CubicSplineInterpolator::Interpolate ( const Real& x, int& i) {
+
         // O(n) search, could do bisection, but if these are sorted, then this is quick
-        while(Spline.x[i] < x && i<Spline.x.size()) {
+        while(Spline.x[i] < x && i<Spline.x.size()-1) {
             ++i;
         }
         --i;
 
 //         if ( x > Spline.x[i] ) {
-//             std::cout << "DOOM\t" << x << "\t" << i << "\t" << Spline.x[i];
-//             throw std::runtime_error("CubicSplineInterpolator::Interpolate ATTEMPT TO INTERPOLATE PAST LAST KNOT");
+//             std::cout << "DOOM in interplate\t x=" <<  x << "\ti=" << i << "\tSpline.x[i]=" << Spline.x[i] << std::endl;
+//             std::cout <<"Spline.x.size()" << Spline.x.size() << std::endl;
+//             std::cout << "Spline.x" << Spline.x.transpose() << std::endl;
+//             //throw std::runtime_error("CubicSplineInterpolator::Interpolate ATTEMPT TO INTERPOLATE PAST LAST KNOT");
 //         }
 
+
         return Spline.a[i] + Spline.b[i]*(x-Spline.x[i]) + Spline.c[i]*((x-Spline.x[i])*(x-Spline.x[i])) +
                Spline.d[i]*((x-Spline.x[i])*(x-Spline.x[i])*(x-Spline.x[i]) );
     }		// -----  end of method CubicSplineInterpolator::Interpolate  -----

Modules/LemmaCore/src/ProgressBar.cpp

@@ -2,6 +2,7 @@
 #include <iostream>
 #include <iomanip>
 #include <sstream>
+#include <algorithm> 
 
 #define LENGTH_OF_PROGRESS_BAR 55
 #define PERCENTAGE_BIN_SIZE (100.0/LENGTH_OF_PROGRESS_BAR)

Modules/LemmaCore/src/helper.cpp

@@ -185,6 +185,8 @@ std::string enum2String( const DIPOLESOURCETYPE& Type ) {
             return std::string("GROUNDEDELECTRICDIPOLE");
         case UNGROUNDEDELECTRICDIPOLE:
             return std::string("UNGROUNDEDELECTRICDIPOLE");
+        case GROUNDINGPOINT:
+            return std::string("GROUNDINGPOINT");
         case MAGNETICDIPOLE:
             return std::string("MAGNETICDIPOLE");
         default:
@@ -258,6 +260,7 @@ DIPOLESOURCETYPE string2Enum<DIPOLESOURCETYPE>( const std::string& str ) {
     if      (str == "NOSOURCETYPE")              return NOSOURCETYPE;
     if      (str == "GROUNDEDELECTRICDIPOLE")    return GROUNDEDELECTRICDIPOLE;
     if      (str == "UNGROUNDEDELECTRICDIPOLE")  return UNGROUNDEDELECTRICDIPOLE;
+    if      (str == "GROUNDINGPOINT")            return GROUNDINGPOINT;
     if      (str == "MAGNETICDIPOLE")            return MAGNETICDIPOLE;
     else {
         throw std::runtime_error("string not recognized as DipoleSource");

python/long.md

@@ -0,0 +1 @@
+Originally, Lemma is an Elecromagnetics modelling API. The scope of the project has expanded somewhat from this, and currently you can find documentation at https://lemmasoftware.org

python/publish.sh

@@ -2,8 +2,8 @@
 rm -rf build dist clean pyLemma.egg.info
 python setup.py build
 python setup.py bdist_wheel
-#auditwheel repair ./dist/pyLemma*.whl -w ./clean
-#twine upload build/* #clean/*
+auditwheel repair ./dist/pyLemma*.whl -w ./clean
+twine upload build/* #clean/*
 
 #rm -rf dist
 #python setup.py build

python/setup.py

@@ -19,13 +19,16 @@ class InstallPlatlib(install):
         if self.distribution.has_ext_modules():
             self.install_lib = self.install_platlib
 
+with open("../README.md", "r") as fh:
+    long_description = fh.read()
+
 setup(
   name             = 'pyLemma',
-  version          = '0.0.13', 
+  version          = '0.4.0', 
   author           = 'Trevor Irons and others',
   author_email     = 'Trevor.Irons@lemmasoftware.org',
-  description      = 'A short description of the app/lib',
-  long_description = 'A longer one',
+  description      = 'PyLemma is a wrapper to Lemma',
+  long_description = long_description,
   classifiers=[
         'Development Status :: 3 - Alpha',
         'Intended Audience :: Developers',

vim/c.vim

@@ -9,7 +9,7 @@ highlight eType ctermfg=Magenta guifg=Magenta
 syn keyword eType  MAGUNITS  TEMPUNITS  TIMEUNITS FREQUENCYUNITS FEMCOILORIENTATION ORIENTATION FIELDTYPE FIELDCOMPONENT SPATIALCOORDINANT HANKELTRANSFORMTYPE FIELDCALCULATIONS WINDOWTYPE DIPOLESOURCETYPE 
 
 highlight eeType ctermfg=Cyan guifg=Cyan
-syn keyword eeType  TESLA NANOTESLA GAUSS CELCIUS KELVIN SEC MILLISEC MICROSEC NANOSEC PICOSEC HZ KHZ MHZ GHZ COAXIAL COPLANAR HFIELDREAL HFIELDIMAG EFIELDREAL EFIELDIMAG XCOMPONENT YCOMPONENT ZCOMPONENT XCOORD YCOORD ZCOORD X Y Z NX  NY  NZ ANDERSON801 CHAVE FHTKEY201 FHTKEY101 FHTKEY51 FHTKONG241 FHTKONG121 FHTKONG61 QWEKEY E H BOTH HAMMING HANNING RECTANGULAR NOSOURCETYPE GROUNDEDELECTRICDIPOLE UNGROUNDEDELECTRICDIPOLE MAGNETICDIPOLE
+syn keyword eeType  TESLA NANOTESLA GAUSS CELCIUS KELVIN SEC MILLISEC MICROSEC NANOSEC PICOSEC HZ KHZ MHZ GHZ COAXIAL COPLANAR HFIELDREAL HFIELDIMAG EFIELDREAL EFIELDIMAG XCOMPONENT YCOMPONENT ZCOMPONENT XCOORD YCOORD ZCOORD X Y Z NX  NY  NZ ANDERSON801 CHAVE FHTKEY201 FHTKEY101 FHTKEY51 FHTKONG241 FHTKONG121 FHTKONG61 QWEKEY E H BOTH HAMMING HANNING RECTANGULAR NOSOURCETYPE GROUNDEDELECTRICDIPOLE GROUNDINGPOINT UNGROUNDEDELECTRICDIPOLE MAGNETICDIPOLE
 
 " Namespaces
 highlight nspace ctermfg=Red guifg=Red