1
=== added file 'AUTHORS'
2
--- AUTHORS	1970-01-01 00:00:00 +0000
3
+++ AUTHORS	2009-12-01 14:23:09 +0000
4
@@ -0,0 +1,2 @@
5
1
Shawkat Nizamov <nizamov.shawkat@gmail.com>
6
2
Emmanuel Lambert <Emmanuel.Lambert@intec.ugent.be>
7
0
3
8
=== renamed file 'AUTHORS' => 'AUTHORS.moved'
9
=== added file 'COPYING'
10
--- COPYING	1970-01-01 00:00:00 +0000
11
+++ COPYING	2009-12-01 14:23:09 +0000
12
@@ -0,0 +1,674 @@
13
1
                    GNU GENERAL PUBLIC LICENSE
14
2
                       Version 3, 29 June 2007
15
3
16
4
 Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
17
5
 Everyone is permitted to copy and distribute verbatim copies
18
6
 of this license document, but changing it is not allowed.
19
7
20
8
                            Preamble
21
9
22
10
  The GNU General Public License is a free, copyleft license for
23
11
software and other kinds of works.
24
12
25
13
  The licenses for most software and other practical works are designed
26
14
to take away your freedom to share and change the works.  By contrast,
27
15
the GNU General Public License is intended to guarantee your freedom to
28
16
share and change all versions of a program--to make sure it remains free
29
17
software for all its users.  We, the Free Software Foundation, use the
30
18
GNU General Public License for most of our software; it applies also to
31
19
any other work released this way by its authors.  You can apply it to
32
20
your programs, too.
33
21
34
22
  When we speak of free software, we are referring to freedom, not
35
23
price.  Our General Public Licenses are designed to make sure that you
36
24
have the freedom to distribute copies of free software (and charge for
37
25
them if you wish), that you receive source code or can get it if you
38
26
want it, that you can change the software or use pieces of it in new
39
27
free programs, and that you know you can do these things.
40
28
41
29
  To protect your rights, we need to prevent others from denying you
42
30
these rights or asking you to surrender the rights.  Therefore, you have
43
31
certain responsibilities if you distribute copies of the software, or if
44
32
you modify it: responsibilities to respect the freedom of others.
45
33
46
34
  For example, if you distribute copies of such a program, whether
47
35
gratis or for a fee, you must pass on to the recipients the same
48
36
freedoms that you received.  You must make sure that they, too, receive
49
37
or can get the source code.  And you must show them these terms so they
50
38
know their rights.
51
39
52
40
  Developers that use the GNU GPL protect your rights with two steps:
53
41
(1) assert copyright on the software, and (2) offer you this License
54
42
giving you legal permission to copy, distribute and/or modify it.
55
43
56
44
  For the developers' and authors' protection, the GPL clearly explains
57
45
that there is no warranty for this free software.  For both users' and
58
46
authors' sake, the GPL requires that modified versions be marked as
59
47
changed, so that their problems will not be attributed erroneously to
60
48
authors of previous versions.
61
49
62
50
  Some devices are designed to deny users access to install or run
63
51
modified versions of the software inside them, although the manufacturer
64
52
can do so.  This is fundamentally incompatible with the aim of
65
53
protecting users' freedom to change the software.  The systematic
66
54
pattern of such abuse occurs in the area of products for individuals to
67
55
use, which is precisely where it is most unacceptable.  Therefore, we
68
56
have designed this version of the GPL to prohibit the practice for those
69
57
products.  If such problems arise substantially in other domains, we
70
58
stand ready to extend this provision to those domains in future versions
71
59
of the GPL, as needed to protect the freedom of users.
72
60
73
61
  Finally, every program is threatened constantly by software patents.
74
62
States should not allow patents to restrict development and use of
75
63
software on general-purpose computers, but in those that do, we wish to
76
64
avoid the special danger that patents applied to a free program could
77
65
make it effectively proprietary.  To prevent this, the GPL assures that
78
66
patents cannot be used to render the program non-free.
79
67
80
68
  The precise terms and conditions for copying, distribution and
81
69
modification follow.
82
70
83
71
                       TERMS AND CONDITIONS
84
72
85
73
  0. Definitions.
86
74
87
75
  "This License" refers to version 3 of the GNU General Public License.
88
76
89
77
  "Copyright" also means copyright-like laws that apply to other kinds of
90
78
works, such as semiconductor masks.
91
79
92
80
  "The Program" refers to any copyrightable work licensed under this
93
81
License.  Each licensee is addressed as "you".  "Licensees" and
94
82
"recipients" may be individuals or organizations.
95
83
96
84
  To "modify" a work means to copy from or adapt all or part of the work
97
85
in a fashion requiring copyright permission, other than the making of an
98
86
exact copy.  The resulting work is called a "modified version" of the
99
87
earlier work or a work "based on" the earlier work.
100
88
101
89
  A "covered work" means either the unmodified Program or a work based
102
90
on the Program.
103
91
104
92
  To "propagate" a work means to do anything with it that, without
105
93
permission, would make you directly or secondarily liable for
106
94
infringement under applicable copyright law, except executing it on a
107
95
computer or modifying a private copy.  Propagation includes copying,
108
96
distribution (with or without modification), making available to the
109
97
public, and in some countries other activities as well.
110
98
111
99
  To "convey" a work means any kind of propagation that enables other
112
100
parties to make or receive copies.  Mere interaction with a user through
113
101
a computer network, with no transfer of a copy, is not conveying.
114
102
115
103
  An interactive user interface displays "Appropriate Legal Notices"
116
104
to the extent that it includes a convenient and prominently visible
117
105
feature that (1) displays an appropriate copyright notice, and (2)
118
106
tells the user that there is no warranty for the work (except to the
119
107
extent that warranties are provided), that licensees may convey the
120
108
work under this License, and how to view a copy of this License.  If
121
109
the interface presents a list of user commands or options, such as a
122
110
menu, a prominent item in the list meets this criterion.
123
111
124
112
  1. Source Code.
125
113
126
114
  The "source code" for a work means the preferred form of the work
127
115
for making modifications to it.  "Object code" means any non-source
128
116
form of a work.
129
117
130
118
  A "Standard Interface" means an interface that either is an official
131
119
standard defined by a recognized standards body, or, in the case of
132
120
interfaces specified for a particular programming language, one that
133
121
is widely used among developers working in that language.
134
122
135
123
  The "System Libraries" of an executable work include anything, other
136
124
than the work as a whole, that (a) is included in the normal form of
137
125
packaging a Major Component, but which is not part of that Major
138
126
Component, and (b) serves only to enable use of the work with that
139
127
Major Component, or to implement a Standard Interface for which an
140
128
implementation is available to the public in source code form.  A
141
129
"Major Component", in this context, means a major essential component
142
130
(kernel, window system, and so on) of the specific operating system
143
131
(if any) on which the executable work runs, or a compiler used to
144
132
produce the work, or an object code interpreter used to run it.
145
133
146
134
  The "Corresponding Source" for a work in object code form means all
147
135
the source code needed to generate, install, and (for an executable
148
136
work) run the object code and to modify the work, including scripts to
149
137
control those activities.  However, it does not include the work's
150
138
System Libraries, or general-purpose tools or generally available free
151
139
programs which are used unmodified in performing those activities but
152
140
which are not part of the work.  For example, Corresponding Source
153
141
includes interface definition files associated with source files for
154
142
the work, and the source code for shared libraries and dynamically
155
143
linked subprograms that the work is specifically designed to require,
156
144
such as by intimate data communication or control flow between those
157
145
subprograms and other parts of the work.
158
146
159
147
  The Corresponding Source need not include anything that users
160
148
can regenerate automatically from other parts of the Corresponding
161
149
Source.
162
150
163
151
  The Corresponding Source for a work in source code form is that
164
152
same work.
165
153
166
154
  2. Basic Permissions.
167
155
168
156
  All rights granted under this License are granted for the term of
169
157
copyright on the Program, and are irrevocable provided the stated
170
158
conditions are met.  This License explicitly affirms your unlimited
171
159
permission to run the unmodified Program.  The output from running a
172
160
covered work is covered by this License only if the output, given its
173
161
content, constitutes a covered work.  This License acknowledges your
174
162
rights of fair use or other equivalent, as provided by copyright law.
175
163
176
164
  You may make, run and propagate covered works that you do not
177
165
convey, without conditions so long as your license otherwise remains
178
166
in force.  You may convey covered works to others for the sole purpose
179
167
of having them make modifications exclusively for you, or provide you
180
168
with facilities for running those works, provided that you comply with
181
169
the terms of this License in conveying all material for which you do
182
170
not control copyright.  Those thus making or running the covered works
183
171
for you must do so exclusively on your behalf, under your direction
184
172
and control, on terms that prohibit them from making any copies of
185
173
your copyrighted material outside their relationship with you.
186
174
187
175
  Conveying under any other circumstances is permitted solely under
188
176
the conditions stated below.  Sublicensing is not allowed; section 10
189
177
makes it unnecessary.
190
178
191
179
  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
192
180
193
181
  No covered work shall be deemed part of an effective technological
194
182
measure under any applicable law fulfilling obligations under article
195
183
11 of the WIPO copyright treaty adopted on 20 December 1996, or
196
184
similar laws prohibiting or restricting circumvention of such
197
185
measures.
198
186
199
187
  When you convey a covered work, you waive any legal power to forbid
200
188
circumvention of technological measures to the extent such circumvention
201
189
is effected by exercising rights under this License with respect to
202
190
the covered work, and you disclaim any intention to limit operation or
203
191
modification of the work as a means of enforcing, against the work's
204
192
users, your or third parties' legal rights to forbid circumvention of
205
193
technological measures.
206
194
207
195
  4. Conveying Verbatim Copies.
208
196
209
197
  You may convey verbatim copies of the Program's source code as you
210
198
receive it, in any medium, provided that you conspicuously and
211
199
appropriately publish on each copy an appropriate copyright notice;
212
200
keep intact all notices stating that this License and any
213
201
non-permissive terms added in accord with section 7 apply to the code;
214
202
keep intact all notices of the absence of any warranty; and give all
215
203
recipients a copy of this License along with the Program.
216
204
217
205
  You may charge any price or no price for each copy that you convey,
218
206
and you may offer support or warranty protection for a fee.
219
207
220
208
  5. Conveying Modified Source Versions.
221
209
222
210
  You may convey a work based on the Program, or the modifications to
223
211
produce it from the Program, in the form of source code under the
224
212
terms of section 4, provided that you also meet all of these conditions:
225
213
226
214
    a) The work must carry prominent notices stating that you modified
227
215
    it, and giving a relevant date.
228
216
229
217
    b) The work must carry prominent notices stating that it is
230
218
    released under this License and any conditions added under section
231
219
    7.  This requirement modifies the requirement in section 4 to
232
220
    "keep intact all notices".
233
221
234
222
    c) You must license the entire work, as a whole, under this
235
223
    License to anyone who comes into possession of a copy.  This
236
224
    License will therefore apply, along with any applicable section 7
237
225
    additional terms, to the whole of the work, and all its parts,
238
226
    regardless of how they are packaged.  This License gives no
239
227
    permission to license the work in any other way, but it does not
240
228
    invalidate such permission if you have separately received it.
241
229
242
230
    d) If the work has interactive user interfaces, each must display
243
231
    Appropriate Legal Notices; however, if the Program has interactive
244
232
    interfaces that do not display Appropriate Legal Notices, your
245
233
    work need not make them do so.
246
234
247
235
  A compilation of a covered work with other separate and independent
248
236
works, which are not by their nature extensions of the covered work,
249
237
and which are not combined with it such as to form a larger program,
250
238
in or on a volume of a storage or distribution medium, is called an
251
239
"aggregate" if the compilation and its resulting copyright are not
252
240
used to limit the access or legal rights of the compilation's users
253
241
beyond what the individual works permit.  Inclusion of a covered work
254
242
in an aggregate does not cause this License to apply to the other
255
243
parts of the aggregate.
256
244
257
245
  6. Conveying Non-Source Forms.
258
246
259
247
  You may convey a covered work in object code form under the terms
260
248
of sections 4 and 5, provided that you also convey the
261
249
machine-readable Corresponding Source under the terms of this License,
262
250
in one of these ways:
263
251
264
252
    a) Convey the object code in, or embodied in, a physical product
265
253
    (including a physical distribution medium), accompanied by the
266
254
    Corresponding Source fixed on a durable physical medium
267
255
    customarily used for software interchange.
268
256
269
257
    b) Convey the object code in, or embodied in, a physical product
270
258
    (including a physical distribution medium), accompanied by a
271
259
    written offer, valid for at least three years and valid for as
272
260
    long as you offer spare parts or customer support for that product
273
261
    model, to give anyone who possesses the object code either (1) a
274
262
    copy of the Corresponding Source for all the software in the
275
263
    product that is covered by this License, on a durable physical
276
264
    medium customarily used for software interchange, for a price no
277
265
    more than your reasonable cost of physically performing this
278
266
    conveying of source, or (2) access to copy the
279
267
    Corresponding Source from a network server at no charge.
280
268
281
269
    c) Convey individual copies of the object code with a copy of the
282
270
    written offer to provide the Corresponding Source.  This
283
271
    alternative is allowed only occasionally and noncommercially, and
284
272
    only if you received the object code with such an offer, in accord
285
273
    with subsection 6b.
286
274
287
275
    d) Convey the object code by offering access from a designated
288
276
    place (gratis or for a charge), and offer equivalent access to the
289
277
    Corresponding Source in the same way through the same place at no
290
278
    further charge.  You need not require recipients to copy the
291
279
    Corresponding Source along with the object code.  If the place to
292
280
    copy the object code is a network server, the Corresponding Source
293
281
    may be on a different server (operated by you or a third party)
294
282
    that supports equivalent copying facilities, provided you maintain
295
283
    clear directions next to the object code saying where to find the
296
284
    Corresponding Source.  Regardless of what server hosts the
297
285
    Corresponding Source, you remain obligated to ensure that it is
298
286
    available for as long as needed to satisfy these requirements.
299
287
300
288
    e) Convey the object code using peer-to-peer transmission, provided
301
289
    you inform other peers where the object code and Corresponding
302
290
    Source of the work are being offered to the general public at no
303
291
    charge under subsection 6d.
304
292
305
293
  A separable portion of the object code, whose source code is excluded
306
294
from the Corresponding Source as a System Library, need not be
307
295
included in conveying the object code work.
308
296
309
297
  A "User Product" is either (1) a "consumer product", which means any
310
298
tangible personal property which is normally used for personal, family,
311
299
or household purposes, or (2) anything designed or sold for incorporation
312
300
into a dwelling.  In determining whether a product is a consumer product,
313
301
doubtful cases shall be resolved in favor of coverage.  For a particular
314
302
product received by a particular user, "normally used" refers to a
315
303
typical or common use of that class of product, regardless of the status
316
304
of the particular user or of the way in which the particular user
317
305
actually uses, or expects or is expected to use, the product.  A product
318
306
is a consumer product regardless of whether the product has substantial
319
307
commercial, industrial or non-consumer uses, unless such uses represent
320
308
the only significant mode of use of the product.
321
309
322
310
  "Installation Information" for a User Product means any methods,
323
311
procedures, authorization keys, or other information required to install
324
312
and execute modified versions of a covered work in that User Product from
325
313
a modified version of its Corresponding Source.  The information must
326
314
suffice to ensure that the continued functioning of the modified object
327
315
code is in no case prevented or interfered with solely because
328
316
modification has been made.
329
317
330
318
  If you convey an object code work under this section in, or with, or
331
319
specifically for use in, a User Product, and the conveying occurs as
332
320
part of a transaction in which the right of possession and use of the
333
321
User Product is transferred to the recipient in perpetuity or for a
334
322
fixed term (regardless of how the transaction is characterized), the
335
323
Corresponding Source conveyed under this section must be accompanied
336
324
by the Installation Information.  But this requirement does not apply
337
325
if neither you nor any third party retains the ability to install
338
326
modified object code on the User Product (for example, the work has
339
327
been installed in ROM).
340
328
341
329
  The requirement to provide Installation Information does not include a
342
330
requirement to continue to provide support service, warranty, or updates
343
331
for a work that has been modified or installed by the recipient, or for
344
332
the User Product in which it has been modified or installed.  Access to a
345
333
network may be denied when the modification itself materially and
346
334
adversely affects the operation of the network or violates the rules and
347
335
protocols for communication across the network.
348
336
349
337
  Corresponding Source conveyed, and Installation Information provided,
350
338
in accord with this section must be in a format that is publicly
351
339
documented (and with an implementation available to the public in
352
340
source code form), and must require no special password or key for
353
341
unpacking, reading or copying.
354
342
355
343
  7. Additional Terms.
356
344
357
345
  "Additional permissions" are terms that supplement the terms of this
358
346
License by making exceptions from one or more of its conditions.
359
347
Additional permissions that are applicable to the entire Program shall
360
348
be treated as though they were included in this License, to the extent
361
349
that they are valid under applicable law.  If additional permissions
362
350
apply only to part of the Program, that part may be used separately
363
351
under those permissions, but the entire Program remains governed by
364
352
this License without regard to the additional permissions.
365
353
366
354
  When you convey a copy of a covered work, you may at your option
367
355
remove any additional permissions from that copy, or from any part of
368
356
it.  (Additional permissions may be written to require their own
369
357
removal in certain cases when you modify the work.)  You may place
370
358
additional permissions on material, added by you to a covered work,
371
359
for which you have or can give appropriate copyright permission.
372
360
373
361
  Notwithstanding any other provision of this License, for material you
374
362
add to a covered work, you may (if authorized by the copyright holders of
375
363
that material) supplement the terms of this License with terms:
376
364
377
365
    a) Disclaiming warranty or limiting liability differently from the
378
366
    terms of sections 15 and 16 of this License; or
379
367
380
368
    b) Requiring preservation of specified reasonable legal notices or
381
369
    author attributions in that material or in the Appropriate Legal
382
370
    Notices displayed by works containing it; or
383
371
384
372
    c) Prohibiting misrepresentation of the origin of that material, or
385
373
    requiring that modified versions of such material be marked in
386
374
    reasonable ways as different from the original version; or
387
375
388
376
    d) Limiting the use for publicity purposes of names of licensors or
389
377
    authors of the material; or
390
378
391
379
    e) Declining to grant rights under trademark law for use of some
392
380
    trade names, trademarks, or service marks; or
393
381
394
382
    f) Requiring indemnification of licensors and authors of that
395
383
    material by anyone who conveys the material (or modified versions of
396
384
    it) with contractual assumptions of liability to the recipient, for
397
385
    any liability that these contractual assumptions directly impose on
398
386
    those licensors and authors.
399
387
400
388
  All other non-permissive additional terms are considered "further
401
389
restrictions" within the meaning of section 10.  If the Program as you
402
390
received it, or any part of it, contains a notice stating that it is
403
391
governed by this License along with a term that is a further
404
392
restriction, you may remove that term.  If a license document contains
405
393
a further restriction but permits relicensing or conveying under this
406
394
License, you may add to a covered work material governed by the terms
407
395
of that license document, provided that the further restriction does
408
396
not survive such relicensing or conveying.
409
397
410
398
  If you add terms to a covered work in accord with this section, you
411
399
must place, in the relevant source files, a statement of the
412
400
additional terms that apply to those files, or a notice indicating
413
401
where to find the applicable terms.
414
402
415
403
  Additional terms, permissive or non-permissive, may be stated in the
416
404
form of a separately written license, or stated as exceptions;
417
405
the above requirements apply either way.
418
406
419
407
  8. Termination.
420
408
421
409
  You may not propagate or modify a covered work except as expressly
422
410
provided under this License.  Any attempt otherwise to propagate or
423
411
modify it is void, and will automatically terminate your rights under
424
412
this License (including any patent licenses granted under the third
425
413
paragraph of section 11).
426
414
427
415
  However, if you cease all violation of this License, then your
428
416
license from a particular copyright holder is reinstated (a)
429
417
provisionally, unless and until the copyright holder explicitly and
430
418
finally terminates your license, and (b) permanently, if the copyright
431
419
holder fails to notify you of the violation by some reasonable means
432
420
prior to 60 days after the cessation.
433
421
434
422
  Moreover, your license from a particular copyright holder is
435
423
reinstated permanently if the copyright holder notifies you of the
436
424
violation by some reasonable means, this is the first time you have
437
425
received notice of violation of this License (for any work) from that
438
426
copyright holder, and you cure the violation prior to 30 days after
439
427
your receipt of the notice.
440
428
441
429
  Termination of your rights under this section does not terminate the
442
430
licenses of parties who have received copies or rights from you under
443
431
this License.  If your rights have been terminated and not permanently
444
432
reinstated, you do not qualify to receive new licenses for the same
445
433
material under section 10.
446
434
447
435
  9. Acceptance Not Required for Having Copies.
448
436
449
437
  You are not required to accept this License in order to receive or
450
438
run a copy of the Program.  Ancillary propagation of a covered work
451
439
occurring solely as a consequence of using peer-to-peer transmission
452
440
to receive a copy likewise does not require acceptance.  However,
453
441
nothing other than this License grants you permission to propagate or
454
442
modify any covered work.  These actions infringe copyright if you do
455
443
not accept this License.  Therefore, by modifying or propagating a
456
444
covered work, you indicate your acceptance of this License to do so.
457
445
458
446
  10. Automatic Licensing of Downstream Recipients.
459
447
460
448
  Each time you convey a covered work, the recipient automatically
461
449
receives a license from the original licensors, to run, modify and
462
450
propagate that work, subject to this License.  You are not responsible
463
451
for enforcing compliance by third parties with this License.
464
452
465
453
  An "entity transaction" is a transaction transferring control of an
466
454
organization, or substantially all assets of one, or subdividing an
467
455
organization, or merging organizations.  If propagation of a covered
468
456
work results from an entity transaction, each party to that
469
457
transaction who receives a copy of the work also receives whatever
470
458
licenses to the work the party's predecessor in interest had or could
471
459
give under the previous paragraph, plus a right to possession of the
472
460
Corresponding Source of the work from the predecessor in interest, if
473
461
the predecessor has it or can get it with reasonable efforts.
474
462
475
463
  You may not impose any further restrictions on the exercise of the
476
464
rights granted or affirmed under this License.  For example, you may
477
465
not impose a license fee, royalty, or other charge for exercise of
478
466
rights granted under this License, and you may not initiate litigation
479
467
(including a cross-claim or counterclaim in a lawsuit) alleging that
480
468
any patent claim is infringed by making, using, selling, offering for
481
469
sale, or importing the Program or any portion of it.
482
470
483
471
  11. Patents.
484
472
485
473
  A "contributor" is a copyright holder who authorizes use under this
486
474
License of the Program or a work on which the Program is based.  The
487
475
work thus licensed is called the contributor's "contributor version".
488
476
489
477
  A contributor's "essential patent claims" are all patent claims
490
478
owned or controlled by the contributor, whether already acquired or
491
479
hereafter acquired, that would be infringed by some manner, permitted
492
480
by this License, of making, using, or selling its contributor version,
493
481
but do not include claims that would be infringed only as a
494
482
consequence of further modification of the contributor version.  For
495
483
purposes of this definition, "control" includes the right to grant
496
484
patent sublicenses in a manner consistent with the requirements of
497
485
this License.
498
486
499
487
  Each contributor grants you a non-exclusive, worldwide, royalty-free
500
488
patent license under the contributor's essential patent claims, to
501
489
make, use, sell, offer for sale, import and otherwise run, modify and
502
490
propagate the contents of its contributor version.
503
491
504
492
  In the following three paragraphs, a "patent license" is any express
505
493
agreement or commitment, however denominated, not to enforce a patent
506
494
(such as an express permission to practice a patent or covenant not to
507
495
sue for patent infringement).  To "grant" such a patent license to a
508
496
party means to make such an agreement or commitment not to enforce a
509
497
patent against the party.
510
498
511
499
  If you convey a covered work, knowingly relying on a patent license,
512
500
and the Corresponding Source of the work is not available for anyone
513
501
to copy, free of charge and under the terms of this License, through a
514
502
publicly available network server or other readily accessible means,
515
503
then you must either (1) cause the Corresponding Source to be so
516
504
available, or (2) arrange to deprive yourself of the benefit of the
517
505
patent license for this particular work, or (3) arrange, in a manner
518
506
consistent with the requirements of this License, to extend the patent
519
507
license to downstream recipients.  "Knowingly relying" means you have
520
508
actual knowledge that, but for the patent license, your conveying the
521
509
covered work in a country, or your recipient's use of the covered work
522
510
in a country, would infringe one or more identifiable patents in that
523
511
country that you have reason to believe are valid.
524
512
525
513
  If, pursuant to or in connection with a single transaction or
526
514
arrangement, you convey, or propagate by procuring conveyance of, a
527
515
covered work, and grant a patent license to some of the parties
528
516
receiving the covered work authorizing them to use, propagate, modify
529
517
or convey a specific copy of the covered work, then the patent license
530
518
you grant is automatically extended to all recipients of the covered
531
519
work and works based on it.
532
520
533
521
  A patent license is "discriminatory" if it does not include within
534
522
the scope of its coverage, prohibits the exercise of, or is
535
523
conditioned on the non-exercise of one or more of the rights that are
536
524
specifically granted under this License.  You may not convey a covered
537
525
work if you are a party to an arrangement with a third party that is
538
526
in the business of distributing software, under which you make payment
539
527
to the third party based on the extent of your activity of conveying
540
528
the work, and under which the third party grants, to any of the
541
529
parties who would receive the covered work from you, a discriminatory
542
530
patent license (a) in connection with copies of the covered work
543
531
conveyed by you (or copies made from those copies), or (b) primarily
544
532
for and in connection with specific products or compilations that
545
533
contain the covered work, unless you entered into that arrangement,
546
534
or that patent license was granted, prior to 28 March 2007.
547
535
548
536
  Nothing in this License shall be construed as excluding or limiting
549
537
any implied license or other defenses to infringement that may
550
538
otherwise be available to you under applicable patent law.
551
539
552
540
  12. No Surrender of Others' Freedom.
553
541
554
542
  If conditions are imposed on you (whether by court order, agreement or
555
543
otherwise) that contradict the conditions of this License, they do not
556
544
excuse you from the conditions of this License.  If you cannot convey a
557
545
covered work so as to satisfy simultaneously your obligations under this
558
546
License and any other pertinent obligations, then as a consequence you may
559
547
not convey it at all.  For example, if you agree to terms that obligate you
560
548
to collect a royalty for further conveying from those to whom you convey
561
549
the Program, the only way you could satisfy both those terms and this
562
550
License would be to refrain entirely from conveying the Program.
563
551
564
552
  13. Use with the GNU Affero General Public License.
565
553
566
554
  Notwithstanding any other provision of this License, you have
567
555
permission to link or combine any covered work with a work licensed
568
556
under version 3 of the GNU Affero General Public License into a single
569
557
combined work, and to convey the resulting work.  The terms of this
570
558
License will continue to apply to the part which is the covered work,
571
559
but the special requirements of the GNU Affero General Public License,
572
560
section 13, concerning interaction through a network will apply to the
573
561
combination as such.
574
562
575
563
  14. Revised Versions of this License.
576
564
577
565
  The Free Software Foundation may publish revised and/or new versions of
578
566
the GNU General Public License from time to time.  Such new versions will
579
567
be similar in spirit to the present version, but may differ in detail to
580
568
address new problems or concerns.
581
569
582
570
  Each version is given a distinguishing version number.  If the
583
571
Program specifies that a certain numbered version of the GNU General
584
572
Public License "or any later version" applies to it, you have the
585
573
option of following the terms and conditions either of that numbered
586
574
version or of any later version published by the Free Software
587
575
Foundation.  If the Program does not specify a version number of the
588
576
GNU General Public License, you may choose any version ever published
589
577
by the Free Software Foundation.
590
578
591
579
  If the Program specifies that a proxy can decide which future
592
580
versions of the GNU General Public License can be used, that proxy's
593
581
public statement of acceptance of a version permanently authorizes you
594
582
to choose that version for the Program.
595
583
596
584
  Later license versions may give you additional or different
597
585
permissions.  However, no additional obligations are imposed on any
598
586
author or copyright holder as a result of your choosing to follow a
599
587
later version.
600
588
601
589
  15. Disclaimer of Warranty.
602
590
603
591
  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
604
592
APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
605
593
HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
606
594
OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
607
595
THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
608
596
PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
609
597
IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
610
598
ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
611
599
612
600
  16. Limitation of Liability.
613
601
614
602
  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
615
603
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
616
604
THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
617
605
GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
618
606
USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
619
607
DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
620
608
PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
621
609
EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
622
610
SUCH DAMAGES.
623
611
624
612
  17. Interpretation of Sections 15 and 16.
625
613
626
614
  If the disclaimer of warranty and limitation of liability provided
627
615
above cannot be given local legal effect according to their terms,
628
616
reviewing courts shall apply local law that most closely approximates
629
617
an absolute waiver of all civil liability in connection with the
630
618
Program, unless a warranty or assumption of liability accompanies a
631
619
copy of the Program in return for a fee.
632
620
633
621
                     END OF TERMS AND CONDITIONS
634
622
635
623
            How to Apply These Terms to Your New Programs
636
624
637
625
  If you develop a new program, and you want it to be of the greatest
638
626
possible use to the public, the best way to achieve this is to make it
639
627
free software which everyone can redistribute and change under these terms.
640
628
641
629
  To do so, attach the following notices to the program.  It is safest
642
630
to attach them to the start of each source file to most effectively
643
631
state the exclusion of warranty; and each file should have at least
644
632
the "copyright" line and a pointer to where the full notice is found.
645
633
646
634
    <one line to give the program's name and a brief idea of what it does.>
647
635
    Copyright (C) <year>  <name of author>
648
636
649
637
    This program is free software: you can redistribute it and/or modify
650
638
    it under the terms of the GNU General Public License as published by
651
639
    the Free Software Foundation, either version 3 of the License, or
652
640
    (at your option) any later version.
653
641
654
642
    This program is distributed in the hope that it will be useful,
655
643
    but WITHOUT ANY WARRANTY; without even the implied warranty of
656
644
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
657
645
    GNU General Public License for more details.
658
646
659
647
    You should have received a copy of the GNU General Public License
660
648
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
661
649
662
650
Also add information on how to contact you by electronic and paper mail.
663
651
664
652
  If the program does terminal interaction, make it output a short
665
653
notice like this when it starts in an interactive mode:
666
654
667
655
    <program>  Copyright (C) <year>  <name of author>
668
656
    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
669
657
    This is free software, and you are welcome to redistribute it
670
658
    under certain conditions; type `show c' for details.
671
659
672
660
The hypothetical commands `show w' and `show c' should show the appropriate
673
661
parts of the General Public License.  Of course, your program's commands
674
662
might be different; for a GUI interface, you would use an "about box".
675
663
676
664
  You should also get your employer (if you work as a programmer) or school,
677
665
if any, to sign a "copyright disclaimer" for the program, if necessary.
678
666
For more information on this, and how to apply and follow the GNU GPL, see
679
667
<http://www.gnu.org/licenses/>.
680
668
681
669
  The GNU General Public License does not permit incorporating your program
682
670
into proprietary programs.  If your program is a subroutine library, you
683
671
may consider it more useful to permit linking proprietary applications with
684
672
the library.  If this is what you want to do, use the GNU Lesser General
685
673
Public License instead of this License.  But first, please read
686
674
<http://www.gnu.org/philosophy/why-not-lgpl.html>.
687
0
675
688
=== renamed file 'COPYING' => 'COPYING.moved'
689
=== added file 'COPYRIGHT'
690
--- COPYRIGHT	1970-01-01 00:00:00 +0000
691
+++ COPYRIGHT	2009-12-01 14:23:09 +0000
692
@@ -0,0 +1,17 @@
693
1
/* Copyright (C) 2008-2009  MSU, Phys, Group of Nanophotonics & Metamaterials
694
2
 * Copyright (C) 2009 Universiteit Gent - INTEC - IMEC
695
3
 *
696
4
 * This program is free software; you can redistribute it and/or modify
697
5
 * it under the terms of the GNU General Public License as published by
698
6
 * the Free Software Foundation; either version 2 of the License, or
699
7
 * (at your option) any later version.
700
8
 *
701
9
 * This program is distributed in the hope that it will be useful,
702
10
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
703
11
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
704
12
 * GNU General Public License for more details.
705
13
 *
706
14
 * You should have received a copy of the GNU General Public License
707
15
 * along with this program; if not, write to the Free Software
708
16
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
709
17
 */
710
0
18
711
=== renamed file 'COPYRIGHT' => 'COPYRIGHT.moved'
712
=== added file 'README'
713
--- README	1970-01-01 00:00:00 +0000
714
+++ README	2009-12-01 14:23:09 +0000
715
@@ -0,0 +1,27 @@
716
1
Meep (or MEEP) is a free finite-difference time-domain (FDTD)
717
2
simulation software package developed at MIT to model electromagnetic
718
3
systems.  You can download Meep and learn more about it at the Meep
719
4
home page:
720
5
721
6
	http://ab-initio.mit.edu/meep/
722
7
723
8
The Meep home page also links to a meep-discuss mailing list for
724
9
discussions about Meep and FDTD simulations, and a meep-announce
725
10
mailing list for announcements of Meep releases.
726
11
727
12
Python-meep is a wrapper around libmeep. It allows the scripting of
728
13
simulation with meep using Python.
729
14
730
15
Prerequisites for installation:
731
16
    - libmeep (or libmeep-mpi)  version 1.1.1
732
17
    - mpi support (if used) - tested with OpenMPI version 1.3.3
733
18
    - swig version 1.3.40 
734
19
    - gcc/g++ (required by swig)
735
20
    - numpy, scipy and company
736
21
737
22
There is a mailinglist available for Python-meep users :
738
23
	python-meep@lists.launchpad.net  
739
24
740
25
A tutorial on how to write Python-meep scripts is available in the doc/html subdirectory:
741
26
	python_meep_documentation.html
742
27
743
0
28
744
=== renamed file 'README' => 'README.moved'
745
=== added directory 'doc'
746
=== renamed directory 'doc' => 'doc.moved'
747
=== added directory 'doc/html'
748
=== added directory 'doc/html-sources'
749
=== added file 'doc/html-sources/Makefile'
750
--- doc/html-sources/Makefile	1970-01-01 00:00:00 +0000
751
+++ doc/html-sources/Makefile	2009-12-01 14:23:10 +0000
752
@@ -0,0 +1,88 @@
753
1
# Makefile for Sphinx documentation
754
2
#
755
3
756
4
# You can set these variables from the command line.
757
5
SPHINXOPTS    =
758
6
SPHINXBUILD   = sphinx-build
759
7
PAPER         =
760
8
761
9
# Internal variables.
762
10
PAPEROPT_a4     = -D latex_paper_size=a4
763
11
PAPEROPT_letter = -D latex_paper_size=letter
764
12
ALLSPHINXOPTS   = -d _build/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
765
13
766
14
.PHONY: help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
767
15
768
16
help:
769
17
	@echo "Please use \`make <target>' where <target> is one of"
770
18
	@echo "  html      to make standalone HTML files"
771
19
	@echo "  dirhtml   to make HTML files named index.html in directories"
772
20
	@echo "  pickle    to make pickle files"
773
21
	@echo "  json      to make JSON files"
774
22
	@echo "  htmlhelp  to make HTML files and a HTML help project"
775
23
	@echo "  qthelp    to make HTML files and a qthelp project"
776
24
	@echo "  latex     to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
777
25
	@echo "  changes   to make an overview of all changed/added/deprecated items"
778
26
	@echo "  linkcheck to check all external links for integrity"
779
27
	@echo "  doctest   to run all doctests embedded in the documentation (if enabled)"
780
28
781
29
clean:
782
30
	-rm -rf _build/*
783
31
784
32
html:
785
33
	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) _build/html
786
34
	@echo
787
35
	@echo "Build finished. The HTML pages are in _build/html."
788
36
789
37
dirhtml:
790
38
	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) _build/dirhtml
791
39
	@echo
792
40
	@echo "Build finished. The HTML pages are in _build/dirhtml."
793
41
794
42
pickle:
795
43
	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) _build/pickle
796
44
	@echo
797
45
	@echo "Build finished; now you can process the pickle files."
798
46
799
47
json:
800
48
	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) _build/json
801
49
	@echo
802
50
	@echo "Build finished; now you can process the JSON files."
803
51
804
52
htmlhelp:
805
53
	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) _build/htmlhelp
806
54
	@echo
807
55
	@echo "Build finished; now you can run HTML Help Workshop with the" \
808
56
	      ".hhp project file in _build/htmlhelp."
809
57
810
58
qthelp:
811
59
	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) _build/qthelp
812
60
	@echo
813
61
	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
814
62
	      ".qhcp project file in _build/qthelp, like this:"
815
63
	@echo "# qcollectiongenerator _build/qthelp/python-meep-doc.qhcp"
816
64
	@echo "To view the help file:"
817
65
	@echo "# assistant -collectionFile _build/qthelp/python-meep-doc.qhc"
818
66
819
67
latex:
820
68
	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) _build/latex
821
69
	@echo
822
70
	@echo "Build finished; the LaTeX files are in _build/latex."
823
71
	@echo "Run \`make all-pdf' or \`make all-ps' in that directory to" \
824
72
	      "run these through (pdf)latex."
825
73
826
74
changes:
827
75
	$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) _build/changes
828
76
	@echo
829
77
	@echo "The overview file is in _build/changes."
830
78
831
79
linkcheck:
832
80
	$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) _build/linkcheck
833
81
	@echo
834
82
	@echo "Link check complete; look for any errors in the above output " \
835
83
	      "or in _build/linkcheck/output.txt."
836
84
837
85
doctest:
838
86
	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) _build/doctest
839
87
	@echo "Testing of doctests in the sources finished, look at the " \
840
88
	      "results in _build/doctest/output.txt."
841
0
89
842
=== added file 'doc/html-sources/conf.py'
843
--- doc/html-sources/conf.py	1970-01-01 00:00:00 +0000
844
+++ doc/html-sources/conf.py	2009-12-01 14:23:10 +0000
845
@@ -0,0 +1,194 @@
846
1
# -*- coding: utf-8 -*-
847
2
#
848
3
# python-meep-doc documentation build configuration file, created by
849
4
# sphinx-quickstart on Fri Aug 21 10:42:04 2009.
850
5
#
851
6
# This file is execfile()d with the current directory set to its containing dir.
852
7
#
853
8
# Note that not all possible configuration values are present in this
854
9
# autogenerated file.
855
10
#
856
11
# All configuration values have a default; values that are commented out
857
12
# serve to show the default.
858
13
859
14
import sys, os
860
15
861
16
# If extensions (or modules to document with autodoc) are in another directory,
862
17
# add these directories to sys.path here. If the directory is relative to the
863
18
# documentation root, use os.path.abspath to make it absolute, like shown here.
864
19
#sys.path.append(os.path.abspath('.'))
865
20
866
21
# -- General configuration -----------------------------------------------------
867
22
868
23
# Add any Sphinx extension module names here, as strings. They can be extensions
869
24
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
870
25
extensions = []
871
26
872
27
# Add any paths that contain templates here, relative to this directory.
873
28
templates_path = ['_templates']
874
29
875
30
# The suffix of source filenames.
876
31
source_suffix = '.txt'
877
32
878
33
# The encoding of source files.
879
34
#source_encoding = 'utf-8'
880
35
881
36
# The master toctree document.
882
37
master_doc = 'python_meep_documentation'
883
38
884
39
# General information about the project.
885
40
project = u'python-meep-doc'
886
41
copyright = u'2009, Emmanuel.Lambert@intec.ugent.be'
887
42
888
43
# The version info for the project you're documenting, acts as replacement for
889
44
# |version| and |release|, also used in various other places throughout the
890
45
# built documents.
891
46
#
892
47
# The short X.Y version.
893
48
version = '0.1'
894
49
# The full version, including alpha/beta/rc tags.
895
50
release = '0.1'
896
51
897
52
# The language for content autogenerated by Sphinx. Refer to documentation
898
53
# for a list of supported languages.
899
54
#language = None
900
55
901
56
# There are two options for replacing |today|: either, you set today to some
902
57
# non-false value, then it is used:
903
58
#today = ''
904
59
# Else, today_fmt is used as the format for a strftime call.
905
60
#today_fmt = '%B %d, %Y'
906
61
907
62
# List of documents that shouldn't be included in the build.
908
63
#unused_docs = []
909
64
910
65
# List of directories, relative to source directory, that shouldn't be searched
911
66
# for source files.
912
67
exclude_trees = ['_build']
913
68
914
69
# The reST default role (used for this markup: `text`) to use for all documents.
915
70
#default_role = None
916
71
917
72
# If true, '()' will be appended to :func: etc. cross-reference text.
918
73
#add_function_parentheses = True
919
74
920
75
# If true, the current module name will be prepended to all description
921
76
# unit titles (such as .. function::).
922
77
#add_module_names = True
923
78
924
79
# If true, sectionauthor and moduleauthor directives will be shown in the
925
80
# output. They are ignored by default.
926
81
#show_authors = False
927
82
928
83
# The name of the Pygments (syntax highlighting) style to use.
929
84
pygments_style = 'sphinx'
930
85
931
86
# A list of ignored prefixes for module index sorting.
932
87
#modindex_common_prefix = []
933
88
934
89
935
90
# -- Options for HTML output ---------------------------------------------------
936
91
937
92
# The theme to use for HTML and HTML Help pages.  Major themes that come with
938
93
# Sphinx are currently 'default' and 'sphinxdoc'.
939
94
html_theme = 'default'
940
95
941
96
# Theme options are theme-specific and customize the look and feel of a theme
942
97
# further.  For a list of options available for each theme, see the
943
98
# documentation.
944
99
#html_theme_options = {}
945
100
946
101
# Add any paths that contain custom themes here, relative to this directory.
947
102
#html_theme_path = []
948
103
949
104
# The name for this set of Sphinx documents.  If None, it defaults to
950
105
# "<project> v<release> documentation".
951
106
#html_title = None
952
107
953
108
# A shorter title for the navigation bar.  Default is the same as html_title.
954
109
#html_short_title = None
955
110
956
111
# The name of an image file (relative to this directory) to place at the top
957
112
# of the sidebar.
958
113
#html_logo = None
959
114
960
115
# The name of an image file (within the static path) to use as favicon of the
961
116
# docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
962
117
# pixels large.
963
118
#html_favicon = None
964
119
965
120
# Add any paths that contain custom static files (such as style sheets) here,
966
121
# relative to this directory. They are copied after the builtin static files,
967
122
# so a file named "default.css" will overwrite the builtin "default.css".
968
123
html_static_path = ['_static']
969
124
970
125
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
971
126
# using the given strftime format.
972
127
#html_last_updated_fmt = '%b %d, %Y'
973
128
974
129
# If true, SmartyPants will be used to convert quotes and dashes to
975
130
# typographically correct entities.
976
131
#html_use_smartypants = True
977
132
978
133
# Custom sidebar templates, maps document names to template names.
979
134
#html_sidebars = {}
980
135
981
136
# Additional templates that should be rendered to pages, maps page names to
982
137
# template names.
983
138
#html_additional_pages = {}
984
139
985
140
# If false, no module index is generated.
986
141
#html_use_modindex = True
987
142
988
143
# If false, no index is generated.
989
144
#html_use_index = True
990
145
991
146
# If true, the index is split into individual pages for each letter.
992
147
#html_split_index = False
993
148
994
149
# If true, links to the reST sources are added to the pages.
995
150
#html_show_sourcelink = True
996
151
997
152
# If true, an OpenSearch description file will be output, and all pages will
998
153
# contain a <link> tag referring to it.  The value of this option must be the
999
154
# base URL from which the finished HTML is served.
1000
155
#html_use_opensearch = ''
1001
156
1002
157
# If nonempty, this is the file name suffix for HTML files (e.g. ".xhtml").
1003
158
#html_file_suffix = ''
1004
159
1005
160
# Output file base name for HTML help builder.
1006
161
htmlhelp_basename = 'python-meep-docdoc'
1007
162
1008
163
1009
164
# -- Options for LaTeX output --------------------------------------------------
1010
165
1011
166
# The paper size ('letter' or 'a4').
1012
167
#latex_paper_size = 'letter'
1013
168
1014
169
# The font size ('10pt', '11pt' or '12pt').
1015
170
#latex_font_size = '10pt'
1016
171
1017
172
# Grouping the document tree into LaTeX files. List of tuples
1018
173
# (source start file, target name, title, author, documentclass [howto/manual]).
1019
174
latex_documents = [
1020
175
  ('python_meep_documentation', 'python-meep-doc.tex', u'python-meep-doc Documentation',
1021
176
   u'Emmanuel.Lambert@intec.ugent.be', 'manual'),
1022
177
]
1023
178
1024
179
# The name of an image file (relative to this directory) to place at the top of
1025
180
# the title page.
1026
181
#latex_logo = None
1027
182
1028
183
# For "manual" documents, if this is true, then toplevel headings are parts,
1029
184
# not chapters.
1030
185
#latex_use_parts = False
1031
186
1032
187
# Additional stuff for the LaTeX preamble.
1033
188
#latex_preamble = ''
1034
189
1035
190
# Documents to append as an appendix to all manuals.
1036
191
#latex_appendices = []
1037
192
1038
193
# If false, no module index is generated.
1039
194
#latex_use_modindex = True
1040
0
195
1041
=== added directory 'doc/html-sources/images'
1042
=== added file 'doc/html-sources/images/bentwgB.gif'
1043
1
Binary files doc/html-sources/images/bentwgB.gif	1970-01-01 00:00:00 +0000 and doc/html-sources/images/bentwgB.gif	2009-12-01 14:23:10 +0000 differ
196
Binary files doc/html-sources/images/bentwgB.gif	1970-01-01 00:00:00 +0000 and doc/html-sources/images/bentwgB.gif	2009-12-01 14:23:10 +0000 differ
1044
=== added file 'doc/html-sources/images/bentwgNB.gif'
1045
2
Binary files doc/html-sources/images/bentwgNB.gif	1970-01-01 00:00:00 +0000 and doc/html-sources/images/bentwgNB.gif	2009-12-01 14:23:10 +0000 differ
197
Binary files doc/html-sources/images/bentwgNB.gif	1970-01-01 00:00:00 +0000 and doc/html-sources/images/bentwgNB.gif	2009-12-01 14:23:10 +0000 differ
1046
=== added file 'doc/html-sources/images/fluxes.png'
1047
3
Binary files doc/html-sources/images/fluxes.png	1970-01-01 00:00:00 +0000 and doc/html-sources/images/fluxes.png	2009-12-01 14:23:10 +0000 differ
198
Binary files doc/html-sources/images/fluxes.png	1970-01-01 00:00:00 +0000 and doc/html-sources/images/fluxes.png	2009-12-01 14:23:10 +0000 differ
1048
=== added file 'doc/html-sources/make.bat'
1049
--- doc/html-sources/make.bat	1970-01-01 00:00:00 +0000
1050
+++ doc/html-sources/make.bat	2009-12-01 14:23:10 +0000
1051
@@ -0,0 +1,112 @@
1052
1
@ECHO OFF
1053
2
1054
3
REM Command file for Sphinx documentation
1055
4
1056
5
set SPHINXBUILD=sphinx-build
1057
6
set ALLSPHINXOPTS=-d _build/doctrees %SPHINXOPTS% .
1058
7
if NOT "%PAPER%" == "" (
1059
8
	set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS%
1060
9
)
1061
10
1062
11
if "%1" == "" goto help
1063
12
1064
13
if "%1" == "help" (
1065
14
	:help
1066
15
	echo.Please use `make ^<target^>` where ^<target^> is one of
1067
16
	echo.  html      to make standalone HTML files
1068
17
	echo.  dirhtml   to make HTML files named index.html in directories
1069
18
	echo.  pickle    to make pickle files
1070
19
	echo.  json      to make JSON files
1071
20
	echo.  htmlhelp  to make HTML files and a HTML help project
1072
21
	echo.  qthelp    to make HTML files and a qthelp project
1073
22
	echo.  latex     to make LaTeX files, you can set PAPER=a4 or PAPER=letter
1074
23
	echo.  changes   to make an overview over all changed/added/deprecated items
1075
24
	echo.  linkcheck to check all external links for integrity
1076
25
	echo.  doctest   to run all doctests embedded in the documentation if enabled
1077
26
	goto end
1078
27
)
1079
28
1080
29
if "%1" == "clean" (
1081
30
	for /d %%i in (_build\*) do rmdir /q /s %%i
1082
31
	del /q /s _build\*
1083
32
	goto end
1084
33
)
1085
34
1086
35
if "%1" == "html" (
1087
36
	%SPHINXBUILD% -b html %ALLSPHINXOPTS% _build/html
1088
37
	echo.
1089
38
	echo.Build finished. The HTML pages are in _build/html.
1090
39
	goto end
1091
40
)
1092
41
1093
42
if "%1" == "dirhtml" (
1094
43
	%SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% _build/dirhtml
1095
44
	echo.
1096
45
	echo.Build finished. The HTML pages are in _build/dirhtml.
1097
46
	goto end
1098
47
)
1099
48
1100
49
if "%1" == "pickle" (
1101
50
	%SPHINXBUILD% -b pickle %ALLSPHINXOPTS% _build/pickle
1102
51
	echo.
1103
52
	echo.Build finished; now you can process the pickle files.
1104
53
	goto end
1105
54
)
1106
55
1107
56
if "%1" == "json" (
1108
57
	%SPHINXBUILD% -b json %ALLSPHINXOPTS% _build/json
1109
58
	echo.
1110
59
	echo.Build finished; now you can process the JSON files.
1111
60
	goto end
1112
61
)
1113
62
1114
63
if "%1" == "htmlhelp" (
1115
64
	%SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% _build/htmlhelp
1116
65
	echo.
1117
66
	echo.Build finished; now you can run HTML Help Workshop with the ^
1118
67
.hhp project file in _build/htmlhelp.
1119
68
	goto end
1120
69
)
1121
70
1122
71
if "%1" == "qthelp" (
1123
72
	%SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% _build/qthelp
1124
73
	echo.
1125
74
	echo.Build finished; now you can run "qcollectiongenerator" with the ^
1126
75
.qhcp project file in _build/qthelp, like this:
1127
76
	echo.^> qcollectiongenerator _build\qthelp\python-meep-doc.qhcp
1128
77
	echo.To view the help file:
1129
78
	echo.^> assistant -collectionFile _build\qthelp\python-meep-doc.ghc
1130
79
	goto end
1131
80
)
1132
81
1133
82
if "%1" == "latex" (
1134
83
	%SPHINXBUILD% -b latex %ALLSPHINXOPTS% _build/latex
1135
84
	echo.
1136
85
	echo.Build finished; the LaTeX files are in _build/latex.
1137
86
	goto end
1138
87
)
1139
88
1140
89
if "%1" == "changes" (
1141
90
	%SPHINXBUILD% -b changes %ALLSPHINXOPTS% _build/changes
1142
91
	echo.
1143
92
	echo.The overview file is in _build/changes.
1144
93
	goto end
1145
94
)
1146
95
1147
96
if "%1" == "linkcheck" (
1148
97
	%SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% _build/linkcheck
1149
98
	echo.
1150
99
	echo.Link check complete; look for any errors in the above output ^
1151
100
or in _build/linkcheck/output.txt.
1152
101
	goto end
1153
102
)
1154
103
1155
104
if "%1" == "doctest" (
1156
105
	%SPHINXBUILD% -b doctest %ALLSPHINXOPTS% _build/doctest
1157
106
	echo.
1158
107
	echo.Testing of doctests in the sources finished, look at the ^
1159
108
results in _build/doctest/output.txt.
1160
109
	goto end
1161
110
)
1162
111
1163
112
:end
1164
0
113
1165
=== added file 'doc/html-sources/python_meep_documentation.txt'
1166
--- doc/html-sources/python_meep_documentation.txt	1970-01-01 00:00:00 +0000
1167
+++ doc/html-sources/python_meep_documentation.txt	2009-12-01 14:23:10 +0000
1168
@@ -0,0 +1,1299 @@
1169
1
PYTHON-MEEP BINDING DOCUMENTATION
1170
2
====================================
1171
3
1172
4
Primary author of this documentation : EL : Emmanuel.Lambert@intec.ugent.be 
1173
5
1174
6
Document history : 
1175
7
1176
8
::
1177
9
1178
10
	* EL-19/20/21-08-2009 : document creation 
1179
11
	* EL-24-08-2009 : small improvements & clarifications.
1180
12
	* EL-25/26-08-2009 : sections 7 & 8 were added.
1181
13
	* EL-03-04/09/2009 : 
1182
14
               -class "structure_eps_pml" (removed again in v0.8).
1183
15
	       -port to Meep 1.1.1 (class 'volume' was renamed to 'grid_volume' and class 'geometric_volume' to 'volume'
1184
16
	       -minor changes in the bent waveguide sample, to make it more consistent with the Scheme version
1185
17
	* EL-07-08/09/2009 : sections 3.2, 8.2, 8.3 : defining a material function with inline C/C++
1186
18
	* EL-10/09/2009 : additions for MPI mode (multiprocessor)
1187
19
	* EL-21/10/2009 : amplitude factor callback function
1188
20
	* EL-22/10/2009 : keyword arguments for runUntilFieldsDecayed
1189
21
	* EL-01/12/2009 : alignment with version 0.8 - III
1190
22
	* EL-01/12/2009 : release 1.0 / added info about environment variables for inline C/C++
1191
23
1192
24
1193
25
1194
26
**1. The general structure of a python-meep program**
1195
27
-----------------------------------------------------
1196
28
1197
29
In general terms, a python-meep program can be structured as follows :
1198
30
1199
31
    * import the python-meep binding :
1200
32
		``from meep import *``
1201
33
      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/``
1202
34
1203
35
      If you are running in MPI mode (multiprocessor, see section 9), then you have to import module ``meep_mpi`` instead :
1204
36
		``from meep_mpi import *``
1205
37
1206
38
    * define a computational grid volume
1207
39
        See section 2 below which explains usage of the ``grid_volume`` class. 
1208
40
1209
41
    * define the waveguide structure (describing the geometry, PML and materials)
1210
42
        See section 3 below which explains usage of the ``structure`` class.
1211
43
		
1212
44
    * create an object which will hold the calculated fields 
1213
45
        See section 4 below which explains usage of the ``field`` class.
1214
46
1215
47
    * define the sources
1216
48
        See section 5 below which explains usage of the ``add_point_source`` and ``add_volume_source`` functions.
1217
49
1218
50
    * run the simulation (iterate over the time-steps)
1219
51
        See section 6 below.
1220
52
1221
53
Section 7 gives details about defining and retrieving fluxes.
1222
54
1223
55
Section 9 gives some complete examples.
1224
56
1225
57
Section 10 outlines some differences between Scheme-Meep and Python-Meep.
1226
58
1227
59
1228
60
**2. Defining the computational grid volume**
1229
61
---------------------------------------------
1230
62
1231
63
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).
1232
64
		* ``volcyl(rsize, zsize, resolution)``
1233
65
			Defines a cyclical computational grid volume.
1234
66
		* ``volone(zsize, resolution)`` *alternatively called* ``vol1d(zsize, resolution)``
1235
67
			Defines a 1-dimensional computational grid volume along the Z-axis.
1236
68
		* ``voltwo(xsize, ysize, resolution)`` *alternatively called* ``vol2d(xsize, ysize, a)``
1237
69
			Defines a 2-dimensional computational grid volumes along the X- and Y-axes
1238
70
		* ``vol3d(xsize, ysize, zsize, resolution)``
1239
71
			Defines a 3-dimensional computational grid volume.
1240
72
1241
73
    e.g.: ``v = volone(6, 10)`` defines a 1-dimensional computational volume of lenght 6, with 10 pixels per distance unit.
1242
74
1243
75
1244
76
**3. Defining the waveguide structure**
1245
77
---------------------------------------
1246
78
	
1247
79
The waveguide structure is defined using the class ``structure``, of which the constructor has the following arguments :
1248
80
1249
81
		* *required* : the computational grid volume (a reference to an object of type ``grid_volume``, see section 2 above)
1250
82
1251
83
		* *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.
1252
84
1253
85
		* *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 :
1254
86
			- ``no_pml()`` describing a conditionless boundary region (no PML)
1255
87
			- ``pml(thickness)`` : decribing a perfectly matching layer (PML) of a certain thickness (double value) on the boundaries in all directions. 
1256
88
			- ``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``). 
1257
89
			- ``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.
1258
90
1259
91
		* *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:
1260
92
			- ``identity`` : no symmetry
1261
93
			- ``rotate4(direction, grid_volume)`` : defines a 90° rotational symmetry with 'direction' the axis of rotation.
1262
94
			- ``rotate2(direction, grid_volume)`` : defines a 180° rotational symmetry with 'direction' the axis of rotation. 
1263
95
			- ``mirror(direction, grid_volume)`` : defines a mirror symmetry plane with 'direction' normal to the mirror plane.
1264
96
			- ``r_to_minus_r_symmetry`` : defines a mirror symmetry in polar coordinates
1265
97
1266
98
		* optional: the number of chunks in which to split up the calculated geometry. If you leave this empty, it is auto-configured. Otherwise, you would set this to a factor which is a multiple of the number of processors in your MPI run (for multiprocessor configuration).
1267
99
1268
100
    e.g. : if ``v`` is a variable pointing to the computational grid volume, then :
1269
101
		 ``s = structure(v, one)`` defines a structure with all air (eps=1), 
1270
102
    which is equivalent to:
1271
103
		 ``s = structure(v, one, no_pml(), identity(), 1)``
1272
104
1273
105
    Another example : ``s = structure(v, EPS, pml(0.1,Y) )`` with EPS a custom material function, which is explained in the note below.
1274
106
1275
107
1276
108
3.1. Defining a material function
1277
109
________________________________________
1278
110
1279
111
In order to describe the geometry of the waveguide, we have to provide a 'material function' returning 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.
1280
112
It is also possible (and faster) to write your material function in inline C/C++ (see 3.3)
1281
113
1282
114
E.g. :
1283
115
1284
116
::
1285
117
1286
118
        class epsilon(Callback):  #inherit from Callback for integration with the meep core library
1287
119
                def __init__(self):
1288
120
                        Callback.__init__(self)
1289
121
                def double_vec(self,vec): #override of function in the Callback class to set the eps function
1290
122
                        self.set_double(self.eps(vec))
1291
123
                        return 
1292
124
                def eps(self,vec):   #return the epsilon value for a certain point (indicated by the vector v)
1293
125
                        v = vec 
1294
126
                        r = v.x()*v.x() + v.y()*v.y()
1295
127
                        dr = sqrt(r)
1296
128
                        while dr>1:
1297
129
                            dr-=1
1298
130
                        if dr > 0.7001:
1299
131
                            return 12.0
1300
132
                        return 1.0
1301
133
1302
134
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.
1303
135
So, one should write : ``vec.x()`` and ``vec.y()``.
1304
136
 
1305
137
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 :
1306
138
1307
139
::
1308
140
1309
141
        set_EPS_Callback(epsilon().__disown__())
1310
142
        s = structure(v, EPS, no_pml(), identity())
1311
143
1312
144
1313
145
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.
1314
146
1315
147
***Important remark*** : at the end of our program, we should call : ``del_EPS_Callback()`` in order to clean up the global variable.
1316
148
1317
149
For better performance, you can define your EPS material function with inline C/C++ : we refer to section 3.3 for details about this.
1318
150
1319
151
3.2 Eps-averaging
1320
152
_________________
1321
153
1322
154
EPS-averaging (anisotrpic averaging) is disabled by default, making this behaviour consistent with the behaviour of the Meep C++ core.
1323
155
1324
156
You can enable EPS-averaging using the function ``use_averaging`` :
1325
157
1326
158
::
1327
159
1328
160
	#enable EPS-averaging
1329
161
	use_averaging(True)
1330
162
	...
1331
163
	#disable EPS-averaging
1332
164
	use_averaging(False)
1333
165
	...
1334
166
1335
167
1336
168
Enabling EPS-averaging results in slower performance, but more accurate results. 
1337
169
1338
170
1339
171
3.3. Defining a material function with inline C/C++
1340
172
_________________________________________________________
1341
173
1342
174
The approach described in 3.1 lets the Meep core library call back to Python for every query of the epsilon-function. This creates a lot of overhead.
1343
175
An approach which has a lot better performance is to define this epsilon-function with an inline C-function or C++ class.
1344
176
1345
177
* If our epsilon-function needs *no other parameters than the position vector (X, Y, Z)*, then we can suffice with an inline C-function (the geometry dependencies are then typically hardcoded).
1346
178
1347
179
* If our epsilon-function needs to base it's calculation on *a more complex set of parameters (e.g. parameters depending on the geometry)*, then we have to write a C++ class. 
1348
180
1349
181
For example, in the bent-waveguide example (section 8.3), we can define a generic C++ class which can return the epsilon-value for both the "bend" and "no bend" case, with variable size parameters.
1350
182
We can also take a simpler approach (section 8.2) and write a function in which the geometry size parameters are hardcoded : we then need 2 inline C-functions : one for the "bend" case and one for the "no bend" case. 
1351
183
1352
184
Make sure the following environment variables are defined :
1353
185
	* MEEP_INCLUDE : should point to the path containing meep.hpp (and a subdirectory 'meep' with vec.hpp and mympi.hpp), e.g.: ``/usr/include``
1354
186
	* MEEP_LIB : should point to the path containing libmeep.so, e.g. : ``/usr/lib``
1355
187
	* PYTHONMEEP_INCLUDE : should point to the path containing custom.hpp, e.g. : ``/usr/local/lib/python2.6/dist-packages``
1356
188
	
1357
189
1358
190
3.3.1 Inline C-function
1359
191
.......................
1360
192
1361
193
First we create a header file, e.g. "eps_function.hpp" which contains our EPS-function.
1362
194
Not that the geometry dependencies are hardcoded (``upperLimitHorizontalWg = 4`` and ``lowerLimitHorizontalWg = 5``).
1363
195
1364
196
::
1365
197
1366
198
1367
199
	namespace meep
1368
200
	{
1369
201
		static double myEps(const vec &v) {
1370
202
			double xCo = v.x();
1371
203
			double yCo = v.y();
1372
204
			double upperLimitHorizontalWg = 4;
1373
205
			double lowerLimitHorizontalWg = 5;
1374
206
			if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
1375
207
				return 1.0;
1376
208
			}
1377
209
			else return 12.0;
1378
210
		}
1379
211
	}
1380
212
1381
213
1382
214
Then, in the Python program, we prepare and set the callback function as shown below.
1383
215
``prepareCallbackCfunction`` returns a pointer to the C-function, which we deliver to the Meep core using ``set_EPS_CallbackInlineFunction``.
1384
216
1385
217
::
1386
218
1387
219
	def initEPS(isWgBent):
1388
220
		funPtr = prepareCallbackCfunction("myEps","eps_function.hpp")  #name of your function / name of header file
1389
221
		set_EPS_CallbackInlineFunction(funPtr)
1390
222
		print "EPS function successfully set."
1391
223
		return funPtr
1392
224
1393
225
We refer to section 8.2 below for a full example.
1394
226
1395
227
1396
228
3.3.2 Inline C++-class
1397
229
......................
1398
230
1399
231
A more complex approach is to have a C++ object that can accept more parameters when it is constructed. 
1400
232
For example this is the case if want to change the parameters of the geometry from inside Python without touching the C++ code.
1401
233
1402
234
We create a header file "eps_class.hpp" which contains the definition of the class (the class must inherit from ``Callback``).
1403
235
In the example below, the parameters ``upperLimitHorizontalWg`` and ``widthWg`` will be communicated from Python upon construction of the object.
1404
236
If these parameters then change (depending on the geometry), the C++ object will follow automatically.
1405
237
1406
238
1407
239
::
1408
240
1409
241
	using namespace meep;
1410
242
1411
243
	namespace meep
1412
244
	{
1413
245
1414
246
	class myEpsCallBack : virtual public Callback {
1415
247
1416
248
	    public:	
1417
249
		   myEpsCallBack() : Callback() { };
1418
250
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
1419
251
1420
252
		   myEpsCallBack(double upperLimitHorizontalWg, double widthWg)  : Callback() {
1421
253
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
1422
254
			_widthWg = widthWg;
1423
255
		   };	
1424
256
1425
257
		   double double_vec(const vec &v) { //return the EPS-value, depending on the position vector
1426
258
			double eps = myEps(v, _upperLimitHorizontalWg, _widthWg);
1427
259
			return eps;
1428
260
		   };
1429
261
		   
1430
262
		   complex<double> complex_vec(const vec &x) { return 0; }; //no need to implement
1431
263
		   complex<double> complex_time(const double &t) { return 0; }; //no need to implement 
1432
264
1433
265
1434
266
	     private:
1435
267
		   double _upperLimitHorizontalWg;
1436
268
		   double _widthWg;
1437
269
1438
270
		   double myEps(const vec &v, double upperLimitHorizontalWg, double widthWg) {
1439
271
			double xCo = v.x();
1440
272
			double yCo = v.y();
1441
273
	    	        if ((yCo < upperLimitHorizontalWg) || (yCo > upperLimitHorizontalWg+widthWg)){
1442
274
				return 1.0;
1443
275
			    }
1444
276
			}
1445
277
			return 12.0;
1446
278
		   }
1447
279
1448
280
	};
1449
281
1450
282
	}
1451
283
1452
284
1453
285
The syntax in Python is a little bit more complex in this case.
1454
286
1455
287
We will need to import the module ``scipy.weave`` : 
1456
288
1457
289
``from scipy.weave import *``
1458
290
1459
291
(this is not required for the previous case of a simple inline function)
1460
292
1461
293
First we create a list with the names of the parameters that we want to pass to the C++ class. These variables must be declared in the scope where the ``inline`` function call happens (see below).
1462
294
1463
295
``c_params = ['upperLimitHorizontalWg','widthWg']``
1464
296
1465
297
Then, we prepare the code snippet, using the function ``prepareCallbackCObjectCode`` and passing the class name and parameter names list.
1466
298
1467
299
``c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)``
1468
300
1469
301
Finally, we call the ``inline`` function, passing :
1470
302
	* the code snippet
1471
303
	* the list of parameter names
1472
304
	* the inline libraries, include directories and headers (helper functions are provided for this, see below). The call to ``getInlineHeaders`` should receive the name of the header file (with the definition of the C++ class) as an argument.
1473
305
1474
306
``funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )``
1475
307
1476
308
::
1477
309
1478
310
1479
311
	def initEPS():
1480
312
		#the set of parameters that we want to pass to the Callback object upon construction
1481
313
		#all these variables must be declared in the scope where the "inline" function call happens
1482
314
		c_params = ['upperLimitHorizontalWg','widthWg'] 
1483
315
		#the C-code snippet for constructing the Callback object
1484
316
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
1485
317
		#do the actual inline C-call and fetch the pointer to the Callback object
1486
318
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
1487
319
		#set the pointer to the callback object in the Python-meep core 
1488
320
		set_EPS_CallbackInlineObject(funPtr)
1489
321
		print "EPS function successfully set."
1490
322
		return
1491
323
 
1492
324
1493
325
We refer to section 8.3 below for a full example.
1494
326
1495
327
1496
328
1497
329
**4. Defining the initial field**
1498
330
---------------------------------
1499
331
1500
332
This is optional.
1501
333
1502
334
We create an object of type ``fields``, which will contain the calculated field.
1503
335
1504
336
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``.
1505
337
1506
338
::
1507
339
1508
340
            class fi(Callback):  #inherit from Callback for integration with the meep core library
1509
341
                def __init__(self):
1510
342
                    Callback.__init__(self)		
1511
343
                def complex_vec(self,v): #override of function in the Callback class to set the field initialization function
1512
344
		     #return the field value for a certain point (indicated by the vector v)
1513
345
                    return complex(1.0,0)
1514
346
1515
347
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 : 
1516
348
1517
349
	``set_INIF_Callback(fi().__disown__())  #link the INIF variable to the fi class``
1518
350
1519
351
We refer to section 3-note1 for more information about the function ``__disown__()``.
1520
352
1521
353
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 :
1522
354
1523
355
::
1524
356
1525
357
    f = fields(s)
1526
358
    comp = Hy
1527
359
    f.initialize_field(comp, INIF)
1528
360
1529
361
1530
362
***Important remark*** : at the end of our program, we should then call : ``del_INIF_Callback()`` in order to clean up the global variable.
1531
363
1532
364
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.
1533
365
1534
366
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. :
1535
367
1536
368
	``f.use_bloch(vec(0.0))``
1537
369
1538
370
*to be further elaborated - what is the exact meaning of the vector argument? (not well understood at this time)*
1539
371
1540
372
1541
373
1542
374
**5. Defining the sources**
1543
375
---------------------------
1544
376
1545
377
The definition of the current sources can be done through 2 functions of the ``fields`` class :
1546
378
        * ``add_point_source(component, src_time, vec, complex)``
1547
379
	* ``add_volume_source(component, src_time, volume, complex)``
1548
380
1549
381
1550
382
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).
1551
383
1552
384
For a point source, we must specify the center point of the current source using a vector (object of type ``vec``).
1553
385
1554
386
For a volume source, we must specify an object of type ``volume`` (*to be elablorated*).
1555
387
1556
388
The last argument is an overall complex amplitude number, multiplying the current source (default 1.0).
1557
389
1558
390
The following variants are available :
1559
391
        * ``add_point_source(component, double, double, double, double, vec centerpoint, complex amplitude, int is_continuous)``
1560
392
                * This is a shortcut function so that no ``src_time`` object must be created. *This function is preferably used for point sources.*
1561
393
                * The four real arguments define the central frequency, spectral width, peaktime and cutoff. 
1562
394
1563
395
        * ``add_volume_source(component, src_time, volume)``
1564
396
1565
397
        * ``add_volume_source(component, src_time, volume, AMPL)``
1566
398
                * AMPL is a built-in reference to a callback function. Such a callback function returns a factor to multiply the source amplitude with (complex value). It receives 1 parameter, i.e. a vector indicating a position RELATIVE to the CENTER of the source. See the example below.
1567
399
1568
400
1569
401
Three classes, inheriting from ``src_time``, are predefined and can be used off the shelf :
1570
402
	* ``gaussian_src_time`` for a Gaussian-pulse source. The constructor demands 2 arguments of type double :
1571
403
                    * the center frequency ω, in units of 2πc 
1572
404
		    * the frequency width w used in the Gaussian
1573
405
	* ``continuous_src_time`` for a continuous-wave source proportional to exp(−iωt). The constructor demands 4 arguments : 
1574
406
                    * the frequency ω, in units 2πc/distance (complex number)
1575
407
		    * the temporal width of smoothing (default 0)
1576
408
                    * the start time (default 0)
1577
409
                    * the end time (default infinity = never turn off)
1578
410
	* ``custom_src_time`` for a user-specified source function f(t) with start/end times. The constructor demands 4 arguments :
1579
411
		    * The function f(t) specifying the time-dependence of the source
1580
412
                    * *...(2nd argument unclear, to be further elaborated)...*
1581
413
                    * the start time of the source (default -infinity)
1582
414
                    * the end time of the source (default +infinity)    
1583
415
1584
416
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 :
1585
417
           
1586
418
::
1587
419
1588
420
	#define a continuous source 
1589
421
	srcFreq = 0.125
1590
422
	srcWidth = 20 
1591
423
	src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
1592
424
	srcComp = Ez
1593
425
	#make it a line source of size 1 starting on position(6,3)
1594
426
	srcGeo = volume(vec(6,3),vec(6,4))
1595
427
	f.add_volume_source(srcComp, src, srcGeo)
1596
428
1597
429
1598
430
Here is an example of the implementation of a callback function for the amplitude factor :
1599
431
1600
432
::
1601
433
1602
434
	class amplitudeFactor(Callback):
1603
435
	    def __init__(self):
1604
436
		Callback.__init__(self)
1605
437
		master_printf("Callback function for amplitude factor activated.\n")
1606
438
1607
439
	    def complex_vec(self,vec):
1608
440
		#BEWARE, these are coordinates RELATIVE to the source center !!!!
1609
441
		x = vec.x()
1610
442
		y = vec.y()
1611
443
		master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
1612
444
		result = complex(1.0,0)
1613
445
		return result
1614
446
1615
447
	...
1616
448
		#define a continuous source 
1617
449
		srcFreq = 0.125
1618
450
		srcWidth = 20 
1619
451
		src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
1620
452
		srcComp = Ez
1621
453
		#make it a line source of size 1 starting on position(6,3)
1622
454
		srcGeo = volume(vec(6,3),vec(6,4))
1623
455
		#create callback object for amplitude factor
1624
456
		af = amplitudeFactor()
1625
457
		set_AMPL_Callback(af.__disown__())
1626
458
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)
1627
459
1628
460
1629
461
**6. Running the simulation, retrieving field values and exporting HDF5 files**
1630
462
-------------------------------------------------------------------------------
1631
463
1632
464
We can now time-step and retrieve various field values along the way.
1633
465
The actual time step value can be retrieved or set through the variable ``f.dt``.
1634
466
1635
467
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.
1636
468
1637
469
To trigger a step in time, you call the function ``f.step()``.
1638
470
1639
471
To step until the source has fully decayed :
1640
472
1641
473
::
1642
474
1643
475
        while (f.time() < f.last_source_time()):
1644
476
             f.step()
1645
477
1646
478
The function ``runUntilFieldsDecayed`` mimicks the behaviour of 'stop-when-fields-decayed' in Meep-Scheme.
1647
479
This will run time steps until the source has decayed to 0.001 of the peak amplitude. After that, by default an additional 50 time steps will be run.
1648
480
The function has 7 arguments, of which 4 are mandatory and 3 are optional keywords arguments :
1649
481
* 4 regular arguments : reference to the field, reference to the computational grid volume, the source component, the monitor point.
1650
482
* keyword argument ``pHDF5OutputFile`` : reference to a HDF5 file (constructed with the function ``prepareHDF5File``); default : None (no ouput to files)
1651
483
* keyword argument ``pH5OutputIntervalSteps`` : step interval for output to HDF5 (default : 50)
1652
484
* keyword argument ``pDecayedStopAfterSteps`` : the number of steps to continue after the source has decayed to 0.001 of the peak amplitude at the probing point (default: 50)
1653
485
 
1654
486
We further refer to section 8 below where this function is applied in an example.
1655
487
1656
488
A rich functionality is available for retrieving field information. Some examples :
1657
489
1658
490
	* ``f.energy_in_box(v.surroundings())``
1659
491
	* ``f.electric_energy_in_box(v.surroundings())``
1660
492
	* ``f.magnetic_energy_in_box(v.surroundings())``
1661
493
	* ``f.thermo_energy_in_box(v.surroundings())``
1662
494
	* ``f.total_energy()``
1663
495
	* ``f.field_energy_in_box(v.surroundings())``
1664
496
	* ``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``)
1665
497
	* ``f.flux_in_box(X, v.surroundings())`` where the first argument is the direction (``X, Y, Z, R`` or ``P``)
1666
498
1667
499
We can probe the field at certain points by defining a *monitor point* as follows : 
1668
500
1669
501
::
1670
502
1671
503
        m = monitor_point()
1672
504
        p = vec(2.10) #vector identifying the point that we want to probe
1673
505
        f.get_point(m, p)
1674
506
        m.get_component(Hx)
1675
507
1676
508
We can export the dielectric function and e.g. the Ex component of the field to HDF5 files as follows :
1677
509
1678
510
::
1679
511
1680
512
        #make sure you start your Python session with 'sudo' or write rights to the current path
1681
513
        feps = prepareHDF5File("eps.h5")
1682
514
        f.output_hdf5(Dielectric, v.surroundings(), feps)   #export the Dielectric structure so that we can visually verify it 
1683
515
	fex = prepareHDF5File("ex.h5")
1684
516
	while (f.time() < f.last_source_time()):
1685
517
             f.step()
1686
518
             f.output_hdf5(Ex, v.surroundings(), fex, 1) #export the Ex component, appending to the file "ex.h5"
1687
519
1688
520
1689
521
1690
522
**7. Defining and retrieving fluxes**
1691
523
--------------------------------------
1692
524
1693
525
First we define a flux plane. 
1694
526
This is done through the creation of an object of type ``volume`` (specifying 2 vectors as arguments).
1695
527
1696
528
Then we apply this flux plane to the field, specifying 4 parameters :
1697
529
	* the reference to the ``volume`` object
1698
530
	* the minimum frequency (in the example below, this is ``srcFreqCenter-(srcPulseWidth/2.0)``)
1699
531
	* the maximum frequency (in the example below this is ``srcFreqCenter+(srcPulseWidth/2.0)`` )
1700
532
	* the number of discrete frequencies that we want to monitor in the flux (in the example below, this is 100).
1701
533
1702
534
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.
1703
535
1704
536
E.g., if ``f`` is the field, then we proceed as follows : 
1705
537
1706
538
::
1707
539
1708
540
	#define the flux plane and flux parameters
1709
541
	fluxplane = volume(vec(1,2),vec(1,3))
1710
542
	flux = f.add_dft_flux_plane(fluxplane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
1711
543
1712
544
	#run the calculation 
1713
545
	while (f.time() < f.last_source_time()):
1714
546
	    f.step()
1715
547
1716
548
	#retrieve the flux data
1717
549
	fluxdata = getFluxData(flux)
1718
550
	frequencies = fluxdata[0]
1719
551
	fluxvalues = fluxdata[1]
1720
552
1721
553
1722
554
1723
555
**8. The "90 degree bent waveguide example in Python**
1724
556
------------------------------------------------------
1725
557
1726
558
We have ported the "90 degree bent waveguide" example from the Meep-Scheme tutorial to Python.
1727
559
1728
560
The original example can be found on the following URL : http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial
1729
561
(section 'A 90° bend').
1730
562
1731
563
You can find the source code below and also in the ``/samples/bent_waveguide`` directory of the Python-Meep distribution. 
1732
564
1733
565
Sections 8.2 and 8.3 contain the same example but with the EPS-callback function as inline C-function and inline C++-class.
1734
566
1735
567
The following animated gifs can be produced from the HDF5-files (see the script included in directory 'samples') :
1736
568
1737
569
.. image:: images/bentwgNB.gif
1738
570
1739
571
.. image:: images/bentwgB.gif
1740
572
1741
573
1742
574
And here is the graph of the transmission, reflection and loss fluxes, showing the same results as the example in Scheme:
1743
575
1744
576
.. image:: images/fluxes.png
1745
577
   :height: 315
1746
578
   :width: 443
1747
579
1748
580
1749
581
8.1 With a Python class as EPS-function
1750
582
________________________________________
1751
583
1752
584
1753
585
A bottleneck in this version is the epsilon-function, which is written in Python. 
1754
586
This means that the Meep core library must do a callback to the Python function, which creates a lot of overhead.
1755
587
An approach which has a much better performance is to write this EPS-function in C : the Meep core library can then directly call back to a C-function.
1756
588
These approaches are described in 8.2 and 8.3.
1757
589
1758
590
::
1759
591
1760
592
	from meep import *
1761
593
	from math import *
1762
594
	from python_meep import *
1763
595
	import numpy 
1764
596
	import matplotlib.pyplot
1765
597
	import sys
1766
598
1767
599
	res = 10.0
1768
600
	gridSizeX = 16.0
1769
601
	gridSizeY = 32.0
1770
602
	padSize = 4.0
1771
603
	wgLengthX = gridSizeX - padSize
1772
604
	wgLengthY = gridSizeY - padSize
1773
605
	wgWidth = 1.0 #width of the waveguide
1774
606
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
1775
607
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
1776
608
	srcFreqCenter = 0.15 #gaussian source center frequency
1777
609
	srcPulseWidth = 0.1 #gaussian source pulse width
1778
610
	srcComp = Ez #gaussian source component
1779
611
1780
612
	#this function plots the waveguide material as a function of a vector(X,Y) 
1781
613
	class epsilon(Callback):
1782
614
	    def __init__(self, pIsWgBent):
1783
615
		Callback.__init__(self)
1784
616
		self.isWgBent = pIsWgBent
1785
617
	    def double_vec(self,vec):
1786
618
		if (self.isWgBent): #there is a bend
1787
619
		    if ((vec.x()<wgLengthX) and (vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
1788
620
		        return 12.0
1789
621
		    elif ((vec.x()>=wgLengthX-wgWidth) and (vec.x()<=wgLengthX) and vec.y()>= padSize ):
1790
622
		        return 12.0
1791
623
		    else:
1792
624
		        return 1.0
1793
625
		else: #there is no bend 
1794
626
		    if ((vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
1795
627
		        return 12.0
1796
628
		    else:
1797
629
		        return 1.0
1798
630
1799
631
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
1800
632
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
1801
633
		#in all directions,
1802
634
		#using the material function EPS
1803
635
		material = epsilon(pIsWgBent)
1804
636
		set_EPS_Callback(material.__disown__())
1805
637
		s = structure(pCompVol, EPS, pml(1.0) )     
1806
638
		f = fields(s)      
1807
639
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
1808
640
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
1809
641
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
1810
642
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
1811
643
		print "Field created..."
1812
644
		return f
1813
645
	       
1814
646
1815
647
	#FIRST WE WORK OUT THE CASE WITH NO BEND
1816
648
	#----------------------------------------------------------------
1817
649
	print "*1* Starting the case with no bend..."
1818
650
	#create the computational grid 
1819
651
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
1820
652
1821
653
	#create the field
1822
654
	wgBent = 0 #no bend
1823
655
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
1824
656
1825
657
	bendFnEps = "./bentwgNB_Eps.h5"
1826
658
	bendFnEz = "./bentwgNB_Ez.h5"
1827
659
	#export the dielectric structure (so that we can visually verify the waveguide structure)
1828
660
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
1829
661
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
1830
662
1831
663
	#create the file for the field components 
1832
664
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
1833
665
1834
666
	#define the flux plane for the reflected flux
1835
667
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
1836
668
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
1837
669
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
1838
670
1839
671
	#define the flux plane for the transmitted flux
1840
672
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
1841
673
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
1842
674
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
1843
675
1844
676
	print "Calculating..."
1845
677
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
1846
678
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
1847
679
	print "Done..!"
1848
680
1849
681
	#construct 2-dimensional array with the flux plane data
1850
682
	#see python_meep.py 
1851
683
	noBendReflFlux = getFluxData(noBendReflectedFlux)
1852
684
	noBendTransmFlux = getFluxData(noBendTransmFlux)
1853
685
1854
686
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
1855
687
	noBendReflectedFlux.scale_dfts(-1);
1856
688
	f = open("minusflux.h5", 'w') #truncate file if already exists
1857
689
	f.close()   
1858
690
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
1859
691
1860
692
	del_EPS_Callback()
1861
693
1862
694
1863
695
	#AND SECONDLY FOR THE CASE WITH BEND
1864
696
	#----------------------------------------------------------------
1865
697
	print "*2* Starting the case with bend..."
1866
698
	#create the computational grid 
1867
699
	bendVol = voltwo(gridSizeX,gridSizeY,res)
1868
700
1869
701
	#create the field
1870
702
	wgBent = 1 #there is a bend
1871
703
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
1872
704
1873
705
	#export the dielectric structure (so that we can visually verify the waveguide structure)
1874
706
	bendFnEps = "./bentwgB_Eps.h5"
1875
707
	bendFnEz = "./bentwgB_Ez.h5"
1876
708
	bendDielectricFile = prepareHDF5File(bendFnEps)
1877
709
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
1878
710
1879
711
	#create the file for the field components 
1880
712
	bendFileOutputEz = prepareHDF5File(bendFnEz)
1881
713
1882
714
	#define the flux plane for the reflected flux
1883
715
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
1884
716
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
1885
717
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
1886
718
1887
719
	#load the minused reflection flux from the "no bend" case as initalization
1888
720
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
1889
721
1890
722
1891
723
	#define the flux plane for the transmitted flux
1892
724
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
1893
725
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
1894
726
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
1895
727
1896
728
	print "Calculating..."
1897
729
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
1898
730
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
1899
731
	print "Done..!"
1900
732
1901
733
	#construct 2-dimensional array with the flux plane data
1902
734
	#see python_meep.py 
1903
735
	bendReflFlux = getFluxData(bendReflectedFlux)
1904
736
	bendTransmFlux = getFluxData(bendTransmFlux)
1905
737
1906
738
	del_EPS_Callback()
1907
739
1908
740
	#SHOW THE RESULTS IN A PLOT
1909
741
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
1910
742
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
1911
743
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
1912
744
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
1913
745
		                  
1914
746
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
1915
747
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
1916
748
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
1917
749
1918
750
	matplotlib.pyplot.show()
1919
751
1920
752
1921
753
8.2 With an inline C-function as EPS-function
1922
754
______________________________________________
1923
755
1924
756
The header file "eps_function.hpp" :
1925
757
1926
758
::
1927
759
1928
760
	using namespace meep;
1929
761
1930
762
	namespace meep
1931
763
	{
1932
764
	static double myEps(const vec &v, bool isWgBent) {
1933
765
		double xCo = v.x();
1934
766
		double yCo = v.y();
1935
767
		double upperLimitHorizontalWg = 4;
1936
768
		double wgLengthX = 12;
1937
769
		double leftLimitVerticalWg = 11;
1938
770
		double lowerLimitHorizontalWg = 5;
1939
771
		if (isWgBent){ //there is a bend
1940
772
		    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
1941
773
			return 1.0;            
1942
774
		    }
1943
775
		    else {
1944
776
			if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
1945
777
			     return 1.0;
1946
778
			}
1947
779
			else {
1948
780
			     return 12.0;
1949
781
			}
1950
782
		    }
1951
783
		}
1952
784
		else { //there is no bend 
1953
785
		    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
1954
786
			return 1.0;
1955
787
		    }
1956
788
		}
1957
789
		return 12.0;
1958
790
	   }
1959
791
1960
792
	static double myEpsBentWg(const vec &v) {
1961
793
		return myEps(v, true);
1962
794
	   }
1963
795
1964
796
	static double myEpsStraightWg(const vec &v) {
1965
797
		return myEps(v, false);
1966
798
	   }
1967
799
	}
1968
800
1969
801
1970
802
1971
803
And the actual Python program :
1972
804
1973
805
1974
806
::
1975
807
1976
808
1977
809
	from meep import *  
1978
810
1979
811
	from math import *
1980
812
	import numpy 
1981
813
	import matplotlib.pyplot
1982
814
	import sys
1983
815
1984
816
	res = 10.0
1985
817
	gridSizeX = 16.0
1986
818
	gridSizeY = 32.0
1987
819
	padSize = 4.0
1988
820
	wgLengthX = gridSizeX - padSize
1989
821
	wgLengthY = gridSizeY - padSize
1990
822
	wgWidth = 1.0 #width of the waveguide
1991
823
	upperLimitHorizontalWg = padSize
1992
824
	lowerLimitHorizontalWg = padSize+wgWidth
1993
825
	leftLimitVerticalWg = wgLengthX-wgWidth
1994
826
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
1995
827
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
1996
828
	srcFreqCenter = 0.15 #gaussian source center frequency
1997
829
	srcPulseWidth = 0.1 #gaussian source pulse width
1998
830
	srcComp = Ez #gaussian source component
1999
831
2000
832
	def initEPS(isWgBent):
2001
833
		if (isWgBent):
2002
834
		    funPtr = prepareCallbackCfunction("myEpsBentWg","eps_function.hpp")
2003
835
		else:
2004
836
		    funPtr = prepareCallbackCfunction("myEpsStraightWg","eps_function.hpp")
2005
837
		set_EPS_CallbackInlineFunction(funPtr)
2006
838
		print "EPS function successfully set."
2007
839
		return funPtr
2008
840
	 
2009
841
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
2010
842
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
2011
843
		#in all directions,
2012
844
		#using the material function EPS
2013
845
		s = structure(pCompVol, EPS, pml(1.0) )     
2014
846
		f = fields(s)      
2015
847
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
2016
848
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
2017
849
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
2018
850
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
2019
851
		print "Field created..."
2020
852
		return f
2021
853
2022
854
2023
855
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C-FUNCTION FOR EPS\n")
2024
856
	       
2025
857
	master_printf("Running on %d processor(s)...\n",count_processors())
2026
858
2027
859
	#FIRST WE WORK OUT THE CASE WITH NO BEND
2028
860
	#----------------------------------------------------------------
2029
861
	master_printf("*1* Starting the case with no bend...")
2030
862
2031
863
	#set EPS material function
2032
864
	initEPS(0)
2033
865
2034
866
	#create the computational grid 
2035
867
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
2036
868
2037
869
	#create the field
2038
870
	wgBent = 0 #no bend
2039
871
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
2040
872
2041
873
	bendFnEps = "./bentwgNB_Eps.h5"
2042
874
	bendFnEz = "./bentwgNB_Ez.h5"
2043
875
	#export the dielectric structure (so that we can visually verify the waveguide structure)
2044
876
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
2045
877
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
2046
878
2047
879
	#create the file for the field components 
2048
880
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
2049
881
2050
882
	#define the flux plane for the reflected flux
2051
883
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
2052
884
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
2053
885
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
2054
886
2055
887
	#define the flux plane for the transmitted flux
2056
888
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
2057
889
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
2058
890
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
2059
891
2060
892
	master_printf("Calculating...")
2061
893
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
2062
894
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
2063
895
	master_printf("Done..!")
2064
896
2065
897
	#construct 2-dimensional array with the flux plane data
2066
898
	#see python_meep.py 
2067
899
	noBendReflFlux = getFluxData(noBendReflectedFlux)
2068
900
	noBendTransmFlux = getFluxData(noBendTransmFlux)
2069
901
2070
902
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
2071
903
	noBendReflectedFlux.scale_dfts(-1);
2072
904
	f = open("minusflux.h5", 'w') #truncate file if already exists
2073
905
	f.close()   
2074
906
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
2075
907
2076
908
	del_EPS_Callback() #destruct the inline-created object
2077
909
2078
910
2079
911
	#AND SECONDLY FOR THE CASE WITH BEND
2080
912
	#----------------------------------------------------------------
2081
913
	master_printf("*2* Starting the case with bend...")
2082
914
2083
915
	#set EPS material function
2084
916
	initEPS(1)
2085
917
2086
918
	#create the computational grid 
2087
919
	bendVol = voltwo(gridSizeX,gridSizeY,res)
2088
920
2089
921
	#create the field
2090
922
	wgBent = 1 #there is a bend
2091
923
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
2092
924
2093
925
	#export the dielectric structure (so that we can visually verify the waveguide structure)
2094
926
	bendFnEps = "./bentwgB_Eps.h5"
2095
927
	bendFnEz = "./bentwgB_Ez.h5"
2096
928
	bendDielectricFile = prepareHDF5File(bendFnEps)
2097
929
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
2098
930
2099
931
	#create the file for the field components 
2100
932
	bendFileOutputEz = prepareHDF5File(bendFnEz)
2101
933
2102
934
	#define the flux plane for the reflected flux
2103
935
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
2104
936
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
2105
937
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
2106
938
2107
939
	#load the minused reflection flux from the "no bend" case as initalization
2108
940
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
2109
941
2110
942
2111
943
	#define the flux plane for the transmitted flux
2112
944
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
2113
945
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
2114
946
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
2115
947
2116
948
	master_printf("Calculating...")
2117
949
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
2118
950
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
2119
951
	master_printf("Done..!")
2120
952
2121
953
	#construct 2-dimensional array with the flux plane data
2122
954
	#see python_meep.py 
2123
955
	bendReflFlux = getFluxData(bendReflectedFlux)
2124
956
	bendTransmFlux = getFluxData(bendTransmFlux)
2125
957
2126
958
	del_EPS_Callback()
2127
959
2128
960
	#SHOW THE RESULTS IN A PLOT
2129
961
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
2130
962
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
2131
963
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
2132
964
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
2133
965
		                  
2134
966
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
2135
967
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
2136
968
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
2137
969
2138
970
	matplotlib.pyplot.show()
2139
971
2140
972
2141
973
2142
974
8.3 With an inline C++ class as EPS-function
2143
975
______________________________________________
2144
976
2145
977
2146
978
The header file "eps_class.hpp" :
2147
979
2148
980
2149
981
::
2150
982
2151
983
2152
984
	using namespace meep;
2153
985
2154
986
	namespace meep
2155
987
	{
2156
988
2157
989
	class myEpsCallBack : virtual public Callback {
2158
990
2159
991
	    public:	
2160
992
		   myEpsCallBack() : Callback() { };
2161
993
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
2162
994
2163
995
		   myEpsCallBack(bool isWgBent,double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX)  : Callback() {
2164
996
		    	_IsWgBent = isWgBent;
2165
997
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
2166
998
			_leftLimitVerticalWg = leftLimitVerticalWg;
2167
999
			_lowerLimitHorizontalWg = lowerLimitHorizontalWg;
2168
1000
			_wgLengthX = wgLengthX;
2169
1001
		   };	
2170
1002
2171
1003
		   double double_vec(const vec &v) {
2172
1004
			double eps = myEps(v, _IsWgBent, _upperLimitHorizontalWg, _leftLimitVerticalWg, _lowerLimitHorizontalWg, _wgLengthX);
2173
1005
			//cout << "X="<<v.x()<<"--Y="<<v.y()<<"--eps="<<eps<<"-"<<_upperLimitHorizontalWg<<"--"<<_leftLimitVerticalWg<<"--"<<_lowerLimitHorizontalWg<<"--"<<_wgLengthX;
2174
1006
			return eps;
2175
1007
		   };
2176
1008
		   
2177
1009
		   complex<double> complex_vec(const vec &x) { return 0; };
2178
1010
		   complex<double> complex_time(const double &t) { return 0; };   
2179
1011
2180
1012
2181
1013
	     private:
2182
1014
		   bool _IsWgBent;;
2183
1015
		   double _upperLimitHorizontalWg;
2184
1016
		   double _leftLimitVerticalWg;
2185
1017
		   double _lowerLimitHorizontalWg;
2186
1018
		   double _wgLengthX;
2187
1019
2188
1020
		   double myEps(const vec &v, bool isWgBent, double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX) {
2189
1021
			double xCo = v.x();
2190
1022
			double yCo = v.y();
2191
1023
			if (isWgBent){ //there is a bend
2192
1024
			    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
2193
1025
				return 1.0;            
2194
1026
			    }
2195
1027
			    else {
2196
1028
				if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
2197
1029
				     return 1.0;
2198
1030
				}
2199
1031
				else {
2200
1032
				     return 12.0;
2201
1033
				}
2202
1034
			    }
2203
1035
			}
2204
1036
			else { //there is no bend 
2205
1037
			    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
2206
1038
				return 1.0;
2207
1039
			    }
2208
1040
			}
2209
1041
			return 12.0;
2210
1042
		   }
2211
1043
2212
1044
	};
2213
1045
2214
1046
	}
2215
1047
2216
1048
2217
1049
The Python program :
2218
1050
2219
1051
2220
1052
::
2221
1053
2222
1054
2223
1055
	from meep import *  
2224
1056
2225
1057
	from math import *
2226
1058
	import numpy 
2227
1059
	import matplotlib.pyplot
2228
1060
	import sys
2229
1061
2230
1062
	from scipy.weave import *
2231
1063
2232
1064
	res = 10.0
2233
1065
	gridSizeX = 16.0
2234
1066
	gridSizeY = 32.0
2235
1067
	padSize = 4.0
2236
1068
	wgLengthX = gridSizeX - padSize
2237
1069
	wgLengthY = gridSizeY - padSize
2238
1070
	wgWidth = 1.0 #width of the waveguide
2239
1071
	upperLimitHorizontalWg = padSize
2240
1072
	lowerLimitHorizontalWg = padSize+wgWidth
2241
1073
	leftLimitVerticalWg = wgLengthX-wgWidth
2242
1074
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
2243
1075
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
2244
1076
	srcFreqCenter = 0.15 #gaussian source center frequency
2245
1077
	srcPulseWidth = 0.1 #gaussian source pulse width
2246
1078
	srcComp = Ez #gaussian source component
2247
1079
2248
1080
2249
1081
	def initEPS():
2250
1082
		#the set of parameters that we want to pass to the Callback object upon construction
2251
1083
		c_params = ['isWgBent','upperLimitHorizontalWg','leftLimitVerticalWg','lowerLimitHorizontalWg','wgLengthX'] #all these variables must be globally declared in the scope where the "inline" function call happens
2252
1084
		#the C-code snippet for constructing the Callback object
2253
1085
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
2254
1086
		#do the actual inline C-call and fetch the pointer to the Callback object
2255
1087
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
2256
1088
		#set the pointer to the callback object in the Python-meep core 
2257
1089
		set_EPS_CallbackInlineObject(funPtr)
2258
1090
		print "EPS function successfully set."
2259
1091
		return
2260
1092
	 
2261
1093
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
2262
1094
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
2263
1095
		#in all directions,
2264
1096
		#using the material function EPS
2265
1097
		s = structure(pCompVol, EPS, pml(1.0) )     
2266
1098
		f = fields(s)      
2267
1099
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
2268
1100
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
2269
1101
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
2270
1102
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
2271
1103
		print "Field created..."
2272
1104
		return f
2273
1105
	       
2274
1106
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C++ CLASS FOR EPS\n")
2275
1107
2276
1108
	master_printf("Running on %d processor(s)...\n",count_processors())
2277
1109
2278
1110
	#FIRST WE WORK OUT THE CASE WITH NO BEND
2279
1111
	#----------------------------------------------------------------
2280
1112
	master_printf("*1* Starting the case with no bend...")
2281
1113
2282
1114
	#set EPS material function
2283
1115
	isWgBent = 0
2284
1116
	initEPS()
2285
1117
2286
1118
	#create the computational grid 
2287
1119
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
2288
1120
2289
1121
	#create the field
2290
1122
	wgBent = 0 #no bend
2291
1123
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
2292
1124
2293
1125
	bendFnEps = "./bentwgNB_Eps.h5"
2294
1126
	bendFnEz = "./bentwgNB_Ez.h5"
2295
1127
	#export the dielectric structure (so that we can visually verify the waveguide structure)
2296
1128
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
2297
1129
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
2298
1130
2299
1131
	#create the file for the field components 
2300
1132
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
2301
1133
2302
1134
	#define the flux plane for the reflected flux
2303
1135
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
2304
1136
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
2305
1137
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
2306
1138
2307
1139
	#define the flux plane for the transmitted flux
2308
1140
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
2309
1141
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
2310
1142
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
2311
1143
2312
1144
	master_printf("Calculating...")
2313
1145
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
2314
1146
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
2315
1147
	master_printf("Done..!")
2316
1148
2317
1149
	#construct 2-dimensional array with the flux plane data
2318
1150
	#see python_meep.py 
2319
1151
	noBendReflFlux = getFluxData(noBendReflectedFlux)
2320
1152
	noBendTransmFlux = getFluxData(noBendTransmFlux)
2321
1153
2322
1154
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
2323
1155
	noBendReflectedFlux.scale_dfts(-1);
2324
1156
	f = open("minusflux.h5", 'w') #truncate file if already exists
2325
1157
	f.close()   
2326
1158
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
2327
1159
2328
1160
	del_EPS_Callback() #destruct the inline-created object
2329
1161
2330
1162
2331
1163
	#AND SECONDLY FOR THE CASE WITH BEND
2332
1164
	#----------------------------------------------------------------
2333
1165
	master_printf("*2* Starting the case with bend...")
2334
1166
2335
1167
	#set EPS material function
2336
1168
	isWgBent = 1
2337
1169
	initEPS()
2338
1170
2339
1171
	#create the computational grid 
2340
1172
	bendVol = voltwo(gridSizeX,gridSizeY,res)
2341
1173
2342
1174
	#create the field
2343
1175
	wgBent = 1 #there is a bend
2344
1176
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
2345
1177
2346
1178
	#export the dielectric structure (so that we can visually verify the waveguide structure)
2347
1179
	bendFnEps = "./bentwgB_Eps.h5"
2348
1180
	bendFnEz = "./bentwgB_Ez.h5"
2349
1181
	bendDielectricFile = prepareHDF5File(bendFnEps)
2350
1182
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
2351
1183
2352
1184
	#create the file for the field components 
2353
1185
	bendFileOutputEz = prepareHDF5File(bendFnEz)
2354
1186
2355
1187
	#define the flux plane for the reflected flux
2356
1188
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
2357
1189
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
2358
1190
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
2359
1191
2360
1192
	#load the minused reflection flux from the "no bend" case as initalization
2361
1193
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
2362
1194
2363
1195
2364
1196
	#define the flux plane for the transmitted flux
2365
1197
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
2366
1198
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
2367
1199
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
2368
1200
2369
1201
	master_printf("Calculating...")
2370
1202
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
2371
1203
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
2372
1204
	master_printf("Done..!")
2373
1205
2374
1206
	#construct 2-dimensional array with the flux plane data
2375
1207
	#see python_meep.py 
2376
1208
	bendReflFlux = getFluxData(bendReflectedFlux)
2377
1209
	bendTransmFlux = getFluxData(bendTransmFlux)
2378
1210
2379
1211
	del_EPS_Callback()
2380
1212
2381
1213
	#SHOW THE RESULTS IN A PLOT
2382
1214
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
2383
1215
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
2384
1216
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
2385
1217
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
2386
1218
		                  
2387
1219
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
2388
1220
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
2389
1221
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
2390
1222
2391
1223
	matplotlib.pyplot.show()
2392
1224
2393
1225
2394
1226
2395
1227
**9. Running in MPI mode (multiprocessor configuration)**
2396
1228
----------------------------------------------------------
2397
1229
2398
1230
* We assume that an MPI implementation is installed on your machine (e.g. OpenMPI).
2399
1231
If you want to run python-meep in MPI mode, then you must must import the ``meep_mpi`` module instead of the ``meep`` module : 
2400
1232
2401
1233
``from meep_mpi import *``
2402
1234
2403
1235
Then start up the Python script as follows :
2404
1236
2405
1237
``mpirun -n 2 ./myscript.py``
2406
1238
2407
1239
The ``-n`` parameter indicates the number of processors requested.
2408
1240
2409
1241
* For printing output to the console, use the ``master_printf`` statement. This will generate output on the master node only (regular Python ``print`` statements will run on all nodes).
2410
1242
2411
1243
* If you output HDF5 files, make sure your HDF5 library is MPI-enable, otherwise your Python script will stall upon writing HDF5 due to a deadlock.
2412
1244
2413
1245
2414
1246
2415
1247
**10. Differences between Python-Meep and Scheme-Meep (libctl)**
2416
1248
------------------------------------------------------------------
2417
1249
2418
1250
**note**: this section does NOT apply to the UGent Intec Photonics Research Group (apart from the coordinate system, the default behaviour for us is made consistent with Scheme-Meep) 
2419
1251
2420
1252
The general rule is that Python-Meep has consistent behaviour with the **C++ core of Meep**.
2421
1253
The default version of Python-Meep (compiled from the LATEST_RELEASE branch) has 3 differences compared with the Scheme version of Meep :
2422
1254
	* in Python-meep, the center of the coordinate system is in the upper left corner (in Scheme-Meep v1.1.1, the center of the coordinate system is in the middle of your computational volume).
2423
1255
	* in Python-meep, eps-averaging is disabled by default (see section 3.2 for details on how to enable eps-averaging)
2424
1256
	* in Python-meep, calculation is done with complex fields by default (in Scheme-Meep v1.1.1, real fields are used by default). You can call function use_real_fields() on your fields-object to enable calculation with real fields only.
2425
1257
2426
1258
On starting your script, Python-Meep will warn you about these differences. You can suppress these warning by setting the global variables ``DISABLE_WARNING_COORDINATESYSTEM``, ``DISABLE_WARNING_EPS_AVERAGING`` and ``DISABLE_WARNING_REAL_FIELDS`` to  ``True``. You can add site-specific customisations to the file ``meep-site-init.py`` : in this script, you can for example suppress the warning, or enable EPS-averaging by default.
2427
1259
2428
1260
Add the following code to ``meep-site-init.py`` if you want *to calculate with real fields only by default* :
2429
1261
2430
1262
::
2431
1263
2432
1264
	#by default enable calculation with real fields only (consistent with Scheme-meep)
2433
1265
	def fields(structure, m = 0, store_pol_energy=0):
2434
1266
	   f = fields_orig(structure, m, store_pol_energy)
2435
1267
	   f.use_real_fields()
2436
1268
	   return f
2437
1269
2438
1270
	#add a new construct 'fields_complex' that you can use to force calculation with complex fields
2439
1271
	def fields_complex(structure, m = 0, store_pol_energy=0):
2440
1272
	   master_printf("Calculation with complex fields enalbed.\n")
2441
1273
	   return fields_orig(structure, m, store_pol_energy)
2442
1274
2443
1275
	global DISABLE_WARNING_REAL_FIELDS
2444
1276
	DISABLE_WARNING_REAL_FIELDS = True
2445
1277
2446
1278
2447
1279
Add the following code to ``meep-site-init.py`` if you want *to enable EPS-averaging by default* :
2448
1280
2449
1281
::
2450
1282
2451
1283
	use_averaging(True)
2452
1284
2453
1285
	global DISABLE_WARNING_EPS_AVERAGING
2454
1286
	DISABLE_WARNING_EPS_AVERAGING = True
2455
1287
2456
1288
2457
1289
2458
1290
2459
1291
2460
1292
2461
1293
2462
1294
2463
1295
2464
1296
2465
1297
2466
1298
2467
1299
2468
0
1300
2469
=== added file 'doc/html/.buildinfo'
2470
--- doc/html/.buildinfo	1970-01-01 00:00:00 +0000
2471
+++ doc/html/.buildinfo	2009-12-01 14:23:09 +0000
2472
@@ -0,0 +1,4 @@
2473
1
# Sphinx build info version 1
2474
2
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
2475
3
config: be01740a8ae8c3dde5cf900e473f3548
2476
4
tags: fbb0d17656682115ca4d033fb2f83ba1
2477
0
5
2478
=== added directory 'doc/html/_images'
2479
=== added file 'doc/html/_images/bentwgB.gif'
2480
1
Binary files doc/html/_images/bentwgB.gif	1970-01-01 00:00:00 +0000 and doc/html/_images/bentwgB.gif	2009-12-01 14:23:10 +0000 differ
6
Binary files doc/html/_images/bentwgB.gif	1970-01-01 00:00:00 +0000 and doc/html/_images/bentwgB.gif	2009-12-01 14:23:10 +0000 differ
2481
=== added file 'doc/html/_images/bentwgNB.gif'
2482
2
Binary files doc/html/_images/bentwgNB.gif	1970-01-01 00:00:00 +0000 and doc/html/_images/bentwgNB.gif	2009-12-01 14:23:10 +0000 differ
7
Binary files doc/html/_images/bentwgNB.gif	1970-01-01 00:00:00 +0000 and doc/html/_images/bentwgNB.gif	2009-12-01 14:23:10 +0000 differ
2483
=== added file 'doc/html/_images/fluxes.png'
2484
3
Binary files doc/html/_images/fluxes.png	1970-01-01 00:00:00 +0000 and doc/html/_images/fluxes.png	2009-12-01 14:23:10 +0000 differ
8
Binary files doc/html/_images/fluxes.png	1970-01-01 00:00:00 +0000 and doc/html/_images/fluxes.png	2009-12-01 14:23:10 +0000 differ
2485
=== added directory 'doc/html/_sources'
2486
=== added directory 'doc/html/_sources/.svn'
2487
=== added file 'doc/html/_sources/.svn/all-wcprops'
2488
--- doc/html/_sources/.svn/all-wcprops	1970-01-01 00:00:00 +0000
2489
+++ doc/html/_sources/.svn/all-wcprops	2009-12-01 14:23:10 +0000
2490
@@ -0,0 +1,11 @@
2491
1
K 25
2492
2
svn:wc:ra_dav:version-url
2493
3
V 59
2494
4
/svn/python-meep/!svn/ver/44/trunk/doc/_build/html/_sources
2495
5
END
2496
6
python_meep_documentation.txt
2497
7
K 25
2498
8
svn:wc:ra_dav:version-url
2499
9
V 89
2500
10
/svn/python-meep/!svn/ver/70/trunk/doc/_build/html/_sources/python_meep_documentation.txt
2501
11
END
2502
0
12
2503
=== added file 'doc/html/_sources/.svn/entries'
2504
--- doc/html/_sources/.svn/entries	1970-01-01 00:00:00 +0000
2505
+++ doc/html/_sources/.svn/entries	2009-12-01 14:23:10 +0000
2506
@@ -0,0 +1,62 @@
2507
1
9
2508
2
2509
3
dir
2510
4
44
2511
5
http://zita.intec.ugent.be/svn/python-meep/trunk/doc/_build/html/_sources
2512
6
http://zita.intec.ugent.be/svn/python-meep
2513
7
2514
8
2515
9
2516
10
2009-10-19T12:09:25.091644Z
2517
11
44
2518
12
elambert
2519
13
2520
14
2521
15
svn:special svn:externals svn:needs-lock
2522
16
2523
17
2524
18
2525
19
2526
20
2527
21
2528
22
2529
23
2530
24
2531
25
2532
26
2533
27
556452f8-113f-44a6-93c3-95c7ba0365e5
2534
28
2535
0
29
2536
30
python_meep_documentation.txt
2537
31
file
2538
32
70
2539
33
2540
34
2541
35
2542
36
2009-12-01T08:16:44.000000Z
2543
37
6d3942cec680bab591d7645d0696ea53
2544
38
2009-12-01T08:17:30.437718Z
2545
39
70
2546
40
elambert
2547
41
2548
42
2549
43
2550
44
2551
45
2552
46
2553
47
2554
48
2555
49
2556
50
2557
51
2558
52
2559
53
2560
54
2561
55
2562
56
2563
57
2564
58
2565
59
2566
60
2567
61
2568
62
55924
2569
63
2570
1
64
2571
2
65
2572
=== added file 'doc/html/_sources/.svn/format'
2573
--- doc/html/_sources/.svn/format	1970-01-01 00:00:00 +0000
2574
+++ doc/html/_sources/.svn/format	2009-12-01 14:23:10 +0000
2575
@@ -0,0 +1,1 @@
2576
1
9
2577
0
2
2578
=== added directory 'doc/html/_sources/.svn/prop-base'
2579
=== added directory 'doc/html/_sources/.svn/props'
2580
=== added directory 'doc/html/_sources/.svn/text-base'
2581
=== added file 'doc/html/_sources/.svn/text-base/python_meep_documentation.txt.svn-base'
2582
--- doc/html/_sources/.svn/text-base/python_meep_documentation.txt.svn-base	1970-01-01 00:00:00 +0000
2583
+++ doc/html/_sources/.svn/text-base/python_meep_documentation.txt.svn-base	2009-12-01 14:23:10 +0000
2584
@@ -0,0 +1,1293 @@
2585
1
PYTHON-MEEP BINDING DOCUMENTATION
2586
2
====================================
2587
3
2588
4
Primary author of this documentation : EL : Emmanuel.Lambert@intec.ugent.be 
2589
5
2590
6
Document history : 
2591
7
2592
8
::
2593
9
2594
10
	* EL-19/20/21-08-2009 : document creation 
2595
11
	* EL-24-08-2009 : small improvements & clarifications.
2596
12
	* EL-25/26-08-2009 : sections 7 & 8 were added.
2597
13
	* EL-03-04/09/2009 : 
2598
14
               -class "structure_eps_pml" (removed again in v0.8).
2599
15
	       -port to Meep 1.1.1 (class 'volume' was renamed to 'grid_volume' and class 'geometric_volume' to 'volume'
2600
16
	       -minor changes in the bent waveguide sample, to make it more consistent with the Scheme version
2601
17
	* EL-07-08/09/2009 : sections 3.2, 8.2, 8.3 : defining a material function with inline C/C++
2602
18
	* EL-10/09/2009 : additions for MPI mode (multiprocessor)
2603
19
	* EL-21/10/2009 : amplitude factor callback function
2604
20
	* EL-22/10/2009 : keyword arguments for runUntilFieldsDecayed
2605
21
	* EL-01/12/2009 : alignment with version 0.8 - III
2606
22
2607
23
2608
24
2609
25
**1. The general structure of a python-meep program**
2610
26
-----------------------------------------------------
2611
27
2612
28
In general terms, a python-meep program can be structured as follows :
2613
29
2614
30
    * import the python-meep binding :
2615
31
		``from meep import *``
2616
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/``
2617
33
2618
34
      If you are running in MPI mode (multiprocessor, see section 9), then you have to import module ``meep_mpi`` instead :
2619
35
		``from meep_mpi import *``
2620
36
2621
37
    * define a computational grid volume
2622
38
        See section 2 below which explains usage of the ``grid_volume`` class. 
2623
39
2624
40
    * define the waveguide structure (describing the geometry, PML and materials)
2625
41
        See section 3 below which explains usage of the ``structure`` class.
2626
42
		
2627
43
    * create an object which will hold the calculated fields 
2628
44
        See section 4 below which explains usage of the ``field`` class.
2629
45
2630
46
    * define the sources
2631
47
        See section 5 below which explains usage of the ``add_point_source`` and ``add_volume_source`` functions.
2632
48
2633
49
    * run the simulation (iterate over the time-steps)
2634
50
        See section 6 below.
2635
51
2636
52
Section 7 gives details about defining and retrieving fluxes.
2637
53
2638
54
Section 9 gives some complete examples.
2639
55
2640
56
Section 10 outlines some differences between Scheme-Meep and Python-Meep.
2641
57
2642
58
2643
59
**2. Defining the computational grid volume**
2644
60
---------------------------------------------
2645
61
2646
62
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).
2647
63
		* ``volcyl(rsize, zsize, resolution)``
2648
64
			Defines a cyclical computational grid volume.
2649
65
		* ``volone(zsize, resolution)`` *alternatively called* ``vol1d(zsize, resolution)``
2650
66
			Defines a 1-dimensional computational grid volume along the Z-axis.
2651
67
		* ``voltwo(xsize, ysize, resolution)`` *alternatively called* ``vol2d(xsize, ysize, a)``
2652
68
			Defines a 2-dimensional computational grid volumes along the X- and Y-axes
2653
69
		* ``vol3d(xsize, ysize, zsize, resolution)``
2654
70
			Defines a 3-dimensional computational grid volume.
2655
71
2656
72
    e.g.: ``v = volone(6, 10)`` defines a 1-dimensional computational volume of lenght 6, with 10 pixels per distance unit.
2657
73
2658
74
2659
75
**3. Defining the waveguide structure**
2660
76
---------------------------------------
2661
77
	
2662
78
The waveguide structure is defined using the class ``structure``, of which the constructor has the following arguments :
2663
79
2664
80
		* *required* : the computational grid volume (a reference to an object of type ``grid_volume``, see section 2 above)
2665
81
2666
82
		* *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.
2667
83
2668
84
		* *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 :
2669
85
			- ``no_pml()`` describing a conditionless boundary region (no PML)
2670
86
			- ``pml(thickness)`` : decribing a perfectly matching layer (PML) of a certain thickness (double value) on the boundaries in all directions. 
2671
87
			- ``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``). 
2672
88
			- ``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.
2673
89
2674
90
		* *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:
2675
91
			- ``identity`` : no symmetry
2676
92
			- ``rotate4(direction, grid_volume)`` : defines a 90° rotational symmetry with 'direction' the axis of rotation.
2677
93
			- ``rotate2(direction, grid_volume)`` : defines a 180° rotational symmetry with 'direction' the axis of rotation. 
2678
94
			- ``mirror(direction, grid_volume)`` : defines a mirror symmetry plane with 'direction' normal to the mirror plane.
2679
95
			- ``r_to_minus_r_symmetry`` : defines a mirror symmetry in polar coordinates
2680
96
2681
97
		* optional: the number of chunks in which to split up the calculated geometry. If you leave this empty, it is auto-configured. Otherwise, you would set this to a factor which is a multiple of the number of processors in your MPI run (for multiprocessor configuration).
2682
98
2683
99
    e.g. : if ``v`` is a variable pointing to the computational grid volume, then :
2684
100
		 ``s = structure(v, one)`` defines a structure with all air (eps=1), 
2685
101
    which is equivalent to:
2686
102
		 ``s = structure(v, one, no_pml(), identity(), 1)``
2687
103
2688
104
    Another example : ``s = structure(v, EPS, pml(0.1,Y) )`` with EPS a custom material function, which is explained in the note below.
2689
105
2690
106
2691
107
3.1. Defining a material function
2692
108
________________________________________
2693
109
2694
110
In order to describe the geometry of the waveguide, we have to provide a 'material function' returning 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.
2695
111
It is also possible (and faster) to write your material function in inline C/C++ (see 3.3)
2696
112
2697
113
E.g. :
2698
114
2699
115
::
2700
116
2701
117
        class epsilon(Callback):  #inherit from Callback for integration with the meep core library
2702
118
                def __init__(self):
2703
119
                        Callback.__init__(self)
2704
120
                def double_vec(self,vec): #override of function in the Callback class to set the eps function
2705
121
                        self.set_double(self.eps(vec))
2706
122
                        return 
2707
123
                def eps(self,vec):   #return the epsilon value for a certain point (indicated by the vector v)
2708
124
                        v = vec 
2709
125
                        r = v.x()*v.x() + v.y()*v.y()
2710
126
                        dr = sqrt(r)
2711
127
                        while dr>1:
2712
128
                            dr-=1
2713
129
                        if dr > 0.7001:
2714
130
                            return 12.0
2715
131
                        return 1.0
2716
132
2717
133
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.
2718
134
So, one should write : ``vec.x()`` and ``vec.y()``.
2719
135
 
2720
136
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 :
2721
137
2722
138
::
2723
139
2724
140
        set_EPS_Callback(epsilon().__disown__())
2725
141
        s = structure(v, EPS, no_pml(), identity())
2726
142
2727
143
2728
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.
2729
145
2730
146
***Important remark*** : at the end of our program, we should call : ``del_EPS_Callback()`` in order to clean up the global variable.
2731
147
2732
148
For better performance, you can define your EPS material function with inline C/C++ : we refer to section 3.3 for details about this.
2733
149
2734
150
3.2 Eps-averaging
2735
151
_________________
2736
152
2737
153
EPS-averaging (anisotrpic averaging) is disabled by default, making this behaviour consistent with the behaviour of the Meep C++ core.
2738
154
2739
155
You can enable EPS-averaging using the function ``use_averaging`` :
2740
156
2741
157
::
2742
158
2743
159
	#enable EPS-averaging
2744
160
	use_averaging(True)
2745
161
	...
2746
162
	#disable EPS-averaging
2747
163
	use_averaging(False)
2748
164
	...
2749
165
2750
166
2751
167
Enabling EPS-averaging results in slower performance, but more accurate results. 
2752
168
2753
169
2754
170
3.3. Defining a material function with inline C/C++
2755
171
_________________________________________________________
2756
172
2757
173
The approach described in 3.1 lets the Meep core library call back to Python for every query of the epsilon-function. This creates a lot of overhead.
2758
174
An approach which has a lot better performance is to define this epsilon-function with an inline C-function or C++ class.
2759
175
2760
176
* If our epsilon-function needs *no other parameters than the position vector (X, Y, Z)*, then we can suffice with an inline C-function (the geometry dependencies are then typically hardcoded).
2761
177
2762
178
* If our epsilon-function needs to base it's calculation on *a more complex set of parameters (e.g. parameters depending on the geometry)*, then we have to write a C++ class. 
2763
179
2764
180
For example, in the bent-waveguide example (section 8.3), we can define a generic C++ class which can return the epsilon-value for both the "bend" and "no bend" case, with variable size parameters.
2765
181
We can also take a simpler approach (section 8.2) and write a function in which the geometry size parameters are hardcoded : we then need 2 inline C-functions : one for the "bend" case and one for the "no bend" case. 
2766
182
	
2767
183
2768
184
3.3.1 Inline C-function
2769
185
.......................
2770
186
2771
187
First we create a header file, e.g. "eps_function.hpp" which contains our EPS-function.
2772
188
Not that the geometry dependencies are hardcoded (``upperLimitHorizontalWg = 4`` and ``lowerLimitHorizontalWg = 5``).
2773
189
2774
190
::
2775
191
2776
192
2777
193
	namespace meep
2778
194
	{
2779
195
		static double myEps(const vec &v) {
2780
196
			double xCo = v.x();
2781
197
			double yCo = v.y();
2782
198
			double upperLimitHorizontalWg = 4;
2783
199
			double lowerLimitHorizontalWg = 5;
2784
200
			if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
2785
201
				return 1.0;
2786
202
			}
2787
203
			else return 12.0;
2788
204
		}
2789
205
	}
2790
206
2791
207
2792
208
Then, in the Python program, we prepare and set the callback function as shown below.
2793
209
``prepareCallbackCfunction`` returns a pointer to the C-function, which we deliver to the Meep core using ``set_EPS_CallbackInlineFunction``.
2794
210
2795
211
::
2796
212
2797
213
	def initEPS(isWgBent):
2798
214
		funPtr = prepareCallbackCfunction("myEps","eps_function.hpp")  #name of your function / name of header file
2799
215
		set_EPS_CallbackInlineFunction(funPtr)
2800
216
		print "EPS function successfully set."
2801
217
		return funPtr
2802
218
2803
219
We refer to section 8.2 below for a full example.
2804
220
2805
221
2806
222
3.3.2 Inline C++-class
2807
223
......................
2808
224
2809
225
A more complex approach is to have a C++ object that can accept more parameters when it is constructed. 
2810
226
For example this is the case if want to change the parameters of the geometry from inside Python without touching the C++ code.
2811
227
2812
228
We create a header file "eps_class.hpp" which contains the definition of the class (the class must inherit from ``Callback``).
2813
229
In the example below, the parameters ``upperLimitHorizontalWg`` and ``widthWg`` will be communicated from Python upon construction of the object.
2814
230
If these parameters then change (depending on the geometry), the C++ object will follow automatically.
2815
231
2816
232
2817
233
::
2818
234
2819
235
	using namespace meep;
2820
236
2821
237
	namespace meep
2822
238
	{
2823
239
2824
240
	class myEpsCallBack : virtual public Callback {
2825
241
2826
242
	    public:	
2827
243
		   myEpsCallBack() : Callback() { };
2828
244
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
2829
245
2830
246
		   myEpsCallBack(double upperLimitHorizontalWg, double widthWg)  : Callback() {
2831
247
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
2832
248
			_widthWg = widthWg;
2833
249
		   };	
2834
250
2835
251
		   double double_vec(const vec &v) { //return the EPS-value, depending on the position vector
2836
252
			double eps = myEps(v, _upperLimitHorizontalWg, _widthWg);
2837
253
			return eps;
2838
254
		   };
2839
255
		   
2840
256
		   complex<double> complex_vec(const vec &x) { return 0; }; //no need to implement
2841
257
		   complex<double> complex_time(double &t) { return 0; }; //no need to implement //-->> SUBJECT TO CHANGE - in Intec branch of v0.8, this is complex_time(double t)
2842
258
2843
259
2844
260
	     private:
2845
261
		   double _upperLimitHorizontalWg;
2846
262
		   double _widthWg;
2847
263
2848
264
		   double myEps(const vec &v, double upperLimitHorizontalWg, double widthWg) {
2849
265
			double xCo = v.x();
2850
266
			double yCo = v.y();
2851
267
	    	        if ((yCo < upperLimitHorizontalWg) || (yCo > upperLimitHorizontalWg+widthWg)){
2852
268
				return 1.0;
2853
269
			    }
2854
270
			}
2855
271
			return 12.0;
2856
272
		   }
2857
273
2858
274
	};
2859
275
2860
276
	}
2861
277
2862
278
2863
279
The syntax in Python is a little bit more complex in this case.
2864
280
2865
281
We will need to import the module ``scipy.weave`` : 
2866
282
2867
283
``from scipy.weave import *``
2868
284
2869
285
(this is not required for the previous case of a simple inline function)
2870
286
2871
287
First we create a list with the names of the parameters that we want to pass to the C++ class. These variables must be declared in the scope where the ``inline`` function call happens (see below).
2872
288
2873
289
``c_params = ['upperLimitHorizontalWg','widthWg']``
2874
290
2875
291
Then, we prepare the code snippet, using the function ``prepareCallbackCObjectCode`` and passing the class name and parameter names list.
2876
292
2877
293
``c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)``
2878
294
2879
295
Finally, we call the ``inline`` function, passing :
2880
296
	* the code snippet
2881
297
	* the list of parameter names
2882
298
	* the inline libraries, include directories and headers (helper functions are provided for this, see below). The call to ``getInlineHeaders`` should receive the name of the header file (with the definition of the C++ class) as an argument.
2883
299
2884
300
``funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )``
2885
301
2886
302
::
2887
303
2888
304
2889
305
	def initEPS():
2890
306
		#the set of parameters that we want to pass to the Callback object upon construction
2891
307
		#all these variables must be declared in the scope where the "inline" function call happens
2892
308
		c_params = ['upperLimitHorizontalWg','widthWg'] 
2893
309
		#the C-code snippet for constructing the Callback object
2894
310
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
2895
311
		#do the actual inline C-call and fetch the pointer to the Callback object
2896
312
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
2897
313
		#set the pointer to the callback object in the Python-meep core 
2898
314
		set_EPS_CallbackInlineObject(funPtr)
2899
315
		print "EPS function successfully set."
2900
316
		return
2901
317
 
2902
318
2903
319
We refer to section 8.3 below for a full example.
2904
320
2905
321
2906
322
2907
323
**4. Defining the initial field**
2908
324
---------------------------------
2909
325
2910
326
This is optional.
2911
327
2912
328
We create an object of type ``fields``, which will contain the calculated field.
2913
329
2914
330
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``.
2915
331
2916
332
::
2917
333
2918
334
            class fi(Callback):  #inherit from Callback for integration with the meep core library
2919
335
                def __init__(self):
2920
336
                    Callback.__init__(self)		
2921
337
                def complex_vec(self,v): #override of function in the Callback class to set the field initialization function
2922
338
		     #return the field value for a certain point (indicated by the vector v)
2923
339
                    return complex(1.0,0)
2924
340
2925
341
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 : 
2926
342
2927
343
	``set_INIF_Callback(fi().__disown__())  #link the INIF variable to the fi class``
2928
344
2929
345
We refer to section 3-note1 for more information about the function ``__disown__()``.
2930
346
2931
347
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 :
2932
348
2933
349
::
2934
350
2935
351
    f = fields(s)
2936
352
    comp = Hy
2937
353
    f.initialize_field(comp, INIF)
2938
354
2939
355
2940
356
***Important remark*** : at the end of our program, we should then call : ``del_INIF_Callback()`` in order to clean up the global variable.
2941
357
2942
358
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.
2943
359
2944
360
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. :
2945
361
2946
362
	``f.use_bloch(vec(0.0))``
2947
363
2948
364
*to be further elaborated - what is the exact meaning of the vector argument? (not well understood at this time)*
2949
365
2950
366
2951
367
2952
368
**5. Defining the sources**
2953
369
---------------------------
2954
370
2955
371
The definition of the current sources can be done through 2 functions of the ``fields`` class :
2956
372
        * ``add_point_source(component, src_time, vec, complex)``
2957
373
	* ``add_volume_source(component, src_time, volume, complex)``
2958
374
2959
375
2960
376
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).
2961
377
2962
378
For a point source, we must specify the center point of the current source using a vector (object of type ``vec``).
2963
379
2964
380
For a volume source, we must specify an object of type ``volume`` (*to be elablorated*).
2965
381
2966
382
The last argument is an overall complex amplitude number, multiplying the current source (default 1.0).
2967
383
2968
384
The following variants are available :
2969
385
        * ``add_point_source(component, double, double, double, double, vec centerpoint, complex amplitude, int is_continuous)``
2970
386
                * This is a shortcut function so that no ``src_time`` object must be created. *This function is preferably used for point sources.*
2971
387
                * The four real arguments define the central frequency, spectral width, peaktime and cutoff. 
2972
388
2973
389
        * ``add_volume_source(component, src_time, volume)``
2974
390
2975
391
        * ``add_volume_source(component, src_time, volume, AMPL)``
2976
392
                * AMPL is a built-in reference to a callback function. Such a callback function returns a factor to multiply the source amplitude with (complex value). It receives 1 parameter, i.e. a vector indicating a position RELATIVE to the CENTER of the source. See the example below.
2977
393
2978
394
2979
395
Three classes, inheriting from ``src_time``, are predefined and can be used off the shelf :
2980
396
	* ``gaussian_src_time`` for a Gaussian-pulse source. The constructor demands 2 arguments of type double :
2981
397
                    * the center frequency ω, in units of 2πc 
2982
398
		    * the frequency width w used in the Gaussian
2983
399
	* ``continuous_src_time`` for a continuous-wave source proportional to exp(−iωt). The constructor demands 4 arguments : 
2984
400
                    * the frequency ω, in units 2πc/distance (complex number)
2985
401
		    * the temporal width of smoothing (default 0)
2986
402
                    * the start time (default 0)
2987
403
                    * the end time (default infinity = never turn off)
2988
404
	* ``custom_src_time`` for a user-specified source function f(t) with start/end times. The constructor demands 4 arguments :
2989
405
		    * The function f(t) specifying the time-dependence of the source
2990
406
                    * *...(2nd argument unclear, to be further elaborated)...*
2991
407
                    * the start time of the source (default -infinity)
2992
408
                    * the end time of the source (default +infinity)    
2993
409
2994
410
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 :
2995
411
           
2996
412
::
2997
413
2998
414
	#define a continuous source 
2999
415
	srcFreq = 0.125
3000
416
	srcWidth = 20 
3001
417
	src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
3002
418
	srcComp = Ez
3003
419
	#make it a line source of size 1 starting on position(6,3)
3004
420
	srcGeo = volume(vec(6,3),vec(6,4))
3005
421
	f.add_volume_source(srcComp, src, srcGeo)
3006
422
3007
423
3008
424
Here is an example of the implementation of a callback function for the amplitude factor :
3009
425
3010
426
::
3011
427
3012
428
	class amplitudeFactor(Callback):
3013
429
	    def __init__(self):
3014
430
		Callback.__init__(self)
3015
431
		master_printf("Callback function for amplitude factor activated.\n")
3016
432
3017
433
	    def complex_vec(self,vec):
3018
434
		#BEWARE, these are coordinates RELATIVE to the source center !!!!
3019
435
		x = vec.x()
3020
436
		y = vec.y()
3021
437
		master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
3022
438
		result = complex(1.0,0)
3023
439
		return result
3024
440
3025
441
	...
3026
442
		#define a continuous source 
3027
443
		srcFreq = 0.125
3028
444
		srcWidth = 20 
3029
445
		src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
3030
446
		srcComp = Ez
3031
447
		#make it a line source of size 1 starting on position(6,3)
3032
448
		srcGeo = volume(vec(6,3),vec(6,4))
3033
449
		#create callback object for amplitude factor
3034
450
		af = amplitudeFactor()
3035
451
		set_AMPL_Callback(af.__disown__())
3036
452
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)
3037
453
3038
454
3039
455
**6. Running the simulation, retrieving field values and exporting HDF5 files**
3040
456
-------------------------------------------------------------------------------
3041
457
3042
458
We can now time-step and retrieve various field values along the way.
3043
459
The actual time step value can be retrieved or set through the variable ``f.dt``.
3044
460
3045
461
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.
3046
462
3047
463
To trigger a step in time, you call the function ``f.step()``.
3048
464
3049
465
To step until the source has fully decayed :
3050
466
3051
467
::
3052
468
3053
469
        while (f.time() < f.last_source_time()):
3054
470
             f.step()
3055
471
3056
472
The function ``runUntilFieldsDecayed`` mimicks the behaviour of 'stop-when-fields-decayed' in Meep-Scheme.
3057
473
This will run time steps until the source has decayed to 0.001 of the peak amplitude. After that, by default an additional 50 time steps will be run.
3058
474
The function has 7 arguments, of which 4 are mandatory and 3 are optional keywords arguments :
3059
475
* 4 regular arguments : reference to the field, reference to the computational grid volume, the source component, the monitor point.
3060
476
* keyword argument ``pHDF5OutputFile`` : reference to a HDF5 file (constructed with the function ``prepareHDF5File``); default : None (no ouput to files)
3061
477
* keyword argument ``pH5OutputIntervalSteps`` : step interval for output to HDF5 (default : 50)
3062
478
* keyword argument ``pDecayedStopAfterSteps`` : the number of steps to continue after the source has decayed to 0.001 of the peak amplitude at the probing point (default: 50)
3063
479
 
3064
480
We further refer to section 8 below where this function is applied in an example.
3065
481
3066
482
A rich functionality is available for retrieving field information. Some examples :
3067
483
3068
484
	* ``f.energy_in_box(v.surroundings())``
3069
485
	* ``f.electric_energy_in_box(v.surroundings())``
3070
486
	* ``f.magnetic_energy_in_box(v.surroundings())``
3071
487
	* ``f.thermo_energy_in_box(v.surroundings())``
3072
488
	* ``f.total_energy()``
3073
489
	* ``f.field_energy_in_box(v.surroundings())``
3074
490
	* ``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``)
3075
491
	* ``f.flux_in_box(X, v.surroundings())`` where the first argument is the direction (``X, Y, Z, R`` or ``P``)
3076
492
3077
493
We can probe the field at certain points by defining a *monitor point* as follows : 
3078
494
3079
495
::
3080
496
3081
497
        m = monitor_point()
3082
498
        p = vec(2.10) #vector identifying the point that we want to probe
3083
499
        f.get_point(m, p)
3084
500
        m.get_component(Hx)
3085
501
3086
502
We can export the dielectric function and e.g. the Ex component of the field to HDF5 files as follows :
3087
503
3088
504
::
3089
505
3090
506
        #make sure you start your Python session with 'sudo' or write rights to the current path
3091
507
        feps = prepareHDF5File("eps.h5")
3092
508
        f.output_hdf5(Dielectric, v.surroundings(), feps)   #export the Dielectric structure so that we can visually verify it 
3093
509
	fex = prepareHDF5File("ex.h5")
3094
510
	while (f.time() < f.last_source_time()):
3095
511
             f.step()
3096
512
             f.output_hdf5(Ex, v.surroundings(), fex, 1) #export the Ex component, appending to the file "ex.h5"
3097
513
3098
514
3099
515
3100
516
**7. Defining and retrieving fluxes**
3101
517
--------------------------------------
3102
518
3103
519
First we define a flux plane. 
3104
520
This is done through the creation of an object of type ``volume`` (specifying 2 vectors as arguments).
3105
521
3106
522
Then we apply this flux plane to the field, specifying 4 parameters :
3107
523
	* the reference to the ``volume`` object
3108
524
	* the minimum frequency (in the example below, this is ``srcFreqCenter-(srcPulseWidth/2.0)``)
3109
525
	* the maximum frequency (in the example below this is ``srcFreqCenter+(srcPulseWidth/2.0)`` )
3110
526
	* the number of discrete frequencies that we want to monitor in the flux (in the example below, this is 100).
3111
527
3112
528
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.
3113
529
3114
530
E.g., if ``f`` is the field, then we proceed as follows : 
3115
531
3116
532
::
3117
533
3118
534
	#define the flux plane and flux parameters
3119
535
	fluxplane = volume(vec(1,2),vec(1,3))
3120
536
	flux = f.add_dft_flux_plane(fluxplane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3121
537
3122
538
	#run the calculation 
3123
539
	while (f.time() < f.last_source_time()):
3124
540
	    f.step()
3125
541
3126
542
	#retrieve the flux data
3127
543
	fluxdata = getFluxData(flux)
3128
544
	frequencies = fluxdata[0]
3129
545
	fluxvalues = fluxdata[1]
3130
546
3131
547
3132
548
3133
549
**8. The "90 degree bent waveguide example in Python**
3134
550
------------------------------------------------------
3135
551
3136
552
We have ported the "90 degree bent waveguide" example from the Meep-Scheme tutorial to Python.
3137
553
3138
554
The original example can be found on the following URL : http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial
3139
555
(section 'A 90° bend').
3140
556
3141
557
You can find the source code below and also in the ``/samples/bent_waveguide`` directory of the Python-Meep distribution. 
3142
558
3143
559
Sections 8.2 and 8.3 contain the same example but with the EPS-callback function as inline C-function and inline C++-class.
3144
560
3145
561
The following animated gifs can be produced from the HDF5-files (see the script included in directory 'samples') :
3146
562
3147
563
.. image:: images/bentwgNB.gif
3148
564
3149
565
.. image:: images/bentwgB.gif
3150
566
3151
567
3152
568
And here is the graph of the transmission, reflection and loss fluxes, showing the same results as the example in Scheme:
3153
569
3154
570
.. image:: images/fluxes.png
3155
571
   :height: 315
3156
572
   :width: 443
3157
573
3158
574
3159
575
8.1 With a Python class as EPS-function
3160
576
________________________________________
3161
577
3162
578
3163
579
A bottleneck in this version is the epsilon-function, which is written in Python. 
3164
580
This means that the Meep core library must do a callback to the Python function, which creates a lot of overhead.
3165
581
An approach which has a much better performance is to write this EPS-function in C : the Meep core library can then directly call back to a C-function.
3166
582
These approaches are described in 8.2 and 8.3.
3167
583
3168
584
::
3169
585
3170
586
	from meep import *
3171
587
	from math import *
3172
588
	from python_meep import *
3173
589
	import numpy 
3174
590
	import matplotlib.pyplot
3175
591
	import sys
3176
592
3177
593
	res = 10.0
3178
594
	gridSizeX = 16.0
3179
595
	gridSizeY = 32.0
3180
596
	padSize = 4.0
3181
597
	wgLengthX = gridSizeX - padSize
3182
598
	wgLengthY = gridSizeY - padSize
3183
599
	wgWidth = 1.0 #width of the waveguide
3184
600
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
3185
601
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
3186
602
	srcFreqCenter = 0.15 #gaussian source center frequency
3187
603
	srcPulseWidth = 0.1 #gaussian source pulse width
3188
604
	srcComp = Ez #gaussian source component
3189
605
3190
606
	#this function plots the waveguide material as a function of a vector(X,Y) 
3191
607
	class epsilon(Callback):
3192
608
	    def __init__(self, pIsWgBent):
3193
609
		Callback.__init__(self)
3194
610
		self.isWgBent = pIsWgBent
3195
611
	    def double_vec(self,vec):
3196
612
		if (self.isWgBent): #there is a bend
3197
613
		    if ((vec.x()<wgLengthX) and (vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
3198
614
		        return 12.0
3199
615
		    elif ((vec.x()>=wgLengthX-wgWidth) and (vec.x()<=wgLengthX) and vec.y()>= padSize ):
3200
616
		        return 12.0
3201
617
		    else:
3202
618
		        return 1.0
3203
619
		else: #there is no bend 
3204
620
		    if ((vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
3205
621
		        return 12.0
3206
622
		    else:
3207
623
		        return 1.0
3208
624
3209
625
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
3210
626
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
3211
627
		#in all directions,
3212
628
		#using the material function EPS
3213
629
		material = epsilon(pIsWgBent)
3214
630
		set_EPS_Callback(material.__disown__())
3215
631
		s = structure(pCompVol, EPS, pml(1.0) )     
3216
632
		f = fields(s)      
3217
633
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
3218
634
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
3219
635
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
3220
636
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
3221
637
		print "Field created..."
3222
638
		return f
3223
639
	       
3224
640
3225
641
	#FIRST WE WORK OUT THE CASE WITH NO BEND
3226
642
	#----------------------------------------------------------------
3227
643
	print "*1* Starting the case with no bend..."
3228
644
	#create the computational grid 
3229
645
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
3230
646
3231
647
	#create the field
3232
648
	wgBent = 0 #no bend
3233
649
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3234
650
3235
651
	bendFnEps = "./bentwgNB_Eps.h5"
3236
652
	bendFnEz = "./bentwgNB_Ez.h5"
3237
653
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3238
654
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
3239
655
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
3240
656
3241
657
	#create the file for the field components 
3242
658
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
3243
659
3244
660
	#define the flux plane for the reflected flux
3245
661
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3246
662
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3247
663
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3248
664
3249
665
	#define the flux plane for the transmitted flux
3250
666
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
3251
667
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
3252
668
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3253
669
3254
670
	print "Calculating..."
3255
671
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
3256
672
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
3257
673
	print "Done..!"
3258
674
3259
675
	#construct 2-dimensional array with the flux plane data
3260
676
	#see python_meep.py 
3261
677
	noBendReflFlux = getFluxData(noBendReflectedFlux)
3262
678
	noBendTransmFlux = getFluxData(noBendTransmFlux)
3263
679
3264
680
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
3265
681
	noBendReflectedFlux.scale_dfts(-1);
3266
682
	f = open("minusflux.h5", 'w') #truncate file if already exists
3267
683
	f.close()   
3268
684
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
3269
685
3270
686
	del_EPS_Callback()
3271
687
3272
688
3273
689
	#AND SECONDLY FOR THE CASE WITH BEND
3274
690
	#----------------------------------------------------------------
3275
691
	print "*2* Starting the case with bend..."
3276
692
	#create the computational grid 
3277
693
	bendVol = voltwo(gridSizeX,gridSizeY,res)
3278
694
3279
695
	#create the field
3280
696
	wgBent = 1 #there is a bend
3281
697
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3282
698
3283
699
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3284
700
	bendFnEps = "./bentwgB_Eps.h5"
3285
701
	bendFnEz = "./bentwgB_Ez.h5"
3286
702
	bendDielectricFile = prepareHDF5File(bendFnEps)
3287
703
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
3288
704
3289
705
	#create the file for the field components 
3290
706
	bendFileOutputEz = prepareHDF5File(bendFnEz)
3291
707
3292
708
	#define the flux plane for the reflected flux
3293
709
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3294
710
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3295
711
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3296
712
3297
713
	#load the minused reflection flux from the "no bend" case as initalization
3298
714
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
3299
715
3300
716
3301
717
	#define the flux plane for the transmitted flux
3302
718
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
3303
719
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
3304
720
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3305
721
3306
722
	print "Calculating..."
3307
723
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
3308
724
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
3309
725
	print "Done..!"
3310
726
3311
727
	#construct 2-dimensional array with the flux plane data
3312
728
	#see python_meep.py 
3313
729
	bendReflFlux = getFluxData(bendReflectedFlux)
3314
730
	bendTransmFlux = getFluxData(bendTransmFlux)
3315
731
3316
732
	del_EPS_Callback()
3317
733
3318
734
	#SHOW THE RESULTS IN A PLOT
3319
735
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
3320
736
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
3321
737
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
3322
738
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
3323
739
		                  
3324
740
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
3325
741
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
3326
742
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
3327
743
3328
744
	matplotlib.pyplot.show()
3329
745
3330
746
3331
747
8.2 With an inline C-function as EPS-function
3332
748
______________________________________________
3333
749
3334
750
The header file "eps_function.hpp" :
3335
751
3336
752
::
3337
753
3338
754
	using namespace meep;
3339
755
3340
756
	namespace meep
3341
757
	{
3342
758
	static double myEps(const vec &v, bool isWgBent) {
3343
759
		double xCo = v.x();
3344
760
		double yCo = v.y();
3345
761
		double upperLimitHorizontalWg = 4;
3346
762
		double wgLengthX = 12;
3347
763
		double leftLimitVerticalWg = 11;
3348
764
		double lowerLimitHorizontalWg = 5;
3349
765
		if (isWgBent){ //there is a bend
3350
766
		    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
3351
767
			return 1.0;            
3352
768
		    }
3353
769
		    else {
3354
770
			if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
3355
771
			     return 1.0;
3356
772
			}
3357
773
			else {
3358
774
			     return 12.0;
3359
775
			}
3360
776
		    }
3361
777
		}
3362
778
		else { //there is no bend 
3363
779
		    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
3364
780
			return 1.0;
3365
781
		    }
3366
782
		}
3367
783
		return 12.0;
3368
784
	   }
3369
785
3370
786
	static double myEpsBentWg(const vec &v) {
3371
787
		return myEps(v, true);
3372
788
	   }
3373
789
3374
790
	static double myEpsStraightWg(const vec &v) {
3375
791
		return myEps(v, false);
3376
792
	   }
3377
793
	}
3378
794
3379
795
3380
796
3381
797
And the actual Python program :
3382
798
3383
799
3384
800
::
3385
801
3386
802
3387
803
	from meep import *  
3388
804
3389
805
	from math import *
3390
806
	import numpy 
3391
807
	import matplotlib.pyplot
3392
808
	import sys
3393
809
3394
810
	res = 10.0
3395
811
	gridSizeX = 16.0
3396
812
	gridSizeY = 32.0
3397
813
	padSize = 4.0
3398
814
	wgLengthX = gridSizeX - padSize
3399
815
	wgLengthY = gridSizeY - padSize
3400
816
	wgWidth = 1.0 #width of the waveguide
3401
817
	upperLimitHorizontalWg = padSize
3402
818
	lowerLimitHorizontalWg = padSize+wgWidth
3403
819
	leftLimitVerticalWg = wgLengthX-wgWidth
3404
820
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
3405
821
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
3406
822
	srcFreqCenter = 0.15 #gaussian source center frequency
3407
823
	srcPulseWidth = 0.1 #gaussian source pulse width
3408
824
	srcComp = Ez #gaussian source component
3409
825
3410
826
	def initEPS(isWgBent):
3411
827
		if (isWgBent):
3412
828
		    funPtr = prepareCallbackCfunction("myEpsBentWg","eps_function.hpp")
3413
829
		else:
3414
830
		    funPtr = prepareCallbackCfunction("myEpsStraightWg","eps_function.hpp")
3415
831
		set_EPS_CallbackInlineFunction(funPtr)
3416
832
		print "EPS function successfully set."
3417
833
		return funPtr
3418
834
	 
3419
835
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
3420
836
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
3421
837
		#in all directions,
3422
838
		#using the material function EPS
3423
839
		s = structure(pCompVol, EPS, pml(1.0) )     
3424
840
		f = fields(s)      
3425
841
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
3426
842
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
3427
843
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
3428
844
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
3429
845
		print "Field created..."
3430
846
		return f
3431
847
3432
848
3433
849
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C-FUNCTION FOR EPS\n")
3434
850
	       
3435
851
	master_printf("Running on %d processor(s)...\n",count_processors())
3436
852
3437
853
	#FIRST WE WORK OUT THE CASE WITH NO BEND
3438
854
	#----------------------------------------------------------------
3439
855
	master_printf("*1* Starting the case with no bend...")
3440
856
3441
857
	#set EPS material function
3442
858
	initEPS(0)
3443
859
3444
860
	#create the computational grid 
3445
861
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
3446
862
3447
863
	#create the field
3448
864
	wgBent = 0 #no bend
3449
865
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3450
866
3451
867
	bendFnEps = "./bentwgNB_Eps.h5"
3452
868
	bendFnEz = "./bentwgNB_Ez.h5"
3453
869
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3454
870
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
3455
871
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
3456
872
3457
873
	#create the file for the field components 
3458
874
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
3459
875
3460
876
	#define the flux plane for the reflected flux
3461
877
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3462
878
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3463
879
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3464
880
3465
881
	#define the flux plane for the transmitted flux
3466
882
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
3467
883
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
3468
884
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3469
885
3470
886
	master_printf("Calculating...")
3471
887
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
3472
888
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
3473
889
	master_printf("Done..!")
3474
890
3475
891
	#construct 2-dimensional array with the flux plane data
3476
892
	#see python_meep.py 
3477
893
	noBendReflFlux = getFluxData(noBendReflectedFlux)
3478
894
	noBendTransmFlux = getFluxData(noBendTransmFlux)
3479
895
3480
896
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
3481
897
	noBendReflectedFlux.scale_dfts(-1);
3482
898
	f = open("minusflux.h5", 'w') #truncate file if already exists
3483
899
	f.close()   
3484
900
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
3485
901
3486
902
	del_EPS_Callback() #destruct the inline-created object
3487
903
3488
904
3489
905
	#AND SECONDLY FOR THE CASE WITH BEND
3490
906
	#----------------------------------------------------------------
3491
907
	master_printf("*2* Starting the case with bend...")
3492
908
3493
909
	#set EPS material function
3494
910
	initEPS(1)
3495
911
3496
912
	#create the computational grid 
3497
913
	bendVol = voltwo(gridSizeX,gridSizeY,res)
3498
914
3499
915
	#create the field
3500
916
	wgBent = 1 #there is a bend
3501
917
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3502
918
3503
919
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3504
920
	bendFnEps = "./bentwgB_Eps.h5"
3505
921
	bendFnEz = "./bentwgB_Ez.h5"
3506
922
	bendDielectricFile = prepareHDF5File(bendFnEps)
3507
923
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
3508
924
3509
925
	#create the file for the field components 
3510
926
	bendFileOutputEz = prepareHDF5File(bendFnEz)
3511
927
3512
928
	#define the flux plane for the reflected flux
3513
929
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3514
930
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3515
931
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3516
932
3517
933
	#load the minused reflection flux from the "no bend" case as initalization
3518
934
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
3519
935
3520
936
3521
937
	#define the flux plane for the transmitted flux
3522
938
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
3523
939
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
3524
940
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3525
941
3526
942
	master_printf("Calculating...")
3527
943
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
3528
944
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
3529
945
	master_printf("Done..!")
3530
946
3531
947
	#construct 2-dimensional array with the flux plane data
3532
948
	#see python_meep.py 
3533
949
	bendReflFlux = getFluxData(bendReflectedFlux)
3534
950
	bendTransmFlux = getFluxData(bendTransmFlux)
3535
951
3536
952
	del_EPS_Callback()
3537
953
3538
954
	#SHOW THE RESULTS IN A PLOT
3539
955
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
3540
956
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
3541
957
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
3542
958
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
3543
959
		                  
3544
960
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
3545
961
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
3546
962
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
3547
963
3548
964
	matplotlib.pyplot.show()
3549
965
3550
966
3551
967
3552
968
8.3 With an inline C++ class as EPS-function
3553
969
______________________________________________
3554
970
3555
971
3556
972
The header file "eps_class.hpp" :
3557
973
3558
974
3559
975
::
3560
976
3561
977
3562
978
	using namespace meep;
3563
979
3564
980
	namespace meep
3565
981
	{
3566
982
3567
983
	class myEpsCallBack : virtual public Callback {
3568
984
3569
985
	    public:	
3570
986
		   myEpsCallBack() : Callback() { };
3571
987
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
3572
988
3573
989
		   myEpsCallBack(bool isWgBent,double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX)  : Callback() {
3574
990
		    	_IsWgBent = isWgBent;
3575
991
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
3576
992
			_leftLimitVerticalWg = leftLimitVerticalWg;
3577
993
			_lowerLimitHorizontalWg = lowerLimitHorizontalWg;
3578
994
			_wgLengthX = wgLengthX;
3579
995
		   };	
3580
996
3581
997
		   double double_vec(const vec &v) {
3582
998
			double eps = myEps(v, _IsWgBent, _upperLimitHorizontalWg, _leftLimitVerticalWg, _lowerLimitHorizontalWg, _wgLengthX);
3583
999
			//cout << "X="<<v.x()<<"--Y="<<v.y()<<"--eps="<<eps<<"-"<<_upperLimitHorizontalWg<<"--"<<_leftLimitVerticalWg<<"--"<<_lowerLimitHorizontalWg<<"--"<<_wgLengthX;
3584
1000
			return eps;
3585
1001
		   };
3586
1002
		   
3587
1003
		   complex<double> complex_vec(const vec &x) { return 0; };
3588
1004
		   complex<double> complex_time(double &t) { return 0; };   //-->> SUBJECT TO CHANGE - in Intec branch of v0.8, this is complex_time(double t)
3589
1005
3590
1006
3591
1007
	     private:
3592
1008
		   bool _IsWgBent;;
3593
1009
		   double _upperLimitHorizontalWg;
3594
1010
		   double _leftLimitVerticalWg;
3595
1011
		   double _lowerLimitHorizontalWg;
3596
1012
		   double _wgLengthX;
3597
1013
3598
1014
		   double myEps(const vec &v, bool isWgBent, double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX) {
3599
1015
			double xCo = v.x();
3600
1016
			double yCo = v.y();
3601
1017
			if (isWgBent){ //there is a bend
3602
1018
			    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
3603
1019
				return 1.0;            
3604
1020
			    }
3605
1021
			    else {
3606
1022
				if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
3607
1023
				     return 1.0;
3608
1024
				}
3609
1025
				else {
3610
1026
				     return 12.0;
3611
1027
				}
3612
1028
			    }
3613
1029
			}
3614
1030
			else { //there is no bend 
3615
1031
			    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
3616
1032
				return 1.0;
3617
1033
			    }
3618
1034
			}
3619
1035
			return 12.0;
3620
1036
		   }
3621
1037
3622
1038
	};
3623
1039
3624
1040
	}
3625
1041
3626
1042
3627
1043
The Python program :
3628
1044
3629
1045
3630
1046
::
3631
1047
3632
1048
3633
1049
	from meep import *  
3634
1050
3635
1051
	from math import *
3636
1052
	import numpy 
3637
1053
	import matplotlib.pyplot
3638
1054
	import sys
3639
1055
3640
1056
	from scipy.weave import *
3641
1057
3642
1058
	res = 10.0
3643
1059
	gridSizeX = 16.0
3644
1060
	gridSizeY = 32.0
3645
1061
	padSize = 4.0
3646
1062
	wgLengthX = gridSizeX - padSize
3647
1063
	wgLengthY = gridSizeY - padSize
3648
1064
	wgWidth = 1.0 #width of the waveguide
3649
1065
	upperLimitHorizontalWg = padSize
3650
1066
	lowerLimitHorizontalWg = padSize+wgWidth
3651
1067
	leftLimitVerticalWg = wgLengthX-wgWidth
3652
1068
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
3653
1069
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
3654
1070
	srcFreqCenter = 0.15 #gaussian source center frequency
3655
1071
	srcPulseWidth = 0.1 #gaussian source pulse width
3656
1072
	srcComp = Ez #gaussian source component
3657
1073
3658
1074
3659
1075
	def initEPS():
3660
1076
		#the set of parameters that we want to pass to the Callback object upon construction
3661
1077
		c_params = ['isWgBent','upperLimitHorizontalWg','leftLimitVerticalWg','lowerLimitHorizontalWg','wgLengthX'] #all these variables must be globally declared in the scope where the "inline" function call happens
3662
1078
		#the C-code snippet for constructing the Callback object
3663
1079
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
3664
1080
		#do the actual inline C-call and fetch the pointer to the Callback object
3665
1081
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
3666
1082
		#set the pointer to the callback object in the Python-meep core 
3667
1083
		set_EPS_CallbackInlineObject(funPtr)
3668
1084
		print "EPS function successfully set."
3669
1085
		return
3670
1086
	 
3671
1087
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
3672
1088
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
3673
1089
		#in all directions,
3674
1090
		#using the material function EPS
3675
1091
		s = structure(pCompVol, EPS, pml(1.0) )     
3676
1092
		f = fields(s)      
3677
1093
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
3678
1094
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
3679
1095
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
3680
1096
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
3681
1097
		print "Field created..."
3682
1098
		return f
3683
1099
	       
3684
1100
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C++ CLASS FOR EPS\n")
3685
1101
3686
1102
	master_printf("Running on %d processor(s)...\n",count_processors())
3687
1103
3688
1104
	#FIRST WE WORK OUT THE CASE WITH NO BEND
3689
1105
	#----------------------------------------------------------------
3690
1106
	master_printf("*1* Starting the case with no bend...")
3691
1107
3692
1108
	#set EPS material function
3693
1109
	isWgBent = 0
3694
1110
	initEPS()
3695
1111
3696
1112
	#create the computational grid 
3697
1113
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
3698
1114
3699
1115
	#create the field
3700
1116
	wgBent = 0 #no bend
3701
1117
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3702
1118
3703
1119
	bendFnEps = "./bentwgNB_Eps.h5"
3704
1120
	bendFnEz = "./bentwgNB_Ez.h5"
3705
1121
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3706
1122
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
3707
1123
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
3708
1124
3709
1125
	#create the file for the field components 
3710
1126
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
3711
1127
3712
1128
	#define the flux plane for the reflected flux
3713
1129
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3714
1130
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3715
1131
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3716
1132
3717
1133
	#define the flux plane for the transmitted flux
3718
1134
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
3719
1135
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
3720
1136
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3721
1137
3722
1138
	master_printf("Calculating...")
3723
1139
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
3724
1140
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
3725
1141
	master_printf("Done..!")
3726
1142
3727
1143
	#construct 2-dimensional array with the flux plane data
3728
1144
	#see python_meep.py 
3729
1145
	noBendReflFlux = getFluxData(noBendReflectedFlux)
3730
1146
	noBendTransmFlux = getFluxData(noBendTransmFlux)
3731
1147
3732
1148
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
3733
1149
	noBendReflectedFlux.scale_dfts(-1);
3734
1150
	f = open("minusflux.h5", 'w') #truncate file if already exists
3735
1151
	f.close()   
3736
1152
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
3737
1153
3738
1154
	del_EPS_Callback() #destruct the inline-created object
3739
1155
3740
1156
3741
1157
	#AND SECONDLY FOR THE CASE WITH BEND
3742
1158
	#----------------------------------------------------------------
3743
1159
	master_printf("*2* Starting the case with bend...")
3744
1160
3745
1161
	#set EPS material function
3746
1162
	isWgBent = 1
3747
1163
	initEPS()
3748
1164
3749
1165
	#create the computational grid 
3750
1166
	bendVol = voltwo(gridSizeX,gridSizeY,res)
3751
1167
3752
1168
	#create the field
3753
1169
	wgBent = 1 #there is a bend
3754
1170
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
3755
1171
3756
1172
	#export the dielectric structure (so that we can visually verify the waveguide structure)
3757
1173
	bendFnEps = "./bentwgB_Eps.h5"
3758
1174
	bendFnEz = "./bentwgB_Ez.h5"
3759
1175
	bendDielectricFile = prepareHDF5File(bendFnEps)
3760
1176
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
3761
1177
3762
1178
	#create the file for the field components 
3763
1179
	bendFileOutputEz = prepareHDF5File(bendFnEz)
3764
1180
3765
1181
	#define the flux plane for the reflected flux
3766
1182
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
3767
1183
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
3768
1184
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
3769
1185
3770
1186
	#load the minused reflection flux from the "no bend" case as initalization
3771
1187
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
3772
1188
3773
1189
3774
1190
	#define the flux plane for the transmitted flux
3775
1191
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
3776
1192
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
3777
1193
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
3778
1194
3779
1195
	master_printf("Calculating...")
3780
1196
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
3781
1197
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
3782
1198
	master_printf("Done..!")
3783
1199
3784
1200
	#construct 2-dimensional array with the flux plane data
3785
1201
	#see python_meep.py 
3786
1202
	bendReflFlux = getFluxData(bendReflectedFlux)
3787
1203
	bendTransmFlux = getFluxData(bendTransmFlux)
3788
1204
3789
1205
	del_EPS_Callback()
3790
1206
3791
1207
	#SHOW THE RESULTS IN A PLOT
3792
1208
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
3793
1209
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
3794
1210
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
3795
1211
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
3796
1212
		                  
3797
1213
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
3798
1214
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
3799
1215
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
3800
1216
3801
1217
	matplotlib.pyplot.show()
3802
1218
3803
1219
3804
1220
3805
1221
**9. Running in MPI mode (multiprocessor configuration)**
3806
1222
----------------------------------------------------------
3807
1223
3808
1224
* We assume that an MPI implementation is installed on your machine (e.g. OpenMPI).
3809
1225
If you want to run python-meep in MPI mode, then you must must import the ``meep_mpi`` module instead of the ``meep`` module : 
3810
1226
3811
1227
``from meep_mpi import *``
3812
1228
3813
1229
Then start up the Python script as follows :
3814
1230
3815
1231
``mpirun -n 2 ./myscript.py``
3816
1232
3817
1233
The ``-n`` parameter indicates the number of processors requested.
3818
1234
3819
1235
* For printing output to the console, use the ``master_printf`` statement. This will generate output on the master node only (regular Python ``print`` statements will run on all nodes).
3820
1236
3821
1237
* If you output HDF5 files, make sure your HDF5 library is MPI-enable, otherwise your Python script will stall upon writing HDF5 due to a deadlock.
3822
1238
3823
1239
3824
1240
3825
1241
**10. Differences between Python-Meep and Scheme-Meep (libctl)**
3826
1242
------------------------------------------------------------------
3827
1243
3828
1244
**note**: this section does NOT apply to the UGent Intec Photonics Research Group (apart from the coordinate system, the default behaviour for us is made consistent with Scheme-Meep) 
3829
1245
3830
1246
The general rule is that Python-Meep has consistent behaviour with the **C++ core of Meep**.
3831
1247
The default version of Python-Meep (compiled from the LATEST_RELEASE branch) has 3 differences compared with the Scheme version of Meep :
3832
1248
	* in Python-meep, the center of the coordinate system is in the upper left corner (in Scheme-Meep v1.1.1, the center of the coordinate system is in the middle of your computational volume).
3833
1249
	* in Python-meep, eps-averaging is disabled by default (see section 3.2 for details on how to enable eps-averaging)
3834
1250
	* in Python-meep, calculation is done with complex fields by default (in Scheme-Meep v1.1.1, real fields are used by default). You can call function use_real_fields() on your fields-object to enable calculation with real fields only.
3835
1251
3836
1252
On starting your script, Python-Meep will warn you about these differences. You can suppress these warning by setting the global variables ``DISABLE_WARNING_COORDINATESYSTEM``, ``DISABLE_WARNING_EPS_AVERAGING`` and ``DISABLE_WARNING_REAL_FIELDS`` to  ``True``. You can add site-specific customisations to the file ``meep-site-init.py`` : in this script, you can for example suppress the warning, or enable EPS-averaging by default.
3837
1253
3838
1254
Add the following code to ``meep-site-init.py`` if you want *to calculate with real fields only by default* :
3839
1255
3840
1256
::
3841
1257
3842
1258
	#by default enable calculation with real fields only (consistent with Scheme-meep)
3843
1259
	def fields(structure, m = 0, store_pol_energy=0):
3844
1260
	   f = fields_orig(structure, m, store_pol_energy)
3845
1261
	   f.use_real_fields()
3846
1262
	   return f
3847
1263
3848
1264
	#add a new construct 'fields_complex' that you can use to force calculation with complex fields
3849
1265
	def fields_complex(structure, m = 0, store_pol_energy=0):
3850
1266
	   master_printf("Calculation with complex fields enalbed.\n")
3851
1267
	   return fields_orig(structure, m, store_pol_energy)
3852
1268
3853
1269
	global DISABLE_WARNING_REAL_FIELDS
3854
1270
	DISABLE_WARNING_REAL_FIELDS = True
3855
1271
3856
1272
3857
1273
Add the following code to ``meep-site-init.py`` if you want *to enable EPS-averaging by default* :
3858
1274
3859
1275
::
3860
1276
3861
1277
	use_averaging(True)
3862
1278
3863
1279
	global DISABLE_WARNING_EPS_AVERAGING
3864
1280
	DISABLE_WARNING_EPS_AVERAGING = True
3865
1281
3866
1282
3867
1283
3868
1284
3869
1285
3870
1286
3871
1287
3872
1288
3873
1289
3874
1290
3875
1291
3876
1292
3877
1293
3878
0
1294
3879
=== added directory 'doc/html/_sources/.svn/tmp'
3880
=== added directory 'doc/html/_sources/.svn/tmp/prop-base'
3881
=== added directory 'doc/html/_sources/.svn/tmp/props'
3882
=== added directory 'doc/html/_sources/.svn/tmp/text-base'
3883
=== added file 'doc/html/_sources/python_meep_documentation.txt'
3884
--- doc/html/_sources/python_meep_documentation.txt	1970-01-01 00:00:00 +0000
3885
+++ doc/html/_sources/python_meep_documentation.txt	2009-12-01 14:23:10 +0000
3886
@@ -0,0 +1,1299 @@
3887
1
PYTHON-MEEP BINDING DOCUMENTATION
3888
2
====================================
3889
3
3890
4
Primary author of this documentation : EL : Emmanuel.Lambert@intec.ugent.be 
3891
5
3892
6
Document history : 
3893
7
3894
8
::
3895
9
3896
10
	* EL-19/20/21-08-2009 : document creation 
3897
11
	* EL-24-08-2009 : small improvements & clarifications.
3898
12
	* EL-25/26-08-2009 : sections 7 & 8 were added.
3899
13
	* EL-03-04/09/2009 : 
3900
14
               -class "structure_eps_pml" (removed again in v0.8).
3901
15
	       -port to Meep 1.1.1 (class 'volume' was renamed to 'grid_volume' and class 'geometric_volume' to 'volume'
3902
16
	       -minor changes in the bent waveguide sample, to make it more consistent with the Scheme version
3903
17
	* EL-07-08/09/2009 : sections 3.2, 8.2, 8.3 : defining a material function with inline C/C++
3904
18
	* EL-10/09/2009 : additions for MPI mode (multiprocessor)
3905
19
	* EL-21/10/2009 : amplitude factor callback function
3906
20
	* EL-22/10/2009 : keyword arguments for runUntilFieldsDecayed
3907
21
	* EL-01/12/2009 : alignment with version 0.8 - III
3908
22
	* EL-01/12/2009 : release 1.0 / added info about environment variables for inline C/C++
3909
23
3910
24
3911
25
3912
26
**1. The general structure of a python-meep program**
3913
27
-----------------------------------------------------
3914
28
3915
29
In general terms, a python-meep program can be structured as follows :
3916
30
3917
31
    * import the python-meep binding :
3918
32
		``from meep import *``
3919
33
      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/``
3920
34
3921
35
      If you are running in MPI mode (multiprocessor, see section 9), then you have to import module ``meep_mpi`` instead :
3922
36
		``from meep_mpi import *``
3923
37
3924
38
    * define a computational grid volume
3925
39
        See section 2 below which explains usage of the ``grid_volume`` class. 
3926
40
3927
41
    * define the waveguide structure (describing the geometry, PML and materials)
3928
42
        See section 3 below which explains usage of the ``structure`` class.
3929
43
		
3930
44
    * create an object which will hold the calculated fields 
3931
45
        See section 4 below which explains usage of the ``field`` class.
3932
46
3933
47
    * define the sources
3934
48
        See section 5 below which explains usage of the ``add_point_source`` and ``add_volume_source`` functions.
3935
49
3936
50
    * run the simulation (iterate over the time-steps)
3937
51
        See section 6 below.
3938
52
3939
53
Section 7 gives details about defining and retrieving fluxes.
3940
54
3941
55
Section 9 gives some complete examples.
3942
56
3943
57
Section 10 outlines some differences between Scheme-Meep and Python-Meep.
3944
58
3945
59
3946
60
**2. Defining the computational grid volume**
3947
61
---------------------------------------------
3948
62
3949
63
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).
3950
64
		* ``volcyl(rsize, zsize, resolution)``
3951
65
			Defines a cyclical computational grid volume.
3952
66
		* ``volone(zsize, resolution)`` *alternatively called* ``vol1d(zsize, resolution)``
3953
67
			Defines a 1-dimensional computational grid volume along the Z-axis.
3954
68
		* ``voltwo(xsize, ysize, resolution)`` *alternatively called* ``vol2d(xsize, ysize, a)``
3955
69
			Defines a 2-dimensional computational grid volumes along the X- and Y-axes
3956
70
		* ``vol3d(xsize, ysize, zsize, resolution)``
3957
71
			Defines a 3-dimensional computational grid volume.
3958
72
3959
73
    e.g.: ``v = volone(6, 10)`` defines a 1-dimensional computational volume of lenght 6, with 10 pixels per distance unit.
3960
74
3961
75
3962
76
**3. Defining the waveguide structure**
3963
77
---------------------------------------
3964
78
	
3965
79
The waveguide structure is defined using the class ``structure``, of which the constructor has the following arguments :
3966
80
3967
81
		* *required* : the computational grid volume (a reference to an object of type ``grid_volume``, see section 2 above)
3968
82
3969
83
		* *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.
3970
84
3971
85
		* *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 :
3972
86
			- ``no_pml()`` describing a conditionless boundary region (no PML)
3973
87
			- ``pml(thickness)`` : decribing a perfectly matching layer (PML) of a certain thickness (double value) on the boundaries in all directions. 
3974
88
			- ``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``). 
3975
89
			- ``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.
3976
90
3977
91
		* *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:
3978
92
			- ``identity`` : no symmetry
3979
93
			- ``rotate4(direction, grid_volume)`` : defines a 90° rotational symmetry with 'direction' the axis of rotation.
3980
94
			- ``rotate2(direction, grid_volume)`` : defines a 180° rotational symmetry with 'direction' the axis of rotation. 
3981
95
			- ``mirror(direction, grid_volume)`` : defines a mirror symmetry plane with 'direction' normal to the mirror plane.
3982
96
			- ``r_to_minus_r_symmetry`` : defines a mirror symmetry in polar coordinates
3983
97
3984
98
		* optional: the number of chunks in which to split up the calculated geometry. If you leave this empty, it is auto-configured. Otherwise, you would set this to a factor which is a multiple of the number of processors in your MPI run (for multiprocessor configuration).
3985
99
3986
100
    e.g. : if ``v`` is a variable pointing to the computational grid volume, then :
3987
101
		 ``s = structure(v, one)`` defines a structure with all air (eps=1), 
3988
102
    which is equivalent to:
3989
103
		 ``s = structure(v, one, no_pml(), identity(), 1)``
3990
104
3991
105
    Another example : ``s = structure(v, EPS, pml(0.1,Y) )`` with EPS a custom material function, which is explained in the note below.
3992
106
3993
107
3994
108
3.1. Defining a material function
3995
109
________________________________________
3996
110
3997
111
In order to describe the geometry of the waveguide, we have to provide a 'material function' returning 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.
3998
112
It is also possible (and faster) to write your material function in inline C/C++ (see 3.3)
3999
113
4000
114
E.g. :
4001
115
4002
116
::
4003
117
4004
118
        class epsilon(Callback):  #inherit from Callback for integration with the meep core library
4005
119
                def __init__(self):
4006
120
                        Callback.__init__(self)
4007
121
                def double_vec(self,vec): #override of function in the Callback class to set the eps function
4008
122
                        self.set_double(self.eps(vec))
4009
123
                        return 
4010
124
                def eps(self,vec):   #return the epsilon value for a certain point (indicated by the vector v)
4011
125
                        v = vec 
4012
126
                        r = v.x()*v.x() + v.y()*v.y()
4013
127
                        dr = sqrt(r)
4014
128
                        while dr>1:
4015
129
                            dr-=1
4016
130
                        if dr > 0.7001:
4017
131
                            return 12.0
4018
132
                        return 1.0
4019
133
4020
134
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.
4021
135
So, one should write : ``vec.x()`` and ``vec.y()``.
4022
136
 
4023
137
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 :
4024
138
4025
139
::
4026
140
4027
141
        set_EPS_Callback(epsilon().__disown__())
4028
142
        s = structure(v, EPS, no_pml(), identity())
4029
143
4030
144
4031
145
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.
4032
146
4033
147
***Important remark*** : at the end of our program, we should call : ``del_EPS_Callback()`` in order to clean up the global variable.
4034
148
4035
149
For better performance, you can define your EPS material function with inline C/C++ : we refer to section 3.3 for details about this.
4036
150
4037
151
3.2 Eps-averaging
4038
152
_________________
4039
153
4040
154
EPS-averaging (anisotrpic averaging) is disabled by default, making this behaviour consistent with the behaviour of the Meep C++ core.
4041
155
4042
156
You can enable EPS-averaging using the function ``use_averaging`` :
4043
157
4044
158
::
4045
159
4046
160
	#enable EPS-averaging
4047
161
	use_averaging(True)
4048
162
	...
4049
163
	#disable EPS-averaging
4050
164
	use_averaging(False)
4051
165
	...
4052
166
4053
167
4054
168
Enabling EPS-averaging results in slower performance, but more accurate results. 
4055
169
4056
170
4057
171
3.3. Defining a material function with inline C/C++
4058
172
_________________________________________________________
4059
173
4060
174
The approach described in 3.1 lets the Meep core library call back to Python for every query of the epsilon-function. This creates a lot of overhead.
4061
175
An approach which has a lot better performance is to define this epsilon-function with an inline C-function or C++ class.
4062
176
4063
177
* If our epsilon-function needs *no other parameters than the position vector (X, Y, Z)*, then we can suffice with an inline C-function (the geometry dependencies are then typically hardcoded).
4064
178
4065
179
* If our epsilon-function needs to base it's calculation on *a more complex set of parameters (e.g. parameters depending on the geometry)*, then we have to write a C++ class. 
4066
180
4067
181
For example, in the bent-waveguide example (section 8.3), we can define a generic C++ class which can return the epsilon-value for both the "bend" and "no bend" case, with variable size parameters.
4068
182
We can also take a simpler approach (section 8.2) and write a function in which the geometry size parameters are hardcoded : we then need 2 inline C-functions : one for the "bend" case and one for the "no bend" case. 
4069
183
4070
184
Make sure the following environment variables are defined :
4071
185
	* MEEP_INCLUDE : should point to the path containing meep.hpp (and a subdirectory 'meep' with vec.hpp and mympi.hpp), e.g.: ``/usr/include``
4072
186
	* MEEP_LIB : should point to the path containing libmeep.so, e.g. : ``/usr/lib``
4073
187
	* PYTHONMEEP_INCLUDE : should point to the path containing custom.hpp, e.g. : ``/usr/local/lib/python2.6/dist-packages``
4074
188
	
4075
189
4076
190
3.3.1 Inline C-function
4077
191
.......................
4078
192
4079
193
First we create a header file, e.g. "eps_function.hpp" which contains our EPS-function.
4080
194
Not that the geometry dependencies are hardcoded (``upperLimitHorizontalWg = 4`` and ``lowerLimitHorizontalWg = 5``).
4081
195
4082
196
::
4083
197
4084
198
4085
199
	namespace meep
4086
200
	{
4087
201
		static double myEps(const vec &v) {
4088
202
			double xCo = v.x();
4089
203
			double yCo = v.y();
4090
204
			double upperLimitHorizontalWg = 4;
4091
205
			double lowerLimitHorizontalWg = 5;
4092
206
			if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
4093
207
				return 1.0;
4094
208
			}
4095
209
			else return 12.0;
4096
210
		}
4097
211
	}
4098
212
4099
213
4100
214
Then, in the Python program, we prepare and set the callback function as shown below.
4101
215
``prepareCallbackCfunction`` returns a pointer to the C-function, which we deliver to the Meep core using ``set_EPS_CallbackInlineFunction``.
4102
216
4103
217
::
4104
218
4105
219
	def initEPS(isWgBent):
4106
220
		funPtr = prepareCallbackCfunction("myEps","eps_function.hpp")  #name of your function / name of header file
4107
221
		set_EPS_CallbackInlineFunction(funPtr)
4108
222
		print "EPS function successfully set."
4109
223
		return funPtr
4110
224
4111
225
We refer to section 8.2 below for a full example.
4112
226
4113
227
4114
228
3.3.2 Inline C++-class
4115
229
......................
4116
230
4117
231
A more complex approach is to have a C++ object that can accept more parameters when it is constructed. 
4118
232
For example this is the case if want to change the parameters of the geometry from inside Python without touching the C++ code.
4119
233
4120
234
We create a header file "eps_class.hpp" which contains the definition of the class (the class must inherit from ``Callback``).
4121
235
In the example below, the parameters ``upperLimitHorizontalWg`` and ``widthWg`` will be communicated from Python upon construction of the object.
4122
236
If these parameters then change (depending on the geometry), the C++ object will follow automatically.
4123
237
4124
238
4125
239
::
4126
240
4127
241
	using namespace meep;
4128
242
4129
243
	namespace meep
4130
244
	{
4131
245
4132
246
	class myEpsCallBack : virtual public Callback {
4133
247
4134
248
	    public:	
4135
249
		   myEpsCallBack() : Callback() { };
4136
250
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
4137
251
4138
252
		   myEpsCallBack(double upperLimitHorizontalWg, double widthWg)  : Callback() {
4139
253
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
4140
254
			_widthWg = widthWg;
4141
255
		   };	
4142
256
4143
257
		   double double_vec(const vec &v) { //return the EPS-value, depending on the position vector
4144
258
			double eps = myEps(v, _upperLimitHorizontalWg, _widthWg);
4145
259
			return eps;
4146
260
		   };
4147
261
		   
4148
262
		   complex<double> complex_vec(const vec &x) { return 0; }; //no need to implement
4149
263
		   complex<double> complex_time(const double &t) { return 0; }; //no need to implement 
4150
264
4151
265
4152
266
	     private:
4153
267
		   double _upperLimitHorizontalWg;
4154
268
		   double _widthWg;
4155
269
4156
270
		   double myEps(const vec &v, double upperLimitHorizontalWg, double widthWg) {
4157
271
			double xCo = v.x();
4158
272
			double yCo = v.y();
4159
273
	    	        if ((yCo < upperLimitHorizontalWg) || (yCo > upperLimitHorizontalWg+widthWg)){
4160
274
				return 1.0;
4161
275
			    }
4162
276
			}
4163
277
			return 12.0;
4164
278
		   }
4165
279
4166
280
	};
4167
281
4168
282
	}
4169
283
4170
284
4171
285
The syntax in Python is a little bit more complex in this case.
4172
286
4173
287
We will need to import the module ``scipy.weave`` : 
4174
288
4175
289
``from scipy.weave import *``
4176
290
4177
291
(this is not required for the previous case of a simple inline function)
4178
292
4179
293
First we create a list with the names of the parameters that we want to pass to the C++ class. These variables must be declared in the scope where the ``inline`` function call happens (see below).
4180
294
4181
295
``c_params = ['upperLimitHorizontalWg','widthWg']``
4182
296
4183
297
Then, we prepare the code snippet, using the function ``prepareCallbackCObjectCode`` and passing the class name and parameter names list.
4184
298
4185
299
``c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)``
4186
300
4187
301
Finally, we call the ``inline`` function, passing :
4188
302
	* the code snippet
4189
303
	* the list of parameter names
4190
304
	* the inline libraries, include directories and headers (helper functions are provided for this, see below). The call to ``getInlineHeaders`` should receive the name of the header file (with the definition of the C++ class) as an argument.
4191
305
4192
306
``funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )``
4193
307
4194
308
::
4195
309
4196
310
4197
311
	def initEPS():
4198
312
		#the set of parameters that we want to pass to the Callback object upon construction
4199
313
		#all these variables must be declared in the scope where the "inline" function call happens
4200
314
		c_params = ['upperLimitHorizontalWg','widthWg'] 
4201
315
		#the C-code snippet for constructing the Callback object
4202
316
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
4203
317
		#do the actual inline C-call and fetch the pointer to the Callback object
4204
318
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
4205
319
		#set the pointer to the callback object in the Python-meep core 
4206
320
		set_EPS_CallbackInlineObject(funPtr)
4207
321
		print "EPS function successfully set."
4208
322
		return
4209
323
 
4210
324
4211
325
We refer to section 8.3 below for a full example.
4212
326
4213
327
4214
328
4215
329
**4. Defining the initial field**
4216
330
---------------------------------
4217
331
4218
332
This is optional.
4219
333
4220
334
We create an object of type ``fields``, which will contain the calculated field.
4221
335
4222
336
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``.
4223
337
4224
338
::
4225
339
4226
340
            class fi(Callback):  #inherit from Callback for integration with the meep core library
4227
341
                def __init__(self):
4228
342
                    Callback.__init__(self)		
4229
343
                def complex_vec(self,v): #override of function in the Callback class to set the field initialization function
4230
344
		     #return the field value for a certain point (indicated by the vector v)
4231
345
                    return complex(1.0,0)
4232
346
4233
347
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 : 
4234
348
4235
349
	``set_INIF_Callback(fi().__disown__())  #link the INIF variable to the fi class``
4236
350
4237
351
We refer to section 3-note1 for more information about the function ``__disown__()``.
4238
352
4239
353
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 :
4240
354
4241
355
::
4242
356
4243
357
    f = fields(s)
4244
358
    comp = Hy
4245
359
    f.initialize_field(comp, INIF)
4246
360
4247
361
4248
362
***Important remark*** : at the end of our program, we should then call : ``del_INIF_Callback()`` in order to clean up the global variable.
4249
363
4250
364
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.
4251
365
4252
366
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. :
4253
367
4254
368
	``f.use_bloch(vec(0.0))``
4255
369
4256
370
*to be further elaborated - what is the exact meaning of the vector argument? (not well understood at this time)*
4257
371
4258
372
4259
373
4260
374
**5. Defining the sources**
4261
375
---------------------------
4262
376
4263
377
The definition of the current sources can be done through 2 functions of the ``fields`` class :
4264
378
        * ``add_point_source(component, src_time, vec, complex)``
4265
379
	* ``add_volume_source(component, src_time, volume, complex)``
4266
380
4267
381
4268
382
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).
4269
383
4270
384
For a point source, we must specify the center point of the current source using a vector (object of type ``vec``).
4271
385
4272
386
For a volume source, we must specify an object of type ``volume`` (*to be elablorated*).
4273
387
4274
388
The last argument is an overall complex amplitude number, multiplying the current source (default 1.0).
4275
389
4276
390
The following variants are available :
4277
391
        * ``add_point_source(component, double, double, double, double, vec centerpoint, complex amplitude, int is_continuous)``
4278
392
                * This is a shortcut function so that no ``src_time`` object must be created. *This function is preferably used for point sources.*
4279
393
                * The four real arguments define the central frequency, spectral width, peaktime and cutoff. 
4280
394
4281
395
        * ``add_volume_source(component, src_time, volume)``
4282
396
4283
397
        * ``add_volume_source(component, src_time, volume, AMPL)``
4284
398
                * AMPL is a built-in reference to a callback function. Such a callback function returns a factor to multiply the source amplitude with (complex value). It receives 1 parameter, i.e. a vector indicating a position RELATIVE to the CENTER of the source. See the example below.
4285
399
4286
400
4287
401
Three classes, inheriting from ``src_time``, are predefined and can be used off the shelf :
4288
402
	* ``gaussian_src_time`` for a Gaussian-pulse source. The constructor demands 2 arguments of type double :
4289
403
                    * the center frequency ω, in units of 2πc 
4290
404
		    * the frequency width w used in the Gaussian
4291
405
	* ``continuous_src_time`` for a continuous-wave source proportional to exp(−iωt). The constructor demands 4 arguments : 
4292
406
                    * the frequency ω, in units 2πc/distance (complex number)
4293
407
		    * the temporal width of smoothing (default 0)
4294
408
                    * the start time (default 0)
4295
409
                    * the end time (default infinity = never turn off)
4296
410
	* ``custom_src_time`` for a user-specified source function f(t) with start/end times. The constructor demands 4 arguments :
4297
411
		    * The function f(t) specifying the time-dependence of the source
4298
412
                    * *...(2nd argument unclear, to be further elaborated)...*
4299
413
                    * the start time of the source (default -infinity)
4300
414
                    * the end time of the source (default +infinity)    
4301
415
4302
416
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 :
4303
417
           
4304
418
::
4305
419
4306
420
	#define a continuous source 
4307
421
	srcFreq = 0.125
4308
422
	srcWidth = 20 
4309
423
	src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
4310
424
	srcComp = Ez
4311
425
	#make it a line source of size 1 starting on position(6,3)
4312
426
	srcGeo = volume(vec(6,3),vec(6,4))
4313
427
	f.add_volume_source(srcComp, src, srcGeo)
4314
428
4315
429
4316
430
Here is an example of the implementation of a callback function for the amplitude factor :
4317
431
4318
432
::
4319
433
4320
434
	class amplitudeFactor(Callback):
4321
435
	    def __init__(self):
4322
436
		Callback.__init__(self)
4323
437
		master_printf("Callback function for amplitude factor activated.\n")
4324
438
4325
439
	    def complex_vec(self,vec):
4326
440
		#BEWARE, these are coordinates RELATIVE to the source center !!!!
4327
441
		x = vec.x()
4328
442
		y = vec.y()
4329
443
		master_printf("Fetching amplitude factor for x=%f - y=%f\n" %(x,y) )
4330
444
		result = complex(1.0,0)
4331
445
		return result
4332
446
4333
447
	...
4334
448
		#define a continuous source 
4335
449
		srcFreq = 0.125
4336
450
		srcWidth = 20 
4337
451
		src = continuous_src_time(srcFreq, srcWidth, 0, infinity)
4338
452
		srcComp = Ez
4339
453
		#make it a line source of size 1 starting on position(6,3)
4340
454
		srcGeo = volume(vec(6,3),vec(6,4))
4341
455
		#create callback object for amplitude factor
4342
456
		af = amplitudeFactor()
4343
457
		set_AMPL_Callback(af.__disown__())
4344
458
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, AMPL)
4345
459
4346
460
4347
461
**6. Running the simulation, retrieving field values and exporting HDF5 files**
4348
462
-------------------------------------------------------------------------------
4349
463
4350
464
We can now time-step and retrieve various field values along the way.
4351
465
The actual time step value can be retrieved or set through the variable ``f.dt``.
4352
466
4353
467
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.
4354
468
4355
469
To trigger a step in time, you call the function ``f.step()``.
4356
470
4357
471
To step until the source has fully decayed :
4358
472
4359
473
::
4360
474
4361
475
        while (f.time() < f.last_source_time()):
4362
476
             f.step()
4363
477
4364
478
The function ``runUntilFieldsDecayed`` mimicks the behaviour of 'stop-when-fields-decayed' in Meep-Scheme.
4365
479
This will run time steps until the source has decayed to 0.001 of the peak amplitude. After that, by default an additional 50 time steps will be run.
4366
480
The function has 7 arguments, of which 4 are mandatory and 3 are optional keywords arguments :
4367
481
* 4 regular arguments : reference to the field, reference to the computational grid volume, the source component, the monitor point.
4368
482
* keyword argument ``pHDF5OutputFile`` : reference to a HDF5 file (constructed with the function ``prepareHDF5File``); default : None (no ouput to files)
4369
483
* keyword argument ``pH5OutputIntervalSteps`` : step interval for output to HDF5 (default : 50)
4370
484
* keyword argument ``pDecayedStopAfterSteps`` : the number of steps to continue after the source has decayed to 0.001 of the peak amplitude at the probing point (default: 50)
4371
485
 
4372
486
We further refer to section 8 below where this function is applied in an example.
4373
487
4374
488
A rich functionality is available for retrieving field information. Some examples :
4375
489
4376
490
	* ``f.energy_in_box(v.surroundings())``
4377
491
	* ``f.electric_energy_in_box(v.surroundings())``
4378
492
	* ``f.magnetic_energy_in_box(v.surroundings())``
4379
493
	* ``f.thermo_energy_in_box(v.surroundings())``
4380
494
	* ``f.total_energy()``
4381
495
	* ``f.field_energy_in_box(v.surroundings())``
4382
496
	* ``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``)
4383
497
	* ``f.flux_in_box(X, v.surroundings())`` where the first argument is the direction (``X, Y, Z, R`` or ``P``)
4384
498
4385
499
We can probe the field at certain points by defining a *monitor point* as follows : 
4386
500
4387
501
::
4388
502
4389
503
        m = monitor_point()
4390
504
        p = vec(2.10) #vector identifying the point that we want to probe
4391
505
        f.get_point(m, p)
4392
506
        m.get_component(Hx)
4393
507
4394
508
We can export the dielectric function and e.g. the Ex component of the field to HDF5 files as follows :
4395
509
4396
510
::
4397
511
4398
512
        #make sure you start your Python session with 'sudo' or write rights to the current path
4399
513
        feps = prepareHDF5File("eps.h5")
4400
514
        f.output_hdf5(Dielectric, v.surroundings(), feps)   #export the Dielectric structure so that we can visually verify it 
4401
515
	fex = prepareHDF5File("ex.h5")
4402
516
	while (f.time() < f.last_source_time()):
4403
517
             f.step()
4404
518
             f.output_hdf5(Ex, v.surroundings(), fex, 1) #export the Ex component, appending to the file "ex.h5"
4405
519
4406
520
4407
521
4408
522
**7. Defining and retrieving fluxes**
4409
523
--------------------------------------
4410
524
4411
525
First we define a flux plane. 
4412
526
This is done through the creation of an object of type ``volume`` (specifying 2 vectors as arguments).
4413
527
4414
528
Then we apply this flux plane to the field, specifying 4 parameters :
4415
529
	* the reference to the ``volume`` object
4416
530
	* the minimum frequency (in the example below, this is ``srcFreqCenter-(srcPulseWidth/2.0)``)
4417
531
	* the maximum frequency (in the example below this is ``srcFreqCenter+(srcPulseWidth/2.0)`` )
4418
532
	* the number of discrete frequencies that we want to monitor in the flux (in the example below, this is 100).
4419
533
4420
534
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.
4421
535
4422
536
E.g., if ``f`` is the field, then we proceed as follows : 
4423
537
4424
538
::
4425
539
4426
540
	#define the flux plane and flux parameters
4427
541
	fluxplane = volume(vec(1,2),vec(1,3))
4428
542
	flux = f.add_dft_flux_plane(fluxplane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
4429
543
4430
544
	#run the calculation 
4431
545
	while (f.time() < f.last_source_time()):
4432
546
	    f.step()
4433
547
4434
548
	#retrieve the flux data
4435
549
	fluxdata = getFluxData(flux)
4436
550
	frequencies = fluxdata[0]
4437
551
	fluxvalues = fluxdata[1]
4438
552
4439
553
4440
554
4441
555
**8. The "90 degree bent waveguide example in Python**
4442
556
------------------------------------------------------
4443
557
4444
558
We have ported the "90 degree bent waveguide" example from the Meep-Scheme tutorial to Python.
4445
559
4446
560
The original example can be found on the following URL : http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial
4447
561
(section 'A 90° bend').
4448
562
4449
563
You can find the source code below and also in the ``/samples/bent_waveguide`` directory of the Python-Meep distribution. 
4450
564
4451
565
Sections 8.2 and 8.3 contain the same example but with the EPS-callback function as inline C-function and inline C++-class.
4452
566
4453
567
The following animated gifs can be produced from the HDF5-files (see the script included in directory 'samples') :
4454
568
4455
569
.. image:: images/bentwgNB.gif
4456
570
4457
571
.. image:: images/bentwgB.gif
4458
572
4459
573
4460
574
And here is the graph of the transmission, reflection and loss fluxes, showing the same results as the example in Scheme:
4461
575
4462
576
.. image:: images/fluxes.png
4463
577
   :height: 315
4464
578
   :width: 443
4465
579
4466
580
4467
581
8.1 With a Python class as EPS-function
4468
582
________________________________________
4469
583
4470
584
4471
585
A bottleneck in this version is the epsilon-function, which is written in Python. 
4472
586
This means that the Meep core library must do a callback to the Python function, which creates a lot of overhead.
4473
587
An approach which has a much better performance is to write this EPS-function in C : the Meep core library can then directly call back to a C-function.
4474
588
These approaches are described in 8.2 and 8.3.
4475
589
4476
590
::
4477
591
4478
592
	from meep import *
4479
593
	from math import *
4480
594
	from python_meep import *
4481
595
	import numpy 
4482
596
	import matplotlib.pyplot
4483
597
	import sys
4484
598
4485
599
	res = 10.0
4486
600
	gridSizeX = 16.0
4487
601
	gridSizeY = 32.0
4488
602
	padSize = 4.0
4489
603
	wgLengthX = gridSizeX - padSize
4490
604
	wgLengthY = gridSizeY - padSize
4491
605
	wgWidth = 1.0 #width of the waveguide
4492
606
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
4493
607
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
4494
608
	srcFreqCenter = 0.15 #gaussian source center frequency
4495
609
	srcPulseWidth = 0.1 #gaussian source pulse width
4496
610
	srcComp = Ez #gaussian source component
4497
611
4498
612
	#this function plots the waveguide material as a function of a vector(X,Y) 
4499
613
	class epsilon(Callback):
4500
614
	    def __init__(self, pIsWgBent):
4501
615
		Callback.__init__(self)
4502
616
		self.isWgBent = pIsWgBent
4503
617
	    def double_vec(self,vec):
4504
618
		if (self.isWgBent): #there is a bend
4505
619
		    if ((vec.x()<wgLengthX) and (vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
4506
620
		        return 12.0
4507
621
		    elif ((vec.x()>=wgLengthX-wgWidth) and (vec.x()<=wgLengthX) and vec.y()>= padSize ):
4508
622
		        return 12.0
4509
623
		    else:
4510
624
		        return 1.0
4511
625
		else: #there is no bend 
4512
626
		    if ((vec.y() >= padSize) and (vec.y() <= padSize+wgWidth)):
4513
627
		        return 12.0
4514
628
		    else:
4515
629
		        return 1.0
4516
630
4517
631
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
4518
632
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
4519
633
		#in all directions,
4520
634
		#using the material function EPS
4521
635
		material = epsilon(pIsWgBent)
4522
636
		set_EPS_Callback(material.__disown__())
4523
637
		s = structure(pCompVol, EPS, pml(1.0) )     
4524
638
		f = fields(s)      
4525
639
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
4526
640
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
4527
641
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
4528
642
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
4529
643
		print "Field created..."
4530
644
		return f
4531
645
	       
4532
646
4533
647
	#FIRST WE WORK OUT THE CASE WITH NO BEND
4534
648
	#----------------------------------------------------------------
4535
649
	print "*1* Starting the case with no bend..."
4536
650
	#create the computational grid 
4537
651
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
4538
652
4539
653
	#create the field
4540
654
	wgBent = 0 #no bend
4541
655
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
4542
656
4543
657
	bendFnEps = "./bentwgNB_Eps.h5"
4544
658
	bendFnEz = "./bentwgNB_Ez.h5"
4545
659
	#export the dielectric structure (so that we can visually verify the waveguide structure)
4546
660
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
4547
661
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
4548
662
4549
663
	#create the file for the field components 
4550
664
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
4551
665
4552
666
	#define the flux plane for the reflected flux
4553
667
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
4554
668
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
4555
669
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
4556
670
4557
671
	#define the flux plane for the transmitted flux
4558
672
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
4559
673
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
4560
674
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
4561
675
4562
676
	print "Calculating..."
4563
677
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
4564
678
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
4565
679
	print "Done..!"
4566
680
4567
681
	#construct 2-dimensional array with the flux plane data
4568
682
	#see python_meep.py 
4569
683
	noBendReflFlux = getFluxData(noBendReflectedFlux)
4570
684
	noBendTransmFlux = getFluxData(noBendTransmFlux)
4571
685
4572
686
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
4573
687
	noBendReflectedFlux.scale_dfts(-1);
4574
688
	f = open("minusflux.h5", 'w') #truncate file if already exists
4575
689
	f.close()   
4576
690
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
4577
691
4578
692
	del_EPS_Callback()
4579
693
4580
694
4581
695
	#AND SECONDLY FOR THE CASE WITH BEND
4582
696
	#----------------------------------------------------------------
4583
697
	print "*2* Starting the case with bend..."
4584
698
	#create the computational grid 
4585
699
	bendVol = voltwo(gridSizeX,gridSizeY,res)
4586
700
4587
701
	#create the field
4588
702
	wgBent = 1 #there is a bend
4589
703
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
4590
704
4591
705
	#export the dielectric structure (so that we can visually verify the waveguide structure)
4592
706
	bendFnEps = "./bentwgB_Eps.h5"
4593
707
	bendFnEz = "./bentwgB_Ez.h5"
4594
708
	bendDielectricFile = prepareHDF5File(bendFnEps)
4595
709
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
4596
710
4597
711
	#create the file for the field components 
4598
712
	bendFileOutputEz = prepareHDF5File(bendFnEz)
4599
713
4600
714
	#define the flux plane for the reflected flux
4601
715
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
4602
716
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
4603
717
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
4604
718
4605
719
	#load the minused reflection flux from the "no bend" case as initalization
4606
720
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
4607
721
4608
722
4609
723
	#define the flux plane for the transmitted flux
4610
724
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
4611
725
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
4612
726
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
4613
727
4614
728
	print "Calculating..."
4615
729
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
4616
730
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
4617
731
	print "Done..!"
4618
732
4619
733
	#construct 2-dimensional array with the flux plane data
4620
734
	#see python_meep.py 
4621
735
	bendReflFlux = getFluxData(bendReflectedFlux)
4622
736
	bendTransmFlux = getFluxData(bendTransmFlux)
4623
737
4624
738
	del_EPS_Callback()
4625
739
4626
740
	#SHOW THE RESULTS IN A PLOT
4627
741
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
4628
742
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
4629
743
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
4630
744
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
4631
745
		                  
4632
746
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
4633
747
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
4634
748
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
4635
749
4636
750
	matplotlib.pyplot.show()
4637
751
4638
752
4639
753
8.2 With an inline C-function as EPS-function
4640
754
______________________________________________
4641
755
4642
756
The header file "eps_function.hpp" :
4643
757
4644
758
::
4645
759
4646
760
	using namespace meep;
4647
761
4648
762
	namespace meep
4649
763
	{
4650
764
	static double myEps(const vec &v, bool isWgBent) {
4651
765
		double xCo = v.x();
4652
766
		double yCo = v.y();
4653
767
		double upperLimitHorizontalWg = 4;
4654
768
		double wgLengthX = 12;
4655
769
		double leftLimitVerticalWg = 11;
4656
770
		double lowerLimitHorizontalWg = 5;
4657
771
		if (isWgBent){ //there is a bend
4658
772
		    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
4659
773
			return 1.0;            
4660
774
		    }
4661
775
		    else {
4662
776
			if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
4663
777
			     return 1.0;
4664
778
			}
4665
779
			else {
4666
780
			     return 12.0;
4667
781
			}
4668
782
		    }
4669
783
		}
4670
784
		else { //there is no bend 
4671
785
		    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
4672
786
			return 1.0;
4673
787
		    }
4674
788
		}
4675
789
		return 12.0;
4676
790
	   }
4677
791
4678
792
	static double myEpsBentWg(const vec &v) {
4679
793
		return myEps(v, true);
4680
794
	   }
4681
795
4682
796
	static double myEpsStraightWg(const vec &v) {
4683
797
		return myEps(v, false);
4684
798
	   }
4685
799
	}
4686
800
4687
801
4688
802
4689
803
And the actual Python program :
4690
804
4691
805
4692
806
::
4693
807
4694
808
4695
809
	from meep import *  
4696
810
4697
811
	from math import *
4698
812
	import numpy 
4699
813
	import matplotlib.pyplot
4700
814
	import sys
4701
815
4702
816
	res = 10.0
4703
817
	gridSizeX = 16.0
4704
818
	gridSizeY = 32.0
4705
819
	padSize = 4.0
4706
820
	wgLengthX = gridSizeX - padSize
4707
821
	wgLengthY = gridSizeY - padSize
4708
822
	wgWidth = 1.0 #width of the waveguide
4709
823
	upperLimitHorizontalWg = padSize
4710
824
	lowerLimitHorizontalWg = padSize+wgWidth
4711
825
	leftLimitVerticalWg = wgLengthX-wgWidth
4712
826
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
4713
827
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
4714
828
	srcFreqCenter = 0.15 #gaussian source center frequency
4715
829
	srcPulseWidth = 0.1 #gaussian source pulse width
4716
830
	srcComp = Ez #gaussian source component
4717
831
4718
832
	def initEPS(isWgBent):
4719
833
		if (isWgBent):
4720
834
		    funPtr = prepareCallbackCfunction("myEpsBentWg","eps_function.hpp")
4721
835
		else:
4722
836
		    funPtr = prepareCallbackCfunction("myEpsStraightWg","eps_function.hpp")
4723
837
		set_EPS_CallbackInlineFunction(funPtr)
4724
838
		print "EPS function successfully set."
4725
839
		return funPtr
4726
840
	 
4727
841
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
4728
842
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
4729
843
		#in all directions,
4730
844
		#using the material function EPS
4731
845
		s = structure(pCompVol, EPS, pml(1.0) )     
4732
846
		f = fields(s)      
4733
847
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
4734
848
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
4735
849
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
4736
850
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
4737
851
		print "Field created..."
4738
852
		return f
4739
853
4740
854
4741
855
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C-FUNCTION FOR EPS\n")
4742
856
	       
4743
857
	master_printf("Running on %d processor(s)...\n",count_processors())
4744
858
4745
859
	#FIRST WE WORK OUT THE CASE WITH NO BEND
4746
860
	#----------------------------------------------------------------
4747
861
	master_printf("*1* Starting the case with no bend...")
4748
862
4749
863
	#set EPS material function
4750
864
	initEPS(0)
4751
865
4752
866
	#create the computational grid 
4753
867
	noBendVol = voltwo(gridSizeX,gridSizeY,res)
4754
868
4755
869
	#create the field
4756
870
	wgBent = 0 #no bend
4757
871
	noBendField = createField(noBendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
4758
872
4759
873
	bendFnEps = "./bentwgNB_Eps.h5"
4760
874
	bendFnEz = "./bentwgNB_Ez.h5"
4761
875
	#export the dielectric structure (so that we can visually verify the waveguide structure)
4762
876
	noBendDielectricFile =  prepareHDF5File(bendFnEps)
4763
877
	noBendField.output_hdf5(Dielectric, noBendVol.surroundings(), noBendDielectricFile)
4764
878
4765
879
	#create the file for the field components 
4766
880
	noBendFileOutputEz = prepareHDF5File(bendFnEz)
4767
881
4768
882
	#define the flux plane for the reflected flux
4769
883
	noBendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
4770
884
	noBendReflectedFluxPlane = volume(vec(noBendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
4771
885
	noBendReflectedFlux = noBendField.add_dft_flux_plane(noBendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
4772
886
4773
887
	#define the flux plane for the transmitted flux
4774
888
	noBendTransmfluxPlaneXPos = gridSizeX - 1.5;  #the X-coordinate of our transmission flux plane
4775
889
	noBendTransmFluxPlane = volume(vec(noBendTransmfluxPlaneXPos,wgHorYCen-wgWidth),vec(noBendTransmfluxPlaneXPos,wgHorYCen+wgWidth))
4776
890
	noBendTransmFlux = noBendField.add_dft_flux_plane(noBendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
4777
891
4778
892
	master_printf("Calculating...")
4779
893
	noBendProbingpoint = vec(noBendTransmfluxPlaneXPos,wgHorYCen) #the point at the end of the waveguide that we want to probe to check if source has decayed
4780
894
	runUntilFieldsDecayed(noBendField, noBendVol, srcComp, noBendProbingpoint, noBendFileOutputEz)
4781
895
	master_printf("Done..!")
4782
896
4783
897
	#construct 2-dimensional array with the flux plane data
4784
898
	#see python_meep.py 
4785
899
	noBendReflFlux = getFluxData(noBendReflectedFlux)
4786
900
	noBendTransmFlux = getFluxData(noBendTransmFlux)
4787
901
4788
902
	#save the reflection flux from the "no bend" case as minus flux in the temporary file 'minusflux.h5'
4789
903
	noBendReflectedFlux.scale_dfts(-1);
4790
904
	f = open("minusflux.h5", 'w') #truncate file if already exists
4791
905
	f.close()   
4792
906
	noBendReflectedFlux.save_hdf5(noBendField, "minusflux")
4793
907
4794
908
	del_EPS_Callback() #destruct the inline-created object
4795
909
4796
910
4797
911
	#AND SECONDLY FOR THE CASE WITH BEND
4798
912
	#----------------------------------------------------------------
4799
913
	master_printf("*2* Starting the case with bend...")
4800
914
4801
915
	#set EPS material function
4802
916
	initEPS(1)
4803
917
4804
918
	#create the computational grid 
4805
919
	bendVol = voltwo(gridSizeX,gridSizeY,res)
4806
920
4807
921
	#create the field
4808
922
	wgBent = 1 #there is a bend
4809
923
	bendField = createField(bendVol, wgLengthX, wgWidth, wgBent, srcFreqCenter, srcPulseWidth, srcComp)
4810
924
4811
925
	#export the dielectric structure (so that we can visually verify the waveguide structure)
4812
926
	bendFnEps = "./bentwgB_Eps.h5"
4813
927
	bendFnEz = "./bentwgB_Ez.h5"
4814
928
	bendDielectricFile = prepareHDF5File(bendFnEps)
4815
929
	bendField.output_hdf5(Dielectric, bendVol.surroundings(), bendDielectricFile)
4816
930
4817
931
	#create the file for the field components 
4818
932
	bendFileOutputEz = prepareHDF5File(bendFnEz)
4819
933
4820
934
	#define the flux plane for the reflected flux
4821
935
	bendReflectedfluxPlaneXPos = 1.5 #the X-coordinate of our reflection flux plane
4822
936
	bendReflectedFluxPlane = volume(vec(bendReflectedfluxPlaneXPos,wgHorYCen-wgWidth),vec(bendReflectedfluxPlaneXPos,wgHorYCen+wgWidth))
4823
937
	bendReflectedFlux = bendField.add_dft_flux_plane(bendReflectedFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100)
4824
938
4825
939
	#load the minused reflection flux from the "no bend" case as initalization
4826
940
	bendReflectedFlux.load_hdf5(bendField, "minusflux")
4827
941
4828
942
4829
943
	#define the flux plane for the transmitted flux
4830
944
	bendTransmfluxPlaneYPos = padSize + wgLengthY - 1.5; #the Y-coordinate of our transmission flux plane
4831
945
	bendTransmFluxPlane = volume(vec(wgVerXCen - wgWidth,bendTransmfluxPlaneYPos),vec(wgVerXCen + wgWidth,bendTransmfluxPlaneYPos))
4832
946
	bendTransmFlux = bendField.add_dft_flux_plane(bendTransmFluxPlane, srcFreqCenter-(srcPulseWidth/2.0), srcFreqCenter+(srcPulseWidth/2.0), 100 )
4833
947
4834
948
	master_printf("Calculating...")
4835
949
	bendProbingpoint = vec(wgVerXCen,bendTransmfluxPlaneYPos) #the point at the end of the waveguide that we want to probe to check if source has decayed
4836
950
	runUntilFieldsDecayed(bendField, bendVol, srcComp, bendProbingpoint, bendFileOutputEz)
4837
951
	master_printf("Done..!")
4838
952
4839
953
	#construct 2-dimensional array with the flux plane data
4840
954
	#see python_meep.py 
4841
955
	bendReflFlux = getFluxData(bendReflectedFlux)
4842
956
	bendTransmFlux = getFluxData(bendTransmFlux)
4843
957
4844
958
	del_EPS_Callback()
4845
959
4846
960
	#SHOW THE RESULTS IN A PLOT
4847
961
	frequencies = bendReflFlux[0] #should be equal to bendTransmFlux.keys() or noBendTransmFlux.keys() or ...
4848
962
	Ptrans = [x / y for x,y in zip(bendTransmFlux[1], noBendTransmFlux[1])] 
4849
963
	Prefl = [ abs(x / y) for x,y in zip(bendReflFlux[1], noBendTransmFlux[1]) ]
4850
964
	Ploss = [ 1-x-y for x,y in zip(Ptrans, Prefl)]
4851
965
		                  
4852
966
	matplotlib.pyplot.plot(frequencies, Ptrans, 'bo')
4853
967
	matplotlib.pyplot.plot(frequencies, Prefl, 'ro')
4854
968
	matplotlib.pyplot.plot(frequencies, Ploss, 'go' )
4855
969
4856
970
	matplotlib.pyplot.show()
4857
971
4858
972
4859
973
4860
974
8.3 With an inline C++ class as EPS-function
4861
975
______________________________________________
4862
976
4863
977
4864
978
The header file "eps_class.hpp" :
4865
979
4866
980
4867
981
::
4868
982
4869
983
4870
984
	using namespace meep;
4871
985
4872
986
	namespace meep
4873
987
	{
4874
988
4875
989
	class myEpsCallBack : virtual public Callback {
4876
990
4877
991
	    public:	
4878
992
		   myEpsCallBack() : Callback() { };
4879
993
		   ~myEpsCallBack() { cout << "Callback object destructed." << endl; };
4880
994
4881
995
		   myEpsCallBack(bool isWgBent,double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX)  : Callback() {
4882
996
		    	_IsWgBent = isWgBent;
4883
997
			_upperLimitHorizontalWg = upperLimitHorizontalWg;
4884
998
			_leftLimitVerticalWg = leftLimitVerticalWg;
4885
999
			_lowerLimitHorizontalWg = lowerLimitHorizontalWg;
4886
1000
			_wgLengthX = wgLengthX;
4887
1001
		   };	
4888
1002
4889
1003
		   double double_vec(const vec &v) {
4890
1004
			double eps = myEps(v, _IsWgBent, _upperLimitHorizontalWg, _leftLimitVerticalWg, _lowerLimitHorizontalWg, _wgLengthX);
4891
1005
			//cout << "X="<<v.x()<<"--Y="<<v.y()<<"--eps="<<eps<<"-"<<_upperLimitHorizontalWg<<"--"<<_leftLimitVerticalWg<<"--"<<_lowerLimitHorizontalWg<<"--"<<_wgLengthX;
4892
1006
			return eps;
4893
1007
		   };
4894
1008
		   
4895
1009
		   complex<double> complex_vec(const vec &x) { return 0; };
4896
1010
		   complex<double> complex_time(const double &t) { return 0; };   
4897
1011
4898
1012
4899
1013
	     private:
4900
1014
		   bool _IsWgBent;;
4901
1015
		   double _upperLimitHorizontalWg;
4902
1016
		   double _leftLimitVerticalWg;
4903
1017
		   double _lowerLimitHorizontalWg;
4904
1018
		   double _wgLengthX;
4905
1019
4906
1020
		   double myEps(const vec &v, bool isWgBent, double upperLimitHorizontalWg, double leftLimitVerticalWg, double lowerLimitHorizontalWg, double wgLengthX) {
4907
1021
			double xCo = v.x();
4908
1022
			double yCo = v.y();
4909
1023
			if (isWgBent){ //there is a bend
4910
1024
			    if ((yCo < upperLimitHorizontalWg) || (xCo>wgLengthX)){
4911
1025
				return 1.0;            
4912
1026
			    }
4913
1027
			    else {
4914
1028
				if ((xCo < leftLimitVerticalWg) && (yCo > lowerLimitHorizontalWg)) {
4915
1029
				     return 1.0;
4916
1030
				}
4917
1031
				else {
4918
1032
				     return 12.0;
4919
1033
				}
4920
1034
			    }
4921
1035
			}
4922
1036
			else { //there is no bend 
4923
1037
			    if ((yCo < upperLimitHorizontalWg) || (yCo > lowerLimitHorizontalWg)){
4924
1038
				return 1.0;
4925
1039
			    }
4926
1040
			}
4927
1041
			return 12.0;
4928
1042
		   }
4929
1043
4930
1044
	};
4931
1045
4932
1046
	}
4933
1047
4934
1048
4935
1049
The Python program :
4936
1050
4937
1051
4938
1052
::
4939
1053
4940
1054
4941
1055
	from meep import *  
4942
1056
4943
1057
	from math import *
4944
1058
	import numpy 
4945
1059
	import matplotlib.pyplot
4946
1060
	import sys
4947
1061
4948
1062
	from scipy.weave import *
4949
1063
4950
1064
	res = 10.0
4951
1065
	gridSizeX = 16.0
4952
1066
	gridSizeY = 32.0
4953
1067
	padSize = 4.0
4954
1068
	wgLengthX = gridSizeX - padSize
4955
1069
	wgLengthY = gridSizeY - padSize
4956
1070
	wgWidth = 1.0 #width of the waveguide
4957
1071
	upperLimitHorizontalWg = padSize
4958
1072
	lowerLimitHorizontalWg = padSize+wgWidth
4959
1073
	leftLimitVerticalWg = wgLengthX-wgWidth
4960
1074
	wgHorYCen = padSize +  wgWidth/2.0 #horizontal waveguide center Y-pos
4961
1075
	wgVerXCen = wgLengthX - wgWidth/2.0 #vertical waveguide center X-pos (in case there is a bend)
4962
1076
	srcFreqCenter = 0.15 #gaussian source center frequency
4963
1077
	srcPulseWidth = 0.1 #gaussian source pulse width
4964
1078
	srcComp = Ez #gaussian source component
4965
1079
4966
1080
4967
1081
	def initEPS():
4968
1082
		#the set of parameters that we want to pass to the Callback object upon construction
4969
1083
		c_params = ['isWgBent','upperLimitHorizontalWg','leftLimitVerticalWg','lowerLimitHorizontalWg','wgLengthX'] #all these variables must be globally declared in the scope where the "inline" function call happens
4970
1084
		#the C-code snippet for constructing the Callback object
4971
1085
		c_code = prepareCallbackCObjectCode("myEpsCallBack", c_params)
4972
1086
		#do the actual inline C-call and fetch the pointer to the Callback object
4973
1087
		funPtr = inline(c_code,c_params, libraries=getInlineLibraries(), include_dirs = getInlineInclude(), headers = getInlineHeaders("eps_class.hpp") )
4974
1088
		#set the pointer to the callback object in the Python-meep core 
4975
1089
		set_EPS_CallbackInlineObject(funPtr)
4976
1090
		print "EPS function successfully set."
4977
1091
		return
4978
1092
	 
4979
1093
	def createField(pCompVol, pWgLengthX, pWgWidth, pIsWgBent, pSrcFreqCenter, pSrcPulseWidth, pSrcComp):               
4980
1094
		#we create a structure with PML of thickness = 1.0 on all boundaries, 
4981
1095
		#in all directions,
4982
1096
		#using the material function EPS
4983
1097
		s = structure(pCompVol, EPS, pml(1.0) )     
4984
1098
		f = fields(s)      
4985
1099
		#define a gaussian line source of length 'wgWidth' at X=wgLength/2, Y=padSize
4986
1100
		srcGaussian = gaussian_src_time(pSrcFreqCenter, pSrcPulseWidth )     
4987
1101
		srcGeo = volume(vec(1.0,padSize),vec(1.0,padSize+pWgWidth))
4988
1102
		f.add_volume_source(pSrcComp, srcGaussian, srcGeo, 1)
4989
1103
		print "Field created..."
4990
1104
		return f
4991
1105
	       
4992
1106
	master_printf("BENT WAVEGUIDE SAMPLE WITH INLINE C++ CLASS FOR EPS\n")
4993
1107
4994
1108
	master_printf("Running on %d processor(s)...\n",count_processors())
4995
1109
4996
1110
	#FIRST WE WORK OUT THE CASE WITH NO BEND
4997
1111
	#----------------------------------------------------------------
4998
1112
	master_printf("*1* Starting the case with no bend...")
4999
1113
5000
1114
	#set EPS material function
Status:	Merged
Merge reported by:	Emmanuel Lambert
Merged at revision:	not available
Proposed branch:	lp:~emmanuel-lambert/python-meep/intec
Merge into:	lp:~nizamov-shawkat/python-meep/devel
Diff against target:	22461 lines (+21840/-0) 106 files modified AUTHORS (+2/-0) COPYING (+674/-0) COPYRIGHT (+17/-0) README (+27/-0) doc/html-sources/Makefile (+88/-0) doc/html-sources/conf.py (+194/-0) doc/html-sources/make.bat (+112/-0) doc/html-sources/python_meep_documentation.txt (+1299/-0) doc/html/.buildinfo (+4/-0) doc/html/_sources/.svn/all-wcprops (+11/-0) doc/html/_sources/.svn/entries (+62/-0) doc/html/_sources/.svn/format (+1/-0) doc/html/_sources/.svn/text-base/python_meep_documentation.txt.svn-base (+1293/-0) doc/html/_sources/python_meep_documentation.txt (+1299/-0) doc/html/_static/basic.css (+405/-0) doc/html/_static/default.css (+210/-0) doc/html/_static/doctools.js (+232/-0) doc/html/_static/jquery.js (+32/-0) doc/html/_static/pygments.css (+61/-0) doc/html/_static/searchtools.js (+467/-0) doc/html/genindex.html (+89/-0) doc/html/objects.inv (+3/-0) doc/html/python_meep_documentation.html (+1355/-0) doc/html/search.html (+91/-0) doc/html/searchindex.js (+1/-0) make-mpi (+10/-0) meep-common.py (+265/-0) meep-site-init.py (+24/-0) meep-site-init.py.intec (+40/-0) meep-site-init.py.standard (+24/-0) meep_common.i (+149/-0) meep_mpi.i (+27/-0) openmpi-bug.py (+17/-0) samples/bent_waveguide/eps_class.hpp (+78/-0) samples/bent_waveguide/eps_function.hpp (+55/-0) samples/bent_waveguide/go (+11/-0) samples/bent_waveguide/go-pictures (+31/-0) samples/bent_waveguide/go_inline_class (+11/-0) samples/bent_waveguide/go_inline_func (+11/-0) samples/bent_waveguide/python_meep_bent_wg.py (+203/-0) samples/bent_waveguide/python_meep_bent_wg_inline_class.py (+188/-0) samples/bent_waveguide/python_meep_bent_wg_inline_function.py (+177/-0) samples/circular_planewave_3d.py (+101/-0) samples/use_averaging.py (+41/-0) setup-mpi.py (+60/-0) setup.py (+8/-0) tests-intec/.svn/all-wcprops (+89/-0) tests-intec/.svn/entries (+511/-0) tests-intec/.svn/format (+1/-0) tests-intec/.svn/prop-base/2D_convergence.py.svn-base (+5/-0) tests-intec/.svn/prop-base/bench.py.svn-base (+5/-0) tests-intec/.svn/prop-base/bragg_transmission.py.svn-base (+5/-0) tests-intec/.svn/prop-base/convergence_cyl_waveguide.py.svn-base (+5/-0) tests-intec/.svn/prop-base/cylindrical.py.svn-base (+5/-0) tests-intec/.svn/prop-base/flux.py.svn-base (+5/-0) tests-intec/.svn/prop-base/harmonics.py.svn-base (+5/-0) tests-intec/.svn/prop-base/known_results.py.svn-base (+5/-0) tests-intec/.svn/prop-base/one_dimensional.py.svn-base (+5/-0) tests-intec/.svn/prop-base/physical.py.svn-base (+5/-0) tests-intec/.svn/prop-base/pml.py.svn-base (+5/-0) tests-intec/.svn/prop-base/symmetry.py.svn-base (+5/-0) tests-intec/.svn/prop-base/three_d.py.svn-base (+5/-0) tests-intec/.svn/prop-base/two_dimensional.py.svn-base (+5/-0) tests-intec/.svn/text-base/2D_convergence.py.svn-base (+132/-0) tests-intec/.svn/text-base/bench.py.svn-base (+280/-0) tests-intec/.svn/text-base/bragg_transmission.py.svn-base (+219/-0) tests-intec/.svn/text-base/convergence_cyl_waveguide.py.svn-base (+184/-0) tests-intec/.svn/text-base/cylindrical.py.svn-base (+307/-0) tests-intec/.svn/text-base/flux.py.svn-base (+320/-0) tests-intec/.svn/text-base/harmonics.py.svn-base (+117/-0) tests-intec/.svn/text-base/known_results.py.svn-base (+192/-0) tests-intec/.svn/text-base/one_dimensional.py.svn-base (+156/-0) tests-intec/.svn/text-base/physical.py.svn-base (+64/-0) tests-intec/.svn/text-base/pml.py.svn-base (+209/-0) tests-intec/.svn/text-base/symmetry.py.svn-base (+1059/-0) tests-intec/.svn/text-base/three_d.py.svn-base (+240/-0) tests-intec/.svn/text-base/two_dimensional.py.svn-base (+379/-0) tests-intec/2D_convergence.py (+132/-0) tests-intec/bench.py (+280/-0) tests-intec/bragg_transmission.py (+219/-0) tests-intec/convergence_cyl_waveguide.py (+183/-0) tests-intec/cylindrical.py (+307/-0) tests-intec/flux.py (+320/-0) tests-intec/harmonics.py (+117/-0) tests-intec/known_results.py (+192/-0) tests-intec/one_dimensional.py (+156/-0) tests-intec/physical.py (+64/-0) tests-intec/pml.py (+209/-0) tests-intec/symmetry.py (+1059/-0) tests-intec/three_d.py (+240/-0) tests-intec/two_dimensional.py (+379/-0) tests/2D_convergence.py (+134/-0) tests/bench.py (+282/-0) tests/bragg_transmission.py (+221/-0) tests/convergence_cyl_waveguide.py (+185/-0) tests/cylindrical.py (+309/-0) tests/flux.py (+321/-0) tests/harmonics.py (+119/-0) tests/known_results.py (+194/-0) tests/one_dimensional.py (+158/-0) tests/physical.py (+66/-0) tests/pml.py (+211/-0) tests/symmetry.py (+1059/-0) tests/three_d.py (+242/-0) tests/two_dimensional.py (+381/-0) weave-bug.py (+12/-0)
To merge this branch:	bzr merge lp:~emmanuel-lambert/python-meep/intec
Related bugs:	Link a bug report
Related blueprints:	callback functions with inline C or C++ (Undefined)
Reviewer	Review Type	Date Requested	Status
Nizamov Shawkat		2009-08-27	Approve on 2009-08-27
Review via email: mp+10791@code.launchpad.net