/
scatter.py
executable file
·123 lines (106 loc) · 7.57 KB
/
scatter.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
#!/usr/bin/env python
#-*- coding: utf-8 -*-
""" An example python-meep simulation of a dielectric sphere scattering a broadband impulse,
illustrating the use of the convenient functions provided by meep_utils.py
(c) 2014 Filip Dominec, see http://f.dominec.eu/meep for more information
The simulated models of structures are stored in the `metamaterial_models.py' module; see its code for the
description of each structure and command-line parameters accepted. Several of these parameters are not
passed to the model; see the definition of `process_param()' in `meep_utils.py'.
"""
import time, sys, os
import numpy as np
from scipy.constants import c, epsilon_0, mu_0
import meep_utils, meep_materials, metamaterial_models
from meep_utils import in_sphere, in_xcyl, in_ycyl, in_zcyl, in_xslab, in_yslab, in_zslab, in_ellipsoid
import meep_mpi as meep
#import meep
# Model selection
model_param = meep_utils.process_param(sys.argv[1:])
model = metamaterial_models.models[model_param.get('model', 'default')](**model_param)
## Initialize volume, structure and the fields according to the model
vol = meep.vol3d(model.size_x, model.size_y, model.size_z, 1./model.resolution)
vol.center_origin()
s = meep_utils.init_structure(model=model, volume=vol, pml_axes=meep.Z)
f = meep.fields(s)
# Define the Bloch-periodic boundaries (any transversal component of k-vector is allowed)
f.use_bloch(meep.X, getattr(model, 'Kx', 0) / (-2*np.pi))
f.use_bloch(meep.Y, getattr(model, 'Ky', 0) / (-2*np.pi))
# Add the field source (see meep_utils for an example of how an arbitrary source waveform is defined)
if not getattr(model, 'frequency', None): ## (temporal source shape)
#src_time_type = meep.band_src_time(model.src_freq/c, model.src_width/c, model.simtime*c/1.1)
src_time_type = meep.gaussian_src_time(model.src_freq/c, model.src_width/c)
else:
src_time_type = meep.continuous_src_time(getattr(model, 'frequency', None)/c)
srcvolume = meep.volume( ## (spatial source shape)
meep.vec(-model.size_x/2, -model.size_y/2, -model.size_z/2+model.pml_thickness),
meep.vec( model.size_x/2, model.size_y/2, -model.size_z/2+model.pml_thickness))
#f.add_volume_source(meep.Ex, src_time_type, srcvolume)
## Replace the f.add_volume_source(meep.Ex, srctype, srcvolume) line with following:
## Option for a custom source (e.g. exciting some waveguide mode)
class SrcAmplitudeFactor(meep.Callback):
## The source amplitude is complex -> phase factor modifies its direction
## todo: implement in MEEP: we should define an AmplitudeVolume() object and reuse it for monitors later
def __init__(self, Kx=0, Ky=0):
meep.Callback.__init__(self)
(self.Kx, self.Ky) = Kx, Ky
def complex_vec(self, vec): ## Note: the 'vec' coordinates are _relative_ to the source center
# (oblique) plane wave source:
return np.exp(-1j*(self.Kx*vec.x() + self.Ky*vec.y()))
# (oblique) Gaussian beam source:
#return np.exp(-1j*(self.Kx*vec.x() + self.Ky*vec.y()) - (vec.x()/100e-6)**2 - (vec.y()/100e-6)**2)
af = SrcAmplitudeFactor(Kx=getattr(model, 'Kx', 0), Ky=getattr(model, 'Ky', 0))
meep.set_AMPL_Callback(af.__disown__())
#f.add_volume_source(meep.Ex, src_time_type, srcvolume, meep.AMPL)
f.add_volume_source(meep.Ex, src_time_type, srcvolume, meep.AMPL)
## Define monitor planes, and the field output for visualisation (controlled by keywords in the 'comment' parameter)
monitor_options = {'size_x':model.size_x, 'size_y':model.size_y, 'resolution':model.resolution, 'Kx':getattr(model, 'Kx', 0), 'Ky':getattr(model, 'Ky', 0)}
monitor1_Ex = meep_utils.AmplitudeMonitorPlane(f, comp=meep.Ex, z_position=model.monitor_z1, **monitor_options)
monitor1_Hy = meep_utils.AmplitudeMonitorPlane(f, comp=meep.Hy, z_position=model.monitor_z1, **monitor_options)
monitor2_Ex = meep_utils.AmplitudeMonitorPlane(f, comp=meep.Ex, z_position=model.monitor_z2, **monitor_options)
monitor2_Hy = meep_utils.AmplitudeMonitorPlane(f, comp=meep.Hy, z_position=model.monitor_z2, **monitor_options)
slices = []
#if not "noepssnapshot" in str(model.comment):
#slices += [meep_utils.Slice(model=model, field=f, components=(meep.Dielectric), at_t=0, name='EPS')]
if "narrowfreq-snapshots" in str(model.comment):
slices += [meep_utils.Slice(model=model, field=f, components=meep.Ex, at_y=0, at_t=np.inf,
name=('At%.3eHz'%getattr(model, 'frequency', None)) if getattr(model, 'frequency', None) else '',
outputpng=True, outputvtk=False)]
if "fieldevolution" in str(model.comment):
slices += [meep_utils.Slice(model=model, field=f, components=(meep.Ex), at_y=0, name='FieldEvolution',
min_timestep=.1/model.src_freq, outputgif=True, outputhdf=True, outputvtk=True)]
if "snapshote" in str(model.comment):
slices += [meep_utils.Slice(model=model, field=f, components=(meep.Ex, meep.Ey, meep.Ez), at_t=np.inf, name='SnapshotE')]
## Run the FDTD simulation or the frequency-domain solver
if not getattr(model, 'frequency', None): ## time-domain computation
f.step(); timer = meep_utils.Timer(simtime=model.simtime); meep.quiet(True) # use custom progress messages
while (f.time()/c < model.simtime): # timestepping cycle
f.step()
timer.print_progress(f.time()/c)
for monitor in (monitor1_Ex, monitor1_Hy, monitor2_Ex, monitor2_Hy): monitor.record(field=f)
for slice_ in slices: slice_.poll(f.time()/c)
for slice_ in slices: slice_.finalize()
meep_utils.notify(model.simulation_name, run_time=timer.get_time())
else: ## frequency-domain computation
f.solve_cw(getattr(model, 'MaxTol',0.001), getattr(model, 'MaxIter', 5000), getattr(model, 'BiCGStab', 8))
for monitor in (monitor1_Ex, monitor1_Hy, monitor2_Ex, monitor2_Hy): monitor.record(field=f)
for slice_ in slices: slice_.finalize()
meep_utils.notify(model.simulation_name)
## Get the reflection and transmission of the structure
if meep.my_rank() == 0:
#t = monitor1_Ex.get_time()
#Ex1, Hy1, Ex2, Hy2 = [mon.get_field_waveform() for mon in (monitor1_Ex, monitor1_Hy, monitor2_Ex, monitor2_Hy)]
freq, s11, s12, columnheaderstring = meep_utils.get_s_parameters(monitor1_Ex, monitor1_Hy, monitor2_Ex, monitor2_Hy,
frequency_domain=True if getattr(model, 'frequency', None) else False,
frequency=getattr(model, 'frequency', None), ## procedure compatible with both FDTD and FDFD
intf=getattr(model, 'interesting_frequencies', [0, model.src_freq+model.src_width]), ## clip the frequency range for plotting
pad_zeros=1.0, ## speed-up FFT, and stabilize eff-param retrieval
Kx=getattr(model, 'Kx', 0), Ky=getattr(model, 'Ky', 0), ## enable oblique incidence (works only if monitors in vacuum)
eps1=getattr(model, 'mon1eps', 1), eps2=getattr(model, 'mon2eps', 1)) ## enable monitors inside dielectrics
if len(s11)>0:
meep_utils.savetxt(fname=model.simulation_name+".dat", fmt="%.6e",
## Save 5 columns: freq, amplitude/phase for reflection/transmission:
X=zip(freq, np.abs(s11), np.angle(s11), np.abs(s12), np.angle(s12)),
header=model.parameterstring+columnheaderstring) ## Export header
with open("./last_simulation_name.dat", "w") as outfile: outfile.write(model.simulation_name)
print "Scattering parameters succesfully saved to "+model.simulation_name+".dat"
meep.all_wait() # Wait until all file operations are finished