{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f4d55ed",
   "metadata": {},
   "source": [
    "# Modeling BLS signal from single magnetic layer\n",
    "This example shows how SpinWaveToolkit (SWT) can be used to calculate the expected BLS signal of a single magnetic layer. \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 toolkit)\n",
    "Nf_common = 51     # 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_min = 1e-6         # minimum k (avoid zero if necessary)\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, KY = np.meshgrid(kx_grid, ky_grid, indexing='ij')"
   ]
  },
  {
   "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 `getBLSsignal` function."
   ]
  },
  {
   "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 Slavin-Kalinikos model...\n",
      "Computing BLS signal from Slavin-Kalinikos model...\n"
     ]
    }
   ],
   "source": [
    "# Now Bloch2D is an array of shape (Nf_common, Nk, Nk), where for each frequency,\n",
    "# you have the Bloch amplitude defined on the Kx,Ky grid.\n",
    "print(\"Preparing focal field...\")\n",
    "objective = SWT.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",
    "\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(2*np.pi*3e9, 2*np.pi*18e9, Nf_common)\n",
    "\n",
    "print(\"Preparing Bloch functions in 3D (f,kx,ky) using Slavin-Kalinikos model...\")\n",
    "# Loop over all grid points in the Kx,Ky plane.\n",
    "for i in range(Nk):\n",
    "    for j in range(Nk):\n",
    "        # Compute polar coordinates for the current (KX, KY) point.\n",
    "        kxi_val = np.sqrt(KX[i, j]**2 + KY[i, j]**2)\n",
    "        # For a nearly zero radius, set phi = 0 (or any default value) to avoid NaN.\n",
    "        phi_val = np.arctan2(KY[i, j], KX[i, j]) if kxi_val > 1e-12 else 0.0\n",
    "\n",
    "        # Create a SingleLayer model for the current kxi and phi.\n",
    "        # Note: We pass kxi as a one-element array.\n",
    "        model = SWT.SingleLayer(Bext=Bext, kxi=np.array([kxi_val]),\n",
    "                                theta=theta, phi=phi_val,\n",
    "                                d=d_layer, material=material)\n",
    "        # Compute the Bloch functions for n=0,1,2.\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)\n",
    "        w1, bf1 = model.GetBlochFunction(n=1, Nf=Nf_common)\n",
    "        w2, bf2 = model.GetBlochFunction(n=2, Nf=Nf_common)\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, bf0[:, 0], left=0, right=0)\n",
    "        bf2_interp = np.interp(w_common, w2, bf0[:, 0], left=0, right=0)\n",
    "\n",
    "        # Sum Bloch functions for n=0,1,2\n",
    "        Bloch2D[:, i, j] = bf0_interp + bf1_interp + bf2_interp\n",
    "\n",
    "\n",
    "Bloch3 = [Bloch2D, np.zeros((Nf_common, Nk, Nk), dtype=complex), Bloch2D*1j]\n",
    "\n",
    "\n",
    "print(\"Computing BLS signal from Slavin-Kalinikos model...\")\n",
    "# Compute the BLS signal using the Bloch functions and the electric field.\n",
    "Ex, Ey = SWT.getBLSsignal(SweepBloch=w_common, KxKyBloch=[kx_grid, ky_grid], Bloch=np.array(Bloch3),\n",
    "                          Exy=Exy, E=E, Nq=50, DF=[1, -8.1653 + 1j*15.348, 17.237 + 1j*0.43004],\n",
    "                          PM=[1,1,1], d=[d_layer], NA=NA, source_layer_index=1, output_layer_index=0,\n",
    "                          wavelength=532e-9, collectionSpot=0.5e-6, focalLength=1e-3)"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEUCAYAAAAyfG1zAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABABklEQVR4nO3deVxU5f4H8M8ZGIYdAVllU1FxF8UNN9TEFbdMKkNcskwz06um997UupXmmuXPrRIs70UrwSyLcmFxFxXUcENFRdkUZdgHmHl+f9CMDgwwB2aF7/v1mlfNmXOe8/UAXx6elWOMMRBCCDE6An0HQAghpGEogRNCiJGiBE4IIUaKEjghhBgpSuCEEGKkKIETQoiRogROCCFGihI4IYQYKUrghBBipCiBE72LjIwEx3FKLycnJwQFBeHXX3+tcT7HcXj33XfrLLOiogI7d+5E79694eDgAEtLS3h7e2PChAmIiYlpdMz37t0Dx3GIjIxsdFmaYGjxEN2gBE4MRkREBM6cOYPTp09j165dMDExQUhICH755RfeZYWFhWHBggUYOnQo9u7di19++QX//ve/YWpqij/++KPRsbq5ueHMmTMYO3Zso8sipKFM9R0AIXJdunRBQECA4v2oUaNgb2+PqKgohISEqF1Oeno69u/fj5UrV+Kjjz5SHB8+fDjmzJkDmUzW6FhFIhH69evX6HIIaQyqgRODZW5uDjMzMwiFQl7X5eXlAaiqJasiENT/bf/jjz+ib9++sLOzg6WlJdq0aYNZs2YpPq+tyeLnn39Gt27dIBKJ0KZNG2zZsgWrV68Gx3FK58mbgb7//nt07NgRlpaW6N69e40mo9u3b2PmzJlo164dLC0t0apVK4SEhODq1avqPArSxFENnBgMqVSKyspKMMaQk5OD9evXo7i4GK+//jqvcjp27IgWLVrgo48+gkAgQHBwMHx8fNS+/syZMwgNDUVoaChWr14Nc3Nz3L9/H8ePH6/zutjYWEyePBmDBw/G/v37UVlZiQ0bNiAnJ0fl+YcPH0ZSUhI+/vhjWFtbY926dZg0aRJu3ryJNm3aAAAyMzPh6OiItWvXwsnJCU+fPsWePXvQt29fJCcno0OHDmr/u0gTxAjRs4iICAagxkskErFt27bVOB8Amz9/fp1lHj58mLVs2VJRlqOjI3vllVfYoUOH6o1nw4YNDADLz8+v9Zz09HQGgEVERCiO9e7dm3l6ejKJRKI4VlhYyBwdHVn1HzUAzMXFhRUUFCiOZWdnM4FAwNasWVPrfSsrK1l5eTlr164dW7RoUZ3xkKavyTShJCYmIiQkBO7u7uA4DgcPHtT7/aqPrJC/1q9fr9XYjNV3332HpKQkJCUl4ffff0d4eDjmz5+PrVu38i5rzJgxePDgAWJiYrBkyRJ07twZBw8exPjx4+sdwdK7d28AwNSpU/HDDz/g0aNH9d6vuLgYFy5cwMSJE2FmZqY4bm1tXWv7/dChQ2FjY6N47+LiAmdnZ9y/f19xrLKyEp999hk6deoEMzMzmJqawszMDGlpabh+/Xq9cZGmrckk8OLiYnTv3r1BP+zaul9WVpbSa/fu3eA4Di+//LJOYjQ2HTt2REBAAAICAjBq1Cjs3LkTwcHBWLZsGfLz83mXZ2FhgYkTJ2L9+vVISEjA7du30alTJ/zf//0fUlNTa71u8ODBOHjwICorKzF9+nR4eHigS5cuiIqKqvWaZ8+egTEGFxeXGp+pOgYAjo6ONY6JRCKUlpYq3i9evBgffvghJk6ciF9++QXnzp1DUlISunfvrnQeaZ6aTAIfPXo0PvnkE0yePFnl5+Xl5Vi2bBlatWoFKysr9O3bF/Hx8Vq7HwC4uroqvX7++WcMHTpU0b5J6tetWzeUlpbi1q1bjS7Ly8sLb731FgDUmcABYMKECTh27BjEYjHi4+Ph4eGB119/HWfOnFF5vr29PTiOU9nenZ2d3eCY9+7di+nTp+Ozzz7DyJEj0adPHwQEBODJkycNLpM0HU0mgddn5syZOHXqFPbt24crV67glVdewahRo5CWlqaT++fk5ODw4cOYPXu2Tu7XVKSkpAAAnJyc1L6msLAQRUVFKj+TNzu4u7urVZZIJMKQIUPw+eefAwCSk5NVnmdlZYWAgAAcPHgQ5eXliuNFRUUqJyOpi+M4iEQipWOHDx9Wq1mHNH3NYhTKnTt3EBUVhYcPHyp+cJcsWYLY2FhERETgs88+03oMe/bsgY2NTZ019ubur7/+QmVlJYCqoYDR0dE4cuQIJk2ahNatWyude+fOHfz00081yujUqRNKSkowcuRIvPrqqxgyZAjc3Nzw7NkzHD58GLt27UJQUBACAwNrjWPlypV4+PAhhg8fDg8PD+Tn52PLli0QCoUYMmRIrdd9/PHHGDt2LEaOHImFCxdCKpVi/fr1sLa2xtOnTxv0TMaNG4fIyEj4+fmhW7duuHjxItavXw8PD48GlUealmaRwC9dugTGGNq3b690XCKRKNoh7927VyNJVNfQDjUA2L17N6ZNmwZzc/MGXd8czJw5U/H/dnZ2aN26NTZt2oR58+bVODc2NhaxsbE1jq9atQrvv/8+Fi9ejOPHj+Pnn3/G48ePIRQK0a5dO3zyySdYvHhxnWPB+/btiwsXLuCDDz7A48eP0aJFCwQEBOD48ePo3LlzrdeNGjUKBw4cwMqVKxEaGgpXV1fMmzcPmZmZ+P7773k+jSryXxxr1qxBUVERevbsiejoaPz73/9uUHmkaeEYa3q70nMch5iYGEycOBEAsH//fkybNg2pqakwMTFROtfa2hqurq6oqKjAnTt36izX3t5eZYdU9ftVd+LECQwePBgpKSno3r17g/5NxDhVVFSgR48eaNWqFf788099h0OamGZRA/f394dUKkVubi4GDRqk8hyhUAg/Pz+t3P/bb79Fr169KHk3A7Nnz8aIESPg5uaG7Oxs7NixA9evX8eWLVv0HRppgppMAi8qKsLt27cV79PT05GSkgIHBwe0b98e06ZNw/Tp07Fx40b4+/vjyZMnOH78OLp27YoxY8Zo9H5eXl6K4wUFBfjxxx+xcePGxv0DiVEoLCzEkiVLFM02PXv2xG+//YaXXnpJ36GRpki/84g0Jy4uTuVsvvDwcMYYY+Xl5WzlypXMx8eHCYVC5urqyiZNmsSuXLmilfvJ7dy5k1lYWNQ5q48QQhqiSbaBE0JIc9BsxoETQkhTQwmcEEKMlFF3YspkMmRmZsLGxqbGesuEEGKMGGMoLCyEu7t7vWvXG3UCz8zMhKenp77DIIQQjcvIyKh3xq1eE/jq1auVtrwCqlZuU3fxH/lSnBkZGbC1tdV4fIQQomsFBQXw9PRUWmq4NnqvgXfu3BlHjx5VvK8+U7Iu8mYTW1tbSuCEkCZFnWZhvSdwU1NTuLq66jsMQggxOnofhZKWlgZ3d3e0bt0ar776Ku7evavvkAghxCjotQbet29ffPfdd2jfvj1ycnLwySefIDAwEKmpqSp3K5FIJJBIJIr3BQUFugyXEEIMikHNxCwuLkbbtm2xbNkyLF68uMbnqjo9AUAsFlMbOCGkSSgoKICdnZ1aeU3vTSgvsrKyQteuXWvdJWfFihUQi8WKV0ZGho4jbJoqpTIY0O9xQoiaDCqBSyQSXL9+HW5ubio/F4lEihEnNPJEM0rKKzF4XRze/v6ivkMhhPCk1zbwJUuWICQkBF5eXsjNzcUnn3yCgoIChIeH6zOsZiX9STEyxWV4XCSBTMYgENCMVkKMhV4T+MOHD/Haa6/hyZMncHJyQr9+/XD27Fl4e3vrM6xmpbRcCgCokDLkFZfDyUZUzxWEEEOh1wS+b98+fd6eACj5O4EDQLa4jBI4IUbEoNrAie4pJfCCMj1GQgjhixJ4M1daUan4/2xxqR4jIYTwRQm8mXuxBp4lpho4IcaEEngzV0pNKIQYLUrgzVz1TkxCiPGgBN7MUQInxHhRAm/mSsufd2JmictoSj0hRoQSeDP3Yg28tEKKgrLKOs4mhBgSSuDNXEmFVOk9NaMQYjwogTdzZeXKCTyLxoITYjQogTdzJdUSeA4NJSTEaFACb+bkTSjyNVBoMg8hxoMSeDMnH4XSpqUVAGoDJ8SYUAJv5uRNKG2crAHQbExCjAkl8GZOPpWeauCEGB9K4M3c8xp4VQKnNnBCjAcl8GZMJmMorVBuQhGXVigtcEUIMVyUwJuxssrnidrZRgQrMxMA1A5OiLGgBN6MvTgG3EJoAlc7cwA0mYcQY0EJvBmTN5WYCwUQCDhFAqfJPIQYB0rgzZi8Bm5pVrW3tautBQDqyCTEWPDalV4sFiMmJgYnTpzAvXv3UFJSAicnJ/j7+2PkyJEIDAzUVpxEC0r+nsRjIaxq+3b7uwZOQwkJMQ5q1cCzsrIwZ84cuLm54eOPP0ZxcTF69OiB4cOHw8PDA3FxcRgxYgQ6deqE/fv3aztmoiGlihp4VQJ3oQROiFFRqwbevXt3TJ8+HefPn0eXLl1UnlNaWoqDBw9i06ZNyMjIwJIlSzQaKNG8kmoJ3M327wRObeCEGAW1EnhqaiqcnJzqPMfCwgKvvfYaXnvtNTx+/FgjwRHtki9kZfF3An8+CoUSOCHGQK0mlPqSd2PPB4A1a9aA4zi8//77vK8lDVNWvRPz7wT+pEiCCqlMb3ERQtSjsVEoWVlZePDgQYOuTUpKwq5du9CtWzdNhUPUoOjE/LsG7mBpBjMTARgDcgsl+gyNEKIGjSXwYcOGoXXr1ryvKyoqwrRp0/D111/D3t5eU+EQNSiaUP4ehSIQcHCxq1oXPJsm8xBi8DSWwL/77jscP36c93Xz58/H2LFj8dJLL9V7rkQiQUFBgdKLNFz1USgA4CrvyBRTDZwQQ8drHHhdevfuzfuaffv24dKlS0hKSlLr/DVr1uCjjz7ifR+imnwUisWLCdzOAsAzmk5PiBHQ20zMjIwMLFy4EHv37oW5ubla16xYsQJisVjxysjI0HKUTZtiGKHw+e9xmsxDiPHgXQMXCATgOK7Wz6VS9ZYivXjxInJzc9GrVy+laxMTE7F161ZIJBKYmJgoXSMSiSASifiGTGoh307txSYUl7+bULJoLDghBo93Ao+JiVF6X1FRgeTkZOzZs4dX88bw4cNx9epVpWMzZ86En58fPvjggxrJm2ieqiYUeQ08h2rghBg83gl8woQJNY5NmTIFnTt3xv79+zF79my1yrGxsakxq9PKygqOjo61zvYkmiXfzEGpE5Mm8xBiNDTWBt63b18cPXpUU8URHag+lR54Pgolp6AMMhnTS1yEEPVoZBRKaWkpvvrqK3h4eDSqnPj4eE2EQ9T0vAnl+beBk40IAg6olDHkFZfDyYb6HAgxVLwTuL29vVInJmMMhYWFsLS0xN69ezUaHNEuVZ2YQhMBnGxEyCmQIFtcRgmcEAPGO4F/8cUXSu8FAgGcnJzQt29fmklpZBQ1cKFyh7GrrTlyCiTIEpeiq4edPkIjhKiBdwIPDw/XRhxED1TNxASqOjIvPxTT1mqEGDjaUq2ZYowp1kKxNFP+Pe5mR1urEWIMNJbAO3bsSGO3jUi5VAbp36NMLFTUwAGajUmIodPYWihr1qyBWCzWVHFEy+TNJ4CKJhTamYcQo6CxBD5x4kRNFUV0QD6JR2jCQWii/IcY1cAJMQ7UBt5M1TYCBXg+nT5LXAbGaDIPIYaqQTXwn376CT/88AMePHiA8vJypc8uXbqkkcCIdpWqWAdFTr6gVWmFFAVllbCzEOo0NkKIenjXwL/88kvMnDkTzs7OSE5ORp8+feDo6Ii7d+9i9OjR2oiRaEFJueoRKABgLjSBvWVV0qZmFEIMF+8Evm3bNuzatQtbt26FmZkZli1bhiNHjuC9996jTkwjotgPU0UTCiDf2AG0sQMhBox3An/w4AECAwMBABYWFigsLAQAhIWFISoqSrPREa2pbRKPnKtt1RR6msxDiOHincBdXV2Rl5cHAPD29sbZs2cBAOnp6dThZURUrQX+InkNPDOfEjghhop3Ah82bBh++eUXAMDs2bOxaNEijBgxAqGhoZg0aZLGAyTaUaJiLfAXedhXJfD7ecU6i4kQwg/vUSi7du2CTCYDAMydOxcODg44efIkQkJCMHfuXI0HSLTj+UqEqr8F2rvYAABu5RTpLCZCCD8N2hNTIHhecZ86dSqmTp2q0aCI9tXXhNLexRoAcOdxEaQyBhNB7fugEkL0Q60mlAcPHvAq9NGjRw0KhuiOohOzllEonvaWMBcKIKmU4cHTEl2GRghRk1oJvHfv3pgzZw7Onz9f6zlisRhff/01unTpgujoaI0FSLRD1XZqLxIIOPg6V9XCb+UU6iwuQoj61GpCuX79Oj777DOMGjUKQqEQAQEBcHd3h7m5OZ49e4Zr164hNTUVAQEBWL9+PU3oMQKqtlOrrr2zDf56VIC0nEKM7Oyqq9AIIWpSqwbu4OCADRs2IDMzE9u3b0f79u3x5MkTpKWlAQCmTZuGixcv4tSpU5S8jURpRc3t1Kpr93dH5k3qyCTEIPHqxDQ3N8fkyZMxefJkbcVDdKS+TkzgeUdmGjWhEGKQaDXCZqq+NnDg+VDCu4+LUSmV6SQuQoj6KIE3U2X1TOQBgFYtLGBpZoJyqQz38mgkCiGGhhJ4M/V8PfDaW9EEAg7tnKkZhRBDRQm8mapvMSu5djQjkxCDxTuBFxdrbm2M7du3o1u3brC1tYWtrS369++P33//XWPlk9qVlNc/CgV43pFJY8EJMTy8E7iLiwtmzZqFkydPNvrmHh4eWLt2LS5cuIALFy5g2LBhmDBhAlJTUxtdNqmbvAnFvJaZmHLPa+CUwAkxNLwTeFRUFMRiMYYPH4727dtj7dq1yMzMbNDNQ0JCMGbMGLRv3x7t27fHp59+Cmtra8UStUQ7pDIGSWXVqJL6auAd/k7g6U+KUV5JI1EIMSS8E3hISAgOHDiAzMxMvPPOO4iKioK3tzfGjRuH6OhoVFZWNigQqVSKffv2obi4GP3791d5jkQiQUFBgdKL8CffkR6ofTVCOTc7c9iITFEpY7hHS8sSYlAa3Inp6OiIRYsW4fLly9i0aROOHj2KKVOmwN3dHStXrkRJiXrDzq5evQpra2uIRCLMnTsXMTEx6NSpk8pz16xZAzs7O8XL09OzoeE3a/L2b44DzIV1fwtwHAdfagcnxCA1OIFnZ2dj3bp16NixI5YvX44pU6bg2LFj2Lx5M2JiYjBx4kS1yunQoQNSUlJw9uxZvPPOOwgPD8e1a9dUnrtixQqIxWLFKyMjo6HhN2uKHemFJuC4+peJbe/8dzt4NiVwQgwJ7/XAo6OjERERgT/++AOdOnXC/Pnz8cYbb6BFixaKc3r06AF/f3+1yjMzM4Ovry8AICAgAElJSdiyZQt27txZ41yRSASRSMQ3ZFKNOrMwX9ROUQOnoYSEGBLeCXzmzJl49dVXcerUKfTu3VvlOW3atMG//vWvBgXEGINEImnQtUQ96qyD8iLF7jy5VAMnxJDwTuBZWVmwtLSs8xwLCwusWrWq3rL++c9/YvTo0fD09ERhYSH27duH+Ph4xMbG8g2L8PB8Mwf1vvwdXKsS+P28EkgqpRCZqpf4CSHaxbsN3MbGBrm5uTWO5+XlwcSE3w92Tk4OwsLC0KFDBwwfPhznzp1DbGwsRowYwTcswoO8E1PdGrizjQi25qaQyhjuPqaRKIQYCt41cMaYyuMSiQRmZma8yvr222/53p5oQKkaC1m9iOM4tHexwYX7z3ArpxAd3Wy1GR4hRE1qJ/Avv/wSQNUP8zfffANra2vFZ1KpFImJifDz89N8hETj+HZiAlUzMuUJnBBiGNRO4Js3bwZQVQPfsWOHUnOJmZkZfHx8sGPHDs1HSDROne3UqmtPI1EIMThq/wSnp6cDAIYOHYro6GjY29trLSiiXYq1wOtZB+VF8pEotKwsIYaDdxt4XFycNuIgOsS3ExN4nsDvPy1BWYW03kWwCCHap1YCX7x4Mf7zn//AysoKixcvrvPcTZs2aSQwoj0NaQNvaW0Ge0shnpVU4HZuEbq0stNWeIQQNamVwJOTk1FRUaH4/9qoMy2b6J+6mzm8iOM4tHOxwfn0p0jLLaQETogBUCuBv9hsQk0oxq8hnZhAVUfm+fSnuJlNHZmEGIJGb6lWUFCAgwcP4saNG5qIh+hAyQuLWfFBHZmEGBbeCXzq1KnYunUrAKC0tBQBAQGYOnUqunbtigMHDmg8QKJ5pRXqbadWXTtnWhOFEEPCO4EnJiZi0KBBAICYmBgwxpCfn48vv/wSn3zyicYDJJrHdzErOflY8IynpYqRLIQQ/eGdwMViMRwcHAAAsbGxePnll2FpaYmxY8ciLS1N4wESzWtIJyYAOFqL0NK6ajnflAf5mg6LEMIT7wTu6emJM2fOoLi4GLGxsQgODgYAPHv2DObm5hoPkGheQ4YRyr3U0RkA8MuVhu2DSgjRHN4J/P3338e0adPg4eEBd3d3BAUFAahqWunataum4yNa8LwTk/c8Lozv4Q4A+O1qNiSV0nrOJoRoE++f4Hnz5qFPnz7IyMjAiBEjIBBU/Q5o06YNtYEbidLyhnViAkDf1o5wsRUhp0CChJuPEdzZVdPhEULU1KBhhAEBAZg0aZLSioRjx47FgAEDNBYY0Q7GGEp4Lif7IhMBh/Hdq2rhP1+mZhRC9Il3DVwqlSIyMhLHjh1Dbm4uZDKZ0ufHjx/XWHBE8ySVMsiXdOc7CkVuQo9W+PpEOo5ey0FhWQVszIUajJAQoi7eCXzhwoWIjIzE2LFj0aVLF5o+b2Tk7d8AYMlzJqZcZ3dbtHGywt3HxfgzNQcv9/LQVHiEEB54/wTv27cPP/zwA8aMGaONeIiWycdvm5kKYCJo2C9fjuMwsUcrbDpyCz9fzqQEToie8G4DNzMzg6+vrzZiITpQ1oj27xfJ28FPpj3G40JJo+MihPDHO4H/4x//wJYtW2rdG5MYNsUY8Eau5+3T0grdPVtAxoDDNCacEL3g3YRy8uRJxMXF4ffff0fnzp0hFCp3YEVHR2ssOKJ5DZ1Gr8rEHu64nJGPgymZmDGgdaPLI4TwwzuBt2jRApMmTdJGLEQHnk+jb1gH5ovGdnPDf369hpSMfNzPK4a3o1WjyySEqI/3T3FERIQ24iA6oskauLONOQb4tsSJtCc4lJKJBcPbNbpMQoj6GjSRp7KyEkePHsXOnTtRWFi1tGhmZiaKimihf0NX0ohZmKpM6NEKAHAw5RH1ixCiY7wT+P3799G1a1dMmDAB8+fPx+PHjwEA69atw5IlS3iVtWbNGvTu3Rs2NjZwdnbGxIkTcfPmTb4hER5KNTQKRW5kZxeITAW487gYqZkFGimTEKIe3gl84cKFCAgIwLNnz2BhYaE4PmnSJBw7doxXWQkJCZg/fz7Onj2LI0eOoLKyEsHBwSguLuYbFlGTvAlFU7vK25gL8VJHFwDAIZpaT4hONWgUyqlTp2BmZqZ03NvbG48ePeJVVmxsrNL7iIgIODs74+LFixg8eDDf0IgaGrOUbG3G93DH4atZ+OFCBt4d5gtbmlpPiE7wroHLZDJIpTWXEX348CFsbGwaFYxYLAYAxYYRRPOer0TY+FEocsP9nNHWyQr5JRXYlXBXY+USQurGO4GPGDECX3zxheI9x3EoKirCqlWrGjW9njGGxYsXY+DAgejSpYvKcyQSCQoKCpRehJ+GbmhcF1MTAZaN8gMAfHPyLnILyjRWNiGkdrwT+ObNm5GQkIBOnTqhrKwMr7/+Onx8fPDo0SN8/vnnDQ7k3XffxZUrVxAVFVXrOWvWrIGdnZ3i5enp2eD7NVcN3U6tPsGdXNDTqwXKKmT44hhtrUeILvBO4O7u7khJScHSpUvx9ttvw9/fH2vXrkVycjKcnZ0bFMSCBQtw6NAhxMXFwcOj9oWRVqxYAbFYrHhlZGQ06H7NmTbawIGqv8Q++LsWvj8pA3ce05BSQrSNd0NoYmIiAgMDMXPmTMycOVNxvLKyEomJibw6HxljWLBgAWJiYhAfH4/Wreueji0SiSASifiGTF4g38zBQoNt4HJ92zhiuJ8zjt3IxYY/bmL7G700fg9CyHO8a+BDhw7F06dPaxwXi8UYOnQor7Lmz5+PvXv34n//+x9sbGyQnZ2N7OxslJaW8g2LqKkx26mpY9koP3Ac8Ptf2Uh+8Ewr9yCEVOGdwBljKjdxyMvLg5UVv7Uwtm/fDrFYjKCgILi5uSle+/fv5xsWUZMmp9Kr0sHVBi/3rGoGW/v7DZqdSYgWqf139OTJkwFUtXXOmDFDqSlDKpXiypUrCAwM5HVz+uHWvVINLSdbl0Uj2uPQ5UycS3+K+FuPMbRDw/pGCCF1U7sGLh/5wRiDjY2N0mgQV1dXvPXWW9i7d682YyUa8HwqvebbwOVatbDAjEAfAMDnv9+AVEa/qAnRBrV/iuWrEPr4+GDJkiW8m0uIYdB2E4rcvKC2iDr/ADeyC/HjhQy82sdLq/cjpDni3Qa+atUqSt5GTFvjwKtrYWmG94ZVLS/76eHryMynjmlCNI13As/JyUFYWBjc3d1hamoKExMTpRcxXJVSGcqlMgDaT+AAMGtga/h7tUChpBIfHLhCfR6EaBjvhtAZM2bgwYMH+PDDD+Hm5qZyRAoxTPIx4ID2m1AAwETAYcMr3TFmywmcSHuC/557gDf6eWv9voQ0Fw1ajfDEiRPo0aOHFsIh2iRvPjERcDAzadBeHry1dbLGslF++M+v1/DZb9cxuJ0TvBwtdXJvQpo63j/Fnp6e9KewkXpxR3pd/uU0M9AHfVo7oKRciqU/XYaMRqUQohG8E/gXX3yB5cuX4969e1oIh2iTfDs1cx00n7xIIOCwYUp3WJqZ4Fz6U0SevqfT+xPSVPFO4KGhoYiPj0fbtm1hY2MDBwcHpRcxXLoagaKKl6Ml/jmmIwBg3R83cJcWuyKk0Xi3gb+4FjgxLtpYC5yPaX298EdqNk6kPcGiHy7ju1l9YGdBu/cQ0lC8E3h4eLg24iA6oK2lZNXFcRw+f7kbRn6RiMsZ+Qj56iS2TeuJLq3s9BIPIcZOrSaUF3e+qb4jDu2QYzxKKzS/nRpf7i0s8L83+8HD3gIPnpZg8vbT2Hf+AXWME9IAaiVwe3t75ObmAgBatGgBe3v7Gi/5cWK4dDWNvj5dPezw64KBGO7njPJKGZZHX8WSH68o2ugJIepRqyp2/PhxRQdlXFycVgMi2qPPTszqWlia4evpAdiecAcb/7yJA5ceIjVTjC9f80d7l8Ztjk1Ic6FWAh8yZIjK/yfGRd9t4NUJBBzmD/WFv1cLvBeVjBvZhRj1RSJe6eWJ90e0g5udhb5DJMSg6WY6HjEIz0eh6K8NXJXAti1x+L1BCO7kAhkD9l/IQND6eKz5/TrEJRX6Do8Qg0UJvBkpqzCsGviLXGzNsWt6AA680x99fBwgqZRhZ8JdDFp3HNvj76BYUqnvEAkxOJTAmxH5TEx9d2LWpZe3A/a/3Q+7ZwSgg4sNCsoq8XnsDQz8/Di+PJZGNXJCXkAJvBkxtDbw2nAch2F+Lvht4SBsfKU7fBwt8aykApuO3ELg2mNY8/t1PC6U6DtMQvSu0Qk8ISEBv/32G549ox3IDZ0hjUJRh4mAw8u9PHB08RB8+Zo//FxtUFwuxc6Euxj4+XGsPpSKJ0WUyEnzpXYCX79+PVatWqV4zxjDqFGjMHToUIwbNw4dO3ZEamqqVoIkmvF8HLhhdWLWx9REgPHd3fH7wkH4NjwA/l4tIKmUIfL0PQxZF4cvjt5CEbWRk2ZI7QQeFRWFTp06Kd7/9NNPSExMxIkTJ/DkyRMEBATgo48+0kqQRDPkbeDa3JFemziOw/COLoh+JxD/e7MvunvYobhcii+OpiFofRy+O3MP5ZUyfYdJiM6oncDT09PRrVs3xfvffvsNL7/8MgYMGAAHBwf8+9//xpkzZ7QSJNGMTHEZAMDZVqTnSBqH4zgE+rbEwfkDsG1aT7RuaYUnReVY+XMqRmxOwM8pjyClNcdJM6B2Aq+oqIBI9PwH/8yZMwgMDFS8d3d3x5MnTzQbHdGYwrIKRcefT8umsSk1x3EY09UNfy4ajP9M7IKW1iLczyvBwn0pGL0lEbF/ZdEaK6RJUzuB+/r6IjExEQDw4MED3Lp1S2lW5sOHD+Ho6Kj5CIlG3M8rAQC0tDaDrXnTWsJVaCJAWD9vJCwNwpLg9rA1N8WtnCLM3XsJ4746ieM3ciiRkyZJ7QT+zjvv4N1338Xs2bMxevRo9O/fX6lN/Pjx4/D39+d188TERISEhMDd3R0cx+HgwYO8rifqS39SDADwcWwatW9VrESmeHdYO5z4YBjeG+YLKzMTpGYWYFbkBby8/TQuZ+TrO0RCNErtBP72229jy5YtePr0KQYPHowDBw4ofZ6ZmYlZs2bxunlxcTG6d++OrVu38rqO8HdPnsCbSPNJXewshFgc3AEnPhiGt4e0gblQgEsP8jFx2ymsiL6Cp8Xl+g6REI3gmIb+tiwuLsbFixcxePDghgXCcYiJicHEiRPVvqagoAB2dnYQi8WwtbVt0H2bi8U/pCD60iMsHdkB84f66jscncotKMPa2BuIvvQIANDCUoilIzvg1d5eMBHobnNnQtTBJ69pbCbm7du3MXToUE0Vp5JEIqENJBroXjNoQqmNs605Nk3tgR/n9oefqw3ySyrwr5i/MGnbKVx5mK/v8AhpMKOaSr9mzRrY2dkpXp6envoOyWgo2sBbWuo5Ev3p7eOAXxcMxKqQTrARmeLKQzEmbzuNb07cpU5OYpSMKoGvWLECYrFY8crIyNB3SEZBXFKBZ38vAtUca+AvMjURYOaA1ji2ZAjGdHVFpYzhk8PX8fb3F2mhLGJ0jCqBi0Qi2NraKr1I/dLzqmrfLrYiWImMaxq9tjjbmOP/Xu+J/0zoDDMTAf68loOxX52gkSrEqKj903zo0KE6P09PT290MEQ7mnP7d104jkNYfx/08LTHvP9dRMbTUkzZcRr/GtMR4YE+4Djq4CSGTe0Ers7oEL7f8EVFRbh9+7bifXp6OlJSUuDg4AAvLy9eZZHaydu/WzeDIYQNUbXJ8iAs++ky/kjNwepfruFGdiE+ndSVRqkQg6Z2E4pMJqv3JZXy21X8woUL8Pf3V0wAWrx4Mfz9/bFy5Up+/wpSp3t5zWcMeEPZWQix441e+HBcJwg4YF9SBt7bl0yLYxGDptcG0aCgIOr91wFqQlEPx3GYPbA13O3M8d6+ZBy+koViSSW2T+tl0LsYkeaLdydmXl6e4v8zMjKwcuVKLF26VLFOCjEsjDHcpSYUXkZ3dcM34b1hLhQg/uZjhEecR2EZjVAhhkftBH716lX4+PjA2dkZfn5+SElJQe/evbF582bs2rULw4YNo7VMDNDT4nIUllWtA+7t2HzHgPM1pL0Tvp/dFzYiU5xPf4pp35yjKfjE4KidwJctW4auXbsiISEBQUFBGDduHMaMGQOxWIxnz57h7bffxtq1a7UZK2kAeft3qxYWMDfSjRz0pbePA6Le6gcHKzNceShG6M4ztBcnMShqJ/CkpCR8+umnGDhwIDZs2IDMzEzMmzcPAoEAAoEACxYswI0bN7QZK2mA9CdVy8g25xmYjdGllR1+eLs/XG3NkZZbhDeoJk4MiNoJ/OnTp3B1dQUAWFtbw8rKCg4ODorP7e3tUVhYqPkISaNQB2bj+TpbY99b/eBiK8LNnEK88c055JdQEif6x6sTs/o4b5roYPjkszCpA7NxfFpa4b9v9kNLaxGuZRVg+u7zKKCOTaJnvIYRzpgxQ7GtWllZGebOnQsrq6rEIJFQ26Ahohq45vg6W+O/b/bFq7vO4MpDMWbsPo/vZveFNS1PQPRE7Rp4eHg4nJ2dFSsBvvHGG3B3d1e8d3Z2xvTp07UZK+GJMfbCKoSUwDWhg6sN9r7ZF3YWQlx6kI9ZkUkoKa/Ud1ikmVK76hAREaHNOIgWPC6UoKRcCgEHeDlQJ6amdHa3w/ez+2Da1+dwPv0p3v7+InbP6A2hiVGtDUeaAPqOa8Lkte9W9hYwM6UvtSZ182iByFl9YCE0wYm0J1h9KJVmFROdo5/qJuyeogPTWs+RNE29vO2x5dUe4Djgv+ceYPepe/oOiTQzlMCbMPkY8NY0A1Nrgju7YsVoPwDAJ4ev4dj1HD1HRJoTSuBNWHPaiV6f5gxqg9f6eIIxYEFUMq5l0l6tRDcogTdhtIysbnAch48ndMEAX0eUlEvx5p4k5BaU6Tss0gxQAm+iZLLnQwhb0xhwrROaCLDt9V5o42SFTHEZ5nx3AaXl/NbHJ4QvSuBNVHZBGSSVMpgKOHjYW+g7nGbBzlKIiBm9YW8pxOWHYiz58TJkMhqZQrSHEngTJW//9nSwhCmNT9YZb0cr7AwLgNCEw+GrWdhyLE3fIZEmjH6ymyj5Gig+NAJF5/q0dsCnE7sCALYcS8Ohy5l6jog0VZTAm6h7T2gMuD5N7e2Jtwa3AQAs/fEyUjLy9RsQaZIogTdRijHgtA643nwwyg/D/ZwhqZRhzncXkCUu1XdIpImhBN5E0RBC/TMRcNjymj86uNjgcaEEb+65QAtfEY2iBN4ESWUMD/L+3omHhhDqlbXIFN+EB8DRygypmQVYuC8FkkoaXkg0gxJ4E5SZX4pyqQxmJgK4t6AhhPrm6WCJnWG9YGYiwJFrOQjdeZaaU4hGUAJvguQTeLwcLWEioF2TDEGAjwO+CQ+AnYUQKRn5CPnqJM7dzdN3WMTIUQJvYlIy8rHm96rNpWkbNcMyuL0Tfnl3IDq62eJJUTmmfXMOEafSaRla0mB6T+Dbtm1D69atYW5ujl69euHEiRP6DskoFZRVYOXPf2HStlO4nlUAOwsh5g5po++wSDVejpaIficQE3q4o1LG8NEv17BofwqKJdS5SfjjmB5//e/fvx9hYWHYtm0bBgwYgJ07d+Kbb77BtWvX4OXlVe/1BQUFsLOzg1gshq2trQ4iNjyMMRy+moWPf7mG3MKqfUkn+bfCv8Z2REtrkZ6jI7VhjCHi1D18+tt1SGUMIlMBBvi2xPCOzhjm5ww3O+q7aK745DW9JvC+ffuiZ8+e2L59u+JYx44dMXHiRKxZs6be6xuawAvKKgx6jQrGgEoZQ6VMhopKhgqZDJVShpLySmSLy5ApLkNWfimyxGW4+6QY17Oqli9t3dIKn0zsggG+LfX8LyDqOns3D8sPXMG9v0cNyXV2t8Wgdk5wshHB1twUthZC2JoLYWthCnOhCUw4DiYCDgIBBxOOg0AAcKjq7+A4gEPVKokAoKoXhKOuEZ2xNRdCwKMvik9e09t22uXl5bh48SKWL1+udDw4OBinT59WeY1EIoFEIlG8Lyho2LrL074+h6uPxA261hCZmQjwTlBbvBPUFuZCE32HQ3jo18YRcUuCcCunCEev5+DY9RwkZ+QjNbMAqbSueJNwdsVwuNqZa6VsvSXwJ0+eQCqVwsXFRem4i4sLsrOzVV6zZs0afPTRR7oIzyAITTiYCgQwNeFgZiKAyFQAFztzuNtZwM3OHG4tLOBuZ47uni1ouKAR4zgOHVxt0MHVBvOH+iKvSIK4m4+RkvEMBaWVKCirQEFpBcSlFSgoq0RZhRSMVY33lzIG2d//Bar+eiPNh94SuBxX7W85xliNY3IrVqzA4sWLFe8LCgrg6enJ+54x8wJ5X6NrJgKu1udAmjZHaxGm9PLAlF4eDS6D1ZHQa8vxNBpGO7Q5lFdvCbxly5YwMTGpUdvOzc2tUSuXE4lEEIka3zFHy6uSpk7R/s0rd1CFwdjoLZOZmZmhV69eOHLkiNLxI0eOIDDQ8GvIhBCib3ptQlm8eDHCwsIQEBCA/v37Y9euXXjw4AHmzp2rz7AIIcQo6DWBh4aGIi8vDx9//DGysrLQpUsX/Pbbb/D29tZnWIQQYhT0Og68sWgiDyGkqeGT16g3jxBCjJTehxE2hvyPh4ZO6CGEEEMjz2fqNI4YdQIvLCwEgAaNBSeEEENWWFgIOzu7Os8x6jZwmUyGzMxM2NjYKMa9yif3ZGRkGF27OMWuHxS7flDsqjHGUFhYCHd3dwgEdbdyG3UNXCAQwMND9Ww1W1tbo/umkKPY9YNi1w+Kvab6at5y1IlJCCFGihI4IYQYqSaXwEUiEVatWqWRNVN0jWLXD4pdPyj2xjPqTkxCCGnOmlwNnBBCmgtK4IQQYqQogRNCiJGiBE4IIUbKqBL4mjVr0Lt3b9jY2MDZ2RkTJ07EzZs367wmPj4eHMfVeN24cUNHUVdZvXp1jRhcXV3rvCYhIQG9evWCubk52rRpgx07dugoWmU+Pj4qn+H8+fNVnq/PZ56YmIiQkBC4u7uD4zgcPHhQ6XPGGFavXg13d3dYWFggKCgIqamp9ZZ74MABdOrUCSKRCJ06dUJMTIxOY6+oqMAHH3yArl27wsrKCu7u7pg+fToyMzPrLDMyMlLl16KsrExnsQPAjBkzasTQr1+/esvV93MHoPL5cRyH9evX11qmrp67USXwhIQEzJ8/H2fPnsWRI0dQWVmJ4OBgFBcX13vtzZs3kZWVpXi1a9dOBxEr69y5s1IMV69erfXc9PR0jBkzBoMGDUJycjL++c9/4r333sOBAwd0GHGVpKQkpbjluyi98sordV6nj2deXFyM7t27Y+vWrSo/X7duHTZt2oStW7ciKSkJrq6uGDFihGJdHVXOnDmD0NBQhIWF4fLlywgLC8PUqVNx7tw5ncVeUlKCS5cu4cMPP8SlS5cQHR2NW7duYfz48fWWa2trq/R1yMrKgrm5ZndJr++5A8CoUaOUYvjtt9/qLNMQnjuAGs9u9+7d4DgOL7/8cp3l6uK5gxmx3NxcBoAlJCTUek5cXBwDwJ49e6a7wFRYtWoV6969u9rnL1u2jPn5+Skde/vtt1m/fv00HBl/CxcuZG3btmUymUzl54byzAGwmJgYxXuZTMZcXV3Z2rVrFcfKysqYnZ0d27FjR63lTJ06lY0aNUrp2MiRI9mrr76q8Zjlqseuyvnz5xkAdv/+/VrPiYiIYHZ2dpoNrh6qYg8PD2cTJkzgVY6hPvcJEyawYcOG1XmOrp67UdXAqxOLxQAABweHes/19/eHm5sbhg8fjri4OG2HplJaWhrc3d3RunVrvPrqq7h7926t5545cwbBwcFKx0aOHIkLFy6goqJC26HWqry8HHv37sWsWbMUC4jVxhCe+YvS09ORnZ2t9FxFIhGGDBmC06dP13pdbV+Luq7RBbFYDI7j0KJFizrPKyoqgre3Nzw8PDBu3DgkJyfrJsBq4uPj4ezsjPbt22POnDnIzc2t83xDfO45OTk4fPgwZs+eXe+5unjuRpvAGWNYvHgxBg4ciC5dutR6npubG3bt2oUDBw4gOjoaHTp0wPDhw5GYmKjDaIG+ffviu+++wx9//IGvv/4a2dnZCAwMRF5ensrzs7Oz4eLionTMxcUFlZWVePLkiS5CVungwYPIz8/HjBkzaj3HUJ55ddnZ2QCg8rnKP6vtOr7XaFtZWRmWL1+O119/vc7FlPz8/BAZGYlDhw4hKioK5ubmGDBgANLS0nQYLTB69Gj897//xfHjx7Fx40YkJSVh2LBhkEgktV5jiM99z549sLGxweTJk+s8T2fPXet1fC2ZN28e8/b2ZhkZGbyvHTduHAsJCdFCVOorKipiLi4ubOPGjSo/b9euHfvss8+Ujp08eZIBYFlZWboIUaXg4GA2btw43tfp45mj2p/Dp06dYgBYZmam0nlvvvkmGzlyZK3lCIVC9r///U/p2N69e5lIJNJovC+qHvuLysvL2YQJE5i/vz8Ti8W8ypVKpax79+5swYIFGohStbpil8vMzGRCoZAdOHCg1nMM7bkzxliHDh3Yu+++y7tcbT13o6yBL1iwAIcOHUJcXFyty8nWpV+/fjqvgVRnZWWFrl271hqHq6trjZpGbm4uTE1N4ejoqIsQa7h//z6OHj2KN998k/e1hvDM5aN+VD3X6jW96tfxvUZbKioqMHXqVKSnp+PIkSO8lzIVCATo3bu33r8Wbm5u8Pb2rjMOQ3ruAHDixAncvHmzQd//2nruRpXAGWN49913ER0djePHj6N169YNKic5ORlubm4ajo4fiUSC69ev1xpH//79FaM95P78808EBARAKBTqIsQaIiIi4OzsjLFjx/K+1hCeeevWreHq6qr0XMvLy5GQkIDAwMBar6vta1HXNdogT95paWk4evRog36RM8aQkpKi969FXl4eMjIy6ozDUJ673LfffotevXqhe/fuvK/V2nPXaH1ey9555x1mZ2fH4uPjWVZWluJVUlKiOGf58uUsLCxM8X7z5s0sJiaG3bp1i/31119s+fLlDECdf7ppwz/+8Q8WHx/P7t69y86ePcvGjRvHbGxs2L1791TGfffuXWZpackWLVrErl27xr799lsmFArZTz/9pNO45aRSKfPy8mIffPBBjc8M6ZkXFhay5ORklpyczACwTZs2seTkZMVIjbVr1zI7OzsWHR3Nrl69yl577TXm5ubGCgoKFGWEhYWx5cuXK96fOnWKmZiYsLVr17Lr16+ztWvXMlNTU3b27FmdxV5RUcHGjx/PPDw8WEpKitL3v0QiqTX21atXs9jYWHbnzh2WnJzMZs6cyUxNTdm5c+d0FnthYSH7xz/+wU6fPs3S09NZXFwc69+/P2vVqpXBP3c5sVjMLC0t2fbt21WWoa/nblQJHIDKV0REhOKc8PBwNmTIEMX7zz//nLVt25aZm5sze3t7NnDgQHb48GGdxx4aGsrc3NyYUChk7u7ubPLkySw1NbXWuBljLD4+nvn7+zMzMzPm4+NT6zePLvzxxx8MALt582aNzwzpmcuHMFZ/hYeHM8aqhhKuWrWKubq6MpFIxAYPHsyuXr2qVMaQIUMU58v9+OOPrEOHDkwoFDI/Pz+t/DKqK/b09PRav//j4uJqjf39999nXl5ezMzMjDk5ObHg4GB2+vRpncZeUlLCgoODmZOTExMKhczLy4uFh4ezBw8eKJVhiM9dbufOnczCwoLl5+erLENfz52WkyWEECNlVG3ghBBCnqMETgghRooSOCGEGClK4IQQYqQogRNCiJGiBE4IIUaKEjghhBgpSuBEYxhjeOutt+Dg4ACO45CSkqLvkAxWeXk5fH19cerUKb3FsGTJErz33nt6uz9pPErgRGNiY2MRGRmJX3/9FVlZWXUu89vc7dq1C97e3hgwYIDS8bi4OIwbNw5OTk4wNzdH27ZtERoaqrQUr3zLuvz8/Brl+vj44IsvvlArhmXLliEiIgLp6emN+acQPaIETjTmzp07cHNzQ2BgIFxdXWFqalrjnPLycj1EZni++uqrGqvabdu2DcOHD4ejoyP279+P69ev4/vvv0dgYCAWLVqk8RicnZ0RHByst71WiQZofHI+aZbCw8OV1pHw9vZmjFWtETF//ny2aNEi5ujoyAYPHswYYyw1NZWNHj2aWVlZMWdnZ/bGG2+wx48fK8orKipiYWFhzMrKirm6urINGzawIUOGsIULFyrOgYq1m+3s7JTWxnn48CGbOnUqa9GiBXNwcGDjx49n6enpSnFPmDCBrV+/nrm6ujIHBwc2b948Vl5erjinrKyMLV26lHl4eDAzMzPm6+vLvvnmGyaTyVjbtm3Z+vXrlWK4evUq4ziO3b59W+WzunjxIhMIBEpred+/f58JhUK2aNEilde8uH1dXVvWeXt7s82bNzPGqrb1goo1PlatWqU4PzIyknl6eqq8JzF8VAMnGrFlyxZ8/PHH8PDwQFZWFpKSkhSf7dmzB6ampjh16hR27tyJrKwsDBkyBD169MCFCxcQGxuLnJwcTJ06VXHN0qVLERcXh5iYGPz555+Ij4/HxYsXecVUUlKCoUOHwtraGomJiTh58iSsra0xatQopb8E4uLicOfOHcTFxWHPnj2IjIxEZGSk4vPp06dj3759+PLLL3H9+nXs2LED1tbW4DgOs2bNQkREhNJ9d+/ejUGDBqFt27Yq40pMTET79u2V1vI+cOAAKioqsGzZMpXX1Ld9nSqhoaFKG+pGRUXB1NRUqdmmT58+yMjIwP3793mXTwyAvn+DkKZj8+bNipq33JAhQ1iPHj2Ujn344YcsODhY6VhGRoZitcPCwkJmZmbG9u3bp/g8Ly+PWVhY8KqBf/vtt6xDhw5KtVeJRMIsLCzYH3/8wRirqoF7e3uzyspKxTmvvPIKCw0NZYwxdvPmTQaAHTlyROW/OTMzk5mYmCiWCS0vL2dOTk4sMjKylqdUtSl09U1x586dy2xtbZWO/fTTT8zKykrxunLlCmPseQ38xc/kL47jFDXwF92+fZs5OjqydevWKR0Xi8UMAIuPj681XmK4ajZSEqJhAQEBSu8vXryIuLg4WFtb1zj3zp07KC0tRXl5Ofr376847uDggA4dOvC678WLF3H79m3Y2NgoHS8rK8OdO3cU7zt37gwTExPFezc3N1y9ehUAkJKSAhMTEwwZMkTlPdzc3DB27Fjs3r0bffr0wa+//oqysjK88sortcZVWloKc3PzGser17JHjhyJlJQUPHr0CEFBQZBKpUqfnzhxosa/LSgoqEa5YrEY48aNw+jRo7F06VKlzywsLABU/bVCjA8lcKJ1VlZWSu9lMhlCQkLw+eef1zjXzc1N7W2nOI4Dq7YackVFhdJ9evXqhf/+9781rnVyclL8f/UdjjiOg0wmA/A8wdXlzTffRFhYGDZv3oyIiAiEhobC0tKy1vNbtmyp+AUh165dO4jFYmRnZyu2frO2toavr6/KzmCgaoeh6jvSVz9XKpUiNDQUtra2+Prrr2uU8fTpUwDKz4MYD2oDJzrXs2dPpKamwsfHB76+vkovKysr+Pr6QigU4uzZs4prnj17hlu3bimV4+TkhKysLMX7tLQ0pZpkz549kZaWBmdn5xr3sbOzUyvWrl27QiaTISEhodZzxowZAysrK2zfvh2///47Zs2aVWeZ/v7+uHHjhtIvnylTpkAoFKr8pdYYixYtwtWrVxETE6Oy1v/XX39BKBSic+fOGr0v0Q1K4ETn5s+fj6dPn+K1117D+fPncffuXfz555+YNWsWpFIprK2tMXv2bCxduhTHjh3DX3/9hRkzZkAgUP52HTZsGLZu3YpLly7hwoULmDt3rlJtetq0aWjZsiUmTJiAEydOID09HQkJCVi4cCEePnyoVqw+Pj4IDw/HrFmzcPDgQaSnpyM+Ph4//PCD4hwTExPMmDEDK1asgK+vr1LTjypDhw5FcXExUlNTFce8vLywceNGbNmyBeHh4YiLi8O9e/dw6dIlfPnll4r78BEREYFt27Zhx44dEAgEyM7ORnZ2NoqKihTnnDhxAoMGDVLrLw1ieCiBE51zd3fHqVOnIJVKMXLkSHTp0gULFy6EnZ2dIkmvX78egwcPxvjx4/HSSy9h4MCB6NWrl1I5GzduhKenJwYPHozXX38dS5YsUWq6sLS0RGJiIry8vDB58mR07NgRs2bNQmlpKa/d3Ldv344pU6Zg3rx58PPzw5w5c1BcXKx0zuzZs1FeXl5v7RsAHB0dMXny5BpNOwsWLMCff/6Jx48fY8qUKWjXrh3GjBmD9PR0xMbGomvXrmrHDAAJCQmQSqUYP3483NzcFK8NGzYozomKisKcOXN4lUsMB22pRoxGUFAQevToofZMQ106deoUgoKC8PDhQ7i4uNR7/tWrV/HSSy+p7GTVlcOHD2Pp0qW4cuVKre3sxLBRDZyQRpBIJLh9+zY+/PBDTJ06Va3kDVS1ra9btw737t3TboB1KC4uRkREBCVvI0ZfOUIaISoqCrNnz0aPHj3w/fff87o2PDxcS1Gp58WJU8Q4URMKIYQYKWpCIYQQI0UJnBBCjBQlcEIIMVKUwAkhxEhRAieEECNFCZwQQowUJXBCCDFSlMAJIcRIUQInhBAj9f9pixCYTpLInQAAAABJRU5ErkJggg==",
      "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(Ey*np.conj(Ey)), label='Slavin-Kalinikos model')\n",
    "plt.xlabel('frequency (GHz)')\n",
    "plt.ylabel('BLS intensity (a.u.)')\n",
    "plt.title('BLS signal')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}