{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "Working in 1D dimensions.\n",
      "Computational cell is 0 x 0 x 12 with resolution 200\n",
      "time for set_epsilon = 0.00153399 s\n",
      "-----------\n",
      "field decay(t = 50.0025): 0.25018932773921454 / 0.25018932773921454 = 1.0\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",
      "Field time usage:\n",
      "        connecting chunks: 0.000319719 s\n",
      "            time stepping: 1.3083 s\n",
      "            communicating: 0.475856 s\n",
      "     Fourier transforming: 0.416269 s\n",
      "          everything else: 0.920692 s\n",
      "\n",
      "-----------\n",
      "Initializing structure...\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.00103307 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 35056 (time=87.64), 0.000114105 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": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "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": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEKCAYAAADJvIhZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XuUVNWZ9/Hv0y1XQRFtu4kgtA5eEBW5E7zgJYJOEu8Tnc54yaghRicTkxhm0IlOxlmJmjUZfQ0GDUYTMsSoUSQqYiKoiYigguIVEbAVGuSmchPo5/1jn6aL6rqcLqq6q7t+n7VqddepfU6dU6uph733c55t7o6IiEhLKGvtExARkdKhoCMiIi1GQUdERFqMgo6IiLQYBR0REWkxCjoiItJiChp0zGycmb1tZkvMbEKK183Mbo9eX2Rmg6Ptnc1snpktNLPFZnZT0n7XRMddbGa3RNv6mdkWM3s1etxVyGsTEZHm26tQBzazcuBO4EtALfCSmU139zcSmp0B9I8eI4BJ0c9twCnu/pmZdQCeN7Mn3H2umZ0MnAUc4+7bzOzAhOO95+6DCnVNIiKyZwrZ0xkOLHH3pe7+OTCNECwSnQXc78FcoIeZ9Yqefxa16RA9Gu5i/RbwE3ffBuDuqwt4DSIikkcF6+kABwEfJDyvJfRisrU5CFgZ9ZQWAH8H3OnuL0ZtDgNOMLObga3A9939pei1ajN7BfgEuN7dn0s+KTO7ErgSYO+99x5yxBFH7MElioiUngULFnzs7hW57FvIoGMptiXX3Enbxt13AoPMrAfwRzMb6O6vE855P2AkMAx4wMwOAVYCB7v7WjMbAjxiZke5+ye7Hdx9MjAZYOjQoT5//vzcr1BEpASZ2fJc9y3k8Fot0CfheW/go+a2cfcNwGxgXMI+D0dDcPOAeuAAd9/m7mujfRYA7xF6RSIiUiQKGXReAvqbWbWZdQQuBKYntZkOXBxlsY0ENrr7SjOriHo4mFkX4DTgrWifR4BTotcOAzoCH0f7lEfbDyEkJywt4PWJiEgzFWx4zd13mNnVwEygHJji7ovNbHz0+l3A48CZwBJgM3BZtHsv4L4oiJQBD7j7jOi1KcAUM3sd+By4xN3dzE4E/tPMdgA7gfHuvq5Q1yciIs1npby0geZ0RErT9u3bqa2tZevWra19KkWtc+fO9O7dmw4dOuy23cwWuPvQXI5ZyEQCEZGiVFtbS/fu3enXrx9mqfKZxN1Zu3YttbW1VFdX5+24JV0GZ8EC6NcPpk5t7TMRkZa0detW9t9/fwWcDMyM/fffP++9wZIOOgDLl8OVVyrwiJQaBZzsCvEZlXzQAdi8GSZObO2zEBFp/xR0IitWtPYZiIi0fwo6kYMPbu0zEJFiNXVqmP8tK9M88J5S0AHM4NNP9QclIk1NnRrmfZcvB/f8zAMvW7aMI444gssvv5yBAwdSU1PD008/zejRo+nfvz/z5s1j06ZNfOMb32DYsGEcd9xxPProowDs3LmTH/zgBwwbNoxjjjmGX/7ylwDMnj2bE088kXPOOYcBAwYwfvx46uvr8/ER5FXJp0x36waffQbrottIG/6gAGpqWu+8RKRl/Ou/wquvpn997lzYtm33bZs3wz//M9x9d+p9Bg2Cn/888/suWbKEP/zhD0yePJlhw4bxu9/9jueff57p06fz3//93wwYMIBTTjmFKVOmsGHDBoYPH85pp53G1KlT2XfffXnppZfYtm0bo0eP5vTTTwdg3rx5vPHGG/Tt25dx48bx8MMPc/755zfj0yi8ku7pDBkCPXs23a7EAhFpkBxwsm2Pq7q6mqOPPpqysjKOOuooTj31VMyMo48+mmXLlvHUU0/xk5/8hEGDBjFmzBi2bt3KihUreOqpp7j//vsZNGgQI0aMYO3atbz77rsADB8+nEMOOYTy8nIuuuginn/++T07yQIo+Z7OBx+k3q7EApHSkK1H0q9fGAFJ1rcvzJ6d+/t26tRp1+9lZWW7npeVlbFjxw7Ky8t56KGHOPzww3fbz9254447GDt27G7bZ8+e3STFuRjTwku6pwPpEwiUWCAiADffDF277r6ta9ewvZDGjh3LHXfcQUOpsldeeWXX9kmTJrF9+3YA3nnnHTZt2gSE4bX333+f+vp6fv/733P88ccX9iRzUPJBJ9UfFIT/2SipQERqamDy5NCzMQs/J08u/JzvDTfcwPbt2znmmGMYOHAgN9xwAwCXX345AwYMYPDgwQwcOJBvfvOb7NixA4BRo0YxYcIEBg4cSHV1Neecc05hTzIHKvg5fz5Tp4Y5nFRd6K5dW+YPTERazptvvsmRRx7Z2qeRV7Nnz+a2225jxowZ2Rs3Q6rPak8KfpZ8TwdCQFm2LPwPJpmSCkRE8qfkEwkSpUseUFKBiBS7MWPGMGbMmNY+jazU00mgpAIRkcIq7aCzYEGYGTSDqqqUSQWqViAikj+lHXQS1dU1yVLp1i2UvVi3Ln/lL0RESpmCTpKGpIL6elUrEBHJNwWdDFStQESK2a9//WuuvvrqrG0++uijXc8vv/xy3njjjUKfWloKOoka5neiOR4lFogIVVW7fzckfEe0BclB55577mHAgAGtdj4KOunU1aVMLCgvD1WplVggUiLq6pq3Paazzz6bIUOGcNRRRzF58mQAunXrxsSJEzn22GMZOXIkddF7PPbYY4wYMYLjjjuO0047bdf2Bp9++inV1dW7SuN88skn9OvXjz/84Q/Mnz+fmpoaBg0axJYtWxgzZgzz588H4Mknn2Tw4MEce+yxnHrqqXt0PXHpPp0MGqoQTJwYhtT22gu2b4e1a8N2LYMg0g5kW9sgk3T3xcRY22DKlCn07NmTLVu2MGzYMM477zw2bdrEyJEjufnmm7nuuuu4++67uf766zn++OOZO3cuZsY999zDLbfcws9+9rNdx+revTtjxozhT3/6E2effTbTpk3jvPPO44ILLuDOO+/ktttuY+jQ3QsIrFmzhiuuuIJnn32W6upq1jWs71Jgpd3TGTIka5PExILKyqavK7FARHJx++237+rRfPDBB7z77rt07NiRL3/5ywAMGTKEZcuWAVBbW8vYsWM5+uijufXWW1m8eHGT411++eXce++9ANx7771cdtllGd9/7ty5nHjiiVRXVwPQM1XmVAGop9MMH36YersSC0TasGxrG2RaHiDHtQ1mz57N008/zQsvvEDXrl13rZfToUOHXcsRlJeX7yrkec0113Dttdfy1a9+ldmzZ3PjjTc2Oebo0aNZtmwZc+bMYefOnQwcODDjObh7qyx9UNo9nWySJgyVWCAi+bBx40b2228/unbtyltvvcXcuXOztj/ooIMAuO+++9K2u/jii7nooot26+V0796dTz/9tEnbUaNGMWfOHN5//30ADa+1mFRjZsmiSTstgyBSgtJ9R8T57khj3Lhx7Nixg2OOOYYbbriBkSNHZmx/4403csEFF3DCCSdwwAEHpG1XU1PD+vXrueiii3Ztu/TSSxk/fvyuRIIGFRUVTJ48mXPPPZdjjz2Wr33tazlfT3NoaYMoiwPI3I2OPictgyDS9rXHpQ0AHnzwQR599FF+85vf5O2YWtqglWkZBBEpRtdccw0TJkzYtdhbsVIiQY60DIKIFJM77rijtU8hFvV04lK1ApF2pZSnFuIqxGekoJMo7sRgmmoFABs2qFqBSLHr3Lkza9euVeDJwN1Zu3YtnTt3zutxNbyWaNWq3Z9nSCxIrlbQowesXw8bN4btqlYgUrx69+5NbW0ta9asae1TKWqdO3emd+/eeT2mstcSs9eSxchma9CvX+qMtr59Q+KBiEh7oey11pA0x6PEAhGR7AoadMxsnJm9bWZLzGxCitfNzG6PXl9kZoOj7Z3NbJ6ZLTSzxWZ2U9J+10THXWxmtyRs/7foWG+b2dhCXttu6uqUWCAiEkPBgo6ZlQN3AmcAA4CLzCx5EYczgP7R40pgUrR9G3CKux8LDALGmdnI6LgnA2cBx7j7UcBt0fYBwIXAUcA44BfROeSuGXccp0ssWL9eiQUiIg0K2dMZDixx96Xu/jkwjRAsEp0F3O/BXKCHmfWKnn8WtekQPRomUb4F/MTdtwG4++qEY01z923u/j6wJDqH3K1aFeZuGh4Z1NSEagR9+4YRt4aCrZ98EnZtSCxQ4BGRUlbIoHMQkLjgc220LVYbMys3s1eB1cAsd38xanMYcIKZvWhmc8xsWDPeDzO70szmm9n8fGeuJC6D0L1709dVsUBESl0hg06q1K/k7kLaNu6+090HAb2B4WbWUKd7L2A/YCTwA+ABC/W547wf7j7Z3Ye6+9CKiop4VxJHUkVqJRaIiDRVyKBTC/RJeN4b+Ki5bdx9AzCbME/TsM/D0RDcPKAeOCDm++2ZZlSkTpdAUFamOR4RKV2FDDovAf3NrNrMOhIm+acntZkOXBxlsY0ENrr7SjOrMLMeAGbWBTgNeCva5xHglOi1w4COwMfRsS40s05mVk1ITpiX1ytKnOPJIl1iwc6dmuMRkdJVsKDj7juAq4GZwJvAA+6+2MzGm9n4qNnjwFLCpP/dwFXR9l7AM2a2iBC8Zrn7jOi1KcAhZvY6ITnhkqjXsxh4AHgDeBL4trvvLNT1ZZOcWFCeIo9OczwiUmpUkSBTRYJM4i7zWlkJq1ZRVpa6g2QWEg9ERNoKVSQoZlnmeHTzqIiUEgWdXDVzqdpUczxm4T4eJRaISKlQ0MlVM24chaZzPPvsE3Zbv16JBSJSOhR0WlDizaP77df0dSUWiEh7p6DTEpIqUoNuHhWR0qSgky/NWHUUlFggIqVJQSdfmjnHk+7m0Y8/VmKBiLRfCjqtJDmxYP/9w/ZNm5RYICLtl4JOa4jmd2q+bizbWkV9PXTr1rSZEgtEpL1R0Glt0RyPEgtEpBQo6BRKM28eTZdAsO++YX5H8zwi0h4o6BRKnhILNmwI8zua5xGR9kBBp0gkJxb07Qs9ejRtp3keEWnL9mrtExB2VayuAWoqK6F+FRCG1FLRPI+ItFXq6bSUZt48CrqBVETaHwWdltLMOR5IP8+zcaMSC0SkbVLQKWLJ8zz77hu2b9igxAIRaZsUdIpRQoHQmu9V7apMrcQCEWnrFHSKXcIcj24gFZG2TkGntTTz5lFIn0DQo4duIBWRtkFBp7XkMbFg/XrdQCoibYOCTluQUCB002ZjdVnVrhtI99mnaXPN84hIsdLNoW1QRX0d9VHnSDeQikhbop5OschhjgfSz/OUlWmOR0SKj4JOsUie49nDeZ6dOzXHIyLFR0GnrUqY5/m4Q9WuG0jLy5s21RyPiBQLzem0A1021rFsQ/hdczwiUszU0ylmebyXR0VCRaQYKOgUszzey1NZGVKslVwgIq1Jw2vtTE1N+DlxYhhS69MHOnWCefMa2zQkFyS2FxFpCerptBeJRUK/bizbWkV9fQgw27Y1ba7kAhFpDQo6bUlz5ngSCoV+8EHqJkouEJGWpqDTluR4L48KhYpIsVDQKQEqFCoixSJW0DGzSjP7lZk9ET0fYGb/HGO/cWb2tpktMbMJKV43M7s9en2RmQ2Otnc2s3lmttDMFpvZTQn73GhmH5rZq9HjzGh7PzPbkrD9rrgfQruV5gbSvn0bVyFNpHkeESm0uNlrvwbuBRq+kt4Bfg/8Kt0OZlYO3Al8CagFXjKz6e7+RkKzM4D+0WMEMCn6uQ04xd0/M7MOwPNm9oS7z432+x93vy3F277n7oNiXlNJSbyBFHQTqYi0jrjDawe4+wNAPYC77wB2ZtlnOLDE3Ze6++fANOCspDZnAfd7MBfoYWa9ouefRW06RI94ExilJs+FQnv21DyPiBRO3KCzycz2J/riN7ORwMYs+xwEJOZN1UbbYrUxs3IzexVYDcxy9xcT2l0dDcdNMbP9ErZXm9krZjbHzE5IdVJmdqWZzTez+WvWrMlyCW1ADjeQQvp5nnXrNM8jIoUTN+hcC0wHDjWzvwL3A9dk2cdSbEv+Vkzbxt13RkNlvYHhZjYwen0ScCgwCFgJ/CzavhI42N2Pi873d2bWZIkzd5/s7kPdfWhFRUWWS2hnku7lSZ7n6dq1adzSPI+I5FOsOR13f9nMTgIOJwSKt919e5bdaoE+Cc97Ax81t427bzCz2cA44HV333UDipndDcyI2m0jzAXh7gvM7D3gMGB+nGssRZrnEZGWFjd77dtAN3df7O6vA93M7Kosu70E9DezajPrCFxI6C0lmg5cHGWxjQQ2uvtKM6swsx7Re3cBTgPeip73Stj/HOD1aHtFlLyAmR1CSE5YGuf62pUc53hA8zwiUnhxh9eucPdd/yd29/XAFZl2iJINrgZmAm8CD7j7YjMbb2bjo2aPEwLDEuBuoCGQ9QKeMbNFhOA1y91nRK/dYmavRa+dDHw32n4isMjMFgIPAuPdfV3M62s/cryBFDTPIyKFZx7jSyn6gj/Wo8ZRj2KRux9V4PMrqKFDh/r8+SUw+mapps5S27JvJUf2WLWrWGhdXerabX37wrJl+TtFEWk7zGyBuw/NZd+49+nMBB6Ibrh0YDzwZC5vKMVN8zwiUkhxh9d+CPwF+BbwbeDPwHWFOinJM83ziEiRiDW81l6VzPBaKpmG3BL+JqZODXM4mzc33T3xT6drV5g8WevziJSCPRlei5u9NtrMZpnZO2a21MzeN7PSywwrFRnu5zn44LAonO7nEZFcxJ3T+RUhS2wB2cvfSDujeR4RyZe4czob3f0Jd1/t7msbHgU9MymsAszz9OmTeruISIO4QecZM7vVzEaZ2eCGR0HPTAqrAPfzQAg8Si4QkXTiDq+NiH4mThw5cEp+T0eKVkLiQQ1w7r6VHNm38X6evn3huecamzfcRApKLhCRRnFrr51c6BORtiV5nqdfv6ZtGpILFHREpEHcng5m9vfAUUDnhm3u/p+FOClpJZWVoQRBDtIlESi5QEQSxU2Zvgv4GmE5AwMuAPoW8LykNTR3nichtbrejZVUpWwyYYJuJBWRIG4iwRfd/WJgvbvfBIxi9yUJRKhi915S586wzz7w05+qYKiIBHGDzpbo52Yz+wKwHaguzClJW5a4KNw990D37k3b6EZSkdIVN+jMiNa3uRV4GVgGTCvUSUkRaeb9PMuWQX19+FlTA7W1qdstX64hN5FSFDeR4JZoZc6HzGwGIZlga+FOS4rGqlVNt2Wq25b0Wl1ZJQd6imMQAk/DT6VXi5SGuD2dFxp+cfdt7r4xcZtIOhX1dWlvJE2kITeR0pAx6JhZlZkNAbqY2XEJ1QjGADG+SqRdauaQ2+TJu8/1pKP0apH2L1tPZyxwG9Ab+FnC47vAvxf21KRoNTO1uqZm97medIGnoiLM7WiuR6T9irtc9Xnu/lALnE+LKun1dPKtGUtiA9RRSRWNcz0N6/PstRfs2NHYTuv0iBSfgq+nAwyJstca3nA/M/uvXN5QBKCSupTp1YkBBzTXI9LexA06Z7j7rkpb7r4eOLMwpyRtUg5LJSQOuX3jG/DZZ6nbaa5HpP2IG3TKzaxTwxMz6wJ0ytBeSk2qeZ5mLoWebp0ezfWItB9x79P5LfBnM7uXsKTBN4D7CnZWUhqS5oHe3LeSA7quYvPm3ZusXg2XXAI7ozVrdV+PSNsVq6fj7rcA/wUcSag0/eNom0jedNlY1yS9+le/CvXbdiYtkq65HpG2KVb2GoCZ9QX6u/vTZtYVKHf3Twt6dgWm7LUWUFXVvOUSUvw9lpWlHqkzC3NCItKyCp69ZmZXAA8Cv4w2HQQ8kssbSonZg+USGh511nTJBIAePeA3v9Fcj0hbEndO59vAcOBFAHd/18wOLNhZiSRoKKWTONdTVgbr18Ollzb2djTXI1L84mavbXP3zxuemNlehIQCkRaRPNdz//1wwAFNh9c01yNS3OIGnTlm9u+EGmxfAv4APFa405J2LYd7epJL6dTUwNq1qdvqvh6R4hU36EwA1gCvAd8EHgeuL9RJSTuXyz09zZjr6dQJbr1Vcz0ixSh29lp7pOy1ItPM+m0Ae3f13eZ6OnQIP7dv372dariJ5E/BstfM7DUzW5TukdvpiqSRw7Bb8lzPvfemPozmekSKQ8aejpkdAWxJ97q7Ly/ESbUU9XTaiEw9IN3XI9LiCnmfzu+iwPJf7r48+ZHLG4rkVTPmetzhiitg0iTN94i0lmz36XQ0s0uAL5rZuckvuvvDhTktkdyluq+nSxc46aRQVueeexq3694ekZaVraczHhgJ9AC+kvT4craDm9k4M3vbzJaY2YQUr5uZ3R69vsjMBkfbO5vZPDNbaGaLzeymhH1uNLMPzezV6HFmwmv/Fh3rbTMbG+cDkDYgD3M9d98NTzwBvXo1bav5HpEW5O5ZH8A/x2mXtE858B5wCNARWAgMSGpzJvAEYITg9mK03YBu0e8dCJUQRkbPbwS+n+L9BkTv0Qmojt67PNM5DhkyxKUNS514nf5RWelmqV8yc//tb9379g2/9+0bnotIU8B8b2ZMaHjEvU9nmpldb2aTAcysv5ll6+kMB5a4+1IP1QymAWcltTkLuD+6jrlADzPrFT1vWNKrQ/TIltt9FjDN3be5+/vAkugcRIK6urRr9riH5ROWLw+/Nwy7ab5HJL/iBp0pwOfAF6PntYSlDjI5CPgg4XlttC1WGzMrN7NXgdXALHd/MaHd1dFw3BQz268Z74eZXWlm881s/po1a7JcgrQ3N98c7tlJ1KULdO6s5RNEWkLcoHOoh/VztgO4+xbCEFgmqV5P7q2kbePuO919ENAbGG5mA6PXJwGHAoOAlcDPmvF+uPtkdx/q7kMrKiqyXIIUtRzL6aSa79m2LXV7ldQRya+4QefzaInqMOFidiiQ5p/pLrVAn4TnvYGPmtvG3TcAs4Fx0fO6KCDVA3fTOIQW5/2kPcmlnE6PHtR83Vi23Kj38LPm6+nTrMvKQsabllAQyY+sQcfMDLgLeBLoY2ZTgT8D12XZ9SWgv5lVm1lH4EJgelKb6cDFURbbSGCju680swoz6xG9fxfgNOCt6Hli/tE5wOsJx7rQzDqZWTXQH5iX7fqkxGzcmHJzQ5p1ok6dQoC5/HLN94jkS9agE2UqfAc4F7gU+D9gqLvPzrLfDuBqYCbwJvCAuy82s/FmNj5q9jiwlDDpfzdwVbS9F/BMVGrnJcKczozotVsayvMAJwPfjd5vMfAA8AYhQH7b3ZNG6aUkpBt2yzIcl2qp7HffhYqKph0ozfeI5CZWwU8zuxP4tbu/VPhTajkqg1OCmltUtLKSstWrMpbVmTo1BKAVK+Dgg0Oygm40lfZsT8rgxF059GTgm2a2HNhEmLR3dz8mlzcVaTPq6ji4bxhSS+YORx8NS5bA1q1hmyociGQWN5HgDELG2Ck0ViP4SqFOSqSYpEuz/sd/hDfeaAw4DTT0JpJerKDjKYp9ugp+SluUxzTrqVPTJ8sp1Voktbg9HZH2IZc0627dUqZZU1WVtsJBp04wb14ITEq1FmkUd05HpHRt2pR6e10dN/82zOEkr15aXg4jRoRg07CGj+Z7RNTTEQlyTLNONfR2772wciX06NF00TjN90ipU9ARgfTDbqtWZd7vwANTDr1171+V7j5UVqxoDEYafpNSo+E1kT2RrmhsllTrY46BE0+E++5rHJrT8JuUAvV0RAokVap1167wrW+F3ydN2n0uCDT8Ju2fgo5INrnO93yvik2bDafxsWmz8YuHq1i0KH1xBKVbS3umoCOSTa7zPXV1abeXlZE23XrvvWHpUs33SPukOR2RVnLzzU3TrcvLw/NDDw2/Nywsp/keaS/U0xHZEzkOvUHqdOv77oMPPoB99tFKptI+KeiI7Ilch97+5V+goiJluvUXBlfx6aepd1u+vHFJIA2/SVuk4TWR1nDXXbB9e+rXMqRbA/TuDaNGwXPPqbq1tD3q6YgUSqaht2XLMu6aLt36xz+Gc86BWbNU3VrapliLuLVXWsRNWlWmBeUqKlLfeFpZCatWUVaWuk5pw8JyIoW0J4u4qacjUowyVDqA9OnWAFddBYsXa85HipOCjkgblGr4rXNnOP54mDIFBg6Eiy8Ocz3ujXM+CjzS2hR0RFpLrunW772XstrBlq3Gs+9UUVurCtdSvBR0RFpLrunW/ftnrHZwwAGkrXC9fDnMmNF4D5CG4KSlKeiItDU/+lHWJunmfMrK4CtfgepqOO88uOIKDcFJy1LQESlGmYbesgWdnTt5c0PVbkNvDY/PulXx0ENwxBHw8MOwZcvuu2oITgpNKdNKmZa2KFO6dZ8+oZZOOtG/eaVdS66UMi0ijY48MlazdENw7jBmDPzmN2Hpbc35SD4p6Ii0RZmG32bOjHWIVGnXXbrABRfAhx+GlOtvfENzPpJfCjoibVGumW8AI0bAL39JzbWVTdKuN28xHni2infeSR3XNm+GH/4w/5cjpUNBR6TUbN4M48fD6tWpX6+rwyz9yx9+CGPHhh7P5s1Ku5bmUZVpkfaosjL1vTyVlbBoESxYAMOGZTzEwQenrnS9777w9tvw9a9Dp06wY4cWm5P41NMRaY8yDb+ZwdAsiUe3385b6w5MmXa9kiqWLoXZs2GvvbTYnDSPgo6INPWd79D509RFR7tsrKOsDE46afelthMtXw633RaG4hpoGE5AQUdEUnnttVjN0qVdd+wIP/hBuGXolFNC5QNVPxBQ0BEpXZnSrgcOzLzvd74DL76YtvLBhi4hA+4//iPcp3rPPap+IIEqEqgigUhqmaoedOoE27Zl3j/6bnGH8vL01Q+2bAmHg9DzmTgRVqwIvaibb1ZCQjEq2ooEZjbOzN42syVmNiHF62Zmt0evLzKzwdH2zmY2z8wWmtliM7spxb7fNzM3swOi5/3MbIuZvRo97irktYmUtLq6UK4gBrPM1Q8OPBAuuQSuuy4MuWkIrn0rWNAxs3LgTuAMYABwkZkNSGp2BtA/elwJTIq2bwNOcfdjgUHAODMbmXDsPsCXgBVJx3vP3QdFj/H5viaRkpJp+G3ffeHSSzPv/9OfwtKlAGmH4T7pWsW558Kjj8KttzZNTNAQXPtTyJ7OcGCJuy9198+BacDRD3FhAAAQNElEQVRZSW3OAu73YC7Qw8x6Rc8/i9p0iB6JnfP/Aa5L2iYi+bQnVQ8AJkyAQw+FYcPosjH1+j/dN9dx772h45RuNG/FCvj88xyvQYpOIYPOQUBiqdvaaFusNmZWbmavAquBWe7+YrT9q8CH7r4wxXtWm9krZjbHzE5IdVJmdqWZzTez+WvSrUMvInvu/fdD9yXT3FCkU6fsQ3D/9E+hR7Rli9Kv27JCBp1Uf2nJPZO0bdx9p7sPAnoDw81soJl1BSYC/5Fiv5XAwe5+HHAt8Dsz26fJwd0nu/tQdx9aUVHRjMsRkSYyDcH16wff/z7Mm5f5GIsXg3vaIbiN0RDcn/4EZ58dluK+5BLN/bRVhQw6tUCfhOe9gY+a28bdNwCzgXHAoUA1sNDMlkXtXzazKnff5u5ro30WAO8Bh+XrYkQkhT0dgoOQnn3kkWmH4PbZXMeUKWEI7qmnwj1AqaogTGiSqiTFqJBB5yWgv5lVm1lH4EJgelKb6cDFURbbSGCju680swoz6wFgZl2A04C33P01dz/Q3fu5ez9C0Brs7quifcqjfQ4hJCcsLeD1iUg+/OIX0Lt31mYdOsCXvgSbNqV+vbY2rAP085+HkT3QMFwxKljQcfcdwNXATOBN4AF3X2xm482sIbPscUJgWALcDVwVbe8FPGNmiwjBa5a7z8jylicCi8xsIfAgMN7d1+X1okQkN5mG4b71LXj66cz7P/LIrmhTZ6mH4VZZFevWwXe/C4ccEuaILrtMw3DFRjeH6uZQkeKQLeGgc+fQ1XnssfRt3HnvvRCjJk5Mff9quurZEl/R3hwqIpIXf/5z6KYsWpS16aGHwve+lz7NesWKsDrq/ffDxx83btdQXMtQT0c9HZHiUFWVfg2ghsQE9xAV0nnqKTjxROjUiTXlVVTUNz1eHZUc12sVK1eGQ40aFaaUHn0Utm5tbNe1K0yerDI8qexJT0dBR0FHpG3JNgzXrVtY2vShh9I2qd/pvPJKGKl77DF4+eXU7fr2hWXLcj/V9krDayIiECJITQ3MnZuxWVkZDBkCN94YFlFNF8eWLw/JdQ3ZcA00FJc7LVctIm1LpqW4v/zl8Mg2DDdtGpx+OvTsCYSMuApvesxVVNLr22Fo74gj4MwzQ/WE//3fxjpxWqK7eTS8puE1kfYp2zBcWRmMGAFnnBEW/knjnbedJ56Axx8PS3SnS1AopaE4Da+JiDTHCy/A9dfDjh3wox9lbHrYYWHNupkzYd26zENxM2bAZ581fU3DcY3U01FPR6R9ipMNB7B6dfqbVwGeeQa++MVQfwfSZsWtopJerKJDBzj++JDLMHYsvP46fPObuy/b0NYz45S9liMFHREBsg/F7b13qLFz+umh25PG07OcmTNDr+i118K2sjKor2/ati0Pxyno5EhBR0SAzEHn0UdDFHnqKViyJPNxEr5PP/oIZs1Kv9adGWzcCN27N/90W5vmdERE9kSm2nBf/SrceSe8+y68917m4zz55K4acV/4QliCYXVZ6lpxH3kVPXvCCSfATTfBX/8K27e3//kf9XTU0xGR5sg2FNehQyhzcOqpcNppMHp02qb/NsF5+ulwr5B7SMfesWP3pRuKcf5Hw2s5UtARkWbLFHRmzgx14v7851DmINv3a/T6unUhX+HSS1Nnv/XsGRLu+vePtRBrwWl4TUSkpWQaijv9dPjpT2H+/FBN9MEHMx/r9dfBnZ494bzz4N3PUg/FLV5XxeGHQ58+YdnuKVN2T0JoS0NyqkggItIccVdFbYgkmRx9NFRUwEknwcknU0Xq1VOrqGPSpNAbmjkTfvvbsL1fv7BUw9y5jTetFnuFBA2vaXhNRAop03jYvfeGSPLMM/DBB5mPE31Xu8MbbzTu9sgjqVOyDzoorKZaCBpeExEpVpmG4y69FO67L3RPsmXG3X47LFyI1e/kqKPg6qtDIe0P61MPyc3/sIq+feHii+FXvwrJd4l9jNYaklNPRz0dESkWcbIEevQIJQ9OPDE8Ro5M2/SC8505c2DNmvC8V68wkte5c6h5muv6Qcpey5GCjogUlUxBZ/lyeO45mDMHnn0W3n47+/HccYe33mrcbc6ccONqKnGX8lbQyZGCjogUlbj14iC0e+65sPZ2On/8Y7hPqKJi1yZ3qCurSpm0sIpKak5ZxQknhM7UyJFhTbwGU6fCxImwfPlQ3OfnlLytoKOgIyJtWZwhucMPD1Fk9Ojw87DD0jYdfJyzcGFITigvh0GDQtWE+nq4+27YsgUg96CjlGkRkfbqr3+F558Pjz/+MWQUZPHyy6Em3Ny5YbfnnoO77tp9/mdPqKejno6ItGVxh+Tq68M80PPPN97Ik8rDD4cyPlVVuzZ9/jms69Q4JDcUmO+u4bXmUtARkZIUZ0iuujqsIzRqVPg5ePCul/Yk6Gh4TUREGv3tb6HQ29/+Bn/5S95v4NHNoSIipSbTDaujRsG114a6cR9+GIq8/e53eXtrBR0RkVKzalXInU5+JKdlm4UlTi+6KG9vraAjIiItRkFHRESySzck10xKJBARkewSht4WmC3I9TDq6YiISItR0BERkRajoCMiIi1GQUdERFqMgo6IiLSYggYdMxtnZm+b2RIzm5DidTOz26PXF5nZ4Gh7ZzObZ2YLzWyxmd2UYt/vm5mb2QEJ2/4tOtbbZja2kNcmIiLNV7CUaTMrB+4EvgTUAi+Z2XR3fyOh2RlA/+gxApgU/dwGnOLun5lZB+B5M3vC3edGx+4THXdFwvsNAC4EjgK+ADxtZoe5+85CXaOIiDRPIXs6w4El7r7U3T8HpgFnJbU5C7jfg7lADzPrFT3/LGrTIXoklsP+H+C6pG1nAdPcfZu7vw8sic5BRESKRCFvDj0I+CDheS2hF5OtzUHAyqintAD4O+BOd38RwMy+Cnzo7gtt9/LcBwFzUxxrN2Z2JdCwmMQ2M3u9mdcl6R0AfNzaJ9GO6PPMH32W+XV4rjsWMuikWmshefGetG2iYbFBZtYD+KOZDQSWAhOB03N8P9x9MjAZwMzmu/vQtFcgzaLPM7/0eeaPPsv8MrOcFyIr5PBaLdAn4Xlv4KPmtnH3DcBsYBxwKFANLDSzZVH7l82sKub7iYhIKypk0HkJ6G9m1WbWkTDJPz2pzXTg4iiLbSSw0d1XmllF1MPBzLoApwFvuftr7n6gu/dz936EQDPY3VdFx7rQzDqZWTUhOWFeAa9PRESaqWDDa+6+w8yuBmYC5cAUd19sZuOj1+8CHgfOJEz6bwYui3bvBdwXzeuUAQ+4+4ws77fYzB4A3gB2AN+Okbk2OberkzT0eeaXPs/80WeZXzl/nubeZNpDRESkIFSRQEREWoyCjoiItJiSCDrZyvEktBtmZjvN7PyWPL+2JkZ5ozFmttHMXo0e/9Ea59kWxPnbjD7PV6OSUHNa+hzbkhh/mz9I+Lt8Pfr33rM1zrUtiPF57mtmjyWULLss1XF24+7t+kFIYngPOAToCCwEBqRp9xdCcsP5rX3exfqI83kCY4AZrX2uxf6I+Vn2ICTHHBw9P7C1z7tYH3H/rSe0/wrwl9Y+72J9xPz7/Hfgp9HvFcA6oGOm45ZCTydOOR6Aa4CHgNUteXJtUNzPU7KL81n+I/Cwu68AcHf9fabX3L/Ni4D/a5Eza5vifJ4OdLdQHqYbIejsyHTQUgg66Urt7GJmBwHnAHe14Hm1VVk/z8ioqMv9hJkd1TKn1ubE+SwPA/Yzs9lmtsDMLm6xs2t74v5tYmZdCTecP9QC59VWxfk8/x9wJOFG/NeA77h7faaDFrIMTrGIUx7n58AP3X1nUj03aSrO5/ky0NdDlfAzgUcIN+vK7uJ8lnsBQ4BTgS7AC2Y2193fKfTJtUGxSmFFvgL81d3XFfB82ro4n+dY4FXgFELFmFlm9py7f5LuoKXQ04lTHmcoMC0qrXM+8AszO7tlTq/NiVO66BOPqoS7++NAh8R1j2SXuKWinnT3Te7+MfAscGwLnV9b05xSWBeiobVs4nyelxGGf93dlwDvA0dkOmgpBJ2s5XjcvdobS+s8CFzl7o+0/Km2CVk/TzOrisZ4MbPhhL+ztS1+psUvTqmoR4ETzGyvaEhoBPBmC59nWxHn88TM9gVOIny2kl6cz3MFoReOmVUSqk8vzXTQdj+85vHK8UhMMT/P84FvmdkOYAtwoUfpLdIozmfp7m+a2ZPAIqAeuMfdtRxHCs34t34O8JS7b2qlU20TYn6ePwZ+bWavEYbjfhj1yNNSGRwREWkxpTC8JiIiRUJBR0REWoyCjoiItBgFHRERaTEKOiIi0mIUdEQKJCpdMzTPx+xhZlclPB9jZhlX1U1o+3MzOzFP5zHNzFRlQppNQUekbekBXJW1VZKofP9Id382T+cxCbguT8eSEqKgIyXFzK4zs3+Jfv8fM/tL9PupZvbb6PdJZjY/Wh/kpmjbGWb2QMJxxpjZY9Hvp5vZC2b2spn9wcy6pXjflG3MbJmZ3RRtf83Mjoi2V5jZrGj7L81seVRK6CfAodF6MLdGh+9mZg+a2VtmNrWhGkSS84EnE85nWUNpIjMbamazo99vNLP7zOypqM25ZnZLdG5PmlmH6BDPAaeZWbu/wVzyS0FHSs2zwAnR70MJX9gdgOMJX6QAE919KHAMcJKZHQPMAkaa2d5Rm68Bv4++uK8HTnP3wcB84NrEN4zR5uNo+yTg+9G2HxHWehkM/BE4ONo+AXjP3Qe5+w+ibccB/woMIKx9MjrFdY8GFsT5gAiFG/+eUMb+t8Az7n40obrE3wNElYSXoDpw0kwKOlJqFgBDzKw7sA14gRB8TqAx6PyDmb0MvAIcRVi4agehp/CV6H/3f0+o3TWS8GX/VzN7FbgE6Jv0ntnaPJxwbv2i348nrF+Cuz8JrM9wTfPcvTYKBK8mHCNRL2BNhmMkesLdtxNK1ZfT2EN6LenYq4EvxDymCFACtddEErn79qia+GXA3wg1zU4m/O/+TTOrJvQ2hrn7ejP7NdA52v33wLcJC1W95O6fRkNZs9z9ogxvm63NtujnThr/TTZnjY1tCb8nHiPRFhqvI/G8ADokbd8GoTdjZtsT6ubVJx27c3RckdjU05FS9CwhsDxL6N2MB16Nvlz3ATYBG6OquWck7DcbGAxcQQhAAHOB0Wb2dxAWBzOzw5LeL06bZM8D/xC1Px3YL9r+KdC9WVcbvAn8XdK2YdHPMYQeTXMdBizOYT8pYQo6UoqeIww3veDudcDWaBvuvpAwrLYYmAL8tWEnd98JzCAEohnRtjXApcD/mdkiQoDZbT2ROG1SuAk4PRrmOwNYCXzq7msJw3SvJyQSxPEnQnBJ9HUzmw9UAn3M7Ny4B4sC8hZ3X9mMcxBRlWmRYmRmnYCdUXn5UcAkdx+0h8d8Hviyu2+IhhiHZitDn+FY3wU+cfdf7ck5SenRnI5IcToYeMDMyoDPCUN6e+p70XE35OFYG4Df5OE4UmLU0xERkRajOR0REWkxCjoiItJiFHRERKTFKOiIiEiLUdAREZEW8/8Bn0H19uIxnVMAAAAASUVORK5CYII=\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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}