{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Effective Vertical Anisotropy of Layered Aquifers \n",
    "\n",
    "_Mark Bakker and Bram Bot._\n",
    "\n",
    "Notebook to run experiments presented in the paper \"The Effective Vertical Anisotropy of\n",
    "Layered Aquifers.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference: M. Bakker and B. Bot (2024) The effective vertical anisotropy of layered aquifers. Groundwater. Available online early: [doi](https://doi.org/10.1111/gwat.13432)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import brentq\n",
    "\n",
    "import timml as tml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to generate hydraulic conductivities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generatek(ksection=20 * [0.1, 1], nsections=10, seed=1):\n",
    "    \"\"\"Generate k.\n",
    "\n",
    "    Input:\n",
    "    ksection: k values in the section\n",
    "    nsection: number of sections\n",
    "    seed: seed of random number generator\n",
    "    \"\"\"\n",
    "    nk = len(ksection)\n",
    "    # nlayers = nk * nsections\n",
    "    kaq = np.zeros((nsections, nk))\n",
    "    rng = np.random.default_rng(seed)\n",
    "    for i in range(nsections):\n",
    "        kaq[i] = rng.choice(ksection, nk, replace=False)\n",
    "    return kaq.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to create a model with a canal given `kx` and `kz`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makemodel(kx, kz, d=4, returnmodel=False):\n",
    "    \"\"\"Creates model with river at center, and water supplied from infinitiy.\n",
    "\n",
    "    d is depth of river.\n",
    "    \"\"\"\n",
    "    H = 40  # thickness of model\n",
    "    # d = 4 # depth of river\n",
    "    naq = len(kx)\n",
    "    ml = tml.Model3D(kaq=kx, z=np.linspace(H, 0, naq + 1), kzoverkh=kz / kx)\n",
    "    tml.LineSink1D(ml, xls=0, sigls=2, layers=np.arange(int(d * 10)))\n",
    "    tml.Constant(ml, xr=2000, yr=0, hr=0, layer=0)\n",
    "    ml.solve(silent=True)\n",
    "    if returnmodel:\n",
    "        return ml\n",
    "    return ml.head(0, 0)[0]\n",
    "\n",
    "\n",
    "def func(kz, kx, h0, d=4, nlayers=400):\n",
    "    \"\"\"Computes head difference.\"\"\"\n",
    "    hnew = makemodel(kx * np.ones(nlayers), kz * np.ones(nlayers), d, returnmodel=False)\n",
    "    return hnew - h0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to create a model with a partially penetrating well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makemodelradial(kx, kz, d=4, rw=0.1, returnmodel=False):\n",
    "    \"\"\"Creates model with river at center, and water supplied from infinitiy.\"\"\"\n",
    "    H = 40  # thickness of model\n",
    "    # d = 4 # depth of river\n",
    "    Qw = 1000\n",
    "    naq = len(kx)\n",
    "    ml = tml.Model3D(kaq=kx, z=np.linspace(H, 0, naq + 1), kzoverkh=kz / kx)\n",
    "    tml.Well(ml, xw=0, yw=0, Qw=Qw, rw=rw, layers=np.arange(int(d * 10)))\n",
    "    tml.Constant(ml, xr=2000, yr=0, hr=0, layer=0)\n",
    "    ml.solve(silent=True)\n",
    "    if returnmodel:\n",
    "        return ml\n",
    "    return ml.head(0, 0)[0]\n",
    "\n",
    "\n",
    "def funcradial(kz, kx, h0, d=4, rw=0.1, nlayers=400):\n",
    "    \"\"\"Computes head difference.\"\"\"\n",
    "    hnew = makemodelradial(\n",
    "        kx * np.ones(nlayers), kz * np.ones(nlayers), d, rw, returnmodel=False\n",
    "    )\n",
    "    return hnew - h0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find effective vertical hydraulic conductivity for one realization of 400 layers and time it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.6 s ± 226 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "kaq = 80 * [1, 3.16, 10, 31.6, 100]  # 400 k values\n",
    "k = generatek(ksection=kaq, nsections=1)  # random order of k values\n",
    "h0 = makemodel(k, k)  # head at canal\n",
    "kx = np.mean(k)  # equivalent horizontal k\n",
    "# vertical hydraulic conductivity:\n",
    "%timeit kz = brentq(func, a=0.001 * kx, b=kx, args=(kx, h0)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the experiment for Figure 3a. This is commented out because it takes a long time to run. The number of the realization is printed to the screen every 10 realizations. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So 1000 realizations takes on the order of 3000 seconds (on this machine), so around 50 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 960 970 980 990 \n",
      " completed\n"
     ]
    }
   ],
   "source": [
    "kaq = np.array(80 * [1, 3.16, 10, 31.6, 100])\n",
    "ntot = 1000\n",
    "aniso = np.zeros(ntot)\n",
    "\n",
    "for i in range(ntot):\n",
    "    k = generatek(kaq, nsections=1, seed=i)\n",
    "    h0 = makemodel(k, k)\n",
    "    kx = np.mean(k)\n",
    "    kz = brentq(func, a=0.001 * kx, b=kx, args=(kx, h0))\n",
    "    aniso[i] = kx / kz\n",
    "    if i % 10 == 0:\n",
    "        print(i, end=\" \")\n",
    "print(\"\\n completed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create figure 3a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import lognorm\n",
    "\n",
    "\n",
    "def create_fig3a():\n",
    "    plt.figure(figsize=(3, 3))\n",
    "    plt.subplot(211)\n",
    "    plt.hist(aniso, bins=np.arange(2, 20, 0.5), density=True)\n",
    "    p5, p50, p95 = np.percentile(aniso, [5, 50, 95])\n",
    "    # print('p5, p50, p95', p5, p50, p95)\n",
    "    plt.axvline(p5, color=\"C1\")\n",
    "    plt.axvline(p95, color=\"C1\")\n",
    "    plt.axvline(p50, color=\"C2\")\n",
    "    kheq = np.mean(kaq)\n",
    "    kveq = len(kaq) / np.sum(1 / kaq)\n",
    "    anisoeq = kheq / kveq\n",
    "    plt.axvline(anisoeq, color=\"k\", linestyle=\":\", linewidth=1)\n",
    "    #\n",
    "    shape, loc, scale = lognorm.fit(aniso)\n",
    "    # print('shape, loc, scale: ', shape, loc, scale)\n",
    "    x = np.linspace(0, 20, 100)\n",
    "    pdf1 = lognorm.pdf(x, shape, loc, scale)\n",
    "    plt.plot(x, pdf1, \"k--\", lw=1)\n",
    "    plt.xlim(0, 20)\n",
    "    plt.ylim(0, 0.25)\n",
    "    plt.xticks(np.arange(0, 21, 4))\n",
    "    plt.xlabel(r\"$\\alpha_{eff}$ for $d=4$ m\")\n",
    "    plt.ylabel(\"pdf\")\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASEAAACqCAYAAADiHfm0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnUElEQVR4nO3deVzU1foH8M/sA8giO4OAKCiIJgKauGGldNU0bVHT1HIpQgPETM28laYY3VwKcQm7divNqz/yalmK5VJpuSCIgKKIQAKyKDuzn98fXOdGICLMzHcYnvfrNa/gzJlznhmZp/Ndzjk8xhgDIYRwhM91AISQro2SECGEU5SECCGcoiRECOEUJSFCCKcoCRFCOEVJiBDCKSHXARibVqtFUVERrK2twePxuA6HkE6JMYaamhrIZDLw+R0by3S5JFRUVAQPDw+uwyDELBQWFqJHjx4daqPLJSFra2sAjR+ejY0Nx9EYkLIO+Khv489LrgJiK6N1Xa+qx6ido1B1tgopa1LQ072n0fruMA4/t86kuroaHh4euu9TR3S5JHTvEMzGxsbMk5AAkPz3cNPGxqhfJqFKCImLBM4TneHj4wNLkaXR+u4wDj+3zkgfpzToxDQxCE2DBrWZtaiuruY6FGLiKAkRg1DeVuLmhzeRm5vLdSjExFESIgYhcZegz4d90K9fP65DISaOkhAxCL6ID7GTGBKJhOtQiImjJEQMQlmhRNG/ilBYUMh1KMTEURIiBqFVaFGfW4+6ujquQyEmjpIQMQipTAqf93zg5+/HdSjExFESIoRwipIQMQh5oRxXoq8g41IG16EQE0dJiBiEwFoA+yfs4ejkyHUoxMRREiIGIbITwXmSM9zc3LgOhZg4SkLEIDRyDeqv16O2tpbrUIiJoyREDEJZosSN92/g2rVrXIdCTBwlIWIQEpkEPmt94OdHl+hJ6zhPQomJifD29oZUKkVwcDB+/vnn+9ZNTk7G2LFj4eTkBBsbG4SGhuLIkSNGjJa0FV/Mh9RdCgsLC65DISaO0/WE9u7di5iYGCQmJmL48OHYvn07xo0bh6ysLHh6ejarf+rUKYwdOxbr1q2DnZ0d/vnPf2LixIn4/fffMWjQIA7egXnrufy7ZmU3109o02tVd1SoOFqBW6NuwdfbV9+hETPC6Uhow4YNmDdvHubPnw9/f39s2rQJHh4e2Lp1a4v1N23ahDfffBODBw+Gr68v1q1bB19fXxw6dMjIkZMH0TRoUJ1WjaqqKq5DISaOsySkVCpx4cIFhIeHNykPDw/H6dOn29SGVqtFTU0N7O3t71tHoVCgurq6yYMYntRdij7r+6BfAC3lQVrHWRIqLy+HRqOBi4tLk3IXFxeUlJS0qY2PPvoIdXV1mDp16n3rxMXFwdbWVvegRe4JMS2cn5j+6xq1jLE2rVu7Z88evPvuu9i7dy+cnZ3vW2/FihWoqqrSPQoLaWkJY5DfkiPnzRxkZWZxHQoxcZydmHZ0dIRAIGg26iktLW02OvqrvXv3Yt68edi3bx/GjBnTal2JREILa3FAYCmAzWAb2NnZcR0KMXGcjYTEYjGCg4ORkpLSpDwlJQXDhg277+v27NmDl156Cbt378aECW27UkOMT9RdBNfnXSFzl3EdCjFxnB6OxcbGIikpCZ999hmys7OxePFiFBQUICIiAkDjodTs2bN19ffs2YPZs2fjo48+wtChQ1FSUoKSkhK6AmNk58+fx8mTJwEAZWVliIiIwOeffw6lUqmro1VqIS+Uo6GhgaswSSfBaRKaNm0aNm3ahNWrVyMwMBCnTp3C4cOH4eXlBQAoLi5GQUGBrv727duhVquxcOFCuLm56R7R0dFcvYUuRVNfhXnz5mHw4MH45JNPADReYDhz5gxefvll+Pj4YPv27WCMQVGkwPVV13HlyhWOoyamjvPNDyMjIxEZGdnic7t27Wry+4kTJwwfEGmR6s4t3N6zAv8n0GDr1q1YsGABAMDf3x/p6enIzMzEunXrEBERgcwrmZD4S9Dr773Qp08fjiMnpo7zJERMH2MM5d9uAE9siayM3yGTNT3Po7uz2mMGnKZ4Y3+NExwkx2HZyxJWVrSDKWkdJSHyQDweD44Tl4AnkmLYxxcBXLxvXcs+oQBPCXmxHIVbCnEp4BKGhgw1XrCk0+H8PiFi2qrP/wea+iqIussg7Hb/O9P/SlWmguKWAjOmzcCdO3cMGCHp7CgJkfuqyzqJuz9+Cnn+pYd+rfUj1vBd54uqyiq88MILUKvVBoiQmANKQqRFt2/fxp2jibD0HwUr/5HtakPiJsHnX32OY8eO4a233tJzhMRcUBIiLVq9ejUAwH5sRLteL78lx7W3r+GV/7sB28fmI+lidYtLgxBCJ6ZJM+Xl5di5cydsQmdAYGHTrjb4Uj6s/KzAF1vAJmSSrlylUkEkEukrVGIGaCREmnF0dERqaiqsg55qdxtiBzFkL8ogtPnflj9VZ/6NJ598EhqNRh9hEjNBSYg0UVFRAbVajX79+oEvav/EX61SC8VtBZj6f1M5JDI/nDhxAhs3btRHqMRMUBIiTSxcuBB/+9vfOtyOokiBa8uuQVn+v6VTpF6PIDo6GqtWraJdOIgOJSGic/PmTezfvx9TpkzpcFtiFzF6LusJUfemd1e///77kMlkmD9/PrRabYf7IZ0fJSGis3nzZtjY2OCll17qcFsCCwG6+XcDX9J0t42ANSdQN2Q+Uiul6Ln0AF0xI3R1jDSqrKxEUlISoqOj9TLfS12lxt1f70LkUgmBZdOVL6Vej0Dq9UiH+yDmgUZCBABw48YN9O7dG4sWLdJLe6oqFcq/K4em9u5969RmHkfFkS166Y90XpSECAAgKCgIaWlpcHV11Ut7Fp4W8N/iD7GL9/0rMS1q077H0aNH9dIn6ZwoCRFcu3YNly9fNnq/VgGPQ+I5AIsWLYJCoTB6/8Q0tDkJ2dvbo7y8HAAwd+5c1NTUGCwoYlxxcXGYNGkSGGN6a1NRrEDumlyoKm7dtw6Px4P9mAjk5eVhw4YNeuubdC5tTkJKpVK3ceDnn38OuVxusKCI8VRXV+t2L2nLVkttxRPxIHWXgidsfYqG2MkLUVFROHbsmF6TIOk82nx1LDQ0FJMnT0ZwcDAYY4iKioKFhUWLdT/77DO9BUgMa8+ePZDL5Xq5LP9nYkcx3Oe6o+aKM/CA3LJ27VqIxWK9JkHSebQ5CX355ZfYuHEjcnNzwePxUFVVRaMhM/DFF1/gySefhLu7u17bZWoGda0aTKMGjy9uta5UKgUA/Prrr+Dz+QgNDdVrLMS0tTkJubi4YP369QAAb29vfPHFF3BwcDBYYMQw/nxzINNqcFfuAKmL/u/Zkf8hR+67uXB9KR8SF/8H1meMYcWKFSgvL0d6ejrNtO9C2nV1LC8vjxKQGeDxBbAf8wosffW/BrTYWQyvGC+I7Np2yZ/H4+Hjjz/G1atXkZCQoPd4iOlq80jo448/bnOjUVFR7QqGGA9jDHUZx2DhMwQCS1u9ty+wFMA60Bo1V6weeE7onsDAQLz66qt499138cILL+jtniVi2tqchP66/EJZWRnq6+t1e41XVlbC0tISzs7OlIQ6AVXpDVR8vxnOz70Li94hem9fXa1G1dkqCOyrILBwarXunw8RNaKRqFV+ibVr1+o2WCTmrc2HY3l5ebrH2rVrERgYiOzsbNy5cwd37txBdnY2goKCsGbNGkPGS/SkLvME+Ja2kHoPMkj7qrsqlHxdAk11xUO9TmBhA6dn36a/oy6kXRNYV61ahf3796Nv3766sr59+2Ljxo147rnnMHPmTL0FSPSPMS3qrvwCy77DweMLAEDvs9ktvCwQkBSAmiu92nw4do+0RwDs7Oxw+/ZtODk5gc+nG/vNWbv+dYuLi6FSqZqVazQa3L59u8NBEcNSFl2FpqYMVn7t20XDGEpLS9GnTx/s3LmT61CIgbUrCT3xxBNYsGABzp8/r7vL9fz583j11VcxZswYvQZI9I8nlKDbI+GQ9OhnsD4UJQrkxedBdaeoXa93dnbG5MmTsWLFCto80cy1Kwl99tlncHd3x5AhQyCVSiGRSDBkyBC4ubkhKSlJ3zESPRO79ILDuCjdoZgh8AQ8CK2FQAf6iI+Ph0qlwsqVK/UYGTE17Ton5OTkhMOHD+PatWvIzs6GWq1G//790adPH33HR/RMWV4AVXkBLPuEGjQJiZ3E8HjNAzVXXB76nBDwv3NUgsHTsG3bdhyo7wuJmy9urp+g50gJ19q9suLOnTuxceNG3YLlvr6+iImJwfz58/UWHNG/2rQfUH/1F1j2HWbQfpiWQavQgmk16MiUMOugp6AqLwBfLNVfcMSktOtwbNWqVYiOjsbEiROxb98+7Nu3DxMnTsTixYvx9ttvP1RbiYmJ8Pb2hlQqRXBwMH7++ef71i0uLsaMGTPQt29f8Pl8xMTEtCf8Lkur1aI+53TjVTGeYa84yQvkyH4tG8rSmx1qh8cXwOFvr0Pk4KGfwIjJaddf4tatW/Hpp5/q1qGZNGkS4uLisGPHDmzbtq3N7ezduxcxMTFYuXIlLl68iJEjR2LcuHEoKChosb5CoYCTkxNWrlyJgQMHtif0Lu3s2bPQ1JTDsu9wg/cldhLDI9IDQjsXvbSnLMtHyRdvoKSkRC/tEdPRriSk0WgQEtL8Ltvg4GCo1eo2t7NhwwbMmzcP8+fPh7+/PzZt2gQPDw9s3bq1xfo9e/bE5s2bMXv2bNja6n+qgbnbv38/+FZ2kLg/eEJpRwmsBLAdYguBtJt+2uvWHaq7RVi8eLFe2iOmo13nhF588UVs3bq12Wp4O3bsaPONikqlEhcuXMDy5cublIeHh+P06dPtCatFCoWiydKh9xZm64oGDRoE26HPG/SE9D3qWjVqLtaA160GAmnHJzsLLGzQ/YkF+Prrj/CTui8seg/WPUcnqzu3Dp2YPnr0KIYObZyB/dtvv6GwsBCzZ89GbGysrt79lu0sLy+HRqOBi0vT4bqLi4teh9xxcXF477339NZeZzZz5kyszLAzSl+qchVu7bwF15dK9ZKEAMCq32jUXf4JFUcTIZuXCL645UX1SOfSriR0+fJlBAUFAQByc3MBNF62d3JyarJgeltWyvtrHcaYXlfYW7FiRZOkWF1dDQ+PrneS89ChQ3Byan0iqT5JvaQI2BmAmpxeemuTx+PB/smFKP33KqgriyF21l/bhDvtSkLHjx/vcMeOjo4QCATNRj2lpaXNRkcdIZFIIJFI9NZeZ8QYQ2xsLMLCwgDHjm/x3BY8Hg8Q/Pe/elw6WmTnCtn8bUY5pCTGwdnMQLFYjODgYKSkpDQpT0lJwbBhhr2Hpau5lHEZ169fx3PPPWe0PhWlCuRvyofqrv6vZvH4AqhrKlBxZAu0Kjl6Lv+u2YN0HpxuAx0bG4tZs2YhJCQEoaGh2LFjBwoKChAREQGg8VDq1q1b+Ne//qV7TVpaGgCgtrYWZWVlSEtLg1gsRr9+hpsH1dntT/4GdnZ2ePzxx4ETKQ9+QSfAlA2ozTgGnlAM+ycWcB0O6QBOk9C0adNQUVGB1atXo7i4GP3798fhw4fh5eUFoPHmxL/eMzRo0P/Wv7lw4QJ2794NLy8v3Lx505ihdxqMMez7v2/w9NNPQyxufcF5fZI4S+AV44WaK656PRy7R+TQA93D5uDuT0mw6D0YFj0D9d8JMQpOkxAAREZGIjIyssXndu3a1ayM9qZ6OGot8PyzUzDmyfFG7ZcxBmj/e6HBQH1Yh0xC/fWzqDi8CbK5CeDr6Z4kYly0WpSZEwl4WPPu3xtPShuRPF+OzHmZUN6+YbA+eDw+HCfEgGlUUNzKNlg/xLAoCZm5xHNKFBQUGr1fkaMI7vPcIbR1Nmg/QhtnuL/yaZObF0nnQknIjF0u1WDhYTnS0i8ZvW9hNyG6j+wOgYW1wfviSyzBNCrcPfk5VOXGT7ikYzg/J0QMZ+9lFeykwNIzWkSdM+5la02dBrWZtWDiWggk9gbvj2k1aLj2Oxqu/w7XWS3fpU9ME42EzBRjDF9nqjHFTwS+0Pi7mSrLlChMLIS60jhrjvNFUjg+vQzqqlJUHPmELmB0IpSEzNTFtHRcv6PF9P7cbKcs9ZTCf6s/xM49jdan2MkLDuOiUJ91kvYs60TocMxMWVlaYtFgER73FgBK4/fP4/MgsBA0Tq8w4qDEyn8UFEVXkZ6ervd5iMQwaCRkpvr27YNPxltAyOfmS6gsU6JwayFURjoc+7Puj81FUlISeDwetFqt0fsnD4dGQmbo4sWLuPD7aczRMIgE3CQhpmFQ16gBrcboffP4AvB4PBw5cgTLly9HSkoKHB0dW5xTRmsRcY9GQmYoISEB6z74BwQc/utKXCXwftMbInsZZzH4+vqiqKgITz31FOrq6jiLg7SORkJm4t7/5bVKOf74cg/sBz8FPu8Ix1Fxq1evXjh8+DDCwsIwdepUsH6vgCegP3lTQyMhM1N/7QyYsgG2/UdzGkdDfgMy52dCWWK4aRttERwcjOTkZKSkpODucdpS2hRREjIzdRk/QuLRH2I97XLRXqLuIrhOd4XARj9Lu3ZEeHg4vvnmG1iHPM11KKQFlITMTLdHxsJ26PNchwGhjRAOYxwgsDSNXVEmTJgAkZ0rtIo6VP78FZim7bvCEMOiA2QzY9Xv3mx5OadxaOo1qM+phxZ14EuMt47RgyiKr6Hqt31QluXBadKbXIdDQCMhs8E0KlR8/zFUFaYxgVNZqmxc3rXStDYrtOgZCOdnVkKel4rb//477ty5w3VIXR6NhMxE/dXTqL10FNaDJ3MdCgBA2kOKvpv6ouGWFyf9t7bOtEXvwXCe9j7Kkt/HyJEjceHCBUiltNc9VygJmYma1O8g8XwEYkdPrkMBAPCEPIjsRJCXCI06baOtpD36wfXFDxH7qIASEMcoCZmBixcvQnErC46TV3Adio6yXImyg2WwDCiF0KYH1+G0SGTvjjXXgDXLv0PVb/sA8ODy6HhcoT0VjYqSkBn4xz/+AYGtCyx9h3Idig5TMchvyWHRV8V1KG2iVcpRfWYvVPkXcftZLVy60elSY6FP2gxs2LABThOXmtSGgBI3CXqv6g2RgzvXobRJ91Gz4Dx1NeSlN9F/ax3+ndk5kqc5oCTUySmVSri4uEDi7sd1KJ2ehXcQvOduwOieAnx9WUULoxkJJaFO7OrVq/Dw8EBqairXoTTTUNCA7IXZUN7O4zqUhyK0ssW+5y2x+1kL8Hg8JCcnIz4+HnI5t/ddmTNKQp3Y0qVLIZFIEBAQwHUozYhsRXCc4AhBt+5ch9IuUiEP/n//AXM3fYNly1egm4sXvvrqK2g0xl+axNxREuqkDh48iEOHDmHjxo2QSCRch9OM0FYIp/FOEFjZcR1Kh9gNmw7ZvC0Qu/TCiy++iAEDBiAvr3ON7kwdXR3rhDyX/B+KkiIh9Q5G7FkJlhh5J4220DRo0HCzAVplA/hG3H7aEEQOHnB+5m0oiq6iMP0IRm/NAI+fhfprvyNnZyzs7Oy4DrFTo5FQJ6Str4bQ1hn2Y1812TWUlbeVuPnBTajuFnEdit5IZH3hMC4KPL4AmoZqlP1nPWQyGebNm4dff/2VTmS3EyWhToYxBqGtM1xnrIeoO3erFj6IRCaB7we+EDt6cB2KQQgsbOAesRMrV67Ejz/+iBEjRmDgwIG0pnU7UBLqRC5fvoygoCCTmxTaEr6YD4mLBDxh5z4Ua42wmz1WrlyJGzdu4Pjx44iMjASfz0d9fT38/f0RGRmJQ4cOobq6mutQTRoloU7i0qVLCA8Ph1arhcDSjutwHkhZoUTRl0VQV5dzHYrB8fl8jB49GhEREQCAuro6jB07FkeOHMGkSZNgb2+PESNG6K6s3b17l8twTQ6dmDZh92aCywsyUJr8PkR2rhCMWgqB2PQnXGrlWtRdqYOkZwPXoRjd4I/OApZPgj0XDlllCeT5acioKoVAIABjDD4+PpBKpRg0aBAGDhyI/v37Y8KECbCxseE6dE5wnoQSExPx4Ycfori4GAEBAdi0aRNGjhx53/onT55EbGwsMjMzIZPJ8Oabb+r+D2SONPJalO5/DxJZHzhNeRt8iSXXIbWJ1F0K3/d9UXPFwyRn0etLa0uG8Hg8iLq7QdTdTVem1WqRlJSEc+fOIS0tDbt27UJRURHy8/NhY2ODJUuW4OLFi+jVqxd69uwJT09PDB8+HE98egWMacHjNT14MYctizhNQnv37kVMTAwSExMxfPhwbN++HePGjUNWVhY8PZsvSZGXl4fx48djwYIF+PLLL/Hrr78iMjISTk5OePbZZzl4B/qnVqtx4cIFfPXVV9AKRkAg7QaX6WshdvUxqblhpH0EAgGmTJmCKVOm6Mru3r2ru8zv6+uLP/74A+np6Th48CDKysqwZcsWAF6oyzqJO0cSILC0A9/SFgJLW7xncR7vvPMOGGP45JNPYGtrC1tbW1hbW8Pa2hqBgYEQi8Wora2FQNC4bImpXVHlMQ6vKz766KMICgrC1q1bdWX+/v6YPHky4uLimtVftmwZDh48iOzsbF1ZREQE0tPTcebMmRb7UCgUUCgUut+rqqrg6emJDRs2wMLif2s2+Pn5ISgoCKWlpTh27FiTNqRSKZ555hkAwIEDB1BfX9/k+cceewxubm5IT09HZmZmk+e8vb0RGhqKqqoqHDp0CFqtFlqtFmq1GlqtFq+88goA4K233sK5c+dw6dIlyOVyuLi4gI2OgsTNt9XP8H6kkOO8dCEAIES+BXIY8RCOp4TI8u/I35IPx6feh9ixt/H67iBOP7cWaFWN00X4IilUFX9Anp8GTV0VtPIaaBqqILRzRfdRc6BVNuDW9vmApunEW7f52yDsZo/yQ/9AQ+7ZxkKBGDyhGHbDpqLbwL+hoeASKk99AZ5AAB5fhEd9nOHj44NNmzYBABYsWNB4LlIggEAggFAoxKJFizBkyBBUVlbC1raD64gzjigUCiYQCFhycnKT8qioKDZq1KgWXzNy5EgWFRXVpCw5OZkJhUKmVCpbfM0777zD0HhAQA960EPPj9zc3A7nAs4Ox8rLy6HRaODi0nRrGhcXF5SUtHwJuqSkpMX6arUa5eXlcHNza/aaFStWIDY2Vvd7ZWUlvLy8UFBQ0PEMbkDV1dXw8PBAYWGhSZ+wpDj1q7PEee+Iwt7evsNtcX5i+q/Hp4yxVo9ZW6rfUvk9EomkxblVtra2Jv2PfI+NjQ3FqUcUp37x+R2/y4ez+4QcHR0hEAiajXpKS0ubjXbucXV1bbG+UCiEgwP3m+wRQh4eZ0lILBYjODgYKSkpTcpTUlIwbNiwFl8TGhrarP7Ro0cREhICkUhksFgJIQbU4bNKHfD1118zkUjEdu7cybKyslhMTAyzsrJiN2/eZIwxtnz5cjZr1ixd/Rs3bjBLS0u2ePFilpWVxXbu3MlEIhHbv39/m/uUy+XsnXfeYXK5XO/vR58oTv2iOPVLn3FymoQYY2zLli3My8uLicViFhQUxE6ePKl7bs6cOSwsLKxJ/RMnTrBBgwYxsVjMevbsybZu3WrkiAkh+sTpfUKEEEITWAkhnKIkRAjhFCUhQginKAkRQjjV5ZJQYmIivL29IZVKERwcjJ9//pnrkFoVFxcHHo+HmJgYrkNpQq1W4+2334a3tzcsLCzQq1cvrF69mvPlTU+dOoWJEydCJpOBx+PhwIEDuudUKhWWLVuGAQMGwMrKCjKZDLNnz0ZRkfHXwW4tznuys7MxadIk3az4oUOHoqCgwGgxxsXFYfDgwbC2toazszMmT56Mq1evNqnDGMO7774LmUwGCwsLjB49utkk7gfpUkno3tIhK1euxMWLFzFy5EiMGzfOqP+wD+PcuXPYsWMHHnnkEa5DaeaDDz7Atm3bkJCQgOzsbMTHx+PDDz/EJ598wmlcdXV1GDhwIBISEpo9V19fj9TUVKxatQqpqalITk5GTk4OJk2aZFJxAkBubi5GjBgBPz8/nDhxAunp6Vi1ahWkUuPN6j958iQWLlyI3377DSkpKVCr1QgPD0ddXZ2uTnx8PDZs2ICEhAScO3cOrq6uGDt2LGpqatreEce3CBjVkCFDWERERJMyPz8/tnz5co4iur+amhrm6+vLUlJSWFhYGIuOjuY6pCYmTJjA5s6d26TsmWeeYS+++CJHETUHgH3zzTet1jl79iwDwPLz840TVAtainPatGkm9VkyxlhpaSkDoLuXT6vVMldXV7Z+/XpdHblczmxtbdm2bdva3G6XGQkplUpcuHAB4eHhTcrDw8Nx+vRpjqK6v4ULF2LChAkYM2YM16G0aMSIEfjxxx+Rk5MDAEhPT8cvv/yC8ePHcxzZw6mqqgKPxzOpvcO0Wi2+++479OnTB08++SScnZ3x6KOPtnjIZkxVVVUAoJs5n5eXh5KSkibfKYlEgrCwsIf6TnE+i95Y2rN0CFe+/vprpKam4ty5c1yHcl/Lli1DVVUV/Pz8IBAIoNFosHbtWrzwwgtch9Zmcrkcy5cvx4wZM0xqxnppaSlqa2uxfv16vP/++/jggw/www8/4JlnnsHx48cRFhZm9JgYY4iNjcWIESPQv39/ANB9b1r6TuXn57e57S6ThO552KVDjK2wsBDR0dE4evSoUY//H9bevXvx5ZdfYvfu3QgICEBaWhpiYmIgk8kwZ84crsN7IJVKhenTp0Or1SIxMZHrcJq4d3L/6aefxuLFiwEAgYGBOH36NLZt28ZJElq0aBEuXbqEX375pdlzHf1OdZkk1J6lQ7hw4cIFlJaWIjg4WFem0Whw6tQpJCQkQKFQQCDgfq3ppUuXYvny5Zg+fToAYMCAAcjPz0dcXJzJJyGVSoWpU6ciLy8PP/30k0mNgoDGv1WhUIh+/fo1Kff3928xCRja66+/joMHD+LUqVPo0aOHrtzV1RVA44jozwsKPux3qsucE2rP0iFceOKJJ5CRkYG0tDTdIyQkBDNnzkRaWppJJCCg8UrTXxe0EggEnF+if5B7CejatWs4duyYSa5DJRaLMXjw4GaXw3NycuDl5WW0OBhjWLRoEZKTk/HTTz/B29u7yfPe3t5wdXVt8p1SKpU4efLkQ32nusxICABiY2Mxa9YshISEIDQ0FDt27EBBQYFJbRlkbW2tO+a+x8rKCg4ODs3KuTRx4kSsXbsWnp6eCAgIwMWLF7FhwwbMnTuX07hqa2tx/fp13e95eXlIS0uDvb09ZDIZnnvuOaSmpuLbb7+FRqPRjYzt7e0hFhtvt9jW4vT09MTSpUsxbdo0jBo1Co899hh++OEHHDp0CCdOnDBajAsXLsTu3bvxn//8B9bW1rrPytbWFhYWFrr719atWwdfX1/4+vpi3bp1sLS0xIwZM9rekV6v4XUCrS0dYqpM8RJ9dXU1i46OZp6enkwqlbJevXqxlStXMoVCwWlcx48fb3FB9jlz5rC8vLz7Lth+/Phxk4nznp07dzIfHx8mlUrZwIED2YEDB4wa4/0+q3/+85+6Olqtlr3zzjvM1dWVSSQSNmrUKJaRkfFQ/dBSHoQQTnWZc0KEENNESYgQwilKQoQQTlESIoRwipIQIYRTlIQIIZyiJEQI4RQlIUIIpygJEUI4RUmItAtjDK+88grs7e3B4/GQlpZmlH7Xr1+P0NBQo/RFjIOSEGmXH374Abt27cK3336L4uJio02uTU9Px8CBAw3StqluKmDuKAmRdsnNzYWbmxuGDRsGV1dXCIXtW5BBqVQ+VP309HQEBga2q6/WmPKmAuaOkpAZOnv2LEaPHg0LCwv4+fnpvmD62lXipZdewuuvv46CggLweDz07NkTAKBQKBAVFQVnZ2dIpVKMGDGi2RK1o0ePxqJFixAbGwtHR0eMHTv2vv1kZ2fr3segQYNw/vx55OTk6H0kVFtbi5kzZ+LTTz9F9+7dH1h/9OjReP311xETE4Pu3bvDxcUFO3bsQF1dHV5++WVYW1ujd+/e+P777/Uap9nS59R/wr0zZ84wqVTK4uLiWE5ODpsyZQobP3488/HxYampqXrpo7Kykq1evZr16NGDFRcXs9LSUsYYY1FRUUwmk7HDhw+zzMxMNmfOHNa9e3dWUVGhe21YWBjr1q0bW7p0Kbty5QrLzs5usY/s7GxmbW3NlixZwq5fv86Sk5OZTCZjfD6f1dbWNqm7du1aZmVl1erj1KlT930/s2fPZjExMbr4HrRsSlhYGLO2tmZr1qxhOTk5bM2aNYzP57Nx48axHTt2sJycHPbaa68xBwcHVldX15aPtEujJGRmQkND2cyZM3W/7927l/H5fDZlypQ2vX7OnDmsX79+bN26dU1+/quNGzcyLy8v3e+1tbVMJBKxr776SlemVCqZTCZj8fHxurKwsDAWGBj4wDgef/zxZlveTJ8+nfXp06dZ3YqKCnbt2rVWH/X19S32s2fPHta/f3/W0NCgi68tSWjEiBG639VqNbOysmKzZs3SlRUXFzMA7MyZMw98r11dl1pZ0dz98ccfOHPmDD788ENdmVgsBmMM77333gNfn56ejtu3byMzMxPp6ek4depUm3fTzM3NhUqlwvDhw3VlIpEIQ4YMQXZ2dpO6ISEhrbaVn5+Pn376CampqU3KRSJRi4di9vb2um1oHkZHNhX487kjgUAABwcHDBgwQFd2b43l0tLSh46rq6FzQmbk3pf9z1/yq1evYsiQIU2+IFevXsX48eMRHByM0aNHo7y8HFlZWRg/fjwyMjLg7Oys+3nUqFFt6pv9d228tuy8YGVl1WpbaWlpEAqFTWIGgNTU1BZPSq9btw7dunVr9dHSdt9/3lRAKBRCKBTi5MmT+PjjjyEUCqHRaO4bo0gkavI7j8drUnbvPZv6mtumgEZCZqSqqqrJQvh37txBfHx8k8vnCoUCCxcuxK5du9CjRw8kJCQgKSkJy5cvx9SpUzFq1ChMmTIFixcv1v3cFj4+PhCLxfjll1906wurVCqcP3/+oS958/l8aLVaKJVK3VW3w4cPIzMzs8UkFBERgalTp7bapru7e7Oye5sK/NnLL78MPz8/LFu2zGQ2FTB3lITMSGBgIDQaDeLj4/H8888jOjoaXl5eyM7ORn5+Pry8vHDgwAFkZWXhqaeeAtCYlObPnw8AyMjIQGRkZLOf28LKygqvvfYali5dqlusPT4+HvX19Zg3b95DvY/g4GCIRCK88cYbeOONN3D58mW89tprAKDXw7HOsqmAuaPDMTPi4+OD1atXY/PmzRg0aBDc3Nxw9OhReHh46LaTzsjIwEcffaTbTig7OxtLliwB0Hhep3fv3s1+bqv169fj2WefxaxZsxAUFITr16/jyJEjbbrs/WcymQxJSUk4dOgQQkJCsHnzZsyZMweOjo4tjmhI50YL3XcxCQkJOH/+PHbt2gWgMSkNGDAAZWVlGDduHM6fP9/kZ0IMjUZCXczLL7+MyspK+Pn5YeDAgdi9ezeAxmQUEBDQ7GdCDI1GQoQQTtFIiBDCKUpChBBOURIihHCKkhAhhFOUhAghnKIkRAjhFCUhQginKAkRQjhFSYgQwilKQoQQTlESIoRw6v8BtcmqKgbCqs8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Only run if anisotropy has been computed\n",
    "create_fig3a()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}