{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflectance Spectrum of Air-Silica Interface\n",
    "\n",
    "We will compute the broadband [reflectance](https://en.wikipedia.org/wiki/Reflectance) spectrum at normal incidence for a planar interface of vacuum/air and [fused quartz](https://en.wikipedia.org/wiki/Fused_quartz). As a validation, we will compare the simulated result with the analytic reflectance computed using the [Fresnel equations](https://en.wikipedia.org/wiki/Fresnel_equations). This is similar to an example in [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](https://meep.readthedocs.io/en/latest/Python_Tutorials/Basics/#angular-reflectance-spectrum-of-a-planar-interface) which involved computing the broadband reflectance spectrum over a range of angles.\n",
    "\n",
    "The wavelength-dependent, lossless permittivity of fused quartz, measured experimentally at 20°C, can be approximated by the [Sellmeier equation](https://en.wikipedia.org/wiki/Sellmeier_equation):\n",
    "\n",
    "$$\\varepsilon(\\lambda) = 1 + \\frac{0.6961663\\lambda^2}{\\lambda^2-0.0684043^2} + \\frac{0.4079426\\lambda^2}{\\lambda^2-0.1162414^2} + \\frac{0.8974794\\lambda^2}{\\lambda^2-9.896161^2}$$\n",
    "\n",
    "The wavelength λ is in units of microns. This equation is valid from 0.21 to 6.7 μm. The Sellmeier form for the permittivity of fused quartz can be imported into Meep as a [Lorentzian susceptibility](https://meep.readthedocs.io/en/latest/Materials/#material-dispersion) via a slight reorganization to convert the wavelength dependence into frequency as described in [Materials](https://meep.readthedocs.io/en/latest/Materials/#sellmeier-coefficients). This is implemented in the [materials library](https://github.com/NanoComp/meep/blob/master/python/materials.py#L153-L173).\n",
    "\n",
    "The simulation involves a 1d cell. A planewave current source with a pulsed profile spanning visible wavelengths of 0.4 to 0.8 μm is normally incident on the fused quartz from air. The reflectance is computed using the convention of two separate runs: (1) an empty cell to obtain the incident power, and (2) with the fused quartz to obtain the reflected power. The details of this type of calculation are described in [Tutorial/Basics/Transmittance Spectrum of a Waveguide Bend](https://meep.readthedocs.io/en/latest/Python_Tutorials/Basics/#transmittance-spectrum-of-a-waveguide-bend). The grid resolution, and by direct extension the time resolution via the [Courant condition](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition), must be made sufficiently fine to obtain agreement with the analytic results and to ensure [numerical stability](https://meep.readthedocs.io/en/latest/Materials/#numerical-stability). Coarse resolutions may lead to field instabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000452042 s\n",
      "Working in 1D dimensions.\n",
      "Computational cell is 0 x 0 x 12 with resolution 200\n",
      "time for set_epsilon = 0.00140619 s\n",
      "-----------\n",
      "field decay(t = 50.0025): 0.25018932773921454 / 0.25018932773921454 = 1.0\n",
      "on time step 32117 (time=80.2925), 0.000124548 s/step\n",
      "field decay(t = 100.0025): 4.358317201613301e-16 / 0.25018932773921454 = 1.742007639173244e-15\n",
      "run 0 finished at t = 100.0025 (40001 timesteps)\n",
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.00031209 s\n",
      "Working in 1D dimensions.\n",
      "Computational cell is 0 x 0 x 12 with resolution 200\n",
      "     block, center = (0,0,3)\n",
      "          size (1e+20,1e+20,6)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (1,1,1)\n",
      "time for set_epsilon = 0.00108504 s\n",
      "lorentzian susceptibility: frequency=0.101049, gamma=0\n",
      "lorentzian susceptibility: frequency=8.60279, gamma=0\n",
      "lorentzian susceptibility: frequency=14.619, gamma=0\n",
      "-----------\n",
      "field decay(t = 50.0025): 0.16530523240803463 / 0.16530523240803463 = 1.0\n",
      "on time step 24293 (time=60.7325), 0.000164657 s/step\n",
      "field decay(t = 100.0025): 1.4363078056736928e-16 / 0.16530523240803463 = 8.68882239691211e-16\n",
      "run 0 finished at t = 100.0025 (40001 timesteps)\n"
     ]
    }
   ],
   "source": [
    "import meep as mp\n",
    "from meep.materials import fused_quartz\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "resolution = 200  # pixels/μm\n",
    "\n",
    "dpml = 1.0\n",
    "sz = 10+2*dpml\n",
    "cell_size = mp.Vector3(z=sz)\n",
    "pml_layers = [mp.PML(dpml)]\n",
    "\n",
    "wvl_min = 0.4\n",
    "wvl_max = 0.8\n",
    "fmin = 1/wvl_max\n",
    "fmax = 1/wvl_min\n",
    "fcen = 0.5*(fmax+fmin)\n",
    "df = fmax-fmin\n",
    "nfreq = 50\n",
    "\n",
    "sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df), component=mp.Ex, center=mp.Vector3(z=-0.5*sz+dpml))]\n",
    "\n",
    "sim = mp.Simulation(cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    sources=sources,\n",
    "                    dimensions=1,\n",
    "                    resolution=resolution)\n",
    "\n",
    "refl_fr = mp.FluxRegion(center=mp.Vector3(z=-0.25*sz))\n",
    "refl = sim.add_flux(fcen, df, nfreq, refl_fr)\n",
    "\n",
    "sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(), 1e-9))\n",
    "\n",
    "empty_flux = mp.get_fluxes(refl)\n",
    "empty_data = sim.get_flux_data(refl)\n",
    "sim.reset_meep()\n",
    "\n",
    "geometry = [mp.Block(mp.Vector3(mp.inf,mp.inf,0.5*sz), center=mp.Vector3(z=0.25*sz), material=fused_quartz)]\n",
    "\n",
    "sim = mp.Simulation(cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    geometry=geometry,\n",
    "                    sources=sources,\n",
    "                    dimensions=1,\n",
    "                    resolution=resolution)\n",
    "\n",
    "refl = sim.add_flux(fcen, df, nfreq, refl_fr)\n",
    "sim.load_minus_flux_data(refl, empty_data)\n",
    "\n",
    "sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mp.Vector3(), 1e-9))\n",
    "\n",
    "refl_flux = mp.get_fluxes(refl)\n",
    "R_meep = -1*np.divide(refl_flux,empty_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A plot of the reflectance spectrum based on the simulated data and the analytic Fresnel equation is generated and shown below. There is agreement between the simulated and analytic results. Note that the reflectance spectra is plotted as a function of wavelength, not frequency from which the Meep data is obtained. Thus, the data points are not equally spaced: the spacing is smaller at low wavelengths (high frequencies) than at high wavelengths (low frequencies)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEKCAYAAADJvIhZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5RU1Zn38e/TLVdBEW27iSC0Dl4QFbkTvOAlgk4S7xOdznjJqCFGJxOTGGbQiU7GWYmaNRl9DQYNRhMyxKhRJCpiIqiJiKCC4hURsBUa5KZyE+jn/WOfpovqupwuqrqru36ftWp116l9Tp1Tq6mHvfdznm3ujoiISEsoa+0TEBGR0qGgIyIiLUZBR0REWoyCjoiItBgFHRERaTEKOiIi0mIKGnTMbJyZvW1mS8xsQorXzcxuj15fZGaDo+2dzWyemS00s8VmdlPSftdEx11sZrdE2/qZ2RYzezV63FXIaxMRkebbq1AHNrNy4E7gS0At8JKZTXf3NxKanQH0jx4jgEnRz23AKe7+mZl1AJ43syfcfa6ZnQycBRzj7tvM7MCE473n7oMKdU0iIrJnCtnTGQ4scfel7v45MI0QLBKdBdzvwVygh5n1ip5/FrXpED0a7mL9FvATd98G4O6rC3gNIiKSRwXr6QAHAR8kPK8l9GKytTkIWBn1lBYAfwfc6e4vRm0OA04ws5uBrcD33f2l6LVqM3sF+AS43t2fSz4pM7sSuBJg7733HnLEEUfswSWKiJSeBQsWfOzuFbnsW8igYym2JdfcSdvG3XcCg8ysB/BHMxvo7q8Tznk/YCQwDHjAzA4BVgIHu/taMxsCPGJmR7n7J7sd3H0yMBlg6NChPn/+/NyvUESkBJnZ8lz3LeTwWi3QJ+F5b+Cj5rZx9w3AbGBcwj4PR0Nw84B64AB33+bua6N9FgDvEXpFIiJSJAoZdF4C+ptZtZl1BC4Epie1mQ5cHGWxjQQ2uvtKM6uIejiYWRfgNOCtaJ9HgFOi1w4DOgIfR/uUR9sPISQnLC3g9YmISDMVbHjN3XeY2dXATKAcmOLui81sfPT6XcDjwJnAEmAzcFm0ey/gviiIlAEPuPuM6LUpwBQzex34HLjE3d3MTgT+08x2ADuB8e6+rlDXJyIizWelvLSB5nREStP27dupra1l69atrX0qRa1z58707t2bDh067LbdzBa4+9BcjlnIRAIRkaJUW1tL9+7d6devH2ap8pnE3Vm7di21tbVUV1fn7bglXQZnwQLo1w+mTm3tMxGRlrR161b2339/BZwMzIz9998/773Bkg46AMuXw5VXKvCIlBoFnOwK8RmVfNAB2LwZJk5s7bMQEWn/FHQiK1a09hmIiLR/CjqRgw9u7TMQkWI1dWqY/y0r0zzwnlLQAczg00/1ByUiTU2dGuZ9ly8H9/zMAy9btowjjjiCyy+/nIEDB1JTU8PTTz/N6NGj6d+/P/PmzWPTpk184xvfYNiwYRx33HE8+uijAOzcuZMf/OAHDBs2jGOOOYZf/vKXAMyePZsTTzyRc845hwEDBjB+/Hjq6+vz8RHkVcmnTHfrBp99Buui20gb/qAAampa77xEpGX867/Cq6+mf33uXNi2bfdtmzfDP/8z3H136n0GDYKf/zzz+y5ZsoQ//OEPTJ48mWHDhvG73/2O559/nunTp/Pf//3fDBgwgFNOOYUpU6awYcMGhg8fzmmnncbUqVPZd999eemll9i2bRujR4/m9NNPB2DevHm88cYb9O3bl3HjxvHwww9z/vnnN+PTKLyS7ukMGQI9ezbdrsQCEWmQHHCybY+rurqao48+mrKyMo466ihOPfVUzIyjjz6aZcuW8dRTT/GTn/yEQYMGMWbMGLZu3cqKFSt46qmnuP/++xk0aBAjRoxg7dq1vPvuuwAMHz6cQw45hPLyci666CKef/75PTvJAij5ns4HH6TersQCkdKQrUfSr18YAUnWty/Mnp37+3bq1GnX72VlZbuel5WVsWPHDsrLy3nooYc4/PDDd9vP3bnjjjsYO3bsbttnz57dJMW5GNPCS7qnA+kTCJRYICIAN98MXbvuvq1r17C9kMaOHcsdd9xBQ6myV155Zdf2SZMmsX37dgDeeecdNm3aBIThtffff5/6+np+//vfc/zxxxf2JHNQ8kEn1R8UhP/ZKKlARGpqYPLk0LMxCz8nTy78nO8NN9zA9u3bOeaYYxg4cCA33HADAJdffjkDBgxg8ODBDBw4kG9+85vs2LEDgFGjRjFhwgQGDhxIdXU155xzTmFPMgcq+Dl/PlOnhjmcVF3orl1b5g9MRFrOm2++yZFHHtnap5FXs2fP5rbbbmPGjBnZGzdDqs9qTwp+lnxPB0JAWbYs/A8mmZIKRETyp+QTCRKlSx5QUoGIFLsxY8YwZsyY1j6NrNTTSaCkAhGRwirtoLNgQZgZNIOqqpRJBapWICKSP6UddBLV1TXJUunWLZS9WLcuf+UvRERKmYJOkoakgvp6VSsQEck3BZ0MVK1ARIrZr3/9a66++uqsbT766KNdzy+//HLeeOONQp9aWgo6iRrmd6I5HiUWiAhVVbt/NyR8R7QFyUHnnnvuYcCAAa12Pgo66dTVpUwsKC8PVamVWCBSIurqmrc9prPPPpshQ4Zw1FFHMXnyZAC6devGxIkTOfbYYxk5ciR10Xs89thjjBgxguOOO47TTjtt1/YGn376KdXV1btK43zyySf069ePP/zhD8yfP5+amhoGDRrEli1bGDNmDPPnzwfgySefZPDgwRx77LGceuqpe3Q9cek+nQwaqhBMnBiG1PbaC7Zvh7Vrw3YtgyDSDmRb2yCTdPfFxFjbYMqUKfTs2ZMtW7YwbNgwzjvvPDZt2sTIkSO5+eabue6667j77ru5/vrrOf7445k7dy5mxj333MMtt9zCz372s13H6t69O2PGjOFPf/oTZ599NtOmTeO8887jggsu4M477+S2225j6NDdCwisWbOGK664gmeffZbq6mrWNazvUmCl3dMZMiRrk8TEgsrKpq8rsUBEcnH77bfv6tF88MEHvPvuu3Ts2JEvf/nLAAwZMoRly5YBUFtby9ixYzn66KO59dZbWbx4cZPjXX755dx7770A3HvvvVx22WUZ33/u3LmceOKJVFdXA9AzVeZUAain0wwffph6uxILRNqwbGsbZFoeIMe1DWbPns3TTz/NCy+8QNeuXXetl9OhQ4ddyxGUl5fvKuR5zTXXcO211/LVr36V2bNnc+ONNzY55ujRo1m2bBlz5sxh586dDBw4MOM5uHurLH1Q2j2dbJImDJVYICL5sHHjRvbbbz+6du3KW2+9xdy5c7O2P+iggwC477770ra7+OKLueiii3br5XTv3p1PP/20SdtRo0YxZ84c3n//fQANr7WYVGNmyaJJOy2DIFKC0n1HxPnuSGPcuHHs2LGDY445hhtuuIGRI0dmbH/jjTdywQUXcMIJJ3DAAQekbVdTU8P69eu56KKLdm279NJLGT9+/K5EggYVFRVMnjyZc889l2OPPZavfe1rOV9Pc2hpgyiLA8jcjY4+Jy2DINL2tcelDQAefPBBHn30UX7zm9/k7Zha2qCVaRkEESlG11xzDRMmTNi12FuxUiJBjrQMgogUkzvuuKO1TyEW9XTiUrUCkXallKcW4irEZ6SgkyjuxGCaagUAGzaoWoFIsevcuTNr165V4MnA3Vm7di2dO3fO63E1vJZo1ardn2dILEiuVtCjB6xfDxs3hu2qViBSvHr37k1tbS1r1qxp7VMpap07d6Z37955Paay1xKz15LFyGZr0K9f6oy2vn1D4oGISHuh7LXWkDTHo8QCEZHsChp0zGycmb1tZkvMbEKK183Mbo9eX2Rmg6Ptnc1snpktNLPFZnZT0n7XRMddbGa3JGz/t+hYb5vZ2EJe227q6pRYICISQ8GCjpmVA3cCZwADgIvMLHkRhzOA/tHjSmBStH0bcIq7HwsMAsaZ2cjouCcDZwHHuPtRwG3R9gHAhcBRwDjgF9E55K4ZdxynSyxYv16JBSIiDQrZ0xkOLHH3pe7+OTCNECwSnQXc78FcoIeZ9Yqefxa16RA9GiZRvgX8xN23Abj76oRjTXP3be7+PrAkOofcrVoV5m4aHhnU1IRqBH37hhG3hoKtn3wSdm1ILFDgEZFSVsigcxCQuOBzbbQtVhszKzezV4HVwCx3fzFqcxhwgpm9aGZzzGxYM94PM7vSzOab2fx8Z64kLoPQvXvT11WxQERKXSGDTqrUr+TuQto27r7T3QcBvYHhZtZQp3svYD9gJPAD4AEL9bnjvB/uPtndh7r70IqKinhXEkdSRWolFoiINFXIoFML9El43hv4qLlt3H0DMJswT9Owz8PRENw8oB44IOb77ZlmVKROl0BQVqY5HhEpXYUMOi8B/c2s2sw6Eib5pye1mQ5cHGWxjQQ2uvtKM6swsx4AZtYFOA14K9rnEeCU6LXDgI7Ax9GxLjSzTmZWTUhOmJfXK0qc48kiXWLBzp2a4xGR0lWwoOPuO4CrgZnAm8AD7r7YzMab2fio2ePAUsKk/93AVdH2XsAzZraIELxmufuM6LUpwCFm9johOeGSqNezGHgAeAN4Evi2u+8s1PVlk5xYUJ4ij05zPCJSalSRIFNFgkziLvNaWQmrVlFWlrqDZBYSD0RE2gpVJChmWeZ4dPOoiJQSBZ1cNXOp2lRzPGbhPh4lFohIqVDQyVUzbhyFpnM8++wTdlu/XokFIlI6FHRaUOLNo/vt1/R1JRaISHunoNMSkipSg24eFZHSpKCTL81YdRSUWCAipUlBJ1+aOceT7ubRjz9WYoGItF8KOq0kObFg//3D9k2blFggIu2Xgk5riOZ3ar5uLNtaRX09dOvWtJkSC0SkvVHQaW3RHI8SC0SkFCjoFEozbx5Nl0Cw775hfkfzPCLSHijoFEqeEgs2bAjzO5rnEZH2QEGnSCQnFvTtCz16NG2neR4Racv2au0TEHZVrK4BaioroX4VEIbUUtE8j4i0VerptJRm3jwKuoFURNofBZ2W0sw5Hkg/z7NxoxILRKRtUtApYsnzPPvuG7Zv2KDEAhFpmxR0ilFCgdCa71XtqkytxAIRaesUdIpdwhyPbiAVkbZOQae1NPPmUUifQNCjh24gFZG2QUGnteQxsWD9et1AKiJtg4JOW5BQIHTTZmN1WdWuG0j32adpc83ziEix0s2hbVBFfR31UedIN5CKSFuink6xyGGOB9LP85SVaY5HRIqPgk6xSJ7j2cN5np07NccjIsVHQaetSpjn+bhD1a4bSMvLmzbVHI+IFAvN6bQDXTbWsWxD+F1zPCJSzNTTKWZ5vJdHRUJFpBgo6BSzPN7LU1kZUqyVXCAirUnDa+1MTU34OXFiGFLr0wc6dYJ58xrbNCQXJLYXEWkJ6um0F4lFQr9uLNtaRX19CDDbtjVtruQCEWkNCjptSXPmeBIKhX7wQeomSi4QkZamoNOW5HgvjwqFikixUNApASoUKiLFIlbQMbNKM/uVmT0RPR9gZv8cY79xZva2mS0xswkpXjczuz16fZGZDY62dzazeWa20MwWm9lNCfvcaGYfmtmr0ePMaHs/M9uSsP2uuB9Cu5XmBtK+fRtXIU2keR4RKbS42Wu/Bu4FGr6S3gF+D/wq3Q5mVg7cCXwJqAVeMrPp7v5GQrMzgP7RYwQwKfq5DTjF3T8zsw7A82b2hLvPjfb7H3e/LcXbvufug2JeU0lJvIEUdBOpiLSOuMNrB7j7A0A9gLvvAHZm2Wc4sMTdl7r758A04KykNmcB93swF+hhZr2i559FbTpEj3gTGKUmz4VCe/bUPI+IFE7coLPJzPYn+uI3s5HAxiz7HAQk5k3VRttitTGzcjN7FVgNzHL3FxPaXR0Nx00xs/0Stleb2StmNsfMTkh1UmZ2pZnNN7P5a9asyXIJbUAON5BC+nmedes0zyMihRM36FwLTAcONbO/AvcD12TZx1JsS/5WTNvG3XdGQ2W9geFmNjB6fRJwKDAIWAn8LNq+EjjY3Y+Lzvd3ZtZkiTN3n+zuQ919aEVFRZZLaGeS7uVJnufp2rVp3NI8j4jkU6w5HXd/2cxOAg4nBIq33X17lt1qgT4Jz3sDHzW3jbtvMLPZwDjgdXffdQOKmd0NzIjabSPMBeHuC8zsPeAwYH6cayxFmucRkZYWN3vt20A3d1/s7q8D3czsqiy7vQT0N7NqM+sIXEjoLSWaDlwcZbGNBDa6+0ozqzCzHtF7dwFOA96KnvdK2P8c4PVoe0WUvICZHUJITlga5/ralRzneEDzPCJSeHGH165w913/J3b39cAVmXaIkg2uBmYCbwIPuPtiMxtvZuOjZo8TAsMS4G6gIZD1Ap4xs0WE4DXL3WdEr91iZq9Fr50MfDfafiKwyMwWAg8C4919Xczraz9yvIEUNM8jIoVnHuNLKfqCP9ajxlGPYpG7H1Xg8yuooUOH+vz5JTD6ZqmmzlLbsm8lR/ZYtatYaF1d6tptffvCsmX5O0URaTvMbIG7D81l37j36cwEHohuuHRgPPBkLm8oxU3zPCJSSHGH134I/AX4FvBt4M/AdYU6KckzzfOISJGINbzWXpXM8FoqmYbcEv4mpk4NczibNzfdPfFPp2tXmDxZ6/OIlII9GV6Lm7022sxmmdk7ZrbUzN43s9LLDCsVGe7nOfjgsCic7ucRkVzEndP5FSFLbAHZy99IO6N5HhHJl7hzOhvd/Ql3X+3uaxseBT0zKawCzPP06ZN6u4hIg7hB5xkzu9XMRpnZ4IZHQc9MCqsA9/NACDxKLhCRdOIOr42IfiZOHDlwSn5PR4pWQuJBDXDuvpUc2bfxfp6+feG55xqbN9xECkouEJFGcWuvnVzoE5G2JXmep1+/pm0akgsUdESkQdyeDmb298BRQOeGbe7+n4U4KWkllZWhBEEO0iURKLlARBLFTZm+C/gaYTkDAy4A+hbwvKQ1NHeeJyG1ut6NlVSlbDJhgm4kFZEgbiLBF939YmC9u98EjGL3JQlEqGL3XlLnzrDPPvDTn6pgqIgEcYPOlujnZjP7ArAdqC7MKUlblrgo3D33QPfuTdvoRlKR0hU36MyI1re5FXgZWAZMK9RJSRFp5v08y5ZBfX34WVMDtbWp2y1friE3kVIUN5HglmhlzofMbAYhmWBr4U5LisaqVU23ZarblvRaXVklB3qKYxACT8NPpVeLlIa4PZ0XGn5x923uvjFxm0g6FfV1aW8kTaQhN5HSkDHomFmVmQ0BupjZcQnVCMYAMb5KpF1q5pDb5Mm7z/Wko/RqkfYvW09nLHAb0Bv4WcLju8C/F/bUpGg1M7W6pmb3uZ50gaeiIsztaK5HpP2Ku1z1ee7+UAucT4sq6fV08q0ZS2ID1FFJFY1zPQ3r8+y1F+zY0dhO6/SIFJ+Cr6cDDImy1xrecD8z+69c3lAEoJK6lOnViQEHNNcj0t7EDTpnuPuuSlvuvh44szCnJG1SDkslJA65feMb8Nlnqdtprkek/YgbdMrNrFPDEzPrAnTK0F5KTap5nmYuhZ5unR7N9Yi0H3Hv0/kt8Gczu5ewpME3gPsKdlZSGpLmgd7ct5IDuq5i8+bdm6xeDZdcAjujNWt1X49I2xWrp+PutwD/BRxJqDT942ibSN502VjXJL36V78K9dt2Ji2SrrkekbYpVvYagJn1Bfq7+9Nm1hUod/dPC3p2BabstRZQVdW85RJS/D2WlaUeqTMLc0Ii0rIKnr1mZlcADwK/jDYdBDySyxtKidmD5RIaHnXWdMkEgB494De/0VyPSFsSd07n28Bw4EUAd3/XzA4s2FmJJGgopZM411NWBuvXw6WXNvZ2NNcjUvziZq9tc/fPG56Y2V6EhAKRFpE813P//XDAAU2H1zTXI1Lc4gadOWb274QabF8C/gA8VrjTknYth3t6kkvp1NTA2rWp2+q+HpHiFTfoTADWAK8B3wQeB64v1ElJO5fLPT3NmOvp1AluvVVzPSLFKHb2Wnuk7LUi08z6bQB7d/Xd5no6dAg/t2/fvZ1quInkT8Gy18zsNTNblO6R2+mKpJHDsFvyXM+996Y+jOZ6RIpDxp6OmR0BbEn3ursvL8RJtRT1dNqITD0g3dcj0uIKeZ/O76LA8l/uvjz5kcsbiuRVM+Z63OGKK2DSJM33iLSWbPfpdDSzS4Avmtm5yS+6+8OFOS2R3KW6r6dLFzjppFBW5557Grfr3h6RlpWtpzMeGAn0AL6S9PhytoOb2Tgze9vMlpjZhBSvm5ndHr2+yMwGR9s7m9k8M1toZovN7KaEfW40sw/N7NXocWbCa/8WHettMxsb5wOQNiAPcz133w1PPAG9ejVtq/kekRbk7lkfwD/HaZe0TznwHnAI0BFYCAxIanMm8ARghOD2YrTdgG7R7x0IlRBGRs9vBL6f4v0GRO/RCaiO3rs80zkOGTLEpQ1LnXid/lFZ6WapXzJz/+1v3fv2Db/37Ruei0hTwHxvZkxoeMS9T2eamV1vZpMBzKy/mWXr6QwHlrj7Ug/VDKYBZyW1OQu4P7qOuUAPM+sVPW9Y0qtD9MiW230WMM3dt7n7+8CS6BxEgrq6tGv2uIflE5YvD783DLtpvkckv+IGnSnA58AXo+e1hKUOMjkI+CDheW20LVYbMys3s1eB1cAsd38xod3V0XDcFDPbrxnvh5ldaWbzzWz+mjVrslyCtDc33xzu2UnUpQt07qzlE0RaQtygc6iH9XO2A7j7FsIQWCapXk/uraRt4+473X0Q0BsYbmYDo9cnAYcCg4CVwM+a8X64+2R3H+ruQysqKrJcghS1HMvppJrv2bYtdXuV1BHJr7hB5/Noieow4WJ2KJDmn+kutUCfhOe9gY+a28bdNwCzgXHR87ooINUDd9M4hBbn/aQ9yaWcTo8e1HzdWLbcqPfws+br6dOsy8pCxpuWUBDJj6xBx8wMuAt4EuhjZlOBPwPXZdn1JaC/mVWbWUfgQmB6UpvpwMVRFttIYKO7rzSzCjPrEb1/F+A04K3oeWL+0TnA6wnHutDMOplZNdAfmJft+qTEbNyYcnNDmnWiTp1CgLn8cs33iORL1qATZSp8BzgXuBT4P2Cou8/Ost8O4GpgJvAm8IC7Lzaz8WY2Pmr2OLCUMOl/N3BVtL0X8ExUauclwpzOjOi1WxrK8wAnA9+N3m8x8ADwBiFAftvdk0bppSSkG3bLMhyXaqnsd9+FioqmHSjN94jkJlbBTzO7E/i1u79U+FNqOSqDU4KaW1S0spKy1asyltWZOjUEoBUr4OCDQ7KCbjSV9mxPyuDEXTn0ZOCbZrYc2ESYtHd3PyaXNxVpM+rqOLhvGFJL5g5HHw1LlsDWrWGbKhyIZBY3keAMQsbYKTRWI/hKoU5KpJikS7P+x3+EN95oDDgNNPQmkl6soOMpin26Cn5KW5THNOupU9MnyynVWiS1uD0dkfYhlzTrbt1SpllTVZW2wkGnTjBvXghMSrUWaRR3TkekdG3alHp7XR03/zbM4SSvXlpeDiNGhGDTsIaP5ntE1NMRCXJMs0419HbvvbByJfTo0XTROM33SKlT0BGB9MNuq1Zl3u/AA1MOvXXvX5XuPlRWrGgMRhp+k1Kj4TWRPZGuaGyWVOtjjoETT4T77mscmtPwm5QC9XRECiRVqnXXrvCtb4XfJ03afS4INPwm7Z+Cjkg2uc73fK+KTZsNp/GxabPxi4erWLQofXEEpVtLe6agI5JNrvM9dXVpt5eVkTbdeu+9YelSzfdI+6Q5HZFWcvPNTdOty8vD80MPDb83LCyn+R5pL9TTEdkTOQ69Qep06/vugw8+gH320Uqm0j4p6IjsiVyH3v7lX6CiImW69RcGV/Hpp6l3W768cUkgDb9JW6ThNZHWcNddsH176tcypFsD9O4No0bBc8+purW0PerpiBRKpqG3Zcsy7pou3frHP4ZzzoFZs1TdWtqmWIu4tVdaxE1aVaYF5SoqUt94WlkJq1ZRVpa6TmnDwnIihbQni7ippyNSjDJUOoD06dYAV10FixdrzkeKk4KOSBuUavitc2c4/niYMgUGDoSLLw5zPe6Ncz4KPNLaFHREWkuu6dbvvZey2sGWrcaz71RRW6sK11K8FHREWkuu6db9+2esdnDAAaStcL18OcyY0XgPkIbgpKUp6Ii0NT/6UdYm6eZ8ysrgK1+B6mo47zy44goNwUnLUtARKUaZht6yBZ2dO3lzQ9VuQ28Nj8+6VfHQQ3DEEfDww7Bly+67aghOCk0p00qZlrYoU7p1nz6hlk460b95pV1LrpQyLSKNjjwyVrN0Q3DuMGYM/OY3YeltzflIPinoiLRFmYbfZs6MdYhUadddusAFF8CHH4aU6298Q3M+kl8KOiJtUa6ZbwAjRsAvf0nNtZVN0q43bzEeeLaKd95JHdc2b4Yf/jD/lyOlQ0FHpNRs3gzjx8Pq1alfr6vDLP3LH34IY8eGHs/mzUq7luZRlWmR9qiyMvW9PJWVsGgRLFgAw4ZlPMTBB6eudL3vvvD22/D1r0OnTrBjhxabk/jU0xFpjzINv5nB0CyJR7ffzlvrDkyZdr2SKpYuhdmzYa+9tNicNI+Cjog09Z3v0PnT1EVHu2yso6wMTjpp96W2Ey1fDrfdFobiGmgYTkBBR0RSee21WM3SpV137Ag/+EG4ZeiUU0LlA1U/EFDQESldmdKuBw7MvO93vgMvvpi28sGGLiED7j/+I9ynes89qn4ggSoSqCKBSGqZqh506gTbtmXeP/pucYfy8vTVD7ZsCYeD0POZOBFWrAi9qJtvVkJCMSraigRmNs7M3jazJWY2IcXrZma3R68vMrPB0fbOZjbPzBaa2WIzuynFvt83MzezA6Ln/cxsi5m9Gj3uKuS1iZS0urpQriAGs8zVDw48EC65BK67Lgy5aQiufStY0DGzcuBO4AxgAHCRmQ1IanYG0D96XAlMirZvA05x92OBQcA4MxuZcOw+wJeAFUnHe8/dB0WP8fm+JpGSkmn4bd994dJLM+//05/C0qUAaYfhPulaxbnnwqOPwq23Nk1M0BBc+1PIns5wYIm7L3X3z4FpwNEPcWEAABA0SURBVFlJbc4C7vdgLtDDzHpFzz+L2nSIHomd8/8BrkvaJiL5tCdVDwAmTIBDD4Vhw+iyMfX6P90313HvvaHjlG40b8UK+PzzHK9Bik4hg85BQGKp29poW6w2ZlZuZq8Cq4FZ7v5itP2rwIfuvjDFe1ab2StmNsfMTkh1UmZ2pZnNN7P5a9KtQy8ie+7990P3JdPcUKRTp+xDcP/0T6FHtGWL0q/bskIGnVR/ack9k7Rt3H2nuw8CegPDzWygmXUFJgL/kWK/lcDB7n4ccC3wOzPbp8nB3Se7+1B3H1pRUdGMyxGRJjINwfXrB9//Psybl/kYixeDe9ohuI3RENyf/gRnnx2W4r7kEs39tFWFDDq1QJ+E572Bj5rbxt03ALOBccChQDWw0MyWRe1fNrMqd9/m7mujfRYA7wGH5etiRCSFPR2Cg5CefeSRaYfg9tlcx5QpYQjuqafCPUCpqiBMaJKqJMWokEHnJaC/mVWbWUfgQmB6UpvpwMVRFttIYKO7rzSzCjPrAWBmXYDTgLfc/TV3P9Dd+7l7P0LQGuzuq6J9yqN9DiEkJywt4PWJSD784hfQu3fWZh06wJe+BJs2pX69tjasA/Tzn4eRPdAwXDEqWNBx9x3A1cBM4E3gAXdfbGbjzawhs+xxQmBYAtwNXBVt7wU8Y2aLCMFrlrvPyPKWJwKLzGwh8CAw3t3X5fWiRCQ3mYbhvvUtePrpzPs/8siuaFNnqYfhVlkV69bBd78LhxwS5oguu0zDcMVGN4fq5lCR4pAt4aBz59DVeeyx9G3cee+9EKMmTkx9/2q66tkSX9HeHCoikhd//nPopixalLXpoYfC976XPs16xYqwOur998PHHzdu11Bcy1BPRz0dkeJQVZV+DaCGxAT3EBXSeeopOPFE6NSJNeVVVNQ3PV4dlRzXaxUrV4ZDjRoVppQefRS2bm1s17UrTJ6sMjyp7ElPR0FHQUekbck2DNetW1ja9KGH0jap3+m88koYqXvsMXj55dTt+vaFZctyP9X2SsNrIiIQIkhNDcydm7FZWRkMGQI33hgWUU0Xx5YvD8l1DdlwDTQUlzstVy0ibUumpbi//OXwyDYMN20anH469OwJhIy4Cm96zFVU0uvbYWjviCPgzDND9YT//d/GOnFaort5NLym4TWR9inbMFxZGYwYAWecERb+SeOdt50nnoDHHw9LdKdLUCiloTgNr4mINMcLL8D118OOHfCjH2VsethhYc26mTNh3brMQ3EzZsBnnzV9TcNxjdTTUU9HpH2Kkw0HsHp1+ptXAZ55Br74xVB/B9Jmxa2ikl6sokMHOP74kMswdiy8/jp885u7L9vQ1jPjlL2WIwUdEQGyD8XtvXeosXP66aHbk8bTs5yZM0Ov6LXXwrayMqivb9q2LQ/HKejkSEFHRIDMQefRR0MUeeopWLIk83ESvk8/+ghmzUq/1p0ZbNwI3bs3/3Rbm+Z0RET2RKbacF/9Ktx5J7z7Lrz3XubjPPnkrhpxX/hCWIJhdVnqWnEfeRU9e8IJJ8BNN8Ff/wrbt7f/+R/1dNTTEZHmyDYU16FDKHNw6qlw2mkwenTapv82wXn66XCvkHtIx96xY/elG4px/kfDazlS0BGRZssUdGbODHXi/vznUOYg2/dr9Pq6dSFf4dJLU2e/9ewZEu7694+1EGvBaXhNRKSlZBqKO/10+OlPYf78UE30wQczH+v118Gdnj3hvPPg3c9SD8UtXlfF4YdDnz5h2e4pU3ZPQmhLQ3KqSCAi0hxxV0VtiCSZHH00VFTASSfBySdTRerVU6uoY9Kk0BuaORN++9uwvV+/sFTD3LmNN60We4UEDa9peE1ECinTeNi994ZI8swz8MEHmY8TfVe7wxtvNO72yCOpU7IPOiisploIGl4TESlWmYbjLr0U7rsvdE+yZcbdfjssXIjV7+Soo+Dqq0Mh7Q/rUw/Jzf+wir594eKL4Ve/Csl3iX2M1hqSU09HPR0RKRZxsgR69AglD048MTxGjkzb9ILznTlzYM2a8LxXrzCS17lzqHma6/pByl7LkYKOiBSVTEFn+XJ47jmYMweefRbefjv78dxxh7featxtzpxw42oqcZfyVtDJkYKOiBSVuPXiILR77rmw9nY6f/xjuE+oomLXJneoK6tKmbSwikpqTlnFCSeEztTIkWFNvAZTp8LEibB8+VDc5+eUvK2go6AjIm1ZnCG5ww8PUWT06PDzsMPSNh18nLNwYUhOKC+HQYNC1YT6erj7btiyBSD3oKOUaRGR9uqvf4Xnnw+PP/4xZBRk8fLLoSbc3Llht+eeg7vu2n3+Z0+op6Oejoi0ZXGH5OrrwzzQ88833siTysMPhzI+VVW7Nn3+Oazr1DgkNxSY767hteZS0BGRkhRnSK66OqwjNGpU+Dl48K6X9iToaHhNREQa/e1vodDb3/4Gf/lL3m/g0c2hIiKlJtMNq6NGwbXXhrpxH34Yirz97nd5e2sFHRGRUrNqVcidTn4kp2WbhSVOL7oob2+toCMiIi1GQUdERLJLNyTXTEokEBGR7BKG3haYLcj1MOrpiIhIi1HQERGRFqOgIyIiLUZBR0REWoyCjoiItJiCBh0zG2dmb5vZEjObkOJ1M7Pbo9cXmdngaHtnM5tnZgvNbLGZ3ZRi3++bmZvZAQnb/i061ttmNraQ1yYiIs1XsJRpMysH7gS+BNQCL5nZdHd/I6HZGUD/6DECmBT93Aac4u6fmVkH4Hkze8Ld50bH7hMdd0XC+w0ALgSOAr4APG1mh7n7zkJdo4iINE8hezrDgSXuvtTdPwemAWcltTkLuN+DuUAPM+sVPf8satMheiSWw/4f4LqkbWcB09x9m7u/DyyJzkFERIpEIW8OPQj4IOF5LaEXk63NQcDKqKe0APg74E53fxHAzL4KfOjuC2338twHAXNTHGs3ZnYl0LCYxDYze72Z1yXpHQB83Non0Y7o88wffZb5dXiuOxYy6KRaayF58Z60baJhsUFm1gP4o5kNBJYCE4HTc3w/3H0yMBnAzOa7+9C0VyDNos8zv/R55o8+y/wys5wXIivk8Fot0CfheW/go+a2cfcNwGxgHHAoUA0sNLNlUfuXzawq5vuJiEgrKmTQeQnob2bVZtaRMMk/PanNdODiKIttJLDR3VeaWUXUw8HMugCnAW+5+2vufqC793P3foRAM9jdV0XHutDMOplZNSE5YV4Br09ERJqpYMNr7r7DzK4GZgLlwBR3X2xm46PX7wIeB84kTPpvBi6Ldu8F3BfN65QBD7j7jCzvt9jMHgDeAHYA346RuTY5t6uTNPR55pc+z/zRZ5lfOX+e5t5k2kNERKQgVJFARERajIKOiIi0mJIIOtnK8SS0G2ZmO83s/JY8v7YmRnmjMWa20cxejR7/0Rrn2RbE+duMPs9Xo5JQc1r6HNuSGH+bP0j4u3w9+vfeszXOtS2I8Xnua2aPJZQsuyzVcXbj7u36QUhieA84BOgILAQGpGn3F0Jyw/mtfd7F+ojzeQJjgBmtfa7F/oj5WfYgJMccHD0/sLXPu1gfcf+tJ7T/CvCX1j7vYn3E/Pv8d+Cn0e8VwDqgY6bjlkJPJ045HoBrgIeA1S15cm1Q3M9TsovzWf4j8LC7rwBwd/19ptfcv82LgP9rkTNrm+J8ng50t1Aephsh6OzIdNBSCDrpSu3sYmYHAecAd7XgebVVWT/PyKioy/2EmR3VMqfW5sT5LA8D9jOz2Wa2wMwubrGza3vi/m1iZl0JN5w/1ALn1VbF+Tz/H3Ak4Ub814DvuHt9poMWsgxOsYhTHufnwA/dfWdSPTdpKs7n+TLQ10OV8DOBRwg368ru4nyWewFDgFOBLsALZjbX3d8p9Mm1QbFKYUW+AvzV3dcV8Hzaujif51jgVeAUQsWYWWb2nLt/ku6gpdDTiVMeZygwLSqtcz7wCzM7u2VOr82JU7roE4+qhLv740CHxHWPZJe4paKedPdN7v4x8CxwbAudX1vTnFJYF6KhtWzifJ6XEYZ/3d2XAO8DR2Q6aCkEnazleNy92htL6zwIXOXuj7T8qbYJWT9PM6uKxngxs+GEv7O1LX6mxS9OqahHgRPMbK9oSGgE8GYLn2dbEefzxMz2BU4ifLaSXpzPcwWhF46ZVRKqTy/NdNB2P7zm8crxSEwxP8/zgW+Z2Q5gC3ChR+kt0ijOZ+nub5rZk8AioB64x921HEcKzfi3fg7wlLtvaqVTbRNifp4/Bn5tZq8RhuN+GPXI01IZHBERaTGlMLwmIiJFQkFHRERajIKOiIi0GAUdERFpMQo6IiLSYhR0RAokKl0zNM/H7GFmVyU8H2NmGVfVTWj7czM7MU/nMc3MVGVCmk1BR6Rt6QFclbVVkqh8/0h3fzZP5zEJuC5Px5ISoqAjJcXMrjOzf4l+/x8z+0v0+6lm9tvo90lmNj9aH+SmaNsZZvZAwnHGmNlj0e+nm9kLZvaymf3BzLqleN+UbcxsmZndFG1/zcyOiLZXmNmsaPsvzWx5VEroJ8Ch0Xowt0aH72ZmD5rZW2Y2taEaRJLzgScTzmdZQ2kiMxtqZrOj3280s/vM7Kmozblmdkt0bk+aWYfoEM8Bp5lZu7/BXPJLQUdKzbPACdHvQwlf2B2A4wlfpAAT3X0ocAxwkpkdA8wCRprZ3lGbrwG/j764rwdOc/fBwHzg2sQ3jNHm42j7JOD70bYfEdZ6GQz8ETg42j4BeM/dB7n7D6JtxwH/CgwgrH0yOsV1jwYWxPmACIUb/55Qxv63wDPufjShusTfA0SVhJegOnDSTAo6UmoWAEPMrDuwDXiBEHxOoDHo/IOZvQy8AhxFWLhqB6Gn8JXof/d/T6jdNZLwZf9XM3sVuATom/Se2do8nHBu/aLfjyesX4K7Pwmsz3BN89y9NgoEryYcI1EvYE2GYyR6wt23E0rVl9PYQ3ot6dirgS/EPKYIUAK110QSufv2qJr4ZcDfCDXNTib87/5NM6sm9DaGuft6M/s10Dna/ffAtwkLVb3k7p9GQ1mz3P2iDG+brc226OdOGv9NNmeNjW0JvyceI9EWGq8j8bwAOiRt3wahN2Nm2xPq5tUnHbtzdFyR2NTTkVL0LCGwPEvo3YwHXo2+XPcBNgEbo6q5ZyTsNxsYDFxBCEAAc4HRZvZ3EBYHM7PDkt4vTptkzwP/ELU/Hdgv2v4p0L1ZVxu8Cfxd0rZh0c8xhB5Ncx0GLM5hPylhCjpSip4jDDe94O51wNZoG+6+kDCsthiYAvy1YSd33wnMIASiGdG2NcClwP+Z2SJCgNltPZE4bVK4CTg9GuY7A1gJfOruawnDdK8nJBLE8SdCcEn0dTObD1QCfczs3LgHiwLyFndf2YxzEFGVaZFiZGadgJ1ReflRwCR3H7SHx3we+LK7b4iGGIdmK0Of4VjfBT5x91/tyTlJ6dGcjkhxOhh4wMzKgM8JQ3p76nvRcTfk4VgbgN/k4ThSYtTTERGRFqM5HRERaTEKOiIi0mIUdEREpMUo6IiISItR0BERkRbz/wGfQfX24jGdUwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "freqs = mp.get_flux_freqs(refl)\n",
    "wvls = np.divide(1,freqs)\n",
    "\n",
    "eps_quartz = lambda l: 1+0.6961663*math.pow(l,2)/(pow(l,2)-pow(0.0684043,2))+0.4079426*pow(l,2)/(pow(l,2)-pow(0.1162414,2))+0.8974794*pow(l,2)/(pow(l,2)-pow(9.896161,2))\n",
    "R_fresnel = lambda l: math.pow(math.fabs(1-math.sqrt(eps_quartz(l)))/(1+math.sqrt(eps_quartz(l))),2)\n",
    "R_analytic = [R_fresnel(i) for i in wvls]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(wvls,R_meep,'bo-',label='meep')\n",
    "plt.plot(wvls,R_analytic,'rs-',label='analytic')\n",
    "plt.xlabel(\"wavelength (μm)\")\n",
    "plt.ylabel(\"reflectance\")\n",
    "plt.axis([0.4, 0.8, 0.0340, 0.0365])\n",
    "plt.xticks([t for t in np.arange(0.4,0.9,0.1)])\n",
    "plt.legend(loc='upper right')\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}