summaryrefslogtreecommitdiff
path: root/openEMS/python/Tutorials/RCS_Sphere.py
diff options
context:
space:
mode:
Diffstat (limited to 'openEMS/python/Tutorials/RCS_Sphere.py')
-rw-r--r--openEMS/python/Tutorials/RCS_Sphere.py126
1 files changed, 126 insertions, 0 deletions
diff --git a/openEMS/python/Tutorials/RCS_Sphere.py b/openEMS/python/Tutorials/RCS_Sphere.py
new file mode 100644
index 0000000..df801ce
--- /dev/null
+++ b/openEMS/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 <thorsten.liebig@gmx.de>
+"""
+
+### 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