vms_radial_burn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1126920a",
   "metadata": {},
   "source": [
    "# Radial Burn\n",
    "In radial burn, combustion starts along the symmetry axis and proceeds radially outwards. Thus, at any instant, the unburned portion resembles a hollow pipe and both the axial and transverse radius of gyration of the system are variable as shown in Figure {numref}`fig-4-12`.\n",
    "\n",
    "```{figure} images/Fig11.png\n",
    ":name: fig-4-12\n",
    ":height: 300px\n",
    "\n",
    "Radially burning cylinder.\n",
    "```\n",
    "\n",
    "For the results presented here, the initial radius and length are $1$ m. The instantaneous moment of inertia scalars for the radially burning cylinder are\n",
    "\n",
    "```{math}\n",
    ":label: eq58\n",
    "I = m\\left(\\frac{R^2 + r^2}{4} + \\frac{h^2}{3}\\right)\n",
    "```\n",
    "\n",
    "and\n",
    "\n",
    "```{math}\n",
    ":label: eq59\n",
    "J = m\\frac{R^2 + r^2}{2}.\n",
    "```\n",
    "\n",
    "Their corresponding time derivatives are\n",
    "\n",
    "```{math}\n",
    "\\begin{aligned}\n",
    "\\dot I &= \\dot m\\left(\\frac{R^2}{2} + \\frac{h^2}{3}\\right) \\\\\n",
    "\\dot J &= \\dot m {r^2}.\n",
    "\\end{aligned}\n",
    "```\n",
    "\n",
    "The angular speeds evaluate to\n",
    "\n",
    "```{math}\n",
    ":label: eq60\n",
    "\\omega_3 = \\omega_{30} \\frac{R^4}{(R^2 + r^2)\\sqrt{(R^2 + r^2)(R^2 - r^2)}} \\\\\n",
    "\\omega_{12} = \\omega_0 \\left(\\frac{R^2 + \\frac{4h^2}{3}}{R^2 + \\frac{4h^2}{3} + r^2}\\right)^{\\frac{3R^2 + \\frac{16h^2}{3}}{2R^2 + \\frac{4h^2}{3}}}\n",
    "\\left(\\frac{R^2 - r^2}{R^2}\\right)^{\\frac{-R^2 + \\frac{8h^2}{3}}{2R^2 + \\frac{4h^2}{3}}}.\n",
    "```\n",
    "\n",
    "Evidently, the spin rate is not constant for a cylinder subject to radial burn. It varies such that it appears to damp out, slowly, for about two-thirds of the burn. Towards the end of the burn, i.e., as $r$ approaches $R$, the spin rate grows without bounds. Analytically, this growth in spin rate has been found to begin at the instant $r/R = 0.707$ {cite}`ekemao3`. The transverse angular speed, on the other hand, can either grow, decay, or fluctuate depending on the shape of the cylinder before the start of the burn. The corresponding body and space surfaces for the radially burning cylinder are shown in Figures {numref}`fig-4-15` and {numref}`fig-4-16`, respectively, for a burn duration of $90$ s; this duration is chosen purely to study a case where the end mass of the system approaches the mass of a payload. \n",
    "\n",
    "```{figure} images/Fig12.png\n",
    ":name: fig-4-15\n",
    ":height: 300px\n",
    "\n",
    "Radial burn: body surface.\n",
    "```\n",
    "\n",
    "```{figure} images/Fig13.png\n",
    ":name: fig-4-16\n",
    ":height: 300px\n",
    "\n",
    "Radial burn: space surface.\n",
    "```\n",
    "\n",
    "The case presented here is for $R/h=2$ and is a nutationally stable configuration; this is in agreement with Mao and Eke's work {cite}`ekemao3` where it was shown that transverse rates are unbounded for oblate radially burning cylinders with a $R/h \\geq \\sqrt{8/3}$. The cone angle does in fact monotonically decrease though this is not as evident from Figure {numref}`fig-4-15` when compared to the uniform and end burns; this is primarily due to the fact that the system's end mass is approximately 315 kg at the end of this radial burn whereas in the other burns it is nearly zero. One approach to verifying the nutational stability of the system is by plotting the time evolution of $\\theta$. However, the stability can also be verified from the angular speeds. For the extreme case of oblateness represented by a radially burning flat disk (i.e., $R \\gg h$), Equation {eq}`eq68` can be reformulated as\n",
    "\n",
    "```{math}\n",
    ":label: eq61\n",
    "\\omega_{12} = \\omega_0 \\frac{\\omega_3}{\\omega_{30}}\n",
    "```\n",
    "\n",
    "by dropping the terms involving $h$ from the solution to the transverse speed. Extending this rationale to the inertia scalars given by Equations {eq}`eq58` and {eq}`eq59`, we also get $I/J = 1/2$. Then, the nutation angle from Equation {eq}`eq39` can be expressed as\n",
    "\n",
    "```{math}\n",
    ":label: eq62\n",
    "\\theta = \\tan^{-1}\\left(\\frac{\\omega_0}{2\\omega_{30}}\\right) = \\text{constant}.\n",
    "```\n",
    "\n",
    "which is clearly bounded. In other words, even the most oblate configuration of a radially burning cylinder results in a constant nutation angle and thus does not show growth in the coning motion. However, unlike the preceding burns, the radial burn would require a nutation damper system to attenuate the cone angle and bring the system into a state of pure spin."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93fc41f9",
   "metadata": {},
   "source": [
    "# Simulation Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea332813",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "976e9ff3",
   "metadata": {},
   "source": [
    "## Define the global physical parameters for the variable mass system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb393ee5",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_ini = 1020\n",
    "rho_exhaust = 0.002 * rho_ini\n",
    "L = 1  # m\n",
    "h = L / 2  # m\n",
    "R = L  # m\n",
    "u = 5  # m/s\n",
    "m_dot = -rho_exhaust * u * np.pi * R**2\n",
    "w0 = 0.2  # rad/s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e93a0c34",
   "metadata": {},
   "source": [
    "## Initial angular rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62120119",
   "metadata": {},
   "outputs": [],
   "source": [
    "w10 = 0  # rad/s\n",
    "w20 = w0  # rad/s\n",
    "w30 = 0.3  # rad/s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d971f79d",
   "metadata": {},
   "source": [
    "## Initial mass properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46f87afe",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf0 = np.pi * R**2 * L * rho_ini  # kg\n",
    "I10 = mf0 * (R**2 / 4 + h**2 / 3)  # kg*m^2\n",
    "I30 = mf0 * R**2 / 2  # kg*m^2\n",
    "chi0 = 0  # Initial precession angle\n",
    "Fd0 = 0  # Initial F value\n",
    "tb = -mf0 / m_dot - 1  # Burn time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31d30c6b",
   "metadata": {},
   "source": [
    "## Quaternion initial conditions (representing no initial rotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b50c2434",
   "metadata": {},
   "outputs": [],
   "source": [
    "q0 = [1, 0, 0, 0]  # (w, x, y, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62839651",
   "metadata": {},
   "source": [
    "## Combine initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc0a3d16",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y0 = [mf0, I10, I30, w10, w20, w30, chi0, Fd0] + q0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "671b0911",
   "metadata": {},
   "source": [
    "## Define the differential equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4384a727",
   "metadata": {},
   "outputs": [],
   "source": [
    "def uniburn(t, w):\n",
    "    m, I1, I3, w1, w2, w3, chi, F, qw, qx, qy, qz = w\n",
    "    # Mass varying terms\n",
    "    md = -rho_exhaust * np.pi * R**2 * u  # mass\n",
    "    I1d = md * (R**2 / 4 + h**2 / 3)  # transverse moment of inertia\n",
    "    I3d = md * R**2 / 2  # spin moment of inertia\n",
    "    # Equations of motion for the axisymmetric cylinder undergoing uniform burn\n",
    "    w1d = (I1 - I3) * w2 * w3 / I1 - (I1d - md * (h**2 + R**2 / 4)) * w1 / I1\n",
    "    w2d = -(I1 - I3) * w1 * w3 / I1 - (I1d - md * (h**2 + R**2 / 4)) * w2 / I1\n",
    "    w3d = -(I3d - md * R**2 / 2) * w3\n",
    "    chid = (1 - I3 / I1) * w3\n",
    "    Fd = -(I1d - md * (h**2 + R**2 / 4)) / I1\n",
    "    # Quaternion derivative\n",
    "    omega_quat = np.array([0, w1, w2, w3])\n",
    "    quat = np.array([qw, qx, qy, qz])\n",
    "    quat_dot = 0.5 * np.array(quat_mult(quat, omega_quat))\n",
    "    wd = [md, I1d, I3d, w1d, w2d, w3d, chid, Fd, quat_dot[0], quat_dot[1], quat_dot[2], quat_dot[3]]\n",
    "    return wd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a87a3f8b",
   "metadata": {},
   "source": [
    "## Function to multiply two quaternions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72da90de",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quat_mult(q, r):\n",
    "    w1, x1, y1, z1 = q\n",
    "    w2, x2, y2, z2 = r\n",
    "    return [\n",
    "        w1*w2 - x1*x2 - y1*y2 - z1*z2,\n",
    "        w1*x2 + x1*w2 + y1*z2 - z1*y2,\n",
    "        w1*y2 - x1*z2 + y1*w2 + z1*x2,\n",
    "        w1*z2 + x1*y2 - y1*x2 + z1*w2\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e87f0eae",
   "metadata": {},
   "source": [
    "## Time span for the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72063701",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 0.01\n",
    "t_eval = np.arange(0, tb, dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95f428ed",
   "metadata": {},
   "source": [
    "## Solve the ODEs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "933acba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(uniburn, [0, tb], Y0, t_eval=t_eval, atol=1e-9, rtol=1e-8)\n",
    "\n",
    "# Extract results if the integration was successful\n",
    "if sol.status == 0:\n",
    "    print(\"Integration successful.\")\n",
    "    t = sol.t\n",
    "    m = sol.y[0]\n",
    "    I1 = sol.y[1]\n",
    "    I3 = sol.y[2]\n",
    "    omega1 = sol.y[3]\n",
    "    omega2 = sol.y[4]\n",
    "    omega3 = sol.y[5]\n",
    "    chi = sol.y[6]\n",
    "    F = sol.y[7]\n",
    "    qw, qx, qy, qz = sol.y[8], sol.y[9], sol.y[10], sol.y[11]\n",
    "else:\n",
    "    print(f\"Integration failed with status {sol.status}: {sol.message}\")\n",
    "\n",
    "# Print the final time to check if it reached the end of the burn time\n",
    "print(f\"Final time: {sol.t[-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "663e3936",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the quaternion norm over time\n",
    "quat_norm = np.sqrt(qw**2 + qx**2 + qy**2 + qz**2)\n",
    "norm_violation_indices = np.where(np.abs(quat_norm - 1) > 1e-6)[0]\n",
    "if len(norm_violation_indices) > 0:\n",
    "    print(f\"Quaternion constraint violated at time steps: {t[norm_violation_indices]}\")\n",
    "else:\n",
    "    print(\"Quaternion constraint satisfied throughout the simulation.\")\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(t, quat_norm, label='Quaternion Norm')\n",
    "plt.axhline(y=1.0, color='r', linestyle='--', label='Ideal Norm (1.0)')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Quaternion Norm')\n",
    "plt.title('Quaternion Norm vs Time')\n",
    "plt.ylim([0.99, 1.01])  # Adjust y-limits for better visualization\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "685c6503",
   "metadata": {},
   "source": [
    "## Function to convert quaternion to Euler angles (Z-X-Z sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00642cab",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quat_to_euler_zxz(q):\n",
    "    \"\"\"\n",
    "    Convert a quaternion to Euler angles (Z-X-Z sequence).\n",
    "    \"\"\"\n",
    "    w, x, y, z = q\n",
    "    # Compute the Euler angles\n",
    "    psi = np.arctan2(2*(w*z + x*y), 1 - 2*(y**2 + z**2))\n",
    "    theta = np.arccos(2*(w*y - z*x))\n",
    "    phi = np.arctan2(2*(w*z + y*x), 1 - 2*(x**2 + y**2))\n",
    "    return psi, theta, phi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a03c8af5",
   "metadata": {},
   "source": [
    "## Extract Euler angles from quaternions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "241ac442",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi, theta, phi = np.zeros(len(qw)), np.zeros(len(qw)), np.zeros(len(qw))\n",
    "for i in range(len(qw)):\n",
    "    psi[i], theta[i], phi[i] = quat_to_euler_zxz([qw[i], qx[i], qy[i], qz[i]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "522a8d4e",
   "metadata": {},
   "source": [
    "## Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc2f3b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "\n",
    "# Plot angular velocities\n",
    "plt.subplot(4, 1, 1)\n",
    "plt.plot(sol.t, omega1, label='omega1')\n",
    "plt.plot(sol.t, omega2, label='omega2')\n",
    "plt.plot(sol.t, omega3, label='omega3')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Angular Velocity (rad/s)')\n",
    "plt.title('Angular Velocity vs Time')\n",
    "plt.legend()\n",
    "\n",
    "# Plot Euler angles\n",
    "plt.subplot(4, 1, 2)\n",
    "plt.plot(sol.t, psi, label='psi')\n",
    "plt.plot(sol.t, theta, label='theta')\n",
    "plt.plot(sol.t, phi, label='phi')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Euler Angles (rad)')\n",
    "plt.title('Euler Angles vs Time')\n",
    "plt.legend()\n",
    "\n",
    "# Plot quaternions\n",
    "plt.subplot(4, 1, 3)\n",
    "plt.plot(sol.t, qw, label='qw')\n",
    "plt.plot(sol.t, qx, label='qx')\n",
    "plt.plot(sol.t, qy, label='qy')\n",
    "plt.plot(sol.t, qz, label='qz')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Quaternion Components')\n",
    "plt.title('Quaternion Components vs Time')\n",
    "plt.legend()\n",
    "\n",
    "# Plot mass and moments of inertia\n",
    "plt.subplot(4, 1, 4)\n",
    "plt.plot(sol.t, m, label='mass')\n",
    "plt.plot(sol.t, I1, label='I1')\n",
    "plt.plot(sol.t, I3, label='I3')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Mass and Inertias')\n",
    "plt.title('Mass and Moments of Inertia vs Time')\n",
    "plt.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec44129f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the T-handle's total mechanical energy over time\n",
    "E = 0.5 * (I1 * omega1**2 + I1 * omega2**2 + I3 * omega3**2)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, E, '-b', linewidth=2)\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Total mechanical energy (J)')\n",
    "plt.ylim([min(E) * 0.8, max(E) * 1.2])  # Set fixed y-axis limits\n",
    "plt.title('Total Mechanical Energy vs Time')\n",
    "plt.show()\n",
    "\n",
    "# Plot the components of the angular momentum about the mass center and the total angular momentum over time\n",
    "H1 = I1 * omega1  # kg-m^2/s\n",
    "H2 = I1 * omega2  # kg-m^2/s\n",
    "H3 = I3 * omega3  # kg-m^2/s\n",
    "H = np.sqrt(H1**2 + H2**2 + H3**2)  # kg-m^2/s\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, H1, label='H \\cdot e1')\n",
    "plt.plot(sol.t, H2, label='H \\cdot e2')\n",
    "plt.plot(sol.t, H3, label='H \\cdot e3')\n",
    "plt.plot(sol.t, H, label='||H||')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Angular momentum (kg-m^2/s)')\n",
    "plt.title('Angular Momentum Components vs Time')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2eee70c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to convert quaternion to rotation matrix\n",
    "def quat_to_rot_matrix(q):\n",
    "    \"\"\"\n",
    "    Convert a quaternion q to a rotation matrix.\n",
    "    \"\"\"\n",
    "    w, x, y, z = q\n",
    "    return np.array([\n",
    "        [1 - 2*(y**2 + z**2), 2*(x*y - z*w), 2*(x*z + y*w)],\n",
    "        [2*(x*y + z*w), 1 - 2*(x**2 + z**2), 2*(y*z - x*w)],\n",
    "        [2*(x*z - y*w), 2*(y*z + x*w), 1 - 2*(x**2 + y**2)]\n",
    "    ])\n",
    "\n",
    "def animate_t_handle_quat(qw, qx, qy, qz, dt):\n",
    "    # Specify dimensions for the T-handle\n",
    "    LAG = 0.5  # cm\n",
    "    LBC = 4  # cm\n",
    "    LAD = 2  # cm\n",
    "\n",
    "    # Initialize arrays to store the T-handle's orientation and key points\n",
    "    e1 = np.zeros((3, len(qw)))\n",
    "    e2 = np.zeros((3, len(qw)))\n",
    "    e3 = np.zeros((3, len(qw)))\n",
    "    xA, yA, zA = np.zeros(len(qw)), np.zeros(len(qw)), np.zeros(len(qw))\n",
    "    xB, yB, zB = np.zeros(len(qw)), np.zeros(len(qw)), np.zeros(len(qw))\n",
    "    xC, yC, zC = np.zeros(len(qw)), np.zeros(len(qw)), np.zeros(len(qw))\n",
    "    xD, yD, zD = np.zeros(len(qw)), np.zeros(len(qw)), np.zeros(len(qw))\n",
    "\n",
    "    # Calculate the orientation of the T-handle over time\n",
    "    for k in range(len(qw)):\n",
    "        q = [qw[k], qx[k], qy[k], qz[k]]\n",
    "        R = quat_to_rot_matrix(q)\n",
    "        e1[:, k] = R @ np.array([1, 0, 0])\n",
    "        e2[:, k] = R @ np.array([0, 1, 0])\n",
    "        e3[:, k] = R @ np.array([0, 0, 1])\n",
    "        xA[k] = -LAG * e2[0, k]\n",
    "        yA[k] = -LAG * e2[1, k]\n",
    "        zA[k] = -LAG * e2[2, k]\n",
    "        xB[k] = xA[k] + LBC / 2 * e1[0, k]\n",
    "        yB[k] = yA[k] + LBC / 2 * e1[1, k]\n",
    "        zB[k] = zA[k] + LBC / 2 * e1[2, k]\n",
    "        xC[k] = xA[k] - LBC / 2 * e1[0, k]\n",
    "        yC[k] = yA[k] - LBC / 2 * e1[1, k]\n",
    "        zC[k] = zA[k] - LBC / 2 * e1[2, k]\n",
    "        xD[k] = xA[k] + LAD * e2[0, k]\n",
    "        yD[k] = yA[k] + LAD * e2[1, k]\n",
    "        zD[k] = zA[k] + LAD * e2[2, k]\n",
    "\n",
    "    # Set up the figure window\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    ax.set_xlabel('X (cm)')\n",
    "    ax.set_ylabel('Y (cm)')\n",
    "    ax.set_zlabel('Z (cm)')\n",
    "    ax.set_xlim([-LBC, LBC])\n",
    "    ax.set_ylim([-LBC, LBC])\n",
    "    ax.set_zlim([-LAD, LAD])\n",
    "    ax.set_title('T-handle Animation')\n",
    "\n",
    "    # Draw the T-handle\n",
    "    AD, = ax.plot([xA[0], xD[0]], [yA[0], yD[0]], [zA[0], zD[0]], 'k-', linewidth=5)\n",
    "    BC, = ax.plot([xB[0], xC[0]], [yB[0], yC[0]], [zB[0], zC[0]], 'k-', linewidth=5)\n",
    "\n",
    "    # Animate the T-handle's motion by updating the figure with its current orientation\n",
    "    def update(k):\n",
    "        AD.set_data([xA[k], xD[k]], [yA[k], yD[k]])\n",
    "        AD.set_3d_properties([zA[k], zD[k]])\n",
    "        BC.set_data([xB[k], xC[k]], [yB[k], yC[k]])\n",
    "        BC.set_3d_properties([zB[k], zC[k]])\n",
    "        return AD, BC,\n",
    "\n",
    "    ani = FuncAnimation(fig, update, frames=len(qw), interval=dt * 1000, blit=True)\n",
    "    plt.show()\n",
    "\n",
    "# Example usage\n",
    "# Assuming `qw`, `qx`, `qy`, `qz` are the quaternion components obtained from the previous solution\n",
    "animate_t_handle_quat(qw, qx, qy, qz, dt)"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}