{
    "info": {
        "author": "Floris laporte",
        "author_email": "floris.laporte@gmail.com",
        "bugtrack_url": null,
        "classifiers": [
            "Development Status :: 3 - Alpha",
            "License :: OSI Approved :: MIT License",
            "Operating System :: OS Independent",
            "Programming Language :: Python :: 3",
            "Topic :: Scientific/Engineering",
            "Topic :: Scientific/Engineering :: Physics"
        ],
        "description": "# Python 3D FDTD Simulator\n\nA 3D electromagnetic FDTD simulator written in Python. The FDTD simulator has\nan optional PyTorch backend, enabling FDTD simulations on a GPU.\n\n**NOTE: This library is under construction. Only some minimal features are implemented and\nthe API might change considerably.**\n\n## Installation\nThe `fdtd`-library can be installed with `pip`:\n```\npip install fdtd\n```\nThe development version can be installed by cloning the repository\n```\ngit clone http://github.com/flaport/fdtd\n```\nand linking it with pip\n```\npip install -e fdtd\n```\n\n## Dependencies\n- python 3.6+\n- numpy\n- matplotlib\n- tqdm\n- pytorch (optional)\n\n## Contributing\nThe library is still in a very early stage of development, but all improvements\nor additions (for example new objects, sources or detectors) are welcome.\nPlease make a pull-request \ud83d\ude0a.\n\n## Introduction to the library\n\n### Imports\n\nThe `fdtd` library is simply imported as follows:\n\n\n```python\nimport fdtd\n```\n\n### Setting the backend\n\nThe `fdtd` library allows to choose a backend. The `\"numpy\"` backend is the\ndefault one, but there are also several additional PyTorch backends:\n- `\"numpy\"` (defaults to float64 arrays)\n- `\"torch\"` (defaults to float64 tensors)\n- `\"torch.float32\"`\n- `\"torch.float64\"`\n- `\"torch.cuda\"` (defaults to float64 tensors)\n- `\"torch.cuda.float32\"`\n- `\"torch.cuda.float64\"`\n\nFor example, this is how to choose the `\"torch\"` backend:\n```python\nfdtd.set_backend(\"torch\")\n```\n\nIn general, the `\"numpy\"` backend is preferred for standard CPU calculations\nwith `\"float64\"` precision. In general, `\"float64\"` precision is always\npreferred over `\"float32\"` for FDTD simulations, however, `\"float32\"` might\ngive a significant performance boost.\n\nThe `\"cuda\"` backends are only available for computers with a GPU.\n\n\n### The FDTD-grid\n\nThe FDTD grid defines the simulation region. \n```python\n# signature\nfdtd.Grid(\n    shape: Tuple[Number, Number, Number],\n    grid_spacing: float = 155e-9,\n    permittivity: float = 1.0,\n    permeability: float = 1.0,\n    courant_number: float = None,\n)\n```\n\nA grid is defined by its `shape`, which is just a 3D tuple of `Number`-types\n(integers or floats). If the shape is given in floats, it denotes the width,\nheight and length of the grid in meters. If the shape is given in integers, it\ndenotes the width, height and length of the grid in terms of the\n`grid_spacing`. Internally, these numbers will be translated to three integers:\n`grid.Nx`, `grid.Ny` and `grid.Nz`.\n\nA `grid_spacing` can be given. For stability reasons, it is recommended to\nchoose a grid spacing that is at least 10 times smaller than the *smallest*\nwavelength in the grid. This means that for a grid containing a source with\nwavelength `1550nm` and a material with refractive index of `3.1`, the\nrecommended minimum `grid_spacing` turns out to be `50pm`\n\nFor the `permittivity` and `permeability` floats or arrays with the following\nshapes\n\n- `(grid.Nx, grid.Ny, grid.Nz)`\n- or `(grid.Nx, grid.Ny, grid.Nz, 1)`\n- or `(grid.Nx, grid.Ny, grid.Nz, 3)`\n\nare expected. In the last case, the shape implies the possibility for different\npermittivity for each of the major axes (so-called *uniaxial* or *biaxial*\nmaterials).  Internally, these variables will be converted (for performance\nreasons) to their inverses `grid.inverse_permittivity` array and a\n`grid.inverse_permeability` array of shape `(grid.Nx, grid.Ny, grid.Nz, 3)`. It\nis possible to change those arrays after making the grid.\n\nFinally, the `courant_number` of the grid determines the relation between the\n`time_step` of the simulation and the `grid_spacing` of the grid. If not given,\nit is chosen to be the maximum number allowed by the [Courant-Friedrichs-Lewy\nCondition](https://en.wikipedia.org/wiki/Courant\u2013Friedrichs\u2013Lewy_condition):\n`1` for `1D` simulations, `1/\u221a2` for `2D` simulations and `1/\u221a3` for `3D`\nsimulations (the dimensionality will be derived by the shape of the grid). For\nstability reasons, it is recommended not to change this value.\n\n\n```python\ngrid = fdtd.Grid(\n    shape = (25e-6, 15e-6, 1), # 25um x 15um x 1 (grid_spacing) --> 2D FDTD\n)\n\nprint(grid)\n```\n\n    Grid(shape=(161,97,1), grid_spacing=1.55e-07, courant_number=0.70)\n\n\n\n### Adding an object to the grid\n\nAn other option to locally change the `permittivity` or `permeability` in the\ngrid is to add an `Object` to the grid.\n```python\n# signature\nfdtd.Object(\n    permittivity: Tensorlike,\n    name: str = None\n)\n```\nAn object defines a part of the grid with modified update equations, allowing\nto introduce for example absorbing materials or biaxial materials for which\nmixing between the axes are present through `Pockels coefficients` or many\nmore. In this case we'll make an object with a different `permittivity` than\nthe grid it is in.\n\nJust like for the grid, the `Object` expects a `permittivity` to be a floats or\nan array of the following possible shapes\n\n- `(obj.Nx, obj.Ny, obj.Nz)`\n- or `(obj.Nx, obj.Ny, obj.Nz, 1)`\n- or `(obj.Nx, obj.Ny, obj.Nz, 3)`\n\nNote that the values `obj.Nx`, `obj.Ny` and `obj.Nz` are not given to the\nobject constructor. They are in stead derived from its placing in the grid:\n\n\n```python\ngrid[11:32, 30:84, 0] = fdtd.Object(permittivity=1.7**2, name=\"object\")\n```\n\nSeveral things happen here. First of all, the object is given the space\n`[11:32, 30:84, 0]` in the grid. Because it is given this space, the object's\n`Nx`, `Ny` and `Nz` are automatically set. Furthermore, by supplying a name to\nthe object, this name will become available in the grid:\n\n\n```python\nprint(grid.object)\n```\n\n        Object(name='object')\n            @ x=11:32, y=30:84, z=0:1\n\n\n\nA second object can be added to the grid:\n\n\n```python\ngrid[13e-6:18e-6, 5e-6:8e-6, 0] = fdtd.Object(permittivity=1.5**2)\n```\n\nHere, a slice with floating point numbers was chosen. These floats will be\nreplaced by integer `Nx`, `Ny` and `Nz` during the registration of the object.\nSince the object did not receive a name, the object won't be available as an\nattribute of the grid. However, it is still available via the `grid.objects`\nlist:\n\n\n```python\nprint(grid.objects)\n```\n\n    [Object(name='object'), Object(name=None)]\n\n\nThis list stores all objects (i.e. of type `fdtd.Object`) in the order that\nthey were added to the grid.\n\n### Adding a source to the grid\n\nSimilarly as to adding an object to the grid, an `fdtd.LineSource` can also be\nadded:\n```python\n# signature\nfdtd.LineSource(\n    period: Number = 15, # timesteps or seconds\n    power: float = 1.0,\n    phase_shift: float = 0.0,\n    name: str = None,\n)\n```\n\n\nAnd also just like an `fdtd.Object`, an `fdtd.Source` size is defined by its\nplacement on the grid:\n\n\n```python\ngrid[7.5e-6:8.0e-6, 11.8e-6:13.0e-6, 0] = fdtd.LineSource(\n    period = 1550e-9 / (3e8), name=\"source\"\n)\n```\n\nHowever, it is important to note that in this case a `LineSource` is added to\nthe grid, i.e. the source spans the diagonal of the cube defined by the slices.\nInternally, these slices will be converted into lists to ensure this behavior:\n\n\n```python\nprint(grid.source)\n```\n\n        LineSource(period=14, power=1.0, phase_shift=0.0, name='source')\n            @ x=[48, ... , 51], y=[76, ... , 83], z=[0, ... , 0]\n\n\n\nNote that one could also have supplied lists to index the grid in the first\nplace. This feature could be useful to create a `LineSource` of arbitrary\nshape.\n\n### Adding a detector to the grid\n\n```python\n# signature\nfdtd.LineDetector(\n    name=None\n)\n```\n\nAdding a detector to the grid works the same as adding a source\n\n\n```python\ngrid[12e-6, :, 0] = fdtd.LineDetector(name=\"detector\")\n```\n\n\n```python\nprint(grid.detector)\n```\n\n        LineDetector(name='detector')\n            @ x=[77, ... , 77], y=[0, ... , 96], z=[0, ... , 0]\n\n\n\n### Adding grid boundaries\n\n```python\n# signature\nfdtd.PML(\n    a: float = 1e-8, # stability factor\n    name: str = None\n)\n```\n\nAlthough, having an object, source and detector to simulate is in principle\nenough to perform an FDTD simulation, One also needs to define a grid boundary\nto prevent the fields to be reflected. One of those boundaries that can be\nadded to the grid is a [Perfectly Matched\nLayer](https://en.wikipedia.org/wiki/Perfectly_matched_layer) or `PML`. These\nare basically absorbing boundaries.\n\n\n\n```python\n# x boundaries\ngrid[0:10, :, :] = fdtd.PML(name=\"pml_xlow\")\ngrid[-10:, :, :] = fdtd.PML(name=\"pml_xhigh\")\n\n# y boundaries\ngrid[:, 0:10, :] = fdtd.PML(name=\"pml_ylow\")\ngrid[:, -10:, :] = fdtd.PML(name=\"pml_yhigh\")\n```\n\n### Grid summary\n\nA simple summary of the grid can be shown by printing out the grid:\n\n\n```python\nprint(grid)\n```\n\n    Grid(shape=(161,97,1), grid_spacing=1.55e-07, courant_number=0.70)\n\n    sources:\n        LineSource(period=14, power=1.0, phase_shift=0.0, name='source')\n            @ x=[48, ... , 51], y=[76, ... , 83], z=[0, ... , 0]\n\n    detectors:\n        LineDetector(name='detector')\n            @ x=[77, ... , 77], y=[0, ... , 96], z=[0, ... , 0]\n\n    boundaries:\n        PML(name='pml_xlow')\n            @ x=0:10, y=:, z=:\n        PML(name='pml_xhigh')\n            @ x=-10:, y=:, z=:\n        PML(name='pml_ylow')\n            @ x=:, y=0:10, z=:\n        PML(name='pml_yhigh')\n            @ x=:, y=-10:, z=:\n\n    objects:\n        Object(name='object')\n            @ x=11:32, y=30:84, z=0:1\n        Object(name=None)\n            @ x=84:116, y=32:52, z=0:1\n\n\n\n### Running a simulation\nRunning a simulation is as simple as using the `grid.run` method.\n```python\ngrid.run(\n    total_time: Number,\n    progress_bar: bool = True\n)\n```\nJust like for the lengths in the grid, the `total_time` of the simulation\ncan be specified as an integer (number of `time_steps`) or as a float (in\nseconds).\n\n\n```python\ngrid.run(total_time=100)\n```\n\n### Grid visualization\n\nLet's visualize the grid. This can be done with the `grid.visualize` method:\n\n```python\n# signature\ngrid.visualize(\n    grid,\n    x=None,\n    y=None,\n    z=None,\n    cmap=\"Blues\",\n    pbcolor=\"C3\",\n    pmlcolor=(0, 0, 0, 0.1),\n    objcolor=(1, 0, 0, 0.1),\n    srccolor=\"C0\",\n    detcolor=\"C2\",\n    show=True,\n)\n```\n\nThis method will by default visualize all objects in the grid, as well as the\npower at the current `time_step` at a certain `x`, `y` **OR** `z`-plane. By\nsetting `show=False`, one can disable the immediate visualization of the\nmatplotlib image.\n\n\n```python\ngrid.visualize(z=0)\n```\n\n\n![png](images/grid.png) \n\n## Background\n\nAn as quick as possible explanation of the FDTD discretization of the Maxwell equations.\n\n### Update Equations\n\n\nAn electromagnetic FDTD solver solves the time-dependent Maxwell Equations\n\n```python\n    curl(H) = \u03b5*\u03b50*dE/dt\n    curl(E) = -\u00b5*\u00b50*dH/dt\n```\n\nThese two equations are called *Ampere's Law* and *Faraday's Law* respectively.\n\nIn these equations, \u03b5 and \u00b5 are the relative permittivity and permeability tensors\nrespectively. \u03b50 and \u00b50 are the vacuum permittivity and permeability and their\nsquare root can be absorbed into E and H respectively, such that `E := \u221a\u03b50*E`\nand `H := \u221a\u00b50*H`.\n\nDoing this, the Maxwell equations can be written as update equations:\n```python\n    E  += c*dt*inv(\u03b5)*curl(H)\n    H  -= c*dt*inv(\u00b5)*curl(E)\n```\nThe electric and magnetic field can then be discretized on a grid with\ninterlaced Yee-coordinates, which in 3D looks like this:\n\n![grid discretization in 3D](images/yee.svg)\n\nAccording to the Yee discretization algorithm, there are inherently two types\nof fields on the grid: `E`-type fields on integer grid locations and `H`-type\nfields on half-integer grid locations.\n\nThe beauty of these interlaced coordinates is that they enable a very natural\nway of writing the curl of the electric and magnetic fields: the curl of an\nH-type field will be an E-type field and vice versa.\n\nThis way, the curl of E can be written as\n```python\n    curl(E)[m,n,p] = (dEz/dy - dEy/dz, dEx/dz - dEz/dx, dEy/dx - dEx/dy)[m,n,p]\n                   =( ((Ez[m,n+1,p]-Ez[m,n,p])/dy - (Ey[m,n,p+1]-Ey[m,n,p])/dz),\n                      ((Ex[m,n,p+1]-Ex[m,n,p])/dz - (Ez[m+1,n,p]-Ez[m,n,p])/dx),\n                      ((Ey[m+1,n,p]-Ey[m,n,p])/dx - (Ex[m,n+1,p]-Ex[m,n,p])/dy) )\n                   =(1/du)*( ((Ez[m,n+1,p]-Ez[m,n,p]) - (Ey[m,n,p+1]-Ey[m,n,p])), [assume dx=dy=dz=du]\n                             ((Ex[m,n,p+1]-Ex[m,n,p]) - (Ez[m+1,n,p]-Ez[m,n,p])),\n                             ((Ey[m+1,n,p]-Ey[m,n,p]) - (Ex[m,n+1,p]-Ex[m,n,p])) )\n\n```\nthis can be written efficiently with array slices (note that the factor\n`(1/du)` was left out):\n\n```python\ndef curl_E(E):\n    curl_E = np.zeros(E.shape)\n    curl_E[:,:-1,:,0] += E[:,1:,:,2] - E[:,:-1,:,2]\n    curl_E[:,:,:-1,0] -= E[:,:,1:,1] - E[:,:,:-1,1]\n\n    curl_E[:,:,:-1,1] += E[:,:,1:,0] - E[:,:,:-1,0]\n    curl_E[:-1,:,:,1] -= E[1:,:,:,2] - E[:-1,:,:,2]\n\n    curl_E[:-1,:,:,2] += E[1:,:,:,1] - E[:-1,:,:,1]\n    curl_E[:,:-1,:,2] -= E[:,1:,:,0] - E[:,:-1,:,0]\n    return curl_E\n```\n\nThe curl for H can be obtained in a similar way (note again that the factor\n`(1/du)` was left out):\n```python\ndef curl_H(H):\n    curl_H = np.zeros(H.shape)\n\n    curl_H[:,1:,:,0] += H[:,1:,:,2] - H[:,:-1,:,2]\n    curl_H[:,:,1:,0] -= H[:,:,1:,1] - H[:,:,:-1,1]\n\n    curl_H[:,:,1:,1] += H[:,:,1:,0] - H[:,:,:-1,0]\n    curl_H[1:,:,:,1] -= H[1:,:,:,2] - H[:-1,:,:,2]\n\n    curl_H[1:,:,:,2] += H[1:,:,:,1] - H[:-1,:,:,1]\n    curl_H[:,1:,:,2] -= H[:,1:,:,0] - H[:,:-1,:,0]\n    return curl_H\n```\n\nThe update equations can now be rewritten as\n```python\n    E  += (c*dt/du)*inv(\u03b5)*curl_H\n    H  -= (c*dt/du)*inv(\u00b5)*curl_E\n```\n\nThe number `(c*dt/du)` is a dimensionless parameter called the *courant number* `sc`. For\nstability reasons, the Courant number should always be smaller than `1/\u221aD`, with `D`\nthe dimension of the simulation. This can be intuitively be understood as the condition\nthat information should always travel slower than the speed of light through the grid.\nin the FDTD method described here, information can only travel to the neighbouring grid\ncells (through application of the curl). It would therefore take `D` time steps to\ntravel over the diagonal of a `D`-dimensional cube (square in `2D`, cube in `3D`), the\nCourant condition follows then automatically from the fact that the length of this\ndiagonal is `1/\u221aD`.\n\nThis yields the final update equations for the FDTD algorithm:\n\n\n```python\n    E  += sc*inv(\u03b5)*curl_H\n    H  -= sc*inv(\u00b5)*curl_E\n```\n\nThis is also how it is implemented:  \n\n```python\nclass Grid:\n    # ... [initialization]\n\n    def step(self):\n        self.update_E()\n        self.update_H()\n\n    def update_E(self):\n        self.E += self.courant_number * self.inverse_permittivity * curl_H(self.H)\n\n    def update_H(self):\n        self.H -= self.courant_number * self.inverse_permeability * curl_E(self.E)\n```\n\n\n### Sources\n\n\nAmpere's Law can be updated to incorporate a current density:\n```python\n    curl(H) = J + \u03b5*\u03b50*dE/dt\n```\nMaking again the usual substitutions `sc := c*dt/du`, `E := \u221a\u03b50*E` and `H := \u221a\u00b50*H`, the\nupdate equations can be modified to include the current density:\n\n```python\n    E += sc*inv(\u03b5)*curl_H - dt*inv(\u03b5)*J/\u221a\u03b50\n```\n\nMaking one final substitution `Es := -dt*inv(\u03b5)*J/\u221a\u03b50` allows us to write this in a\nvery clean way:\n```python\n    E += sc*inv(\u03b5)*curl_H + Es\n```\n\nWhere we defined Es as the *electric field source term*.\n\nIt is often useful to also define a *magnetic field source term* `Hs`, which would be\nderived from the *magnetic current density* if it were to exist. In the same way,\nFaraday's update equation can be rewritten as\n```python\n    H  -= sc*inv(\u00b5)*curl_E + Hs\n```\n\n```python\nclass Source:\n    # ... [initialization]\n    def update_E(self):\n        # electric source function here\n\n    def update_H(self):\n        # magnetic source function here\n\nclass Grid:\n    # ... [initialization]\n    def update_E(self):\n        # ... [electric field update equation]\n        for source in self.sources:\n            source.update_E()\n\n    def update_H(self):\n        # ... [magnetic field update equation]\n        for source in self.sources:\n            source.update_H()\n\n```\n\n### Lossy Medium\n\nWhen a material has a *electric conductivity* \u03c3e, a conduction-current will ensure that the medium is lossy. Ampere's law with a conduction current becomes\n```python\n    curl(H) = \u03c3e*E + \u03b5*\u03b50*dE/dt\n```\n\nmaking the usual substitutions, this becomes:\n```python\n    E(t+dt) - E(t) = sc*inv(\u03b5)*curl_H(t+dt/2) - dt*inv(\u03b5)*\u03c3e*E(t+dt/2)/\u03b50\n```\n\nThis update equation depends on the electric field on a half-integer time step (a *magnetic field timestep*). We need to make a substitution to interpolate the electric field to this time step:\n```python\n    (1 + 0.5*dt*inv(\u03b5)*\u03c3/\u221a\u03b50)*E(t+dt) = sc*inv(\u03b5)*curl_H(t+dt/2) + (1 - 0.5*dt*inv(\u03b5)*\u03c3e/\u03b50)*E(t)\n```\n\nWhich, after substitution `\u03c3 := inv(\u03b5)*\u03c3e/\u03b50` yield the new update equations:\n```python\n    f = 0.5*dt*\u03c3\n    E *= inv(1 + f) * (1 - f)\n    E += inv(1 + f)*sc*inv(\u03b5)*curl_H\n```\n\nIf we want to keep track of the absorbed energy:\n\nNote that the more complicated the permittivity tensor \u03b5 is, the more time consuming this\nalgorithm will be. It is therefore sometimes the right decision to transfer the absorption to the magnetic domain by introducing a (*nonphysical*) magnetic conductivity, because\nthe permeability tensor \u00b5 is usually just equal to one.\n\nWhich, after substitution `\u03c3 := inv(\u00b5)*\u03c3m/\u00b50`, we get the magnetic field update equations:\n```python\n    f = 0.5*dt*\u03c3\n    H *= inv(1 + f) * (1 - f)\n    H += inv(1 + f)*sc*inv(\u00b5)*curl_E\n```\n\n### Energy Density and Poynting Vector\nThe electromagnetic energy density can be given by\n```python\n    e = (1/2)*\u03b5*\u03b50*E**2 + (1/2)*\u00b5*\u00b50*H**2\n```\nmaking the above substitutions, this becomes in simulation units:\n```python\n    e = (1/2)*\u03b5*E**2 + (1/2)*\u00b5*H**2\n```\nThe Poynting vector is given by\n```python\n    P = E\u00d7H\n```\nWhich in simulation units becomes\n```python\n    P = c*E\u00d7H\n```\nThe energy introduced by a source `Es` can be derived from tracking the change in energy density\n```python\n    de = \u03b5*Es\u00b7E + (1/2)*\u03b5*Es**2\n```\nThis could also be derived from Poyntings energy conservation law:\n```python\n    de/dt = -grad(S) - J\u00b7E\n```\nwhere the first term just describes the redistribution of energy in a volume and the second term describes\nthe energy introduced by a current density.\n\nNote: although it is unphysical, one could also have introduced a magnetic source. This source would have introduced the following energy:\n```python\n    de = \u03b5*Hs\u00b7H + (1/2)*\u00b5*Hs**2\n```\nSince the \u00b5-tensor is usually just equal to one, using a magnetic source term is often more efficient.\n\nSimilarly, one can also keep track of the absorbed energy due to an electric conductivity in the following way:\n```python\n    f = 0.5*dt*\u03c3\n    Enoabs = E + sc*inv(\u03b5)*curl_H\n    E *= inv(1 + f) * (1 - f)\n    E += inv(1 + f)*sc*inv(\u03b5)*curl_H\n    dE = Enoabs - E\n    e_abs += \u03b5*E*dE + 0.5*\u03b5*dE**2\n```\n\nor if we want to keep track of the absorbed energy by magnetic a magnetic conductivity:\n```python\n    f = 0.5*dt*inv(\u00b5)*\u03c3\n    Hnoabs = E + sc*inv(\u00b5)*curl_E\n    H *= inv(1 + f) * (1 - f)\n    H += inv(1 + f)*sc*inv(\u00b5)*curl_E\n    dH = Hnoabs - H\n    e_abs += \u00b5*H*dH + 0.5*\u00b5*dH**2\n```\n\nThe electric term and magnetic term in the energy density are usually of the same size. Therefore, the same amount of energy will be absorbed by introducing a *magnetic conductivity* \u03c3m\nas by introducing a *electric conductivity* \u03c3e if:\n```python\n    inv(\u00b5)*\u03c3m/\u00b50 = inv(\u03b5)*\u03c3e/\u03b50\n```\n\n### Boundary Conditions\n\n#### Periodic Boundary Conditions\n\nAssuming we want periodic boundary conditions along the `X`-direction, then we have to\nmake sure that the fields at `Xlow` and `Xhigh` are the same. This has to be enforced\nafter performing the update equations:\n\nNote that the electric field `E` is dependent on `curl_H`, which means that\nthe first indices of `E` will not be updated through the update equations.\nIt's those indices that need to be set through the periodic boundary condition.\nConcretely: `E[0]` needs to be set to equal `E[-1]`. For the magnetic field, the\ninverse is true: `H` is dependent on `curl_E`, which means that its last indices\nwill not be set. This has to be done by the boundary condition: `H[-1]` needs to\nbe set equal to `H[0]`:\n\n```python\nclass PeriodicBoundaryX:\n    # ... [initialization]\n    def update_E(self):\n        self.grid.E[0, :, :, :] = self.grid.E[-1, :, :, :]\n\n    def update_H(self):\n        self.grid.H[-1, :, :, :] = self.grid.H[0, :, :, :]\n\nclass Grid:\n    # ... [initialization]\n    def update_E(self):\n        # ... [electric field update equation]\n        # ... [electric field source update equations]\n        for boundary in self.boundaries:\n            boundary.update_E()\n\n    def update_H(self):\n        # ... [magnetic field update equation]\n        # ... [magnetic field source update equations]\n        for boundary in self.boundaries:\n            boundary.update_H()\n```\n\n#### Perfectly Matched Layer\na Perfectly Matched Layer (PML) is the state of the art for\nintroducing absorbing boundary conditions in an FDTD grid.\nA PML is an impedance-matched absorbing area in the grid. It turns out that\nfor a impedance-matching condition to hold, the PML can only be absorbing in\na single direction. This is what makes a PML in fact a nonphysical material.\n\n\nConsider Ampere's law for the `Ez` component, where the usual substitutions\n`E := \u221a\u03b50*E`, `H := \u221a\u00b50*H` and `\u03c3 := inv(\u03b5)*\u03c3e/\u03b50` are\nalready introduced:\n```python\n    \u03b5*dEz/dt + \u03b5*\u03c3*Ez = c*dHy/dx - c*dHx/dy\n```\nThis becomes in the frequency domain:\n```python\n    i\u03c9*\u03b5*Ez + \u03b5*\u03c3*Ez = c*dHy/dx - c*dHx/dy\n```\nWe can split this equation in a x-propagating wave and a y-propagating wave:\n```python\n    i\u03c9*\u03b5*Ezx + \u03b5*\u03c3x*Ezx = i\u03c9*\u03b5*(1 + \u03c3x/i\u03c9)*Ezx = c*dHy/dx\n    i\u03c9*\u03b5*Ezy + \u03b5*\u03c3y*Ezy = i\u03c9*\u03b5*(1 + \u03c3y/i\u03c9)*Ezy = -c*dHx/dy\n```\n\nWe can define the `S`-operators as follows\n```python\n    Su = 1 + \u03c3u/i\u03c9          with u in {x, y, z}\n```\nIn general, we prefer to add a stability factor `au` and a scaling factor `ku` to `Su`:\n```python\n    Su = ku + \u03c3u/(i\u03c9+au)    with u in {x, y, z}\n```\nSumming the two equations for `Ez` back together after dividing by the respective `S`-operator gives\n```python\n    i\u03c9*\u03b5*Ez = (c/Sx)*dHy/dx - (c/Sy)*dHx/dy\n```\nConverting this back to the time domain gives\n```python\n    \u03b5*dEz/dt = c*sx[*]dHy/dx - c*sx[*]dHx/dy\n```\nwhere `sx` denotes the inverse Fourier transform of `(1/Sx)` and `[*]` denotes a convolution.\nThe expression for `su` can be proven [after some derivation] to look as follows:\n```python\n    su = (1/ku)*\u03b4(t) + Cu(t)    with u in {x, y, z}\n```\nwhere `\u03b4(t)` denotes the Dirac delta function and `C(t)` an exponentially\ndecaying function given by:\n```python\n    Cu(t) = -(\u03c3u/ku**2)*exp(-(au+\u03c3u/ku)*t)     for all t > 0 and u in {x, y, z}\n```\nPlugging this in gives:\n```python\n    dEz/dt = (c/kx)*inv(\u03b5)*dHy/dx - (c/ky)*inv(\u03b5)*dHx/dy + c*inv(\u03b5)*Cx[*]dHy/dx - c*inv(\u03b5)*Cx[*]dHx/dy\n           = (c/kx)*inv(\u03b5)*dHy/dx - (c/ky)*inv(\u03b5)*dHx/dy + c*inv(\u03b5)*\u0424ez/du      with du=dx=dy=dz\n```\nThis can be written as an update equation:\n```python\n    Ez += (1/kx)*sc*inv(\u03b5)*dHy - (1/ky)*sc*inv(\u03b5)*dHx + sc*inv(\u03b5)*\u0424ez\n```\nWhere we defined `\u0424eu` as\n```python\n    \u0424eu = \u03a8euv - \u03a8ezw           with u, v, w in {x, y, z}\n```\nand `\u03a8euv` as the convolution updating the component `Eu` by taking the derivative of `Hw` in the `v` direction:\n```python\n    \u03a8euv = dv*Cv[*]dHw/dv     with u, v, w in {x, y, z}\n```\nThis can be rewritten [after some derivation] as an update equation in itself:\n```python\n     \u03a8euv = bv*\u03a8euv + cv*dv*(dHw/dv)\n          = bv*\u03a8euv + cv*dHw            with u, v, w in {x, y, z}\n```\nWhere the constants `bu` and `cu` are derived to be:\n```python\n    bu = exp(-(au + \u03c3u/ku)*dt)              with u in {x, y, z}\n    cu = \u03c3u*(bu - 1)/(\u03c3u*ku + au*ku**2)     with u in {x, y, z}\n```\nThe final PML algorithm for the electric field now becomes:\n\n1. Update `\u0424e=[\u0424ex, \u0424ey, \u0424ez]` by using the update equation for the `\u03a8`-components.\n2. Update the electric fields the normal way\n3. Add `\u0424e` to the electric fields.\n\nor as python code:\n```python\nclass PML(Boundary):\n    # ... [initialization]\n    def update_phi_E(self): # update convolution\n        self.psi_Ex *= self.bE\n        self.psi_Ey *= self.bE\n        self.psi_Ez *= self.bE\n\n        c = self.cE\n        Hx = self.grid.H[self.locx]\n        Hy = self.grid.H[self.locy]\n        Hz = self.grid.H[self.locz]\n\n        self.psi_Ex[:, 1:, :, 1] += (Hz[:, 1:, :] - Hz[:, :-1, :]) * c[:, 1:, :, 1]\n        self.psi_Ex[:, :, 1:, 2] += (Hy[:, :, 1:] - Hy[:, :, :-1]) * c[:, :, 1:, 2]\n\n        self.psi_Ey[:, :, 1:, 2] += (Hx[:, :, 1:] - Hx[:, :, :-1]) * c[:, :, 1:, 2]\n        self.psi_Ey[1:, :, :, 0] += (Hz[1:, :, :] - Hz[:-1, :, :]) * c[1:, :, :, 0]\n\n        self.psi_Ez[1:, :, :, 0] += (Hy[1:, :, :] - Hy[:-1, :, :]) * c[1:, :, :, 0]\n        self.psi_Ez[:, 1:, :, 1] += (Hx[:, 1:, :] - Hx[:, :-1, :]) * c[:, 1:, :, 1]\n\n        self.phi_E[..., 0] = self.psi_Ex[..., 1] - self.psi_Ex[..., 2]\n        self.phi_E[..., 1] = self.psi_Ey[..., 2] - self.psi_Ey[..., 0]\n        self.phi_E[..., 2] = self.psi_Ez[..., 0] - self.psi_Ez[..., 1]\n\n    def update_E(self): # update PML located at self.loc\n        self.grid.E[self.loc] += (\n            self.grid.courant_number\n            * self.grid.inverse_permittivity[self.loc]\n            * self.phi_E\n        )\n\nclass Grid:\n    # ... [initialization]\n    def update_E(self):\n        for boundary in self.boundaries:\n            boundary.update_phi_E()\n        # ... [electric field update equation]\n        # ... [electric field source update equations]\n        for boundary in self.boundaries:\n            boundary.update_E()\n```\n\nThe same has to be applied for the magnetic field.\n\nThese update equations for the PML were based on [Schneider, Chap. 11](https://www.eecs.wsu.edu/~schneidj/ufdtd).\n\n## License\n\u00a9 Floris laporte - [MIT License](license)\n\n\n",
        "description_content_type": "text/markdown",
        "docs_url": null,
        "download_url": "",
        "downloads": {
            "last_day": -1,
            "last_month": -1,
            "last_week": -1
        },
        "home_page": "http://github.com/flaport/fdtd",
        "keywords": "",
        "license": "",
        "maintainer": "",
        "maintainer_email": "",
        "name": "fdtd",
        "package_url": "https://pypi.org/project/fdtd/",
        "platform": "",
        "project_url": "https://pypi.org/project/fdtd/",
        "project_urls": {
            "Homepage": "http://github.com/flaport/fdtd"
        },
        "release_url": "https://pypi.org/project/fdtd/0.0.0/",
        "requires_dist": null,
        "requires_python": "",
        "summary": "a 3D electromagnetic FDTD simulator written in Python",
        "version": "0.0.0"
    },
    "last_serial": 4926552,
    "releases": {
        "0.0.0": [
            {
                "comment_text": "",
                "digests": {
                    "md5": "6c45cd7bb7b96ca76bc78ee7174e556e",
                    "sha256": "c865379a6e5e3a993d8af050d4a6355356cd9ce5b150360a079fbd8e31e322fd"
                },
                "downloads": -1,
                "filename": "fdtd-0.0.0-py3-none-any.whl",
                "has_sig": false,
                "md5_digest": "6c45cd7bb7b96ca76bc78ee7174e556e",
                "packagetype": "bdist_wheel",
                "python_version": "py3",
                "requires_python": null,
                "size": 26546,
                "upload_time": "2019-03-11T19:04:14",
                "url": "https://files.pythonhosted.org/packages/ba/87/a2c62e92087cc39b1b694f870a2812577110d3bc1e46bec2c0e6d51f542d/fdtd-0.0.0-py3-none-any.whl"
            },
            {
                "comment_text": "",
                "digests": {
                    "md5": "5a54ced69cd9bd4b96d05b61f6d0f550",
                    "sha256": "398935001570663d969001b291ab4c2e9d256d9c0eb6c0f46aacb2f876f13e9a"
                },
                "downloads": -1,
                "filename": "fdtd-0.0.0.tar.gz",
                "has_sig": false,
                "md5_digest": "5a54ced69cd9bd4b96d05b61f6d0f550",
                "packagetype": "sdist",
                "python_version": "source",
                "requires_python": null,
                "size": 30925,
                "upload_time": "2019-03-11T19:04:17",
                "url": "https://files.pythonhosted.org/packages/00/95/495ce07e2c294c5f7d031215977ba0048a56aa06bb0f67744f0c65c58c95/fdtd-0.0.0.tar.gz"
            }
        ]
    },
    "urls": [
        {
            "comment_text": "",
            "digests": {
                "md5": "6c45cd7bb7b96ca76bc78ee7174e556e",
                "sha256": "c865379a6e5e3a993d8af050d4a6355356cd9ce5b150360a079fbd8e31e322fd"
            },
            "downloads": -1,
            "filename": "fdtd-0.0.0-py3-none-any.whl",
            "has_sig": false,
            "md5_digest": "6c45cd7bb7b96ca76bc78ee7174e556e",
            "packagetype": "bdist_wheel",
            "python_version": "py3",
            "requires_python": null,
            "size": 26546,
            "upload_time": "2019-03-11T19:04:14",
            "url": "https://files.pythonhosted.org/packages/ba/87/a2c62e92087cc39b1b694f870a2812577110d3bc1e46bec2c0e6d51f542d/fdtd-0.0.0-py3-none-any.whl"
        },
        {
            "comment_text": "",
            "digests": {
                "md5": "5a54ced69cd9bd4b96d05b61f6d0f550",
                "sha256": "398935001570663d969001b291ab4c2e9d256d9c0eb6c0f46aacb2f876f13e9a"
            },
            "downloads": -1,
            "filename": "fdtd-0.0.0.tar.gz",
            "has_sig": false,
            "md5_digest": "5a54ced69cd9bd4b96d05b61f6d0f550",
            "packagetype": "sdist",
            "python_version": "source",
            "requires_python": null,
            "size": 30925,
            "upload_time": "2019-03-11T19:04:17",
            "url": "https://files.pythonhosted.org/packages/00/95/495ce07e2c294c5f7d031215977ba0048a56aa06bb0f67744f0c65c58c95/fdtd-0.0.0.tar.gz"
        }
    ]
}