{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test line-sink discharge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (4, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of one line-sink with uniform strength to a row of wells with constant strength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=20)\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.LineSinkBase(ml1, x1=-10, y1=-10, x2=10, y2=10, Qls=1000)\n",
    "ml1.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"head at center of line-sink:\", ml1.head(ls1.xc[0], ls1.yc[0]))\n",
    "print(\"discharge of line-sink:\", ls1.discharge())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml2 = tfs.ModelMaq(kaq=20)\n",
    "rf2 = tfs.Constant(ml2, xr=0, yr=20, hr=30)\n",
    "N = 20\n",
    "d = 20 / N\n",
    "xw = np.arange(-10 + d / 2, 10, d)\n",
    "yw = np.arange(-10 + d / 2, 10, d)\n",
    "for i in range(N):\n",
    "    tfs.Well(ml2, xw[i], yw[i], Qw=1000 / N)\n",
    "ml2.solve(silent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.plots.contour(\n",
    "    [-20, 20, -20, 20], 50, [0], np.arange(20, 31, 1), color=\"C0\", labels=False\n",
    ")\n",
    "ml2.plots.contour(\n",
    "    [-20, 20, -20, 20],\n",
    "    50,\n",
    "    [0],\n",
    "    np.arange(20, 31, 1),\n",
    "    color=\"C1\",\n",
    "    labels=False,\n",
    "    newfig=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Head line-sink of higher order with control point on the line-sink"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(ml1, -10, 0, 10, 0, 20, order=7, layers=0)  # set control point at y=0.5\n",
    "ml1.solve()\n",
    "\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Head line-sink of higher order with control point $\\Delta y$ off line-sink\n",
    "Useful when line-sink is used to simulate horizontal well with finite radius or ditch with finite width.\n",
    "Note that this means that the head inside this distance $\\Delta y$ is essentially meaningless (but a value is returned)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(\n",
    "    ml1, -10, 0, 10, 0, 20, order=7, layers=0, dely=0.5\n",
    ")  # set control point at 0.5 from line-sink\n",
    "ml1.solve()\n",
    "\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of head line-sink with row of well with specified head. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(ml1, -10, -10, 10, 10, 20, order=7, layers=0)\n",
    "ml1.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml2 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf2 = tfs.Constant(ml2, xr=0, yr=20, hr=30)\n",
    "N = 50\n",
    "d = 20 / N\n",
    "xw = np.arange(-10 + d / 2, 10, d)\n",
    "yw = np.arange(-10 + d / 2, 10, d)\n",
    "for i in range(N):\n",
    "    tfs.HeadWell(ml2, xw[i], yw[i], 20, layers=0)\n",
    "ml2.solve(silent=True)\n",
    "Qwell = 0\n",
    "for i in range(N):\n",
    "    Qwell += ml2.elementlist[i + 1].discharge()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"discharge of line-sink:\", ls1.discharge())\n",
    "print(\"discharge of wells:\", Qwell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = ml1.plots.contour(\n",
    "    [-20, 20, -20, 20], 50, [0], np.arange(20, 31, 1), color=\"C0\", labels=False\n",
    ")\n",
    "ml2.plots.contour(\n",
    "    [-20, 20, -20, 20],\n",
    "    50,\n",
    "    [0],\n",
    "    np.arange(20, 31, 1),\n",
    "    color=\"C1\",\n",
    "    labels=False,\n",
    "    ax=ax,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-100, 100, 100)\n",
    "h1 = ml1.headalongline(x, 0)\n",
    "h2 = ml2.headalongline(x, 0)\n",
    "plt.figure(figsize=(8, 2))\n",
    "for ilay in range(2):\n",
    "    plt.plot(x, h1[ilay], label=f\"line-sink layer {ilay}\")\n",
    "    plt.plot(x, h2[ilay], label=f\"row of wells layer {ilay}\")\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher-order head line-sink with resistance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=3, z=[10, 0])\n",
    "ls1 = tfs.River(ml, -10, 0, 10, 0, hls=2, wh=1, res=0.2, order=2)\n",
    "rf = tfs.Constant(ml, 0, 20, 4)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"inflow at each control-point\")\n",
    "for i in range(3):\n",
    "    print(\n",
    "        i,\n",
    "        (ml.head(ls1.xc[i], ls1.yc[i]) - ls1.hc[i]) * ls1.wh / ls1.res,\n",
    "        np.sum(ls1.strengthinf[i] * ls1.parameters[:, 0]),\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"head inside at each control point\")\n",
    "for i in range(3):\n",
    "    print(\"icp, head: \", i, ls1.headinside(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Higher order head line-sink with resistance in one layer of multiple layers system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[20, 10, 10, 0], c=[1000])\n",
    "lslayer = 0\n",
    "order = 2\n",
    "ls1 = tfs.River(ml, -10, 0, 10, 0, hls=-2, order=order, wh=1, res=2, layers=[lslayer])\n",
    "rf = tfs.Constant(ml, 0, 20, 2)\n",
    "ml.solve()\n",
    "\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[1, 2], z=[20, 12, 10, 0], c=[1000])\n",
    "order = 2\n",
    "ls1 = tfs.River(ml, -10, 0, 10, 0, hls=-2, order=order, wh=1, res=2, layers=[0, 1])\n",
    "rf = tfs.Constant(ml, 0, 2000, 2)\n",
    "ml.solve()\n",
    "\n",
    "print(\"inflow at each control-point\")\n",
    "for i in range(order + 1):\n",
    "    print(\n",
    "        i,\n",
    "        (ml.head(ls1.xc[i], ls1.yc[i]) - ls1.hc[i])[lslayer] * ls1.wh / ls1.res,\n",
    "        np.sum(ls1.strengthinf[i] * ls1.parameters[:, 0]),\n",
    "    )\n",
    "print(\"head inside at each control point\")\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying heads along line-sinks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Give one value that is applied at all control points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(ml1, -10, 0, 10, 0, hls=20, order=2, layers=[0])\n",
    "ml1.solve()\n",
    "print(ml1.headalongline(ls1.xc, ls1.yc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Give `order + 1` values, which is applied at the `order + 1` control points. This may not be so useful, as the user needs to know where those control points are. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(ml1, -10, 0, 10, 0, hls=[20, 19, 18], order=2, layers=[0])\n",
    "ml1.solve()\n",
    "print(ml1.headalongline(ls1.xc, ls1.yc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Give two values, the head at the beginning and end of the line-sink. Linear interpolation in between. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.River(ml1, -10, 0, 10, 0, hls=[19, 20], order=2, layers=[0])\n",
    "ml1.solve()\n",
    "print(ml1.headalongline(ls1.xc, ls1.yc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ditch\n",
    "Specify total discharge of ditch. Uniform head inside ditch is computed. \n",
    "\n",
    "Single ditch in one layer of two-aquifer system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.Ditch(ml1, -10, -10, 10, 10, Qls=1000, order=2, layers=[0])\n",
    "ml1.solve()\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single ditch with resistance, in both layers of two-aquifer system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.Ditch(ml1, -10, -10, 10, 10, Qls=1000, order=2, res=0.1, layers=[0, 1])\n",
    "ml1.solve()\n",
    "for icp in range(ls1.ncp):\n",
    "    print(\"icp, head inside: \", icp, ls1.headinside(icp))\n",
    "print(\"discharge: \", ls1.discharge())\n",
    "print(\"total discharge: \", np.sum(ls1.discharge()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Head line-sink string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.RiverString(\n",
    "    ml1, xy=[(-10, 0), (0, 0), (10, 0), (10, 10)], hls=20, order=5, layers=[0]\n",
    ")\n",
    "ml1.solve()\n",
    "ml1.plots.contour([-20, 20, -20, 20], 41, [0], 40, labels=False)\n",
    "for ils in range(ls1.nls):\n",
    "    print(\"line-sink, head inside: \", ils, ls1.lslist[ils].headinside())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Head line-sink string with linearly varying head."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# linearly varying head along line-sink\n",
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=30)\n",
    "ls1 = tfs.RiverString(\n",
    "    ml1, xy=[(-10, 0), (0, 0), (10, 0), (10, 10)], hls=[20, 22], order=5, layers=[0]\n",
    ")\n",
    "ml1.solve()\n",
    "ml1.plots.contour([-20, 20, -20, 20], 41, [0], 40)\n",
    "for ils in range(ls1.nls):\n",
    "    print(\"line-sink, spec head, head inside: \", ils, ls1.lslist[ils].hc)\n",
    "    for icp in range(ls1.lslist[ils].ncp):\n",
    "        print(\"icp, head \", ls1.lslist[ils].headinside(icp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xls1 = np.linspace(-10, 10, 50)\n",
    "yls1 = np.linspace(0, 0, 50)\n",
    "hls1 = ml1.headalongline(xls1, yls1)\n",
    "plt.figure()\n",
    "plt.plot(xls1, hls1[0])\n",
    "xls2 = np.linspace(10, 10, 50)\n",
    "yls2 = np.linspace(0, 10, 50)\n",
    "hls2 = ml1.headalongline(xls2, yls2)\n",
    "plt.plot(10 + yls2, hls2[0])\n",
    "plt.xlabel(\"length along line-sink string (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Head line-sink string with resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=200, hr=2)\n",
    "ls1 = tfs.RiverString(\n",
    "    ml1,\n",
    "    xy=[(-10, 0), (0, 0), (10, 0), (10, 10)],\n",
    "    hls=1,\n",
    "    res=2,\n",
    "    wh=5,\n",
    "    order=5,\n",
    "    layers=[0],\n",
    ")\n",
    "ml1.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ils in range(ls1.nls):\n",
    "    for icp in range(ls1.lslist[ils].ncp):\n",
    "        print(\"ils, icp, head inside:\", ils, icp, ls1.lslist[ils].headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate total discharge and per linesink (using two different methods), check they\n",
    "are all equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qtot = ls1.discharge()\n",
    "Qls_sum = np.sum(ls1.discharge_per_linesink(), axis=1)\n",
    "Qper_linesink_sum = np.sum([e.discharge() for e in ls1.lslist], axis=0)\n",
    "\n",
    "\n",
    "assert np.allclose(Qtot, Qls_sum)\n",
    "assert np.allclose(Qtot, Qper_linesink_sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot head in aquifer along line-sink string\n",
    "# (note there is a resistance, so this is not constant)\n",
    "xls1 = np.linspace(-10, 10, 50)\n",
    "yls1 = np.linspace(0, 0, 50)\n",
    "hls1 = ml1.headalongline(xls1, yls1)\n",
    "plt.figure()\n",
    "plt.plot(xls1, hls1[0])\n",
    "xls2 = np.linspace(10, 10, 50)\n",
    "yls2 = np.linspace(0, 10, 50)\n",
    "hls2 = ml1.headalongline(xls2, yls2)\n",
    "plt.plot(10 + yls2, hls2[0])\n",
    "plt.xlabel(\"length along string (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ditch string\n",
    "Ditch string without resistance, screened in one layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=1)\n",
    "ls1 = tfs.DitchString(\n",
    "    ml1, xy=[(-10, 0), (0, 0), (10, 0)], Qls=100, res=0, order=5, layers=[0]\n",
    ")\n",
    "\n",
    "ml1.solve()\n",
    "print(\"discharge:\", ls1.discharge())\n",
    "\n",
    "for ils in range(ls1.nls):\n",
    "    for icp in range(ls1.lslist[ils].ncp):\n",
    "        print(\"ils, icp, head inside:\", ils, icp, ls1.lslist[ils].headinside(icp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.plots.contour([-20, 20, -20, 20], 81, [0], 20, labels=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ditch in multiple layers with resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=1)\n",
    "ls1 = tfs.DitchString(\n",
    "    ml1, xy=[(-10, 0), (0, 0), (10, 0)], Qls=100, res=0.2, order=5, layers=[0, 1]\n",
    ")\n",
    "\n",
    "ml1.solve()\n",
    "print(\"discharge:\", ls1.discharge())\n",
    "\n",
    "for ils in range(ls1.nls):\n",
    "    for icp in range(ls1.lslist[ils].ncp):\n",
    "        print(\"ils, icp, head inside:\", ils, icp, ls1.lslist[ils].headinside(icp))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ditch in different layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tfs.ModelMaq(kaq=[20, 10], z=[20, 12, 10, 0], c=[100])\n",
    "rf1 = tfs.Constant(ml1, xr=0, yr=20, hr=1)\n",
    "ls1 = tfs.DitchString(\n",
    "    ml1,\n",
    "    xy=[(-10, -10), (-5, -5), (5, 5), (10, 10)],\n",
    "    Qls=100,\n",
    "    wh=2,\n",
    "    res=0,\n",
    "    order=5,\n",
    "    layers=[0, 1, 0],\n",
    ")\n",
    "ml1.solve()\n",
    "\n",
    "for ils in range(ls1.nls):\n",
    "    for icp in range(ls1.lslist[ils].ncp):\n",
    "        print(\"ils, icp, head inside:\", ils, icp, ls1.lslist[ils].headinside(icp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.plots.contour([-20, 20, -20, 20], 81, layers=[0], levels=20, labels=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.plots.contour([-20, 20, -20, 20], 81, layers=[1], levels=20, labels=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Angle well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(\n",
    "    kaq=1,\n",
    "    z=np.arange(10, -0.1, -0.2),\n",
    "    kzoverkh=0.1,\n",
    "    topboundary=\"semi\",\n",
    "    topres=0,\n",
    "    topthick=2,\n",
    "    hstar=7,\n",
    ")\n",
    "xy = list(zip(np.linspace(-10, 10, 21), np.zeros(21), strict=False))\n",
    "ls1 = tfs.DitchString(\n",
    "    ml, xy=xy, Qls=100, wh=2, res=5, order=2, layers=np.arange(10, 30, 1)\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.vcontour(win=[-20, 20, 0, 0], n=100, levels=20, layout=True, horizontal_axis=\"x\")\n",
    "plt.xlabel(\"x (m)\");"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}