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_fields_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()
```

**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 @@
|