{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ad954d5f-faff-45be-9291-d7875b453349",
   "metadata": {},
   "source": [
    "# Simulate balloon angioplasty of a coronary artery "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3da803f6-96cb-4c5f-b926-a06f4700c8b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymecht as pmt\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb7cf411-e727-4d66-b876-73b17ff2b829",
   "metadata": {},
   "source": [
    "Balloon angioplasty is a procedure wherein a balloon is used to inflate an artery. If we assume that both artery and balloon are perfect cylinders, the procedure can be simulated using pymecht using a `LayeredTube`. Thus, we start by creating artery and balloon separately (both made of single layer) of neo-Hookean material. The parameters of artery and balloon are different -- both material and geometric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6df3c335-fd14-4f97-8073-4ea00675f8a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "An object of type TubeInflationwith radius as input, pressure as output, and the following material\n",
      "Material model with 1 component:\n",
      "Component1: NH\n",
      " An object of type TubeInflationwith radius as input, pressure as output, and the following material\n",
      "Material model with 1 component:\n",
      "Component1: NH\n",
      "\n",
      "------------------------------------------------------------------\n",
      "Keys              Value       Fixed?      Lower bound Upper bound \n",
      "------------------------------------------------------------------\n",
      "Ri                1.50        No          0.50        1.50        \n",
      "thick             1.00        No          0.00        1.00        \n",
      "omega             0.00        No          0.00        0.00        \n",
      "L0                1.00        No          1.00        1.00        \n",
      "lambdaZ           1.00        No          1.00        1.00        \n",
      "mu_0              0.23        No          1.00e-04    1.00e+02    \n",
      "------------------------------------------------------------------\n",
      " ------------------------------------------------------------------\n",
      "Keys              Value       Fixed?      Lower bound Upper bound \n",
      "------------------------------------------------------------------\n",
      "Ri                0.58        No          0.50        1.50        \n",
      "thick             0.10        No          0.00        1.00        \n",
      "omega             0.00        No          0.00        0.00        \n",
      "L0                1.00        No          1.00        1.00        \n",
      "lambdaZ           1.00        No          1.00        1.00        \n",
      "mu_0              1.72        No          1.00e-04    1.00e+02    \n",
      "------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "artery_mat = pmt.MatModel('nh')\n",
    "balloon_mat = pmt.MatModel('nh')\n",
    "artery = pmt.TubeInflation(artery_mat,force_measure='pressure')\n",
    "balloon = pmt.TubeInflation(balloon_mat,force_measure='pressure')\n",
    "print(artery,balloon)\n",
    "artery_params = artery.parameters\n",
    "balloon_params = balloon.parameters\n",
    "artery_params.set('Ri',1.5)\n",
    "artery_params.set('thick',1)\n",
    "balloon_params.set('Ri',0.58)\n",
    "balloon_params.set('thick',0.1)\n",
    "artery_params.set('mu_0',0.227)\n",
    "balloon_params.set('mu_0',1.72)\n",
    "print(artery_params,balloon_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d14249fe-8cca-436a-98a8-26f5fb0acbb2",
   "metadata": {},
   "source": [
    "Next, we combine them into a `LayeredTube` with balloon on the inside and artery on the outside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5f37d484-ec9b-4f76-9373-a1e6558843f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "An object of type LayeredTubewith 2 layers:\n",
      "Layer1: An object of type TubeInflationwith radius as input, pressure as output, and the following material\n",
      "Material model with 1 component:\n",
      "Component1: NH\n",
      "Layer2: An object of type TubeInflationwith radius as input, pressure as output, and the following material\n",
      "Material model with 1 component:\n",
      "Component1: NH\n",
      " ------------------------------------------------------------------\n",
      "Keys              Value       Fixed?      Lower bound Upper bound \n",
      "------------------------------------------------------------------\n",
      "Ri_layer0         0.58        No          0.50        1.50        \n",
      "thick_layer0      0.10        No          0.00        1.00        \n",
      "omega_layer0      0.00        No          0.00        0.00        \n",
      "L0_layer0         1.00        No          1.00        1.00        \n",
      "lambdaZ_layer0    1.00        No          1.00        1.00        \n",
      "mu_0_layer0       1.72        No          1.00e-04    1.00e+02    \n",
      "Ri_layer1         1.50        No          0.50        1.50        \n",
      "thick_layer1      1.00        No          0.00        1.00        \n",
      "omega_layer1      0.00        No          0.00        0.00        \n",
      "L0_layer1         1.00        No          1.00        1.00        \n",
      "lambdaZ_layer1    1.00        No          1.00        1.00        \n",
      "mu_0_layer1       0.23        No          1.00e-04    1.00e+02    \n",
      "------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "combined = pmt.LayeredTube(balloon,artery)\n",
    "combined_params = combined.parameters\n",
    "print(combined,combined_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb1961de-a98f-40ea-83e2-c19cbe77d705",
   "metadata": {},
   "source": [
    "In order to simulate their interaction, one consideration is to define when they are in contact (i.e., they bear the load together) and when they are not. For this, we have to calculate the outer radius of the balloon and see if it lower than the reference inner radius of the artery. If yes, then the balloon and artery are separate and only balloon bears the applied pressure, otherwise they are in contact and must be simulated together. Thus, we vary the balloon radius over a range and keep track of whether they are in contact or not by using the `outer_radius` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "271d1f42-de8e-46e3-a261-2685857c34ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "r_balloon = np.linspace(0.58,2,100)\n",
    "p_all = []\n",
    "contact = []\n",
    "r_all = []\n",
    "p_art = []\n",
    "for ri in r_balloon:\n",
    "    rio = balloon.outer_radius([ri],balloon_params)[0]\n",
    "    r_all.append(rio)\n",
    "    if rio >= artery_params['Ri'].value:\n",
    "        contact.append(True)\n",
    "        p_all.append(combined.disp_controlled([ri])[0])\n",
    "        p_art.append(artery.disp_controlled([rio])[0])\n",
    "    else:\n",
    "        contact.append(False)\n",
    "        p_all.append(balloon.disp_controlled([ri])[0])\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdb9807f-2706-4aa0-8a2a-672cfeec6f7f",
   "metadata": {},
   "source": [
    "Lastly, we plot the results separating the points when the two are in contact from the points when they are not. Contact pressure is simply the pressure supported by the artery (i.e., the blue curve in the plot below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "897ebd39-24ce-4b91-9043-20b21109a790",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAADTCAYAAABa1gxPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsDElEQVR4nO2deXwUVbbHv8cYiAIim+OCQBhZxATCKrLJIggDA6KOICDghqiDo/Pc9Sni6PMN6CgiuKKCElAERESFJyCCG4lEkX0LCqKyh0WQJOf9cSuhCVk6pCvdSc7386lPV9+6t+p0defk3lP3np+oKoZhGH5wSrgNMAyj9GIOxjAM3zAHYxiGb5iDMQzDN8zBGIbhG+ZgDMPwDV8djIh0F5G1IrJBRO7P5fhwEVkhIikiskREGnnldUTkd688RURe9NNOwzD8QfyaByMiUcA6oCuwFVgGXKuqqwLqnKGqad5+b+A2Ve0uInWAOaoa54txhmEUC372YFoBG1R1k6r+AUwF+gRWyHIuHhUAm/VnGKWIU30893nATwHvtwIX56wkIrcD/wTKAZ0DDsWKyHIgDXhYVT/P72LVq1fXOnXqFNVmwzBOguTk5J2qWiNnuZ8OJihU9QXgBREZADwMDAG2A7VUdZeINAdmichFOXo8iMgwYBhArVq1SEpKKmbrDcMAEJEtuZX7OUTaBpwf8L6mV5YXU4ErAFT1iKru8vaTgY1A/ZwNVPVlVW2hqi1q1DjBeRqGEWb8dDDLgHoiEisi5YD+wOzACiJSL+BtT2C9V17DCxIjInWBesAmH201DMMHfBsiqWq6iPwd+ASIAiaq6koRGQUkqeps4O8ichlwFNiDGx4BdABGichRIBMYrqq7/bLVMAx/8O0xdXHTokULtRiMYYQHEUlW1RY5y20mr2EY2WRmZvL9998zYcIEBg0axJYtucZugybsT5EMwwgfR48eJSkpicWLF7N48WKWLl3Kvn37ADj77LMZNmwYtWvXPunzm4MxjDJElkNZuHAhixYtYunSpRw6dAiABg0acM0119CuXTvatWtHbGwsIlKk65mDMYxSjKqyatUq5s+fz//93//x2WefceDAAQDi4uK48cYb6dChAx06dOCss84K+fXNwRhGKWPnzp3Mmzcve9u+fTsA9erVY9CgQXTu3JmOHTtSHHPHzMEYRgknMzOT5ORk5s6dy9y5c1m2bBmqStWqVenatSvdunXjsssuo1atWsVumzkYwyiBHDx4kPnz5/PBBx/w4Ycf8uuvvyIiXHzxxYwcOZLu3bvTvHlzoqKiwmqnORjDKCHs2LGD2bNnM2vWLObPn8+RI0eoXLkyPXr0oGfPnnTv3p3q1auH28zjMAdjGBHMzz//zIwZM5g+fTqff/45mZmZ1K5dm1tuuYU+ffrQvn17oqOjw21mnpiDMYwI45dffmH69OlMmzaNpUuXoqo0atSIhx56iL59+5KQkFDkx8fFhTkYw4gA9u7dy3vvvUdiYiILFy4kMzOTuLg4HnvsMa6++mouvPDCcJt4UpiDMYww8ccff/Dxxx8zadIk5syZw5EjR7jgggt46KGH6N+/P40aNQq3iUXGHIxhFDMpKSm8/vrrTJkyhZ07d3LWWWdxyy23MHDgQFq2bFlihj/B4KuDEZHuwHO4dA2vqupTOY4PB24HMoADwLCspOAi8gBwo3fsDlX9xE9bDcNP9u7dy5QpU3j11VdZvnw55cqVo0+fPgwZMoRu3bpFdKC2SKiqLxvOqWwE6uLy7X4HNMpR54yA/d7Ax95+I69+eSDWO09Uftdr3ry5GkYkkZmZqV988YUOGTJETzvtNAU0ISFBn3/+ed21a1e4zQspuBxPJ/xd+tmDyVYVABCRLFWBbNkSzVtVoA8wVVWPAJtFZIN3vi99tNcwQsLBgweZMmUK48ePJyUlhYoVKzJ48GBuvvlmmjdvHm7zipVIVRU4D/gqR9vz/DHTMELDpk2bGDduHBMnTmTfvn00btyYF198kQEDBlCpUqVwmxcWwh7k1dxVBYIip6qAYRQ3qsqSJUt45plneP/994mKiuKqq65ixIgRtGnTplQFbE8GPx3MyagKTChMW1V9GXgZXMrMohhrGIUhPT2dGTNmMHr0aJKSkqhatSoPPPAAt912G+edZ53tLCJSVcCr119EyotILE5V4BsfbTWMoDh8+DAvvvgiDRo0oF+/fuzbt48JEybw008/8cQTT5hzyUFEqgp49d7BBYTTgdtVNcMvWw2jIA4cOMCLL77I008/zS+//EKrVq0YPXo0ffr0CfuK5UjGVAUMIx/279/PuHHjePrpp9m1axeXXXYZDz74IB07dizz8ZVA8lIVCHuQ1zAikYMHDzJu3DhGjx7Nrl276NGjB4888gitW7cOt2klCnMwhhHAkSNHeOmll3jyySf59ddf6dGjByNHjqRVq1bhNq1EYg7GMICMjAzefvttHnnkEbZs2ULHjh2ZMWMGbdq0CbdpJRoTXjPKPPPnz6d58+YMGTKEatWqMW/ePBYsWGDOJQSYgzHKLKtXr6Znz55069aNtLQ0EhMTWbZsGV27drUAbogwB2OUOfbs2cOdd95JfHw8S5YsYfTo0axevZr+/ftzyin2JxFKLAZjlBkyMzOZOHEiDzzwALt372bYsGGMGjWqWPSByirmro0ywfLly2nTpg0333wzDRs2JDk5mQkTJphz8RlzMEapZv/+/dx11120aNGCzZs3M3nyZBYvXkxCQkK4TSsTmIMxSg+q8Pvvx97v3cunn37Kc889x/Dhw1m7di2DBg2yAG4xYjEYI7JRhV274OefYft2yMyEHj3csfvvh6Qk+O032LEDdu6E9u1hwQJ3fNUq+vTpw8qVK0tsVv6SjjkYI7zs3w+pqbBli9sOHoR773XHBg2Cd9+FP/44Vr9+/WMOZts2OHQI6taFiy+GGjUg0JHExSEi5lzCiDkYw3/27IG1a2H9eti8Gf77v0EEbr8dxo8/vm7lynDPPe54+/Zw7rlQsyacc47bP/fcY3UnT87/umecEfrPYhSKcKsK/BO4CZeSYQdwg6pu8Y5lACu8qj+qam8/bTVCwJ49sGIFtGoFMTHw0kvwyCNuCJPFKafAbbdB9erwl79ArVpQpw7Uru1ezzrLOReAW24Jx6cwQohvDkZEooAXgK64nLrLRGS2erIkHsuBFqp6SERuBf4N9POO/a6qCX7ZZ4SADRvg7bfh228hJQV+/NGVL1sGLVo4p9Grlxu21K8PDRpAbCyUK+fq9ezpNqPUEm5VgYUB9b8CBvloj3GyHD0Ky5fDkiXw1VeuB9KxI/z0Ezz2GDRsCG3buvLGjZ0zAeje3W1GmSXsqgIB3Ah8FPA+RkSScMOnp1R1VsgtNPJnxw4YOBCWLnXBVHC9kquvdvvt2kFaGlSsGD4bjYgmIoK8IjIIaAFcGlBcW1W3iUhdYIGIrFDVjTnamapAqNi2DebOhY8/dj2Q//kfqFrVPdW54Qbo0MH1UgKDrNHRbjOMPAi7qoCXk/ch4FJPaA0AVd3mvW4SkUVAU5zCIwF1TFWgqIwdC2++6eIoAOef74Y5AFFRrvdiGCdJ0DN5RaSKiFwkInVFJJh2wagKNAVeAnqr6m8B5VVEpLy3Xx1oS0DsxigC69fDf/7jJrABfPcdlC/veiwrVri5KI8+Gl4bjVJDvj0YEamME6e/Fqe8uAOIAf4kIl8B43MEarMJUlVgNFAReNebvp31OPpC4CURycQ5wadyPH0yCsO+fZCYCG+8AV9/7R4D9+oF9erByy+7noph+EC+qgIiMh+YBHygqntzHGsOXAesUNXX/DQyGExVIA+WLoWuXd0anbg4GDIErr0WTL/HCCEnpSqgql3zOZYMJIfANiOUpKfDzJmul3L11dCsGdx0E1x3nZubYgv9jGIk6CCviFTBKSzGZJWp6mI/jDJOgiNH3BDo3/+GTZtcr+Xqq+G001wg1zDCQFBBXhG5CViMi6c85r2O9M8so1DMnAkXXADDh7sp+DNmwEcfFdzOMHwm2KdI/wBaAltUtRPukfFev4wygiAz0/VawM1FqVUL5s1zM2379rXArRERBOtgDqvqYQARKa+qa4AG/pll5MuyZW7S26hR7n3Pnm4af9euFmMxIopgHcxWETkTmAXMF5H3gS1+GWXkwZ49boXxxRe7tAeNGrlyEXMsRkQSbJD3Zu8x9UgRWQhUBj72zSrjRObPd0+CduyAO++EkSMt34kR8RQ00e6vwETgqDfp7RpV/axYLDOO5+yzXea2uXPdo+cyyNGjR9m6dSuHDx8OtylllpiYGGrWrEl0kGvQCurBPAG0V9U1InIxLl/LpQW0MULFZ585h/K//wvx8W7SXBkeCm3dupVKlSpRp04dS9wdBlSVXbt2sXXrVmJjY4NqU1AMJt0L6KKqXwOVimijEQyZmfD449C5M8ya5WIvUKadC8Dhw4epVq2aOZcwISJUq1atUD3IgnowZ3lpLXN9r6rPFNJGoyD27nWxljlzXNLrCRMs30oA5lzCS2Hvf0EO5hWO77XkfG+EElXo1s1ljxs3zmWIsz8oowRT0FqkxwBEpKqq7g48JiLBDcKM4BFxKSgrVHAJnoyIIyoqivj4eFSVqKgoxo0bR5s2bfJtU7FiRQ4cOEBqaiq9evXihx9+KCZrC+aNN94gKSmJcePG+XL+YB9TfyAiPVQ1DUBELgTeBeJ8saqsMWeOExYbNuyY5o8RkZx22mmkpKQA8Mknn/DAAw/w2Wf2YDUvgp1o9yTOyVT00jRMJ4gE3SLSXUTWisgGEbk/l+P/FJFVIvK9iHwqIrUDjg0RkfXeNiTYD1TiePttuOIKeO01txLaKDGkpaVRpUoVAA4cOECXLl1o1qwZ8fHxvP/++/m2PXz4MNdffz3x8fE0bdqUhQsX5lv+xhtvcOWVV9K9e3fq1avHvVnidDn49NNPadq0KfHx8dxwww0c8ZaT1KlTh0cffTTbvjVr1hzXbv/+/cTGxnL06NHszxb4/mQJqgejqh+KSDQwDxeD6auq6/JrUxTZEhGpCjyKy9OrQLLXdk8hP19kM2kSDB3qMvS//z6cGhEpkksOHTueWHbNNS52deiQ013KydChbtu581jy8iwWLSrwkr///jsJCQkcPnyY7du3s8CTqY2JiWHmzJmcccYZ7Ny5k9atW9O7d+88g6IvvPACIsKKFStYs2YN3bp1Y926dXmWA6SkpLB8+XLKly9PgwYNGDFiBOeffywr7eHDhxk6dCiffvop9evXZ/DgwUyYMIE777wTgOrVq/Ptt98yfvx4xowZw6uvvprdtlKlSnTs2JEPP/yQK664gqlTp3LllVcGPd8lL/LtwYjI8yIyVkTGAp1xM3g3A3/3yvIjW7ZEVf8AsmRLslHVharqpavnK1zeXoDLgfmquttzKvOB0qV/8dZb7ofeubMbIlWy2HlJIGuItGbNGj7++GMGDx6MqqKqPPjggzRu3JjLLruMbdu28euvv+Z5niVLljBokBsENGzYkNq1a7Nu3bo8ywG6dOlC5cqViYmJoVGjRmzZcvxqnbVr1xIbG0t9TzZmyJAhLF58LKPKlVdeCUDz5s1JTU09waabbrqJ119/HYDXX3+d66+//iTv0jEK+peZM0VcYRJMFUW2JLe2J6RgK9GqAj//DJ06wezZcPrp4bamZJJfj+P00/M/Xr16UD2W/LjkkkvYuXMnO3bsYO7cuezYsYPk5GSio6OpU6dOyGccly9fPns/KiqK9EIOqbPa59W2bdu2pKamsmjRIjIyMoiLK3qINd8ejKq+md9W5Kt7BMiWjC5MO1V9WVVbqGqLGjVqhMocf8n6Yu+9Fz75xJxLCWbNmjVkZGRQrVo19u3bx1lnnUV0dDQLFy48oXeRk/bt2/P2228DsG7dOn788UcaNGiQZ3kwNGjQgNTUVDZs2ADA5MmTufTSwk28Hzx4MAMGDAhJ7wUKHiJ9ICJ/9eIvOY/VFZFRInJDHs0LK1vSO0C2JKi2JY41a5x86hdfuPcWcylxZMVgEhIS6NevH2+++SZRUVEMHDiQpKQk4uPjmTRpEg0bNsz3PLfddhuZmZnEx8fTr18/3njjDcqXL59neTDExMTw+uuv87e//Y34+HhOOeUUhg8fXqjPN3DgQPbs2cO1115bqHZ5kjV+zG0DzsYFXjfiZEjmAgtwcZj5QJ982p4KbAJicYoE3wEX5aiTpXVUL0d5Ve8aVbxtM1A1P1ubN2+uEc1vv6nWratao4bqpk3htqZEsmrVqnCbUOp59913ddCgQfnWye17wCmFnPB3WdBEu1+Ae4F7RaQOcA7wO7BOjwVn82p70rIlqrpbRB73nBrAKM0x0a9EkZ4Of/ubi7ssXOgE4A0jwhgxYgQfffQRc+fODdk5g+6jq2oqkFqYk6vqXFyvJ7DskYD9y/JpOxGXKqLkc999bmX05MnQunW4rTGMXHn++edDfs6glR2NkyQjw+k+//3vbvGiYZQhLMroN1FRTlUxIyPclhhGsVMYberTRMQSfQdLerpbW7RunVvEaE+MjDJIsLpIfwVS8PLwikiCiMzOt1FZ54kn4JVX4Ntvw22JYYSNYHswI3FT//cCqGoK7vGzkRspKfCvf8GAAdC/f7itMULMrFmzEJETFgwGsnfvXsaPH1+MVuXPyJEjGTNmTLFfN1gHc1RV9+Uo01AbUyo4etStMapWzSRbSymJiYm0a9eOxMTEXI+np6eflINRVTIzM0NhYsQQrINZKSIDgCgRqScizwNf+GhXyWX8ePjuO3jxRedkjFLFgQMHWLJkCa+99hpTp07NLl+0aBHt27end+/eNGrUiPvvv5+NGzeSkJDAPffcA8Do0aNp2bIljRs35tFHHwUgNTWVBg0aMHjwYOLi4nj88cezVz8DvPLKK9x1110n2JGYmEh8fDxxcXHcd9992eUVK1bkoYceokmTJrRu3fqEBZcbN26kWYAqxfr16497H2qCjTyOwE3nPwJMwU2e+5dfRpVobr4ZatRwOV4MX8ktW0NOevWCu+8+Vr+I2Rp4//336d69O/Xr16datWokJyfTvHlzAL799lt++OEHYmNjSU1N5YcffshOTjVv3jzWr1/PN998g6rSu3dvFi9eTK1atVi/fj1vvvkmrVu35sCBAzRp0oTRo0cTHR3N66+/zksvvXScDT///DP33XcfycnJVKlShW7dujFr1iyuuOIKDh48SOvWrXniiSe49957eeWVV3j44Yez2/75z3+mcuXKpKSkkJCQELJV03lRYA/Gy+vyoao+pKotve1h9aRkDQ9VpxV9+uku9mKUShITE+nvxdX69+9/3DCpVatWecp5zJs3j3nz5tG0aVOaNWvGmjVrWL9+PQC1a9emtTcBs2LFinTu3Jk5c+awZs0ajh49Snx8/HHnWrZsGR07dqRGjRqceuqpDBw4MDstQ7ly5ejVqxdQcFqGjIwMpk2bxgAff68F9mBUNUNEMkWkci5xGCOLOXPgH/9wCox//nO4rSkTFDbbQmD9k8nWsHv3bhYsWMCKFSsQETIyMhARRo92SQAqVKiQZ1tV5YEHHuCWW245rjw1NfWEdjfddBNPPvkkDRs2LHTvIjo6OjvJVV5pGa666ioee+wxOnfuTPPmzanm41A+2BjMAWCFiLyWlYAqiIRTZYejR10/vFw5KGl5aYygmT59Otdddx1btmwhNTWVn376idjYWD7//PMT6laqVIn9+/dnv7/88suZOHEiBw4cAGDbtm389ttvuV7n4osv5qeffmLKlCm5rmpu1aoVn332GTt37iQjI4PExMRCpWWIiYnh8ssv59Zbb/V1eATBx2BmeJuRGy+95CbUzZ4NRUwxaEQuiYmJxwVUwfUGEhMT6dev33Hl1apVo23btsTFxdGjRw9Gjx7N6tWrueSSSwA3FHrrrbeIiorK9VrXXHMNKSkp2Tl/AznnnHN46qmn6NSpE6pKz5496dOnTy5nyZuBAwcyc+ZMunXrVqh2hUXcSuuST4sWLTQpKWcCvmJg71644AJo3Bg+/dR0jHxk9erVXHjhheE2o1jo1asXd911F126dPHl/GPGjGHfvn08/vjjhW6b2/cgIsmq2iJn3aB6MCKymVzmvahq3QLadQeew6VreFVVn8pxvAPwLNAY6K+q0wOOZQArvLc/qmrvYGwtdt58E3btgjFjzLkYRWbv3r20atWKJk2a+OZc+vbty8aNG7MTlvtJsEOkQM8UA/wNlxQqT4JUFfgRGArcncspflfVhCDtCx8jRkDLluDjXAKj7HDmmWdmJ/n2i5kzZ/p6/kCCCvKq6q6AbZuqPgv0LKBZMKoCqar6PVAypy8ePQqnnAIFKPsZRlkl2MWOzQK2FiIynIJ7P0EpA+RDjIgkichXInJFHnYN8+ok7dixoxCnDgG//eYy082aVbzXNYwSRLBDpKcD9tNxme2uCbk1x1NbVbeJSF1ggYisUNWNgRVU9WXgZXBBXp/tOZ5nn3UpMMtI0NEwToZglR07ncS5i6QMoKrbvNdNIrKIYwnCw09amltzdNVVTiXAMIxcCXaI9A8ROUMcr4rItyJS0AP0ZUA9EYkVkXJAfyCoHDIiUkVEynv71YG2wKr8WxUjL78M+/a5XLtGmaJixYphue6zzz7LoUP55tnPk5SUlJAm8i4Mwc7kvUFV04BuQDXgOuCp/BqoajqQpSqwGngnS1VARHoDiEhLEdmKeyr1kois9JpfCCSJyHfAQuCpHE+fwsfRo/DMM9ClC7Q44bG/YfhCSXUw+eoi6TGdou+91+dwwvcAy4NpW1xbseoiffWVanJy8V3PUNXI0EWqUKGCqqouXLhQL730Ur3qqqu0QYMGOmDAAM3MzDyh/vr167VLly7auHFjbdq0qW7YsEEzMzP17rvv1osuukjj4uJ06tSp+Z7zueee0+joaI2Li9OOHTuqqurw4cO1efPm2qhRI33kkUeyr/fNN9/oJZdcoo0bN9aWLVvq3r179fzzz9fq1atrkyZNsq9VFEKmixRAsojMw2Wxe0BEKlFSHy2Hgovzk9g2ioM777wzOxVCqEhISODZZ58Nuv7y5ctZuXIl5557Lm3btmXp0qW0a9fuuDoDBw7k/vvvp2/fvhw+fJjMzExmzJhBSkoK3333HTt37qRly5Z06NAhz3PecccdPPPMMyxcuJDq1asD8MQTT1C1alUyMjLo0qUL33//PQ0bNqRfv35MmzaNli1bkpaWxumnn86oUaNISkpi3LhxIbtXwRLsEOlG4H6gpTrBtWjA31VSkciyZS6Rd44kPkbZpFWrVtSsWZNTTjmFhISEE1Ij7N+/n23bttG3b1/ALTI8/fTTWbJkCddeey1RUVH86U9/4tJLL2XZsmVBnTOLd955h2bNmtG0aVNWrlzJqlWrWLt2Leeccw4tW7YE4IwzzuDUMCebD/bqlwApqnrQE6pvhhsulS2ee84taHz66YLrGr5SmJ6GXwRqRueVGsGPc27evJkxY8awbNkyqlSpwtChQzl8ODLTMwXbg5kAHBKRJsB/4R4XT/LNqkjkt9/gnXdcOrRKlcJtjVECqFSpEjVr1mSWNxnzyJEjHDp0iPbt2zNt2jQyMjLYsWMHixcvplWrVgWeKyv9Q1paGhUqVKBy5cr8+uuvfPTRRwA0aNCA7du3Z/eG9u/fT3p6+gmpI4qTYB1MuhfI6QOMU9UXgLL1VzZ5snuCNHx4uC0xShCTJ09m7NixNG7cmDZt2vDLL7/Qt29fGjduTJMmTejcuTP//ve/Ofvss/M9z7Bhw+jevTudOnWiSZMmNG3alIYNGzJgwADatm0LuGx206ZNY8SIETRp0oSuXbty+PBhOnXqxKpVq0hISGDatGnF8bGzCSpdg4h8htNEugFoD/wGfKeq8fk2LEZ8TdegChddBJUrw5df+nMNo0DKUrqGSCbk6RqAfsAA3HyYX0SkFjC6yJaWFA4fhs6dbVGjYRSSYJcK/CIi7wH1vKKdQPGt+Q43p50GYXjEZxglnWCXCtwMTAey9BPOA2b5ZFNkcegQLF0KpUwQq6QSzJDe8I/C3v9gg7y349YDpXkXWQ+cVagrlVQ++ADatXNOxggrMTEx7Nq1y5xMmFBVdu3aRUxMTNBtgo3BHFHVP7LkEETkVMqKdOyUKXDeeeBF6o3wUbNmTbZu3Uqx5/4xsomJiaFmzZpB1w/WwXwmIg8Cp4lIV+A24IOTsK9ksXs3fPQR3HGHy1xnhJXo6Og8hc2MyCTYv5r7gB24JNy3AHOBh/NtURp47z0398WUGg3jpAhWOna1qr6iqn9T1au9/QKHSCLSXUTWisgGEbk/l+MdvNwy6SJydY5jQ0RkvbcNKdSnChWzZ0P9+tC0aVgubxglnWClY9eKSC1V/THYExdFVUBEqgKP4tQMFLeae7aq7gn2+iFh+nTYssXkSAzjJAk2BlMFWCki3wAHswo1f62ibFUBABHJUhXIdjCqmuody/kM+HJgvqru9o7PB7oDiRQn5cu7HoxhGCdFsA7mv0/i3LmpCgSbSCUoRQIRGQYMA6gVak3o2293+XbvuCO05zWMMkS+MRgRiRGRO3EpLRsCS1X1s6ytOAzMD1V9WVVbqGqLGjVqhO7E+/bBK6/ATz8VXNcwjDwpKMj7Ji4OsgLowfHyJQVRFFWBIikSFJkPP3RPj7xEQYZhnBwFDZEaZa2YFpHXgG8Kce5sVQGcc+iPWzAZDJ8AT4pIFe99N+CBQly7aMyYAWefDa1bF9slDaM0UlAP5mjWjjqVgKDRIqgKeMHdx3FOahkwKivg6zt//AGffAK9e9vkOsMoIgX1YJqISJq3L7iZvGnevqrqGfk1VtW5uEl5gWWPBOwvww1/cms7EZhYgH2hZ88e6N4drryy2C9tGKWNfB2MqkYVlyERw5/+BO++G24rDKNUYGOAnGzdGm4LDKPUYA4mkNRUOP98eO21cFtiGKUCczCBeNnZySGeZRjGyWEOJpCPP4bYWFseYBghwhxMFunpsGgRdOtmixsNI0SYg8kiKQnS0px6gGEYIcEcTBb168OkSXDZZeG2xDBKDeFVxo4kqlaF664LtxWGUaqwHgw4YbVx42Bb8a2nNIxIIyMDvv4annvOZYmdMqXo57QeDMAXX8CIEe4J0nknpJ0xjFLLhx+61549XQKB9u3d63nnhWatrzkYgAULICrK3V3DKIWounmkixfD9u1wv5ch+1//gpgY52BiYmDuXGjYEAqhTJIv5mDA3fVmzeCMfNduGkaJQRU2b4aFC93si0WLjq2COfdcuOce9z916lS3/C6LUD/j8DUGE4SqQHkRmeYd/1pE6njldUTkdxFJ8bYXfTPyyBH45hvrvRilgiVL4IYb3Gj/z3+Gm26CefPgkktcmHHFCpeoMcpbxly7tuu5+IVvPZggVQVuBPao6gUi0h/4X6Cfd2yjqib4ZV82q1a5HDC2PMAogaSmOsfxj3+4ZXSrVsGsWdCpk+uldOoEF14Yvrmjfg6RClQV8N6P9PanA+NEivlWNG3qFBzLly/WyxrGybBxo4uTNG4Ml14Khw7B889Dx47OwQwZAjfeeKyHEm78dDDBqApk11HVdBHZB1TzjsWKyHIgDXhYVT/PeYGQqQqceebJtzUMH0lPd8OeOXPctnatK/+v/3IO5sIL3f/HChVceaT9n4zUIO92oJaq7hKR5sAsEblIVdMCK6nqy8DLAC1atChQafIEMjPhmmtg6FDo1SsEZhtG0dm/32VtnTXL9Vb27IFy5Vwv5bbb4C9/gQsucHVFjjmXSMRPBxOMMkBWna0icipQGdjlydIeAVDVZBHZCNQHkkJq4dq1Tn+6Z8+QntYwCsuhQ3D66W7/8svhyy+hWjWXGrp3b+jaFSpVCq+NJ4OfDiYYVYHZwBDgS+BqYIGqqojUAHZ7srV1gXrAppBbuGSJe7UArxFGxo2Dhx5yE8krVoRRo1yPpU0bODVSxxhB4ttj6mBUBYDXgGoisgH4J5D1KLsD8L2IpOCCv8N9URX4+mv3byKrv2kYPrN/P0ye7IY5X37pylq2hFtvdTMmwM1F6dCh5DsX8DkGE4SqwGGcZEnOdu8B7/lpG+Dmv7RqZflfDF85etTlMnv7bZg9G37/HWrVgh073PGLL3ZbaaQU+MiTJCMDatRwoXjDCDGqkJICb7wBiYnOmVSvDtdf7xYSXnJJ2ZDdKrsOJioKPv003FYYpZQePdyToHLloE8fGDzYBW+jo8NtWfFSdh2Mqg2NjJCRnAyvvuoCtlFRzqn07g3XXgtVqhTcvrRSBjppeXDddSZubxSJXbtg7163v2kTTJsG69a597fe6uaslGXnAmXZwSxZEnnTHo0SQXKyi6XUrOmm6QNccQX8/LObWWsco2w6mF9/hS1b3BMkwwiC9HR45x1o2xZatHDqwkOHHpMwj472d1VySaVsxmCWLXOvLVuG1w4j4klLg1degbFj4ccfoW5d+M9/XA+mcuVwWxf5lE0Hs3y5C/AmJITbEiOCmTTJZVJNS3OzGZ5/3q0qiZSVyiWBsulg4uLcL6ckLu4wfGX9ercm6Lzz3ATvHj3g7rvdsMgoPOLWFZZ8WrRooUlJoV0LaZQt9u1z6SSHDIHx48NtTclCRJJV9QQ3XPaCvEeOwC+/hNsKI0JYsQIee8ztV64Mb70Fjz4aXptKE2XPwXz1FZxzDsyfH25LjDCyejX06+cywz399LGE2H37Hp8E2ygaZc/BpKS41/j4sJphhIfUVDcEiotzyZwefNCVhUqmwzieiFQV8I494JWvFZHLQ2bU8uVw9tluM8oMO3fCXXdBgwZuPstdd7nZt0884VSDDX/wzcEEqAr0ABoB14pIoxzVslUFgP/gVAXw6vUHLgK6A+O98xWdlBR7PF3GGDvWSXiMHesWHa5fD2PGuMX0hr/42YPJVhVQ1T+ALFWBQPoAb3r704EunqpAH2Cqqh5R1c3ABu98RePIEVi50ikJGGWGrHksP/zgJs3ZcKj4iFRVgfOAr3K0PUE0utCqAqrw2msWfyljPPhg2ci9EomU6Il2hVYViIlxfWSjTGHOJXz4eesLoypAoKpAkG0Nw4hw/HQw2aoCIlIOF7SdnaNOlqoABKgKeOX9vadMsThVgW98tNUwDB/wbYjkxVSyVAWigIlZqgJAkqrOxqkKTPZUBXbjnBBevXdwMrPpwO2qmuGXrYZh+IOtRTIMo8jYWiTDMIodczCGYfhGqRkiicgOYEuIT1sd2Bnic4aCSLQrEm2CyLQrEm2CotlVW1VPmBtdahyMH4hIUm7jynATiXZFok0QmXZFok3gj102RDIMwzfMwRiG4RvmYPLn5XAbkAeRaFck2gSRaVck2gQ+2GUxGMMwfMN6MIZh+EaZdTBBZNv7j4ikeNs6EdkbcCwj4FjO9VVFsWmiiPwmIj/kcVxEZKxn8/ci0izg2BARWe9tQ3Jr75NNAz1bVojIFyLSJOBYqleeIiIhnWYdhF0dRWRfwPf0SMCxfL97H226J8CeH7zfUVXvmC/3SkTOF5GFIrJKRFaKyD9yqePf70pVy9yGWxu1EagLlAO+AxrlU38Ebi1V1vsDPtnVAWgG/JDH8b8AHwECtAa+9sqrApu81yrefpVisqlN1rVw2Qu/DjiWClQP073qCMwp6ncfSpty1P0rbnGvr/cKOAdo5u1XAtbl/Lx+/q7Kag8mmGx7gVwLJPptlKouxi36zIs+wCR1fAWcKSLnAJcD81V1t6ruAebjUo36bpOqfuFdE1ySsGLJFxfEvcqLwn73ftlUXL+p7ar6rbe/H1jNicnbfPtdlVUHk1u2vRMy5gGISG0gFlgQUBwjIkki8pWIXOGblSeSl91Bfx6fuRH3nzALBeaJSLKXfbC4uUREvhORj0TkIq8s7PdKRE7H/aG+F1Ds+70Sl1S/KfB1jkO+/a5KdEa7YqI/MF2PTxdRW1W3iUhdYIGIrFDVjWGyLyIQkU44B9MuoLidd5/OAuaLyBrvv3xx8C3uezogIn8BZuHyCkUCfwWWqmpgb8fXeyUiFXEO7U5VTQvVeQuirPZgCpMxrz85urKqus173QQswv1XKA7ysjusGQBFpDHwKtBHVXdllQfcp9+AmYQicXuQqGqaqh7w9ucC0SJSncjIlpjfbyrk90pEonHO5W1VnZFLFf9+V6EOKpWEDddz24Qb+mQF+i7KpV5DXPBNAsqqAOW9/erAekIUJPTOWYe8A5c9OT4Y940eC8Zt9myr4u1XLSabauFUH9rkKK8AVArY/wLoHuLvMT+7zubYPK9WwI/efQvqu/fDJu94ZVycpkJx3CvvM08Cns2njm+/qzI5RNLgsu2B+08zVb277XEh8JKIZOJ6gE+p6qpQ2CUiibinH9VFZCvwKBDt2fwiMBcX8d8AHAKu947tFpHHcWlKAUbp8d1vP216BKcEMV5EANLVLZj7EzDTKzsVmKKqH4fCpiDtuhq4VUTSgd+B/t73mOt3X0w2AfQF5qnqwYCmft6rtsB1wAoRSfHKHsT9Y/D9d2UzeQ3D8I2yGoMxDKMYMAdjGIZvmIMxDMM3zMEYhuEb5mAMw/ANczBGSAhYYf6DiHwgImcWsv0iEWnh7c8tbHsjMjEHY4SK31U1QVXjcBPJbj/ZE6nqX1R1b8gsM8KGORjDD77EWxQnIq1E5EsRWe7li2nglZ8mIlNFZLWIzAROy2rs5UapLiJ1AnOriMjdIjLS27/Dy3HyvYhMLdZPZwRNmZzJa/iHiEQBXXC64wBrgPbe7OnLgCeBq4BbgUOqeqG3lunbQl7qfiBWVY/YcCpyMQdjhIrTvKno5+Fyjsz3yisDb4pIPVxKgmivvAMwFkBVvxeR7wt5ve+Bt0VkFm6ltBGB2BDJCBW/q2oCUBu3aC4rBvM4sNCLzfwViCnEOdM5/jca2LYn8AIug9wyEbF/lhGIORgjpKjqIeAO4L+8P/rKHFviPzSg6mJgAICIxAGNczndr8BZIlJNRMoDvbz6pwDnq+pC4D7vGhVD/2mMomIOxgg5qrocN4S5Fvg38D8ispzjh+QTgIoishoYBSTncp6j3rFvcEOuNd6hKOAtEVkBLAfG2lOnyMRWUxuG4RvWgzEMwzfMwRiG4RvmYAzD8A1zMIZh+IY5GMMwfMMcjGEYvmEOxjAM3zAHYxiGb/w/CRzA0Db6SzYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_all = np.array(p_all)\n",
    "p_art = np.array(p_art)\n",
    "r_all = np.array(r_all)\n",
    "contact = np.array(contact)\n",
    "fig,ax = plt.subplots(1,1,figsize=(4,3))\n",
    "\n",
    "ax.plot(r_all[~contact],p_all[~contact],'--r',label='Balloon only')\n",
    "ax.plot(r_all[contact],p_art,'-.b',label='Artery only')\n",
    "ax.plot(r_all[contact],p_all[contact],'-k',label='In contact')\n",
    "ax.set_xlabel('Radius')\n",
    "ax.set_ylabel('Pressure (kPa)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}