~nizamov-shawkat/python-meep/devel : revision 15

1

PYTHON-MEEP BINDING DOCUMENTATION

2

====================================

3

4

Author of this documentation : Emmanuel.Lambert@intec.ugent.be

5

(reverse-engineered from samples and source code).

6

7

Document history :

8

9

::

10

11

* 19/20/21-08-2009 : document creation

12

* 24-08-2009 : small improvements & clarifications.

13

* 25/26-08-2009 : sections 7 & 8 were added.

14

* 03-04/09/2009 : -a new class "structure_eps_pml" was introduced to make eps averaging the default behaviour (consistent with the Scheme behaviour).

15

-port to Meep 1.1.1 (class 'volume' was renamed to 'grid_volume' and class 'geometric_volume' to 'volume'

16

-minor changes in the bent waveguide sample, to make it more consistent with the Scheme version

17

18

19

**0. Installing the python-meep binding**

20

-----------------------------------------

21

22

Currently the following platform is supported : Ubuntu 9.04 - 32bit.

23

A debian package is under construction. For the moment, contact Manu and we'll go through the manual installation process.

24

25

**1. The general structure of a python-meep program**

26

-----------------------------------------------------

27

28

In general terms, a python-meep program can be structured as follows :

29

30

* import the python-meep binding :

31

``from meep import *``

32

This will load the library ``_meep.so`` and Python-files ``meep.py`` and ``python_meep.py`` from path ``/usr/local/lib/python2.6/dist-packages/``

33

34

* initialize Meep

35

We create an object of type MPI, giving the program arguments (``sys.argv``) as initialisation parameter (by default these are empty, which means that we will run on a single processor system -- *to be further elaborated*).

36

37

``meep = MPI(sys.argsv)``

38

39

* define a computational grid volume

40

See section 2 below which explains usage of the ``grid_volume`` class.

41

42

* define the waveguide structure (describing the geometry, PML and materials)

43

See section 3 below which explains usage of the ``structure_eps_pml`` class.

44

45

* create an object which will hold the calculated fields

46

See section 4 below which explains usage of the ``field`` class.

47

48

* define the sources

49

See section 5 below which explains usage of the ``add_point_source`` and ``add_volume_source`` functions.

50

51

* run the simulation (iterate over the time-steps)

52

See section 6 below.

53

54

55

**2. Defining the computational grid volume**

56

---------------------------------------------

57

58

The following set of 'factory functions' is provided, aimed at creating a ``grid_volume`` object. The first arguments define the size of the computational volume, the last argument is the computational grid resolution (in pixels per distance unit).

59

* ``volcyl(rsize, zsize, resolution)``

60

Defines a cyclical computational grid volume.

61

* ``volone(zsize, resolution)`` *alternatively called* ``vol1d(zsize, resolution)``

62

Defines a 1-dimensional computational grid volume along the Z-axis.

63

* ``voltwo(xsize, ysize, resolution)`` *alternatively called* ``vol2d(xsize, ysize, a)``

64

Defines a 2-dimensional computational grid volumes along the X- and Y-axes

65

* ``vol3d(xsize, ysize, zsize, resolution)``

66

Defines a 3-dimensional computational grid volume.

67

68

e.g.: ``v = volone(6, 10)`` defines a 1-dimensional computational volume of lenght 6, with 10 pixels per distance unit.

69

70

71

**3. Defining the waveguide structure**

72

---------------------------------------

73

74

The waveguide structure is defined using the class ``structure_eps_pml``, of which the constructor has the following arguments :

75

76

* *required* : the computational grid volume (a reference to an object of type ``grid_volume``, see section 2 above)

77

78

* *required* : a function defining the dielectric properties of the materials in the computational grid volume (thus describing the actual waveguide structure). For all-air, the predefined function 'one' can be used (epsilon = constant value 1). See note 3.1 below for more information about defining your own custom material function.

79

80

* *optional* : a boundary region: this is a reference to an object of type ``boundary_region``. There are a number of predefined functions that can be used to create such an object :

81

- ``no_pml()`` describing a conditionless boundary region (no PML)

82

- ``pml(thickness)`` : decribing a perfectly matching layer (PML) of a certain thickness (double value) on the boundaries in all directions.

83

- ``pml(thickness, direction)`` : decribing a perfectly matching layer (PML) of a certain thickness (double value) in a certain direction (``X, Y, Z, R or P``).

84

- ``pml(thickness, direction, boundary_side)`` : describing a PML of a certain thickness (double value), in a certain direction (``X, Y, Z, R or P``) and on the ``High`` or ``Low`` side. E.g. if boundary_side is ``Low`` and direction is ``X``, then a PML layer is added to the −x boundary. The default puts PML layers on both sides of the direction.

85

86

* *optional* : a parameter to disable (value 0) or enable (value 1) EPS-averaging (anisotrpic averaging). EPS-averaging is enabled by default, making this consistent with the behaviour in Scheme.

87

88

* *optional* : a function defining a symmetry to exploit, in order to speed up the FDTD calculation (reference to an object of type ``symmetry``). The following predefined functions can be used to create a ``symmetry`` object:

89

- identity : no symmetry

90

- rotate4(direction, grid_volume) : defines a 90° rotational symmetry with 'direction' the axis of rotation.

91

- rotate2(direction, grid_volume) : defines a 180° rotational symmetry with 'direction' the axis of rotation.

92

- mirror(direction, grid_volume) : defines a mirror symmetry plane with 'direction' normal to the mirror plane.

93

- r_to_minus_r_symmetry : defines a mirror symmetry in polar coordinates

94

95

* optional : the number of chunks in which to split up the calculated geometry (used for multi-processor configurations, only if you have set this up in the ``MPI`` initialization, mentioned in section 1 above)

96

97

e.g. : if ``v`` is a variable pointing to the computational grid volume, then :

98

``s = structure_eps_pml(v, one) `` defines a structure with all air (eps=1),

99

which is equivalent to:

100

``s = structure_eps_pml(v, one, no_pml(), 1, identity(), 1)``

101

102

Another example : ``s = structure_eps_pml(v, EPS, pml(0.1,Y) )`` with EPS a custom material function, which is explained in the note below.

103

104

To disable EPS-averaging : ``s = structure_eps_pml(v, EPS, pml(0.1,Y), 0 )``.

105

106

107

note 3.1. Defining a material function

108

________________________________________

109

110

In order to describe the geometry of the waveguide, we have to provide a 'material function' returning the the dielectric variable epsilon as a function of the position (identified by a vector). In python-meep, we can do this by defining a class that inherits from class ``Callback``. Through this inheritance, the core meep library (written in C++) will be able to call back to the Python function which describes the material properties.

111

112

E.g. :

113

114

::

115

116

class epsilon(Callback): #inherit from Callback for integration with the meep core library

117

def __init__(self):

118

Callback.__init__(self)

119

def double_vec(self,vec): #override of function in the Callback class to set the eps function

120

self.set_double(self.eps(vec))

121

return

122

def eps(self,vec): #return the epsilon value for a certain point (indicated by the vector v)

123

v = vec

124

r = v.x()*v.x() + v.y()*v.y()

125

dr = sqrt(r)

126

while dr>1:

127

dr-=1

128

if dr > 0.7001:

129

return 12.0

130

return 1.0

131

132

Please note the brackets when referring to the x- and y-components of the vector ``vec``. These are crucial : no error message will be thrown if you refer to it as vec.x or vec.y, but the value will always be zero.

133

So, one should write : ``vec.x()`` and ``vec.y()``.

134

135

The meep-python library has a 'global' variable EPS, which is used as a reference for communication between the Meep core library and the Python code. We assign our epsilon-function as follows to the global EPS variable :

136

137

::

138

139

set_EPS_Callback(epsilon().__disown__())

140

s = structure(v, EPS, no_pml(), identity(), 1)

141

142

143

144

The call to function ``__disown__()`` is for memory management purposes and is *absolutely required*. An improvement of the python-meep binding could consist of making this call transparant for the end user. But for now, the user must manually provide it.

145

146

***Important remark*** : at the end of our program, we should call : ``del_EPS_Callback()`` in order to clean up the global variable.

147

148

**4. Defining the initial field**

149

---------------------------------

150

151

We create an object of type ``fields``, which will contain the calculated field.

152

153

We must first create a Python class that inherits from class ``Callback`` and that will define the function for initialization of the field. Inheritance from class ``Callback`` is required, because the core meep library (written in C++) will have to call back to the Python function. For example, let's call our initialization class ``fi``.

154

155

::

156

157

class fi(Callback): #inherit from Callback for integration with the meep core library

158

def __init__(self):

159

Callback.__init__(self)

160

def complex_vec(self,v): #override of function in the Callback class to set the field initialization function

161

self.set_complex(self.f(v))

162

return

163

def f(self,v): #return the field value for a certain point (indicated by the vector v)

164

return 1.0

165

166

The meep-python library has a 'global' variable INIF, that is used to bind the meep core library to our Python field initialization class. To set INIF, we have to use the following statement :

167

168

``set_INIF_Callback(fi().__disown__()) #link the INIF variable to the fi class``

169

170

We refer to section 3-note1 for more information about the function ``__disown__()``.

171

172

E.g.: If ``s`` is a variable pointing to the structure and ``comp`` denotes the component which we are initializing, then the complete field initialization code looks as follows :

173

174

::

175

176

f = fields(s)

177

comp = Hy

178

f.initialize_field(comp, INIF)

179

180

181

***Important remark*** : at the end of our program, we should then call : ``del_INIF_Callback()`` in order to clean up the global variable.

182

183

The call to ``initialize_field`` is not mandatory. If the initial conditions are zeros for all components, then we can rely on the automatic initialization at creation of the object.

184

185

We can additionally define **Bloch-periodic boundary conditions** over the field. This is done with the ``use_bloch`` function of the field class, e.g. :

186

187

``f.use_bloch(vec(0.0))``

188

189

*to be further elaborated - what is the exact meaning of the vector argument? (not well understood at this time)*

190

191

192

193

**5. Defining the sources**

194

---------------------------

195

196

The definition of the current sources can be done through 2 functions of the ``fields`` class :

197

* ``add_point_source(component, src_time, vec, complex)``

198

* ``add_volume_source(component, src_time, volume, complex)``

199

200

201

Each require as arguments an electromagnetic component (e.g. ``Ex, Ey, ...`` and ``Hx, Hy, ...``) and an object of type ``src_time``, which specifies the time dependence of the source (see below).

202

203

For a point source, we must specify the center point of the current source using a vector (object of type ``vec``).

204

205

For a volume source, we must specify an object of type ``volume`` (*to be elablorated*).

206

207

The last argument is an overall complex amplitude number, multiplying the current source (default 1.0).

208

209

The following variants are available :

210

* ``add_point_source(component, double, double, double, double, vec centerpoint, complex amplitude, int is_continuous)``

211

* This is a shortcut function so that no ``src_time`` object must be created. **This function is preferably used for point sources.

212

* The four real arguments define the central frequency, spectral width, peaktime and cutoff.

213

214

* ``add_volume_source(component, src_time, volume, complex A(vec), complex amp)``

215

* A is a function with a single argument (a vector) describing the current amplitude for a certain point. The position argument is relative to the center of the current source.

216

* *to be further elaborated : what is then the difference between A and amp ? (not well understood)*

217

218

If you want to have full control over the time dependence of the current source, you can create a subclass of class ``src_time`` to specify the time dependence of the current source. This subclass should implement/override the following methodes :

219

220

* ``complex dipole(double time)`` : returns the dipole moment of the source as a function of time

221

* ``double last_time()`` : returns the end time of the source

222

* ``complex frequency()`` : returns the frequency of the source

223

* ``src_time clone() `` (to create a deep copy) and ``bool is_equal(src_time)`` (to compare 2 ``src_time`` objects)

224

225

For example :

226

227

::

228

229

class mySource(src_time):

230

def __init__(self):

231

src_time.__init__(self)

232

233

def dipole(self, t):

234

...

235

return ...

236

237

def last_time(self):

238

.... #your functionality

239

return ...

240

241

def frequency(self):

242

.... #your functionality

243

return ...

244

245

246

247

Give an instance of this class as second argument to the ``add_point_source`` or ``add_volume_source`` functions.

248

249

Three classes, inheriting from ``src_time``, are predefined and can be used off the shelf :

250

* ``gaussian_src_time`` for a Gaussian-pulse source. The constructor demands 2 arguments of type double :

251

* the center frequency ω, in units of 2πc

252

* the frequency width w used in the Gaussian

253

* ``continuous_src_time`` for a continuous-wave source proportional to exp(−iωt). The constructor demands 4 arguments :

254

* the frequency ω, in units 2πc/distance (complex number)

255

* the temporal width of smoothing (default 0)

256

* the start time (default 0)

257

* the end time (default infinity = never turn off)

258

* ``custom_src_time`` for a user-specified source function f(t) with start/end times. The constructor demands 4 arguments :

259

* The function f(t) specifying the time-dependence of the source

260

* *...(2nd argument unclear, to be further elaborated)...*

261

* the start time of the source (default -infinity)

262

* the end time of the source (default +infinity)

263

264

For example, in order to define a continuous line source of length 1, from point (6,3) to point (6,4) in 2-dimensional geometry :

265

266

::

267

268

#define a continuous source

269

srcFreq = 0.125

270

srcWidth = 20

271

src = continuous_src_time(srcFreq, srcWidth, 0, infinity)

272

srcComp = Ez

273

#make it a line source of size 1 starting on position(6,3)

274

srcGeo = volume(vec(6,3),vec(6,4))

275

f.add_volume_source(srcComp, src, srcGeo, 1)

276

277

278

**6. Running the simulation, retrieving field values and exporting HDF5 files**

279

-------------------------------------------------------------------------------

280

281

We can now time-step and retrieve various field values along the way.

282

The actual time step value can be retrieved or set through the variable ``f.dt``.

283

284

The default time step in Meep is : ``Courant factor / resolution`` (in FDTD, the Courant factor relates the time step size to the spatial discretization: cΔt = SΔx. Default for S is 0.5). If no further parametrization is done, then this default value is used.

285

286

To trigger a step in time, you call the function ``f.step()``.

287

288

To step until the source has fully decayed :

289

290

::

291

292

while (f.time() < f.last_source_time()):

293

f.step()

294

295

The function ``runUntilFieldsDecayed`` mimicks the behaviour of 'stop-when-fields-decayed' in Meep-Scheme.

296

This will run time steps until the source has decayed to 0.001 of the peak amplitude. After that, an additional 50 time steps will be run.

297

We refer to section 8 below where this function is applied in an example.

298

299

A rich functionality is available for retrieving field information. Some examples :

300

301

* ``f.energy_in_box(v.surroundings())``

302

* ``f.electric_energy_in_box(v.surroundings())``

303

* ``f.magnetic_energy_in_box(v.surroundings())``

304

* ``f.thermo_energy_in_box(v.surroundings())``

305

* ``f.total_energy()``

306

* ``f.field_energy_in_box(v.surroundings())``

307

* ``f.field_energy_in_box(component, v.surroundings())`` where the first argument is the electromagnetic component (``Ex, Ey, Ez, Er, Ep, Hx, Hy, Hz, Hr, Hp, Dx, Dy, Dz, Dp, Dr, Bx, By, Bz, Bp, Br, Dielectric`` or ``Permeability``)

308

* ``f.flux_in_box(X, v.surroundings())`` where the first argument is the direction (``X, Y, Z, R`` or ``P``)

309

310

We can probe the field at certain points by defining a *monitor point* as follows :

311

312

::

313

314

m = monitor_point()

315

p = vec(2.10) #vector identifying the point that we want to probe

316

f.get_point(m, p)

317

m.get_component(Hx)

318

319

We can export the dielectric function and e.g. the Ex component of the field to HDF5 files as follows :

320

321

::

322

323

#make sure you start your Python session with 'sudo' or write rights to the current path

324

feps = prepareHDF5File("eps.h5")

325

f.output_hdf5(Dielectric, v.surroundings(), feps) #export the Dielectric structure so that we can visually verify it

326

fex = prepareHDF5File("ex.h5")

327

while (f.time() < f.last_source_time()):

328

f.step()

329

f.output_hdf5(Ex, v.surroundings(), fex, 1) #export the Ex component, appending to the file "ex.h5"

330

331

332

333

**7. Defining and retrieving fluxes**

334

--------------------------------------

335

336

First we define a flux plane.

337

This is done through the creation of an object of type ``volume`` (specifying 2 vectors as arguments).

338

339

Then we apply this flux plane to the field, specifying 4 parameters :

340

* the reference to the ``volume`` object

341

* the minimum frequency (in the example below, this is ``srcFreqCenter-(srcPulseWidth/2.0)``)

342

* the maximum frequency (in the example below this is ``srcFreqCenter+(srcPulseWidth/2.0)`` )

343

* the number of discrete frequencies that we want to monitor in the flux (in the example below, this is 100).

344

345

After running the simulation, we can retrieve the flux values through the function ``getFluxData()`` : this returns a 2-dimensional array with the frequencies and actual flux values.

346

347

E.g., if ``f`` is the field, then we proceed as follows :

348

349

::

350

351

#define the flux plane and flux parameters

352

fluxplane = volume(vec(1,2),vec(1,3))

353

flux = f.add_dft_flux_plane(fluxplane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

354

355

#run the calculation

356

while (f.time() < f.last_source_time()):

357

f.step()

358

359

#retrieve the flux data

360

fluxdata = getFluxData(flux)

361

frequencies = fluxdata[0]

362

fluxvalues = fluxdata[1]

363

364

365

366

**8. The "90 degree bent waveguide example in Python**

367

------------------------------------------------------

368

369

We have ported the "90 degree bent waveguide" example from the Meep-Scheme tutorial to Python.

370

371

The original example can be found on the following URL : http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial

372

(section 'A 90° bend').

373

374

You can find the source code below.

375

376

The following animated gifs can be produced from the HDF5-files (see the script included in directory 'samples') :

377

378

.. image:: images/bentwgNB.gif

379

380

.. image:: images/bentwgB.gif

381

382

383

And here is the graph of the transmission, reflection and loss fluxes, showing the same results as the example in Scheme:

384

385

.. image:: images/fluxes.png

386

:height: 315

387

:width: 443

388

389

390

::

391

392

from meep import *

393

from math import *

394

from python_meep import *

395

import numpy

396

import matplotlib.pyplot

397

import sys

398

399

res = 10.0

400

gridSizeX = 16.0

401

gridSizeY = 32.0

402

padSize = 4.0

403

wgLengthX = gridSizeX - padSize

404

wgLengthY = gridSizeY - padSize

405

wgWidth = 1.0 #width of the waveguide

406

wgHorYCen = padSize + wgWidth/2.0 #horizontal waveguide center Y-pos

407

wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)

408

srcFreqCenter = 0.15 #gaussian source center frequency

409

srcPulseWidth = 0.1 #gaussian source pulse width

410

srcComp = Ez #gaussian source component

411

412

#this function plots the waveguide material as a function of a vector(X,Y)

413

class epsilon(Callback):

414

def __init__(self, pIsWgBent):

415

Callback.__init__(self)

416

self.isWgBent = pIsWgBent

417

def double_vec(self,vec):

418

self.set_double(self.eps(vec))

419

return

420

def eps(self,vec):

421

if (self.isWgBent): #there is a bend

422

if ((vec.x()<wgLengthX) and (vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):

423

return 12.0

424

elif ((vec.x()>=wgLengthX-wgWidth) and (vec.x()<=wgLengthX) and vec.y()>= padSize ):

425

return 12.0

426

else:

427

return 1.0

428

else: #there is no bend

429

if ((vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):

430

return 12.0

431

else:

432

return 1.0

433

434

def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):

435

#we create a structure with PML of thickness = 1.0 on all boundaries,

436

#in all directions,

437

#using the material function EPS

438

material = epsilon(pIsWgBent)

439

set_EPS_Callback(material.__disown__())

440

s = structure_eps_pml(pCompVol, EPS, pml(1.0) )

441

f = fields(s)

442

#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize

443

srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )

444

srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))

445

f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)

446

print "Field created..."

447

return f

448

449

450

meep = MPI(sys.argv)

451

452

#FIRST WE WORK OUT THE CASE WITH NO BEND

453

#----------------------------------------------------------------

454

print "*1* Starting the case with no bend..."

455

#create the computational grid

456

noBendVol = voltwo(gridSizeX,gridSizeY,res)

457

458

#create the field

459

wgBent = 0 #no bend

460

noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

461

462

bendFnEps = "./bentwgNB_Eps.h5"

463

bendFnEz = "./bentwgNB_Ez.h5"

464

#export the dielectric structure (so that we can visually verify the waveguide structure)

465

noBendDielectricFile = prepareHDF5File(bendFnEps)

466

noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)

467

468

#create the file for the field components

469

noBendFileOutputEz = prepareHDF5File(bendFnEz)

470

471

#define the flux plane for the reflected flux

472

noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane

473

noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))

474

noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

475

476

#define the flux plane for the transmitted flux

477

noBendTransmfluxPlaneXPos = gridSizeX - 1.5; #the X-coordinate of our transmission flux plane

478

noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))

479

noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

480

481

print "Calculating..."

482

noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed

483

runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)

484

print "Done..!"

485

486

#construct 2-dimensional array with the flux plane data

487

#see python_meep.py

488

noBendReflFlux = getFluxData(noBendReflectedFlux)

489

noBendTransmFlux = getFluxData(noBendTransmFlux)

490

491

#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'

492

noBendReflectedFlux.scale_dfts(-1);

493

f = open("minusflux.h5", 'w') #truncate file if already exists

494

f.close()

495

noBendReflectedFlux.save_hdf5(noBendField, "minusflux")

496

497

del_EPS_Callback()

498

499

500

#AND SECONDLY FOR THE CASE WITH BEND

501

#----------------------------------------------------------------

502

print "*2* Starting the case with bend..."

503

#create the computational grid

504

bendVol = voltwo(gridSizeX,gridSizeY,res)

505

506

#create the field

507

wgBent = 1 #there is a bend

508

bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)

509

510

#export the dielectric structure (so that we can visually verify the waveguide structure)

511

bendFnEps = "./bentwgB_Eps.h5"

512

bendFnEz = "./bentwgB_Ez.h5"

513

bendDielectricFile = prepareHDF5File(bendFnEps)

514

bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)

515

516

#create the file for the field components

517

bendFileOutputEz = prepareHDF5File(bendFnEz)

518

519

#define the flux plane for the reflected flux

520

bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane

521

bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))

522

bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)

523

524

#load the minused reflection flux from the "no bend" case as initalization

525

bendReflectedFlux.load_hdf5(bendField, "minusflux")

526

527

528

#define the flux plane for the transmitted flux

529

bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane

530

bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))

531

bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )

532

533

print "Calculating..."

534

bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed

535

runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)

536

print "Done..!"

537

538

#construct 2-dimensional array with the flux plane data

539

#see python_meep.py

540

bendReflFlux = getFluxData(bendReflectedFlux)

541

bendTransmFlux = getFluxData(bendTransmFlux)

542

543

del_EPS_Callback()

544

545

#SHOW THE RESULTS IN A PLOT

546

frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...

547

Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])]

548

Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]

549

Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]

550

551

matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')

552

matplotlib.pyplot.plot(frequencies, Prefl, 'ro')

553

matplotlib.pyplot.plot(frequencies, Ploss, 'go' )

554

555

matplotlib.pyplot.show()

556

557

558