From 6086ee7714f80dd3a8cd5c40e2fab5c009660e8b Mon Sep 17 00:00:00 2001 From: Mingde Yin Date: Sun, 24 Jul 2022 22:52:40 -0400 Subject: [PATCH 1/3] Demo --- python_code/plotter.py | 71 ++++++++++++++++++++++++++++++++++------ python_code/supernova.py | 40 +++++++++++++++++++--- 2 files changed, 97 insertions(+), 14 deletions(-) diff --git a/python_code/plotter.py b/python_code/plotter.py index bf138db..7823071 100644 --- a/python_code/plotter.py +++ b/python_code/plotter.py @@ -2,6 +2,30 @@ import numpy as np +def set_axes_equal(ax: plt.Axes): + """Set 3D plot axes to equal scale. + + Make axes of 3D plot have equal scale so that spheres appear as + spheres and cubes as cubes. Required since `ax.axis('equal')` + and `ax.set_aspect('equal')` don't work on 3D. + """ + limits = np.array([ + ax.get_xlim3d(), + ax.get_ylim3d(), + ax.get_zlim3d(), + ]) + origin = np.mean(limits, axis=1) + radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0])) + _set_axes_radius(ax, origin, radius) + + +def _set_axes_radius(ax, origin, radius): + x, y, z = origin + ax.set_xlim3d([x - radius, x + radius]) + ax.set_ylim3d([y - radius, y + radius]) + ax.set_zlim3d([z - radius, z + radius]) + + def plot_trajectory(*filenames: str): fig = plt.figure(figsize=(18, 7), tight_layout=True) @@ -43,7 +67,7 @@ def plot_from_array(t, arr): axes = [] for i in range(3): - axes.append(fig.add_subplot(231 + i)) + axes.append(fig.add_subplot(131 + i)) axes[i].add_patch(earth[i]) axes[i].set_aspect("equal") axes[i].grid() @@ -51,25 +75,52 @@ def plot_from_array(t, arr): x = arr[:, 0] y = arr[:, 1] z = arr[:, 2] - axes[0].plot(x, y) + axes[0].plot(x, y, color="green") axes[0].set_ylabel("y position (m)") axes[0].set_xlabel("x position (m)") - axes[1].plot(y, z) + axes[1].plot(y, z, color="green") axes[1].set_ylabel("z position (m)") axes[1].set_xlabel("y position (m)") - axes[2].plot(x, z) + axes[2].plot(x, z, color="green") axes[2].set_ylabel("z position (m)") axes[2].set_xlabel("x position (m)") axes[0].set_title("xy") axes[1].set_title("yz") axes[2].set_title("xz") - axes.append(fig.add_subplot(212)) - # Plot altitude above each - altitude = np.sqrt(x**2 + y**2 + z**2) - 6.378e6 - axes[3].plot(t/86400, altitude) - axes[3].set_ylabel("Altitude (m)") - axes[3].set_xlabel("Mission Time (d)") + plt.suptitle("Spacecraft Position, Inertial Coordinates") + + plt.show() + + +def plot_3d_from_array(t, arr): + fig = plt.figure(figsize=(6, 6), tight_layout=True) + ax = plt.axes(projection='3d') + + x = arr[:, 0] + y = arr[:, 1] + z = arr[:, 2] + + ax.plot(x, y, z, color="green") + + # Set the labels + ax.set_xlabel('x (m)') + ax.set_ylabel('y (m)') + ax.set_zlabel('z (m)') + + u = np.linspace(0, 2*np.pi, 100) + v = np.linspace(0, np.pi, 100) + x = np.outer(np.cos(u), np.sin(v)) # np.outer() -> outer vector product + y = np.outer(np.sin(u), np.sin(v)) + z = np.outer(np.ones(np.size(u)), np.cos(v)) + ax.plot_surface(x*6378e3, y*6378e3, z*6378e3) + + # Set the title + ax.set_title("Spacecraft Position, Inertial Coordinates") + + ax.set_box_aspect([1, 1, 1]) # IMPORTANT - this is the new, key line + # ax.set_proj_type('ortho') # OPTIONAL - default is perspective (shown in image above) + set_axes_equal(ax) # IMPORTANT - this is also required plt.show() diff --git a/python_code/supernova.py b/python_code/supernova.py index a2208c8..62c0cab 100644 --- a/python_code/supernova.py +++ b/python_code/supernova.py @@ -1,6 +1,8 @@ from ctypes import CDLL, POINTER, c_double, c_int, c_char_p, Structure import numpy as np -from plotter import plot_from_array +from plotter import plot_from_array, plot_3d_from_array +from celest.satellite import Time, Coordinate, Satellite +from celest.encounter import GroundPosition, windows ''' Code Responsible for Loading C Library @@ -47,10 +49,12 @@ def wrap_func(lib, funcname, restype, argtypes): tSpan = (c_double * len(py_tSpan))(*py_tSpan) state: POINTER(c_double) = state(6378e3+270e3, 0, 97 * np.pi/180, 0, 0, 0) - solution: POINTER(AdaptiveSolution) = orbit("RK810".encode("utf-8"), "advanced".encode("utf-8"), tSpan, state, 1e-6) - # py_y0 = [-4749231.102296294, -4975106.82469687, 0.0, -719.6538503589323, 686.980716081442, 7561.282735263496] - # y0 = (c_double * len(py_y0))(*py_y0) + py_y0 = [-4749231.102296294, -4975106.82469687, 0.0, - + 719.6538503589323, 686.980716081442, 7561.282735263496] + y0 = (c_double * len(py_y0))(*py_y0) + solution: POINTER(AdaptiveSolution) = orbit("RK810".encode( + "utf-8"), "simplified".encode("utf-8"), tSpan, y0, 1e-6) n = solution[0].n # Dereference pointer print(f"Steps taken: {n}") @@ -58,3 +62,31 @@ def wrap_func(lib, funcname, restype, argtypes): y = np.ctypeslib.as_array(solution[0].y, shape=(n, 6)) initial_state = np.ctypeslib.as_array(state, shape=(6,)) plot_from_array(t, y) + plot_3d_from_array(t, y) + + # Celest stuff + julian = Time(t/86400) + position = Coordinate(y[:, :3], "gcrs", t/86400) + + toronto = GroundPosition(latitude=43.6532, longitude=-79.3832) + southern_manitoba = GroundPosition(latitude=53.7167, longitude=-98.8167) + + # Get ITRS data. + itrs_x, itrs_y, itrs_z = position.itrs() + + # plot_from_array(t, np.vstack((itrs_x, itrs_y, itrs_z)).T) + + alt, az = position.horizontal(southern_manitoba) + + satellite = Satellite(position=y[:, :3], frame='gcrs', julian=t/86400, + offset=2414900) + + satellite.itrs() + + # Generate ground location windows. + IMG_windows = windows.generate( + satellite=satellite, location=southern_manitoba, enc="image", ang=30, lighting=1) + GL_windows = windows.generate( + satellite=satellite, location=toronto, enc="data_link", ang=10, lighting=0) + + print("hello") From c8bec4db751d31dc1676468eebbd6d717521e2f0 Mon Sep 17 00:00:00 2001 From: Mingde Yin Date: Thu, 22 Jun 2023 22:51:22 -0400 Subject: [PATCH 2/3] overhaul into a viable python package --- .gitignore | 162 ++++++++++++++++++ {python_code => old_python_code}/nn_api.py | 0 python_code/reference_propagator.py | 42 ----- python_code/supernova.py | 92 ---------- setup.py | 31 ++++ {new_python => standalone_python}/JGM3.py | 0 .../RK1210_weights.py | 0 .../RK4_weights.py | 0 {new_python => standalone_python}/__init__.py | 0 .../ode_solver.py | 0 .../propagate_orbit.py | 0 .../propagator_solution.py | 0 {new_python => standalone_python}/rk1210.txt | 0 .../simplified_physics.py | 0 .../stability_region.py | 0 {new_python => standalone_python}/test_RK.py | 0 supernova/__init__.py | 0 supernova/_builder.py | 37 ++++ supernova/api.py | 82 +++++++++ supernova/demos/__init__.py | 0 supernova/demos/integrated_mission.py | 69 ++++++++ {python_code => supernova}/plotter.py | 0 .../test_against_poliastro.py | 89 ++++++++-- 23 files changed, 452 insertions(+), 152 deletions(-) rename {python_code => old_python_code}/nn_api.py (100%) delete mode 100644 python_code/reference_propagator.py delete mode 100644 python_code/supernova.py create mode 100644 setup.py rename {new_python => standalone_python}/JGM3.py (100%) rename {new_python => standalone_python}/RK1210_weights.py (100%) rename {new_python => standalone_python}/RK4_weights.py (100%) rename {new_python => standalone_python}/__init__.py (100%) rename {new_python => standalone_python}/ode_solver.py (100%) rename {new_python => standalone_python}/propagate_orbit.py (100%) rename {new_python => standalone_python}/propagator_solution.py (100%) rename {new_python => standalone_python}/rk1210.txt (100%) rename {new_python => standalone_python}/simplified_physics.py (100%) rename {new_python => standalone_python}/stability_region.py (100%) rename {new_python => standalone_python}/test_RK.py (100%) create mode 100644 supernova/__init__.py create mode 100644 supernova/_builder.py create mode 100644 supernova/api.py create mode 100644 supernova/demos/__init__.py create mode 100644 supernova/demos/integrated_mission.py rename {python_code => supernova}/plotter.py (100%) rename python_code/comparison.py => tests/test_against_poliastro.py (55%) diff --git a/.gitignore b/.gitignore index 738c089..067278d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,165 @@ *.exe *.csv *.so + + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/#use-with-ide +.pdm.toml + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ \ No newline at end of file diff --git a/python_code/nn_api.py b/old_python_code/nn_api.py similarity index 100% rename from python_code/nn_api.py rename to old_python_code/nn_api.py diff --git a/python_code/reference_propagator.py b/python_code/reference_propagator.py deleted file mode 100644 index 10807d9..0000000 --- a/python_code/reference_propagator.py +++ /dev/null @@ -1,42 +0,0 @@ -from numba import njit -from poliastro.twobody.propagation import func_twobody, propagate, cowell -from poliastro.core.perturbations import J2_perturbation -from astropy import units as u -from poliastro.bodies import Earth -import numpy as np - - -#@title Aerodynamic Characteristics - -#@markdown C_D: Drag coefficeint -- dim ensionless - -#@markdown FRONTAL_AREA: Frontal area of S/C in m^2 - -#@markdown MASS: Spacecraft mass in kg - -@njit -def a_d(t0, state, k, J2, R): - return J2_perturbation(t0, state, k, J2, R) - -# Force model -def f(t0, state, k): - du_kep = func_twobody(t0, state, k) - # standard twobody - - ax, ay, az = a_d( - t0, - state, - k, - R=Earth.R.to(u.km).value, - J2=Earth.J2.value, - ) - du_ad = np.array([0, 0, 0, ax, ay, az]) - - return du_kep + du_ad - - -def prop(orb, times, fname): - rr = propagate(orb, times * u.s, method=cowell, f=f) - arr = rr.xyz.value.T * 1000 - times = np.reshape(times, (times.shape[0], 1)) - np.savetxt(fname, np.concatenate((times, arr), axis=-1), delimiter=",") diff --git a/python_code/supernova.py b/python_code/supernova.py deleted file mode 100644 index 62c0cab..0000000 --- a/python_code/supernova.py +++ /dev/null @@ -1,92 +0,0 @@ -from ctypes import CDLL, POINTER, c_double, c_int, c_char_p, Structure -import numpy as np -from plotter import plot_from_array, plot_3d_from_array -from celest.satellite import Time, Coordinate, Satellite -from celest.encounter import GroundPosition, windows - -''' -Code Responsible for Loading C Library -''' -supernova = CDLL(("./python_code/supernova.so")) -# Load library - - -class AdaptiveSolution(Structure): - _fields_ = [("t", POINTER(c_double)), - ("y", POINTER(c_double)), - ("n", c_int)] - ''' - Structure: - t: Times array of shape (n,) - y: position/velocity array of shape (n, 6) with units [m, s] - n: number of steps used - ''' - - -def wrap_func(lib, funcname, restype, argtypes): - ''' Referenced from - https://dbader.org/blog/python-ctypes-tutorial-part-2 - Taken from ESC190 - ''' - func = lib.__getattr__(funcname) - func.restype = restype - func.argtypes = argtypes - return func - - -orbit = wrap_func(supernova, "orbitPTR", POINTER(AdaptiveSolution), - [c_char_p, c_char_p, POINTER(c_double), POINTER(c_double), - c_double]) -state = wrap_func(supernova, "stateFromKeplerian", POINTER(c_double), - [c_double, c_double, c_double, c_double, c_double, c_double]) -# Define functions from supernova C in Python - - -if __name__ == "__main__": - days = 50 - - py_tSpan = [0, 86400 * days] - tSpan = (c_double * len(py_tSpan))(*py_tSpan) - - state: POINTER(c_double) = state(6378e3+270e3, 0, 97 * np.pi/180, 0, 0, 0) - - py_y0 = [-4749231.102296294, -4975106.82469687, 0.0, - - 719.6538503589323, 686.980716081442, 7561.282735263496] - y0 = (c_double * len(py_y0))(*py_y0) - solution: POINTER(AdaptiveSolution) = orbit("RK810".encode( - "utf-8"), "simplified".encode("utf-8"), tSpan, y0, 1e-6) - - n = solution[0].n # Dereference pointer - print(f"Steps taken: {n}") - t = np.ctypeslib.as_array(solution[0].t, shape=(n,)) - y = np.ctypeslib.as_array(solution[0].y, shape=(n, 6)) - initial_state = np.ctypeslib.as_array(state, shape=(6,)) - plot_from_array(t, y) - plot_3d_from_array(t, y) - - # Celest stuff - julian = Time(t/86400) - position = Coordinate(y[:, :3], "gcrs", t/86400) - - toronto = GroundPosition(latitude=43.6532, longitude=-79.3832) - southern_manitoba = GroundPosition(latitude=53.7167, longitude=-98.8167) - - # Get ITRS data. - itrs_x, itrs_y, itrs_z = position.itrs() - - # plot_from_array(t, np.vstack((itrs_x, itrs_y, itrs_z)).T) - - alt, az = position.horizontal(southern_manitoba) - - satellite = Satellite(position=y[:, :3], frame='gcrs', julian=t/86400, - offset=2414900) - - satellite.itrs() - - # Generate ground location windows. - IMG_windows = windows.generate( - satellite=satellite, location=southern_manitoba, enc="image", ang=30, lighting=1) - GL_windows = windows.generate( - satellite=satellite, location=toronto, enc="data_link", ang=10, lighting=0) - - print("hello") diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..6480f5e --- /dev/null +++ b/setup.py @@ -0,0 +1,31 @@ +from setuptools import setup + +with open("README.md", "r", encoding="utf-8") as fh: + long_description = fh.read() + +setup( + name="supernova-cffi", + version="1.0.0", + author="Mingde Yin", + description="Orbit propagation and analysis tool.", + long_description=long_description, + long_description_content_type="text/markdown", + url="https://github.com/spacesys-finch/supernova/", + license="LGPLv3", + packages=["supernova"], + setup_requires=["cffi>=1.15.0"], + cffi_modules=["supernova/_builder.py:ffibuilder"], + install_requires=["cffi>=1.15.0", "numpy", "matplotlib"], + extras_require={"test": ["pytest", "poliastro"], "demo": ["celest"]}, + options={"build_ext": {"inplace": True}}, + classifiers=[ + "Development Status :: 3 - Alpha", + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Programming Language :: C", + "Topic :: Scientific/Engineering :: Physics", + "Intended Audience :: Science/Research", + ], + python_requires=">=3.8", +) diff --git a/new_python/JGM3.py b/standalone_python/JGM3.py similarity index 100% rename from new_python/JGM3.py rename to standalone_python/JGM3.py diff --git a/new_python/RK1210_weights.py b/standalone_python/RK1210_weights.py similarity index 100% rename from new_python/RK1210_weights.py rename to standalone_python/RK1210_weights.py diff --git a/new_python/RK4_weights.py b/standalone_python/RK4_weights.py similarity index 100% rename from new_python/RK4_weights.py rename to standalone_python/RK4_weights.py diff --git a/new_python/__init__.py b/standalone_python/__init__.py similarity index 100% rename from new_python/__init__.py rename to standalone_python/__init__.py diff --git a/new_python/ode_solver.py b/standalone_python/ode_solver.py similarity index 100% rename from new_python/ode_solver.py rename to standalone_python/ode_solver.py diff --git a/new_python/propagate_orbit.py b/standalone_python/propagate_orbit.py similarity index 100% rename from new_python/propagate_orbit.py rename to standalone_python/propagate_orbit.py diff --git a/new_python/propagator_solution.py b/standalone_python/propagator_solution.py similarity index 100% rename from new_python/propagator_solution.py rename to standalone_python/propagator_solution.py diff --git a/new_python/rk1210.txt b/standalone_python/rk1210.txt similarity index 100% rename from new_python/rk1210.txt rename to standalone_python/rk1210.txt diff --git a/new_python/simplified_physics.py b/standalone_python/simplified_physics.py similarity index 100% rename from new_python/simplified_physics.py rename to standalone_python/simplified_physics.py diff --git a/new_python/stability_region.py b/standalone_python/stability_region.py similarity index 100% rename from new_python/stability_region.py rename to standalone_python/stability_region.py diff --git a/new_python/test_RK.py b/standalone_python/test_RK.py similarity index 100% rename from new_python/test_RK.py rename to standalone_python/test_RK.py diff --git a/supernova/__init__.py b/supernova/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/supernova/_builder.py b/supernova/_builder.py new file mode 100644 index 0000000..8bfd4ae --- /dev/null +++ b/supernova/_builder.py @@ -0,0 +1,37 @@ +from cffi import FFI + +ffibuilder = FFI() + +# two public functions: orbit and stateFromKeplerian +ffibuilder.cdef( + """ +typedef struct AdaptiveSolution { + double* t; // timesteps + double (*y)[6]; // function value + int n; // number of steps taken +} solution; + +solution* orbitPTR(char solver[], char model[], double tSpan[], double y0[], double ATOL); +double* stateFromKeplerian(double sma, double ecc, double inc, double raan, double aop, double m); +void free(void* ptr); // free memory allocated by C code +""" +) + + +ffibuilder.set_source( + "supernova.backend", + """ + #include "constants.h" + #include "forces.h" + #include "gravity.h" + #include "numericaldiff.h" + #include "orbit.h" + #include "solvers.h" + #include "vecmath.h" + """, + include_dirs=["./c_code"], + extra_compile_args=["-std=c99", "-O3", "-march=native", "-ffast-math", "-lm"], +) + +if __name__ == "__main__": + ffibuilder.compile(verbose=False) diff --git a/supernova/api.py b/supernova/api.py new file mode 100644 index 0000000..4f45479 --- /dev/null +++ b/supernova/api.py @@ -0,0 +1,82 @@ +from __future__ import annotations + +import numpy as np +import numpy.typing as npt + +from supernova.backend import ffi +from supernova.backend.lib import orbitPTR, stateFromKeplerian, free + + +def propagate_orbit( + solver: str, model: str, t_span: tuple[float], y0: tuple[float], atol: float +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Propagate orbit using given solver and model. + + Parameters + ---------- + solver : str + Name of the solver to use, either "RK810" or "RK1012" + RK810 is a 10th order, 17 stage Runge-Kutta method + RK1012 is a 12th order, 25 stage Runge-Kutta method + model : str + Name of the model to use, either "simplified" or "advanced" + simplified: only J2 perturbation is considered + advanced: JGM-3 20th order spherical harmonic model is used + t_span : tuple[float] + Time span to propagate the orbit over, in seconds + y0 : tuple[float] + Initial state vector, in meters and meters per second (ECI) + (x, y, z, vx, vy, vz) + atol : float + Absolute tolerance for the solver + + Returns + ------- + t : npt.NDArray[np.float64] + Time steps, in seconds + y : npt.NDArray[np.float64] + State vector, in meters and meters per second (ECI) + """ + solution = orbitPTR(solver.encode("utf-8"), model.encode("utf-8"), t_span, y0, atol) + + t_buf = ffi.buffer(solution.t, solution.n * 8) + t = np.frombuffer(t_buf, dtype=np.float64).copy() + t.shape = (solution.n,) + + y_buf = ffi.buffer(solution.y, solution.n * 8 * 6) + y = np.frombuffer(y_buf, dtype=np.float64).copy() + y.shape = (solution.n, 6) + + free(solution.t) + free(solution.y) + + return t, y + + +def state_from_keplerian( + sma: float, ecc: float, inc: float, raan: float, aop: float, m: float +) -> tuple[float]: + """Calculate state vector from Keplerian elements. + + Parameters + ---------- + sma : float + Semi-major axis, in meters + ecc : float + Eccentricity + inc : float + Inclination, in radians + raan : float + Right ascension of the ascending node, in radians + aop : float + Argument of periapsis, in radians + m : float + Mean anomaly, in radians + """ + d_ptr = stateFromKeplerian(sma, ecc, inc, raan, aop, m) + + d_buf = ffi.buffer(d_ptr, 6 * 8) + d = np.frombuffer(d_buf, dtype=np.float64).copy() + + free(d_ptr) + return tuple(*d) diff --git a/supernova/demos/__init__.py b/supernova/demos/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/supernova/demos/integrated_mission.py b/supernova/demos/integrated_mission.py new file mode 100644 index 0000000..31ed938 --- /dev/null +++ b/supernova/demos/integrated_mission.py @@ -0,0 +1,69 @@ +from supernova.api import propagate_orbit +from supernova.plotter import plot_from_array, plot_3d_from_array +from celest.satellite import Satellite +from celest.time import Time +from celest.coordinates import Coordinate, GroundLocation, GCRS, ITRS +from celest.encounter import generate_vtws +from celest import units as u + +from datetime import datetime + + +def jd2000_to_datetime(jd2000: u.Quantity) -> datetime: + return Time(jd2000).datetime + + +if __name__ == "__main__": + days = 50 + + t_span = [0, 86400 * days] + + y0 = [ + -4749231.102296294, + -4975106.82469687, + 0.0, + -719.6538503589323, + 686.980716081442, + 7561.282735263496, + ] + + t, y = propagate_orbit("RK810", "simplified", t_span, y0, 1e-6) + + print(f"Steps taken: {len(t)}") + # plot_from_array(t, y) + # plot_3d_from_array(t, y) + + # Celest stuff + julian = Time(t / 86400, offset=2414900) + + position = GCRS(julian.julian.data, y[:, 0], y[:, 1], y[:, 2], u.m) + velocity = GCRS(julian.julian.data, y[:, 3], y[:, 4], y[:, 5], u.m / u.s) + + toronto = GroundLocation( + latitude=43.6532, + longitude=-79.3832, + height=76, + angular_unit=u.deg, + length_unit=u.m, + ) + + # pos_aa = Coordinate(position).convert_to("AzEl", toronto) + # pos_itrs = Coordinate(position).convert_to(ITRS) + + # plot_from_array(t, pos_itrs.data) + + satellite = Satellite(position=position, velocity=velocity) + + # # Generate ground location windows. + downlinking_windows = generate_vtws( + satellite=satellite, location=toronto, vis_threshold=10 + ) + + for window in downlinking_windows: + rise_time = window.rise_time + + print(jd2000_to_datetime(rise_time).strftime("%Y-%m-%d %H:%M:%S")) + + set_time = window.set_time + + print(jd2000_to_datetime(set_time).strftime("%Y-%m-%d %H:%M:%S")) diff --git a/python_code/plotter.py b/supernova/plotter.py similarity index 100% rename from python_code/plotter.py rename to supernova/plotter.py diff --git a/python_code/comparison.py b/tests/test_against_poliastro.py similarity index 55% rename from python_code/comparison.py rename to tests/test_against_poliastro.py index b90f943..c4ea6a5 100644 --- a/python_code/comparison.py +++ b/tests/test_against_poliastro.py @@ -6,15 +6,61 @@ from astropy import units as u from astropy.time import Time import time -import plotter +import supernova.plotter as plotter from reference_propagator import prop from ctypes import c_double +from numba import njit +from poliastro.twobody.propagation import func_twobody, propagate, cowell +from poliastro.core.perturbations import J2_perturbation +from astropy import units as u +from poliastro.bodies import Earth +import numpy as np + + +# @title Aerodynamic Characteristics + +# @markdown C_D: Drag coefficeint -- dim ensionless + +# @markdown FRONTAL_AREA: Frontal area of S/C in m^2 + +# @markdown MASS: Spacecraft mass in kg + + +@njit +def a_d(t0, state, k, J2, R): + return J2_perturbation(t0, state, k, J2, R) + + +# Force model +def f(t0, state, k): + du_kep = func_twobody(t0, state, k) + # standard twobody + + ax, ay, az = a_d( + t0, + state, + k, + R=Earth.R.to(u.km).value, + J2=Earth.J2.value, + ) + du_ad = np.array([0, 0, 0, ax, ay, az]) + + return du_kep + du_ad + + +def prop(orb, times, fname): + rr = propagate(orb, times * u.s, method=cowell, f=f) + arr = rr.xyz.value.T * 1000 + times = np.reshape(times, (times.shape[0], 1)) + np.savetxt(fname, np.concatenate((times, arr), axis=-1), delimiter=",") + + def get_orbit(): - ''' + """ Get SSO from poliastro - ''' + """ SMA = 6903 * u.km LTAN = 22.5 * u.hourangle Eccentricity = 0.003621614 * u.one @@ -24,16 +70,19 @@ def get_orbit(): epoch_time = "00:00" # Store time as UTC - epoch = Time(f"{epoch_date} {epoch_time}", format='iso', scale="utc") + epoch = Time(f"{epoch_date} {epoch_time}", format="iso", scale="utc") print("Determining intial orbit from parameters...") # Define orbital parameters, along with units - orb = Orbit.heliosynchronous(attractor=Earth, - a=SMA, ecc=Eccentricity, - ltan=LTAN, - argp=argument_of_perigee, - nu=total_anomaly, - epoch=epoch) + orb = Orbit.heliosynchronous( + attractor=Earth, + a=SMA, + ecc=Eccentricity, + ltan=LTAN, + argp=argument_of_perigee, + nu=total_anomaly, + epoch=epoch, + ) # orb = Orbit.from_classical(attractor=Earth, # a=SMA, ecc = Eccentricity, @@ -55,14 +104,14 @@ def compare_propagators(orbit, days: int): t_start = time.perf_counter() orb = get_orbit() - t_orb = time.perf_counter()-t_start + t_orb = time.perf_counter() - t_start print(f"parameter fetch completion in {t_orb} seconds.") - py_y0 = [*(orb.r.value*1000), *(orb.v.value*1000)] + py_y0 = [*(orb.r.value * 1000), *(orb.v.value * 1000)] y0 = (c_double * len(py_y0))(*py_y0) t_start = time.perf_counter() orbit("RK810".encode("utf-8"), tSpan, y0, "FINCH.csv".encode("utf-8"), 1e-4) - t_C = time.perf_counter()-t_start + t_C = time.perf_counter() - t_start print(f"Supernova completion in {t_C} seconds.") arr = np.loadtxt("FINCH.csv", delimiter=",") # Load times @@ -70,15 +119,19 @@ def compare_propagators(orbit, days: int): t_start = time.perf_counter() prop(orb, ref_times, "poliastro.csv") - t_Poli = time.perf_counter()-t_start + t_Poli = time.perf_counter() - t_start print(f"Poliastro propagator completion in {t_Poli} seconds.") print(f"Supernova is {t_Poli/t_C}x faster than Poliastro") # Error analysis arrPoli = np.loadtxt("poliastro.csv", delimiter=",") diff = arr[:, 1:4] - arrPoli[:, 1:] - print(f"RMS error between Supernova and Poliastro is {np.sqrt(np.average(diff**2))} metres across {len(ref_times)} steps.") - print(f"Max error between Supernova and Poliastro is {np.max(np.abs(diff))} metres across {len(ref_times)} steps.") + print( + f"RMS error between Supernova and Poliastro is {np.sqrt(np.average(diff**2))} metres across {len(ref_times)} steps." + ) + print( + f"Max error between Supernova and Poliastro is {np.max(np.abs(diff))} metres across {len(ref_times)} steps." + ) plotter.plot_trajectory("FINCH.csv", "poliastro.csv") plotter.plot_error(np.abs(diff)) @@ -90,9 +143,9 @@ def propagate_orbit(orbit, days: int): tSpan = (c_double * len(py_tSpan))(*py_tSpan) orb = get_orbit() - py_y0 = [*(orb.r.value*1000), *(orb.v.value*1000)] + py_y0 = [*(orb.r.value * 1000), *(orb.v.value * 1000)] y0 = (c_double * len(py_y0))(*py_y0) orbit("RK810".encode("utf-8"), tSpan, y0, "FINCH.csv".encode("utf-8"), 1e-12) - plotter.plot_trajectory("FINCH.csv") \ No newline at end of file + plotter.plot_trajectory("FINCH.csv") From 23bc61345ef8630adc675ff826c72e4133c1a23f Mon Sep 17 00:00:00 2001 From: Mingde Yin Date: Thu, 13 Jul 2023 22:20:57 -0400 Subject: [PATCH 3/3] refactor mission demo into another repo --- setup.py | 6 +- supernova/demos/__init__.py | 0 supernova/demos/integrated_mission.py | 69 --------------------- supernova/plotter.py | 88 +++++++++------------------ 4 files changed, 32 insertions(+), 131 deletions(-) delete mode 100644 supernova/demos/__init__.py delete mode 100644 supernova/demos/integrated_mission.py diff --git a/setup.py b/setup.py index 6480f5e..d471c99 100644 --- a/setup.py +++ b/setup.py @@ -15,8 +15,8 @@ packages=["supernova"], setup_requires=["cffi>=1.15.0"], cffi_modules=["supernova/_builder.py:ffibuilder"], - install_requires=["cffi>=1.15.0", "numpy", "matplotlib"], - extras_require={"test": ["pytest", "poliastro"], "demo": ["celest"]}, + install_requires=["cffi>=1.15.0", "numpy>=1.25", "matplotlib>=3.7"], + extras_require={"test": ["pytest", "poliastro"]}, options={"build_ext": {"inplace": True}}, classifiers=[ "Development Status :: 3 - Alpha", @@ -27,5 +27,5 @@ "Topic :: Scientific/Engineering :: Physics", "Intended Audience :: Science/Research", ], - python_requires=">=3.8", + python_requires=">=3.9", ) diff --git a/supernova/demos/__init__.py b/supernova/demos/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/supernova/demos/integrated_mission.py b/supernova/demos/integrated_mission.py deleted file mode 100644 index 31ed938..0000000 --- a/supernova/demos/integrated_mission.py +++ /dev/null @@ -1,69 +0,0 @@ -from supernova.api import propagate_orbit -from supernova.plotter import plot_from_array, plot_3d_from_array -from celest.satellite import Satellite -from celest.time import Time -from celest.coordinates import Coordinate, GroundLocation, GCRS, ITRS -from celest.encounter import generate_vtws -from celest import units as u - -from datetime import datetime - - -def jd2000_to_datetime(jd2000: u.Quantity) -> datetime: - return Time(jd2000).datetime - - -if __name__ == "__main__": - days = 50 - - t_span = [0, 86400 * days] - - y0 = [ - -4749231.102296294, - -4975106.82469687, - 0.0, - -719.6538503589323, - 686.980716081442, - 7561.282735263496, - ] - - t, y = propagate_orbit("RK810", "simplified", t_span, y0, 1e-6) - - print(f"Steps taken: {len(t)}") - # plot_from_array(t, y) - # plot_3d_from_array(t, y) - - # Celest stuff - julian = Time(t / 86400, offset=2414900) - - position = GCRS(julian.julian.data, y[:, 0], y[:, 1], y[:, 2], u.m) - velocity = GCRS(julian.julian.data, y[:, 3], y[:, 4], y[:, 5], u.m / u.s) - - toronto = GroundLocation( - latitude=43.6532, - longitude=-79.3832, - height=76, - angular_unit=u.deg, - length_unit=u.m, - ) - - # pos_aa = Coordinate(position).convert_to("AzEl", toronto) - # pos_itrs = Coordinate(position).convert_to(ITRS) - - # plot_from_array(t, pos_itrs.data) - - satellite = Satellite(position=position, velocity=velocity) - - # # Generate ground location windows. - downlinking_windows = generate_vtws( - satellite=satellite, location=toronto, vis_threshold=10 - ) - - for window in downlinking_windows: - rise_time = window.rise_time - - print(jd2000_to_datetime(rise_time).strftime("%Y-%m-%d %H:%M:%S")) - - set_time = window.set_time - - print(jd2000_to_datetime(set_time).strftime("%Y-%m-%d %H:%M:%S")) diff --git a/supernova/plotter.py b/supernova/plotter.py index 7823071..77f4668 100644 --- a/supernova/plotter.py +++ b/supernova/plotter.py @@ -1,70 +1,42 @@ +from __future__ import annotations + import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + import numpy as np -def set_axes_equal(ax: plt.Axes): +def set_axes_equal(ax: Axes3D): """Set 3D plot axes to equal scale. Make axes of 3D plot have equal scale so that spheres appear as spheres and cubes as cubes. Required since `ax.axis('equal')` and `ax.set_aspect('equal')` don't work on 3D. """ - limits = np.array([ - ax.get_xlim3d(), - ax.get_ylim3d(), - ax.get_zlim3d(), - ]) + limits = np.array( + [ + ax.get_xlim3d(), + ax.get_ylim3d(), + ax.get_zlim3d(), + ] + ) origin = np.mean(limits, axis=1) radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0])) _set_axes_radius(ax, origin, radius) -def _set_axes_radius(ax, origin, radius): +def _set_axes_radius(ax: Axes3D, origin: np.ndarray, radius: np.ndarray): x, y, z = origin ax.set_xlim3d([x - radius, x + radius]) ax.set_ylim3d([y - radius, y + radius]) ax.set_zlim3d([z - radius, z + radius]) -def plot_trajectory(*filenames: str): - fig = plt.figure(figsize=(18, 7), tight_layout=True) - - earth = [plt.Circle((0, 0), 6.378e6, color='lightblue') for i in range(3)] - axes = [] - - for i in range(3): - axes.append(fig.add_subplot(131 + i)) - axes[i].add_patch(earth[i]) - axes[i].set_aspect("equal") - axes[i].grid() - - for filename in filenames: - arr = np.loadtxt(filename, delimiter=",") - x = arr[:, 1] - y = arr[:, 2] - z = arr[:, 3] - axes[0].plot(x, y, label=filename) - axes[0].set_ylabel("y position (m)") - axes[0].set_xlabel("x position (m)") - axes[1].plot(y, z, label=filename) - axes[1].set_ylabel("z position (m)") - axes[1].set_xlabel("y position (m)") - axes[2].plot(x, z, label=filename) - axes[2].set_ylabel("z position (m)") - axes[2].set_xlabel("x position (m)") - axes[0].set_title("xy") - axes[1].set_title("yz") - axes[2].set_title("xz") - - plt.legend() - plt.show() - - -def plot_from_array(t, arr): +def plot_from_array(t: np.ndarray, arr: np.ndarray) -> None: fig = plt.figure(figsize=(12, 5), tight_layout=True) - earth = [plt.Circle((0, 0), 6.378e6, color='lightblue') for i in range(3)] - axes = [] + earth = [plt.Circle((0, 0), 6.378e6, color="lightblue") for i in range(3)] + axes: list[plt.Axes] = [] for i in range(3): axes.append(fig.add_subplot(131 + i)) @@ -93,27 +65,29 @@ def plot_from_array(t, arr): plt.show() -def plot_3d_from_array(t, arr): +def plot_3d_from_array( + t: np.ndarray, arr: np.ndarray, array_desampling_factor: int = 1, earth_n: int = 20 +): fig = plt.figure(figsize=(6, 6), tight_layout=True) - ax = plt.axes(projection='3d') + ax = plt.axes(projection="3d") - x = arr[:, 0] - y = arr[:, 1] - z = arr[:, 2] + x = arr[::array_desampling_factor, 0] + y = arr[::array_desampling_factor, 1] + z = arr[::array_desampling_factor, 2] ax.plot(x, y, z, color="green") # Set the labels - ax.set_xlabel('x (m)') - ax.set_ylabel('y (m)') - ax.set_zlabel('z (m)') + ax.set_xlabel("x (m)") + ax.set_ylabel("y (m)") + ax.set_zlabel("z (m)") - u = np.linspace(0, 2*np.pi, 100) - v = np.linspace(0, np.pi, 100) + u = np.linspace(0, 2 * np.pi, earth_n) + v = np.linspace(0, np.pi, earth_n) x = np.outer(np.cos(u), np.sin(v)) # np.outer() -> outer vector product y = np.outer(np.sin(u), np.sin(v)) z = np.outer(np.ones(np.size(u)), np.cos(v)) - ax.plot_surface(x*6378e3, y*6378e3, z*6378e3) + ax.plot_surface(x * 6378e3, y * 6378e3, z * 6378e3, alpha=0.5, color="lightblue") # Set the title ax.set_title("Spacecraft Position, Inertial Coordinates") @@ -144,7 +118,3 @@ def plot_error(diff: np.ndarray): axes[2].set_title("z error") plt.show() - - -if __name__ == "__main__": - plot_trajectory("c_code/simplified.csv", "c_code/advanced.csv")