PoC_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#move back to the root directory\n",
    "import os\n",
    "os.chdir('..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "!pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading file from: https://gitlab.orekit.org/orekit/orekit-data/-/archive/master/orekit-data-master.zip\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OpenJDK 64-Bit Server VM warning: Attempt to protect stack guard pages failed.\n",
      "OpenJDK 64-Bit Server VM warning: Attempt to deallocate stack guard pages failed.\n"
     ]
    }
   ],
   "source": [
    "import orekit\n",
    "from orekit.pyhelpers import setup_orekit_curdir\n",
    "orekit.pyhelpers.download_orekit_data_curdir()\n",
    "vm = orekit.initVM()\n",
    "setup_orekit_curdir()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from orekit.pyhelpers import datetime_to_absolutedate\n",
    "from org.hipparchus.geometry.euclidean.threed import Vector3D\n",
    "from org.hipparchus.linear import RealMatrix\n",
    "from org.orekit.propagation import StateCovariance\n",
    "from org.orekit.frames import FramesFactory\n",
    "from org.orekit.orbits import PositionAngleType, CartesianOrbit, KeplerianOrbit, CircularOrbit, EquinoctialOrbit, OrbitType\n",
    "from org.orekit.ssa.metrics import ProbabilityOfCollision \n",
    "from org.orekit.ssa.collision.shorttermencounter.probability.twod import Patera2005\n",
    "from org.orekit.utils import PVCoordinates\n",
    "from org.orekit.utils import Constants\n",
    "import datetime\n",
    "import random\n",
    "from org.hipparchus.linear import MatrixUtils\n",
    "from org.orekit.utils import IERSConventions, PVCoordinates\n",
    "from org.hipparchus.linear import MatrixUtils\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# # Define the covariance matrix",
    "# covariance_values1 = [\n",
    "#     [127.5401258877026, -216.5469301895778, 23.41184551274494, 0.2482739312599897, -0.1375882369099475, 0.008141795305978281],\n",
    "#     [-216.5469301895778, 9738.762752029239, -5.346171445289681, -10.55496325674788, 0.1222189661319367, -0.05578838789977865],\n",
    "#     [23.41184551274494, -5.346171445289681, 24.60870138594973, 0.005950144116086071, -0.02576925085939057, -0.0003344474012710994],\n",
    "#     [0.2482739312599897, -10.55496325674788, 0.005950144116086071, 0.01147849769311488, -0.0001458466399925159, 0.00005801069817622688],\n",
    "#     [-0.1375882369099475, 0.1222189661319367, -0.02576925085939057, -0.0001458466399925159, 0.0001499094267380701, -0.000008757429553163563],\n",
    "#     [0.008141795305978281, -0.05578838789977865, -0.0003344474012710994, 0.00005801069817622688, -0.000008757429553163563, 0.00009287778622823746]\n",
    "# ]\n",
    "\n",
    "# covariance_values2 = [\n",
    "#     [4.7440894789163000000000, -1.2583279067770000000000, -1.2583279067770000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000],\n",
    "#     [-1.2583279067770000000000,6.1279552605419000000000, 2.1279552605419000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000],\n",
    "#     [-1.2583279067770000000000, 2.1279552605419000000000, 6.1279552605419000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000],\n",
    "#     [0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000010000000000000000, 0.0000000000000000000000, 0.0000000000000000000000],\n",
    "#     [0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000001, 0.0000010000000000000000, -0.0000000000000000000001],\n",
    "#     [0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000000, 0.0000000000000000000001, -0.0000000000000000000001, 0.0000010000000000000000]\n",
    "# ]\n",
    "\n",
    "covariance_values3 = [\n",
    "    [-5.826e-1, 2.538e-2, -2.476e-6, -1.628e-4, -1.782e-4, 1.605e-4],\n",
    "    [2.538e-2, 7.537e-1, -8.935e-2, -2.343e-4, -7.149e-4, 5.660e-4],\n",
    "    [-2.476e-6, -8.935e-2, 9.269e-1, 2.591e-4, 6.95e-4, -7.503e-4],\n",
    "    [-1.628e-4, -2.343e-4, 2.591e-4,2.591e-7, 4.042e-7, -3.707e-7],\n",
    "    [-1.782e-4, -7.149e-4, 6.95e-4, 4.042e-7, 1.198e-6, -9.648e-7],\n",
    "    [1.605e-4, 5.660e-4, -7.503e-4, -3.707e-7, -9.648e-7, 1.066e-6]]\n",
    "\n",
    "# covariance_values3 = [\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],\n",
    "#     [1e-1, 1e-1 ,1e-1, 1e-1,1e-1, 1e-1],]\n",
    "\n",
    "# Convert the Python list to a Java double[][]\n",
    "jarray = MatrixUtils.createRealMatrix(len(covariance_values3), len(covariance_values3[0]))\n",
    "for i in range(len(covariance_values3)):\n",
    "    for j in range(len(covariance_values3[i])):\n",
    "        try:\n",
    "            jarray.setEntry(i, j, covariance_values3[i][j])\n",
    "        except:\n",
    "            print(i, j, covariance_values3[i][j])\n",
    "            print(f\"value: {covariance_values3[i][j]}\")\n",
    "            raise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use MatrixUtils to create a RealMatrix object from the 2D array\n",
    "cov_mat1 = jarray\n",
    "cov_mat2  = jarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def symmetrize(matrix):\n",
    "    nrows, ncols = matrix.getRowDimension(), matrix.getColumnDimension()\n",
    "    for i in range(nrows):\n",
    "        for j in range(i+1, ncols):\n",
    "            value = (matrix.getEntry(i, j) + matrix.getEntry(j, i)) / 2.0\n",
    "            matrix.setEntry(i, j, value)\n",
    "            matrix.setEntry(j, i, value)\n",
    "\n",
    "# Apply symmetrization covariance\n",
    "symmetrize(cov_mat1)\n",
    "symmetrize(cov_mat2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of collision: 0.1940657070465181\n"
     ]
    }
   ],
   "source": [
    "def generate_symmetric_covariance_matrix(size):\n",
    "    # Create an empty matrix of the desired size\n",
    "    matrix = MatrixUtils.createRealMatrix(size, size)\n",
    "    \n",
    "    # Fill the upper triangle and diagonal with random values\n",
    "    for i in range(size):\n",
    "        for j in range(i, size):\n",
    "            value = random.random()\n",
    "            matrix.setEntry(i, j, value)\n",
    "            matrix.setEntry(j, i, value)  # Mirror to make it symmetric\n",
    "    \n",
    "    return matrix\n",
    "\n",
    "state1 = [700001.0, 0.0, 0.0, 0.0, 7000.24, 0.0]\n",
    "state2 = [700000.0, 0.0, 0.0, 0.0, 7000.1, 0.0]\n",
    "tca = datetime.datetime(2020, 1, 1, 0, 0, 0)\n",
    "\n",
    "orbit1 = CartesianOrbit(PVCoordinates(Vector3D(float(state1[0]), float(state1[1]), float(state1[2])),\n",
    "                        Vector3D(float(state1[3]), float(state1[4]), float(state1[5]))),\n",
    "                        FramesFactory.getEME2000(),\n",
    "                        datetime_to_absolutedate(tca),\n",
    "                        Constants.WGS84_EARTH_MU)\n",
    "\n",
    "orbit2 = CartesianOrbit(PVCoordinates(Vector3D(float(state2[0]), float(state2[1]), float(state2[2])),\n",
    "                        Vector3D(float(state2[3]), float(state2[4]), float(state2[5]))),\n",
    "                        FramesFactory.getEME2000(),\n",
    "                        datetime_to_absolutedate(tca),\n",
    "                        Constants.WGS84_EARTH_MU)\n",
    "\n",
    "radius1 = 1.0\n",
    "radius2 = 1.0\n",
    "\n",
    "covariance1 = StateCovariance(cov_mat1, datetime_to_absolutedate(tca), FramesFactory.getITRF(IERSConventions.IERS_2010, False), OrbitType.CARTESIAN, PositionAngleType.TRUE)\n",
    "covariance2 = StateCovariance(cov_mat2, datetime_to_absolutedate(tca), FramesFactory.getITRF(IERSConventions.IERS_2010, False), OrbitType.CARTESIAN, PositionAngleType.TRUE)\n",
    "\n",
    "# Patera2005.compute(Orbit primaryAtTCA, StateCovariance primaryCovariance, double primaryRadius, Orbit secondaryAtTCA, StateCovariance secondaryCovariance, double secondaryRadius)\n",
    "patera2005 = Patera2005() \n",
    "poc_result = patera2005.compute(orbit1, covariance1, orbit2, covariance2, radius2, 1e-10)\n",
    "print(f\"Probability of collision: {poc_result.getValue()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#state covariance constructor takes in:\n",
    "#   a 6x6 matrix (RealMatrix)\n",
    "#   an epoch\n",
    "#   a local orbital frame\n",
    "\n",
    "# org.hipparchus.geometry.euclidean.threed.Rotation\trotationFromInertial(AbsoluteDate date, PVCoordinates pv)\n",
    "# Get the rotation from inertial frame to local orbital frame."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "erp_tools_env",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}