Initial commit. Updated material and funcs

2016-11-11 23:22:44 +01:00 · 2016-11-11 23:22:44 +01:00 · 83890dda85
commit 83890dda85
36 changed files with 1858 additions and 0 deletions
--- a/.gitignore
+++ b/.gitignore
@ -0,0 +1,14 @@
+CMakeCache.txt
+CMakeFiles
+Makefile
+*.so
+nonlinPYTHON_wrap.cxx
+*.a
+*.dll
+*.pyd
+__pycache__
+cmake_install.cmake
+doc/html
+doc/latex
+doc/usg.pdf
+*.pyc
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -0,0 +1,139 @@
+# CMakeList.txt for TaSMET
+cmake_minimum_required (VERSION 2.8)
+project(TaSMET)
+set(PACKAGE_VERSION 0.1)
+message("Running Cmake for TaSMET version ${PACKAGE_VERSION}")
+
+# Set the Python version to link against. If this is not set, Cmake tries to find it automatically.
+# set(TaSMET_PY_VERSION "2.7")
+set(TaSMET_PY_VERSION "3.5m")
+
+# Tracer name (name of the variable)
+add_definitions(-DTRACER=1)
+# add_definitions(-DTRACER_IN_COMMON)
+
+add_definitions(-DTASMET_FLOAT=64)
+
+#====================================================
+# Compiler settings ************************************************
+#====================================================
+
+# Always required make flags in case of both compilers
+set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pipe -fPIC -Wall -Wextra -Wno-unused-parameter -Wno-unused-variable  ")
+
+
+# Stop letting Numpy complain about its API
+add_definitions(-DNPY_NO_DEPRECATED_API=NPY_1_4_API_VERSION)
+
+#==================================================
+# Optimized code flags *******************************************
+#==================================================
+
+#Debug mode
+# set(CMAKE_CXX_FLAGS  "${CMAKE_CXX)FLAGS} -g -ggdb")
+
+# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -march=native -mtune=native")
+# set(CMAKE_CLANG "${CMAKE_GCC} -march=native -mtune=native -fopenmp")
+
+# To disable bound checking on std::vector, and to disable assertions
+# add_definitions(-DNDEBUG)
+
+# To increase speed on Armadillo
+# add_definitions(-DARMA_NO_DEBUG)
+
+# Disable traces
+# add_definitions(-DTRACER=0)
+
+# Pre-allocation size for matrices. Very important setting to tune the
+# code in case we know that we are going to run with certain sizes of matrices and vectors. For Nf=6, set this to (2*Nf+1)^2=169
+# For Nf=1: 9
+# For Nf=2: 25
+# For Nf=3: 49
+# For Nf=4: 81
+# For Nf=5: 121
+# For Nf=6: 169
+# For Nf=7: 225
+# For Nf=8: 289
+# For Nf=9: 361
+# For Nf=10: 441
+# For Nf=11: 529
+# For Nf=12: 625
+# Watch out! Setting prealloc too high can give too much overhead for smaller Nf's
+# add_definitions(-DARMA_MAT_PREALLOC=625)
+
+
+# ##########################
+# Finding the presence of the prerequisites
+# ##########################
+
+
+# For importing find directives for Cmake
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${PROJECT_SOURCE_DIR}/common/cmake_tools)
+
+# ##########################
+# Python #####################
+# ##########################
+
+if(TaSMET_PY_VERSION)
+  # Find major version from version string
+  set(PYTHON_LIBRARY "/usr/lib/libpython${TaSMET_PY_VERSION}.so")
+  set(PYTHON_INCLUDE_DIR "/usr/include/python${TaSMET_PY_VERSION}")  
+  set(PYTHON_INCLUDE_DIRS "/usr/include/python${TaSMET_PY_VERSION}")  
+
+endif(TaSMET_PY_VERSION)
+message("Python include dirs: ${PYTHON_INCLUDE_DIRS}")
+
+find_package(PythonLibs REQUIRED)
+string(REGEX MATCH "^." TaSMET_PY_MAJOR_VERSION ${PYTHONLIBS_VERSION_STRING})
+MESSAGE("Python major version: ${TaSMET_PY_MAJOR_VERSION}")  
+
+# Find the site_packages directory of python
+execute_process(COMMAND python${TaSMET_PY_MAJOR_VERSION} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib())" OUTPUT_VARIABLE PYTHON_SITE_PACKAGES OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+# ################################
+# Initialize swig
+# ################################
+
+find_package(SWIG REQUIRED)
+include(${SWIG_USE_FILE})
+
+SET(CMAKE_SWIG_FLAGS  -Wall -DSWIG_PYTHON)
+if(${TaSMET_PY_MAJOR_VERSION}=="3")
+  SET(CMAKE_SWIG_FLAGS  ${CMAKE_SWIG_FLAGS} -py${TaSMET_PY_MAJOR_VERSION})  
+endif(${TaSMET_PY_MAJOR_VERSION}=="3")
+
+include_directories(common/src)
+
+# Armadillo
+find_package(Armadillo REQUIRED)
+add_definitions(-DARMA_USE_SUPERLU -DARMA_USE_CXX11)
+
+# ==================== Compile the code in common and src
+
+# This is the common code (gas and solid libs, etc)
+# link_directories(common)
+
+aux_source_directory(common/src/gas gas)
+include_directories(
+  ${PYTHON_INCLUDE_DIRS}  
+  common/src/swig
+  common/src/gas
+  )
+# Add the code subdirectory
+add_subdirectory(src)
+
+
+# set_source_files_properties( ${swig_generated_file_fullname}
+  # PROPERTIES COMPILE_FLAGS "${SWIG_COMMON_COMPILE_FLAGS} ")
+
+
+# ================================== Installation
+
+# Install common files
+install(FILES ${PROJECT_SOURCE_DIR}/common/__init__.py
+  DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common)
+install(FILES ${PROJECT_SOURCE_DIR}/common/common.py
+  DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME}/common)
+
+# Rest of the files is installed from src/CMakeLists.txt
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@ -0,0 +1,69 @@
+#================
+# The TaSMET code
+#================
+
+include_directories(
+  ${PYTHON_INCLUDE_PATH}
+  ${ARMADILLO_INCLUDE_DIRS}
+  .
+  # duct
+  # duct/cell
+  # duct/connectors
+  # duct/eq
+  # duct/geom
+  # duct/geom/grid
+  # duct/drag
+  # duct/heat
+  # mech
+  # python
+  # seg
+  # sol
+  # sys
+  # var
+  # volume
+  )
+
+add_subdirectory(funcs)
+add_subdirectory(material)
+# add_subdirectory(duct)
+# add_subdirectory(mech)
+# add_subdirectory(seg)
+# add_subdirectory(sol)
+# add_subdirectory(sys)
+# add_subdirectory(var)
+
+add_library(tasmet_src tasmet_tracer.cpp tasmet_exception.cpp)
+target_link_libraries(tasmet_src
+  # duct
+  # mech
+  # seg
+  # sol
+  # sys
+  # var
+  )
+ 
+# set_source_files_properties(swig/nonlin.i
+#   PROPERTIES CPLUSPLUS ON)
+# swig_add_module(TaSMET python swig/nonlin.i)
+
+# set(SWIG_MODULE_TaSMET_FILE_nonlin_EXTRA_DEPS
+#   ${CMAKE_CURRENT_SOURCE_DIR}/common/swig_files/arma_numpy.i
+# )
+# swig_link_libraries(TaSMET tasmet_src common ${PYTHON_LIBRARIES}
+#   ${ARMADILLO_LIBRARIES})
+
+# set_source_files_properties( ${swig_generated_file_fullname}
+  # PROPERTIES COMPILE_FLAGS "${SWIG_COMMON_COMPILE_FLAGS} ")
+  
+# Install swig files to the right place
+# install(TARGETS _TaSMET
+#   DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME})
+# install(FILES ${CMAKE_BINARY_DIR}/src/TaSMET.py
+#   DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME})
+# install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/python/__init__.py
+#   DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME})
+# install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/python/post
+#   DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME})
+# install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/python/gui
+#   DESTINATION ${PYTHON_SITE_PACKAGES}/${PROJECT_NAME})
+
--- a/src/funcs/CMakeLists.txt
+++ b/src/funcs/CMakeLists.txt
@ -0,0 +1,7 @@
+add_library(funcs
+  bessel.cpp
+  cbessj.cpp
+  rottfuncs.cpp
+  skewsine.cpp
+  sph_j.cpp
+  )
--- a/src/funcs/bessel.cpp
+++ b/src/funcs/bessel.cpp
@ -0,0 +1,102 @@
+#include "bessel.h"
+#include "cbessj.h"
+#include "skewsine.h"
+
+namespace {  
+    c besselj0_smallarg(c x){
+        //Abramowich p. 369, polynomial approximation
+        //For complex numbers, unfortunately we have to check
+        //for both imaginary and real part.
+        /* return 1.0		 \
+        //   -2.2499997*pow(x/3.0,2) \
+        //   +1.2656208*pow(x/3.0,4) \
+        //   -0.3163866*pow(x/3.0,6) \
+        //   +0.0444479*pow(x/3.0,8) \
+        //   -0.0039444*pow(x/3.0,10) \
+        //   +0.0002100*pow(x/3.0,12); */
+
+        // External implementation
+        d z[2];
+        d res[2];
+        z[0]=x.real();
+        z[1]=x.imag();
+        CBESSJ(z,0,res);
+        return c(res[0],res[1]);
+    }
+    c besselj0_largearg(c x){
+        // larger than 3" << endl;
+        // Abramovich, p. 370
+        c f0=0.79788456-0.00000077*(3.0/x)-
+            0.00552740*pow(3.0/x,2)-
+            0.00009512*pow(3.0/x,3)+
+            0.00137237*pow(3.0/x,4)-
+            0.00072805*pow(3.0/x,5)+
+            0.00014476*pow(3.0/x,6);
+        c th0=x-0.78539816-
+            0.04166397*(3.0/x)-
+            0.00003954*pow(3.0/x,2)+
+            0.00262573*pow(3.0/x,3)-
+            0.00054125*pow(3.0/x,4)-
+            0.00029333*pow(3.0/x,5)+
+            0.00013558*pow(3.0/x,6);
+        return pow(x,-0.5)*f0*cos(th0);
+    }
+    c besselj1_over_x_smallarg(c x){
+        // return 0.5-0.56249985*pow(x/3.0,2)+0.21093573*pow(x/3.0,4)-0.03954289*pow(x/3.0,6)+0.00443319*pow(x/3.0,8)-0.00031761*pow(x/3.0,10)+0.00001109*pow(x/3.0,12);
+
+        // External implementation
+        d z[2];
+        d res[2];
+        z[0]=x.real();
+        z[1]=x.imag();
+        CBESSJ(z,1,res);
+        return c(res[0],res[1])/x;
+    }
+    c besselj1_over_x_largearg(c x){
+        c f1=
+            0.79788456
+            +0.00000156*(3.0/x)
+            +0.01659667*pow(3.0/x,2)
+            +0.00017105*pow(3.0/x,3)
+            -0.00249511*pow(3.0/x,4)
+            +0.00113653*pow(3.0/x,5)
+            -0.00020033*pow(3.0/x,6);
+
+        c th1=x
+            -2.35619449
+            +0.12499612*(3.0/x)
+            +0.00005650*pow(3.0/x,2)
+            -0.00637879*pow(3.0/x,3)
+            +0.00074348*pow(3.0/x,4)
+            +0.00079824*pow(3.0/x,5)
+            -0.00029166*pow(3.0/x,6);
+        return pow(x,-1.5)*f1*cos(th1);
+
+    }
+
+}
+
+c besselj0(c x){
+    d trstart=11.0;		// Transition start point
+    d trdelta=1.0;		// Transition deltax
+    if (abs(x)<trstart)
+        return besselj0_smallarg(x);
+    else if(abs(x)<(trstart+trdelta))
+        return skewsine((abs(x)-trstart)/trdelta)*besselj0_largearg(x)	\
+            +(1.0-skewsine((abs(x)-trstart)/trdelta))*besselj0_smallarg(x);
+    else
+        return besselj0_largearg(x);
+}
+c besselj1_over_x(c x) {
+    d trstart=8.0;
+    d trdelta=1.0;
+    if (abs(x)<trstart) 
+        return besselj1_over_x_smallarg(x);
+    else if(abs(x)<trstart+trdelta)
+        return skewsine((abs(x)-trstart)/trdelta)*besselj1_over_x_largearg(x)+(1-skewsine((abs(x)-trstart)/trdelta))*besselj1_over_x_smallarg(x);      
+    else {
+        return besselj1_over_x_largearg(x);
+    }
+
+}
+
--- a/src/funcs/bessel.h
+++ b/src/funcs/bessel.h
@ -0,0 +1,23 @@
+// bessel.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef BESSEL_H
+#define BESSEL_H
+#include "tasmet_types.h"
+
+// Bessel function of first kind and zero'th order for complex numbers
+// j_0(x)
+c besselj0(c x);
+
+// Bessel function of first kind and first order for complex numbers:
+// j1(x)/x
+c besselj1_over_x(c x);
+
+
+#endif // BESSEL_H
+//////////////////////////////////////////////////////////////////////
--- a/src/funcs/cbessj.cpp
+++ b/src/funcs/cbessj.cpp
@ -0,0 +1,187 @@
+// Copied from Jean-Pierre Moreau's Home Page
+// I would like to thank him for sharing his code.
+/*****************************************************************
+ *    Complex Bessel Function of the 1st Kind of integer order    *
+ * -------------------------------------------------------------- *
+ * SAMPLE RUN:                                                    *
+ *                                                                *
+ * Complex Bessel Function of the 1st Kind of integer order       *
+ *                                                                * 
+ * Input complex argument (real imaginary): 1 2                   *
+ * Input integer order: 1                                         *
+ *                                                                *
+ * Function value: 1.291848  1.010488                             *
+ *                                                                *
+ *                                                                *
+ *                          C++ Release 1.2 By J-P Moreau, Paris. *
+ *                                   (www.jpmoreau.fr)            *
+ * -------------------------------------------------------------- *
+ * Release 1.1: Corrected bug in Function ZLn (atan replaced by   *
+ *              function ATAN2) 11/10/2005.                       *
+ * Release 1.2: ZPower replaced by IZpower (integer exponant).    *
+ *              Limitations increased in CDIV, MAXK=20            * 
+ *****************************************************************/
+#include <stdio.h>
+#include <math.h>
+
+#define  MAXK  20    //09/21/2009 (15 before)
+#define  FALSE  0
+#define  TRUE   1
+
+double HALF = 0.5, ONE = 1.0, FPF = 5.5;
+double PI=4.0*atan(1.0);
+
+int nu;              // order of complex Bessel
+double z[2], z1[2];  // Complex numbers
+double x,y;
+
+// Z=Z1/Z2
+void CDIV(double *Z1, double *Z2, double *Z) {
+  double D;
+  D=Z2[0]*Z2[0]+Z2[1]*Z2[1];
+  if (D > 1e30) {
+    Z[0]=0.0; Z[1]=0.0; return;
+  }
+  if (D == 0) return;
+  Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D;
+  Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D;
+}
+
+// Z=Z1*Z2
+void CMUL(double *Z1, double *Z2, double *Z) {
+  Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1];
+  Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0];
+}
+
+// compute Z^N
+void IZPower(double *z, int n, double *z1) {
+  double temp[2],temp1[2]; int i;
+  if (n==0) {
+    z1[0]=1.0;
+    z1[1]=0.0;
+  }
+  else if (n==1) {
+    z1[0]=z[0];
+    z1[1]=z[1];
+  }
+  else {
+    temp1[0]=z[0]; temp1[1]=z[1];
+    for (i=2; i<=n; i++) {
+      CMUL(temp1,z,temp);
+      temp1[0]=temp[0];
+      temp1[1]=temp[1];
+    }
+    z1[0]=temp[0];
+    z1[1]=temp[1];
+  }
+}
+
+double Fact(int k) {
+  int i; double f;
+  f=ONE;
+  for (i=2; i<=k; i++)  f *= ONE*i;
+  return f;
+}
+
+/******************************************
+ *           FUNCTION  GAMMA(X)            *
+ * --------------------------------------- *
+ * Returns the value of Gamma(x) in double *
+ * precision as EXP(LN(GAMMA(X))) for X>0. *
+ ******************************************/
+double Gamma(double xx) {
+  double cof[6];
+  double stp,x,tmp,ser;
+  int j;
+  cof[0]=76.18009173;
+  cof[1]=-86.50532033;
+  cof[2]=24.01409822;
+  cof[3]=-1.231739516;
+  cof[4]=0.120858003e-2;
+  cof[5]=-0.536382e-5;
+  stp=2.50662827465;
+  x=xx-ONE;
+  tmp=x+FPF;
+  tmp=(x+HALF)*log(tmp)-tmp;
+  ser=ONE;
+  for (j=0; j<6; j++) {
+    x += ONE;
+    ser += cof[j]/x;
+  }
+  return (exp(tmp+log(stp*ser)));
+}
+
+// main subroutine
+void CBESSJ(double *z, int nu, double *z1)  {
+  /*--------------------------------------------------
+    inf.     (-z^2/4)^k
+    Jnu(z) = (z/2)^nu x Sum  ------------------
+    k=0  k! x Gamma(nu+k+1)
+    (nu must be >= 0). Here k=15.
+    ---------------------------------------------------*/
+  int k;
+  double sum[2],tmp[2],tmp1[2],tmp2[2],tmp3[2];
+  // calculate (z/2)^nu in tmp3
+  tmp[0]=2.0; tmp[1]=0.0;
+  CDIV(z,tmp,tmp1);
+  IZPower(tmp1,nu,tmp3);
+  sum[0]=0.0; sum[1]=0.0;
+  //calculate Sum
+  for (k=0; k<=MAXK; k++) {
+    // calculate (-z^2/4)^k
+    IZPower(z,2,tmp);
+    tmp[0]=-tmp[0]; tmp[1]=-tmp[1];
+    tmp1[0]=4.0; tmp1[1]=0.0;
+    CDIV(tmp,tmp1,tmp2);
+
+    if (z[1]==0) {             //case real number
+      tmp[0]=pow(tmp2[0],k);
+      tmp[1]=0.0;
+    }
+    else                       //case complex number
+      IZPower(tmp2,k,tmp);
+
+    // divide by k!
+    tmp1[0]=Fact(k); tmp1[1]=0.0;
+    if (z[1]==0) {
+      tmp2[0]=tmp[0]/tmp1[0];
+      tmp2[1]=0.0;
+    }
+    else
+      CDIV(tmp,tmp1,tmp2);
+    // divide by Gamma(nu+k+1)
+    tmp1[0]=Gamma(ONE*(nu+k+1)); tmp1[1]=0.0;
+
+    if (z[1]==0) {
+      tmp[0]=tmp2[0]/tmp1[0];
+      tmp[1]=0.0;
+    }
+    else
+      CDIV(tmp2,tmp1,tmp);
+    // actualize sum
+    sum[0] += tmp[0];
+    sum[1] += tmp[1];
+
+  }
+  // multiply (z/2)^nu by sum
+  CMUL(tmp3,sum,z1);
+}
+
+
+// int main() {
+
+//   double z1[2];
+
+//   printf("\n Complex Bessel Function of the 1st Kind of integer order\n\n");
+//   printf(" Input complex argument (real imaginary): ");
+//   scanf("%lf %lf", &x, &y);
+//   z[0]=x; z[1]=y;
+//   printf(" Input integer order: "); scanf("%d", &nu);
+
+//   CBESSJ(z,nu,z1);
+
+//   printf("\n Function value: %f  %f\n\n", z1[0], z1[1]);
+
+// }
+
+// end of file cbessj.cpp
--- a/src/funcs/cbessj.h
+++ b/src/funcs/cbessj.h
@ -0,0 +1,22 @@
+// cbessj.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Implementation of the besselJ function for complex numbers
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef CBESSJ_H
+#define CBESSJ_H
+
+// z: input, z1: result;
+
+void CBESSJ(double *z, int nu, double *z1);
+
+
+#endif // CBESSJ_H
+//////////////////////////////////////////////////////////////////////
+
+
+
+
--- a/src/funcs/rottfuncs.cpp
+++ b/src/funcs/rottfuncs.cpp
@ -0,0 +1,126 @@
+// rottfuncs.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+// Wat doet dit bestand?
+//////////////////////////////////////////////////////////////////////
+#include "rottfuncs.h"
+#include "bessel.h"
+
+namespace {
+    template <class array_type>
+    array_type f_inviscid(const array_type& rh_over_delta){
+        TRACE(0,"f_inviscid");
+        return 0.0*rh_over_delta;
+    }
+    
+    template <class array_type>
+    array_type f_vert(const array_type& rh_over_delta){ //Vertical plate f_nu/f_k
+        TRACE(0,"f_vert");
+        c a=I+1.0;
+        return tanh(a*rh_over_delta)/(a*rh_over_delta);
+    }
+
+    template <class array_type>
+    array_type f_blapprox(const array_type& rh_over_delta){
+        TRACE(0,"f_blapprox()");
+        return (1.0-I)/rh_over_delta/2.0;
+    }
+  
+    c f_circ(const c& rh_over_delta) {
+        TRACE(0,"f_circ(c rh_ov_delta)");
+        // vc s=sq2*rh_over_delta;
+        c a(-1,1.0);
+        c R_over_delta=2.0*rh_over_delta;
+        c jvarg=a*R_over_delta;
+        c j0=besselj0(jvarg);
+        c j1_over_jxarg=besselj1_over_x(jvarg);;
+        return 2.0*j1_over_jxarg/j0;
+    }
+    vc f_circ(const vc& rh_over_delta) {
+        TRACE(0,"f_circ()");
+        return vc(rh_over_delta).transform([] (c x) {return f_circ(x);});
+    }
+
+    c f_square(const c& rh_over_delta){
+        TRACE(0,"f_square()");
+        // Adjusted sligthly for improved speed.
+        // Needs checking if rh is defined right.
+        // For square pore with dimensions AxA : rh=S/PI = A^2 / 4A = A/4
+        c F(0,0);
+        c C(0,0);
+        us n,m,msq,nsq;
+        d pisq=pow(number_pi,2);
+        for(n=1; n<11; n=n+2)
+            {
+                for(m = 1; m<11; m=m+2)
+                    {
+                        msq=m*m;
+                        nsq=n*n;
+                        C = 1.0-I*pisq/(32.0*pow(rh_over_delta,2))*((double)(msq+nsq));
+                        F+=(1.0/msq/nsq)/C;
+                    }
+            }
+        return 1.0-64.0/(pisq*pisq)*F;
+    }
+    // To transform this one would introduce an unallowable loop overhead
+    vc f_square(const vc& rh_over_delta){
+        TRACE(0,"f_square()");
+        // Adjusted sligthly for improved speed.
+        // Needs checking if rh is defined right.
+        // For square pore with dimensions AxA : rh=S/PI = A^2 / 4A = A/4
+        vc F(rh_over_delta.size(),fillwith::zeros);
+        vc C(rh_over_delta.size(),fillwith::zeros);
+        us n,m,msq,nsq;
+        d pisq=pow(number_pi,2);
+        for(n=1; n<11; n=n+2)
+            {
+                for(m = 1; m<11; m=m+2)
+                    {
+                        msq=m*m;
+                        nsq=n*n;
+                        C = 1.0-I*pisq/(32.0*pow(rh_over_delta,2))*(msq+nsq);
+                        F+=(1.0/msq/nsq)/C;
+                    }
+            }
+        return 1.0-64.0/(pisq*pisq)*F;
+    }
+}
+RottFuncs::RottFuncs(const cshape shape):_shape(shape){
+    const string shapestr = cshapeToString(shape);
+    TRACE(8,"RottFuncs::RottFuncs(" << shapestr << ")");
+
+    switch (shape) {
+    case CIRC: 
+        TRACE(10,"Fptr set to circ")
+            f_ptr=&f_circ;        
+        f_ptrc=&f_circ;
+        break;
+    case VERT:
+        TRACE(10,"Fptr set to vert");
+        f_ptr=&f_vert<vc>;
+        f_ptrc=&f_vert<c>;
+        break;
+    case SQUARE:
+        TRACE(10,"Fptr set to square");
+        f_ptr=&f_square;
+        f_ptrc=&f_square;
+        break;
+    case BLAPPROX:
+        TRACE(10,"Fptr set to blapprox");
+        f_ptr=&f_blapprox<vc>;
+        f_ptrc=&f_blapprox<c>;
+        break;
+    case INVISCID:
+        TRACE(10,"Fptr set to inviscid");
+        f_ptr=&f_inviscid<vc>;
+        f_ptrc=&f_inviscid<c>;
+        break;
+    default:
+        FATAL("Unhandled shape");
+        break;
+    }
+}
+
+//////////////////////////////////////////////////////////////////////
--- a/src/funcs/rottfuncs.h
+++ b/src/funcs/rottfuncs.h
@ -0,0 +1,43 @@
+// rottfuncs.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Implementation of thermoviscous (Rott's) functions
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef ROTTFUNCS_H
+#define ROTTFUNCS_H
+#include "tasmet_types.h"
+#include "tasmet_tracer.h"
+#include "tasmet_enum.h"
+
+DECLARE_ENUM(cshape,INVISCID,BLAPPROX,VERT,SQUARE,CIRC)
+
+class RottFuncs {		// Thermoviscous Rott functions
+    vc (*f_ptr)(const vc& rh_over_delta);
+    c (*f_ptrc)(const c& rh_over_delta);
+    cshape _shape;
+public:
+    RottFuncs(const cshape shape=cshape::INVISCID);
+    RottFuncs(const RottFuncs& other) =delete;
+    RottFuncs& operator=(const RottFuncs&) =delete;
+    ~RottFuncs(){}
+
+    vc operator()(const vc& rh_over_delta) const {
+        TRACE(5,"fx vc");
+        return f_ptr(rh_over_delta);
+    }
+    vc operator()(const vd& rh_over_delta) const {
+        TRACE(5,"fx vd");
+        return f_ptr((1.0+0.0*I)*rh_over_delta);
+    }
+
+    c operator()(const c& rh_over_delta) const {return f_ptrc(rh_over_delta);}
+    c operator()(const d& rh_over_delta) const {c n(rh_over_delta,0); return f_ptrc(n);}
+    cshape getCshape() const {return _shape;}
+};
+
+
+#endif // ROTTFUNCS_H
+//////////////////////////////////////////////////////////////////////
--- a/src/funcs/skewsine.cpp
+++ b/src/funcs/skewsine.cpp
@ -0,0 +1,27 @@
+// skewsine.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+// Implementation of a skew sine
+//////////////////////////////////////////////////////////////////////
+
+#include "skewsine.h"
+
+// Equation for a skew sine between x=0 and x=1
+d skewsine(d x){
+    if (x<0.0)
+        return 0.0;
+    else if (x<1.0)
+        return x-(1.0/(2.0*number_pi))*sin(x*number_pi*2.0);
+    else
+        return 1.0;
+}
+
+vd skewsine(const vd& x){
+    vd result=x;
+    result.transform( [](d val) {return skewsine(val);});
+    return result;
+}
+
+//////////////////////////////////////////////////////////////////////
--- a/src/funcs/skewsine.h
+++ b/src/funcs/skewsine.h
@ -0,0 +1,13 @@
+#pragma once
+
+#ifndef _SKEWSINE_H_
+#define _SKEWSINE_H_
+#include "tasmet_types.h"
+
+
+d skewsine(d x); 		// Skew sine between x=0 and x=1
+vd skewsine(const vd& x); 		// Skew sine between x=0 and x=1  
+
+
+
+#endif /* _SKEWSINE_H_ */
--- a/src/funcs/sph_j.cpp
+++ b/src/funcs/sph_j.cpp
@ -0,0 +1,32 @@
+// sph_j.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#include <gsl/gsl_sf_bessel.h>
+#include "sph_j.h"
+#include "tasmet_exception.h"
+
+
+vd sph_j(int n,d x){
+    if (n < 1)
+        throw TaSMETError("n needs to be > 1");
+    vd result(n);
+    gsl_sf_bessel_jl_steed_array(n,x,result.memptr());
+    return result;
+}
+
+dmat sph_j(int n, const vd& x){
+    if (n < 1)
+        throw TaSMETError("n needs to be > 1");
+
+    dmat result(x.size(),n);
+    for (us i = 0; i < x.size(); i++)
+        result.row(i) = sph_j(n, x(i)).t();
+    return result;
+}
+
+
+//////////////////////////////////////////////////////////////////////
--- a/src/funcs/sph_j.h
+++ b/src/funcs/sph_j.h
@ -0,0 +1,32 @@
+// sph_j.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Implementation of spherical bessel functions
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef SPH_J_H
+#define SPH_J_H
+#include "tasmet_types.h"
+
+
+// Returns the first n spherical Bessel functions for a given double value
+vd sph_j(int n,d x);
+// Returns a 2D array. The first axis is the value of x and the
+// second axis (say the column number) corresponds to the n'th
+// spherical Bessel function value. Hence the result vector looks like:
+
+// [ j_0(x_0) j_1(x_0) ... j_n(x_0) ]
+// [ j_0(x_1) j_1(x_1_ ... j_n(x_1) ]
+//     .       .
+//     .       .
+//     .       .
+// [ j_0(x_n) j_1(x_n) ... j_n(x_n) ]
+  
+dmat sph_j(int n, const vd& x);
+
+
+
+#endif // SPH_J_H
+//////////////////////////////////////////////////////////////////////
--- a/src/material/CMakeLists.txt
+++ b/src/material/CMakeLists.txt
@ -0,0 +1,7 @@
+add_library(material
+  gas.cpp
+  air.cpp
+  helium.cpp
+  nitrogen.cpp
+  solid.cpp
+  )
--- a/src/material/air.cpp
+++ b/src/material/air.cpp
@ -0,0 +1,54 @@
+// air.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "air.h"
+
+namespace  {
+    // Air-specific data
+    d kappac[]={-0.00227583562,1.15480022E-4,   \
+                -7.90252856E-8,4.11702505E-11,  \
+                -7.43864331E-15};
+
+    d cpc[]={1047.63657,-0.372589265,           \
+             9.45304214E-4,-6.02409443E-7,      \
+             1.2858961E-10};
+    d muc[]={-8.38278E-7,8.35717342E-8,         \
+             -7.69429583E-11,4.6437266E-14,     \
+             -1.06585607E-17};
+} // namespace
+  
+d Air::h(d T,d p) const {
+    return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5);
+}
+vd Air::h(const vd& T,const vd& p) const {
+    return cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5);
+}
+vd Air::cp(const vd& T,const vd& p) const {
+    return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4);
+}
+d Air::cp(d T,d p) const {
+    return cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4);
+}
+
+vd Air::kappa(const vd& T,const vd& p) const {
+    return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4);
+}
+d Air::kappa(d T,d p) const {
+    return kappac[0]+kappac[1]*T+kappac[2]*pow(T,2)+kappac[3]*pow(T,3)+kappac[4]*pow(T,4);
+}
+vd Air::mu(const vd& T,const vd& p) const {
+    return muc[1]*T+muc[2]*pow(T,2)+muc[3]*pow(T,3)+muc[4]*pow(T,4)+muc[0];
+}
+d Air::mu(d T,d p) const {
+    return muc[1]*T+muc[2]*pow(T,2)+muc[3]*pow(T,3)+muc[4]*pow(T,4)+muc[0];
+}
+
+  
+
+
+//////////////////////////////////////////////////////////////////////
--- a/src/material/air.h
+++ b/src/material/air.h
@ -0,0 +1,34 @@
+// air.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Implementation of the gas air
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef AIR_H
+#define AIR_H
+#include "perfectgas.h"
+
+
+class Air : public PerfectGas {
+protected:
+    d Rs() const {return 287;}
+public:
+    Air():PerfectGas(air){}
+    d cp(d T,d p) const;
+    vd cp(const vd& T,const vd& p) const;
+    d h(d T,d p) const;
+    vd h(const vd& T,const vd& p) const;
+    d mu(d T,d p) const;
+    vd mu(const vd& T,const vd& p) const;
+    d kappa(d T,d p) const;
+    vd kappa(const vd& T,const vd& p) const;
+    ~Air(){}
+};
+
+  
+
+
+#endif // AIR_H
+//////////////////////////////////////////////////////////////////////
--- a/src/material/gas.cpp
+++ b/src/material/gas.cpp
@ -0,0 +1,34 @@
+// gas.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#include "tasmet_tracer.h"
+#include "tasmet_exception.h"
+#include "gas.h"
+#include "helium.h"
+#include "air.h"
+#include "nitrogen.h"
+
+Gas* Gas::newGas(const GasType gastype) {
+
+    switch (gastype) {
+    case helium:
+        return new Helium();
+        break;
+    case air:
+        return new Air();
+    case nitrogen:
+        return new Nitrogen();
+    default:
+        FATAL("Gas type not known");
+        break;
+    }
+    return nullptr;
+}
+
+
+//////////////////////////////////////////////////////////////////////
+
--- a/src/material/gas.h
+++ b/src/material/gas.h
@ -0,0 +1,96 @@
+// gas.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Gas properties interface
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef GAS_H
+#define GAS_H
+#include <string>
+#include <armadillo>
+#include "tasmet_enum.h"
+#include "tasmet_types.h"
+
+
+
+class Gas{
+public:
+    DECLARE_ENUM(GasType,air,helium,nitrogen);
+private:    
+    GasType _gastype;
+protected:
+    Gas(GasType gastype):_gastype(gastype){}
+public:
+
+
+
+    Gas(const Gas& other) =delete;
+    Gas& operator=(const Gas&)=delete;
+    virtual ~Gas(){}
+
+    // Static method to generate a Gas
+    static Gas* newGas(const GasType gastype);
+    operator GasType() { return _gastype;}
+
+    // Dimensionless numbers:
+
+    // Specific heat ratio
+    d gamma(d T,d p) const {return cp(T,p)/cv(T,p);}
+    vd gamma(const vd& T,const vd& p) const {return cp(T,p)/cv(T,p);}
+
+    // Prandtl number
+    vd pr(const vd& T,const vd& p) const {return mu(T,p)%cp(T,p)/kappa(T,p);}
+    d pr(d T,d p) const {return mu(T,p)*cp(T,p)/kappa(T,p);}
+
+    // Virtuals that are part of the interface
+    // Density [kg/m^3]
+    virtual d rho(d T,d p) const=0;
+    virtual vd rho(const vd& T,const vd& p) const=0;
+
+    // Internal energy [J/kg]
+    virtual d e(d T,d p) const=0;
+    virtual vd e(const vd& T,const vd& p) const=0;
+
+    // Static enthalpy [J/kg]
+    virtual d h(d T,d p) const=0;
+    virtual vd h(const vd& T,const vd& p) const=0;
+
+    // Specific heat at constant pressure
+    virtual d cp(d T,d p) const=0;
+    virtual vd cp(const vd& T,const vd& p) const=0;
+
+    // Specific heat at constant density
+    virtual d cv(d T,d p) const=0;
+    virtual vd cv(const vd& T,const vd& p) const=0;
+
+    virtual d beta(d T,d p) const=0;
+    virtual vd beta(const vd& T,const vd& p) const=0;
+
+    // Adiabatic speed of sound [m/s]
+    virtual d cm(d T,d p) const=0;
+    virtual vd cm(const vd& T,const vd& p) const=0;
+
+    // Dynamic viscosity [Pa*s]
+    virtual d mu(d T,d p) const=0;
+    virtual vd mu(const vd& T,const vd& p) const=0;
+
+    // Thermal conductivity [W/mK]
+    virtual d kappa(d T,d p) const=0;
+    virtual vd kappa(const vd& T,const vd& p) const=0;
+
+
+
+
+
+
+
+
+
+};
+
+
+
+#endif // GAS_H
+//////////////////////////////////////////////////////////////////////
--- a/src/material/helium.cpp
+++ b/src/material/helium.cpp
@ -0,0 +1,33 @@
+// helium.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "helium.h"
+
+  
+vd Helium::cp(const vd& T,const vd& p) const {
+    return cp(0.0,0.0)*arma::ones(T.size());
+}
+vd Helium::h(const vd& T,const vd& p) const {
+    return cp(T,p)%T;
+}
+d Helium::mu(d T,d p) const {
+    return 0.412e-6*pow(T,0.68014);
+}
+d Helium::kappa(d T,d p) const {
+    return 0.0025672*pow(T,0.716);
+}
+vd Helium::mu(const vd& T,const vd& p) const {
+    return 0.412e-6*pow(T,0.68014);
+}
+vd Helium::kappa(const vd& T,const vd& p) const {
+    return 0.0025672*pow(T,0.716);
+}
+
+  
+
+//////////////////////////////////////////////////////////////////////
--- a/src/material/helium.h
+++ b/src/material/helium.h
@ -0,0 +1,39 @@
+// helium.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef HELIUM_H
+#define HELIUM_H
+#include "perfectgas.h"
+
+class Helium :public PerfectGas {
+protected:
+    d Rs() const {return 2070;}
+public:
+    Helium():PerfectGas(helium){}
+
+    d cp(d T,d p) const { return 5195;}
+    vd cp(const vd& T,const vd& p) const;
+
+    d h(d T,d p) const {return cp(0.0,0.0)*T;}
+    vd h(const vd& T,const vd& p) const;
+
+    d mu(d T,d p) const;
+    vd mu(const vd& T,const vd& p) const;
+
+    d kappa(d T,d p) const;
+    vd kappa(const vd& T,const vd& p) const;
+
+    virtual ~Helium(){}
+};
+
+  
+
+
+
+#endif // HELIUM_H
+//////////////////////////////////////////////////////////////////////
--- a/src/material/nitrogen.cpp
+++ b/src/material/nitrogen.cpp
@ -0,0 +1,49 @@
+// nitrogen.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "nitrogen.h"
+  
+namespace  {
+    // Nitrogen-specific data
+    d cpc[]={3.29868,0.00140824,                \
+             -3.96322e-6,5.64152e-9,            \
+             -2.44486e-12};
+
+    d kappavals[]={6.995e-3,0.0631e-3};
+} // namespace
+  
+d Nitrogen::h(d T,d p) const {
+    return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5));
+}
+vd Nitrogen::h(const vd& T,const vd& p) const {
+    return Rs()*(cpc[0]*T+0.5*cpc[1]*pow(T,2)+(1/3.0)*cpc[2]*pow(T,3)+cpc[3]*0.25*pow(T,4)+cpc[4]*(0.2)*pow(T,5));
+}
+vd Nitrogen::cp(const vd& T,const vd& p) const {
+    return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4));
+}
+d Nitrogen::cp(d T,d p) const {
+    return Rs()*(cpc[0]+cpc[1]*T+cpc[2]*pow(T,2)+cpc[3]*pow(T,3)+cpc[4]*pow(T,4));
+}
+// Document Mina
+vd Nitrogen::kappa(const vd& T,const vd& p) const {
+    return kappavals[1]*T+kappavals[0];
+}
+d Nitrogen::kappa(d T,d p) const {
+    return kappavals[1]*T+kappavals[0];
+}
+// http://www.lmnoeng.com/Flow/GasViscosity.php
+vd Nitrogen::mu(const vd& T,const vd& p) const {
+    return (0.01781/1000)*(411.55/(T+111))%pow(T/300.55,1.5);
+}
+d Nitrogen::mu(d T,d p) const {
+    return (0.01781/1000)*(411.55/(T+111))*pow(T/300.55,1.5);
+}
+
+  
+
+//////////////////////////////////////////////////////////////////////
--- a/src/material/nitrogen.h
+++ b/src/material/nitrogen.h
@ -0,0 +1,36 @@
+// nitrogen.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef NITROGEN_H
+#define NITROGEN_H
+
+#include "perfectgas.h"
+
+  
+class Nitrogen : public PerfectGas {
+protected:
+    d Rs() const {return 297;}
+public:
+    Nitrogen():PerfectGas(nitrogen){}
+
+    d cp(d T,d p) const;
+    vd cp(const vd& T,const vd& p) const;
+    d h(d T,d p) const;
+    vd h(const vd& T,const vd& p) const;
+    d mu(d T,d p) const;
+    vd mu(const vd& T,const vd& p) const;
+    d kappa(d T,d p) const;
+    vd kappa(const vd& T,const vd& p) const;
+    ~Nitrogen(){}
+};
+
+  
+
+#endif // NITROGEN_H
+
+//////////////////////////////////////////////////////////////////////
--- a/src/material/perfectgas.cpp
+++ b/src/material/perfectgas.cpp
@ -0,0 +1,7 @@
+// perfectgas.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
--- a/src/material/perfectgas.h
+++ b/src/material/perfectgas.h
@ -0,0 +1,78 @@
+// perfectgas.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef PERFECTGAS_H
+#define PERFECTGAS_H
+#include "gas.h"
+#ifdef DEBUG_GAS
+#define checkzero(x)                                    \
+    TRACE(-1,"Divide by zero testing entered.");	\
+    try {if(min(abs(x))<1e-13) throw 0; }               \
+    catch(int a){					\
+        TRACE(0,"Divide by zero encountered.");         \
+    }
+#else
+#define checkzero(x)
+#endif
+
+  
+class PerfectGas: public Gas {
+protected:
+    PerfectGas(GasType gastype): Gas(gastype){}
+
+    // Not implemented for a perfect gas:
+    // mu, kappa, h, cp, Rs
+
+    virtual d Rs() const=0;
+
+public:
+    ~PerfectGas(){}
+
+    vd rho(const vd&T,const vd&p) const {
+        checkzero(T);
+        return p/Rs()/T;
+    }
+    vd p(const vd& T,const vd& rho) const {
+        return rho%(Rs()*T);
+    }
+    vd cv(const vd& T,const vd& p) const {
+        return cp(T,p)-Rs();
+    }
+    vd e(const vd& T,const vd& p) const {
+        return h(T,p)-Rs()*T;
+    }
+    d beta(d T,d p) const {
+        checkzero(T);
+        return 1/T;
+    }
+    vd beta(const vd& T,const vd& p) const {
+        return 1/T;
+    }
+    vd cm(const vd& T,const vd& p) const {
+        return sqrt(gamma(T,p)*Rs()%T);
+    }
+    d rho(d T,d p) const {
+        checkzero(T);
+        return p/Rs()/T;
+    }
+    d cv(d T,d p) const {
+        return cp(T,p)-Rs();
+    }
+    d e(d T,d p) const {
+        return h(T,p)-Rs()*T;
+    }
+    d cm(d T,d p) const {
+        d csq=gamma(T,p)*Rs()*T;
+        return sqrt(csq);
+    }
+    
+};
+
+#endif // PERFECTGAS_H
+//////////////////////////////////////////////////////////////////////
+
--- a/src/material/solid.cpp
+++ b/src/material/solid.cpp
@ -0,0 +1,66 @@
+#include "solid.h"
+#include "tasmet_tracer.h"
+
+//Stainless steel
+//Container class
+Solid* Solid::newSolid(const SolidType solidtype){
+    switch (solidtype) {
+    case stainless: 
+        return new Stainless();
+        break;
+    case kapton:
+        return new Kapton();
+        break;
+    case copper:
+        return new Copper();
+        break;
+    case stainless_hopkins:
+        return new StainlessHopkins();
+        break;
+    default:
+        FATAL("Undefined solid type");
+        break;
+    }
+    return nullptr;
+}
+
+vd StainlessHopkins::c(const vd& T) const{return 490*pow(T,0);}
+d StainlessHopkins::c(const d& T) const{return 490;}   
+vd StainlessHopkins::rho(const vd& T) const{return 7900*pow(T,0);}
+d StainlessHopkins::rho(const d& T) const{return 7900;}  
+vd StainlessHopkins::kappa(const vd& T) const{return 14.9*pow(T,0);}
+d StainlessHopkins::kappa(const d& T) const{return 14.9;}   
+
+   
+vd Stainless::c(const vd& T) const {
+    vd arg=1.7054e-6*pow(T,-0.88962)+22324.0/pow(T,6);
+    return pow(arg,-1.0/3.0)+15/T;
+}
+d Stainless::c(const d& T) const {
+    d arg=1.7054e-6*pow(T,-0.88962)+22324.0/pow(T,6);
+    return pow(arg,-1.0/3.0)+15/T;
+}
+vd Stainless::rho(const vd& T) const {
+    return 8274.55-1055.23*exp(-1.0*pow((T-273.15-2171.05)/2058.08,2));
+}
+d Stainless::rho(const d& T) const {
+    return 8274.55-1055.23*exp(-1.0*pow((T-273.15-2171.05)/2058.08,2));
+}
+vd Stainless::kappa(const vd& T) const {
+    return pow(266800.0*pow(T,-5.2)+0.21416*pow(T,-1.6),-0.25);
+}
+d Stainless::kappa(const d& T) const {
+    return pow(266800.0*pow(T,-5.2)+0.21416*pow(T,-1.6),-0.25);
+}
+
+//Copper
+vd Copper::kappa(const vd& T) const {
+    return 398.0-0.567*(T-300.0);
+}
+d Copper::kappa(const d& T) const { return 398.0-0.567*(T-300.0); }
+vd Copper::c(const vd& T) const {return 420.0*pow(T,0);}
+d Copper::c(const d& T) const {return 420.0*pow(T,0);}
+vd Copper::rho(const vd& T) const {return 9000.0*pow(T,0);}
+d Copper::rho(const d& T) const {return 9000.0*pow(T,0);}
+
+
--- a/src/material/solid.h
+++ b/src/material/solid.h
@ -0,0 +1,78 @@
+#pragma once
+#ifndef SOLID_H
+#define SOLID_H
+
+#include "tasmet_types.h"
+#include "tasmet_enum.h"
+
+
+class Solid{
+public:
+    DECLARE_ENUM(SolidType,stainless,stainless_hopkins,copper,kapton);
+private:
+    SolidType _solidtype;
+protected:
+    Solid(const SolidType);
+public:
+    Solid(const Solid&) = delete;
+    Solid& operator=(const Solid&) = delete;
+    operator SolidType() const {return _solidtype;}
+
+    Solid* newSolid(const SolidType);
+
+    virtual vd kappa(const vd& T) const=0;
+    virtual d kappa(const d& T) const=0;
+    virtual vd c(const vd& T) const=0;
+    virtual d c(const d& T) const=0;
+    virtual vd rho(const vd& T) const=0;
+    virtual d rho(const d& T) const=0;
+    virtual ~Solid();
+
+};
+
+class StainlessHopkins:public Solid
+{
+public:
+    StainlessHopkins(): Solid(stainless_hopkins){}
+    vd kappa(const vd& T) const;
+    d kappa(const d& T) const;
+    vd c(const vd& T) const;
+    d c(const d& T) const;
+    vd rho(const vd& T) const;
+    d rho(const d& T) const;
+};
+
+class Stainless: public Solid{
+public:
+    Stainless(): Solid(stainless){}
+    vd kappa(const vd& T) const;
+    d kappa(const d& T) const;
+    vd c(const vd& T) const;
+    d c(const d& T) const;
+    vd rho(const vd& T) const;
+    d rho(const d& T) const;
+
+};
+class Copper: public Solid{
+public:
+    Copper(): Solid(copper){}
+    vd kappa(const vd& T) const;
+    d kappa(const d& T) const;
+    vd c(const vd& T) const;
+    d c(const d& T) const;
+    vd rho(const vd& T) const;
+    d rho(const d& T) const;
+};
+
+class Kapton: public Solid {
+public:
+    Kapton():Solid(kapton){}
+    vd kappa(const vd& T) const { return 0.2*(1.0-exp(-T/100.0));}
+    d kappa(const d& T) const { return 0.2*(1.0-exp(-T/100.0));}
+    vd c(const vd& T) const {return 3.64*T;}
+    d c(const d& T) const {return 3.64*T;}
+    vd rho(const vd& T) const {return 1445.0-0.085*T;}
+    d rho(const d& T) const {return 1445.0-0.085*T;}
+};
+
+#endif /* SOLID_H */
--- a/src/tasmet_consolecolors.h
+++ b/src/tasmet_consolecolors.h
@ -0,0 +1,27 @@
+// consolecolors.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Print text from C++ to stdout in color
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef CONSOLECOLORS_H
+#define CONSOLECOLORS_H
+#include <iostream>
+
+#ifndef SWIG
+#define red "\e[31m"
+#define green "\e[32m"
+#define def " \e[39m"
+
+#endif  // SWIG
+
+// Command to clear the content of a console
+inline void clearConsole(){
+  std::cout << "\033c" << std::endl;
+}
+
+#endif // CONSOLECOLORS_H
+//////////////////////////////////////////////////////////////////////
+
--- a/src/tasmet_enum.h
+++ b/src/tasmet_enum.h
@ -0,0 +1,61 @@
+// tasmet_enum.h
+//
+// Author: J.A. de Jong 
+//
+// Description: Copied and adapted code from
+// http://stackoverflow.com/questions/5093460/how-to-
+// convert-an-enum-type-variable-to-a-string
+
+// Used to automatically couple a enum to the describing string.
+// Use it as follows:
+// 
+// #include "tasmet_enum.h"
+// DECLARE_ENUM(OS_type, Windows, Linux, Apple)
+
+// void main() {
+//     cout << OS_typeToString(Windows) << endl;
+//     if (StringToOS_type("Linux") == Linux) { cout << "yay!" << endl; }
+//     cout << "Max: " << MAX_NUMBER_OF_OS_type << endl;
+// }
+
+#pragma once
+#ifndef TASMET_ENUM_H
+#define TASMET_ENUM_H
+#include <string>
+#include <sstream>
+
+// Search and remove whitespace from both ends of the string
+inline std::string TrimEnumString(const std::string &s)
+{
+    std::string::const_iterator it = s.begin();
+    while (it != s.end() && isspace(*it)) { it++; }
+    std::string::const_reverse_iterator rit = s.rbegin();
+    while (rit.base() != it && isspace(*rit)) { rit++; }
+    return std::string(it, rit.base());
+}
+inline void SplitEnumArgs(const char* szArgs, std::string Array[], int nMax)
+{
+    std::stringstream ss(szArgs);
+    std::string strSub;
+    int nIdx = 0;
+    while (ss.good() && (nIdx < nMax)) {
+        getline(ss, strSub, ',');
+        Array[nIdx] = TrimEnumString(strSub);
+        nIdx++;
+    }
+};
+#define DECLARE_ENUM(ename, ...) \
+    enum ename { __VA_ARGS__, MAX_NUMBER_OF_##ename }; \
+    static std::string ename##Strings[MAX_NUMBER_OF_##ename]; \
+    inline const char* ename##ToString(ename e) { \
+        if (ename##Strings[0].empty()) { SplitEnumArgs(#__VA_ARGS__, ename##Strings, MAX_NUMBER_OF_##ename); } \
+        return ename##Strings[e].c_str(); \
+    } \
+    inline ename StringTo##ename(const char* szEnum) { \
+        for (int i = 0; i < MAX_NUMBER_OF_##ename; i++) { if (ename##Strings[i] == szEnum) { return (ename)i; } } \
+        return MAX_NUMBER_OF_##ename; \
+    }
+
+
+#endif // TASMET_ENUM_H
+//////////////////////////////////////////////////////////////////////
--- a/src/tasmet_exception.cpp
+++ b/src/tasmet_exception.cpp
@ -0,0 +1,18 @@
+// tasmet_exception.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+// Implementation of exception TasMETBadAlloc
+//////////////////////////////////////////////////////////////////////
+
+#include "tasmet_exception.h"
+
+const char*  TasMETBadAlloc::what() const throw() {
+	return "Error: memory allocation failed. "
+   		"Please make sure enough memory is available and restart the application";
+}
+
+
+//////////////////////////////////////////////////////////////////////
+
--- a/src/tasmet_exception.h
+++ b/src/tasmet_exception.h
@ -0,0 +1,26 @@
+// tasmet_exception.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_EXCEPTION_H
+#define TASMET_EXCEPTION_H
+
+#include <string>
+#include <stdexcept>
+
+class TaSMETError : public std::runtime_error {
+public:
+    TaSMETError(const std::string& msg = "") : std::runtime_error(msg) {}
+};
+
+class TasMETBadAlloc: public std::bad_alloc {
+	virtual const char* what() const throw();
+};
+
+#endif // TASMET_EXCEPTION_H
+//////////////////////////////////////////////////////////////////////
+
--- a/src/tasmet_io.h
+++ b/src/tasmet_io.h
@ -0,0 +1,25 @@
+// tasmet_io.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_IO_H
+#define TASMET_IO_H
+#include <iostream>
+#include <iomanip>
+
+//Spoiling global namespace with often used functions and variables 
+using std::cout; /* Output to stdout */				
+using std::endl; 
+using std::setprecision; 
+using std::setiosflags; 
+using std::cin;		/* Input from stdin */ 
+using std::ostream; 
+using std::cerr; 
+using std::ios; 
+
+#endif // TASMET_IO_H
+//////////////////////////////////////////////////////////////////////
--- a/src/tasmet_tracer.cpp
+++ b/src/tasmet_tracer.cpp
@ -0,0 +1,13 @@
+// tasmet_tracer.cpp
+//
+// last-edit-by: J.A. de Jong 
+// 
+// Description:
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "tasmet_tracer.h"
+
+
+
+//////////////////////////////////////////////////////////////////////
--- a/src/tasmet_tracer.h
+++ b/src/tasmet_tracer.h
@ -0,0 +1,89 @@
+// tasmet_tracer.h
+//
+// Author: J.A. de Jong 
+// Description:
+// When this file is included, we are able to use the following macro's:
+// TRACE(10, "Test");
+// 
+// double x=2;
+// VARTRACE(2,x);
+//
+// WARN("Warning message");
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_TRACER_H
+#define TASMET_TRACER_H
+#include <type_traits>
+
+// **************************************** Tracer code
+// If PP variable TRACER is not defined, we automatically set it on.
+#ifndef TRACER
+static_assert(false,"TRACER macro not defined");
+#endif
+
+// Not so interesting part
+#define rawstr(x) #x
+#define namestr(x) rawstr(x)
+#define annestr(x) namestr(x)
+#define FILEWITHOUTPATH ( strrchr(__FILE__, '/') ? strrchr(__FILE__, '/') + 1 : strrchr(__FILE__, '\\') ? strrchr(__FILE__, '\\') + 1 : __FILE__ )
+#define POS FILEWITHOUTPATH << ":" << __LINE__ <<  ": "
+// End not so interesting part
+
+#define RAWWARNING(a) std::cout << red << a << def << "\n";
+#define WARN(a) RAWWARNING(POS << "WARNING: " << a)
+
+// SHOULD NOT BE USED IN A LIBRARY!!
+#define FATAL(a) WARN(a); abort(); 
+
+#if TRACER == 1
+#include "tasmet_consolecolors.h"
+#include <iostream>
+
+// This variable can be used to increase the trace level of one
+// specific file at compile time
+#ifndef TRACERPLUS
+#define TRACERPLUS (0)
+#endif
+
+
+// Initialize MAXTRACELEVEL and BUILDINTRACELEVEL
+#ifndef MAXTRACELEVEL
+#define MAXTRACELEVEL (5000) 	// To which the TRACELEVEL initially is set
+#endif
+
+// Define this preprocessor definition to overwrite
+// Use -O flag for compiler to remove the dead functions!
+// In that case all cout's for TRACE() are removed from code
+#ifndef BUILDINTRACELEVEL
+#define BUILDINTRACELEVEL (-10)
+#endif
+
+extern int tasmet_tracer_level;
+// Use this preprocessor command to introduce one tasmet_tracer_level integer per unit
+/* Introduce one static logger */
+// We trust that the compiler will eliminate 'dead code', which means
+// that if variable BUILDINTRACERLEVEL is set, the inner if statement
+// will not be reached.
+#define WRITETRACE(level,a) \
+    if(level>=BUILDINTRACELEVEL && level >=tasmet_tracer_level) {       \
+        std::cout << a << "\n";                                         \
+    }
+
+
+#define TRACE(l,a) WRITETRACE(l+TRACERPLUS,annestr(tasmet_tracer_level) << "-" << (l)<< "-" << POS << a)
+#define VARTRACE(l,a) TRACE(l,#a " = " << a)
+
+#define INITTRACE(ll)                                                   \
+    std::cout << "INITTRACE with tracelevel " << ll << " for " << annestr(tasmet_tracer_level) << "\n"; \
+    tasmet_tracer_level=ll;
+#else  // TRACER !=1
+
+#define VARTRACE(l,a) 
+#define TRACE(l,a)
+#define INITTRACE(a)
+#define WARN(a)
+
+
+#endif	// ######################################## TRACER ==1
+#endif // TASMET_TRACER_H
+//////////////////////////////////////////////////////////////////////
--- a/src/tasmet_types.h
+++ b/src/tasmet_types.h
@ -0,0 +1,81 @@
+// tasmet_types.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Typedefs and namespace pollution for stuff that is almost always
+// needed.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef TASMET_TYPES_H
+#define TASMET_TYPES_H
+#include <type_traits>
+
+// We require C++
+#ifndef __cplusplus
+#error The c++ compiler should be used.
+#endif
+
+// Some header files I (almost) always need
+#include <vector>		// C++ std vector
+#include <string>		// C++ string class
+#include <complex>
+#include <armadillo>		// Our linear algebra package
+
+typedef size_t us;		  /* Size type I always use */
+
+#define fillwith arma::fill	    // Fill vector or matrix with ...
+
+// To change the whole code to 32-bit floating points, change this to
+// float.
+#if TASMET_FLOAT == 32
+// typedef float d;		/* Shortcut for double */			
+#elif TASMET_FLOAT == 64
+typedef double d;		/* Shortcut for double */			
+#else
+static_assert(false,"TASMET_FLOAT should be either 32 or 64");
+#endif
+
+typedef std::complex<d> c; /* Shortcut for c++ complex number */
+
+using arma::abs;
+using std::string;		/* C++ Strings */				
+using std::vector;              // C++ Vectors
+
+/* Armadillo functions and objects we often need */				
+using arma::sin;
+using arma::cos;
+using arma::exp;
+using arma::ones;								
+using arma::zeros;								
+using arma::eye;								
+using arma::inv;								
+using arma::det;								
+using arma::log_det;								
+using arma::norm;								
+using arma::span;								
+using arma::linspace;								
+using arma::diagmat;								
+
+// Constants
+// I need these numbers so often, that they can be in the global namespace
+
+const c I(0,1); //Complex unity
+const c sqI=sqrt(I);
+const d sq2=sqrt(2.0);
+const c sqmI=sqrt(-1.0*I);
+const c minI(0,-1); //Minus complex unity
+const d number_pi=arma::datum::pi; // The number 3.14159265359..
+
+
+typedef arma::Col<d> vd;		  /* Column vector of doubles */		
+typedef arma::Col<c> vc;	  /* Column vector of complex numbers */	
+typedef arma::Mat<d> dmat; /* (Dense) Matrix of doubles */		
+typedef arma::Mat<c> cmat;	  /* (Dense) matrix of complex numbers */	
+
+typedef arma::SpMat<d> sdmat;	  /* Sparse matrix of doubles */		
+
+
+#endif // TASMET_TYPES_H
+//////////////////////////////////////////////////////////////////////
+
--- a/src/tasmet_utils.h
+++ b/src/tasmet_utils.h
@ -0,0 +1,71 @@
+// utils.h
+//
+// Author: J.A. de Jong 
+//
+// Description:
+// Some generic utils.
+//////////////////////////////////////////////////////////////////////
+#pragma once
+#ifndef UTILS_H
+#define UTILS_H
+#include <vector>
+#include <map>
+#include "tracer.h"
+#include "vtypes.h"
+#include <typeinfo>
+#include <exception>
+#include <cmath>
+
+
+// Purge a vector of components
+template<typename T>
+void purge(std::vector<T>& vec){
+    TRACE(10,"purge(vector)");
+    for (T& it: vec){
+        delete it;
+        it=nullptr;
+    }
+    vec.clear();
+}
+
+// Purge a vector of components
+template<typename Key,typename T>
+void purge(std::map<Key,T>& map){
+    TRACE(10,"purge(map)");
+    for (auto& it: map){
+        delete it.second;
+        it.second=nullptr;
+    }
+    map.clear();
+}
+
+template<typename T>
+const T& min(const T& x,const T& y)  {
+    return x<=y? x : y;
+}
+
+template<typename T>
+const T& max(const T& x,const T& y)  {
+    return x<=y? y : x;
+}
+
+template<typename SegType,typename Sys>
+SegType* copySeg(const SegType& t,const Sys& sys) {
+    SegType* newt=new SegType(t);
+    if(!newt){
+        WARN("Copying " << typeid(t).name() << "failed!");
+    }
+    return newt;
+} // copySeg
+template<typename T>
+void makeNormal(T& c) {
+    for(auto& val: c)
+        if(!std::isnormal(val))
+            val=0;
+}
+  
+} // namespace utils
+
+
+//////////////////////////////////////////////////////////////////////
+