{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Dipole Emission in Light Emitting Diodes\n",
    "\n",
    "In addition to the two source time profiles of a [continuous wave](https://meep.readthedocs.io/en/latest/Python_User_Interface/#continuoussource) (CW) and finite-bandwidth [pulse](https://meep.readthedocs.io/en/latest/Python_User_Interface/#gaussiansource), Meep supports a [custom source](https://meep.readthedocs.io/en/latest/Python_User_Interface/#customsource) for defining an arbitrary time profile. This feature can be used to model **spatially incoherent** random (i.e., [white noise](https://en.wikipedia.org/wiki/White_noise)) dipole emission in a [light-emitting diode](https://en.wikipedia.org/wiki/Light-emitting_diode) (LED), spontaneously recombining excitonic emission in an [organic light-emitting diode](https://en.wikipedia.org/wiki/OLED) (OLED), as well as near-field heat transfer.   Such incoherent emission processes are very different from *coherent* emission by a single collection of sources with a fixed phase relation, and more complicated modeling techniques are required.\n",
    "\n",
    "This tutorial example involves computing the radiated [flux](https://meep.readthedocs.io/en/latest/Introduction/#transmittancereflectance-spectra) from $N=10$ dipole emitters of a 2d LED-like periodic structure with a thin (1d) light-emitting layer. A schematic of the unit cell geometry and simulation layout is shown below. A silver (Ag) back reflector is used to direct nearly all the flux upwards into the $+y$ direction. PML is used to terminate the air region at the top of the cell. (Note: PML is not necessary at the bottom of the cell due to the Ag layer which is effectively a lossless mirror with [skin depth](https://en.wikipedia.org/wiki/Skin_effect) of a few nanometers at a wavelength of 1 μm.) The emitting layer is placed within a lossless dielectric substrate with wavelength-independent refractive index of 3.45.\n",
    "\n",
    "![](https://meep.readthedocs.io/en/latest/images/LED_layout.png)\n",
    "\n",
    "One can take two different approaches to computing the radiated flux based on the type of emitter: (1) random or (2) deterministic. In Method 1 (brute-force Monte-Carlo), each emitter is a white-noise dipole: every timestep for every dipole is an independent random number. A single run involves all $N$ dipoles which are modeled using a `CustomSource`. The stochastic results for the radiated flux are averaged over multiple trials/iterations via [Monte-Carlo sampling](https://en.wikipedia.org/wiki/Monte_Carlo_method). Method 2 exploits the property of [linear time-invariance](https://en.wikipedia.org/wiki/Linear_time-invariant_system) of the materials/geometry and involves a sequence of $N$ separate runs each with a single deterministic dipole (i.e., pulse time profile, `GaussianSource`) at different positions in the emitting layer. Because dipoles at different positions are uncorrelated, the radiated flux from the ensemble is simply the average of all the individual iterations. The two approaches converge towards identical results, but Method 1 is more computationally expensive than Method 2 due to the much larger number of trials/iterations ($\\gg  N$) required to attain low noise variance.  (Even more sophisticated deterministic methods exist to reduce the number of separate simulations, especially at high resolutions; for example, replacing the point-dipole sources with a [rapidly converging set of smooth basis functions](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.81.012119), or fancier methods that exploit [trace-estimation methods](http://doi.org/10.1103/PhysRevB.92.134202) and/or transform volumetric sources to [surface sources](http://doi.org/10.1103/PhysRevB.88.054305).)\n",
    "\n",
    "*Note regarding normalization:* To directly compare the results for the radiated flux from the two methods, one might scale the spectrum from Method 2 in post processing to correct for the difference in spectrum between a Gaussian pulse and white noise. However, it is usually more convenient to *nondimensionalize* the results (for both methods) by dividing the flux spectrum for the textured surface with a reference spectrum computed by the *same method*, for example emission from a flat surface or in a homogeneous medium. This way, the details of the source spectra cancel automatically, and the nondimensionalized results can be compared as-is without any tricky scaling calculations. This kind of nondimensionalized comparison is useful to determine the *emission enhancement* (or suppression) of one structure relative to another as well as the *light-extraction efficiency* (the ratio of the radiated flux to the total flux emitted by the dipoles).  In order to compute the *absolute* (not relative) light emission by a particular structure, using either Method 1 or Method 2, one would need to rescale the output ([thanks to linearity](http://doi.org/10.1103/PhysRevLett.107.114302)) to convert the input spectrum (from white noise or Gaussian) to the actual emission spectrum (e.g. determined from the gain spectrum of a light-emitting diode)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import meep as mp\n",
    "from meep.materials import Ag\n",
    "import numpy as np\n",
    "\n",
    "import argparse\n",
    "\n",
    "parser = argparse.ArgumentParser()\n",
    "parser.add_argument(\"-res\", type=int, default=50, help=\"resolution (pixels/um)\")\n",
    "parser.add_argument(\n",
    "    \"-nr\", type=int, default=20, help=\"number of random trials (method 1)\"\n",
    ")\n",
    "parser.add_argument(\"-nd\", type=int, default=10, help=\"number of dipoles\")\n",
    "parser.add_argument(\"-nf\", type=int, default=500, help=\"number of frequencies\")\n",
    "parser.add_argument(\n",
    "    \"-textured\",\n",
    "    action=\"store_true\",\n",
    "    default=False,\n",
    "    help=\"flat (default) or textured surface\",\n",
    ")\n",
    "parser.add_argument(\n",
    "    \"-method\",\n",
    "    type=int,\n",
    "    choices=[1, 2],\n",
    "    default=1,\n",
    "    help=\"type of method: (1) random dipoles with nr trials or (2) single dipole with 1 run per dipole\",\n",
    ")\n",
    "args = parser.parse_args()\n",
    "\n",
    "resolution = args.res\n",
    "\n",
    "dpml = 1.0\n",
    "dair = 1.0\n",
    "hrod = 0.7\n",
    "wrod = 0.5\n",
    "dsub = 5.0\n",
    "dAg = 0.5\n",
    "\n",
    "sx = 1.1\n",
    "sy = dpml + dair + hrod + dsub + dAg\n",
    "\n",
    "cell_size = mp.Vector3(sx, sy)\n",
    "\n",
    "pml_layers = [mp.PML(direction=mp.Y, thickness=dpml, side=mp.High)]\n",
    "\n",
    "fcen = 1.0\n",
    "df = 0.2\n",
    "nfreq = args.nf\n",
    "ndipole = args.nd\n",
    "ntrial = args.nr\n",
    "run_time = 2 * nfreq / df\n",
    "\n",
    "geometry = [\n",
    "    mp.Block(\n",
    "        material=mp.Medium(index=3.45),\n",
    "        center=mp.Vector3(0, 0.5 * sy - dpml - dair - hrod - 0.5 * dsub),\n",
    "        size=mp.Vector3(mp.inf, dsub, mp.inf),\n",
    "    ),\n",
    "    mp.Block(\n",
    "        material=Ag,\n",
    "        center=mp.Vector3(0, -0.5 * sy + 0.5 * dAg),\n",
    "        size=mp.Vector3(mp.inf, dAg, mp.inf),\n",
    "    ),\n",
    "]\n",
    "\n",
    "if args.textured:\n",
    "    geometry.append(\n",
    "        mp.Block(\n",
    "            material=mp.Medium(index=3.45),\n",
    "            center=mp.Vector3(0, 0.5 * sy - dpml - dair - 0.5 * hrod),\n",
    "            size=mp.Vector3(wrod, hrod, mp.inf),\n",
    "        )\n",
    "    )\n",
    "\n",
    "\n",
    "def compute_flux(m=1, n=0):\n",
    "    if m == 1:\n",
    "        sources = []\n",
    "        for n in range(ndipole):\n",
    "            sources.append(\n",
    "                mp.Source(\n",
    "                    mp.CustomSource(src_func=lambda t: np.random.randn()),\n",
    "                    component=mp.Ez,\n",
    "                    center=mp.Vector3(\n",
    "                        sx * (-0.5 + n / ndipole), -0.5 * sy + dAg + 0.5 * dsub\n",
    "                    ),\n",
    "                )\n",
    "            )\n",
    "    else:\n",
    "        sources = [\n",
    "            mp.Source(\n",
    "                mp.GaussianSource(fcen, fwidth=df),\n",
    "                component=mp.Ez,\n",
    "                center=mp.Vector3(\n",
    "                    sx * (-0.5 + n / ndipole), -0.5 * sy + dAg + 0.5 * dsub\n",
    "                ),\n",
    "            )\n",
    "        ]\n",
    "\n",
    "    sim = mp.Simulation(\n",
    "        cell_size=cell_size,\n",
    "        resolution=resolution,\n",
    "        k_point=mp.Vector3(),\n",
    "        boundary_layers=pml_layers,\n",
    "        geometry=geometry,\n",
    "        sources=sources,\n",
    "    )\n",
    "\n",
    "    flux_mon = sim.add_flux(\n",
    "        fcen,\n",
    "        df,\n",
    "        nfreq,\n",
    "        mp.FluxRegion(center=mp.Vector3(0, 0.5 * sy - dpml), size=mp.Vector3(sx)),\n",
    "    )\n",
    "\n",
    "    sim.run(until=run_time)\n",
    "\n",
    "    flux = mp.get_fluxes(flux_mon)\n",
    "    freqs = mp.get_flux_freqs(flux_mon)\n",
    "\n",
    "    return freqs, flux\n",
    "\n",
    "\n",
    "if args.method == 1:\n",
    "    fluxes = np.zeros((nfreq, ntrial))\n",
    "    for t in range(ntrial):\n",
    "        freqs, fluxes[:, t] = compute_flux(m=1)\n",
    "else:\n",
    "    fluxes = np.zeros((nfreq, ndipole))\n",
    "    for d in range(ndipole):\n",
    "        freqs, fluxes[:, d] = compute_flux(m=2, n=d)\n",
    "\n",
    "\n",
    "if mp.am_master():\n",
    "    with open(\n",
    "        \"method{}_{}_res{}_nfreq{}_ndipole{}.npz\".format(\n",
    "            args.method,\n",
    "            \"textured\" if args.textured else \"flat\",\n",
    "            resolution,\n",
    "            nfreq,\n",
    "            ndipole,\n",
    "        ),\n",
    "        \"wb\",\n",
    "    ) as f:\n",
    "        np.savez(f, freqs=freqs, fluxes=fluxes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are five items to note in this script. (1) The frequency discretization of the flux spectrum must be sufficiently fine to resolve noisy features. In this example, the frequency range is 0.9 to 1.1 with spacing of 0.0004. (2) The runtime must be long enough for the DFT spectrum to resolve these oscillations. Due to the Fourier Uncertainty Principle, the runtime should be at least ~1/frequency resolution. Here, we found that a larger runtime of 2/frequency resolution was sufficient to [converge](https://meep.readthedocs.io/en/latest/FAQ/#checking-convergence) to the desired accuracy.  Technically, what we are doing is [spectral density estimation](https://en.wikipedia.org/wiki/Spectral_density_estimation) of the [periodogram](https://en.wikipedia.org/wiki/Periodogram) by an ensemble average with a [rectangular window](https://en.wikipedia.org/wiki/Window_function#Rectangular_window), but many variations on this general idea exist.  (3) The material properties for Ag are imported from the [materials library](https://meep.readthedocs.io/en/latest/Materials/#materials-library). (4) At the end of the run, the flux spectra from all iterations are saved to a file using NumPy's [uncompressed `.npz` format](https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html#module-numpy.lib.format). This data is used to plot the results in post processing. (5) For Method 1, independent random numbers can be used for the white-noise dipole emitters in the trial runs for the two cases of a flat and textured surface, since all that matters is that they average to the same power spectrum.\n",
    "\n",
    "Results for Methods 1 and 2 for the two cases of a flat and textured surface are generated using the following shell script:\n",
    "```\n",
    "#!/bin/bash\n",
    "\n",
    "# Method 1: flat surface\n",
    "python stochastic_emitter.py -method 1 -res 50 -nf 500 -nd 10 -nr 500\n",
    "\n",
    "# Method 1: textured surface\n",
    "python stochastic_emitter.py -method 1 -res 50 -nf 500 -nd 10 -nr 500 -textured\n",
    "\n",
    "# Method 2: flat surface\n",
    "python stochastic_emitter.py -method 2 -res 50 -nf 500 -nd 10\n",
    "\n",
    "# Method 2: textured surface\n",
    "python stochastic_emitter.py -method 2 -res 50 -nf 500 -nd 10 -textured\n",
    "```\n",
    "\n",
    "Afterwards, the four NumPy files produced by each run are used to plot the normalized flux for each method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "method1_f0 = np.load(\"method1_flat_res50_nfreq500_ndipole10.npz\")\n",
    "method1_f1 = np.load(\"method1_textured_res50_nfreq500_ndipole10.npz\")\n",
    "\n",
    "method1_freqs = method1_f0[\"freqs\"]\n",
    "method1_f0_mean = np.mean(method1_f0[\"fluxes\"], axis=1)\n",
    "method1_f1_mean = np.mean(method1_f1[\"fluxes\"], axis=1)\n",
    "\n",
    "method2_f0 = np.load(\"method2_flat_res50_nfreq500_ndipole10.npz\")\n",
    "method2_f1 = np.load(\"method2_textured_res50_nfreq500_ndipole10.npz\")\n",
    "\n",
    "method2_freqs = method2_f0[\"freqs\"]\n",
    "method2_f0_mean = np.mean(method2_f0[\"fluxes\"], axis=1)\n",
    "method2_f1_mean = np.mean(method2_f1[\"fluxes\"], axis=1)\n",
    "\n",
    "plt.semilogy(method1_freqs, method1_f1_mean / method1_f0_mean, \"b-\", label=\"Method 1\")\n",
    "plt.semilogy(method2_freqs, method2_f1_mean / method2_f0_mean, \"r-\", label=\"Method 2\")\n",
    "plt.xlabel(\"frequency\")\n",
    "plt.ylabel(\"normalized flux\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Results for Method 1 for three different numbers of trials/iterations (10, 50, and 500) are shown in the following three figures. Each trial/iteration involves two runs: one each for the flat and textured surface. As the number of trials/iterations is increased, the \"noisiness\" in the plot is gradually reduced. However, the total runtime increases significantly.\n",
    "\n",
    "![](https://meep.readthedocs.io/en/latest/images/stochastic_emitter_trials.png)\n",
    "\n",
    "The next figure shows a comparison of the normalized radiated flux for Method 1 (500 trials) and 2 (20 runs; 10 runs each for the flat and textured surface). The results show good agreement over the entire bandwidth spectrum. The Method 1 results required almost *four days* whereas the Method 2 results were obtained in less than forty minutes. In general, deterministic approaches tend to be more efficient than brute-force Monte-Carlo.\n",
    "\n",
    "![](https://meep.readthedocs.io/en/latest/images/stochastic_emitter_normalized_flux_comparison.png)\n",
    "\n",
    "One situation in which you may need to perform brute-force Monte-Carlo simulations is that of nonlinear or time-varying media, because the equivalence between random and deterministic simulations above relied on linearity and time-invariance.   However, in such media one also cannot directly employ white-noise sources, but you must instead input the noise with the correct spectrum for your desired emission process.   For example, to [model thermal emission in a nonlinear medium](http://doi.org/10.1103/PhysRevB.91.115406) one must have a noise spectrum consistent with the [fluctuation-dissipation theorem](https://en.wikipedia.org/wiki/Fluctuation-dissipation_theorem), which can be achieved using the `NoisyLorentzianSusceptibility` feature in Meep."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}