diff --git a/.gitignore b/.gitignore index 8ccadff..151f2da 100644 --- a/.gitignore +++ b/.gitignore @@ -12,7 +12,6 @@ Makefile* *.pro.user* *.user *.orig -openEMS localPaths.pri .directory @@ -22,3 +21,11 @@ CMakeCache.txt cmake_install.cmake install_manifest.txt localConfig.cmake + +#python +*.pyc +*.pyo +python/**/*.cpp +!python/doc +python/doc/_build +python/doc/Tutorials/__* diff --git a/python/README.md b/python/README.md new file mode 100644 index 0000000..72f64be --- /dev/null +++ b/python/README.md @@ -0,0 +1,14 @@ +# openEMS python interface + +## Install +* Simple version: +```python +python setup.py install +``` + +* Extended options, e.g. for custom install path at */opt/openEMS*: +```python +python setup.py build_ext -I/opt/openEMS/include -L/opt/openEMS/lib -R/opt/openEMS/lib" +pyhton setup.py install +``` +**Note:** The install command may require root on Linux, or add --user to install to ~/.local diff --git a/python/Tutorials/Bent_Patch_Antenna.py b/python/Tutorials/Bent_Patch_Antenna.py new file mode 100644 index 0000000..ef2cb0e --- /dev/null +++ b/python/Tutorials/Bent_Patch_Antenna.py @@ -0,0 +1,198 @@ +# -*- coding: utf-8 -*- +""" + Bent Patch Antenna Tutorial + + Tested with + - python 3.4 + - openEMS v0.0.33+ + + (C) 2016 Thorsten Liebig + +""" + +### Import Libraries +import os, tempfile +from pylab import * +from mpl_toolkits.mplot3d import Axes3D + +from CSXCAD import CSXCAD + +from openEMS.openEMS import openEMS +from openEMS.physical_constants import * + + +### Setup the simulation +Sim_Path = os.path.join(tempfile.gettempdir(), 'Bent_Patch') + +post_proc_only = False + +unit = 1e-3 # all length in mm + +f0 = 2.4e9 # center frequency, frequency of interest! +lambda0 = round(C0/f0/unit) # wavelength in mm +fc = 0.5e9 # 20 dB corner frequency + +# patch width in alpha-direction +patch_width = 32 # resonant length in alpha-direction +patch_radius = 50 # radius +patch_length = 40 # patch length in z-direction + +#substrate setup +substrate_epsR = 3.38 +substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epsR +substrate_width = 80 +substrate_length = 90 +substrate_thickness = 1.524 +substrate_cells = 4 + +#setup feeding +feed_pos = -5.5 #feeding position in x-direction +feed_width = 2 #feeding port width +feed_R = 50 #feed resistance + +# size of the simulation box +SimBox_rad = 2*100 +SimBox_height = 1.5*200 + +### Setup FDTD parameter & excitation function +FDTD = openEMS(CoordSystem=1) # init a cylindrical FDTD +f0 = 2e9 # center frequency +fc = 1e9 # 20 dB corner frequency +FDTD.SetGaussExcite(f0, fc) +FDTD.SetBoundaryCond(['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR']) # boundary conditions + +### Setup the Geometry & Mesh +# init a cylindrical mesh +CSX = CSXCAD.ContinuousStructure(CoordSystem=1) +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) + +### Setup the geometry using cylindrical coordinates +# calculate some width as an angle in radiant +patch_ang_width = patch_width/(patch_radius+substrate_thickness) +substr_ang_width = substrate_width/patch_radius +feed_angle = feed_pos/patch_radius + +# create patch +patch = CSX.AddMetal('patch') # create a perfect electric conductor (PEC) +start = [patch_radius+substrate_thickness, -patch_ang_width/2, -patch_length/2 ] +stop = [patch_radius+substrate_thickness, patch_ang_width/2, patch_length/2 ] +CSX.AddBox(patch, priority=10, start=start, stop=stop, edges2grid='all') # add a box-primitive to the metal property 'patch' + +# create substrate +substrate = CSX.AddMaterial('substrate', epsilon=substrate_epsR, kappa=substrate_kappa ) +start = [patch_radius , -substr_ang_width/2, -substrate_length/2] +stop = [patch_radius+substrate_thickness, substr_ang_width/2, substrate_length/2] +substrate.AddBox(start=start, stop=stop, edges2grid='all') + +# save current density oon the patch +jt_patch = CSX.AddDump('Jt_patch', dump_type=3, file_type=1) +start = [patch_radius+substrate_thickness, -substr_ang_width/2, -substrate_length/2] +stop = [patch_radius+substrate_thickness, +substr_ang_width/2, substrate_length/2] +jt_patch.AddBox(start=start, stop=stop) + +# create ground +gnd = CSX.AddMetal('gnd') # create a perfect electric conductor (PEC) +start = [patch_radius, -substr_ang_width/2, -substrate_length/2] +stop = [patch_radius, +substr_ang_width/2, +substrate_length/2] +gnd.AddBox(priority=10, start=start, stop=stop, edges2grid='all') + +# apply the excitation & resist as a current source +start = [patch_radius , feed_angle, 0] +stop = [patch_radius+substrate_thickness, feed_angle, 0] +port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'r', 1.0, priority=50, edges2grid='all') + +### Finalize the Mesh +# add the simulation domain size +mesh.AddLine('r', patch_radius+np.array([-20, SimBox_rad])) +mesh.AddLine('a', [-0.75*pi, 0.75*pi]) +mesh.AddLine('z', [-SimBox_height/2, SimBox_height/2]) + +# add some lines for the substrate +mesh.AddLine('r', patch_radius+np.linspace(0,substrate_thickness,substrate_cells)) + +# generate a smooth mesh with max. cell size: lambda_min / 20 +max_res = C0 / (f0+fc) / unit / 20 +max_ang = max_res/(SimBox_rad+patch_radius) # max res in radiant +mesh.SmoothMeshLines(0, max_res, 1.4) +mesh.SmoothMeshLines(1, max_ang, 1.4) +mesh.SmoothMeshLines(2, max_res, 1.4) + +## Add the nf2ff recording box +nf2ff = FDTD.CreateNF2FFBox() + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'bent_patch.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + +### Postprocessing & plotting +f = np.linspace(max(1e9,f0-fc),f0+fc,401) +port.CalcPort(Sim_Path, f) +Zin = port.uf_tot / port.if_tot +s11 = port.uf_ref/port.uf_inc +s11_dB = 20.0*np.log10(np.abs(s11)) + +figure() +plot(f/1e9, s11_dB) +grid() +ylabel('s11 (dB)') +xlabel('frequency (GHz)') + +P_in = 0.5*np.real(port.uf_tot * np.conj(port.if_tot)) # antenna feed power + +# plot feed point impedance +figure() +plot( f/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' ) +grid() +plot( f/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' ) +title( 'feed point impedance' ) +xlabel( 'frequency (MHz)' ) +ylabel( 'impedance ($\Omega$)' ) +legend( ) + + +idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0] +if not len(idx)==1: + print('No resonance frequency found for far-field calulation') +else: + f_res = f[idx[0]] + theta = np.arange(-180.0, 180.0, 2.0) + print("Calculate NF2FF") + nf2ff_res_phi0 = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, 0, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xz.h5') + + figure(figsize=(15, 7)) + ax = subplot(121, polar=True) + E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + nf2ff_res_phi0.Dmax + ax.plot(np.deg2rad(theta), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xz-plane') + ax.grid(True) + ax.set_xlabel('theta (deg)') + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) + ax.legend(loc=3) + + phi = theta + nf2ff_res_theta90 = nf2ff.CalcNF2FF(Sim_Path, f_res, 90, phi, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xy.h5') + + ax = subplot(122, polar=True) + E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + nf2ff_res_theta90.Dmax + ax.plot(np.deg2rad(phi), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xy-plane') + ax.grid(True) + ax.set_xlabel('phi (deg)') + suptitle('Bent Patch Anteanna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14) + ax.legend(loc=3) + + print( 'radiated power: Prad = {:.2e} Watt'.format(nf2ff_res_theta90.Prad[0])) + print( 'directivity: Dmax = {:.1f} ({:.1f} dBi)'.format(nf2ff_res_theta90.Dmax[0], 10*np.log10(nf2ff_res_theta90.Dmax[0]))) + print( 'efficiency: nu_rad = {:.1f} %'.format(100*nf2ff_res_theta90.Prad[0]/real(P_in[idx[0]]))) + +show() + diff --git a/python/Tutorials/CRLH_Extraction.py b/python/Tutorials/CRLH_Extraction.py new file mode 100644 index 0000000..4e1850e --- /dev/null +++ b/python/Tutorials/CRLH_Extraction.py @@ -0,0 +1,239 @@ +# -*- coding: utf-8 -*- +""" + Tutorials / CRLH_Extraction + + Describtion at: + http://openems.de/index.php/Tutorial:_CRLH_Extraction + + Tested with + - python 3.4 + - openEMS v0.0.34+ + + (C) 2016 Thorsten Liebig +""" + + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import * + +### Class to represent single CRLH unit cells +class CRLH_Cells: + def __init__(self, LL, LW, Top, Bot, GLT, GLB, SL, SW, VR): + self.LL = LL # Line length + self.LW = LW # Line width + self.Top = Top # top signal height + self.Bot = Bot # bottom signal height + self.GLT = GLT # gap length top + self.GLB = GLB # gap length bottom + self.SL = SL # stub length + self.SW = SW # stub width + self.VR = VR # via radius + self.props = dict() # property dictionary + self.edge_resolution = None + + def createProperties(self, CSX): + for p in ['metal_top', 'metal_bot', 'via']: + self.props[p] = CSX.AddMetal(p) + + def setEdgeResolution(self, res): + self.edge_resolution = res + + def createCell(self, translate = [0,0,0]): + def append_mesh(mesh1, mesh2): + for n in range(3): + if mesh1[n] is None: + mesh1[n] = mesh2[n] + elif mesh2[n] is None: + continue + else: + mesh1[n] += mesh2[n] + return mesh1 + translate = array(translate) + start = [-self.LL/2 , -self.LW/2, self.Top] + translate + stop = [-self.GLT/2, self.LW/2, self.Top] + translate + box = self.props['metal_top'].AddBox(start, stop, priority=10) + mesh = box.GetGridHint('x', metal_edge_res=self.edge_resolution, down_dir=False) + append_mesh(mesh, box.GetGridHint('y', metal_edge_res=self.edge_resolution) ) + + start = [+self.LL/2 , -self.LW/2, self.Top] + translate + stop = [+self.GLT/2, self.LW/2, self.Top] + translate + box = self.props['metal_top'].AddBox(start, stop, priority=10) + append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution, up_dir=False) ) + + start = [-(self.LL-self.GLB)/2, -self.LW/2, self.Bot] + translate + stop = [+(self.LL-self.GLB)/2, self.LW/2, self.Bot] + translate + box = self.props['metal_bot'].AddBox(start, stop, priority=10) + append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution) ) + + start = [-self.SW/2, -self.LW/2-self.SL, self.Bot] + translate + stop = [+self.SW/2, self.LW/2+self.SL, self.Bot] + translate + box = self.props['metal_bot'].AddBox(start, stop, priority=10) + append_mesh(mesh, box.GetGridHint('xy', metal_edge_res=self.edge_resolution) ) + + start = [0, -self.LW/2-self.SL+self.SW/2, 0 ] + translate + stop = [0, -self.LW/2-self.SL+self.SW/2, self.Bot] + translate + + self.props['via'].AddCylinder(start, stop, radius=self.VR, priority=10) + + start[1] *= -1 + stop [1] *= -1 + self.props['via'].AddCylinder(start, stop, radius=self.VR, priority=10) + + return mesh + + +if __name__ == '__main__': + ### Setup the simulation + Sim_Path = os.path.join(tempfile.gettempdir(), 'CRLH_Extraction') + post_proc_only = False + + unit = 1e-6 # specify everything in um + + feed_length = 30000 + + substrate_thickness = [1524, 101 , 254 ] + substrate_epsr = [3.48, 3.48, 3.48] + + CRLH = CRLH_Cells(LL = 14e3, LW = 4e3, GLB = 1950, GLT = 4700, SL = 7800, SW = 1000, VR = 250 , \ + Top = sum(substrate_thickness), \ + Bot = sum(substrate_thickness[:-1])) + + # frequency range of interest + f_start = 0.8e9 + f_stop = 6e9 + + ### Setup FDTD parameters & excitation function + CSX = ContinuousStructure() + FDTD = openEMS(EndCriteria=1e-5) + FDTD.SetCSX(CSX) + mesh = CSX.GetGrid() + mesh.SetDeltaUnit(unit) + + CRLH.createProperties(CSX) + + FDTD.SetGaussExcite((f_start+f_stop)/2, (f_stop-f_start)/2 ) + BC = {'PML_8' 'PML_8' 'MUR' 'MUR' 'PEC' 'PML_8'} + FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'PML_8'] ) + + ### Setup a basic mesh and create the CRLH unit cell + resolution = C0/(f_stop*sqrt(max(substrate_epsr)))/unit /30 # resolution of lambda/30 + CRLH.setEdgeResolution(resolution/4) + + mesh.SetLines('x', [-feed_length-CRLH.LL/2, 0, feed_length+CRLH.LL/2]) + mesh.SetLines('y', [-30000, 0, 30000]) + + substratelines = cumsum(substrate_thickness) + mesh.SetLines('z', [0, 20000]) + mesh.AddLine('z', cumsum(substrate_thickness)) + mesh.AddLine('z', linspace(substratelines[-2],substratelines[-1],4)) + + # create the CRLH unit cell (will define additional fixed mesh lines) + mesh_hint = CRLH.createCell() + mesh.AddLine('x', mesh_hint[0]) + mesh.AddLine('y', mesh_hint[1]) + + # Smooth the given mesh + mesh.SmoothMeshLines('all', resolution, 1.2) + + ### Setup the substrate layer + substratelines = [0] + substratelines.tolist() + start, stop = mesh.GetSimArea() + + for n in range(len(substrate_thickness)): + sub = CSX.AddMaterial( 'substrate_{}'.format(n), epsilon=substrate_epsr[n] ) + start[2] = substratelines[n] + stop [2] = substratelines[n+1] + + sub.AddBox( start, stop ) + + ### Add the feeding MSL ports + pec = CSX.AddMetal( 'PEC' ) + port = [None, None] + x_lines = mesh.GetLines('x') + portstart = [ x_lines[0], -CRLH.LW/2, substratelines[-1]] + portstop = [ -CRLH.LL/2, CRLH.LW/2, 0] + port[0] = FDTD.AddMSLPort( 1, pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=feed_length/2, priority=10) + + + portstart = [ x_lines[-1], -CRLH.LW/2, substratelines[-1]] + portstop = [ +CRLH.LL/2 , CRLH.LW/2, 0] + port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=feed_length/2, priority=10) + + ### Run the simulation + if 1: # debugging only + CSX_file = os.path.join(Sim_Path, 'CRLH_Extraction.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + + if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + + ### Post-Processing + f = linspace( f_start, f_stop, 1601 ) + for p in port: + p.CalcPort( Sim_Path, f, ref_impedance = 50, ref_plane_shift = feed_length) + + # calculate and plot scattering parameter + s11 = port[0].uf_ref / port[0].uf_inc + s21 = port[1].uf_ref / port[0].uf_inc + + plot(f/1e9,20*log10(abs(s11)),'k-' , linewidth=2, label='$S_{11}$') + plot(f/1e9,20*log10(abs(s21)),'r--', linewidth=2, label='$S_{21}$') + grid() + legend(loc=3) + ylabel('S-Parameter (dB)') + xlabel('frequency (GHz)') + ylim([-40, 2]) + + ### Extract CRLH parameter form ABCD matrix + A = ((1+s11)*(1-s11) + s21*s21)/(2*s21) + C = ((1-s11)*(1-s11) - s21*s21)/(2*s21) / port[1].Z_ref + + Y = C + Z = 2*(A-1)/C + + iZ = imag(Z) + iY = imag(Y) + + fse = interp(0, iZ, f) + fsh = interp(0, iY, f) + + df = f[1]-f[0] + fse_idx = np.where(f>fse)[0][0] + fsh_idx = np.where(f>fsh)[0][0] + + LR = 0.5*(iZ[fse_idx]-iZ[fse_idx-1])/(2*pi*df) + CL = 1/(2*pi*fse)**2/LR + + CR = 0.5*(iY[fsh_idx]-iY[fsh_idx-1])/(2*pi*df) + LL = 1/(2*pi*fsh)**2/CR + + print(' Series tank: CL = {:.2f} pF, LR = {:.2f} nH -> f_se = {:.2f} GHz '.format(CL*1e12, LR*1e9, fse*1e-9)) + print(' Shunt tank: CR = {:.2f} pF, LL = {:.2f} nH -> f_sh = {:.2f} GHz '.format(CR*1e12, LL*1e9, fsh*1e-9)) + + ### Calculate analytical wave-number of an inf-array of cells + w = 2*pi*f + wse = 2*pi*fse + wsh = 2*pi*fsh + beta_calc = real(arccos(1-(w**2-wse**2)*(w**2-wsh**2)/(2*w**2/CR/LR))) + + # plot + figure() + beta = -angle(s21)/CRLH.LL/unit + plot(abs(beta)*CRLH.LL*unit/pi,f*1e-9,'k-', linewidth=2, label=r'$\beta_{CRLH,\ 1\ cell}$' ) + grid() + plot(beta_calc/pi,f*1e-9,'c--', linewidth=2, label=r'$\beta_{CRLH,\ \infty\ cells}$') + plot(real(port[1].beta)*CRLH.LL*unit/pi,f*1e-9,'g-', linewidth=2, label=r'$\beta_{MSL}$') + ylim([1, 6]) + xlabel(r'$|\beta| p / \pi$') + ylabel('frequency (GHz)') + legend(loc=2) + + show() \ No newline at end of file diff --git a/python/Tutorials/Helical_Antenna.py b/python/Tutorials/Helical_Antenna.py new file mode 100644 index 0000000..3211ec8 --- /dev/null +++ b/python/Tutorials/Helical_Antenna.py @@ -0,0 +1,191 @@ +# -*- coding: utf-8 -*- +""" + Helical Antenna Tutorial + + Tested with + - python 3.4 + - openEMS v0.0.33+ + + (C) 2015-2016 Thorsten Liebig + +""" + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import CSXCAD + +from openEMS import openEMS +from openEMS.physical_constants import * + + +### Setup the simulation +Sim_Path = os.path.join(tempfile.gettempdir(), 'Helical_Ant') +post_proc_only = False + +unit = 1e-3 # all length in mm + +f0 = 2.4e9 # center frequency, frequency of interest! +lambda0 = round(C0/f0/unit) # wavelength in mm +fc = 0.5e9 # 20 dB corner frequency + +Helix_radius = 20 # --> diameter is ~ lambda/pi +Helix_turns = 10 # --> expected gain is G ~ 4 * 10 = 40 (16dBi) +Helix_pitch = 30 # --> pitch is ~ lambda/4 +Helix_mesh_res = 3 + +gnd_radius = lambda0/2 + +# feeding +feed_heigth = 3 +feed_R = 120 #feed impedance + +# size of the simulation box +SimBox = array([1, 1, 1.5])*2.0*lambda0 + +### Setup FDTD parameter & excitation function +FDTD = openEMS(EndCriteria=1e-4) +FDTD.SetGaussExcite( f0, fc ) +FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'PML_8'] ) + +### Setup Geometry & Mesh +CSX = CSXCAD.ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) + +max_res = floor(C0 / (f0+fc) / unit / 20) # cell size: lambda/20 + +# create helix mesh +mesh.AddLine('x', [-Helix_radius, 0, Helix_radius]) +mesh.SmoothMeshLines('x', Helix_mesh_res) +# add the air-box +mesh.AddLine('x', [-SimBox[0]/2-gnd_radius, SimBox[0]/2+gnd_radius]) +# create a smooth mesh between specified fixed mesh lines +mesh.SmoothMeshLines('x', max_res, ratio=1.4) + +# copy x-mesh to y-direction +mesh.SetLines('y', mesh.GetLines('x')) + +# create helix mesh in z-direction +mesh.AddLine('z', [0, feed_heigth, Helix_turns*Helix_pitch+feed_heigth]) +mesh.SmoothMeshLines('z', Helix_mesh_res) + +# add the air-box +mesh.AddLine('z', [-SimBox[2]/2, max(mesh.GetLines('z'))+SimBox[2]/2 ]) +# create a smooth mesh between specified fixed mesh lines +mesh.SmoothMeshLines('z', max_res, ratio=1.4) + +### Create the Geometry +## * Create the metal helix using the wire primitive. +## * Create a metal gorund plane as cylinder. +# create a perfect electric conductor (PEC) +helix_metal = CSX.AddMetal('helix' ) + +ang = linspace(0,2*pi,21) +coil_x = Helix_radius*cos(ang) +coil_y = Helix_radius*sin(ang) +coil_z = ang/2/pi*Helix_pitch + +Helix_x=np.array([]) +Helix_y=np.array([]) +Helix_z=np.array([]) +zpos = feed_heigth +for n in range(Helix_turns-1): + Helix_x = r_[Helix_x, coil_x] + Helix_y = r_[Helix_y, coil_y] + Helix_z = r_[Helix_z ,coil_z+zpos] + zpos = zpos + Helix_pitch + +p = np.array([Helix_x, Helix_y, Helix_z]) +helix_metal.AddCurve(p) + +# create ground circular ground +gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC) + +# add a box using cylindrical coordinates +start = [0, 0, -0.1] +stop = [0, 0, 0.1] +gnd.AddCylinder(start, stop, radius=gnd_radius) + +# apply the excitation & resist as a current source +start = [Helix_radius, 0, 0] +stop = [Helix_radius, 0, feed_heigth] +port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'z', 1.0, priority=5) + +# nf2ff calc +nf2ff = FDTD.CreateNF2FFBox(opt_resolution=[lambda0/15]*3) + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'helix.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + +### Postprocessing & plotting +freq = linspace( f0-fc, f0+fc, 501 ) +port.CalcPort(Sim_Path, freq) + +Zin = port.uf_tot / port.if_tot +s11 = port.uf_ref / port.uf_inc + +## Plot the feed point impedance +figure() +plot( freq/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' ) +grid() +plot( freq/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' ) +title( 'feed point impedance' ) +xlabel( 'frequency (MHz)' ) +ylabel( 'impedance ($\Omega$)' ) +legend( ) + +## Plot reflection coefficient S11 +figure() +plot( freq/1e6, 20*log10(abs(s11)), 'k-', linewidth=2 ) +grid() +title( 'reflection coefficient $S_{11}$' ) +xlabel( 'frequency (MHz)' ) +ylabel( 'reflection coefficient $|S_{11}|$' ) + +### Create the NFFF contour +## * calculate the far field at phi=0 degrees and at phi=90 degrees +theta = arange(0.,180.,1.) +phi = arange(-180,180,2) +disp( 'calculating the 3D far field...' ) + +nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f0, theta, phi, read_cached=True, verbose=True ) + +Dmax_dB = 10*log10(nf2ff_res.Dmax[0]) +E_norm = 20.0*log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0]) + +theta_HPBW = theta[ np.where(squeeze(E_norm[:,phi==0]) + +""" + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import * + + +### Setup the simulation +Sim_Path = os.path.join(tempfile.gettempdir(), 'NotchFilter') +post_proc_only = False + +unit = 1e-6 # specify everything in um +MSL_length = 50000 +MSL_width = 600 +substrate_thickness = 254 +substrate_epr = 3.66 +stub_length = 12e3 +f_max = 7e9 + +### Setup FDTD parameters & excitation function +FDTD = openEMS() +FDTD.SetGaussExcite( f_max/2, f_max/2 ) +FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] ) + +### Setup Geometry & Mesh +CSX = ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) + +resolution = C0/(f_max*sqrt(substrate_epr))/unit/50 # resolution of lambda/50 +third_mesh = array([2*resolution/3, -resolution/3])/4 + +## Do manual meshing +mesh.AddLine('x', 0) +mesh.AddLine('x', MSL_width/2+third_mesh) +mesh.AddLine('x', -MSL_width/2-third_mesh) +mesh.SmoothMeshLines('x', resolution/4) + +mesh.AddLine('x', [-MSL_length, MSL_length]) +mesh.SmoothMeshLines('x', resolution) + +mesh.AddLine('y', 0) +mesh.AddLine('y', MSL_width/2+third_mesh) +mesh.AddLine('y', -MSL_width/2-third_mesh) +mesh.SmoothMeshLines('y', resolution/4) + +mesh.AddLine('y', [-15*MSL_width, 15*MSL_width+stub_length]) +mesh.AddLine('y', (MSL_width/2+stub_length)+third_mesh) +mesh.SmoothMeshLines('y', resolution) + +mesh.AddLine('z', linspace(0,substrate_thickness,5)) +mesh.AddLine('z', 3000) +mesh.SmoothMeshLines('z', resolution) + +## Add the substrate +substrate = CSX.AddMaterial( 'RO4350B', epsilon=substrate_epr) +start = [-MSL_length, -15*MSL_width, 0] +stop = [+MSL_length, +15*MSL_width+stub_length, substrate_thickness] +substrate.AddBox(start, stop ) + +## MSL port setup +port = [None, None] +pec = CSX.AddMetal( 'PEC' ) +portstart = [ -MSL_length, -MSL_width/2, substrate_thickness] +portstop = [ 0, MSL_width/2, 0] +port[0] = FDTD.AddMSLPort( 1, pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=MSL_length/3, priority=10) + +portstart = [MSL_length, -MSL_width/2, substrate_thickness] +portstop = [0 , MSL_width/2, 0] +port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=MSL_length/3, priority=10 ) + +## Filter-Stub Definition +start = [-MSL_width/2, MSL_width/2, substrate_thickness] +stop = [ MSL_width/2, MSL_width/2+stub_length, substrate_thickness] +pec.AddBox(start, stop, priority=10 ) + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'notch.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + +### Post-processing and plotting +f = linspace( 1e6, f_max, 1601 ) +for p in port: + p.CalcPort( Sim_Path, f, ref_impedance = 50) + +s11 = port[0].uf_ref / port[0].uf_inc +s21 = port[1].uf_ref / port[0].uf_inc + +plot(f/1e9,20*log10(abs(s11)),'k-',linewidth=2 , label='$S_{11}$') +grid() +plot(f/1e9,20*log10(abs(s21)),'r--',linewidth=2 , label='$S_{21}$') +legend() +ylabel('S-Parameter (dB)') +xlabel('frequency (GHz)') +ylim([-40, 2]) + +show() diff --git a/python/Tutorials/RCS_Sphere.py b/python/Tutorials/RCS_Sphere.py new file mode 100644 index 0000000..df801ce --- /dev/null +++ b/python/Tutorials/RCS_Sphere.py @@ -0,0 +1,126 @@ +# -*- coding: utf-8 -*- +""" + Tutorials / radar cross section of a metal sphere + + Tested with + - python 3.4 + - openEMS v0.0.34+ + + (C) 2016 Thorsten Liebig +""" + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import * +from openEMS.ports import UI_data + +### Setup the simulation +Sim_Path = os.path.join(tempfile.gettempdir(), 'RCS_Sphere') +post_proc_only = False + +unit = 1e-3 # all length in mm + +sphere_rad = 200 + +inc_angle = 0 #incident angle (to x-axis) in deg + +# size of the simulation box +SimBox = 1200 +PW_Box = 750 + +### Setup FDTD parameters & excitation function +FDTD = openEMS(EndCriteria=1e-5) + +f_start = 50e6 # start frequency +f_stop = 1000e6 # stop frequency +f0 = 500e6 +FDTD.SetGaussExcite( 0.5*(f_start+f_stop), 0.5*(f_stop-f_start) ) + +FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'PML_8', 'PML_8', 'PML_8', 'PML_8'] ) + +### Setup Geometry & Mesh +CSX = ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) + +#create mesh +mesh.SetLines('x', [-SimBox/2, 0, SimBox/2]) +mesh.SmoothMeshLines('x', C0 / f_stop / unit / 20) # cell size: lambda/20 +mesh.SetLines('y', mesh.GetLines('x')) +mesh.SetLines('z', mesh.GetLines('x')) + +### Create a metal sphere and plane wave source +sphere_metal = CSX.AddMetal( 'sphere' ) # create a perfect electric conductor (PEC) +sphere_metal.AddSphere(priority=10, center=[0, 0, 0], radius=sphere_rad) + +# plane wave excitation +k_dir = [cos(inc_angle), sin(inc_angle), 0] # plane wave direction +E_dir = [0, 0, 1] # plane wave polarization --> E_z + +pw_exc = CSX.AddExcitation('plane_wave', exc_type=10, exc_val=E_dir) +pw_exc.SetPropagationDir(k_dir) +pw_exc.SetFrequency(f0) + +start = np.array([-PW_Box/2, -PW_Box/2, -PW_Box/2]) +stop = -start +pw_exc.AddBox(start, stop) + +# nf2ff calc +nf2ff = FDTD.CreateNF2FFBox() + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'RCS_Sphere.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + +### Postprocessing & plotting +# get Gaussian pulse stength at frequency f0 +ef = UI_data('et', Sim_Path, freq=f0) + +Pin = 0.5*norm(E_dir)**2/Z0 * abs(ef.ui_f_val[0])**2 +# +nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f0, 90, arange(-180, 180.1, 2)) +RCS = 4*pi/Pin[0]*nf2ff_res.P_rad[0] + +fig = figure() +ax = fig.add_subplot(111, polar=True) +ax.plot( nf2ff_res.phi, RCS[0], 'k-', linewidth=2 ) +ax.grid(True) + +# calculate RCS over frequency +freq = linspace(f_start,f_stop,100) +ef = UI_data( 'et', Sim_Path, freq ) # time domain/freq domain voltage +Pin = 0.5*norm(E_dir)**2/Z0 * abs(np.array(ef.ui_f_val[0]))**2 + +nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, freq, 90, 180+inc_angle, outfile='back_nf2ff.h5') + +back_scat = np.array([4*pi/Pin[fn]*nf2ff_res.P_rad[fn][0][0] for fn in range(len(freq))]) + +figure() +plot(freq/1e6,back_scat, linewidth=2) +grid() +xlabel('frequency (MHz)') +ylabel('RCS ($m^2$)') +title('radar cross section') + +figure() +semilogy(sphere_rad*unit/C0*freq,back_scat/(pi*sphere_rad*unit*sphere_rad*unit), linewidth=2) +ylim([10^-2, 10^1]) +grid() +xlabel('sphere radius / wavelength') +ylabel('RCS / ($\pi a^2$)') +title('normalized radar cross section') + +show() \ No newline at end of file diff --git a/python/Tutorials/Rect_Waveguide.py b/python/Tutorials/Rect_Waveguide.py new file mode 100644 index 0000000..5d38115 --- /dev/null +++ b/python/Tutorials/Rect_Waveguide.py @@ -0,0 +1,125 @@ +# -*- coding: utf-8 -*- +""" + Rectangular Waveguide Tutorial + + Describtion at: + http://openems.de/doc/openEMS/Tutorials.html#rectangular-waveguide + + Tested with + - python 3.4 + - openEMS v0.0.34+ + + (C) 2015-2016 Thorsten Liebig + +""" + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import * + +### Setup the simulation +Sim_Path = os.path.join(tempfile.gettempdir(), 'Rect_WG') + +post_proc_only = False +unit = 1e-6; #drawing unit in um + +# waveguide dimensions +# WR42 +a = 10700; #waveguide width +b = 4300; #waveguide heigth +length = 50000; + +# frequency range of interest +f_start = 20e9; +f_0 = 24e9; +f_stop = 26e9; +lambda0 = C0/f_0/unit; + +#waveguide TE-mode definition +TE_mode = 'TE10'; + +#targeted mesh resolution +mesh_res = lambda0/30 + +### Setup FDTD parameter & excitation function +FDTD = openEMS(NrTS=1e4); +FDTD.SetGaussExcite(0.5*(f_start+f_stop),0.5*(f_stop-f_start)); + +# boundary conditions +FDTD.SetBoundaryCond([0, 0, 0, 0, 3, 3]); + +### Setup geometry & mesh +CSX = ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) + +mesh.AddLine('x', [0, a]) +mesh.AddLine('y', [0, b]) +mesh.AddLine('z', [0, length]) + +## Apply the waveguide port +ports = [] +start=[0, 0, 10*mesh_res]; +stop =[a, b, 15*mesh_res]; +mesh.AddLine('z', [start[2], stop[2]]) +ports.append(FDTD.AddRectWaveGuidePort( 0, start, stop, 'z', a*unit, b*unit, TE_mode, 1)) + +start=[0, 0, length-10*mesh_res]; +stop =[a, b, length-15*mesh_res]; +mesh.AddLine('z', [start[2], stop[2]]) +ports.append(FDTD.AddRectWaveGuidePort( 1, start, stop, 'z', a*unit, b*unit, TE_mode)) + +mesh.SmoothMeshLines('all', mesh_res, ratio=1.4) + +### Define dump box... +Et = CSX.AddDump('Et', file_type=0, sub_sampling=[2,2,2]) +start = [0, 0, 0]; +stop = [a, b, length]; +Et.AddBox(start, stop); + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'rect_wg.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + +### Postprocessing & plotting +freq = linspace(f_start,f_stop,201) +for port in ports: + port.CalcPort(Sim_Path, freq) + +s11 = ports[0].uf_ref / ports[0].uf_inc +s21 = ports[1].uf_ref / ports[0].uf_inc +ZL = ports[0].uf_tot / ports[0].if_tot +ZL_a = ports[0].ZL # analytic waveguide impedance + +## Plot s-parameter +figure() +plot(freq*1e-6,20*log10(abs(s11)),'k-',linewidth=2, label='$S_{11}$') +grid() +plot(freq*1e-6,20*log10(abs(s21)),'r--',linewidth=2, label='$S_{21}$') +legend(); +ylabel('S-Parameter (dB)') +xlabel(r'frequency (MHz) $\rightarrow$') + +## Compare analytic and numerical wave-impedance +figure() +plot(freq*1e-6,real(ZL), linewidth=2, label='$\Re\{Z_L\}$') +grid() +plot(freq*1e-6,imag(ZL),'r--', linewidth=2, label='$\Im\{Z_L\}$') +plot(freq*1e-6,ZL_a,'g-.',linewidth=2, label='$Z_{L, analytic}$') +ylabel('ZL $(\Omega)$') +xlabel(r'frequency (MHz) $\rightarrow$') +legend() + +show() diff --git a/python/Tutorials/Simple_Patch_Antenna.py b/python/Tutorials/Simple_Patch_Antenna.py new file mode 100644 index 0000000..cd80f78 --- /dev/null +++ b/python/Tutorials/Simple_Patch_Antenna.py @@ -0,0 +1,151 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Dec 18 20:56:53 2015 + +@author: thorsten +""" + +### Import Libraries +import os, tempfile +from pylab import * + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import * + +### General parameter setup +Sim_Path = os.path.join(tempfile.gettempdir(), 'Simp_Patch') + +post_proc_only = False + +# patch width (resonant length) in x-direction +patch_width = 32 # +# patch length in y-direction +patch_length = 40 + +#substrate setup +substrate_epsR = 3.38 +substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epsR +substrate_width = 60 +substrate_length = 60 +substrate_thickness = 1.524 +substrate_cells = 4 + +#setup feeding +feed_pos = -6 #feeding position in x-direction +feed_R = 50 #feed resistance + +# size of the simulation box +SimBox = np.array([200, 200, 150]) + +# setup FDTD parameter & excitation function +f0 = 2e9 # center frequency +fc = 1e9 # 20 dB corner frequency + +### FDTD setup +## * Limit the simulation to 30k timesteps +## * Define a reduced end criteria of -40dB +FDTD = openEMS(NrTS=30000, EndCriteria=1e-4) +FDTD.SetGaussExcite( f0, fc ) +FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR'] ) + + +CSX = ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(1e-3) +mesh_res = C0/(f0+fc)/1e-3/20 + +### Generate properties, primitives and mesh-grid +#initialize the mesh with the "air-box" dimensions +mesh.AddLine('x', [-SimBox[0]/2, SimBox[0]/2]) +mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2] ) +mesh.AddLine('z', [-SimBox[2]/3, SimBox[2]*2/3] ) + +# create patch +patch = CSX.AddMetal( 'patch' ) # create a perfect electric conductor (PEC) +start = [-patch_width/2, -patch_length/2, substrate_thickness] +stop = [ patch_width/2 , patch_length/2, substrate_thickness] +patch.AddBox(priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch' +FDTD.AddEdges2Grid(dirs='xy', properties=patch, metal_edge_res=mesh_res/2) + +# create substrate +substrate = CSX.AddMaterial( 'substrate', epsilon=substrate_epsR, kappa=substrate_kappa) +start = [-substrate_width/2, -substrate_length/2, 0] +stop = [ substrate_width/2, substrate_length/2, substrate_thickness] +substrate.AddBox( priority=0, start=start, stop=stop ) + +# add extra cells to discretize the substrate thickness +mesh.AddLine('z', linspace(0,substrate_thickness,substrate_cells+1)) + +# create ground (same size as substrate) +gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC) +start[2]=0 +stop[2] =0 +gnd.AddBox(start, stop, priority=10) + +FDTD.AddEdges2Grid(dirs='xy', properties=gnd) + +# apply the excitation & resist as a current source +start = [feed_pos, 0, 0] +stop = [feed_pos, 0, substrate_thickness] +port = FDTD.AddLumpedPort(1, feed_R, start, stop, 'z', 1.0, priority=5, edges2grid='xy') + +mesh.SmoothMeshLines('all', mesh_res, 1.4) + +# Add the nf2ff recording box +nf2ff = FDTD.CreateNF2FFBox() + +### Run the simulation +if 0: # debugging only + CSX_file = os.path.join(Sim_Path, 'simp_patch.xml') + if not os.path.exists(Sim_Path): + os.mkdir(Sim_Path) + CSX.Write2XML(CSX_file) + os.system(r'AppCSXCAD "{}"'.format(CSX_file)) + +if not post_proc_only: + FDTD.Run(Sim_Path, verbose=3, cleanup=True) + + +### Post-processing and plotting +f = np.linspace(max(1e9,f0-fc),f0+fc,401) +port.CalcPort(Sim_Path, f) +s11 = port.uf_ref/port.uf_inc +s11_dB = 20.0*np.log10(np.abs(s11)) +figure() +plot(f/1e9, s11_dB, 'k-', linewidth=2, label='$S_{11}$') +grid() +legend() +ylabel('S-Parameter (dB)') +xlabel('Frequency (GHz)') + +idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0] +if not len(idx)==1: + print('No resonance frequency found for far-field calulation') +else: + f_res = f[idx[0]] + theta = np.arange(-180.0, 180.0, 2.0) + phi = [0., 90.] + nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3]) + + figure() + E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + nf2ff_res.Dmax[0] + plot(theta, np.squeeze(E_norm[:,0]), 'k-', linewidth=2, label='xz-plane') + plot(theta, np.squeeze(E_norm[:,1]), 'r--', linewidth=2, label='yz-plane') + grid() + ylabel('Directivity (dBi)') + xlabel('Theta (deg)') + title('Frequency: {} GHz'.format(f_res/1e9)) + legend() + +Zin = port.uf_tot/port.if_tot +figure() +plot(f/1e9, np.real(Zin), 'k-', linewidth=2, label='$\Re\{Z_{in}\}$') +plot(f/1e9, np.imag(Zin), 'r--', linewidth=2, label='$\Im\{Z_{in}\}$') +grid() +legend() +ylabel('Zin (Ohm)') +xlabel('Frequency (GHz)') + +show() diff --git a/python/doc/Tutorials/Antenna_Tutorials.rst b/python/doc/Tutorials/Antenna_Tutorials.rst new file mode 100644 index 0000000..4a5927a --- /dev/null +++ b/python/doc/Tutorials/Antenna_Tutorials.rst @@ -0,0 +1,9 @@ +Antennas +-------- + +.. toctree:: + :maxdepth: 1 + + Simple_Patch_Antenna + Helical_Antenna + Bent_Patch_Antenna diff --git a/python/doc/Tutorials/Bent_Patch_Antenna.rst b/python/doc/Tutorials/Bent_Patch_Antenna.rst new file mode 100644 index 0000000..9582563 --- /dev/null +++ b/python/doc/Tutorials/Bent_Patch_Antenna.rst @@ -0,0 +1,35 @@ +Bent Patch Antenna +================== + +* Setup & Simulate a bent patch antenna using a cylindrical mesh + +Introduction +------------- +**This tutorial covers:** + +* Setup of a Bent Patch Antenna (see for comparison: :ref:`simple_patch_antenna`) +* setup of a *cylindrical FDTD mesh*. +* Calculate the S-Parameter and input impedance +* Calculate far-field pattern 2D/3D + + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__Bent_Patch_Antenna.txt + +Images +------------- +.. figure:: images/Bent_Patch.png + :width: 49% + :alt: alternate text + + 3D view of the Bent Patch Antenna (AppCSXCAD) + +.. figure:: images/Bent_Patch_Pattern.png + :width: 80% + :alt: Farfield pattern + + Farfield pattern on an xy- and xz-plane + diff --git a/python/doc/Tutorials/CRLH_Extraction.rst b/python/doc/Tutorials/CRLH_Extraction.rst new file mode 100644 index 0000000..156cf73 --- /dev/null +++ b/python/doc/Tutorials/CRLH_Extraction.rst @@ -0,0 +1,40 @@ +CRLH Parameter Extraction +========================= + +* Setup a composite-right/left-handed (CRLH) unit cell and extract the equivalent circuit parameter. + +Introduction +------------- +**This tutorial covers:** + +* Setup a feeding mircostrip line & port +* Apply an inhomogeneous mesh used for improved accuracy and simulation speed +* Use an internal clss to setup a CRLH unit cell +* Use the port voltages and currents to extract the unit cell equivalent circuit parameter + +.. figure:: images/CRLH_cell.png + :width: 80% + :alt: CRLH unit cell with feeding MSL. + + CRLH unit cell with feeding MSL. + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__CRLH_Extraction.txt + +Images +------------- + +.. figure:: images/CRLH_Spara.png + :width: 80% + :alt: CRLH cell S-parameter + + CRLH cell S-parameter + +.. figure:: images/CRLH_dispersion.png + :width: 80% + :alt: CRLH unit cell dispersion diagram + + CRLH unit cell dispersion diagram diff --git a/python/doc/Tutorials/Helical_Antenna.rst b/python/doc/Tutorials/Helical_Antenna.rst new file mode 100644 index 0000000..6a94081 --- /dev/null +++ b/python/doc/Tutorials/Helical_Antenna.rst @@ -0,0 +1,32 @@ +Helical Antenna +=============== + +Introduction +------------- +**This tutorial covers:** + +* setup of a helix using the wire primitive +* setup a lumped feeding port (R_in = 120 Ohms) +* adding a near-field to far-field (nf2ff) box using an efficient subsampling +* calculate the S-Parameter of the antenna +* calculate and plot the far-field pattern + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__Helical_Antenna.txt + +Images +------------- +.. figure:: images/Helix_Ant.png + :width: 49% + :alt: alternate text + + 3D view of the Helical Antenna (AppCSXCAD) + +.. figure:: images/Helix_Ant_Pattern.png + :width: 49% + :alt: alternate text + + Far-Field pattern showing a right-handed circular polarization. diff --git a/python/doc/Tutorials/Intro_Tutorials.rst b/python/doc/Tutorials/Intro_Tutorials.rst new file mode 100644 index 0000000..8126f33 --- /dev/null +++ b/python/doc/Tutorials/Intro_Tutorials.rst @@ -0,0 +1,10 @@ +.. _intro_tutorials: + +Introductional Tutorials +------------------------ + + +.. toctree:: + + Rect_Waveguide + RCS_Sphere diff --git a/python/doc/Tutorials/MSL_NotchFilter.rst b/python/doc/Tutorials/MSL_NotchFilter.rst new file mode 100644 index 0000000..c84883e --- /dev/null +++ b/python/doc/Tutorials/MSL_NotchFilter.rst @@ -0,0 +1,27 @@ +Microstrip Notch Filter +======================= + + * A straight MSL line with a open-ended stub to create a simple microwave filter. + +Introduction +------------- +**This tutorial covers:** + + +* Setup a mircostrip line (MSL) and MSL port +* Apply an inhomogeneous mesh used for improved accuracy and simulation speed +* Calculate the S-Parameter of the filter + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__MSL_NotchFilter.txt + +Images +------------- +.. figure:: images/Notch_Filter_SPara.png + :width: 49% + :alt: S-Parameter over Frequency + + S-Parameter over Frequency diff --git a/python/doc/Tutorials/MicroWave_Tutorials.rst b/python/doc/Tutorials/MicroWave_Tutorials.rst new file mode 100644 index 0000000..92fc0ce --- /dev/null +++ b/python/doc/Tutorials/MicroWave_Tutorials.rst @@ -0,0 +1,10 @@ +.. _microwave_tutorials: + +Micro Wave Tutorials +-------------------- + + +.. toctree:: + + MSL_NotchFilter + CRLH_Extraction diff --git a/python/doc/Tutorials/RCS_Sphere.rst b/python/doc/Tutorials/RCS_Sphere.rst new file mode 100644 index 0000000..1cbced8 --- /dev/null +++ b/python/doc/Tutorials/RCS_Sphere.rst @@ -0,0 +1,32 @@ +Metal Sphere Radar Cross Section +================================ + + * A 3D simulation demonstrating a the total-field/scattered-field approach on a metallic sphere with a RCS (radar cross section) calculation. + +Introduction +------------- +**This tutorial covers:** + + +* The total-field/scattered-field approach +* Calculation of a radar cross section (RCS) + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__RCS_Sphere.txt + +Images +------------- +.. figure:: images/RCS_pattern.png + :width: 49% + :alt: Radar cross section pattern + + Radar cross section pattern + +.. figure:: images/RCS_norm.png + :width: 49% + :alt: normalized radar cross section + + Normalized radar cross Section over normalized wavelength diff --git a/python/doc/Tutorials/Rect_Waveguide.rst b/python/doc/Tutorials/Rect_Waveguide.rst new file mode 100644 index 0000000..8b955ac --- /dev/null +++ b/python/doc/Tutorials/Rect_Waveguide.rst @@ -0,0 +1,27 @@ +Rectangular Waveguide +===================== + + * A simple rectangular waveguide, showing the openEMS mode profile capabilities. + +Introduction +------------- +**This tutorial covers:** + +* Setup a mode profile excitation +* Create voltage and current probes using the mode profile +* Calculate the waveguide impedance and S-Parameter + + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__Rect_Waveguide.txt + +Images +------------- +.. figure:: images/Rect_WG_SPara.png + :width: 49% + :alt: S-Parameter over Frequency + + S-Parameter over Frequency diff --git a/python/doc/Tutorials/Simple_Patch_Antenna.rst b/python/doc/Tutorials/Simple_Patch_Antenna.rst new file mode 100644 index 0000000..4345c56 --- /dev/null +++ b/python/doc/Tutorials/Simple_Patch_Antenna.rst @@ -0,0 +1,42 @@ +.. _simple_patch_antenna: + +Simple Patch Antenna +==================== + +Introduction +------------ +A simple patch antenna for 2.4 GHz. + +**This tutorial covers:** + +* Setup a patch, substrate and ground. +* Setup of a lumped feeding port. +* Adding a near-field to far-field (nf2ff) recording box. +* Calculate the S-Parameter of the antenna. +* Calculate and plot the far-field pattern + +Python Script +------------- +Get the latest version `from git `_. + +.. include:: ./__Simple_Patch_Antenna.txt + +Images +------ +.. figure:: images/Simp_Patch_S11.png + :width: 49% + :alt: S11 over Frequency + + S-Parameter over Frequency + +.. figure:: images/Simp_Patch_Zin.png + :width: 49% + :alt: Input Impedance + + Antenna Input Impedance + +.. figure:: images/Simp_Patch_Pattern.png + :width: 49% + :alt: Farfield pattern + + Farfield pattern for the xy- and yz-plane. diff --git a/python/doc/Tutorials/images/Bent_Patch.png b/python/doc/Tutorials/images/Bent_Patch.png new file mode 100644 index 0000000..0beddf9 Binary files /dev/null and b/python/doc/Tutorials/images/Bent_Patch.png differ diff --git a/python/doc/Tutorials/images/Bent_Patch_Pattern.png b/python/doc/Tutorials/images/Bent_Patch_Pattern.png new file mode 100644 index 0000000..4fc3211 Binary files /dev/null and b/python/doc/Tutorials/images/Bent_Patch_Pattern.png differ diff --git a/python/doc/Tutorials/images/Bent_Patch_SPara.png b/python/doc/Tutorials/images/Bent_Patch_SPara.png new file mode 100644 index 0000000..d3970fb Binary files /dev/null and b/python/doc/Tutorials/images/Bent_Patch_SPara.png differ diff --git a/python/doc/Tutorials/images/CRLH_Spara.png b/python/doc/Tutorials/images/CRLH_Spara.png new file mode 100644 index 0000000..d0f4a3b Binary files /dev/null and b/python/doc/Tutorials/images/CRLH_Spara.png differ diff --git a/python/doc/Tutorials/images/CRLH_cell.png b/python/doc/Tutorials/images/CRLH_cell.png new file mode 100644 index 0000000..85ae9b7 Binary files /dev/null and b/python/doc/Tutorials/images/CRLH_cell.png differ diff --git a/python/doc/Tutorials/images/CRLH_dispersion.png b/python/doc/Tutorials/images/CRLH_dispersion.png new file mode 100644 index 0000000..65fc163 Binary files /dev/null and b/python/doc/Tutorials/images/CRLH_dispersion.png differ diff --git a/python/doc/Tutorials/images/Helix_Ant.png b/python/doc/Tutorials/images/Helix_Ant.png new file mode 100644 index 0000000..16c532c Binary files /dev/null and b/python/doc/Tutorials/images/Helix_Ant.png differ diff --git a/python/doc/Tutorials/images/Helix_Ant_Pattern.png b/python/doc/Tutorials/images/Helix_Ant_Pattern.png new file mode 100644 index 0000000..020aae1 Binary files /dev/null and b/python/doc/Tutorials/images/Helix_Ant_Pattern.png differ diff --git a/python/doc/Tutorials/images/Notch_Filter_SPara.png b/python/doc/Tutorials/images/Notch_Filter_SPara.png new file mode 100644 index 0000000..9d1d7be Binary files /dev/null and b/python/doc/Tutorials/images/Notch_Filter_SPara.png differ diff --git a/python/doc/Tutorials/images/RCS_norm.png b/python/doc/Tutorials/images/RCS_norm.png new file mode 100644 index 0000000..57712ba Binary files /dev/null and b/python/doc/Tutorials/images/RCS_norm.png differ diff --git a/python/doc/Tutorials/images/RCS_pattern.png b/python/doc/Tutorials/images/RCS_pattern.png new file mode 100644 index 0000000..b3d9d23 Binary files /dev/null and b/python/doc/Tutorials/images/RCS_pattern.png differ diff --git a/python/doc/Tutorials/images/Rect_WG_SPara.png b/python/doc/Tutorials/images/Rect_WG_SPara.png new file mode 100644 index 0000000..0c0ff27 Binary files /dev/null and b/python/doc/Tutorials/images/Rect_WG_SPara.png differ diff --git a/python/doc/Tutorials/images/Simp_Patch_Pattern.png b/python/doc/Tutorials/images/Simp_Patch_Pattern.png new file mode 100644 index 0000000..2cd90d4 Binary files /dev/null and b/python/doc/Tutorials/images/Simp_Patch_Pattern.png differ diff --git a/python/doc/Tutorials/images/Simp_Patch_S11.png b/python/doc/Tutorials/images/Simp_Patch_S11.png new file mode 100644 index 0000000..381bb3c Binary files /dev/null and b/python/doc/Tutorials/images/Simp_Patch_S11.png differ diff --git a/python/doc/Tutorials/images/Simp_Patch_Zin.png b/python/doc/Tutorials/images/Simp_Patch_Zin.png new file mode 100644 index 0000000..bbabe54 Binary files /dev/null and b/python/doc/Tutorials/images/Simp_Patch_Zin.png differ diff --git a/python/doc/Tutorials/index.rst b/python/doc/Tutorials/index.rst new file mode 100644 index 0000000..fd5c5dc --- /dev/null +++ b/python/doc/Tutorials/index.rst @@ -0,0 +1,12 @@ +.. _tutorials: + +######### +Tutorials +######### + +.. toctree:: + :maxdepth: 2 + + Intro_Tutorials + MicroWave_Tutorials + Antenna_Tutorials diff --git a/python/doc/conf.py b/python/doc/conf.py new file mode 100644 index 0000000..b5c7b76 --- /dev/null +++ b/python/doc/conf.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# openEMS documentation build configuration file, created by +# sphinx-quickstart on Thu Sep 8 20:42:18 2016. +# +# This file is execfile()d with the current directory set to its +# containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +import sys +import os +import sphinx_rtd_theme + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +#sys.path.insert(0, os.path.abspath('.')) + +# -- General configuration ------------------------------------------------ + +# If your documentation needs a minimal Sphinx version, state it here. +#needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + 'sphinx.ext.intersphinx', + 'sphinx.ext.todo', + 'sphinx.ext.mathjax', + 'numpydoc', + 'sphinx.ext.autosummary', +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The encoding of source files. +#source_encoding = 'utf-8-sig' + +# The master toctree document. +master_doc = 'index' + +# General information about the project. +project = 'openEMS' +copyright = '2016, Thorsten Liebig' +author = 'Thorsten Liebig' + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +version = '0.0.34' +# The full version, including alpha/beta/rc tags. +release = '0.0.34' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# There are two options for replacing |today|: either, you set today to some +# non-false value, then it is used: +#today = '' +# Else, today_fmt is used as the format for a strftime call. +#today_fmt = '%B %d, %Y' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +exclude_patterns = ['_build'] + +# The reST default role (used for this markup: `text`) to use for all +# documents. +#default_role = None + +# If true, '()' will be appended to :func: etc. cross-reference text. +#add_function_parentheses = True + +# If true, the current module name will be prepended to all description +# unit titles (such as .. function::). +#add_module_names = True + +# If true, sectionauthor and moduleauthor directives will be shown in the +# output. They are ignored by default. +#show_authors = False + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# A list of ignored prefixes for module index sorting. +#modindex_common_prefix = [] + +# If true, keep warnings as "system message" paragraphs in the built documents. +#keep_warnings = False + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = True + + +# -- Options for HTML output ---------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +html_theme = 'sphinx_rtd_theme' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +#html_theme_options = {} + +# Add any paths that contain custom themes here, relative to this directory. +#html_theme_path = [] + +# The name for this set of Sphinx documents. If None, it defaults to +# " v documentation". +#html_title = None + +# A shorter title for the navigation bar. Default is the same as html_title. +#html_short_title = None + +# The name of an image file (relative to this directory) to place at the top +# of the sidebar. +#html_logo = None + +# The name of an image file (relative to this directory) to use as a favicon of +# the docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 +# pixels large. +#html_favicon = None + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + +# Add any extra paths that contain custom files (such as robots.txt or +# .htaccess) here, relative to this directory. These files are copied +# directly to the root of the documentation. +#html_extra_path = [] + +# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, +# using the given strftime format. +#html_last_updated_fmt = '%b %d, %Y' + +# If true, SmartyPants will be used to convert quotes and dashes to +# typographically correct entities. +#html_use_smartypants = True + +# Custom sidebar templates, maps document names to template names. +#html_sidebars = {} + +# Additional templates that should be rendered to pages, maps page names to +# template names. +#html_additional_pages = {} + +# If false, no module index is generated. +#html_domain_indices = True + +# If false, no index is generated. +#html_use_index = True + +# If true, the index is split into individual pages for each letter. +#html_split_index = False + +# If true, links to the reST sources are added to the pages. +#html_show_sourcelink = True + +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. +#html_show_sphinx = True + +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. +#html_show_copyright = True + +# If true, an OpenSearch description file will be output, and all pages will +# contain a tag referring to it. The value of this option must be the +# base URL from which the finished HTML is served. +#html_use_opensearch = '' + +# This is the file name suffix for HTML files (e.g. ".xhtml"). +#html_file_suffix = None + +# Language to be used for generating the HTML full-text search index. +# Sphinx supports the following languages: +# 'da', 'de', 'en', 'es', 'fi', 'fr', 'h', 'it', 'ja' +# 'nl', 'no', 'pt', 'ro', 'r', 'sv', 'tr' +#html_search_language = 'en' + +# A dictionary with options for the search language support, empty by default. +# Now only 'ja' uses this config value +#html_search_options = {'type': 'default'} + +# The name of a javascript file (relative to the configuration directory) that +# implements a search results scorer. If empty, the default will be used. +#html_search_scorer = 'scorer.js' + +# Output file base name for HTML help builder. +htmlhelp_basename = 'openEMSdoc' + +# -- Options for LaTeX output --------------------------------------------- + +latex_elements = { +# The paper size ('letterpaper' or 'a4paper'). +#'papersize': 'letterpaper', + +# The font size ('10pt', '11pt' or '12pt'). +#'pointsize': '10pt', + +# Additional stuff for the LaTeX preamble. +#'preamble': '', + +# Latex figure (float) alignment +#'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'openEMS.tex', 'openEMS Documentation', + 'Thorsten Liebig', 'manual'), +] + +# The name of an image file (relative to this directory) to place at the top of +# the title page. +#latex_logo = None + +# For "manual" documents, if this is true, then toplevel headings are parts, +# not chapters. +#latex_use_parts = False + +# If true, show page references after internal links. +#latex_show_pagerefs = False + +# If true, show URL addresses after external links. +#latex_show_urls = False + +# Documents to append as an appendix to all manuals. +#latex_appendices = [] + +# If false, no module index is generated. +#latex_domain_indices = True + + +# -- Options for manual page output --------------------------------------- + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'openems', 'openEMS Documentation', + [author], 1) +] + +# If true, show URL addresses after external links. +#man_show_urls = False + + +# -- Options for Texinfo output ------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'openEMS', 'openEMS Documentation', + author, 'openEMS', 'One line description of project.', + 'Miscellaneous'), +] + +# Documents to append as an appendix to all manuals. +#texinfo_appendices = [] + +# If false, no module index is generated. +#texinfo_domain_indices = True + +# How to display URL addresses: 'footnote', 'no', or 'inline'. +#texinfo_show_urls = 'footnote' + +# If true, do not generate a @detailmenu in the "Top" node's menu. +#texinfo_no_detailmenu = False + +numpydoc_show_class_members = False + +# Example configuration for intersphinx: refer to the Python standard library. +intersphinx_mapping = {'CSXCAD': ('http://openems.de/doc/CSXCAD/', None)} diff --git a/python/doc/convert_tutorials.py b/python/doc/convert_tutorials.py new file mode 100755 index 0000000..009937b --- /dev/null +++ b/python/doc/convert_tutorials.py @@ -0,0 +1,63 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- +""" +Created on Sat Sep 10 17:12:53 2016 + +@author: thorsten +""" + +import os +import glob + +DOC_DIR = os.path.dirname(__file__) +ROOT_DIR = os.path.join(DOC_DIR, '..') + +def main(): + in_path = os.path.join(ROOT_DIR, 'Tutorials') + + fns = glob.glob(os.path.join(in_path, '*.py')) + + for fn in fns: + bn = os.path.basename(fn) + out_fn = os.path.join(DOC_DIR, 'Tutorials', '__' + bn.replace('.py', '.txt')) + + in_code_block = False + in_ignore_block = False + out_fh = open(out_fn, 'w') + for line in open(fn, 'r'): + if in_ignore_block==False and line.startswith('"""'): + in_ignore_block = True + in_code_block = False + continue + elif in_ignore_block==True and line.startswith('"""'): + in_ignore_block = False + in_code_block = False + continue + elif in_ignore_block==True: + in_code_block = False + continue + elif line.startswith('# -*-'): + continue + elif not line.startswith('##'): + if not in_code_block: + if len(line.strip())==0: + continue + out_fh.write('\n.. code-block:: python\n\n') + in_code_block = True + out_fh.write(' ' + line) + elif line.startswith('###'): + if in_code_block: + out_fh.write('\n') + in_code_block = False + line = line.replace('#','').strip() + out_fh.write('**' + line + '**\n\n') +# out_fh.write('"'*len(line) + '\n') + elif line.startswith('##'): + if in_code_block: + out_fh.write('\n') + in_code_block = False + out_fh.write(line.replace('#','').strip() + '\n') + out_fh.close() + +if __name__ == '__main__': + main() diff --git a/python/doc/index.rst b/python/doc/index.rst new file mode 100644 index 0000000..f4cb9ea --- /dev/null +++ b/python/doc/index.rst @@ -0,0 +1,24 @@ +.. openEMS documentation master file, created by + sphinx-quickstart on Thu Sep 8 20:42:18 2016. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to openEMS's documentation! +=================================== + +Contents: + +.. toctree:: + :maxdepth: 3 + + Tutorials/index + openEMS_API + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + diff --git a/python/doc/nf2ff.rst b/python/doc/nf2ff.rst new file mode 100644 index 0000000..c4fbb01 --- /dev/null +++ b/python/doc/nf2ff.rst @@ -0,0 +1,16 @@ +.. _nf2ff: + +NF2FF +----- + +.. automodule:: openEMS.nf2ff + + NF2FF + ----- + .. autoclass:: nf2ff + :members: + + NF2FF Results + ----------------- + .. autoclass:: nf2ff_results + :members: diff --git a/python/doc/openEMS.rst b/python/doc/openEMS.rst new file mode 100644 index 0000000..c73dece --- /dev/null +++ b/python/doc/openEMS.rst @@ -0,0 +1,8 @@ +.. _openems: + +openEMS +------- + +.. automodule:: openEMS + :members: openEMS + :undoc-members: diff --git a/python/doc/openEMS_API.rst b/python/doc/openEMS_API.rst new file mode 100644 index 0000000..ae7e381 --- /dev/null +++ b/python/doc/openEMS_API.rst @@ -0,0 +1,11 @@ +.. _openems_api: + +openEMS Python Interface +======================== + +.. toctree:: + + openEMS + ports + nf2ff + diff --git a/python/doc/ports.rst b/python/doc/ports.rst new file mode 100644 index 0000000..79971e6 --- /dev/null +++ b/python/doc/ports.rst @@ -0,0 +1,37 @@ +.. _ports: + +Ports +----- + +.. automodule:: openEMS.ports + + Port (Base Class) + ----------------- + .. autoclass:: Port + :members: + :show-inheritance: + + Lumped Port + ----------- + .. autoclass:: LumpedPort + :members: + :show-inheritance: + + MSL Port + -------- + .. autoclass:: MSLPort + :members: + :show-inheritance: + + Waveguide Port + -------------- + .. autoclass:: WaveguidePort + :members: + :show-inheritance: + + Rect Waveguide Port + ------------------- + .. autoclass:: RectWGPort + :members: + :show-inheritance: + diff --git a/python/openEMS/__init__.py b/python/openEMS/__init__.py new file mode 100644 index 0000000..fb1ddb9 --- /dev/null +++ b/python/openEMS/__init__.py @@ -0,0 +1,4 @@ +# -*- coding: utf-8 -*- +# +# Shortcut openEMS import +from openEMS.openEMS import openEMS diff --git a/python/openEMS/_nf2ff.pxd b/python/openEMS/_nf2ff.pxd new file mode 100644 index 0000000..7e1b764 --- /dev/null +++ b/python/openEMS/_nf2ff.pxd @@ -0,0 +1,49 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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..8016694 --- /dev/null +++ b/python/openEMS/_nf2ff.pyx @@ -0,0 +1,59 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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 SetRadius(self, radius): + self.thisptr.SetRadius(radius) + + 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/automesh.py b/python/openEMS/automesh.py new file mode 100644 index 0000000..93472fa --- /dev/null +++ b/python/openEMS/automesh.py @@ -0,0 +1,77 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun Feb 19 20:29:25 2017 + +@author: thorsten +""" + +import sys +import numpy as np + +from CSXCAD import CSPrimitives +from CSXCAD.Utilities import CheckNyDir, GetMultiDirs + +def mesh_hint_from_primitive(primitive, dirs, **kw): + if primitive.GetType() is CSPrimitives.POINT: + return mesh_hint_from_point(primitive, dirs, **kw) + if primitive.GetType() is CSPrimitives.BOX: + return mesh_hint_from_box(primitive, dirs, **kw) + else: + return None + +def mesh_hint_from_point(point, dirs, **kw): + """ mesh_hint_from_point(point, dirs) + + Get a grid hint for the coordinates of the point. + + :param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all' + :returns: (3,) list of mesh hints + """ + hint = [None, None, None] + coord = point.GetCoord() + for ny in GetMultiDirs(dirs): + hint[ny] = [coord[ny],] + return hint + +def mesh_hint_from_box(box, dirs, **kw): + """ mesh_hint_from_box(box, dirs, metal_edge_res=None, **kw) + + Get a grid hint for the edges of the given box with an an optional 2D metal + edge resolution. + + :param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all' + :param metal_edge_res: float -- 2D flat edge resolution + :returns: (3,) list of mesh hints + """ + metal_edge_res = kw.get('metal_edge_res', None) + up_dir = kw.get('up_dir' , True) + down_dir = kw.get('down_dir', True) + + if metal_edge_res is None: + mer = 0 + else: + mer = np.array([-1.0, 2.0])/3 * metal_edge_res + if box.HasTransform(): + sys.stderr.write('FDTD::automesh: Warning, cannot add edges to grid with transformations enabled\n') + return + hint = [None, None, None] + start = np.fmin(box.GetStart(), box.GetStop()) + stop = np.fmax(box.GetStart(), box.GetStop()) + for ny in GetMultiDirs(dirs): + hint[ny] = [] + if metal_edge_res is not None and stop[ny]-start[ny]>metal_edge_res: + if down_dir: + hint[ny].append(start[ny]-mer[0]) + hint[ny].append(start[ny]-mer[1]) + if up_dir: + hint[ny].append(stop[ny]+mer[0]) + hint[ny].append(stop[ny]+mer[1]) + elif stop[ny]-start[ny]: + if down_dir: + hint[ny].append(start[ny]) + if up_dir: + hint[ny].append(stop[ny]) + else: + hint[ny].append(start[ny]) + return hint + diff --git a/python/openEMS/nf2ff.py b/python/openEMS/nf2ff.py new file mode 100644 index 0000000..b6d38d0 --- /dev/null +++ b/python/openEMS/nf2ff.py @@ -0,0 +1,210 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import os +import numpy as np +import h5py +from openEMS import _nf2ff +from openEMS import utilities + +class nf2ff: + """ + Create an nf2ff recording box. The nf2ff can either record in time-domain + or frequency-domain. Further more certain directions and boundary condition + mirroring can be enabled/disabled. + + :param name: str -- Name for this recording box. + :param start/stop: (3,) array -- Box start/stop coordinates. + :param directions: (6,) bool array -- Enable/Disables directions. + :param mirror: (6,) int array -- 0 (Off), 1 (PEC) or 2 (PMC) boundary mirroring + :param frequency: array like -- List of frequencies (FD-domain recording) + """ + 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 + + if np.isscalar(self.freq): + self.freq = [self.freq] + + 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] + self.e_dump.AddBox(l_start, l_stop) + self.h_dump.AddBox(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] + self.e_dump.AddBox(l_start, l_stop) + self.h_dump.AddBox(l_start, l_stop) + + def CalcNF2FF(self, sim_path, freq, theta, phi, radius=1, center=[0,0,0], outfile=None, read_cached=False, verbose=0): + """ CalcNF2FF(sim_path, freq, theta, phi, center=[0,0,0], outfile=None, read_cached=True, verbose=0): + + Calculate the far-field after the simulation is done. + + :param sim_path: str -- Simulation path + :param freq: array like -- list of frequency for transformation + :param theta/phi: array like -- Theta/Phi angles to calculate the far-field + :param radius: float -- Radius to calculate the far-field (default is 1m) + :param center: (3,) array -- phase center, must be inside the recording box + :param outfile: str -- File to save results in. (defaults to recording name) + :param read_cached: bool -- enable/disable read already existing results (default off) + :param verbose: int -- set verbose level (default 0) + + :returns: nf2ff_results class instance + """ + if np.isscalar(freq): + freq = [freq] + self.freq = freq + if np.isscalar(theta): + theta = [theta] + self.theta = theta + if np.isscalar(phi): + phi = [phi] + self.phi = phi + self.center = center + + if outfile is None: + fn = os.path.join(sim_path, self.name + '.h5') + else: + fn = os.path.join(sim_path, outfile) + if not read_cached or not os.path.exists(fn): + 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]) + + nfc.SetRadius(radius) + + 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) + + result = nf2ff_results(fn) + if result.phi is not None: + assert np.abs((result.r-radius)/radius)<1e-6, 'Radius does not match. Did you read an invalid chached result? Try "read_cached=False"' + assert utilities.Check_Array_Equal(np.rad2deg(result.theta), self.theta, 1e-4), 'Theta array does not match. Did you read an invalid chached result? Try "read_cached=False"' + assert utilities.Check_Array_Equal(np.rad2deg(result.phi), self.phi, 1e-4), 'Phi array does not match. Did you read an invalid chached result? Try "read_cached=False"' + assert utilities.Check_Array_Equal(result.freq, self.freq, 1e-6, relative=True), 'Frequency array does not match. Did you read an invalid chached result? Try "read_cached=False"' + return result + +class nf2ff_results: + """ + nf2ff result class containing all results obtained by the nf2ff calculation. + Usueally returned from nf2ff.CalcNF2FF + + Available attributes: + + * `fn` : file name + * `theta`: theta angles + * `phi` : phi angles + * `r` : radius + * `freq` : frequencies + * `Dmax` : directivity over frequency + * `Prad` : total radiated power over frequency + + * `E_theta` : theta component of electric field over frequency/theta/phi + * `E_phi` : phi component of electric field over frequency/theta/phi + * `E_norm` : abs component of electric field over frequency/theta/phi + * `E_cprh` : theta component of electric field over frequency/theta/phi + * `E_cplh` : theta component of electric field over frequency/theta/phi + * `P_rad` : radiated power (S) over frequency/theta/phi + """ + def __init__(self, fn): + self.fn = fn + h5_file = h5py.File(fn, 'r') + mesh_grp = h5_file['Mesh'] + self.phi = np.array(mesh_grp['phi']) + self.theta = np.array(mesh_grp['theta']) + self.r = np.array(mesh_grp['r']) + + data = h5_file['nf2ff'] + self.freq = np.array(data.attrs['Frequency']) + + self.Dmax = np.array(data.attrs['Dmax']) + self.Prad = np.array(data.attrs['Prad']) + + THETA, PHI = np.meshgrid(self.theta, self.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(self.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..c531dc1 --- /dev/null +++ b/python/openEMS/openEMS.pxd @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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(string 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..60abbbb --- /dev/null +++ b/python/openEMS/openEMS.pyx @@ -0,0 +1,440 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import os, sys, shutil +import numpy as np +cimport openEMS +from . import ports, nf2ff, automesh + +from CSXCAD.Utilities import GetMultiDirs + +cdef class openEMS: + """ openEMS + + This class is the main control class for the FDTD options and setup and + to run the final simulation. + + Examples + -------- + + >>> CSX = CSXCAD.ContinuousStructure() + >>> + >>> grid = CSX.GetGrid() + >>> grid.SetLines('x', np.arange(-50,50,1)) + >>> grid.SetLines('y', np.arange(-50,50,1)) + >>> grid.SetLines('z', np.arange(-2,2.1,1)) + >>> grid.SetDeltaUnit(1e-3) + >>> + >>> FDTD = openEMS(NrTS=1e4, EndCriteria=1e-4) + >>> + >>> FDTD.SetCSX(CSX) + >>> FDTD.SetBoundaryCond(['PML_8', 'PML_8', 'PML_8', 'PML_8', 'PEC', 'PEC']) + >>> FDTD.SetGaussExcite(0, 10e9) + >>> + >>> FDTD.AddLumpedPort(port_nr=1, R=50, start=[10, 0, -2], stop=[10, 0, 2], p_dir='z', excite=1) + >>> + >>> FDTD.Run(sim_path='/tmp/test') + + :param NrTS: max. number of timesteps to simulate (e.g. default=1e9) + :param EndCriteria: end criteria, e.g. 1e-5, simulations stops if energy has decayed by this value (<1e-4 is recommended, default=1e-5) + :param MaxTime: max. real time in seconds to simulate + :param OverSampling: nyquist oversampling of time domain dumps + :param CoordSystem: choose coordinate system (0 Cartesian, 1 Cylindrical) + :param MultiGrid: define a cylindrical sub-grid radius + :param TimeStep: force to use a given timestep (dangerous!) + :param TimeStepFactor: reduce the timestep by a given factor (>0 to <=1) + :param TimeStepMethod: 1 or 3 chose timestep method (1=CFL, 3=Rennigs (default)) + :param CellConstantMaterial: set to 1 to assume a material is constant inside a cell (material probing in cell center) + """ + @staticmethod + def WelcomeScreen(): + """ + Show the openEMS welcome screen. + """ + _openEMS.WelcomeScreen() + + def __cinit__(self, *args, **kw): + self.thisptr = new _openEMS() + self.__CSX = None + + if 'NrTS' in kw: + self.SetNumberOfTimeSteps(kw['NrTS']) + del kw['NrTS'] + else: + self.SetNumberOfTimeSteps(1e9) + if 'EndCriteria' in kw: + self.SetEndCriteria(kw['EndCriteria']) + del kw['EndCriteria'] + if 'MaxTime' in kw: + self.SetMaxTime(kw['MaxTime']) + del kw['MaxTime'] + if 'OverSampling' in kw: + self.SetOverSampling(kw['OverSampling']) + del kw['OverSampling'] + if 'CoordSystem' in kw: + self.SetCoordSystem(kw['CoordSystem']) + del kw['CoordSystem'] + if 'TimeStep' in kw: + self.SetTimeStep(kw['TimeStep']) + del kw['TimeStep'] + if 'TimeStepFactor' in kw: + self.SetTimeStepFactor(kw['TimeStepFactor']) + del kw['TimeStepFactor'] + if 'TimeStepMethod' in kw: + self.SetTimeStepMethod(kw['TimeStepMethod']) + del kw['TimeStepMethod'] + if 'CellConstantMaterial' in kw: + self.SetCellConstantMaterial(kw['CellConstantMaterial']) + del kw['CellConstantMaterial'] + if 'MultiGrid' in kw: + self.SetMultiGrid(kw['MultiGrid']) + del kw['MultiGrid'] + + assert len(kw)==0, 'Unknown keyword arguments: "{}"'.format(kw) + + def __dealloc__(self): + del self.thisptr + if self.__CSX is not None: + self.__CSX.thisptr = NULL + + def SetNumberOfTimeSteps(self, val): + """ SetNumberOfTimeSteps(val) + + Set the number of timesteps. E.g. 5e4 (default is 1e9) + """ + self.thisptr.SetNumberOfTimeSteps(val) + + def SetEndCriteria(self, val): + """ SetEndCriteria(val) + + Set the end critera value. E.g. 1e-6 for -60dB + """ + self.thisptr.SetEndCriteria(val) + + def SetOverSampling(self, val): + """ SetOverSampling(val) + + Set the time domain signal oversampling as multiple of the Nyquist-rate. + """ + self.thisptr.SetOverSampling(val) + + def SetCellConstantMaterial(self, val): + """ SetCellConstantMaterial(val) + + Set cell material averaging to assume constant material inside each primary cell. (Advanced option) + + :param val: bool -- Enable or Disable (default disabled) + """ + self.thisptr.SetCellConstantMaterial(val) + + def SetCoordSystem(self, val): + """ SetCoordSystem(val) + + Set the coordinate system. 0 --> Cartesian (default), 1 --> cylindrical + """ + assert (val==0 or val==1), 'SetCoordSystem: Invalid coordinate system' + if val==0: + self.thisptr.SetCylinderCoords(False) + elif val==1: + self.thisptr.SetCylinderCoords(True) + + def SetMultiGrid(self, radii): + """ SetMultiGrid(radii) + + Define radii at which a cylindrical multi grid should be defined. + + :param radii: array like, multigrid radii + + See Also + -------- + openEMS.SetCylinderCoords + """ + assert len(radii)>0, 'SetMultiGrid: invalid multi grid definition' + + grid_str = ','.join(['{}'.format(x) for x in radii]) + self.thisptr.SetupCylinderMultiGrid(grid_str.encode('UTF-8')) + + def SetCylinderCoords(self): + """ SetCylinderCoords() + + Enable use of cylindircal coordinates. + + See Also + -------- + openEMS.SetMultiGrid + """ + self.thisptr.SetCylinderCoords(True) + + def SetTimeStepMethod(self, val): + """ SetTimeStepMethod(val) + + Set the time step calculation method. (Advanced option) + + Options: + + * 1: CFL criteria + * 3: Advanced Rennings criteria (default) + + :param val: int -- 1 or 3 (See above) + """ + self.thisptr.SetTimeStepMethod(val) + + def SetTimeStep(self, val): + """ SetTimeStep(val) + + Set/force the timestep. (Advanced option) + + It is highly recommended to not use this method! You may use the + SetTimeStepFactor instead to reduce the time step if necessary! + """ + self.thisptr.SetTimeStep(val) + + def SetTimeStepFactor(self, val): + """ SetTimeStepFactor(val) + + Set a time step factor (>0..1) to increase FDTD stability. + + :param val: float -- >0..1 + """ + self.thisptr.SetTimeStepFactor(val) + + def SetMaxTime(self, val): + """ SetMaxTime(val) + + Set max simulation time for a max. number of timesteps. + """ + self.thisptr.SetMaxTime(val) + + def SetGaussExcite(self, f0, fc): + """ SetGaussExcite(f0, fc) + + Set a Gaussian pulse as excitation signal. + + :param f0: float -- Center frequency in Hz. + :param fc: float -- -20dB bandwidth in Hz. + """ + self.thisptr.SetGaussExcite(f0, fc) + + + def SetBoundaryCond(self, BC): + """ SetBoundaryCond(BC) + + Set the boundary conditions for all six FDTD directions. + + Options: + + * 0 or 'PEC' : perfect electric conductor (default) + * 1 or 'PMC' : perfect magnetic conductor, useful for symmetries + * 2 or 'MUR' : simple MUR absorbing boundary conditions + * 3 or 'PML-8' : PML absorbing boundary conditions + + :param BC: (8,) array or list -- see options above + """ + 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): + """ AddLumpedPort(port_nr, R, start, stop, p_dir, excite=0, **kw) + + Add a lumped port wit the given values and location. + + See Also + -------- + openEMS.ports.LumpedPort + """ + assert self.__CSX is not None, 'AddLumpedPort: CSX is not set!' + port = ports.LumpedPort(self.__CSX, port_nr, R, start, stop, p_dir, excite, **kw) + edges2grid = kw.get('edges2grid', None) + if edges2grid is not None: + grid = self.__CSX.GetGrid() + for n in GetMultiDirs(edges2grid): + grid.AddLine(n, start[n]) + if start[n] != stop[n]: + grid.AddLine(n, stop[n]) + return port + + def AddWaveGuidePort(self, port_nr, start, stop, p_dir, E_func, H_func, kc, excite=0, **kw): + """ AddWaveGuidePort(self, port_nr, start, stop, p_dir, E_func, H_func, kc, excite=0, **kw) + + Add a arbitrary waveguide port. + + See Also + -------- + openEMS.ports.WaveguidePort + """ + assert self.__CSX is not None, 'AddWaveGuidePort: CSX is not set!' + return ports.WaveguidePort(self.__CSX, port_nr, start, stop, p_dir, E_func, H_func, kc, excite, **kw) + + def AddRectWaveGuidePort(self, port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw): + """ AddRectWaveGuidePort(port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw) + + Add a rectilinear waveguide port. + + See Also + -------- + openEMS.ports.RectWGPort + """ + assert self.__CSX is not None, 'AddRectWaveGuidePort: CSX is not set!' + 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): + """ AddMSLPort(port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw) + + Add a microstrip transmission line port. + + See Also + -------- + openEMS.ports.MSLPort + """ + assert self.__CSX is not None, 'AddMSLPort: CSX is not set!' + 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): + """ CreateNF2FFBox(name='nf2ff', start=None, stop=None, **kw) + + Create a near-field to far-field box. + + This method will automatically adept the recording box to the current + FDTD grid and boundary conditions. + + Notes + ----- + * Make sure the mesh grid and all boundary conditions are finially defined. + + See Also + -------- + openEMS.nf2ff.nf2ff + """ + assert self.__CSX is not None, 'CreateNF2FFBox: CSX is not set!' + 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): + """ SetCSX(CSX) + + Set the CSXCAD Continuous Structure for CAD data handling. + + See Also + -------- + CSXCAD.ContinuousStructure + """ + self.__CSX = CSX + self.thisptr.SetCSX(CSX.thisptr) + + def GetCSX(self): + return self.__CSX + + def AddEdges2Grid(self, dirs, primitives=None, properties=None, **kw): + """ AddEdges2Grid(primitives, dirs, **kw) + + Add the edges of the given primitives to the FDTD grid. + + :param dirs: primitives -- one or more primitives + :param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all' + """ + csx = self.GetCSX() + if csx is None: + raise Exception('AddEdges2Grid: Unable to access CSX!') + prim_list = [] + if primitives is not None and type(primitives) is not list: + prim_list.append(primitives) + elif primitives is not None: + prim_list += primitives + + if properties is not None and type(properties) is not list: + prim_list += properties.GetAllPrimitives() + elif primitives is not None: + for prop in properties: + prim_list += prop.GetAllPrimitives() + + grid = csx.GetGrid() + for prim in prim_list: + hint = automesh.mesh_hint_from_primitive(prim, dirs, **kw) + if hint is None: + continue + for n in range(3): + if hint[n] is None: + continue + grid.AddLine(n, hint[n]) + + def Run(self, sim_path, cleanup=False, setup_only=False, verbose=None): + """ Run(sim_path, cleanup=False, setup_only=False, verbose=None) + + Run the openEMS FDTD simulation. + + :param sim_path: str -- path to run in and create result data + :param cleanup: bool -- remove exisiting sim_path to cleanup old results + :param setup_only: bool -- only perform FDTD setup, do not run simulation + :param verbose: int -- set the openEMS verbosity level 0..3 + """ + 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..8f014b5 --- /dev/null +++ b/python/openEMS/physical_constants.py @@ -0,0 +1,26 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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..1f0ea2a --- /dev/null +++ b/python/openEMS/ports.py @@ -0,0 +1,433 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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 + if type(fns)==str: + fns = [fns] + self.fns = fns + + if np.isscalar(freq): + freq = [freq] + 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: + """ + The port base class. + + :param CSX: Continuous Structure + :param port_nr: int -- port number + :param R: float -- port reference impedance, e.g. 50 (Ohms) + :param start, stop: (3,) array -- Start/Stop box coordinates + :param p_dir: int -- port direction + :param excite: float -- port excitation amplitude + :param priority: int -- priority of all contained primtives + :param PortNamePrefix: str -- a prefix for all ports-names + :param delay: float -- a positiv delay value to e.g. emulate a phase shift + """ + 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)) + +class LumpedPort(Port): + """ + The lumped port. + + See Also + -------- + 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], 'LumpedPort: start and stop may not be identical in excitation direction' + + 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')) + + lumped_R.AddBox(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) + exc.AddBox(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) + u_probe.AddBox(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) + i_probe.AddBox(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): + """ + The microstrip transmission line port. + + :param prop_dir: int/str -- direction of propagation + + See Also + -------- + 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] + metal_prop.AddBox(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) + u_probe.AddBox(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) + i_probe.AddBox(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) + exc.AddBox(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: + metal_prop.AddBox(R_start, R_stop) + else: + lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.feed_R) + lumped_R.AddBox(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): + """ + Base class for any waveguide port. + + See Also + -------- + Port, RectWGPort + + """ + def __init__(self, CSX, port_nr, start, stop, exc_dir, E_WG_func, H_WG_func, kc, 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.kc = kc + self.E_func = E_WG_func + self.H_func = H_WG_func + + 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) + exc.SetWeightFunction([str(x) for x in self.E_func]) + exc.AddBox(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) + u_probe.AddBox(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) + i_probe.AddBox(m_start, m_stop) + + + 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): + """ + Rectangular waveguide port. + + :param a,b: float -- Width/Height of rectangular waveguide port + + See Also + -------- + Port, WaveguidePort + + """ + def __init__(self, CSX, port_nr, start, stop, exc_dir, a, b, mode_name, excite=0, **kw): + Port.__init__(self, CSX, port_nr, start, stop, excite=0, **kw) + self.exc_ny = CheckNyDir(exc_dir) + self.ny_P = (self.exc_ny+1)%3 + self.ny_PP = (self.exc_ny+2)%3 + self.WG_size = [a, b] + + self.WG_mode = mode_name + 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]) + + 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] + + 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] + + kc = np.sqrt((self.M*np.pi/a)**2 + (self.N*np.pi/b)**2) + + a /= self.unit + b /= self.unit + E_func = [0,0,0] + H_func = [0,0,0] + if self.N>0: + 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: + 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: + 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: + 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) + + super(RectWGPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, exc_dir=exc_dir, E_WG_func=E_func, H_WG_func=H_func, kc=kc, excite=excite, **kw) + diff --git a/python/openEMS/utilities.py b/python/openEMS/utilities.py new file mode 100644 index 0000000..22bfbf6 --- /dev/null +++ b/python/openEMS/utilities.py @@ -0,0 +1,66 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de) +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +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)