{
 "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 timflow.steady as tfs"
   ]
  },
  {
   "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 = tfs.Model3D(kaq=kx, z=np.linspace(H, 0, naq + 1), kzoverkh=kz / kx)\n",
    "    tfs.LineSink1D(ml, xls=0, sigls=2, layers=np.arange(int(d * 10)))\n",
    "    tfs.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 = tfs.Model3D(kaq=kx, z=np.linspace(H, 0, naq + 1), kzoverkh=kz / kx)\n",
    "    tfs.Well(ml, xw=0, yw=0, Qw=Qw, rw=rw, layers=np.arange(int(d * 10)))\n",
    "    tfs.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": "markdown",
   "metadata": {},
   "source": [
    "Find effective vertical hydraulic conductivity for one realization of 400 layers. This takes significant time, so it is commented out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "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",
    "kz = brentq(func, a=0.001 * kx, b=kx, args=(kx, h0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the experiment for Figure 3a. 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": null,
   "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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import lognorm\n",
    "\n",
    "\n",
    "def create_fig3a(kaq, aniso):\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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAACqCAYAAAAJKkK3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJtVJREFUeJzt3XlclNX+B/DPMys7siPKIqGipiAIrqklSZaaueY1MTXrGi5IkstVqdRANPG6pBfU9FYmUhdL82JqLmUYBpIbgqKCyg4Css4wc35/+HOuE6CAM/PMwPf9es3rNfNw5jkflvnyrOdwjDEGQgjhkYDvAIQQQoWIEMI7KkSEEN5RISKE8I4KESGEd1SICCG8o0JECOGdiO8AuqZUKpGbmwtzc3NwHMd3HEIMEmMMDx48gJOTEwSCZ9+eaXeFKDc3F87OznzHIKRNuHPnDjp37vzM62l3hcjc3BzAwx+ghYUFz2m0TFYFfNb94fMPMgCJqVa7q5ZX46X4lwAA+4ftx+GEw5g8eTLs7e212q/O6Pjnqc8qKirg7Oys+jw9q3ZXiB7tjllYWLSDQiQEpP+/+2lhofUPjkgugtBYCADw8PDA8uXLtdqfzun452kINHV4gw5WE62oqKjA8ePHUVFRwXcUYgCoEBGtyMrKwssvv4wbN27wHYUYgHa3a0Z0o2fPnrh16xY6duzIdxRiAKgQEa2QSqVwc3PjOwYxELRrRrTiTs4dBAcHIycnh+8oxABQISJaUVVVhaSkJFRVVfEdhRgA2jUjWuHZwxOpqal8xyAGgraICCG8o0JEtOLSxUvo2LEjLl68yHcUYgCoEBGtsLWzRXBwMOzs7PiOQgwAHSMiWtGxY0esWLGC7xjEQNAWEdGKyspKJCUlobKyku8oxABQISJacf36dQwaNAiZmZl8RyEGgHbNiFZ4enri8uXLcHd35zsKMQB6sUW0bds2uLm5wcjICP3790dycnKTbWNjY/HCCy/AysoKVlZWCAgIeGJ7wg9jY2P06tULxsbGfEchBoD3QhQXF4fQ0FCEh4cjNTUVXl5eCAwMRGFhYaPtT506halTp+LkyZNISkqCs7MzRo4ciXv37uk4OXmSe3fvISwsDHfv3uU7CjEAvBeijRs3Ys6cOZg5cyZ69uyJHTt2wMTEBLt37260/ddff433338f3t7e8PT0xM6dO6FUKnHixAkdJydPUl5ejh9++AHl5eV8RyEGgNdjRDKZDCkpKVi2bJlqmUAgQEBAAJKSkpq1jurqasjlclhbWzf69bq6OtTV1ale00BdutGzV09kZGTwHYMYCF63iIqLi6FQKODg4KC23MHBAfn5+c1ax5IlS+Dk5ISAgIBGvx4REQFLS0vVgwbOJ0T/8L5r9iwiIyOxf/9+JCQkwMjIqNE2y5YtQ3l5uepx584dHadsn65euQoPDw9cuXKF7yjEAPC6a2ZrawuhUIiCggK15QUFBXB0dHziezds2IDIyEgcP34cffr0abKdVCqFVCrVSF7SfB06dMCkSZPQoUMHvqMQA8DrFpFEIoGvr6/ageZHB54HDhzY5PuioqKwevVqJCYmol+/frqISlrIqZMTIiIi0KlTJ76jEAPA+65ZaGgoYmNjsXfvXqSnp2Pu3LmoqqrCzJkzAQBBQUFqB7PXrVuHlStXYvfu3XBzc0N+fj7y8/PpVgI9U1NTg4sXL6KmpobvKMQA8H5l9ZQpU1BUVIRVq1YhPz8f3t7eSExMVB3AzsnJUZvSdvv27ZDJZJg4caLaesLDw/HRRx/pMjp5gmvXrmFI/yFISUmBj48P33GInuO9EAHAvHnzMG/evEa/durUKbXXt2/f1n4g0iSZTIbk5GRUVlbCysoK3bt3b/Q4ULdu3ZCcnIzu3bvrPiQxOHpRiAj/3Jb+2GDZ7cjX1F7HxsYiNDRUbTf4hRdewJkzZxq819TUFH5+fpoPStokKkSkSW5LfwRjDPVl+RBbdYSi2hTLly/HyJEj4eDggJKSElXblJQUbP/Xdih8FBAaC5GXl4cvd3+J9957j+Y2I09FhYg8UfnZb/Ag9TA6/X0XhCaWaicOOnfurHp+8+ZNxH0Th7qEOriFuqGkuAQ7d+7EhAkTqBCRp6JCRJpUfi4e5Wf3ocPQIAgkD++ib2oXbtKkSXje+3n0fbEvbq69iaohVXTDK2k23k/fE/1UlXEWZaf3wnLQm7AcOLlZ73F1c4X7cndIO0kRNC1I7R4/Qp6EtohIA0xRj/sndsKk2yBYDpn21PaqrSROBnNPIdwWu6Ho2Hj4+flh//796Nmzp5YTE0NHW0SkAU4ogsOba2D9yjxwHNfi9wskAkjsXDBo0CCEh4fTWFHkqagQETV1+TeglNVCbN0JQmOLVq9HZGGLjz76COfOncPEiRNpN408ERUioqKU1aDou09w/+SuZ14Xq5ehsrIS+/btQ2pqKhYsWKCBhKStokJEVCqSE6CoqYDlgIlPb/wUsuI76Nq1K0xNTbF9+3bExMRg7969GkhJ2iIqRAQAUF9Ziork72DhOxYiS4env+EpxFZOcJj6KcZ9mYVPMh1g1mckQvacevagpE2is2YEAFBxLh6cUNLsU/VPI5Aaw8jlf+NEWb8yv1UHvkn7QFtEBABg7N4PViPegcDITCPrU1SVofz376Coug8AqiIUGRmJlStXaqQP0nZQISIAAGN3X5g9P0Jj61NU3kfFuXhVIXpEqVTi008/xblz5zTWFzF8VIjauaqqKkyaNAnyYs2O5S1x6ALnhfshsVef6fXDDz9Ev379MGPGDFRXV2u0T2K4qBC1c7t370ZCQgI4sUQn/XmsOIq7fWbi+s3bcBw+rdF710j70+xCZG1tjeLiYgDArFmz8ODBA62FIrrBGMPnn3+OCRMmaORM2ePkJfeQ9+UHkJc0vPFVbOMMyyHToCgvBGNKjfZLDFOzz5rJZDJUVFTA1tYWe/fuxbp162Bubq7NbEQLHt8Cqc25iIJr13C/bxAan4yp9TiRGBJbV3Cixre0LPzfAMDRmTQCoAWFaODAgRg3bhx8fX3BGMOCBQtgbGzcaNumposm+uXBn0chsu4MqXNvja9bZGkPm1FNX03NcQ83xqszkxAdnYlFixZpPAMxHM0uRF999RWio6ORlZUFjuNQXl6O2tpabWYjWmYd8HcoKgq1slXCFPVQVFdCaGwBTtj0n5ms4CZCF0di/WUJJHZuquV/HaaWtG3NLkQODg6IjIwEAHTp0gVffvklbGxstBaMaBdjDEJjcwiNtbN7LSvKRv6eMDjO2ASpo0eT7SwHTkbVtV9QkrgFjtOiwAmEWslD9FurzprdunWLipABY4wh/8vFqLx84umNW0ncwRF2E8MhtnryMLGcSAybUfMhy83AgwtHtJaH6LdmbxFt3ry52SulO631W929a5DlZUBoHqS1PgRGpjB5rnmzeBh17gUz71GounIS5j6vqY4fkfaj2YUoOjpa7XVRURGqq6tVc1qVlZXBxMQE9vb2VIj0XNWVExBa2MHIRfMHqR9RVJejOv13mHi+AKGJ5VPbW704C5xAREWonWr2b/3WrVuqx9q1a+Ht7Y309HSUlpaitLQU6enp8PHxwerVq7WZlzwjVi9DdfovMO31klY/9IqKEpT+vBOKB8XNai+QGIMTiSErykZN9p9ay0X0U6v+EleuXIktW7aozeLZvXt3REdHY8WKFRoLRzSvNucSlHVVMHv+Ja32I3F0h+vig5A4PNei95Wf3YeSQxtQVlamnWBEL7WqEOXl5aG+vr7BcoVCgYKCgmcORbTH2N0XTu/thNi6E99RGmX10hwo5bVYsmQJ31GIDrWqEI0YMQLvvfceUlNTVctSUlIwd+5cBAQEaCwc0az6+nowxiDu4Kj1vuSluSjY/w/IS1s2cL7IwhZWw2ciJiYGp0+f1lI6om9aVYh2794NR0dH9OvXD1KpFFKpFH5+fnBwcMDOnTs1nZFoSGxsLPJ2BYMp5NrvTCCEwMQSaMV1QWber2Dw4MH48MMPwRjTQjiib1o1QqOdnR2OHDmC69evIz09HQDg6emJbt26aTQc0axvvvkGQgtbcEKx1vsSd3CA3dgPW/VejhPg3//+N6RSKd2L1k60eqjYXbt2ITo6GtevXwcAdO3aFSEhIXjnnXc0Fo5ozp07d/HLL7/A5lXd3NPFlAowWTU4sbRVV0u7uz8cx6isrAzFxcXw8Gj66mxi+Fq1a7Zq1SosXLgQY8aMQXx8POLj4zFmzBgsWrQIq1atatG6tm3bBjc3NxgZGaF///5ITk5usu2VK1cwYcIEuLm5geM4bNq0qTXx26W4+O8glUph0m2gTvqTFd7GnU2TISu89Uzr+dvf/oaJEydCLtfB7iThTasK0fbt2xEbG4uIiAiMHTsWY8eORUREBGJiYvD55583ez1xcXEIDQ1FeHg4UlNT4eXlhcDAQBQWFjbavrq6Gu7u7oiMjISjo/YPuLYlP586jdGjR0MgNdFJf6IODrB9fSlEz3hgfM2aNbh8+TIiIiI0lIzoo1YVIrlcjn79+jVY7uvr2+hp/aZs3LgRc+bMwcyZM9GzZ0/s2LEDJiYmTQ4j4ufnh/Xr1+PNN9+EVCptTfR26/DBbxETE6Oz/oRGZjD1HALhMw7G7+Pjg2XLlmH16tX480+60LGtatUxounTp2P79u3YuHGj2vKYmBhMmzatWeuQyWRISUnBsmXLVMsEAgECAgKQlJTUmliNqqurU5vuuKKiQmPrNhS19QxGAgGsra111qei5gFqrl+AsYd/q6aufnwAN6bwBWfVGf6BE1CTmwmBgG4DaWue6WD1Tz/9hAEDBgAAfv/9d+Tk5CAoKAihoaGqdn8tVo8UFxdDoVDAwUF9iFIHBwdcu3attbEaiIiIwMcff6yx9RmSHqsSkW4EDNhZhXunZ8N08HSd9V1fXoiSI5vgOGNTqwrR4zihGLajQ6GU1VARaqNaVYguX74MHx8fAEBWVhYAwNbWFra2trh8+bKqnT6cel22bJlaYayoqICzszOPiXQrvUiBPwuU6DTYVaf9Shzc4RL2PaCh+9kezQaiUCiQm5vbrn6H7UGrCtHJkyefuWNbW1sIhcIGt4QUFBRo9ED0owsu26u4K3JYSAHTLn1R9/TmGsNxXKsuZnyakJAQJCYmIi0tDaamphpfP+EHb9u5EokEvr6+OHHif4NzKZVKnDhxAgMH6uYUc1vHGEPclXqM8xRDINL+RYyPk9/PR+F3n0B+P0+j612wYAFyc3OxcOFCja6X8IvXHe7Q0FDExsZi7969SE9Px9y5c1FVVYWZM2cCAIKCgtQOZstkMqSlpSEtLQ0ymQz37t1DWloabty4wde3oNfqH5Qip1yJKb1afShQ73Tt2hWbN2/Grl27EB8fz3ccoiG8/oVOmTIFRUVFWLVqFfLz8+Ht7Y3ExETVAeycnBy1g5O5ubno27ev6vWGDRuwYcMGDBs2DKdOndJ1fL0ntrBBUZg5JEIAMh33beUI+wktu7i1uWbNmoWjR49izpw5eOGFF+iasjaA93+V8+bNw7x58xr92l+Li5ubG90E2UyMMShltTCx4OeEAWMMUCoATqDxkxYcxyEmJgYJCQkNzroSw0TnQtuoP/74A9e3zMKlAgUv/csKbiJn/euQFWRpZf0dOnTAzJkzwXGc2plaYpioELVRcXFxEEiM0cOOn1+xyNIeNq+GQGRpr9V+fvnlF/Tu3RvfffedVvsh2kWFqA1SKpWIi4uDefcBEAn42TUTGpvDrHfAM1/M+DRDhgzB5MmT8fbbb+Pq1ata7YtoDxWiNujs2bO4e/cuLHoM5i2DorYSVdd+haK2Uqv9cByHpE6TUGdkA68hAXAJiVO7PYQYBipEbdD58+fh5uYG486evGWoLytA8feRqC/L13pfAokx7MavgLLmAcrOfqP1/ojmUSFqg0JDQ3HlyhVe5wiT2LvBOeQAJPZddNKf2Koj7KesQYcXdHc/HdEc3k/fE80qLS2FpaUlTEx0M+5QUziBEJxUovH1Pmm3S+r4cBRHefEdfPHFF6oLY4n+oy2iNmbu3Ll49dVX+Y4BeVkBin6IglwHu2Z/VZV+GrNnz8aBAwd03jdpHSpEbUhFRQV++OEHjBgxgu8ogFIBZXX5w4sadcxyyN8wbdo0vPXWWzh+/LjO+yctR4WoDYmPj0ddXR2mTp3KdxSIrZ3g8OZaXiZy5DgBdu/ejYCAAIwbNw5nz57VeQbSMlSI2pDdu3fj5ZdfprF6AIjFYnz77bcYPnw4DbxvAKgQtREPHjxAaWkpZs+ezXcUAIAs/yayN4zT2i0ezWFiYoLDhw9j+PDhqK+vpwse9RgVojbC3NwcV69excSJE/mOAgAQWtjA+qV3IDS35TsKAGD9+vXw9/enURr0FJ2+bwPkcjlc5+6ExNZFbbkxT3kAQGhiCXOf0bz1/9fT/EpZNwwaNAijRo3C/v378frrr/OUjDSGtojagMOHDyNv1/uQl97jO4qKsrYK1Vnnoayr4jsKAEAgMcKhQ4cwevRojB8/Htu2beM7EnkMFaI2YMuWLZA4deflDFVT5GX5KPr2Y40PFfssuocfR7J7EEx9xmDpzsNwXXKY7kvTE7RrZuAuXryIkydPwnZMGN9R1EjsXNEp+N9av/u+pThOAOsRc8CYEhzHoeZWKkpKBsDGxobvaO0abREZuM2bN6NTp04w6c7fnfaN4YQiiMyswQn1838dxwnA6mUo+e9m+Pj4aHRST9JyVIgMnJ2dHcLCwvTuA19fXoiS/25GfXkh31GaxIkkcHxrPTp37oyhQ4ciMjISCsXDK8Hdlv7Y4EG0hwqRgYuIiNDLqXVYvRyy4myweh2P2t9CIgs7nDp1CosXL8by5csRFqZfu7jtBRUiA1VWVob169ejqko/zkr9ldimEzpO/wxim858R3kqsViMiIgInDlzRlXU5fdzwRT1PCdrP/Rre540W3R0NKKiovDWW2/RjKfPSH23qxxM+ScKD4SDE0lgNWIOjN28AQA9ViUi3Qiq5zUwwu3I13Sety2iLSIDVFJSgujoaAQHB6Njx458x2mUrOAW7vzzTcgKb/IdpcU4gRB245ZCIDVFYdwKFP5nDWRF2XzHatOoEBkg91ffRWWtHHEyH709kCo0s4LFgEkQmlrxHaVVJA7PwWHaOtiOCYOs8BaKDkaAMSXfsdos2jUzMNnZ2aj443tYDnoTQhNLvuM0SWjaAZb9J/Ad45lwHAfTnsNg0n0Q6ssLwXECXClU4KPTdajxyQQ69eE7YptBhcjAuLq6wn7CKhi56PeHQFlXA1n+NUgcPCCQ8jts7bPihOL/v2q9FoVVDGn5SmR/tRySjt1ge+k4TLoPgUBipPYeOnbUMrRrZkBu3rwJxhiM3X3BicR8x3ki+f1cFHyzHPL7uXxH0agXu4iQMc8UnSc8PIZUcmQTys/uAwCwehntvrUSbREZiOTkZAwdOhSbNm0CoP8Dn0lsneH0bgxEejIMiCYJOA5mHv0g9BgCeVk+OOHDfwoVKYfw4PxBGHcbCIdbqTBy7q32D4O2kppGW0QG4O7du3j99dfh4+NjMDNTcCIJxFZO4ESan8lDn4g7OEJk/vA+NWO3vjDpMRQ1WX+g8MAq3Nk8FQ9SH55IYIp6KJW0tdQU2iLSc3l5eRg9ejTEYjESEhIglUr5jtQs9RXFqPj9ECz8x0NkYcd3HJ2QOLjD2sEdVi+9A3lxNmqyzkNs5woAqLpyEra2MzF48GD0798f/v7+8PHxga1t29tibA0qRHpu/fr1KCoqwtGjR+Hg4MB3nGZTympQm3MJZt6j+I6icxzHQWLnBomdm2qZpGM3vLdwIc6ePYvPPvsMZWVlmDx5MuLi4lBcXIyoqCj06NEDnp6e6Nq1K2xsbMBxHH/fhI7pRSHatm0b1q9fj/z8fHh5eWHLli3w9/dvsn18fDxWrlyJ27dvo2vXrli3bp1ezOWlKbW1tUhJScHgwYMRERGBxYsXw8nJie9YLSKxdYbTrK18x9AbEjtXhIe/DwBQKpW4ceOGalft3r17OHDgALKz/3fRpIuLi+p1REQEJBIJIs8UQmhmDaGZNUQWduCE4jZz3In3QhQXF4fQ0FDs2LED/fv3x6ZNmxAYGIiMjAzY29s3aP/bb79h6tSpiIiIwOjRo7Fv3z6MGzcOqampeP7553n4DjSDMYaMjAzs378fe/bsQVlZGe7evQszMzODK0KkcY1deHo70hNeXl64ffs2qqurkZmZiRs3bqCiokLVJiEhAVevXlW7r9DxrfWQduqBTz/9FN999x1sbGxgbW0NKysrBAYGYty4cSguLsbRo0dhbm4OMzMzmJmZwcLCAp6engCA6upqSKVSCIVC7X/zT8ExxhifAfr37w8/Pz9s3frwv6dSqYSzszPmz5+PpUuXNmg/ZcoUVFVV4fDhw6plAwYMgLe3N3bs2NGgfV1dHerq6lSvy8vL4eLigo0bN8LY+H+jOvfr1w/dunXDzZs3ce7cObV12NraYuTIkQCAffv2Nehj9OjRsLCwwG+//Ybbt28DeFhYAODjXx9AYt8F9Q9KUJudBjAllo70QG1tLUxMTPDuu++ivr4eds7PQVldBoiNYOLhDwu/cRBbP9sNo0aoxR9GwQ+/v9ptqIXRU97xjDgZzLqtBQCUng1CUUIU7N5YDomtq3b71RGd/zwfwxgDq6uCoqoM9dX3IbV/DgKpCdb61OLnn3/G/fv3UVpaioqKCtw2fx4W/V5H7b2rKIr/SG09nIklOr8bCwC4tzsYyooigBMAIgk4oRi2Yz7AjZ2LsGPHDnz55ZeQSCQQiUQQi8UIDAzEwoULkZubi0WLFiExMRFlZWWwtNTAhbWMR3V1dUwoFLKEhAS15UFBQWzs2LGNvsfZ2ZlFR0erLVu1ahXr06dPo+3Dw8MZAHrQgx5aeGRlZWmiFDBed82Ki4uhUCgaHIR1cHDAtWvXGn1Pfn5+o+3z8xufY33ZsmUIDQ1VvS4rK4OrqytycnI0U8m1pKKiAs7Ozrhz5w4sLPRruNXHUU7NMpScj/YsrK2tNbI+3o8RaZtUKm30lLelpaVe/6IfsbCwoJwaRDk1SyDQzKWIvF7QaGtrC6FQiIKCArXlBQUFcHR0bPQ9jo6OLWpPCNF/vBYiiUQCX19fnDhxQrVMqVTixIkTGDhwYKPvGThwoFp7ADh27FiT7QkhBkAjR5qewf79+5lUKmV79uxhV69eZe+++y7r0KEDy8/PZ4wxNn36dLZ06VJV+7NnzzKRSMQ2bNjA0tPTWXh4OBOLxezSpUvN6q+2tpaFh4ez2tparXw/mkI5NYtyapamc/JeiBhjbMuWLczFxYVJJBLm7+/Pzp07p/rasGHD2IwZM9TaHzhwgHXr1o1JJBLWq1cv9uOPP+o4MSFEk3i/jogQQujue0II76gQEUJ4R4WIEMI7KkSEEN61u0K0bds2uLm5wcjICP3790dycjLfkZ4oMjISHMchJCSE7yhqFAoFVq5ciS5dusDY2BjPPfccVq9eDb7PfZw5cwZjxoyBk5MTOI7DwYMHVV+Ty+VYsmQJevfuDVNTUzg5OSEoKAi5ubofV/tJOR9JT0/H2LFjYWlpCVNTU/j5+SEnJ0enOSMiIuDn5wdzc3PY29tj3LhxyMjIUGtTW1uL4OBg2NjYwMzMDBMmTGhw0fHTtKtC9GjIkfDwcKSmpsLLywuBgYEoLCzkO1qjzp8/j3/961/o00f/ZuxYt24dtm/fjq1btyI9PR3r1q1DVFQUtmzZwmuuqqoqeHl5Ydu2bQ2+Vl1djdTUVKxcuRKpqan4z3/+g4yMDIwdO1avcgJAVlYWhgwZAk9PT5w6dQoXL17EypUrYWSkuzv+AeD06dMIDg7GuXPncOzYMcjlcowcOVJtSJJFixbh0KFDiI+Px+nTp5Gbm4vx48e3rCOeLx/QKX9/fxYcHKx6rVAomJOTE4uIiOAxVeMePHjAunbtyo4dO8aGDRvGFi5cyHckNa+99hqbNWuW2rLx48ezadOm8ZSoIQANRnb4q+TkZAaAZWdn6yZUIxrLOWXKFPbWW2/xE+gJCgsLGQB2+vRpxhhjZWVlTCwWs/j4eFWb9PR0BoAlJSU1e73tZotIJpMhJSUFAQEBqmUCgQABAQFISkriMVnjgoOD8dprr6nl1SeDBg3CiRMnkJmZCQD4888/8euvv2LUKMMaGra8vBwcx6FDhw58R1FRKpX48ccf0a1bNwQGBsLe3h79+/dvdPdN18rLywFAddd9SkoK5HK52t+pp6cnXFxcWvS5avN33z/SmiFH+LJ//36kpqbi/PnzfEdp0tKlS1FRUQFPT08IhUIoFAqsXbsW06ZN4ztas9XW1mLJkiWYOnWqXt3pXlhYiMrKSkRGRmLNmjVYt24dEhMTMX78eJw8eRLDhg3jJZdSqURISAgGDx6sGg01Pz8fEomkQSF/0tA8jWk3hchQ3LlzBwsXLsSxY8d0fjygJQ4cOICvv/4a+/btQ69evZCWloaQkBA4OTlhxowZfMd7KrlcjsmTJ4Mxhu3bt/MdR82jsaxff/11LFq0CADg7e2N3377DTt27OCtEAUHB+Py5cv49ddfNb7udlOIWjPkCB9SUlJQWFgIHx8f1TKFQoEzZ85g69atqKur04sxhsPCwrB06VK8+eabAIDevXsjOzsbERERel+IHhWh7Oxs/Pzzz3q1NQQ8/FsViUTo2bOn2vIePXpopQg0x7x583D48GGcOXMGnTv/bwhjR0dHyGQylJWVqW0VtfRz1W6OEbVmyBE+jBgxApcuXUJaWprq0a9fP0ybNg1paWl6UYSAh2eg/joollAo1PtJBB8VoevXr+P48eOwsbHhO1IDEokEfn5+DU6TZ2ZmwtXVVadZGGOYN28eEhIS8PPPP6NLly5qX/f19YVYLFb7XGVkZCAnJ6dFn6t2s0UEAKGhoZgxYwb69esHf39/bNq0CVVVVXo1e6q5uXmD2UhMTU1hY2OjV7OUjBkzBmvXroWLiwt69eqFCxcuYOPGjZg1axavuSorK3Hjxg3V61u3biEtLQ3W1tbo2LEjJk6ciNTUVBw+fBgKhUJ1HMPa2hoSie5mpX1SThcXF4SFhWHKlCkYOnQoXnzxRSQmJuLQoUM4deqUzjICD3fH9u3bh++//x7m5uaqn5elpSWMjY1haWmJ2bNnIzQ0FNbW1rCwsMD8+fMxcOBADBgwoPkdafjsnt570pAj+kofT99XVFSwhQsXMhcXF2ZkZMTc3d3ZP/7xD1ZXV8drrpMnTzY6yPuMGTPYrVu3mhwE/uTJk3qT85Fdu3YxDw8PZmRkxLy8vNjBgwd1mpEx1uTP64svvlC1qampYe+//z6zsrJiJiYm7I033mB5eXkt6oeGASGE8K7dHCMihOgvKkSEEN5RISKE8I4KESGEd1SICCG8o0JECOEdFSJCCO+oEBFCeEeFiBDCOypEpFUYY3j33XdhbW0NjuOQlpamk37XrFnTsnuYiEGgQkRaJTExEXv27MHhw4eRl5ensxty//zzT3h7e2tl3fo6UUF7QIWItEpWVhY6duyIQYMGwdHRESJR6wZykMlkLWqvrUKkzxMVtAdUiNqg5ORkDB8+HMbGxvD09MQff/yBmJgYjc1W8fbbb2P+/PnIyckBx3Fwc3MDANTV1WHBggWwt7eHkZERhgwZ0mC42+HDh2PevHkICQmBra0tAgMDm+wnJSUFQ4cOhbGxMfr27Yvff/8dWVlZGi9ElZWVmDZtGmJjY2FlZfXU9sOHD8f8+fMREhICKysrODg4IDY2VjWkjLm5OTw8PPDf//5XoznbNE0OGUD4l5SUxIyMjFhUVBTLzMxk48aNY2PGjGHu7u4sNTVVI32UlZWxTz75hHXu3Jnl5eWxwsJCxhhjCxYsYE5OTuzIkSPsypUrbMaMGczKyoqVlJSo3jts2DBmZmbGwsLC2LVr19i1a9ca7SM9PZ2Zm5uzFStWsBs3brBvv/2WOTo6MoFAwKqqqtTarl27lpmamj7x8aRZOoKCglhISIgq39OGXBk2bBgzNzdnq1evZpmZmWz16tVMKBSyUaNGsZiYGJaZmcnmzp3LbGxsGmQljaNC1MYMHDiQTZ8+XfU6Li6OCQQC9sYbbzTr/XPmzGHe3t5szZo1as//Kjo6mrm6uqpeV1ZWMrFYzL7++mvVMplMxpycnFhUVJRq2bBhw1jfvn2fmuOll15S+z4YY2zixImse/fuDdqWlJSw69evP/Ehl8sb7eebb75hzz//PKupqVHla04hGjJkiOp1fX09MzU1Vcubl5fX4il12rN2NUJjW3f37l0kJSVhw4YNqmUikQiMMXz88cdPff/FixeRk5ODCxcu4OLFi/jwww9x4cKFZvWdlZUFuVyOwYMHq5aJxWL4+/sjPT1dra2vr+8T1/VoLOnU1FS15WKxuNHdMmtra9X0Ni3xLBMVPH4sSSgUwsbGBr1791YtezRbjL5O3qlv6BhRG/LoA//4wPsZGRnw9/dX+5Bcv34dr732Gnx9fTF06FAUFhbi6tWrGDVqFC5fvgxra2vV80GDBmk8p6mp6RO/npaWBpFIpJYZAC5cuNBoIfr0009hZmb2xEdjUzU/PlGBSCSCSCTC6dOnsXnzZohEIigUiiYzisVitdccx6kt4zgOAPR+DG99QVtEbUh5eTmEQqHqQ1BaWooNGzbAy8tL1aaurg7vv/8+vvjiC3Tu3Bk7duxATEwMVqxYgalTp2LAgAGYOHEiFi9erHreHM899xwkEgnOnj2rGuBdLpfj/PnzLT4dLhAIoFQqIZPJVGfjjhw5gmvXrjVaiP7+979j8uTJT1ynk5NTg2WPJip43MyZM+Hp6YklS5bozUQF7QEVojbE29sbCoUCUVFRmDRpEhYuXAg3NzdcvXoV2dnZcHV1xcGDB3HlyhWMHj0awMPC9PbbbwMALl26hHfeeafB8+YwNTXF3LlzERYWphoAPioqCtXV1Zg9e3aLvo9HM0OEhYXhgw8+wOXLlzF37lzV9/hXrd01M5SJCtoD2jVrQzw8PPDJJ5/gn//8J/r27QsnJyf89NNP6NSpE1555RUADwvMZ599ppqqKD09HUuWLAHwcJeta9euDZ43V2RkJCZMmIDp06fDx8cHN27cwNGjR5t1SvxxTk5O2LlzJ3744Qf06tULn332GYKCguDg4KBXc9ARzaHB89uZrVu34o8//sCePXsAPDxA3adPHxQXF2PkyJFITU1Ve06ILtAWUTszc+ZMlJWVwdPTE15eXvjqq68APNxSerQ78vhzQnSBtogIIbyjLSJCCO+oEBFCeEeFiBDCOypEhBDeUSEihPCOChEhhHdUiAghvKNCRAjhHRUiQgjvqBARQnhHhYgQwrv/A/LXgr5Ab2WoAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 300x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Only run if anisotropy has been computed\n",
    "create_fig3a(kaq, aniso)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "timflow (3.13.11)",
   "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.11"
  },
  "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
}