python: initial code for python interface using cython

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

.gitignore

@@ -12,7 +12,6 @@ Makefile*
 *.pro.user*
 *.user
 *.orig
-openEMS
 localPaths.pri
 .directory
 
@@ -22,3 +21,8 @@ CMakeCache.txt
 cmake_install.cmake
 install_manifest.txt
 localConfig.cmake
+
+#python
+*.pyc
+*.pyo
+python/**/*.cpp

python/openEMS/_nf2ff.pxd

@@ -0,0 +1,38 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Dec 20 22:43:35 2015
+
+@author: thorsten
+"""
+
+from libcpp.string cimport string
+from libcpp.vector cimport vector
+from libcpp.complex cimport complex
+from libcpp cimport bool
+cimport cython.numeric
+
+cdef extern from "openEMS/nf2ff.h":
+    cdef cppclass cpp_nf2ff "nf2ff":
+        cpp_nf2ff(vector[float] freq, vector[float] theta, vector[float] phi, vector[float] center, unsigned int numThreads) except +
+
+        bool AnalyseFile(string E_Field_file, string H_Field_file)
+
+        void SetRadius(float radius)
+        void SetPermittivity(vector[float] permittivity);
+        void SetPermeability(vector[float] permeability);
+
+        void SetMirror(int _type, int _dir, float pos);
+
+        double GetTotalRadPower(size_t f_idx)
+        double GetMaxDirectivity(size_t f_idx)
+
+        complex[double]** GetETheta(size_t f_idx)
+        complex[double]** GetEPhi(size_t f_idx)
+        double** GetRadPower(size_t f_idx)
+
+        bool Write2HDF5(string filename)
+
+        void SetVerboseLevel(int level)
+
+cdef class _nf2ff:
+    cdef  cpp_nf2ff *thisptr

python/openEMS/_nf2ff.pyx

@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Dec 20 22:42:19 2015
+
+@author: thorsten
+"""
+
+cimport _nf2ff
+import numpy as np
+import os
+from CSXCAD.Utilities import CheckNyDir
+
+cdef class _nf2ff:
+    def __cinit__(self, freq, theta, phi, center, numThreads=0, **kw):
+        if type(freq) in [float, int]:
+            freq = list(float(freq))
+        if type(theta) in [float, int]:
+            theta = list(float(theta))
+        if type(phi) in [float, int]:
+            phi = list(float(phi))
+        self.thisptr = new cpp_nf2ff(freq, theta, phi, center, numThreads)
+
+        if 'verbose' in kw:
+            self.SetVerboseLevel(kw['verbose'])
+            del kw['verbose']
+
+        assert len(kw)==0, 'Unknown keyword(s): {}'.format(kw)
+
+    def AnalyseFile(self, e_file, h_file):
+        assert os.path.exists(e_file)
+        assert os.path.exists(h_file)
+        return self.thisptr.AnalyseFile(e_file.encode('UTF-8'), h_file.encode('UTF-8'))
+
+    def SetMirror(self, mirr_type, ny, pos):
+        if mirr_type<=0:
+            return
+        assert mirr_type<3
+        ny = CheckNyDir(ny)
+        self.thisptr.SetMirror(mirr_type, ny, pos)
+
+    def Write2HDF5(self, filename):
+        return self.thisptr.Write2HDF5(filename.encode('UTF-8'))
+
+    def SetVerboseLevel(self, level):
+        self.thisptr.SetVerboseLevel(level)

python/openEMS/nf2ff.py

@@ -0,0 +1,139 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Dec 20 20:53:12 2015
+
+@author: thorsten
+"""
+
+import os
+import numpy as np
+import h5py
+from openEMS import _nf2ff
+from openEMS import utilities
+
+class nf2ff:
+    def __init__(self, CSX, name, start, stop, **kw):
+        self.CSX   = CSX
+        self.name  = name
+        self.start = start
+        self.stop  = stop
+
+        self.freq  = None
+        self.theta = None
+        self.phi   = None
+        self.center = None
+
+        self.directions = [True]*6 # all directions by default
+        if 'directions' in kw:
+            self.directions = kw['directions']
+            del kw['directions']
+            assert len(self.directions)==6
+
+        self.mirror = [0]*6
+        if 'mirror' in kw:
+            self.mirror = kw['mirror']
+            del kw['mirror']
+            assert len(self.mirror)==6
+
+        self.dump_type = 0  # default Et/Ht
+        self.dump_mode = 1  # default cell interpolated
+
+        self.freq = None   # broadband recording by defualt
+        if 'frequency' in kw:
+            self.freq = kw['frequency']
+            del kw['frequency']
+            self.dump_type = 10  # Ef/Hf
+
+        self.e_file = '{}_E'.format(self.name)
+        self.h_file = '{}_H'.format(self.name)
+
+        self.e_dump = CSX.AddDump(self.e_file, dump_type=self.dump_type  , dump_mode=self.dump_mode, file_type=1, **kw)
+        self.h_dump = CSX.AddDump(self.h_file, dump_type=self.dump_type+1, dump_mode=self.dump_mode, file_type=1, **kw)
+        if self.freq is not None:
+            self.e_dump.SetFrequency(self.freq)
+            self.h_dump.SetFrequency(self.freq)
+
+#        print(self.directions)
+        for ny in range(3):
+            pos = 2*ny
+            if self.directions[pos]:
+                l_start = np.array(start)
+                l_stop  = np.array(stop)
+                l_stop[ny] = l_start[ny]
+                CSX.AddBox(self.e_dump, l_start, l_stop)
+                CSX.AddBox(self.h_dump, l_start, l_stop)
+            if self.directions[pos+1]:
+                l_start = np.array(start)
+                l_stop  = np.array(stop)
+                l_start[ny] = l_stop[ny]
+                CSX.AddBox(self.e_dump, l_start, l_stop)
+                CSX.AddBox(self.h_dump, l_start, l_stop)
+
+    def CalcNF2FF(self, sim_path, freq, theta, phi, center=[0,0,0], read_cached=True, verbose=0):
+        if not hasattr(freq, "__iter__"):
+            freq = [freq]
+        self.freq  = freq
+        self.theta = theta
+        self.phi   = phi
+        self.center = center
+
+        fn = os.path.join(sim_path, self.name + '.h5')
+        if os.path.exists(fn) and read_cached:
+            self.ReadNF2FF(sim_path)
+            return
+
+        nfc = _nf2ff._nf2ff(self.freq, np.deg2rad(theta), np.deg2rad(phi), center, verbose=verbose)
+
+        for ny in range(3):
+            nfc.SetMirror(self.mirror[2*ny]  , ny, self.start[ny])
+            nfc.SetMirror(self.mirror[2*ny+1], ny, self.stop[ny])
+
+        for n in range(6):
+            fn_e = os.path.join(sim_path, self.e_file + '_{}.h5'.format(n))
+            fn_h = os.path.join(sim_path, self.h_file + '_{}.h5'.format(n))
+            if os.path.exists(fn_e) and os.path.exists(fn_h):
+                assert nfc.AnalyseFile(fn_e, fn_h)
+
+        nfc.Write2HDF5(fn)
+
+        self.ReadNF2FF(sim_path)
+
+    def ReadNF2FF(self, sim_path):
+        h5_file = h5py.File(os.path.join(sim_path, self.name + '.h5'), 'r')
+        mesh_grp = h5_file['Mesh']
+        phi   = np.array(mesh_grp['phi'])
+        theta = np.array(mesh_grp['theta'])
+        self.r     = np.array(mesh_grp['r'])
+
+        data  = h5_file['nf2ff']
+        freq = np.array(data.attrs['Frequency'])
+
+        if self.phi is not None:
+            assert utilities.Check_Array_Equal(np.rad2deg(theta), self.theta, 1e-4)
+            assert utilities.Check_Array_Equal(np.rad2deg(phi), self.phi, 1e-4)
+            assert utilities.Check_Array_Equal(freq, self.freq, 1e-6, relative=True)
+        self.Dmax = np.array(data.attrs['Dmax'])
+        self.Prad = np.array(data.attrs['Prad'])
+
+        THETA, PHI = np.meshgrid(theta, phi, indexing='ij')
+        cos_phi = np.cos(PHI)
+        sin_phi = np.sin(PHI)
+
+        self.E_theta = []
+        self.E_phi   = []
+        self.P_rad   = []
+        self.E_norm  = []
+        self.E_cprh  = []
+        self.E_cplh  = []
+        for n in range(len(freq)):
+            E_theta = np.array(h5_file['/nf2ff/E_theta/FD/f{}_real'.format(n)]) + 1j*np.array(h5_file['/nf2ff/E_theta/FD/f{}_imag'.format(n)])
+            E_theta = np.swapaxes(E_theta, 0, 1)
+            E_phi   = np.array(h5_file['/nf2ff/E_phi/FD/f{}_real'.format(n)])   + 1j*np.array(h5_file['/nf2ff/E_phi/FD/f{}_imag'.format(n)])
+            E_phi   = np.swapaxes(E_phi, 0, 1)
+            self.P_rad  .append(np.swapaxes(np.array(h5_file['/nf2ff/P_rad/FD/f{}'.format(n)]), 0, 1))
+
+            self.E_theta.append(E_theta)
+            self.E_phi  .append(E_phi)
+            self.E_norm .append(np.sqrt(np.abs(E_theta)**2 + np.abs(E_phi)**2))
+            self.E_cprh .append((cos_phi+1j*sin_phi) * (E_theta+1j*E_phi)/np.sqrt(2.0))
+            self.E_cplh .append((cos_phi-1j*sin_phi) * (E_theta-1j*E_phi)/np.sqrt(2.0))

python/openEMS/openEMS.pxd

@@ -0,0 +1,49 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Dec 14 00:12:43 2015
+
+@author: thorsten
+"""
+
+from libcpp.string cimport string
+from libcpp cimport bool
+
+from CSXCAD.CSXCAD cimport _ContinuousStructure, ContinuousStructure
+
+cdef extern from "openEMS/openems.h":
+    cdef cppclass _openEMS "openEMS":
+        _openEMS() except +
+        void SetNumberOfTimeSteps(unsigned int val)
+        void SetCSX(_ContinuousStructure* csx)
+
+        void SetEndCriteria(double val)
+        void SetOverSampling(int val)
+        void SetCellConstantMaterial(bool val)
+
+        void SetCylinderCoords(bool val)
+        #void SetupCylinderMultiGrid(std::vector<double> val)
+
+        void SetTimeStepMethod(int val)
+        void SetTimeStep(double val)
+        void SetTimeStepFactor(double val)
+        void SetMaxTime(double val)
+
+        void Set_BC_Type(int idx, int _type)
+        int Get_BC_Type(int idx)
+        void Set_BC_PML(int idx, unsigned int size)
+        int Get_PML_Size(int idx)
+        void Set_Mur_PhaseVel(int idx, double val)
+
+        void SetGaussExcite(double f0, double fc)
+
+        void SetVerboseLevel(int level)
+
+        int SetupFDTD()
+        void RunFDTD()
+
+        @staticmethod
+        void WelcomeScreen()
+
+cdef class openEMS:
+    cdef  _openEMS *thisptr
+    cdef readonly ContinuousStructure CSX      # hold a C++ instance which we're wrapping

python/openEMS/openEMS.pyx

@@ -0,0 +1,176 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Dec 13 23:50:24 2015
+
+@author: thorsten
+"""
+
+import os, sys, shutil
+import numpy as np
+cimport openEMS
+from . import ports, nf2ff
+
+cdef class openEMS:
+    @staticmethod
+    def WelcomeScreen():
+        _openEMS.WelcomeScreen()
+
+    """ __cinit__
+    NrTS:           max. number of timesteps to simulate (e.g. default=1e9)
+    EndCriteria:    end criteria, e.g. 1e-5, simulations stops if energy has decayed by this value (<1e-4 is recommended, default=1e-5)
+    MaxTime:        max. real time in seconds to simulate
+    OverSampling:   nyquist oversampling of time domain dumps
+    CoordSystem:    choose coordinate system (0 Cartesian, 1 Cylindrical)
+    #MultiGrid:      define a cylindrical sub-grid radius
+    TimeStep:       force to use a given timestep (dangerous!)
+    TimeStepFactor: reduce the timestep by a given factor (>0 to <=1)
+    TimeStepMethod: 1 or 3 chose timestep method (1=CFL, 3=Rennigs (default))
+    CellConstantMaterial: set to 1 to assume a material is constant inside a cell (material probing in cell center)
+    """
+    def __cinit__(self, *args, **kw):
+        self.thisptr = new _openEMS()
+        self.CSX = None
+
+        if 'NrTS' in kw:
+            self.SetNumberOfTimeSteps(kw['NrTS'])
+        else:
+            self.SetNumberOfTimeSteps(1e9)
+        if 'EndCriteria' in kw:
+            self.SetEndCriteria(kw['EndCriteria'])
+        if 'MaxTime' in kw:
+            self.SetMaxTime(kw['MaxTime'])
+        if 'OverSampling' in kw:
+            self.SetOverSampling(kw['OverSampling'])
+        if 'CoordSystem' in kw:
+            self.SetCoordSystem(kw['CoordSystem'])
+        if 'TimeStep' in kw:
+            self.SetTimeStep(kw['TimeStep'])
+        if 'TimeStepFactor' in kw:
+            self.SetTimeStepFactor(kw['TimeStepFactor'])
+        if 'TimeStepMethod' in kw:
+            self.SetTimeStepMethod(kw['TimeStepMethod'])
+        if 'CellConstantMaterial' in kw:
+            self.SetCellConstantMaterial(kw['CellConstantMaterial'])
+
+    def __dealloc__(self):
+        del self.thisptr
+        if self.CSX is not None:
+            self.CSX.thisptr = NULL
+
+    def SetNumberOfTimeSteps(self, val):
+        self.thisptr.SetNumberOfTimeSteps(val)
+
+    def SetEndCriteria(self, val):
+        self.thisptr.SetEndCriteria(val)
+
+    def SetOverSampling(self, val):
+        self.thisptr.SetOverSampling(val)
+
+    def SetCellConstantMaterial(self, val):
+        self.thisptr.SetCellConstantMaterial(val)
+
+    def SetCoordSystem(self, val):
+        assert (val==0 or val==1)
+        if val==0:
+            pass
+        elif val==1:
+            self.SetCylinderCoords()
+
+    def SetCylinderCoords(self):
+        self.thisptr.SetCylinderCoords(True)
+
+    def SetTimeStepMethod(self, val):
+        self.thisptr.SetTimeStepMethod(val)
+    def SetTimeStep(self, val):
+        self.thisptr.SetTimeStep(val)
+    def SetTimeStepFactor(self, val):
+        self.thisptr.SetTimeStepFactor(val)
+    def SetMaxTime(self, val):
+        self.thisptr.SetMaxTime(val)
+
+    def SetGaussExcite(self, f0, fc):
+        self.thisptr.SetGaussExcite(f0, fc)
+
+
+    def SetBoundaryCond(self, BC):
+        assert len(BC)==6
+        for n in range(len(BC)):
+            if type(BC[n])==int:
+                self.thisptr.Set_BC_Type(n, BC[n])
+                continue
+            if BC[n] in ['PEC', 'PMC', 'MUR']:
+                self.thisptr.Set_BC_Type(n, ['PEC', 'PMC', 'MUR'].index(BC[n]))
+                continue
+            if BC[n].startswith('PML_'):
+                size = int(BC[n].strip('PML_'))
+                self.thisptr.Set_BC_PML(n, size)
+                continue
+            raise Exception('Unknown boundary condition')
+
+    def AddLumpedPort(self, port_nr, R, start, stop, p_dir, excite=0, **kw):
+        assert self.CSX is not None
+        return ports.LumpedPort(self.CSX, port_nr, R, start, stop, p_dir, excite, **kw)
+
+    def AddRectWaveGuidePort(self, port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw):
+        assert self.CSX is not None
+        return ports.RectWGPort(self.CSX, port_nr, start, stop, p_dir, a, b, mode_name, excite, **kw)
+
+    def AddMSLPort(self, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw):
+        assert self.CSX is not None
+        return ports.MSLPort(self.CSX, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite, **kw)
+
+    def CreateNF2FFBox(self, name='nf2ff', start=None, stop=None, **kw):
+        assert self.CSX is not None
+        directions = [True]*6
+        mirror     = [0]*6
+        BC_size = [0]*6
+        BC_type = [0]*6
+        for n in range(6):
+            BC_type[n] = self.thisptr.Get_BC_Type(n)
+            if BC_type[n]==0:
+                directions[n]= False
+                mirror[n]    = 1  # PEC mirror
+            elif BC_type[n]==1:
+                directions[n]= False
+                mirror[n]    = 2  # PMC mirror
+            elif BC_type[n]==2:
+                BC_size[n] = 2
+            elif BC_type[n]==3:
+                BC_size[n] = self.thisptr.Get_PML_Size(n)+1
+
+        if start is None or stop is None:
+            grid = self.CSX.GetGrid()
+            assert grid.IsValid(), 'Error::CreateNF2FFBox: Grid is invalid'
+            start = np.zeros(3)
+            stop  = np.zeros(3)
+            for n in range(3):
+                l = grid.GetLines(n)
+                BC_type = self.thisptr.Get_BC_Type(2*n)
+                assert len(l)>(BC_size[2*n]+BC_size[2*n+1]), 'Error::CreateNF2FFBox: not enough lines in some direction'
+                start[n] = l[BC_size[2*n]]
+                stop[n]  = l[-1*BC_size[2*n+1]-1]
+        return nf2ff.nf2ff(self.CSX, name, start, stop, directions=directions, mirror=mirror, **kw)
+
+    def SetCSX(self, ContinuousStructure CSX):
+        self.CSX = CSX
+        self.thisptr.SetCSX(CSX.thisptr)
+
+    def Run(self, sim_path, cleanup=False, setup_only=False, verbose=None):
+        if cleanup and os.path.exists(sim_path):
+            shutil.rmtree(sim_path)
+            os.mkdir(sim_path)
+        if not os.path.exists(sim_path):
+            os.mkdir(sim_path)
+        cwd = os.getcwd()
+        os.chdir(sim_path)
+        if verbose is not None:
+            self.thisptr.SetVerboseLevel(verbose)
+        assert os.getcwd() == sim_path
+        _openEMS.WelcomeScreen()
+        cdef int EC
+        EC = self.thisptr.SetupFDTD()
+        if EC!=0:
+            print('Run: Setup failed, error code: {}'.format(EC))
+        if setup_only or EC!=0:
+            return EC
+        self.thisptr.RunFDTD()

python/openEMS/physical_constants.py

@@ -0,0 +1,15 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec 18 20:58:01 2015
+
+@author: thorsten
+"""
+
+import numpy as np
+
+C0 = 299792458 # m/s
+MUE0 = 4e-7*np.pi # N/A^2
+EPS0 = 1/(MUE0*C0**2) # F/m
+
+# free space wave impedance
+Z0 = np.sqrt(MUE0/EPS0) # Ohm

python/openEMS/ports.py

@@ -0,0 +1,364 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Dec 17 22:53:39 2015
+
+@author: thorsten
+"""
+
+import os
+import numpy as np
+from CSXCAD.Utilities import CheckNyDir
+from openEMS import utilities
+
+from openEMS.physical_constants import *
+
+class UI_data:
+    def __init__(self, fns, path, freq, signal_type='pulse', **kw):
+        self.path = path
+        self.fns  = fns
+        self.freq = freq
+
+        self.ui_time = []
+        self.ui_val  = []
+        self.ui_f_val = []
+
+        for fn in fns:
+            tmp = np.loadtxt(os.path.join(path, fn),comments='%')
+            self.ui_time.append(tmp[:,0])
+            self.ui_val.append(tmp[:,1])
+            self.ui_f_val.append(utilities.DFT_time2freq(tmp[:,0], tmp[:,1], freq, signal_type=signal_type))
+
+# Port Base-Class
+class Port:
+    def __init__(self, CSX, port_nr, start, stop, excite, **kw):
+        self.CSX      = CSX
+        self.number   = port_nr
+        self.excite   = excite
+        self.start    = np.array(start, np.float)
+        self.stop     = np.array(stop, np.float)
+        self.Z_ref    = None
+        self.U_filenames = []
+        self.I_filenames = []
+
+        self.priority = 0
+        if 'priority' in kw:
+            self.priority = kw['priority']
+
+        self.prefix = ''
+        if 'PortNamePrefix' in kw:
+            self.prefix = kw['PortNamePrefix']
+        self.delay = 0
+
+        if 'delay' in kw:
+            self.delay = kw['delay']
+
+        self.lbl_temp = self.prefix + 'port_{}' +  '_{}'.format(self.number)
+
+    def ReadUIData(self, sim_path, freq, signal_type ='pulse'):
+        self.u_data = UI_data(self.U_filenames, sim_path, freq, signal_type )
+        self.uf_tot = 0
+        self.ut_tot = 0
+        for n in range(len(self.U_filenames)):
+            self.uf_tot += self.u_data.ui_f_val[n]
+            self.ut_tot += self.u_data.ui_val[n]
+
+        self.i_data = UI_data(self.I_filenames, sim_path, freq, signal_type )
+        self.if_tot = 0
+        self.it_tot = 0
+        for n in range(len(self.U_filenames)):
+            self.if_tot += self.i_data.ui_f_val[n]
+            self.it_tot += self.i_data.ui_val[n]
+
+
+    def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
+        self.ReadUIData(sim_path, freq, signal_type)
+
+        if ref_impedance is not None:
+            self.Z_ref = ref_impedance
+        assert self.Z_ref  is not None
+
+        if ref_plane_shift is not None:
+            assert hasattr(self, 'beta')
+            shift = ref_plane_shift
+            if self.measplane_shift:
+                shift -= self.measplane_shift
+            shift *= self.CSX.GetGrid().GetDeltaUnit()
+            phase = np.real(self.beta)*shift
+            uf_tot = self.uf_tot * np.cos(-phase) + 1j * self.if_tot * self.Z_ref * np.sin(-phase)
+            if_tot = self.if_tot * np.cos(-phase) + 1j * self.uf_tot / self.Z_ref * np.sin(-phase)
+            self.uf_tot = uf_tot
+            self.if_tot = if_tot
+
+        self.uf_inc = 0.5 * ( self.uf_tot + self.if_tot * self.Z_ref )
+        self.if_inc = 0.5 * ( self.if_tot + self.uf_tot / self.Z_ref )
+        self.uf_ref = self.uf_tot - self.uf_inc
+        self.if_ref = self.if_inc - self.if_tot
+
+        if type(self.Z_ref) == float:
+            self.ut_inc = 0.5 * ( self.ut_tot + self.it_tot * self.Z_ref )
+            self.it_inc = 0.5 * ( self.it_tot + self.ut_tot / self.Z_ref )
+            self.ut_ref = self.ut_tot - self.ut_inc
+            self.it_ref = self.it_inc - self.it_tot
+
+        # calc some more port parameter
+        # incoming power
+        self.P_inc = 0.5*np.real(self.uf_inc*np.conj(self.if_inc))
+        # reflected power
+        self.P_ref = 0.5*np.real(self.uf_ref*np.conj(self.if_ref))
+        # accepted power (incoming - reflected)
+        self.P_acc = 0.5*np.real(self.uf_tot*np.conj(self.if_tot))
+
+# Lumped-Port
+class LumpedPort(Port):
+    def __init__(self, CSX,  port_nr, R, start, stop, exc_dir, excite=0, **kw):
+        super(LumpedPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
+        self.R = R
+        self.exc_ny  = CheckNyDir(exc_dir)
+
+        self.direction = np.sign(self.stop[self.exc_ny]-self.start[self.exc_ny])
+        assert self.start[self.exc_ny]!=self.stop[self.exc_ny]
+
+        if self.R > 0:
+            lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.R)
+        elif self.R==0:
+            lumped_R = CSX.AddMetal(self.lbl_temp.format('resist'))
+
+        CSX.AddBox(lumped_R, self.start, self.stop, priority=self.priority)
+
+        if excite!=0:
+            exc_vec = np.zeros(3)
+            exc_vec[self.exc_ny] = -1*self.direction*excite
+            exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=exc_vec, delay=self.delay)
+            CSX.AddBox(exc, self.start, self.stop, priority=self.priority)
+
+        self.U_filenames = [self.lbl_temp.format('ut'), ]
+        u_start = 0.5*(self.start+self.stop)
+        u_start[self.exc_ny] = self.start[self.exc_ny]
+        u_stop  = 0.5*(self.start+self.stop)
+        u_stop[self.exc_ny]  = self.stop[self.exc_ny]
+        u_probe = CSX.AddProbe(self.U_filenames[0], p_type=0, weight=-1*self.direction)
+        CSX.AddBox(u_probe, u_start, u_stop)
+
+        self.I_filenames = [self.lbl_temp.format('it'), ]
+        i_start = np.array(self.start)
+        i_start[self.exc_ny] = 0.5*(self.start[self.exc_ny]+self.stop[self.exc_ny])
+        i_stop  = np.array(self.stop)
+        i_stop[self.exc_ny]  = 0.5*(self.start[self.exc_ny]+self.stop[self.exc_ny])
+        i_probe = CSX.AddProbe(self.I_filenames[0], p_type=1, weight=self.direction, norm_dir=self.exc_ny)
+        CSX.AddBox(i_probe, i_start, i_stop)
+
+    def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
+        if ref_impedance is None:
+            self.Z_ref = self.R
+        if ref_plane_shift is not None:
+            Warning('A lumped port does not support a reference plane shift! Ignoring...')
+        super(LumpedPort, self).CalcPort(sim_path, freq, ref_impedance, ref_plane_shift, signal_type)
+
+class MSLPort(Port):
+    def __init__(self, CSX, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw):
+        super(MSLPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
+        self.exc_ny  = CheckNyDir(exc_dir)
+        self.prop_ny = CheckNyDir(prop_dir)
+        self.direction   = np.sign(stop[self.prop_ny]-start[self.prop_ny])
+        self.upside_down = np.sign(stop[self.exc_ny]  -start[self.exc_ny])
+        assert (self.start!=self.stop).all()
+#        assert stop[self.prop_ny]!=start[self.prop_ny], 'port length in propergation direction may not be zero!'
+#        assert stop[self.exc_ny] !=start[self.exc_ny], 'port length in propergation direction may not be zero!'
+        assert self.exc_ny!=self.prop_ny
+
+        self.feed_shift = 0
+        if 'FeedShift' in kw:
+            self.feed_shift = kw['FeedShift']
+        self.measplane_shift = 0.5*np.abs(self.start[self.prop_ny]-self.stop[self.prop_ny])
+        if 'MeasPlaneShift' in kw:
+            self.measplane_shift =  kw['MeasPlaneShift']
+        self.measplane_pos = self.start[self.prop_ny] + self.measplane_shift*self.direction
+        self.feed_R = np.inf
+        if 'Feed_R' in kw:
+            self.feed_R = kw['Feed_R']
+
+        # add metal msl-plane
+        MSL_start = np.array(self.start)
+        MSL_stop  = np.array(self.stop)
+        MSL_stop[self.exc_ny] = MSL_start[self.exc_ny]
+        CSX.AddBox( metal_prop, MSL_start, MSL_stop, priority=self.priority )
+
+        mesh = CSX.GetGrid()
+        prop_lines = mesh.GetLines(self.prop_ny)
+        assert len(prop_lines)>5, 'At least 5 lines in propagation direction required!'
+        meas_pos_idx = np.argmin(np.abs(prop_lines-self.measplane_pos))
+        if meas_pos_idx==0:
+            meas_pos_idx=1
+        if meas_pos_idx>=len(prop_lines)-1:
+            meas_pos_idx=len(prop_lines)-2
+        self.measplane_shift = np.abs(self.start[self.prop_ny]-prop_lines[meas_pos_idx])
+        prope_idx = np.array([meas_pos_idx-1, meas_pos_idx, meas_pos_idx+1], np.int)
+        if self.direction<0:
+            prope_idx = np.flipud(prope_idx)
+        u_prope_pos = prop_lines[prope_idx]
+        self.U_filenames = []
+        self.U_delta = np.diff(u_prope_pos)
+        suffix = ['A', 'B', 'C']
+        for n in range(len(prope_idx)):
+            u_start = 0.5*(self.start+self.stop)
+            u_stop  = 0.5*(self.start+self.stop)
+            u_start[self.prop_ny] = u_prope_pos[n]
+            u_stop[self.prop_ny]  = u_prope_pos[n]
+            u_start[self.exc_ny]  = self.start[self.exc_ny]
+            u_stop[self.exc_ny]   = self.stop [self.exc_ny]
+            u_name = self.lbl_temp.format('ut') + suffix[n]
+            self.U_filenames.append(u_name)
+            u_probe = CSX.AddProbe(u_name, p_type=0, weight=self.upside_down)
+            CSX.AddBox(u_probe, u_start, u_stop)
+
+        i_prope_pos = u_prope_pos[0:2] + np.diff(u_prope_pos)/2.0
+        self.I_filenames = []
+        self.I_delta = np.diff(i_prope_pos)
+        i_start = np.array(self.start)
+        i_stop  = np.array(self.stop)
+        i_stop[self.exc_ny] = self.start[self.exc_ny]
+        for n in range(len(i_prope_pos)):
+            i_start[self.prop_ny] = i_prope_pos[n]
+            i_stop[self.prop_ny]  = i_prope_pos[n]
+            i_name = self.lbl_temp.format('it') + suffix[n]
+            self.I_filenames.append(i_name)
+            i_probe = CSX.AddProbe(i_name, p_type=1, weight=self.direction, norm_dir=self.prop_ny)
+            CSX.AddBox(i_probe, i_start, i_stop)
+
+        if excite!=0:
+            excide_pos_idx = np.argmin(np.abs(prop_lines-(self.start[self.prop_ny] + self.feed_shift*self.direction)))
+            exc_start = np.array(self.start)
+            exc_stop  = np.array(self.stop)
+            exc_start[self.prop_ny] = prop_lines[excide_pos_idx]
+            exc_stop [self.prop_ny] = prop_lines[excide_pos_idx]
+            exc_vec = np.zeros(3)
+            exc_vec[self.exc_ny] = -1*self.upside_down*excite
+            exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=exc_vec, delay=self.delay)
+            CSX.AddBox(exc, exc_start, exc_stop, priority=self.priority)
+
+        if self.feed_R>=0 and not np.isinf(self.feed_R):
+            R_start = np.array(self.start)
+            R_stop  = np.array(self.stop)
+            R_stop [self.prop_ny] = R_start[self.prop_ny]
+            if self.feed_R==0:
+                CSX.AddBox(metal_prop, R_start, R_stop)
+            else:
+                lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.feed_R)
+                CSX.AddBox(lumped_R, R_start, R_stop)
+
+    def ReadUIData(self, sim_path, freq, signal_type ='pulse'):
+        self.u_data = UI_data(self.U_filenames, sim_path, freq, signal_type )
+        self.uf_tot = self.u_data.ui_f_val[1]
+
+        self.i_data = UI_data(self.I_filenames, sim_path, freq, signal_type )
+        self.if_tot = 0.5*(self.i_data.ui_f_val[0]+self.i_data.ui_f_val[1])
+
+        unit = self.CSX.GetGrid().GetDeltaUnit()
+        Et = self.u_data.ui_f_val[1]
+        dEt = (self.u_data.ui_f_val[2] - self.u_data.ui_f_val[0]) / (np.sum(np.abs(self.U_delta)) * unit)
+        Ht = self.if_tot # space averaging: Ht is now defined at the same pos as Et
+        dHt = (self.i_data.ui_f_val[1] - self.i_data.ui_f_val[0]) / (np.abs(self.I_delta[0]) * unit)
+
+        beta = np.sqrt( - dEt * dHt / (Ht * Et) )
+        beta[np.real(beta) < 0] *= -1 # determine correct sign (unlike the paper)
+        self.beta = beta
+
+        # determine ZL
+        self.Z_ref = np.sqrt(Et * dEt / (Ht * dHt))
+
+class WaveguidePort(Port):
+    def __init__(self, CSX, port_nr, start, stop, exc_dir, mode_name, excite=0, **kw):
+        super(WaveguidePort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
+        self.exc_ny  = CheckNyDir(exc_dir)
+        self.ny_P  = (self.exc_ny+1)%3
+        self.ny_PP = (self.exc_ny+2)%3
+        self.direction = np.sign(stop[self.exc_ny]-start[self.exc_ny])
+        self.ref_index = 1
+
+        assert not (self.excite!=0 and stop[self.exc_ny]==start[self.exc_ny]), 'port length in excitation direction may not be zero if port is excited!'
+
+        self.InitMode(mode_name)
+
+        if excite!=0:
+            e_start = np.array(start)
+            e_stop  = np.array(stop)
+            e_stop[self.exc_ny] = e_start[self.exc_ny]
+            e_vec = np.ones(3)
+            e_vec[self.exc_ny]=0
+            exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=e_vec, delay=self.delay)
+            CSX.AddBox(exc, e_start, e_stop, priority=self.priority)
+
+        # voltage/current planes
+        m_start = np.array(start)
+        m_stop  = np.array(stop)
+        m_start[self.exc_ny] = m_stop[self.exc_ny]
+        self.measplane_shift = np.abs(stop[self.exc_ny] - start[self.exc_ny])
+
+        self.U_filenames = [self.lbl_temp.format('ut'), ]
+
+        u_probe = CSX.AddProbe(self.U_filenames[0], p_type=10, mode_function=self.E_func)
+        CSX.AddBox(u_probe, m_start, m_stop)
+
+        self.I_filenames = [self.lbl_temp.format('it'), ]
+        i_probe = CSX.AddProbe(self.I_filenames[0], p_type=11, weight=self.direction, mode_function=self.H_func)
+        CSX.AddBox(i_probe, m_start, m_stop)
+
+    def InitMode(self, wg_mode):
+        self.WG_mode = wg_mode
+        assert len(self.WG_mode)==4, 'Invalid mode definition'
+        self.unit = self.CSX.GetGrid().GetDeltaUnit()
+        if self.WG_mode.startswith('TE'):
+            self.TE = True
+            self.TM = False
+        else:
+            self.TE = False
+            self.TM = True
+        self.M = float(self.WG_mode[2])
+        self.N = float(self.WG_mode[3])
+        self.kc = None
+        self.E_func = [0,0,0]
+        self.H_func = [0,0,0]
+
+    def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
+        k = 2.0*np.pi*freq/C0*self.ref_index
+        self.beta = np.sqrt(k**2 - self.kc**2)
+        self.ZL = k * Z0 / self.beta    #analytic waveguide impedance
+        if ref_impedance is None:
+            self.Z_ref = self.ZL
+        super(WaveguidePort, self).CalcPort(sim_path, freq, ref_impedance, ref_plane_shift, signal_type)
+
+class RectWGPort(WaveguidePort):
+    def __init__(self, CSX, port_nr, start, stop, exc_dir, a, b, mode_name, excite=0, **kw):
+        self.WG_size = [a, b]
+        super(RectWGPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, exc_dir=exc_dir, mode_name=mode_name, excite=excite, **kw)
+
+    def InitMode(self, wg_mode):
+        super(RectWGPort, self).InitMode(wg_mode)
+        assert self.TE, 'Currently only TE-modes are supported! Mode found: {}'.format(self.WG_mode)
+
+        # values by David M. Pozar, Microwave Engineering, third edition
+        a = self.WG_size[0]
+        b = self.WG_size[1]
+        self.kc = np.sqrt((self.M*np.pi/a)**2 + (self.N*np.pi/b)**2)
+
+        xyz = 'xyz'
+        if self.start[self.ny_P]!=0:
+            name_P = '({}-{})'.format(xyz[self.ny_P], self.start[self.ny_P])
+        else:
+            name_P = xyz[self.ny_P]
+        if self.start[self.ny_PP]!=0:
+            name_PP = '({}-{})'.format(xyz[self.ny_P], self.start[self.ny_P])
+        else:
+            name_PP = xyz[self.ny_P]
+
+        a /= self.unit
+        b /= self.unit
+        if self.N>0:
+            self.E_func[self.ny_P]  = '{}*cos({}*{})*sin({}*{})'.format(self.N/b   , self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
+        if self.M>0:
+            self.E_func[self.ny_PP] = '{}*sin({}*{})*cos({}*{})'.format(-1*self.M/a, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
+
+        if self.M>0:
+            self.H_func[self.ny_P]  = '{}*sin({}*{})*cos({}*{})'.format(self.M/a, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
+        if self.N>0:
+            self.H_func[self.ny_PP] = '{}*cos({}*{})*sin({}*{})'.format(self.N/b, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)

python/openEMS/utilities.py

@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec 18 19:21:19 2015
+
+@author: thorsten
+"""
+
+import numpy as np
+
+def DFT_time2freq( t, val, freq, signal_type='pulse'):
+    assert len(t)==len(val)
+    assert len(freq)>0
+    f_val = np.zeros(len(freq))*1j
+    for n_f in range(len(freq)):
+        f_val[n_f] = np.sum( val*np.exp( -1j*2*np.pi*freq[n_f] * t ) )
+
+    if signal_type == 'pulse':
+        f_val *= t[1]-t[0]
+    elif signal_type == 'periodic':
+        f_val /= len(t)
+    else:
+        raise Exception('Unknown signal type: "{}"'.format(signal_type))
+
+    return 2*f_val  # single-sided spectrum
+
+def Check_Array_Equal(a,b, tol, relative=False):
+    a = np.array(a)
+    b = np.array(b)
+    if a.shape!=b.shape:
+        return False
+    if tol==0:
+        return (a==b).all()
+    if relative:
+        d = np.abs((a-b)/a)
+    else:
+        d = np.abs((a-b))
+    return np.max(d)<tol
+
+if __name__=="__main__":
+    import pylab as plt
+
+    t = np.linspace(0,2,201)
+
+    s = np.sin(2*np.pi*2*t)
+    plt.plot(t,s)
+
+    f = np.linspace(0,3,101)
+    sf = DFT_time2freq(t, s, f, 'periodic')
+
+    plt.figure()
+    plt.plot(f, np.abs(sf))
+
+    plt.show()
+
+

python/setup.py

@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Dec 13 23:48:22 2015
+
+@author: thorsten
+"""
+
+from distutils.core import setup
+from distutils.extension import Extension
+from Cython.Build import cythonize
+
+import os, sys
+ROOT_DIR = os.path.dirname(__file__)
+
+sys.path.append(os.path.join(ROOT_DIR,'..','..','CSXCAD','python'))
+
+extensions = [
+    Extension("*", [os.path.join(os.path.dirname(__file__), "openEMS","*.pyx")],
+        language="c++",             # generate C++ code
+        libraries    = ['CSXCAD','openEMS', 'nf2ff']),
+]
+
+setup(
+  name="openEMS",
+  version = '0.0.33',
+  description = "Python interface for the openEMS FDTD library",
+  classifiers = [
+      'Development Status :: 3 - Alpha',
+      'Intended Audience :: Developers',
+      'Intended Audience :: Information Technology',
+      'Intended Audience :: Science/Research',
+      'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
+      'Programming Language :: Python',
+      'Topic :: Scientific/Engineering',
+      'Topic :: Software Development :: Libraries :: Python Modules',
+      'Operating System :: POSIX :: Linux',
+      'Operating System :: Microsoft :: Windows',
+  ],
+  author = 'Thorsten Liebig',
+  author_email = 'Thorsten.Liebig@gmx.de',
+  maintainer = 'Thorsten Liebig',
+  maintainer_email = 'Thorsten.Liebig@gmx.de',
+  url = 'http://openEMS.de',
+  packages=["openEMS", ],
+  package_data={'openEMS': ['*.pxd']},
+  ext_modules = cythonize(extensions)
+ )