{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6eeef840-c338-4ec9-ae10-a9e75b658053",
   "metadata": {},
   "source": [
    "# Energy and angular dispersive analysis\n",
    "\n",
    "`refnx` able to deal with reflectivity from systems containing energy dispersive materials, i.e. those whose optical properties change as a function of wavelength. For neutrons this mainly corresponds to elements with strong absorption effects. The treatment below is also able to deal with energy dispersive X-ray measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f2d9a637-f1a9-47c4-8ce7-dd9475032b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from refnx.reflect import ReflectModelTL, ReflectModel, SLD, MaterialSLD\n",
    "from refnx.util import q, xray_wavelength"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8ad07e9-8658-433a-995c-742ec780b91d",
   "metadata": {},
   "source": [
    "Start off by creating a `MaterialSLD`. This is a variant of `Scatterer`, whose optical properties are controlled by a formula, mass density (g/cc), and whether the material is being used for neutron or X-ray calculation.\n",
    "\n",
    "Here we'll calculate the SLD of the material at two different wavelengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9b54ac43-b365-4862-b22e-6d6392d5544c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SLD: (6.544686245235056+0.7318962182191739j) at 2.8 Angstrom\n",
      "SLD: (6.471928664642279+0.648390380281317j) at 18.0 Angstrom\n"
     ]
    }
   ],
   "source": [
    "gdgao_disp = MaterialSLD(\"GdGa5O12\", 7, probe='neutron') # can be 'x-ray'\n",
    "\n",
    "gdgao_disp.wavelength = 2.8   # Angstrom\n",
    "print(f\"SLD: {complex(gdgao_disp)} at {gdgao_disp.wavelength} Angstrom\")\n",
    "gdgao_disp.wavelength = 18.\n",
    "print(f\"SLD: {complex(gdgao_disp)} at {gdgao_disp.wavelength} Angstrom\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38228996-5edc-43c5-9019-54f92838e3b3",
   "metadata": {},
   "source": [
    "For comparison let's create a non-dispersive version. By non-dispersive we mean that the optical properties don't change as a function of wavelength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d976c1fe-a739-4446-8e84-5b66edd87fbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdgao_nondisp = SLD(6.5 + 0.68j)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d5078ff-1f36-4d03-a2dd-0029af54563c",
   "metadata": {},
   "source": [
    "The `MaterialSLD.density` attribute can be allowed to vary during a fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8d37781c-c504-446d-b40c-92fa9817ae00",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:   'density'   , value=7  (fixed) , bounds=[-inf, inf]>\n"
     ]
    }
   ],
   "source": [
    "print(gdgao_disp.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2ac79a0-0a43-441d-a782-b2e1cb419268",
   "metadata": {},
   "source": [
    "In comparison, with the non-dispersive analogue one can allow the real and imaginary part of the SLD to vary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0d240af4-fb86-4f68-9564-42c38da109f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "Parameters:       ''       \n",
      "<Parameter:   ' - sld'    , value=6.5  (fixed) , bounds=[-inf, inf]>\n",
      "<Parameter:   ' - isld'   , value=0.68  (fixed) , bounds=[-inf, inf]>\n"
     ]
    }
   ],
   "source": [
    "print(gdgao_nondisp.parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a66c1e4-8a88-4480-934b-b50696d3fc6c",
   "metadata": {},
   "source": [
    "Now we create two `Structure`s that are ostensibly the same, but one has a dispersive material in it, the other a non-dispersive analogue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "42637c4e-1c87-480d-b6ea-5754bc9d58e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "air = SLD(0.0)\n",
    "si = SLD(2.07)\n",
    "\n",
    "s_disp = air | gdgao_disp(300, 5) | si(0, 3)\n",
    "s_nondisp = air | gdgao_nondisp(300, 5) | si(0, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "851691b9-7e3e-421b-89db-84260f0efc5c",
   "metadata": {},
   "source": [
    "Now we generate `theta`/`wavelength` arrays, with a corresponding Q value.\n",
    "Subsequently we create a `ReflectModelTL` and a `ReflectModel`. `ReflectModelTL` is a variant of `ReflectModel`. Instead of calculating reflectivity as a function of Q (a. la. `ReflectModel`), it calculates as a function of incident angle and wavelength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f816cf11-9180-4f89-b1b3-cfbd137a59a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "npnts = 201\n",
    "theta = np.ones(npnts) * 0.65\n",
    "wavelength = np.geomspace(2.8, 18, npnts)\n",
    "qq = q(theta, wavelength)\n",
    "\n",
    "model_disp = ReflectModelTL(s_disp)\n",
    "model_nondisp = ReflectModel(s_nondisp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b0a42b5-0a46-4875-bad0-f7fb63019651",
   "metadata": {},
   "source": [
    "Now let's compare the reflectivity from the dispersive and non-dispersive analogues. The reflectivities are almost identical, the energy dispersive absorption effect has little effect in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a199c625-2073-4ffa-956d-8fdab09ce99d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGhCAYAAACphlRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPkElEQVR4nO3deXxU9aH+8c+ZSTLZE5KQjSSEHcKSsCUiKiBBREXRaq22ira1vW2wWqyt2l9V2tvqra1Sr7m1dSkutaJVccEFRAFFdmQNO4EEQhKyTfZt5vz+iEYpiwSSnJnJ8369zkuZnDnfZ2bQPHPO95xjmKZpIiIiIuIlbFYHEBEREekIlRcRERHxKiovIiIi4lVUXkRERMSrqLyIiIiIV1F5EREREa+i8iIiIiJexc/qAJ3N7XZTVFREWFgYhmFYHUdERETOgGma1NTUkJiYiM12+n0rPldeioqKSE5OtjqGiIiInIXCwkKSkpJOu47PlZewsDCg7cWHh4dbnEZERETORHV1NcnJye2/x0/H58rLl4eKwsPDVV5ERES8zJlM+dCEXREREfEqKi8iIiLiVVReRERExKv43JwXERERl8tFS0uL1THka/z9/bHb7Z2yLY8sL++88w533XUXbrebX/3qV/zwhz+0OpKIiHgB0zQpLi6mqqrK6ihyEpGRkcTHx5/zddg8rry0trYyd+5cPv74YyIiIhg7dixXX3010dHRVkcTEREP92VxiY2NJTg4WBcr9RCmaVJfX09paSkACQkJ57Q9jysv69atY/jw4fTp0weAGTNmsGTJEm644QaLk4mIiCdzuVztxUVfeD1PUFAQAKWlpcTGxp7TIaROn7C7cuVKZs6cSWJiIoZhsGjRohPWyc3NJTU1lcDAQLKysli3bl37z4qKitqLC0CfPn04cuRIZ8cUEREf8+Ucl+DgYIuTyKl8+dmc63ykTi8vdXV1pKenk5ube9KfL1y4kLlz5/LAAw+wadMm0tPTmT59evuupI5qamqiurr6uEVERHouHSryXJ312XR6eZkxYwb//d//zdVXX33Snz/66KPcdttt3HrrraSlpfHkk08SHBzMs88+C0BiYuJxe1qOHDlCYmLiKcd76KGHiIiIaF90XyMRERHf1q3XeWlubmbjxo1kZ2d/FcBmIzs7m9WrVwOQmZnJ9u3bOXLkCLW1tbz33ntMnz79lNu89957cTqd7UthYWGXvw4RERFfMHnyZO688872P6empjJ//nzL8pypbp2wW1ZWhsvlIi4u7rjH4+Li2LVrV1sgPz/+/Oc/M2XKFNxuN7/85S9PO/HK4XDgcDi6NLeIiEhPsH79ekJCQqyO8Y087mwjgCuvvJIrr7zS6hjHKa9toryumQC7DX8/G/52gwC7jQA/G/52G342Q8dZRUTEq/Xu3dvqCGekW8tLTEwMdrudkpKS4x4vKSkhPj6+O6N02I635tN/999pMv1pwI9m/GnCn2az7d+b8afF8KfVCMBl88dlBOC2++OyOXDbAjDtAZg2B6ZfAKbdgeHnALsDwz8Aw8+B4R+I4ReEzRGEzS8IuyMIuyMYe0AQAQEOAuw2HP62L/5pJ+iLJTDARuAX/+5v190eRES80eTJkxk1ahSBgYE8/fTTBAQE8F//9V88+OCDABQUFHD77bezbNkybDYbl156Kf/7v//bfiTjwQcfZNGiRdx111385je/obKykhkzZvDUU08RFhYGtJ1Q85Of/ITXX3+dsLAwfvGLX5yQIzU1lTvvvJM777wT0zSZN28ezz77LCUlJURHR3Pttdfy+OOPt6/7gx/8gLy8PN566y0iIyO57777yMnJ6fL3q1vLS0BAAGPHjmXZsmXMmjULALfbzbJly5gzZ845bTs3N5fc3FxcLlcnJD1RiMtJklEGZ7JzxfxicXfO2C2mnUYCaMSfJgJoNANoJIAaAmg0/b/4WQDNBNBqc9BqD/zin8G0+gXT6heC279tISAUW2AIQSGRhIZFEBYeSWR4OFFhDqJCAogKDsBPJUhEfIRpmjS0dM3vhW8S5G/v0B755557jrlz57J27VpWr17NLbfcwsSJE5k6dSpXXXUVoaGhrFixgtbWVnJycrj++utZvnx5+/P379/PokWLeOedd6isrOTb3/42Dz/8ML///e8BuPvuu1mxYgVvvvkmsbGx3HfffWzatImMjIyT5nnttdd47LHHePnllxk+fDjFxcVs2bLluHUeeeQR7rvvPubNm8cHH3zAHXfcweDBg5k2bVqH36+O6PTyUltby759+9r/nJ+fz+bNm4mKiiIlJYW5c+cye/Zsxo0bR2ZmJvPnz6euro5bb731nMbNyckhJyeH6upqIiIizvVlnGDsrJ+B8xpobQJXE2ZLI+6WJlpbGnE1N+FqacTV0oS7pQF3S1Pb0tqI2dKM2dqE2dr4xXObMVxNGK5mDFczNlcTNnczdnczdlcjfu4m/M0m/M2vzoH3N1z400AYDW0PfNN/C+4vllag6Ztfm8s0qCSMMjOCvWY41fZI6v2jaHRE4wqKxgyJxR4Wh394LMFRCfQKDycqJIDo0AB6BQcQ4KeyIyKeqaHFRdr9H1gydt5vpxMccOa/ZkeNGsUDDzwAwKBBg3jiiSdYtmwZANu2bSM/P7/9jNrnn3+e4cOHs379esaPHw+07QxYsGBB+56Wm266iWXLlvH73/+e2tpannnmGV588UWmTp0KtJWlpKSkU+YpKCggPj6e7Oxs/P39SUlJITMz87h1Jk6cyD333APA4MGDWbVqFY899pj3lZcNGzYwZcqU9j/PnTsXgNmzZ7NgwQKuv/56jh07xv33309xcTEZGRm8//77J0zi9Thh8W3LFwzA/sXSJdxuaG1sW1oaTvlPs6WB1uZ6WhsbaG2qp7W5HldzA+7mesymWmiuw2iuw2ipw95Sh721DntrPf6uehzutjJkN0xiqCbG+No1clq+WGqBY8dHqzGDKDPDKSCCz80InPZIGvyjaHZE4wqOwQhtKzuOiDjCIqOICg0kOiSgbc9OSACB/l32romIeK1Ro0Yd9+eEhARKS0vZuXMnycnJx10KJC0tjcjISHbu3NleXlJTU9uLy9efD217ZZqbm8nKymr/eVRUFEOGDDllnuuuu4758+fTv39/Lr30Ui677DJmzpyJn99X1WHChAnHPWfChAndcrZSp5eXyZMnY5rmadeZM2fOOR8m8nk2GwQEty2nYQD+Xywd5nZDSz001UB9Oa7aY9RXFNFYVUyzswSzthTqyvBvLCOwqZzg1kr8zRbCjAbCjAb68bW5S18vO1+73mCT6Uc54ZSZEew22/7ptPeiMSCKlqAYCI7BHt5WdEJ7xRETEUJMqIPYMAcxoQ6CAlR0ROTsBfnbyfvtqS+30dVjd4S///H/JzcMA7f7zOcfnOvz/1NycjK7d+/mww8/ZOnSpfz0pz/lkUceYcWKFSeM1d088mwj6SY2GzhC25bwBOxA2BfLSZkmNDqhrgzqSnHXltJYWUxjVTEt1aWYtaXYGsoIaCgjqLkch7seh9FKIhUkGhXHb+vLslMNFLc95DYNKgij3AznkBnBOnpRZYuiwdGblqBY3GFx+IUn4OiVSGREJHERgSREBJIQHkR4kJ/O9hKRExiG0aFDN55o2LBhFBYWUlhY2L73JS8vj6qqKtLS0s5oGwMGDMDf35+1a9eSkpICQGVlJXv27GHSpEmnfF5QUBAzZ85k5syZ5OTkMHToULZt28aYMWMAWLNmzXHrr1mzhmHDhp3Ny+wQ7/5Ev6arJ+wKYBgQFNm2xAzEBgR/sZxUSwPUHWtbao9h1pXSVFVCk7OY1upS+KLsOJrKCWypwva1w1dDOPzVdpq/WJxfPVRtBnHMjKTU7MVuIik3omgMjKUlJB4ikvCPSiY8pg9xkaEkRAQSHxFIVHAANpsKjoh4l+zsbEaOHMl3v/td5s+fT2trKz/96U+ZNGkS48aNO6NthIaG8oMf/IC7776b6OhoYmNj+fWvf43Nduo5iwsWLMDlcpGVlUVwcDAvvvgiQUFB9O3bt32dVatW8cc//pFZs2axdOlSXn31VRYvXnzOr/mb+Ex56eoJu3IW/IMgMqVtoe0QV+AXywncLqgvh9pSqCvFrCmhqaqIpoojtDqLobYY//pSApuOEeBuJNxoINxoYABHv9rGlyWnEjjYdpZWCb0oMqNZZUZTQjR1gfG4QhPxi0omKHYAcbFxJEcFkdwrmN5hDu29ERGPYxgGb775JrfffjsXXXTRcadKd8QjjzxCbW0tM2fOJCwsjLvuugun03nK9SMjI3n44YeZO3cuLpeLkSNH8vbbbx934di77rqLDRs2MG/ePMLDw3n00UdPe1X8zmKY3zRBxct8WV6cTifh4eFWx5HOZpptc3RqS6CmGGpLaKk6Qn3FEVoqj2BUFxFQV0Rw0zHsfPNeuCozhAIzlgIzjiNGHLVBSbRGpOAX3Y+w2FT6xkbQPyaElOhgHH6afyPiyRobG8nPz6dfv34EBp70a5J0oq9fE+ZMne4z6sjvb5/Z8yI9hGFAYHjbEjMIaJusfMK+NrerrdxUHwHnYVorC2koP0RLxWEM52Ec9UUEt1QSadQRaeQzivy25zXRNuG4FFrzbBwxYygwY1lLPBWBfWmKHIB/7BCi+vSnX+8w+vcOJSE8UIejRES6kcqL+CabHSL6tC3JmfhxkonITbVQdQgqD+Iqz6euZB+tZfnYnYcIqT+CH830NUrpSykXsr1tgvGxtqVhewAHzXg+NxN5y9aHhvD+2GKHEdF3BIMSohkcH0rvUB2GEhHpCiov0nM5QiFuOMQNxw4ct5PS7YbaYqg8iFlxgPqju2ks3o1fxT5C6goIoplhRgHDKGhbv6Ztad1n44CZwBozhUN+qdRHDsGeMJKkvoMYkRTJ4LgwXdRPRHzCwYMHLRvbZ8qLzjaSTmWzQXgihCdi9D2fkNHQfp9VV2vbHpvyfbhKd1N7ZCeu0t0EO/cS2FrNYOMIgzkC5uq2ycOVUL0jiN1mMq+YqRwLHw59xhDffyQjkqIYHB+q+TQiIh2gCbsincU0oboISvNoKdpGXeFWjNIdhNYcwG62nrB6rRnIDjOVbeYAysKHE5Aylv6DRjA2NYqkXkE65CTSQZqw6/k0YVfE0xhG+zwb/0HTiPzy8dZmKN+LWbydmvwNtBRuIqxyB6HuBrKMXWSxC+oWw06ozAtls3sAb/kPpzEhk6hB55HeP56RfSJ013ARkS+ovIh0Nb8AiBuOETec8PTr2x5zu+DYbswjG6nN34Dr8EZCq3bRi1qm2Lcwxb0FjrxE02E/tn3Un+eMYdTEjiNy2EWMH9KPYQnh2HWGk4j0UCovIlaw2SEuDSMujbAxN7U91toEJTtoObSGmt0rcRxdT0hzGeOMPYxjDxx7k9ZSG1uX9+dp2yhqEi4gbviFTBqWREr06e+BJSLiS1ReRDyFnwP6jMG/zxiizv9p2xyaynzch1bj3LUSW8FnRDQUMMbYxxj2wdHXqS9ysO6DobwbPAZz4CWkZ4xnXL9ondEkIj7NZ8qLzjYSn2MYENUfW1R/eo3+bttjVYW49i/HuWMpgYWfENxSwWT7FiY3bYEd/yB/Wxz/MsZR3udiUjIuZtqIZCKCrb37q4icvcmTJ5ORkcH8+fPP6oq2vspnyovubSQ9QmQy9rE3ETX2prY9M6V5NOxeRu329+l1bC39bCX0YzEULcZ5JJhP3kmnMHYKcWNncXF6PyKDA6x+BSJyltavX09ISMg3r9gD+Ex5EelxDAPihhMUN5ygi34GTTW49y2javPbBOZ/SERrFVcYq6FsNQ3v/4mP3x3NvrjpJIydySXp/bRHRsTL9O7d2+oIHkMHxkV8hSMM2/BZRH33GYLvOwDfX0LlmByqApMJMpq5zL6Wn5X9lhnvX8jyh2fx17/nsiLvMC63T13qScRr1dXVcfPNNxMaGkpCQgJ//vOfj/t5amoq8+fPB8A0TR588EFSUlJwOBwkJibys5/97Lh1f/e733HDDTcQEhJCnz59yM3N7c6X06W050XEF9nskJJFr5QsmPl7OLqZqnUvY9+5iLCmYq4yPoWiTylf+BAL7ZNpGHEjUy68iP69Q61OLtK5TBNa6q0Z2z+4bQ/pGbr77rtZsWIFb775JrGxsdx3331s2rSJjIyME9Z97bXXeOyxx3j55ZcZPnw4xcXFbNmy5bh1HnnkEe677z7mzZvHBx98wB133MHgwYOZNm3aub4yy6m8iPg6w4DE0UTOGg1XPoR5eB3la1/GsfstolvLudH9Nmx9mw2bB/Nu5OUkX3Ajl44ZoFsWiG9oqYc/JFoz9n1FEHBmc1Rqa2t55plnePHFF5k6dSoAzz33HElJSSddv6CggPj4eLKzs/H39yclJYXMzMzj1pk4cSL33HMPAIMHD2bVqlU89thjPlFedNhIpCex2TBSziPmuvmE3buHluv/RXHCVFzYGGfbw5zqx7h48UW88/vree71tznqbLA6sUiPsH//fpqbm8nKymp/LCoqiiFDhpx0/euuu46Ghgb69+/PbbfdxhtvvEFr6/G3IZkwYcIJf965c2fnh7eA9ryI9FR2P/yHXUb8sMugppiatc/j2vAckY2H+Za5FLYuZc3mYbzV50YmzPguo1KirU4s0nH+wW17QKwau4skJyeze/duPvzwQ5YuXcpPf/pTHnnkEVasWIG/v+9PxveZPS+5ubmkpaUxfvx4q6OIeJ+weMKyf0nkL7fRetPbFCVdRit2zrPt5MdHf0Pk01k8/9gv+WTHQXzsXq7i6wyj7dCNFUsH5rsMGDAAf39/1q5d2/5YZWUle/bsOeVzgoKCmDlzJo8//jjLly9n9erVbNu2rf3na9asOW79NWvWMGzYsA68eZ7LZ/a86DovIp3AZsNvwEUkDrgInEco+ziXoG0vkMIxbnb+jYpX/skLQVeRMO1nTB09BJvuryTSKUJDQ/nBD37A3XffTXR0NLGxsfz617/GZjv5PoYFCxbgcrnIysoiODiYF198kaCgIPr27du+zqpVq/jjH//IrFmzWLp0Ka+++iqLFy/urpfUpXymvIhIJ4voQ8ysP8Bl/4/KNS/gXvU40U2Hubnxn9S+9RqvLbmc3pfcxaQxaRgd+IYpIif3yCOPUFtby8yZMwkLC+Ouu+7C6XSedN3IyEgefvhh5s6di8vlYuTIkbz99ttER391ePeuu+5iw4YNzJs3j/DwcB599FGmT5/eXS+nSxmmj+0D/nLPi9PpJDw83Oo4Ir7D7aLu839Tv+yP9K7fB0Cd6eCd4Fn0ueyXTBwxQCVGLNXY2Eh+fj79+vUjMDDQ6jiW8tRbCZzuM+rI72+fmfMiIl3MZidk7PX0vnsDtdf8k6MhaYQYTVzfsJAR/76IhY/dxc6CEqtTikgPoPIiIh1jGISOuoKEX3xGzVX/4FhgPyKNOr5T/QzRz2Sy8O8PU+K06KJgItIjqLyIyNkxDMJGX0PvX26kfNpfKPOLJ9ao4vqihyh9dCKvLXqN5la31SlFeqSDBw963CGjzqTyIiLnxmYneuItxNyzlcPj76XeCGKkcYBvbf4+nzw8i41bt1udUER8jMqLiHQOPwdJl99D0NzNHEz5Fm4MprauYPBr2bz+1/spq9ahJOkePnYeik/prM/GZ8qLLlIn4hmMsHhSv/8sdTd/SEHwcMKMBq4p+QtHHr2Ij1Z8rF8s0mW+vLJsfb2Ksqf68rM516sA61RpEek6bjeFS58gevUfCKaBFtPOh1HfIXP2w0RH6r9P6XxHjx6lqqqK2NhYgoODdfq+hzBNk/r6ekpLS4mMjCQhIeGEdTry+1vlRUS6XEvlYQ69mMPA8uUA7KEvpdMe54KJky3NJb7HNE2Ki4upqqqyOoqcRGRkJPHx8SctlSovKi8iHqng05cJX/ZLIk0nTaYfS+J/xNTvzyPYEWB1NPExLpeLlpYWq2PI1/j7+2O320/5c5UXlRcRj9VUdZTCBT9kYNWnAGy2jyT8hqfpP3CoxclExEq6wq6IeCxHZAID73iH/ef9gXoCyXBtI+qFi1m5+EVN5hWRM6LyIiLdzzAYcGkOjT9Yzv6AIUQadVy0PoelT9xOQ2Oz1elExMOpvIiIZaKSh9HvFyvZlvhtAC4pf4Hdf86m5GiBxclExJOpvIiIpWwBgYz80VPsu+AxGnCQ0bIF42+T2LnhY6ujiYiHUnkREY8wMPv7VH33Aw7ZkomlgtS3v81nbz1jdSwR8UAqLyLiMRIGjSbmzpVsD84kyGjm/E1zWf70PZhu3eBRRL6i8iIiHiUkPIq0ue+yKeE7AEw+/Fc2/OU7tDQ1WJxMRDyFyouIeBybnz9jfvw3No74f7SaNsY7P2Dfo9Opr6mwOpqIeACfKS+6MaOI7xl77d1sn/w0tWYQw5q2UPSXaVQeK7I6lohYTFfYFRGPt3PTJ8S9dQNR1HDIlkTwD96id58BVscSkU6kK+yKiE8ZNuZCam54h2Ji6Os+jOvp6Rw7uMPqWCJiEZUXEfEKfYdk4L71PQqMROLNY/gtmMGx/ZusjiUiFlB5ERGvkdh3MPYffMBuoz+9cGJ/YRal+VutjiUi3UzlRUS8Sp+kFEJ/tJg9Rj+icGJ7/kpK8nUISaQnUXkREa/TJyGRsB+9w34jhRizEuP5mTiP7LU6loh0E5UXEfFKCQlJBP9wMQeNPsSa5TQ+czkNxw5aHUtEuoHKi4h4rYQ+KbhvepMC4olzl1D1t8tpqSmzOpaIdDGVFxHxav37D6LqutcpMmNIaD3M4f+7EndTvdWxRKQLqbyIiNcbNXw4hZc/j9MMoV/DDnb/9TvgdlkdS0S6iMqLiPiErMyJbLng/2gy/RhWtYJdC3LAty4gLiJfUHkREZ9x0bRZLE/7HQBDC/7FwbcftjiRiHQFlRcR8SmXfPsnLIr9KQCpmx6mdO2rFicSkc6m8iIiPsUwDC697b95O+hKAMLem0Nd4TaLU4lIZ1J5ERGfE+hvJ+vHf2W9MZIgGql97npcdZVWxxKRTqLyIiI+KTYylODvPs9hM4a41iMcfOoGnYEk4iN8przk5uaSlpbG+PHjrY4iIh5i+MD+7JvyJA1mAAOqVnPo37+2OpKIdALDNH3rXMLq6moiIiJwOp2Eh4dbHUdEPMCr/3iU6w7NA6DyimfoNe5aixOJyH/qyO9vn9nzIiJyKjO/dwevB84CwLH4dlrLDlgbSETOicqLiPi8QH87Y77/OBvNoQSb9ZQu+B64WqyOJSJnSeVFRHqE1NgIKmfk4jSDSazdQcFrmv8i4q1UXkSkx8g+bxyL+90HQFLe36navsTiRCJyNlReRKRHuea7P+XdgEuxYcIbP8asPWZ1JBHpIJUXEelRAv3t9P/e4+wz+xDpqqD4+e/rBo4iXkblRUR6nKEpcWwY92eaTH8SSlfi/ORJqyOJSAeovIhIj3TtZdN5PvRWAAI+nodZecjiRCJyplReRKRH8rPbmHLzb9hoDiHIbODoP/9Lh49EvITKi4j0WAPjwtl//v/QaPqTWPYZ5Z8+a3UkETkDKi8i0qNdO20yr4bfDEDgR7/BdB6xOJGIfBOVFxHp0Ww2gwtvepAt5gBCzDpKXvqpDh+JeDiVFxHp8VJjw9k29g80m3biS5ZTv+llqyOJyGmovIiIANddNo0XHd8BwHz3V1BfYXEiETkVlRcREcDhZ2fotb9hlzuZEJeT8rf+n9WRROQUVF5ERL5w/uAElqT+AoBeu17CVbjR4kQicjIqLyIiX3PDdTfwNhdiw6Ty1dvB7bI6koj8B5UXEZGv6R3moHHKPKrNIGKqd1Cz7kWrI4nIf1B5ERH5D9dcOIZ/h9wAgPnhPGiusziRiHydyouIyH+w2wyGX303Be7ehLeWU770T1ZHEpGvUXkRETmJrEGJvBv/EwBCN/wfVB+1OJGIfEnlRUTkFC659kdsdA/GYTZS8qZOnRbxFB5ZXq6++mp69erFtddea3UUEenB+seGsXFo26nTvfe/hrtoi8WJRAQ8tLzccccdPP/881bHEBHhW1fO4l1zIjZMyl6/W/c9EvEAHlleJk+eTFhYmNUxRESIDnVQcf69NJn+xJatpSnvXasjifR4HS4vK1euZObMmSQmJmIYBosWLTphndzcXFJTUwkMDCQrK4t169Z1RlYREUtce/H5vOI3E4Dad38DbrfFiUR6tg6Xl7q6OtLT08nNzT3pzxcuXMjcuXN54IEH2LRpE+np6UyfPp3S0tL2dTIyMhgxYsQJS1FRUYdfQFNTE9XV1cctIiKdKdDfTtQld1NtBhNdt5/6La9ZHUmkR/Pr6BNmzJjBjBkzTvnzRx99lNtuu41bb70VgCeffJLFixfz7LPPcs899wCwefPms0t7Eg899BDz5s3rtO2JiJzMpeOH8cJHs7il6SUal/ye4PRrwGa3OpZIj9Spc16am5vZuHEj2dnZXw1gs5Gdnc3q1as7c6h29957L06ns30pLCzsknFEpGez2wxip91JlRlCVEM+9Z+/YnUkkR6rU8tLWVkZLpeLuLi44x6Pi4ujuLj4jLeTnZ3Nddddx7vvvktSUtJpi4/D4SA8PPy4RUSkK0wfM5jXA68GoPHDP4Cr1eJEIj1Thw8bdYcPP/zQ6ggiIiew2wwSLrmTircWEdVQQP2mlwke/z2rY4n0OJ265yUmJga73U5JSclxj5eUlBAfH9+ZQ4mIWGL66IG8FvQtAJqXPQSuFosTifQ8nVpeAgICGDt2LMuWLWt/zO12s2zZMiZMmNCZQ50gNzeXtLQ0xo8f36XjiEjPZrMZJF9yB2VmOJGNh6lf/0+rI4n0OB0uL7W1tWzevLn9jKH8/Hw2b95MQUEBAHPnzuWpp57iueeeY+fOnfzkJz+hrq6u/eyjrpKTk0NeXh7r16/v0nFERC7J6M/rQdcB0PLxw9DabHEikZ6lw3NeNmzYwJQpU9r/PHfuXABmz57NggULuP766zl27Bj3338/xcXFZGRk8P77758wiVdExFvZbAZ9p8/h2KJ/07vpKPWbXiE4U3NfRLqLYZq+daOO6upqIiIicDqdOvNIRLqM222y4H9+xvebnqc8ZCDRv9gAhmF1LBGv1ZHf3x55byMREU9nsxnETP4v6kwH0XX7aNmz1OpIIj2Gz5QXTdgVke526fhhvGWfBkDFkj9ZnEak5/CZ8qIJuyLS3QL8bLiy/otW00Zc+VrMos+tjiTSI/hMeRERscLMi87jfdouBVH6/iMWpxHpGVReRETOQUSQP0eH/wiAmIL3oPKQxYlEfJ/Ki4jIObrskul84h6JHTdlHz5mdRwRn6fyIiJyjvpEBrGt72wAwvL+BfUVFicS8W0+U150tpGIWOmi6dexw90Xh9mI89O/WR1HxKf5THnR2UYiYqURSZGsiL4eAGP9M+BqtTiRiO/ymfIiImK1Ydk3U2aGE95yjOad71odR8RnqbyIiHSSi4Yl8Z5/NgAVy//P4jQivkvlRUSkk9htBubYW3GbBvFlq6F8v9WRRHySyouISCe67MIsVpgZAJQt/6u1YUR8lM+UF51tJCKeICbUwe7kbwMQnPcytDRYnEjE9/hMedHZRiLiKTKmXEuhuzfBrhoaPn/V6jgiPsdnyouIiKfIGtCbD4JmAFD36ZMWpxHxPSovIiKdzDAMQs67hSbTj5jqHZhHNlkdScSnqLyIiHSBy84bxRIzC4CyjzVxV6QzqbyIiHSBiCB/Cgfc2Pbv+9+EhiprA4n4EJUXEZEucsHFl7PLnUyA2UTtplesjiPiM1ReRES6yKjkXqwKnQZA/boXLE4j4jt8przoOi8i4olCxt2AyzSIdW6Fsn1WxxHxCT5TXnSdFxHxRNmZ6XxipgNQvvp5i9OI+AafKS8iIp4oJtTB7rjLAfDbthDcbosTiXg/lRcRkS6WOvE6qs1gIpqLceV/YnUcEa+n8iIi0sUmj0hhiTEBgLJVC6wNI+IDVF5ERLqYw89O+cBvARCZ/x401VqcSMS7qbyIiHSDcRfM4KA7DofZQOO2N62OI+LVVF5ERLrBmL69+DjwYgCca3TWkci5UHkREekGhmFgpN8AQO+yteA8bHEiEe/lM+VFF6kTEU83dcJ41riHYcPEufZFq+OIeC2fKS+6SJ2IeLrkqGA+73UpAK2bda8jkbPlM+VFRMQbxGddR7NpJ7p+P5TkWR1HxCupvIiIdKOLMwZ/dbuAdS9bnEbEO6m8iIh0o4hgf/bHTgfAtuMNME2LE4l4H5UXEZFulpB5NY2mP70aCzCPbrY6jojXUXkREelmk0f152NzDADlaxdanEbE+6i8iIh0s7BAfw4ltJ115L9zkQ4diXSQyouIiAWSs2ZRbzqIaD6KWbTZ6jgiXkXlRUTEAlNGpLDSzACgfMOr1oYR8TIqLyIiFggO8ONQXDYA9l1v69CRSAeovIiIWKT3mCtpMv3o1VAApTutjiPiNVReREQsMmlkPz51jwTAuek1i9OIeA+fKS+6MaOIeJvoUAe7oqYA4NrxlsVpRLyHz5QX3ZhRRLxRePqVtJh2omr3QPl+q+OIeAWfKS8iIt5oUvoQ1rmHAlC/fbHFaUS8g8qLiIiFUqKD2RYyAYDabe9YnEbEO6i8iIhYLCDtcgCiyzZAo9PiNCKeT+VFRMRiWePGss+diB0XzbuWWB1HxOOpvIiIWCwtIZy1/pkAlG3SWUci30TlRUTEYoZh0DzgEgAijiwHV6u1gUQ8nMqLiIgHGDo+myozhBBXNa0Fa6yOI+LRVF5ERDzA+P69WWWMAaBk/ZsWpxHxbCovIiIewM9uo6LPxQAEHNCkXZHTUXkREfEQvUdfRotpp3fjQag4YHUcEY+l8iIi4iHOS+vPenMIAFWb37Y4jYjnUnkREfEQkcEB7A6/AIDGHbpVgMipqLyIiHgQ29BLAehdrqvtipyKyouIiAcZnTGW/e4E7Lho3fOh1XFEPJLKi4iIBxmRGMEq+3gAKjXvReSkfKa85ObmkpaWxvjx462OIiJy1mw2g5rkKQAEFa4A07Q4kYjn8ZnykpOTQ15eHuvXr7c6iojIOUkeNYV600FoSwWU7LA6jojH8ZnyIiLiKyYOTWSNexgANXm6YJ3If1J5ERHxMNGhDg6Etx0Cr8tbanEaEc+j8iIi4oH8BmcDEFW+AVoaLE4j4llUXkREPNDI9EyOmlEEmM24Dn5mdRwRj6LyIiLigTJSerHWSAfg2Ob3LE4j4llUXkREPJDdZlAZ33arAFv+cmvDiHgYlRcREQ8VPeoSAGLr90JNicVpRDyHyouIiIc6b+RgtrlTAajbqVsFiHxJ5UVExEPFhgWyPXAsAFXbP7A4jYjnUHkREfFgjSmTAAgv+lS3ChD5gp/VAURE5NQSR06mfq+DsNZyKM2DuOFWR+o2FTX17Nu1FeehLbgqDmJrqsVorcftH4LpiMAenUpkajoDB4+gV2ig1XGlG6m8iIh4sMyBCax1D2WKfQs1eUsI8/HyUlpeyZalLxCUv4SRjRvJNOpPvfIRYCtUmKGsCBxHQ2o2I7JvJKl3dLflFWuovIiIeLBeIQHsCc1kSsMWGnYuJWzKz62O1CX27N/HoXceYVzFO0wzatseNKCBQEocqdSHpWI6IiAgCJrrsDVWElpzkLjmg0QZtUxqWg67l1O16yGWhF9CTPbPGZOebulrkq6j8iIi4uHM/lNgx1P0OrYeWhrB33cOkZSUl7Pl5d8xsfQlBhtNYECpLZbSfleRMO4qooecT6rNfuoNuFoo37WKko1v0vvQYnq7Srik5g2aX3+LpR9cSvyV9zNy6NDue0HSLQzT9K0ZYNXV1UREROB0OgkPD7c6jojIOVuWV0zawgkkGBVw85vQf7LVkc6ZaZp8+uFb9F91F304BkB+4DDsF91FynnXwOkKy6m4XRR//i61H89nYO0GAGrMIFYk/ZgLb7yXiBDfKX2+qCO/v3W2kYiIh8vsH80ad9tcl+pdH1uc5tzVNTTy4eM/ZuKns+nDMUpsseRPyaXfr1aTcv51Z1dcAGx24sfOZOAvllF67escCkwjzGjgiiPzKfzTBWzavKlzX4hYRuVFRMTDhQX6UxjZdr2X5r3LrQ1zjo4cLWLXn6czrXIhNsNkR9xVRP9iPf0mfQ8Mo9PGiR0xlb6/XMX+zN9SSzAjzL0MfOMy3vlXLi63Tx1w6JFUXkREvID/gLbrvfSq3AZNtRanOTu7d22n5W9TGdu6mQYc7J/yfwz/yfP4BUd2zYA2GwMuuwP7nDXkB48k3Gjgit33seQvP6ausblrxpRuofIiIuIFRgwfSYG7N3ZcmAWrrY7TYXk7thD28lWkUkSprTc1Ny5mwKTvdsvYQTF96XfXcnYP/CEAM5wL+fzPV1FSXtEt40vnU3kREfEC4/pGsZYv5r3s/MjiNB2Tt2Mz0a/MIpEyjtiTCP7JR8QOHt+9Iex+DPnen8m/6DFa8OOCls84mnsFRaXHujeHdAqVFxERLxAUYKe4V9sv/Nb9KyxOc+YOFRwk/NXriDMqKPTrS2TOUkJ7p1iWp9/F36fimleoJZgM9w7KnryCw0d1x25vo/IiIuIlHIMmA9DLuRMaqizNcibKKsqpX/AtkijlqC2eqJ++T0hUotWxiBs1lcYbXqOGEEa5d1H595mUlJVbHUs6wOPKS2FhIZMnTyYtLY1Ro0bx6quvWh1JRMQjjEobxn53AjbcmIdWWR3ntJpbWtn/5I0Mc++jinActyzyiOLypZgh59PyvbdwEsZIczeHnryOqpo6q2PJGfK48uLn58f8+fPJy8tjyZIl3HnnndTV6S+UiMjolEjWfzHvpSrPs+e9rPzHr8lqXkMT/tR960WiUoZZHekEUQPH0XjdSzQSQGbrRrbkfo/G5larY8kZ8LjykpCQQEZGBgDx8fHExMRQUaEZ4SIiDj87x2KyADDzP7E4zamt+vANphz5GwD7xz9In5GTLE50anHDL+LYjKdoxcakxo/46Klf4WMXnvdJHS4vK1euZObMmSQmJmIYBosWLTphndzcXFJTUwkMDCQrK4t169adVbiNGzficrlITk4+q+eLiPiakCGTAYiq2Q11njdP43BhAUM+uQO7YbKt90zSLp9jdaRvlJw1i4NZvwPgsmNPs/T1Zy1OJN+kw+Wlrq6O9PR0cnNzT/rzhQsXMnfuXB544AE2bdpEeno606dPp7S0tH2djIwMRowYccJSVFTUvk5FRQU333wzf//730+bp6mpierq6uMWERFflTFsELvcbV/o3Ac9a++Ly+XmyIs/JsZwcsgvlbQfnv7/355k4Iw57Ey+HoDzt97Hxg3edy2dnuScbsxoGAZvvPEGs2bNan8sKyuL8ePH88QTTwDgdrtJTk7m9ttv55577jmj7TY1NTFt2jRuu+02brrpptOu++CDDzJv3rwTHteNGUXEF7W43Cz87Q18z3if8rTZRH/7casjtft44V+YsvN+Wkw75Td+QPyQbr6WyzkyW5vZ/+glDKz/nH2kEPGzT+gdFWl1rB7DshszNjc3s3HjRrKzs78awGYjOzub1avPrMWapsktt9zCxRdf/I3FBeDee+/F6XS2L4WFhWedX0TE0/nbbVTGts17sR1caXGar+Qf2MvYvIcB2Dkkx+uKC4DhF0CfH75EpRHBQArY+kwObt0HySN1ankpKyvD5XIRFxd33ONxcXEUFxef0TZWrVrFwoULWbRoERkZGWRkZLBt27ZTru9wOAgPDz9uERHxZRFDp+A2DXrV50ON9RdYM02Tkld+TrhRz/6AoYy8/n6rI521oKhE6i//K27TYGrdO3z4mvcc+upJPO5sowsuuAC3283mzZvbl5EjR1odS0TEY4wZ2p88sy8ArgPWX2139QcLOa/xE1pNGyHX5mLY/a2OdE76jLucXQN/AMD47b9jf/4BixPJf+rU8hITE4Pdbqek5PhvAiUlJcTHx3fmUCIiPdawhHA22UYAULnD2uu9OGtqSF7zAADbk24gfvA4S/N0lmE3/IEC/wH0Mmoo+VcOra0uqyPJ13RqeQkICGDs2LEsW7as/TG3282yZcuYMGFCZw51gtzcXNLS0hg/3vuOs4qIdITdZuCMb/t/qn+BtWccbfrXPJIppsyIYtiNf7A0S2cy/ByEfOcpWrBzfvNnLHvtSasjydd0uLzU1ta2H84ByM/PZ/PmzRQUFAAwd+5cnnrqKZ577jl27tzJT37yE+rq6rj11ls7Nfh/ysnJIS8vj/Xr13fpOCIiniBq2CRaTRsRjYehqsCSDIUF+WQeeR6AYxPuxxESaUmOrhI9YCx7h/wXAJl5D3FIJ4R4jA6Xlw0bNjB69GhGjx4NtJWV0aNHc//9bRO0rr/+ev70pz9x//33k5GRwebNm3n//fdPmMQrIiJnb/zQVLaYAwBo2bfckgwH//0bQowm9gUMZdi0WyzJ0NWGXfcghf6p9DJq2L9QV9/1FB0uL5MnT8Y0zROWBQsWtK8zZ84cDh06RFNTE2vXriUrK6szM4uI9HiDYkP53D4KAGfeh90+ft6WdZzvfAcAv0t/D4bR7Rm6g+EXgH3mYwBMrnmX1SvetziRgAeebSQiIt/MMAzq+kwEIPDwKujGPQKmaVL/7v/DbphsD7uQ1DHZ3/wkL5Y46mK2x16BzTCJXnEvdQ1NVkfq8XymvGjCroj0NHFpF9Fo+hPaXAZle7pt3C1rPmJc01pcpkHsNQ9127hWGnjjo1QTyhAzn1X//ovVcXo8nykvmrArIj1N1uBENrgHA9C8d3m3jGmaJu6P284q2hY9g9h+PeM6XIGRcRwedTsAGftyKTnmeTfF7El8pryIiPQ0qdHBbAvIAKB6Z/fMe9n82RLGNG+g1bSRPOuBbhnTUwy78ucctScQa1Sx5ZXfWh2nR1N5ERHxUoZh0JR8AQChRavB3bUXUjNNE1a03b9oW8xlRKcM7dLxPI3h56BhUtuZtReWvsTuvbstTtRzqbyIiHixpLQJVJtBBLpqoHhrl461bd3HjG7eRItpJ/lq771/0bnof+EN7A8cTpDRTNFb2vtiFZ8pL5qwKyI90YTB8ax1DwOgcU/X3iqgacWjAGyPmkZM0pAuHctjGQaBM9pKy8Tq98jLO/WNg6Xr+Ex50YRdEemJ+kQGsTtoDAC1uz7usnEO7N7K2LpPAYi79O4uG8cb9EnPZk/IOAIMF6WL/9vqOD2Sz5QXEZGeyux3IQDhpeuhtblLxih+/0/YDJPtwVkkDvGNmy+ei/AZbZOVL6hdwratmyxO0/OovIiIeLlBIzIpM8MJcDfCkQ2dvv3SogLGVrwLgGPyXZ2+fW8UP+Iidoedh5/hpvL9nnGtG0+i8iIi4uUmDOjNancaALU7l3X69ve+8ygOo4W9/kMZNP6STt++t4qc0TZpeULdMnbv3mlxmp5F5UVExMtFBPtzKKLtZIWGPZ0776W6upIRR14BoDEzx2fvYXQ24tImsid4NP6GiyPv/cnqOD2KyouIiA8IGDgZgF4VW6C5rtO2u+2tJ4gw6jhsS2T4lBs7bbu+InBK22G08yrf5lBhocVpeg6fKS86VVpEerIRw9M5bMbgRytm/spO2WZzUxP99y0AoHj4bdj8/Dplu74kZdwVHAoYSLDRxN53HrU6To/hM+VFp0qLSE82JjWKleZoAKq2vNMp29z8/rMkUEYFEYy8/Medsk2fYxi0TLgDgLHFr1JeWWVtnh7CZ8qLiEhPFuhv52jcZAD89y8F0zyn7ZluNzFbngRgb7/v4ggMOdeIPmvARTdSYoujl1HD5veetjpOj6DyIiLiI+LTp1FvOghtKoGS7ee0rW0r36C/+yD1poOhV87tpIS+ybD7UTrsJgCS975AS2vX3mNKVF5ERHzG5BEprHKPAKB+++Jz2pZ99eMAbI27mohevc85m68bfOlPacDBYPMg61Z0zmE7OTWVFxERH9EnMoi8sAkANOx496y3c2DLJwxv2kyraSPlcl2U7kw4wqLZG385ALZ1T1qcxvepvIiI+JDAYZcC0KtyK9QeO6ttOD/8MwCbwqeS2Hdwp2XzdX2mt03czWxcze7deRan8W0qLyIiPiQrYyTb3KnYMGnd9X6Hn196aBejqpcDEJmtvS4dEd0vgz3Bo7EbJoc/+rvVcXyaz5QXXedFRARG9YngU7+2Q0c1G17q8PML330Eu2HyuWMcg9MndHY83zf6ZgDSit+irqHJ4jC+y2fKi67zIiICNptB9aCrAYgoXg3Ow2f83NqKo6QVvwWAe8LPuiSfrxs46QachJFglLPxo39bHcdn+Ux5ERGRNpOyxrHalYYNk+ZNZ773Zc/bjxFkNLPLNpDRF87swoS+yxYQxMGktvcuYMsLFqfxXSovIiI+JqtfFMuDpgLQtOGfZ3TBupZ6J/3z24pO2agfY7Pr18PZSs7+CQDjmtayd/9ei9P4Jv3tFBHxMYZh0GvstdSbDsLqDsLhbz6cvuvfvyOSGg6RwLgZs7s+pA+LSh3F/sDh+BluCj561uo4PknlRUTEB12ZNYT33W0nMNSu/sdp160/dohBB54D4EDG3QQ6HF2ez9c1jbgBgNQjb9OqK+52OpUXEREflBgZxI6Etom7wXkL4dieU6576NV7CaSZLbY0Jl5+Szcl9G0DJ3+PJvwZQCFbNnxidRyfo/IiIuKjMiZexlLXGGy4aHrv3pOuU7VrBUNK2q7GW3XhAwT427szos8KCO3F3sgLAahZ+6LFaXyPyouIiI+6bGQCr0f/mBbTjuPAh7D/o+N+7q4swHjlZmyGyYcBF3PhpOkWJfVNQeNuBGB4xRLqGhotTuNbfKa86CJ1IiLHs9sMfvytS3nRlQ1A46I7oPiLu0031VD+9LeIcFex0+xL0vf+D5vNsDCt7+l/3lVUEU5vw8nny9+wOo5P8ZnyoovUiYicKCM5kkMjfsZRM4rAmgLcf5+M65/foeWPg+ldt4cyM5wDU59iaEqC1VF9juEXwKGEtntNGdtesTiNb/GZ8iIiIid3+xXjuTXgTyx1jcXmbsG+9z38XfUcdMfx8oD/4fKLsqyO6LN6n/89AEbVrabSWW1xGt9hmOYZXL3Ii1RXVxMREYHT6SQ8PNzqOCIiHqGirpnHluymdsNLpBkHWR90ARMmzeCmCan46YJ0XcftpvR3g4k1j/HJ2L9w4cxbrE7ksTry+9uvmzKJiIiFokIC+N3VI9l7/n3kl9Vxy9BY/FVaup7NRlHiJcQe+Sf2nYtA5aVT6G+uiEgPMigujEuGx6u4dKO4Cd8B2g4dVVQ5LU7jG/S3V0REpAslDL+QUltvQo1Gtq983eo4PkHlRUREpCsZBkWJlwBg3/mmxWF8g8qLiIhIF4s7r+3QUXr9apzVtRan8X4qLyIiIl0sIe0Cyo0oQo1GdqxebHUcr6fyIiIi0tVsNg73vgiA1rx3LA7j/VReREREukFo+lUADK76lKaWVovTeDefKS+6t5GIiHiyfuNnUE8g8UYFOzZ+YnUcr+Yz5UX3NhIREU9mCwjiQHjbrRicn+uso3PhM+VFRETE09mGzQCgT+nH+NjdebqVyouIiEg36X/+t3CZBoPNg+zZnWd1HK+l8iIiItJNAiNiORCYBsDRTTpl+mypvIiIiHSjuuTJAAQf+tjaIF5M5UVERKQbxY+5AoC0xs+pqau3OI13UnkRERHpRvFDz6OKcEKNBnauX2Z1HK+k8iIiItKdbDYKep0HQH3eBxaH8U4qLyIiIt3Mb0jbXaYTjq3SKdNnQeVFRESkm6Vmtc17GWIe4OChfIvTeB+VFxERkW4W3CuBA/6DAChYpxs1dpTKi4iIiAWcCRcA4F+w0uIk3kflRURExAIRw7MB6F+zidZWl8VpvIvKi4iIiAX6ZkyhGT/ijXL27NxsdRyvovIiIiJiAbsjhPzA4QCUbl1qcRrvovIiIiJikfo+EwEIOrzK4iTeReVFRETEIjEjpwEwqP5zmlpaLE7jPXymvOTm5pKWlsb48eOtjiIiInJGkkZMpJ5Aoowadm1ea3Ucr+Ez5SUnJ4e8vDzWr19vdRQREZEzYvg5OBiSDkDlDs17OVM+U15ERES8UVNy2/Vewo6utjiJ91B5ERERsVDsiCkADGjcrnkvZ0jlRURExEKJQ7NowEGkUcfe7Zr6cCZUXkRERCxk+AVwKKjtei/lO5ZbG8ZLqLyIiIhYrD4hEwBHkc44OhMqLyIiIhbrNXQSAP3qtuJyuS1Oc2rOukbeXLaSVoszqryIiIhYLCV9Ei2mnTijgv378qyOc0orFz7KZStn8V7uzy3NofIiIiJiMbsjhEOOQQAUb/3Y4jQnt+9QARMP5eJvuEjr18fSLCovIiIiHqC6d9sV4m2Fnne9F9M0OfDKfUQZtRzxT2XAZdrzIiIi0uOFDL4QgD7VmzFN0+I0x1v16XKm1r4DgP8VfwK7v6V5VF5EREQ8QN+MiwHoxxGKiossTvOVFpeboI9/g90w2R09ldj0aVZHUnkRERHxBIERvTlsb5tLcnjrJxan+cqnHy5irHsbzfiR8p0/Wx0HUHkRERHxGGURIwFoPLjG4iRtWl1uItY+CsDexKsJ6t3P4kRtVF5EREQ8hJHcdrG6iLLN1gb5wqoP32SMexst+NH/mt9YHaedyouIiIiHiB3WNml3QPMuGputvUmjy20Stq7tMNHuxKsJiulraZ6vU3kRERHxEPGDRlNPIGFGA/vzNlmaZcNnHzLGtY0W7PSf9f8szfKfVF5EREQ8hGH3pyBwCADluz+1NEvzZ08CsDt6GsGxqZZm+U8qLyIiIh6krvdoAOxHNliW4UD+ATLrlgPQO/sOy3KcisqLiIiIBwkecB4ACTXbLMtw4INcHEYr+x3DiBt2vmU5TkXlRURExIOkjGq7w3Sq+zAlpaXdPn51XR2jjr4GgGvcj7p9/DOh8iIiIuJBQqISKbbFYjNMCrav6vbxP1/yT2KNSsqNXgya8t1uH/9MqLyIiIh4mJLQ4QA0HOr+eS8hO18BoLDvNRh+jm4f/0x4XHmpqqpi3LhxZGRkMGLECJ566imrI4mIiHQrV9woAAKPbe3WcfMPHiCjaSMAKVN+2K1jd4Sf1QH+U1hYGCtXriQ4OJi6ujpGjBjBNddcQ3R0tNXRREREukXYgEzYC4n1uzBNE8MwumXcgx/9g36Gm32ONAb2TeuWMc+Gx+15sdvtBAcHA9DU1IRpmh53a3AREZGulJw2AYAkSjlacrRbxnS53KQULgKgIe3b3TLm2epweVm5ciUzZ84kMTERwzBYtGjRCevk5uaSmppKYGAgWVlZrFu3rkNjVFVVkZ6eTlJSEnfffTcxMTEdjSkiIuK1AsOjKbIlAHBkx+puGXPrhpUMMAtowp/BU2d3y5hnq8Plpa6ujvT0dHJzc0/684ULFzJ37lweeOABNm3aRHp6OtOnT6f0a6d7fTmf5T+XoqIiACIjI9myZQv5+fm89NJLlJSUnDJPU1MT1dXVxy0iIiLe7ljYMAAaDm3slvFq1r4IwO7Ii3CERnXLmGerw3NeZsyYwYwZM07580cffZTbbruNW2+9FYAnn3ySxYsX8+yzz3LPPfcAsHnz5jMaKy4ujvT0dD755BOuvfbak67z0EMPMW/evI69CBEREQ/XGp8Ozo+6ZdJuc0srgys+AsAx+vouH+9cdeqcl+bmZjZu3Eh2dvZXA9hsZGdns3r1me32KikpoaamBgCn08nKlSsZMmTIKde/9957cTqd7UthYeG5vQgREREPENE/E4Ckhl1dPvdz+/qPiKecWoIYOOHKLh2rM3Tq2UZlZWW4XC7i4uKOezwuLo5du3ad0TYOHTrEj370o/aJurfffjsjR4485foOhwOHwzPPQxcRETlbSWkT4D1I5BhFRw+TmJjcZWPVbmq7ou7+XheQHhDUZeN0Fo87VTozM/OMDyuJiIj4qsCwXhyxJdLHXURR3uouKy+trS4Gli0DIGDk1V0yRmfr1MNGMTEx2O32EybYlpSUEB8f35lDiYiI+LxjYW3XWqk/tKnLxtixcQWJHKMeB4POn9Vl43SmTi0vAQEBjB07lmXLlrU/5na7WbZsGRMmTOjMoU6Qm5tLWloa48eP79JxREREuosrIR2AoLKum7Tr3Nh2yGhv+ET8AkO6bJzO1OHDRrW1tezbt6/9z/n5+WzevJmoqChSUlKYO3cus2fPZty4cWRmZjJ//nzq6urazz7qKjk5OeTk5FBdXU1ERESXjiUiItIdIgdkwq6vJu129pV2TdMktbRth4N95FWduu2u1OHysmHDBqZMmdL+57lz5wIwe/ZsFixYwPXXX8+xY8e4//77KS4uJiMjg/fff/+ESbwiIiJyeklpWbjfMUgwyjlypJA+SSmduv29eZsZzFFaTLvXHDKCsygvkydP/sZTtubMmcOcOXPOOpSIiIiAI6QXhfY+JLsPU7RzdaeXl5KNbzIY2BucTlpIZKduuyt53L2NRERE5Ctl4W2TdrviSrsRhW0XpmtIndbp2+5KPlNeNGFXRER8kTu+bdJucCdP2q0oL2NY83YAUs7zjlOkv+Qz5SUnJ4e8vDzWr19vdRQREZFO02tg25V2kxt3d+qVdvd89ib+hovDtj707jus07bbHXymvIiIiPiipLTzcJsGcVRwpPBQp23X3LsEgOK4SZ22ze6i8iIiIuLBAoLDOeyXBMDRnWd2n8Bv0tLayuDqtm2Fp1/RKdvsTiovIiIiHu7LSbuNBRs6ZXt7N68iGie1BDFgbPY3P8HD+Ex50YRdERHxWfEZAASXb++UzVVsfR+A/SFjsPt7382Nfaa8aMKuiIj4qshOnrQbefQTAFr7TfmGNT2Tz5QXERERX5U0LAuXaRBLJUcKDpzTtpyVFQxuzmvb7jjvm+8CKi8iIiIeLyA4jMN+bVfXLdq55py2tXf9+wQYLoqMeOJSvesU6S+pvIiIiHiB8i8m7TYXnNuVdpt3t92IsSh6wjlnsorKi4iIiDdI+OJKu+XbznoTpmmSVPEZAI6h3neW0ZdUXkRERLxAr0HnAZDStBu3y31W2ziSv5sUs4hW08aAzMs6M1638pnyolOlRUTElyUNG0+raSMGJ4cLz27SbuGGxQDscwwlODyqM+N1K58pLzpVWkREfJl/YGj7pN3is7zSrv+h5QA4Ey7srFiW8JnyIiIi4uvKI4YD0HwWV9ptaWlmUG3bZN/o9Bmdmqu7qbyIiIh4CSMxA4CQih0dfu7ezz8hwqijmhD6j7qgk5N1L5UXERERL9FrYBYAfRt34ergpN2qbV/cEiB0HDY//07P1p1UXkRERLxE8rDxtJh2oowaDu3r2N6XXsWfAuDy0lsCfJ3Ki4iIiJfwcwRzMGAQACV5K8/4ec7KcgY17wIgZbz3niL9JZUXERERL+KMGQ2AWbDujJ+zf+1i/Aw3hbY+xKYM6apo3cZnyouu8yIiIj2Bo1/bxep6O7ec8XOa97TdEuCoF98S4Ot8przoOi8iItITJI2aDEB/10Gqqiq+cX3TNEmubLsuTODQaV0Zrdv4THkRERHpCXrFp1JixGA3TPK3fPKN6x8+kEcfs4Rm087ArEu7IWHXU3kRERHxMkVhIwGo2//ZN657pP2WAMMJDo3syljdRuVFRETEy7Qmts3vDCnd9I3rOg59DEB1H+++JcDXqbyIiIh4meihbVfI7deQd9qL1bU0NzGo7nMAYkb5xiEjUHkRERHxOn2Hn0ej6U+kUcuBnafe+7Lv8+WEGg1UEkb/URO7MWHXUnkRERHxMnZ/BweC2ua9FG9865TrObd9AMCBsPHY7PZuydYdVF5ERES8UMOAtsNA0YVLT7lOdHHbVXhd/b3/lgBf5zPlRRepExGRnqTfxOsAGNqyk+IjBSf8/NCujQxq3UuraaPfhFndnK5r+Ux50UXqRESkJ4lK7M8+v0HYDJMDn/37hJ8fXf4MANtCzqN3fEp3x+tSPlNeREREepqK5LYr5gbtf/+4x1uamxhU/E7bH0Z/r7tjdTmVFxERES+VOOHbAKQ1bML5tVsF7Fjxb6JxUkYkIyZda1W8LqPyIiIi4qWSBmVw2JaIw2hh50f//OoHn78IwN74y/EPcFiUruuovIiIiHgrw6AoZSYAw7f8noO7NvH5+/9gZF3bjRgTJv/QynRdRuVFRETEi4258bfsDBhBmNFAyMvXMGr1z7EbJmsjLiV16Bir43UJlRcREREv5hcQSPwPX6GY3vSmErthsqbXVYyd86LV0bqMyouIiIiX6xXbh9pvvcBuv8F8lvJjsm5fgJ+/v9Wxuoyf1QFERETk3A0cOQFG9oxrnWnPi4iIiHgVlRcRERHxKj5TXnRvIxERkZ7BME3TtDpEZ6quriYiIgKn00l4eLjVcUREROQMdOT3t8/seREREZGeQeVFREREvIrKi4iIiHgVlRcRERHxKiovIiIi4lVUXkRERMSrqLyIiIiIV1F5EREREa+i8iIiIiJeReVFREREvIqf1QE625d3O6iurrY4iYiIiJypL39vn8ldi3yuvNTU1ACQnJxscRIRERHpqJqaGiIiIk67js/dmNHtdlNUVERYWBiGYXTZOOPHj2f9+vVdtn0rxqyuriY5OZnCwkLd1FI8khX/3cmJ9DmcnK++L931ukzTpKamhsTERGy2089q8bk9LzabjaSkpC4fx263d/sv+O4aMzw8XOVFPJIV/93JifQ5nJyvvi/d+bq+aY/LlzRh9yzl5OT0iDFFPIn+G/AM+hxOzlffF098XT532EjOXnV1NRERETidTp/89iAiIr5Be16kncPh4IEHHsDhcFgdRURE5JS050VERES8iva8iIiIiFdReRERERGvovIiIiIiXkXlRUQEKCwsZPLkyaSlpTFq1CheffVVqyP1SPocTk7vy/E0YVdEBDh69CglJSVkZGRQXFzM2LFj2bNnDyEhIVZH61H0OZyc3pfjac+LnBV9CxBfk5CQQEZGBgDx8fHExMRQUVFhbageSJ/Dyel9OZ7Ki5wVPz8/5s+fT15eHkuWLOHOO++krq7O6lhigYceeojx48cTFhZGbGwss2bNYvfu3Z06xsqVK5k5cyaJiYkYhsGiRYtOul5ubi6pqakEBgaSlZXFunXrzmq8jRs34nK5vOoGr3/9618ZNWpU++09JkyYwHvvvdepY3j75/Dwww9jGAZ33nlnp2zvS97+vngjlRc5K/oWIF9asWIFOTk5rFmzhqVLl9LS0sIll1xyyjK7atUqWlpaTng8Ly+PkpKSkz6nrq6O9PR0cnNzT5lj4cKFzJ07lwceeIBNmzaRnp7O9OnTKS0tbV8nIyODESNGnLAUFRW1r1NRUcHNN9/M3//+9zN9CzxCUlISDz/8MBs3bmTDhg1cfPHFXHXVVezYseOk6/e0z2H9+vX87W9/Y9SoUaddr6e9L17LFJ+0YsUK84orrjATEhJMwHzjjTdOWOeJJ54w+/btazocDjMzM9Ncu3btWY21YcMGc/jw4eeYWHxFaWmpCZgrVqw44Wcul8tMT083r732WrO1tbX98V27dplxcXHm//zP/3zj9k/19zkzM9PMyck5bqzExETzoYceOuPsjY2N5oUXXmg+//zzZ/wcT9arVy/z6aefPuHxnvY51NTUmIMGDTKXLl1qTpo0ybzjjjtOul5Pe1+8mfa8+Khv+iagbwHSVZxOJwBRUVEn/Mxms/Huu+/y+eefc/PNN+N2u9m/fz8XX3wxs2bN4pe//OVZjdnc3MzGjRvJzs4+bqzs7GxWr159RtswTZNbbrmFiy++mJtuuumscngKl8vFyy+/TF1dHRMmTDjh5z3tc8jJyeHyyy8/LtfJ9LT3xatZXJ6kG3CSbwL6FiBdweVymZdffrk5ceLE06536NAhMyUlxbz++uvNlJQU8+abbzbdbvcZjXGyv89HjhwxAfOzzz477vG7777bzMzMPKPtfvLJJ6ZhGGZ6enr7snXr1jN6rqfYunWrGRISYtrtdjMiIsJcvHjxadfvCZ/Dv/71L3PEiBFmQ0ODaZrmafe8fKknvC/ezs/C3iQW+fJbwL333tv+mL4FSGfIyclh+/btfPrpp6ddLyUlhRdeeIFJkybRv39/nnnmGQzD6KaUJ3fBBRfgdrstzXCuhgwZwubNm3E6nfz73/9m9uzZrFixgrS0tJOu7+ufQ2FhIXfccQdLly4lMDDwjJ/n6++LL9Bhox6orKwMl8tFXFzccY/HxcVRXFx8RttYtWoVCxcuZNGiRWRkZJCRkcG2bdu6Iq54iTlz5vDOO+/w8ccfk5SUdNp1S0pK+NGPfsTMmTOpr6/n5z//+TmNHRMTg91uP2FCZUlJCfHx8ee0bW8SEBDAwIEDGTt2LA899BDp6en85S9/OeX6vv45bNy4kdLSUsaMGYOfnx9+fn6sWLGCxx9/HD8/P1wu10mf5+vviy/Qnhc5K/oWIF8yTZPbb7+dN954g+XLl9OvX7/Trl9WVsbUqVMZNmwYr776Knv27GHy5Mk4HA7+9Kc/nVWGgIAAxo4dy7Jly5g1axYAbrebZcuWMWfOnLPapi9wu900NTWd9Gc94XOYOnXqCV+qbr31VoYOHcqvfvUr7Hb7Cc/pCe+LL1B56YH0LUA6U05ODi+99BJvvvkmYWFh7XvvIiIiCAoKOm5dt9vNjBkz6Nu3LwsXLsTPz4+0tDSWLl3KxRdfTJ8+fU76Lbe2tpZ9+/a1/zk/P5/NmzcTFRVFSkoKAHPnzmX27NmMGzeOzMxM5s+fT11dHbfeemsXvnrPce+99zJjxgxSUlKoqanhpZdeYvny5XzwwQcnrNtTPoewsDBGjBhx3GMhISFER0ef8Dj0nPfFJ1g96Ua6HqeYsDtnzpz2P7tcLrNPnz4dmrArYpptf79OtvzjH/846fpLlixpnzz5dZs2bTILCwtP+pyPP/74pGPMnj37uPX+93//10xJSTEDAgLMzMxMc82aNef68rzG97//fbNv375mQECA2bt3b3Pq1KnmkiVLTrl+T/0cvmnCbk99X7yN7m3ko77+TWD06NE8+uijTJkypf2bwMKFC5k9ezZ/+9vf2r8FvPLKK+zateuEuTAiIiKeROXFRy1fvpwpU6ac8Pjs2bNZsGABAE888QSPPPIIxcXFZGRk8Pjjj5OVldXNSUVERDpG5UVERES8ik6VFhEREa+i8iIiIiJeReVFREREvIrKi4iIiHgVlRcRERHxKiovIiIi4lVUXkRERMSrqLyIiIiIV1F5EREREa+i8iIiIiJeReVFREREvIrKi4iIiHiV/w/DzVBmiaVfxwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(qq, model_nondisp(qq), label='nondisp')\n",
    "\n",
    "# note how we provide theta and wavelength to the ReflectModelTL object.\n",
    "plt.plot(qq, model_disp(np.c_[theta, wavelength]), label='disp')\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25f47802-b6c0-4b41-8748-74f17d940c36",
   "metadata": {},
   "source": [
    "To reassure ourselves let's loko at the slab representation of the the dispersive `Structure` at two different wavelengths. We can see that the real and imaginary components of the SLD (second column) do change, just not by much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "274ac51a-006e-404a-8d16-f20bc7575051",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  0.           0.           0.           0.           0.        ]\n",
      " [300.           6.54468625   0.73189622   5.           0.        ]\n",
      " [  0.           2.07         0.           3.           0.        ]]\n",
      "\n",
      "[[  0.           0.           0.           0.           0.        ]\n",
      " [300.           6.47192866   0.64839038   5.           0.        ]\n",
      " [  0.           2.07         0.           3.           0.        ]]\n"
     ]
    }
   ],
   "source": [
    "print(s_disp.slabs(wavelength=2.8))\n",
    "print()\n",
    "print(s_disp.slabs(wavelength=18.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0da8e77-1d80-4cd2-b172-43cc3f47bc7e",
   "metadata": {},
   "source": [
    "For comparison here is the non-dispersive system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "76fb6a62-3892-45a6-b3b4-170178edadf1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  0.     0.     0.     0.     0.  ]\n",
      " [300.     6.5    0.68   5.     0.  ]\n",
      " [  0.     2.07   0.     3.     0.  ]]\n",
      "\n",
      "[[  0.     0.     0.     0.     0.  ]\n",
      " [300.     6.5    0.68   5.     0.  ]\n",
      " [  0.     2.07   0.     3.     0.  ]]\n"
     ]
    }
   ],
   "source": [
    "print(s_nondisp.slabs(wavelength=2.8))\n",
    "print()\n",
    "print(s_nondisp.slabs(wavelength=18.0))"
   ]
  }
 ],
 "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.12.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}