/
Python_User_Interface.md.in
970 lines (613 loc) · 54.8 KB
/
Python_User_Interface.md.in
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
---
# Python User Interface
---
This page is a listing of the functions exposed by the Python interface. For a gentler introduction, see [Tutorial/Basics](Python_Tutorials/Basics.md). Note that this page is not a complete listing of all functions. In particular, because of the [SWIG wrappers](#swig-wrappers), every function in the C++ interface is accessible from the Python module, but not all of these functions are documented or intended for end users. See also the instructions for [Parallel Meep](Parallel_Meep.md).
The Python API functions and classes can be found in the `meep` module, which should be installed in your Python system by Meep's `make install` script. If you installed into a nonstandard location (e.g. your home directory), you may need to set the `PYTHONPATH` environment variable as documented in [Building From Source](Build_From_Source.md#building-from-source). You typically import the `meep` module in Python via `import meep as mp`.
[TOC]
Predefined Variables
--------------------
These are available directly via the `meep` package.
**`air`, `vacuum` [`Medium` class ]**
—
Two aliases for a predefined material type with a dielectric constant of 1.
**`perfect_electric_conductor` or `metal` [`Medium` class ]**
—
A predefined material type corresponding to a perfect electric conductor at the boundary of which the parallel electric field is zero. Technically, $\varepsilon = -\infty$.
**`perfect_magnetic_conductor` [`Medium` class ]**
—
A predefined material type corresponding to a perfect magnetic conductor at the boundary of which the parallel magnetic field is zero. Technically, $\mu = -\infty$.
**`inf` [`number`]**
—
A big number (10<sup>20</sup>) to use for "infinite" dimensions of objects.
Constants (Enumerated Types)
----------------------------
Several of the functions/classes in Meep ask you to specify e.g. a field component or a direction in the grid. These should be one of the following constants (which are available directly via the `meep` package):
**`direction` constants**
—
Specify a direction in the grid. One of `X`, `Y`, `Z`, `R`, `P` for $x$, $y$, $z$, $r$, $\phi$, respectively.
**`side` constants**
—
Specify particular boundary in the positive `High` (e.g., +`X`) or negative `Low` (e.g., -`X`) direction.
**`boundary_condition` constants**
—
`Metallic` (i.e., zero electric field) or `Magnetic` (i.e., zero magnetic field).
**`component` constants**
—
Specify a particular field or other component. One of `Ex`, `Ey`, `Ez`, `Er`, `Ep`, `Hx`, `Hy`, `Hz`, `Hy`, `Hp`, `Hz`, `Bx`, `By`, `Bz`, `By`, `Bp`, `Bz`, `Dx`, `Dy`, `Dz`, `Dr`, `Dp`, `Dielectric`, `Permeability`, for $E_x$, $E_y$, $E_z$, $E_r$, $E_\phi$, $H_x$, $H_y$, $H_z$, $H_r$, $H_\phi$, $B_x$, $B_y$, $B_z$, $B_r$, $B_\phi$, $D_x$, $D_y$, $D_z$, $D_r$, $D_\phi$, ε, μ, respectively.
**`derived_component` constants**
—
These are additional components which are not actually stored by Meep but are computed as needed, mainly for use in output functions. One of `Sx`, `Sy`, `Sz`, `Sr`, `Sp`, `EnergyDensity`, `D_EnergyDensity`, `H_EnergyDensity` for $S_x$, $S_y$, $S_z$, $S_r$, $S_\phi$ (components of the Poynting vector $\mathrm{Re}\,\mathbf{E}^* \times \mathbf{H}$), $(\mathbf{E}^* \cdot \mathbf{D} + \mathbf{H}^* \cdot \mathbf{B})/2$, $\mathbf{E}^* \cdot \mathbf{D}/2$, $\mathbf{H}^* \cdot \mathbf{B}/2$, respectively.
The Simulation Class
---------------------
The `Simulation` class is the primary abstraction of the high-level interface. Minimally, a simulation script amounts to passing the desired keyword arguments to the `Simulation` constructor and calling the `run` method on the resulting instance.
@@ Simulation @@
@@ Simulation.__init__ @@
@@ Simulation.run @@
### Output File Names
The output filenames used by Meep, e.g. for HDF5 files, are automatically prefixed by the
input variable `filename_prefix`. If `filename_prefix` is `None` (the default), however,
then Meep constructs a default prefix based on the current Python filename with `".py"`
replaced by `"-"`: e.g. `test.py` implies a prefix of `"test-"`. You can get this prefix,
or set the output folder, with these methods of the `Simulation` class:
@@ Simulation.get_filename_prefix @@
@@ Simulation.use_output_directory @@
### Simulation Time
The `Simulation` class provides the following time-related methods:
@@ Simulation.meep_time @@
@@ Simulation.print_times @@
@@ Simulation.time_spent_on @@
@@ Simulation.mean_time_spent_on @@
@@ Simulation.output_times @@
### Field Computations
Meep supports a large number of functions to perform computations on the fields. Most of them are accessed via the lower-level C++/SWIG interface. Some of them are based on the following simpler, higher-level versions. They are accessible as methods of a `Simulation` instance.
@@ Simulation.set_boundary @@
@@ Simulation.phase_in_material @@
@@ Simulation.get_field_point @@
@@ Simulation.get_epsilon_point @@
@@ Simulation.get_mu_point @@
@@ Simulation.get_epsilon_grid @@
@@ Simulation.initialize_field @@
@@ Simulation.add_dft_fields @@
@@ Simulation.flux_in_box @@
@@ Simulation.electric_energy_in_box @@
@@ Simulation.magnetic_energy_in_box @@
@@ Simulation.field_energy_in_box @@
@@ Simulation.modal_volume_in_box @@
@@ Simulation.integrate_field_function @@
@@ Simulation.max_abs_field_function @@
@@ Simulation.integrate2_field_function @@
### Reloading Parameters
Once the fields/simulation have been initialized, you can change the values of various parameters by using the following functions (which are members of the `Simulation` class):
@@ Simulation.reset_meep @@
@@ Simulation.restart_fields @@
@@ Simulation.change_k_point @@
@@ Simulation.change_sources @@
@@ Simulation.set_materials @@
### Flux Spectra
Given a bunch of [`FluxRegion`](#fluxregion) objects, you can tell Meep to accumulate the Fourier transforms of the fields in those regions in order to compute the Poynting flux spectra. (Note: as a matter of convention, the "intensity" of the electromagnetic fields refers to the Poynting flux, *not* to the [energy density](#energy-density-spectra).) See also [Introduction/Transmittance/Reflectance Spectra](Introduction.md#transmittancereflectance-spectra) and [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend). These are attributes of the `Simulation` class. The most important function is:
@@ Simulation.add_flux @@
As described in the tutorial, you normally use `add_flux` via statements like:
```python
transmission = sim.add_flux(...)
```
to store the flux object in a variable. You can create as many flux objects as you want, e.g. to look at powers flowing in different regions or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the flux region multiplied by the number of electric and magnetic field components required to get the Poynting vector multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
Once you have called `add_flux`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the [run functions](#run-functions). At any time, you can ask for Meep to print out the current flux spectrum via the `display_fluxes` method.
@@ Simulation.display_fluxes @@
You might have to do something lower-level if you have multiple flux regions corresponding to *different* frequency ranges, or have other special needs. `display_fluxes(f1, f2, f3)` is actually equivalent to `meep.display_csv("flux", meep.get_flux_freqs(f1), meep.get_fluxes(f1), meep.get_fluxes(f2), meep.get_fluxes(f3))`, where `display_csv` takes a bunch of lists of numbers and prints them as a comma-separated table; this involves calling two lower-level functions:
@@ get_flux_freqs @@
@@ get_fluxes @@
As described in [Introduction/Transmittance/Reflectance Spectra](Introduction.md#transmittancereflectance-spectra) and [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend), for a reflection spectrum you often want to save the Fourier-transformed fields from a "normalization" run and then load them into another run to be subtracted. This can be done via:
@@ Simulation.save_flux @@
@@ Simulation.load_flux @@
@@ Simulation.load_minus_flux @@
Sometimes it is more convenient to keep the Fourier-transformed fields in memory rather than writing them to a file and immediately loading them back again. To that end, the `Simulation` class exposes the following three methods:
@@ Simulation.get_flux_data @@
@@ Simulation.load_flux_data @@
@@ Simulation.load_minus_flux_data @@
The `Simulation` class also provides some aliases for the corresponding "flux" methods.
* **`save_mode`**
* **`load_mode`**
* **`load_minus_mode`**
* **`get_mode_data`**
* **`load_mode_data`**
* **`load_minus_mode_data`**
### Mode Decomposition
Given a structure, Meep can decompose the Fourier-transformed fields into a superposition of its harmonic modes. For a theoretical background, see [Features/Mode Decomposition](Mode_Decomposition.md).
@@ Simulation.get_eigenmode_coefficients @@
The flux object should be created using `add_mode_monitor`. (You could also use `add_flux`, but with `add_flux` you need to be more careful about symmetries that bisect the flux plane: the `add_flux` object should only be used with `get_eigenmode_coefficients` for modes of the same symmetry, e.g. constrained via `eig_parity`. On the other hand, the performance of `add_flux` planes benefits more from symmetry.) `eig_vol` is the volume passed to [MPB](https://mpb.readthedocs.io) for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequency-domain fields are tabulated, which is the default (i.e. `flux.where`). `eig_parity` should be one of [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has $z$ and/or $y$ mirror symmetry *in the source region*, just as for `EigenModeSource` above. If the structure has both $y$ and $z$ mirror symmetry, you can combine more than one of these, e.g. `EVEN_Z+ODD_Y`. Default is `NO_PARITY`, in which case MPB computes all of the bands which will still be even or odd if the structure has mirror symmetry, of course. This is especially useful in 2d simulations to restrict yourself to a desired polarization. `eig_resolution` is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep `resolution` in which case the structure is linearly interpolated from the Meep pixels. `eig_tolerance` is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to `1e-12`. (Note that this is the tolerance for the frequency eigenvalue $\omega$; the tolerance for the mode profile is effectively the square root of this.) For examples, see [Tutorial/Mode Decomposition](Python_Tutorials/Mode_Decomposition.md).
Technically, MPB computes $\omega_n(\mathbf{k})$ and then inverts it with Newton's method to find the wavevector $\mathbf{k}$ normal to `eig_vol` and mode for a given frequency; in rare cases (primarily waveguides with *nonmonotonic* dispersion relations, which doesn't usually happen in simple dielectric waveguides), MPB may need you to supply an initial "guess" for $\mathbf{k}$ in order for this Newton iteration to converge. You can supply this initial guess with `kpoint_func`, which is a function `kpoint_func(f, n)` that supplies a rough initial guess for the $\mathbf{k}$ of band number $n$ at frequency $f=\omega/2\pi$. (By default, the $\mathbf{k}$ components in the plane of the `eig_vol` region are zero. However, if this region spans the *entire* cell in some directions, and the cell has Bloch-periodic boundary conditions via the `k_point` parameter, then the mode's $\mathbf{k}$ components in those directions will match `k_point` so that the mode satisfies the Meep boundary conditions, regardless of `kpoint_func`.) If `direction` is set to `mp.NO_DIRECTION`, then `kpoint_func` is not only the initial guess and the search direction of the $\mathbf{k}$ vectors, but is also taken to be the direction of the waveguide, allowing you to [detect modes in oblique waveguides](Python_Tutorials/Eigenmode_Source.md#oblique-waveguides) (not perpendicular to the flux plane).
**Note:** for planewaves in homogeneous media, the `kpoints` may *not* necessarily be equivalent to the actual wavevector of the mode. This quantity is given by `kdom`.
Note that Meep's MPB interface only supports dispersionless non-magnetic materials but it does support anisotropic $\varepsilon$. Any nonlinearities, magnetic responses $\mu$ conductivities $\sigma$, or dispersive polarizations in your materials will be *ignored* when computing the mode decomposition. PML will also be ignored.
@@ Simulation.add_mode_monitor @@
`add_mode_monitor` works properly with arbitrary symmetries, but may be suboptimal because the Fourier-transformed region does not exploit the symmetry. As an optimization, if you have a mirror plane that bisects the mode monitor, you can instead use `add_flux` to gain a factor of two, but in that case you *must* also pass the corresponding `eig_parity` to `get_eigenmode_coefficients` in order to only compute eigenmodes with the corresponding mirror symmetry.
@@ Simulation.get_eigenmode @@
The following top-level function is also available:
@@ get_eigenmode_freqs @@
@@ DiffractedPlanewave[methods-with-docstrings] @@
### Energy Density Spectra
Very similar to flux spectra, you can also compute **energy density spectra**: the energy density of the electromagnetic fields as a function of frequency, computed by Fourier transforming the fields and integrating the energy density:
$$ \frac{1}{2}ε|\mathbf{E}|^2 + \frac{1}{2}μ|\mathbf{H}|^2 $$
The usage is similar to the flux spectra: you define a set of [`EnergyRegion`](#EnergyRegion) objects telling Meep where it should compute the Fourier-transformed fields and energy densities, and call `add_energy` to add these regions to the current simulation over a specified frequency bandwidth, and then use `display_electric_energy`, `display_magnetic_energy`, or `display_total_energy` to display the energy density spectra at the end. There are also `save_energy`, `load_energy`, and `load_minus_energy` functions that you can use to subtract the fields from two simulation, e.g. in order to compute just the energy from scattered fields, similar to the flux spectra. The function used to add an [`EnergyRegion`](#EnergyRegion) is as follows:
@@ Simulation.add_energy @@
As for flux regions, you normally use `add_energy` via statements like:
```py
En = sim.add_energy(...)
```
to store the energy object in a variable. You can create as many energy objects as you want, e.g. to look at the energy densities in different objects or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the energy region multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
Once you have called `add_energy`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the `run` functions. At any time, you can ask for Meep to print out the current energy density spectrum via:
@@ Simulation.display_electric_energy @@
@@ Simulation.display_magnetic_energy @@
@@ Simulation.display_total_energy @@
You might have to do something lower-level if you have multiple energy regions corresponding to *different* frequency ranges, or have other special needs. `display_electric_energy(e1, e2, e3)` is actually equivalent to `meep.display_csv("electric_energy", meep.get_energy_freqs(e1), meep.get_electric_energy(e1), meep.get_electric_energy(e2), meep.get_electric_energy(e3))`, where `display_csv` takes a bunch of lists of numbers and prints them as a comma-separated table; this involves calling lower-level functions:
@@ get_energy_freqs @@
@@ get_electric_energy @@
@@ get_magnetic_energy @@
@@ get_total_energy @@
As described in [Introduction/Transmittance/Reflectance Spectra](Introduction.md#transmittancereflectance-spectra) and [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend) for flux computations, to compute the energy density from the scattered fields you often want to save the Fourier-transformed fields from a "normalization" run and then load them into another run to be subtracted. This can be done via:
@@ Simulation.save_energy @@
@@ Simulation.load_energy @@
@@ Simulation.load_minus_energy @@
### Force Spectra
Very similar to flux spectra, you can also compute **force spectra**: forces on an object as a function of frequency, computed by Fourier transforming the fields and integrating the vacuum [Maxwell stress tensor](https://en.wikipedia.org/wiki/Maxwell_stress_tensor):
$$\sigma_{ij} = E_i^*E_j + H_i^*H_j - \frac{1}{2} δ_{ij} \left( |\mathbf{E}|^2 + |\mathbf{H}|^2 \right)$$
over a surface $S$ via $\mathbf{F} = \int_S \sigma d\mathbf{A}$. You should normally **only evaluate the stress tensor over a surface lying in vacuum**, as the interpretation and definition of the stress tensor in arbitrary media is often problematic (the subject of extensive and controversial literature). It is fine if the surface *encloses* an object made of arbitrary materials, as long as the surface itself is in vacuum.
See also [Tutorial/Optical Forces](Python_Tutorials/Optical_Forces.md).
Most commonly, you will want to **normalize** the force spectrum in some way, just as for flux spectra. Most simply, you could divide two different force spectra to compute the ratio of forces on two objects. Often, you will divide a force spectrum by a flux spectrum, to divide the force $F$ by the incident power $P$ on an object, in order to compute the useful dimensionless ratio $Fc$/$P$ where $c=1$ in Meep units. For example, it is a simple exercise to show that the force $F$ on a perfectly reflecting mirror with normal-incident power $P$ satisfies $Fc$/$P=2$, and for a perfectly absorbing (black) surface $Fc$/$P=1$.
The usage is similar to the [flux spectra](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend): you define a set of [`ForceRegion`](#ForceRegion) objects telling Meep where it should compute the Fourier-transformed fields and stress tensors, and call `add_force` to add these regions to the current simulation over a specified frequency bandwidth, and then use `display_forces` to display the force spectra at the end. There are also `save_force`, `load_force`, and `load_minus_force` functions that you can use to subtract the fields from two simulation, e.g. in order to compute just the force from scattered fields, similar to the flux spectra. The function used to add a [`ForceRegion`](#ForceRegion) object is defined as follows:
@@ Simulation.add_force @@
As for flux regions, you normally use `add_force` via statements like:
```py
Fx = sim.add_force(...)
```
to store the force object in a variable. You can create as many force objects as you want, e.g. to look at forces on different objects, in different directions, or in different frequency ranges. Note, however, that Meep has to store (and update at every time step) a number of Fourier components equal to the number of grid points intersecting the force region, multiplied by the number of electric and magnetic field components required to get the stress vector, multiplied by `nfreq`, so this can get quite expensive (in both memory and time) if you want a lot of frequency points over large regions of space.
Once you have called `add_force`, the Fourier transforms of the fields are accumulated automatically during time-stepping by the `run` functions. At any time, you can ask for Meep to print out the current force spectrum via:
@@ Simulation.display_forces @@
You might have to do something lower-level if you have multiple force regions corresponding to *different* frequency ranges, or have other special needs. `display_forces(f1, f2, f3)` is actually equivalent to `meep.display_csv("force", meep.get_force_freqs(f1), meep.get_forces(f1), meep.get_forces(f2), meep.get_forces(f3))`, where `display_csv` takes a bunch of lists of numbers and prints them as a comma-separated table; this involves calling two lower-level functions:
@@ get_force_freqs @@
@@ get_forces @@
As described in [Introduction/Transmittance/Reflectance Spectra](Introduction.md#transmittancereflectance-spectra) and [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend) for flux computations, to compute the force from the scattered fields often requires saving the Fourier-transformed fields from a "normalization" run and then loading them into another run to be subtracted. This can be done via these `Simulation` methods:
@@ Simulation.save_force @@
@@ Simulation.load_force @@
@@ Simulation.load_minus_force @@
To keep the fields in memory and avoid writing to and reading from a file, use the following three `Simulation` methods:
@@ Simulation.get_force_data @@
@@ Simulation.load_force_data @@
@@ Simulation.load_minus_force_data @@
### LDOS spectra
Meep can also calculate the LDOS (local density of states) spectrum, as described in [Tutorial/Local Density of States](Python_Tutorials/Local_Density_of_States.md). To do this, you simply pass the following step function to your `run` command:
@@ Ldos @@
@@ get_ldos_freqs @@
@@ dft_ldos @@
Analytically, the per-polarization LDOS is exactly proportional to the power radiated by an $\ell$-oriented point-dipole current, $p(t)$, at a given position in space. For a more mathematical treatment of the theory behind the LDOS, refer to the relevant discussion in Section 4.4 ("Currents and Fields: The Local Density of States") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of the book [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707), but for now it is defined as:
$$\operatorname{LDOS}_{\ell}(\vec{x}_0,\omega)=-\frac{2}{\pi}\varepsilon(\vec{x}_0)\frac{\operatorname{Re}[\hat{E}_{\ell}(\vec{x}_0,\omega)\hat{p}(\omega)^*]}{|\hat{p}(\omega)|^2}$$
where the $|\hat{p}(\omega)|^2$ normalization is necessary for obtaining the power exerted by a unit-amplitude dipole (assuming linear materials), and hats denote Fourier transforms. It is this quantity that is computed by the `dft_ldos` command for a single dipole source. For a volumetric source, the numerator and denominator are both integrated over the current volume, but "LDOS" computation is less meaningful in this case.
### Near-to-Far-Field Spectra
Meep can compute a near-to-far-field transformation in the frequency domain as described in [Tutorial/Near-to-Far Field Spectra](Python_Tutorials/Near_to_Far_Field_Spectra.md): given the fields on a "near" bounding surface inside the cell, it can compute the fields arbitrarily far away using an analytical transformation, assuming that the "near" surface and the "far" region lie in a single homogeneous non-periodic 2d, 3d, or cylindrical region. That is, in a simulation *surrounded by PML* that absorbs outgoing waves, the near-to-far-field feature can compute the fields outside the cell as if the outgoing waves had not been absorbed (i.e. in the fictitious infinite open volume). Moreover, this operation is performed on the Fourier-transformed fields: like the flux and force spectra above, you specify a set of desired frequencies, Meep accumulates the Fourier transforms, and then Meep computes the fields at *each frequency* for the desired far-field points.
This is based on the principle of equivalence: given the Fourier-transformed tangential fields on the "near" surface, Meep computes equivalent currents and convolves them with the analytical Green's functions in order to compute the fields at any desired point in the "far" region. For details, see Section 4.2.1 ("The Principle of Equivalence") in [Chapter 4](http://arxiv.org/abs/arXiv:1301.5366) ("Electromagnetic Wave Source Conditions") of the book [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707). Since the "far" fields are computed using the full Green's functions, they should be able to be computed *anywhere* outside of the near-field surface monitor. The only limiting factor should be discretization errors but for any given distannce, the "far" fields should converge to the actual DFT fields at that location with resolution (assuming the distance separation is >> resolution).
There are three steps to using the near-to-far-field feature: first, define the "near" surface(s) as a set of surfaces capturing *all* outgoing radiation in the desired direction(s); second, run the simulation, typically with a pulsed source, to allow Meep to accumulate the Fourier transforms on the near surface(s); third, tell Meep to compute the far fields at any desired points (optionally saving the far fields from a grid of points to an HDF5 file). To define the near surfaces, use this `Simulation` method:
@@ Simulation.add_near2far @@
Each `Near2FarRegion` is identical to `FluxRegion` except for the name: in 3d, these give a set of planes (**important:** all these "near surfaces" must lie in a single *homogeneous* material with *isotropic* ε and μ — and they should *not* lie in the PML regions) surrounding the source(s) of outgoing radiation that you want to capture and convert to a far field. Ideally, these should form a closed surface, but in practice it is sufficient for the `Near2FarRegion`s to capture all of the radiation in the direction of the far-field points. **Important:** as for flux computations, each `Near2FarRegion` should be assigned a `weight` of ±1 indicating the direction of the outward normal relative to the +coordinate direction. So, for example, if you have six regions defining the six faces of a cube, i.e. the faces in the $+x$, $-x$, $+y$, $-y$, $+z$, and $-z$ directions, then they should have weights +1, -1, +1, -1, +1, and -1 respectively. Note that, neglecting discretization errors, all near-field surfaces that enclose the same outgoing fields are equivalent and will yield the same far fields with a discretization-induced difference that vanishes with increasing resolution etc.
After the simulation run is complete, you can compute the far fields. This is usually for a pulsed source so that the fields have decayed away and the Fourier transforms have finished accumulating.
If you have Bloch-periodic boundary conditions, then the corresponding near-to-far transformation actually needs to perform a "lattice sum" of infinitely many periodic copies of the near fields. This doesn't happen by default, which means the default `near2far` calculation may not be what you want for periodic boundary conditions. However, if the `Near2FarRegion` spans the entire cell along the periodic directions, you can turn on an approximate lattice sum by passing `nperiods > 1`. In particular, it then sums `2*nperiods+1` Bloch-periodic copies of the near fields whenever a far field is requested. You can repeatedly double `nperiods` until the answer converges to your satisfaction; in general, if the far field is at a distance $d$, and the period is $a$, then you want `nperiods` to be much larger than $d/a$. (Future versions of Meep may use fancier techniques like [Ewald summation](https://en.wikipedia.org/wiki/Ewald_summation) to compute the lattice sum more rapidly at large distances.)
@@ Simulation.get_farfield @@
@@ Simulation.output_farfields @@
@@ Simulation.get_farfields @@
This lower-level function is also available:
@@ get_near2far_freqs @@
(Multi-frequency `get_farfields` and `output_farfields` can be accelerated by
[compiling Meep](Build_From_Source.md#meep) with `--with-openmp` and using the
`OMP_NUM_THREADS` environment variable to specify multiple threads.)
For a scattered-field computation, you often want to separate the scattered and incident fields. As described in [Introduction/Transmittance/Reflectance Spectra](Introduction.md#transmittancereflectance-spectra) and [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](Scheme_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend) for flux computations, you can do this by saving the Fourier-transformed incident from a "normalization" run and then load them into another run to be subtracted. This can be done via these `Simulation` methods:
@@ Simulation.save_near2far @@
@@ Simulation.load_near2far @@
@@ Simulation.load_minus_near2far @@
To keep the fields in memory and avoid writing to and reading from a file, use the following three methods:
@@ Simulation.get_near2far_data @@
@@ Simulation.load_near2far_data @@
@@ Simulation.load_minus_near2far_data @@
See also this lower-level function:
@@ scale_near2far_fields @@
And this [`DftNear2Far`](#DftNear2Far) method:
@@ DftNear2Far.flux @@
### Load and Dump Simulation State
These functions add support for saving and restoring parts of the
`Simulation` state.
For all functions listed below, when dumping/loading state to/from a distributed filesystem
(using say, parallel HDF5) and running in an MPI environment, setting
`single_parallel_file=True` (the default) will result in all
processes writing/reading to/from the same/single file after computing their respective offsets into this file.
When set to `False`, each process writes/reads data for the chunks it owns
to/from a separate, process unique file.
#### Load and Dump Structure
These functions dump the raw ε and μ data to disk and load it back for doing multiple simulations with the same materials but different sources etc. The only prerequisite is that the dump/load simulations have the same [chunks](Chunks_and_Symmetry.md) (i.e. the same grid, number of processors, symmetries, and PML). When using `split_chunks_evenly=False`, you must also dump the original chunk layout using `dump_chunk_layout` and load it into the new `Simulation` using the `chunk_layout` parameter. Currently only stores dispersive and non-dispersive $\varepsilon$ and $\mu$ but not nonlinearities. Note that loading data from a file in this way overwrites any `geometry` data passed to the `Simulation` constructor.
@@ Simulation.dump_structure @@
@@ Simulation.load_structure @@
#### Load and Dump Chunk Layout
@@ Simulation.dump_chunk_layout @@
To load a chunk layout into a `Simulation`, use the `chunk_layout` argument to the constructor, passing either a file obtained from `dump_chunk_layout` or another `Simulation` instance. Note that when using `split_chunks_evenly=False` this parameter is required when saving and loading flux spectra, force spectra, or near-to-far spectra so that the two runs have the same chunk layout. Just pass the `Simulation` object from the first run to the second run:
```python
# Split chunks based on amount of work instead of size
sim1 = mp.Simulation(..., split_chunks_evenly=False)
norm_flux = sim1.add_flux(...)
sim1.run(...)
sim1.save_flux(...)
# Make sure the second run uses the same chunk layout as the first
sim2 = mp.Simulation(..., chunk_layout=sim1)
flux = sim2.add_flux(...)
sim2.load_minus_flux(...)
sim2.run(...)
```
#### Load and Dump Fields
These functions can be used to dump (and later load) the time-domain fields, auxiliary
fields for PMLs, polarization fields (for dispersive materials), and the DFT fields
at a certain timestamp. The timestamp at which the dump happens is also saved so that
the simulation can continue from where it was saved. The one pre-requisite of this
feature is that it needs the `Simulation` object to have been setup *exactly* the
same as the one it was dumped from.
@@ Simulation.dump_fields @@
@@ Simulation.load_fields @@
#### Checkpoint and Restore
The above dump/load related functions can be used to implement a
checkpoint/restore *like* feature. The caveat of this feature is that
it does *not* store all the state required to re-create the `Simulation` object
itself. Instead, it expects the user to create and set up the new `Simulation`
object to be *exactly* the same as the one the state was dumped from.
@@ Simulation.dump @@
@@ Simulation.load @@
Example usage:
```python
# Setup, run and dump the simulation.
sim1 = mp.Simulation(...)
sim1.run(..., until=xx)
sim1.dump(dirname, dump_structure=True, dump_fields=True, ...)
...
# Later in the same or a new process/run
sim2 = mp.Simulation(<same setup as sim1>)
sim2.load(dirname, load_structure=True, load_fields=True, ...)
sim2.run(...) # continues from t=xx
```
### Frequency-Domain Solver
Meep contains a frequency-domain solver that computes the fields produced in a geometry in response to a [continuous-wave (CW) source](https://en.wikipedia.org/wiki/Continuous_wave). This is based on an [iterative linear solver](https://en.wikipedia.org/wiki/Iterative_method) instead of time-stepping. For details, see Section 5.3 ("Frequency-domain solver") of [Computer Physics Communications, Vol. 181, pp. 687-702 (2010)](http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf). Benchmarking results have shown that in many instances, such as cavities (e.g., [ring resonators](Python_Tutorials/Frequency_Domain_Solver.md)) with long-lived resonant modes, this solver converges much faster than simply running an equivalent time-domain simulation with a CW source (using the default `width` of zero for no transient turn-on), time-stepping until all transient effects from the source turn-on have disappeared, especially if the fields are desired to a very high accuracy.
To use the frequency-domain solver, simply define a `ContinuousSource` with the desired frequency and [initialize the fields and geometry](#initializing-the-structure-and-fields) via `init_sim()`:
```py
sim = mp.Simulation(...)
sim.init_sim()
sim.solve_cw(tol, maxiters, L)
```
The first two parameters to the frequency-domain solver are the tolerance `tol` for the iterative solver (10<sup>−8</sup>, by default) and a maximum number of iterations `maxiters` (10<sup>4</sup>, by default). Finally, there is a parameter $L$ that determines a tradeoff between memory and work per step and convergence rate of the iterative algorithm, biconjugate gradient stabilized ([BiCGSTAB-L](https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method)), that is used; larger values of $L$ will often lead to faster convergence at the expense of more memory and more work per iteration. Default is $L=2$, and normally a value ≥ 2 should be used.
The frequency-domain solver supports arbitrary geometries, PML, boundary conditions, symmetries, parallelism, conductors, and arbitrary nondispersive materials. Lorentz-Drude dispersive materials are not currently supported in the frequency-domain solver, but since you are solving at a known fixed frequency rather than timestepping, you should be able to pick conductivities etcetera in order to obtain any desired complex ε and μ at that frequency.
The frequency-domain solver requires you to use complex-valued fields, via `force_complex_fields=True`.
After `solve_cw` completes, it should be as if you had just run the simulation for an infinite time with the source at that frequency. You can call the various field-output functions and so on as usual at this point. For examples, see [Tutorial/Frequency Domain Solver](Python_Tutorials/Frequency_Domain_Solver.md) and [Tutorial/Mode Decomposition/Reflectance and Transmittance Spectra for Planewave at Oblique Incidence](Python_Tutorials/Mode_Decomposition.md#reflectance-and-transmittance-spectra-for-planewave-at-oblique-incidence).
**Note:** The convergence of the iterative solver can sometimes encounter difficulties. For example, increasing the diameter of a ring resonator relative to the wavelength increases the [condition number](https://en.wikipedia.org/wiki/Condition_number), which worsens the convergence of iterative solvers. The general way to improve this is to implement a more sophisticated iterative solver that employs [preconditioners](https://en.wikipedia.org/wiki/Preconditioner). Preconditioning wave equations (Helmholtz-like equations) is notoriously difficult to do well, but some possible strategies are discussed in [Issue #548](https://github.com/NanoComp/meep/issues/548). In the meantime, a simpler way improving convergence (at the expense of computational cost) is to increase the $L$ parameter and the number of iterations.
### Frequency-Domain Eigensolver
Building on the frequency-domain solver above, Meep also includes a frequency-domain eigensolver that computes resonant frequencies and modes in the frequency domain. The usage is very similar to `solve_cw`:
```py
sim = mp.Simulation(...)
sim.init_sim()
eigfreq = sim.solve_eigfreq(tol, maxiters, guessfreq, cwtol, cwmaxiters, L)
```
The `solve_eig` routine performs repeated calls to `solve_cw` in a way that converges to the resonant mode whose frequency is *closest* to the source frequency. The complex resonant-mode frequency is returned, and the mode Q can be computed from `eigfreq.real / (-2*eigfreq.imag)`. Upon return, the fields should be the corresponding resonant mode (with an arbitrary scaling).
The resonant mode is converged to a relative error of roughly `tol`, which defaults to `1e-7`. A maximum of `maxiters` (defaults to `100`) calls to `solve_cw` are performed. The tolerance for each `solve_cw` call is `cwtol` (defaults to `tol*1e-3`) and the maximum iterations is `cwmaxiters` (10<sup>4</sup>, by default); the `L` parameter (defaults to `10`) is also passed through to `solve_cw`.
The closer the input frequency is to the resonant-mode frequency, the faster `solve_eig` should converge. Instead of using the source frequency, you can instead pass a `guessfreq` argument to `solve_eigfreq` specifying an input frequency (which may even be complex).
Technically, `solve_eig` is using a [shift-and-invert power iteration](https://en.wikipedia.org/wiki/Inverse_iteration) to compute the resonant mode, as reviewed in [Frequency-Domain Eigensolver](Eigensolver_Math.md).
As for `solve_cw` above, you are required to set `force_complex_fields=True` to use `solve_eigfreq`.
### GDSII Support
This feature is only available if Meep is built with [libGDSII](Build_From_Source.md#libgdsii). It so, then the following functions are available:
@@ GDSII_layers @@
@@ GDSII_prisms @@
@@ GDSII_vol @@
### Data Visualization
This module provides basic visualization functionality for the simulation domain. The intent of the module is to provide functions that can be called with *no customization options whatsoever* and will do useful relevant things by default, but which can also be customized in cases where you *do* want to take the time to spruce up the output. The `Simulation` class provides the following methods:
@@ Simulation.plot2D @@
@@ Simulation.plot3D @@
@@ Simulation.visualize_chunks @@
An animated visualization is also possible via the [Animate2D](#Animate2D) class.
<a id="run-functions"></a>
### Run and Step Functions
The actual work in Meep is performed by `run` functions, which time-step the simulation for a given amount of time or until a given condition is satisfied. These are attributes of the `Simulation` class.
The run functions, in turn, can be modified by use of [step functions](#predefined-step-functions): these are called at every time step and can perform any arbitrary computation on the fields, do outputs and I/O, or even modify the simulation. The step functions can be transformed by many [modifier functions](#step-function-modifiers), like `at_beginning`, `during_sources`, etcetera which cause them to only be called at certain times, etcetera, instead of at every time step.
A common point of confusion is described in [The Run Function Is Not A Loop](The_Run_Function_Is_Not_A_Loop.md). Read this article if you want to make Meep do some customized action on each time step, as many users make the same mistake. What you really want to in that case is to write a step function, as described below.
@@ Simulation.run @@
In particular, a useful value for `until_after_sources` or `until` is often `stop_when_field_decayed`, which is demonstrated in [Tutorial/Basics](Python_Tutorials/Basics.md#transmittance-spectrum-of-a-waveguide-bend). These top-level functions are available:
@@ stop_when_fields_decayed @@
@@ stop_when_energy_decayed @@
@@ stop_when_dft_decayed @@
@@ stop_after_walltime @@
@@ stop_on_interrupt @@
Finally, another run function, useful for computing $\omega(\mathbf{k})$ band diagrams, is available via these `Simulation` methods:
@@ Simulation.run_k_points @@
@@ Simulation.run_k_point @@
### Predefined Step Functions
Several useful step functions are predefined by Meep. These are available directly via the `meep` package but require a `Simulation` instance as an argument.
#### Output Functions
The most common step function is an output function, which outputs some field component to an [HDF5](https://en.wikipedia.org/wiki/HDF5) file. Normally, you will want to modify this by one of the `at_*` functions, below, as outputting a field at *every* time step can get quite time- and storage-consuming.
Note that although the various field components are stored at different places in the [Yee lattice](Yee_Lattice.md), when they are outputted they are all linearly interpolated to the same grid: to the points at the *centers* of the Yee cells, i.e. $(i+0.5,j+0.5,k+0.5)\cdotΔ$ in 3d.
@@ Simulation.output_dft @@
@@ output_epsilon @@
@@ output_mu @@
@@ output_poynting @@
@@ output_hpwr @@
@@ output_dpwr @@
@@ output_tot_pwr @@
@@ output_png @@
@@ output_hfield @@
@@ output_hfield_x @@
@@ output_hfield_y @@
@@ output_hfield_z @@
@@ output_hfield_r @@
@@ output_hfield_p @@
@@ output_bfield @@
@@ output_bfield_x @@
@@ output_bfield_y @@
@@ output_bfield_z @@
@@ output_bfield_r @@
@@ output_bfield_p @@
@@ output_efield @@
@@ output_efield_x @@
@@ output_efield_y @@
@@ output_efield_z @@
@@ output_efield_r @@
@@ output_efield_p @@
@@ output_dfield @@
@@ output_dfield_x @@
@@ output_dfield_y @@
@@ output_dfield_z @@
@@ output_dfield_r @@
@@ output_dfield_p @@
@@ output_sfield @@
@@ output_sfield_x @@
@@ output_sfield_y @@
@@ output_sfield_z @@
@@ output_sfield_r @@
@@ output_sfield_p @@
More generally, it is possible to output an arbitrary function of position and zero or more field components, similar to the `Simulation.integrate_field_function` method, described above. This is done by:
@@ Simulation.output_field_function @@
See also [Field Functions](Field_Functions.md), and [Synchronizing the Magnetic and Electric Fields](Synchronizing_the_Magnetic_and_Electric_Fields.md) if you want to do computations combining the electric and magnetic fields.
#### Array Slices
The output functions described above write the data for the fields and materials for the entire cell to an HDF5 file. This is useful for post-processing large datasets which may not fit into memory as you can later read in the HDF5 file to obtain field/material data as a NumPy array. However, in some cases it is convenient to bypass the disk altogether to obtain the data *directly* in the form of a NumPy array without writing/reading HDF5 files. Additionally, you may want the field/material data on just a subregion (or slice) of the entire volume. This functionality is provided by the `get_array` method which takes as input a subregion of the cell and the field/material component. The method returns a NumPy array containing values of the field/material at the current simulation time.
@@ Simulation.get_array @@
@@ Simulation.get_dft_array @@
#### Array Metadata
@@ Simulation.get_array_metadata @@
The following are some examples of how array metadata can be used.
**Labeling Axes in Plots of Grid Quantities**
```python
# using the geometry from the bend-flux tutorial example
import matplotlib.pyplot as plt
import numpy as np
eps_array=sim.get_epsilon()
(x,y,z,w)=sim.get_array_metadata()
plt.figure()
ax = plt.subplot(111)
plt.pcolormesh(x,y,np.transpose(eps_array),shading='gouraud')
ax.set_aspect('equal')
plt.show()
```
![](images/PermittivityWithLabeledAxes.png)
**Computing Quantities Defined by Integrals of Field-Dependent Functions Over Grid Regions**
+ energy stored in the $\mathbf{E}$-field in a region $\mathcal{V}$:
$$ \mathcal{E}=
\frac{1}{2}\int_{\mathcal V} \epsilon |\mathbf{E}|^2\,dV
$$
+ Poynting flux through a surface $\mathcal{S}$:
$$\mathcal{S}=\frac{1}{2}\text{Re }\int_{\mathcal S}
\Big(\mathbf{E}^*\times \mathbf{H}\Big)\times d\mathbf{A}
$$
```python
import numpy as np
# E-field modal volume in box from time-domain fields
box = mp.Volume(center=box_center, size=box_size)
(Ex,Ey,Ez) = [sim.get_array(vol=box, component=c, cmplx=True) for c in [mp.Ex, mp.Ey, mp.Ez]]
eps = sim.get_array(vol=box, component=mp.Dielectric)
(x,y,z,w) = sim.get_array_metadata(vol=box)
energy_density = np.real(eps*(np.conj(Ex)*Ex + np.conj(Ey)*Ey + np.conj(Ez)*Ez)) # array
energy = np.sum(w*energy_density) # scalar
# x-directed Poynting flux through monitor from frequency-domain fields
monitor = mp.FluxRegion(center=mon_center, size=mon_size)
dft_cell = sim.add_flux(freq, 0, 1, monitor)
sim.run(...) # timestep until DFTs converged
(Ey,Ez,Hy,Hz) = [sim.get_dft_array(dft_cell,c,0) for c in [mp.Ey, mp.Ez, mp.Hy, mp.Hz]]
(x,y,z,w) = sim.get_array_metadata(dft=dft_cell)
flux_density = np.real( np.conj(Ey)*Hz - np.conj(Ez)*Hy ) # array
flux = np.sum(w*flux_density) # scalar
```
#### Source Slices
@@ Simulation.get_source @@
#### Harminv Step Function
The following step function collects field data from a given point and runs [Harminv](https://github.com/NanoComp/harminv) on that data to extract the frequencies, decay rates, and other information.
* [Harminv class](#Harminv)
### Step-Function Modifiers
Rather than writing a brand-new step function every time something a bit different is required, the following "modifier" functions take a bunch of step functions and produce *new* step functions with modified behavior. See also [Tutorial/Basics](Python_Tutorials/Basics.md) for examples.
#### Miscellaneous Step-Function Modifiers
@@ combine_step_funcs @@
@@ synchronized_magnetic @@
#### Controlling When a Step Function Executes
@@ when_true @@
@@ when_false @@
@@ at_every @@
@@ after_time @@
@@ before_time @@
@@ at_time @@
@@ after_sources @@
@@ after_sources_and_time @@
@@ during_sources @@
@@ at_beginning @@
@@ at_end @@
#### Modifying HDF5 Output
@@ in_volume @@
@@ in_point @@
@@ to_appended @@
@@ with_prefix @@
### Writing Your Own Step Functions
A step function can take two forms. The simplest is just a function with one argument (the simulation instance), which is called at every time step unless modified by one of the modifier functions above. e.g.
```py
def my_step(sim):
print("Hello world!")
```
If one then does `sim.run(my_step, until=100)`, Meep will run for 100 time units and print "Hello world!" at every time step.
This suffices for most purposes. However, sometimes you need a step function that opens a file, or accumulates some computation, and you need to clean up (e.g. close the file or print the results) at the end of the run. For this case, you can write a step function of two arguments: the second argument will either be `step` when it is called during time-stepping, or `finish` when it is called at the end of the run:
```py
def my_step(sim, todo):
if todo == 'step':
# do something
elif todo == 'finish':
# do something else
# access simulation attributes
sim.fields ...etc.
```
Low-Level Functions
-------------------
By default, Meep initializes C++ objects like `meep::structure` and `meep::fields` in the `Simulation` object based on attributes like `sources` and `geometry`. Theses objects are then accessible via `simulation_instance.structure` and `simulation_instance.fields`. Given these, you can then call essentially any function in the C++ interface, because all of the C++ functions are automatically made accessible to Python by the wrapper-generator program [SWIG](https://en.wikipedia.org/wiki/SWIG).
### Initializing the Structure and Fields
The `structure` and `fields` variables are automatically initialized when any of the run functions is called, or by various other functions such as `add_flux`. To initialize them separately, you can call `Simulation.init_sim()` manually, or `Simulation._init_structure(k_point)` to just initialize the structure.
If you want to time step more than one field simultaneously, the easiest way is probably to do something like:
```py
sim = Simulation(cell_size, resolution).init_sim()
my_fields = sim.fields
sim.fields = None
sim.reset_meep()
```
and then change the geometry etc. and re-run `sim.init_sim()`. Then you'll have two field objects in memory.
### SWIG Wrappers
If you look at a function in the C++ interface, then there are a few simple rules to infer the name of the corresponding Python function.
- First, all functions in the `meep::` namespace are available in the Meep Python module from the top-level `meep` package.
- Second, any method of a class is accessible via the standard Python class interface. For example, `meep::fields::step`, which is the function that performs a time-step, is exposed to Python as `fields_instance.step()` where a fields instance is usually accessible from Simulation.fields.
- C++ constructors are called using the normal Python class instantiation. E.g., `fields = mp.fields(...)` returns a new `meep::fields` object. Calling destructors is not necessary because objects are automatically garbage collected.
Some argument type conversion is performed automatically, e.g. types like complex numbers are converted to `complex<double>`, etcetera. `Vector3` vectors are converted to `meep::vec`, but to do this it is necessary to know the dimensionality of the problem in C++. The problem dimensions are automatically initialized by `Simulation._init_structure`, but if you want to pass vector arguments to C++ before that time you should call `Simulation.require_dimensions()`, which infers the dimensions from the `cell_size`, `k_point`, and `dimensions` variables.
<a id="classes"></a>
Class Reference
---------------
Classes are complex datatypes with various properties which may have default values. Classes can be "subclasses" of other classes. Subclasses inherit all the properties of their superclass and can be used in any place the superclass is expected.
The `meep` package defines several types of classes. The most important of these is the `Simulation` class. Classes which are available directly from the `meep` package are constructed with:
```py
mp.ClassName(prop1=val1, prop2=val2, ...)
```
The most numerous are the geometric object classes which are the same as those used in [MPB](https://mpb.readthedocs.io). You can get a list of the available classes (and constants) in the Python interpreter with:
```py
import meep
[x for x in dir(meep) if x[0].isupper()]
```
More information, including their property types and default values, is available with the standard python `help` function: `help(mp.ClassName)`.
The following classes are available directly via the `meep` package.
@@ Medium[methods-with-docstrings] @@
@@ MaterialGrid[methods-with-docstrings] @@
@@ Susceptibility[methods-with-docstrings] @@
@@ LorentzianSusceptibility[methods-with-docstrings] @@
@@ DrudeSusceptibility[methods-with-docstrings] @@
@@ MultilevelAtom[methods-with-docstrings] @@
@@ Transition[methods-with-docstrings] @@
@@ NoisyLorentzianSusceptibility[methods-with-docstrings] @@
@@ NoisyDrudeSusceptibility[methods-with-docstrings] @@
@@ GyrotropicLorentzianSusceptibility[methods-with-docstrings] @@
@@ GyrotropicDrudeSusceptibility[methods-with-docstrings] @@
@@ GyrotropicSaturatedSusceptibility[methods-with-docstrings] @@
@@ Vector3[methods-with-docstrings] @@
@@ GeometricObject[methods-with-docstrings] @@
@@ Sphere[methods-with-docstrings] @@
@@ Cylinder[methods-with-docstrings] @@
@@ Wedge[methods-with-docstrings] @@
@@ Cone[methods-with-docstrings] @@
@@ Block[methods-with-docstrings] @@
@@ Ellipsoid[methods-with-docstrings] @@
@@ Prism[methods-with-docstrings] @@
@@ Matrix[methods-with-docstrings] @@
**Related function:**
@@ get_rotation_matrix @@
@@ Symmetry[methods-with-docstrings] @@
@@ Rotate2[methods-with-docstrings] @@
@@ Rotate4[methods-with-docstrings] @@
@@ Mirror[methods-with-docstrings] @@
@@ Identity[methods-with-docstrings] @@
@@ PML[methods-with-docstrings] @@
@@ Absorber[methods-with-docstrings] @@
@@ Source[methods-with-docstrings] @@
@@ SourceTime @@
@@ EigenModeSource[methods-with-docstrings] @@
@@ GaussianBeamSource[methods-with-docstrings] @@
@@ ContinuousSource[methods-with-docstrings] @@
@@ GaussianSource[methods-with-docstrings] @@
@@ CustomSource[methods-with-docstrings] @@
@@ FluxRegion[methods-with-docstrings] @@
@@ EnergyRegion[methods-with-docstrings] @@
@@ ForceRegion[methods-with-docstrings] @@
@@ Volume[methods-with-docstrings] @@
** Related function:**
@@ get_center_and_size @@
@@ Animate2D[methods-with-docstrings] @@
@@ Harminv[methods-with-docstrings] @@
@@ Verbosity @@
@@ Verbosity.__init__ @@
@@ Verbosity.__call__ @@
@@ Verbosity.add_verbosity_var @@
@@ Verbosity.get @@
@@ Verbosity.set @@
@@ BinaryPartition[methods-with-docstrings] @@
Miscellaneous Functions Reference
---------------------------------
@@ interpolate @@
#### Flux functions
@@ get_flux_freqs @@
@@ get_fluxes @@
@@ scale_flux_fields @@
@@ get_eigenmode_freqs @@
#### Energy Functions
@@ get_energy_freqs @@
@@ get_electric_energy @@
@@ get_magnetic_energy @@
@@ get_total_energy @@
#### Force Functions
@@ get_force_freqs @@
@@ get_forces @@
#### LDOS Functions
@@ Ldos @@
@@ get_ldos_freqs @@
@@ dft_ldos @@
#### Near2Far Functions
@@ get_near2far_freqs @@
@@ scale_near2far_fields @@
#### GDSII Functions
@@ GDSII_layers @@
@@ GDSII_prisms @@
@@ GDSII_vol @@
#### Run and Step Functions
@@ stop_when_fields_decayed @@
@@ stop_after_walltime @@
@@ stop_on_interrupt @@
#### Output Functions
@@ output_epsilon @@
@@ output_mu @@
@@ output_poynting @@
@@ output_hpwr @@
@@ output_dpwr @@
@@ output_tot_pwr @@
@@ output_png @@
@@ output_hfield @@
@@ output_hfield_x @@
@@ output_hfield_y @@
@@ output_hfield_z @@
@@ output_hfield_r @@
@@ output_hfield_p @@
@@ output_bfield @@
@@ output_bfield_x @@
@@ output_bfield_y @@
@@ output_bfield_z @@
@@ output_bfield_r @@
@@ output_bfield_p @@
@@ output_efield @@
@@ output_efield_x @@
@@ output_efield_y @@
@@ output_efield_z @@
@@ output_efield_r @@
@@ output_efield_p @@
@@ output_dfield @@
@@ output_dfield_x @@
@@ output_dfield_y @@
@@ output_dfield_z @@
@@ output_dfield_r @@
@@ output_dfield_p @@
@@ output_sfield @@
@@ output_sfield_x @@
@@ output_sfield_y @@
@@ output_sfield_z @@
@@ output_sfield_r @@
@@ output_sfield_p @@