{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#1a\n",
        "# Explanation of the Physics in the Code\n",
        "\n",
        "# 1️Newton's Law of Universal Gravitation\n",
        "# -----------------------------------------\n",
        "# The force between two masses due to gravity is given by:\n",
        "\n",
        "# F = (G * M * m) / r**2\n",
        "\n",
        "# Where:\n",
        "# F = gravitational force (N)\n",
        "# G = gravitational constant (6.67430 × 10⁻¹¹ m³/kg/s²)\n",
        "# M = mass of Earth (kg)\n",
        "# m = mass of the satellite (kg)\n",
        "# r = distance between Earth and the satellite (m)\n",
        "\n",
        "# This force acts toward Earth's center, keeping the satellite in orbit.\n",
        "\n",
        "# Centripetal Force & Orbital Velocity\n",
        "# ----------------------------------------\n",
        "# A satellite remains in orbit because its velocity is just right to counteract gravity:\n",
        "\n",
        "# v = sqrt(G * M / r)\n",
        "\n",
        "# Where:\n",
        "# v = orbital velocity (m/s)\n",
        "# r = distance from the center of Earth (m)\n",
        "\n",
        "# This equation ensures that the satellite is in continuous free-fall,\n",
        "# meaning it \"falls\" around the Earth without ever hitting it.\n",
        "\n",
        "#  Kepler's Laws of Motion\n",
        "# ---------------------------\n",
        "# First Law: The satellite follows an elliptical orbit with Earth at one focus.\n",
        "# Second Law: The satellite moves faster when closer to Earth and slower when farther.\n",
        "# Third Law: The time to complete one orbit (T) depends on its distance:\n",
        "\n",
        "# T^2 = (4 * pi**2 * r^3) / (G * M)\n",
        "\n",
        "# Where:\n",
        "# T = orbital period (seconds)\n",
        "# r = average distance from Earth (m)\n",
        "\n",
        "# What You’ll See When You Run the Code:\n",
        "# --------------------------------------\n",
        "# A realistic orbital path around Earth.\n",
        "# If you lower the altitude, the satellite will orbit faster.\n",
        "# If you increase the velocity, the satellite's orbit will stretch into an ellipse.\n",
        "\n",
        "# Try modifying the altitude or velocity in the simulation to see different orbits!\n"
      ],
      "metadata": {
        "id": "5hxMXxnJBfzl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#1b\n",
        "# Constants\n",
        "altitude = 500e3  # Orbital altitude in meters (500 km above Earth)\n",
        "# Try different altitude values:\n",
        "# altitude = 2000e3   # 2,000 km (Higher Low Earth Orbit)\n",
        "# altitude = 20000e3  # 20,000 km (Medium Earth Orbit)\n",
        "# altitude = 35786e3  # 35,786 km (Geostationary Orbit)\n",
        "\n",
        "# Initial Conditions for the Satellite\n",
        "r0 = R_earth + altitude  # Initial distance from Earth's center (m)\n",
        "\n",
        "# Compute the velocity needed for a circular orbit\n",
        "v0 = np.sqrt(G * M_earth / r0)  # Circular orbital speed (m/s)\n",
        "\n",
        "# Modify velocity for different orbit types:\n",
        "# v0 = np.sqrt(G * M_earth / r0) * 1.1  # 10% faster → Elliptical orbit\n",
        "# v0 = np.sqrt(G * M_earth / r0) * 0.9  # 10% slower → Elliptical orbit (closer to Earth)\n",
        "# v0 = np.sqrt(G * M_earth / r0) * 1.414  # Escape velocity → Satellite leaves Earth's orbit\n",
        "# v0 = np.sqrt(G * M_earth / r0) * 0.8  # Too slow → Satellite will crash back to Earth\n",
        "\n",
        "# Set initial conditions [x, y, vx, vy]\n",
        "initial_state = [r0, 0, 0, v0]  # Starts at (r0,0) moving in the y-direction\n",
        "\n"
      ],
      "metadata": {
        "id": "dGje5TqwDXWX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#2\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import solve_ivp\n",
        "from matplotlib.animation import FuncAnimation\n",
        "from IPython.display import HTML\n",
        "\n",
        "# Constants\n",
        "G = 6.67430e-11  # Universal Gravitational Constant (m³/kg/s²)\n",
        "M_earth = 5.972e24  # Mass of Earth (kg)\n",
        "R_earth = 6.378e6  # Radius of Earth (m)\n",
        "altitude = 500e3  # Orbital altitude in meters (500 km above Earth)\n",
        "\n",
        "# Initial Conditions for the Satellite\n",
        "r0 = R_earth + altitude  # Initial distance from Earth's center (m)\n",
        "v0 = np.sqrt(G * M_earth / r0) # Circular orbital speed (m/s)\n",
        "initial_state = [r0, 0, 0, v0]  # [x, y, vx, vy] -> Start at (r0,0) with velocity along y-axis\n",
        "\n",
        "# Simulation Time\n",
        "t_max = 6000  # Simulate for ~1 full orbit\n",
        "t_eval = np.linspace(0, t_max, num=500)  # Fewer points for smoother animation\n",
        "\n",
        "# Equations of Motion: Newton's Laws & Gravity\n",
        "def orbital_dynamics(t, state):\n",
        "    \"\"\"\n",
        "    Computes acceleration due to gravity at a given position.\n",
        "    \"\"\"\n",
        "    x, y, vx, vy = state  # Unpack position and velocity\n",
        "    r = np.sqrt(x**2 + y**2)  # Compute radial distance from Earth's center\n",
        "\n",
        "    # Compute acceleration due to gravity: a = -GM/r^2 (directed towards center)\n",
        "    ax = -G * M_earth * x / r**3  # Acceleration in x-direction\n",
        "    ay = -G * M_earth * y / r**3  # Acceleration in y-direction\n",
        "\n",
        "    return [vx, vy, ax, ay]  # Return derivatives for solver\n",
        "\n",
        "# Solve the Equations of Motion\n",
        "solution = solve_ivp(orbital_dynamics, [0, t_max], initial_state, t_eval=t_eval, method='RK45')\n",
        "\n",
        "# Extract the position values\n",
        "x_vals, y_vals = solution.y[0], solution.y[1]\n",
        "\n",
        "# Setup Figure for Animation\n",
        "fig, ax = plt.subplots(figsize=(6,6))\n",
        "ax.set_xlim(-r0 / 1e6, r0 / 1e6)  # Convert to km\n",
        "ax.set_ylim(-r0 / 1e6, r0 / 1e6)\n",
        "\n",
        "# Plot Earth\n",
        "earth = plt.Circle((0, 0), R_earth / 1e6, color='blue', alpha=0.5)  # Convert to km\n",
        "ax.add_patch(earth)\n",
        "\n",
        "# Plot the orbit path\n",
        "ax.plot(x_vals / 1e6, y_vals / 1e6, 'gray', linestyle='--', alpha=0.5, label=\"Orbit Path\")\n",
        "\n",
        "# Satellite point\n",
        "satellite, = ax.plot([], [], 'ro', markersize=5, label=\"Satellite\")\n",
        "\n",
        "ax.set_xlabel(\"X Position (km)\")\n",
        "ax.set_ylabel(\"Y Position (km)\")\n",
        "ax.set_title(\"Animated Satellite Orbit\")\n",
        "ax.legend()\n",
        "ax.set_aspect('equal')\n",
        "\n",
        "# Update function for animation\n",
        "def update(frame):\n",
        "    # Wrap single values in a list to ensure FuncAnimation works\n",
        "    satellite.set_data([x_vals[frame] / 1e6], [y_vals[frame] / 1e6])\n",
        "    return satellite,\n",
        "\n",
        "plt.rcParams['animation.embed_limit'] = 50  # Increase to 50MB or more if needed\n",
        "# Run animation\n",
        "ani = FuncAnimation(fig, update, frames=len(t_eval), interval=20, blit=True)\n",
        "\n",
        "# Display animation\n",
        "HTML(ani.to_jshtml())\n"
      ],
      "metadata": {
        "id": "ahj2etxCBnzJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "upce2OfEQ5br"
      },
      "outputs": [],
      "source": [
        "#3a\n",
        "#satellite tle data\n",
        "#https://celestrak.org/NORAD/elements/satnogs.txt\n",
        "#http://mstl.atl.calpoly.edu/~ops/keps/kepler.txt\n",
        "#https://www.amsat.org/amsat/ftp/keps/current/nasabare.txt\n",
        "#https://db.satnogs.org/satellites/\n",
        "#https://celestrak.org/NORAD/elements/gp.php?GROUP=active&FORMAT=tle\n",
        "#space debris tle data\n",
        "#https://celestrak.org/NORAD/elements/gp.php?GROUP=cosmos-1408-debris&FORMAT=tle\n",
        "#https://celestrak.org/NORAD/elements/gp.php?GROUP=fengyun-1c-debris&FORMAT=tle\n",
        "#https://celestrak.org/NORAD/elements/gp.php?GROUP=iridium-33-debris&FORMAT=tle\n",
        "#https://celestrak.org/NORAD/elements/gp.php?GROUP=cosmos-2251-debris&FORMAT=tle\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#3b\n",
        "!pip install skyfield\n",
        "#The command pip install skyfield installs the Skyfield library, which is a powerful Python package used for astronomy and satellite orbit calculations. It's the tool you need to work with TLE data and simulate the motion of real satellites."
      ],
      "metadata": {
        "id": "fsYmIfash1td"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#3c\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import datetime as dt\n",
        "import math as m\n",
        "from matplotlib.animation import FuncAnimation\n",
        "from IPython.display import HTML\n",
        "#matplotlib.pyplot for plotting the orbits\n",
        "#numpy for math operations and array handling\n",
        "#datetime for handling TLE timestamps\n",
        "#math for trigonometric and orbital calculations\n",
        "\n",
        "mu = 398600.4418\n",
        "r = 6781\n",
        "D = 24 * 0.997269\n",
        "\n",
        "def plot_tle(data):\n",
        "    fig = plt.figure()\n",
        "    ax = plt.axes(projection='3d', computed_zorder=False)\n",
        "\n",
        "    # Plot Earth\n",
        "    u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]\n",
        "    ax.plot_wireframe(r * np.cos(u) * np.sin(v),\n",
        "                     r * np.sin(u) * np.sin(v),\n",
        "                     r * np.cos(v), color=\"b\", alpha=0.5, lw=0.5, zorder=0)\n",
        "\n",
        "    # Set plot labels and title\n",
        "    ax.set_xlabel(\"X-axis (km)\")\n",
        "    ax.set_ylabel(\"Y-axis (km)\")\n",
        "    ax.set_zlabel(\"Z-axis (km)\")\n",
        "    ax.xaxis.set_tick_params(labelsize=7)\n",
        "    ax.yaxis.set_tick_params(labelsize=7)\n",
        "    ax.zaxis.set_tick_params(labelsize=7)\n",
        "    ax.set_aspect('equal', adjustable='box')\n",
        "\n",
        "    lines = []\n",
        "    for sat_name, tle in data.items():\n",
        "        line, = ax.plot([], [], [], label=sat_name)\n",
        "        lines.append(line)\n",
        "\n",
        "    def init():\n",
        "        for line in lines:\n",
        "            line.set_data([], [])\n",
        "            line.set_3d_properties([])\n",
        "        return lines\n",
        "\n",
        "    def update(frame):\n",
        "        for line, (sat_name, tle) in zip(lines, data.items()):\n",
        "            tle1, tle2 = tle\n",
        "            if tle1[0] != \"1\":\n",
        "                continue\n",
        "\n",
        "            year_str = tle1[18:20]\n",
        "            if int(year_str) > int(dt.date.today().year % 100):\n",
        "                year_prefix = \"19\"\n",
        "            else:\n",
        "                year_prefix = \"20\"\n",
        "\n",
        "            orb = {\"t\": dt.datetime.strptime(\n",
        "                year_prefix + year_str + \" \" + tle1[20:23] + \" \" +\n",
        "                str(int(24 * float(tle1[23:33]) // 1)) + \" \" +\n",
        "                str(int(((24 * float(tle1[23:33]) % 1) * 60) // 1)) + \" \" +\n",
        "                str(int((((24 * float(tle1[23:33]) % 1) * 60) % 1) // 1)),\n",
        "                \"%Y %j %H %M %S\"\n",
        "            )}\n",
        "\n",
        "            orb.update({\n",
        "                \"name\": tle2[2:7],\n",
        "                \"e\": float(\".\" + tle2[26:34]),\n",
        "                \"a\": (mu / ((2 * m.pi * float(tle2[52:63]) / (D * 3600)) ** 2)) ** (1. / 3),\n",
        "                \"i\": float(tle2[9:17]) * m.pi / 180,\n",
        "                \"RAAN\": float(tle2[17:26]) * m.pi / 180,\n",
        "                \"omega\": float(tle2[34:43]) * m.pi / 180\n",
        "            })\n",
        "\n",
        "            orb.update({\"b\": orb[\"a\"] * m.sqrt(1 - orb[\"e\"] ** 2),\n",
        "                        \"c\": orb[\"a\"] * orb[\"e\"]})\n",
        "\n",
        "            R = np.matmul(np.array([[m.cos(orb[\"RAAN\"]), -m.sin(orb[\"RAAN\"]), 0],\n",
        "                                     [m.sin(orb[\"RAAN\"]), m.cos(orb[\"RAAN\"]), 0],\n",
        "                                     [0, 0, 1]]),\n",
        "                          (np.array([[1, 0, 0],\n",
        "                                     [0, m.cos(orb[\"i\"]), -m.sin(orb[\"i\"])],\n",
        "                                     [0, m.sin(orb[\"i\"]), m.cos(orb[\"i\"])]])))\n",
        "            R = np.matmul(R, np.array([[m.cos(orb[\"omega\"]), -m.sin(orb[\"omega\"]), 0],\n",
        "                                     [m.sin(orb[\"omega\"]), m.cos(orb[\"omega\"]), 0],\n",
        "                                     [0, 0, 1]]))\n",
        "\n",
        "            # Ensure x, y, z are NumPy arrays and not nested lists\n",
        "            x = []\n",
        "            y = []\n",
        "            z = []\n",
        "            for i in np.linspace(0, 2 * m.pi, 100):\n",
        "                P = np.matmul(R, np.array([[orb[\"a\"] * m.cos(i)],\n",
        "                                         [orb[\"b\"] * m.sin(i)],\n",
        "                                         [0]])) - np.matmul(R, np.array([[orb[\"c\"]],\n",
        "                                                                         [0],\n",
        "                                                                         [0]]))\n",
        "                x.append(P[0][0])  # Extract the value from the nested array\n",
        "                y.append(P[1][0])  # Extract the value from the nested array\n",
        "                z.append(P[2][0])  # Extract the value from the nested array\n",
        "\n",
        "            x = np.array(x)\n",
        "            y = np.array(y)\n",
        "            z = np.array(z)\n",
        "\n",
        "\n",
        "            line.set_data(x[:frame], y[:frame])\n",
        "            line.set_3d_properties(z[:frame])\n",
        "\n",
        "        return lines\n",
        "\n",
        "    ani = FuncAnimation(fig, update, frames=100, init_func=init, blit=True)\n",
        "\n",
        "    plt.title(\"Orbits plotted in the ECE frame\")\n",
        "    if len(data) < 5:\n",
        "        ax.legend()\n",
        "    else:\n",
        "        fig.subplots_adjust(right=0.8)\n",
        "        ax.legend(loc='center left', bbox_to_anchor=(1.07, 0.5), fontsize=7)\n",
        "        #A legend is added so we differentiate between the different satellite orbits\n",
        "\n",
        "    return ani\n",
        "\n",
        "# Define the TLE data as a dictionary\n",
        "tle_data = {\n",
        "    \"LILACSAT-2\": (\n",
        "        '1 40908U 15049K   25040.75362640  .00014926  00000-0  44552-3 0  9999',\n",
        "        '2 40908  97.5165  58.4844 0008714 220.7406 139.3186 15.34689321519828'\n",
        "    ),\n",
        "    \"IO-86\": (\n",
        "        '1 40931U 15052B   25038.81889873  .00001817  00000-0  15548-3 0  9993',\n",
        "        '2 40931   6.0001 108.1159 0012838 187.7041 172.2889 14.78568669506492'\n",
        "    ),\n",
        "    \"Horyu-4\": (\n",
        "        '1 41340U 16012D   25040.50964053  .00024423  00000-0  83068-3 0  9998',\n",
        "        '2 41340  30.9968 215.0775 0005942 258.1882 101.8095 15.29227906494341'\n",
        "    ),\n",
        "    \"Lapan A3\": (\n",
        "        '1 41603U 16040E   25040.82089363  .00005724  00000-0  20899-3 0  9994',\n",
        "        '2 41603  97.1391  50.2222 0010651 346.8368  13.2593 15.28446454479488'\n",
        "    ),\n",
        "    \"CAS-2T\": (\n",
        "        '1 41847U 16066G   25040.81371518  .00003589  00000-0  52043-3 0  9990',\n",
        "        '2 41847  98.4339 116.3024 0348372 180.5013 179.5879 14.42893172433742'\n",
        "    ),\n",
        "    \"CAS-4B\": (\n",
        "        '1 42759U 17034B   25040.80858593  .00043921  00000-0  84809-3 0  9990',\n",
        "        '2 42759  43.0123 306.9566 0016345 180.6949 179.3890 15.46942498423955'\n",
        "    ),\n",
        "    \"CAS-4A\": (\n",
        "        '1 42761U 17034D   25040.83991301  .00043338  00000-0  82365-3 0  9998',\n",
        "        '2 42761  43.0161 304.7417 0017725 187.1675 172.8933 15.47379287423971'\n",
        "    ),\n",
        "    \"TechnoSat\": (\n",
        "        '1 42829U 17042E   25040.39008121  .00004774  00000-0  38520-3 0  9997',\n",
        "        '2 42829  97.3894 214.9223 0012152 342.6431  17.4377 15.00254609412697'\n",
        "    ),\n",
        "    \"AO-91\": (\n",
        "        '1 43017U 17073E   25040.58219161  .00013379  00000-0  72834-3 0  9997',\n",
        "        '2 43017  97.5498 274.4797 0181459 214.7102 144.2175 15.01471107391187'\n",
        "    ),\n",
        "    \"S-Net D\": (\n",
        "        '1 43186U 18014G   25040.29591181  .00007102  00000-0  44445-3 0  9997',\n",
        "        '2 43186  97.5738 304.5247 0007708 281.1198  78.9164 15.09684555384211'\n",
        "    ),\n",
        "    \"S-Net B\": (\n",
        "        '1 43187U 18014H   25040.29990714  .00006884  00000-0  43107-3 0  9993',\n",
        "        '2 43187  97.5744 304.4942 0008350 275.0534  84.9742 15.09670862384215'\n",
        "    ),\n",
        "    \"S-Net A\": (\n",
        "        '1 43188U 18014J   25040.69435521  .00008315  00000-0  51936-3 0  9998',\n",
        "        '2 43188  97.5751 305.0192 0008334 274.3627  85.6650 15.09697017384278'\n",
        "    ),\n",
        "    \"S-Net C\": (\n",
        "        '1 43189U 18014K   25040.63344299  .00007613  00000-0  47604-3 0  9991',\n",
        "        '2 43189  97.5769 304.9524 0009830 265.3769  94.6337 15.09675364382815'\n",
        "    ),\n",
        "    \"Ten-Koh\": (\n",
        "        '1 43677U 18084G   25040.72854077  .00005254  00000-0  43457-3 0  9991',\n",
        "        '2 43677  98.0827 208.2460 0011283 307.3159  52.7033 14.99330140342270'\n",
        "    )\n",
        "}\n",
        "\n",
        "# Call the plotting function with the TLE data\n",
        "ani = plot_tle(tle_data)\n",
        "\n",
        "# Display the animation\n",
        "HTML(ani.to_jshtml())"
      ],
      "metadata": {
        "id": "HvPKmRUNkbeg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#4\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from matplotlib.animation import FuncAnimation\n",
        "from IPython.display import HTML\n",
        "\n",
        "# Constants\n",
        "mu = 398600.4418  # Earth's gravitational parameter (km^3/s^2)\n",
        "r_earth = 6378  # Radius of Earth (km)\n",
        "r_orbit = 700  # Approximate altitude of satellites (km)\n",
        "r = r_earth + r_orbit  # Total orbital radius\n",
        "\n",
        "# Define satellite orbits\n",
        "num_satellites = 3\n",
        "angles = np.linspace(0, 2 * np.pi, 100)\n",
        "orbits = [\n",
        "    (r * np.cos(angles), r * np.sin(angles), np.zeros_like(angles)),  # Circular orbit\n",
        "    (r * np.cos(angles), np.zeros_like(angles), r * np.sin(angles)),  # Perpendicular orbit\n",
        "    (r * np.cos(angles), r * np.sin(angles) * np.cos(np.pi / 4), r * np.sin(angles) * np.sin(np.pi / 4))  # Tilted orbit\n",
        "]\n",
        "\n",
        "# **Find the Collision Point** - The position of satellites at the collision frame\n",
        "collision_frame = 50  # Frame where collision happens\n",
        "collision_positions = [\n",
        "    (orbits[i][0][collision_frame], orbits[i][1][collision_frame], orbits[i][2][collision_frame])\n",
        "    for i in range(num_satellites)\n",
        "]\n",
        "collision_point = np.mean(collision_positions, axis=0)  # Average the positions to approximate the collision point\n",
        "#The collision point is approximated by averaging the positions of the three satellites at that frame.\n",
        "#This point becomes the origin of the explosion and debris.\n",
        "\n",
        "# Explosion and debris parameters\n",
        "debris_spread = 50  # Speed of debris spread\n",
        "fireball_lifetime = 15  # How long the explosion fire lasts\n",
        "num_debris = 200  # Number of debris particles\n",
        "#A collision is simulated at a specific animation frame (collision_frame = 50).\n",
        "\n",
        "# Initialize figure\n",
        "fig = plt.figure()\n",
        "ax = plt.axes(projection='3d')\n",
        "\n",
        "# Plot Earth\n",
        "u, v = np.mgrid[0:2 * np.pi:20j, 0:np.pi:10j]\n",
        "ax.plot_wireframe(r_earth * np.cos(u) * np.sin(v),\n",
        "                   r_earth * np.sin(u) * np.sin(v),\n",
        "                   r_earth * np.cos(v), color=\"b\", alpha=0.3, lw=0.5)\n",
        "\n",
        "# Initialize satellite paths\n",
        "lines = [ax.plot([], [], [], label=f\"Satellite {i+1}\")[0] for i in range(num_satellites)]\n",
        "debris_points, = ax.plot([], [], [], 'ro', markersize=3, alpha=0.8)  # Debris particles\n",
        "fireball, = ax.plot([], [], [], 'yo', markersize=6, alpha=0.9)  # Fire explosion effect\n",
        "\n",
        "# Pre-generate random explosion directions\n",
        "debris_directions = np.random.normal(size=(num_debris, 3))\n",
        "debris_directions /= np.linalg.norm(debris_directions, axis=1).reshape(-1, 1)  # Normalize to unit vectors\n",
        "\n",
        "# Animation functions\n",
        "def init():\n",
        "    for line in lines:\n",
        "        line.set_data([], [])\n",
        "        line.set_3d_properties([])\n",
        "    debris_points.set_data([], [])\n",
        "    debris_points.set_3d_properties([])\n",
        "    fireball.set_data([], [])\n",
        "    fireball.set_3d_properties([])\n",
        "    return lines + [debris_points, fireball]\n",
        "\n",
        "#To animate it, first the 3D plot is initizaled with Earth's wireframe\n",
        "#Each satellite gets its own line plotted and updated frame-by-frame.\n",
        "#Explosion and debris are updated conditionally depending on the frame:\n",
        "#Before collision: satellites are shown in orbit.\n",
        "#At and after collision:, Satellites disappear, Fireball expands, Debris scatters from the collision point.\n",
        "#update(frame): Animates satellite orbits, Triggers explosion at frame 50, Expands fireball effect and debris cloud, Ensures smooth transition from orbiting satellites to chaos.\n",
        "\n",
        "def update(frame):\n",
        "    for i, line in enumerate(lines):\n",
        "        x, y, z = orbits[i]\n",
        "        if frame < collision_frame:\n",
        "            line.set_data(x[:frame], y[:frame])\n",
        "            line.set_3d_properties(z[:frame])\n",
        "        else:\n",
        "            line.set_data(x[:collision_frame], y[:collision_frame])\n",
        "            line.set_3d_properties(z[:collision_frame])\n",
        "\n",
        "\n",
        "    if collision_frame <= frame < collision_frame + fireball_lifetime:\n",
        "        # Expanding fireball effect centered at the actual collision point\n",
        "        num_flames = 100\n",
        "        explosion_radius = 100 * (frame - collision_frame) / fireball_lifetime\n",
        "        theta = np.random.uniform(0, 2 * np.pi, num_flames)\n",
        "        phi = np.random.uniform(0, np.pi, num_flames)\n",
        "        x_fire = collision_point[0] + explosion_radius * np.sin(phi) * np.cos(theta)\n",
        "        y_fire = collision_point[1] + explosion_radius * np.sin(phi) * np.sin(theta)\n",
        "        z_fire = collision_point[2] + explosion_radius * np.cos(phi)\n",
        "        fireball.set_data(x_fire, y_fire)\n",
        "        fireball.set_3d_properties(z_fire)\n",
        "    else:\n",
        "        fireball.set_data([], [])\n",
        "        fireball.set_3d_properties([])\n",
        "\n",
        "    if frame >= collision_frame:\n",
        "        # Debris spreading from the **correct collision point**\n",
        "        debris_distance = debris_spread * (frame - collision_frame)\n",
        "        x_debris = collision_point[0] + debris_distance * debris_directions[:, 0]\n",
        "        y_debris = collision_point[1] + debris_distance * debris_directions[:, 1]\n",
        "        z_debris = collision_point[2] + debris_distance * debris_directions[:, 2]\n",
        "        debris_points.set_data(x_debris, y_debris)\n",
        "        debris_points.set_3d_properties(z_debris)\n",
        "\n",
        "    return lines + [debris_points, fireball]\n",
        "\n",
        "# Run animation\n",
        "ani = FuncAnimation(fig, update, frames=100, init_func=init, blit=True)\n",
        "\n",
        "# Display animation\n",
        "plt.title(\"Satellite Collision with Fire and Debris\")\n",
        "ax.legend()\n",
        "HTML(ani.to_jshtml())\n"
      ],
      "metadata": {
        "id": "dOXhDxbW-fLl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#5a\n",
        "!pip install matplotlib numpy sgp4\n",
        "#The SGP4 orbital model is the industry-standard algorithm used to compute a satellite's position and velocity from a TLE at a given time."
      ],
      "metadata": {
        "id": "NeJFkmQOt5wg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#5b\n",
        "from sgp4.api import Satrec, jday\n",
        "from sgp4.conveniences import sat_epoch_datetime\n",
        "import numpy as np\n",
        "from datetime import datetime\n",
        "import matplotlib.pyplot as plt\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "from matplotlib.animation import FuncAnimation\n",
        "from IPython.display import HTML\n",
        "#sgp4.api, sat_epoch_datetime: Used to load and propagate TLEs.\n",
        "#datetime, numpy: For date/time and numerical operations.\n",
        "\n",
        "# ===================\n",
        "# UTILITY FUNCTIONS\n",
        "# ===================\n",
        "\n",
        "def calculate_satellite_positions(sat, jd, fr):\n",
        "    e, r, v = sat.sgp4(jd, fr)\n",
        "    return e, np.array(r), np.array(v)\n",
        "\n",
        "def calculate_relative_position_velocity(r1, v1, r2, v2):\n",
        "    return np.subtract(r1, r2), np.subtract(v1, v2)\n",
        "\n",
        "def calculate_distance(relative_position):\n",
        "    return np.linalg.norm(relative_position)\n",
        "\n",
        "def estimate_collision_risk(distance, relative_velocity, safe_distance=10.0, max_velocity=3.0):\n",
        "    if distance < safe_distance:\n",
        "        velocity_magnitude = np.linalg.norm(relative_velocity)\n",
        "        return min(1.0, velocity_magnitude / max_velocity)\n",
        "    return 0.0\n",
        "\n",
        "# ===================\n",
        "# TLE DATA\n",
        "# ===================\n",
        "\n",
        "tle_cas4a = [\n",
        "    '1 42761U 17034D   25040.83991301  .00043338  00000-0  82365-3 0  9998',\n",
        "    '2 42761  43.0161 304.7417 0017725 187.1675 172.8933 15.47379287423971'\n",
        "]\n",
        "\n",
        "# We'll scan MA only (Mean Anomaly), keeping everything else aligned\n",
        "tle_cas4b_template = [\n",
        "    '1 42759U 17034B   25040.83991301  .00043921  00000-0  84809-3 0  9990',\n",
        "    '2 42759  43.0161 304.7417 0017725 187.1675 {:8.4f} 15.47379287423971'\n",
        "]\n",
        "\n",
        "# Load CAS-4A and get shared epoch time\n",
        "sat1 = Satrec.twoline2rv(*tle_cas4a)\n",
        "epoch = sat_epoch_datetime(sat1)\n",
        "jd, fr = jday(epoch.year, epoch.month, epoch.day, epoch.hour, epoch.minute, epoch.second)\n",
        "\n",
        "# ===================\n",
        "# BRUTE FORCE SCAN\n",
        "# ===================\n",
        "\n",
        "center_ma = 172.8933\n",
        "scan_range = np.linspace(center_ma - 0.1, center_ma + 0.1, 1000)\n",
        "\n",
        "distances = []\n",
        "best_ma = None\n",
        "min_distance = float('inf')\n",
        "best_sat2 = None\n",
        "best_rel_pos = None\n",
        "\n",
        "for ma in scan_range:\n",
        "    tle2_line2 = tle_cas4b_template[1].format(ma)\n",
        "    sat2 = Satrec.twoline2rv(tle_cas4b_template[0], tle2_line2)\n",
        "\n",
        "    e1, r1, v1 = calculate_satellite_positions(sat1, jd, fr)\n",
        "    e2, r2, v2 = calculate_satellite_positions(sat2, jd, fr)\n",
        "\n",
        "    if e1 == 0 and e2 == 0:\n",
        "        rel_pos, rel_vel = calculate_relative_position_velocity(r1, v1, r2, v2)\n",
        "        dist = calculate_distance(rel_pos)\n",
        "        distances.append(dist)\n",
        "\n",
        "        if dist < min_distance:\n",
        "            min_distance = dist\n",
        "            best_ma = ma\n",
        "            best_sat2 = sat2\n",
        "            best_rel_pos = rel_pos\n",
        "    else:\n",
        "        distances.append(np.nan)\n",
        "\n",
        "# ===================\n",
        "# REPORT RESULTS\n",
        "# ===================\n",
        "\n",
        "print(f\"\\n Closest Mean Anomaly for CAS-4B: {best_ma:.4f}°\")\n",
        "print(f\" Minimum Distance to CAS-4A: {min_distance:.4f} km\")\n",
        "risk = estimate_collision_risk(min_distance, rel_vel)\n",
        "print(f\" Collision Risk: {risk:.4f}\")\n",
        "if risk > 0.5:\n",
        "    print(\" High risk of collision detected!\")\n",
        "elif risk > 0:\n",
        "    print(\" Moderate risk of collision.\")\n",
        "else:\n",
        "    print(\" No significant collision risk.\")\n",
        "\n",
        "# ===================\n",
        "# PLOT DISTANCE VS. MEAN ANOMALY\n",
        "# ===================\n",
        "\n",
        "plt.figure(figsize=(10, 5))\n",
        "plt.plot(scan_range, distances, label=\"Distance (km)\")\n",
        "plt.axvline(best_ma, color='r', linestyle='--', label=\"Closest Approach\")\n",
        "plt.title(\"Minimum Distance vs. Mean Anomaly of CAS-4B\")\n",
        "plt.xlabel(\"Mean Anomaly (degrees)\")\n",
        "plt.ylabel(\"Distance to CAS-4A (km)\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "kUG51JvIvjo6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7fenTkjZa0oh"
      },
      "outputs": [],
      "source": [
        "#5c\n",
        "# Import necessary modules\n",
        "from sgp4.api import Satrec, jday\n",
        "from datetime import datetime\n",
        "import numpy as np  # Import numpy for array operations\n",
        "\n",
        "# Function to ensure positive values are formatted correctly\n",
        "def ensure_positive_sign(tle1, tle2):\n",
        "    \"\"\"\n",
        "    This function ensures that the inclination, RAAN, argument of perigee, and mean anomaly\n",
        "    have a + sign where appropriate by replacing any '-' with '+'.\n",
        "    \"\"\"\n",
        "    # Split TLE strings into lines, ensuring they have two elements\n",
        "    tle1_lines = tle1.splitlines() if isinstance(tle1, str) else tle1 # Handle cases where tle1 is already a list\n",
        "    tle2_lines = tle2.splitlines() if isinstance(tle2, str) else tle2 # Handle cases where tle2 is already a list\n",
        "\n",
        "    # Ensure tle2_lines has at least two elements before accessing index 1\n",
        "    if len(tle2_lines) > 1:\n",
        "        tle2_lines[1] = tle2_lines[1].replace('-', '+', 1)\n",
        "\n",
        "    return tle1_lines, tle2_lines\n",
        "\n",
        "# Define the TLE data for the satellites you want to test\n",
        "tle_data = {\n",
        "    \"LILACSAT-2\": (\n",
        "        '1 40908U 15049K   25040.75362640  .00014926  00000-0  44552-3 0  9999',\n",
        "        '2 40908  97.5165  58.4844 0008714 220.7406 139.3186 15.34689321519828'\n",
        "    ),\n",
        "    \"IO-86\": (\n",
        "        '1 40931U 15052B   25038.81889873  .00001817  00000-0  15548-3 0  9993',\n",
        "        '2 40931   6.0001 108.1159 0012838 187.7041 172.2889 14.78568669506492'\n",
        "    ),\n",
        "    \"Horyu-4\": (\n",
        "        '1 41340U 16012D   25040.50964053  .00024423  00000-0  83068-3 0  9998',\n",
        "        '2 41340  30.9968 215.0775 0005942 258.1882 101.8095 15.29227906494341'\n",
        "    ),\n",
        "    \"Lapan A3\": (\n",
        "        '1 41603U 16040E   25040.82089363  .00005724  00000-0  20899-3 0  9994',\n",
        "        '2 41603  97.1391  50.2222 0010651 346.8368  13.2593 15.28446454479488'\n",
        "    ),\n",
        "    \"CAS-2T\": (\n",
        "        '1 41847U 16066G   25040.81371518  .00003589  00000-0  52043-3 0  9990',\n",
        "        '2 41847  98.4339 116.3024 0348372 180.5013 179.5879 14.42893172433742'\n",
        "    ),\n",
        "    \"TechnoSat\": (\n",
        "        '1 42829U 17042E   25040.39008121  .00004774  00000-0  38520-3 0  9997',\n",
        "        '2 42829  97.3894 214.9223 0012152 342.6431  17.4377 15.00254609412697'\n",
        "    ),\n",
        "    \"AO-91\": (\n",
        "        '1 43017U 17073E   25040.58219161  .00013379  00000-0  72834-3 0  9997',\n",
        "        '2 43017  97.5498 274.4797 0181459 214.7102 144.2175 15.01471107391187'\n",
        "    ),\n",
        "    \"S-Net D\": (\n",
        "        '1 43186U 18014G   25040.29591181  .00007102  00000-0  44445-3 0  9997',\n",
        "        '2 43186  97.5738 304.5247 0007708 281.1198  78.9164 15.09684555384211'\n",
        "    ),\n",
        "    \"S-Net B\": (\n",
        "        '1 43187U 18014H   25040.29990714  .00006884  00000-0  43107-3 0  9993',\n",
        "        '2 43187  97.5744 304.4942 0008350 275.0534  84.9742 15.09670862384215'\n",
        "    ),\n",
        "    \"S-Net A\": (\n",
        "        '1 43188U 18014J   25040.69435521  .00008315  00000-0  51936-3 0  9998',\n",
        "        '2 43188  97.5751 305.0192 0008334 274.3627  85.6650 15.09697017384278'\n",
        "    ),\n",
        "    \"S-Net C\": (\n",
        "        '1 43189U 18014K   25040.63344299  .00007613  00000-0  47604-3 0  9991',\n",
        "        '2 43189  97.5769 304.9524 0009830 265.3769  94.6337 15.09675364382815'\n",
        "    ),\n",
        "    \"Ten-Koh\": (\n",
        "        '1 43677U 18084G   25040.72854077  .00005254  00000-0  43457-3 0  9991',\n",
        "        '2 43677  98.0827 208.2460 0011283 307.3159  52.7033 14.99330140342270'\n",
        "    )\n",
        "}\n",
        "\n",
        "# Function to calculate and propagate satellite positions\n",
        "def calculate_satellite_positions(tle1, tle2):\n",
        "    tle1_lines, tle2_lines = ensure_positive_sign(tle1, tle2)\n",
        "\n",
        "    # Convert back the lines to TLE format\n",
        "    tle1 = \"\\n\".join(tle1_lines)\n",
        "    tle2 = \"\\n\".join(tle2_lines)\n",
        "\n",
        "    sat = Satrec.twoline2rv(tle1, tle2)\n",
        "\n",
        "    current_time = datetime.utcnow()\n",
        "    jd, fr = jday(current_time.year, current_time.month, current_time.day,\n",
        "                  current_time.hour, current_time.minute, current_time.second)\n",
        "\n",
        "    # Propagate the satellite position and velocity\n",
        "    e, r, v = sat.sgp4(jd, fr)\n",
        "    return e, r, v\n",
        "\n",
        "# Function to calculate relative position and velocity\n",
        "def calculate_relative_position_velocity(r1, v1, r2, v2):\n",
        "    relative_position = np.subtract(r1, r2)\n",
        "    relative_velocity = np.subtract(v1, v2)\n",
        "    return relative_position, relative_velocity\n",
        "\n",
        "def calculate_distance(relative_position):\n",
        "    return np.linalg.norm(relative_position)\n",
        "\n",
        "def estimate_collision_risk(distance, relative_velocity, safe_distance=10.0, max_velocity=3.0):\n",
        "    \"\"\"Estimates the risk of collision based on distance and relative velocity.\n",
        "\n",
        "    Args:\n",
        "        distance (float): Distance between the two satellites in km.\n",
        "        relative_velocity (np.ndarray): Relative velocity between the two satellites in km/s.\n",
        "        safe_distance (float): Distance below which a collision is considered possible (in km).\n",
        "        max_velocity (float): Relative velocity above which a collision is considered certain (in km/s).\n",
        "\n",
        "    Returns:\n",
        "        float: Collision risk (0.0 to 1.0).\n",
        "    \"\"\"\n",
        "    collision_risk = 0.0\n",
        "    if distance < safe_distance:\n",
        "        velocity_magnitude = np.linalg.norm(relative_velocity)\n",
        "        if velocity_magnitude > max_velocity:\n",
        "            collision_risk = 1.0\n",
        "        else:\n",
        "            collision_risk = velocity_magnitude / max_velocity\n",
        "    return collision_risk\n",
        "\n",
        "# Function to display results for each satellite\n",
        "def display_satellite_results(sat_name, tle_data):\n",
        "    e, r, v = calculate_satellite_positions(tle_data[0], tle_data[1])\n",
        "    print(f\"Results for {sat_name}:\")\n",
        "    print(f\"Propagation status: {e}\")\n",
        "    print(f\"Position (r): {r} km\")\n",
        "    print(f\"Velocity (v): {v} km/s\")\n",
        "    return r, v\n",
        "\n",
        "# Test the satellites with the provided TLEs\n",
        "satellite_pairs = [\n",
        "    (\"LILACSAT-2\", tle_data[\"LILACSAT-2\"]),\n",
        "    (\"IO-86\", tle_data[\"IO-86\"]),\n",
        "    (\"Horyu-4\", tle_data[\"Horyu-4\"]),\n",
        "    (\"Lapan A3\", tle_data[\"Lapan A3\"]),\n",
        "    (\"CAS-2T\", tle_data[\"CAS-2T\"]),\n",
        "    (\"TechnoSat\", tle_data[\"TechnoSat\"]),\n",
        "    (\"AO-91\", tle_data[\"AO-91\"]),\n",
        "    (\"S-Net D\", tle_data[\"S-Net D\"]),\n",
        "    (\"S-Net B\", tle_data[\"S-Net B\"]),\n",
        "    (\"S-Net A\", tle_data[\"S-Net A\"]),\n",
        "    (\"S-Net C\", tle_data[\"S-Net C\"]),\n",
        "    (\"Ten-Koh\", tle_data[\"Ten-Koh\"])\n",
        "]\n",
        "\n",
        "# Loop through the satellite pairs and calculate collision risk\n",
        "for i in range(len(satellite_pairs)):\n",
        "    for j in range(i+1, len(satellite_pairs)):\n",
        "        sat1_name, tle1 = satellite_pairs[i]\n",
        "        sat2_name, tle2 = satellite_pairs[j]\n",
        "\n",
        "        # Calculate positions and velocities for both satellites\n",
        "        r1, v1 = display_satellite_results(sat1_name, tle1)\n",
        "        r2, v2 = display_satellite_results(sat2_name, tle2)\n",
        "\n",
        "        # Calculate relative position and velocity\n",
        "        relative_position, relative_velocity = calculate_relative_position_velocity(r1, v1, r2, v2)\n",
        "\n",
        "        # Calculate distance and collision risk\n",
        "        distance = calculate_distance(relative_position)\n",
        "        collision_risk = estimate_collision_risk(distance, relative_velocity)\n",
        "\n",
        "        print(f\"Collision risk between {sat1_name} and {sat2_name}: {collision_risk:.4f}\")\n",
        "        print(\"-\" * 50)\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#5d\n",
        "from sgp4.api import Satrec, jday\n",
        "from datetime import datetime\n",
        "import numpy as np\n",
        "\n",
        "# Realistic collision-prone pair with a slight mean anomaly offset\n",
        "tle_data = {\n",
        "    \"S-Net B\": (\n",
        "        '1 43187U 18014H   25040.29990714  .00006884  00000-0  43107-3 0  9993',\n",
        "        '2 43187  97.5744 304.4942 0008350 275.0534  84.9742 15.09670862384215'\n",
        "    ),\n",
        "    \"S-Net A\": (\n",
        "        '1 43188U 18014J   25040.29990714  .00006884  00000-0  43107-3 0  9993',\n",
        "        '2 43188  97.5744 304.4942 0008350 275.0534  84.9842 15.09670862384215'\n",
        "    )\n",
        "}\n",
        "\n",
        "def calculate_satellite_positions(tle1, tle2):\n",
        "    sat = Satrec.twoline2rv(tle1, tle2)\n",
        "    now = datetime.utcnow()\n",
        "    jd, fr = jday(now.year, now.month, now.day, now.hour, now.minute, now.second)\n",
        "    e, r, v = sat.sgp4(jd, fr)\n",
        "    return e, np.array(r), np.array(v)\n",
        "\n",
        "def calculate_relative_position_velocity(r1, v1, r2, v2):\n",
        "    return r1 - r2, v1 - v2\n",
        "\n",
        "def calculate_distance(relative_position):\n",
        "    return np.linalg.norm(relative_position)\n",
        "\n",
        "def estimate_collision_risk(distance, relative_velocity, safe_distance=10.0, max_velocity=3.0):\n",
        "    if distance == 0.0:\n",
        "        return 1.0\n",
        "    if distance < safe_distance:\n",
        "        velocity_magnitude = np.linalg.norm(relative_velocity)\n",
        "        return 1.0 if velocity_magnitude > max_velocity else velocity_magnitude / max_velocity\n",
        "    return 0.0\n",
        "\n",
        "def apply_delta_v(velocity, delta_v_vector):\n",
        "    return velocity + delta_v_vector\n",
        "\n",
        "# === Main Simulation ===\n",
        "sat1_name, tle1 = \"S-Net A\", tle_data[\"S-Net A\"]\n",
        "sat2_name, tle2 = \"S-Net B\", tle_data[\"S-Net B\"]\n",
        "\n",
        "e1, r1, v1 = calculate_satellite_positions(tle1[0], tle1[1])\n",
        "e2, r2, v2 = calculate_satellite_positions(tle2[0], tle2[1])\n",
        "\n",
        "print(f\"\\n--- {sat1_name} ---\\nPosition: {r1}\\nVelocity: {v1}\")\n",
        "print(f\"\\n--- {sat2_name} ---\\nPosition: {r2}\\nVelocity: {v2}\")\n",
        "\n",
        "relative_position, relative_velocity = calculate_relative_position_velocity(r1, v1, r2, v2)\n",
        "distance = calculate_distance(relative_position)\n",
        "collision_risk = estimate_collision_risk(distance, relative_velocity)\n",
        "\n",
        "print(\"\\n--- Relative Stats ---\")\n",
        "print(f\"Relative Position Vector (km): {relative_position}\")\n",
        "print(f\"Relative Velocity Vector (km/s): {relative_velocity}\")\n",
        "print(f\"Distance Between Satellites: {distance:.6f} km\")\n",
        "print(f\"Relative Speed: {np.linalg.norm(relative_velocity)*1000:.2f} m/s\")\n",
        "print(f\"Collision Risk Score: {collision_risk:.4f}\")\n",
        "\n",
        "# === Apply small delta-v to mitigate collision ===\n",
        "if collision_risk > 0:\n",
        "    print(\"Applying small avoidance maneuver (∆v)...\")\n",
        "    delta_v = np.array([0.0, 0.01, 0.0])  # 10 m/s cross-track bump\n",
        "    v1_adjusted = apply_delta_v(v1, delta_v)\n",
        "    relative_position_new, relative_velocity_new = calculate_relative_position_velocity(r1, v1_adjusted, r2, v2)\n",
        "    distance_new = calculate_distance(relative_position_new)\n",
        "    collision_risk_new = estimate_collision_risk(distance_new, relative_velocity_new)\n",
        "\n",
        "    print(\"\\n--- After Avoidance Maneuver ---\")\n",
        "    print(f\"New Relative Velocity Vector: {relative_velocity_new}\")\n",
        "    print(f\"New Distance Between Satellites: {distance_new:.6f} km\")\n",
        "    print(f\"New Collision Risk Score: {collision_risk_new:.4f}\")\n",
        "\n",
        "    if collision_risk_new < collision_risk:\n",
        "        print(\"aneuver Successful: Collision Avoided\")\n",
        "    else:\n",
        "        print(\"Maneuver Ineffective: Further Planning Required\")\n",
        "\n",
        "else:\n",
        "    print(\"No Avoidance Needed\")\n"
      ],
      "metadata": {
        "id": "VTp9GuASPReA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#6a\n",
        "from sgp4.api import Satrec, jday\n",
        "from datetime import datetime, timedelta\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "#sgp4.api, jday Used to load and propagate satellite orbits from TLEs.\n",
        "#datetime, timedelta For defining the time range of the simulation.\n",
        "#numpy For calculating vector distances between satellites.\n",
        "#matplotlib.pyplot To plot the distance between the satellites over time.\n",
        "\n",
        "# Use mid TLE for TIMED and Cosmos 2221\n",
        "tle_TIMED = (\n",
        "    \"1 26998U 01055B   24059.55288321  .00001515  00000-0  15772-3 0  9996\",\n",
        "    \"2 26998  74.0701 120.4821 0001965 249.2297 110.8675 14.90837151205485\"\n",
        ")\n",
        "\n",
        "tle_COSMOS = (\n",
        "    \"1 22236U 92080A   24059.56507143  .00002036  00000-0  20268-3 0  9997\",\n",
        "    \"2 22236  82.5028 265.2326 0016263 147.5211 212.7013 14.91072594670299\"\n",
        ")\n",
        "\n",
        "# Create satellite models\n",
        "sat_TIMED = Satrec.twoline2rv(*tle_TIMED)\n",
        "sat_COSMOS = Satrec.twoline2rv(*tle_COSMOS)\n",
        "\n",
        "# Propagation time window: 3 hours, every 10 seconds\n",
        "start_time = datetime(2024, 2, 28, 0, 0, 0)\n",
        "seconds = range(0, 3 * 60 * 60, 1)  # 3 hours in 10-second steps\n",
        "timestamps = [start_time + timedelta(seconds=i) for i in seconds]\n",
        "jd_fr_list = [jday(t.year, t.month, t.day, t.hour, t.minute, t.second + t.microsecond * 1e-6) for t in timestamps]\n",
        "\n",
        "#For each second in the time window:\n",
        "#Use sat.sgp4(jd, fr) to calculate the satellite’s 3D position.\n",
        "#If both satellites propagate successfully (e == 0), calculate the distance between them.\n",
        "#If propagation fails, store NaN to avoid errors in analysis.\n",
        "#After the full 3-hour simulation: Identify the shortest distance between the two satellites, Record the exact time it occurred.\n",
        "#Print both values in the output.\n",
        "\n",
        "positions_TIMED = []\n",
        "positions_COSMOS = []\n",
        "distances = []\n",
        "\n",
        "for jd, fr in jd_fr_list:\n",
        "    e1, r1, _ = sat_TIMED.sgp4(jd, fr)\n",
        "    e2, r2, _ = sat_COSMOS.sgp4(jd, fr)\n",
        "    if e1 == 0 and e2 == 0:\n",
        "        positions_TIMED.append(r1)\n",
        "        positions_COSMOS.append(r2)\n",
        "        distances.append(np.linalg.norm(np.array(r1) - np.array(r2)))\n",
        "    else:\n",
        "        positions_TIMED.append([np.nan] * 3)\n",
        "        positions_COSMOS.append([np.nan] * 3)\n",
        "        distances.append(np.nan)\n",
        "\n",
        "# Results\n",
        "min_distance = np.nanmin(distances)\n",
        "min_index = distances.index(min_distance)\n",
        "closest_time = timestamps[min_index]\n",
        "\n",
        "print(f\" Closest approach on: {closest_time}\")\n",
        "print(f\" Minimum distance: {min_distance:.4f} km\")\n",
        "\n",
        "# Plot distance over time\n",
        "plt.figure(figsize=(10, 5))\n",
        "plt.plot(timestamps, distances, label=\"Separation Distance\")\n",
        "plt.axvline(closest_time, color='red', linestyle='--', label='Closest Approach')\n",
        "plt.title(\"Distance Between TIMED and Cosmos 2221 (Feb 28, 2024)\")\n",
        "plt.ylabel(\"Distance (km)\")\n",
        "plt.xlabel(\"Time (UTC)\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.tight_layout()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "id": "-tu5t58JGj3c"
      },
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}