diff --git a/docs/examples/Tutorial_Wave_Propagation_ASM_Development.ipynb b/docs/examples/Tutorial_Wave_Propagation_ASM_Development.ipynb new file mode 100644 index 00000000..a51c5534 --- /dev/null +++ b/docs/examples/Tutorial_Wave_Propagation_ASM_Development.ipynb @@ -0,0 +1,549 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "title", + "metadata": {}, + "source": [ + "# Tutorial - Wave Propagation (Angular Spectrum Method)\n", + "\n", + "This tutorial demonstrates a first implementation of wave propagation in Optiland using the **Angular Spectrum Method (ASM)**.\n", + "\n", + "In this notebook we will:\n", + "- Build a simple lens system\n", + "- Configure fields and wavelengths\n", + "- Compute a propagated complex field at multiple `z` planes\n", + "- Visualize normalized intensity along propagation\n", + "\n", + "> Note: this notebook is intended as a development tutorial. The API and output tensor layout may change." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "imports", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import ipywidgets as widgets\n", + "\n", + "from optiland.optic import Optic\n", + "from optiland.physical_apertures import PolygonAperture\n", + "from optiland.solves import QuickFocusSolve\n", + "import optiland.backend as be" + ] + }, + { + "cell_type": "markdown", + "id": "backend", + "metadata": {}, + "source": [ + "## 1. Backend selection\n", + "\n", + "Optiland supports multiple backends. Here we use **PyTorch** to exercise GPU and complex FFT paths." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "set-backend", + "metadata": {}, + "outputs": [], + "source": [ + "be.set_backend(\"torch\")" + ] + }, + { + "cell_type": "markdown", + "id": "build-lens", + "metadata": {}, + "source": [ + "## 2. Build an example lens\n", + "\n", + "We build a multi-surface lens and define:\n", + "- An entrance pupil (EPD)\n", + "- Two fields (angles)\n", + "- Two wavelengths (one primary)\n", + "\n", + "Then we run a quick focus solve and visualize the layout." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "lens-definition", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "╒════╤═════════════════╤═══════════╤═══════════╤═════════════╤════════════╤═════════╤═════════════════╕\n", + "│ │ Type │ Comment │ Radius │ Thickness │ Material │ Conic │ Semi-aperture │\n", + "╞════╪═════════════════╪═══════════╪═══════════╪═════════════╪════════════╪═════════╪═════════════════╡\n", + "│ 0 │ Planar │ │ inf │ inf │ Air │ 0 │ 2.30501 │\n", + "│ 1 │ Standard │ │ 22.0136 │ 3.25896 │ SK16 │ 0 │ 2.30501 │\n", + "│ 2 │ Standard │ │ -435.76 │ 6.00755 │ Air │ 0 │ 2.03363 │\n", + "│ 3 │ Standard │ │ -22.2133 │ 0.99997 │ F2 │ 0 │ 1.20443 │\n", + "│ 4 │ Stop - Standard │ │ 20.2919 │ 4.75041 │ Air │ 0 │ 1.14025 │\n", + "│ 5 │ Standard │ │ 79.6836 │ 2.95208 │ SK16 │ 0 │ 1.68566 │\n", + "│ 6 │ Standard │ │ -18.3953 │ 42.3963 │ Air │ 0 │ 1.87058 │\n", + "│ 7 │ Planar │ │ inf │ nan │ Air │ 0 │ 3.49561 │\n", + "╘════╧═════════════════╧═══════════╧═══════════╧═════════════╧════════════╧═════════╧═════════════════╛\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "(
, )" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAACUCAYAAABhn+PEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMAlJREFUeJzt3Xl4U1XCP/DvzZ423ReogEWkgAiUgSoDCOqIDjh0FH3dwHEGRYUBnlcHeXmAV2CccUE2fYGRTRlZBJHl9+g4wui8isv8RinooIAL/ASBAaErTZv9nt8fadKkSdoE0t572+/nefIkd8nNSU4D95tz7jmSEEKAiIiIiIiIkkKndAGIiIiIiIjaE4YsIiIiIiKiJGLIIiIiIiIiSiKGLCIiIiIioiRiyCIiIiIiIkoihiwiIiIiIqIkYsgiIiIiIiJKIoPSBVA7WZZht9thMpkgSZLSxSEiIiIiIoUIIeB2u2Gz2aDTxW6vYshqgd1ux7Jly5QuBhERERERqcTjjz+O9PT0mNsZslpgMpkA+D9Is9mscGmUJ4SA3W6HzWZjy54Gsf60i3WnXaw77WLdaRfrTrvUXnculwvLli0LZoRYGLJaEKhcs9nMkAX/H77H44HZbFblHz41j/WnXaw77WLdaRfrTrtYd9rVtO4CywBgNBpVU58tlYMhi4iIiIiIlCEE4KkHHNWAowpwVMFQdQaAE3BUw1NXhWf/r3/X2bNnt9iCpBYMWUREREREdGl8nsag5Gy4j7YcbZvPHTyMBMAaWDBnAJZsAOPa8I0kB0MWEREREREBsgy4LsQXjBzV4ctue/RjGiyANQuwZPrvrZlATs/w5ZDtwpIBu2yCLacAkt4IuN3As8+2+ltPNoYsIiIiIqL2QgjA42i+NSnWsrMGEHLkMSW9PwyFBqO0zkBenyhBqcljozXyeC2V324HdNqOKdouPRERERFRe+Tz+ENPokHJURXW/S6MOT0kCGX6H2d0i2hNilg2pwEqGXBCKxiyiIiIiIhagxCN3e/ivl6p4eaujX5MvTkyCGX3aDkoWTIAPU/92wo/aSIiIiKi5ngciQ3kELotavc7XWTXOlsnIK9380HpYrrfkSIYsoiIiIio/fN5G1uKEr1eyeeKfkxTWnjXO2sWkNElxvVJIaHJlAbodK38hklJmglZ69evx/vvv4/jx4/DbDZjwIABmD59Orp37x7zOW+99RZ+//vfh60zmUz4xz/+0cqlJSIiIqKkEwJw1cYVjKx1FYD7AuBouK4pZvc7U2QQyuoOFAyM3ZoU7H5nbIM3TVqkmZB14MAB3HXXXejbty98Ph9WrlyJadOm4Y033oDVGrvZNDU1FTt27Aguq2WWaCIiIqIOy+NMbCCH0OuVhC/KAaWI0e+ENRfIjzb6XZNlo5WDOlDSaSZkLV++PGx5wYIFuPnmm3HkyBEMGjQo5vMkSUJubm5rF4+IiIioY/F5/aPfJTrxrKMK8DqjH9NkaxKEMoG0ghaGCc/yj5oX2v1OCDjtdthsNgYoUoRmQlZTdrt/wrP09PRm93M4HBg7diyEEOjduzemTp2KK6+8Mub+brcbbnfjsJcul78PrhACQogklFzbAp8DPwttYv1pF+tOu1h32tUh6k4Ifze6aJPMhoWlyIloJdeF6IfUGRtCUVZj61JmN6Bz/zi635ku7b0EH3aAumunmtadiFKvSor39TUZsmRZxpIlS1BcXIyePXvG3K+wsBBPPvkkioqKYLfbsWnTJjz44IPYtm0bOnXqFPU569evx9q1a4PLBoMBY8aMgd1uh8fjSfp70RohBBwOBwB2vdQi1p92se60i3WnXZqqO68TkrMGkqsGkqMacNVAclb71zmrGu6r/dsb1iOwPUr3OwEJsKRDmDMhLJkQlgz/ffrl/nuzfxmB9cH7TMBwEd3vBACHG0CM+Z0SpKm6ozBN6y70/Luuri6sMUQJgQaYlmgyZC1cuBDHjh3DunXrmt1vwIABGDBgQHC5uLgY//Ef/4GdO3diypQpUZ8zceJETJgwIbjscrmwYsUK2Gw2mM3m5LwBDQuk9/LycsV/SQhITU2NGZopXKDObDYb/9PRGNaddrHutKvN6072xdH9LnR940S1ktcR/T0YU2OPftfQgiRCW5SCcyqlA5IOEgAt/tXye6ddTesuNFSlpqbCZLqE1s4kMBrjG+xEcyFr4cKF+Pjjj7FmzZqET6wNBgN69+6NkydPxtzHZDKFVV7gg5QkiV/SBrIsY9OmTaoJWZIkYcqUKcjLy1O6KJoQ+Fvm37P2sO60i3WnXQnXnRCA257AQA5Vjd3zXDXRj6kzRnaxy7w8sstdlEEdJIOyJ6RK4vdOu0LrLrT+1FCf8b6+ZkKWEALPP/88PvjgA6xevRpdunRJ+Bg+nw9Hjx7F8OHDW6GEHUegCfenP/3pRdVDMtntduzZswf19fWKloOIiNoZrwtwVEFXcRqodCc2v5LsjXJAf/e7iIEbsq5oefQ7UyoHbyDSGM2ErIULF2L37t1YsmQJUlJSUF5eDsDflGixWAAA8+bNQ35+PqZNmwYAWLt2Lfr374+uXbvCbrdjw4YNOHv2LG6//Xal3ka7EAhZeXl5SEtLU7QsshxlFnUiIiIgvu53sYKTpx4SgNSmxzSmRAah3N6NXfJiBSVLBqDTt8W7JiIViCtkLV26NOEDP/TQQ8jIyEj4ebFs374dAPDoo4+GrZ8/fz5KS0sBAGfPnoUuZPjOCxcu4I9//CMqKiqQnp6OPn364OWXX0aPHj2SVq6OKBCympufjIiIKCmEANx1ic2nFFh2XoB/RIUmdIbIIJTeBejULywYCWsmHMIMa3YBJGu2f5uB12cTUcviCllbtmxB//79477Q64svvsDdd9+d1JBVVlbW4j5r1qwJW54xYwZmzJiRtDKQn9Ppn9si0IKoBmq5PoyUdejQIVRWVir2+hkZGWGD7RBRCK87saAU+liOMbqvOSOyBSmre+zWpMCyKc65k4SAz24HONcSESUo7u6CixcvRnZ2dlz7jhw58qILROrncDhgsVjCWg2VovTFj6QeFRUV2L59O0wmkyJ/m06nE7m5uQxZ1L7Jsn9whngHcghd9tRFP6bBGhmEcosir12yZgKWJnMqsfsdEalUXCFr/vz5/hmz4zRnzhzk5ORcdKFI3Vwul+qGs2dLFgWuzxs1ahTy8/Pb/PUPHjyIw4cPt/nrEiVMCMBTn1i3u2CIqkHU7neSPjIopRcAnfo236JkyQSM6ukVQUSULHGFrLFjxyZ00NGjR19UYUgbnE6nakIWW7ISc+HCBezduxdWqzXssyspKUFmZqZyBUsCpYO21WqF0+mEz+eDXs9f16kN+DwXF5QcVc10v0uPHAo8qzC8RSlaaDKnsTsdEVGISxpdsL6+PmJ0t0RavCgxH3xzDqeqok822HYEvjjpgdORAembKhh1Eox6HUx6CQa9BKNOgkkvwdiwbNLrGvYJuemSP8eB0ifYWmG1WvGf//mfOH78eHBdt27d8N133ylXqCTx+XwAoFg31sA1inV1dUhPT1ekDKRBsgy4LsQZlGrCt7nt0Y9psER2tcu+suVhwi0ZgF4zgw4TEalawv+anj59Gs8//zz2798fNgOzEAKSJOGzzz5LagGp0ZbPfsB7R84pXQwI2QAZ6fjLj6cv+hiGhuBl0kthj416HQyhQS3ksTFKYBNeD77xdMaoC25wzMiWGY1GzJw5E1OnTg2umzNnjmpaJi+F0iErMNpmfX09Q1ZHIwTgcTQTlKphvnAO8NVFGf2uBhBRpqKQdJFByNYJyOvTfFCyZgJGjvxKRKS0hEPWk08+CcA/J1V2dja7a7Wh1b8qUboIEEJgw4YNSElJxU8GD4bHJ+D2CXh8Al458FiGR/avC94alt0+ucm+At7A+hjP8fhkODwy3D4vvMHjCHhlGS6vDLu3E87ZY3R9oQj33nsvFi9ejO+//x7dunXDxIkTlS5SUgRa1ZUKWYGgyomxNSzQ/S4w6Wwi3fB87ujHNKUB1kzoTemALccfhDK6xb4+KTj6XRqggsGFiIjo4iQcsr777jts3LgR3bt3b4XikBZ4PB6YzSZ/y5IeSFGwLHa7Hdu3b8eAgv4KlkJbjEYj5s6di0mTJrWbViyAIYsayDLgrk3g+qTqxuVY3e/05sgWo+wrAOugFgZ1yAD0RkAI1Nvt/u70/GGSiKhDSDhk9e3bFz/++CNDVgfm8XjinjON2ta//vUvuFyuFve76qqrMHDgQIwbN64NStU2AiFLqdZ1o9EISZKC88jRJQp0v0uoRanafx+z+11GeBCy5QN5vZsPStYsdr8jIqKEXVR3wWeeeQbnz5/HlVdeCYMh/BBFRUVJKxypj8/ng8/nU03ICpxQc+ALv7179+LChQsxt0uShKysLFRVVWHgwIHYuXMnpkyZ0oYlbD1Kt2RJkgSTyQSHQ+nBaVTE542j612Mbb4YPxaYbE0GdcgEMrq0PEy4OZ3d74iIqM0kHLKqqqpw+vRp/P73vw+ukySJA190EIHBTtQSsijcnXfeCa/XG3O7ECLYErl06dJ2dU1lIGgr+Z5MJlNcLYmaIkTD6HfVLQej4HLDZLXu2ujH1Jsig1BWd6BgYMuDOuj5bw8REalfwiHrqaeeQq9evfD0009z4IsOKBCymrZgkvaoqUUyGZTuLgj4f3xQbcjyOMKvP4q7G14NIHxRDiiFDNaQ6Q9CKblATlELQamh+x3/7yAionYs4TPlM2fOYOnSpejWrVtrlIdUzuPxj+J3ySFLCEB4IfnckHwuSN6Ge5+rcZ3siVwX3M9/b3PaMUYchrFmGICel/4GOwghBLxeb7sKy2oJWaFTWySd7AXqKxrmS0rweiVvjGvFjKlNglAGkN+35aDE7ndEREQxJXyGdc011+Dbb79lyGrvhADcgTldGm/Gs8fxk/p/oNuxb5F6ytAYggKBKCQswRceiBrDk9v/ONrF6S0VS2+C0JshdEYIvRk+nRGXox4eZ3XyP4N2TJZlCCHaVchSQ3dBvV4f/CEiJiEAV238AzmEPE5zxbjeTmeMDEKZhUBBcfPDhFsyAYMpGW+diIiIQiR8hjVixAgsXboUR48eRc+ePSNO0q6//vqkFY6SQJYBV02TsFQdEZ6irpcjTxYzIGGwZIbuBxsko8UfeALBR28CdP7HsjktYlvwXtd0nRnC0GSdzhTxXOiN/hHCQtTV1eGNN97A+PyBbfJxtheB67b0er3CJUmeNgtZQobkdUDntkPy1EHnsUPn9t/3ce6H1eUG/ve75kNUrO53lozwoJSSDeRcGQxGTskCc2ZnSNasJqPfpbD7HRERkYokHLKeffZZAMC6desitnHgi1ZWXwnUnW8hKDW5OWsARBl5T2f0n8CF/sqd3aPhxC3k5K3J8vEzlXj7r+9g5MiRyMjIaNO3Hw2vCbw4ga51HaolS/ZC56mH5LFD57ZD56lreOwPSOGBKcq24HPqIEX7TgHIA+CRTMAXeeEtRvl9Wh4m3JwO6JoJvULAY7fDzLmWiIiIVC/hM6x9+/a1RjkoHv9nCvDt7vB1wespQk7WMro2E5ayLumXb59cAUC5YbIptjNnzuCf//wnBg0a1GL9+Hz+lhRN1qMs+69Lsv8I1J0D7P5bp+OHcJs4iM57D0DvrYeuIRgFw1Osa5IACEgQxlTIxlTIJlvDYxtkUyp8qZ0gG1MgGxvWm2z+/Yw2CFPDfsZUCJMNH332BeqdbkyaNKkNPxAiIiJSm/bzM3ZHcOMc4LrHG4OSAtdTaPrkvJ2rqKjAV199heLi4hbrJ9CSpZrugkL4W17t50KC04/BAOVf17BcVx7Z3c6cjnRjOhyQIcnZkC3Z8KZdDtmUCtEQgpqGp9DQJAzWiK6oF0NnMMHn4zxZREREHd1FhaxDhw6hrKwMVVVVwZO1gN/97ndJKVgs27Ztw8aNG1FRUYGioiLMnDkT/fr1i7n/e++9h5deeglnzpxBt27dMH36dFx33XWtWsZWU1CsdAkUn/CVYgvUSdPvZDRtUo+BAR6ahqTQABXSEhVxDaDBCtjyAVsn/33Xa/yPU/Ma1nUCbHlAaj5gSsE3Bw9i165duP9n9yvWDVKSpOAPEURERNRxJXwm8sorr+Cll15CYWFhxDxZrX19zN/+9jcsW7YMs2fPRr9+/bBlyxZMnz4dO3bsQHZ2dsT+//rXvzB37lxMnToVI0aMwO7du/HEE09g06ZN6NmTw31fDDUMk03RBeokcG1Scy5pkAh3XfQWpmjrmnbR05vCg1Ln/o2BKRie8v03U2LXHsXzvlubTqeLK+QSERFR+5ZwyNq6dSvmzZuH0tLS1ihPszZv3ozbb78dv/zlLwEAs2fPxscff4w333wTv/nNbyL237p1K4YOHYoHHngAADBlyhR8+umn2LZtG+bMmdOWRW831HAiG41ay9WWLiZkBXldzQSnJq1ObnuTF9b7Q1EgJOX2BrpfF73VyZLZ6oM2KPkDgCRJ/FskIiKixEOWJEkoLm77bmsejwdff/01Jk6cGFyn0+lw7bXX4uDBg1Gfc/DgQUyYMCFs3dChQ/HBBx/EfB232x02majL5QLgPylV+uTJ4fbh2Hl7yzu2oqMVLvzoNuG7cidMJuV/sXc6nagQKfj2vAOulGqli6OoH85dgCw8OPv9IaTo3DB46mDw1sHgtftvHv9jq9uOFFcNRoujSH/3ffgOToHeHT7/koAEryUHXmsuvCl58FovgzdzALzWPHhT8uCx5vkfW3Phs2S1fD2Tq+GGmtZ6+zha4UKFSME35x3QK9Sd9VSdhHK3CV+eqk76sYUQqHc4kGL1siVZY1h32sW60y7WnXY1rTtvyPyT9S4vjEajgqWL/4f9hEPW+PHj8cYbb2DGjBkJF+pSVFdXw+fzRXQLzM7OxvHjx6M+p6KiIur+FRUVMV9n/fr1WLt2bXDZYDBgzJgxsNvtLU8y2soOn6nFPa98rmgZ/Lpg41unlC5EiL54861TANRUpuTQQUYOLiBPqkaeVOO/R8N9YB389/2levwCAD4NP4ZPSKhFCuywolakoLbh/pjcCQdqU/EPWw+cF5koRwbOiwycF5moRBp8Tj1Q3VzpyhtuatIXb+46ruDrWwB0xaYVnyhYBiIiovbDAB9+ZfU/PnTyPAZ0i7xEqC0FGmBaknDI+tWvfoXHHnsMt912G3r06BFxgfmiRYsSPaSqTJw4Maz1y+VyYcWKFbDZbDCbzQqWDOhfaMVb04YrWoajx47hwP79uHbIEJhNbTuyYTROpxPv/f093HLzLSgsLFS6OPERAjr3BRgd52GoPw+Do+FWf96/LnS9swKSCG8x9JqzGlqX8uCx9oY3JQ911jwcrZPwz0P/D/2vuR5SSia8Bht8hlT49JZgFz0hBLxeL+rr61H2+msYOXIkRt44SolPIem+/fZb7P1wL8aMuVWxlqzDhw+jvLwc99xzT9KP3fjLnpW/ymoM6067WHfaxbrTrqZ15/V48OYGfyPD1d3yYLNZFS1fvC1pCYesRYsWoaysDCUlJcjIyGizP9zMzEzo9XpUVlaGra+srEROTk7U5+Tk5CS0PwCYTCaYQsJD4IOUJEnxL2mK2YD+XTMVLYOuxoLTJjd65VpgsVgULQsA1NcDn0v16J2fgiKFPxt4HE2uY/ox/Jqm0Htfk19BjCmN1y6l5QMFvUIGgQi5T82DwWCK+sWtPnIE5w9vQ37PQTHrRggBj8cDu92HDMmJwiyL4n9TySJXmPGVVI9euRbFRhesOS7DYPa0ymcqhIDdboDNZlP83yJKDOtOu1h32sW6066mded2u/Fmw7YUs0Hx+oz39RM+E3n77bfx/PPPt/kw6EajEX369MFnn32GG264AYB/pLt9+/bh7rvvjvqcAQMGYN++fRg/fnxw3aeffor+/fu3RZHbpcCQ30pfn9ZmZJ9/XqawoBQjPLmaXG+kM/iHFw8EpPyrgB43RA9PZtslFzWREQMD+7SnkfCU/kcX8H+eqpl7jIiIiBSTcMhKT09H165dW6MsLZowYQIWLFiAvn374uqrr8Zrr70Gh8MRHOlw3rx5yM/Px7Rp0wAA9957Lx555BFs2rQJ1113Hfbs2YPDhw9zZMFLkMhcTKolBOCsia/Fqb4caNJdD9asxoCUXgBcNjAyNNk6AdZsoA27rSUy91W7qEcVkmWZc8gRERFR4iHrkUcewerVqzF//vw27y52yy23oKqqCqtWrUJFRQV69eqF5cuXB7v/nT17NuwEp7i4GE8//TT+9Kc/YeXKlejWrRsWL17MObIuQeBXelWenHuc8bU4ReuuZ7ACaZ0ag9LlQ2J214NB2WvzYglMgttRQ1YiQ9i3FrZkEREREXCR82SdPn0at9xyCwoKCiKufdi8eXPSChfNPffcE/Oi8jVr1kSsGzVqFEaNah8X9qtBq4UsISB5HdB57JDcdug8dujcDbfgulro3HVh+8B5AVPFGWRtWwN4YszfFOyu1wfocX1ki9NFTHyrRomELFWH5Yukhu6CPp9PsevBiIiISD0SPhsIXA9FHVPgBDJwQg8AkD0R4UfnroXOU9ckMIWuaxqY6iAJX4xXBWRjqv9mSoNstEGYUiGbMuA25+PbSiuK+g1HXverFe2up7RevXph7NixcYWNQMjyer2tXaw2o4aWLJ/PFzZwDhEREXVMF9VdkNohrxtwVgOOqia38HV51WdxV/X3yHjvJRhkpz88+ZwxDyt0BsjGNMim1IZwlAbZmAqvrQtkUxqEMRWyyRayT5o/QBltDYEqFcKYCuiid8Gqq6vDu8ffQF7f8cgrKmqdz0YjbDYbOnfuHFdwCoSssLCscWoIWV6vF6mpqYq9PhEREakD+7W0J0IAnvqYASlqgAoEK7c9+jGNqf6BHqxZgDUTkikV5frOkPOvgDUzPzwMmRofB1qcoDe1STc8NXQV0xKdTuefe4ItWUnl9XrZkkVERETxhayf/exn2LlzJzIzM+M66C9+8QusW7cOBQUFl1I2aurrt4EfD/tDUdRWpyrA547+XEtGSFjK8nery+vduGzJDN/eEKqaDvLgrKnB+xs34ic9f6LYKJN06SRJgsFggMfjUbooSaOG6QW8Xm/ckxQSERFR+xVXyKqtrcUnn3wCmy2+uXxqamraVTck1fhyO3D8o/AglFkIFAyMEpAaQpI1yx+wYnS3S1TgV/r21ALSnvh8vmbrRggBn88XnFybISu53G43zGZ1jj5JREREbSfu7oILFixoxWJQXO5ar3QJgiFLLSfnHWZS5Djt3LkTNTU1MbdLkoTs7GxUVlaitra2Xc3ppIbugh6Pp82ntiAiIiL1iStk7du3r7XLQRqh0+mg1+tVE7ICeE2W35gxY+ByuWJuF0JACAG3243f/OY3eP7559uwdK1L6bm/Ap8rQxYRERFx4AtKmNFohNsd49ovUlTv3r2b3S6EgN1ux6uvvoojR45gy5YteOyxx9qmcK1M6e6CHo8HQghYrVZFXp+IiIjUgyFLQ5weH3yyst3jhBCQdSbUODyorPfC45Ph8Qm4fQJe2X/v8ckhjxtucuCxHNw3+LyGe48sN+7f5DmR6/yP3V4Zbt9g9D1dhyuvVPSj0Qyn04nnnnsOALBw4UI8+uij7SIYKN2S5XT6pzJoD58lERERXRqGLA2ZvuVzvHv4R6WLASDHf/fF1wk9Sy8BRr0Eo14Ho05qeCyFP9brgutSjDoYLYF9dDDpJRj0EkwN+xp0EmSvG1/96wtcntmx58hKxIYNG3D69GkAwNmzZ7F69ep20ZoVmPtLqZAV6KaZkpKiyOsTERGResQdss6fP4+8vLzWLAu14JGRPfDL4ssULoXAvn1lcNTZMWjggOihqSEIGXS6kMcS9LrkXzdVW1sL76HzyLNx2Ox4OJ1OLF26NGxde2nNUktLFicjJiIiorhD1t13341Zs2Zh9OjRrVkeasY13bOVLgKEENCfNuKHH2oxrHu60sUJXn/DgS/i4/V68dprryElJSXsM2sPQ/Ir3ZLlcDgAsCWLiIiIEghZv/3tb/HMM8/g/fffx5w5c5CRkdGa5SIVs1gswV/tSVtsNhsGDRoEm83W7oKp0u/H4XDAarUGwx4RERF1XHGHrLvuugvDhg3DH/7wB9x9992YO3cuRo4c2ZplI5Uym83NDhOuBKVPsEl5BoP/n7N33nlHkb8HWZbZpZqIiIgAJDjwRZcuXbBq1Sq8/vrrmDlzJq644oqIX203b96c1AKS+lgsFrhcLsiyrPhktsLnQaqoA3wcUr6jy8rKwoQJE1BVVaVYGdLS0hR7bSIiIlKPhEcXPHPmDN5//32kp6fj+uuvZ9eYDigwQILT6Uz+9SeyFzpnNfTOSuidldA5q0Luq6B3VELnqoTeWQWdoxKF7gsYAOBsxVDgyl7JLQtpTs+ePZUuAhEREVFiIWvXrl144YUXcO2112Lbtm3IyspqrXKRigVClsPhaDlkJRia9O4LEYcQehN8lmz4LFmQLdnw2rrAndsfPks27MKCjw8cwXWZnCSLiIiIiNQh7pA1ffp0HDp0CDNnzsTYsWNbs0wR/v3vf2PdunUoKytDRUUFcnNzceutt+LBBx+E0Rh76O5HHnkEBw4cCFt3xx13YM6cOa1d5PZJCMBTj1RfDfLFeTgP74GcYYTBVQWDuwYGVzWMrmoY3A03VzUMntqIw8g6E7zmTHhNGfCYs+A15cKb2hNecxa85gx4TZnwmDPhNWXCa86ErLcCMa6xsdvt+FZyY5iZA7EQERERkTrEHbJ8Ph+2bNmCTp06tWZ5ojp+/DiEEJgzZw66du2KY8eO4emnn4bD4WhxEtVx48bh0UcfDS5bLJZWLq0GyDLgugA4qvw3Z3XjY0cV4Gi63LhO8rmQD2AyABz1H84DPeqRgjpYUY0U1MOKeuSjDoWoa1iua9hejxS4ZSPglICoAxT6AFQ03OKj1+s5NxERERERqUbcIetPf/pTa5ajWcOGDcOwYcOCy127dsWJEyewY8eOFkOWxWJBbm5uK5dQIT5vlIDUUlhqCFUiylxCOgNgzQq/ZRYCBQODy8KaCQcsqJONkC3ZkC3ZEMaUYEtTasOtLVmtVnZdJSIiIiLVSHjgC7Ww2+1IT295Mtx33nkHf/3rX5GTk4ORI0di0qRJzbZmud1uuN2NI9UFhioXQgQnvlXMewuAY//bEJRqILkir18CAGGwNoSizMawlH9V42NLyPrQfUy2mN3ygscWAl67HTkqm2dJ8brRiMDfMT8v7WHdaRfrTrtYd9rFutOupnUXWodqqNN4X1+TIevkyZN4/fXXW2zFGj16NAoKCpCXl4fvvvsOy5cvx4kTJ7Bo0aKYz1m/fj3Wrl0bXDYYDBgzZgzsdjs8Hk+y3sJFMZoyoetUDGHJDN4QfJzhvzdnAEZr4gf3APDUtbibEAIOhwMA56bSItafdrHutIt1p12sO+1i3WlX07oLPf+uq6sLawxRQrxzxSoaspYvX45XX3212X22b9+O7t27B5fPnTuH6dOnY9SoURg3blyzz73jjjuCj3v27Inc3FxMmTIFp06dQteuXaM+Z+LEiZgwYUJw2eVyYcWKFbDZbDCbzXG8q1Y08jFlXx+N6d2mspYsig/rT7tYd9rFutMu1p12se60q2ndhYaq1NRUmEwmpYoGAM0OuhdK0ZB1//33o7S0tNl9unTpEnx8/vx5TJ48GQMGDMDcuXMTfr1+/foB8LeExQpZJpMprPICH6QkSfySNgh8Fvw8tIn1p12sO+1i3WkX6067WHfaFVp3ofWnhvqM9/UVDVlZWVlxD1hw7tw5TJ48GX369MH8+fOh0+kSfr1vvvkGANrvQBhERERERKS4xJOKAs6dO4dHH30UnTt3xmOPPYaqqiqUl5ejvLw8bJ8777wTX331FQDg1KlTWLduHY4cOYJ///vf2Lt3L+bPn49BgwahqKhIqbdCRERERETtnCYGvvj0009x8uRJnDx5ErfeemvYtrKyMgCA1+vFiRMn4HT6J18yGAz47LPPsGXLFjgcDnTq1Ak/+9nP8NBDD7V5+YmIiIiIqOPQRMgqLS1t8dqtyy67LBi4AKBz585Ys2ZNaxeNiIiIiIgojCa6CxIREREREWkFQxYREREREVESMWQRERERERElEUMWERERERFREjFkERERERERJZEmRhckIiIiIqKOx2g0Yvbs2cHHWsGQRUREREREqiRJEkwmk9LFSBhDVguEEAAAl8ulcEnUQQgBl8sFo9EISZKULg4liPWnXaw77WLdaRfrTrtYd9ql9roLZIJARoiFIasFbrcbALBs2TKFS0JERERERGrgdrthsVhibpdESzGsg5NlGXa7HSaTSZVpuq3Z7Xb84he/wNtvvw2bzaZ0cShBrD/tYt1pF+tOu1h32sW60y61150QAm63GzabDTpd7DEE2ZLVAp1Oh/T0dKWLoRoejwderxdmsxlms1np4lCCWH/axbrTLtaddrHutIt1p11aqLvmWrACOIQ7ERERERFREjFkERERERERJRFDFiXEZDLh4Ycf1uRQmsT60zLWnXax7rSLdaddrDvtai91x4EviIiIiIiIkogtWUREREREREnEkEVERERERJREDFlERERERERJxJBFRERERESURAxZlJBt27ahtLQUw4YNw69//Wt89dVXSheJmjhw4AAef/xxjB49GiUlJfjggw/CtgshsGrVKvz85z/H8OHD8dvf/hY//PCDMoWlMOvXr8cDDzyAkSNH4uabb8aMGTNw/PjxsH1cLhcWLlyIm266CSNGjMDMmTNRUVGhTIEpaPv27bj33ntx/fXX4/rrr8fEiRPxySefBLez3rTjz3/+M0pKSrBkyZLgOtafeq1evRolJSVhtzvvvDO4nXWnbufOncOTTz6Jm266CcOHD8c999yDw4cPB7dr+ZyFIYvi9re//Q3Lli3Dww8/jE2bNqFXr16YPn06KisrlS4ahXA4HCgqKsKsWbOibn/11VexdetWzJ49G3/+859hsVgwffp0uFyuNi4pNXXgwAHcddddWL9+PVauXAmv14tp06bB4XAE91m6dCk+/PBDPPfcc1izZg3Ky8sxc+ZMBUtNAJCfn49p06Zh48aN2LBhA0pKSjBjxgwcO3YMAOtNKw4dOoSdO3eiqKgobD3rT9169OiB3bt3B28vv/xycBvrTr0uXLiAhx56CAaDAS+++CK2bduGxx9/HOnp6cF9NH3OIoji9MADD4jnnnsuuOzz+cTo0aPF+vXrlSsUNWvw4MHi/fffDy7LsixuueUWsWHDhuC62tpaMXToULF7924FSkjNqaysFIMHDxb79+8XQvjrasiQIeLdd98N7vP999+LwYMHi4MHDypVTIrhxhtvFLt27WK9aURdXZ0YN26c+Oc//ykefvhhsXjxYiEEv3dqt2rVKnHfffdF3ca6U7f/+Z//EQ899FDM7Vo/Z2FLFsXF4/Hg66+/xpAhQ4LrdDodrr32Whw8eFDBklEiTp8+jYqKClx77bXBdTabDf369cOXX36pYMkoGrvdDgDBX/WOHDkCr9cb9j3s3r07OnfuzO+hivh8PuzZswcOhwMDBgxgvWnEwoULMXz48LB6Avi904IffvgBo0ePxm233Yb//u//xtmzZwGw7tTuww8/xFVXXYVZs2bh5ptvxvjx47Fr167gdq2fsxiULgBpQ3V1NXw+H7Kzs8PWZ2dnR1wzQuoV6Ieek5MTtj47O5t91FVGlmUsWbIExcXF6NmzJwB//RmNRqSlpYXty/pTh6NHj2LixIlwu92wWq1YtGgRevTogW+//Zb1pnJ79uzB119/jQ0bNkRs4/dO3fr164cFCxagsLAQ5eXlWLt2LSZNmoTXX3+ddadyp0+fxo4dOzBhwgRMnDgRhw8fxuLFi2E0GjF27FjNn7MwZBERqdDChQtx7NgxrFu3TumiUJwKCwvx2muvwW634+9//zsWLFiANWvWKF0sasHZs2exZMkSrFy5EmazWeniUIKGDx8efFxUVIR+/fph7NixePfdd2GxWBQsGbVElmX07dsXU6dOBQD06dMHx44dw44dOzB27FiFS3fp2F2Q4pKZmQm9Xh8xyEVlZWXELwykXoG6avoLEOtRXRYuXIiPP/4Yq1atQqdOnYLrc3Jy4PF4UFtbG7Y/608djEYjunXrhquuugrTpk1Dr169sGXLFtabyn399deorKzE/fffjyFDhmDIkCE4cOAAtm7diiFDhiA7O5v1pyFpaWkoLCzEqVOn+N1TudzcXFxxxRVh66644opgd0+tn7MwZFFcjEYj+vTpg88++yy4TpZl7Nu3DwMGDFCwZJSILl26ICcnB/v27Quus9vt+Oqrr9C/f38FS0aAf6jahQsX4oMPPsBLL72ELl26hG2/6qqrYDAYwr6Hx48fx9mzZ/k9VCFZluHxeFhvKnfNNddg69at2Lx5c/DWt29fjB49OviY9acd9fX1OHXqFHJzc/ndU7ni4mKcOHEibN2JEydQUFAAQPvnLOwuSHGbMGECFixYgL59++Lqq6/Ga6+9BofDgdLSUqWLRiHq6+tx8uTJ4PLp06fxzTffICMjA507d8Z9992Hl19+Gd26dUOXLl3w0ksvIS8vDzfccINyhSYA/has3bt3Y8mSJUhJSUF5eTkA/4W+FosFNpsNt912G5YtW4aMjAykpqZi0aJFGDBggCb+w2nPVqxYgWHDhqFz586or6/H7t27sX//fixfvpz1pnKpqanB6x4DLBYLMjMzg+tZf+r1wgsvYMSIESgoKMD58+exevVq6HQ6/PznP+d3T+XGjx+PBx98EK+88gpuvvlmHDp0CLt27cLcuXMBAJIkafqcRRJCCKULQdrx+uuvY+PGjaioqECvXr0wc+ZM9OvXT+liUYiysjJMnjw5Yv3YsWOxYMECCCGwevVq7Nq1C7W1tRg4cCBmzZqFwsJCBUpLoUpKSqKunz9/fvDHDJfLhRdeeAF79uyB2+3G0KFDMWvWLOTm5rZlUamJp556Cvv27UN5eTlsNhuKiorwwAMP4Kc//SkA1pvWPPLII+jduzdmzJgBgPWnZrNnz8bnn3+OmpoaZGVlobi4GFOnTkXXrl0BsO7U7qOPPsKKFStw8uRJXHbZZZgwYQLGjRsX3K7lcxaGLCIiIiIioiTiNVlERERERERJxJBFRERERESURAxZREREREREScSQRURERERElEQMWUREREREREnEkEVERERERJREDFlERERERERJxJBFRERERESURAxZRETUYZWWlqKkpAQlJSWora1t89cvKysLvv6MGTPa/PWJiKh1GJQuABER0aUoKyvD5MmTY24fPHgwVq9eHXP75MmTcfvtt8Nms7VG8ZpVXFyM3bt3Y8mSJXC73W3++kRE1DoYsoiISNMCQaWpDz/8EM8++yzuuuuuZp+fkpKC3Nzc1ipes4xGI3Jzc2E2mxmyiIjaEXYXJCIiTQsEldBbbW0tXnzxRUycOBGjRo1K6HhvvfUWbrjhBnz00Ue44447MHz4cPzXf/0XnE4n/vKXv6C0tBQ33ngjFi1aBJ/PF3xeaWkp1q1bh3nz5mHEiBEYO3Ys9u7di6qqKvzud7/DiBEjcO+99+Lw4cPJ/giIiEhlGLKIiKhdqa2txYwZMzBo0CBMmTLloo7hdDqxdetWPPPMM1i+fDn279+PJ554Ap988glefPFFPPXUU9i5cyf+/ve/hz1vy5YtKC4uxubNm3Hddddh3rx5mD9/Pm699VZs2rQJXbt2xfz58yGESMZbJSIilWLIIiKidkOWZcydOxd6vR5//OMfIUnSRR3H6/Vi9uzZ6NOnDwYNGoSbbroJX3zxBZ588kn06NEDI0aMQElJCcrKysKeN2zYMNx55524/PLLMWnSJNTV1aFv374YNWoUCgsL8etf/xrff/89KioqkvF2iYhIpRiyiIio3Vi5ciW+/PJLLFmyBKmpqRd9HIvFgq5duwaXc3JycNlllyElJSW4Ljs7G1VVVWHPKyoqCnsOAPTs2TPsOQAinkdERO0LB74gIqJ2Yc+ePdi0aRNeeOEFXH755Zd0LIMh8r/HpuskSYIsyzH3CbSiRVvX9HlERNS+sCWLiIg075tvvsEf/vAHTJs2DUOHDlW6OERE1MGxJYuIiDSturoaTzzxBAYPHowxY8agvLw8bLter0dWVpZCpSMioo6IIYuIiDTt448/xpkzZ3DmzBmMHj06YntBQQHeeustBUpGREQdlSQ4jiwREXVQpaWluO+++zB+/HhFy7FgwQLU1tZiyZIlipaDiIiSg9dkERFRh7Z8+XKMGDECdru9zV/7888/x4gRI/DOO++0+WsTEVHrYUsWERF1WGfOnIHX6wUAdOnSBTpd2/726HQ6cf78eQCA1WpFbm5um74+ERG1DoYsIiIiIiKiJGJ3QSIiIiIioiRiyCIiIiIiIkoihiwiIiIiIqIkYsgiIiIiIiJKIoYsIiIiIiKiJGLIIiIiIiIiSiKGLCIiIiIioiRiyCIiIiIiIkqi/w9twStYWD7cZwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAGwCAYAAACw8T6LAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMfRJREFUeJzt3XtY1HW+B/D3DJdBVAZYkIETiZdWIS+4mBPWqVznEYrHlbNuiaGoS6CWuoalcjalpCLNYx6VZ+1RET3Zau2qqbUUBzVPSZAoeQl9woN3Bi8cZgITufzOH22/GrnNwNy+M+/X8/ye5Dvf33e+v2l4z4fv7zczCkmSJBARkXCUjp4AERF1DwOciEhQDHAiIkExwImIBMUAJyISFAOciEhQDHAiIkF5OnoCrqC1tRXXrl1D3759oVAoHD0dIhKcJEn4/vvvERYWBqWy4zqbAW4F165dQ3h4uKOnQUQu5vLly7jvvvs6vJ0BbgV9+/YFAMTE/xmeXj4Ong0Ria656Q7KCt6Qs6UjDHAr+GnZxNPLhwFORFbT1ZIsT2ISEQmKAU5EJCgGOBGRoBjgRESCYoATEQmKAU5EJCgGOBGRoBjgRESCYoATEQmKAU5EJCihAvzIkSOYOHEiwsLCoFAosHfv3i73OXz4MH7zm99ApVJh8ODByM/Pb9MnNzcXERER8PHxgVarRWlpqfUnT0RkZUIFeENDA0aOHInc3Fyz+ldVVSEhIQHjxo1DeXk5Fi5ciOeeew6ffvqp3GfXrl3IyMhAVlYWjh8/jpEjRyIuLg7Xr1+31WEQEVmFQpIkydGT6A6FQoE9e/YgMTGxwz5LlizBxx9/jNOnT8ttSUlJqKurQ0FBAQBAq9XioYcewoYNGwD8+Nne4eHhmD9/PpYuXWrWXIxGI9RqNbQTs/lhVkTUY81Nd1CyfxkMBgP8/Pw67CdUBW6p4uJi6HQ6k7a4uDgUFxcDAO7evYuysjKTPkqlEjqdTu7TnsbGRhiNRpONSFSDllQ4egrUTS4d4Hq9HiEhISZtISEhMBqN+OGHH3Dz5k20tLS020ev13c4bk5ODtRqtbzxyxxIVD+FN0NcTC4d4LaSmZkJg8Egb5cvX3b0lIjIDbn0FzpoNBrU1NSYtNXU1MDPzw+9evWCh4cHPDw82u2j0Wg6HFelUkGlUtlkzkT2wqpbfC5dgcfGxqKoqMikrbCwELGxsQAAb29vxMTEmPRpbW1FUVGR3IfIFbUX3gx08QhVgdfX16OyslL+uaqqCuXl5QgMDMT999+PzMxMXL16Fdu3bwcAzJkzBxs2bMDixYvxxz/+EQcPHsQHH3yAjz/+WB4jIyMDM2bMwOjRozFmzBisXbsWDQ0NmDVrlt2Pj8jWGNKuRagAP3bsGMaNGyf/nJGRAQCYMWMG8vPzUV1djUuXLsm3DxgwAB9//DFefPFF/Od//ifuu+8+bN68GXFxcXKfKVOm4MaNG1i+fDn0ej2io6NRUFDQ5sQmEZGzEfY6cGfC68BJBOZU3+dXRtphJtQVXgdORDJzl064xCIWoZZQiMgyDGTXxgqciEhQDHAiF8Xq2/UxwIlcUE/Cm8EvDq6BE7kQhq97YQVO5CIY3u6HAU5EJCgGOJELsHb1zWpeDAxwIsExbN0XT2ISCYrBTazAiYgExQAnEpA9qm9W+M6PAU4kGAYr/YRr4ESCYHDTvViBEwmA4U3tYYATEQmKAU7k5BxZfbPyd25cAydyUgxP6gorcCInxPAmczDAiYgExQAncjLOVn0723zoZwxwIifCsCRL8CQmkRNgcFN3sAInIhIUA5zIwUSovkWYoztigBM5EIOReoJr4EQOwOAma2AFTmRnDG+yFgY4EZmFLzzOhwFOZEcMQbIm4QI8NzcXERER8PHxgVarRWlpaYd9n3jiCSgUijZbQkKC3GfmzJltbo+Pj7fHoZAbGbSkwiXC2xWOwZUIdRJz165dyMjIwMaNG6HVarF27VrExcXh3Llz6NevX5v+u3fvxt27d+Wfb926hZEjR+Lpp5826RcfH4+tW7fKP6tUKtsdBLkdhh7ZilAV+Jo1a5CWloZZs2YhKioKGzduhK+vL/Ly8trtHxgYCI1GI2+FhYXw9fVtE+AqlcqkX0BAgD0Oh4ioR4QJ8Lt376KsrAw6nU5uUyqV0Ol0KC4uNmuMLVu2ICkpCb179zZpP3z4MPr164chQ4Zg7ty5uHXrVqfjNDY2wmg0mmxE7WH1TbYkTIDfvHkTLS0tCAkJMWkPCQmBXq/vcv/S0lKcPn0azz33nEl7fHw8tm/fjqKiIqxcuRKff/45nnzySbS0tHQ4Vk5ODtRqtbyFh4d376DIpblqeLvqcYlIqDXwntiyZQuGDx+OMWPGmLQnJSXJ/x4+fDhGjBiBQYMG4fDhwxg/fny7Y2VmZiIjI0P+2Wg0MsRJxoAjexGmAg8KCoKHhwdqampM2mtqaqDRaDrdt6GhATt37kRqamqX9zNw4EAEBQWhsrKywz4qlQp+fn4mGxGRvQkT4N7e3oiJiUFRUZHc1traiqKiIsTGxna674cffojGxkZMmzaty/u5cuUKbt26hdDQ0B7PmdwPq2+yJ2ECHAAyMjKwadMmbNu2DRUVFZg7dy4aGhowa9YsAEBKSgoyMzPb7LdlyxYkJibiV7/6lUl7fX09Xn75ZXz11Ve4cOECioqKMGnSJAwePBhxcXF2OSZyHe4U3u50rM5MqDXwKVOm4MaNG1i+fDn0ej2io6NRUFAgn9i8dOkSlErT16Rz587hiy++wGeffdZmPA8PD5w8eRLbtm1DXV0dwsLCMGHCBGRnZ/NacDIbw4wcRSFJkuToSYjOaDRCrVZDOzEbnl4+jp4O2ZE7h/f5lZGOnoLLam66g5L9y2AwGDo9xybUEgoREf2MAU7UTe5cfQM8fmcg1Bo4kTNgcJGzYAVOZAGGNzkTBjgRkaAY4ERmYvVNzoYBTmQGhnf7+Lg4Fk9iEnWCAUXOjBU4EZGgGOBEHWD1Tc6OAU7UDoa3+fhYOQ7XwIl+gWFEImEFTvRPDG8SDQOciEhQDHAisPruKT5+jsEAJ7fH8CFR8SQmuS0GN4mOFTgRkaAY4OSWWH1bHx9T+2OAk9th0JCr4Bo4uQ0GN7kaVuBERIJigJNbYPVNrogBTi6P4W0/fKzti2vg5LIYJuTqWIGTS2J4kztggBMRCYoBTi6H1bdj8fG3HwY4uRSGB7kTnsQkl8DgJnfECpyISFDCBXhubi4iIiLg4+MDrVaL0tLSDvvm5+dDoVCYbD4+PiZ9JEnC8uXLERoail69ekGn0+G7776z9WGQFbH6dj78f2IfQgX4rl27kJGRgaysLBw/fhwjR45EXFwcrl+/3uE+fn5+qK6ulreLFy+a3L5q1SqsW7cOGzduRElJCXr37o24uDjcuXPH1odDVsCgIHcmVICvWbMGaWlpmDVrFqKiorBx40b4+voiLy+vw30UCgU0Go28hYSEyLdJkoS1a9filVdewaRJkzBixAhs374d165dw969e+1wRNRdg5ZUMLzJ7QkT4Hfv3kVZWRl0Op3cplQqodPpUFxc3OF+9fX16N+/P8LDwzFp0iScOXNGvq2qqgp6vd5kTLVaDa1W2+mYjY2NMBqNJhsRkb0JE+A3b95ES0uLSQUNACEhIdDr9e3uM2TIEOTl5eGjjz7Ce++9h9bWVowdOxZXrlwBAHk/S8YEgJycHKjVankLDw/vyaGRhVh5i4H/n2xPmADvjtjYWKSkpCA6OhqPP/44du/ejeDgYLz77rs9GjczMxMGg0HeLl++bKUZExGZT5gADwoKgoeHB2pqakzaa2pqoNFozBrDy8sLo0aNQmVlJQDI+1k6pkqlgp+fn8lG9sGqTiz8/2VbwgS4t7c3YmJiUFRUJLe1traiqKgIsbGxZo3R0tKCU6dOITQ0FAAwYMAAaDQakzGNRiNKSkrMHpOIyFGEeidmRkYGZsyYgdGjR2PMmDFYu3YtGhoaMGvWLABASkoK/uVf/gU5OTkAgBUrVuDhhx/G4MGDUVdXh7fffhsXL17Ec889B+DHK1QWLlyI119/HQ888AAGDBiAZcuWISwsDImJiY46TOoAqzkiU0IF+JQpU3Djxg0sX74cer0e0dHRKCgokE9CXrp0CUrlz39U/N///R/S0tKg1+sREBCAmJgYHD16FFFRUXKfxYsXo6GhAenp6airq8Ojjz6KgoKCNm/4ISJyNgpJkiRHT0J0RqPxx8sPJ2bD04vBbwusvsV2fmWko6cglOamOyjZvwwGg6HTc2zCrIETEZEpBjg5PVbfRO1jgBMRCYoBTkQ2x7+ibIMBTk6Nv/hEHWOAk9NieBN1jgFOTonhTdQ1BjgR2QVflK2PAU5EJCgGODkdVmpE5mGAk1NheBOZjwFORCQoBjgR2Q3/wrIuBjg5Df5yE1mGAU5OgeFNZDkGOBGRoBjgRGRX/GvLehjg5HD8hSbqHgY4EZGgGODkUKy+ibqPAU5EdscXbutggJPD8JeYqGcY4EREgmKAk0Ow+ibqOQY4ETkEX8R7jgFOdsdfXCLrYIATEQmKAU52xeqbyHoY4ETkMHxB7xkGOBGRoBjgZDestoisS7gAz83NRUREBHx8fKDValFaWtph302bNuFf//VfERAQgICAAOh0ujb9Z86cCYVCYbLFx8fb+jDcDsObyPqECvBdu3YhIyMDWVlZOH78OEaOHIm4uDhcv3693f6HDx/G1KlTcejQIRQXFyM8PBwTJkzA1atXTfrFx8ejurpa3v7617/a43DcBsObyDaECvA1a9YgLS0Ns2bNQlRUFDZu3AhfX1/k5eW123/Hjh14/vnnER0djaFDh2Lz5s1obW1FUVGRST+VSgWNRiNvAQEB9jgcIgJf4HtCmAC/e/cuysrKoNPp5DalUgmdTofi4mKzxrh9+zaampoQGBho0n748GH069cPQ4YMwdy5c3Hr1q1Ox2lsbITRaDTZiIjsTZgAv3nzJlpaWhASEmLSHhISAr1eb9YYS5YsQVhYmMmLQHx8PLZv346ioiKsXLkSn3/+OZ588km0tLR0OE5OTg7UarW8hYeHd++g3ACrKyLb8XT0BOzlrbfews6dO3H48GH4+PjI7UlJSfK/hw8fjhEjRmDQoEE4fPgwxo8f3+5YmZmZyMjIkH82Go0M8XYwvIlsS5gKPCgoCB4eHqipqTFpr6mpgUaj6XTf1atX46233sJnn32GESNGdNp34MCBCAoKQmVlZYd9VCoV/Pz8TDYi6j6+2HePMAHu7e2NmJgYkxOQP52QjI2N7XC/VatWITs7GwUFBRg9enSX93PlyhXcunULoaGhVpk3EZGtCBPgAJCRkYFNmzZh27ZtqKiowNy5c9HQ0IBZs2YBAFJSUpCZmSn3X7lyJZYtW4a8vDxERERAr9dDr9ejvr4eAFBfX4+XX34ZX331FS5cuICioiJMmjQJgwcPRlxcnEOO0VWwoiKyPaHWwKdMmYIbN25g+fLl0Ov1iI6ORkFBgXxi89KlS1Aqf35N+stf/oK7d+/iD3/4g8k4WVlZePXVV+Hh4YGTJ09i27ZtqKurQ1hYGCZMmIDs7GyoVCq7HpsrYXgT2YdCkiTJ0ZMQndFohFqthnZiNjy9fLrewcUxwKm7zq+MdPQUnEJz0x2U7F8Gg8HQ6Tk2oZZQiIjoZwxwsipW30T2wwAnIhIUA5yshtU39RSfQ5ZhgBMRCYoBTlbByomshc8l8zHAiYgExQCnHmPFROQYDHAiIkExwKlHWH2TLfB5ZR4GOBGRoBjg1G2skogcy6xPI/zlt8+Y65VXXmnz3ZNERGQ9ZgX42rVrERsbC29vb7MG/eKLLzBv3jwGOBF126AlFfx0wi6Y/Xnge/bsQb9+/czq27dv325PiMTA5RMixzNrDXzr1q1Qq9VmD/ruu++2+fZ4ch0MbyLnYFYFPmPGDIsGffbZZ7s1GXJ+DG8i59Gjq1Dq6+thNBpNNiIia2HB0DmLA7yqqgoJCQno3bs31Go1AgICEBAQAH9/fwQEBNhijkRE1A6Lv9R42rRpkCQJeXl5CAkJgUKhsMW8yAmxGiJyLhYH+DfffIOysjIMGTLEFvMhJ8XwJnI+Fi+hPPTQQ7h8+bIt5kJERBawuALfvHkz5syZg6tXr2LYsGHw8vIyuX3EiBFWmxwREd/Q0zGLA/zGjRs4f/48Zs2aJbcpFApIkgSFQoGWlharTpAcj8snRM7J4gD/4x//iFGjRuGvf/0rT2ISETmQxQF+8eJF7Nu3D4MHD7bFfMjJsPomcl4Wn8T87W9/i2+++cYWcyEiahcLifZZXIFPnDgRL774Ik6dOoXhw4e3OYn5u9/9zmqTI8fiLw2Rc7M4wOfMmQMAWLFiRZvbeBKTiMh+LA7w1tZWW8yDnAyrbyLnx69UIyIhsKhoq1sB/vXXX2PVqlV46aWXkJGRYbLZWm5uLiIiIuDj4wOtVovS0tJO+3/44YcYOnQofHx8MHz4cHzyyScmt0uShOXLlyM0NBS9evWCTqfDd999Z8tDcHr8RSESg8UB/uabb0Kr1WLr1q04duwYTpw4IW/l5eU2mOLPdu3ahYyMDGRlZeH48eMYOXIk4uLicP369Xb7Hz16FFOnTkVqaipOnDiBxMREJCYm4vTp03KfVatWYd26ddi4cSNKSkrQu3dvxMXF4c6dOzY9FiKinlJIkiRZskNISAhWrlyJmTNn2mhKHdNqtXjooYewYcMGAD+ux4eHh2P+/PlYunRpm/5TpkxBQ0MDDhw4ILc9/PDDiI6OxsaNGyFJEsLCwrBo0SK89NJLAACDwYCQkBDk5+cjKSnJrHkZjUao1WpoJ2bD08vHCkfqOKy+yZm5y1vqm5vuoGT/MhgMBvj5+XXYz+IKXKlU4pFHHunR5Lrj7t27KCsrg06nM5mLTqdDcXFxu/sUFxeb9AeAuLg4uX9VVRX0er1JH7VaDa1W2+GYANDY2MgvsiByABYYpiwO8BdffBG5ubm2mEunbt68iZaWljbftRkSEgK9Xt/uPnq9vtP+P/3XkjEBICcnB2q1Wt7Cw8MtPh4iop6y+DLCl156CQkJCRg0aBCioqLavJFn9+7dVpucs8rMzDQ5YWs0Gl0ixFndEInF4gBfsGABDh06hHHjxuFXv/qV3T7MKigoCB4eHqipqTFpr6mpgUajaXcfjUbTaf+f/ltTU4PQ0FCTPtHR0R3ORaVSQaVSdecwnBbDm0g8Fgf4tm3b8Pe//x0JCQm2mE+HvL29ERMTg6KiIiQmJgL48SRmUVER5s2b1+4+sbGxKCoqwsKFC+W2wsJCxMbGAgAGDBgAjUaDoqIiObCNRiNKSkowd+5cWx4OEVGPWbwGHhgYiEGDBtliLl3KyMjApk2bsG3bNlRUVGDu3LloaGiQP5s8JSUFmZmZcv8//elPKCgowH/8x3/g7NmzePXVV3Hs2DE58BUKBRYuXIjXX38d+/btw6lTp5CSkoKwsDD5RYKInAv/WvyZxRX4q6++iqysLGzduhW+vr62mFOHpkyZghs3bmD58uXQ6/WIjo5GQUGBfBLy0qVLUCp/fk0aO3Ys3n//fbzyyiv493//dzzwwAPYu3cvhg0bJvdZvHgxGhoakJ6ejrq6Ojz66KMoKCiAj4/YlwNagr8QRGKy+DrwUaNG4fz585AkCREREW1OYh4/ftyqExSByNeBM7xJRK5+Pbi514FbXIFzacF1MLyJxGZxgGdlZdliHkREZuMXHf+In0ZIRCQoswI8MDAQN2/eNHvQ+++/HxcvXuz2pMj2uHxCJD6zllDq6urwj3/8A2q12qxBb926xW/mcWIMbyLXYPYa+IwZM2w5DyIii3Ad3MwA59eoERE5H57EdDNcPiFyHQxwIiJBMcDdCKtvcjXu/pw2O8CvXbtmy3kQEZGFzA7wBx98EO+//74t50I25O6VCrkud35umx3gb7zxBmbPno2nn34atbW1tpwTERGZwewAf/7553Hy5EncunULUVFR2L9/vy3nRVbkzhUKkSuz6MOsBgwYgIMHD2LDhg34/e9/j8jISHh6mg7hjh8nS0TkCBZ/GuHFixexe/duBAQEYNKkSW0CnJwLq29yB+76rkyL0nfTpk1YtGgRdDodzpw5g+DgYFvNi4iIumB2gMfHx6O0tBQbNmxASkqKLedEVsLqm8i1mR3gLS0tOHnyJO677z5bzoeIiMxk9lUohYWFDG8iclru+Bcn30rvotzxyUzkbhjgLojhTeQeGOAuhuFN5D4Y4ETkMtytgGGAExEJigHuQtyt+iBydwxwF8HwJnI/DHAiIkExwInIpbjTX6MMcBfgTk9YIvoZA1xwDG8i9yVMgNfW1iI5ORl+fn7w9/dHamoq6uvrO+0/f/58DBkyBL169cL999+PBQsWwGAwmPRTKBRttp07d9r6cIiIekyYAE9OTsaZM2dQWFiIAwcO4MiRI0hPT++w/7Vr13Dt2jWsXr0ap0+fRn5+PgoKCpCamtqm79atW1FdXS1viYmJNjwSIrI1d/nLVIiv06moqEBBQQG+/vprjB49GgCwfv16PPXUU1i9ejXCwsLa7DNs2DD8/e9/l38eNGgQ3njjDUybNg3Nzc0m3yTk7+8PjUZj+wOxMnd5khJR+4SowIuLi+Hv7y+HNwDodDoolUqUlJSYPY7BYICfn1+br4F74YUXEBQUhDFjxiAvLw+SJHU6TmNjI4xGo8lGRGRvQlTger0e/fr1M2nz9PREYGAg9Hq9WWPcvHkT2dnZbZZdVqxYgd/+9rfw9fXFZ599hueffx719fVYsGBBh2Pl5OTgtddes/xArIjVNxE5tAJfunRpuycRf7mdPXu2x/djNBqRkJCAqKgovPrqqya3LVu2DI888ghGjRqFJUuWYPHixXj77bc7HS8zMxMGg0HeLl++3OM5EpF1uUOR49AKfNGiRZg5c2anfQYOHAiNRoPr16+btDc3N6O2trbLtevvv/8e8fHx6Nu3L/bs2QMvL69O+2u1WmRnZ6OxsREqlardPiqVqsPb7MEdnphE1DWHBnhwcLBZ32wfGxuLuro6lJWVISYmBgBw8OBBtLa2QqvVdrif0WhEXFwcVCoV9u3bBx8fny7vq7y8HAEBAQ4NaCIicwixBh4ZGYn4+HikpaVh48aNaGpqwrx585CUlCRfgXL16lWMHz8e27dvx5gxY2A0GjFhwgTcvn0b7733nsnJxuDgYHh4eGD//v2oqanBww8/DB8fHxQWFuLNN9/ESy+95MjD7RSrbyL6iRABDgA7duzAvHnzMH78eCiVSkyePBnr1q2Tb29qasK5c+dw+/ZtAMDx48flK1QGDx5sMlZVVRUiIiLg5eWF3NxcvPjii5AkCYMHD8aaNWuQlpZmvwMjIpsZtKQC51dGOnoaNqOQurpmjrpkNBqhVquhnZgNT6+ul2m6i9U3keVEDPDmpjso2b9MvvS5I0JcB05ERG0xwAXB6puI7sUAJyKX5srFDwOciEhQDHABuHIFQUTdxwB3cgxvIuoIA9yJMbyJqDMMcCJyea5aDDHAiYgExQB3Uq5aMRCR9TDAnRDDm4jMwQAnIrfgioURA5yISFAMcCfjilUCEdkGA9yJMLyJyBIMcCJyG65WJDHAiYgExQB3Eq5WGRCR7THAiYgExQB3Aqy+iezHlX7fGOBERIJigDuYK1UDRKJwld87BjgRkaAY4A7kKlUAETkGA5yISFAMcAdh9U3kWK7wO8gAJyISFAPcAVzhlZ+IHI8BTkQkKAY4Ebkt0f8aZoDbmehPGCJyHsIEeG1tLZKTk+Hn5wd/f3+kpqaivr6+032eeOIJKBQKk23OnDkmfS5duoSEhAT4+vqiX79+ePnll9Hc3GyTY2B4E5E1eTp6AuZKTk5GdXU1CgsL0dTUhFmzZiE9PR3vv/9+p/ulpaVhxYoV8s++vr7yv1taWpCQkACNRoOjR4+iuroaKSkp8PLywptvvmnV+TO8icjahKjAKyoqUFBQgM2bN0Or1eLRRx/F+vXrsXPnTly7dq3TfX19faHRaOTNz89Pvu2zzz7Dt99+i/feew/R0dF48sknkZ2djdzcXNy9e9fWh0VETkDk4kqIAC8uLoa/vz9Gjx4tt+l0OiiVSpSUlHS6744dOxAUFIRhw4YhMzMTt2/fNhl3+PDhCAkJkdvi4uJgNBpx5syZDsdsbGyE0Wg02YiI7E2IJRS9Xo9+/fqZtHl6eiIwMBB6vb7D/Z599ln0798fYWFhOHnyJJYsWYJz585h9+7d8ri/DG8A8s+djZuTk4PXXnvNomM4vzJS6Fd6InI+Dq3Aly5d2uYk473b2bNnuz1+eno64uLiMHz4cCQnJ2P79u3Ys2cPzp8/36N5Z2ZmwmAwyNvly5fN2u/8ysge3S8R0S85tAJftGgRZs6c2WmfgQMHQqPR4Pr16ybtzc3NqK2thUajMfv+tFotAKCyshKDBg2CRqNBaWmpSZ+amhoA6HRclUoFlUpl9v3+0r0hzqqciLrLoQEeHByM4ODgLvvFxsairq4OZWVliImJAQAcPHgQra2tciibo7y8HAAQGhoqj/vGG2/g+vXr8hJNYWEh/Pz8EBUVZeHRdA8DncjxBi2pEPIvZCFOYkZGRiI+Ph5paWkoLS3Fl19+iXnz5iEpKQlhYWEAgKtXr2Lo0KFyRX3+/HlkZ2ejrKwMFy5cwL59+5CSkoLHHnsMI0aMAABMmDABUVFRmD59Or755ht8+umneOWVV/DCCy90u8LuKRGfRETkGEKcxAR+vJpk3rx5GD9+PJRKJSZPnox169bJtzc1NeHcuXPyVSbe3t747//+b6xduxYNDQ0IDw/H5MmT8corr8j7eHh44MCBA5g7dy5iY2PRu3dvzJgxw+S6cUdgVU5E5lBIkiQ5ehKiMxqNUKvV0E7MhqeXj83vj4FOZH3O9Ndvc9MdlOxfBoPBYPLelXsJsYRCRGRrIhZGwiyh0M+4xEJEAAPcJTDQidwTl1BckDOt5RGR7bACd1GsyoksJ9r14AxwN8FAJ3I9XEJxUyJVGUTUPlbgboxVOZHYGOAkY6ATibUOzgCnDjHQiZwb18DJbKJUJUTughU4WYRVOZHzYIBTjzDQiRyHSyhERPcQpRBhBU5WxYqcyH4Y4GRTDHQi2+ESCtkVr2Qhsh5W4GR3rMpJBCK8oYcBTg7HQCfqHi6hkNNx9qqHyFmwAienxKqcqGsMcBICA50cwdnXwRngJCQGOhHXwMlFOHOVRGQrrMDJZbAqJ3fDACeXxUAna3DmdXAuoZDbcNZfQqLuYgVOboVVOXWHs1bhDHByawx0EhmXUIiIBMUKnOgXWJGTSISpwGtra5GcnAw/Pz/4+/sjNTUV9fX1Hfa/cOECFApFu9uHH34o92vv9p07d9rjkEgA51dGmmzkvpzxxVyYAE9OTsaZM2dQWFiIAwcO4MiRI0hPT++wf3h4OKqrq0221157DX369MGTTz5p0nfr1q0m/RITE218NCQqhjg5EyGWUCoqKlBQUICvv/4ao0ePBgCsX78eTz31FFavXo2wsLA2+3h4eECj0Zi07dmzB8888wz69Olj0u7v79+mL1FHuMxCzkKICry4uBj+/v5yeAOATqeDUqlESUmJWWOUlZWhvLwcqampbW574YUXEBQUhDFjxiAvLw+SJHU6VmNjI4xGo8lG7ovLLOQoQgS4Xq9Hv379TNo8PT0RGBgIvV5v1hhbtmxBZGQkxo4da9K+YsUKfPDBBygsLMTkyZPx/PPPY/369Z2OlZOTA7VaLW/h4eGWHRC5NIa463K2v7YcuoSydOlSrFy5stM+FRU9f8B++OEHvP/++1i2bFmb237ZNmrUKDQ0NODtt9/GggULOhwvMzMTGRkZ8s9Go5EhTia4zEL24NAAX7RoEWbOnNlpn4EDB0Kj0eD69esm7c3NzaitrTVr7fpvf/sbbt++jZSUlC77arVaZGdno7GxESqVqt0+KpWqw9uI2sNAJ1twaIAHBwcjODi4y36xsbGoq6tDWVkZYmJiAAAHDx5Ea2srtFptl/tv2bIFv/vd78y6r/LycgQEBDCgyaYY6GQNQqyBR0ZGIj4+HmlpaSgtLcWXX36JefPmISkpSb4C5erVqxg6dChKS0tN9q2srMSRI0fw3HPPtRl3//792Lx5M06fPo3Kykr85S9/wZtvvon58+fb5biIfsJ1c3E404utEJcRAsCOHTswb948jB8/HkqlEpMnT8a6devk25uamnDu3Dncvn3bZL+8vDzcd999mDBhQpsxvby8kJubixdffBGSJGHw4MFYs2YN0tLSbH48RPdiVU6WUkhdXTNHXTIajVCr1dBOzIanl4+jp0MuioHuPGz9F1Nz0x2U7F8Gg8EAPz+/DvsJsYRCRERtCbOEQuTuuMRC92KAEwmKge44zvIFD1xCIXIRzhAoZF+swIlcCKty98IAJ3JhDHTXxiUUIjfCZRbrcYYXQ1bgRG6GVbnrYIATuTkGurgY4ERkgoEuDq6BE1GnuG7eMUe/uLECJ6IusSp3TgxwIrIYA905cAmFiEhQrMCJqMfaWyd3l6rckZ+LwgAnIpvgMovtcQmFiEhQrMCJyC5YkVsfA5yIHMKVAt1R6+BcQiEip8A3DFmOFTgROQ1XqsrtgQFORE6Lgd45LqEQkTC4zGKKFTgRCcVZq3JHnMhkgBOR0Jw10O2BAU5ELsWdAp1r4ETk0lx53ZwVOBG5PHtV5fZeB2eAE5HbcZVlFi6hEBEJihU4Ebk9UStyYSrwN954A2PHjoWvry/8/f3N2keSJCxfvhyhoaHo1asXdDodvvvuO5M+tbW1SE5Ohp+fH/z9/ZGamor6+nobHAERieL8ykiTzRL2DH9hAvzu3bt4+umnMXfuXLP3WbVqFdatW4eNGzeipKQEvXv3RlxcHO7cuSP3SU5OxpkzZ1BYWIgDBw7gyJEjSE9Pt8UhEJGgnPVKFoUkSZKjJ2GJ/Px8LFy4EHV1dZ32kyQJYWFhWLRoEV566SUAgMFgQEhICPLz85GUlISKigpERUXh66+/xujRowEABQUFeOqpp3DlyhWEhYWZNSej0Qi1Wg3txGx4evn06PiISAxdVdo9Cf3mpjso2b8MBoMBfn5+HfYTpgK3VFVVFfR6PXQ6ndymVquh1WpRXFwMACguLoa/v78c3gCg0+mgVCpRUlLS4diNjY0wGo0mGxG5l86WWexVsbtsgOv1egBASEiISXtISIh8m16vR79+/Uxu9/T0RGBgoNynPTk5OVCr1fIWHh5u5dkTkWjc7gsdli5dCoVC0el29uxZR06xXZmZmTAYDPJ2+fJlR0+JiJxAd0569oRDLyNctGgRZs6c2WmfgQMHdmtsjUYDAKipqUFoaKjcXlNTg+joaLnP9evXTfZrbm5GbW2tvH97VCoVVCpVt+ZFRGQtDg3w4OBgBAcH22TsAQMGQKPRoKioSA5so9GIkpIS+UqW2NhY1NXVoaysDDExMQCAgwcPorW1FVqt1ibzIiKyFmHWwC9duoTy8nJcunQJLS0tKC8vR3l5uck120OHDsWePXsAAAqFAgsXLsTrr7+Offv24dSpU0hJSUFYWBgSExMBAJGRkYiPj0daWhpKS0vx5ZdfYt68eUhKSjL7ChQiIkcR5p2Yy5cvx7Zt2+SfR40aBQA4dOgQnnjiCQDAuXPnYDAY5D6LFy9GQ0MD0tPTUVdXh0cffRQFBQXw8fn5Ur8dO3Zg3rx5GD9+PJRKJSZPnox169bZ56CIiHpAuOvAnRGvAycia3L768CJiFwdA5yISFAMcCIiQTHAiYgExQAnIhIUA5yISFAMcCIiQTHAiYgExQAnIhKUMG+ld2Y/vZm1uelOFz2JiLr2U5Z09UZ5vpXeCq5cucIvdSAiq7t8+TLuu+++Dm9ngFtBa2srrl27hr59+0KhULS53Wg0Ijw8HJcvX+70cw3IPHw8rYuPp3VZ4/GUJAnff/89wsLCoFR2vNLNJRQrUCqVnb5K/sTPz4+/IFbEx9O6+HhaV08fT7Va3WUfnsQkIhIUA5yISFAMcDtQqVTIysri92haCR9P6+LjaV32fDx5EpOISFCswImIBMUAJyISFAOciEhQDHAiIkExwG2ktrYWycnJ8PPzg7+/P1JTU1FfX9/pPk888QQUCoXJNmfOHDvN2Lnk5uYiIiICPj4+0Gq1KC0t7bT/hx9+iKFDh8LHxwfDhw/HJ598YqeZisGSxzM/P7/N89DHx8eOs3VuR44cwcSJExEWFgaFQoG9e/d2uc/hw4fxm9/8BiqVCoMHD0Z+fr5V5sIAt5Hk5GScOXMGhYWFOHDgAI4cOYL09PQu90tLS0N1dbW8rVq1yg6zdS67du1CRkYGsrKycPz4cYwcORJxcXG4fv16u/2PHj2KqVOnIjU1FSdOnEBiYiISExNx+vRpO8/cOVn6eAI/vovwl8/Dixcv2nHGzq2hoQEjR45Ebm6uWf2rqqqQkJCAcePGoby8HAsXLsRzzz2HTz/9tOeTkcjqvv32WwmA9PXXX8tt//jHPySFQiFdvXq1w/0ef/xx6U9/+pMdZujcxowZI73wwgvyzy0tLVJYWJiUk5PTbv9nnnlGSkhIMGnTarXS7NmzbTpPUVj6eG7dulVSq9V2mp3YAEh79uzptM/ixYulBx980KRtypQpUlxcXI/vnxW4DRQXF8Pf3x+jR4+W23Q6HZRKJUpKSjrdd8eOHQgKCsKwYcOQmZmJ27dv23q6TuXu3bsoKyuDTqeT25RKJXQ6HYqLi9vdp7i42KQ/AMTFxXXY35105/EEgPr6evTv3x/h4eGYNGkSzpw5Y4/puiRbPj/5YVY2oNfr0a9fP5M2T09PBAYGQq/Xd7jfs88+i/79+yMsLAwnT57EkiVLcO7cOezevdvWU3YaN2/eREtLC0JCQkzaQ0JCcPbs2Xb30ev17fbv7LF2F915PIcMGYK8vDyMGDECBoMBq1evxtixY3HmzBmzPrSNTHX0/DQajfjhhx/Qq1evbo/NALfA0qVLsXLlyk77VFRUdHv8X66RDx8+HKGhoRg/fjzOnz+PQYMGdXtcIkvExsYiNjZW/nns2LGIjIzEu+++i+zsbAfOjO7FALfAokWLMHPmzE77DBw4EBqNps0JoubmZtTW1kKj0Zh9f1qtFgBQWVnpNgEeFBQEDw8P1NTUmLTX1NR0+NhpNBqL+ruT7jye9/Ly8sKoUaNQWVlpiym6vI6en35+fj2qvgFehWKR4OBgDB06tNPN29sbsbGxqKurQ1lZmbzvwYMH0draKoeyOcrLywEAoaGh1j4Up+Xt7Y2YmBgUFRXJba2trSgqKjKpCn8pNjbWpD8AFBYWdtjfnXTn8bxXS0sLTp065VbPQ2uy6fOzx6dBqV3x8fHSqFGjpJKSEumLL76QHnjgAWnq1Kny7VeuXJGGDBkilZSUSJIkSZWVldKKFSukY8eOSVVVVdJHH30kDRw4UHrsscccdQgOs3PnTkmlUkn5+fnSt99+K6Wnp0v+/v6SXq+XJEmSpk+fLi1dulTu/+WXX0qenp7S6tWrpYqKCikrK0vy8vKSTp065ahDcCqWPp6vvfaa9Omnn0rnz5+XysrKpKSkJMnHx0c6c+aMow7BqXz//ffSiRMnpBMnTkgApDVr1kgnTpyQLl68KEmSJC1dulSaPn263P9///d/JV9fX+nll1+WKioqpNzcXMnDw0MqKCjo8VwY4DZy69YtaerUqVKfPn0kPz8/adasWdL3338v315VVSUBkA4dOiRJkiRdunRJeuyxx6TAwEBJpVJJgwcPll5++WXJYDA46Agca/369dL9998veXt7S2PGjJG++uor+bbHH39cmjFjhkn/Dz74QPr1r38teXt7Sw8++KD08ccf23nGzs2Sx3PhwoVy35CQEOmpp56Sjh8/7oBZO6dDhw5JANpsPz2GM2bMkB5//PE2+0RHR0ve3t7SwIEDpa1bt1plLvw4WSIiQXENnIhIUAxwIiJBMcCJiATFACciEhQDnIhIUAxwIiJBMcCJiATFACciEhQDnKiHIiIi5K8eq6urs/v9Hz58WL7/xMREu98/OQ4DnAg/fmDT2LFj8fvf/96k3WAwIDw8HH/+85873X/FihWorq6GWq225TTbNXbsWFRXV+OZZ56x+32TYzHAiQB4eHggPz8fBQUF2LFjh9w+f/58BAYGIisrq9P9+/btC41GA4VCYeuptuHt7Q2NRtPjjyYl8TDAif7p17/+Nd566y3Mnz8f1dXV+Oijj7Bz505s374d3t7eFo2Vn58Pf39/HDhwAEOGDIGvry/+8Ic/4Pbt29i2bRsiIiIQEBCABQsWoKWlRd4vIiICr7/+OlJSUtCnTx/0798f+/btw40bNzBp0iT06dMHI0aMwLFjx6x9+CQgBjjRL8yfPx8jR47E9OnTkZ6ejuXLl2PkyJHdGuv27dtYt24ddu7ciYKCAhw+fBj/9m//hk8++QSffPIJ/uu//gvvvvsu/va3v5ns98477+CRRx7BiRMnkJCQgOnTpyMlJQXTpk3D8ePHMWjQIKSkpICfQ0f8OFmie1RUVEgApOHDh0tNTU1d9u/fv7/0zjvvmLRt3bpVAiBVVlbKbbNnz5Z8fX1NPlY4Li5Omj17tslY06ZNk3+urq6WAEjLli2T24qLiyUAUnV1tcl9zpgxQ5o0aZK5h0kugBU40T3y8vLg6+uLqqoqXLlypdvj+Pr6mnwVXkhICCIiItCnTx+Ttnu/fm/EiBEmtwM/fkfqvW337kfuhwFO9AtHjx7FO++8gwMHDmDMmDFITU3t9lKFl5eXyc8KhaLdttbW1g73++mkaHtt9+5H7ocBTvRPt2/fxsyZMzF37lyMGzcOW7ZsQWlpKTZu3OjoqRG1iwFO9E+ZmZmQJAlvvfUWgB+vCFm9ejUWL16MCxcuOHZyRO1ggBMB+Pzzz5Gbm4utW7fC19dXbp89ezbGjh3bo6UUIlvhd2IS9VBERAQWLlyIhQsXOnQeM2fORF1dHfbu3evQeZD9sAInsoIlS5agT58+MBgMdr/v//mf/0GfPn1M3kFK7oEVOFEPXbx4EU1NTQCAgQMHQqm0b130ww8/4OrVqwCAPn36QKPR2PX+yXEY4EREguISChGRoBjgRESCYoATEQmKAU5EJCgGOBGRoBjgRESCYoATEQmKAU5EJKj/B9nufxXfDpJYAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "lens = Optic()\n", + "\n", + "ap = PolygonAperture(x=[-0.5, 0, 0.5, 1], y=[-0.5, 0.5, 1, -1])\n", + "\n", + "lens.add_surface(index=0, radius=np.inf, thickness=np.inf)\n", + "lens.add_surface(index=1, radius=22.01359, thickness=3.25896, material=\"SK16\")\n", + "lens.add_surface(index=2, radius=-435.76044, thickness=6.00755)\n", + "lens.add_surface(index=3, radius=-22.21328, thickness=0.99997, material=(\"F2\", \"schott\"), aperture=ap)\n", + "lens.add_surface(index=4, radius=20.29192, thickness=4.75041, is_stop=True)\n", + "lens.add_surface(index=5, radius=79.68360, thickness=2.95208, material=\"SK16\")\n", + "lens.add_surface(index=6, radius=-18.39533, thickness=42.20778)\n", + "lens.add_surface(index=7)\n", + "\n", + "lens.set_aperture(aperture_type=\"EPD\", value=3)\n", + "\n", + "lens.set_field_type(field_type=\"angle\")\n", + "lens.add_field(y=0)\n", + "lens.add_field(y=4)\n", + "\n", + "lens.add_wavelength(value=0.48)\n", + "lens.add_wavelength(value=0.55, is_primary=True)\n", + "\n", + "lens.update_paraxial()\n", + "\n", + "solver = QuickFocusSolve(lens)\n", + "solver.apply()\n", + "\n", + "lens.info()\n", + "lens.draw()\n", + "ap.view()" + ] + }, + { + "cell_type": "markdown", + "id": "asm", + "metadata": {}, + "source": [ + "## 3. Compute the propagated field\n", + "\n", + "We compute the complex field at multiple planes `z`.\n", + "\n", + "Output convention used in this development notebook:\n", + "- If the propagator returns `(W, N, N)`, we treat it as a single field and store as `(F=1, W, N, N)`.\n", + "- If the propagator returns `(F, W, N, N)`, we store it directly.\n", + "\n", + "The final stacked output is shaped as:\n", + "- `field_stack`: `(S, F, W, N, N)` where `S = num_slices`." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "helpers", + "metadata": {}, + "outputs": [], + "source": [ + "def _is_torch(x):\n", + " return x.__class__.__module__.startswith(\"torch\")\n", + "\n", + "\n", + "def _to_numpy(x):\n", + " if _is_torch(x):\n", + " x = x.detach()\n", + " if getattr(x, \"is_cuda\", False):\n", + " x = x.cpu()\n", + " x = x.numpy()\n", + " return np.asarray(x)\n", + "\n", + "\n", + "def compute_field_stack(\n", + " lens,\n", + " distance,\n", + " wave_resolution,\n", + " pitch,\n", + " field=(0.0, 0.0),\n", + " wavelengths=\"all\",\n", + " num_slices=50,\n", + " beam_waist=\"auto\",\n", + "):\n", + " import time\n", + " import numpy as np\n", + "\n", + " z = np.linspace(0.0, float(distance) * 1.1, num_slices)\n", + "\n", + " t0 = time.perf_counter()\n", + " U0 = lens.wave_propagator.compute_field(\n", + " z_target=float(z[0]),\n", + " num_points=wave_resolution,\n", + " dx=pitch,\n", + " field=field,\n", + " wavelengths=wavelengths,\n", + " beam_waist=beam_waist,\n", + " )\n", + " dt = time.perf_counter() - t0\n", + " print(f\"slice 1/{num_slices} | z={z[0]:.6g} | compute_field: {dt:.4f}s\")\n", + "\n", + " nd = getattr(U0, \"ndim\", None)\n", + " is_torch = U0.__class__.__module__.startswith(\"torch\")\n", + "\n", + " if nd == 3:\n", + " F, W, N = 1, int(U0.shape[0]), int(U0.shape[1])\n", + " elif nd == 4:\n", + " F, W, N = int(U0.shape[0]), int(U0.shape[1]), int(U0.shape[2])\n", + " else:\n", + " raise ValueError(f\"Expected (W,N,N) or (F,W,N,N), got {getattr(U0,'shape',None)}\")\n", + "\n", + " if is_torch:\n", + " import torch\n", + "\n", + " stack = torch.empty((num_slices, F, W, N, N), dtype=U0.dtype, device=U0.device)\n", + " if nd == 3:\n", + " stack[0, 0].copy_(U0)\n", + " else:\n", + " stack[0].copy_(U0)\n", + "\n", + " for i in range(1, num_slices):\n", + " zi = float(z[i])\n", + " t0 = time.perf_counter()\n", + " U = lens.wave_propagator.compute_field(\n", + " z_target=zi,\n", + " num_points=wave_resolution,\n", + " dx=pitch,\n", + " field=field,\n", + " wavelengths=wavelengths,\n", + " beam_waist=beam_waist,\n", + " )\n", + " dt = time.perf_counter() - t0\n", + " print(f\"slice {i+1}/{num_slices} | z={zi:.6g} | compute_field: {dt:.4f}s\")\n", + "\n", + " if nd == 3:\n", + " stack[i, 0].copy_(U)\n", + " else:\n", + " stack[i].copy_(U)\n", + "\n", + " return {\"z\": z, \"field\": stack}\n", + "\n", + " stack = np.empty((num_slices, F, W, N, N), dtype=U0.dtype)\n", + " stack[0, 0] = U0 if nd == 3 else stack.__setitem__((0,), U0)\n", + "\n", + " if nd == 3:\n", + " stack[0, 0] = U0\n", + " else:\n", + " stack[0] = U0\n", + "\n", + " for i in range(1, num_slices):\n", + " zi = float(z[i])\n", + " t0 = time.perf_counter()\n", + " U = lens.wave_propagator.compute_field(\n", + " z_target=zi,\n", + " num_points=wave_resolution,\n", + " dx=pitch,\n", + " field=field,\n", + " wavelengths=wavelengths,\n", + " beam_waist=beam_waist,\n", + " )\n", + " dt = time.perf_counter() - t0\n", + " print(f\"slice {i+1}/{num_slices} | z={zi:.6g} | compute_field: {dt:.4f}s\")\n", + "\n", + " if nd == 3:\n", + " stack[i, 0] = U\n", + " else:\n", + " stack[i] = U\n", + "\n", + " return {\"z\": z, \"field\": stack}" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "run-stack", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "wave_resolution: 600\n", + "slice 1/30 | z=0 | compute_field: 0.0319s\n", + "slice 2/30 | z=3.7931 | compute_field: 0.3158s\n", + "slice 3/30 | z=7.58621 | compute_field: 0.3276s\n", + "slice 4/30 | z=11.3793 | compute_field: 0.6429s\n", + "slice 5/30 | z=15.1724 | compute_field: 0.7725s\n", + "slice 6/30 | z=18.9655 | compute_field: 0.9692s\n", + "slice 7/30 | z=22.7586 | compute_field: 0.9294s\n", + "slice 8/30 | z=26.5517 | compute_field: 0.9284s\n", + "slice 9/30 | z=30.3448 | compute_field: 0.9453s\n", + "slice 10/30 | z=34.1379 | compute_field: 0.9404s\n", + "slice 11/30 | z=37.931 | compute_field: 0.9349s\n", + "slice 12/30 | z=41.7241 | compute_field: 0.9331s\n", + "slice 13/30 | z=45.5172 | compute_field: 0.9466s\n", + "slice 14/30 | z=49.3103 | compute_field: 0.9376s\n", + "slice 15/30 | z=53.1034 | compute_field: 0.9239s\n", + "slice 16/30 | z=56.8966 | compute_field: 0.9456s\n", + "slice 17/30 | z=60.6897 | compute_field: 1.1196s\n", + "slice 18/30 | z=64.4828 | compute_field: 1.1215s\n", + "slice 19/30 | z=68.2759 | compute_field: 1.1048s\n", + "slice 20/30 | z=72.069 | compute_field: 1.1066s\n", + "slice 21/30 | z=75.8621 | compute_field: 1.0994s\n", + "slice 22/30 | z=79.6552 | compute_field: 1.1782s\n", + "slice 23/30 | z=83.4483 | compute_field: 1.0989s\n", + "slice 24/30 | z=87.2414 | compute_field: 1.0968s\n", + "slice 25/30 | z=91.0345 | compute_field: 1.0945s\n", + "slice 26/30 | z=94.8276 | compute_field: 1.0739s\n", + "slice 27/30 | z=98.6207 | compute_field: 1.1010s\n", + "slice 28/30 | z=102.414 | compute_field: 1.0952s\n", + "slice 29/30 | z=106.207 | compute_field: 1.1070s\n", + "slice 30/30 | z=110 | compute_field: 1.0909s\n" + ] + } + ], + "source": [ + "pitch = 5e-3\n", + "wave_resolution = int(lens.aperture.value / pitch)\n", + "print(\"wave_resolution:\", wave_resolution)\n", + "\n", + "stack = compute_field_stack(\n", + " lens=lens,\n", + " distance=100,\n", + " wave_resolution=wave_resolution,\n", + " pitch=pitch,\n", + " field=[(0, 0), (0.1, 0.2)],\n", + " wavelengths=\"all\",\n", + " num_slices=30,\n", + " beam_waist=\"auto\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "viz", + "metadata": {}, + "source": [ + "## 4. Visualizing Intensity Along Propagation\n", + "\n", + "The propagated complex field can be inspected as normalized intensity at different axial positions `z`.\n", + "\n", + "The interactive viewer allows you to:\n", + "\n", + "- Scroll through propagation slices.\n", + "- Select field index and wavelength.\n", + "- Control zoom around the optical axis.\n", + "- Choose the normalization strategy.\n", + "\n", + "Normalization modes:\n", + "\n", + "- `norm=\"global\"` uses a single maximum over all slices for the selected field and wavelength. This preserves relative power variations along propagation.\n", + "- `norm=\"per_frame\"` normalizes each slice independently. This enhances local structure but removes axial power comparison.\n", + "\n", + "The display is cropped around the center of the grid to emphasize the main lobe and focal evolution." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "plot-slices", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "c01636cc46fc4987b2000a8e28309dd0", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(Play(value=0, interval=80, max=29), IntSlider(value=0, description='slice', max=29)))" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "b3619738b55c4fa9bfe422bfc2ddfb38", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "HBox(children=(IntSlider(value=0, description='field', max=1), IntSlider(value=0, description='wl', max=1), In…" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "2853ad323e4e4013b30aca0cd0295e8f", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Output()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def visualize_field_stack(stack, pitch, cmap=\"inferno\"):\n", + " z = stack[\"z\"]\n", + " U = stack[\"field\"]\n", + "\n", + " is_torch = U.__class__.__module__.startswith(\"torch\")\n", + " S, F, W, H, W2 = U.shape\n", + "\n", + " def to_numpy(x):\n", + " if x.__class__.__module__.startswith(\"torch\"):\n", + " x = x.detach()\n", + " if getattr(x, \"is_cuda\", False):\n", + " x = x.cpu()\n", + " x = x.numpy()\n", + " return np.asarray(x)\n", + "\n", + " def intensity(slice_idx, field_idx, wl_idx, norm):\n", + " f = U[slice_idx, field_idx, wl_idx]\n", + " if is_torch:\n", + " I = f.abs() ** 2\n", + " if norm == \"global\":\n", + " denom = (U[:, field_idx, wl_idx].abs() ** 2).max().clamp_min(1e-12)\n", + " I = I / denom\n", + " else:\n", + " I = I / I.max().clamp_min(1e-12)\n", + " return to_numpy(I).astype(np.float64, copy=False)\n", + "\n", + " I = (np.abs(f) ** 2).astype(np.float64, copy=False)\n", + " if norm == \"global\":\n", + " I /= float((np.abs(U[:, field_idx, wl_idx]) ** 2).max() + 1e-12)\n", + " else:\n", + " I /= float(I.max() + 1e-12)\n", + " return I\n", + "\n", + " out = widgets.Output()\n", + "\n", + " play = widgets.Play(interval=80, value=0, min=0, max=S - 1, step=1)\n", + " slice_slider = widgets.IntSlider(value=0, min=0, max=S - 1, step=1, description=\"slice\")\n", + " widgets.jslink((play, \"value\"), (slice_slider, \"value\"))\n", + "\n", + " field_slider = widgets.IntSlider(value=0, min=0, max=F - 1, step=1, description=\"field\")\n", + " wl_slider = widgets.IntSlider(value=0, min=0, max=W - 1, step=1, description=\"wl\")\n", + " zoom_slider = widgets.IntSlider(value=0, min=1, max=15, step=1, description=\"zoom\")\n", + " norm = widgets.ToggleButtons(options=[\"global\", \"per_frame\"], value=\"per_frame\", description=\"norm\")\n", + "\n", + " def render(*_):\n", + " s = int(slice_slider.value)\n", + " f = int(field_slider.value)\n", + " w = int(wl_slider.value)\n", + " zoom = int(zoom_slider.value)\n", + " nrm = str(norm.value)\n", + "\n", + " I = intensity(s, f, w, nrm)\n", + "\n", + " cy, cx = H // 2, W2 // 2\n", + " half = max(int(H // (2 * zoom)), 5)\n", + " y0, y1 = max(cy - half, 0), min(cy + half, H)\n", + " x0, x1 = max(cx - half, 0), min(cx + half, W2)\n", + " crop = I[y0:y1, x0:x1]\n", + " extent = half * pitch\n", + "\n", + " with out:\n", + " out.clear_output(wait=True)\n", + " fig, ax = plt.subplots(figsize=(7, 7))\n", + " im = ax.imshow(\n", + " crop,\n", + " cmap=cmap,\n", + " origin=\"upper\",\n", + " extent=[-extent, extent, -extent, extent],\n", + " vmin=0.0,\n", + " vmax=1.0,\n", + " )\n", + " fig.colorbar(im, ax=ax, label=\"intensity (normalized)\")\n", + " ax.set_title(f\"z={z[s]:.6g} | slice={s} | field={f} | wl={w}\")\n", + " ax.set_xlabel(\"x\")\n", + " ax.set_ylabel(\"y\")\n", + " plt.show()\n", + "\n", + " for wdg in (slice_slider, field_slider, wl_slider, zoom_slider, norm):\n", + " wdg.observe(render, names=\"value\")\n", + "\n", + " display(widgets.HBox([play, slice_slider]))\n", + " display(widgets.HBox([field_slider, wl_slider, zoom_slider, norm]))\n", + " display(out)\n", + " render()\n", + "\n", + "visualize_field_stack(stack, pitch)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2ee2e77-386c-4b0e-b04c-ecdc5907fe74", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/optiland/backend/numpy_backend.py b/optiland/backend/numpy_backend.py index a61829ab..c1bffaec 100644 --- a/optiland/backend/numpy_backend.py +++ b/optiland/backend/numpy_backend.py @@ -187,6 +187,18 @@ def fftconvolve( return _fftconvolve(a, b, mode=mode) +def fftfreq(n: int, d: float = 1.0) -> NDArray: + return np.fft.fftfreq(n, d=d) + + +def fft2(x: ArrayLike) -> NDArray: + return np.fft.fft2(array(x)) + + +def ifft2(x: ArrayLike) -> NDArray: + return np.fft.ifft2(array(x)) + + def grid_sample( input, grid, mode="bilinear", padding_mode="zeros", align_corners=False ): @@ -265,3 +277,32 @@ def polyval(coeffs: ArrayLike, x: ArrayLike) -> NDArray: x: A number or an array of numbers at which to evaluate. """ return np.polyval(coeffs, x) + + +def meshgrid(*arrays: ArrayLike, indexing: str = "xy"): + return np.meshgrid(*arrays, indexing=indexing) + + +def pad( + x: NDArray, + pad_width: tuple[tuple[int, int], tuple[int, int]], + mode: Literal["constant"] = "constant", + constant_values: float = 0.0, +) -> NDArray: + (pt, pb), (pl, pr) = pad_width + pads = [(0, 0)] * x.ndim + pads[-2] = (pt, pb) + pads[-1] = (pl, pr) + return np.pad(x, pads, mode=mode, constant_values=constant_values) + + +def clamp(x: ArrayLike, min: float | None = None, max: float | None = None) -> NDArray: + return np.clip(array(x), a_min=min, a_max=max) + + +def broadcast_to(x: ArrayLike, shape) -> NDArray: + return np.broadcast_to(array(x), shape) + + +def zeros(shape) -> NDArray: + return np.zeros(shape, dtype=complex) diff --git a/optiland/backend/torch_backend.py b/optiland/backend/torch_backend.py index ed7d71e4..ec247ac6 100644 --- a/optiland/backend/torch_backend.py +++ b/optiland/backend/torch_backend.py @@ -309,6 +309,10 @@ def cast(x: ArrayLike) -> Tensor: return x.to(device=get_device(), dtype=get_precision()) +def clamp(x: Tensor, min: float | None = None, max: float | None = None) -> Tensor: + return torch.clamp(x, min=min, max=max) + + def copy(x: Tensor) -> Tensor: return x.clone() @@ -357,8 +361,8 @@ def flip(x: Tensor) -> Tensor: return torch.flip(x, dims=(0,)) -def meshgrid(*arrays: Tensor) -> tuple[Tensor, ...]: - return torch.meshgrid(*arrays, indexing="xy") +def meshgrid(*arrays: Tensor, indexing: str = "xy") -> tuple[Tensor, ...]: + return torch.meshgrid(*arrays, indexing=indexing) def roll( @@ -865,6 +869,18 @@ def eye(n: int) -> Tensor: # -------------------------- # Signal Processing # -------------------------- +def fftfreq(n: int, d: float = 1.0) -> Tensor: + return torch.fft.fftfreq(n, d=d, device=get_device(), dtype=get_precision()) + + +def fft2(x: Tensor) -> Tensor: + return torch.fft.fft2(x) + + +def ifft2(x: Tensor) -> Tensor: + return torch.fft.ifft2(x) + + def fftconvolve( in1: Tensor, in2: Tensor, mode: Literal["full", "valid", "same"] = "full" ) -> Tensor: @@ -997,6 +1013,7 @@ def grid_sample( "load", # Utilities "cast", + "clamp", "copy", "is_array_like", "size", @@ -1063,4 +1080,8 @@ def grid_sample( # Simulation "fftconvolve", "grid_sample", + # FFT + "fftfreq", + "fft2", + "ifft2", ] diff --git a/optiland/optic/optic.py b/optiland/optic/optic.py index 31713e48..1e441807 100644 --- a/optiland/optic/optic.py +++ b/optiland/optic/optic.py @@ -42,6 +42,7 @@ SurfaceSagViewer, ) from optiland.wavelength import WavelengthGroup +from optiland.wavepropagation import BaseWavePropagator if TYPE_CHECKING: from matplotlib.axes import Axes @@ -129,6 +130,7 @@ def _initialize_attributes(self): self.fields: FieldGroup = FieldGroup() self.wavelengths: WavelengthGroup = WavelengthGroup() + self.wave_propagator: BaseWavePropagator = BaseWavePropagator(self) self.paraxial: Paraxial = Paraxial(self) self.aberrations: Aberrations = Aberrations(self) self.ray_tracer: RealRayTracer = RealRayTracer(self) diff --git a/optiland/wavepropagation/__init__.py b/optiland/wavepropagation/__init__.py new file mode 100644 index 00000000..1d8bd223 --- /dev/null +++ b/optiland/wavepropagation/__init__.py @@ -0,0 +1,4 @@ +# flake8: noqa + +from .base import BaseWavePropagator +from .asm import AngularSpectrumPropagator \ No newline at end of file diff --git a/optiland/wavepropagation/asm.py b/optiland/wavepropagation/asm.py new file mode 100644 index 00000000..6f9e09ee --- /dev/null +++ b/optiland/wavepropagation/asm.py @@ -0,0 +1,67 @@ +import optiland.backend as be +from math import pi + + +class AngularSpectrumPropagator: + def __init__( + self, + num_points: int, + dx: float, + evanescent: str = "clamp", + ): + self.dx = float(dx) + self.evanescent = evanescent + M_orig, N_orig = (num_points, num_points) + self.Mpad = int(round(M_orig * 0.5)) + self.Npad = int(round(N_orig * 0.5)) + M = M_orig + 2 * self.Mpad + N = N_orig + 2 * self.Npad + fx = be.fftfreq(N, d=self.dx) + fy = be.fftfreq(M, d=self.dx) + FY, FX = be.meshgrid(fy, fx, indexing="ij") + self.fx2_fy2 = (FX * FX + FY * FY)[None, None] + + def __call__(self, input_field, distance, wavelengths): + return self.forward(input_field, distance, wavelengths) + + def forward(self, input_field, distance, wavelengths): + if getattr(input_field, "ndim", None) != 4: + raise ValueError("input_field must have shape (B, C, M, N).") + + B = input_field.shape[0] + distance = be.array(distance).reshape(-1, 1, 1, 1) + if distance.shape[0] == 1 and B > 1: + distance = be.broadcast_to(distance, (B, 1, 1, 1)) + wavelengths = be.array(wavelengths).reshape(-1, 1, 1, 1) + if wavelengths.shape[0] == 1 and B > 1: + wavelengths = be.broadcast_to(wavelengths, (B, 1, 1, 1)) + + k = 2.0 * pi / (wavelengths * 1e-3) + + if self.Mpad > 0 or self.Npad > 0: + padded = be.pad( + input_field, ((self.Mpad, self.Mpad), (self.Npad, self.Npad)) + ) + else: + padded = input_field + + spectrum = be.fft2(padded) + argument = k * k - (2.0 * pi) ** 2 * self.fx2_fy2 + + if self.evanescent == "clamp": + kz = be.sqrt(be.clamp(argument, min=0.0)) + H = be.exp(1j * be.to_complex(kz) * be.to_complex(distance)) + elif self.evanescent == "decay": + kz = be.sqrt(be.to_complex(argument)) + H = be.exp(1j * kz * be.to_complex(distance)) + else: + raise ValueError('evanescent must be "clamp" or "decay".') + + out = be.ifft2(spectrum * H) + + if self.Mpad > 0: + out = out[:, :, self.Mpad : -self.Mpad, :] + if self.Npad > 0: + out = out[:, :, :, self.Npad : -self.Npad] + + return out diff --git a/optiland/wavepropagation/base.py b/optiland/wavepropagation/base.py new file mode 100644 index 00000000..b1aabd5c --- /dev/null +++ b/optiland/wavepropagation/base.py @@ -0,0 +1,154 @@ +from __future__ import annotations + +import warnings +from math import pi +from typing import TYPE_CHECKING + +import optiland.backend as be +from optiland.utils import resolve_wavelengths + +if TYPE_CHECKING: + from optiland.optic import Optic + + +class BaseWavePropagator: + def __init__(self, optic: Optic): + self.optic = optic + + def _build_propagator(self, num_points: int, dx: float): + from optiland.wavepropagation.asm import AngularSpectrumPropagator as Propagator + + return Propagator(num_points, dx) + + def compute_surface_phase(self, surf, X, Y, wl_arr): + interaction = getattr(surf, "interaction_model", None) + if interaction is None: + raise ValueError("Surface has no interaction_model") + + interaction_type = getattr(interaction, "interaction_type", None) + if interaction_type not in ("phase", "refractive_reflective"): + raise NotImplementedError( + f"Unsupported interaction_type: {interaction_type}" + ) + + if interaction_type == "phase": + phase_profile = getattr(interaction, "phase_profile", None) + if phase_profile is None: + raise ValueError("Phase surface has no phase_profile") + phi = be.array(phase_profile.get_phase(X, Y, wl_arr)) + else: + n = surf.material_post.n(wl_arr) if hasattr(surf, "material_post") else 1.0 + sag = be.array(surf.geometry.sag(X, Y)) + phi = ( + (2 * pi / (wl_arr * 1e-3))[:, None, None] + * (n - 1.0)[:, None, None] + * sag[None, :, :] + ) + + if getattr(phi, "ndim", None) == 2: + phi = phi[None] + + dphi_x = be.max(be.abs(phi[:, 1:] - phi[:, :-1])) + dphi_y = be.max(be.abs(phi[:, :, 1:] - phi[:, :, :-1])) + dphi_max = be.max(be.array([dphi_x, dphi_y])) + + if dphi_max > be.pi: + warnings.warn( + f"Phase under-sampled: phase jump > π between pixels. Max phase jump = {dphi_max:.2f} rad.", + RuntimeWarning, + stacklevel=2, + ) + + phase = be.exp(1j * phi) + if interaction_type != "phase": + phase = phase.conj() + return phase + + def create_input_field(self, X, Y, wl_arr, field, w0: float | str | None = "auto"): + angle_x, angle_y = field + + k = (2 * pi / (wl_arr * 1e-3))[:, None, None] + ax = be.array(float(angle_x)) + ay = be.array(float(angle_y)) + + phase = k * (X[None] * be.sin(ax) + Y[None] * be.sin(ay)) + field = be.exp(1j * phase) + + if w0 is not None: + if w0 == "auto": + dx = float(X[0, 1] - X[0, 0]) + w0 = ((X.shape[0] * dx) / 2) / 2 + r2 = X**2 + Y**2 + field = field * be.exp(-(r2 / (w0**2)))[None] + + return field + + def compute_field( + self, + z_target: float, + num_points: int, + dx: float, + field: list[tuple[float, float]] | tuple[float, float] | None = None, + wavelengths: str | float | list = "primary", + beam_waist: float | str | None = "auto", + ): + if field is None: + fields = [(0.0, 0.0)] + elif isinstance(field, tuple) and len(field) == 2: + fields = [field] + else: + fields = list(field) + + F = len(fields) + if F == 0: + raise ValueError("No fields provided.") + + wavelengths_resolved = resolve_wavelengths(self.optic, wavelengths) + wl_arr = be.array([float(w) for w in wavelengths_resolved]) + W = int(wl_arr.shape[0]) + if W == 0: + raise ValueError("No wavelengths resolved.") + + x = be.linspace(-(num_points // 2) * dx, (num_points // 2) * dx, num_points) + y = be.copy(x) + Y, X = be.meshgrid(y, x, indexing="ij") + + field_arr = be.zeros((F, W, num_points, num_points)) + 0j + for i, f in enumerate(fields): + field_arr[i] = self.create_input_field( + X=X, Y=Y, wl_arr=wl_arr, field=f, w0=beam_waist + ) + + propagator = self._build_propagator(num_points, dx) + current_z = 0.0 + wl_batch = wl_arr.repeat(F) + + def propagate(field_arr, dist: float): + flat = field_arr.reshape(F * W, num_points, num_points) + flat = propagator(flat[:, None], dist, wl_batch)[:, 0] + return flat.reshape(F, W, num_points, num_points) + + for surf in self.optic.surface_group.surfaces: + phase = self.compute_surface_phase(surf, X, Y, wl_arr) + field_arr = field_arr * phase[None] + + if getattr(surf, "aperture", None): + aperture = be.array(surf.aperture.contains(X, Y)) + field_arr = field_arr * aperture[None, None] + + t = surf.thickness if surf.thickness != float("inf") else 0.0 + + if current_z + t >= z_target: + remaining = z_target - current_z + if remaining > 0: + field_arr = propagate(field_arr, remaining) + return field_arr + + if t > 0: + field_arr = propagate(field_arr, t) + current_z += t + + if z_target > current_z: + field_arr = propagate(field_arr, z_target - current_z) + + return field_arr