{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f4d55ed",
   "metadata": {},
   "source": [
    "# Modeling BLS signal from single magnetic layer using Green's functions\n",
    "This example shows how SpinWaveToolkit (SWT) can be used to calculate the expected BLS signal of a single magnetic layer using the Green's function formalism and laser electric fields in real space.\n",
    "\n",
    "Source publication: Wojewoda et al., *Physical Review B* **110**, 224428 (2024). DOI: [10.1103/PhysRevB.110.224428](https://doi.org/10.1103/PhysRevB.110.224428)\n",
    "\n",
    "We start by impoting the modules we need and by defining the parameters of the system and model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b1ba96bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import numpy as np\n",
    "import SpinWaveToolkit as SWT\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8805b570",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define parameters\n",
    "Bext = 50e-3         # External field [T]\n",
    "theta = np.pi / 2    # Out-of-plane angle (fixed)\n",
    "d_layer = 30e-9      # Magnetic layer thickness [m]\n",
    "material = SWT.NiFe  # Material (from built-ins)\n",
    "Nf_common = 101       # Number of frequency points for Bloch function\n",
    "\n",
    "# Define the Kx,Ky grid limits and resolution.\n",
    "Nk = 100             # resolution in Kx and Ky\n",
    "k_max = 15e6         # maximum k (as in original kxi range)\n",
    "NA = 0.75            # Numerical Aperture (NA) for the lens\n",
    "\n",
    "# Create a regular Kx,Ky grid\n",
    "kx_grid = np.linspace(-k_max, k_max, Nk)\n",
    "ky_grid = np.linspace(-k_max, k_max, Nk)\n",
    "kx_grid[np.abs(kx_grid) < 1e-6] = 1e-6  # avoid division by zero at k=0\n",
    "kx_grid[np.abs(kx_grid) < 1e-6] = 1e-6\n",
    "KX, KY = np.meshgrid(kx_grid, ky_grid, indexing='ij')\n",
    "K = np.sqrt(KX**2 + KY**2)\n",
    "PHI = np.arctan2(KY, KX)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73175d39",
   "metadata": {},
   "source": [
    "Now we calculate the electric field incident on the magnetic layer using the `ObjectiveLens` class, get the Bloch functions in 3D and compute the BLS signal using the `get_signal_GF_focal()` function.\n",
    "\n",
    "In the current state, the `get_signal_GF_focal()` function implements only linear magneto-optical response, whereas the reciprocity theorem-related functions (`get_signal_RT_focal()` and `get_signal_RT_pupil()`) can use also quadratic magneto-optical effect.\n",
    "\n",
    "Since our input field is linearly polarized along *x* and only the linear MO response is assumed, the maximum BLS signal is expected to be at 90° rotated polarization (along *y*). We can define more output polarizer types/orientations, but here we show only the most common case for BLS on spin waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8770378b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing focal field...\n",
      "Preparing Bloch functions in 3D (f,kx,ky) using Kalinikos-Slavin model...\n",
      "Computing BLS signal from Kalinikos-Slavin model...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\228634\\AppData\\Local\\Temp\\ipykernel_23984\\293609646.py:43: UserWarning: `get_signal_GF_focal` is an experimental function and may be subject to change. Please verify results carefully.\n",
      "  sigma = SWT.bls.get_signal_GF_focal(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done\n"
     ]
    }
   ],
   "source": [
    "print(\"Preparing focal field...\")\n",
    "objective = SWT.bls.ObjectiveLens(NA=NA, wavelength=532e-9, f0=10, f=1e-3)\n",
    "x, y, Ex, Ey, Ez = objective.getFocalField(z=0, rho_max=10e-6, N=400)\n",
    "E = [Ex, Ey, Ez]\n",
    "Exy = [x, y]\n",
    "\n",
    "# Preallocate the output array for the Bloch function amplitudes.\n",
    "# The result will be complex amplitudes defined on a common frequency axis.\n",
    "# Shape: (Nf_common, Nk, Nk)\n",
    "Bloch2D = np.zeros((Nf_common, Nk, Nk), dtype=complex)\n",
    "\n",
    "# We will also store the common frequency axis.\n",
    "w_common = np.linspace(4, 27, Nf_common)*2e9*np.pi  # (rad/s)\n",
    "\n",
    "print(\"Preparing Bloch functions in 3D (f,kx,ky) using Kalinikos-Slavin model...\")\n",
    "# Create a SingleLayer model for the current kxi and phi.\n",
    "# Note: We pass kxi and phi as arbitrary values for now.\n",
    "model = SWT.SingleLayer(Bext=Bext, kxi=(0,), theta=theta,\n",
    "                        phi=0, d=d_layer, material=material)\n",
    "# Loop over all grid points in the Kx,Ky plane.\n",
    "for i in range(Nk):\n",
    "    for j in range(Nk):\n",
    "        # Set correct kxi and phi for the current grid point.\n",
    "        model.kxi = np.array((K[i, j],))  # must be an array\n",
    "        model.phi = PHI[i, j]\n",
    "        # Compute the Bloch functions for n=0,1; higher n are outside `w_common`\n",
    "        # The returned w has shape (Nf_common,) and bf has shape (Nf_common, 1),\n",
    "        # since we are passing kxi as length-1 array.\n",
    "        # We pass also parameters for Bose-Einstein distribution (optional).\n",
    "        w0, bf0 = model.GetBlochFunction(n=0, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "        w1, bf1 = model.GetBlochFunction(n=1, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "        # Interpolate the Bloch functions to the common frequency axis.\n",
    "        bf0_interp = np.interp(w_common, w0, bf0[:, 0], left=0, right=0)\n",
    "        bf1_interp = np.interp(w_common, w1, bf1[:, 0], left=0, right=0)\n",
    "\n",
    "        # Sum Bloch functions for n=0,1\n",
    "        Bloch2D[:, i, j] = bf0_interp + bf1_interp\n",
    "\n",
    "Bloch3 = [Bloch2D, np.zeros((Nf_common, Nk, Nk), dtype=complex), Bloch2D*1j]\n",
    "\n",
    "print(\"Computing BLS signal from Kalinikos-Slavin model...\")\n",
    "# Compute the BLS signal using the Bloch functions and the electric field.\n",
    "sigma = SWT.bls.get_signal_GF_focal(\n",
    "    # The Bloch functions and related axis vectors\n",
    "    SweepBloch=w_common, KxKyBloch=[kx_grid, ky_grid], Bloch=np.array(Bloch3),\n",
    "    Exy=Exy, E=E,  # the incoming electric field from the objective\n",
    "    Nq=50,  # number of points to use in the k-space integration\n",
    "    # dielectric functions of the multilayer system: air, NiFe, substrate\n",
    "    DF=[1, -8.1653 + 1j*15.348, 17.237 + 1j*0.43004],\n",
    "    PM=[1,1,1],  # same for permeabilities (response of NiFe is included in the Bloch functions)\n",
    "    d=[d_layer], # thicknesses of the inner layers, here just NiFe\n",
    "    NA=NA, # numerical aperture of the collecting objective (here the same as focusing)\n",
    "    source_layer_index=1,  # layer with the magnetic film\n",
    "    output_layer_index=0,  # layer in which we observe the scattered light (air)\n",
    "    # other parameters of the optical setup\n",
    "    wavelength=532e-9, collectionSpot=0.5e-6, focalLength=1e-3,\n",
    "    output_analyzer=\"linear\", output_analyzer_angle_deg=90,\n",
    ")\n",
    "print(\"Done\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c53533fe",
   "metadata": {},
   "source": [
    "Finally, we plot the resulting BLS signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c35a5e8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjUsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvWftoOwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPixJREFUeJztnQd8FOXWh096BQKk0TuE3qv0LqAUCyBNRIoiIIgIXhFU/BAuSBGQi15Br9IRREUQ6b33EgKETkggEEJ6me93TjLLbrJJdsOG7M7+Hx12dubd2Xd2N3PmdAdFURQCAAAATMTR1IEAAAAAA8EBAADALCA4AAAAmAUEBwAAALOA4AAAAGAWEBwAAADMAoIDAACAWUBwAAAAMAsIDgAAAGYBwQGAjbFz505ycHCQR2tg6tSpMh9gP0BwALtm2bJlctHTX/z9/alNmzb0119/ZRrP+997771sj5mamko//fQTNW7cmIoUKUIFChSgypUr08CBA+ngwYN5eDYAPB+cn9P7AGDVfP7551SuXDni0m337t0TgdKlSxf6/fffqVu3bmYda/To0bRw4ULq3r079evXj5ydnSk4OFgEUfny5alJkybPNNeWLVtSXFwcubq6PtNxAMgtEBwAENGLL75IDRo00D0fMmQIBQQE0IoVK8wSHCx0Fi1aREOHDqUlS5YY7Js7dy5FREQ881wdHR3J3d39mY8DQG6BqQoAI/j4+JCHh4doC+YQGhoqWssLL7yQaZ9qBsuJlStXUv369cXEVbBgQapZsybNmzcvRx8Hazms0fC8GzVqRHv27KHWrVvLkvG1q1evpi+//JJKliwpQqhdu3Z0+fJlg+Px61977TUqXbo0ubm5UalSpWjs2LGi7QD7BhoHAEQUFRVF9+/fl4t+eHg4ffPNN/TkyRPq37+/WccpU6aMPK5Zs0Yuup6enma9fuvWrdS3b1+5kM+YMUO2Xbhwgfbt20djxozJ8nXffvut+F5atGghF/dr165Rjx49qHDhwiIcMvLVV1+J5jJ+/Hg595kzZ4pZ7dChQ7oxfA6xsbH0zjvvUNGiRenw4cPyudy6dUv2AfsFggMAImrfvr3Bc77D/uGHH6hDhw5mHadYsWLiBGfnOF+w+W6ftY+uXbtSUFBQjq//888/RcvYsmULOTk5mfSeiYmJNHnyZGrYsCFt375dpyXVqlWL3nzzTaOCIz4+nk6ePKnzk7CAYcF09uxZqlGjhmxjwcXai8qwYcOoYsWK9PHHH9ONGzdEEwH2CUxVAKSbefhun5eff/5Zoqrefvtt+vXXX80+1tKlS2nBggXibF+/fr3c1VetWlW0iNu3b+doIouJiZF5mMrRo0fpwYMH4lfRN62xBsECwRiDBw82cK6zpsJcvXpVt01faPCcWCNr1qyZaGUnTpwweX5Ae0BwAEAkPgHWOnjhCy7f+VerVk3MP3xHbw5sAho5ciQdO3ZMLra//fabON9ZG+jTp0+2r3333XcldJfHs6bw1ltv0ebNm7N9zfXr1+WRtQF9WIiULVvW6GsyaguqgHn48KFuG2sVrLFwSLG3tzf5+flRq1atZB+bt4D9AsEBQBYXf9Y67t69SyEhIbk+DvsGXn75Zdq0aZNcdPfu3au70BuDnedsQtq4caO8bseOHSJEBg0aRJYkKzOY2kk6JSVFzHQsQD/66CPasGGDaEEcpqzmqgD7BYIDgCxITk6WR3aSWwI13JeFUXawCemll16SsN4rV67Q8OHDxWeSMeopo0M+436ePzvJc8OZM2fo0qVLNHv2bBEcnJPC2ljx4sVzdTygLSA4ADBCUlIS/f3333IRZ/+EqYSFhdH58+czbWdz17Zt20STyWhS0od9FfrweHZyMwkJCVkKJNZsvvvuO52wY3755RcD01NuNBJVA1HX9cOCgf2CqCoAiCSr++LFi7LO4bjLly8XE9XEiRMlyimjM3ratGmZjsERVJwTwf6Stm3bijM8MDBQjseJhKdOnaL333+ffH19s5wHO+QjIyPl9ezjYLMWh8DWqVMnSwHGwo3rRY0aNUpe9/rrr4umwWalChUq5KqOFEeA8WvZsc8Off4M1q1bl2tBBDSGAoAds3TpUr6lNljc3d2VOnXqKN9++62SmppqMD7jWP3liy++UB4/fqzMmzdP6dSpk1KyZEnFxcVFKVCggNK0aVPlu+++y3S8jKxdu1bp2LGj4u/vr7i6uiqlS5dWhg8frty9e1c3ZseOHfJ+/KjP/PnzlTJlyihubm5Ko0aNlH379in169dXOnfunOm1a9asMXhtaGiobOfPQ+X8+fNK+/btFW9vb8XX11cZOnSocurUqUzjpkyZItuA/eDA/+S38AIAWB52YHMkVK9evcSMBYClgI8DAA3ACX0Z7wHZoc5mL/2SIwBYAmgcAGgArkHFpUa4zAk7yo8fP07//e9/xS/C+SSopAssCZzjAGgATvTjIoTz588XLYOT9rj0CdekgtAAlgYaBwAAALOAjwMAAIBZQHAAAAAwC/g4niHU8c6dO9JsJzcJVgAAYG2w5yI6OlpKy3DVgqyA4MglLDTYGQkAAFrj5s2bRvu4qEBw5BLWNNQPOGNJCgAAsEUeP34sN8Tq9S0rIDhyiWqeYqEBwQEA0BI5md/hHAcAAGAWEBwAAADMAoIDAACAWUBwAAAAMAsIDgAAyGOSU7TVox2CAwAA8pDNZ8OoxtQttOlM9r3mbQkIDgAAyEMOh0ZSfFIqHbpq2E/eloHgAACAPCQuKVkeYxJTSCtAcAAAQB4Sly4wYhLSBIgWgOAAAIA8JC4pXXBA4wAAAGAKcUlpEVWx0DgAAACYQny6pvEEggMAAIA5pqpYmKoAAACYQmxissGjFoDgAACAPCQ+3ccBUxUAAACzTFUsQFJSFdIC+S44Fi5cSGXLliV3d3dq3LgxHT58ONvxa9asoaCgIBlfs2ZN2rRpk8H+X3/9lTp27EhFixaVZiQnT57MdIzWrVvLPv1lxIgRFj83AACI0/NtaMVcla+CY9WqVTRu3DiaMmUKHT9+nGrXrk2dOnWi8PBwo+P3799Pffv2pSFDhtCJEyeoR48espw9e1Y3JiYmhpo3b04zZszI9r2HDh1Kd+/e1S0zZ860+PkBAOwbRVF0GgcTk6ANB7mDwmeWT7CG0bBhQ1qwYIE8T01NlX63o0aNookTJ2Ya37t3bxEMf/zxh25bkyZNqE6dOrR48WKDsdeuXaNy5cqJgOH9GTUO3jZ37txn6s1bqFAhioqKQutYAIBR4pNSKGjyZt3zbR+0ogp+3mStmHpdyzeNIzExkY4dO0bt27d/OhlHR3l+4MABo6/h7frjGdZQshqfHb/88gv5+vpSjRo1aNKkSRQbG5vt+ISEBPlQ9RcAAMiOjCG4sRrROJzz643v379PKSkpFBAQYLCdn1+8eNHoa8LCwoyO5+3m8MYbb1CZMmWoePHidPr0afroo48oODhY/CNZMX36dPrss8/Meh8AgH0Tp2em0lJkVb4Jjvxk2LBhunV2sBcrVozatWtHV65coQoVKhh9DWsl7I9RYY2DzWoAAGCKY1xLzvF8ExxsJnJycqJ79+4ZbOfngYGBRl/D280Zb46vhbl8+XKWgsPNzU0WAAAwx8ehRY0j33wcrq6uVL9+fdq2bZtuGzvH+XnTpk2Nvoa3649ntm7dmuV4U1FDdlnzAACAvDJVxWqk7Ei+mqrY9DNo0CBq0KABNWrUSKKcOGpq8ODBsn/gwIFUokQJ8S8wY8aMoVatWtHs2bOpa9eutHLlSjp69CgtWbJEd8zIyEi6ceMG3blzR56z74JhrYQXNkctX76cunTpIrke7OMYO3YstWzZkmrVqpUvnwMAwD5MVTEa0TjyVXBweG1ERAR9+umn4uDmENnNmzfrHOAsADjSSqVZs2Zy0f/kk0/o448/pkqVKtGGDRskMkpl48aNOsHD9OnTRx45V2Tq1Kmi6fzzzz86IcV+ildeeUWOCQAAliQ2k+DQhsaRr3kctgzyOAAAObHhxG16f9XT6hXDW5anSV2qkrVi9XkcAABgbz6OGI1EVUFwAADAc/NxpJAWgOAAAIDnpXEk2KFz/NGjR7R+/Xras2cPXb9+Xcp0+Pn5Ud26daX0BzuvAQAAGOZx+Hi60KPYJM2E45qkcXBo69tvvy15DtOmTaO4uDiJgOJs65IlS9KOHTuoQ4cOVK1aNal4CwAAgHSCoqiXq6YSAE3SOFij4HwLLkrIwsEYLEw4NJbDXG/evEnjx4+39FwBAMAmTVVFvd3oSkSMfZUcOX/+vCTLZYeHh4f0yuDlwYMHlpofAADYLPHpGoeft5v9OcdzEhrPOh4AALStcbjKI8JxM8Bd9DjTGwAAgKHg8E3XOLTSj8NigqNt27bScQ8AAIBhHocqOBJTUikxOZVsHYvVqvrpp59y7KIHAAD2bKpi2EHu6vz0uV0LDu4dDgAAILPGUcDdmVydHUXbiElMIR9PsmmQOQ4AAHmscXi4OJGXq5NmssfN1ji4zLmDg0OW+7mPOAAAANJljnu4OpGnqzM9jE2yT8HBJUf0SUpKohMnTtCPP/5In332mSXnBgAAmjBVebo4k7db2uVWC2VHzBYc3bt3z7Tt1VdfperVq0u5kSFDhlhqbgAAYLMoiqIzVbm7OpKnm5Nmyo5YzMfRpEmTTP3AAQDAXklITqXU9DZ5aT4OVeOA4NDVqZo/f770BwcAAEA6/wbjzoIjXePQQtkRs01VhQsXNnCOszoWHR1Nnp6e9PPPP1t6fgAAYJPEpQsOFycHcnFy1Gkcdukc5+q3GaOsuCdH48aNRagAAAAgnWOctQ1G9XFwHofdCQ4urw4AAMD0HA7GS42q0oDGgQRAAADIQx+HZ3rin5dqqoJz/ClVq1YlJ6e0DwgAAOyd2AymKlXjsEvneFZMnz6doqKiLHU4AADQhI/DQ6dx2HHJkazo0aOHpQ4FAACa83F4qhoHTFUAAACyrVOVLji806Oq7LLkCLN27VpavXq1dPxLTEw02Hf8+HFLzQ0AAGw/HNc1XeNId47bZckRzhAfPHgwBQQESHHDRo0aSY/xq1ev0osvvpg3swQAABsjLimt05+n6hxXS45owDlutuBYtGgRLVmyhL755htydXWlCRMm0NatW2n06NFwjgMAQEYfh+oc1yUA2qHGweapZs2aybqHh4eUG2EGDBhAK1assPwMAQDABolLFxAZEwA5qopLNdmV4AgMDKTIyEhZL126NB08eFDWQ0NDbf7DAAAASxGnllTPIDi4Yi5XzrUrwdG2bVvauHGjrLOvY+zYsdShQwfq3bs39ezZMy/mCAAANkdcYqqBqUrVPLSQy2F2VBX7N1JT0z6QkSNHimN8//799PLLL9Pw4cPzYo4AAGDz4bhOjg6yzpoIZ48X9Sb76jnOi0qfPn1kAQAAkHUCoGquEsFh4w5yR1Md4uZw+/bt3M4HAAA0WXKEedrMyQ4ER8OGDcUMdeTIkSzHcCjud999RzVq1KB169ZZco4AAGBzxBrRONQkQFvvyWGSqer8+fP05ZdfihPc3d2d6tevT8WLF5f1hw8fyv5z585RvXr1aObMmdSlS5e8nzkAAFgx8UY0Dl3ZEXvQONgB/vXXX9Pdu3dpwYIFVKlSJbp//z6FhITI/n79+tGxY8fowIEDEBoAAECZw3G1VHbELOc4J/y9+uqrsgAAADDPOa6ux9tbHgcAAIDcmarUdXWfrQLBAQAAeahxeOoJDtVspe6zVSA4AADAwiQmp1Iy1xbJ4ONQTVUQHAAAAAzQFwwGPg5XR4McD7sRHDExMRadwMKFC6ls2bIS2tu4cWM6fPhwtuPXrFlDQUFBMr5mzZq0adMmg/2//vordezYUSLBHBwc6OTJk5mOER8fryuX4u3tTa+88grdu3fPoucFALBf4tMFB5cZcXFyyOwctzeNgxs4vfXWW7R3795nfvNVq1bRuHHjaMqUKdI5sHbt2tSpUycKDw83Op5rYvXt25eGDBkiTaS4zzkvZ8+eNRBszZs3pxkzZmT5vlyY8ffffxchtGvXLrpz5w716tXrmc8HAAAMssZdnOQGVms+Di6Fbhbr169Xunfvrri4uCiVKlVSpk+frty+fVvJDY0aNVJGjhype56SkqIUL15cjmmM119/XenatavBtsaNGyvDhw/PNDY0NJQNjMqJEycMtj969EjmvmbNGt22CxcuyNgDBw6YPPeoqCh5DT8CAIA+5+9EKWU++kOp/8VWg+0/H7wm24f+eESxRky9rpmtcfAd/oYNG6Qe1YgRI2j58uVUpkwZ6tatm5iJkpNNS2zhXuWcNNi+fXvdNi6eyM85kdAYvF1/PMMaSlbjjcHvmZSUZHAcNn1xb5HsjpOQkECPHz82WAAAwNSIKsbuneN+fn5iZjp9+rRklf/zzz+SGMilSD799FOKjY3N9vWceZ6SkiKmL334eVhYmNHX8HZzxmd1DG556+PjY9Zxpk+fToUKFdItpUqVMvk9AQD2a6rSx259HCrsTOa6VNWqVaOJEyeK0Ni2bRvNnj1bNA/WTLTEpEmTpJCjuty8eTO/pwQAsHLB4Z5B41Cf27rGYXY/DhYKS5cupS1btojQePfdd6l///4Gd/Dck7xq1arZHsfX15ecnJwyRTPxc25Pawzebs74rI7BZrJHjx4ZzDmn47i5uckCAACmlxtxNG6qsrdwXG4Xy+aoffv2Sajre++9l8nsw/v/9a9/ZXscNhdxlV3WUlS4syA/b9q0qdHX8Hb98czWrVuzHG8Mfk8XFxeD4wQHB0vPEXOOAwAAuTdVpdqXxsEVcj09PXMshsghtjnBPpJBgwZRgwYNqFGjRjR37lwJp2XhxAwcOJBKlCgh/gVmzJgx1KpVKzGHde3alVauXElHjx6VdrYqkZGRIgQ4xFYVCgxrE7ywf4LDefm9ixQpQgULFqRRo0aJ0GjSpIm5HwcAAGTjHDe8xHraq6mqQIECIjz8/f0Ntj948EC2scPbVHr37k0RERHiTGfHdJ06dWjz5s06BzgLAP02tWwC4yiuTz75hD7++GMp784RXtw8SmXjxo06wcOobW1ZkE2dOlXW58yZI8flxD+OluLIrEWLFpn7UQAAgFFiVR9HBo3DXSOmKgeOyTXnBXzB5Yt8RsHBd/gVKlSguLg4sgc4HJe1F3aUs9YCAAAqX2+9RPO3hdCAJmXoix5Pb2zvP0mgBtP+kfXQ6V0MkgNt6bpmssYxf/58eeQT/f7776VUhwprGbt375Z8CAAAsHfiEpMzlVSX53oaSEJyaiaNxFYwWXCweYdhBWXx4sUSEaXv6OZ6U7wdAADsnTgjTZwYfUHB5izNC47Q0FB5bNOmjYTkFi5cOC/nBQAANu/j8MigcXDRQ1dnRym7bssOcrOd4zt27MibmQAAgEaIz6LkiKqFiOBI1Ljg4NDVL774gry8vGQ9O7j8CAAA2DOxWeRxqNui4pJsuuyISYKDS5hzYUB1PSusLUIAAACsyVSlv03zpip98xRMVQAAkHtTlRZyORwtEffLSXgXL160zIwAAECjCYD69atsWeMwW3C8/vrrtGDBAlnnZD8uF8LbuI3runXr8mKOAABgU8QlGi85om+qircnwcGJfi1atJD19evXS14HV5rlBMFp06blxRwBAEATeRxaqZBrtuDgVHQuDshwXSmu98RFD7noYEhISF7MEQAAbFTjcMq0Twt9x80WHNz5jluschVbFhwdO3aU7Q8fPiR3d/e8mCMAANgMqanKU40jizwOWxccZicAvv/++9SvXz+pVcW9xlu3bq0zYbGfAwAA7JmE5Ke9NoyaqlQfR6IdCQ7u+Me9M7h1aocOHXRlz8uXLw8fBwDA7olNL3CYo4/DnjQOhiOpeNGHfRwAAGDvxKZrEm7OjuTo6KBJH4fZgoNLqC9btkxar4aHh0u7V322b99uyfkBAIBmkv8MMscTU+1HcHD7VhYcrGFw5z2UGQEAANPqVOlvt+U8DrMFB/f5Xr16NXXp0iVvZgQAADZMXDYRVbLdHsNxuWlTxYoV82Y2AACgkRwOjywEh7vOVGVHguODDz6gefPmScY4AAAAQ1RNwtPFWbMah9mmqr1790qF3L/++ouqV69OLi4uBvu5OyAAANgrsTloHKrT3K58HD4+PtSzZ8+8mQ0AAGi4TpXdhuMuXbo0b2YC8pSklFR6a9kRqhJQgD7pVi2/pwOAZolLTwDMMhzXHoscMsnJyfTPP//Qf/7zH4qOjpZtd+7coSdPnlh6fsBCBIdF056Q+/T93lC6GoHvCYA878XhmkMeR5IdCY7r169LTaru3bvTyJEjKSIiQrbPmDGDxo8fnxdzBBbgScLTMggrDt/I17kAYB/OcSfN5nE45iYBkMuNcDVcDw8P3Xb2e3A2ObD++jlrj92ihGTb/dECYMvhuB7pgiMpRRETsl34OPbs2UP79++XfA59ypYtS7dv37bk3IAFeZLwVFA8jE2izWfDqHudEvk6JwDsM4/DUbfOWoeL0zN38H7umD1jrk3F9aoycuvWLSpQoICl5gUsTKyeqYqBuQqAvCE2h6gqVydHUmsf2qqfw2zBwY2b5s6dq3vOtarYKT5lyhSUIbEBH0ejskXkR3vwaiRdgZMcAIsTn033P/WaqfNz2GihQ7MFx+zZs2nfvn1UrVo1io+PpzfeeENnpmIHObDuSI8K/t7Upoq/rK+E1gFAHiYAOmc5RjVjxSYZWgI06+MoWbIknTp1ilatWiWPrG0MGTJEugLqO8uBdRGTrnF4uTpRz3olaNvFcNE6AADP11RlkARoo7kcZgsObhHbrFkzERS86Od28L6WLVtaeo7AAsSkR1V5uTlTUGBBWWdTFfdHNtZsBgCQN6YqLdSrMttU1aZNG4qMzHynGhUVJfuAdRKTHlXl5eZEZYp6krOjg6jUdx/H5/fUANAUsenmJ1WrMIau77i9CA6uimusedODBw/Iy8vLUvMCeWWqcnOW8L+yvmnf1eVwOMgBsCRx6Q7v7DSOp6aqVG2bqnr16iWPLDTefPNNcnNz0+3j8NzTp0+LCQtYuakq3WFX0c9bhAYvrSr75e/kANBgrSoPF+2aqkwWHIUKFdJpHJyvoe8I52TAJk2a0NChQ/NmlsCCpqp0weHvTXQOGgcAlkRRlKclRzTs43A2tyouh95yTSqYpWw3qkonONId5AAAy5CQnEqp6T3usipyaODjsJeoKk70A7YbW65qHBX80gUHNA4ALEa8ngaRVZFDLfTkMNs5fu/ePRowYAAVL16cnJ2dycnJyWAB1p05zlFVTAX/NI3xQUwiPYxJzNe5AaC1GzRXJ0dyzqYGld2YqlTYMX7jxg2aPHkyFStWzGiEFbDe6riqxuHp6kwlfDzo9qM4uhzxhBp6FcnnGQKgoV4cLtnfk3ukFzq0mwRA7jnOFXLr1KmTNzMCFodLqHMJZ1VgqHD5EREc4U+oYVkIDgAsZaryzKbciBZ6cphtqipVqpREDgDbIVavpLrqHFdDchlEVgFg6TpVTtmOszsfB1fGnThxIl27ds1ik1i4cKFEa7m7u1Pjxo3p8OHD2Y5fs2YNBQUFyXjuRrhp0yaD/SzYPv30UzGlcdhw+/btKSQkxGAMvx+b2fSXr776irTs33BzNrS7qpFVEBwAWIY4E+pUGbSPTbQTwdG7d2/auXMnVahQQfI5ihQpYrCYCxdLHDdunERrHT9+nGrXrk2dOnWi8PBwo+O5iVTfvn2lsOKJEyeoR48espw9e1Y3ZubMmTR//nxavHgxHTp0SEKH+ZhczVefzz//nO7evatbRo0aRVq+C/JO92+oQHAAkEfJf67ZCw41x8NWNQ6zfRz6vTgswddffy2Jg4MHD5bnfLH/888/6YcffhDNJiPz5s2jzp0704cffijPv/jiC9q6dSstWLBAXsvaBs/xk08+kb7ozE8//UQBAQG0YcMG6tOnj+5YLPgCAwNJ66gah2d6RFVGwcF+Dnae52SXBQBkjynJf1rwcZh9pRg0aJDF3jwxMZGOHTtGkyZN0m1zdHQU09KBAweMvoa3s4aiD2sTLBSY0NBQCgsLk2PoZ72zCYxfqy842DTFgqd06dLSV2Ts2LESYmyMhIQEWVQeP35MNhdRlUEwFPFylSUyJpGuRsRQjRJp1QEAAM/o43DRto/DJMHBF8mCBQuadMFUx5nC/fv3pc4VawP68POLFy8afQ0LBWPjebu6X92W1Rhm9OjRVK9ePTGvsfmLhRebq1gDMsb06dPps88+Iy2UG9GHHeSHY9K6AUJwAJC3/cYz5XEkalhwFC5cWC6q/v7+5OPjYzR3Q62aa6wfuTWir7XUqlVL6m0NHz5cBIR+AUcVFiz6r2EByhFmtlYZNyMcknv4WiRdDIumNMMeACC3xJnQi8OwrLqGq+Nu375d5/jesWOHxd7c19dXss05G10ffp6V74G3ZzdefeRtHFWlPya73BM2ZXEzKo4Wq1KlSqb9LEyMCRTbMlVl/jHXK+1DKw7foH2X7+fDzADQZvc/95yiquzBVNWqVSuj688K3+XXr1+ftm3bJpFRTGpqqjx/7733jL6madOmsv/999/XbWPnOG9nypUrJ8KDx6iCgrUDjq565513spzLyZMnxb/CWpXWeJKNqapVlbSS6qdvRVFEdAL5FbBN4QiALWkc7vZgqspL2PzDDvcGDRpQo0aNJCIqJiZGF2U1cOBAKlGihJiQmDFjxojwmj17NnXt2pVWrlxJR48epSVLlsh+NpexUJk2bRpVqlRJBAmXR+HaWqpwYic5CxLuWMiRVfycHeP9+/cXs5zWyE7j8C/gTjVLFKIzt6No16UIerV+yXyYIQAa83G4mJjHkZSSZXM8aybfBQfnhUREREjCHjuvWUvYvHmzzrnNdbFYE1DhZlHLly+XcNuPP/5YhANHVNWoUUM3ZsKECSJ8hg0bRo8ePaLmzZvLMTlhkGGTEwucqVOnSqQUCxcWHBmjtbQXjmv8625TxU8Ex46L4RAcAFgiAdDVtJIjain2nExb1ka+Cw6GzVJZmaY42TAjr732mixZwdKbk/t4MQZHUx08eJDsreRIxgRAldZB/jR/+2XaHRJBSSmp0loWAJD34biqlmJrggNXCDvgSbqpKiu7a+2SPpLPER2fTMevP3zOswNAi0UOnbId5+ToQK7OjjbrIH9mwbFr1y6pFfXwIS441kpsNuG46o9Y7Tu+Izjiuc4NAC36Ez1yEBy2HlllsuCYMWOGOJlV2KHDpT/YwdytWzeqWrUqnTt3Lq/mCSyRAJiN3bV1enTVzmDjNcIAADkTl56XkZOpSt90HBWXRJoVHFyMUN8BvXbtWtq9e7f05uAMcI6KstXMaq0To2vilPWPuWUlP3J0IEkEvPMo7jnODgDtFTn0NEHjUEPfOQxes4KDa0BxhrUKm6deffVVeuGFFyQ5kKOcsqovBaw3c1ylsJcr1S2dFoq8/sTt5zY3ALTZAdDJZMERrmXBwVnV+pnTLCQ4NFaF8yRY8wDWR0z6j9krhxDBfo1Ly+MPe0NtNjEJAFuojsv424PGwf032DSl5lZcunSJWrZsqdt/69YtKlq0aN7MElhI48j+x/xy7eJUqogHPYhJlDIkAIC8KXKoJt8yEdGGfYI0JThGjhwpuRbcQOnFF1+UEh/VqlUzqGdVt27dvJonyCWpqYpOfc7OVMVwd8B3WlWU9f/sviK9ygEAppGUkkrJqWlttT1dck6RswsfBzdb4q56kZGRommsW7fOYP+dO3forbfeyos5AgsUXWO8TGjU9Er9EhRY0J3uPU6gtcdu5e3kANAQsXrmXXdXR5NNVZr2cTAsGNavX0/ffvttpuq1s2bNytQDA1hPDgdHTLm75Px1uzk70fBW5WX9251XKDHZNss+A5BfyX9OnNxnQvUFu9A4ciIkJIRatGhhqcMBC9ep8nJ1NrmQWp+GpcnX25VuPYyjmZuNN9QCABjXODxdnEz6W/Mv+FRwsEnZlkDJEY1jqn9DH3bs/V/PmrL+/d5Q+vvc086JAADjPE5P5DP1b62oV5rgYL/IIxtLAoTg0DhPK+OaV0StY/VAGtK8nKyPX3OKbkbG5sn8ANAKt9MTZ4v7pEVL5QTXquIacUy4jUVWQXDYSe2crCrjZsdHnYOoTikfehyfTO+tOKGz4QIAMqNWXCju40Gm4udtm34Ok68mGzduzDGzHFhv9z9TEpKM3REteKMudZ2/l07dfEQT1p6meX3q2FzTGQCeB7cepgmOEoVNFxzs5wi+F03hjzUqONTuedmBC4r1RlXlRuNgShb2pG/716OB/z1MG0/dofJ+XvR++8oWniUA2jFVlcyFxmFrIbkmm6q4F3hOS0oKTBlW6+MwIYcjK5pV8KVpPdIKXM79J4TWIb8DgEzczoXG4acXWWVLwMehcXITVWWMPo1K07CWafkdH6w5RbP/DqYUGwshBCAvuROVex+HrTnHzb6aPHjwQFeT6ubNm/Tdd99RXFwcvfTSSwa1q4CV1anKhY/DmLOchcV/94bSN9sv09nbUfT163Wksi4A9v539ig2LaS2hBmCw7+gu7Y1jjNnzlDZsmXJ39+fgoKC6OTJk9SwYUOaM2cOLVmyhNq2bUsbNmzI29mCZ+jF8ezt5TkjdnK3ajSnd21yc3aUboEtZ+6gef+EUHS8bcWhA5AX/o2C7s5UwN1F81FVJguOCRMmUM2aNaVCbuvWraXrX9euXSkqKkraxg4fPpy++uqrvJ0tyH33PzPzOLKjZ92StO6dZhQUWICiE5Jpzj+XRID836YLdDn8icXeBwDb8294mvU6/exxW8Lk29AjR45IBVxu5lS7dm3RMt59911ydEyTPaNGjaImTZrk5VxBHjVxyg01ShSiTaNb0J9n7tKcrZfo6v0YWrL7qiz1SvtQl5rFqFP1QCpVxLw/JABskVvpGoc5Zir9elV8A8Yl2U0px24NmHw14aq4amFDb29v8vLyosKF0zrGMbweHR2dN7MEz26qeoaoqqxwdHSgl2oXpxdrBNL2i+G0+uhNeTx+45Es0/68QNWKFaSutYpR15rFqKyvl8XnAIA1Jf+VMDFrXKWAm7MUH41PShWto3RR27jRMutqkjFPA3kbtmSqsrzg0O/jwSVKeLn3OJ42nblLW86F0eHQSDp/97Es/94SDCECNMvtXITiqtdQ1jpuRsZJZJUmBcebb76pax8bHx9PI0aMEM2DSUiwLRudvWDJqCpTCCjoToNfKCdLZEyiFEhkc9b+Kw8MhAj7R9pXDaB2Vf2pdkkf0V4AsHXneAkf8y/83AmQBYct+TlMFhyDBg0yeN6/f/9MYwYOHGiZWQGry+PIDVzAjfM/eMkoRC6GRcuyYMdlKurlSq2q+FGbKv7UsrIfFfIwPSoFAFvWOGw1e9zkq8nSpUvzdiYgb/txWDCq6lmFyMOYRNp5KZz+OR9Ouy5FSI/zX4/floVDfuuXKUwdqgbQy3WKiwYDgDWTmJxK99IT+EytjGssssqWkgCf/20oeG4oCvcbf/aSI5aGEwY5pJcX/qM7dv0h7QwOF8d6SPgT8Y3wMv2vC/RCRV96pV5J6lwjkNxdbCPiBNgX9x7Hk6KkFQX1Te+xYQ62mMthPVcTYHHuRMVTUooid/FFva0zu5v/2JpWKCrLpC5Vpe8HC5A/Tt+hI9ce0p6Q+7IU2uhCPeuWoNcblKJqxQvm97QByFwV18cjV766pxoHBAewAoLDHstjeV8v6SVuC3Dex6BmZWW58SCWfj1xi9YcvSXOx2X7r8nCjnXWQrrXKa4r2QBA/jvGPXL1elvsPQ7BoWGCw9KyuKsEFiBbhEMTuYT7qLaVaE9IBK06cpO2XQgXp/qXmy7oTFk96pSgF2sGWpU5DtifY7x4LvwbalQVA40DWJXGwXfotgyb2lpX8ZflUWyiRGaxI519I6opa+rGc9Sjbgnq26g0TFnguXL7UWyuQ3H1TVUPniRI+Hx+RECaC8qqaxi+M2eqBGrnQurj6Ur9GpeRWlm7P2xD4zpUptJFPKVkw/8OXqcu8/dQ3yUHafvFe5SKsu/gOXDnUXyuQ3FV53iZop7EP9e9l++TLQDBoVGSUlLpSsQTTWgc2ZmyRrerRDvHt6afhzSmLjUDRTs5cPUBvbXsKHWYs4t+OXRdagABYK0+DgcHB2ob5C/r2y7cI1sAgkOjhN6PkYgqzhjP7Q/aVuBIluaVfGlRv/q0e0IbGtqinNQAuhIRQ/9af5aafbWNvv47WMxcAFiS+KSUpy1jc6lxMO2CAuRx+8UIm9CUITg0bqaqHFjArsp5sJD8V9dqtH9SW/q0WzUqVcSDHsYm0fztl+mFr7bTjM0XJQERAEuw+WyY5CKx0HiWG7RG5YqQt5sz3X+SQGduR5G1A8GhUbTiGM8t3EznreblaOf4NrSoXz2qWqwgxSSm0Lc7r1Cb2TvFhIXWt+BZWX30pjy+Vr/UM92gcT5Ty8q+sr7tYjhZOxAcWg/FDbBPwaHCPg/uDbJpdHP6bmADEaTc4pNNWD0X7ZP2twDkhpuRsVJ3jYuEv1K/xDMfr63OXGX9fg4IDo0SfO+x5iKqngV2QHaoFkB/jGpOU16qJj6Q07eiqPvCfTRrSzAlJMOBDsxjTbq20byiL5U0s/OfMVpX8RMhdPb2YwqLsu66VRAcGi1syGWabTn5Ly97h3DJ923jW0lvEDZXcYXebvP30plb0D6AaaSkKrT22C1Zf61BKYsc09fbjeqW8pF1LrtjzUBwaJBL96J1pQy4Ki0wnq278I16tLh/PfmD5eKKbLpasD2EklNS83t6wMrZd/m+1ILjFgAdq6WZmCxBu6ppx/rt5G2rjq6C4NAgwekRVfbqGDeHzjWK0daxLaUrYXKqQrP+vkS9lxwU+zUAWeVIfbfnqqz3qFPcolWb2R/n4uRAh0IjadbfwWStQHBoWHDYu2PcnDLvC96oS1+/Xlt8H1zKpMu8PfT7qTv5PTVghXkbI/53TMrcODs6UL8mZSx6/HK+XvRVr1qyvmjnFVp15AZZI1YhOBYuXEhly5Yld3d3aty4MR0+fDjb8WvWrKGgoCAZX7NmTdq0aVOmPhSffvopFStWjDw8PKh9+/YUEhJiMCYyMpL69etHBQsWJB8fHxoyZAg9eZIWiWTLcL2b3SERsg7/hnnO8171StKmMS2kkRSXMBm14gSNW31SPlMAuIvloB8OS7ism7OjROlVzoObs1fql5SKCAxH/3EpHWsL3sh3wbFq1SoaN24cTZkyhY4fP061a9emTp06UXi4cefQ/v37qW/fvnKhP3HiBPXo0UOWs2fP6sbMnDmT5s+fT4sXL6ZDhw5JX3Q+JvdJV2Ghce7cOdq6dSv98ccftHv3bho2bBjZMhxa+vKCfXQ1IkbunLkNKzC/rPuqYU1odNuKEuHCxRTbzNpJ/ztwDXkfdkh0fJI0GXtv+XFq8n/bxIRUwM2Z/jekMbVJLxOSF4xtX0naBrD5dPKGs9R8xg5auOOy/I2zqSy/cVD49jwfYQ2jYcOGtGDBAnmemppKpUqVolGjRtHEiRMzje/duzfFxMTIxV6lSZMmVKdOHREUfDrFixenDz74gMaPHy/7o6KiKCAggJYtW0Z9+vShCxcuULVq1ejIkSPUoEEDGbN582bq0qUL3bp1S16fE48fP6ZChQrJsVlrMQV2uqqtXM2Fv6VURaEUXlK5s18KRccn08PYRAq5F00X70bTprN3KT4pVdTd7wbWp4r+0DieBTZZ8R/t+btpoc3cG53LuKeFX3pIkqG3u7PYpJ0dHcnRkciB/3PgxzQthtFPC0vfBJ4B/b8FdZ0f+UKmOpR5G6/KY+rTv5vkFIWSU1OlHE9CcqpkfcclpUhVWl7474n7YnCJ88vhT3RNmlSqFy9IM16pRTVKFMrz82QB8eP+a/T9nlAKe/z0ppe1HfZfFivkQYGF3OV3WcDdmbzdXcjT1YncXRyl/w43bwsyMxzf1OtavtbvTUxMpGPHjtGkSZN02xwdHcW0dODAAaOv4e2soejD2sSGDRtkPTQ0lMLCwuQYKvxBsIDi17Lg4Ec2T6lCg+Hx/N6sofTs2TPT+yYkJMii/wGbC9dO6jR3N+UlrSr70fy+dSXaAzwbbLL6fVRzyTKf/fcl6Y2+8dQdWYD9UKyQO7WvGkC9G5Z6LgJDxcXJkd5uUZ4GNi0rUVa/nbxDp249khvGU7eiZMmOZhWK0vKhTSgvyFfBcf/+fUpJSRFtQB9+fvHiRaOvYaFgbDxvV/er27Ib4+9vqGY6OztTkSJFdGMyMn36dPrss88ov+E7VnbKebg4yR0v32mU9/OSO4uaJQtRy0p+ki0NLAN/lvyH26dhaTp+4yHtDblPh69FSr2rx/FJ9CQ+WcwJfDer3gGD5/e34OjgQPxzZ02P/1fX+XtTNT9e54XH8sWY152dHOSunO/e+Q6d60RxHwwfDxcJY+cQ7bK+/HdVQEr55yeuzo6SK8ILa0+hD2Io5N4T6XXOmgj/FlmY8O+RnfesSfFjKQskJWaF9XcMsRJYK9LXdFjjYJOaOVQO8KbLX76Y6znIHwmEQr798TYpX1QWU1D0zCemjAXmkSYY7O9vwdHRgSr4ecuSn+Sr4PD19SUnJye6d8+wNgs/DwwMNPoa3p7dePWRt3FUlf4Y9oOoYzI635OTkyXSKqv3dXNzk+VZ4B863+kA7cPftenXNfwmgG2Rr1FVrq6uVL9+fdq2bZtuGzvH+XnTpk2Nvoa3649nODJKHV+uXDm5+OuPYe2AfRfqGH589OiR+FdUtm/fLu/NvhAAAADZoOQzK1euVNzc3JRly5Yp58+fV4YNG6b4+PgoYWFhsn/AgAHKxIkTdeP37dunODs7K7NmzVIuXLigTJkyRXFxcVHOnDmjG/PVV1/JMX777Tfl9OnTSvfu3ZVy5copcXFxujGdO3dW6tatqxw6dEjZu3evUqlSJaVv374mzzsqKortC/IIAABawNTrWr4LDuabb75RSpcurbi6uiqNGjVSDh48qNvXqlUrZdCgQQbjV69erVSuXFnGV69eXfnzzz8N9qempiqTJ09WAgICRCi1a9dOCQ4ONhjz4MEDERTe3t5KwYIFlcGDByvR0dEmzxmCAwCgNUy9ruV7Hoetkps8DgAA0MJ1Ld8zxwEAANgWCMfNJaqilptEQAAAsEbU61lOhigIjlwSHZ1WgdbcXA4AALCF6xubrLICPo5cwqG7d+7coQIFCuRpIpKaaHjz5k278aXY4znb63njnAuSNcHigIUG1+vjEkxZAY0jl/CHWrJkyef2fvwDs7YfWV5jj+dsr+eNc7YestM0VOAcBwAAYBYQHAAAAMwCgsPK4fpY3OTqWetk2RL2eM72et44Z9sEznEAAABmAY0DAACAWUBwAAAAMAsIDgAAAGYBwQEAAMAsIDislKlTp6Z3kXu6BAUFkZbYvXs3vfTSS5Klyue3YcMGg/0ct/Hpp59KJ0cPDw9q3749hYSEkJbP+c0338z0vXfu3JlsmenTp1PDhg2lyoK/vz/16NGDgoODDcbEx8fTyJEjqWjRouTt7U2vvPJKpk6fWjvn1q1bZ/quR4wYQbYABIcVU716dbp7965u2bt3L2mJmJgYql27Ni1cuNDo/pkzZ9L8+fNp8eLF0sHRy8uLOnXqJBcZrZ4zw4JC/3tfsWIF2TK7du0SoXDw4EHp1pmUlEQdO3aUz0Jl7Nix9Pvvv9OaNWtkPJfz6dWrF2n5nJmhQ4cafNf8m7cJnk97EGAu3Nmwdu3air3AP8X169cbNOMKDAxU/v3vf+u2PXr0SBpzrVixQtHiOTPctIw7VmqZ8PBwOfddu3bpvlfu4rlmzRrdGO7uyWMOHDigaPGc1SZ1Y8aMUWwRaBxWDJtl2KRRvnx56tevH924cYPshdDQUAoLCxPzlH4NHe4Jf+DAAdIyO3fuFPNGlSpV6J133qEHDx6QluAmQUyRIkXk8dixY3JHrv9ds1m2dOnSmvmuozKcs8ovv/xCvr6+VKNGDZo0aRLFxsaSLYAih1YKXyCXLVsmFw9WYT/77DNq0aIFnT17VuymWoeFBhMQEGCwnZ+r+7QIm6nYRFOuXDm6cuUKffzxx/Tiiy/KBdTJyYm0UFX6/fffpxdeeEEulgx/n66uruTj46PJ7zrVyDkzb7zxBpUpU0ZuDk+fPk0fffSR+EF+/fVXsnYgOKwUvlio1KpVSwQJ/8hWr15NQ4YMyde5gbyjT58+uvWaNWvKd1+hQgXRQtq1a0e2Dtv9+eZHa/663JzzsGHDDL5rDgLh75hvGPg7t2ZgqrIR+G6scuXKdPnyZbIHAgMD5TFjZA0/V/fZA2ymZFOGFr739957j/744w/asWOHQUsC/j4TExPp0aNHmvuu38vinI3BN4eMLXzXEBw2wpMnT+ROhO9K7AE21fBFY9u2bQYNcDi6qmnTpmQv3Lp1S3wctvy9cxwAX0DXr19P27dvl+9Wn/r165OLi4vBd80mG/bp2ep3reRwzsY4efKkPNrCdw1TlZUyfvx4ifdn8xSHJnI1TbZx9+3bl7QkDPXvrtghzn887EBkxyjbhadNm0aVKlWSP7zJkyeLPZhj4rV4zrywL4tzGFho8o3ChAkTqGLFihKGbMummuXLl9Nvv/0m/jnVb8HBDpyfw49sfh03bpx8BtzcaNSoUSI0mjRpQlo85ytXrsj+Ll26SO4K+zg4JLlly5ZinrR68jusCxind+/eSrFixRRXV1elRIkS8vzy5cuKltixY4eEKGZcOCRVDcmdPHmyEhAQIGG47dq1U4KDgxWtnnNsbKzSsWNHxc/PT8JTy5QpowwdOlQJCwtTbBlj58vL0qVLdWPi4uKUd999VylcuLDi6emp9OzZU7l7966i1XO+ceOG0rJlS6VIkSLy265YsaLy4YcfKlFRUYotgLLqAAAAzAI+DgAAAGYBwQEAAMAsIDgAAACYBQQHAAAAs4DgAAAAYBYQHAAAAMwCggMAAIBZQHAATcDpSFw0jjOPuZOaWr4BZIbrQnE2+v79+/NtDhMnTpTscGCbQHAATbB582YpQ88F5bgMvX75amAId1TkEi7NmjUz2M6F+Lp160Z+fn7k7u4uFVp79+4t7W5VuEovC+aMBQmZsmXL0ty5c00uqfPjjz/S1atXLXBG4HkDwQE0gVoAki+GXOfJ2dnZ6J22vcOa2YIFCzKV5l+0aJGU9Oa6SatWrZIig1ygjz9PrqFkabjiL9ff+vbbby1+bPAcyO+aJwA8K1znSb8eENd4Ultzjhw5UtpzFi1aVGndurVsP3PmjNK5c2fFy8tL8ff3V/r3769ERETojvfkyRNlwIABsp/b186aNStTm09jbV8LFSpkUH+J6xG99tprsp1rML388stKaGhopjax3B6X34frFnG9psTERN2Y+Ph4ZcKECUrJkiWlblmFChWU77//Xup48bp+a13mxIkTMreQkBCjn9WRI0cUR0dH5fHjx7pt169fl9pYY8eONfoafq+MtbYePnyYaRx/7nPmzJF1/hyM1WrilsgqP/74o5wXsD2gcQCbZ968efT5559LvwM2Ux05ckS3j80h3F1u3759YqJhE0vbtm2pbt26dPToUTFxcd+H119/XfeaDz/8kHbt2iWVTf/++28xzxw/ftysOXErVL6j5sqoe/bskff39vaWDn/6mg+bh1hb4keeK5vbeFEZOHAgrVixgubPn08XLlyg//znP3IcNhe99dZbtHTpUoP35edcYZV9GMbguXBfF/0ukuvWrZP5ciVeY/B7mQubuPi7UBc+B9YCuQueSqNGjaRs/LVr18w+Pshn8ltyAWAJ+E5X1TRUWEuoW7euwbYvvvhCKtDqc/PmTbkb5sq70dHRcme/evVq3f4HDx4oHh4eZmkc//vf/5QqVaoY3K0nJCTIcbZs2aLTOHjOycnJujGsoXAlZIbnw++zdetWo+d8+/ZtxcnJSTl06JA8Z03F19dXWbZsWZafE59D27ZtDbaNGDFCKViwoMG2tWvXisalLqdPnzbQOPT3qYuDg4NO49CHqzqzNjVz5kyD7VwJlo+1c+fOLOcLrBP04wCahpsE6XPq1Cm5u+e79ozwnX9cXJxoBGo3NoYjtbj3uznw+3DfjYz94ePj4+V9VKpXr27QS5z9NGfOnJF1jgzjfa1atTL6HtybpGvXrvTDDz/I3fvvv/9OCQkJ9Nprr2U5Lz4/dnznpFWwtsTvf/v2bWrdujWlpKRk0lwynhuPy0hUVJQ43HmerMnpw30pmNjY2CznC6wTCA6gaby8vDI1UuIGWTNmzMg0li/aprbt5Attxo4EbO7Rfx8WWr/88kum13LUkgp3vst43NTUVIMLa3a8/fbbNGDAAJozZ46YqdhE5Onpma1TWhVMKtwoiy/w3GxIbdXKgpXNXcaCDBiOyuJ2xvpkHMvChufDjZmWLFmS6RiRkZGZPg9gG8DHAeyKevXq0blz5yR0lC+M+gsLGQ5B5Ys5t6hVefjwIV26dMngOHyxY9u9SkhIiMGdM78Pb/P398/0PtwFzhRq1qwpQoT9LVnBHeR43hydxP4a9ntkB/t2Ll68aCD0Xn31VTlnY8L0WeBoLBZSGzZsMKrlnD17Vt6XtS5gW0BwALuCW3rynS634GUnOpuNtmzZQoMHD5Y7ZL7T5lBVNqtwr2i+uL355pvk6Gj4p8IOdg5rPXHihDjZR4wYYaA99OvXT+7uu3fvLmYdbhHLTvbRo0eLQ9gUWLgNGjRIhAFffNVjrF69WjeGTVk8v0mTJonmkFOP7jZt2og2xMJThdv0zp49W4IM+P3YlMcOaw4IYKe8+j7mwNoPh/hyQAJrUazN8MLvrcKfS4sWLUzSrIB1AcEB7Ar2C3CEEwuJjh07yl099zZns4sqHP7973/LBY1NWu3bt6fmzZtn8pXwhbZUqVIy7o033pCENn0TEa9z4hxflHv16kVVq1YVgcQ+DjbdmAprEqwRvPvuuxQUFERDhw6lmJgYgzF8XPbLsPDLCc7T6NmzZyYTGmdxcwRZRESEvB8LIdZmWFixJsOfkzmwlsSf8csvvywmQHWZNWuWbszKlSvlfIDtgdaxAJgAO37r1Kljcmb084Tv3Dl57+bNmxQQEJDj+NOnT1OHDh1E2zIWJPA8+Ouvv+iDDz6QuWTlRwHWCzQOAGwUjqBis9fUqVMlksoUocHUqlVL/BmsTeQXrDWxOQtCwzbBtwaAjcJJdWymYk3op59+Muu17BfJT9gcBmwXmKoAAACYBUxVAAAAzAKCAwAAgFlAcAAAADALCA4AAABmAcEBAADALCA4AAAAmAUEBwAAALOA4AAAAGAWEBwAAADIHP4fSrT3q/gPH5oAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the BLS signal\n",
    "plt.figure(figsize=(4, 2.5))\n",
    "plt.plot(w_common/2/np.pi/1e9, np.abs(sigma), label='Kalinikos-Slavin model')\n",
    "plt.xlabel('frequency (GHz)')\n",
    "plt.ylabel('BLS intensity (a.u.)')\n",
    "plt.title('BLS signal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca7f418a",
   "metadata": {},
   "source": [
    "*(Note: Sometimes the relative amplitude of PSSW peaks does not match the experiment, since it may depend on a vast number of parameters. If that is the case, calculate the BLS signal separately and normalize to each respective peak.)*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b0acfd9",
   "metadata": {},
   "source": [
    "### Tip on speeding things up\n",
    "Calculations of 2D dispersion and Bloch functions using some dispersion models can be greatly accelerated when we pass the `K` and `PHI` as flattened arrays. What happens is that if the model does not iterate over wavenumbers, it can calculate the dispersion for the whole 1d array of wavenumbers at once (therefore there is no need for the double for loop). If we supply a parameter of the same shape, it may be broadcasted together. The only thing we need to do then is to reshape the acquired data back to our square reciprocal space grid.\n",
    "\n",
    "Unfortunately, this is currently supported only by the `SingleLayer` class at the moment. Maybe we will try to implement it in a limited form also in other classes.\n",
    "\n",
    "For example this code from above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92408f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ...\n",
    "# Create a SingleLayer model for the current kxi and phi.\n",
    "# Note: We pass kxi and phi as arbitrary values for now.\n",
    "model = SWT.SingleLayer(Bext=Bext, kxi=(0,), theta=theta,\n",
    "                        phi=0, d=d_layer, material=material)\n",
    "# Loop over all grid points in the Kx,Ky plane.\n",
    "for i in range(Nk):\n",
    "    for j in range(Nk):\n",
    "        # Set correct kxi and phi for the current grid point.\n",
    "        model.kxi = np.array((K[i, j],))  # must be an array\n",
    "        model.phi = PHI[i, j]\n",
    "        # Compute the Bloch functions for n=0,1; higher n are outside `w_common`\n",
    "        # The returned w has shape (Nf_common,) and bf has shape (Nf_common, 1),\n",
    "        # since we are passing kxi as length-1 array.\n",
    "        # We pass also parameters for Bose-Einstein distribution (optional).\n",
    "        w0, bf0 = model.GetBlochFunction(n=0, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "        w1, bf1 = model.GetBlochFunction(n=1, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "        # Interpolate the Bloch functions to the common frequency axis.\n",
    "        bf0_interp = np.interp(w_common, w0, bf0[:, 0], left=0, right=0)\n",
    "        bf1_interp = np.interp(w_common, w1, bf1[:, 0], left=0, right=0)\n",
    "\n",
    "        # Sum Bloch functions for n=0,1\n",
    "        Bloch2D[:, i, j] = bf0_interp + bf1_interp\n",
    "\n",
    "Bloch3 = [Bloch2D, np.zeros((Nf_common, Nk, Nk), dtype=complex), Bloch2D*1j]\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7da22565",
   "metadata": {},
   "source": [
    "can be replaced by this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c36dc222",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ...\n",
    "# Create a SingleLayer model for the current kxi and phi.\n",
    "# Note: We pass kxi and phi as flattened arrays of same shape.\n",
    "model = SWT.SingleLayer(Bext=Bext, kxi=K.flatten(), theta=theta,\n",
    "                        phi=PHI.flatten(), d=d_layer, material=material)\n",
    "# Compute the Bloch functions for n=0,1.\n",
    "# The returned w has shape (Nf_common,) and bf has shape (Nf_common, len(kxi))\n",
    "w0, bf0 = model.GetBlochFunction(n=0, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "w1, bf1 = model.GetBlochFunction(n=1, Nf=Nf_common, temp=300, mu=-1e12*SWT.H)\n",
    "# reshape Bloch functions to match the Kx,Ky grid shape\n",
    "bf0 = bf0.reshape((Nf_common, Nk, Nk))\n",
    "bf1 = bf1.reshape((Nf_common, Nk, Nk))\n",
    "# Loop over all grid points in the Kx,Ky plane.\n",
    "for i in range(Nk):\n",
    "    for j in range(Nk):\n",
    "        # Interpolate the Bloch functions to the common frequency axis.\n",
    "        bf0_interp = np.interp(w_common, w0, bf0[:, i, j], left=0, right=0)\n",
    "        bf1_interp = np.interp(w_common, w1, bf1[:, i, j], left=0, right=0)\n",
    "\n",
    "        # Sum Bloch functions for n=0,1,2\n",
    "        Bloch2D[:, i, j] = bf0_interp + bf1_interp\n",
    "\n",
    "Bloch3 = [Bloch2D, np.zeros((Nf_common, Nk, Nk), dtype=complex), Bloch2D*1j]\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6846a26",
   "metadata": {},
   "source": [
    "which is much faster in most cases. More specifically, here it is about 100× faster, but is depends greatly on number of points in the grid (more convenient with large number of points).\n",
    "\n",
    "*(Note that the double for loop was not totally eliminated, since we need to perform the interpolation of the Bloch functions at each point.)*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv (3.13.5)",
   "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.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}