{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "Meep: using complex fields.\n",
      "Meep progress: 10.11/200.0 = 5.1% done in 4.0s, 75.2s to go\n",
      "Meep progress: 19.91/200.0 = 10.0% done in 8.0s, 72.4s to go\n",
      "Meep progress: 30.03/200.0 = 15.0% done in 12.0s, 67.9s to go\n",
      "Meep progress: 40.33/200.0 = 20.2% done in 16.0s, 63.4s to go\n",
      "Meep progress: 50.46/200.0 = 25.2% done in 20.0s, 59.3s to go\n",
      "Meep progress: 60.7/200.0 = 30.4% done in 24.0s, 55.1s to go\n",
      "Meep progress: 70.97/200.0 = 35.5% done in 28.0s, 50.9s to go\n",
      "Meep progress: 81.13/200.0 = 40.6% done in 32.0s, 46.9s to go\n",
      "Meep progress: 91.17/200.0 = 45.6% done in 36.0s, 43.0s to go\n",
      "Meep progress: 101.35000000000001/200.0 = 50.7% done in 40.0s, 39.0s to go\n",
      "Meep progress: 111.68/200.0 = 55.8% done in 44.0s, 34.8s to go\n",
      "Meep progress: 121.94/200.0 = 61.0% done in 48.0s, 30.7s to go\n",
      "Meep progress: 132.32/200.0 = 66.2% done in 52.0s, 26.6s to go\n",
      "Meep progress: 142.6/200.0 = 71.3% done in 56.0s, 22.6s to go\n",
      "Meep progress: 152.43/200.0 = 76.2% done in 60.0s, 18.7s to go\n",
      "Meep progress: 162.74/200.0 = 81.4% done in 64.0s, 14.7s to go\n",
      "Meep progress: 172.87/200.0 = 86.4% done in 68.0s, 10.7s to go\n",
      "Meep progress: 183.22/200.0 = 91.6% done in 72.0s, 6.6s to go\n",
      "Meep progress: 193.47/200.0 = 96.7% done in 76.0s, 2.6s to go\n",
      "run 0 finished at t = 200.0 (20000 timesteps)\n",
      "-----------\n",
      "Initializing structure...\n",
      "Meep: using complex fields.\n",
      "Meep progress: 9.91/300.0 = 3.3% done in 4.0s, 117.2s to go\n",
      "Meep progress: 19.97/300.0 = 6.7% done in 8.0s, 112.3s to go\n",
      "Meep progress: 29.560000000000002/300.0 = 9.9% done in 12.0s, 109.9s to go\n",
      "Meep progress: 39.53/300.0 = 13.2% done in 16.0s, 105.5s to go\n",
      "Meep progress: 49.58/300.0 = 16.5% done in 20.0s, 101.1s to go\n",
      "Meep progress: 59.56/300.0 = 19.9% done in 24.0s, 97.0s to go\n",
      "Meep progress: 69.63/300.0 = 23.2% done in 28.0s, 92.7s to go\n",
      "Meep progress: 79.63/300.0 = 26.5% done in 32.0s, 88.6s to go\n",
      "Meep progress: 89.79/300.0 = 29.9% done in 36.0s, 84.3s to go\n",
      "Meep progress: 99.65/300.0 = 33.2% done in 40.0s, 80.5s to go\n",
      "Meep progress: 109.44/300.0 = 36.5% done in 44.0s, 76.7s to go\n",
      "Meep progress: 119.5/300.0 = 39.8% done in 48.0s, 72.5s to go\n",
      "Meep progress: 129.28/300.0 = 43.1% done in 52.0s, 68.7s to go\n",
      "Meep progress: 138.88/300.0 = 46.3% done in 56.0s, 65.0s to go\n",
      "Meep progress: 148.73/300.0 = 49.6% done in 60.0s, 61.1s to go\n",
      "Meep progress: 158.18/300.0 = 52.7% done in 64.0s, 57.4s to go\n",
      "Meep progress: 168.13/300.0 = 56.0% done in 68.0s, 53.4s to go\n",
      "Meep progress: 178.14000000000001/300.0 = 59.4% done in 72.0s, 49.3s to go\n",
      "Meep progress: 188.20000000000002/300.0 = 62.7% done in 76.0s, 45.2s to go\n",
      "Meep progress: 198.17000000000002/300.0 = 66.1% done in 80.0s, 41.1s to go\n",
      "Meep progress: 208.09/300.0 = 69.4% done in 84.0s, 37.1s to go\n",
      "Meep progress: 218.17000000000002/300.0 = 72.7% done in 88.0s, 33.0s to go\n",
      "Meep progress: 228.17000000000002/300.0 = 76.1% done in 92.1s, 29.0s to go\n",
      "Meep progress: 238.24/300.0 = 79.4% done in 96.1s, 24.9s to go\n",
      "Meep progress: 248.07/300.0 = 82.7% done in 100.1s, 20.9s to go\n",
      "Meep progress: 258.12/300.0 = 86.0% done in 104.1s, 16.9s to go\n",
      "Meep progress: 268.3/300.0 = 89.4% done in 108.1s, 12.8s to go\n",
      "Meep progress: 278.44/300.0 = 92.8% done in 112.1s, 8.7s to go\n",
      "Meep progress: 288.27/300.0 = 96.1% done in 116.1s, 4.7s to go\n",
      "Meep progress: 298.33/300.0 = 99.4% done in 120.1s, 0.7s to go\n",
      "run 0 finished at t = 300.0 (30000 timesteps)\n",
      "refl:, 0, -0.82, 0.00005706\n",
      "refl:, 1, 1.09, 0.00000404\n",
      "refl:, 2, -2.73, 0.00000529\n",
      "refl:, 3, 3.00, 0.00006039\n",
      "refl:, 4, -4.65, 0.00005601\n",
      "refl:, 5, 4.91, 0.00000845\n",
      "refl:, 6, -6.57, 0.00000815\n",
      "refl:, 7, 6.83, 0.00006645\n",
      "refl:, 8, -8.49, 0.00005697\n",
      "refl:, 9, 8.76, 0.00015744\n",
      "refl:, 10, -10.43, 0.00001317\n",
      "refl:, 11, 10.70, 0.04414669\n",
      "refl:, 12, -12.38, 0.00005969\n",
      "refl:, 13, 12.65, 0.00041534\n",
      "refl:, 14, -14.34, 0.00001986\n",
      "refl:, 15, 14.62, 0.00009002\n",
      "refl:, 16, -16.32, 0.00006408\n",
      "refl:, 17, 16.60, 0.00012194\n",
      "refl:, 18, -18.32, 0.00003124\n",
      "refl:, 19, 18.60, 0.00011016\n",
      "refl:, 20, -20.34, 0.00006994\n",
      "refl:, 21, 20.63, 0.00011303\n",
      "refl:, 22, -22.40, 0.00004960\n",
      "refl:, 23, 22.69, 0.00013856\n",
      "refl:, 24, -24.48, 0.00007689\n",
      "refl:, 25, 24.77, 0.00014294\n",
      "refl:, 26, -26.59, 0.00007941\n",
      "refl:, 27, 26.89, 0.00017844\n",
      "refl:, 28, -28.75, 0.00008474\n",
      "refl:, 29, 29.06, 0.00020622\n",
      "refl:, 30, -30.95, 0.00012869\n",
      "refl:, 31, 31.27, 0.00023543\n",
      "refl:, 32, -33.21, 0.00009471\n",
      "refl:, 33, 33.53, 0.00032429\n",
      "refl:, 34, -35.52, 0.00021572\n",
      "refl:, 35, 35.85, 0.00032455\n",
      "refl:, 36, -37.90, 0.00011932\n",
      "refl:, 37, 38.24, 0.00060212\n",
      "refl:, 38, -40.37, 0.00042988\n",
      "refl:, 39, 40.72, 0.00054810\n",
      "refl:, 40, -42.92, 0.00031748\n",
      "refl:, 41, 43.29, 0.00174005\n",
      "refl:, 42, -45.59, 0.00064620\n",
      "refl:, 43, 45.97, 0.00100453\n",
      "refl:, 44, -48.39, 0.00034986\n",
      "refl:, 45, 48.79, 0.00105225\n",
      "refl:, 46, -51.35, 0.00032074\n",
      "refl:, 47, 51.78, 0.00098224\n",
      "refl:, 48, -54.52, 0.00044285\n",
      "refl:, 49, 54.98, 0.00071074\n",
      "refl:, 50, -57.96, 0.00023208\n",
      "refl:, 51, 58.47, 0.00109322\n",
      "refl:, 52, -61.76, 0.00035650\n",
      "refl:, 53, 62.33, 0.00079058\n",
      "refl:, 54, -66.11, 0.00026145\n",
      "refl:, 55, 66.78, 0.00071195\n",
      "refl:, 56, -71.38, 0.00002471\n",
      "refl:, 57, 72.24, 0.00028160\n",
      "refl:, 58, -78.81, 0.00001352\n",
      "tran:, 0, -1.23, 0.00092964\n",
      "tran:, 1, 1.63, 0.01528639\n",
      "tran:, 2, -4.10, 0.00767314\n",
      "tran:, 3, 4.50, 0.00094525\n",
      "tran:, 4, -6.98, 0.00092643\n",
      "tran:, 5, 7.38, 0.04285359\n",
      "tran:, 6, -9.88, 0.00453034\n",
      "tran:, 7, 10.28, 0.00097440\n",
      "tran:, 8, -12.80, 0.00093868\n",
      "tran:, 9, 13.21, 0.38618480\n",
      "tran:, 10, -15.75, 0.00292543\n",
      "tran:, 11, 16.17, 0.00214517\n",
      "tran:, 12, -18.75, 0.00095546\n",
      "tran:, 13, 19.18, 0.38260810\n",
      "tran:, 14, -21.81, 0.00198520\n",
      "tran:, 15, 22.24, 0.00107209\n",
      "tran:, 16, -24.93, 0.00098420\n",
      "tran:, 17, 25.37, 0.04148387\n",
      "tran:, 18, -28.13, 0.00137337\n",
      "tran:, 19, 28.59, 0.00113879\n",
      "tran:, 20, -31.43, 0.00101836\n",
      "tran:, 21, 31.90, 0.01425444\n",
      "tran:, 22, -34.85, 0.00094045\n",
      "tran:, 23, 35.35, 0.00121503\n",
      "tran:, 24, -38.43, 0.00105001\n",
      "tran:, 25, 38.94, 0.00671291\n",
      "tran:, 26, -42.18, 0.00061482\n",
      "tran:, 27, 42.73, 0.00129876\n",
      "tran:, 28, -46.18, 0.00106158\n",
      "tran:, 29, 46.76, 0.00353975\n",
      "tran:, 30, -50.49, 0.00037006\n",
      "tran:, 31, 51.12, 0.00137924\n",
      "tran:, 32, -55.24, 0.00100766\n",
      "tran:, 33, 55.94, 0.00185785\n",
      "tran:, 34, -60.63, 0.00022601\n",
      "tran:, 35, 61.46, 0.00134905\n",
      "tran:, 36, -67.15, 0.00077375\n",
      "tran:, 37, 68.20, 0.00106048\n",
      "tran:, 38, -76.29, 0.00022316\n",
      "mode-coeff:, 0.061048, 0.937868, 0.998915\n",
      "poynting-flux:, 0.061102, 0.938344, 0.999447\n"
     ]
    }
   ],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "import meep as mp\n",
    "import math\n",
    "import cmath\n",
    "import numpy as np\n",
    "\n",
    "resolution = 50        # pixels/μm\n",
    "\n",
    "dpml = 1.0             # PML thickness\n",
    "dsub = 3.0             # substrate thickness\n",
    "dpad = 3.0             # length of padding between grating and PML\n",
    "gp = 10.0              # grating period\n",
    "gh = 0.5               # grating height\n",
    "gdc = 0.5              # grating duty cycle\n",
    "\n",
    "sx = dpml+dsub+gh+dpad+dpml\n",
    "sy = gp\n",
    "\n",
    "cell_size = mp.Vector3(sx,sy,0)\n",
    "pml_layers = [mp.PML(thickness=dpml,direction=mp.X)] \n",
    "\n",
    "wvl = 0.5              # center wavelength\n",
    "fcen = 1/wvl           # center frequency\n",
    "df = 0.05*fcen         # frequency width\n",
    "\n",
    "ng = 1.5\n",
    "glass = mp.Medium(index=ng)\n",
    "\n",
    "use_cw_solver = False  # CW solver or time stepping?\n",
    "tol = 1e-6             # CW solver tolerance\n",
    "max_iters = 2000       # CW solver max iterations\n",
    "L = 10                 # CW solver L\n",
    "\n",
    "# rotation angle of incident planewave; counter clockwise (CCW) about Z axis, 0 degrees along +X axis\n",
    "theta_in = math.radians(10.7)\n",
    "\n",
    "# k (in source medium) with correct length (plane of incidence: XY)\n",
    "k = mp.Vector3(fcen*ng).rotate(mp.Vector3(z=1), theta_in)\n",
    "\n",
    "symmetries = []\n",
    "eig_parity = mp.ODD_Z\n",
    "if theta_in == 0:\n",
    "  k = mp.Vector3(0,0,0)\n",
    "  symmetries = [mp.Mirror(mp.Y)]\n",
    "  eig_parity += mp.EVEN_Y\n",
    "\n",
    "def pw_amp(k,x0):\n",
    "  def _pw_amp(x):\n",
    "    return cmath.exp(1j*2*math.pi*k.dot(x+x0))\n",
    "  return _pw_amp\n",
    "\n",
    "src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)\n",
    "sources = [mp.Source(mp.ContinuousSource(fcen,fwidth=df) if use_cw_solver else mp.GaussianSource(fcen,fwidth=df),\n",
    "                     component=mp.Ez,\n",
    "                     center=src_pt,\n",
    "                     size=mp.Vector3(0,sy,0),\n",
    "                     amp_func=pw_amp(k,src_pt))]\n",
    "\n",
    "sim = mp.Simulation(resolution=resolution,\n",
    "                    cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    k_point=k,\n",
    "                    default_material=glass,\n",
    "                    sources=sources,\n",
    "                    symmetries=symmetries)\n",
    "\n",
    "refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)\n",
    "refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))\n",
    "\n",
    "if use_cw_solver:\n",
    "  sim.init_sim()\n",
    "  sim.solve_cw(tol, max_iters, L)\n",
    "else:\n",
    "  sim.run(until_after_sources=100)\n",
    "\n",
    "input_flux = mp.get_fluxes(refl_flux)\n",
    "input_flux_data = sim.get_flux_data(refl_flux)\n",
    "\n",
    "sim.reset_meep()\n",
    "\n",
    "geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),\n",
    "            mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]\n",
    "\n",
    "sim = mp.Simulation(resolution=resolution,\n",
    "                    cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    geometry=geometry,\n",
    "                    k_point=k,\n",
    "                    sources=sources,\n",
    "                    symmetries=symmetries)\n",
    "\n",
    "refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))\n",
    "sim.load_minus_flux_data(refl_flux,input_flux_data)\n",
    "\n",
    "tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)\n",
    "tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))\n",
    "\n",
    "if use_cw_solver:\n",
    "  sim.init_sim()\n",
    "  sim.solve_cw(tol, max_iters, L)\n",
    "else:\n",
    "  sim.run(until_after_sources=200)\n",
    "\n",
    "nm_r = np.floor((fcen*ng-k.y)*gp)-np.ceil((-fcen*ng-k.y)*gp) # number of reflected orders\n",
    "if theta_in == 0:\n",
    "  nm_r = nm_r/2 # since eig_parity removes degeneracy in y-direction\n",
    "nm_r = int(nm_r)\n",
    "\n",
    "res = sim.get_eigenmode_coefficients(refl_flux, range(1,nm_r+1), eig_parity=eig_parity)\n",
    "r_coeffs = res.alpha\n",
    "\n",
    "Rsum = 0\n",
    "for nm in range(nm_r):\n",
    "  r_kdom = res.kdom[nm]\n",
    "  Rmode = abs(r_coeffs[nm,0,1])**2/input_flux[0]\n",
    "  r_angle = np.sign(r_kdom.y)*math.acos(r_kdom.x/(ng*fcen))\n",
    "  print(\"refl:, {}, {:.2f}, {:.8f}\".format(nm,math.degrees(r_angle),Rmode))\n",
    "  Rsum += Rmode\n",
    "\n",
    "nm_t = np.floor((fcen-k.y)*gp)-np.ceil((-fcen-k.y)*gp)       # number of transmitted orders\n",
    "if theta_in == 0:\n",
    "  nm_t = nm_t/2 # since eig_parity removes degeneracy in y-direction\n",
    "nm_t = int(nm_t)\n",
    "\n",
    "res = sim.get_eigenmode_coefficients(tran_flux, range(1,nm_t+1), eig_parity=eig_parity)\n",
    "t_coeffs = res.alpha\n",
    "\n",
    "Tsum = 0\n",
    "for nm in range(nm_t):\n",
    "  t_kdom = res.kdom[nm]\n",
    "  Tmode = abs(t_coeffs[nm,0,0])**2/input_flux[0]\n",
    "  t_angle = np.sign(t_kdom.y)*math.acos(t_kdom.x/fcen)\n",
    "  print(\"tran:, {}, {:.2f}, {:.8f}\".format(nm,math.degrees(t_angle),Tmode))\n",
    "  Tsum += Tmode\n",
    "\n",
    "print(\"mode-coeff:, {:.6f}, {:.6f}, {:.6f}\".format(Rsum,Tsum,Rsum+Tsum))\n",
    "\n",
    "r_flux = mp.get_fluxes(refl_flux)\n",
    "t_flux = mp.get_fluxes(tran_flux)\n",
    "Rflux = -r_flux[0]/input_flux[0]\n",
    "Tflux =  t_flux[0]/input_flux[0]\n",
    "print(\"poynting-flux:, {:.6f}, {:.6f}, {:.6f}\".format(Rflux,Tflux,Rflux+Tflux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}