From 4ebe163aeb5b91054b20ec900283664a8a49292f Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Sun, 28 Aug 2016 21:32:48 +0200 Subject: [PATCH] python: initial code for python interface using cython Signed-off-by: Thorsten Liebig --- .gitignore | 6 +- python/openEMS/__init__.py | 0 python/openEMS/_nf2ff.pxd | 38 +++ python/openEMS/_nf2ff.pyx | 45 ++++ python/openEMS/nf2ff.py | 139 ++++++++++ python/openEMS/openEMS.pxd | 49 ++++ python/openEMS/openEMS.pyx | 176 +++++++++++++ python/openEMS/physical_constants.py | 15 ++ python/openEMS/ports.py | 364 +++++++++++++++++++++++++++ python/openEMS/utilities.py | 55 ++++ python/setup.py | 47 ++++ 11 files changed, 933 insertions(+), 1 deletion(-) create mode 100644 python/openEMS/__init__.py create mode 100644 python/openEMS/_nf2ff.pxd create mode 100644 python/openEMS/_nf2ff.pyx create mode 100644 python/openEMS/nf2ff.py create mode 100644 python/openEMS/openEMS.pxd create mode 100644 python/openEMS/openEMS.pyx create mode 100644 python/openEMS/physical_constants.py create mode 100644 python/openEMS/ports.py create mode 100644 python/openEMS/utilities.py create mode 100644 python/setup.py diff --git a/.gitignore b/.gitignore index 8ccadff..a0fc011 100644 --- a/.gitignore +++ b/.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 diff --git a/python/openEMS/__init__.py b/python/openEMS/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/openEMS/_nf2ff.pxd b/python/openEMS/_nf2ff.pxd new file mode 100644 index 0000000..bb00c34 --- /dev/null +++ b/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 diff --git a/python/openEMS/_nf2ff.pyx b/python/openEMS/_nf2ff.pyx new file mode 100644 index 0000000..8fc5b09 --- /dev/null +++ b/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) diff --git a/python/openEMS/nf2ff.py b/python/openEMS/nf2ff.py new file mode 100644 index 0000000..a414740 --- /dev/null +++ b/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)) diff --git a/python/openEMS/openEMS.pxd b/python/openEMS/openEMS.pxd new file mode 100644 index 0000000..f47da82 --- /dev/null +++ b/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 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 diff --git a/python/openEMS/openEMS.pyx b/python/openEMS/openEMS.pyx new file mode 100644 index 0000000..1d4f37f --- /dev/null +++ b/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() diff --git a/python/openEMS/physical_constants.py b/python/openEMS/physical_constants.py new file mode 100644 index 0000000..b99a965 --- /dev/null +++ b/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 diff --git a/python/openEMS/ports.py b/python/openEMS/ports.py new file mode 100644 index 0000000..99ae1d5 --- /dev/null +++ b/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) diff --git a/python/openEMS/utilities.py b/python/openEMS/utilities.py new file mode 100644 index 0000000..ea2950f --- /dev/null +++ b/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)