Moved LemmaCore into Modules, cleaned up build process and documentation.

CMakeLists.txt

@@ -70,7 +70,7 @@ ExternalProject_Add(YAML-CPP
 	UPDATE_COMMAND ""
 	PATCH_COMMAND ""
     PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/yaml-cpp
-    CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} -DBUILD_SHARED_LIBS=ON -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} /
+    CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} /
 		-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
 )
 add_compile_options(-DHAVE_YAMLCPP) 
@@ -89,6 +89,8 @@ else()
 		BUILD_COMMAND ""
 		INSTALL_COMMAND ""
     	PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/CxxTest
+    	CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} /
+			-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
 	)
 	set (CXXTEST_INCLUDE_DIR  "${CMAKE_CURRENT_BINARY_DIR}/external/CxxTest/src/CxxTest/")
 	set (CXXTEST_PYTHON_TESTGEN_EXECUTABLE "${CMAKE_CURRENT_BINARY_DIR}/external/CxxTest/src/CxxTest/bin/cxxtestgen" )
@@ -104,8 +106,8 @@ if ( LEMMA_TINYXML_SUPPORT )
 		UPDATE_COMMAND ""
 		PATCH_COMMAND ""
     	PREFIX ${CMAKE_CURRENT_BINARY_DIR}/external/tinyxml2
-    	CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} -DBUILD_SHARED_LIBS=ON -DBUILD_STATIC_LIBS=ON 
-		#-DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS}
+    	CMAKE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${CMAKE_INSTALL_PREFIX} -DBUILD_SHARED_LIBS=${BUILD_SHARED_LIBS} /
+			-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
 	)
 	add_compile_options(-DTINYXMLSUPPORT) 
 endif()
@@ -138,7 +140,7 @@ if ( LEMMA_VTK6_SUPPORT )
 	include(${VTK_USE_FILE}) 
 	add_compile_options(-DLEMMAUSEVTK) 
 	# Compile Matplot_vtk if VTK is present
-	add_subdirectory(Matplot_vtk)	
+	#add_subdirectory(Matplot_vtk)	
 	include_directories ("${PROJECT_SOURCE_DIR}/Matplot_vtk")
 endif()
 
@@ -152,7 +154,7 @@ if ( LEMMA_VTK7_SUPPORT )
 	include(${VTK_USE_FILE}) 
 	add_compile_options(-DLEMMAUSEVTK) 
 	# Compile Matplot_vtk if VTK is present
-	add_subdirectory(Matplot_vtk)	
+	#add_subdirectory(Matplot_vtk)	
 	include_directories ("${PROJECT_SOURCE_DIR}/Matplot_vtk")
 endif()
 
@@ -178,6 +180,8 @@ else()
 	message(FATAL_ERROR "The compiler ${CMAKE_CXX_COMPILER} has no C++11 or C++14 support. Please use a different C++ compiler.")
 endif()
 
+option( LEMMA_BUILD_EXAMPLES "Compile example Lemma applications" OFF )
+
 option( LEMMA_USE_OPENMP  "Use OpenMP in Lemma" OFF)
 CHECK_CXX_COMPILER_FLAG(-fopenmp COMPILER_SUPPORTS_OPENMP)
 if(COMPILER_SUPPORTS_OPENMP AND LEMMA_USE_OPENMP )
@@ -190,42 +194,29 @@ set(LEMMA_VERSION_PATCH "1")
 set(LEMMA_VERSION "${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}.${LEMMA_VERSION_PATCH}")
 
 configure_file (
-	"${PROJECT_SOURCE_DIR}/include/LemmaConfig.h.in"
+	"${PROJECT_SOURCE_DIR}/config/LemmaConfig.h.in"
 	"${PROJECT_BINARY_DIR}/include/LemmaConfig.h"
 )
-include_directories("${PROJECT_BINARY_DIR}/include")
-
-include_directories ("${PROJECT_SOURCE_DIR}/LemmaCore/include")
 
+include_directories("${PROJECT_BINARY_DIR}/include")
 include_directories ("${CMAKE_INSTALL_PREFIX}/include/")
 link_directories ("${CMAKE_INSTALL_PREFIX}/lib/")
 
-
+#######################################################################
+# CXXTEST Framework
 if(CXXTEST_FOUND)
 	option (LEMMA_ENABLE_TESTING  "Turn on unit testing" ON)
 	if (LEMMA_ENABLE_TESTING)
-    include_directories(${CXXTEST_INCLUDE_DIR})
-    enable_testing()
-    #CXXTEST_ADD_TEST(tester_executable runner.cpp
-    #                    ${CMAKE_CURRENT_SOURCE_DIR}/src/myTestsuite.h)
-	#CXXTEST_ADD_TEST(unittest_foo foo_test.cc
-    #                  ${CMAKE_CURRENT_SOURCE_DIR}/src/myTestSuite.h)
-    #target_link_libraries(unittest_foo foo) # as needed
-#	else()
-#	message(FATAL_ERROR "CxxTest not found.")
+    	#include_directories(${CXXTEST_INCLUDE_DIR}) # Add this in testing CMakeLists.txt instead, cleaner
+    	enable_testing()
 	endif()
 endif()
 
-add_subdirectory (LemmaCore) 
-
 ########################################################################
 # Build Modules
 add_subdirectory (Modules) 
 
 ########################################################################
-# 
-include_directories(${CMAKE_INSTALL_PREFIX}/include)
-
 # add a target to generate API documentation with Doxygen
 # ALL make documentation build by default if possible
 find_package(Doxygen)

Documentation/Doxyfile.in

@@ -586,8 +586,9 @@ USE_MATHJAX = YES
 # with spaces.
 #INPUT                  = ./include ./Documentation/dox ./modules/Wavelet/include ./modules/Merlin/include ./modules/FEM4EllipticPDE/include ./modules/EMSchur3D/include
 #INPUT                  = @CMAKE_CURRENT_SOURCE_DIR@/LemmaCore/include @CMAKE_CURRENT_SOURCE_DIR@/Modules/FEM4EllipticPDE.git/include
-INPUT   = @CMAKE_CURRENT_SOURCE_DIR@/LemmaCore/include @CMAKE_CURRENT_SOURCE_DIR@/Modules/*/include
+INPUT   = @CMAKE_CURRENT_SOURCE_DIR@/Modules/*/include
 INPUT  += @CMAKE_CURRENT_SOURCE_DIR@/Modules/
+INPUT  += @CMAKE_CURRENT_SOURCE_DIR@/Documentation/dox
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is

Modules/FDEM1D/CMakeLists.txt

@@ -1,24 +1,26 @@
-include_directories(${CMAKE_INSTALL_PREFIX}/include)
-
-include_directories( "${CMAKE_CURRENT_SOURCE_DIR}/include" )
 add_subdirectory("src")
-
 add_library( fdem1d ${FEM1DSOURCE} )  
-target_link_libraries(fdem1d "lemmacore")
+target_include_directories( fdem1d PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include" )
+
 set_target_properties(fdem1d PROPERTIES 
 	VERSION  "${LEMMA_VERSION}"
 	SOVERSION "${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}"
-	PROJECT_LABEL "Lemma ${LABEL_SUFFIX}"
+	PROJECT_LABEL "FDEM1D ${LABEL_SUFFIX}"
 )
 
-install ( TARGETS fdem1d DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
+# Linking
+target_link_libraries(fdem1d "lemmacore")
 
-if (LEMMA_BUILD_EXAMPLES)
-	add_subdirectory(examples)
-endif()
-	
+# Testing
 if (LEMMA_ENABLE_TESTING)
 	add_subdirectory(testing)
 endif()
 
+# Install
 install ( TARGETS fdem1d DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
+
+# Examples
+if (LEMMA_BUILD_EXAMPLES)
+	add_subdirectory(examples)
+endif()
+	

Modules/FDEM1D/testing/CMakeLists.txt

@@ -1,3 +1,3 @@
-
+include_directories(${CXXTEST_INCLUDE_DIR})
 CXXTEST_ADD_TEST(unittestFEM1D_GetNameCheck GetNameCheck.cc ${CMAKE_CURRENT_SOURCE_DIR}/GetNameCheck.h)
 target_link_libraries(unittestFEM1D_GetNameCheck "lemmacore" "fdem1d")

Modules/LemmaCore/CMakeLists.txt

@@ -0,0 +1,39 @@
+add_subdirectory("src")
+add_library( lemmacore ${SOURCE} )  
+target_include_directories( lemmacore PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include" )
+
+set_target_properties(lemmacore PROPERTIES 
+	VERSION  "${LEMMA_VERSION}"
+	SOVERSION "${LEMMA_VERSION_MAJOR}.${LEMMA_VERSION_MINOR}"
+	PROJECT_LABEL "LemmaCore ${LABEL_SUFFIX}"
+)
+
+# External project dependencies 
+add_dependencies(lemmacore YAML-CPP)
+if (LEMMA_TINYXML_SUPPORT)
+	add_dependencies(lemmacore TINYXML2)
+endif()
+if (LEMMA_MATIO_SUPPORT)
+	add_dependencies(lemmacore MATIO)
+endif()
+
+# Linking
+if ( LEMMA_VTK6_SUPPORT OR LEMMA_VTK7_SUPPORT ) 
+	target_link_libraries(lemmacore ${VTK_LIBRARIES})
+	target_link_libraries(lemmacore "matplot")
+endif()
+
+target_link_libraries(lemmacore "yaml-cpp")
+
+# Testing
+if (LEMMA_ENABLE_TESTING)
+	add_subdirectory(testing)
+endif()
+
+# Install
+install ( TARGETS lemmacore DESTINATION ${CMAKE_INSTALL_PREFIX}/lib )
+
+# Examples
+if (LEMMA_BUILD_EXAMPLES)
+	add_subdirectory(examples)
+endif()

Modules/LemmaCore/LemmaCore.dox

@@ -0,0 +1,24 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/15/2016 04:28:16 PM
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Lemma Software, LLC
+ */
+
+/**
+    \defgroup LemmaCore LemmaCore
+    @brief    LemmaCore provides the general framework and skeleton for applications.
+    @details  This module is the glue holding Lemma together.
+*/
+

Modules/LemmaCore/examples/CMakeLists.txt

@@ -0,0 +1,33 @@
+#FILE(GLOB  examples  "*.cpp")
+
+#add_executable( datatem datatem.cpp )
+#target_link_libraries(datatem "lemmacore")
+
+#add_executable( wireantenna wireantenna.cpp )
+#target_link_libraries(wireantenna "lemmacore")
+
+#add_executable( DCIPElectrode  DCIPElectrode.cpp )
+#target_link_libraries(DCIPElectrode "lemmacore")
+
+#add_executable( digitalfiltercostrans  digitalfiltercostrans.cpp )
+#target_link_libraries(digitalfiltercostrans "lemmacore")
+
+#add_executable( dipolefreqsweep  dipolefreqsweep.cpp )
+#target_link_libraries(dipolefreqsweep "lemmacore")
+
+#add_executable( dipolesource  dipolesource.cpp )
+#target_link_libraries(dipolesource "lemmacore")
+
+#add_executable( edipole  edipole.cpp )
+#target_link_libraries(edipole "lemmacore")
+
+#add_executable( emdipearth1d  emdipearth1d.cpp )
+#target_link_libraries(emdipearth1d "lemmacore")
+
+#add_executable( femforward  femforward.cpp )
+#target_link_libraries(femforward "lemmacore")
+
+add_executable( filter  filter.cpp )
+target_link_libraries(filter "lemmacore")
+
+

Modules/LemmaCore/examples/filter.cpp

@@ -0,0 +1,52 @@
+// ===========================================================================
+//
+//       Filename:  utfilter.cpp
+//
+//        Created:  07/20/2010 10:12:07 AM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @author   Trevor Irons
+  @date     07/20/2010
+  @version   0.0
+ **/
+
+#include "WindowFilter.h"
+using namespace Lemma;
+
+int main() {
+
+    auto Hamming = WindowFilter::NewSP();
+
+        Hamming->SetWindowType( HANNING );
+        Hamming->SetSamplingInterval(1e-4);
+        Hamming->SetBandwidth(300);
+        Hamming->SetNumberOfSamples(1./1e-4);
+
+    std::cout << Hamming->GetFilterCoefficients() //.transpose()
+              << std::endl;
+
+
+}

Modules/LemmaCore/examples/utORS.cpp

@@ -0,0 +1,342 @@
+// ===========================================================================
+//
+//       Filename:  utdipolesource.cpp
+//
+//    Description:
+//
+//        Version:  0.0
+//        Created:  12/02/2009 11:57:14 AM
+//       Revision:  none
+//       Compiler:  g++ (c++)
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include <iostream>
+#include <fstream>
+
+#include "dipolesource.h"
+#include "layeredearth.h"
+#include "receiverpoints.h"
+#include "emearth1d.h"
+#include "WireAntenna.h"
+#include "PolygonalWireAntenna.h"
+
+#if  LEMMAUSEVTK
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderLargeImage.h"
+#include "vtkPNGWriter.h"
+#include "vtkRectilinearGrid.h"
+#include "vtkXMLRectilinearGridWriter.h"
+#include "vtkDoubleArray.h"
+#include "vtkFieldData.h"
+#include "vtkCellData.h"
+#endif     // -----  not LEMMA_USE_VTK  -----
+
+// For testing purposes disable VTK and run scale.sh
+//#undef LEMMAUSEVTK
+
+#include "timer.h"
+
+using namespace Lemma;
+
+double randDouble(double low, double high) {
+	//srand(time(0));
+	double temp;
+	/* swap low & high around if the user makes no sense */
+	if (low > high)	{
+	temp = low;
+	low = high;
+	high = temp;
+	}
+	/* calculate the random number & return it */
+	temp = (rand() / (static_cast<double>(RAND_MAX) + 1.0))
+	* (high - low) + low;
+	return temp;
+}
+
+
+int main() {
+
+	// Keep track of time
+	jsw_timer timer;
+
+
+	srand(time(0));
+
+	PolygonalWireAntenna *pa = PolygonalWireAntenna::New();
+	pa->SetNumberOfFrequencies(1);
+	pa->SetFrequency(0, 540000);
+    Real Depth  =   370;  // nominal depth
+    Real Width  = 0.14;  // transmitter height
+    Real Height = 0.70;//  .14;  //.014;  // transmitter width
+
+
+	pa->SetNumberOfPoints(5);
+    /*
+    pa->SetPoint(0, Vector3r(       0,       0, Depth));
+	pa->SetPoint(1, Vector3r(       0,       0, Depth+Height));
+	pa->SetPoint(2, Vector3r(       0,   Width, Depth+Height));
+	pa->SetPoint(3, Vector3r(       0,   Width, Depth));
+	pa->SetPoint(4, Vector3r(       0,       0, Depth));
+    */
+
+    pa->SetPoint(0, Vector3r(   -Width/2., -Height/2., Depth));
+	pa->SetPoint(1, Vector3r(    Width/2., -Height/2., Depth));
+	pa->SetPoint(2, Vector3r(    Width/2.,  Height/2., Depth));
+	pa->SetPoint(3, Vector3r(   -Width/2.,  Height/2., Depth));
+	pa->SetPoint(4, Vector3r(   -Width/2., -Height/2., Depth));
+
+	pa->SetCurrent(1.);
+    pa->SetNumberOfTurns(6);
+	//Vector3r rp = Vector3r::Random(3);
+	//rp << 150., 10, 0.;
+	//rp << -27.1456, 15.2350, -1e-3;
+	//rp << randDouble(-35,35), randDouble(-35,35), randDouble(-35,35);
+	//rp << 3.22806, -13.1548, 14.9695;
+	//rp.setRandom(3);
+	//std::cout <<  "rp " << rp.transpose() << std::endl;
+	//pa->ApproximateWithElectricDipoles(rp);
+
+
+	// Define model
+    Real Sigma = 1.0/20. ; // .05;
+	VectorXcr sigma(2);
+		sigma << Complex(0.,0), Complex(Sigma,0);
+	VectorXr  thick(1);
+		thick << 10;
+	LayeredEarthEM *earth = LayeredEarthEM::New();
+		earth->SetNumberOfLayers(2);
+		earth->SetLayerConductivity(sigma);
+		//earth->SetLayerThickness(thick);
+
+	// Receivers
+	ReceiverPoints *receivers = ReceiverPoints::New();
+		Vector3r loc;
+		Real ox    =    -5.*Width - .003373;
+		Real oy    =    -3.*Height - .003373;
+		Real oz    =    Depth - Height/2. - .003373;
+		int nx     =   120; // 60;
+		int ny     =   180; //  50;
+		int nz     =   100; // 40;
+		Real hx    =    11.*Width/nx;
+		Real hy    =     7.*Height/ny;
+		Real hz    =     1.*Height/nz;
+
+        receivers->SetNumberOfReceivers(nx*ny*nz);
+		int ir = 0;
+		for (int iz=0; iz<nz; ++iz) {
+		for (int iy=0; iy<ny; ++iy) {
+		for (int ix=0; ix<nx; ++ix) {
+			loc << ox+ix*hx, oy+iy*hy, oz+iz*hz;
+			receivers->SetLocation(ir, loc);
+			++ ir;
+		}
+        }
+        }
+
+	// EmEarth
+	EMEarth1D  *EmEarth = EMEarth1D::New();
+		//EmEarth->AttachWireAntenna(wire);
+		EmEarth->AttachWireAntenna(pa);
+		EmEarth->AttachLayeredEarthEM(earth);
+		EmEarth->AttachReceiverPoints(receivers);
+		EmEarth->SetFieldsToCalculate(E);
+        //EmEarth->SetHankelTransformMethod(GAUSSIANQUADRATURE);
+        EmEarth->SetHankelTransformMethod(ANDERSON801);
+
+	// Do calculation
+	timer.begin();
+	EmEarth->CalculateWireAntennaFields();
+	Real paTime = timer.end();
+	std::cout << "Polygonal wire antennae time: " << paTime << "\n";
+
+	//EmEarth->AttachWireAntenna(wire);
+	//timer.begin();
+	//EmEarth->CalculateWireAntennaFields();
+	//Real waTime = timer.end();
+	//std::cout << "Fixed wire antennae time: " << waTime << "\n";
+/*
+	depth = depth2;
+	std::fstream real("reale_lay.dat", std::ios::out);
+	std::fstream imag("image_lay.dat", std::ios::out);
+	for (int iz=0; iz<nz; ++iz) {
+		Vector3cr temp = receivers->GetEfield(0,iz);
+		real << ox << "\t" << oy << "\t" << depth << "\t"
+				<< temp(0).real() << "\t" << temp(1).real()
+				<< "\t" << temp(2).real() << std::endl;
+		imag << ox << "\t" << oy << "\t" << depth << "\t"
+				<< std::imag(temp(0)) << "\t" << std::imag(temp(1))
+				<< "\t" << std::imag(temp(2)) << std::endl;
+		depth += dx;
+	}
+	real.close();
+	imag.close();
+*/
+
+	//wire->Delete();
+#if  LEMMAUSEVTK
+
+    // Set Coordinates
+    vtkDoubleArray *xCoords = vtkDoubleArray::New();
+    xCoords->InsertNextValue(ox-hx/2.);
+    double xm1 = ox-hx/2.;
+    for (int ix=0; ix<nx; ix++) {
+        xCoords->InsertNextValue(xm1 + hx);
+        xm1 += hx;
+    }
+
+    vtkDoubleArray *yCoords = vtkDoubleArray::New();
+    yCoords->InsertNextValue(oy-hy/2.);
+    double ym1 = oy-hy/2.;
+    for (int iy=0; iy<ny; iy++) {
+        yCoords->InsertNextValue(ym1 + hy);
+        ym1 += hy;
+    }
+
+    vtkDoubleArray *zCoords = vtkDoubleArray::New();
+    zCoords->InsertNextValue(oz-hz/2.);
+    double zm1 = oz-hz/2.;
+    for (int iz=0; iz<nz; iz++) {
+        zCoords->InsertNextValue(zm1 + hz);
+        zm1 += hz;
+    }
+
+    vtkDoubleArray *EReal = vtkDoubleArray::New();
+    vtkDoubleArray *EImag = vtkDoubleArray::New();
+    vtkDoubleArray *Watts = vtkDoubleArray::New();
+    EReal->SetNumberOfComponents(3);
+    EImag->SetNumberOfComponents(3);
+    Watts->SetNumberOfComponents(1);
+
+    ir = 0;
+    Real WattsTotal(0);
+    Real WattsInterior(0);
+    for (int iz=0; iz<nz; ++iz) {
+    for (int iy=0; iy<ny; ++iy) {
+    for (int ix=0; ix<nx; ++ix) {
+        //sigmaArray->InsertTuple1(i,  sigma[ix][iy][iz] );
+		Vector3cr E = receivers->GetEfield(0, ir);
+        EReal-> InsertTuple3(ir, real(E(0)), real(E(1)), real(E(2)));
+        EImag-> InsertTuple3(ir, imag(E(0)), imag(E(1)), imag(E(2)));
+        //std::cout << std::abs(ox+ix*hx) << "\t"  <<  Width/2. << endl; // && std::abs(oy+iy*hy) > Height/2. ) { // && std::abs(oz+iz*hz - Depth) > Width/2. ) {
+        //if ( std::abs( ox+ix*hx ) < Width/2. && std::abs(oy+iy*hy) < Height/2. && std::abs(oz+iz*hz - Depth) < Width/2. ) {
+        if ( std::sqrt( std::pow(ox+ix*hx,2)  + std::pow(oz+iz*hz - Depth, 2) ) < .085 && std::abs(oy+iy*hy) < Height/2. ) {
+            Watts-> InsertTuple1(ir,  1e-20 );
+            WattsInterior += .5*  (( pow((std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
+        } else {
+            Watts-> InsertTuple1(ir, .5* (pow( (std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
+            //Watts-> InsertTuple1(ir,  1e-20 );
+            WattsTotal += .5*  (( pow((std::abs(E(0)) + std::abs(E(1)) + std::abs(E(2))), 2)*Sigma)*hx*hy*hz );
+        }
+		++ ir;
+    }
+    }
+    }
+    std::cout << "Total Power: " << Sigma << "\t" << WattsTotal << "\t" << WattsInterior << endl;
+
+    EReal->SetName("E_real");
+    EImag->SetName("E_imag");
+    Watts->SetName("Power");
+
+    vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
+        rgrid->SetDimensions(nx+1,ny+1,nz+1);
+        rgrid->SetXCoordinates(xCoords);
+        rgrid->SetYCoordinates(yCoords);
+        rgrid->SetZCoordinates(zCoords);
+
+        rgrid->GetCellData()->AddArray(EReal);
+        rgrid->GetCellData()->AddArray(EImag);
+        rgrid->GetCellData()->AddArray(Watts);
+        //rgrid->Update();
+
+    vtkXMLRectilinearGridWriter *gridWrite = vtkXMLRectilinearGridWriter::New();
+        gridWrite->SetInputData(rgrid);
+        gridWrite->SetFileName("ors.vtr");
+        gridWrite->Write();
+        //gridWrite->Update();
+
+#endif
+	EmEarth->Delete();
+	receivers->Delete();
+	earth->Delete();
+
+#if  LEMMAUSEVTKX
+	// Create the usual rendering stuff.
+	vtkRenderer           *renderer = vtkRenderer::New();
+	vtkRenderWindow         *renWin = vtkRenderWindow::New();
+	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+
+	std::cout << "wire antennae approximating " << std::endl;
+	loc << 0, .5*Width, Depth+.5*Height;
+	pa->ApproximateWithElectricDipoles(loc);
+    std::cout <<  "Wire approximated with " << pa->GetNumberOfDipoles() << std::endl;
+
+    vtkActor **pdipActors = new vtkActor*[pa->GetNumberOfDipoles()];
+	for (int id=0; id<pa->GetNumberOfDipoles(); ++id) {
+		pdipActors[id] = pa->GetVtkActor(id);
+		renderer->AddActor(pdipActors[id]);
+	}
+
+    /*
+    vtkActor **dipActors = new vtkActor*[wire->GetNumberOfDipoles()];
+	for (int id=0; id<wire->GetNumberOfDipoles(); ++id) {
+		dipActors[id] = wire->GetVtkActor(id);
+		renderer->AddActor(dipActors[id]);
+	}
+    */
+
+	renderer->SetBackground(1,1,1);
+
+	// Render the window
+	renWin->AddRenderer(renderer);
+	renWin->SetWindowName("Wire antennae");
+
+	iren->SetRenderWindow(renWin);
+	iren->Initialize();
+	iren->Start();
+	iren->Render();
+
+    #if 0
+	cout << "Enter File name?: ";
+	std::string pngName;
+	std::cin >> pngName;
+	vtkPNGWriter *pngwrite = vtkPNGWriter::New();
+	vtkRenderLargeImage *renlarge = vtkRenderLargeImage::New();
+	renlarge->SetInput(renderer);
+	renlarge->SetMagnification(2);
+	pngwrite->SetInputConnection(renlarge->GetOutputPort());
+	pngName.append(".png");
+	pngwrite->SetFileName(pngName.c_str());
+	pngwrite->Write();
+    #endif
+
+#endif     // -----  not LEMMA_USE_VTK  -----
+
+
+	//std::cout << *pa << std::endl;
+	//pa->Delete();
+
+	return 0;
+}

Modules/LemmaCore/examples/utQWEKey.cpp

@@ -0,0 +1,34 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      02/12/2014 01:01:29 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "QWEKey.h"
+using namespace Lemma;
+
+int main() {
+
+    // TODO chech
+
+    QWEKey* QWE = QWEKey::New();
+    QWE->TestPrivate( 5 );
+    QWE->Delete();
+
+    return 0;
+}
+
+

Modules/LemmaCore/examples/utSerialize.cpp

@@ -0,0 +1,86 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/29/2014 01:31:26 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "Lemma"
+
+using namespace Lemma;
+
+#ifdef  HAVE_YAMLCPP
+#include "yaml-cpp/yaml.h"
+
+int main() {
+
+    std::ofstream ofstr("output.yaml");
+    //ofstr << "%YAML 1.2\n" << "---\n";
+    //ofstr << "LEMMA_CLASSES:\n";
+
+    // This works for individual classes or even groups of unrelated classes. But
+    // How should references work? One option would be to hack YAML and allow for references
+    // across 'Docs'. This is maybe undesireable, as then it's not valid YAML anymore.
+    // OR we could have some other kind of IDENTIFIER that could be used and Connections could be
+    // made AFTER the whole file is processed?
+
+    // Alternative is each class is written to it's own file. And then implicit on the Serialize function, the file tags are written
+    // into the YAML instead of the pointer. Alternatively, tags could do this too.
+
+    // Or file could be saved as pointer.yaml
+
+    // *
+    // * Or a file contains all the subclasses it point to. So that it's self-enclosed. This is nice as it's quite
+    // * flexible. The only downwside is possible repeat of classes. Basically if the API is careful though, only 1
+    // * instance will need to be serialized.
+    // *
+
+
+    TEMTransmitter* Seq = TEMTransmitter::New();
+        Seq->SetRepFrequency( 20, KHZ );
+        VectorXr Times (18);
+        VectorXr Amps (18);
+        Times << 0.0, 0.03051, 0.10267, 0.19408, 0.19889, 0.21332, 0.74249, 1.3775, 1.83452, 2.52245, \
+                3.191132, 3.9031135, 4.0, 4.00484486, 4.123904, 4.200182, 4.20732, 4.212946;
+        Amps << 0.0, 14.71872, 62.34372, 114.84372, 117.84372, 118.96872, 118.96872, 118.96872, 118.96872,\
+                118.59372, 119.34372, 120.0, 120.0, 117.94176, 47.60364, 0.8905848, 0.1203888, 0.0;
+        Seq->SetWaveform(Times, Amps, MILLISEC);
+    // OK dump to disk
+    ofstr << *Seq ;
+    ofstr.close();
+    Seq->Delete();
+
+
+    // OK now load
+    std::ifstream ifstr("output.yaml");
+    std::vector<YAML::Node> nodes = YAML::LoadAll(ifstr);
+
+    TEMTransmitter* Seq2 = TEMTransmitter::DeSerialize( nodes[0] );
+    std::cout << "Loaded\n" << *Seq2;
+    Seq2->Delete();
+    ifstr.close();
+
+    //YAML::Node nodes = YAML::LoadFile("output.yaml");
+    //std::cout << nodes[0] << std::endl;
+    //std::cout << nodes[1] << std::endl;
+
+    exit(EXIT_SUCCESS);
+}
+
+#else
+int main() {
+
+}
+#endif

Modules/LemmaCore/examples/utSkyTEM.cpp

@@ -0,0 +1,31 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      02/27/2015 11:02:43 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2015, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2015, Trevor Irons
+ */
+
+// SkyTEM forward modelling programme
+
+#include "TEMSurvey.h"
+
+int main() {
+
+    std::cout << "SkyTEM modelling\n";
+
+
+
+}
+

Modules/LemmaCore/examples/utTEMSurvey.cpp

@@ -0,0 +1,197 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      10/06/2014 10:37:30 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "Lemma"
+using namespace Lemma;
+
+int main() {
+
+    // Each pulse sequence functions as an autonomous transmitter. A record may contain multiple transmitters (moments).
+    // There is some small overhead for wire loop locations, but generality is retained doing this, as well as monitoring of
+    // slight changes in geometry for both pulses.
+    TEMTransmitter* TxHM = TEMTransmitter::New();
+        TxHM->SetRepFrequency( 20, KHZ );
+        VectorXr Times (18);
+        VectorXr Amps (18);
+        Times << 0.0, 0.03051, 0.10267, 0.19408, 0.19889, 0.21332, 0.74249, 1.3775, 1.83452, 2.52245, \
+                3.191132, 3.9031135, 4.0, 4.00484486, 4.123904, 4.200182, 4.20732, 4.212946;
+        Amps << 0.0, 14.71872, 62.34372, 114.84372, 117.84372, 118.96872, 118.96872, 118.96872, 118.96872,\
+                118.59372, 119.34372, 120.0, 120.0, 117.94176, 47.60364, 0.8905848, 0.1203888, 0.0;
+        TxHM->SetWaveform( Times, Amps, MILLISEC );
+        // Define wire loop
+        TxHM->SetNumberOfPoints(8);
+        TxHM->SetPoint(0, Vector3r( -16.10,    2.13,  -34));
+        TxHM->SetPoint(1, Vector3r(  -7.51,   10.72,  -34));
+        TxHM->SetPoint(2, Vector3r(   7.51,   10.72,  -34));
+        TxHM->SetPoint(3, Vector3r(  14.92,    3.31,  -34));
+        TxHM->SetPoint(4, Vector3r(  14.92,   -3.31,  -34));
+        TxHM->SetPoint(5, Vector3r(   7.51,  -10.72,  -34));
+        TxHM->SetPoint(6, Vector3r(  -7.51,  -10.72,  -34));
+        TxHM->SetPoint(7, Vector3r( -16.10,   -2.13,  -34));
+        TxHM->SetNumberOfTurns(8);
+
+    // Each Transmitter needs its own `Receiver(s)', the only difference may be the time gates, but that's OK.
+    //     It's a minor amount of overhead for greatly improved genearality. As sometimes different gates are
+    //     masked.
+    TEMInductiveReceiver* RxHM = TEMInductiveReceiver::New();
+        RxHM->SetComponent( ZCOMPONENT );      // What about Overloaded to take (X), (X,Y), or (X,Y,Z) or what about tilt?
+        RxHM->SetMoment( 1 );                  // Normalized
+        RxHM->SetReferenceTime( 4., MILLISEC );
+        RxHM->SetRxLocation( (Vector3r() << -16.80, 0, -36.00).finished() );
+
+        // Gate Centres  ms
+        VectorXr centres (37);
+        centres << 7.15000000e-04,   2.21500000e-03,   4.21500000e-03,   6.21500000e-03,
+                   8.21500000e-03,   1.02150000e-02,   1.22150000e-02,   1.47150000e-02,
+                   1.82150000e-02,   2.27150000e-02,   2.82150000e-02,   3.52150000e-02,
+                   4.42150000e-02,   5.57150000e-02,   7.02150000e-02,   8.82150000e-02,
+                   1.10715000e-01,   1.38715000e-01,   1.74215000e-01,   2.19715000e-01,
+                   2.76715000e-01,   3.48715000e-01,   4.39715000e-01,   5.53715000e-01,
+                   6.97715000e-01,   8.79215000e-01,   1.10771500e+00,   1.39621500e+00,
+                   1.76021500e+00,   2.21871500e+00,   2.79671500e+00,   3.52571500e+00,
+                   4.44471500e+00,   5.60321500e+00,   7.06321500e+00,   8.90421500e+00,
+                   1.10667200e+01;
+        //centres.array() += 4.;
+
+        // Gate Widths  ms
+        VectorXr widths (37);
+        widths <<   4.30000000e-04,   1.43000000e-03,   3.43000000e-03,   5.43000000e-03,
+                    7.43000000e-03,   9.43000000e-03,   1.14300000e-02,   1.34300000e-02,
+                    1.64300000e-02,   2.04300000e-02,   2.54300000e-02,   3.14300000e-02,
+                    3.94300000e-02,   4.94300000e-02,   6.24300000e-02,   7.84300000e-02,
+                    9.84300000e-02,   1.23430000e-01,   1.54430000e-01,   1.94430000e-01,
+                    2.45430000e-01,   3.08430000e-01,   3.89430000e-01,   4.90430000e-01,
+                    6.17430000e-01,   7.78430000e-01,   9.80430000e-01,   1.23543000e+00,
+                    1.55743000e+00,   1.96343000e+00,   2.47443000e+00,   3.11943000e+00,
+                    3.93243000e+00,   4.95743000e+00,   6.24943000e+00,   7.87743000e+00,
+                    9.93143000e+00;
+
+        RxHM->SetWindows(centres, widths, MILLISEC);
+
+    TEMTransmitter* TxLM = TEMTransmitter::New();
+        TxLM->SetRepFrequency( 20, KHZ );
+        VectorXr TimesLM (18);
+        VectorXr AmpsLM (18);
+
+        TimesLM << -8.00000E-004, -7.86965E-004, -7.66493E-004, -7.23688E-004,
+                   -6.39938E-004, -5.16174E-004, -3.93340E-004, -2.63993E-004,
+                   -1.43952E-004, -7.15990E-006, -2.50712E-006,  0.00000E+000,
+                    2.19597E-007,  1.47193E-006,  3.34398E-006,  4.68669E-006,
+                    5.96484E-006,  7.04934E-006;
+        TimesLM.array() = TimesLM.array() + 8e-4;  // Valgrind Hack += yields error. Correct for SkyTEM convention
+
+
+        AmpsLM <<  0.00000E+000,  3.67188E-002,  6.17188E-002,  1.17969E-001,
+                   2.14844E-001,  3.28906E-001,	 4.75781E-001,  6.30469E-001,
+                   7.82031E-001,  9.92969E-001,	 1.00000E+000,  1.00000E+000,
+                   9.63459E-001,  6.01030E-001,  2.29652E-001,  8.64702E-002,
+                   2.53196E-002,  0.00000E+000;
+        AmpsLM.array() = AmpsLM.array() * 8.;   // Also correct for SkyTEM convention in .geo file
+
+        TxLM->SetWaveform( TimesLM, AmpsLM, SEC );
+        // Define wire loop
+        TxLM->SetNumberOfPoints(8);
+        TxLM->SetPoint(0, Vector3r( -16.10,    2.13,  -34));
+        TxLM->SetPoint(1, Vector3r(  -7.51,   10.72,  -34));
+        TxLM->SetPoint(2, Vector3r(   7.51,   10.72,  -34));
+        TxLM->SetPoint(3, Vector3r(  14.92,    3.31,  -34));
+        TxLM->SetPoint(4, Vector3r(  14.92,   -3.31,  -34));
+        TxLM->SetPoint(5, Vector3r(   7.51,  -10.72,  -34));
+        TxLM->SetPoint(6, Vector3r(  -7.51,  -10.72,  -34));
+        TxLM->SetPoint(7, Vector3r( -16.10,   -2.13,  -34));
+        TxLM->SetNumberOfTurns(8);
+
+
+    TEMInductiveReceiver* RxLM = TEMInductiveReceiver::New();
+        RxLM->SetComponent( ZCOMPONENT );      // What about Overloaded to take (X), (X,Y), or (X,Y,Z) or what about tilt?
+        RxLM->SetMoment( 1 );                  // Normalized
+        RxLM->SetReferenceTime( 8e-4, SEC );
+        RxLM->SetRxLocation( (Vector3r() << -16.80, 0, -36.00).finished() );
+
+        // Gate Centres  ms
+        VectorXr centresLM (26);
+            centresLM << 7.150000E-07, 2.215000E-06, 4.215000E-06, 6.215000E-06, 8.215000E-06,
+                         1.021500E-05, 1.221500E-05, 1.471500E-05, 1.821500E-05, 2.271500E-05,
+                         2.821500E-05, 3.521500E-05, 4.421500E-05, 5.571500E-05, 7.021500E-05,
+                         8.821500E-05, 1.107150E-04, 1.387150E-04, 1.742150E-04, 2.197150E-04,
+                         2.767150E-04, 3.487150E-04, 4.397150E-04, 5.537150E-04, 6.977150E-04,
+                         8.792150E-04;
+
+        VectorXr widthsLM (26);
+            widthsLM << 5.700000E-07, 1.570000E-06, 1.570000E-06, 1.570000E-06, 1.570000E-06,
+                        1.570000E-06, 1.570000E-06, 2.570000E-06, 3.570000E-06, 4.570000E-06,
+                        5.570000E-06, 7.570000E-06, 9.570000E-06, 1.257000E-05, 1.557000E-05,
+                        1.957000E-05, 2.457000E-05, 3.057000E-05, 3.957000E-05, 5.057000E-05,
+                        6.257000E-05, 8.057000E-05, 1.005700E-04, 1.265700E-04, 1.605700E-04,
+                        2.015700E-04;
+
+        RxLM->SetWindows(centresLM, widthsLM, SEC);
+
+
+
+    // Specifies survey, this is the Glue Class, the top of the structure, etc.
+    TEMSurvey* Survey = TEMSurvey::New();
+        Survey->SetNumberOfLines(1);  // Flight lines or
+                                      // Internally each line is a class? But that's sort of hidden to the
+                                      // end user. Having each line seperate is nice for constrained inversion, where
+                                      // each line is a nice thing to deal with.
+        Survey->GetLine(0)->SetNumberOfRecords(1);    // Each Record then contains everything needed for modelling response curve(s)
+        Survey->GetLine(0)->GetRecord(0)->SetNumberOfPulseSequences( 1 );
+        Survey->GetLine(0)->GetRecord(0)->SetTransmitterReceiverPair( 0, TxHM, RxHM );
+        //Survey->GetLine(0)->GetRecord(0)->SetTransmitterReceiverPair( 1, TxHM, RxHM );
+
+
+    LayeredEarthEM* Earth = LayeredEarthEM::New();
+        Earth->SetNumberOfLayers(31);
+        Earth->SetLayerThickness( (VectorXr(29) <<  5.00E+00,  5.40E+00,
+             5.80E+00, 6.30E+00, 6.80E+00, 7.30E+00, 7.80E+00, 8.50E+00,
+             9.10E+00, 9.80E+00, 1.06E+01, 1.14E+01, 1.23E+01, 1.33E+01,
+             1.43E+01, 1.54E+01, 1.66E+01, 1.79E+01, 1.93E+01, 2.08E+01,
+             2.25E+01, 2.42E+01, 2.61E+01, 2.81E+01, 3.03E+01, 3.27E+01,
+             3.52E+01, 3.80E+01, 4.10E+01).finished() );
+
+        VectorXcr rho = ( (VectorXcr(31) << 0, 1.54E+01, 4.21E+01,
+            9.25E+01, 1.26E+02, 1.17E+02, 7.57E+01, 3.21E+01, 1.40E+01, 1.96E+01,
+            2.67E+01, 2.56E+01, 2.00E+01, 1.86E+01, 2.31E+01, 2.97E+01, 3.50E+01,
+            3.79E+01, 3.75E+01, 3.32E+01, 2.52E+01, 1.57E+01, 8.38E+00, 5.38E+00,
+            5.49E+00, 6.34E+00, 7.07E+00, 7.81E+00, 8.67E+00, 9.59E+00, 1.05E+01
+        ).finished() ); // 200 Ohm-m half space
+
+        Earth->SetLayerConductivity( 1./rho.array() );
+        // ALL SET UP
+        std::cout << *Survey << std::endl;
+
+        // OK do your thinc
+        //DataTEM* data =
+        Survey->ForwardModel( Earth );
+
+
+    // Clean up
+
+
+
+//     //Survey->Delete();
+//     Earth->Delete();
+//
+//     TxHM->Delete();
+//     RxHM->Delete();
+//     TxLM->Delete();
+//     RxLM->Delete();
+
+}

Modules/LemmaCore/examples/utXMLParse.cpp

@@ -0,0 +1,168 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      04/03/2014 01:08:28 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "Lemma"
+#ifdef TINYXMLSUPPORT
+#include "tinyxml2.h"
+
+#include <ctime>
+using namespace tinyxml2;
+using namespace Lemma;
+
+void ParseTEM( XMLElement* element );
+
+template <typename T>
+std::vector<T> ReadLine( const int& nr, const std::string& line ) ;
+
+int main(int argc, char **argv) {
+
+    if ( argc > 1 ) {
+		XMLDocument* doc = new XMLDocument();
+		clock_t startTime = clock();
+		doc->LoadFile( argv[1] );
+		clock_t loadTime = clock();
+		int errorID = doc->ErrorID();
+
+        // Get some info
+        //static const char* xml = "<element/>";
+        //doc->Parse( xml );
+        //XMLNode* titleElement = doc->FirstChild(); //Element()->FirstChildElement( "PLAY" )->FirstChildElement( "TITLE" );
+        XMLElement* titleElement = doc->FirstChildElement( ); //->FirstChildElement( "TITLE" );
+        //const char* title = titleElement->Value(); //GetText();
+        const char* title = titleElement->Value(); //GetText();
+
+        printf( "Test file '%s' loaded. ErrorID=%d\n", argv[1], errorID );
+        printf( "Name of FirstChildElement: %s\n", title );
+
+        if (strcmp ( titleElement->Value(), "NGXFile" ) != 0) {
+            std::cout <<"NGX file was not detected!" << std::endl;
+            exit(EXIT_FAILURE);
+        } else {
+            //std::cout <<
+            //titleElement->FirstAttribute() ;//"type");// << std::endl;;
+            if (strcmp(  titleElement->FindAttribute("type")->Value(),"tem") == 0) {
+                // Parse TEM file
+                ParseTEM( titleElement );
+            } else {
+                std::cout << "No parser for NGX filetype " << titleElement->FindAttribute("type")->Value() << std::endl;
+                std::cout << "YOU should write one!\n";
+                exit(EXIT_FAILURE);
+            }
+
+        }
+
+		delete doc; doc = 0;
+
+        clock_t deleteTime = clock();
+		if ( !errorID ) {
+			printf( "========================================================\n");
+			printf( "Load time=%u\n",   (unsigned)(loadTime - startTime) );
+			printf( "Delete time=%u\n", (unsigned)(deleteTime - loadTime) );
+			printf( "Total time=%u\n",  (unsigned)(deleteTime - startTime) );
+		}
+
+
+		exit(0);
+	} else {
+        std::cout << "Enter NGX.x file to parse\n";
+    }
+}
+
+void ParseTEM( XMLElement* element ) {
+
+    std::cout << "This is a NGX.t file" << std::endl;
+
+
+    if (strcmp( element->FindAttribute("version")->Value(),"0.1") == 0) {
+        std::cout << "Version 0.1 detected" << std::endl;
+    } else {
+        std::cout << "Unsupported NGX.t version" << std::endl;
+    }
+
+    // system description
+    std::cout << "System description" << std::endl;
+    std::cout << element->FirstChildElement("system")->FirstChildElement("name")->GetText() << std::endl;
+    XMLElement* system = element->FirstChildElement("system");
+        int ntx = atoi(system->FirstChildElement("transmitters")->FindAttribute("num")->Value());
+        std::cout << "numTx = " << ntx << std::endl;
+        // OK Parse Transmitter(s)
+        XMLElement* tx = system->FirstChildElement("transmitters")->FirstChildElement("transmitter");
+        //    std::cout << "Tx id = " << tx->FindAttribute("id")->Value() << std::endl;
+        //std::cout << tx->FirstChildElement("loop.coordinates")->FindAttribute("npoint")->Value();//->FirstChildElement("northingPoints")->Value() << std::endl;
+        //std::cout << tx->FirstChildElement("loop.coordinates")->FirstChildElement("northing.points")->GetText();//->FirstChildElement("northingPoints")->Value() << std::endl;
+        for (int it=0; it<ntx; ++it) {
+            std::cout << "\nTransmitter " << it << std::endl;
+            if (it > 0) tx = tx->NextSiblingElement();
+            int npoints = atoi(tx->FirstChildElement("loop.coordinates")->FindAttribute("npoint")->Value());
+
+            std::cout << "points= " << npoints; // tx->FirstChildElement("loop.coordinates")->FindAttribute("npoint")->Value();//->FirstChildElement("northingPoints")->Value() << std::endl;
+            std::vector<Real> LP = ReadLine<Real>( npoints, std::string(tx->FirstChildElement("loop.coordinates")->FirstChildElement("northing.points")->GetText()) );
+            for (int ii=0; ii<npoints; ++ii) std::cout << "\t"<< LP[ii] ;
+            //std::cout << tx->FirstChildElement("loop.coordinates")->FirstChildElement("northing.points")->GetText();//->FirstChildElement("northingPoints")->Value() << std::endl;
+            //std::cout << tx->FirstChildElement("loop.coordinates")->FirstChildElement("easting.points")->GetText();//->FirstChildElement("northingPoints")->Value() << std::endl;
+        }
+    std::cout << "\nEND system description\n";
+
+}
+
+
+template <typename T>
+std::vector<T> ReadLine( const int& nr, const std::string& line ) {
+
+    std::vector<T>   lineData;
+    std::stringstream  lineStream(line);
+
+    T value;
+    while(lineStream >> value) {
+        lineData.push_back(value);
+    }
+    return lineData;
+}
+#else
+
+int main() {
+    std::cout << "you have to compile lemma with external tinyxml library to use this" << std::endl;
+}
+
+#endif
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ReadInts
+    //--------------------------------------------------------------------------------------
+//     std::vector<int> ReadInts ( const int& nr ) {
+//         std::string buf;
+//         char *dump = new char[BufferSize];
+//         std::vector<int> vals(0);
+//         while (input >> buf) {
+//             if (buf.substr(0, CommentString.size()) == CommentString) {
+//                 input.getline(dump, BufferSize);
+//             } else {
+//                 vals.push_back( atoi(buf.c_str() ));
+//             }
+//             if (static_cast<int>(vals.size()) == nr) {
+//                 delete [] dump;
+//                 return vals;
+//             }
+//
+//         }
+//         delete [] dump;
+//         return vals;
+//     }		// -----  end of method ASCIIParser::ReadInts  -----
+

Modules/LemmaCore/examples/utYAMLconfig.cpp

@@ -0,0 +1,67 @@
+#ifdef HAVE_YAMLCPP
+
+#include "yaml-cpp/yaml.h"
+#include "yaml-cpp/eventhandler.h"
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+struct Params {
+	bool hasFile;
+	std::string fileName;
+};
+
+Params ParseArgs(int argc, char **argv) {
+	Params p;
+
+	std::vector<std::string> args(argv + 1, argv + argc);
+
+	return p;
+}
+
+class NullEventHandler: public YAML::EventHandler
+{
+public:
+	virtual void OnDocumentStart(const YAML::Mark&) {}
+	virtual void OnDocumentEnd() {}
+
+	virtual void OnNull(const YAML::Mark&, YAML::anchor_t) {}
+	virtual void OnAlias(const YAML::Mark&, YAML::anchor_t) {}
+	virtual void OnScalar(const YAML::Mark&, const std::string&, YAML::anchor_t, const std::string&) {}
+
+	virtual void OnSequenceStart(const YAML::Mark&, const std::string&, YAML::anchor_t) {}
+	virtual void OnSequenceEnd() {}
+
+	virtual void OnMapStart(const YAML::Mark&, const std::string&, YAML::anchor_t) {}
+	virtual void OnMapEnd() {}
+};
+
+void parse(std::istream& input)
+{
+	try {
+		YAML::Node doc = YAML::Load(input);
+		std::cout << doc << "\n";
+	} catch(const YAML::Exception& e) {
+		std::cerr << e.what() << "\n";
+	}
+}
+
+int main(int argc, char **argv)
+{
+	Params p = ParseArgs(argc, argv);
+
+	if(argc > 1) {
+		std::ifstream fin;
+		fin.open(argv[1]);
+		parse(fin);
+	} else {
+		parse(std::cin);
+	}
+
+	return 0;
+}
+#else
+int main() {
+}
+#endif
+

Modules/LemmaCore/examples/utgaussquad.cpp

@@ -0,0 +1,84 @@
+// ===========================================================================
+// 
+//       Filename:  utqaussquad.cpp
+// 
+//    Description:  Test for gaussian quadrature algorithm
+// 
+//        Version:  0.0
+//        Created:  10/05/2010 08:54:26 AM
+//       Revision:  none
+//       Compiler:  Tested with g++
+// 
+//         Author:  M. Andy Kass (MAK)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  Broken Spoke Development, LLC
+//
+//          Email:  mkass@numericalgeo.com
+// 
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "Lemma"
+using namespace Lemma;
+using namespace std;
+
+#ifdef LEMMAUSEVTK
+#include "matplot.h"
+using namespace matplot;
+#endif
+
+int main() {
+
+	int N = 600;
+	int upb = 5;
+	int lowb = -1;
+	VectorXr func(N);
+	VectorXr wx(N);
+	Real integrationres,actres;
+
+	gaussianquadrature *lgqw = gaussianquadrature::New();
+	lgqw->SetFreqs(N,upb,lowb);
+	lgqw->CalcAW();
+
+	wx=lgqw->GetAbscissae();
+	func=cos(wx.array());
+//	for (int ii=0;ii<N;++ii) {
+//		cout << wx(ii) << "  " << func(ii) << endl;
+//	}
+
+	lgqw->SetFunc(func);
+	lgqw->Integrate();
+	integrationres=lgqw->GetResult();
+
+	actres = -sin((Real)upb)+sin((Real)lowb);
+	
+	cout << "Gaussian Quadrature Result: " << integrationres << endl;
+	cout << "Actual Result: " << actres << endl;
+	
+		
+	//Plot things
+	#ifdef LEMMAUSEVTK
+	double colour1[3] = {0.0,0.0,1.0};
+
+	Plot2D_VTK p1("X", "Y", 800,600);
+	p1.plot(wx,func,colour1,".-");
+	p1.show();
+
+	#endif
+
+	return EXIT_SUCCESS;
+}
+

Modules/LemmaCore/examples/utgmrprint.cpp

@@ -0,0 +1,258 @@
+// ===========================================================================
+//   This file is distributed with Lemma,
+//
+//       Filename:  utsnmrinversion1d.cpp
+//
+//        Created:  10/07/2010 08:57:04 AM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     10/07/2010
+  @version   0.0
+ **/
+
+#ifdef LEMMAUSEVTK
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkPNGWriter.h"
+#include "vtkRenderLargeImage.h"
+#endif
+
+#include "receiverpoints.h"
+#include "emearth1d.h"
+#include "PolygonalWireAntenna.h"
+
+using namespace Lemma;
+
+std::vector<Real>  readinpfile(const std::string& fname);
+
+int main(int argc, char** argv) {
+
+    std::cout <<
+    "\n"
+    << "hantenna - a programme for computing the h field from polygonal wire\n"
+    << "loop sources \n\n"
+    << "The following copyrights apply to this application:\n"
+    << "Copyright (C) 2009, 2010, 2011, 2012 Colorado School of Mines\n"
+    << "Copyright (C) 2009, 2010, 2011, 2012 Trevor Irons\n"
+    << "Copyright (C) 2011 Broken Spoke Development, LLC\n\n"
+    << "hantenna was built using Lemma (Lemma is an Electromagnetics Modelling API)\n"
+    << "More information may be found at https://lemmasoftware.org\n\n"
+    << "This program is free software: you can redistribute it and/or modify\n"
+    << "it under the terms of the GNU General Public License as published by\n"
+    << "the Free Software Foundation, either version 3 of the License, or\n"
+    << "(at your option) any later version.\n\n"
+
+    << "This program is distributed in the hope that it will be useful,\n"
+    << "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
+    << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
+    << "GNU General Public License for more details.\n\n"
+
+    << "You should have received a copy of the GNU General Public License\n"
+    << "along with this program.  If not, see <http://www.gnu.org/licenses/>.\n\n";
+
+    if (argc < 4) {
+        std::cout << "usage: hantenna.exe  trans.inp cond.inp points.inp \n";
+        exit(0);
+    }
+
+    std::vector<Real> Trans   = readinpfile(std::string(argv[1]));
+    std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
+    std::vector<Real> Points  = readinpfile(std::string(argv[3]));
+
+    //////////////////////////////////////
+    // Define transmitter
+    PolygonalWireAntenna* trans = PolygonalWireAntenna::New();
+        trans->SetNumberOfPoints((int)(Trans[0]));
+        int ip=1;
+        for ( ; ip<=(int)(Trans[0])*2; ip+=2) {
+            trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], -1e-3));
+        }
+ 	    trans->SetNumberOfFrequencies(1);
+ 	    trans->SetFrequency(0, Trans[ip]);
+        trans->SetCurrent(Trans[ip+1]);
+
+	// Receivers
+ 	ReceiverPoints *receivers = ReceiverPoints::New();
+        int nx = (int)Points[0];
+        int ny = (int)Points[1];
+        int nz = (int)Points[2];
+        Real ox = Points[3];
+        Real oy = Points[4];
+        Real oz = Points[5];
+     	Vector3r loc;
+        VectorXr dx(nx-1); // TODO map the dx, dy, dz vectors
+        VectorXr dy(ny-1);
+        VectorXr dz(nz-1);
+        ip = 6;
+        int ir = 0;
+        for ( ; ip <6+nx-1; ++ip) {
+            dx[ir] = Points[ip];
+            ++ir;
+        }
+        ir = 0;
+        for ( ; ip <6+ny-1+nx-1; ++ip) {
+            dy[ir] = Points[ip];
+            ++ir;
+        }
+        ir = 0;
+        for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
+            dz[ir] = Points[ip];
+            ++ir;
+        }
+ 		receivers->SetNumberOfReceivers(nx*ny*nz);
+ 		ir = 0;
+        Real pz  = oz;
+ 		for (int iz=0; iz<nz; ++iz) {
+ 		    Real py    =  oy;
+ 		    for (int iy=0; iy<ny; ++iy) {
+ 		        Real px    =  ox;
+ 		        for (int ix=0; ix<nx; ++ix) {
+ 			        loc << px, py, pz;
+ 			        receivers->SetLocation(ir, loc);
+ 			        if (ix < nx-1) px += dx[ix];
+ 			        ++ ir;
+     	        }
+                if (iy<ny-1) py += dy[iy];
+            }
+            if (iz<nz-1) pz += dz[iz];
+        }
+
+    ////////////////////////////////////
+    // Define model
+    LayeredEarthEM *earth = LayeredEarthEM::New();
+    VectorXcr sigma;
+    VectorXr  thick;
+ 	earth->SetNumberOfLayers(CondMod[0]+1);
+ 	sigma.resize(CondMod[0]+1); sigma(0) = 0; // airlayer
+    thick.resize(CondMod[0]-1);
+    int ilay=1;
+    for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
+        sigma(ilay/2+1) =  1./CondMod[ilay];
+        thick(ilay/2) =  CondMod[ilay+1];
+    }
+    sigma(ilay/2+1) = 1./ CondMod[ilay];
+	earth->SetLayerConductivity(sigma);
+    if (thick.size() > 0) earth->SetLayerThickness(thick);
+
+	EMEarth1D  *EmEarth = EMEarth1D::New();
+		EmEarth->AttachWireAntenna(trans);
+		EmEarth->AttachLayeredEarthEM(earth);
+		EmEarth->AttachReceiverPoints(receivers);
+		EmEarth->SetFieldsToCalculate(H);
+        EmEarth->CalculateWireAntennaFields();
+
+    #ifdef LEMMAUSEVTK
+	vtkRenderer           *renderer = vtkRenderer::New();
+	vtkRenderWindow         *renWin = vtkRenderWindow::New();
+	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+    trans->ApproximateWithElectricDipoles( (VectorXr (3) << 0,0,0).finished() );
+	vtkActor **dipActors = new vtkActor*[trans->GetNumberOfDipoles()];
+	for (int id=0; id<trans->GetNumberOfDipoles(); ++id) {
+		dipActors[id] = trans->GetVtkActor(id);
+		renderer->AddActor(dipActors[id]);
+	}
+    vtkActor *hfield  = receivers->GetVtkGlyphActor (HFIELDREAL, 3, 1e2, 0);
+ 	vtkActor *hfieldi = receivers->GetVtkGlyphActor (HFIELDIMAG, 3, 1e2, 0);
+    hfieldi->GetProperty()->SetColor(1,0,1);
+ 	hfield->GetProperty()->SetColor(0,1,0);
+ 	renderer->AddActor(hfield);
+ 	renderer->AddActor(hfieldi);
+	renderer->SetBackground(0,0,0);
+	// Render the window
+	renWin->AddRenderer(renderer);
+	iren->SetRenderWindow(renWin);
+	iren->Initialize();
+	iren->Start();
+	iren->Render();
+	cout << "Enter File name?: ";
+	std::string pngName;
+	std::cin >> pngName;
+	vtkPNGWriter *pngwrite = vtkPNGWriter::New();
+	vtkRenderLargeImage *renlarge = vtkRenderLargeImage::New();
+	renlarge->SetInput(renderer);
+	renlarge->SetMagnification(2);
+	pngwrite->SetInputConnection(renlarge->GetOutputPort());
+	pngName.append(".png");
+	pngwrite->SetFileName(pngName.c_str());
+	pngwrite->Write();
+	iren->Delete();
+	renWin->Delete();
+	renderer->Delete();
+    #endif
+
+    ////////////////////////////////////
+    // Report
+ 	std::fstream hreal("hfield.dat", std::ios::out);
+    hreal << *trans << std::endl;
+    hreal << *earth << std::endl;
+    hreal << "Right hand coordinate system, z is positive down\n";
+    hreal << "x[m]\ty[m]\tz[m]\tHx[A/m]\tHy[A/m]\tHz[A/m]\n";
+    hreal.precision(8);
+    int i=0;
+	for (int iz=0; iz<nz; ++iz) {
+	for (int iy=0; iy<ny; ++iy) {
+	for (int ix=0; ix<nx; ++ix) {
+        hreal << receivers->GetLocation(i).transpose() << "\t";
+ 		hreal << std::real(receivers->GetHfield(0, i).transpose()[0]) << "\t";
+ 		hreal << std::imag(receivers->GetHfield(0, i).transpose()[0]) << "\t";
+ 		hreal << std::real(receivers->GetHfield(0, i).transpose()[1]) << "\t";
+ 		hreal << std::imag(receivers->GetHfield(0, i).transpose()[1]) << "\t";
+ 		hreal << std::real(receivers->GetHfield(0, i).transpose()[2]) << "\t";
+ 		hreal << std::imag(receivers->GetHfield(0, i).transpose()[2]) << "\n";
+        ++i;
+    }
+    }
+    }
+    hreal.close();
+
+    // Clean up
+    EmEarth->Delete();
+    earth->Delete();
+    receivers->Delete();
+    trans->Delete();
+}
+
+std::vector<Real>  readinpfile(const std::string& fname) {
+    std::string buf;
+    char dump[255];
+    std::vector<Real> vals;
+    std::fstream input(fname.c_str(), std::ios::in);
+    if (input.fail()) {
+        std::cerr << "Input file " << fname << " failed to open\n";
+        exit(EXIT_FAILURE);
+    }
+    while (input >> buf) {
+        if (buf.substr(0,2) == "//") {
+            input.getline(dump, 255);
+        } else {
+            vals.push_back( atof(buf.c_str() ));
+        }
+    }
+    return vals;
+}

Modules/LemmaCore/examples/uthankel2.cpp

@@ -0,0 +1,93 @@
+// ===========================================================================
+//
+//       Filename:  uthankel2.cpp
+//
+//    Description:
+//
+//        Version:  0.0
+//        Created:  05/14/2010 10:26:59 AM
+//       Revision:  none
+//       Compiler:  Tested with g++, icpc, and MSVC 2000
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "hankeltransformhankel2.h"
+#include "hankeltransformgaussianquadrature.h"
+#include "kernelem1d.h"
+#include "layeredearthem.h"
+
+using namespace Lemma;
+
+int main() {
+
+	// Define model
+	LayeredEarthEM *Earth = LayeredEarthEM::New();
+		VectorXcr sigma(3);
+		sigma << 0., .1, 0.01;
+		VectorXr  thick(1);
+		thick << 10;
+		Earth->SetNumberOfLayers(3);
+		Earth->SetLayerConductivity(sigma);
+		Earth->SetLayerThickness(thick);
+
+	// Source
+	// Test with a single dipole
+	DipoleSource *dipole = DipoleSource::New();
+		dipole->SetMoment(2);
+		dipole->SetType(MAGNETICDIPOLE);
+		dipole->SetPolarisation(ZPOLARISATION);
+		dipole->SetLocation(0,0,0);
+		dipole->SetNumberOfFrequencies(1);
+		dipole->SetFrequency(0,1000);
+		dipole->SetPhase(0);
+
+	KernelEm1D *Kernel = KernelEm1D::New();
+		Kernel->Initialise(Earth);
+		Kernel->SetUpSource(dipole, 0);
+		Kernel->SetUpReceiver(-0.4);
+    //std::vector< KernelEm1DBase* > KernelVec;
+    //    KernelVec.push_back( )
+
+	std::cout.precision(10);
+
+	Hankel2 *dig = Hankel2::New();
+//         dig->ComputeRelated(.23 , Kernel);
+//  		std::cout << "  dig->Zgauss "
+//  				  << dig->Zgauss(0, TE, 0, 100.,
+//  								  1000.*std::sqrt(MU0*EPSILON0), Kernel)
+//  			      << std::endl;
+
+// 	HankelTransformGaussianQuadrature *gauss =
+// 			HankelTransformGaussianQuadrature::New();
+// 		std::cout << "gauss->Zgauss "
+// 				  << gauss->Zgauss(0, TE, 0, 100.,
+// 								  1000.*std::sqrt(MU0*EPSILON0), Kernel)
+// 			      << std::endl;
+
+	Earth->Delete();
+	dig->Delete();
+	//gauss->Delete();
+    dipole->Delete();
+	Kernel->Delete();
+	return EXIT_SUCCESS;
+}

Modules/LemmaCore/examples/uthantenna.cpp

@@ -0,0 +1,257 @@
+// ===========================================================================
+//
+//       Filename:  hantenna.cpp
+//
+//        Created:  10/07/2010 08:57:04 AM
+//       Modified:  14 Aug 2012
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//       Copyright (C) 2011    Broken Spoke Development, LLC
+//       Copyright (C) 2012    Trevor Irons
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//
+//          Email:  tirons@mines.edu
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     10/07/2010
+  @version  $Revision: 208 $
+  $Id$
+ **/
+
+#include "receiverpoints.h"
+#include "emearth1d.h"
+#include "PolygonalWireAntenna.h"
+#ifdef LEMMAUSEOMP
+#include "omp.h"
+#endif
+
+#include "timer.h"
+
+using namespace Lemma;
+
+std::vector<Real>  readinpfile(const std::string& fname);
+
+std::vector<std::string>  readinpfile2(const std::string& fname);
+
+int main(int argc, char** argv) {
+
+    std::cout <<
+    "\n"
+    << "hantenna- $Date: 2015-01-26 10:11:24 -0700 (Mon, 26 Jan 2015) $ $Revision: 208 $\n\n"
+    << "hantenna is a programme for computing the h field from polygonal wire\n"
+    << "loop sources \n\n"
+    << "hantenna was built using Lemma (Lemma is an Electromagnetics Modelling API)\n"
+    << "Lemma is Free and Open Source Software (FOSS) and is released under\n"
+    << "the MPL, it is covered by the following copyrights:\n"
+    << "Copyright (C) 2009, 2010, 2011, 2012      Colorado School of Mines\n"
+    << "Copyright (C) 2009, 2010, 2011, 2012      Trevor P. Irons\n"
+    << "Copyright (C) 2011, 2012                  M. Andy Kass\n\n"
+    << "More information may be found at: https://lemmasoftware.org\n"
+    << "                                  info@lemmasoftware.org\n\n"
+    << "=====================================================================\n"
+    << "This programme is part of Lemma, a geophysical modelling and inversion API \n"
+    << "This Source Code Form is subject to the terms of the Mozilla Public\n"
+    << "License, v. 2.0. If a copy of the MPL was not distributed with this\n"
+    << "file, You can obtain one at http://mozilla.org/MPL/2.0/. \n"
+    << "=====================================================================\n\n\n";
+
+    if (argc < 5) {
+        std::cout << "usage: hantenna.exe  trans.inp cond.inp points.inp config.inp \n";
+        exit(0);
+    }
+
+    #ifdef LEMMAUSEOMP
+    std::cout << "OpenMP is using " << omp_get_max_threads() << " threads" << std::endl;
+    #endif
+
+    std::vector<Real> Trans   = readinpfile(std::string(argv[1]));
+    std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
+    std::vector<Real> Points  = readinpfile(std::string(argv[3]));
+    std::vector<std::string> config  = readinpfile2(std::string(argv[4]));
+
+    //////////////////////////////////////
+    // Define transmitter
+    PolygonalWireAntenna* trans = PolygonalWireAntenna::New();
+        trans->SetNumberOfPoints((int)(Trans[0]));
+        int ip=1;
+        for ( ; ip<=(int)(Trans[0])*2; ip+=2) {
+            trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], -1e-3));
+        }
+ 	    trans->SetNumberOfFrequencies(1);
+ 	    trans->SetFrequency(0, Trans[ip]);
+        trans->SetCurrent(Trans[ip+1]);
+        trans->SetMinDipoleRatio(atof(config[1].c_str()));
+        trans->SetMinDipoleMoment(atof(config[2].c_str()));
+        trans->SetMaxDipoleMoment(atof(config[3].c_str()));
+
+	// Receivers
+ 	ReceiverPoints *receivers = ReceiverPoints::New();
+        int nx = (int)Points[0];
+        int ny = (int)Points[1];
+        int nz = (int)Points[2];
+        Real ox = Points[3];
+        Real oy = Points[4];
+        Real oz = Points[5];
+     	Vector3r loc;
+        VectorXr dx(nx-1);
+        VectorXr dy(ny-1);
+        VectorXr dz(nz-1);
+        ip = 6;
+        int ir = 0;
+        for ( ; ip <6+nx-1; ++ip) {
+            dx[ir] = Points[ip];
+            ++ir;
+        }
+        ir = 0;
+        for ( ; ip <6+ny-1+nx-1; ++ip) {
+            dy[ir] = Points[ip];
+            ++ir;
+        }
+        ir = 0;
+        for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
+            dz[ir] = Points[ip];
+            ++ir;
+        }
+ 		receivers->SetNumberOfReceivers(nx*ny*nz);
+ 		ir = 0;
+        Real pz  = oz;
+ 		for (int iz=0; iz<nz; ++iz) {
+ 		    Real py    =  oy;
+ 		    for (int iy=0; iy<ny; ++iy) {
+ 		        Real px    =  ox;
+ 		        for (int ix=0; ix<nx; ++ix) {
+ 			        loc << px, py, pz;
+ 			        receivers->SetLocation(ir, loc);
+ 			        if (ix < nx-1) px += dx[ix];
+ 			        ++ ir;
+     	        }
+                if (iy<ny-1) py += dy[iy];
+            }
+            if (iz<nz-1) pz += dz[iz];
+        }
+
+    ////////////////////////////////////
+    // Define model
+    LayeredEarthEM *earth = LayeredEarthEM::New();
+    VectorXcr sigma;
+    VectorXr  thick;
+ 	earth->SetNumberOfLayers(CondMod[0]+1);
+ 	sigma.resize(CondMod[0]+1); sigma(0) = 0; // airlayer
+    thick.resize(CondMod[0]-1);
+    int ilay=1;
+    for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
+        sigma(ilay/2+1) =  1./CondMod[ilay];
+        thick(ilay/2) =  CondMod[ilay+1];
+    }
+    sigma(ilay/2+1) = 1./ CondMod[ilay];
+	earth->SetLayerConductivity(sigma);
+    if (thick.size() > 0) earth->SetLayerThickness(thick);
+
+	EMEarth1D  *EmEarth = EMEarth1D::New();
+		EmEarth->AttachWireAntenna(trans);
+		EmEarth->AttachLayeredEarthEM(earth);
+		EmEarth->AttachReceiverPoints(receivers);
+		EmEarth->SetFieldsToCalculate(H);
+        EmEarth->SetHankelTransformMethod(string2Enum<HANKELTRANSFORMTYPE>(config[0]));
+        EmEarth->SetHankelTransformMethod(FHTKEY201);
+
+	// Keep track of time
+	jsw_timer timer;
+  	timer.begin();
+    clock_t launch = clock();
+    EmEarth->CalculateWireAntennaFields(true);    // true=status bar
+	Real paTime = timer.end();
+
+    std::cout << "\n\n===========================================\ncalc. real time: " << paTime/60. << "\t[m]\n";
+
+    std::cout << "calc. user time: " <<  (clock()-launch)/CLOCKS_PER_SEC/60.   << "\t[CPU m]"
+		 << std::endl;
+
+    ////////////////////////////////////
+    // Report
+ 	std::fstream hrep("hfield.yaml", std::ios::out);
+    std::fstream hreal("hfield.dat", std::ios::out);
+
+    hrep << *EmEarth << std::endl;
+    hrep.close();
+    //hreal << *trans << std::endl;
+    //hreal << *earth << std::endl;
+
+    hreal << "// Right hand coordinate system, z is positive down\n";
+    hreal << "// x[m]\ty[m]\tz[m]\tHx[A/m]\tHy[A/m]\tHz[A/m]\n";
+    hreal.precision(8);
+    int i=0;
+	for (int iz=0; iz<nz; ++iz) {
+	for (int iy=0; iy<ny; ++iy) {
+	for (int ix=0; ix<nx; ++ix) {
+        hreal << receivers->GetLocation(i).transpose() << "\t";
+ 		hreal << receivers->GetHfield(0, i).transpose() << "\n";
+        ++i;
+    }
+    }
+    }
+    hreal.close();
+    // Clean up
+    EmEarth->Delete();
+    earth->Delete();
+    receivers->Delete();
+    trans->Delete();
+}
+
+std::vector<Real>  readinpfile(const std::string& fname) {
+    std::string buf;
+    char dump[255];
+    std::vector<Real> vals;
+    std::fstream input(fname.c_str(), std::ios::in);
+    if (input.fail()) {
+        std::cerr << "Input file " << fname << " failed to open\n";
+        exit(EXIT_FAILURE);
+    }
+    while (input >> buf) {
+        if (buf.substr(0,2) == "//") {
+            input.getline(dump, 255);
+        } else {
+            vals.push_back( atof(buf.c_str() ));
+        }
+    }
+    return vals;
+}
+
+std::vector<std::string>  readinpfile2(const std::string& fname) {
+    std::string buf;
+    char dump[255];
+    std::vector<std::string> vals;
+    std::fstream input(fname.c_str(), std::ios::in);
+    if (input.fail()) {
+        std::cerr << "Input file " << fname << " failed to open\n";
+        exit(EXIT_FAILURE);
+    }
+    while (input >> buf) {
+        if (buf.substr(0,2) == "//") {
+            input.getline(dump, 255);
+        } else {
+            vals.push_back( std::string(buf.c_str() ));
+        }
+    }
+    return vals;
+}

Modules/LemmaCore/examples/utlogbarriercg.cpp

@@ -0,0 +1,218 @@
+// ===========================================================================
+//
+//       Filename:  utlogbarriercg.cpp
+//
+//        Created:  12/29/2010 09:40:23 AM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     12/29/2010
+  @version   0.0
+ **/
+
+#include "LemmaObject.h"
+//#include "logbarriercg.h"
+//#include "datasnmr.h"
+#include <ctime>
+
+#ifdef HAVEBOOSTRANDOM
+#include <boost/random.hpp>
+#include <boost/random/normal_distribution.hpp>
+#endif
+
+using namespace std;
+using namespace Lemma;
+
+MatrixXr BuildG(const Real& T2, const Real& Larmor, const Real& T,
+        const Real& dt, const VectorXr& nT2Bins, const Real& scale);
+
+MatrixXr BuildGw(const Real& T2, const Real& Larmor, const Real& T,
+        const Real& dt, const VectorXr& nT2Bins, const Real& scale,
+        const VectorXr& freqs);
+
+/// test driver program for log-barrier optimization algorithm
+int main() {
+    return 0;
+}
+
+#ifdef AA
+    // Parameters
+    Real T2(.200);
+    Real Larmor(2000); // Hz
+    Real wL = 2.*PI * Larmor;
+    Real T(.6);
+    Real dt(1e-4);
+    int  nT2(20);
+    Real T2Low(.020);
+    Real T2Hi(.220);
+
+    // T2 bins
+    VectorXr T2Bins = VectorXr::LinSpaced(nT2, T2Low, T2Hi);
+
+    // Build true model
+    VectorXr m = VectorXr::Zero(T2Bins.size());
+    //m[4]  = .15;
+    //m[5]  = .15;
+    //m[6]  = .15;
+    m[7]  = .15;
+    m[8]  = .15;
+    m[11] = .15;
+    m[12] = .15;
+    m[13] = .15;
+    m[14] = .15;
+
+    // Build G
+    srand ( std::time(NULL) );
+    Real scale = 1e-4*std::abs(VectorXr::Random(1)[0]);
+    MatrixXr G = BuildG(T2, wL, T, dt, T2Bins, scale);
+
+    #ifdef HAVEBOOSTRANDOM
+    Real mean = 0;
+    Real variance = scale*.05*(.5+.15);
+    boost::mt19937 randgen(static_cast<unsigned int> (std::time(0)));
+    boost::normal_distribution<Real> noise(mean, variance);
+    boost::variate_generator<boost::mt19937, boost::normal_distribution<Real> > nD(randgen, noise);
+    VectorXr noisy = G*m;
+//     for (int i=0; i < noisy.size(); ++i) {
+//        noisy(i) += nD();
+//     }
+    #else
+    VectorXr d = G * m;
+    VectorXr noisy = d.array();// + .05*( .5 + .15 ) *Eigen::ArrayXd::Random(d.size());
+    #endif
+
+    // Clean Data, and 10% random noise data
+    VectorXr RefMod = VectorXr::Zero(m.size());
+
+    // todo let alpha's be set.
+    Real alpha_s  =  .0;
+    Real alpha_t2 = 1.0;
+    Eigen::SparseMatrix<Real> Wd (G.rows(), G.rows());
+    Eigen::SparseMatrix<Real> Wm (G.cols(), G.cols());
+
+    /////////////////
+    // Model objective function
+    int ii = 0;
+    for (int iT2=0; iT2<nT2; ++iT2) {
+        Wm.coeffRef(ii, ii) += alpha_s - alpha_t2;
+        if (iT2 < nT2-1) {
+            Wm.coeffRef(ii, ii+1) += alpha_t2;
+        }
+        ++ii;
+    }
+    Wm.finalize();
+    std::fstream  WmOut( "wmMatrix.dat" , std::ios::out);
+    WmOut << Wm << std::endl;
+    WmOut.close();
+
+    // Data objective function
+    /// @todo add measure of suspected std dev to this.
+    for (int id=0; id<G.rows(); ++id) {
+        Wd.coeffRef(id, id) = 1.; // 1/sigma(i);
+    }
+    Wd.finalize();
+
+    // invert time-domain sinusoid exponent
+    //LogBarrierCG(G, RefMod, noisy, (Real)0., (Real)1e5, nT2, Wd, Wm, scale*1e1);
+    return EXIT_SUCCESS;
+
+//     #if 0
+//     // invert FD
+//     Eigen::FFT<Real> fft;
+//     fft.SetFlag(fft.HalfSpectrum);
+//     VectorXcr dw;
+//     fft.fwd(dw, noisy);
+//
+//     // Compute freqs Hz
+//     VectorXr freqs = DataSNMR::GetNMRSampledFrequencies(T/dt, dt);
+//     Real df = freqs[1];
+//     Real window = 300; // Hz
+//     int iL = (int)(Larmor/df);
+//     MatrixXr Gw = BuildGw(T2, wL, T, dt, T2Bins, scale, freqs.segment(iL-window/2, window));
+//
+//     VectorXr dwu = VectorXr::Zero(2.*window);
+//     dwu.head(window) = dw.segment(iL-window/2, window).real();
+//     dwu.tail(window) = dw.segment(iL-window/2, window).imag();
+//
+//     Eigen::DynamicSparseMatrix<Real> Wdw (Gw.rows(), Gw.rows());
+//     for (int id=0; id<Gw.rows(); ++id) {
+//         Wdw.coeffRef(id, id) = 1.; // 1/sigma(i);
+//     }
+//     Wdw.finalize();
+//     std::cout << "psize " << (Gw*RefMod).size() << std::endl;
+//     VectorXr dw2 = Gw*m;
+//     //std::cout << dw2 << std::endl;
+//     std::cout << "Gw " << Gw.rows() << "\t" << Gw.cols() << std::endl;
+//     std::cout << "dwu " << dwu.size() << std::endl;
+//     //RefMod.array() += .1;
+//     VectorXr mod = LogBarrierCG(Gw, RefMod, dwu, 0., 3., nT2, Wdw, Wm,scale*1e5);
+//     //VectorXr mod = LogBarrierCG(Gw, dwu, 0., 1e5);
+//
+//     std::fstream Gws("gw.out", std::ios::out);
+//     Gws << Gw << std::endl;
+//     Gws.close();
+//     #endif
+}
+
+#endif
+
+MatrixXr BuildG(const Real& T2, const Real& wL, const Real& T,
+            const Real& dt, const VectorXr& T2Bins, const Real& scale) {
+    std::cout << "scale " << scale << std::endl;
+    int nt = T/dt;
+    int nT2 = T2Bins.size();
+    MatrixXr G = MatrixXr::Zero(nt,nT2);
+    for (int iT2=0; iT2<nT2; ++iT2) {
+        Real t=0;
+        for (int it=0; it<nt; ++it) {
+            G(it, iT2) = scale * exp(-t/T2Bins(iT2)) * sin(wL * t);
+            t += dt;
+        }
+    }
+    return G;
+}
+
+
+MatrixXr BuildGw(const Real& T2, const Real& wL, const Real& T,
+            const Real& dt, const VectorXr& T2Bins, const Real& scale,
+            const VectorXr &freqs) {
+    VectorXr omegas = freqs.array() * 2.*PI;
+    std::cout << "scale " << scale << std::endl;
+    int nT2 = T2Bins.size();
+    MatrixXr Gw = MatrixXr::Zero(2*freqs.size(), nT2);
+    int nW = freqs.size();
+    Real sc = 1./dt;
+    for (int iT2=0; iT2<nT2; ++iT2) {
+        for (int iw=0; iw<freqs.size(); ++iw) {
+            Real a = 1./T2Bins(iT2);
+            Complex T = wL / ( wL*wL + (a+Complex(0, omegas[iw]))*(a+Complex(0,omegas[iw])) ) ;
+            Gw(iw, iT2)     = sc * scale * std::real(T);
+            Gw(iw+nW, iT2) = sc * scale * std::imag(T);
+        }
+    }
+    return Gw;
+}

Modules/LemmaCore/examples/utpolygonalantenna.cpp

@@ -0,0 +1,67 @@
+// ===========================================================================
+//
+//       Filename:  utpolygonalantennae.cpp
+//
+//    Description:  Unit test of PolygonalAntenna
+//
+//        Version:  0.0
+//        Created:  05/26/2010 04:30:56 PM
+//       Revision:  none
+//       Compiler:  Tested with g++, icpc, and MSVC 2000
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "Lemma"
+
+using namespace Lemma;
+
+int main() {
+
+    PolygonalWireAntenna *pa = PolygonalWireAntenna::New();
+        pa->SetNumberOfPoints(5);
+        pa->SetPoint(0, Vector3r(  0,  0, -1e-3));
+        pa->SetPoint(1, Vector3r( 10,  0, -1e-3));
+        pa->SetPoint(2, Vector3r( 10, 10, -1e-3));
+        pa->SetPoint(3, Vector3r(  0, 10, -1e-3));
+        pa->SetPoint(4, Vector3r(  0,  0, -1e-3));
+        pa->SetNumberOfFrequencies(1);
+        pa->SetFrequency(0,1000.);
+        pa->SetCurrent(10.);
+        pa->SetNumberOfTurns(1);
+
+    Vector3r pos;
+    pos << 23,23,23;
+    pa->ApproximateWithElectricDipoles(pos);
+    std::cout << "pos "<< pos.transpose() << "  ndipoles="
+              << pa->GetNumberOfDipoles() << std::endl;
+    pos << 2,2,2;
+    pa->ApproximateWithElectricDipoles(pos);
+    std::cout << "pos "<< pos.transpose() << "  ndipoles="
+              << pa->GetNumberOfDipoles() << std::endl;
+    pos << 1,1,1;
+    pa->ApproximateWithElectricDipoles(pos);
+    std::cout << "pos "<< pos.transpose() << "  ndipoles="
+              << pa->GetNumberOfDipoles() << std::endl;
+
+    pa->Delete();
+}

Modules/LemmaCore/examples/utreadfemubc.cpp

@@ -0,0 +1,68 @@
+// ===========================================================================
+// 
+//       Filename:  utreadfemubc.cpp
+// 
+//    Description:  
+// 
+//        Version:  0.0
+//        Created:  01/03/2013 02:47:24 PM
+//       Revision:  none
+//       Compiler:  Tested with g++ 
+// 
+//         Author:  M. Andy Kass (MAK)
+//
+//   Organisation:  Broken Spoke Development, LLC
+//                  
+//
+//          Email:  mkass@numericalgeo.com
+// 
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "Lemma"
+
+using namespace Lemma;
+
+#ifdef LEMMAUSEVTK
+#include "matplot.h"
+using namespace matplot;
+#endif
+
+int main() {
+
+	std::string datfile;
+
+	datfile = "126.obs";
+	
+	DataReaderFemUBC* Reader = DataReaderFemUBC::New();
+	DataFEM* thedata = DataFEM::New();
+
+	Reader->SetDataFEM(thedata);
+
+	try {
+	Reader->ReadData(datfile,1);
+	} catch(std::exception& e) {
+		exit(EXIT_FAILURE);
+	}
+
+	std::cout << *Reader << std::endl;
+	std::cout << *thedata << std::endl;
+
+
+	thedata->Delete();
+	Reader->Delete();
+
+	return EXIT_SUCCESS;
+}

Modules/LemmaCore/examples/utsandbox.cpp

@@ -0,0 +1,54 @@
+// ===========================================================================
+//
+//       Filename:  utsandbox.cpp
+//
+//        Created:  07/29/2010 10:34:23 AM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file test sandbox helm1rameters
+  @author   Trevor Irons
+  @date     07/29/2010
+  @version   0.0
+ **/
+
+#include "Lemma"
+
+using namespace Lemma;
+
+int main() {
+
+    PolygonalWireAntenna *helm1 = PolygonalWireAntenna::New();
+	helm1->SetNumberOfFrequencies(1);
+	helm1->SetFrequency(0, 1000);
+	helm1->SetNumberOfPoints(5);
+	helm1->SetPoint(0, Vector3r(   0,   0, 0));
+	helm1->SetPoint(1, Vector3r(   0,   0, 0));
+	helm1->SetPoint(2, Vector3r(   0,   0, 0));
+	helm1->SetPoint(3, Vector3r(   0,   0, 0));
+	helm1->SetPoint(4, Vector3r(   0,   0, 0));
+	helm1->SetCurrent(1.);
+    helm1->Delete();
+}

Modules/LemmaCore/examples/uttemforward.cpp

@@ -0,0 +1,199 @@
+// ===========================================================================
+//
+//       Filename:  uttemforward.cpp
+//
+//    Description:  TEM Forward Modeller
+//
+//        Version:  0.0
+//        Created:  02/10/2011 03:32:18 PM
+//       Revision:  none
+//       Compiler:  Tested with g++, icpc, and MSVC 2000
+//
+//         Author:  M. Andy Kass (MAK)
+//
+//		Copyright:  2011 Trevor Irons and M. Andy Kass
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  Broken Spoke Development, LLC
+//
+//          Email:  mkass@numericalgeo.com
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include <time.h>
+#include "Lemma"
+#include "banner.h"
+
+using namespace Lemma;
+
+#ifdef LEMMAUSEVTK
+#include "matplot.h"
+using namespace matplot;
+#endif
+
+int main(int argv, char** argc) {
+
+	// Just for timing
+	clock_t launch =clock();
+
+	// Banner display
+	std::string name;
+	std::string version;
+	std::string usage;
+	name = "TEM Forward Modeller - 1D";
+	version = "1.0beta";
+	usage = "temforward1d [inputfile]";
+	banner(name,version,usage);
+    if ( argv < 2 ) {
+        std::cout << "No input file specified, should be of format\n";
+        std::cout <<
+"origparabk.obsy        ! Observations-style file.\n\
+origpara.con           ! Starting model.\n\
+100                    ! Number of Hankel kernel evaluations.\n\
+100 1. 1.0E+10         ! Number of Fourier kernel evaluations, & min & max freq.\n\
+N                      ! Is noise to be added?\n\
+5.  1.0E-08  30        ! Noise to be added: percentage, threshold, seed.\n\
+" << std::endl;
+
+        exit(0);
+    }
+
+	time_t curr=time(0);
+	std::cout << std::endl;
+	std::cout << "  Start time: " << ctime(&curr) << std::endl;
+
+	// Define objects
+	PolygonalWireAntenna* Trans = PolygonalWireAntenna::New();
+//     PolygonalWireAntenna* pa = PolygonalWireAntenna::New();
+//         pa->SetNumberOfPoints(5);
+// 	    pa->SetPoint(0, Vector3r(   0,   0, -30.001));
+// 	    pa->SetPoint(1, Vector3r(   0,  20, -35.001));
+// 	    pa->SetPoint(2, Vector3r(  20 , 20, -35.001));
+// 	    pa->SetPoint(3, Vector3r(  20,   0, -30.001));
+// 	    pa->SetPoint(4, Vector3r(   0,   0, -30.001));
+// 	    pa->SetCurrent(100.);
+//         pa->SetNumberOfTurns(1);
+
+
+	//DipoleSource* Trans = DipoleSource::New();
+	ReceiverPoints *Receivers = ReceiverPoints::New();
+	LayeredEarthEM *Earth = LayeredEarthEM::New();
+	VectorXr maintimes;
+
+    /////////////////////////////////////
+	// Using model reader
+	ModelReaderTem1DUBC* Reader = ModelReaderTem1DUBC::New();
+	    Reader->SetEMEarth1D(Earth);
+	    Reader->SetTransmitter(Trans);
+	    Reader->SetReceiver(Receivers);
+	    Reader->ReadParameters( argc[1] );
+	    maintimes = Reader->GetTimes();
+
+	// Attach it to instrumentTEM
+	InstrumentTem *instrument = InstrumentTem::New();
+	    instrument->EMEarthModel(Earth);
+	    instrument->SetTransmitLoop(Trans);
+	    //instrument->SetTransmitLoop(pa);
+	    //instrument->SetDipoleSource(Trans);
+	    instrument->SetReceiver(Receivers);
+
+
+
+        // TODO need to input these
+        VectorXr widths = VectorXr::Ones(maintimes.size()); // * 1e-5;
+        /*
+        widths.segment(0,5).array()  *= 5e-6;
+        widths.segment(5,5).array()  *= 1e-5;
+        widths.segment(10,5).array() *= 5e-5;
+        widths.segment(15,5).array() *= 2e-4;
+        widths.segment(20,5).array() *= 4e-4;
+        //widths.segment(25,5).array() *= 4e-4;
+        */
+        /*
+        widths << .48,4.,8.,8.,8.,8.,8.,12.,12.,12.,16.,16.,16.,20.,24.,24.,28.,32.,32.,36.,40.,44.,52.,
+  				56.,64.,68.,76.,88.,96.,104.,120.,132.,144.,160.,180.,200.,220.,248.,272.,304.,336.,376.,
+  				416.,464.,516.,572.,636.,708.,784.,872.,968.,1076.,1196.,1328.,1476.,1640.,1824.,2024.,
+  				2248.,2500. ;
+        */
+        widths *= 1e-6;
+        std::cout << widths.transpose() << std::endl;
+
+        instrument->SetReferenceTime(.025);
+        //instrument->SetReferenceTime(.025);
+	    instrument->SetTimeGates(maintimes, widths);
+
+        //instrument->SetReceiverType(MAGNETOMETER);
+        instrument->SetReceiverType(INDUCTIVE);
+        // Set the pulse
+        VectorXr Amp(4); Amp << 0,1,1,0;
+        VectorXr Times(4); Times << 0., .0005, .02497, .025;
+
+        //VectorXr Amp(4); Amp << 0,1,1,0;
+        //VectorXr Times(4); Times << -0.025, -0.0245, -0.00003, 0.0;
+        //instrument->SetPulse(Amp, Times);
+
+        //VectorXr Amp(3); Amp << 0,1,0;
+        //VectorXr Times(3); Times << 0, .015, .030;
+        //instrument->SetPulse(Amp, Times);
+
+        //VectorXr Amp(3); Amp << 0,1,0;
+        //VectorXr Times(3); Times << 0, .01, .020;
+        instrument->SetPulse(Amp, Times);
+
+    std::cout << *Earth << "\n" << *Trans << "\n" << *Receivers << std::endl;
+
+	// Perform the forward model
+	//instrument->MakeDirectCalculation( FHTKEY201  );
+	//instrument->MakeLaggedCalculation( ANDERSON801 );
+	//instrument->MakeLaggedCalculation( QWEKEY );
+	//instrument->MakeLaggedCalculation( CHAVE );
+	instrument->MakeLaggedCalculation( FHTKEY201 );
+	//instrument->MakeLaggedCalculation( FHTKEY101 );
+
+    int nlag = instrument->GetMeasurements().rows();
+	// Output results to screen
+	for (int ii=0; ii<nlag; ii++) {
+		std::cout<<"  "<<instrument->GetMeasurements()(ii,0)<<"  "<<instrument->GetMeasurements()(ii,1)<<std::endl;
+	}
+
+	// Output results to file
+	std::ofstream outfile1;
+	outfile1.open("solution.out");
+    //outfile1 << "//timegate [ms]\t dB/dt [nT/s]" << std::endl;
+    outfile1 << "//timegate [s]\t dB/dt [T/s]" << std::endl;
+	for (int ii=0; ii<nlag; ii++) {
+		//outfile1 << 1e3*instrument->GetMeasurements()(ii,0)<< "\t" << (MU0*1e9)*instrument->GetMeasurements()(ii,1)<<std::endl;
+		outfile1 << instrument->GetMeasurements()(ii,0)<< "\t" << MU0*instrument->GetMeasurements()(ii,1) << std::endl;
+		//outfile1 << instrument->GetMeasurements()(ii,0)<< "\t" << instrument->GetMeasurements()(ii,1)<<std::endl;
+	}
+	outfile1.close();
+
+
+	// Timing stuff
+	std::cout << "  Time for execution: " <<  (clock()-launch)/CLOCKS_PER_SEC   << " [CPU] seconds."
+		 << std::endl;
+
+
+	instrument->Delete();
+	Trans->Delete();
+	Earth->Delete();
+	Receivers->Delete();
+	Reader->Delete();
+
+	return EXIT_SUCCESS;
+}
+
+

Modules/LemmaCore/examples/utteminv1d.cpp

@@ -0,0 +1,131 @@
+// ===========================================================================
+//
+//       Filename:  utteminv1d.cpp
+//
+//    Description:  Invert for the difference in water table
+//
+//        Version:  0.0
+//        Created:  03/23/2011 02:07:19 PM
+//       Revision:  none
+//       Compiler:  Tested with g++, icpc, and MSVC 2000
+//
+//         Author:  M. Andy Kass (MAK)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  Broken Spoke Development, LLC
+//
+//          Email:  mkass@numericalgeo.com
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "Lemma"
+#include "banner.h"
+
+using namespace Lemma;
+
+int main () {
+
+	VectorXr maintimes;
+	int ntimes;
+
+	//Banner
+	std::string name;
+	std::string version;
+	std::string usage;
+	name = "TEM Inversion - 1D";
+	version = "1.0beta";
+	usage = "teminv1d [inputfile]";
+	banner(name,version,usage);
+
+	PolygonalWireAntenna* Trans = PolygonalWireAntenna::New();
+	ReceiverPoints* Receivers = ReceiverPoints::New();
+	LayeredEarthEM* StartingMod = LayeredEarthEM::New();
+	//LayeredEarthEM* RecoveredMod = LayeredEarthEM::New();
+	InverseSolverTEM1D* Inverse = InverseSolverTEM1D::New();
+
+
+	Trans->SetNumberOfPoints(5);
+    Trans->SetPoint(0, Vector3r(   0,   0, -1e-3));
+    Trans->SetPoint(1, Vector3r( 100,   0, -1e-3));
+    Trans->SetPoint(2, Vector3r( 100, 100, -1e-3));
+    Trans->SetPoint(3, Vector3r(   0, 100, -1e-3));
+    Trans->SetPoint(4, Vector3r(   0,   0, -1e-3));
+	Trans->SetCurrent(1);
+	Trans->SetNumberOfTurns(1);
+    Trans->SetMinDipoleRatio(1/2.);
+
+	Vector3r loc;
+	Real ox = 50.;
+	Real oy = 50.;
+	Real depth = -1.e-3;
+	Receivers->SetNumberOfReceivers(1);
+	loc << ox,oy,depth;
+	Receivers->SetLocation(0,loc);
+
+	StartingMod->SetNumberOfLayers(22);
+	StartingMod->SetLayerConductivity( (VectorXcr(22) << 0.,1.e-3,
+		1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,
+		1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3).finished() );
+	StartingMod->SetLayerThickness((VectorXr(20)<<4,4,4,4,4,4,4,4,4,4,4,4,4,
+		4,4,4,4,4,4,4).finished());
+	//EM-47 gate times
+	ntimes = 30;
+	maintimes.resize(ntimes);
+	maintimes << 36.,45.25,57.,72.25,92.,117.,148.,186.5,234.,290.,
+		352.5,427.5,525.,647.5,802.5,1002.5,1257.5,1582.5,1997.5,2525.,
+		3197.5,4055.,5147.5,6542.5,8322.5,10592.,13490.,17187.,21902.,
+		27915.;
+	maintimes = maintimes.array() /1000000;
+
+	Inverse->AttachStartMod(StartingMod);
+	Inverse->AttachWireAntenna(Trans);
+	Inverse->AttachReceiver(Receivers);
+	Inverse->SetTimes(maintimes);
+	LayeredEarthEM* RecoveredMod = StartingMod->Clone();
+	Inverse->AttachRecMod(RecoveredMod);
+
+	//Read input data - temp for SEG!
+	/// TODO write a TEM data reader class.  and a functional TEM data class
+	std::fstream infile("rotateddat.dat",std::ios::in);
+	if (infile.fail()) {
+		std::cout << "Shit!" << std::endl;
+	}
+	MatrixXr rotinpdata;
+	rotinpdata.resize(30,1);
+	for (int ii=0;ii<30;ii++) {
+		infile >> rotinpdata(ii);
+	}
+	infile.close();
+
+
+	Inverse->AttachMeasuredData(rotinpdata);
+	//Inverse->SetFreeParams( (VectorXi(2)<<0,1).finished(),
+	//	(VectorXi(3)<<1,2,3).finished());
+
+	std::cout << *Inverse << std::endl;
+
+	Inverse->Calculate();
+	//Inverse->ShowSoln();
+	//Inverse->WriteSoln("temp.dat");
+
+	Inverse->Delete();
+	RecoveredMod->Delete();
+	StartingMod->Delete();
+	Receivers->Delete();
+	Trans->Delete();
+
+	return EXIT_SUCCESS;
+}

Modules/LemmaCore/examples/uttemreader.cpp

@@ -0,0 +1,84 @@
+// ===========================================================================
+// 
+//       Filename:  uttemreader.cpp
+// 
+//    Description:  
+// 
+//        Version:  0.0
+//        Created:  02/19/2011 12:30:59 PM
+//       Revision:  none
+//       Compiler:  Tested with g++, icpc, and MSVC 2000 
+// 
+//         Author:  M. Andy Kass (MAK)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  Broken Spoke Brewery, LLC
+//
+//          Email:  mkass@numericalgeo.com
+// 
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include "Lemma"
+
+using namespace Lemma;
+
+
+int main() {
+
+	/*
+	ModelReaderTem1DUBC* Reader = ModelReaderTem1DUBC::New();
+	PolygonalWireAntenna* Trans = PolygonalWireAntenna::New();
+	ReceiverPoints *Receivers = ReceiverPoints::New();
+	LayeredEarthEM *Earth = LayeredEarthEM::New();
+	Reader->SetEMEarth1D(Earth);
+	Reader->SetTransmitter(Trans);
+	Reader->SetReceiver(Receivers);
+	Reader->ReadParameters("em1dtmfwd.in");
+
+	Reader->Delete();
+	Earth->Delete();
+	Receivers->Delete();
+	Trans->Delete();
+	*/
+
+	int rnObs;
+	int rnGates;
+	std::string inputfile;
+	MatrixXr somedata;
+	Vector3Xr somepositions;
+
+	DataTEM* TheData = DataTEM::New();
+	DataReaderTem* Reader = DataReaderTem::New();
+
+	inputfile = "temreadertest.txt";
+
+	Reader->SetDataTEM(TheData);
+	Reader->ReadData(inputfile);
+
+	rnObs = TheData->GetnObs();
+	rnGates = TheData->GetnGates();
+	somedata = TheData->GetData();
+	somepositions = TheData->GetPositions();
+
+	std::cout << rnObs << " " << rnGates << std::endl;
+	std::cout << somedata << std::endl;
+	std::cout << somepositions << std::endl;
+
+	Reader->Delete();
+	TheData->Delete();
+
+	return EXIT_SUCCESS;
+}

Modules/LemmaCore/examples/uttemsandbox.cpp

@@ -0,0 +1,207 @@
+// ===========================================================================
+//
+//       Filename:  uttemsandbox.cpp
+//
+//        Created:  09/29/2010 07:38:26 AM
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti), M. Andy Kass
+//
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//					Broken Spoke Development, LLC
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//					mkass@numericalgeo.com
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     09/29/2010
+  @version   0.0
+ **/
+
+#include "Lemma"
+using namespace Lemma;
+
+#ifdef LEMMAUSEVTK
+#include "matplot.h"
+using namespace matplot;
+#endif
+
+int main() {
+
+
+    int nfreq  =    6000;  // Number of frequencies
+
+	// Specify receiver times
+	int ntimes	=	30;
+	VectorXr times(ntimes);
+	times <<
+	36.0, 45.25, 57.0, 72.25, 92.0,
+	117.0, 148.0, 186.5, 234.0, 290.0,
+	352.5, 427.5, 525.0, 647.5, 802.5,
+	1002.5, 1257.5, 1582.5, 1997.5, 2525.0,
+	3197.5, 4055.0, 5147.5, 6542.5, 8322.5,
+	10592.0, 13490.0, 17187.0, 21902.0, 27915.0;
+
+	times = times.array()*1.0e-6;
+
+	//CALCULATE ABSCISSA
+	//define bounds
+
+	gaussianquadrature *lgqw = gaussianquadrature::New();
+
+	Real lowb;
+	Real upb;
+	lowb = 1.0;
+	upb = 6000000.0;
+	lgqw->SetFreqs(nfreq,upb,lowb);
+	lgqw->CalcAW();
+
+	VectorXr xu(nfreq);
+	xu=lgqw->GetAbscissae();
+	DipoleSource* Trans=DipoleSource::New();
+	Trans->SetType(MAGNETICDIPOLE);
+	Trans->SetPolarisation(ZPOLARISATION);
+	Trans->SetLocation(0.0,0.0,-1e-4);
+
+
+    // Specify Transmitter
+    //PolygonalWireAntenna* Trans  = PolygonalWireAntenna::New();
+    //Trans->SetNumberOfPoints(5);
+    //Trans->SetPoint(0, Vector3r(   0,   0, -1e-3));
+    //Trans->SetPoint(1, Vector3r( 100,   0, -1e-3));
+    //Trans->SetPoint(2, Vector3r( 100, 100, -1e-3));
+    //Trans->SetPoint(3, Vector3r(   0, 100, -1e-3));
+    //Trans->SetPoint(4, Vector3r(   0,   0, -1e-3));
+    Trans->SetNumberOfFrequencies(nfreq);
+    VectorXr  f(nfreq);
+    for (int ifreq=0; ifreq<nfreq; ++ifreq) {
+    //    Trans->SetFrequency(ifreq, 1e-3 + dfreq*(Real)(ifreq));
+		Trans->SetFrequency(ifreq,xu(ifreq));
+        f(ifreq) = xu(ifreq);
+    }
+	//cout << Trans->GetFrequencies();
+    //Trans->SetCurrent(1);
+    //Trans->SetNumberOfTurns(1);
+
+    // Earth properties, top layer is air layer 0 conductivity is fine
+    LayeredEarthEM *Earth = LayeredEarthEM::New();
+        Earth->SetNumberOfLayers(4);
+        Earth->SetLayerConductivity( (VectorXcr(4) << 0.,1.e-6,1.e-2,1.e-6 ).finished() );
+        Earth->SetLayerThickness( (VectorXr(2) << 50, 20).finished() );
+
+    // Receivers, just 1 in the centre for now
+   	ReceiverPoints *Receivers = ReceiverPoints::New();
+		Vector3r loc;
+		Real ox     =    50.;
+		Real oy     =    50.;
+		Real depth  =  -1e-2;
+		Receivers->SetNumberOfReceivers(1);
+		loc << ox, oy, depth;
+		Receivers->SetLocation(0, loc);
+
+	// EmEarth
+	EMEarth1D  *EmEarth = EMEarth1D::New();
+		//EmEarth->AttachWireAntenna(Trans);
+		EmEarth->AttachDipoleSource(Trans);
+		EmEarth->AttachLayeredEarthEM(Earth);
+		EmEarth->AttachReceiverPoints(Receivers);
+		EmEarth->SetFieldsToCalculate(H);
+        // slower but may be more accurate, depending on frequencies
+        EmEarth->SetHankelTransformMethod(CHAVE);
+        //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
+
+    // Do calculation
+   	//EmEarth->CalculateWireAntennaFields();
+	//EmEarth->MakeCalc();
+	EmEarth->MakeCalc3();
+
+    // Grab Z component
+    // A little painful, I'll try to clean up the API
+    VectorXcr Hw(nfreq);
+    VectorXr  e(nfreq);
+    Hw(0) = 0.; // DC component
+    for (int ifreq=0; ifreq<nfreq; ++ifreq) {
+        Hw(ifreq) = (Receivers->GetHfield(ifreq, 0))(2);
+        e(ifreq) = ifreq;
+    }
+
+    VectorXr Ht(ntimes);                   // Time domain H field
+
+	//Sine transform
+	VectorXr func(nfreq);
+	for (int ii=0;ii<ntimes;++ii) {
+		func = (Hw.imag().array()*((f.array()*times(ii)).sin()));
+		lgqw->SetFunc(func);
+		lgqw->Integrate();
+		Ht(ii)=lgqw->GetResult();
+	}
+	Ht = Ht.array()*(-2.0/PI);
+
+	std::ofstream myfile1;
+	myfile1.open("solution.txt");
+	for (int ii=0;ii<ntimes;++ii) {
+		myfile1 << times(ii) << "   " << Ht(ii) << std::endl;
+	}
+	myfile1.close();
+
+	std::ofstream myfile2;
+	myfile2.open("freq_domain.txt");
+	for (int ii=0;ii<nfreq;++ii) {
+		myfile2 << f(ii) << "  " << Hw(ii).real() << "  " << Hw(ii).imag()
+			<< std::endl;
+	}
+	myfile2.close();
+
+    //std::cout << Ht << std::endl;
+    //std::cout << Hw.real() << std::endl;
+
+    // Quick and dirty plot
+    #ifdef  LEMMAUSEVTK
+    double colour1[3] = {0.0,0.0,1.0};
+    double colour2[3] = {1.0,0.0,0.0};
+
+    Plot2D_VTK p1("Hz", "Re(Hw)", 800, 600);
+    p1.plot(f, Hw.real().eval(), colour1, ".-");
+    p1.show();
+
+    Plot2D_VTK p2("Hz", "Im(Hw)", 800, 600);
+    p2.plot(f, Hw.imag().eval(), colour2, ".-");
+    p2.show();
+
+    Plot2D_VTK p3("Time (seconds)", "Ht", 800, 600);
+    p3.plot(times, Ht, colour2, ".-");
+    p3.show();
+
+	//Plot2D_VTK p4("index", "frequency",800,600);
+	//p4.plot(e,f,colour1,".-");
+	//p4.show();
+
+    #endif
+
+    // Clean up
+    Trans->Delete();
+    Earth->Delete();
+    Receivers->Delete();
+    EmEarth->Delete();
+
+    return EXIT_SUCCESS;
+}

Modules/LemmaCore/examples/wireantenna.cpp

@@ -0,0 +1,319 @@
+// ===========================================================================
+//
+//       Filename:  utdipolesource.cpp
+//
+//    Description:
+//
+//        Version:  0.0
+//        Created:  12/02/2009 11:57:14 AM
+//       Revision:  none
+//       Compiler:  g++ (c++)
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+#include <iostream>
+#include <fstream>
+
+#include "dipolesource.h"
+#include "layeredearth.h"
+#include "receiverpoints.h"
+#include "emearth1d.h"
+#include "WireAntenna.h"
+#include "PolygonalWireAntenna.h"
+
+#if  LEMMAUSEVTK
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderLargeImage.h"
+#include "vtkPNGWriter.h"
+#endif     // -----  not LEMMA_USE_VTK  -----
+
+// For testing purposes disable VTK and run scale.sh
+//#undef LEMMAUSEVTK
+
+#include "timer.h"
+
+using namespace Lemma;
+
+double randDouble(double low, double high) {
+	//srand(time(0));
+	double temp;
+	/* swap low & high around if the user makes no sense */
+	if (low > high)	{
+	temp = low;
+	low = high;
+	high = temp;
+	}
+	/* calculate the random number & return it */
+	temp = (rand() / (static_cast<double>(RAND_MAX) + 1.0))
+	* (high - low) + low;
+	return temp;
+}
+
+
+int main() {
+
+	// Keep track of time
+	jsw_timer timer;
+
+
+	srand(time(0));
+
+	PolygonalWireAntenna *pa = PolygonalWireAntenna::New();
+	pa->SetNumberOfFrequencies(1);
+	pa->SetFrequency(0, 1000);
+	pa->SetNumberOfPoints(5);
+	pa->SetPoint(0, Vector3r(   0,   0, -.001));
+	pa->SetPoint(1, Vector3r( 100,   0, -.001));
+	pa->SetPoint(2, Vector3r( 100, 100, -.001));
+	pa->SetPoint(3, Vector3r(   0, 100, -.001));
+	pa->SetPoint(4, Vector3r(   0,   0, -.001));
+	pa->SetCurrent(1.);
+    pa->SetNumberOfTurns(1);
+	//Vector3r rp = Vector3r::Random(3);
+	//rp << 150., 10, 0.;
+	//rp << -27.1456, 15.2350, -1e-3;
+	//rp << randDouble(-35,35), randDouble(-35,35), randDouble(-35,35);
+	//rp << 3.22806, -13.1548, 14.9695;
+	//rp.setRandom(3);
+	//std::cout <<  "rp " << rp.transpose() << std::endl;
+	//pa->ApproximateWithElectricDipoles(rp);
+
+	WireAntenna *wire = WireAntenna::New();
+		wire->SetNumberOfPoints(5);
+		wire->SetPoint(0, Vector3r(   0,   0, -1e-3));
+		wire->SetPoint(1, Vector3r( 10,   0, -1e-3));
+		wire->SetPoint(2, Vector3r( 10, 10, -1e-3));
+		wire->SetPoint(3, Vector3r(   0, 10, -1e-3));
+		wire->SetPoint(4, Vector3r(   0,   0, -1e-3));
+		// TODO change wire antennae to use my class
+		//wire->SetNumberOfFrequencies(1);
+		wire->SetCurrent(1.);
+		wire->SetNumberOfFrequencies(1);
+		wire->SetFrequency(0, 1000);
+        wire->SetNumberOfTurns(1);
+		//wire->ApproximateWithElectricDipoles(5);
+
+	// Define model
+	VectorXcr sigma(2);
+		sigma << Complex(0.,0), Complex(.1,0);
+	VectorXr  thick(1);
+		thick << 10;
+	LayeredEarthEM *earth = LayeredEarthEM::New();
+		earth->SetNumberOfLayers(2);
+		earth->SetLayerConductivity(sigma);
+		//earth->SetLayerThickness(thick);
+
+	// Receivers
+	ReceiverPoints *receivers = ReceiverPoints::New();
+		Vector3r loc;
+		Real ox    =    50.561 ;
+		Real oy    =   105.235 ;
+		Real depth =    -3.75e1;
+		Real depth2 =   depth;
+		Real dx    =    1.;
+		int nz     =    1;
+		receivers->SetNumberOfReceivers(nz);
+		int ir = 0;
+		for (int iz=0; iz<nz; ++iz) {
+			loc << ox, oy, depth;
+			receivers->SetLocation(ir, loc);
+			depth += dx;
+			++ ir;
+		}
+
+	// EmEarth
+	EMEarth1D  *EmEarth = EMEarth1D::New();
+		//EmEarth->AttachWireAntenna(wire);
+		EmEarth->AttachWireAntenna(pa);
+		EmEarth->AttachLayeredEarthEM(earth);
+		EmEarth->AttachReceiverPoints(receivers);
+		EmEarth->SetFieldsToCalculate(H);
+        //EmEarth->SetHankelTransformMethod(GAUSSIANQUADRATURE);
+
+	// Do calculation
+	timer.begin();
+	EmEarth->CalculateWireAntennaFields();
+	Real paTime = timer.end();
+	std::cout << "Polygonal wire antennae time: " << paTime << "\n";
+
+	//EmEarth->AttachWireAntenna(wire);
+	//timer.begin();
+	//EmEarth->CalculateWireAntennaFields();
+	//Real waTime = timer.end();
+	//std::cout << "Fixed wire antennae time: " << waTime << "\n";
+
+	depth = depth2;
+	std::fstream real("reale_lay.dat", std::ios::out);
+	std::fstream imag("image_lay.dat", std::ios::out);
+	std::fstream hreal("real_lay.dat", std::ios::out);
+	std::fstream himag("imag_lay.dat", std::ios::out);
+	for (int iz=0; iz<nz; ++iz) {
+		Vector3cr temp = receivers->GetEfield(0,iz);
+		real << ox << "\t" << oy << "\t" << depth << "\t"
+				<< temp(0).real() << "\t" << temp(1).real()
+				<< "\t" << temp(2).real() << std::endl;
+		imag << ox << "\t" << oy << "\t" << depth << "\t"
+				<< std::imag(temp(0)) << "\t" << std::imag(temp(1))
+				<< "\t" << std::imag(temp(2)) << std::endl;
+		temp = receivers->GetHfield(0, iz);
+		hreal << ox << "\t" << oy << "\t" << depth << "\t"
+				<< std::real(temp(0)) << "\t" << std::real(temp(1))
+				<< "\t" << std::real(temp(2)) << std::endl;
+		himag << ox << "\t" << oy << "\t" << depth << "\t"
+				<< std::imag(temp(0)) << "\t" << std::imag(temp(1))
+				<< "\t" << std::imag(temp(2)) << std::endl;
+		depth += dx;
+	}
+	real.close();
+	imag.close();
+	hreal.close();
+	himag.close();
+	EmEarth->Delete();
+	receivers->Delete();
+	earth->Delete();
+	//wire->Delete();
+
+#if  LEMMAUSEVTK
+	// Create the usual rendering stuff.
+	vtkRenderer           *renderer = vtkRenderer::New();
+	vtkRenderWindow         *renWin = vtkRenderWindow::New();
+	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+
+	// Line of tx
+	vtkLineSource       *vline = vtkLineSource::New();
+	vtkTubeFilter       *vTube = vtkTubeFilter::New();
+ 	vtkPolyDataMapper   *vMapper = vtkPolyDataMapper::New();
+	vtkActor             *vActor = vtkActor::New();
+	vline->SetPoint1(0,0,0);
+	vline->SetPoint2(10,0,0);
+
+	vTube->SetInputConnection(vline->GetOutputPort());
+	vTube->SetRadius(.2);
+	vTube->SetNumberOfSides(6);
+	vMapper->SetInputConnection(vTube->GetOutputPort());
+	vActor->SetMapper(vMapper);
+	vActor->GetProperty()->SetColor(0.0, .0, 1.0);
+	vActor->GetProperty()->SetOpacity(.15);
+	renderer->AddActor(vActor);
+
+	vtkLineSource       *vline2 = vtkLineSource::New();
+	vtkTubeFilter       *vTube2 = vtkTubeFilter::New();
+ 	vtkPolyDataMapper   *vMapper2 = vtkPolyDataMapper::New();
+	vtkActor             *vActor2 = vtkActor::New();
+	vline2->SetPoint1(10,0,0);
+	vline2->SetPoint2(10,10,0);
+
+	vTube2->SetInputConnection(vline2->GetOutputPort());
+	vTube2->SetRadius(.2);
+	vTube2->SetNumberOfSides(6);
+	vMapper2->SetInputConnection(vTube2->GetOutputPort());
+	vActor2->SetMapper(vMapper2);
+	vActor2->GetProperty()->SetColor(0.0, .0, 1.0);
+	vActor2->GetProperty()->SetOpacity(.15);
+	renderer->AddActor(vActor2);
+
+	vtkLineSource       *vline3 = vtkLineSource::New();
+	vtkTubeFilter       *vTube3 = vtkTubeFilter::New();
+ 	vtkPolyDataMapper   *vMapper3 = vtkPolyDataMapper::New();
+	vtkActor             *vActor3 = vtkActor::New();
+	vline3->SetPoint1(10,10,0);
+	vline3->SetPoint2(0,10,0);
+
+	vTube3->SetInputConnection(vline3->GetOutputPort());
+	vTube3->SetRadius(.2);
+	vTube3->SetNumberOfSides(6);
+	vMapper3->SetInputConnection(vTube3->GetOutputPort());
+	vActor3->SetMapper(vMapper3);
+	vActor3->GetProperty()->SetColor(0.0, .0, 1.0);
+	vActor3->GetProperty()->SetOpacity(.15);
+	renderer->AddActor(vActor3);
+
+	vtkLineSource       *vline4 = vtkLineSource::New();
+	vtkTubeFilter       *vTube4 = vtkTubeFilter::New();
+ 	vtkPolyDataMapper   *vMapper4 = vtkPolyDataMapper::New();
+	vtkActor             *vActor4 = vtkActor::New();
+	vline4->SetPoint1(0,10,0);
+	vline4->SetPoint2(0,0,0);
+
+	vTube4->SetInputConnection(vline4->GetOutputPort());
+	vTube4->SetRadius(.2);
+	vTube4->SetNumberOfSides(6);
+	vMapper4->SetInputConnection(vTube4->GetOutputPort());
+	vActor4->SetMapper(vMapper4);
+	vActor4->GetProperty()->SetColor(0.0, .0, 1.0);
+	vActor4->GetProperty()->SetOpacity(.15);
+	renderer->AddActor(vActor4);
+
+    loc << 50, 50, -1e-3;
+	pa->ApproximateWithElectricDipoles(loc);
+
+    vtkActor **pdipActors = new vtkActor*[pa->GetNumberOfDipoles()];
+    std::cout <<  "Wire approximated with " << pa->GetNumberOfDipoles() << std::endl;
+	for (int id=0; id<pa->GetNumberOfDipoles(); ++id) {
+		pdipActors[id] = pa->GetVtkActor(id);
+		renderer->AddActor(pdipActors[id]);
+	}
+
+    /*
+    vtkActor **dipActors = new vtkActor*[wire->GetNumberOfDipoles()];
+	for (int id=0; id<wire->GetNumberOfDipoles(); ++id) {
+		dipActors[id] = wire->GetVtkActor(id);
+		renderer->AddActor(dipActors[id]);
+	}
+    */
+
+	renderer->SetBackground(1,1,1);
+
+	// Render the window
+	renWin->AddRenderer(renderer);
+	renWin->SetWindowName("Wire antennae");
+
+	iren->SetRenderWindow(renWin);
+	iren->Initialize();
+	iren->Start();
+	iren->Render();
+
+    #if 0
+	cout << "Enter File name?: ";
+	std::string pngName;
+	std::cin >> pngName;
+	vtkPNGWriter *pngwrite = vtkPNGWriter::New();
+	vtkRenderLargeImage *renlarge = vtkRenderLargeImage::New();
+	renlarge->SetInput(renderer);
+	renlarge->SetMagnification(2);
+	pngwrite->SetInputConnection(renlarge->GetOutputPort());
+	pngName.append(".png");
+	pngwrite->SetFileName(pngName.c_str());
+	pngwrite->Write();
+    #endif
+
+#endif     // -----  not LEMMA_USE_VTK  -----
+
+
+	//std::cout << *pa << std::endl;
+	//pa->Delete();
+
+	return 0;
+}

Modules/LemmaCore/include/ASCIIParser.h

@@ -0,0 +1,161 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/23/2013 02:31:24 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2013, Trevor Irons
+ */
+
+#ifndef  ASCIIPARSER_INC
+#define  ASCIIPARSER_INC
+
+#include "LemmaObject.h"
+#include <fstream>
+
+namespace Lemma {
+
+/**
+ * \ingroup LemmaCore
+ * \brief   Provides ASCII input file parsing
+ * \details Able to accomodate various inputs and comment styles
+ */
+class ASCIIParser : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const ASCIIParser &ob);
+
+    public:
+
+    // ====================  LIFECYCLE     =======================
+
+    /**
+     *  Factory method for generating concrete class.
+     *  @return a std::shared_ptr of type ASCIIParser
+     */
+    static std::shared_ptr< ASCIIParser >  NewSP();
+
+    /**
+     *  Constructs an object from a YAML serialization
+     *  @return a std::shared_ptr of type ASCIIParser
+     */
+    static std::shared_ptr< ASCIIParser >  DeSerialize( const YAML::Node& node );
+
+    /**
+     *  Uses YAML to serialize this object.
+     *  @return a YAML::Node
+     */
+    YAML::Node Serialize() const;
+
+    // ====================  OPERATORS     =======================
+
+    // ====================  OPERATIONS    =======================
+
+    /** Opens the file up for reading. Throws an exception if something bad happens.
+        @param[in] fname is the filename to be parsed.
+     */
+    void Open (const std::string& fname);
+
+
+    /** Closes the file. Throws an exception if something bad happens.
+     */
+    void Close ( );
+
+    /** Reads a series of Reals.
+     *  @param[in] nr is the number of reals to read. Enter -1 for the entire file
+     *  @todo throw exception if no file is open
+     */
+    std::vector<Real> ReadReals( const int& nr);
+
+    /** Reads a series of Reals.
+     *  @param[in] nr is the number of ints to read. Enter -1 for the entire file
+     *  @todo throw exception if no file is open
+     */
+    std::vector<int> ReadInts( const int& nr);
+
+    /** Reads a series of space delimited strings
+     *  @param[in] nr is the number of strings to read. Enter -1 for the entire file
+     *  @todo throw exception if no file is open
+     */
+    std::vector< std::string > ReadStrings( const int& nr);
+
+    /**
+     * @param[in] loc is the point in the file to jump to. Uses seekg
+     */
+    void JumpToLocation(const int& loc);
+
+    // ====================  ACCESS        =======================
+
+    /** Sets the comment identifier key.
+     * @param[in] key is a string identifying comments. All text after the key will be
+     *  ignored by the parser. Default is //
+     */
+    void SetCommentString( const std::string& key );
+
+    /** Sets the buffer size. This affects the maximum number of column in a line. Defaults
+     *  is 255.
+     *  @param[in] BufferSize is the size of the buffer to use
+     */
+    void SetBufferSize( const int& BufferSize);
+
+    /**
+     * @return the current position in the file, as reported by istream::tellg
+     */
+    int GetFileLocation();
+
+    /** Returns the name of the underlying class, similiar to Python's type */
+    virtual inline std::string GetName() const {
+        return CName;
+    }
+
+    // ====================  INQUIRY       =======================
+
+    protected:
+
+    // ====================  LIFECYCLE     =======================
+
+    /** Default protected constructor, use New */
+    ASCIIParser ( );
+
+    /** Constructor using YAML::Node */
+    ASCIIParser ( const YAML::Node& node );
+
+    /** Default protected destructor, use Delete */
+    ~ASCIIParser ();
+
+    /**
+     *  @copybrief   LemmaObject::Release()
+     *  @copydetails LemmaObject::Release()
+     */
+    void Release();
+
+    private:
+
+    // ====================  DATA MEMBERS  =========================
+
+    /** ASCII string representation of the class name */
+    static constexpr auto CName = "ASCIIParser";
+
+    /** c++ style file IO */
+    std::fstream input;
+
+    /** comment string, defaults to c++ style // */
+    std::string  CommentString;
+
+    /** Buffer size, max line width supported, defaults to 255 */
+    int BufferSize;
+
+}; // -----  end of class  ASCIIParser  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef ASCIIPARSER_INC  -----

Modules/LemmaCore/include/CubicSplineInterpolator.h

@@ -0,0 +1,195 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/25/2013 08:20:14 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2013, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2013, Trevor Irons
+ */
+
+#ifndef  CUBICSPLINEINTERPOLATOR_INC
+#define  CUBICSPLINEINTERPOLATOR_INC
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+// Simple struct to hold spline terms
+struct SplineSet{
+    VectorXr a;
+    VectorXr b;
+    VectorXr c;
+    VectorXr d;
+    VectorXr x;
+
+    SplineSet( ) {
+    }
+
+    SplineSet(const int&n) {
+        a = VectorXr::Zero(n+1);
+        b = VectorXr::Zero(n);
+        c = VectorXr::Zero(n+1);
+        d = VectorXr::Zero(n);
+        x = VectorXr::Zero(n+1);
+    }
+};
+
+/**
+ * \ingroup LemmaCore
+ * \brief   Real 1D Natural cubic spline interpolator.
+ * \details Splines are fit between knots \f$j\f$ according to the forulae
+ *          \f[ S_j(x) =  a_j + b_j(x - x_j) + c_j(x-x_j)^2 + d_j(x-x_y)^3  \f]
+ *          The spline must satisfy the following conditions
+ *           \f{eqnarray} {
+ *               S_i(x_i) & = & y_i = S_{i-1}(x_i), i = 1,..., n-1 \\
+ *               S'_i(x_i) & = &  S'_{i-1}(x_i), i = 1,..., n-1 \\
+ *               S''_i(x_i) & = & S''_{i-1}(x_i), i = 1,..., n-1 \\
+ *               S''_0(x_0) & = & S''_{n-1}(x_n) = 0
+ *           \f}
+ */
+class CubicSplineInterpolator : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const CubicSplineInterpolator &ob);
+
+    public:
+
+    // ====================  LIFECYCLE     =======================
+
+    /**
+     *  Factory method for generating concrete class.
+     *  @return a std::shared_ptr of type CubicSplineInterpolator
+     */
+    static std::shared_ptr<CubicSplineInterpolator> NewSP();
+
+    /**
+     *  Uses YAML to serialize this object.
+     *  @return a YAML::Node
+     */
+    YAML::Node Serialize() const;
+
+    /**
+     *   Constructs an object from a YAML::Node.
+     */
+    static std::shared_ptr< CubicSplineInterpolator > DeSerialize(const YAML::Node& node);
+
+    // ====================  OPERATORS     =======================
+
+    // ====================  OPERATIONS    =======================
+
+    /** Sets the knots to use for interpolation.
+        @param[in] x are the absissa values
+        @param[in] y are the ordinate values
+     */
+    void SetKnots(const VectorXr& x, const VectorXr& y);
+
+    /** Resets the knots to use for interpolation, when abscissa values haven't changed.
+        @param[in] y are the ordinate values
+     */
+    void ResetKnotOrdinate( const VectorXr& y );
+
+    /** Interpolate a monotonically increasing ordered set.
+        @param[in] x are the interpolation abscissa points
+        @return the ordinate values at x
+     */
+    VectorXr InterpolateOrderedSet(const VectorXr& x);
+
+    /** integrates the spline from x0 to x1. Uses composite Simpson's rule and n is the number of segments
+     *  @param[in] x0 is left argument
+     *  @param[in] x1 is right argument
+     *  @param[in] n is the number of points, must be even
+     */
+    Real Integrate(const Real& x0, const Real& x1, const int& n);
+
+    /** integrates using cubic spline values. Taken from AMRIRA P223F project code Leroi, which in turn was based on
+        This is a modification of the FUNCTION PPVALU in the book
+        "A PRACTICAL GUIDE TO SPLINES"  by C. DE BOOR
+     */
+    Real Integrate(const Real& x0, const Real& x1);
+
+    /** @returns the know abscissa values
+     */
+    VectorXr GetKnotAbscissa();
+
+    /** @returns the know abscissa values
+     */
+    VectorXr GetKnotOrdinate();
+
+    /** Interpolation at a single point.
+        @param[in] x is the interpolation abscissa point
+        @param[in] i is an optional index to start searching at. Defaults to zero
+        @return the ordinate value at x
+     */
+    Real Interpolate(const Real& x, int& i);
+
+    /** Interpolation at a single point.
+        @param[in] x is the interpolation abscissa point
+        @return the ordinate value at x
+     */
+    Real Interpolate(const Real& x);
+
+    // ====================  ACCESS        =======================
+
+    // ====================  INQUIRY       =======================
+
+    /** Returns the name of the underlying class, similiar to Python's type */
+    virtual inline std::string GetName() const {
+        return CName;
+    }
+
+    protected:
+
+    // ====================  LIFECYCLE     =======================
+
+    /** Default protected constructor, use New */
+    CubicSplineInterpolator ( );
+
+    /** Protected DeDerializing constructor, use factory DeSerialize  method*/
+    CubicSplineInterpolator (const YAML::Node& node);
+
+    /** Default protected destructor, smart pointers auto delete */
+    ~CubicSplineInterpolator ();
+
+    /**
+     *  @copybrief   LemmaObject::Release()
+     *  @copydetails LemmaObject::Release()
+     */
+    void Release();
+
+    // ====================  OPERATIONS    =======================
+
+    /** Finds the interval of knots in spline to use for integration.
+     */
+    int Interval(const Real& x);
+
+    private:
+
+    /** ASCII string representation of the class name */
+    static constexpr auto CName = "CubicSplineInterpolator";
+
+    SplineSet Spline;
+
+    int ilo;
+
+    int mflag;
+
+    // ====================  DATA MEMBERS  =========================
+
+}; // -----  end of class  CubicSplineInterpolator  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef CUBICSPLINEINTERPOLATOR_INC  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/Data.h

@@ -0,0 +1,94 @@
+/* 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     03/01/2010
+  @version  $Id: data.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef  _DATA_H_INC
+#define  _DATA_H_INC
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+    // ===================================================================
+    //        Class:  Data
+    /**
+     *  \ingroup LemmaCore
+     *  \brief Abstract class for data from an instrument
+     *  \details Currently provides no functionality, but this will likely
+     *          change in the near future.
+     */
+    // ===================================================================
+    class Data : public LemmaObject {
+
+        friend std::ostream &operator<<(std::ostream &stream,
+                const Data &ob);
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+            /** Copies the *structure* of this class to a new object and
+             * returns a pointer to this object. Memory management is
+             * the responsibility of the receiver!. The data values are
+             * not copied.
+             */
+            virtual Data* Clone()=0;
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            /** Clears all data arrays and containers
+             */
+            virtual void Zero()=0;
+
+            /** Data2 needs to be compatible with the calling class. Data2 will
+             * be cast into the calling class type. Most of the time this
+             * cast will not perform anything, as Data2 *is* the same type
+             * as the calling class.
+             */
+            virtual Real Norm(Data* Data2)=0;
+
+            // ====================  ACCESS        =======================
+
+            // ====================  INQUIRY       =======================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /** Default protected constructor. */
+            Data (  );
+
+            /** Default protected constructor. */
+            ~Data ();
+
+            // ====================  DATA MEMBERS  =========================
+
+        private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "Data";
+
+    }; // -----  end of class  Data  -----
+
+}       // -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef _DATA_H_INC  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/DataReader.h

@@ -0,0 +1,75 @@
+/* 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     03/01/2010
+  @version  $Id: datareader.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef  DATAREADER_H_INC
+#define  DATAREADER_H_INC
+
+#include "LemmaObject.h"
+#include "Data.h"
+
+namespace Lemma {
+
+// ===================================================================
+//        Class:  DataReader
+/// \ingroup LemmaCore
+/// \brief  Abstract class whose derived members read data files and
+/// converts them to Data Classes
+/// \details
+// ===================================================================
+class DataReader : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const DataReader &ob);
+
+    public:
+
+        // ====================  LIFECYCLE     =======================
+
+        // ====================  OPERATORS     =======================
+
+        // ====================  OPERATIONS    =======================
+
+        // ====================  ACCESS        =======================
+        virtual Data* GetData()=0;
+
+        // ====================  INQUIRY       =======================
+        /** Returns the name of the underlying class, similiar to Python's type */
+        virtual inline std::string GetName() const {
+            return CName;
+        }
+
+    protected:
+
+        // ====================  LIFECYCLE     =======================
+
+        /// Default protected constructor.
+        DataReader ( );
+
+        /// Default protected constructor.
+        ~DataReader ( );
+
+        // ====================  DATA MEMBERS  =========================
+
+    private:
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "DataReader";
+
+}; // -----  end of class  DataReader  -----
+
+}       // -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef DATAREADER_H_INC  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/EarthModel.h

@@ -0,0 +1,149 @@
+/* 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     03/23/2010
+  @version  $Id: earthmodel.h 198 2014-12-27 06:14:05Z tirons $
+ **/
+
+#ifndef  EARTHMODEL_INC
+#define  EARTHMODEL_INC
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+    // ===================================================================
+    //        Class:  EarthModel
+    /// \ingroup LemmaGroup
+    /// \brief abstract class for Earth models
+    /// \details
+    // ===================================================================
+    class EarthModel : public LemmaObject {
+
+        // ====================    FRIENDS     ===========================
+        friend std::ostream &operator<<(std::ostream &stream,
+                        const EarthModel &ob);
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+            /** YAML Serializing method
+             */
+            YAML::Node Serialize() const;
+            //static LayeredEarth* DeSerialize(const YAML::Node& node);
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            // ====================  ACCESS        =======================
+
+            /** Sets the earth's magnetic field, in Tesla
+                @param[in] bfield contains the x,y, and z components of
+                the field. EarthModel::Inc, EarthModel::Dec, and
+                EarthModel::BMag are computed:
+                \f{eqnarray*}{ \mathrm{Inc} = &
+                \cos^{-1}( \mathbf{B} \cdot \hat{\mathbf{z}} / B )
+                \\ \mathrm{Dec} =& \cos^{-1}
+                         (\mathbf{B} \cdot \hat{\mathbf{x}} / B )
+                   \\           \mathrm{Mag} =& \Vert \mathbf{B} \Vert  \f}
+               @param[in] unit specifies the units of bfield.
+              */
+            void SetMagneticFieldComponents(const Vector3r &bfield,
+                            const MAGUNITS &unit);
+
+            /// Sets the magnetic field using inclination, declination,
+            /// and magnitude, units of Tesla
+            /// @param[in] inc is the field's inclination
+            /// @param[in] dec is the field's declination
+            /// @param[in] Mag is the field's magnitude
+            /// @param[in] unit specifies the the units used.
+            /// The components are set according to
+            /// \f{eqnarray*} { B_x =& B \cos(inc (\pi/180)) \cos(dec (\pi/180))
+            ///            \\   B_y =& B \cos(inc (\pi/180)) \sin(dec (\pi/180))
+            ///            \\   B_z =& B \sin(inc (\pi/180))
+            /// \f}
+            void SetMagneticFieldIncDecMag(const Real& inc, const Real&dec,
+                            const Real& Mag, const MAGUNITS &unit);
+
+            // ====================  INQUIRY       =======================
+
+            /// Gets the magnetic field of the earth, in T
+            Vector3r GetMagneticField( );
+
+            /// Gets the magnetic field of the earth, in T
+            Vector3r GetMagneticFieldInGauss(  );
+
+            /// Gets the magnetic field unit vector of the earth
+            Vector3r GetMagneticFieldUnitVector( );
+
+            /// Returns the magnitude of the magnetic field \f$ \Vert
+            /// \mathbf{B}_0 \Vert \f$
+            Real GetMagneticFieldMagnitude();
+
+            /// Returns the magnitude of the magnetic field \f$ \Vert
+            /// \mathbf{B}_0 \Vert \f$ in Gauss units
+            Real GetMagneticFieldMagnitudeInGauss( );
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /// Default protected constructor.
+            EarthModel ( );
+
+            /** Deserialize constructor */
+            EarthModel (const YAML::Node& node);
+
+            /// Default protected constructor.
+            ~EarthModel ();
+
+            // ====================  DATA MEMBERS  =========================
+
+        private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "EarthModel";
+
+            /// Magnetic field, in units of nT
+            Vector3r          BField;
+
+            /// Magnetic field unit vector \f$  \hat{\mathbf{B}_0} \f$
+            Vector3r          BFieldUnit;
+
+            /// Inclination of the magnetic field
+            Real              BInc;
+
+            /// Declination of the magnetic field
+            Real              BDec;
+
+            /// Declination of the magnetic field
+            Real              BMag;
+
+    }; // -----  end of class  EarthModel  -----
+
+    //////////////////////////////////////////////////////////////
+    // Error Classes
+
+    /** If the Earth class is NULL valued, throw this
+     */
+    class NullEarth : public std::runtime_error {
+            /** Thrown when an Earth pointer is NULL
+             */
+            public: NullEarth ();
+    };
+
+}       // -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef EARTHMODEL_INC  -----

Modules/LemmaCore/include/Filter.h

@@ -0,0 +1,71 @@
+/* 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     07/20/2010
+  @version  $Id: filter.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+
+    // ===================================================================
+    //  Class:  Filter
+    /**
+      @class  Filter
+      \ingroup LemmaCore
+      \brief  Abstract class for filtering.
+      \details Derived classes include fourier domain and time domain filters.
+     */
+    // ===================================================================
+
+    class Filter : public LemmaObject {
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            // ====================  ACCESS        =======================
+
+            // ====================  INQUIRY       =======================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /// Default protected constructor.
+            Filter ( );
+
+            /// Default protected constructor.
+            ~Filter ();
+
+            // ====================  DATA MEMBERS  =========================
+
+        private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "Filter";
+
+    }; // -----  end of class  Filter  -----
+
+
+}		// -----  end of Lemma  name  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/Grid.h

@@ -0,0 +1,78 @@
+/* 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     10/25/2009
+  @version  $Id: grid.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef __GRID_H
+#define __GRID_H
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+// ===============================================================================
+//        Class:  Grid
+/// \ingroup LemmaCore
+/// \brief Abstract base class for grids.
+/// \details Grids define values in space, they may either be an integral part
+/// of models, or simply the output of a forward modelling procedure.
+// ===============================================================================
+
+class Grid : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const Grid &ob);
+
+    public:
+
+        // ====================  LIFECYCLE     ===================================
+
+        YAML::Node Serialize() const;
+
+        // ====================  OPERATORS     ===================================
+
+        // ====================  OPERATIONS    ===================================
+
+        // ====================  ACCESS        ===================================
+
+        // ====================  INQUIRY       ===================================
+
+        /** Returns the name of the underlying class, similiar to Python's type */
+        virtual inline std::string GetName() const {
+            return this->CName;
+        }
+
+    protected:
+
+        // ====================  LIFECYCLE     ===================================
+
+        /** Protected DeDerializing constructor, use factory DeSerialize  method on non abstract classes*/
+        Grid (const YAML::Node& node);
+
+        /// Default protected constructor.
+        Grid ( );
+
+        /// Default protected constructor.
+        ~Grid ();
+
+        // ====================  DATA MEMBERS  ===================================
+
+    private:
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "Grid";
+
+}; // -----  end of class  Grid  -----
+
+} // namespace Lemma
+
+#endif // __GRID_H
+

Modules/LemmaCore/include/GridReader.h

@@ -0,0 +1,84 @@
+/* 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
+  @date      09/23/2013 11:25:11 AM
+  @version   $Id$
+  @author    Trevor Irons (ti)
+  @email     Trevor.Irons@xri-geo.com
+  @copyright Copyright (c) 2013, XRI Geophysics, LLC
+  @copyright Copyright (c) 2013, Trevor Irons
+ **/
+
+#ifndef  GRIDREADER_INC
+#define  GRIDREADER_INC
+
+#include "LemmaObject.h"
+#include "Grid.h"
+
+namespace Lemma {
+
+    /**
+      @class
+      \ingroup LemmaCore
+      \brief   Abstract class for grid readers.
+      \details Provides a consistent interface for reading grid files in Lemma.
+     */
+    class GridReader : public LemmaObject {
+
+        friend std::ostream &operator << (std::ostream &stream, const GridReader &ob) {
+            stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+            return stream;
+        }
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            virtual void ReadASCIIGridFile( const std::string& name )=0;
+
+            virtual std::shared_ptr< Grid > GetGrid( )=0;
+
+            // ====================  ACCESS        =======================
+
+            // ====================  INQUIRY       =======================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /** Default protected constructor, use New */
+            GridReader ( ) : LemmaObject( ) {
+            }
+
+            /** Default protected constructor, use Delete */
+            ~GridReader () {
+            }
+
+        private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "GridReader";
+
+            // ====================  DATA MEMBERS  =========================
+
+    }; // -----  end of class  GridReader  -----
+
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef GRIDREADER_INC  -----

Modules/LemmaCore/include/Instrument.h

@@ -0,0 +1,71 @@
+/* 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     02/05/2010
+  @version  $Id: instrument.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef __INSTRUMENT_H
+#define __INSTRUMENT_H
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+	// =======================================================================
+	//        Class:  Instrument
+    /// \ingroup LemmaCore
+	/// \brief    Abstract class for instruments
+	/// \details
+	// =======================================================================
+	class Instrument : public LemmaObject {
+
+		friend std::ostream &operator<<(std::ostream &stream,
+			const Instrument &ob);
+
+		public:
+
+			// ====================  LIFECYCLE     ===========================
+
+			// ====================  OPERATORS     ===========================
+
+			// ====================  OPERATIONS    ===========================
+            //virtual void ForwardModelOnOctreeMesh()=0;
+
+			// ====================  ACCESS        ===========================
+
+			// ====================  INQUIRY       ===========================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+			// ====================  DATA MEMBERS  ===========================
+
+		protected:
+
+			// ====================  LIFECYCLE     ===========================
+
+			/// Default protected constructor.
+			Instrument ( );
+
+			/// Default protected constructor.
+			~Instrument ();
+
+		private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "Instrument";
+
+	}; // -----  end of class  Instrument  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif // __INSTRUMENT_H

Modules/LemmaCore/include/LayeredEarth.h

@@ -0,0 +1,174 @@
+/* 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     06/24/2009
+  @version   0.0
+ **/
+
+#ifndef __LAYEREDEARTH_H
+#define __LAYEREDEARTH_H
+
+#include "EarthModel.h"
+
+namespace Lemma {
+
+	const int MAXLAYERS	= 10;
+
+	// =======================================================================
+	//        Class:  LayeredEarth
+    /// \ingroup LemmaCore
+	/// \brief Abstract 1D layered earth.
+	// =======================================================================
+	class LayeredEarth : public EarthModel {
+
+        friend class LayeredEarthEM;
+
+		public:
+
+			// ====================    FRIENDS     ===========================
+			/** Recursively prints out information about this class
+             */
+            friend std::ostream &operator<<(std::ostream &stream,
+						const LayeredEarth &ob);
+
+			// ====================  LIFECYCLE     ===========================
+
+			// ====================  OPERATORS     ===========================
+
+			// ====================  OPERATIONS    ===========================
+
+			// ====================  ACCESS        ===========================
+
+			/** Sets the number of layers and resizes all model parameters to
+			 *  be the correct size, finally initialises all values to free
+			 *  space.
+             */
+			virtual void SetNumberOfLayers(const int& nlay)=0;
+
+			/** Sets the thickness of the layers
+             */
+			void SetLayerThickness(const VectorXr &thick);
+
+			// ====================  INQUIRY       ===========================
+
+			/** Returns the number of layers in the model */
+			int GetNumberOfLayers( );
+
+			/** Returns the number of layers in the model, excluding the air
+			 * layer
+             */
+			int GetNumberOfNonAirLayers( );
+
+			/** Returns the layer that a depth is in.
+			 * @param[in] depth is the depth you desire to know the layer of
+			 * @return the indice of the layer that depth is in.
+             */
+			int GetLayerAtThisDepth(const Real& depth);
+
+			/** Returns the thickness of a layer
+             *  @param[in] lay is the layer indice
+             *  @return the thickness of a layer
+             */
+			Real GetLayerThickness(const int & lay);
+
+			/** Returns the depth of the bottom interface of a layer
+             *  @param[in] lay is the layer indice
+             *  @return the depth of the bottom interface of a layer
+             */
+			Real GetLayerDepth(const int & lay);
+
+            /** YAML Serializing method
+             */
+            YAML::Node Serialize() const;
+            //static LayeredEarth* DeSerialize(const YAML::Node& node);
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+		protected:
+
+			// ====================  LIFECYCLE     ===========================
+
+			/** Default protected constructor. */
+			LayeredEarth ( );
+
+            /** Default protected constructor. */
+			LayeredEarth (const YAML::Node& node);
+
+			/** Default protected constructor. */
+			~LayeredEarth ();
+
+			// ====================  DATA MEMBERS  ===========================
+
+        private:
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "LayeredEarth";
+
+			/** Number of layers in the model, including the air layer,
+			 * and counting from 0
+             */
+			int               NumberOfLayers;
+
+			/** Number of interfaces in the model. Equal to NumberOfLayers - 1
+			 */
+            int               NumberOfInterfaces;
+
+			/** Vector of layer thicknesses */
+			VectorXr          LayerThickness;
+
+	}; // -----  end of class  LayeredEarth  -----
+
+
+    ///////////////////////////////////////////////////////////////////
+    // Error classes
+
+    /** If a model with less than two layers is specified, throw error. Not set up
+     * for whole space solutions.
+     */
+    class EarthModelWithLessThanTwoLayers : public std::runtime_error {
+        public:
+            EarthModelWithLessThanTwoLayers( );
+    };
+
+    /** If a solver has a max number of layers, and this is exceeded, throw this
+     * error.
+     */
+    class EarthModelWithMoreThanMaxLayers : public std::runtime_error {
+        public:
+            /** Thrown when an Earth model with more than the maximum number of
+             *  Layers is given
+             */
+            EarthModelWithMoreThanMaxLayers( );
+    };
+
+    /** If the model parameters have different lengths than the earth model, throw
+     * this error.
+     */
+    class EarthModelParametersDoNotMatchNumberOfLayers : public std::runtime_error {
+        public:
+            /** Thrown when the parameters do not match the number of Layers
+             */
+            EarthModelParametersDoNotMatchNumberOfLayers ( );
+    };
+
+    /** If a request if made for a non-valid earth model parameter, throw this
+     * error.
+     */
+    class RequestForNonValidEarthModelParameter : public std::runtime_error {
+        public:
+            /** Thrown when a request is made for a non-valid earth parameter
+             */
+            RequestForNonValidEarthModelParameter ( );
+    };
+
+} // namespace Lemma
+
+#endif // __LAYEREDEARTH_H

Modules/LemmaCore/include/LemmaCore

@@ -0,0 +1,40 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API */
+
+/* Copyright (C) 2009-2013 Trevor Irons <trevorirons@gmail.com> */
+
+/* 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
+  @version  $Id$
+ **/
+
+// Convenience include all of LemmaCore
+#ifndef  LEMMAINCLUDE_INC
+#define  LEMMAINCLUDE_INC
+#include "ASCIIParser.h"
+#include "CubicSplineInterpolator.h"
+//#include "Data.h"
+//#include "DataReader.h"
+//#include "EarthModel.h"
+//#include "Filter.h"
+//#include "Grid.h"
+//#include "GridReader.h"
+//#include "Instrument.h"
+//#include "LayeredEarth.h"
+//#include "LemmaObject.h"
+#include "RectilinearGrid.h"
+#include "RectilinearGridReader.h"
+#include "RectilinearGridVTKExporter.h"
+#include "WindowFilter.h"
+// helpers and utilities
+#include "banner.h"
+#include "helper.h"
+#include "timer.h"
+#endif   // ----- #ifndef LEMMAINCLUDE_INC  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/LemmaObject.h

@@ -0,0 +1,168 @@
+/* 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     06-13-2016
+ **/
+
+#ifndef __LEMMAOBJECT_H
+#define __LEMMAOBJECT_H
+
+#include "helper.h"
+#include "lemma.h"
+#include "yaml-cpp/yaml.h"
+
+#include <chrono>
+#include <memory>
+
+namespace Lemma {
+
+/**
+  * \ingroup  LemmaCore
+  * \brief    Abstract class providing common interface for Lemma Objects.
+  * \details  Lemma objects can be members of other Lemma, and may be members
+  *           of multiple objects. Since updates should be atomic, and objects
+  *           can be large, it becomes useful to count the number of
+  *           Classes an object is a member of.
+  *           Before C++-11, this was done internally in Lemma, with the inclusion of
+  *           more sophisticated smart pointers, this logic has been offloaded to the
+  *           standard. All Lemma objects should be created as C++-11 Smart pointers, using
+  *           the supplied New method. Calls to Delete are no longer necessary or available.
+  */
+class LemmaObject {
+
+    /**
+     *  Streams class information as YAML::Node
+     */
+    friend YAML::Emitter& operator << (YAML::Emitter& out, const LemmaObject &ob) ;
+
+    friend class LemmaObjectDeleter;
+
+    public:
+
+        // ====================  LIFECYCLE     ==============================
+
+        // Needed because many derived classes have Eigen vectors as members,
+        // causing alignment issues when vectorisation is enabled.
+        EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+
+        // ====================  OPERATORS     ==============================
+
+        // ====================  OPERATIONS    ==============================
+
+        // ====================  ACCESS        ==============================
+
+        // ====================  INQUIRY       ==============================
+
+        /** Returns the name of the underlying class; Run-time type information (RTTI). This approach
+            Was chosen over typeid due to name mangling among various compilers, and the need for consistency
+            in Serialized objects.
+         */
+        virtual inline std::string GetName() const {
+            return this->CName;
+        }
+
+        /**
+         *  Uses YAML to serialize this object.
+         *  @return a YAML::Node
+         *  FOR NOT LemmaObject does not write out any YAML info,
+         *      in the future the reference count could be logged? But the utility
+         *      of that is minimal.
+         *  @note Not every Lemma class needs to be Serializable, for instance HankelTransform
+         *        classes will never need to be Serialized. There may be a need to differentiate these
+         *        two families in to a LemmaInternalClass without serializing and perhaps this class for
+         *        all external classes that might need to be serialized.
+         */
+        virtual YAML::Node Serialize() const {
+            YAML::Node node = YAML::Node();
+            node.SetTag( GetName() );
+            std::time_t now = std::chrono::system_clock::to_time_t( std::chrono::system_clock::now() );
+            node["Serialized"] = std::ctime(&now);
+            return node;
+        };
+
+    protected:
+
+        // ====================  LIFECYCLE     ==============================
+
+        /** Protected default constructor. This is an abstract class and
+         *  cannot be instantiated.
+         */
+        LemmaObject ( );
+
+        /** Protected DeDerializing constructor, use factory DeSerialize  method*/
+        LemmaObject (const YAML::Node& node);
+
+        /** Protected default destructor. This is an abstract class and
+         *  cannot be instantiated. Virtual is necessary so that if base class destructor is
+         *  called, we get the right behaviour.
+         */
+        virtual ~LemmaObject();
+
+    private:
+
+        // ====================  DATA MEMBERS  ==============================
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "LemmaObject";
+
+}; // -----  end of class  LemmaObject  -----
+
+class LemmaObjectDeleter
+{
+    public:
+        void operator()(LemmaObject* p) { delete p; }
+};
+
+    /////////////////////////////////////////////////////////////////
+    // Error Classes
+
+    /** Error called when DeSerializing breaks. If the node type is not the expected one
+     *  this error is thown.
+     */
+    class DeSerializeTypeMismatch : public std::runtime_error {
+        public:
+            DeSerializeTypeMismatch(const std::string& expected, const std::string& got);
+
+    };
+
+    /** If an assignment is made that is out of bounts, throw this.
+     */
+    class AssignmentOutOfBounds : public std::runtime_error {
+        public:
+            /** Throw when an assignment is out of bounds.
+             *  @param[in] ptr is a pointer to the class throwing the exception.
+             */
+            AssignmentOutOfBounds(LemmaObject *ptr);
+    };
+
+    /** If a pointer to a class is requested, but it is NULL valued, throw this
+     */
+    class RequestToReturnNullPointer : public std::runtime_error {
+        public:
+            /** Thrown when the pointer is NULL
+             *  @param[in] ptr is a pointer to the class throwing the exception.
+             */
+            RequestToReturnNullPointer(LemmaObject *ptr);
+    };
+
+    /** If an error in opening a .mat file is encountered, throw this.
+     */
+    class MatFileCannotBeOpened : public std::runtime_error {
+        /** thown when a mat file fails to be opened.
+         */
+            public: MatFileCannotBeOpened();
+    };
+
+	/** Generic file I/O error. */
+	class GenericFileIOError : public std::runtime_error {
+			public: GenericFileIOError(LemmaObject *ptr, const std::string &filename);
+	};
+
+}
+#endif // __LEMMAOBJECT_H

Modules/LemmaCore/include/RectilinearGrid.h

@@ -0,0 +1,216 @@
+/* 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     10/28/2010
+  @version  $Id: rectilineargrid.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef  RECTILINEARGRID_INC
+#define  RECTILINEARGRID_INC
+
+#include "Grid.h"
+
+namespace Lemma {
+
+    // ===================================================================
+    //  Class:  RectilinearGrid
+    /**
+      @class
+      \ingroup LemmaCore
+      \brief   Impliments a rectilinear grid.
+      \details A rectilinear grid can vary regularily in space but must be
+               constant variation across each dimension. In this way three
+               vectors can define the entire grid, and each cell is right
+               aligned with all its neighbours.
+     */
+    // ===================================================================
+    class RectilinearGrid : public Grid {
+
+        friend std::ostream &operator<<(std::ostream &stream,
+			const RectilinearGrid &ob);
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+            /**
+             *  Factory method for generating concrete class.
+             *  @return a std::shared_ptr of type RectilinearGrid
+             */
+            static std::shared_ptr<RectilinearGrid> NewSP();
+
+
+            /**
+             *  Uses YAML to serialize this object.
+             *  @return a YAML::Node
+             */
+            YAML::Node Serialize() const;
+
+            /**
+             *   Constructs an object from a YAML::Node.
+             */
+            static std::shared_ptr< RectilinearGrid > DeSerialize(const YAML::Node& node);
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            /**
+             *  Sets the dimensions in x, y and z
+             *  @param[in] nx is the number of cells in x
+             *  @param[in] ny is the number of cells in y
+             *  @param[in] nz is the number of cells in z
+             */
+            void SetDimensions (const int &nx, const int &ny, const int &nz);
+
+            /**
+             *  Sets the offset in x, y and z. This sets the corner position of the top SW corner.
+             *  @param[in] ox is the number of cells in x
+             *  @param[in] oy is the number of cells in y
+             *  @param[in] oz is the number of cells in z
+             */
+            void SetOffset (const Real &ox, const Real &oy, const Real &oz);
+
+            /**
+             *  Sets the spacing in x, y and z
+             *  @param[in] dx is a Vector of cell spacing in x, must be same
+             *  size ad nx
+             *  @param[in] dy is a Vector of cell spacing in x, must be same
+             *  size ad ny
+             *  @param[in] dz is a Vector of cell spacing in x, must be same
+             *  size ad nz
+             */
+            void SetSpacing (const VectorXr &dx, const VectorXr &dy,
+                const VectorXr &dz);
+
+
+            // ====================  ACCESS        =======================
+
+            /** Returns the number of cells in x
+             *  @return number of cells in x
+             */
+            int GetNx();
+
+            /** Returns the number of cells in y
+             *  @return number of cells in y
+             */
+            int GetNy();
+
+            /** Returns the number of cells in z
+             *  @return number of cells in z
+             */
+            int GetNz();
+
+            /** Returns the offset of cells in x
+             *  @return offset of cells in x
+             */
+            Real GetOx();
+
+            /** Returns the number of cells in y
+             *  @return number of cells in y
+             */
+            Real GetOy();
+
+            /** Returns the offset of cells in z
+             *  @return offset of cells in z
+             */
+            Real GetOz();
+
+            /** Returns Spacing Vector
+             *  @return dx
+             */
+            VectorXr GetDx();
+
+            /** Returns a spacing
+             *  @return dx[ix]
+             */
+            Real GetDx(const int& ix);
+
+            /** Returns Spacing Vector
+             *  @return dy
+             */
+            VectorXr GetDy();
+
+            /** Returns a spacing
+             *  @return dy[iy]
+             */
+            Real GetDy(const int& iy);
+
+            /** Returns Spacing Vector
+             *  @return dz
+             */
+            VectorXr GetDz();
+
+            /** Returns Spacing Vector
+             *  @return dz[iz]
+             */
+            Real GetDz(const int& iz);
+
+            // ====================  INQUIRY       =======================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /** Default protected constructor. */
+            RectilinearGrid ( );
+
+            /** Protected DeDerializing constructor, use factory DeSerialize  method*/
+            RectilinearGrid (const YAML::Node& node);
+
+            /** Default protected constructor. */
+            ~RectilinearGrid ();
+
+            // ====================  DATA MEMBERS  =========================
+
+            /// Number of cells in the x dimension
+            int nx;
+
+            /// Number of cells in the y dimension
+            int ny;
+
+            /// Number of cells in the z dimension
+            int nz;
+
+            /// Grid offset in x dimension
+            Real ox;
+
+            /// Grid offset in y dimension
+            Real oy;
+
+            /// Grid offset in z dimension
+            Real oz;
+
+            /// Cell spacing in the x dimension
+            VectorXr dx;
+
+            /// Cell spacing in the y dimension
+            VectorXr dy;
+
+            /// Cell spacing in the z dimension
+            VectorXr dz;
+
+        private:
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "RectilinearGrid";
+
+    }; // -----  end of class  RectilinearGrid  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef RECTILINEARGRID_INC  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/RectilinearGridReader.h

@@ -0,0 +1,110 @@
+/* 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
+  @date      09/23/2013 11:05:11 AM
+  @version   $Id$
+  @author    Trevor Irons (ti)
+  @email     Trevor.Irons@xri-geo.com
+  @copyright Copyright (c) 2013, XRI Geophysics, LLC
+  @copyright Copyright (c) 2013, Trevor Irons
+ **/
+
+#ifndef  RECTILINEARGRIDREADER_INC
+#define  RECTILINEARGRIDREADER_INC
+
+#include "GridReader.h"
+#include "RectilinearGrid.h"
+#include "ASCIIParser.h"
+
+namespace Lemma {
+
+    /**
+      @class
+      \ingroup LemmaCore
+      \brief   Base Class for rectilinear grid reader data files
+      \details
+     \todo this set of classes is deprecated, use serialize methods!
+     */
+    class RectilinearGridReader : public GridReader {
+
+        friend std::ostream &operator<<(std::ostream &stream,
+                const RectilinearGridReader &ob);
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+            /**
+             *  Factory method for generating concrete class.
+             *  @return a std::shared_ptr of type RectilinearGridReader
+             */
+            static std::shared_ptr< RectilinearGridReader > NewSP();
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            /**
+            Reads a .mod file, a simple representation of the grid. These files
+            have the following format:
+            @verbatim
+            nx ny nz  // number of cells in x, y, and z
+            ox oy oz  // grid offset of top SW corner
+
+            hx[0] hx[1] ... hx[nx-1]  // grid cell spacing in x, from 0 to nx-1
+            hy[0] hy[1] ... hy[ny-1]  // grid cell spacing in y, from 0 to ny-1
+            hz[0] hz[1] ... hz[nz-1]  // grid cell spacing in z, from 0 to nz-1
+
+            // c++ style (//) comments are allowed anywhere, c-style are not
+            @endverbatim
+            */
+            void ReadASCIIGridFile( const std::string& name );
+
+            // ====================  ACCESS        =======================
+            /**
+             * Accessor method for the underlying RectilinearGrid class constructed
+             * from the input file.
+             */
+            std::shared_ptr<Grid> GetGrid();
+
+            // ====================  INQUIRY       =======================
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /** Default protected constructor, use New */
+            RectilinearGridReader ( );
+
+            /** Default protected constructor, use Delete */
+            ~RectilinearGridReader ();
+
+        private:
+
+            // ====================  DATA MEMBERS  =========================
+
+            /** ASCII string representation of the class name */
+            static constexpr auto CName = "RectilinearGridReader";
+
+            /** Object holding constructed Grid */
+            std::shared_ptr<RectilinearGrid> rGrid;
+
+            /** Performs actual file parsing */
+            std::shared_ptr<ASCIIParser> Parser;
+
+    }; // -----  end of class  RectilinearGridReader  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef RECTILINEARGRIDREADER_INC  -----

Modules/LemmaCore/include/RectilinearGridVTKExporter.h

@@ -0,0 +1,124 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/25/2013 08:20:01 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2013, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2013, Trevor Irons
+ */
+
+
+#ifndef  RECTILINEARGRIDVTKEXPORTER_INC
+#define  RECTILINEARGRIDVTKEXPORTER_INC
+
+#ifdef  LEMMAUSEVTK
+
+#include "LemmaObject.h"
+#include "RectilinearGrid.h"
+
+#include <vtkSmartPointer.h>
+#include "vtkXMLRectilinearGridWriter.h"
+#include "vtkRectilinearGrid.h"
+#include "vtkDoubleArray.h"
+#include "vtkFieldData.h"
+#include "vtkCellData.h"
+
+namespace Lemma {
+
+/**
+  \ingroup LemmaCore
+  \brief   Exports a LemmaRectilinearGrid into a VTKRectilinearGrid.
+  \details Results can then be dumped to a VTK file or used directly.
+ */
+class RectilinearGridVTKExporter : public LemmaObject {
+
+    friend std::ostream &operator<<(std::ostream &stream,
+            const RectilinearGridVTKExporter &ob);
+
+    public:
+
+    // ====================  LIFECYCLE     =======================
+
+    /**
+     * @copybrief LemmaObject::New()
+     * @copydetails LemmaObject::New()
+     */
+    static RectilinearGridVTKExporter* New();
+
+    /**
+     *  Factory method for generating concrete class.
+     *  @return a std::shared_ptr of type RectilinearGridVTKExporter
+     */
+    static std::shared_ptr< RectilinearGridVTKExporter > NewSP();
+
+    // ====================  OPERATORS     =======================
+
+    // ====================  OPERATIONS    =======================
+
+    // ====================  ACCESS        =======================
+    /**
+     * Sets the Grid to be exported
+     */
+    void SetGrid( std::shared_ptr<RectilinearGrid> Grid);
+
+    // ====================  INQUIRY       =======================
+
+    /**
+     *  @return Returns a pointer to a vtkRectilinearGrid containing same information as the Lemma Grid.
+     */
+    vtkSmartPointer<vtkRectilinearGrid>  GetVTKGrid();
+
+    /**
+     *  Writes the exported VTKGrid out to a .vtr XML file
+     *  @param[in] name is the filename to write to. The .vtr extension will be added
+     */
+    void WriteVTKGrid( const std::string& name );
+
+    /** Returns the name of the underlying class, similiar to Python's type */
+    virtual inline std::string GetName() const {
+        return this->CName;
+    }
+
+    protected:
+
+    // ====================  LIFECYCLE     =======================
+
+    /** Default protected constructor, use New */
+    RectilinearGridVTKExporter ( );
+
+    /** Default protected destructor, use Delete */
+    ~RectilinearGridVTKExporter ();
+
+    private:
+
+    /** ASCII string representation of the class name */
+    static constexpr auto CName = "RectilinearGridVTKExporter";
+
+    void BuildVTKRectilinearGrid();
+
+    // ====================  DATA MEMBERS  =========================
+
+    /** The Grid to convert from */
+    std::shared_ptr<RectilinearGrid>    Grid;
+
+    /** The Grid to convert to */
+    vtkSmartPointer<vtkRectilinearGrid> VTKGrid;
+
+}; // -----  end of class  RectilinearGridVTKExporter  -----
+
+}		// -----  end of Lemma  name  -----
+
+#endif     // -----  not LEMMAUSEVTK  -----
+
+#endif   // ----- #ifndef RECTILINEARGRIDVTKEXPORTER_INC  -----
+

Modules/LemmaCore/include/WindowFilter.h

@@ -0,0 +1,127 @@
+/* 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     07/20/2010
+  @version  $Id: windowfilter.h 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#ifndef  WINDOWFILTER_INC
+#define  WINDOWFILTER_INC
+#include "Filter.h"
+
+namespace Lemma {
+    /** Types of filter window that are supported.
+     */
+    enum WINDOWTYPE { HAMMING, /*!< A hamming window */
+                      HANNING, /*!< A hanning window */
+                      RECTANGULAR /*!< Rectangular window */
+                    };
+
+    // ===================================================================
+    //  Class:  WindowFilter
+    /**
+      @class   WindowFilter
+      \ingroup LemmaCore
+      \brief   Fourier domain window filter.
+      \details Current types are Hamming and Hanning. Others to be added.
+     */
+    // ===================================================================
+    class WindowFilter : public Filter {
+
+        public:
+
+            // ====================  LIFECYCLE     =======================
+
+           /**
+            *  Factory method for generating concrete class.
+            *  @return a std::shared_ptr of type WindowFilter
+            */
+            static std::shared_ptr< WindowFilter > NewSP();
+
+            // ====================  OPERATORS     =======================
+
+            // ====================  OPERATIONS    =======================
+
+            /// Sets the bandwith of the window.
+            /// @param[in] width Bandwidth in Hz
+            void SetBandwidth(const Real& width);
+
+            /// Sets the sampling rate of an input time domain signal.
+            /// @param[in] dt is the sampling rate, in seconds.
+            void SetSamplingInterval(const Real& dt);
+
+            /// Specifies the length of the time series to apply filter to.
+            void SetNumberOfSamples(const int& nt);
+
+            /// Sets the type of window to compute.
+            void SetWindowType(const WINDOWTYPE& type);
+
+            /** Returns the number of frequency bins for whole data record.
+             *  this is different then the length of the filter. This is
+             *  simply the length of the real to complex FFT, or
+             *  Nw = Nt%2 ? (Nt-1)/2+1 : (Nt)/2+1;
+             *  @return the number of frequency bins.
+             */
+            int GetNumberOfFrequencyBins();
+
+            // ====================  ACCESS        =======================
+
+            // ====================  INQUIRY       =======================
+
+            /// Returns a Vector of the filter coefficients.
+            VectorXr GetFilterCoefficients( );
+
+            /** Returns the name of the underlying class, similiar to Python's type */
+            virtual inline std::string GetName() const {
+                return this->CName;
+            }
+
+        protected:
+
+            // ====================  LIFECYCLE     =======================
+
+            /// Default protected constructor.
+            WindowFilter ( );
+
+            /// Default protected constructor.
+            ~WindowFilter ( );
+
+            // ====================  DATA MEMBERS  =========================
+
+            /// Width of the window
+            int         Width;
+
+            /// Length of time series
+            int         Nt;
+
+            /// Length of the frequency series representation
+            int         Nw;
+
+            /// Sampling rate of a time domain signal
+            Real        SamplingRate;
+
+            /// Bandwith of filter
+            Real        Bandwidth;
+
+            /// The filter coefficients
+            VectorXr    Coefficients;
+
+            /// The type of filter to use
+            WINDOWTYPE  Type;
+
+        private:
+
+        /** ASCII string representation of the class name */
+        static constexpr auto CName = "WindowFilter";
+
+    }; // -----  end of class  WindowFilter  -----
+
+}       // -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef WINDOWFILTER_INC  -----

Modules/LemmaCore/include/banner.h

@@ -0,0 +1,24 @@
+/* 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   M. Andy Kass
+  @date     05/03/2010
+  @version  $Id: banner.h 87 2013-09-05 22:44:05Z tirons $
+ **/
+
+#ifndef __BANNER_H
+#define __BANNER_H
+#include <iostream>
+#include <ostream>
+void banner(const std::string &name, const std::string &version,
+		const std::string &usage);
+#endif
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/include/helper.h

@@ -0,0 +1,288 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      10/02/2014 02:49:55 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#pragma once
+
+#ifndef  HELPER_INC
+#define  HELPER_INC
+
+#include "lemma.h"
+#include "yaml-cpp/yaml.h"
+
+namespace Lemma {
+
+/** \addtogroup LemmaCore
+ * @{
+ */
+
+    /**
+     *  Convenience function for string conversion
+     *  @param[in] t input value to be converted to string
+     *  @return string representation of input value
+     */
+    template <class T>
+    inline std::string to_string (const T& t) {
+        std::stringstream ss;
+        ss << t;
+        return ss.str();
+    }
+
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const FREQUENCYUNITS& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const TIMEUNITS& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const MAGUNITS& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const TEMPUNITS& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const FEMCOILORIENTATION& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const ORIENTATION& Units);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const FIELDCOMPONENT& Comp);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const HANKELTRANSFORMTYPE& Htype);
+    /// convert enums to string saves repeated code useful for YAML serializing
+    std::string enum2String(const FIELDCALCULATIONS& Htype);
+
+    // other way around is a template, where template argument lets us know
+    // which specialisation to use.
+    template <typename T>
+    T string2Enum( const std::string& str );
+
+    // Handy little class that indents a stream.
+    // Based on solution provided here, todo may need to add to some managing class which keeps
+    // track of nesting levels? But perhaps not. A Lemma class will contain pointers to other Lemma
+    // classes. But those are not specifically listed out.
+    // http://stackoverflow.com/questions/9599807/how-to-add-indention-to-the-stream-operator
+    class IndentingOStreambuf : public std::streambuf {
+        std::streambuf*     myDest;
+        bool                myIsAtStartOfLine;
+        std::string         myIndent;
+        std::ostream*       myOwner;
+    protected:
+        virtual int         overflow( int ch )
+        {
+            if ( myIsAtStartOfLine && ch != '\n' ) {
+                myDest->sputn( myIndent.data(), myIndent.size() );
+            }
+            myIsAtStartOfLine = ch == '\n';
+            return myDest->sputc( ch );
+        }
+    public:
+        explicit            IndentingOStreambuf(
+                                std::streambuf* dest, int indent = 4 )
+            : myDest( dest )
+            , myIsAtStartOfLine( true )
+            , myIndent( indent, ' ' )
+            , myOwner( NULL )
+        {
+        }
+        explicit            IndentingOStreambuf(
+                                std::ostream& dest, int indent = 4 )
+            : myDest( dest.rdbuf() )
+            , myIsAtStartOfLine( true )
+            , myIndent( indent, ' ' )
+            , myOwner( &dest )
+        {
+            myOwner->rdbuf( this );
+        }
+        virtual             ~IndentingOStreambuf()
+        {
+            if ( myOwner != NULL ) {
+                myOwner->rdbuf( myDest );
+            }
+        }
+    };
+
+/** @}*/
+
+} // end namespace Lemma
+
+///////////////////////////////////////////////////////
+// YAML Serializing helper functions. Can we move this into helper.h
+namespace YAML {
+
+template<>
+struct convert<Lemma::Complex> {
+  static Node encode(const Lemma::Complex& rhs) {
+    Node node;
+    node["real"] = rhs.real();
+    node["imag"] = rhs.imag();
+
+    // No labels
+    //node.push_back(rhs.real());
+    //node.push_back(rhs.imag());
+
+    node.SetTag( "Complex" ); // too verbose?
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::Complex& rhs) {
+    // Disabled due to overly verbose output. Just believe...
+    if( node.Tag() != "Complex" ) {
+        return false;
+    }
+    rhs = Lemma::Complex( node["real"].as<Lemma::Real>(), node["imag"].as<Lemma::Real>()  );
+    // no label style
+    //rhs = Lemma::Complex( node[0].as<Lemma::Real>(), node[1].as<Lemma::Real>()  );
+    return true;
+  }
+
+};
+
+template<>
+struct convert<Lemma::Vector3Xr> {
+  static Node encode(const Lemma::Vector3Xr& rhs) {
+    Node node;
+    node["size"] = rhs.cols();
+    //node["rows"] = rhs.rows(); // == 3
+    for (int ic=0; ic<rhs.cols(); ++ic) {
+        node[ic].push_back( rhs(0, ic) );
+        node[ic].push_back( rhs(1, ic) );
+        node[ic].push_back( rhs(2, ic) );
+    }
+    node.SetTag( "Vector3Xr" );
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::Vector3Xr& rhs) {
+    if( node.Tag() != "Vector3Xr" ) {
+        return false;
+    }
+    rhs.resize( Eigen::NoChange, node["size"].as<int>() );
+    for (unsigned int ic=0; ic<node.size(); ++ic) {
+        int ir=0;
+        for(YAML::const_iterator it=node[ic].begin();it!=node[ic].end();++it) {
+            rhs(ir, ic) = it->as<Lemma::Real>();
+            ++ir;
+        }
+    }
+    return true;
+  }
+};
+
+template<>
+struct convert<Lemma::VectorXr> {
+  static Node encode(const Lemma::VectorXr& rhs) {
+    Node node;
+    node["size"] = rhs.size();
+    for (int ic=0; ic<rhs.size(); ++ic) {
+        node["data"].push_back( rhs(ic) );
+    }
+    node.SetTag( "VectorXr" );
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::VectorXr& rhs) {
+    if( node.Tag() != "VectorXr" ) {
+        return false;
+    }
+    rhs.resize( node["size"].as<int>() );
+    int ir=0;
+    for(YAML::const_iterator it=node["data"].begin(); it!=node["data"].end(); ++it) {
+        rhs(ir) = it->as<Lemma::Real>();
+        ++ir;
+    }
+    return true;
+  }
+};
+
+template<>
+struct convert<Lemma::VectorXcr> {
+  static Node encode(const Lemma::VectorXcr& rhs) {
+    Node node;
+    node["size"] = rhs.size();
+    for (int ic=0; ic<rhs.size(); ++ic) {
+        node["data"].push_back( rhs(ic) );
+    }
+    node.SetTag( "VectorXcr" );
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::VectorXcr& rhs) {
+    if( node.Tag() != "VectorXcr" ) {
+        return false;
+    }
+    rhs.resize( node["size"].as<int>() );
+    int ir=0;
+    for(YAML::const_iterator it=node["data"].begin(); it!=node["data"].end(); ++it) {
+        rhs(ir) = it->as<Lemma::Complex>();
+        ++ir;
+    }
+    return true;
+  }
+};
+
+
+template<>
+struct convert<Lemma::VectorXi> {
+  static Node encode(const Lemma::VectorXi& rhs) {
+    Node node;
+    node["size"] = rhs.size();
+    for (int ic=0; ic<rhs.size(); ++ic) {
+        node["data"].push_back( rhs(ic) );
+    }
+    node.SetTag( "VectorXi" );
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::VectorXi& rhs) {
+    if( node.Tag() != "VectorXi" ) {
+        return false;
+    }
+    rhs.resize( node["size"].as<int>() );
+    int ir=0;
+    for(YAML::const_iterator it=node["data"].begin(); it!=node["data"].end(); ++it) {
+        rhs(ir) = it->as<int>();
+        ++ir;
+    }
+    return true;
+  }
+};
+
+template<>
+struct convert<Lemma::Vector3r> {
+  static Node encode(const Lemma::Vector3r& rhs) {
+    Node node;
+    for (int ic=0; ic<rhs.size(); ++ic) {
+        node[0].push_back( rhs(ic) );
+    }
+    node.SetTag( "Vector3r" );
+    return node;
+  }
+
+  static bool decode(const Node& node, Lemma::Vector3r& rhs) {
+    if( node.Tag() != "Vector3r" ) {
+        return false;
+    }
+    int ir=0;
+    for(YAML::const_iterator it=node[0].begin(); it!=node[0].end(); ++it) {
+        rhs(ir) = it->as<Lemma::Real>();
+        ++ir;
+    }
+    return true;
+  }
+};
+
+}
+
+#endif   // ----- #ifndef HELPER_INC  -----
+

Modules/LemmaCore/include/lemma.h

@@ -0,0 +1,333 @@
+/* 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     06/19/2009 09:12:20 AM  The Birth of Lemma!
+  @version  $Id: lemma.h 203 2015-01-09 21:19:04Z tirons $
+ **/
+
+/** \mainpage Lemma is an ElectroMagnetics Modelling API
+    \image html lemma.png
+
+    \authors Trevor Irons and M. Andrew Kass and others
+
+    Lemma is a recursive acronym that stands for <B>L</B>emma is an
+    <B>E</B>lectro<B>M</B>agnetics <B>M</B>odelling <B>A</B>PI.
+    Lemma is a cross-platform library delivering an expressive API that
+    can be used to easily create versatile programs.
+    We are building a powerful, flexible,
+    expresive framework that allows for the straightforward creation of
+    geophysical EM applications. Lemma is not itself a program, instead
+    it is a collection of building blocks to make applications.
+    We chose this name because:
+
+    - In mathematics a Lemma is a proven proposition which is used as a
+    stepping stone to a larger result rather than as a statement in-and-of
+    itself.
+    - In addition to the electromagnetic modelling, some other facilities are
+    provided such as numerical optimization and inversion capabilities. These
+    tools are also considered stepping stones to final products.
+
+    We feel that this is a partucularily approprate name, as Lemma's
+    API can be leveraged create powerful applications such as forward
+    modelling and inverting frequency and time-domain
+    surveys of arbitrary survey design, SNMR surveys, CSAMT and more.
+
+    \section Motivation
+    Why another Geophysical EM project? For starters, there
+    aren't that many quality open source packages out there. Those that do
+    exist are generally specialized to perform a single task and extending
+    them is a major undertaking. Lemma's approach is much different, by
+    providing a set of general tools users can easily assemble applications
+    that suite their needs.
+
+    While most geophysical EM modeling packages are concerned with forward
+    modeling <B>DATA</B> between specific transmitters and receivers located
+    on the surface of the earth; Lemma is fundamentally concerned with
+    modelling <B>FIELDS</B> due to sources, the field values can be calculated
+    everywhere within the domain. This allows for solving multiphysics EM
+    problems. Tradional EM data may be forward modelled using the Lemma::Instrument
+    classes.
+
+    \section Capabilities Capabilities
+    In the long term, we have many goals for this software project. Due to its
+    design, Lemma can be built upon and extended easily. The initial aim is to
+    provide flexible 1D and 3D EM modelling in the time and frequency domains.
+    The project is still in beta, but we have made a lot of progress already.
+    We will release our first non-beta release as soon as the following are
+    supported.
+
+    \subsection FDM Frequency-domain forward modelling
+    Lemma was initially called EMMODFD: Electromagnetic Modelling in the Frequency
+    Domain. As such this is the most mature area of Lemma.
+
+    \par 1D
+    Frequency domain solutions to electrical and magnetic dipoles can be computed
+    quasi-analytically in 1D. Calculations can be made in or above the layered
+    media, and complex electrical conductivity and magnetic susceptibility are
+    supported according to the Cole-Cole model. Sources may be embdedded in the
+    media or in the resisitive air layer. Lemma can also can compute fields due
+    to arbitrarily shaped ungrounded wire loops, topography of the loops is also
+    supported. Two separate approaches to solving the Hankel transform, one
+    based on Anderson's digitial filtering technique, and another based on Gaussian
+    quadrature.
+
+    \par 3D
+    A fast 3D solver that can modify the 1D results based
+    on arbitrary electrical conductivity model is nearing completion.
+
+    \par future work
+    We are also planning on supporting grounded wires in the near future.
+
+    \subsection TDM Time-domain forward modelling
+	A 1D time-domain solution has been implemented that utilises both a
+	dipole source as well as a wire loop.  Currently, only one receiver is
+	modelled at a time, but will be generalised.  In addition, utilities
+	to read in data files for modelling have been implemented.
+    We would like to offer 3D time domain support, but this will not be
+    provided before our first stable release.
+
+    \subsection DataFormats Data Formats
+    The EM community is plagued with myriad data formats. Often each equiptment
+    manufacturer provides their own data format and interoperability is a
+    constant struggle. We are working on a flexible data format based on the XML
+    format that can be adapted to many types of data. The template for this
+    format will be publically released and we hope it catches on in the community.
+    At the least, it will provide a mechanism to compare datasets and datatypes
+    within Lemma.
+
+    \section Modules Modules
+    Due to Lemma's design, it is easy to extend the platform. In some cases this
+    extension results in adding functionality that is not directly related to
+    ElectroMagnetics. The following modules utilise parts of Lemma to provide
+    their functionality.
+
+    \subsection Merlin Merlin.
+    This module provides nuclear magnetic resonance functionality.
+
+
+    \section Tutorials
+        - \ref Tutorial - Basic intruduction to Lemma, including aquiring and
+          compiling the code, class structure, and building your own
+          applications.
+        - \ref Extending Tutorial on how to extend Lemma.
+
+
+    \section Development Development and design
+    Ths package is primarily being developed by the Center for Gravity,
+    Electrical, and Magnetic Studies (CGEM) at the Colorado School of Mines (CSM), the United
+    States Geological Survey (USGS), and Broken Spoke Development, LLC. It draws on work by
+    many others including Ki Ha Lee, and Walt Anderson. All new work and
+    interfaces are written entirely in C++. Several small external projects are
+    included, which are written in standard C, and FORTRAN 77. We adapt a
+    modern, test driven, object oriented, C++ framework.
+
+    \section Legalities
+
+    \subsection Copyrights
+    The following copyrights apply to the source.
+    Most of the code was developed either by Trevor
+    Copyright (C) 2008-2010 Trevor Irons <tirons@mines.edu> or
+    M. Andrew Kass Copyright (C) 2010 <mkass@numericalgeo.com>.
+
+    The 1D EM solver was derived (but updated heavily) from a fortran
+    programme written by Ki Ha Lee in 1984. We have communicated with Ki Ha,
+    and he assured us that this code is in the public domian.
+
+    A Gaussian quadrature hankel transform originally written by Alan Chave was
+    ported to C++. This code is in the public domain, and was published in Geophysics.
+
+    A digital filtering approach to the Hankel transform written by Walt
+    Anderson was also rewritten for Lemma. The origional fortran code is also in
+    the public domain. Using this Hankel transform works
+    arround the copyright issue mentioned above, but may not produce accurate
+    results for high frequencies.
+
+    Please note that Ki Ha Lee and Walt Anderson had no part in this work, and
+    the above should not be interpreted as any sort of endorsement by those
+    parties.
+
+    \subsection License
+
+    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/.
+
+    \section Contributing Suggestions and contributions
+        We welcome contributions and suggestions. Feel free to email the development
+        team at info@lemmasoftware.org.
+        Under the terms of the MPL, if you modify a Lemma file, you are obligated to
+        share those contributions back with the community.
+
+    \section Useful Useful links
+        - Home page https://lemmasoftware.org
+        - SVN repository https://svn.lemmasoftware.org
+        - Broken Spoke Develpment http://numericalgeo.com
+        - CGEM at the Coloroado School of Mines http://geophysics.mines.edu/cgem/
+**/
+
+#pragma once
+
+#ifndef __LEMMA_H
+#define __LEMMA_H
+
+    // Include some basic stuff that will always be needed
+    #include <iostream>
+    #include <iomanip>
+    #include <complex>
+    #include <fstream>
+    #include <string>
+    #include <vector>
+    #include <stdexcept>
+    #include <sstream>
+
+    #include <Eigen/Core>
+    #include <cstddef>
+    #include <Eigen/StdVector>
+    #include <Eigen/Sparse>
+    #include <unsupported/Eigen/FFT>
+    //#include <unsupported/Eigen/SparseExtra>
+    #include <Eigen/Geometry>
+
+    /** \brief The only namespace used by Lemma
+     *
+     *  \details The rational behind this namespace is that built-in
+     *  types should be used wherever possible, but not
+     *  not built-in names. This allows for code that is better
+     *  enacsulated and easier to modify. The typedefs and constants
+     *  specified here are defined so that
+     *  precision/inplimentation can easily be changed.
+     *  All floating precision types should be typedefed in this file
+     *  and should not be used natively within any code.
+     *  Lemma uses
+     *  the Eigen Matrix/Vector/Linear Algebra library.
+     *  <http://eigen.tuxfamily.org> and a lot of the namespece typedefs
+     *  are specifying Eigen types.
+     */
+    namespace Lemma {
+
+        /// Real defines precision for the whole API, default is double
+#ifdef  LEMMA_SINGLE_PRECISION
+        typedef float Real;
+#else      // -----      LEMMA_SINGLE_PRECISION  -----
+        typedef double Real;
+#endif     // -----  not LEMMA_SINGLE_PRECISION  -----
+
+        /// Complex version of Real.
+        typedef std::complex<Real> Complex;
+
+        /// A 3 component Eigen vector of Reals
+        typedef Eigen::Matrix<Real, 3, 1> Vector3r;
+
+        /// A 3 X Dynamic Component Eigen matrix of Reals
+        typedef Eigen::Matrix<Real, 3, Eigen::Dynamic> Vector3Xr;
+
+        /// Variable length Eigen vector of Reals
+        typedef Eigen::Matrix<Real, Eigen::Dynamic, 1> VectorXr;
+
+        /// Variable length Eigen vector of integers (int)
+        typedef Eigen::Matrix<int, Eigen::Dynamic, 1> VectorXi;
+
+        /// Variable length Eigen vector of Complexes
+        typedef Eigen::Matrix<Complex, Eigen::Dynamic, 1> VectorXcr;
+
+        /// A 3 Component Eigen vector of Complexes
+        typedef Eigen::Matrix<Complex, 3, 1> Vector3cr;
+
+        /// A 3 X Dynamic Component Eigen matrix of Complexes
+        typedef Eigen::Matrix<Complex, 3, Eigen::Dynamic> Vector3Xcr;
+
+        /// Variable length Eigen Matrix of Reals
+        typedef Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic> MatrixXr;
+
+        /// Variable length Eigen Matrix of ints
+        typedef Eigen::Matrix<int, Eigen::Dynamic, Eigen::Dynamic> MatrixXi;
+
+        /// Variable length Eigen vector of Complexes
+        typedef Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic> MatrixXcr;
+
+        ////////////////////////////////////////
+        // Constants used across the programmes
+
+        /// Restating the obvious, this is pi
+        const Real PI = 4.0*atan(1.0);
+
+        /// Permitivity of Free Space
+        //const Real EPSILON0 = 8.854187817e-12;
+        const Real EPSILON0 = 8.854187817e-12;
+
+        /// Permeability of free space
+        const Real MU0 = 4.*PI*1e-7;
+
+        /// 1/4 of \f$ \pi\f$
+        const Real QPI = .25/PI;
+
+        /// Some functions will convert units from SI (standard) to Gauss
+        /// This is because NMR calculations are much more natural in Gauss
+        enum MAGUNITS  {TESLA, NANOTESLA, GAUSS};
+
+        /// Unit of temperature entered
+        enum TEMPUNITS {CELCIUS, KELVIN};
+
+        /// Unit of time entered
+        enum TIMEUNITS {SEC, MILLISEC, MICROSEC, NANOSEC, PICOSEC};
+
+        /// Unit of time entered
+        enum FREQUENCYUNITS {HZ, KHZ, MHZ, GHZ};
+
+		/// FEM coil relative orientations
+		enum FEMCOILORIENTATION {COAXIAL, COPLANAR};
+
+		/// General orientation relative to coordinate system
+		enum ORIENTATION {X, Y, Z, NX, NY, NZ};
+
+        /// Type of field
+        enum FIELDTYPE         {HFIELDREAL, HFIELDIMAG, EFIELDREAL, EFIELDIMAG};
+
+        /// Compenent of vector field
+        enum FIELDCOMPONENT    {XCOMPONENT=0, YCOMPONENT=1, ZCOMPONENT=2};
+
+        /// Spatial component of vector
+        enum SPATIALCOORDINANT  {XCOORD=0, YCOORD=1, ZCOORD=2};
+
+        /** Evaluation method for Hankel integrals.
+         *  ANDERSON801     Walt Anderson's 801 point filter
+         *  CHAVE           Alan Chave's gaussian quadrature integration method
+         *  FHTKEY201       Key's 201 point filter
+         *  FHTKEY201       Key's 101 point filter
+         *  FHTKEY51        Key's 51 point filter
+         *  QWEKEY          Key's Gaussian quadrature integration method
+         */
+        enum HANKELTRANSFORMTYPE { ANDERSON801, CHAVE, FHTKEY201, FHTKEY101, FHTKEY51, QWEKEY };
+
+        /** Enum is OK because these are the only physically possible sources.
+         @param NOSOURCETYPE is default.
+         @param ELECTRICDIPOLE is an electric dipole
+         @param MAGNETICDIPOLE is a magnetic dipole
+        */
+        enum DipoleSourceType {NOSOURCETYPE, GROUNDEDELECTRICDIPOLE, UNGROUNDEDELECTRICDIPOLE, MAGNETICDIPOLE};
+
+        /// Only three polarizations are supported. They may be summed to
+        /// approximate others
+        /// @param NOPOLARISATION is uninitialized, default value
+        /// @param XPOLARISATION is a dipole oriented in the x direction
+        /// @param YPOLARISATION is a dipole oriented in the y direction
+        /// @param ZPOLARISATION is a dipole oriented in the z direction
+        enum DipoleSourcePolarisation {NOPOLARISATION, XPOLARISATION,
+            YPOLARISATION,  ZPOLARISATION};
+
+        /// The polarity may be either negative or positinve
+        enum DipoleSourcePolarity {NEGATIVE, POSITIVE};
+
+        /** The fields to make calculations on
+        */
+        enum FIELDCALCULATIONS {E, H, BOTH};
+
+    }
+
+#endif // __Lemma_H

Modules/LemmaCore/include/timer.h

@@ -0,0 +1,82 @@
+/* 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     01/27/2010
+  @version  $Id: timer.h 87 2013-09-05 22:44:05Z tirons $
+ **/
+
+#ifndef  TIMER_H_INC
+#define  TIMER_H_INC
+
+#include <iostream>
+#include <fstream>
+#include <ctime>
+#ifdef LEMMAUSEOMP
+#include "omp.h"
+#endif
+
+#include "lemma.h"
+
+namespace Lemma {
+
+/** \brief Simple timer class.
+  * \details Can be used either with multithreaded applications, using
+  * OpenMP defined functions, or using system wall time. System wall time
+  * is not good for parallel programs as it reports cpu core execution time.
+  */
+class jsw_timer {
+
+    public:
+
+        /**
+          * Constructor that automatically starts the timer.
+          */
+        jsw_timer(): start ( std::clock() ), elapsed ( 0.0 ) {
+
+        }
+
+        /** begins the timer.
+          */
+        void begin() {
+            #ifdef LEMMAUSEOMP
+            start =  omp_get_wtime();
+            #else
+            start = std::clock();
+            #endif
+            elapsed = 0.0;
+        }
+
+        /** Ends the timer.
+          * @return the time that the stopwatch ended at.
+          */
+        Real end() {
+
+            #ifdef LEMMAUSEOMP
+            elapsed = static_cast<Real> (omp_get_wtime()) - start;
+            return elapsed;
+            #else
+            elapsed = static_cast<Real> ( std::clock() ) - start;
+            return elapsed /= CLOCKS_PER_SEC;
+            #endif
+
+        }
+
+        /**
+          * @return the elapsed time between start and end.
+          */
+        Real last() const { return elapsed; }
+
+    private:
+        std::clock_t start;
+        Real    elapsed;
+};
+
+}       // -----  end of Lemma  name  -----
+
+#endif   // ----- #ifndef TIMER_H_INC  -----

Modules/LemmaCore/src/ASCIIParser.cpp

@@ -0,0 +1,238 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/23/2013 02:33:41 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2013, Trevor Irons
+ */
+
+#include "ASCIIParser.h"
+
+namespace Lemma {
+
+
+    // ====================  FRIEND METHODS  =====================
+
+    std::ostream &operator << (std::ostream &stream, const ASCIIParser &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     =======================
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ASCIIParser
+    // Description:  constructor (protected)
+    //--------------------------------------------------------------------------------------
+    ASCIIParser::ASCIIParser ( ) : LemmaObject( ),
+            CommentString("//"), BufferSize(255) {
+
+    }  // -----  end of method ASCIIParser::ASCIIParser  (constructor)  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  NewSP()
+    // Description:  public smart pointer factory constructor
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr< ASCIIParser > ASCIIParser::NewSP() {
+        std::shared_ptr<ASCIIParser> sp(new  ASCIIParser( ), LemmaObjectDeleter() );
+        return sp;
+        //return  std::make_shared<ASCIIParser>();
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  Serialize
+    //--------------------------------------------------------------------------------------
+    YAML::Node  ASCIIParser::Serialize (  ) const {
+        YAML::Node node = LemmaObject::Serialize();;
+        node.SetTag( GetName() );
+        node["CommentString"] = CommentString;
+        node["BufferSize"] = BufferSize;
+        return node;
+    }		// -----  end of method ASCIIParser::Serialize  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  DeSerialize
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<ASCIIParser> ASCIIParser::DeSerialize ( const YAML::Node& node  ) {
+        if (node.Tag() != "ASCIIParser") {
+            throw  DeSerializeTypeMismatch( "ASCIIParser", node.Tag());
+        }
+        std::shared_ptr<ASCIIParser> Object(new  ASCIIParser(node), LemmaObjectDeleter() );
+        return Object ;
+    }		// -----  end of method ASCIIParser::DeSerialize  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ~ASCIIParser
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    ASCIIParser::~ASCIIParser () {
+
+    }  // -----  end of method ASCIIParser::~ASCIIParser  (destructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ASCIIParser
+    // Description:  DeSerializing constructor (protected)
+    //--------------------------------------------------------------------------------------
+    ASCIIParser::ASCIIParser (const YAML::Node& node) : LemmaObject(node) {
+        this->CommentString = node["CommentString"].as<std::string>();
+        this->BufferSize = node["BufferSize"].as<int>();
+    }  // -----  end of method ASCIIParser::ASCIIParser  (constructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  Release
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::Release() {
+        delete this;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  Open
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::Open ( const std::string& fname ) {
+        input.open(fname.c_str(), std::ios::in);
+        if (input.fail()) {
+           throw  GenericFileIOError(this, fname);
+        }
+        return ;
+    }		// -----  end of method ASCIIParser::Open  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  Close
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::Close (  ) {
+        input.close();
+        return ;
+    }		// -----  end of method ASCIIParser::Close  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ReadReals
+    //--------------------------------------------------------------------------------------
+    std::vector<Real> ASCIIParser::ReadReals ( const int& nr ) {
+        std::string buf;
+        char *dump = new char[BufferSize];
+        std::vector<Real> vals(0);
+        while (input >> buf) {
+            if (buf.substr(0, CommentString.size()) == CommentString) {
+                input.getline(dump, BufferSize);
+            } else {
+                vals.push_back( atof(buf.c_str() ));
+            }
+            if (static_cast<int>(vals.size()) == nr) {
+                delete [] dump;
+                return vals;
+            }
+        }
+        delete [] dump;
+        return vals;
+    }		// -----  end of method ASCIIParser::ReadReals  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ReadInts
+    //--------------------------------------------------------------------------------------
+    std::vector<int> ASCIIParser::ReadInts ( const int& nr ) {
+        std::string buf;
+        char *dump = new char[BufferSize];
+        std::vector<int> vals(0);
+        while (input >> buf) {
+            if (buf.substr(0, CommentString.size()) == CommentString) {
+                input.getline(dump, BufferSize);
+            } else {
+                vals.push_back( atoi(buf.c_str() ));
+            }
+            if (static_cast<int>(vals.size()) == nr) {
+                delete [] dump;
+                return vals;
+            }
+
+        }
+        delete [] dump;
+        return vals;
+    }		// -----  end of method ASCIIParser::ReadInts  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  ReadStings
+    //--------------------------------------------------------------------------------------
+    std::vector<std::string> ASCIIParser::ReadStrings ( const int& nr ) {
+        std::string buf;
+        char *dump = new char[BufferSize];
+        std::vector<std::string> vals(0);
+        while (input >> buf) {
+            if (buf.substr(0, CommentString.size()) == CommentString) {
+                input.getline(dump, BufferSize);
+            } else {
+                vals.push_back( buf );
+            }
+            if (static_cast<int>(vals.size()) == nr) {
+                delete [] dump;
+                return vals;
+            }
+
+        }
+        delete [] dump;
+        return vals;
+    }		// -----  end of method ASCIIParser::ReadInts  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  SetCommentString
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::SetCommentString ( const std::string& key ) {
+        CommentString = key;
+    }		// -----  end of method ASCIIParser::SetCommentString  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  SetBufferSize
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::SetBufferSize ( const int& size ) {
+        BufferSize = size;
+    }		// -----  end of method ASCIIParser::SetCommentString  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  GetFileLocation
+    //--------------------------------------------------------------------------------------
+    int ASCIIParser::GetFileLocation (  ) {
+        return input.tellg();
+    }		// -----  end of method ASCIIParser::GetFileLocation  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  ASCIIParser
+    //      Method:  JumpToLocation
+    //--------------------------------------------------------------------------------------
+    void ASCIIParser::JumpToLocation ( const int& loc ) {
+        input.seekg( loc );
+        return ;
+    }		// -----  end of method ASCIIParser::JumpToLocation  -----
+
+
+}		// -----  end of Lemma  name  -----
+
+

Modules/LemmaCore/src/CMakeLists.txt

@@ -0,0 +1,20 @@
+set (SOURCE 
+	${SOURCE}
+	${CMAKE_CURRENT_SOURCE_DIR}/helper.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/banner.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/ASCIIParser.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/CubicSplineInterpolator.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/Data.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/DataReader.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/EarthModel.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/Instrument.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/LayeredEarth.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/LemmaObject.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/Grid.cpp       
+	${CMAKE_CURRENT_SOURCE_DIR}/RectilinearGrid.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/RectilinearGridReader.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/RectilinearGridVTKExporter.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/Filter.cpp
+	${CMAKE_CURRENT_SOURCE_DIR}/WindowFilter.cpp
+	PARENT_SCOPE
+)

Modules/LemmaCore/src/CubicSplineInterpolator.cpp

@@ -0,0 +1,355 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      02/07/2014 12:50:52 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "CubicSplineInterpolator.h"
+
+#include<iostream>
+#include<fstream>
+#include<vector>
+#include<algorithm>
+#include<cmath>
+
+namespace Lemma {
+
+    // ====================  FRIEND METHODS  =====================
+
+    std::ostream &operator << (std::ostream &stream, const CubicSplineInterpolator &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     =======================
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  CubicSplineInterpolator
+    // Description:  constructor (protected)
+    //--------------------------------------------------------------------------------------
+    CubicSplineInterpolator::CubicSplineInterpolator ( ) : LemmaObject( ) {
+
+    }  // -----  end of method CubicSplineInterpolator::CubicSplineInterpolator  (constructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  CubicSplineInterpolator
+    // Description:  DeSerializing constructor (protected)
+    //--------------------------------------------------------------------------------------
+    CubicSplineInterpolator::CubicSplineInterpolator (const YAML::Node& node) : LemmaObject(node) {
+    }  // -----  end of method CubicSplineInterpolator::CubicSplineInterpolator  (constructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  NewSP()
+    // Description:  public constructor
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<CubicSplineInterpolator> CubicSplineInterpolator::NewSP() {
+        std::shared_ptr<CubicSplineInterpolator> sp(new  CubicSplineInterpolator( ), LemmaObjectDeleter() );
+        return sp;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Serialize
+    //--------------------------------------------------------------------------------------
+    YAML::Node  CubicSplineInterpolator::Serialize (  ) const {
+        YAML::Node node = LemmaObject::Serialize();;
+        node.SetTag( GetName() );
+        // FILL IN CLASS SPECIFICS HERE
+        return node;
+    }		// -----  end of method CubicSplineInterpolator::Serialize  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  DeSerialize
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<CubicSplineInterpolator> CubicSplineInterpolator::DeSerialize ( const YAML::Node& node  ) {
+        if (node.Tag() != "CubicSplineInterpolator") {
+            throw  DeSerializeTypeMismatch( "CubicSplineInterpolator", node.Tag());
+        }
+        std::shared_ptr<CubicSplineInterpolator> Object(new  CubicSplineInterpolator(node), LemmaObjectDeleter() );
+        return Object ;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  ~CubicSplineInterpolator
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    CubicSplineInterpolator::~CubicSplineInterpolator () {
+
+    }  // -----  end of method CubicSplineInterpolator::~CubicSplineInterpolator  (destructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Release
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    void CubicSplineInterpolator::Release() {
+        delete this;
+    }
+
+    // ====================  OPERATIONS    =======================
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  SetKnots
+    //--------------------------------------------------------------------------------------
+    void CubicSplineInterpolator::SetKnots ( const VectorXr& x, const VectorXr& y ) {
+
+        int n = x.size()-1;
+
+        Spline = SplineSet(n);
+
+        Spline.a = y;
+        Spline.x = x;
+        VectorXr h = VectorXr::Zero(n);
+        for(int i=0; i<n; ++i)
+            h(i) = Spline.x(i+1)-Spline.x(i);
+
+        VectorXr alpha(n-1);
+        for(int i=1; i<n-1; ++i)
+            alpha(i) =  3.*(Spline.a[i+1]-Spline.a[i])/h[i] - 3.*(Spline.a[i]-Spline.a[i-1])/h[i-1] ;
+
+        VectorXr l = VectorXr::Zero(n+1);
+        VectorXr mu  = VectorXr::Zero(n+1);
+        VectorXr z = VectorXr::Zero(n+1);
+        l[0] = 1;
+        mu[0] = 0;
+        z[0] = 0;
+
+        for(int i = 1; i < n-1; ++i) {
+            l[i] = 2 *(Spline.x[i+1]-Spline.x[i-1])-h[i-1]*mu[i-1];
+            mu[i] = h[i]/l[i];
+            z[i] = (alpha[i]-h[i-1]*z[i-1])/l[i];
+        }
+        l[n] = 1;
+        z[n] = 0;
+
+        for(int j = n-1; j >= 0; --j) {
+            Spline.c[j] = z[j] - mu[j] * Spline.c[j+1];
+            Spline.b[j] = (Spline.a[j+1]-Spline.a[j])/h[j]-h[j]*(Spline.c[j+1]+2*Spline.c[j])/3;
+            Spline.d[j] = (Spline.c[j+1]-Spline.c[j])/3/h[j];
+        }
+        // On OSX, this causes a strange bug 'sometimes', alignment?
+        //Spline.c = Spline.c.head(n);
+
+        return;
+    }		// -----  end of method CubicSplineInterpolator::SetKnots  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterplator
+    //      Method:  ResetKnotOrdinate
+    //--------------------------------------------------------------------------------------
+    void CubicSplineInterpolator::ResetKnotOrdinate ( const VectorXr& y ) {
+        VectorXr x = Spline.x;
+        SetKnots(x, y);
+        return ;
+    }		// -----  end of method CubicSplineInterplator::ResetKnotOrdinate  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  InterpolateOrderedSet
+    //--------------------------------------------------------------------------------------
+    VectorXr CubicSplineInterpolator::InterpolateOrderedSet ( const VectorXr& x ) {
+        VectorXr y = VectorXr::Zero(x.size());
+        int ii = 0;
+        for (int iy=0; iy<y.size(); ++iy) {
+            y[iy] = Interpolate(x[iy], ii);
+        }
+        return y;
+    }		// -----  end of method CubicSplineInterpolator::InterpolateOrderedSet  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      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()) {
+            ++i;
+        }
+        --i;
+
+        //if ( x > Spline.x[i] ) {
+        //    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  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Interpolate
+    //--------------------------------------------------------------------------------------
+    Real CubicSplineInterpolator::Interpolate ( const Real& x ) {
+        int ii(0);
+        return Interpolate(x, ii);
+    }
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Integrate
+    //--------------------------------------------------------------------------------------
+    Real CubicSplineInterpolator::Integrate ( const Real& x0, const Real& x1, const int& n ) {
+
+        assert(n > 0);
+        // force n to be even?
+        //n += (n % 2);
+
+        Real S = Interpolate(x0) + Interpolate(x1);
+        Real h = (x1 - x0) / Real(n);
+
+        int ik = 0;
+        for (int i=1; i<n; i+=2) {
+            S += 4. * Interpolate(x0 + (Real)(i)*h, ik);
+        }
+
+        ik = 0;
+        for (int i=2; i<n-1; i+=2) {
+            S += 2. * Interpolate(x0 + (Real)(i)*h, ik);
+        }
+
+        return  S * h / 3.;
+    }
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Integrate
+    //--------------------------------------------------------------------------------------
+    Real CubicSplineInterpolator::Integrate ( const Real& x0, const Real& x1 ) {
+
+
+        int i0 = Interval(x0);
+        int i1 = Interval(x1);
+
+        Real h0 = x0 - Spline.x(i0);
+        if (mflag == -1) h0 = 0;
+        Real h1 = x1 - Spline.x(i1);
+
+        Real cubint = (((Spline.d(i1)*h1/4.0 + Spline.c(i1) )*h1/3.0 +
+                         Spline.b(i1) )*h1/2.0 + Spline.a(i1) )*h1
+                    - (((Spline.d(i0)*h0/4.0 + Spline.c(i0) )*h0/3.0 +
+                         Spline.b(i0) )*h0/2.0 + Spline.a(i0) )*h0;
+
+        // Include integrals over intervening intervals.
+        if (i1 > i0) {
+            for (int i=i0; i<i1-1; ++i) {
+                Real h = Spline.x(i+1) - Spline.x(i);
+                cubint +=  (((Spline.d(i)*h/4.0 + Spline.c(i) )*h/3.0 +
+                              Spline.b(i))*h/2.0 + Spline.a(i) )*h;
+            }
+        }
+        return cubint;
+    }		// -----  end of method CubicSplineInterpolator::Integrate  -----
+
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  Interval
+    //--------------------------------------------------------------------------------------
+    int CubicSplineInterpolator::Interval ( const Real& x ) {
+
+        std::cerr << "ENTERING CubicSplineInterpolator::Inverval. Known bugs here"  << std::endl;
+        int nx = Spline.x.size() - 2; // TODO check if this is correct or just -1
+        // when x not in range
+        if (x <= Spline.x(0) || nx <= 1 ) {
+            mflag = -1;
+            return 1;
+        }
+
+        if (x >= Spline.x(nx)) {
+            mflag = 1;
+            return nx;
+        }
+
+        mflag = 0;
+        if (ilo >= nx) ilo = nx-1;
+        int ihi = ilo+1;
+
+        // if x is already in the interval
+        if ( x<Spline.x(ihi) && x >= Spline.x(ilo) ) {
+            //std::cout << "TRIVIAL INTERVAL " << Spline.x(ilo) << "\t" << x << "\t" << Spline.x(ihi) << std::endl;
+            return ilo;
+        }
+
+        if (x <= Spline.x(ilo)) { // decrease ilo to capture
+            int istep = 1;
+            for (int ix=1; ix<nx; ++ix) {
+                ihi = ilo;
+                ilo = ihi - istep;
+                ilo = std::max(1, ilo);
+                if (x >= Spline.x(ilo) || ilo == 1) break;
+                istep *= 2;
+            }
+        } else if (x >= Spline.x(ihi)) { // increase ihi to capture
+            int istep = 1;
+            for (int ix=1; ix<nx; ++ix) {
+                ilo = ihi;
+                ihi = ilo + istep;
+                ihi = std::min(ihi, nx);
+                if (x <= Spline.x(ihi) || ihi == nx) break;
+                istep *= 2;
+            }
+        }
+
+        // Now Spline.x(ilo) <= x < Spline.x(ihi) --> Narrow the interval.
+        //std::cout << "WIDE INTERVAL " << Spline.x(ilo) << "\t" << x << "\t" << Spline.x(ihi) << std::endl;
+        for (int ix=1; ix<nx; ++ix) {
+            int middle = (ilo+ihi) / 2;
+            if (middle == ilo) break;
+            if (x < Spline.x(middle)) {
+                ihi = middle;
+            } else {
+                ilo = middle;
+            }
+        }
+        assert ( Spline.x(ilo) < x && x < Spline.x(ihi) );
+        return ilo;
+    }		// -----  end of method CubicSplineInterpolator::Inverval  -----
+
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterplator
+    //      Method:  GetAbscissa
+    //--------------------------------------------------------------------------------------
+    VectorXr CubicSplineInterpolator::GetKnotAbscissa (  ) {
+        return Spline.x;
+    }		// -----  end of method CubicSplineInterplator::get_GetAbscissa  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  CubicSplineInterpolator
+    //      Method:  GetKnotOrdinate
+    //--------------------------------------------------------------------------------------
+    VectorXr CubicSplineInterpolator::GetKnotOrdinate (  ) {
+        return Spline.a;
+    }		// -----  end of method CubicSplineInterpolator::GetKnotOrdinate  -----
+
+
+}		// -----  end of Lemma  name  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/src/Data.cpp

@@ -0,0 +1,30 @@
+/* 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     03/01/2010
+  @version  $Id: data.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "Data.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const Data &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+	Data::Data(  ) : LemmaObject( ) {
+	}
+
+
+	Data::~Data() {
+	}
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/DataReader.cpp

@@ -0,0 +1,29 @@
+/* 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     03/01/2010
+  @version  $Id: datareader.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "DataReader.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const DataReader &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+	DataReader::DataReader( ) : LemmaObject( ) {
+	}
+
+	DataReader::~DataReader() {
+	}
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/EarthModel.cpp

@@ -0,0 +1,129 @@
+/* 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     03/23/2010
+  @version  $Id: earthmodel.cpp 210 2015-02-25 02:57:03Z tirons $
+ **/
+
+#include "EarthModel.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const EarthModel &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+	// ====================  LIFECYCLE     =======================
+
+	EarthModel::EarthModel() : LemmaObject( ),
+        BField(0,0,0), BFieldUnit(0,0,0), BInc(0), BDec(0), BMag(0) {
+	}
+
+    EarthModel::EarthModel(const YAML::Node& node) : LemmaObject(node) {
+        BInc = node["BInc"].as<double>();
+        BDec = node["BDec"].as<double>();
+        BMag = node["BMag"].as<double>();
+        BField = node["BField"].as<Vector3r>();
+    }
+
+
+    YAML::Node EarthModel::Serialize() const {
+        YAML::Node node = LemmaObject::Serialize();
+        node.SetTag( this->GetName() );
+        node["BField"] = BField;
+        node["BInc"] = BInc;
+        node["BDec"] = BDec;
+        node["BMag"] = BMag;
+        return node;
+    }
+
+	EarthModel::~EarthModel() {
+	}
+
+ 	// ====================  ACCESS        ===================================
+	void EarthModel::SetMagneticFieldComponents(const Vector3r &bfield,
+					const MAGUNITS &unit) {
+
+ 		this->BField = bfield;
+
+		switch (unit) {
+			case TESLA:
+				break;
+			case NANOTESLA:
+				BField *= 1e-9;
+				break;
+		  	case GAUSS:
+				BField *= 1e-4;
+				break;
+			default:
+				throw "MAGUNITS UNDEFINED\n";
+		}
+
+		BMag = bfield.norm(  );
+	    BInc = std::acos (bfield.dot(Vector3r(0,0,1)) / BMag) ;
+	    BDec = std::acos (bfield.dot(Vector3r(1,0,0)) / BMag) ;
+		BFieldUnit = BField.array() / BMag;
+ 	}
+
+	void EarthModel::SetMagneticFieldIncDecMag(const Real &inc,
+			const Real &dec, const Real &mag, const MAGUNITS &unit) {
+		BMag = mag;
+		BInc = inc;
+		BDec = dec;
+
+		switch (unit) {
+			case TESLA:
+				break;
+			case NANOTESLA:
+				BMag *= 1e-9;
+				break;
+		  	case GAUSS:
+				BMag *= 1e-4;
+				break;
+			default:
+				throw "MAGUNITS UNDEFINED\n";
+		}
+
+		BField(0) = BMag * cos(BInc*(PI/180.)) * cos(BDec*(PI/180.));
+		BField(1) = BMag * cos(BInc*(PI/180.)) * sin(BDec*(PI/180.));
+		BField(2) = BMag * sin(BInc*(PI/180.))  ;
+		BFieldUnit = BField.array() / BMag;
+	}
+
+	// ====================  INQUIRY       =======================
+
+	Vector3r EarthModel::GetMagneticField( ) {
+		return this->BField;
+	}
+
+    Vector3r EarthModel::GetMagneticFieldInGauss( ) {
+		return this->BField*1e4;
+	}
+
+	Vector3r EarthModel::GetMagneticFieldUnitVector() {
+		return this->BFieldUnit;
+	}
+
+	Real  EarthModel::GetMagneticFieldMagnitude() {
+		return this->BMag;
+	}
+
+	Real  EarthModel::GetMagneticFieldMagnitudeInGauss() {
+		return this->BMag*1e4;
+	}
+
+    ////////////////////////////////////////////////////////
+    // Error Classes
+
+    NullEarth::NullEarth() :
+        runtime_error( "NULL VALUED LAYERED EARTH MODEL") {}
+
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/Filter.cpp

@@ -0,0 +1,22 @@
+/* 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     07/20/2010
+  @version  $Id: filter.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "Filter.h"
+
+namespace Lemma {
+
+    Filter::Filter( ) : LemmaObject( ) { }
+
+    Filter::~Filter( ) { }
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/Grid.cpp

@@ -0,0 +1,43 @@
+/* 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     10/25/2009
+  @version  $Id: grid.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "Grid.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const Grid &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     ===================================
+
+	Grid::Grid( ) : LemmaObject( ) {
+	}
+
+    Grid::Grid( const YAML::Node& node ) : LemmaObject(node) {
+
+    }
+
+	Grid::~Grid  ( ) {
+	}
+
+    YAML::Node Grid::Serialize() const {
+
+        YAML::Node node = LemmaObject::Serialize();
+        node.SetTag( this->GetName() );
+        return node;
+
+    }
+
+}

Modules/LemmaCore/src/Instrument.cpp

@@ -0,0 +1,29 @@
+/* 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     02/05/2010
+  @version  $Id: instrument.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "Instrument.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const Instrument &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+	Instrument::Instrument( ) : LemmaObject() {
+	}
+
+	Instrument::~Instrument() {
+	}
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/LayeredEarth.cpp

@@ -0,0 +1,130 @@
+/* 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     06/24/2009
+  @version  $Id: layeredearth.cpp 210 2015-02-25 02:57:03Z tirons $
+ **/
+
+#include "LayeredEarth.h"
+
+namespace Lemma {
+
+	// ====================    FRIENDS     ======================
+
+    std::ostream &operator << (std::ostream &stream, const LayeredEarth &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+	// ====================  LIFECYCLE     ===================================
+
+	LayeredEarth::LayeredEarth( ) : EarthModel( ),
+		NumberOfLayers(0), 	NumberOfInterfaces(0) {
+	}
+
+
+
+	LayeredEarth::~LayeredEarth() {
+	}
+
+    LayeredEarth::LayeredEarth(const YAML::Node& node) : EarthModel(node)
+    {
+        NumberOfLayers = node["NumberOfLayers"].as<int>();
+        NumberOfInterfaces = node["NumberOfInterfaces"].as<int>();
+        LayerThickness = node["LayerThickness"].as<VectorXr>();
+	}
+
+    YAML::Node LayeredEarth::Serialize() const {
+        YAML::Node node = EarthModel::Serialize();
+        node["NumberOfLayers"] = NumberOfLayers;
+        node["NumberOfInterfaces"] = NumberOfInterfaces;
+        node["LayerThickness"] = LayerThickness;
+        node.SetTag( this->GetName() );
+        return node;
+    }
+
+	// ====================  OPERATIONS    ===================================
+
+	// ====================  ACCESS        ===================================
+
+	void LayeredEarth::SetLayerThickness(const VectorXr &thick) {
+		if (thick.size() != this->NumberOfLayers - 2)
+			throw EarthModelParametersDoNotMatchNumberOfLayers( );
+		LayerThickness = thick;
+	}
+
+	// ====================  INQUIRY       ===================================
+
+	int LayeredEarth::GetNumberOfLayers () {
+		return this->NumberOfLayers;
+	}
+
+	int LayeredEarth::GetNumberOfNonAirLayers () {
+		return this->NumberOfLayers - 1;
+	}
+
+	Real LayeredEarth::GetLayerThickness(const int &ilay) {
+		// Take into account infinite top and bottom layers
+		// estimate infinity by 1000 m
+		if (ilay < 0 || ilay > NumberOfLayers - 1) {
+			throw RequestForNonValidEarthModelParameter( );
+		} else if (ilay == 0) {
+			return 1000.;
+		} else if (ilay == NumberOfLayers - 1) {
+			return 1000.;
+		} else {
+			return this->LayerThickness(ilay-1);
+		}
+	}
+
+	Real LayeredEarth::GetLayerDepth(const int &ilay) {
+
+		Real depth = 0;
+
+
+		if (ilay == 0) {
+			return depth;
+		} else {
+			for (int i=1; i<=ilay; ++i) {
+				depth += GetLayerThickness(i);
+			}
+		}
+		return depth;
+	}
+
+	int LayeredEarth::GetLayerAtThisDepth(const Real& depth) {
+
+		if (depth <= 0 || NumberOfLayers < 2) {
+			return 0;
+		}
+
+		Real laydep = 0;
+		for (int ilay=0; ilay<NumberOfLayers-2; ++ilay) {
+			laydep += LayerThickness[ilay];
+			if (laydep >= depth) { return ilay+1; }
+		}
+
+		return NumberOfLayers-1;
+	}
+
+    EarthModelWithLessThanTwoLayers::EarthModelWithLessThanTwoLayers() :
+        runtime_error( "EARTH MODEL WITH LESS THAN TWO LAYERS") { }
+
+    EarthModelWithMoreThanMaxLayers::EarthModelWithMoreThanMaxLayers() :
+        runtime_error( "EARTH MODEL WITH MORE THAN MAX LAYERS") { }
+
+    EarthModelParametersDoNotMatchNumberOfLayers::
+        EarthModelParametersDoNotMatchNumberOfLayers( ) :
+        runtime_error( "EARTH MODEL PARAMETERS DO NOT MATCH NUMBER OF LAYERS")
+        {}
+
+    RequestForNonValidEarthModelParameter::
+        RequestForNonValidEarthModelParameter() :
+        runtime_error( "REQUEST FOR NON VALID EARTH MODEL PARAMETER") {}
+}

Modules/LemmaCore/src/LemmaObject.cpp

@@ -0,0 +1,86 @@
+/* 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     06/25/2009
+  @version  $Id: LemmaObject.cpp 266 2015-04-01 03:24:00Z tirons $
+ **/
+
+#include "LemmaObject.h"
+
+namespace Lemma {
+
+    YAML::Emitter& operator << (YAML::Emitter& out, const LemmaObject& ob) {
+        out << YAML::BeginMap;
+        out << YAML::Key <<"Class Name"    << YAML::Value << ob.GetName();
+        return out;
+    }
+
+    // ====================  LIFECYCLE     ==============================
+
+    // Constructor
+    LemmaObject::LemmaObject(  ) {
+
+    }
+
+    LemmaObject::LemmaObject(const YAML::Node &node) {
+
+    }
+
+    // Destructor
+    LemmaObject::~LemmaObject() {
+    }
+
+    // ====================  OPERATIONS    ==============================
+
+    // ====================  INQUIRY       ==============================
+
+    // ====================  ACCESS        ==============================
+
+    // ====================  ACCESS        ==============================
+
+    // ====================  OPERATORS     ==============================
+
+    //////////////////////////////////////////////////////////////////////
+    //////////////////////////////////////////////////////////////////////
+
+    DeSerializeTypeMismatch::DeSerializeTypeMismatch(const std::string& expected, const std::string& got) :
+        runtime_error("DESERIALIZE TYPE MISMATCH") {
+            std::cerr << "Expected " << expected << " got " << got << std::endl;
+        }
+
+    RequestToReturnNullPointer::
+        RequestToReturnNullPointer(LemmaObject *ptr) :
+        runtime_error("REQUEST TO RETURN NULL POINTER"){
+            std::cerr << "Thrown by instance of "
+                      << ptr->GetName() << std::endl;
+        }
+
+    MatFileCannotBeOpened::
+        MatFileCannotBeOpened() :
+        runtime_error("MATFILE CANNOT BE OPENED"){}
+
+    AssignmentOutOfBounds::
+        AssignmentOutOfBounds(LemmaObject *ptr) :
+        runtime_error("ASSIGNMENT OUT OF BOUNDS"){
+            std::cerr << "Thrown by instance of "
+                      << ptr->GetName() << std::endl;
+       }
+
+	GenericFileIOError::
+		GenericFileIOError(LemmaObject *ptr, const std::string &filename) : runtime_error("FILE I/O ERROR"){
+			std::cerr << std::endl;
+			std::cerr << "FILE I/O ERROR" << std::endl;
+			std::cerr << std::endl;
+			std::cerr << "Thrown by instance of "
+				<< ptr->GetName() << std::endl;
+			std::cerr << "  while trying to access " << filename << std::endl;
+			std::cerr << std::endl;
+		}
+
+} // end of namespace Lemma

Modules/LemmaCore/src/RectilinearGrid.cpp

@@ -0,0 +1,215 @@
+/* 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     10/28/2010
+  @version  $Id: rectilineargrid.cpp 193 2014-11-10 23:51:41Z tirons $
+ **/
+
+#include "RectilinearGrid.h"
+
+namespace Lemma {
+
+    std::ostream &operator << (std::ostream &stream, const RectilinearGrid &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     =======================
+
+    RectilinearGrid::RectilinearGrid( ) : Grid( ), nx(0), ny(0), nz(0) {
+
+    }
+
+    RectilinearGrid::RectilinearGrid( const YAML::Node& node ) : Grid(node) {
+
+        nx = node["nx"].as<int>( );
+        ny = node["ny"].as<int>( );
+        nz = node["nz"].as<int>( );
+
+        ox = node["ox"].as<Real>( );
+        oy = node["oy"].as<Real>( );
+        oz = node["oz"].as<Real>( );
+
+        dx = node["dx"].as< VectorXr >();
+        dy = node["dy"].as< VectorXr >();
+        dz = node["dz"].as< VectorXr >();
+
+    }
+
+    RectilinearGrid::~RectilinearGrid() {
+
+    }
+
+    std::shared_ptr< RectilinearGrid > RectilinearGrid::NewSP() {
+        std::shared_ptr<RectilinearGrid> sp(new  RectilinearGrid( ), LemmaObjectDeleter() );
+        return sp;
+    }
+
+    YAML::Node RectilinearGrid::Serialize() const {
+        YAML::Node node = Grid::Serialize();
+
+        node["nx"] = nx;
+        node["ny"] = ny;
+        node["nz"] = nz;
+
+        node["ox"] = ox;
+        node["oy"] = oy;
+        node["oz"] = oz;
+
+        node["dx"] = dx;
+        node["dy"] = dy;
+        node["dz"] = dz;
+
+        node.SetTag( this->GetName() );
+        return node;
+    }
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  DeSerialize
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<RectilinearGrid> RectilinearGrid::DeSerialize ( const YAML::Node& node  ) {
+        if (node.Tag() != "RectilinearGrid") {
+            throw  DeSerializeTypeMismatch( "RectilinearGrid", node.Tag());
+        }
+        std::shared_ptr<RectilinearGrid> Object(new  RectilinearGrid(node), LemmaObjectDeleter() );
+        return Object ;
+    }
+
+
+    // ====================  OPERATIONS    =======================
+
+    void RectilinearGrid::SetDimensions (const int &inx, const int &iny,
+            const int &inz) {
+        dx.resize(inx);
+        dy.resize(iny);
+        dz.resize(inz);
+        nx = inx;
+        ny = iny;
+        nz = inz;
+    }
+
+
+    void RectilinearGrid::SetOffset (const Real& iox, const Real& ioy, const Real& ioz) {
+        ox = iox;
+        oy = ioy;
+        oz = ioz;
+    }
+
+    void RectilinearGrid::SetSpacing (const VectorXr &idx, const VectorXr &idy,
+                const VectorXr &idz) {
+        nx = idx.size();
+        ny = idy.size();
+        nz = idz.size();
+        dx = idx;
+        dy = idy;
+        dz = idz;
+    }
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetNx
+    //--------------------------------------------------------------------------------------
+    int RectilinearGrid::GetNx (  ) {
+        return nx;
+    }		// -----  end of method RectilinearGrid::GetNx  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetNy
+    //--------------------------------------------------------------------------------------
+    int RectilinearGrid::GetNy (  ) {
+        return ny;
+    }		// -----  end of method RectilinearGrid::GetNy  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetNz
+    //--------------------------------------------------------------------------------------
+    int RectilinearGrid::GetNz (  ) {
+        return nz;
+    }		// -----  end of method RectilinearGrid::GetNz  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetOx
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetOx (  ) {
+        return ox;
+    }		// -----  end of method RectilinearGrid::GetOx  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetOy
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetOy (  ) {
+        return oy;
+    }		// -----  end of method RectilinearGrid::GetOy  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetOz
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetOz (  ) {
+        return oz;
+    }		// -----  end of method RectilinearGrid::GetOz  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDx
+    //--------------------------------------------------------------------------------------
+    VectorXr RectilinearGrid::GetDx (  ) {
+        return dx;
+    }		// -----  end of method RectilinearGrid::GetDx  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDx
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetDx ( const int& ix  ) {
+        return dx[ix];
+    }		// -----  end of method RectilinearGrid::GetDx  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDy
+    //--------------------------------------------------------------------------------------
+    VectorXr RectilinearGrid::GetDy (  ) {
+        return dy;
+    }		// -----  end of method RectilinearGrid::GetDy  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDy
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetDy ( const int& iy  ) {
+        return dy[iy];
+    }		// -----  end of method RectilinearGrid::GetDy  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDz
+    //--------------------------------------------------------------------------------------
+    VectorXr RectilinearGrid::GetDz (  ) {
+        return dz;
+    }		// -----  end of method RectilinearGrid::GetDz  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGrid
+    //      Method:  GetDz
+    //--------------------------------------------------------------------------------------
+    Real RectilinearGrid::GetDz ( const int& iz  ) {
+        return dz[iz];
+    }		// -----  end of method RectilinearGrid::GetDz  -----
+
+}		// -----  end of Lemma  name  -----
+
+/* vim: set tabstop=4 expandtab: */
+/* vim: set filetype=cpp: */

Modules/LemmaCore/src/RectilinearGridReader.cpp

@@ -0,0 +1,113 @@
+/* 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
+  @date      09/23/2013 11:05:48 AM
+  @version   $Id$
+  @author    Trevor Irons (ti)
+  @email     Trevor.Irons@xri-geo.com
+  @copyright Copyright (c) 2013, XRI Geophysics, LLC
+  @copyright Copyright (c) 2013, Trevor Irons
+ **/
+
+#include "RectilinearGridReader.h"
+
+namespace Lemma {
+
+    // ====================  FRIEND METHODS  =====================
+
+    std::ostream &operator << (std::ostream &stream, const RectilinearGridReader &ob) {
+        stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+        return stream;
+    }
+
+    // ====================  LIFECYCLE     =======================
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGridReader
+    //      Method:  RectilinearGridReader
+    // Description:  constructor (protected)
+    //--------------------------------------------------------------------------------------
+    RectilinearGridReader::RectilinearGridReader ( ) : GridReader( ),
+        rGrid( nullptr ), Parser( nullptr ) {
+
+    }  // -----  end of method RectilinearGridReader::RectilinearGridReader  (constructor)  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGridReader
+    //      Method:  NewSP()
+    // Description:  constructor
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr< RectilinearGridReader > RectilinearGridReader::NewSP() {
+        std::shared_ptr<RectilinearGridReader> sp(new  RectilinearGridReader( ), LemmaObjectDeleter() );
+        return sp;
+    }
+
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGridReader
+    //      Method:  ~RectilinearGridReader
+    // Description:  destructor (protected)
+    //--------------------------------------------------------------------------------------
+    RectilinearGridReader::~RectilinearGridReader () {
+    }  // -----  end of method RectilinearGridReader::~RectilinearGridReader  (destructor)  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGridReader
+    //      Method:  ReadASCIIGridFile
+    //--------------------------------------------------------------------------------------
+    void RectilinearGridReader::ReadASCIIGridFile ( const std::string& name  ) {
+
+        rGrid = RectilinearGrid::NewSP();
+        Parser = ASCIIParser::NewSP();
+
+        Parser->SetCommentString("//");
+
+        Parser->Open(name);
+
+        std::vector<int> vals = Parser->ReadInts(3);
+        rGrid->SetDimensions(vals[0], vals[1], vals[2]);
+
+        std::vector<Real> rvals = Parser->ReadReals(3);
+        rGrid->SetOffset(rvals[0], rvals[1], rvals[2]);
+        rvals.clear();
+
+        rvals = Parser->ReadReals( vals[0] );
+        VectorXr temp = VectorXr::Map(&rvals[0], vals[0]);
+        VectorXr hx = temp;
+
+        rvals.clear();
+        rvals = Parser->ReadReals( vals[1] );
+        VectorXr hy = VectorXr::Map(&rvals[0], vals[1]);
+
+        rvals.clear();
+        rvals = Parser->ReadReals( vals[2] );
+        VectorXr hz = VectorXr::Map(&rvals[0], vals[2]);
+
+        rGrid->SetSpacing(hx, hy, hz);
+
+        // Read in model(s)/data? Where should this be done?
+        Parser->Close();
+
+        return;
+
+    }		// -----  end of method RectilinearGridReader::ReadSCIIGridFile  -----
+
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  RectilinearGridReader
+    //      Method:  GetGrid
+    //--------------------------------------------------------------------------------------
+    std::shared_ptr<Grid> RectilinearGridReader::GetGrid (  ) {
+        return std::static_pointer_cast<Grid> (rGrid);
+    }		// -----  end of method RectilinearGridReader::GetGrid  -----
+
+
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/RectilinearGridVTKExporter.cpp

@@ -0,0 +1,147 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      09/25/2013 08:20:14 AM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2013, Trevor Irons
+ */
+
+#if  LEMMAUSEVTK
+
+#include "RectilinearGridVTKExporter.h"
+
+namespace Lemma {
+
+// ====================  FRIEND METHODS  =====================
+
+std::ostream &operator << (std::ostream &stream, const RectilinearGridVTKExporter &ob) {
+    stream << ob.Serialize()  << "\n---\n"; // End of doc --- as a direct stream should encapulste thingy
+    return stream;
+}
+
+
+// ====================  LIFECYCLE     =======================
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  RectilinearGridVTKExporter
+// Description:  constructor (protected)
+//--------------------------------------------------------------------------------------
+RectilinearGridVTKExporter::RectilinearGridVTKExporter ( ) : LemmaObject( ), Grid(nullptr), VTKGrid(nullptr) {
+
+}  // -----  end of method RectilinearGridVTKExporter::RectilinearGridVTKExporter  (constructor)  -----
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  NewSP()
+// Description:  public constructor
+//--------------------------------------------------------------------------------------
+std::shared_ptr< RectilinearGridVTKExporter > RectilinearGridVTKExporter::NewSP() {
+    std::shared_ptr<RectilinearGridVTKExporter> sp(new  RectilinearGridVTKExporter( ), LemmaObjectDeleter() );
+    return sp;
+}
+
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  ~RectilinearGridVTKExporter
+// Description:  destructor (protected)
+//--------------------------------------------------------------------------------------
+RectilinearGridVTKExporter::~RectilinearGridVTKExporter () {
+}  // -----  end of method RectilinearGridVTKExporter::~RectilinearGridVTKExporter  (destructor)  -----
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  SetGrid
+//--------------------------------------------------------------------------------------
+void RectilinearGridVTKExporter::SetGrid ( std::shared_ptr<RectilinearGrid> GridIn ) {
+    Grid = GridIn;
+    // I see no harm in just doing this now.
+    BuildVTKRectilinearGrid();
+    return ;
+}		// -----  end of method RectilinearGridVTKExporter::SetGrid  -----
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  GetVTKRectilinearGrid
+//--------------------------------------------------------------------------------------
+vtkSmartPointer<vtkRectilinearGrid> RectilinearGridVTKExporter::GetVTKGrid (  ) {
+    return VTKGrid;
+}		// -----  end of method RectilinearGridVTKExporter::GetVTKRectilinearGrid  -----
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  WriteVTKGrid
+//--------------------------------------------------------------------------------------
+void RectilinearGridVTKExporter::WriteVTKGrid ( const std::string& fname  ) {
+    vtkXMLRectilinearGridWriter *gridWrite = vtkXMLRectilinearGridWriter::New();
+        gridWrite->SetInputData(VTKGrid);
+        gridWrite->SetFileName( (fname+std::string(".vtr")).c_str() );
+        gridWrite->Update();
+        gridWrite->Write();
+        gridWrite->Delete();
+    return ;
+}		// -----  end of method RectilinearGridVTKExporter::WriteVTKGrid  -----
+
+//--------------------------------------------------------------------------------------
+//       Class:  RectilinearGridVTKExporter
+//      Method:  BuildVTKRectilinearGrid
+//--------------------------------------------------------------------------------------
+void RectilinearGridVTKExporter::BuildVTKRectilinearGrid (  ) {
+
+    // Set Coordinate>s
+    vtkDoubleArray *xCoords = vtkDoubleArray::New();
+    xCoords->InsertNextValue(Grid->GetOx()-Grid->GetDx(0)/2.);
+    double xm1 = Grid->GetOx() - Grid->GetDx(0)/2.;
+    for (int ix=0; ix<Grid->GetNx(); ix++) {
+        xCoords->InsertNextValue(xm1 + Grid->GetDx(ix));
+        xm1 += Grid->GetDx(ix);
+    }
+
+    vtkDoubleArray *yCoords = vtkDoubleArray::New();
+    yCoords->InsertNextValue(Grid->GetOy()-Grid->GetDy(0)/2.);
+    double ym1 = Grid->GetOy()-Grid->GetDy(0)/2.;
+    for (int iy=0; iy<Grid->GetNy(); iy++) {
+        yCoords->InsertNextValue(ym1 + Grid->GetDy(iy));
+        ym1 += Grid->GetDy(iy);
+    }
+
+    vtkDoubleArray *zCoords = vtkDoubleArray::New();
+    zCoords->InsertNextValue(Grid->GetOz()-Grid->GetDz(0)/2.);
+    double zm1 = Grid->GetOz()-Grid->GetDz(0)/2.;
+    for (int iz=0; iz<Grid->GetNz(); iz++) {
+        zCoords->InsertNextValue(zm1 + Grid->GetDz(iz));
+        zm1 += Grid->GetDz(iz);
+    }
+
+    VTKGrid = vtkSmartPointer<vtkRectilinearGrid>::New();
+        VTKGrid->SetDimensions(Grid->GetNx()+1, Grid->GetNy()+1, Grid->GetNz()+1);
+        VTKGrid->SetXCoordinates(xCoords);
+        VTKGrid->SetYCoordinates(yCoords);
+        VTKGrid->SetZCoordinates(zCoords);
+
+    // Clean up
+    xCoords->Delete();
+    yCoords->Delete();
+    zCoords->Delete();
+
+    return ;
+}		// -----  end of method RectilinearGridVTKExporter::BuildVTKRectilinearGrid  -----
+
+
+
+}		// -----  end of Lemma  name  -----
+
+#endif     // -----  not LEMMAUSEVTK  -----

Modules/LemmaCore/src/WindowFilter.cpp

@@ -0,0 +1,83 @@
+/* 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     07/20/2010
+  @version  $Id: windowfilter.cpp 87 2013-09-05 22:44:05Z tirons $
+ **/
+
+#include "WindowFilter.h"
+
+namespace Lemma {
+
+    // ====================  LIFECYCLE     =======================
+
+    WindowFilter::WindowFilter( ) : Filter( ),
+        Width(0), SamplingRate(0), Bandwidth(0), Coefficients(),
+        Type(RECTANGULAR) {}
+
+    std::shared_ptr< WindowFilter > WindowFilter::NewSP() {
+        std::shared_ptr<WindowFilter> sp(new  WindowFilter( ), LemmaObjectDeleter() );
+        return sp;
+    }
+
+    WindowFilter::~WindowFilter() {
+    }
+
+
+    // ====================  OPERATIONS    =======================
+
+    void WindowFilter::SetWindowType(const WINDOWTYPE &type) {
+        Type = type;
+    }
+
+    void WindowFilter::SetBandwidth(const Real& width) {
+        Bandwidth = width;
+    }
+
+    void WindowFilter::SetSamplingInterval(const Real& rate) {
+        SamplingRate = rate;
+    }
+
+    int WindowFilter::GetNumberOfFrequencyBins() {
+        return this->Nw;
+    }
+
+    void WindowFilter::SetNumberOfSamples(const int& nt) {
+        Nt = nt;
+        Nw = Nt%2 ? (Nt-1)/2+1 : (Nt)/2+1;  // odd , even
+    }
+
+    VectorXr WindowFilter::GetFilterCoefficients() {
+        Real  Nyquist = .5/SamplingRate;
+        std::cout << "Window nyquist " << Nyquist << std::endl;
+        std::cout << "Window Nw  " << Nw << std::endl;
+        Real  df      = (Nyquist)/((Real)(Nw)); // df Hz per bin
+        Width         = Bandwidth/df;
+        Coefficients.resize(Width);
+        switch (Type) {
+            case HAMMING:
+           	    for (int n=0; n<Width; ++n) {
+		            Coefficients(n) = 0.5 * (1.- std::cos( (2.*PI*(Real)(n)) /
+                                       ((Real)(Width)-1.) ) );
+	            }
+                break;
+            case HANNING:
+	            for (int n=0; n<Width; ++n) {
+		            Coefficients(n) = 0.54 - 0.46 *
+                                std::cos( (2.*PI*(Real)(n)) /
+                                    ((Real)(Width)-1.) ) ;
+	            }
+                break;
+            case RECTANGULAR:
+                Coefficients.setOnes();
+                break;
+        }
+        return Coefficients;
+    }
+}		// -----  end of Lemma  name  -----

Modules/LemmaCore/src/banner.cpp

@@ -0,0 +1,36 @@
+/* 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   M. Andy Kass
+  @date     05/03/2010
+  @version  $Id: banner.cpp 87 2013-09-05 22:44:05Z tirons $
+ **/
+
+#include "banner.h"
+
+///Prints banner information to the screen
+void banner(const std::string &name, const std::string &version,
+		const std::string &usage) {
+
+	std::cout << std::endl;
+	std::cout << "  " << name.c_str() << version.c_str() << std::endl;
+	std::cout << std::endl;
+	std::cout << "  A program built from Lemma v 1.0b" << std::endl;
+	std::cout << "  Usage: " << usage.c_str() << std::endl;
+	std::cout << std::endl;
+	std::cout << "  Developed jointly by the Center for Gravity, " <<
+		"Electrical, and Magnetic" <<std::endl;
+	std::cout << "   Studies (CGEM), the United States Geologic " <<
+		"Survey (USGS), and Broken" << std::endl;
+	std::cout << "   Spoke Development, LLC." <<std::endl;
+	std::cout << std::endl;
+	std::cout << "  (C) Copyright 2010, Trevor Irons and M. Andy Kass.  " <<
+		"Some rights reserved." << std::endl;
+	std::cout << "  Please see license for details." << std::endl;
+	std::cout << std::endl;
+}

Modules/LemmaCore/src/helper.cpp

@@ -0,0 +1,201 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      10/02/2014 03:21:07 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include "helper.h"
+
+namespace Lemma {
+
+std::string enum2String(const FREQUENCYUNITS& FreqUnits) {
+    std::string t;
+    switch (FreqUnits) {
+        case HZ:
+            t = std::string("HZ");
+            break;
+        case KHZ:
+            t = std::string("KHZ");
+            break;
+        case MHZ:
+            t = std::string("MHZ");
+            break;
+        case GHZ:
+            t = std::string("GHZ");
+            break;
+    }
+    return t;
+}
+
+std::string enum2String(const TIMEUNITS& Units) {
+    std::string t;
+    switch (Units) {
+        case SEC:
+            t = std::string("SEC");
+            break;
+        case MILLISEC:
+            t = std::string("MILLISEC");
+            break;
+        case MICROSEC:
+            t = std::string("MICROSEC");
+            break;
+        case NANOSEC:
+            t = std::string("NANOSEC");
+            break;
+        case PICOSEC:
+            t = std::string("PICOSEC");
+            break;
+    }
+    return t;
+}
+
+
+std::string enum2String(const MAGUNITS& Units) {
+    std::string t;
+    switch (Units) {
+        case TESLA:
+            t = std::string("TESLA");
+            break;
+        case NANOTESLA:
+            t = std::string("NANOTESLA");
+            break;
+        case GAUSS:
+            t = std::string("GAUSS");
+            break;
+    }
+    return t;
+}
+
+std::string enum2String(const FIELDCALCULATIONS& Field) {
+    std::string t;
+    switch (Field) {
+        case E:
+            t = std::string("E");
+            break;
+        case H:
+            t = std::string("H");
+            break;
+        case BOTH:
+            t = std::string("BOTH");
+            break;
+    }
+    return t;
+}
+
+
+std::string enum2String(const TEMPUNITS& Units) {
+    std::string t;
+    switch (Units) {
+        case CELCIUS:
+            t = std::string("CELCIUS");
+            break;
+        case KELVIN:
+            t = std::string("KELVIN");
+            break;
+    }
+    return t;
+}
+
+std::string enum2String(const FIELDCOMPONENT& Comp) {
+    std::string t;
+    switch (Comp) {
+        case XCOMPONENT:
+            t = std::string("XCOMPONENT");
+            break;
+        case YCOMPONENT:
+            t = std::string("YCOMPONENT");
+            break;
+        case ZCOMPONENT:
+            t = std::string("ZCOMPONENT");
+            break;
+    }
+    return t;
+}
+
+std::string enum2String(const HANKELTRANSFORMTYPE& Type) {
+    std::string t;
+    switch (Type) {
+        case ANDERSON801:
+            t = std::string("ANDERSON801");
+            break;
+        case CHAVE:
+            t = std::string("CHAVE");
+            break;
+        case QWEKEY:
+            t = std::string("QWEKEY");
+            break;
+        case FHTKEY201:
+            t = std::string("FHTKEY201");
+            break;
+        case FHTKEY101:
+            t = std::string("FHTKEY101");
+            break;
+        case FHTKEY51:
+            t = std::string("FHTKEY51");
+            break;
+    }
+    return t;
+}
+
+template<>
+FREQUENCYUNITS string2Enum<FREQUENCYUNITS>( const std::string& str ) {
+    if       (str ==  "HZ") return   HZ;
+    else if  (str == "KHZ") return  KHZ;
+    else if  (str == "MHZ") return  MHZ;
+    else if  (str == "GHZ") return  GHZ;
+    else {
+        throw std::runtime_error("string not recognized as FREUENCYUNIT");
+    }
+}
+
+template<>
+HANKELTRANSFORMTYPE string2Enum<HANKELTRANSFORMTYPE>( const std::string& str ) {
+    if       (str ==  "ANDERSON801") return   ANDERSON801;
+    else if  (str == "CHAVE") return  CHAVE;
+    else if  (str == "QWEKEY") return  QWEKEY;
+    else if  (str == "FHTKEY201") return  FHTKEY201;
+    else if  (str == "FHTKEY51") return  FHTKEY51;
+    else if  (str == "FHTKEY101") return  FHTKEY101;
+    else {
+        throw std::runtime_error("string not recognized as HANKELTRANSFORMTYPE");
+    }
+}
+
+
+template<>
+TIMEUNITS string2Enum<TIMEUNITS>( const std::string& str ) {
+    if       (str == "SEC")      return      SEC;
+    else if  (str == "MILLISEC") return MILLISEC;
+    else if  (str == "MICROSEC") return MICROSEC;
+    else if  (str == "NANOSEC")  return  NANOSEC;
+    else if  (str == "PICOSEC")  return  PICOSEC;
+    else {
+        throw std::runtime_error("string not recognized as TIMEUNIT");
+    }
+}
+
+template<>
+FIELDCOMPONENT string2Enum<FIELDCOMPONENT>( const std::string& str) {
+    if       (str == "XCOMPONENT")      return XCOMPONENT;
+    else if  (str == "YCOMPONENT")      return YCOMPONENT;
+    else if  (str == "ZCOMPONENT")      return ZCOMPONENT;
+    else {
+        throw std::runtime_error("string not recognized as FieldComponent");
+    }
+}
+
+}		// -----  end of Lemma  name  -----
+

Modules/LemmaCore/testing/CMakeLists.txt

@@ -0,0 +1,12 @@
+#add_executable(testOne testOne.cpp)
+#add_executable(testTwo testTwo.cpp)
+#add_test( testLemmaCoreOne testOne )
+#add_test( testLemmaCoreTwo testTwo )
+    	
+include_directories(${CXXTEST_INCLUDE_DIR})
+CXXTEST_ADD_TEST(unittest_GetNameCheck GetNameCheck.cc ${CMAKE_CURRENT_SOURCE_DIR}/GetNameCheck.h)
+target_link_libraries(unittest_GetNameCheck "lemmacore")
+
+CXXTEST_ADD_TEST(unittest_SerializeCheck SerializeCheck.cc ${CMAKE_CURRENT_SOURCE_DIR}/SerializeCheck.h)
+target_link_libraries(unittest_SerializeCheck "lemmacore")
+ 

Modules/LemmaCore/testing/GetNameCheck.h

@@ -0,0 +1,69 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      06/16/2016 09:12:46 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Trevor Irons & Lemma Software, LLC
+ */
+
+#include <cxxtest/TestSuite.h>
+#include <LemmaCore>
+
+using namespace Lemma;
+
+class MyTestSuite : public CxxTest::TestSuite
+{
+    public:
+
+    void testASCIIParser( void )
+    {
+        auto Obj = ASCIIParser::NewSP();
+        TS_ASSERT_EQUALS( Obj->GetName(), std::string("ASCIIParser") );
+    }
+
+    void testCubicSplineInterpolator(void)
+    {
+        auto Obj = CubicSplineInterpolator::NewSP();
+        TS_ASSERT_EQUALS( Obj->GetName(), std::string("CubicSplineInterpolator") );
+    }
+
+    void testRectilinearGrid( void )
+    {
+        auto Obj = RectilinearGrid::NewSP();
+        TS_ASSERT_EQUALS( Obj->GetName(), std::string("RectilinearGrid") );
+    }
+
+    void testRectilinearGridReader( void )
+    {
+        auto Obj = RectilinearGridReader::NewSP();
+        TS_ASSERT_EQUALS( Obj->GetName(), std::string("RectilinearGridReader") );
+    }
+
+// INCLUDE GUARDS DON'T WORK WITH CXXTEST, NEED TO MOVE VTK TESTS OUT
+// #ifdef LEMMAUSEVTK
+//     void testRectilinearGridVTKExporter( void )
+//     {
+//         auto Obj = RectilinearGridVTKExporter::NewSP();
+//         TS_ASSERT_EQUALS( Obj->GetName(), std::string("RectilinearGridVTKExporter") );
+//     }
+// #endif
+
+    void testWindowFilter( void )
+    {
+        auto Obj = WindowFilter::NewSP();
+        TS_ASSERT_EQUALS( Obj->GetName(), std::string("WindowFilter") );
+    }
+
+};
+

Modules/LemmaCore/testing/SerializeCheck.h

@@ -0,0 +1,75 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      06/16/2016 09:12:46 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Trevor Irons & Lemma Software, LLC
+ */
+
+#include <cxxtest/TestSuite.h>
+#include <LemmaCore>
+
+using namespace Lemma;
+
+class MyTestSuite : public CxxTest::TestSuite
+{
+    public:
+
+    void testASCIIParser( void )
+    {
+        auto Obj = ASCIIParser::NewSP();
+        YAML::Node node = Obj->Serialize();
+        auto Obj2 = ASCIIParser::DeSerialize(node);
+        TS_ASSERT_EQUALS( Obj->GetName(), Obj2->GetName() );
+    }
+
+
+    void testCubicSplineInterpolator(void)
+    {
+        auto Obj = CubicSplineInterpolator::NewSP();
+        YAML::Node node = Obj->Serialize();
+        auto Obj2 = CubicSplineInterpolator::DeSerialize(node);
+        TS_ASSERT_EQUALS( Obj->GetName(), Obj2->GetName() );
+    }
+
+    void testRectilinearGrid( void )
+    {
+        auto Obj = RectilinearGrid::NewSP();
+        YAML::Node node = Obj->Serialize();
+        auto Obj2 = RectilinearGrid::DeSerialize(node);
+        TS_ASSERT_EQUALS( Obj->GetName(), Obj2->GetName() );
+    }
+
+// /*
+//     void testRectilinearGridReader( void )
+//     {
+//         auto Obj = RectilinearGridReader::NewSP();
+//         TS_ASSERT_EQUALS( Obj->GetName(), std::string("RectilinearGridReader") );
+//     }
+//
+//     void testRectilinearGridVTKExporter( void )
+//     {
+//         auto Obj = RectilinearGridVTKExporter::NewSP();
+//         TS_ASSERT_EQUALS( Obj->GetName(), std::string("RectilinearGridVTKExporter") );
+//     }
+//
+//     void testWindowFilter( void )
+//     {
+//         auto Obj = WindowFilter::NewSP();
+//         TS_ASSERT_EQUALS( Obj->GetName(), std::string("WindowFilter") );
+//     }
+// */
+
+};
+

Modules/LemmaCore/testing/testOne.cpp

@@ -0,0 +1,23 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      06/14/2016 09:30:38 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Trevor Irons & Lemma Software, LLC
+ */
+
+int main() {
+    return 0;
+}
+

Modules/LemmaCore/testing/testTwo.cpp

@@ -0,0 +1,23 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      06/14/2016 09:43:43 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Trevor Irons & Lemma Software, LLC
+ */
+
+int main() {
+    return 0;
+}
+

config/LemmaConfig.h.in

@@ -0,0 +1,25 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* 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
+ * @date      01/15/2016 09:52:10 PM
+ * @author    Trevor Irons (ti)
+ * @email     tirons@egi.utah.edu
+ * @copyright Copyright (c) 2016, University of Utah
+ * @copyright Copyright (c) 2016, Trevor Irons & Lemma Software, LLC
+ */
+
+// the configured options and settings for Tutorial
+#define LEMMA_VERSION_MAJOR @LEMMA_VERSION_MAJOR@
+#define LEMMA_VERSION_MINOR @LEMMA_VERSION_MINOR@
+#define LEMMA_VERSION_PATCH @LEMMA_VERSION_PATCH@
+//#define LEMMA_VERSION @LEMMA_VERSION@
+
+