diff --git a/data/.gitattributes b/data/.gitattributes index 140ae139..74d1bf13 100644 --- a/data/.gitattributes +++ b/data/.gitattributes @@ -1,5 +1,4 @@ -pos.csv filter=lfs diff=lfs merge=lfs -text -v.csv filter=lfs diff=lfs merge=lfs -text +initial_population.csv filter=lfs diff=lfs merge=lfs -text */pos_*.csv filter=lfs diff=lfs merge=lfs -text */v_*.csv filter=lfs diff=lfs merge=lfs -text */shells_*.txt filter=lfs diff=lfs merge=lfs -text diff --git a/data/initial_population.csv b/data/initial_population.csv new file mode 100644 index 00000000..0e244658 --- /dev/null +++ b/data/initial_population.csv @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0aa44614a0ed365b60a32ca236b1ebd45857b97e332c9d23df7c9371fb192e69 +size 2993390 diff --git a/data/pos.csv b/data/pos.csv deleted file mode 100644 index acbbc993..00000000 --- a/data/pos.csv +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:4ccad2d454df1de6116942000582de7665cd43efa865541ba9f770ba13e53e29 -size 1329304 diff --git a/data/v.csv b/data/v.csv deleted file mode 100644 index 9b30043a..00000000 --- a/data/v.csv +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:5dcf21b8eb05ce8fff43db896d2934f64e56ff188434036b1b6329c666097b58 -size 1329518 diff --git a/notebooks/Data Processing.ipynb b/notebooks/Data Processing.ipynb deleted file mode 100644 index 06bb449d..00000000 --- a/notebooks/Data Processing.ipynb +++ /dev/null @@ -1,256 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "0a48e2c9", - "metadata": {}, - "outputs": [], - "source": [ - "### Imports\n", - "%load_ext autoreload\n", - "%autoreload 2\n", - "\n", - "# Append main folder\n", - "import sys\n", - "sys.path.append(\"../\")\n", - "\n", - "import pykep as pk\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "\n", - "starting_t = pk.epoch_from_string('2022-01-01 00:00:00.000')\n", - "lower_cutoff_in_km = 6371 + 200 # Earth radius + ...\n", - "higher_cutoff_in_km = 6371 + 2000" - ] - }, - { - "cell_type": "markdown", - "id": "de523087", - "metadata": {}, - "source": [ - "## Read TLE data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cf98e917", - "metadata": {}, - "outputs": [], - "source": [ - "planet_list = pk.util.read_tle(\"../data/tle.txt\")\n", - "#satcat = pykep.util.read_satcat(\"../data/satcat.txt\")\n", - "print(\"Loaded \",len(planet_list),\" planets.\")" - ] - }, - { - "cell_type": "markdown", - "id": "8cd629fd", - "metadata": {}, - "source": [ - "## Plot some examples" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9f3f3a16", - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure(figsize=(6,6),dpi=100)\n", - "ax = plt.axes(projection='3d');\n", - "for i in range (10):\n", - " pk.orbit_plots.plot_planet(planet_list[i],axes=ax)" - ] - }, - { - "cell_type": "markdown", - "id": "dbd01df1", - "metadata": {}, - "source": [ - "## Propagate all objects to t and discard too low and high ones" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3106348d", - "metadata": {}, - "outputs": [], - "source": [ - "objects = []\n", - "count_too_low = 0\n", - "count_too_high = 0\n", - "for planet in planet_list:\n", - " try:\n", - " pos,v = planet.eph(starting_t)\n", - " \n", - " # convert to km and numpy\n", - " pos = np.asarray(pos) / 1000.0 \n", - " v = np.asarray(v) / 1000.0\n", - " sma,_,_,_,_,_ = pk.ic2par(pos * 1000,v *1000,mu=pk.MU_EARTH)\n", - " \n", - " altitude = np.linalg.norm(pos)\n", - " if altitude < lower_cutoff_in_km:\n", - " count_too_low += 1\n", - " continue\n", - " if sma / 1000. > higher_cutoff_in_km or altitude > higher_cutoff_in_km:\n", - " count_too_high += 1\n", - " continue\n", - " \n", - " objects.append((pos,v))\n", - " except RuntimeError as e:\n", - " print(e, \" propagating \",planet.name)\n", - " \n", - "print(\"Successfully propagated \",len(objects),\" objects.\")\n", - "print(count_too_low,\" had a too small altitude\")\n", - "print(count_too_high,\" had a too high altitude\")" - ] - }, - { - "cell_type": "markdown", - "id": "5c77e85b", - "metadata": {}, - "source": [ - "## Plot and store results" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8aae75ff", - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure(figsize=(6,6),dpi=100)\n", - "ax = plt.axes(projection='3d');\n", - "\n", - "positions = np.array([pos for pos,_ in objects])\n", - "velocities = np.array([v for _,v in objects])\n", - "\n", - "ax.scatter(positions[:,0],positions[:,1],positions[:,2],\".\",alpha=0.25)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "dd111eb3", - "metadata": {}, - "outputs": [], - "source": [ - "np.savetxt(\"../data/pos.csv\",positions,delimiter=\",\")\n", - "np.savetxt(\"../data/v.csv\",velocities,delimiter=\",\")" - ] - }, - { - "cell_type": "markdown", - "id": "d23e9dba", - "metadata": {}, - "source": [ - "# Propagate test set by some time" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cf6eb7d5", - "metadata": {}, - "outputs": [], - "source": [ - "# Load test data\n", - "pos = np.loadtxt(\"../data/pos.csv\",delimiter=\",\")\n", - "v = np.loadtxt(\"../data/v.csv\",delimiter=\",\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d16a0034", - "metadata": {}, - "outputs": [], - "source": [ - "# Propagate by t seconds\n", - "t = 10\n", - "objects = []\n", - "t_end = pk.epoch(starting_t.mjd + t * pk.SEC2DAY,\"mjd\")\n", - "for pos_i,v_i in zip(pos,v):\n", - " try:\n", - " p = pk.planet.keplerian(starting_t,pos_i * 1000.0,v_i * 1000.0,pk.MU_EARTH,1.,1.,1.)\n", - " r,v = p.eph(t_end)\n", - " \n", - " objects.append((np.array(r) / 1000., np.array(v) / 1000.))\n", - " \n", - " except RuntimeError as e:\n", - " print(e, \" propagating \",p.name)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6db2f74d", - "metadata": {}, - "outputs": [], - "source": [ - "# unpack\n", - "positions = np.array([r for r,_ in objects])\n", - "velocities = np.array([v_i for _,v_i in objects])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ce544999", - "metadata": {}, - "outputs": [], - "source": [ - "#look at them\n", - "pos" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6d303ca5", - "metadata": {}, - "outputs": [], - "source": [ - "positions" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4c65416f", - "metadata": {}, - "outputs": [], - "source": [ - "# Save\n", - "np.savetxt(\"../data/pos_test_10s.csv\",positions,delimiter=\",\")\n", - "np.savetxt(\"../data/v_test_10s.csv\",velocities,delimiter=\",\")" - ] - } - ], - "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.7" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/notebooks/Input Dataset Processing.ipynb b/notebooks/Input Dataset Processing.ipynb new file mode 100644 index 00000000..81bdf5d9 --- /dev/null +++ b/notebooks/Input Dataset Processing.ipynb @@ -0,0 +1,448 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "id": "b442b78f", + "metadata": {}, + "source": [ + "# Input Data Processing\n", + "\n", + "This notebook serves the following purposes:\n", + "\n", + "1. Read XML data from space-track.org\n", + "2. Read the CSV SATCAT catalog from celestrak\n", + "3. Based on above data estimate mass and radius (characteristic length), get activity state\n", + "4. Propagate all the satellites to the same point in time\n", + "5. Investigate data and clean-up unwanted data. Then store.\n", + "6. Old code used for generating test data.\n", + "\n", + "## Input files\n", + "- CSV data from the [CelesTrak SATCAT catalog](https://celestrak.com/pub/satcat.csv) following this [format](https://celestrak.com/satcat/satcat-format.php)\n", + "- XML 3LE data from the [Space-Track.org catalog](https://www.space-track.org/) following their [format](https://www.space-track.org/documentation#/tle)\n", + "\n", + "## Output files\n", + "\n", + "- Satellite data in CSV format with data on Satellite ID, Position, Velocity, Mass, Radius (characteristic length) and Activity State" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a48e2c9", + "metadata": {}, + "outputs": [], + "source": [ + "### Imports\n", + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "# Append main folder\n", + "import sys\n", + "sys.path.append(\"../\")\n", + "\n", + "import pykep as pk\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "from tqdm import tqdm \n", + "\n", + "starting_t = pk.epoch_from_string('2022-01-01 00:00:00.000')\n", + "lower_cutoff_in_km = 6371 + 200 # Earth radius + ...\n", + "higher_cutoff_in_km = 6371 + 2000" + ] + }, + { + "cell_type": "markdown", + "id": "de523087", + "metadata": {}, + "source": [ + "## 1. Read XML data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf98e917", + "metadata": {}, + "outputs": [], + "source": [ + "import xml.etree.ElementTree as ET\n", + "def parse_xml(file):\n", + " \"\"\"Parse spacetrack xml. \n", + " This function is inspired by https://github.com/brandon-rhodes/python-sgp4/blob/master/sgp4/omm.py , MIT Licensed\n", + " \"\"\"\n", + " root = ET.parse(file).getroot()\n", + " for segment in root.findall('.//segment'):\n", + " metadata = segment.find('metadata')\n", + " data = segment.find('data')\n", + " meanElements = data.find('meanElements')\n", + " tleParameters = data.find('tleParameters')\n", + " userParameters = data.find('userDefinedParameters')\n", + " fields = {}\n", + " for element in metadata, meanElements, tleParameters, userParameters:\n", + " fields.update((field.items()[0][1], field.text) if len(field.items()) > 0 else (field.tag, field.text) for field in element)\n", + " yield fields" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6d36cf3", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a generator to iterate over the data\n", + "fields = parse_xml(\"../data/spacetrack.xml\")" + ] + }, + { + "cell_type": "markdown", + "id": "afd3af8a", + "metadata": {}, + "source": [ + "## 2. Read SATCAT data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f50d2589", + "metadata": {}, + "outputs": [], + "source": [ + "# Load all the xml data from space-track\n", + "satellites = []\n", + "while True:\n", + " try:\n", + " satellites.append(next(fields))\n", + " if len(satellites) % 5000 == 0:\n", + " print(\"Loaded \",len(satellites), \"sats...\")\n", + " except StopIteration:\n", + " print(\"Loaded \",len(satellites), \"sats...Done\")\n", + " break" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb53b6d1", + "metadata": {}, + "outputs": [], + "source": [ + "# Read satcat data from celestrak\n", + "satcat = pd.read_csv(\"../data/satcat.csv\")\n", + "satcat" + ] + }, + { + "cell_type": "markdown", + "id": "ad908ccc", + "metadata": {}, + "source": [ + "## 3. Compute mass and radius (characteristic length) and get status\n", + "\n", + "This follows the formulas from \n", + "Nicholas L Johnson, Paula H Krisko, J-C Liou, and Phillip D Anz-Meador.\n", + "Nasa’s new breakup model of evolve 4.0. Advances in Space Research, 28(9):1377–\n", + "1384, 2001.\n", + "\n", + "According to space-track , RCS small, medium and large are, respectively < 0.1 , 0.1 < RCS < 1.0 and 1.0 < RCS. For simplicity using above formula we convert this to 15cm, 55cm, 200cm\n", + "\n", + "We get activity status from the celestrak data following https://celestrak.com/satcat/status.php" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a79d2bf", + "metadata": {}, + "outputs": [], + "source": [ + "sats_with_info = []\n", + "for sat in tqdm(satellites):\n", + " \n", + " satcat_sat = satcat[satcat[\"OBJECT_ID\"] == sat[\"OBJECT_ID\"]]\n", + " \n", + " # Skip decayed ones or ones not in celestrak\n", + " if len(satcat_sat) == 0 or satcat_sat[\"OPS_STATUS_CODE\"].values == \"D\":\n", + " continue\n", + " \n", + " # Determine L_C\n", + " if not np.isnan(satcat_sat[\"RCS\"].values):\n", + " sat[\"RADIUS\"] = 2. * np.sqrt(float(satcat_sat[\"RCS\"].values) / np.pi)\n", + " else:\n", + " if sat[\"RCS_SIZE\"] == \"SMALL\":\n", + " sat[\"RADIUS\"] = 0.15\n", + " elif sat[\"RCS_SIZE\"] == \"MEDIUM\":\n", + " sat[\"RADIUS\"] = 0.55\n", + " elif sat[\"RCS_SIZE\"] == \"LARGE\":\n", + " sat[\"RADIUS\"] = 2.0\n", + " else:\n", + " # skip if no info was found\n", + " continue\n", + " \n", + " # Determine Mass\n", + " if sat[\"RADIUS\"] > 0.01:\n", + " sat[\"MASS\"] = 4 / 3 * np.pi *(sat[\"L_C\"] / 2)**3 * 92.937 * sat[\"L_C\"]**(-0.74)\n", + " else:\n", + " sat[\"MASS\"] = 4 / 3 * np.pi *(sat[\"L_C\"] / 2)**3 * 2698.9\n", + " \n", + " \n", + " # Determine if active satellite\n", + " if satcat_sat[\"OPS_STATUS_CODE\"].values in [\"+\",\"P\",\"B\",\"S\",\"X\"]:\n", + " sat[\"TYPE\"] = \"evasive\"\n", + " else:\n", + " sat[\"TYPE\"] = \"passive\"\n", + " \n", + " # Add planet\n", + " t0 = pk.epoch_from_string(sat[\"EPOCH\"].replace(\"T\",\" \"))\n", + " elements = [float(sat[\"SEMIMAJOR_AXIS\"]) * 1000.,\n", + " float(sat[\"ECCENTRICITY\"]),\n", + " float(sat[\"INCLINATION\"]) * pk.DEG2RAD,\n", + " float(sat[\"RA_OF_ASC_NODE\"]) * pk.DEG2RAD,\n", + " float(sat[\"ARG_OF_PERICENTER\"]) * pk.DEG2RAD,\n", + " float(sat[\"MEAN_ANOMALY\"]) * pk.DEG2RAD,\n", + " ]\n", + " planet = pk.planet.keplerian(t0,elements,pk.MU_EARTH,6.67430e-11*sat[\"MASS\"],sat[\"RADIUS\"] / 2,sat[\"RADIUS\"] / 2)\n", + " sat[\"PLANET\"] = planet\n", + " \n", + " sats_with_info.append(sat)\n", + " \n", + "print(\"Now we have a total of \",len(sats_with_info), \"sats.\")" + ] + }, + { + "cell_type": "markdown", + "id": "8cd629fd", + "metadata": {}, + "source": [ + "### Plot some examples" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f3f3a16", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(6,6),dpi=100)\n", + "ax = plt.axes(projection='3d');\n", + "for i in range (10):\n", + " pk.orbit_plots.plot_planet(sats_with_info[i][\"PLANET\"],axes=ax)" + ] + }, + { + "cell_type": "markdown", + "id": "dbd01df1", + "metadata": {}, + "source": [ + "## 4. Propagate all objects to t and discard too low and high ones" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3106348d", + "metadata": {}, + "outputs": [], + "source": [ + "objects = []\n", + "count_too_low = 0\n", + "count_too_high = 0\n", + "\n", + "for sat in sats_with_info:\n", + " try:\n", + " planet = sat[\"PLANET\"]\n", + " pos,v = planet.eph(starting_t)\n", + " \n", + " # convert to km and numpy\n", + " pos = np.asarray(pos) / 1000.0 \n", + " v = np.asarray(v) / 1000.0\n", + " sma,_,_,_,_,_ = pk.ic2par(pos * 1000,v *1000,mu=pk.MU_EARTH)\n", + " \n", + " altitude = np.linalg.norm(pos)\n", + " if altitude < lower_cutoff_in_km:\n", + " count_too_low += 1\n", + " continue\n", + " if sma / 1000. > higher_cutoff_in_km or altitude > higher_cutoff_in_km:\n", + " count_too_high += 1\n", + " continue\n", + " \n", + " objects.append({\"ID\": sat[\"OBJECT_NAME\"],\n", + " \"R\": tuple(pos),\n", + " \"V\": tuple(v),\n", + " \"M\": sat[\"MASS\"],\n", + " \"RADIUS\": sat[\"RADIUS\"],\n", + " \"TYPE\": sat[\"TYPE\"]\n", + " })\n", + " except RuntimeError as e:\n", + " print(e, \" propagating \",planet.name)\n", + " \n", + "print(\"Successfully propagated \",len(objects),\" objects.\")\n", + "print(count_too_low,\" had a too small altitude\")\n", + "print(count_too_high,\" had a too high altitude\")" + ] + }, + { + "cell_type": "markdown", + "id": "5c77e85b", + "metadata": {}, + "source": [ + "## 5. Plot, clean up and store results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8aae75ff", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(6,6),dpi=100)\n", + "ax = plt.axes(projection='3d');\n", + "\n", + "positions = np.array([obj[\"R\"] for obj in objects])\n", + "velocities = np.array([obj[\"V\"] for obj in objects])\n", + "\n", + "ax.scatter(positions[:,0],positions[:,1],positions[:,2],\".\",alpha=0.25)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd111eb3", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert to pandas dataframe and drop ISS and any duplicate entries.\n", + "df = pd.DataFrame(objects)\n", + "df = df.drop(np.argmax(df[\"ID\"] == \"ISS (ZARYA)\"))\n", + "df = df.drop(df[df.ID.str.startswith('STARLINK')].index)\n", + "df = df.drop(df[df.ID.str.startswith('ONEWEB')].index)\n", + "df = df.drop_duplicates(subset=['R'])\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "055b4b52", + "metadata": {}, + "outputs": [], + "source": [ + "# Write to csv\n", + "df.to_csv(\"../data/initial_population.csv\")" + ] + }, + { + "cell_type": "markdown", + "id": "d23e9dba", + "metadata": {}, + "source": [ + "# 6. (deprecated) Propagate test set by some time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf6eb7d5", + "metadata": {}, + "outputs": [], + "source": [ + "# Load test data\n", + "pos = np.loadtxt(\"../data/pos.csv\",delimiter=\",\")\n", + "v = np.loadtxt(\"../data/v.csv\",delimiter=\",\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d16a0034", + "metadata": {}, + "outputs": [], + "source": [ + "# Propagate by t seconds\n", + "t = 10\n", + "objects = []\n", + "t_end = pk.epoch(starting_t.mjd + t * pk.SEC2DAY,\"mjd\")\n", + "for pos_i,v_i in zip(pos,v):\n", + " try:\n", + " p = pk.planet.keplerian(starting_t,pos_i * 1000.0,v_i * 1000.0,pk.MU_EARTH,1.,1.,1.)\n", + " r,v = p.eph(t_end)\n", + " \n", + " objects.append((np.array(r) / 1000., np.array(v) / 1000.))\n", + " \n", + " except RuntimeError as e:\n", + " print(e, \" propagating \",p.name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6db2f74d", + "metadata": {}, + "outputs": [], + "source": [ + "# unpack\n", + "positions = np.array([r for r,_ in objects])\n", + "velocities = np.array([v_i for _,v_i in objects])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce544999", + "metadata": {}, + "outputs": [], + "source": [ + "#look at them\n", + "pos" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d303ca5", + "metadata": {}, + "outputs": [], + "source": [ + "positions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c65416f", + "metadata": {}, + "outputs": [], + "source": [ + "# Save\n", + "np.savetxt(\"../data/pos_test_10s.csv\",positions,delimiter=\",\")\n", + "np.savetxt(\"../data/v_test_10s.csv\",velocities,delimiter=\",\")" + ] + } + ], + "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.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/Output Analysis.ipynb b/notebooks/Output Analysis.ipynb deleted file mode 100644 index 5cec63ec..00000000 --- a/notebooks/Output Analysis.ipynb +++ /dev/null @@ -1,202 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "id": "b5fe1f2e", - "metadata": {}, - "outputs": [], - "source": [ - "### Imports\n", - "%load_ext autoreload\n", - "%autoreload 2\n", - "\n", - "# Append main folder\n", - "import sys\n", - "sys.path.append(\"../\")\n", - "\n", - "from glob import glob\n", - "\n", - "from tqdm import tqdm\n", - "import pykep as pk\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", - "\n", - "dt = 0.01\n", - "starting_t = pk.epoch_from_string('2022-01-01 00:00:00.000')\n", - "iterations = 100\n", - "end_t = pk.epoch(starting_t.mjd2000 + iterations * dt * pk.SEC2DAY)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ba281ee2", - "metadata": {}, - "outputs": [], - "source": [ - "# Load test data\n", - "init_r = np.loadtxt(\"../data/pos.csv\",delimiter=\",\")\n", - "init_v = np.loadtxt(\"../data/v.csv\",delimiter=\",\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3c6adf1e", - "metadata": {}, - "outputs": [], - "source": [ - "files = glob(\"../build/*.vtu\")\n", - "file = \"../build/output_\" + str(iterations) + \".vtu\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bf9e8f72", - "metadata": {}, - "outputs": [], - "source": [ - "from vtk import vtkXMLUnstructuredGridReader\n", - "from vtk.util import numpy_support as VN\n", - "\n", - "reader = vtkXMLUnstructuredGridReader()\n", - "reader.SetFileName(file)\n", - "reader.Update()\n", - "\n", - "data = reader.GetOutput()\n", - "\n", - "v = VN.vtk_to_numpy(data.GetPointData().GetArray('velocity')).astype(\"double\")\n", - "r = VN.vtk_to_numpy(data.GetPoints().GetData()).astype(\"double\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0cf78a93", - "metadata": {}, - "outputs": [], - "source": [ - "init_v" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "93beb687", - "metadata": {}, - "outputs": [], - "source": [ - "v / 1000" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "20fe73f7", - "metadata": {}, - "outputs": [], - "source": [ - "r" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ce842a55", - "metadata": {}, - "outputs": [], - "source": [ - "errs = [[],[],[],[],[],[]]\n", - "\n", - "for v_p,r_p, init_v_p,init_r_p in tqdm(zip(v,r,init_v,init_r)):\n", - " a,e,i,W,w,E = pk.ic2par(r_p,v_p, pk.MU_EARTH)\n", - " oa,oe,oi,oW,ow,oE = pk.ic2par(init_r_p * 1000.0,init_v_p* 1000.0, pk.MU_EARTH)\n", - " errs[0].append(abs(a-oa))\n", - " errs[1].append(abs(e-oe))\n", - " errs[2].append(abs(i-oi))\n", - " errs[3].append(abs(W-oW))\n", - " errs[4].append(abs(w-ow))\n", - " errs[5].append(abs(E-oE))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "06e35f59", - "metadata": {}, - "outputs": [], - "source": [ - "for err,name in zip(errs, [\"a\",\"e\",\"i\",\"W\",\"w\",\"E\"]):\n", - " print(f\"{name} Error - Mean. {np.mean(err):<10.6f} Min. {np.min(err):<8.6f} Max. {np.max(err):<8.6f}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9106a5ec", - "metadata": {}, - "outputs": [], - "source": [ - "for err,name in zip(errs, [\"a\",\"e\",\"i\",\"W\",\"w\",\"E\"]):\n", - " df = pd.DataFrame(data=err)\n", - " df.plot.hist(title=name,bins=100,density=True, cumulative=True)\n", - " plt.xlim(0,100)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c62f0d3f", - "metadata": {}, - "outputs": [], - "source": [ - "idx = 1000\n", - "p_init = pk.planet.keplerian(starting_t, init_r[idx] * 1000.0,init_v[idx] * 1000.0,pk.MU_EARTH,1.,1.,1.)\n", - "p = pk.planet.keplerian(end_t,r[idx], v[idx],pk.MU_EARTH,1.,1.,1.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9686a1c4", - "metadata": {}, - "outputs": [], - "source": [ - "ax = pk.orbit_plots.plot_planet(p_init)\n", - "pk.orbit_plots.plot_planet(p,axes=ax)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7c75bcf5", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "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.7" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/notebooks/Results Analysis.ipynb b/notebooks/Results Analysis.ipynb index 61aef617..71626df1 100644 --- a/notebooks/Results Analysis.ipynb +++ b/notebooks/Results Analysis.ipynb @@ -9,7 +9,21 @@ "\n", "This notebook serves analysing LADDS run results, the following packages need to be installed\n", "\n", - "`conda install tqdm pykep pandas scikit-learn h5py -c conda-forge`" + "`conda install tqdm pykep pandas scikit-learn h5py -c conda-forge`\n", + "\n", + "In general, it will:\n", + "\n", + "1. Load the HDF data generated by a simulation and extract position, velocities, conjunctions.\n", + "2. Create a bunch of plots related to conjunctions\n", + "3. Investigate orbital dynamics (elements, distances between particles) over the simulation\n", + "\n", + "## Input \n", + "\n", + "This notebook uses the input from a simulation run in HDF5 format. For more details refer to the [docs](https://github.com/esa/LADDS#hdf5)\n", + "\n", + "## Output\n", + "\n", + "N/A at the moment but the generated plots may be saved as desired." ] }, { @@ -49,7 +63,7 @@ "id": "ca306edd", "metadata": {}, "source": [ - "### Load the run data\n", + "## 1. Load the run data\n", "Adapt the path below to the h5 file you want to load." ] }, @@ -225,7 +239,7 @@ "id": "16ae0da2", "metadata": {}, "source": [ - "## Let the plotting begin.\n", + "## 2. Conjunction Plots\n", "We start with counting the number of conjunctions for different thresholds over the simulation time and then threshold vs total # of conjunctions." ] }, @@ -389,6 +403,7 @@ "id": "a8d7fac8", "metadata": {}, "source": [ + "## 3. Dynamics\n", "### Now let's orbital elements of all particles over the simulation as a sanity check" ] }, @@ -539,14 +554,6 @@ "ax.set_xlim([0,max_dist])\n", "plt.tight_layout()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "27a29b9e", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": {