{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Focusing Properties of a Binary-Phase Zone Plate\n",
"\n",
"It is also possible to compute a [near-to-far field transformation](https://meep.readthedocs.io/en/latest/Python_User_Interface/#near-to-far-field-spectra) in cylindrical coordinates. This is demonstrated in this example for a binary-phase [zone plate](https://en.wikipedia.org/wiki/Zone_plate) which is a rotationally-symmetric diffractive lens used to focus a normally-incident planewave to a single spot.\n",
"\n",
"Using [scalar theory](http://zoneplate.lbl.gov/theory), the radius of the $n$th zone can be computed as:\n",
"\n",
"
\n",
"\n",
"where $n$ is the zone index (1,2,3,...,$N$), $f$ is the focal length, and $\\lambda$ is the operating wavelength. The main design variable is the number of zones $N$. The design specifications of the zone plate are similar to the binary-phase grating in [Tutorial/Mode Decomposition/Diffraction Spectrum of a Binary Grating](https://meep.readthedocs.io/en/latest/Python_Tutorials/Mode_Decomposition/#diffraction-spectrum-of-a-binary-grating) with refractive index of 1.5 (glass), $\\lambda$ of 0.5 μm, and height of 0.5 μm. The focusing property of the zone plate is verified by the concentration of the electric-field energy density at the focal length of 0.2 mm (which lies *outside* the cell). The planewave is incident from within a glass substrate and spans the entire length of the cell in the radial direction. The cell is surrounded on all sides by PML. A schematic of the simulation geometry for a design with 25 zones and flat-surface termination is shown below. The near-field line monitor is positioned at the edge of the PML.\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using MPI version 3.1, 1 processes\n"
]
}
],
"source": [
"import math\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import meep as mp\n",
"import numpy as np\n",
"\n",
"\n",
"resolution_um = 25\n",
"\n",
"pml_um = 1.0\n",
"substrate_um = 2.0\n",
"padding_um = 2.0\n",
"height_um = 0.5\n",
"focal_length_um = 200\n",
"scan_length_z_um = 100\n",
"farfield_resolution_um = 10\n",
"\n",
"pml_layers = [mp.PML(thickness=pml_um)]\n",
"\n",
"wavelength_um = 0.5\n",
"frequency = 1 / wavelength_um\n",
"frequench_width = 0.2 * frequency\n",
"\n",
"# The number of zones in the zone plate. \n",
"# Odd-numbered zones impart a π phase shift and \n",
"# even-numbered zones impart no phase shift. \n",
"num_zones = 25\n",
"\n",
"# Specify the radius of each zone using the equation \n",
"# from https://en.wikipedia.org/wiki/Zone_plate. \n",
"zone_radius_um = np.zeros(num_zones)\n",
"for n in range(1, num_zones + 1):\n",
" zone_radius_um[n - 1] = math.sqrt(\n",
" n * wavelength_um * (focal_length_um + n * wavelength_um / 4)\n",
" )\n",
"\n",
"size_r_um = zone_radius_um[-1] + padding_um + pml_um\n",
"size_z_um = pml_um + substrate_um + height_um + padding_um + pml_um\n",
"cell_size = mp.Vector3(size_r_um, 0, size_z_um)\n",
"\n",
"# Specify a (linearly polarized) planewave at normal incidence. \n",
"sources = [\n",
" mp.Source(\n",
" mp.GaussianSource(frequency, fwidth=frequench_width, is_integrated=True),\n",
" component=mp.Er,\n",
" center=mp.Vector3(0.5 * size_r_um, 0, -0.5 * size_z_um + pml_um),\n",
" size=mp.Vector3(size_r_um),\n",
" ),\n",
" mp.Source(\n",
" mp.GaussianSource(frequency, fwidth=frequench_width, is_integrated=True),\n",
" component=mp.Ep,\n",
" center=mp.Vector3(0.5 * size_r_um, 0, -0.5 * size_z_um + pml_um),\n",
" size=mp.Vector3(size_r_um),\n",
"\tamplitude=-1j,\n",
" ),\n",
"]\n",
"\n",
"glass = mp.Medium(index=1.5)\n",
"\n",
"# Add the substrate. \n",
"geometry = [\n",
" mp.Block(\n",
" material=glass,\n",
" size=mp.Vector3(size_r_um, 0, pml_um + substrate_um),\n",
" center=mp.Vector3(\n",
" 0.5 * size_r_um, 0, -0.5 * size_z_um + 0.5 * (pml_um + substrate_um)\n",
" ),\n",
" )\n",
"]\n",
"\n",
"# Add the zone plates starting with the ones with largest radius. \n",
"for n in range(num_zones - 1, -1, -1):\n",
" geometry.append(\n",
" mp.Block(\n",
" material=glass if n % 2 == 0 else mp.vacuum,\n",
" size=mp.Vector3(zone_radius_um[n], 0, height_um),\n",
" center=mp.Vector3(\n",
" 0.5 * zone_radius_um[n],\n",
" 0,\n",
" -0.5 * size_z_um + pml_um + substrate_um + 0.5 * height_um,\n",
" ),\n",
" )\n",
" )\n",
"\n",
"sim = mp.Simulation(\n",
" cell_size=cell_size,\n",
" boundary_layers=pml_layers,\n",
" resolution=resolution_um,\n",
" sources=sources,\n",
" geometry=geometry,\n",
" dimensions=mp.CYLINDRICAL,\n",
" m=-1,\n",
")\n",
"\n",
"# Add the near-field monitor (must be entirely in air). \n",
"n2f_monitor = sim.add_near2far(\n",
" frequency,\n",
" 0,\n",
" 1,\n",
" mp.Near2FarRegion(\n",
" center=mp.Vector3(0.5 * (size_r_um - pml_um), 0, 0.5 * size_z_um - pml_um),\n",
" size=mp.Vector3(size_r_um - pml_um, 0, 0),\n",
" ),\n",
" mp.Near2FarRegion(\n",
" center=mp.Vector3(\n",
" size_r_um - pml_um,\n",
" 0,\n",
" 0.5 * size_z_um - pml_um - 0.5 * (height_um + padding_um),\n",
" ),\n",
" size=mp.Vector3(0, 0, height_um + padding_um),\n",
" ),\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" block, center = (26.6946,0,-1.75)\n",
" size (53.3891,0,3)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (25.1946,0,0)\n",
" size (50.3891,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (24.6779,0,0)\n",
" size (49.3559,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (24.1509,0,0)\n",
" size (48.3018,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (23.6128,0,0)\n",
" size (47.2255,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (23.0628,0,0)\n",
" size (46.1255,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (22.5,0,0)\n",
" size (45,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (21.9235,0,0)\n",
" size (43.847,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (21.3322,0,0)\n",
" size (42.6644,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (20.7248,0,0)\n",
" size (41.4495,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" ...(+ 16 objects not shown)...\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAABWCAYAAADR9dHwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAOl0lEQVR4nO2df4wc51nHP8+u73zG8cVuktqXxL7zofiSNCpuzwqpQpAplNppSUCqpVaiQPkjgIIUJApqaolcEREIpIpKrZombUQt2qJGaUkaFbmh2AVES+06P9zUvfhH1sGJYyvBF87u+e529+GPmbms92Zm5867O3e734/0ambeed/v+zwzs/vs+87OO+buCCGEEEkU8jZACCHE0kaBQgghRCoKFEIIIVJRoBBCCJGKAoUQQohUFCiEEEKksiJvA1pB8YqiV/ureZshhBDLhqtmruL111+32J3unksCNgL7gCPAC8B9MWW2A28Cz4bpLzJpD+CN2PfSPr/6b6/2fS/ta1h2MUhf+tKX/lLUZyz++3F0dNQ96Ts1aUerEzAAvDtcXwO8CNxcV2Y78NSCtRsEik6+CKQvfelLP01/WQWKeYbAE8D76vKaHijyPknSl770pZ+n/rINFMAQ8DLQX5e/HXgDeA74F+AdmfQSAsVSOEnSl770pZ+n/mIChbnnO9eTmV0BfA940N2/UbevH6i6+3kzuxP4jLvfkKBzD3BPsDU6CgdbarcQQixX4r72t23bxsGDB2NvZucaKMysB3gK2Ovun85QvgRsc/fX08ttcwUKIYRIYGx+PBj91mhioMjt77FmZsCXgCNJQcLMNgBn3N3N7FaC5z7eaKQ9OgoHFSeEEOISLAwD/sD8DsK2b21LrJfncxS3Ax8FDpvZs2HeJ4FNAO7+EPAh4I/MrAxMAR/2vMfKhBCiy8gtULj7fwLxD3e8VeazwGfbY5EQQog4NIWHEEKIVBQohBBCpKJAIYQQIhUFCiGE6HL2l/an7legEEKILmZ/aT+7HtuVWkaBQgghupQoSDy267HUch35PgqAarX73kdhZpil/uO4aTTz+DbL7sXatNj2/a05yTJRKCzsd1lWfxZifxbNrHpZ/M/icyObGmmk2ZHFl7T209pOa3cx9RrZmmTnQq+riNogsX1oe2rZXAOFme0APgMUgS+6+9/U7bdw/53Az4Dfc/dDWbQXe/BENpbi8W23Ta0OzK3wp5mazfL/cm26XDsW2/5i211svWaeu4UECWgQKMzs/e6+N2HfLndP76+kaxeBzwHvA04BB8zsSXf/SU2xncANYfpF4PPhMpWpqSkOHz68WNOWJZVKhbVr1zI0NNTytqampjh69GjTegF9fX2MjIxcls7k5CTHjx+nWCwuqF65XObaa69l/fr1C6rn7oyPjzM9PZ3pA1ypVBgaGmLt2rWZ9GdmZhgfHw9m7kw5zpVKhXXr1jE4OJjJhiNHjlCtVhM1K5UK/f39DA8Pp2pVq1XGx8eZmZmJ9d/dKRQKjIyM0NPTk6hz7tw5Tp48GXveol/dIyMjrFy5MlHj2LFjXLhwYZ4d1WqV3t5eRkZGEs/R+fPnOX78+Lz9kf033XRT4jVVKpWYmJiYt79SqTA8PEx/f39svfHxcS5evHhJm9VqlVWrVrFly5bYOhMTE5RKpXltuTtbtmyhr68vtl4SCwkScw0lJaBC8Ba662L2HUqr2ygB7yGYDDDavh+4v67MF4CP1GyPAwMZtL0b044dO9zdvVKpeLVaTZpleNGUy2V3dz9w4EBT7R4eHvbZ2Vl39wXbHdn0xBNPLLhdM3PAd+/e7e4+Z0MakX0XLlzwgYGBBbW3Z8+eS2yOo1KpuLv70aNHvVAoZNK96667UnUjzTNnzviaNWsa6t1xxx1z9erPR7Q9OTnpGzZsSNXp6enxEydOXGJDRHSsH3nkkYb2PP/88/M0au3aunVrYt0NGzb45OTkvDrRsdq7d29i3f7+fj979uwlbVer1bn1nTt3JtZ9/PHHL2knant6etqHhoZi69x4441z5aJlVH/Pnj2JbR06dCj1/NcCQYqbijxtmvFGP4WeB74K/MDM6m+LX+7PyeuA/6nZPhXmLbRMYIzZPWZ20My6bjrA6NdhtCwUCi0ZEol+zTR7SKRQKMxpLtTuqN5ibKo9Xlk1autkbTPu/CRRa0sj/ay6tZppxzfLdRRtZ+m5mVnisa2/lrLYVatRW77RNRPXRpZzHmd/Ul6SrfVlFntO6/Pj/FtIbzpzTyKk0T0Kd/dHzOx7wFfCd0Lc6+4/I4hml0Pc2a3XzFImyHR/GHgYIPyl2DWY2SU3yGZnZ1syvl2tVunp6aFSqTRFL7IbgiGgnp6eBd+QrlQq9Pb2ztlUq9mIQqFAtVqda7NcLjf8sLk7K1asoFwuX/LhTWsz7vwktVN7jBt9Ada3Wy6XE4eB6m3OqldfPtKanZ1t6L+ZUS6XY7Wi8xYd+7RjGJWJu65rj2OtRv21VW9D/bUc136S/VHdONuj9Ug3OtfuTrFYnNOrLVu7jM57pBcdpzi/ao9l0vGZT/GSOllJfR+FmR1y93eH6yuAvwJ+C/gd4PPRvsVgZu8Bxtz9/eH2/QDu/tc1Zb4A7Hf3r4Xb48B2dz/dQLurAkXEqlWr2LRpU+oY9OXg4bjt1NQUL7/8ctN0e3p6GBoaWvQ/j4rFIm+++SanT6deFomsW7eO9evXZ/pyjqhWq5RKpUs++I3YsGEDa9euTW0nOsbT09OUSqVMuqtXr2bjxo2JupHmzMwMpVKpYSDt6+tjcHAwtVy1WuXEiRMNA/vQ0BArV66cd01G5+3cuXOcOXMmVWNwcJC+vr7E6/rkyZNMT0/H1i0UCgwPD8fegygWi0xOTvLKK6/E1jUzNm/ezIoVK+bVLRQKnDp1igsXLsTWHRgY4Morr5x3TtydUqnE7OzsvDq9vb3z7jFGdk5MTPDaa6/FttXo+NTy4ovjAGzZMv+e4MmTJ7l48WLyhZmUgGdi8rYDJ4DJtLqNEkFv5gSwGegleN3pO+rKfIDgFagG3Ab8MKP2gserlZSUlLoijSXvS/pObTT09Kn6DHffb2ajwB80qJuKu5fN7I+BvQT9oUfd/QUz+8Nw/0PAtwn+GnuM4O+xH7ucNruBhQy9LEX9xdAMm9rlV5Z2FmNL1uGvZukthLShqaxttOL8ZNVMGpZqxXlcrE7LP/dL7UPfDLp16EkIIRoyFqYY3D126GnpPTUlhBBiSaFAIYQQIhUFCiGEEKkoUAghhEhFgUIIIUQqChRCCCFSyWWacTP7O+A3gBngOPAxd5+IKVcCJgkmJyy7+7Y2mimEEIL8ehRPA7e4+zuBFwlmjk3iV9x9q4KEEELkQy49Cnf/Ts3mD4APNVN/FOi6KWSFECILY/HZabNELYV7FL9PMJ9THA58x8x+ZGb3tNEmIYQQEZczsV+Difn+FfhxTLq7psxu4JuEU4nEaFwbLt9OMGngL6e0dw9BR+IgV7ZoMq0hnD8Ll9KXvvSlvxz1x5L3JX6/tipQZAgkvwt8H/i5jOXHgI9nKjuwhE+S9KUvfennqT8Wk3dbsFxSgQLYAfwEuCalzGpgTc36fwE7cgkUy+kikL70pS/9tDRWt30bzgPB+lILFMcIXnH6bJgeCvOvBb4drg8TDDc9B7wA7M6s38xAsdwuAulLX/rST0tjNetRkFiKPYqWB6JmBYrleBFIX/rSl35aGguXdUECkgNFp76PYhIYz9uOHLgaeD1vI3JCvncf3eo3tMb3QXe/Jm5HLs9RtIFx78IH9MzsYDf6DfK9G33vVr+h/b4vhecohBBCLGEUKIQQQqTSqYHi4bwNyIlu9RvkezfSrX5Dm33vyJvZQgghmken9iiEEEI0iY4KFGa2w8zGzeyYmX0ib3taiZk9amZnzezHNXlvM7OnzexouFyXp42twMw2mtk+MztiZi+Y2X1hfjf43mdmPzSz50LfPxXmd7zvAGZWNLNnzOypcLsr/Ibg3TxmdtjMnjWzg2Fe2/zvmEBhZkXgc8BO4GbgI2Z2c75WtZR/IJgKpZZPAN919xuA74bbnUYZ+FN3vwm4Dbg3PM/d4Ps08F53/wVgK7DDzG6jO3wHuA84UrPdLX5H1L+bp23+d0ygAG4Fjrn7CXefAf4JuDtnm1qGu/878L912XcDXw7Xvwz8ZjttagfuftrdD4XrkwRfHNfRHb67u58PN3vC5HSB72Z2PfAB4Is12R3vdwPa5n8nBYrrCOaPijgV5nUT6939NARfqATTs3csZjYEvAv4b7rE93D45VngLPC0u3eL738P/DlQrcnrBr8jnPnv5mmb/530ZHbcC5r0l64OxcyuAB4H/sTd/88s7f1cnYO7V4CtZrYW+KaZ3ZKzSS3HzD4InHX3H5nZ9pzNyYvb3f1VM3s78LSZ/bSdjXdSj+IUsLFm+3rg1ZxsyYszZjYAEC7P5mxPSzCzHoIg8RV3/0aY3RW+R7j7BLCf4D5Vp/t+O3CXmZUIhpTfa2b/SOf7PYe7vxouzxK87O1W2uh/JwWKA8ANZrbZzHqBDwNP5mxTu3mS4IVQhMsncrSlJVjQdfgScMTdP12zqxt8vybsSWBmq4BfA35Kh/vu7ve7+/XuPkTwuf43d/9tOtzvCDNbbWZronXg1wneFto2/zvqgTszu5NgLLMIPOruD+ZrUesws68B2wlmkTwDPAD8M/B1YBPwMrDL3etveC9rzOyXgP8ADvPWePUnCe5TdLrv7yS4aVkk+JH3dXf/SzO7ig73PSIcevq4u3+wW/w2s2GCXgQEtwu+6u4PttP/jgoUQgghmk8nDT0JIYRoAQoUQgghUlGgEEIIkYoChRBCiFQUKIQQQqTSSU9mC7FkMbMKwV96VwAvAR8NH5oTYsmjHoUQ7WEqnPnzFoLJHO/N2yAhsqJAIUT7+T7dN2GlWMYoUAjRRsL3pvwq3Te9jFjGKFAI0R5WhdODvwG8DXg6X3OEyI4ChRDtYcrdtwKDQC+6RyGWEZrrSYg2YGbn3f2KcP1dBDN9/ry7z+ZrmRCNUY9CiDbj7s8AzxFMmS3Ekkc9CiGEEKmoRyGEECIVBQohhBCpKFAIIYRIRYFCCCFEKgoUQgghUlGgEEIIkYoChRBCiFQUKIQQQqTy/wovnEoJdLeGAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"sim.plot2D(ax=ax)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"-----------\n",
"Initializing structure...\n",
"time for choose_chunkdivision = 0.00028924 s\n",
"Working in Cylindrical dimensions.\n",
"Computational cell is 53.4 x 0 x 6.52 with resolution 25\n",
" block, center = (26.6946,0,-1.75)\n",
" size (53.3891,0,3)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (25.1946,0,0)\n",
" size (50.3891,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (24.6779,0,0)\n",
" size (49.3559,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (24.1509,0,0)\n",
" size (48.3018,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (23.6128,0,0)\n",
" size (47.2255,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (23.0628,0,0)\n",
" size (46.1255,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (22.5,0,0)\n",
" size (45,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (21.9235,0,0)\n",
" size (43.847,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" block, center = (21.3322,0,0)\n",
" size (42.6644,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (1,1,1)\n",
" block, center = (20.7248,0,0)\n",
" size (41.4495,0,0.5)\n",
" axes (1,0,0), (0,1,0), (0,0,1)\n",
" dielectric constant epsilon diagonal = (2.25,2.25,2.25)\n",
" ...(+ 16 objects not shown)...\n",
"time for set_epsilon = 0.536792 s\n",
"-----------\n",
"Meep: using complex fields.\n",
"on time step 309 (time=6.18), 0.0129541 s/step\n",
"on time step 630 (time=12.6), 0.012487 s/step\n",
"on time step 950 (time=19), 0.0125075 s/step\n",
"on time step 1270 (time=25.4), 0.0125189 s/step\n",
"on time step 1592 (time=31.84), 0.0124597 s/step\n",
"on time step 1914 (time=38.28), 0.0124481 s/step\n",
"on time step 2241 (time=44.82), 0.0122668 s/step\n",
"field decay(t = 50.02): 0.11926764409048278 / 0.11926764409048278 = 1.0\n",
"on time step 2565 (time=51.3), 0.0123561 s/step\n",
"on time step 2891 (time=57.82), 0.0122783 s/step\n",
"on time step 3211 (time=64.22), 0.012518 s/step\n",
"on time step 3537 (time=70.74), 0.0122997 s/step\n",
"on time step 3862 (time=77.24), 0.012323 s/step\n",
"on time step 4188 (time=83.76), 0.0122724 s/step\n",
"on time step 4514 (time=90.28), 0.0122736 s/step\n",
"on time step 4839 (time=96.78), 0.0123397 s/step\n",
"field decay(t = 100.04): 6.406243098076629e-08 / 0.11926764409048278 = 5.371316878881679e-07\n",
"run 0 finished at t = 100.04 (5002 timesteps)\n"
]
}
],
"source": [
"# Timestep the fields until they have sufficiently decayed away. \n",
"sim.run(\n",
" until_after_sources=mp.stop_when_fields_decayed(\n",
" 50.0, mp.Er, mp.Vector3(0.5 * size_r_um, 0, 0), 1e-6\n",
" )\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"get_farfields_array working on point 44 of 523 (8% done), 0.0923484 s/point\n",
"get_farfields_array working on point 75 of 523 (14% done), 0.133575 s/point\n",
"get_farfields_array working on point 100 of 523 (19% done), 0.16357 s/point\n",
"get_farfields_array working on point 120 of 523 (22% done), 0.207044 s/point\n",
"get_farfields_array working on point 138 of 523 (26% done), 0.229834 s/point\n",
"get_farfields_array working on point 155 of 523 (29% done), 0.246436 s/point\n",
"get_farfields_array working on point 171 of 523 (32% done), 0.258188 s/point\n",
"get_farfields_array working on point 186 of 523 (35% done), 0.267175 s/point\n",
"get_farfields_array working on point 201 of 523 (38% done), 0.279909 s/point\n",
"get_farfields_array working on point 214 of 523 (40% done), 0.312383 s/point\n",
"get_farfields_array working on point 226 of 523 (43% done), 0.336642 s/point\n",
"get_farfields_array working on point 238 of 523 (45% done), 0.360156 s/point\n",
"get_farfields_array working on point 249 of 523 (47% done), 0.375252 s/point\n",
"get_farfields_array working on point 260 of 523 (49% done), 0.390756 s/point\n",
"get_farfields_array working on point 270 of 523 (51% done), 0.403964 s/point\n",
"get_farfields_array working on point 280 of 523 (53% done), 0.418169 s/point\n",
"get_farfields_array working on point 290 of 523 (55% done), 0.427388 s/point\n",
"get_farfields_array working on point 300 of 523 (57% done), 0.437372 s/point\n",
"get_farfields_array working on point 309 of 523 (59% done), 0.44643 s/point\n",
"get_farfields_array working on point 318 of 523 (60% done), 0.453257 s/point\n",
"get_farfields_array working on point 327 of 523 (62% done), 0.464375 s/point\n",
"get_farfields_array working on point 336 of 523 (64% done), 0.471166 s/point\n",
"get_farfields_array working on point 345 of 523 (65% done), 0.475754 s/point\n",
"get_farfields_array working on point 354 of 523 (67% done), 0.481679 s/point\n",
"get_farfields_array working on point 363 of 523 (69% done), 0.487806 s/point\n",
"get_farfields_array working on point 372 of 523 (71% done), 0.494046 s/point\n",
"get_farfields_array working on point 381 of 523 (72% done), 0.498369 s/point\n",
"get_farfields_array working on point 389 of 523 (74% done), 0.508082 s/point\n",
"get_farfields_array working on point 397 of 523 (75% done), 0.507507 s/point\n",
"get_farfields_array working on point 405 of 523 (77% done), 0.514216 s/point\n",
"get_farfields_array working on point 413 of 523 (78% done), 0.517538 s/point\n",
"get_farfields_array working on point 421 of 523 (80% done), 0.519612 s/point\n",
"get_farfields_array working on point 429 of 523 (82% done), 0.522912 s/point\n",
"get_farfields_array working on point 437 of 523 (83% done), 0.527134 s/point\n",
"get_farfields_array working on point 445 of 523 (85% done), 0.546583 s/point\n",
"get_farfields_array working on point 452 of 523 (86% done), 0.575292 s/point\n",
"get_farfields_array working on point 459 of 523 (87% done), 0.590524 s/point\n",
"get_farfields_array working on point 466 of 523 (89% done), 0.60685 s/point\n",
"get_farfields_array working on point 473 of 523 (90% done), 0.617983 s/point\n",
"get_farfields_array working on point 480 of 523 (91% done), 0.635441 s/point\n",
"get_farfields_array working on point 487 of 523 (93% done), 0.63931 s/point\n",
"get_farfields_array working on point 494 of 523 (94% done), 0.655253 s/point\n",
"get_farfields_array working on point 501 of 523 (95% done), 0.663373 s/point\n",
"get_farfields_array working on point 507 of 523 (96% done), 0.673751 s/point\n",
"get_farfields_array working on point 513 of 523 (98% done), 0.681454 s/point\n",
"get_farfields_array working on point 519 of 523 (99% done), 0.690325 s/point\n",
"get_farfields_array working on point 46 of 1000 (4% done), 0.0886262 s/point\n",
"get_farfields_array working on point 91 of 1000 (9% done), 0.0901391 s/point\n",
"get_farfields_array working on point 136 of 1000 (13% done), 0.090729 s/point\n",
"get_farfields_array working on point 181 of 1000 (18% done), 0.0905702 s/point\n",
"get_farfields_array working on point 225 of 1000 (22% done), 0.0910369 s/point\n",
"get_farfields_array working on point 270 of 1000 (27% done), 0.090324 s/point\n",
"get_farfields_array working on point 315 of 1000 (31% done), 0.0903349 s/point\n",
"get_farfields_array working on point 360 of 1000 (36% done), 0.0901996 s/point\n",
"get_farfields_array working on point 405 of 1000 (40% done), 0.0904465 s/point\n",
"get_farfields_array working on point 450 of 1000 (45% done), 0.0904624 s/point\n",
"get_farfields_array working on point 495 of 1000 (49% done), 0.0901669 s/point\n",
"get_farfields_array working on point 540 of 1000 (54% done), 0.0907005 s/point\n",
"get_farfields_array working on point 585 of 1000 (58% done), 0.0901638 s/point\n",
"get_farfields_array working on point 630 of 1000 (63% done), 0.0902798 s/point\n",
"get_farfields_array working on point 675 of 1000 (67% done), 0.0905745 s/point\n",
"get_farfields_array working on point 720 of 1000 (72% done), 0.0904627 s/point\n",
"get_farfields_array working on point 765 of 1000 (76% done), 0.0908991 s/point\n",
"get_farfields_array working on point 810 of 1000 (81% done), 0.090244 s/point\n",
"get_farfields_array working on point 855 of 1000 (85% done), 0.0902669 s/point\n",
"get_farfields_array working on point 900 of 1000 (90% done), 0.090408 s/point\n",
"get_farfields_array working on point 945 of 1000 (94% done), 0.0902008 s/point\n",
"get_farfields_array working on point 990 of 1000 (99% done), 0.0900602 s/point\n"
]
}
],
"source": [
"farfields_r = sim.get_farfields(\n",
" n2f_monitor,\n",
" farfield_resolution_um,\n",
" center=mp.Vector3(\n",
" 0.5 * (size_r_um - pml_um),\n",
" 0,\n",
" -0.5 * size_z_um + pml_um + substrate_um + height_um + focal_length_um,\n",
" ),\n",
" size=mp.Vector3(size_r_um - pml_um, 0, 0),\n",
")\n",
"\n",
"farfields_z = sim.get_farfields(\n",
" n2f_monitor,\n",
" farfield_resolution_um,\n",
" center=mp.Vector3(\n",
" 0, 0, -0.5 * size_z_um + pml_um + substrate_um + height_um + focal_length_um\n",
" ),\n",
" size=mp.Vector3(0, 0, scan_length_z_um),\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 0.98, 'binary-phase zone plate with focal length $z$ = 200 μm')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEnCAYAAABsR64CAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABhsElEQVR4nO2deXhU1fn4P28ChDWshmqUxK2CigtYa1u1UFHRal1q/dYS1LpQoNalqK3FX4ttqdalxaVC0aLWpLXWqqh1a9SIdrEq7iItUkFRQUEFWUPy/v44d8hkcteZOzN3kvN5nvvM3OWc+86dc++5533f876iqlgsFovFkgtlxRbAYrFYLKWP7UwsFovFkjO2M7FYLBZLztjOxGKxWCw5YzsTi8ViseSM7UwsFovFkjO2M7FYLBZLztjOxGKxWCw50+U6ExF5S0TGeex7TUTGFFai6Pj9BoshadcoqG1FlVdE9hCRF0RknYicG4eMUWQq5vVN2n9rMXS5zsQPVd1LVZuKLYelsBTi4ZTZtmI458VAk6r2U9XrchYwoSS54xCRChH5nYgsczr1F0TkqIxjBonIPSKy3jnuW1H2lxK2M8kzItKt2DJYOiU1wGvFFqKL0w14G/gy0B/4f8CdIlKbdsxvgC3AUGACMFtE9oqwv3RQ1S61AG8BlwCvAx8BtwA90/aNS/t+IfAy8Anwp9Rxzv4fAm8C65y6Tsg4xw+cspuBi4C/ZMhxPTArBzmzle0HwApn32LgsLR9OwB/AT4A/gec6yHb/wGfpi2bMW/JACOAJuBjzMPuaxm/yU/uUOcPuj4u/6Xf9bgdaAU2Or/l4gjX4dvA/WnrS4A709bfBvZzkcftnL7XJuO8jwMtwCan/GcDrvtOwN3O71kN3BChHY/zuf7j0tY9r1mI/30U8IIjw5+d/T93u05h6suQ807at1UFzsnj8+Vl4OvO9z6YjuKzGe3tijD7XepWYLe09VuBn2dc54scGdYDv8N0Ug8517YRGJi3356vipO6OBf8VecGGwT8PfWH0LEz+bdzkwwCFgGT0+r5hrOvDPNwXQ9sn1b2ReccvYDtnf0DnP3dgFXA6BzkjCwbsAfmAbeDc1wtsKvzvQx4Hvgx0APYBVgKHBlwPSud838H6I55oP7IqeMrTiPeI0juqOf3uz4u/6Xnf+VybGg5nH0fO2W2B5YBK9L2fQSUZZ7DZ93zP3U5dxNwlvPd87oD5cBLwK8xD6+ewMER2nFgZxJ0zQL+9x7OdTvP+R0nYh6wHe7JbK9VWrnvYzqtQR77H3D+T7flgRD1D8V08MOd9f2BjRnHXIjzAhK036X+MJ3Jvxw5qjHPmIXOeSowLyE/ifuZmlq6qprrBlV9W1XXADOBUzyOu05V33WOux/YL7VDVf/s7GtV1T8B/wUOzCj7tqpuVNX3gAWYGxdgPPChqj6fg5zZyNaCaVR7ikh3VX1LVd90in0O2E5Vf6qqW1R1KXAT8E0v4USkDPgDZlTyW+AgoC/mzWqLqj6OuUHDyB35/AHXZxsh/qt0Qsvh7Fvn/IYvA48AK0RkuLP+lKq2+sified/GoDfdT8Q89C9SFXXq+omVX067TdEuTZehLlmXr/tIMzL1XWq2qyqd2M6iiAiXSsROQ84FdMxrXE7RlWPUdUBHssxAfV3BxqA21T1DWdzX8zIKZ1PgH4h92fD9aq6UlVXAE8Bz6jqC6q6GbgH07Hkha6qz3877fsyzM3mxvtp3zekHycip2LedGqdTX2BIR7nALgNmIK5yeoww1lEZALwW+eYp1Q13YDnJ2dk2VR1iYicD8wA9hKRR4Dvq+q7GB38DiLycVq95ZgG6cVMTMNPeRPtALyd8QBdhnlLCpI7m/OH+h9D/FfpRJXjSWAMsJvz/WNMR/IFZz0Knv9pAH7XfSdgmapudSsY8dp4Eeaaef22HTCjufRcGJn3jhuhr5WInAOcCXxFVVeHqDsSzkvV7ZgR1Tlpuz7FjNzTqcS8gITZnw0r075vdFnvm0PdvnTVkclOad+HAe9GKSwiNZhO4RxgsKoOwKhcJO2wzEQx9wL7iMjewDGYtxhUtUFV+zrLURllIssZJJuq/kFVD8Y8ABT4pVP0beB/GW9j/VT1aI/zfBPz5nuSqjY7m98FdnJurnS5VwTJHfX8DoHXJ4v/Kqocqc7kEOf7k5jO5Mv4dyZxJhLyu+5vA8PcHEFCXpswZPPfpXgPqBaR9HOm/685XScRmYJRwR6mqh8GHPuQiHzqsTzkUUZos018Pe1eAPgP0E1Edk/bti9tjhNB+93omfZ9gN/vKTRdtTP5rojsKCKDMHrmP0Us3wfTyD8AEJFvA3v7FVDVTcBdGLXQv1V1eZ7k9JTNmZvwFRGpwOh2N2JUX2BUC2tF5Aci0ktEykVkbxH5XOYJRGR/jAPB8ar6QdquZzA694tFpLszr+JY4I4Qcoc+fxphrk+Y/2olRs+fjRxPAmOBXqr6DuZtfDwwGKOf9yL9nLnid93/jXlgXyEifUSkp4h8ySkXuR17kM1/l+KfmDZ4joh0E5HjaK9my/o6icgkTEc5LqOduqKqR6W92GUumS96KWZjnB+OVdWNGfWtxzg+/NS59l8CjsPRSgTt9+DbzvXdDzgM6Oeo2IpOV+1M/gA8ijESLsV4joRGVV8HrsHcCCuBkRgDcBC3Ocf6NZac5AyQrQK4AvgQoyaowjyEUdUWzANoP4w3zofAzRiXx0yOAwYCT6e/uanqFuBrwFFO+RuBU9N0yH5yRzl/isDrE/K/uhy41FHTXBBFDlX9D0Zd8ZSzvtaR5e/Ob/Ji2zlF5EKf4wLxu+5p13U3YDnwDsbQnks7zjx/Nv9duuwnYtRQH2NUwA9gPAQht+t0JbAr8GZaO50YsQ5PnJHddzC/+/20c0xIO2wqxglnFfBHYIqqvhZhfya9MS8HN2EcHk7DOFwUHWmvqrTkExEZBrwBfMZ56Pgd+xbGW6exELKVGvb6dF5E5BlgjqreUmxZkoSIKLC7qi4ptixudNWRScFx9NnfB+4I6kgslq6EiHxZRD7jqLlOA/YBHi62XJZodFVvroIiIn0waoRlGH26xWJpYw/M5MK+mAmUJ6lxp7eUEFbNZbFYLJacsWoui8ViseSM7UwsFovFkjO2M7FYLBZLztjOxGKxWCw5YzsTi8ViseSM7UwsFovFkjO2M7FYLBZLztjOxGKxWCw5YzsTi8ViseSM7UwsFovFkjO2M7FYLBZLztjOxGKxWCw502WiBg8ZMkRra2uLLYalE/P8889/qKrbFfq8tm1b8knYdt3pOxMRORY4trq6mquvvrrY4lg6MWPHjl1WyPPZtm0pBGHbdafvTCyWzoqq3g/cf8ABB5w9ZsyYYotj6eJ0+s7E3nCWzkr6yKSpqanY4li6OJ2+M7E3nKWzYl+ULEmi03cm9oazWCyW/NPpOxOLpbNiR92WJNHpO5Nsb7jGxipuvnkXVq2qoKpqM2edtZRx41blTU6LJSp21B2e3r1h48aO20Xg9tthwoTCy9TZEFUttgwF4YADDtDnnnsu1LENDTBpEmzY0LZNBCZPhhtvzJOAlpJHRJ5X1QMKfd4obburUV0N774bfFyvXu3vd0sbYdu1HZm4MG3aQWzY0LPdNlWYPVt5990VnH/+kvgFtVgiYtVc/owdeyggzuLPxo2KCEyfvshqILLEjkxcKCsznYcX9fV2WGzpiB2ZJAcJ7j88sfd3e8K2axuby4Vhw/z3n3deYeSwWCzRKS/PrXxdXTxydDVsZ+LCzJn+bzarVxdOFovFEp7qamhtzb2eXEY2XRVrM3Ghuhq+9rXdmD+/Gnd9q9LU9GSMUlos0bE2k/Y0Nlbx7rsjCLaRaKhjams/5dZbn49HuC6AtZn40K8ffPqp+z6rV7VkYm0mxaVbN2hpibfOLvJ49MXaTGJgzhzvfdOnF04Oi8XiT0NDuI7ksMNMB1FfH67e6urc5OpK2M7EB7+Rx7KCBhu3WCx+nH568DGHHQaNjeb7hAnhRh3vvms6Kksw1mYSwNChB7FyZU+XPcqll1qfdIul2DQ0wNat/sfssENbR5JOfX2w99bEiValHQZrMwmgocE0JrfLVFMDb72Vu2yWzkGhbSZpL0pn14fV23RCjjjiEJqb/fyBlSee8HaYOemkz7N6dU+8jfLapSczjh07NlS7tp1JCLzcBEXicUO0dA6sAb44BLnxpqu3sq2jRw/YvDmaXJ0Fa4CPkZoa9+1BkxstFkt+mTo1+JigjgRMh+PHli3h5OnK2M4kBDNnQvfu7bd17262WyyW4jF7tv/+oE4iRZgOZ9y4cHV1VWxnEpLMYbCdIWuxJJ8wnUSKKVP89z/2WG6ydHasN1cIpk07iC1b2nt0bdkC06Ztorr6X7kLabFYIhOk4urTJ1p9N94YPNJpaLCeXV5YA3wIvKIIWwO8JR1rgC8s+YjuPW6c/wikTx/vqBidFZvPJEaGDXOfpGgN8JZi0tVjc6l+GW933laqqxcQ9bJceik89ph3vevX27h8XtjOJAQzZ7pnXjz66OLJZLF05bS9QSquKVPKyM81EVasGGNVXS5YA3wIJkyA005rb3RXhdtus6EWLJZi4Bc3D3JLrx1kiP/Od7KvuzNjO5OQPPhgR/3shg024KPFUgzyaeoN6ojWr8/fuUuZku1MRGSEiMwRkbtEJOBdIneWL4+23WKxFIewc0v86OkWjs/iS6I6ExGZJyKrROTVjO3jRWSxiCwRkR8CqOoiVZ0MnAzk3YPGy9hujfAWS2EJUi1HmVvixc03++8PM/O+q5GozgS4FRifvkFEyoHfAEcBewKniMiezr6vAU8DeZ9ONHMm9O7dflvv3nYWvMVSaCZPzv85ggzsQfNRuiKBnYmIHC4iN4nIfs76pHwJo6oLgDUZmw8ElqjqUlXdAtwBHOccf5+qfhHIu2/FhAkwdy5UVZn1oUPNuvXqKF3+9re/cfbZZ/Piiy8CMHfu3OIKZAlFoeZ52CgX0QjjGjwV+DZwqYgMAvbLq0QdqQbeTlt/B/i8iIwBTgQqgAfdCjod3ySAoUOH5uyLX10Nxx+/I3Pn7sbKlcq0aZtZtGhplw1NXer89Kc/5Qc/+AGXX345n//853nyySf57Gc/W2yxLDkQdda7H5Mn+49A7Gz49oRRc32gqh+r6oXAEcDn8ixTJm7vB6qqTap6rqp+R1V/41ZQVecClwELu3XLfUpNY2MVt9228zaxVq7sydVX70FjY1XOdVsKT//+/enbty9TpkzhueeeY/HixcUWqeCOJZ2N3/42vrqCvLqsi3AGquq7AMdlrH8vqEwuC1ALvJq2/gXgkbT1S4BLotY7evRozZWaGlXjlNh+qanJuWpLEbj33nvbrV933XU51Qc8p+5teh6wKr1dO9vHA4uBJcAPM/aVAb9zqy9ziaNtlwpTprjfg6klbvzOlY/zJRGvdp25BL6uq+r8jPXrc+u+IvMssLuI7AysAL4JfCts4ThDTixf7h5mYflyG2KhFOnfv3+7NjFy5Mh8hSW5FbgB+H1qQ5pjyeEY1e2zInKfqr7uOJb80CljSSPOkUcY+vbterG4siWU7kdEwjrAfqyqa7MVRkT+CIwBhojIO8BPVPV3InIO8AhQDsxT1dfC1qkxhpzwjtEleQrdYMk3y0NOFBowYACVlZVZnUNVF4hIbcbmbY4lACKScix5XVXvA+4Tkb8Cf3CrM257YKnQ2uodN6tnz600NT0d6/nOO6+KmTNHeJxTufTSrpvON5NQUYNF5AlA8Y6qhrP/VlX9vc8xBSfOPNmNjVVcffUebN7clm+6oqKFCy9cbBtUiXLBBRcgIvjdByLCkUceyZFHHulbl1+ubKczeUBV93bWTwLGq+pZzvpE4PPAXbQ5lrysHvbAdLpS1GA/D6tsogTnes6uEEU41qjBqjo2d5GKQ5wjkzFjYMQIY3hbv96k8505s5wJE/bETIGxlBovvPBCsU7t6VgCNIWqoEtGDfYamSjV1U9GjhKc2zltFOF0wqq5LlbVK53v31DVP6ft+4Wq/ihfAuZK3DdcdTUcddSuPPDADtx661MAeWrAlkLwxz/+kVNOOQWApqamdurKm266ibPPPjtfp34H2CltfUfg3SgVxPmiVAr4z3zPn6p58GBYvbrw5y01wqq5FqrqqMzvbutJJU5VwC9+YQI8btxoY/iUOqNGjWLhwoUdvrutB+GnDnBRc3UD/gMchnEseRb4VhR7YJwq3FLg6KMPZuNGr/df5Ykn8jNCaGz0t5vk67xJwU99m07YyRfi8d1tvdPz5pvms3dvY5SfOdNOXipV0l+mMl+swrxohSEfjiWOfF1qZLJxo/e+Pn3yN0IYM8YvbJLNb5IibGeiHt/d1hNF3GquxsYqbr99OFCGqvHuOvPMFhYtskb4UmT9+vXb2kX6d7f1bFHVUzy2P4hH9AZLNArtMpzOeefZl0kIr+ZqAdZjRiG9gFTOQQF6qmr3vEkYE3GpuWpr3d2Da2rgrbdyrt5SYMrLy+nTpw+qysaNG+ntRPNUVTZt2kRzc3PougqdA76rqbnGjvU2vudb1eR97sKcv5iEVXOF6kw6A3F1JmVl7ol5RKC1NefqLSVMoTuTFF3FNdjPRTffj7EhQ/yM8Pk/fzGJ1TW4lIlbzVVVdRArV3a0uldVbaKp6V8512+xhKUruQab+HfeRvB8u+dOnuxvhLfuweHVXAWZAZ9P4np7a2iAs89ubwzs3duGoy9V4pwBb0cm+SMJI4NijoyKSdwjk9sIOQOetPhDnZFUh1FXZz7NxEXbkZQqp512WqgZ8KeffjqnnnpqASWzpOPXkQweXDg5vJg6NTjKcGen08+AzwcTJsA558DEiXDddcWWxpILTzzxRLFFyJqupObym/k+efIimpoK4UnpbYSfPbuVk09eUAAZkou1mWRJRcVBLF78MU1Nb8RWp8USha42z8Qd4ec/L0w4I/+Z8GVdfiZ8p+9M8nXDVVVBr16fYcyYz8RWp8ViSS7XXtum3rZ0JEymRYsL/fvDJ58UWwqLxVIorF3Un5xGJiKyPbBGVTfHJE/J0L8/vBspLJ+llHjvvfcYNGgQFRUVxRbFE2szgcK75SZFjuSRq5rrdmBXEfmLmhzxXYKGBnjiCdiwwcyIt95cnY+JEyfy5ptv8vWvf52rr7662OK40lVsJsWKFhwNG6Mrp85EVceJiNCFknk0NMCkSaYjARNaZdIk870rN6TORmNjI6rK66+/XmxRujzTpxdbgnB09RhdkToTEfkG8LCqrhORS4FRwM9UtWgZhoKIWxUwbdpBbNjQfgb8hg0wbdomqqvtDPhSpampiQMPPJDevXtz++2389///peJEyey++67dwEVUrJxi4WXotBzTPw8uvzmwnQJVDX0gkkjCnAw8BQmZ/UzUeoo1jJ69GiNAxFVM9+1/SISS/WWIjFy5EhVVX3qqaf04IMP1nvvvVcPPPDASHUAz2kJt+2kUlbmfs+Ban19YWWpr/eWxSTK7HyEbddR1VwtzudXgdmqOl9EZsTSq5UIw4a5vykNCxtwxpJIysvLAfjrX//KlClTOO6445gxY0ZxhQqgqxjgW1uLkarXnepqsEZ4dyJFDRaRBzBZ4cYBo4GNwL9Vdd/8iBcfccbmSreZgI3N1Rk45phjqK6uprGxkeeff55evXpx4IEH8tJLL4Wuw8bmyg9Ji4mVqzxTp8Ls2d77e/Vq/3wpNmHbddR5JidjMsONV9WPgUHARdHFK10mTDAdx5AhZn2HHWxH0hm48847OfLII3n44YcZMGAAa9as4aqrriq2WJZORlBHAiaIrF+HlVRCqblE5PsZm2ql/a99NDaJSoAJE6BvXzj+eLj/fhg1qtgSWbLlV7/6Vbv1tzIynB1xxBEFlKZr0KMHuOUcO+wwaGwsvDxx0tDg/2IZ1JGkI1Ja0YjDjkz6OcsBwBSg2lkm04XcgtNxEvIlajhqic66detYt24dzz33HLNnz2bFihWsWLGCOXPmWLfgPCDi3pEAPPYYjBtXWHni5rzzvPcNHBi9vtRzphQIGzX4MgAReRQYparrnPUZwJ/zJp0PInI8xhGgCviNqhZ0dNSnj/lcv76QZ7XEzU9+8hPAjEAWLlxIv379AJgxYwbf+MY3iilap8MYr/157LH8y5Er2bgHNzTAxx9HP9fGjaUT3j6qzWQYsCVtfQtQG5cwIjJPRFaJyKsZ28eLyGIRWSIiPwRQ1XtV9WzgdOD/4pIhLHZk0rlYvnw5PXr02Lbeo0ePDiovS/Y0NIQPP5Q+OvGf/V4crr02epnTT8/+fFFUY8Ukqmvw7cC/ReQeTDKsEzCJs+LiVuAG0hJsiUg58BvgcOAd4FkRuU9VUzqIS539BcWOTDoXEydO5MADD+SEE05ARLjnnns47bTTii2WL6XkGnz22QcT9nHz2GNtLrZTp34R6OF6XGXlFpqa/hGThOHJxj1461bvXChh8g4ef/wKzj9/SURJC0sk12AAERmNmbQIsEBjnv0uIrXAA6q6t7P+BWCGqh7prF/iHHqFs/xNVQPNdnG7T65YATvuCL/9relYzjuvbYg7eLB5e7EeXqXF888/z9NPPw3AoYceyv777x+pvHUN9iaqd1J9vbl//MqljikGUdyDx42LR31XLGN83Gl7t6GqzwPPZyVVdlQDb6etvwN8HvgeZr5LfxHZTVXnZBYUkUnAJIChQ4fG+vb26afdgIO56673efzxKlpa2jSGq1fD6ae3smjRG4wbV4gMcJa42HdfM2Xqk08+SfzbfqmQjarqO98J7ihK5WUtqCOZMgX+85/g44I8xYpNWNfgdZixGJjxWLvvqlqZB9m2nd5lm6rqdYBv0lxVnSsi7wHHduvWbXScQlVUmGAA//jH4HYdSYqtW8u4/vrdbGeScI4++mhSbu6q2uH7X//612KK1ymYPDl6ma6kPk4Z18vK/EcfZ5zRCToTVe2Xb0F8eAfYKW19RyB0JhHNY5jubt1g/frunvvXru3R5cNSJ50N1oMi73z6abElKB5Tp/rvnzKl7fvtt/tnctyyxXtfEogaNViACcDOqvozEdkJ2F5V/50X6QzPAruLyM6YUC7fBL4VtnA+jZQVFQdTViZs2VLueYyNJlwaqCqNjY289957nHrqqaxatYrVq1czYsSIYovWZUmiJ1dU5nRQvrcn3eV3woTgtMCJVnWFiQaZWoDZGM+pRc76QODZKHUE1P9H4D2gGTMiOdPZfjTwH+BNYHo2dccdWbW+3j+aqY0mXFpMnjxZp06dqsOHD1dV1TVr1ugBBxwQqQ5s1OAOBEXZ9VsGD05ulG4/udMjGQf9xkymTPE/vk+fwv3Gtt+Qn6jBn1fVUSLygtMRfSQi7n57WaCqp3hsfxB4MJs68zEyaWys4uqr96C11XtEkqKqahNNTXZkknQee+wx5s6dy6OPPrqtnaxZs8Ya4XPEb0Z4EH75QbTIYUb8Ji6GTZKVruJKceON/vNKkmxLitqZNDvzPhRARLYDWmOXKkY0DzaT00+HzSGy3ovAiSf2TEhaUYsfAwcO5JBDDqFfv36MGTOGDz74gMrKSvvf5UhQwii/hzJAeTm0tLhvLybXXuutkkr9niB7ides9p49YdMm73JJVXVF7UyuA+4BqkRkJnASZtJgYsnHyGT5cr8JSG2owrx5LQwatNh6dSWccePGceihh/L2229TV1fHggULOOOMM4oyMil2qKBCkZqP5WcncOtI/LYXijD2jd/+Nru6b77Zv+4wbtPFIJtJi8OBwzBP08dUdVE+BIubOCd21db6pxLNpKYGbGSO5PPGG2/w2GOPoaocdthhkY3vfpO7RGQecAywSp0Juc728cC1QDlws6pekbZvIHC1qp7pd94kT1oMM+nQ75iyMmh10X2Ul8PWrbnLlwtBExf99vfp4+/lFjTJMxc1n18YfLdcKvmctPgG8EbUcp2JmTM7JsjyY/ny/MpjiYfhw4czfPjwfFV/KyUSKqhQhHm7dutIoPgjk1wJGrUEqbqyZeBA/4CTGzeauIPZeMyHnbT4tKoenDF5EQozaTEn8qHmqq6GCy6o4qqr9mDLljIn74BQVtZKa2vHCYzWCJ9cvve973H99de3m7wI8U9aVNUFTqigdA4ElqjqUgARuQM4TkQWYUIFPaSqC2MRIAa8woLssIMJL1RIamoKe764CepIg1Rd2dhNgjqSFBs3Rqs3RdhJiwc7n8WcvJgV+TDAA4wZA48+Cs891zbkdOtIeveGa66xRvik8sorrwBFm7yYyFBBbkybNpKFCwfhZit8912loqKFRx55ut32xsYqYIRrmfYBEb2DJnptr6tbRFNTse2QQcEes88V7x9MEs46ayvV1U+77nNj1qzd+Pjjas/6osrnRqgQ9CJyu/OZg6Nf5+O11/x1lzU1NqVv0pk4cSIA12YTVzx3PEMFqepoVZ3s1pE4B80FLgMWdusWWVsdGa+OxGAm7s6atVu7rddfv5tPmTB4R9lNukOL6Ujzx6ZN0dzZ5s8P25FkT9hWOFpEaoAzROT3ZEilqmtilywm8jkDfsMG7zePo456n4svXgyAnaqQXJ566inuuOMOrrvuOnbZZRcyHVIqK/OqwU1sqKB0glxcDcL8+Tty7707btuydq330YMHy7bRepB7cEfKEj7SF+bM8UtAK6Hk79vXz0gfrg4Il5QsnV69wtfdjjAzG4FzgUXAZmAp8L+0ZWmYOoq95GOWcN++7rNUu3dXPeaY2E9nyQPXXnutDh8+XHv06KE777yz1tbWblt23nnnSHURMFMYk0ju1bT1bs79tDMmacdLwF5+dbgt+Z4B7zULPWhGd9hZ4lFnyZeX5/XnhiZMBAyvmf1hCLoucdSRufTq1bGOoHadWsLaTK4DrhOR2arqMm+za/KlL8Ejj7Tf1t2J+/jAA8aFeOZMq+ZKMueeey7nnnsuU6ZMYXYeU9qJyB+BMcAQEXkH+Imq/k5EzgEewbgGz1PV1yLUWZDkWKrh5lXRIYmT98i9uvrJbSP2IPtAJi0t2en046a11c/Wg+e+yZPD2XuCknCFSZj17W8fgmlaXigiyuOPL9i2JdumFHmeSamSD1/8iy+GX//afN+61QzX161rH92zd29rN+kqdMbkWA0NwZPzMkk9UqIkkIqSPGvwYPjww2gy5YOo881SRHnk+l0XEW/X6TDlwczjCXKzzts8k1Ijn29v779fS0tLDdttt5lRoz7ihRcGsmVLz3bHbNhgIwdb8kMhRiZ+aXPdCeel1XFkEX5k0txcnHS9mdTVVTFzpre3WraeXOn06nUwGze6P6ZV/euaNm0k4Oc4oTz22JPx2XTD6MKc0YsAO4U9PmlLPvTKP/uZ0TNWVqp+97v+uuWaGrO/pqa9vthSfFpbW3X58uU510MnjBqcjU1ANbq+P5tzJIFs7CZRyMVuEiTHYYeFkyFsuw7lGux0OgrcG1Mf1imoqDCfa9ea8AjDhrkfJ2KGw6rmc+LEsB4ylkIgIhx//PHFFqPT0NAA06fnr/5iB3lMJ0jNlCv5VI83NsZbX1Q1179E5HOq+my8YpQmqc4EjG3EK8yKasf1OXOMAd/aUpLBQQcdxLPPPsvnPve5YosSmsIY4P2MzO7qk6lTt7B2bXfP/ZWVHdVUIl9GtbQM8ABlZYe6TlaOS81liKIuNJg5P15zS5RRo9bQ1PRKRDkCCDN8SS3A68BWTJKql4FXgJej1FGsJR+qgDlz2oaMV11lttXXqw4ZYrb17+8/zKypiV0kS5aMGDFCy8vLdZdddtGRI0fq3nvvrSNHjoxUByWg5vJKvjRgQMdjc0ls5af+cVPzRqk7rGttIYh6XbKR3a++KVPcywS5c0c7f36SYx0VZ0dW6qSPTPr0MZ8TJhhvk6OOMsHaPvnEu7wNAJkcHnrooWKLkHf8osV+/HHHAH+5JLbyU/+4jcZrarLzjCo2XvlWvIg70MKcOe55UVQ7bss3kToTVV3mhMXeHUh3W0psM8inKuDNN6sAM9N1+fJFNDWtBGDJkkpgFCtXeqsCAPr1S4ZXisWwbt063nnnHbak+Xbvu+++RZTIn6hte/ZsL3WJYeNG5dJLF20LVbJ6td/xrc6+6PMs3FQz/p5R7VmzJjlqrpYW/2vanvbza8JSWflF1q5196hT9boW3nL17LmVpqbwcb1CE2b4klqAszCqrY+AJ4CNwONR6ijWkg8111/+0jZsvOuutu0vv2y2pdRdpTBc7+rcdNNNuvfee+uAAQN0zJgx2rNnTx07dmykOki4miuMGiY9x3iQeiVbFZgXYT2jknTflJfH89v9iOrRFfTfRPUmDduuQ3tzOZwHfA5Ypqpjgf2BD+Lq2EoNNzUXQCqc0/HHG9WBF2sSG9Gs63Httdfy7LPPUlNTwxNPPMELL7zAdtttV2yxYiOs92DYHONeKWdzId+eUfmgEHlVojrpzHENDZp9fWGJ2plsUtVNACJSoSZR1h7xi1UaBHUme+0Fp53mXd7LldhSeHr27EnPnkZzu3nzZoYPH87ixYuLLFV8BD1gkkBYl98kvYQlyU05hfrYS9KfU3ET1QD/jogMwMw3+ZuIfESEKKedjSfTVJUnnwxXX216/X5O1pennwYvu27KldiSDHbccUc+/vhjjj/+eA4//HAGDhzIDjvsUGyxfIliM8kuvpZXmeB8HY6EPmU7Etb+kKRkc1FtJtnberyvdbqdK+jY88/PYx6YMLowt8WR+GtAj2zrKOQSt82kvl61Z8/2usjevc32+vpg17zBg+1M+KTS1NSk8+fP182bN0cqR4JtJlH0+iLh9PSDB0erN8huEMb+kJItKUS5BiLZnyesDSmOSMMdzx2uXYdqrMDtzud5YY5P4hJ3Z1JT4/3H9u4drnGlOh9L8airq1NV1VmzZuVcV1I7k2zmi3i17/SHUrbzULyvX27li0HUDjVf50nhlRajEJ1JWJtJenKsgSIyKH2Je7QUBhHZRUR+JyJ3FeP8XnNEVq/uOAPeiw0b8ht2whLM888/z7Jly5g3bx4fffQRa9asabd0BrKZL+I352PwYPOZjSE3VdaNMPYHv/LFIEoTycW+EnZ+incyrWiRmbMhrM1kDvAwsAvwPO0VcupszxkRmQccA6xS1b3Tto8HrsUE5r9ZVa9Q1aXAmcXqTAYNipodzh07cbG4TJ48mfHjx7N06VJGjx6NeREziAhLly4tonTxkE07LSvz9q7KZeKdX9lCeEbFzbBh4Sdb5vL7JkyIngogk8mTcysfRNKSY90K3AD8PrVBRMqB3wCHY9KcPisi96nq63mUI2v8bkI3rEdXcSlUcqx8EN4A/2Wihkk3bTg4sVWU0PHQSnX1As9Je36T81IkacIiRJtsOXRoro4D3v9jmLD/J58cY7h5FxKXHEtEaoEHUiMTEfkCMENVj3TWLwFQ1cud9btU9SSPuiYBkwCGDh06+o477ohNzq985cu4B6ZTKipa2by5/Zj2sMPe5+mnt2u3vaKihfHj3+OJJ6qcwHhQWdnM9763JMM7w1IKjB07NpHJseJWb6Q/MqLW7fe4GTIkeBRVUwNvvRXtnPmmb99w83Pq63Ob4xGUbCwokVm2j/rQSd/CGFYKudAxT/ZJGNVWan0iZvQyGKN+exO4JKjeQhngU/lKUp5eAweazw8+aB8EcvvtzUzV7t071tGjhzXMlyIk1ACfjZE8rBE3Srmg3O1hcs17BTYsJmGM8Ll4cqXwq7++PryRPvp5Ywz0KCK3q+pEETlPVWMOVRZ8epdtqqqrgUAtYL5ic9XVVXH11Xt0GGnU1S2munoVo0fvzXvv9eTQQz/g1lt35qWXnqS6Wpk+fQAXXLAfF174IldeOZzm5p4d6t6yxWZnLBS/+MUv+NGPfsRdd93FSSe5DnAtMRFkMwhjh3zwwfjkiYswRnjNswJo+nT/a5dv4zuEnwFfTG+ud4Cd0tZ3JAETJceNW8WFFy5m6NBNiChDh27iwgsXb1NP9ejRSnNzGZ9+2o1evbZSXm5aU//+zQB88kl3Vq6s8Kx/1SrvfZb4+M9//sP777/PQw89xLp161i7dm27pdRpaCi2BG3EMVs8iQ4rg0I8AeP47X6ebEFOAPk2vgPh1FzAucAiYDOwFPhf2rI0TB1hFzqqubo559wZk4z6JWCvqPXmM7WpG6efrjpokAmcl67+uuGGcCqBJAWz68xce+21Onz4cO3Ro4fuvPPOWltbu23ZeeedI9VFAtVc2U4szIeaK0jVEkbNlcQcQGGvca74zeuJM39JJmHbdaK8uUTkj8AYYIiIvAP8RFV/JyLnAI9gXIPnqeprEeosQDa6jrz88j6sWTOQlJZu2TI47bRWZ7gZPCBsbrbh6QvBPvvsw+zZs/n1r3/NBRdc0GF/IdtMPvBTfQweHN1tOFNdEqWOoDkiYdxsjz463LkKSRg1VxwjEz/3YM2zGi0Mkb25RGRf4BBndYGqvhy7VHkgyOMlbiorYd267MuLGBfjVD7t5cvNzTZzpk31my9eeuklnnrqKQAOPfRQ9tlnn0jlQ3u9xETai9LZ9fX1rseMHftlvFxFp09f5OPWqp7bn3iizTW3sTG8a2xl5Rbmz/d+QQpT19Chm7jjjmTZEo87LtilOfO6ZYvf/2kI/s+inzOcl2KkzkREzsW42t7tbDoBmKuq12clZQEpdGeSq8GrpsY9p3zv3jB3ru1Q4ua6665j7ty5nHjiiQDcc889TJo0ie9973uh6yh0Z5LCr20HuZNGbafl5bB1a8dtYedWBT1uguRJvWQliTAuzYMHw4cf5n6ubJ8ruYxcwrbrqFGDzwI+r6rrnZP8EvgnkNjOpFhqrj59vsT69d2zKpvyCps6dTc2bGj/xrNhg/X0ygfXXnstN9xwA7169QLgC1/4Aueccw4jR44ssmTJws0jK+zDPYyqJyh9bxIn+XaSqDs5E7UzESC9ObUQfvprUVDV+4H7DzjggLPHjBlTsPMeeyxkzpHs3t28WaRlhXWlb99y1qzZEy9nopUre1LI39IV6N27N2PHjt2W02TTpk307t3bXucMamo6bgubBz3MMUcf7Z2nPrU/aYRxaY4j9FK2FCqeWdTO5BbgGRG5x1k/HvhdrBLFTLFGJtttVw3sTkqXOXToZs46y8R5uuqqPdiyJWWE79gXr14Ns2d76ayhrKyVpqYF8QvdhTnkkEPYa6+9OOQQYw58+umnGT9+fMkb4IOIaoR3y8ETNuZUmIda0DySJM4zCUMxk2jlEkstCtkY4EcBB2OedAtU9YV8CBY3hbaZ3HSTsXcAnHMOXJ+mCDz5ZLjnno665ygkwXujs7Fw4UKefvppVJVDDz2U/fffP1L5JNpMysrc20q6g0eUAIJudYWxGUA4u4GXvCmSaDMJkjlFHPdsNjaTXM+bL5sJqroQWJiVVF0IR/UOdMwDX1mZW0ciYh4C1ggfL6NGjWLUqFHFFiNWvB4kqe1RotHmqi4JY1sIcg9Oos0kjJorLlVTNu7chSJqDnhLSHqmRUlJ71jAdCa5eHup2jwolnB4qVfiVLuEfbiF6Qj8bCI9ethU14VSWWVD5JFJqVEsm8l//zsIMPMU3nvvTZqa3t62b82aWlRr6d7dhFzJhuXLkxWK25JMvOwZ2eTW8Oo0whjgRcJ1BH42kaSqdsOMuOLy+Iojr0m+6PSdSbG8udJvrr322pUxY3bdtr7QURJ+61tl3HabudEGDYLNm/0zpaUzbJhYT6NOiIjsAkwH+qtHaoUoeOXXSR+ZhM3B4zWaCdMxqYZTy/rF3mpuNiPypKl3w8zcT6J6Lm4ivRaLyDkiMjBfwnQmFqQ5W112WfuAe5WV5nPoUPO5apUxTP7eSQm2/fb+dYsk00WylLnhhhv46KOP8lK3iMwTkVUi8mrG9vEislhElojIDwFUdamqnhnHeRsavDuJ9A4grEHbq9MIYw8IazMIeugmMdBj0L0YdlRW6kTVsXwGk+nwTudGSPQck2LR0ABXXtm2vmaN8exKdSgvvWQ+U8c89JD5TN1wp5/e0WifjircfLNJyiNiliFDkhUhttR4//33+dznPsfJJ5/Mww8/TFQvxwBuBcanb0jLIHoUsCdwiojsGedJ/exq6fNFwj7oC+HeOnOmmY/lRRLf8IPclcOOyvJBIZ/QkdRcqnqpiPw/4Ajg28ANInIn8DtVfTMfAuZKMWwm06YdxKZN7fOUpGauL1q0lN/+djjp/fjZZ7ewePFiJ+z8rlx+uVJZ2YxId1ShrExpbW3f7zc3myXF6tUwcWIrixa9wbhxq5g1azfuv38HWluFsjLl2GPf5fzzl+TxV5c248aN47DDDuPZZ5/lyiuv5Mwzz2TMmDEcffTRVFdX51S3qi5wMoimcyCwRFWXAojIHcBxQKh01BlZRF3b9rJl3nGc6uoW0dRk0iU0N38RE5Dbn5YWdzvdmjVe50k/JpyNb9GiKlTb3x8pUpEhUnInheXL/X9/ZWW8QVtFvox7lteOqBbQthomtHDmAuwLzALeAGYDLwBXZlNXoZZChqD3Cgct4p2hsU8f1V69woWyDgpdP2WK+77DDivYJShZXnzxRT3vvPN0jz320MmTJ+t+++2nF110Uaiy+ITqJk8ZRNWnbZeXe7fD9nKHb1tuhAnBHjalgtf9AcnMsqjqL3OU3x6WKM+DOEL2+7Xr9CXSyMQJ9Hga8CFwM3CRqjaLSBnwX+Di7Lu1zoOXQW7YMG+db5gc0mFYvRp++1v3fY89ZueneHHddddx2223MWTIEM466yyuuuoqunfvTmtrK7vvvjtXpust4yGnDKIQPOpuaXF/Y1Zt/7ZaVnZoh5GvG15pEcKMbMKmVPB7y583r4VBg9oS0CWFuroqrrhiOC0t7tcw7KgsPIcSzkLRfgSad8L0OKkF+CVQk7nN+RwRpa5CL4UcmdTXq/bu3f4NoXdvsz3oLSbfSxKTCyWBiy++WN96660O21RVX3/99VB1EG1k8gXgkbT1Swg5Eslcoo5MMnOxh207XnnMwyS1CpsDPej+SGr79RudxS1zlPs9nvPlYWQCHK6qP8jYdhTwA1VdlE1n1hlJvfmn/MG33x6uuqrj9nwgYpqRF0EujF2Vv/3tb/zyl79st+2hhx7il7/8JSNGjMjHKZ8FdheRnYEVwDeBb0WpINuRSabtI+zIpF8/99FFVdVBrFzZ06VE+jGbaGoKjnQd9Jaf1PlV3nYjZb/9VtDUFKe90nv0lnnugl6rMD0OMAV4BVgPvJy2/A+oD1NHsZdCp+1NT7G5ww5mPYXXW0zQG16YN8Awb4jpsnR1brzxRt177721d+/eOnLkyG1LbW2tTpgwIVJdeLzBAX8E3gOagXeAM53tRwP/wdhHpruVDbMUamTipfuvrw9um2HtHfX1qj16eNeT1JGJ34gqbpmjpGKOA692nbmECvQoIv2BgcDlwA/Tdq1T1URH8w+TjS5uGhuruPrqPdi8uc2XsqKihQsvNPpet4xyFRUtjB//Hn/96w5s3Zranv72oYR7GwmmrKwVVaGqykQyTpoOupB8+umnfPrpp9x0002cffbZ27b37t2bytSEoJCEzUgXF0Ft2y8rX3rmvXCZAkFEefxx9zdd73MZwmZI/OY3vUc56fdQ0pg1azfmz6/G7Rr4XbdsCJ/dMq7sjiHbdZgepzMshRyZeL2lpL+hVFS0354aLXzxi6rjxqlOnarbRhJeb5hxLUn1kik1CPkGF/fi1ba93mAzRxhh33T93rCDbB1hbSZ+I5wkj6gLOTJRDX9vx3OucO061KRFEXna+VwnImudZV1qPadurxPi5bHltv3BB+Gtt9rsKX37mpAq++3XVibfIbdnz26b8NjQALW1JsRGbW3nnwh58MEHA9CvXz8qKyuprKykX79+29a7AmHiRgXN4g6aBR52sqHXcTU1yfZC9JuZ31WiVYQywKvqwc5nv/yK0znwcw0G84DevNl8P/PM9sb5vn3h3Xdh3Tqz3q9fuNg/uTJ9Ovz97zBnjnmnAXPOiRPN9htvzO/53WhogPPOawswOHiwiZoa50Pl6aefBmBd6oKXEEEGeC+jcKarahgDuqpSXf0kXvN+7777IMC9jm7dWqmreyOUi2pdXUcVMXkxYseL3zW8++5NnHxy3Gm2/dWKhgIn0QszfEktwDeAfs73S4G7gf2j1FGsJSmuwX77VFVPPVW1tlb1ssvMvuZm9zL5WPwmW06Z0l4dMnhweLVDfX30svX1qt27d5SlvDw/6o4777xT165dq6qqP/vZz/SEE07QhQsXRqqDElVzhTGgB0288yvfvXu0/8xt0m36PZJEvCYKR1HxRaGsLNw9HQdh23WkRgu87HweDDyFCf/wTJQ6irUUw5urpqZt1nvqRgiyp3z3u+bGvfBC1Z4929cXpNuurg7XwOJaevQIvsG9brKgsn6/Ne4ZxaqqI0eOVFXVp556Sg8++GC999579cADD4xUR6l2JkZ2/yXomsc5P8RL7qR6cqn6//58tNcw92dc1ytsu446zyQVN/SrwGxVnS8iM3IZGWWLiPQBbgS2AE2qmijt/oQJ7uqYIHtKymby6adGxZXOWg/rVGpuyYoVwfNM4mTLFqOG8lI7NTQYtZlXWa9w4g0N/gmX8pFprtyJYvjXv/6VKVOmcNxxxzFjxoz4TxQjcam5ACor/T26gmZx77ffbixb5u7NBOHnhzQ2VrF6tbunUlLnmID/zP2ws/+jETQLvsCz3yHyyOQB4LcYv/gBQAXwUpQ6AuqfB6wibaaws308sBhYAvzQ2TYRONb5/qegugs9MvEiaGTys5+Z9f/7P9Wddw4u57fk2wsstXiNMMLInG25uFUeX/3qV3XSpEm6yy676EcffaSbNm3SffbZJ1IdlPDIJGjUG/SWG9fIpNBeUXHhJ3c+1Fxh7sv4zhWjN1caJwOPAONV9WNgEHBRth2ZC7cSPlT3jkAqfWEWeeOKw8yZHcPL9+7d5inTt6/5vPde+N//2jyqouZx6N+fbYm38o1XqPMgmb1CmodxNjjvvOBjonDnnXdy5JFH8vDDDzNgwADWrFnDVVddFe9JEkyQR1eQR1LQfx3Wo8mvniTnBPGTLYlh8/NBpM5EVTeo6t2q+l9n/T1VfTQuYVR1AZDZrLeF6lbVLUAqVPc7mA4FSiiX/YQJMHeucXUUMZ9z57ape1510ielvL2WLTO5UAYNinaeTz4xn1oAlZfXAyBIZrdkSw0N4TrA1avjdVvu3bs3J554IrvvvjsA22+/PUcccUR8JygCXupAt44j6IEXlLMj1/JB9QwenGzXYPBut9Y12AURqQC+jglat62sqv40XrHaUU3bCARMJ/J54DpMPpWvAve7FQyT86EYVFfDrbe239bUZPTF8+Z11Bdv2ADdum2hvLybS8wi7XB8ijPPbKGyssVHF+5dNgpecZeCosm65XmYNu0gVP3dVNuO3UR1dTwul1u2bGHBggW8//77tKT1cqeddlos9ecDP5tJY2MV4G57cPu/crV5xBVTy801uEePFiZPTl4ek3T82m1xXIMLb1+KaoCfD3wCPA9sjl8cV7xCda/HJOjyRFXnish7wLHdunUbnRfpYiIVgsUr6c26dd350Y8Wcf31u7F2rUlFV1nZzNixq3j44e0z/PINmzeXU1HRgnenEYcOTDnooA89ZfZjw4ZuNDZWtQuPYRKEhSPKsUFceuml9OnTh89+9rP06BEcWiQJqOr9wP0HHHDA2WPGjGm37/TT3cuIwDXX9CTs8SmGDZMOZdJZscKoLb1S+waVTzFmDIwYAWefDRs3mm29e5czYsSejBkTayLKWFnl08+tWtXxeuefcNc7VsIYVlILGYbxfCzkKVR3UgzwXuRiwEwPKulm/PPzgY9j8ZItTJiOzLJRHA3iNMjutddeOddBggzwfvM+3Mg1UKPf/xZ1jojbHKOkzzMpZAj6oPMVywAfdWTyDxEZqaqvxNSXhSGnUN3FSNubDX6uhUHpSquroXv3Q2lu7qhiqKoyQ+xBg6q4+eZdWLWqgqqqzWzcWBYquF842TsOqRsbq/j4Y/f0q35l3dUtKcNP+21xzoqura1l3rx57LLLLrHUV2y8oiak534Pc3yKIJuHn+H8tNOi2TvOO699Smowql4vV/Kkkw+bybXX5jeVRVaE6XFSCyY/dTPGTfdlTFj6l6PUEVB/3kJ1l+rIJOyM73337fh26fc2F+eserc3r7AjjKSMTEaMGKHdunXTz372szpy5Ejde++9t01kDAsJGpl4jUa9RhhBo9cg99a4RiZBo+yk4jeyy5dLs9//FedEybDtOurI5Kicey8fVPUUj+0PAiH9QdpTKiMTN8NjeXkrP/zhG1RXr/KMiZRiu+1GUFk5kObmMjZsKGfoUBNe3qtsdTVccEEVv/jFCFTBK1R5kF3Fa9TkN9JKrz9zdBGunGHZsuhGxlmzduP++3egtVUoK1OOPfZdzj9/CT/+8Y9dj09ym/Fr216xsryMwX6xtSA4uZV7TC3Dhg3hnSWmTfOWI2yCrWLQr5/3pM9s2mk4vO4VZfLkAk9YhMgjEwHqgB8768OAA6PUUawl6SMTVfNWtuOObW8XdXXhy55zjuqAAaqTJqkOHRq+3N57q44e7R4v7LDDOr5xde/epq/t1y+3CYtub61RRiZRE315vX0fdphqa2ur3n777XrZZZepquqyZcv0mWee2VbWKzxOOiRoZOIXZ82NOJJbxRGfqlRD0Pfp4y13ZjKyuPCalBz3CC5su446MrkRaAW+AvwUWAf8BfhcPF1b/JTKyATMaOGWW4TDD/8yANtv/wZNTe8HlmtsrOKWWz7L+vXl3HZbK336bKWp6Z+hyi1evAfNzWX069eMcePVbaOaceNWcdBB7W0tZ521FIArrhjOunXCtGmbWbSoY4KturoqrrpqD7ZsaR/91c3tOf2tta6uiiuv3IPm5vYjtJYW6VBWNZp78Jw57m9yjz2mHHDAKey00yZeeOEFDj30UNatW8dFF13EnDlzaGys4pe/HM7Wrcb+s2wZ1NUpf/7zCs4/P5mRbIMiV4c9PkWQzaShwUyS9ZMnDF5yJH2eyfr13vu8PNxy5bbb3O0mt9+en/MFEqbHSS3AQufzhbRtsYVTyedSCiOTFKnEWX/6U/CxbraPsrJwUXndbCZ9+/qXDYp6nH5cZWXbMX7eJ5lvUuee27Y9NQrwe2sOQ1Ad5eX7q6rqfvvtt61MKpyKn+zpv5sEjUxK2WaS2b569Ur2qETV/9rlMwxMmBFzroRt15EaLfAMUJ7WqWyX3rEkeSmlziT18PrrX4OPDZPVMUq5oJs/zPm8Opww0WDr61UHDjTbd9wxONpyWFVXsPrsQP3977fq/vubTmXVqlXbOha/cumGziR1JlHbRa6xteJUT9XXqw4Z0lZ+7txo5YtB2BeOUiRsu46q5roOuAcYKiIzgZMweU0SSympuVJ062aMkEuWvEBT0ye+x3oZrJctUy69dJFnvmw/Q7efwdSrXLqL77RpB7FhQ3sjamoWf0VFeTsjbboBPzVxM7X/nXfMLP5FixZTV4dr3mvVcKquYMP+uZx11nFUVr5NXV0dCxYs4IwzzuDSS1/HayY5wOrV/te5WETJ9um3HaBHj+C4WHGqpyZMgF694OtfN+unuLrlJItrr4UzzjDRsNOZMiXZ6rlYCdPjpC/AcOC7zjIiavliLaU0MtlrL/NGEyY3U9wjjCC1Rpg33qAkW6n95eXt1S5BdfvJG0S4POeL9IYbbtAbbrhBX3/99VDXqb18nXNkEia5lVcyszA5b9y45562OjZujF6+GBRC5VQMwrbrUCMTEfm+x66jROQoVf1VHB2bxZCKHJz69GPmTBMIcsOGjvv8Jnr5lQNvg6lbud69zcSs2lrzhltW5m50HDTIGA1VzXpLi1n/0peMjEFvzV7hOjID7E2daoJntrSYMmPGeOeCgbamO3BgW4DNhx56iIceeohly7yafkf5Co3fqNtr8qfXRE+/2FrNzcGjv+pq6NXrizQ3t3eP3bIluxhqr78+GBgJwN///iTl5RqpfDHwirnXVQir5kqladoD47l1n7N+LFDAJMPRKUU1V3PzvsBAXnnlH6xYscX32NR8ETcVEHgH2EuVmzVrd9avbx9Dy2/GfarcDTfsxief9GDQoM0ccsgHzJvXFh/MPPC1nTwVFS00N7ewYUP7h026Ss0rj3ZqfkFLi7uqSrVN1TRr1m7Mn9/2EG1pMd5a3iouk/td5A3gH/zrX58H4J///CeDB38euMCnbHv5Co1Gjs0lvPjijowZs2OHPUGxtcLEl/LqsLOJTfXAA23fzzjjy/ziF11IXVSqhBm+pBbgUZwc8M56P+DhKHUUaykVNVd9vfFeAdWddgo/VN5pp2hqjRTvvdfx+DDnfPxxc/zjj4dLrOSXZzylUgvyFAuTGjVsbuyOMhyuN920dtvvW7t2rfbseWSosilVHQlSc0WdZ5KrAd7v/43qzVRf3+bR6NYOLIUlbLuO1GiBN4CKtPUK4I0odRRrKYXOJKzbrVfZzJs4TNnf/759mcGDw53vpZfM8Sk3Xu+HdHBnkOnNlXooZXZsQe69uQW03EOHDdu07VybNm1S2CNU2c5gM/HzxgrTjnL1tstFdkt+yVdnMh14CZgB/AR4kSwi+BZjKYXOJNebKP0GDjPCyMVoev313g8fr4dRmM5y/Xqz/fLL3c/rNwrKLU3xzxX20Z/85Cc6Y8YM3XfffXXAgF+EKpt6209SZxJ1nonfyCTM7PeoUYr98LvWlsKTl87E1Mso4Dxn2T9q+WItpdCZRFVNpJP+1h5WPZZtvu1sgkSm6quvVx02zGzr37+jnMuXm31ecwuCRie5LJ/5zPM6a9YsnTVrli5cuDD07+wMIxO/35rLyCSb0YTXS0G+wpJY/AnbrqPOM0FVFwILIxtnikQpGeCDDNBepOZnmPmk8PbbbfMz/OY/+M298MuM5zaPJIhUfdXVMHGiCTD5ySe0C8fS2FjF7Nm7AhVcdNEWli1b0kH+6moQ+TLuScTU8/cEUVHRwpln9mTfffcF4JNPPmHRotcpL9+NVJiZnj1b2Lq1bFtYldQ54wyFHxdR55NMmAB//zvMnt1xX5jw717egdmEX/dyAshXWBJLTITpcTrDUgojk2xtJvmYBe9XNigoYNDIxO03TpkS/rfHNRJJ/x2ZtiIvOd2CX6bkpJOOTFLXKgg31Vo2hnNrM0kWYdt1wRt+sZZS6ExUs5v4lK16LFubSZTIvqklpXf3y9sS9gGSzfndztejh/eDLxs5k9SZ1Nd3/H1+/2uu3lx+dVhvrtImL50JcA4wMEqZpCyl0plkQy43cX19e6N2GG+uXGwmUUc1bh2iVycYZenoQny9wpqc5IyjMwH6ALcBNwETwpTx6kwyr5HfTPZcvblUo/2HQVx2Wfu2YzuS4pGvzuTnwBLgTmA8IFHKF3PpzJ1JLi7FuZyzZ89oD1vVeEYmqmHDo0RZpivsqvANfeihh3TYsNbIcnrddMA8YBXwasb28ZispUuAHzrbJgLHOt//5FZf5pLPcCphIlCn2kNc80xUVZ97rq28pbiE7UwiGeBV9VIR+X/AEcC3gRtE5E7gd6r6ZlR7TSEoJQN8tqRmpWfmHQmToTGXc37+83vy9NNDXENwZJJyInDLyFdR0cL48e/x8MPbewaBzGTNmi+TrbEdoKysldbWdLl/DvyMAQMe4Morr2TdukmUl3+LlpazgV1DyekTDPFW4Abg96kNIlIO/AY4HJOi+lkRuQ/YEZMOGyBrk3NUA/zMmfDtb3fMvd4t5BNi+nTz6M9EJDhIpBtlwU3KkjCy8eZSEXkfeB/YCgwE7hKRv6nqxXELmCvqE3KiMzFmDPz856m1nsCezpI/dtsNXnrJePBkRktNp0cPuOYaE1JjzBgYMQIuuAA++AA+8xm4+upyJkzYkYYGmDwZPv0Uampg5sxyJkxw/x1ByZwGD4bVq733f+c7Zdx2W2aMMeGii4axatU+vPvuu4wa9QmPP34SqodTU3OlI4+R83vfg48+gh13hCuuMHJ6PTRVdYGI1GZsPhBYoqpLAUTkDuA4TMeyI2YOl+cjVUQmAZMAhg4d2uFFKapnYK6xtbw8A1WV6uonI7/UvPlmH1I59zrrS2CnI8zwJbUA5wLPA48A3wC6O9vLgDej1FXopTOruYrF979vkmmdeKJuU2UNHtw+hamXDeanP20rk64T/9a3VHfdNfjcbqq97t3bbCE1NUY2N7VLKvRKur2osvJara0dpUcccYTeeeedumXLFsfu0KKwSwe9/TnndJQfH3UAUEuamguTvuHmtPWJmNFLH+AWYDYx20yCHCv81IBBxO2B9corVs2VFPzadfoSdWQyBDhRVdu9E6pqq4gcE7kns5Q0lZVmFFFVZSLurlkTrlxDA1x+ufmuakYYkyaZ9U8+gf79g+tIzXmYOrUtwGC6imbZMuje3YyK0kdNvXub3BOpOoYMgfHj4cQTP2TGjLupqanZJuOkSdDcXAY80E5GgJtucpN/yKBwVwBw19Gpqq7HqJCDK/BR4S5aVEVr63DSBzctLa0sWvSGq9oQoKzs0AzVX2p7K01N/vFcvdSXXmrKIJYt640ZvNmRSckQpsfpDIsdmcTPr39t3hxThviwXjd+b7EHH6w6dmy486cHxfRaBg/2d7VOGXrnzw8vo7cb7T6b1aP90XFk8gXgkbT1S8gyNFEcBnhV/+sYhswMiTffHK6cG2+8Ee3clvxBPkYmHnlNPgGeV9UXs+/SLKXI/Pnmc9Mm85n+9u43W9rPOFxZCbvuGu7806fDxo3+x6xZAx9+6L3/6afN53HH/YqBA+Goo2D06JQ9pj8wGtgvUHZD9x5+ezN4FthdRHYGVgDfBL4VobzvyGTZMncbhl9kg6FD3e0sQ4eGC7FfXQ3nnDOEGTP2dtYX0NTUGljOjRUregEmHYAdmZQIYXqc1AL8AfgPcI2zvAHcjrkxLo5SV6EXOzKJF78YWdnOvK+pMXG7TjstnAxh5oIExRhrG9mcorC7duv2fT3qqO9rt257KNQpHKDwy5xGJsAfgfeAZoyB/Uxn+9HO/fQmMN2tbJgls21n+99EnZvixvnnt5WNkkIhk9SoN8qI15IfyNM8k0eAvmnrfYGHgV7A61HqynUBdgF+B9wV5njbmcSL34zpoElqXhFt00OVhHmABM3aDppr0778EQrrtp375pvXaVnZkQobFEa0q89rXg8MWaoFvAdSS2bb9puD43c96utVu3Vrf3y3buEf5HHNXHdTX9oZ8MUjbGcS1QA/DEh3Am0GalR1o4hsDluJiMwDjgFWqereadvHA9diIhberKpXeNWhxqXyTBG5K+JvsMSAn7rHK+VvigcfdN/+2GNt38OozPxSDxvX4ijqtuWYgI5me11ddy6+eBlr1vQCKlzrO+MMY9xP7aur+zCkC0I8eKm5Vq92V3GBv5vu1KlfZOvW9pq6rVth6tQtVFf/I1CeadMOYvPm9mqy9EyaYZk27SA2bsy9HkuBCdPjpBbg/2EiBv/EWZ4DfoxxZ2yIUM+hmFD26QbJcsxwfxfMXf0SZoLBSOCBjKUqrZwdmRSBXJIhRQlVEibDX8oFOMqopuNv+KnC/goztH//GTp69Gj9+tcvU/hUjznmW67lv/Y11X32aVsnIbG5sjWk52qAzyWFQpxyWOIlbLsOPTIREcHM5H0QOBjz6jNZVZ9zDgmdoVkjTOJS1csxoxhLgnAbFYiYSYdBubqDJhym42/wNueaNQuee848biC8I0Dbb1DgdOBound/mmOPVc47bw7/+McB/OUvcMstDa7lhw6FZ54J9zvygbcB3ntk4mV8z62cIdsUCpnk4qJsKR6hOxNVVRG5V1VHYyYuxk018Hba+juk3DlcEJHBwExgfxG5xOl0Mo/xnSVsyZ5UCJeZM83s9IEDtzB1qsk/EnSZ6+qqmDlzBGHCoYTJ5fLccx3rCqMWSf2Gq67agy1bjmfo0L9z1lm9GDduFffe25tZs7YC5ey112bOPntpu9wqjY1V/OEPn2X9+nKGDjX7C41Gju4g+B3nFTVg8GD/cilOPNE9H8qJJ/YMVT5Fq4cDWGtrWaR6LIUlagScf4nI5/IiidcrkQequlpVJ6vqrm4diXPMXOAyYGG3sEGGLFlRVub5V3Vg3LhVVFY2Bx5XUdHCWWf5P6RvvnkXvDqlVasqQsnyhS+spm/f0cyYceu2JF1XX70H69d3A4RVq3py9dV70NhYBeC5P+KkxbxRXh5te4qTT462PRMvW5jXdi+ceaOht1uSgaiGfwiIyOvAHsBbwHrMXayquk/kExs11wPqGOBF5AvADFU90lm/BFO5a0cRlQMOOECfe+654AMtoUjNEG8f2wrmzg1Wc6XKT5zYpprKJIwBHUxAQL863norWJazz4Zbb92T1tbF1NbW8vbbfWhuVkzzfrlDfbW1Xmq6fbeovhTcg8VEmprr7Pr6+m3bZ83ajfnzq2nfySrHHbeC88/3zgh53HFfZO3ajlNlKiu3MH9+sAH+K1/5Mm4ZMEWUxx8PVpOlSHXWmbPpL7zQP3OoJT+MHTv2eVU9IOi4qJ2J67uBZoRXCVlXLe07k24Yn/vDMJO4ngW+paqvRa074zyuN5wlN775Te8JbnfcEU4/Pnast47+iSfCPXy85ABl+vRFoR4+N964K/fd18ytt/4bgFNOOShNrrYmn3ooej004QBUn8s+lHGWuL0oTZ1qOvaWFjMimTQJbrzRvx7xkTzMY8Krkw3bqaczdy585ztt5cO8WFjyg4iE6kyienMJUAf82FkfBhwYpQ6nXF4ncbkt1psrXuLw3PHyCBswIHwdfnNWwmISMbXqLbfcrpdddpkj1zKFZ1w9y7IJp5LPJa62nasXVZwZEj/6yHpwJQXyNM/kRqAV+ArwU2Ad8BdSsaJDoqqneGx/EOMtFhtdIZ9JMYjDc6euroorrhjeIR/Kp5+2cumlb4QaVdx990GYkPvtefXV8HKsWlUN/Io771zJa689z8SJJ3HVVVVs2fJdzAC5fdBCr6CGmze/uyLUCROKiPsIxG/Eks6ECbBkCcyYYdZzGVGEPaclOURVcy1U1VEi8oKq7u9se0lV982bhDFhbSbxkqvNJMWQIe4eRGFVI142ExFvr6BMvvMdmDt3FCIL6dZtf2655QUA6ur2BV5yfSg2NMCPftQWT+zGG6GuLqQ6ICbiVuHGoXZ8441+TJkyGoAnnmjKWpYNG8r56lcPybkeS+6EtZlEHZk0OxnijM5LZDvMSCWx2JFJfogru6NXxkS/gITp5DpCamys4pZbhgPdUW2huVk488wWJk36J1DGGWf8j4kTjSEg/XdVV8PEiVVcfvkI1q6FadM2U2hvLo058VtNjZfNI5xrMJiONUUuMn36aTz1WApIGF1YasFMTLwPYyCficlf/Y0odRRrsTaTZOJlf0glsArCK05WWD192/nrFY5V2EHhR9qt22cV7tQbbwx/XhjVoiXctnO9lqqqL74Yj61j/XprM0kKhLSZRJpnoqoNwMXAL4B3geNV9c9xdm6WrsXMmSaJVSbr1hlVUhATJhjVWk2NUW3V1ERTtbXNsJ8AXAn8CNiBrVvvBb7Bj3/sLsf06W4xwaSkM5fnei0hvtzt1mZSekTNZ1KBianV3yn7DRFBVX+aD+HiwKq5kk11NXTrdjDNze2bYtjc46k6br21/bawf3WbmmwzJuzcJ8BWwLwjffjhjznzzBYWLWo/x8Er53mpM2FCbi64Dz3U9r221hrguxJRbSbzcZJhYe6+xKMx65Ut8dLQ4J3gatWqaGE4suGaa8ykxY0bj6MtGVb7eYebN5dTX78nP//5ntu2RYkvli+S9qLU2GhC05iYreb6uHXEYWhuFkysMJscq1SI6s31qqaFjC8lrDdXMvGeTZ7dZLdsuOkmmDRpb+BVz2MyvcPcvNmSNGmxGMQ5abG5GXo4k/EjPKIseSDspMWoI5N/iMhIVX0lS7kKTtLe3izt8VYXKXV1i2hqyn/4jJ13FuCLHH/8Q/zzn2NDeYdVV8Phh7uFLem6+KVjjopVc5Ue2cTm2h1YilFzZR2bq9Ak5e3N0h6vt9nBg/1zt8dNWdmeiPyXqqpdWLmywvFeNLG5vObPdJwjY0cmcY1MWlvbAlPakUlxCTsyiep7cRSwG3AEcCwmz8ix0cWzWAwzZ5rJjun07g3XXltYOfr3f4gJE5bwr389yq9/fT8VFQ8A9/t6NLlNtuzKuHnm9ehhtkfFjkxKj6idyXLgEOA0NcEdFRgau1SWLkMc7qhxUFk5jP/97yluu+02zjuvhuHDhQMPXMlbb9kAg7mQ7ajCdialR1FicxUSazNJPrm49sbFRx/V09zcyvLlj7Fly0m8/PJQVL/LZz7zFGedtdTVG6my0j1ke1dl+nRjOE+nudlstx1y58fG5rJYgD59RnHwwQtZvHh/PvjgBcdLy8Tm8rKZNDRAXV36lq5tM4kjTlpmObA2k2KTL2+ukovNZbGEoVu37qxb18I77wgtLQAfkNICb9jg/nY9YUJmZ1JYkjbqjisHfBtjADvPpFSI2plcB9wDVInITOAk4NLYpbJYCoiZOHku//znCcAqYDpwF/Dzbcd4ubd6BUcsBEmbkNs2AbRtW+/ecM01uU0+TcJvswSTbWyuyzHJrWxsLktJk5p82Nycis11CbA9cC/wjW3HDRvmXt7NG62rMmECzJrVtp6LM0V6PLTa2nBx2izFJerIBFV9A3gjD7JYLAWnfcDG4c7SES/31tSDcto0WLkybulKj5NPbku3m230glQHn2LZsrZ1a8hPLpE7k1IjaXplS7IIE7CxsnIL1dX/8PQwq66GE07YkTlzmrfELmCJEUfUYLeIzF52K0ty6PSdSdL0ypZkESZg44039vDV2zc0wG23AXTv8n7CcXQmXv9HsQNrWvwp6fwLFkuuBNk8Bg8OfhuePt078nFXI47JhqkwKmG3W5JBpx+ZWCx+pDqKM8+EzS5JFU4+ObiObAIZdlbiGJkY1+zw2y3JwI5MLF2eCRNgwAD3fQ8+GFzey9MrW0RkFxH5nYjcFW/N+SeOkUlNTbTtlmRgOxOLBVjlEek+zKgjXVUmIvNEZJWItEuOIiLjRWSxiCwRkR/61aeqS1X1zHCSJ4s4RiZewT+zCRhpKRy2M7FYgIED3beHGXWkglVC8xbgVmB8+n4nasRvMFG39wROEZE9RWSkiDyQsVTl8DOKThydSVKCf1qiUbI2ExE5HvgqUAX8RlUfLa5EllKloQHWreu4PUr4dBNa5eVXVHWBiNRm7D4QWKKqSwFE5A7gOFW9HJPGodMQV7TfXHPRWwpPUToTEZmHuYlWpacBFpHxwLWYJNI3q+oVXnWo6r3AvSIyELgasJ2JJSvcot0C9OsX2wOtGng7bf0d4PNeB4vIYGAmsL+IXOJ0OpnHTAImAQwdOjQxc6j+9rcqzOALPvOZTZ4Rly2dj2KNTG4FbgB+n9qQpgo4HHOzPSsi92E6lsyb6QxVTbXQS51yFktWeNlF1qyJ7RTueYk9UNXVwGS/ClV1LjAXTNTgJMyhamhoH05l5cqe/PrXezJixJ52lNEFKEpnEocqQEQEuAJ4SFUXup0nqW9vlmQRX7TbIYM8drwD7JS2viPwboSKXUladIdp0w5iw4b213HDBpg2bRPV1dlEDbaUEkmymURSBQDfA8YB/UVkN1Wdk3lAEt/eLMnjmmtM7Kf0EB5Ro92aQITDvJxXnwV2F5GdgRXAN4Fv5SAykLzoDl4ecatW5RY12FIaJKkziaoKuA4TEt+/0oS9vVmSR3U1XHBBFTffvAurVlVQVbWZs85aSnX1qtAZH6dNOwiQMhH5IyYRxxAReQf4iar+TkTOAR7BqG3nqeprucqdtLYdfz4TSykRKdNirCc2aq4HUgZ4EfkCMENVj3TWLwFwMz5mQ1Ky0Vk6JybLYNfOtJiK9ps5wrNuvaVNvjIt5pO8qAKS9vZm6ZyYt/LCnjNpbbu6Gg4/fDfmz98RgLKyVg4//F2qq5eEHuFZSpeijEzSVQHAStpUAUcDs2hTBcQ25zUpb2+WzsnUqTB7th2Z2JFJ5yPRIxNVPcVj+4NAiGhI4Una25ulc3L33QcVW4SiY/OQdG2SpObKC0nzeLF0Trw8mfJJ0l6UvBKNLV+uNDU9WXiBLAWl03cmSbvhLJ2TYthMkvai5JVobNgwsa7BXYBO35kk7YazdE6uuQbq6rS12HIUk5kz3W0mNtpv16DTdyZ2ZGIpBNXVAMsLmlg2aW07jvk6ltKlaPNMCk1SPF4snZewXi9xY9u2JZ+Ebdc2n4nFYrFYcsZ2JhaLxWLJGWszsVhKFNu2LUnC2kwslpiwNhNLZ8TaTCwWi8VSMLrMyEREPgCycd0cAnwYszjZkhRZkiIHJEuWPVS1X6FPGtC2k3R9kiJLUuSA5MjiJ0eNqm4XVEGnt5mkCHMx3BCR54qhunAjKbIkRQ5InizFOK9f207a9UmCLEmRA5IjSxxyWDWXxWKxWHLGdiYWi8ViyRnbmQQzt9gCpJEUWZIiB1hZgkiSTEmRJSlyQHJkyVmOLmOAt1gsFkv+sCMTi8ViseSM7Ux8EJHxIrJYRJaIyA+LLMtbIvKKiLxYSK8hEZknIqtE5NW0bYNE5G8i8l/nc2ARZZkhIiuc6/Kik/o533LsJCJPiMgiEXlNRM5zthf8ukS9JiJyidOeF4vIkXmW409pMrwlIi8622tFZGPavjkxyhH5v8njNfGS5SoReUNEXhaRe0RkgLO9GNclvraiqnZxWTB56N8EdgF6AC8BexZRnreAIUU476HAKODVtG1XAj90vv8Q+GURZZkBXFjga7I9MMr53g/4D7BnMa5LlGviyPgSUAHs7LTv8nzJkbH/GuDHzvdar+MK/d/k+Zp4yXIE0M3Z/ss0WYpxXWJrK3Zk4s2BwBJVXaqqW4A7gOOKLFPBUdUFwJqMzccBtznfbwOOL6IsBUdV31PVhc73dcAioJoiXJeI1+Q44A5V3ayq/wOWYNp5XuUQEQFOBv4Yx7kC5Ij63+TzmrjKoqqPqupW57B/ATvGcb5sZPEpEvm62M7Em2rg7bT1d/C/+PlGgUdF5HkRmVREOQCGqup7YBopUFVkec5xVAbzCqVySyEitcD+wDMk67q4XZNitelDgJWq+t+0bTuLyAsi8qSIHJKPk4b8bwpyTTJkSecM4KG09UJfF4iprdjOxBtx2VZM17cvqeoo4CjguyJyaBFlSRKzgV2B/YD3MOqUgiAifYG/AOer6tpCnTcEXtekWG36FNqPSt4Dhqnq/sD3gT+ISGWcJ4zw3+T9mnjJIiLTga1Ag7OpGNcltrZiOxNv3gF2SlvfEXi3SLKgqu86n6uAe4hpKJ4lK0VkewDnc1WxBFHVlaraoqqtwE0U6LqISHfMTdmgqnc7mxNxXXyuScHbtIh0A04E/pQm32ZVXe18fx6jj/9sjOeM8t/k9Zp4yIKInAYcA0xQx0hRjOsSZ1uxnYk3zwK7i8jOItID+CZwXzEEEZE+ItIv9R1jwHvVv1ReuQ84zfl+GjC/WIKkHhAOJ1CA6+LYAH4HLFLVX6XtSsR18bkm9wHfFJEKEdkZ2B34d57FGQe8oarvpMm3nYiUO993ceRYGsfJsvhv8nZNvGQRkfHAD4CvqeqGtO0Fvy6xtpV8eA50lgU4GuP18CYwvYhy7ILxrHgJeK2QsmDUE+8BzZi3lTOBwcBjwH+dz0FFlOV24BXgZecG2L4AchyMGfK/DLzoLEcX47pEvSbAdKc9LwaOyqcczvZbgckZx37daccvAQuBY4v53+TxmnjJsgRjj0htm1PE6xJbW7Ez4C0Wi8WSM1bNZbFYLJacsZ2JxWKxWHLGdiYWi8ViyRnbmVgsFoslZ2xnYrFYLJacsZ2JxWKxWHLGdiYWi8ViyRnbmXQxnPwFFzrf/5FlHQNEZGoW5Xo5AezKszlvRl09RGSBE67DYilq23bKdun2bTuTIhJHowuoX0TE8z9W1S9mWfUAIJsb7gzgblVtyfK821CTFuAx4P9yrctSeiSwbUMXb9+2MykwIvJnEfmViDwBXJKx71QnFPRLInK7s+37IvKqs5yfdqzX9lox2dRuxIRk2ElEpovJltYI7JF27KcZZW4Sk4XtURHp5ey7V0zY+9ekLfT9FcCuYjKzXeUcVyci/3a2/dajo5yAExPJOWd6Vr4LRWRG2r43RORm5/c1iMg4Efm7mEx5qWB09zp1WoqMiDwubdn6NonINzL2h27bXvsS3rYhRPuO0Lah1Np3vmMH2aVDjJw3gJ+6bN8LEwNniLM+CBiNiZvTB+iLiduzv9d2p1wt0Aoc5Kynju0NVGLiAl3o7Ps0rcxWYD9n/U6gLiWH89kLEwRuMBkZ4YARwP1Ad2f9RuDUjN/XA3g/bT2zjguBGRnyjMS88DwPzMOExT4OuNc5rhz4oNj/qV3a/c9TnPZTnrYtdNvOaLOZ7T6RbTtK+w7btkuxfZeMPq4zICI9MTfST112fwW4S1U/BFDVNSIyEbhHVdc75e/GJBoSj+0vOHUtU9V/Od8PcY7d4BzrFfn4f6r6ovP9eUyjBzhXRE5wvu+EiR76fkbZwzA39rMiAubmzAy/PgT42OPcXvK84sj8GvCYqqqIvJKSTVVbRGSLiPRTkz3OUkRE5FRMvp2va3tVT5S2/QImKKHbvvtIZtuGaO07sG1D6bVv25kUlr2AZ7QtZWc6QsfkM24Javy2p1ifsR4mmufmtO8tQC8RGYMJIf4FVd0gIk1ATw95blPVS1z2pdjoUjb9d3T3kac1bb2V9u22Atjkc15LAXDUWhOA41S1OXM34dt20L4ktm2I1r7Dtm0oofZtbSaFZSQm1LMbjwEni8hgABEZBCwAjheR3mLymJwAPOWz3Y0FwAliPE36AcdGkLc/8JFzsw0HDnK2rwP6Zch+kohUpWQXkZr0ilT1I6DcGZ2lqBGTw6EMOBQzrA+Nc60+cHl4WQqIiByDMVqfqKpuD74obZuAfekkom2Dbd9gO5NC49mZqOprwEzgSRF5CfiVqi7E5IP4NyZf882q+oLXdo96F2Ky3L2IybLm1em48TDQTUReBn4G/MupczXwd8eAeJWqvg5cislR/zLwN2B7l/oexagwUqwGfo9RPbwKnCoiu0aQbyzwYITjLfnhNkwmvr87Ruoz03dGadvO8aHad8LaNoRo30TrUEqqfdt8JpaCISL7A99X1YkiUgs8oKp751Df3cAlqro4Lhktlmzp6u3bjkwsBcN5u3zCx7UyNGJSKd9bKjeapfPT1du3HZlYLBaLJWfsyMRisVgsOWM7E4vFYrHkjO1MLBaLxZIztjOxWCwWS87YzsRisVgsOWM7E4vFYrHkjO1MLBaLxZIztjOxWCwWS878f2IPCNW4Vt6IAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"intensity_r = (\n",
" np.absolute(farfields_r[\"Ex\"]) ** 2\n",
" + np.absolute(farfields_r[\"Ey\"]) ** 2\n",
" + np.absolute(farfields_r[\"Ez\"]) ** 2\n",
")\n",
"intensity_z = (\n",
" np.absolute(farfields_z[\"Ex\"]) ** 2\n",
" + np.absolute(farfields_z[\"Ey\"]) ** 2\n",
" + np.absolute(farfields_z[\"Ez\"]) ** 2\n",
")\n",
"\n",
"# Plot the intensity data and save the result to disk. \n",
"fig, ax = plt.subplots(ncols=2)\n",
"\n",
"ax[0].semilogy(np.linspace(0, size_r_um - pml_um, intensity_r.size), intensity_r, \"bo-\")\n",
"ax[0].set_xlim(-2, 20)\n",
"ax[0].set_xticks(np.arange(0, 25, 5))\n",
"ax[0].grid(True, axis=\"y\", which=\"both\", ls=\"-\")\n",
"ax[0].set_xlabel(r\"$r$ coordinate (μm)\")\n",
"ax[0].set_ylabel(r\"energy density of far fields, |E|$^2$\")\n",
"\n",
"ax[1].semilogy(\n",
" np.linspace(\n",
" focal_length_um - 0.5 * scan_length_z_um,\n",
" focal_length_um + 0.5 * scan_length_z_um,\n",
" intensity_z.size,\n",
" ),\n",
" intensity_z,\n",
" \"bo-\",\n",
")\n",
"ax[1].grid(True, axis=\"y\", which=\"both\", ls=\"-\")\n",
"ax[1].set_xlabel(r\"$z$ coordinate (μm)\")\n",
"ax[1].set_ylabel(r\"energy density of far fields, |E|$^2$\")\n",
"\n",
"fig.suptitle(f\"binary-phase zone plate with focal length $z$ = {focal_length_um} μm\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the volume specified in `get_farfields` via `center` and `size` is in cylindrical coordinates. These points must therefore lie in the $\\phi = 0$ ($rz = xz$) plane. The fields $E$ and $H$ returned by `get_farfields` can be thought of as either cylindrical ($r$,$\\phi$,$z$) or Cartesian ($x$,$y$,$z$) coordinates since these are the same in the $\\phi = 0$ plane (i.e., $E_r=E_x$ and $E_\\phi=E_y$). Also, `get_farfields` tends to gradually *slow down* as the far-field point gets closer to the near-field monitor. This performance degradation is unavoidable and is due to the larger number of $\\phi$ integration points required for accurate convergence of the integral involving the Green's function which diverges as the evaluation point approaches the source point.\n",
"\n",
"Shown below is the far-field energy-density profile around the focal length for both the *r* and *z* coordinate directions for three lens designs with $N$ of 25, 50, and 100. The focus becomes sharper with increasing $N$ due to the enhanced constructive interference of the diffracted beam. As the number of zones $N$ increases, the size of the focal spot (full width at half maximum) at $z = 200$ μm decreases as $1/\\sqrt{N}$ (see eq. 17 of the [reference](http://zoneplate.lbl.gov/theory)). This means that doubling the resolution (halving the spot width) requires quadrupling the number of zones.\n",
"\n",
""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}