python-meep

Merge lp:~emmanuel-lambert/python-meep/intec into lp:~nizamov-shawkat/python-meep/devel

  1. intec
  2. Merge into devel
Proposed by Emmanuel Lambert
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
Reviewer Review Type Date Requested Status
Nizamov Shawkat Approve
Review via email: mp+10791@code.launchpad.net
To post a comment you must log in.
Revision history for this message
Emmanuel Lambert (emmanuel-lambert) wrote :

proposal for merge into Shawkat's branch

Revision history for this message
Nizamov Shawkat (nizamov-shawkat) :
review: Approve
lp:~emmanuel-lambert/python-meep/intec updated
8. By Emmanuel Lambert <email address hidden>

*Introduction of new class 'structure_eps_pml' that always defautls to EPS-averaging enabled (consistent with Scheme behaviour)
*Changes to the bent_waveguide sample to make it more consistent with the Scheme sample
*Conversion to meep-1.1.1

9. By Emmanuel Lambert <email address hidden>

file meep_common.i

10. By Emmanuel Lambert <email address hidden>

Declare EPS-functions through inline C-function or C++-class (improved performance).
Update of sources, sample and documentation.

11. By Emmanuel Lambert <email address hidden>

Examples for EPS-function through inline C or C++

12. By Emmanuel Lambert <email address hidden>

update of setup.py script

13. By Emmanuel Lambert <email address hidden>

changes in make and setup.py to support -I parameter

14. By Emmanuel Lambert <email address hidden>

deletion of .svn files

15. By Emmanuel Lambert <email address hidden>

changes in configuration files

16. By Emmanuel Lambert <email address hidden>

configuration files for MPI

17. By Emmanuel Lambert <email address hidden>

bugfix for MPI version

18. By Emmanuel Lambert <email address hidden>

various improvements; resolution of bug 447309

19. By Emmanuel Lambert <email address hidden>

sources of doc files

20. By Emmanuel Lambert <email address hidden>

sources of documentation files

21. By Emmanuel Lambert <email address hidden>

* callback function for source amplitude factor
* merge of NS code for custom source
* update of documentation, fix issues 457159+457156

22. By Emmanuel Lambert <email address hidden>

merging differences with NS branch - not yet completely functional and tested

23. By Emmanuel Lambert <email address hidden>

forgot to add meep-common.py

24. By Emmanuel Lambert <email address hidden>

removed a redundant file

25. By Emmanuel Lambert <email address hidden>

added 3 warnings at startup / complex double callback now working fine / implemented site-specific initialisations

26. By Emmanuel Lambert <email address hidden>

intermediate commit with update of unit tests (remaining: cylindrical, symmetry, three)

27. By Emmanuel Lambert <email address hidden>

v0.7beta - **all unit tests pass** -> MPI-version remains to be tested

28. By Emmanuel Lambert <email address hidden>

updates to setup script, more specifically for MPI

29. By Emmanuel Lambert <email address hidden>

minor update to bent waveguide sample

30. By Emmanuel Lambert <email address hidden>

further merging with NS / replacing tests dir with default tests and creating new 'tests-intec' subdir

31. By Emmanuel Lambert <email address hidden>

further merging with NS branch / update of tests subdir / added test-intec subdir

32. By Emmanuel Lambert <email address hidden>

reverted structure class to be in line with c++ core / imported extra circular_planewave example / fixed complex_time callback issue

33. By Emmanuel Lambert <email address hidden>

paths for incline C were hardcoded / replaced by environment variables

34. By Emmanuel Lambert <email address hidden>

small error in eps_class.hpp (inline C example)

35. By Emmanuel Lambert <email address hidden>

update of documentation - consistency with version 0.8

36. By Emmanuel Lambert <email address hidden>

minor changes in merging with NS branch

37. By Emmanuel Lambert <email address hidden>

minor change for merge

38. By Emmanuel Lambert <email address hidden>

added authors & GPL stuff

39. By Emmanuel Lambert <email address hidden>

merging

40. By Emmanuel Lambert <email address hidden>

version 1.0

Preview Diff

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

Subscribers

People subscribed via source and target branches