How to plot/visualize (elliptical) orbits of celestial objects in heliocentric reference frame?

50 次查看（过去 30 天）
Giuseppe
Giuseppe 2022-2-3
回答： Karan Singh 2024-1-30
Hi guys!
I'm solving numerical problems of celestial mechanics in Fortran, the data pre-processing and post-processing is performed in Matlab instead.
My problem is to find "design of low-energy trajectories to Near Earth Objects". So it is an interplanetary transfer to small bodies (asteroids and comets) of the Solar system.
Firstly I need to plot an heliocentric reference frame, the Earth orbit and some asteroids orbit, then I want to plot the trajectory of the spacecraft.
Can you help me?
  4 个评论
显示 2更早的评论
James Tursa
James Tursa 2022-2-3
Does your current code generate the 3D position data of the bodies? I.e., do you have the data but are just asking about plotting in 3D in MATLAB? Or are you also asking how to generate the 3D data given orbital elements?
Giuseppe
Giuseppe 2022-2-3
Hi @James Tursa, for now I have the six orbital parameters of 55 asteroids, and I would like to plot (3D and 2D) their orbits whitin an heliocentric reference frame.
Something like this (from NASA JPL solar system dynamics website):

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Karan Singh
Karan Singh 2024-1-30
Ran in:
Hi Giuseppe,
As you wish to plot the 2D and the 3D, the task you've described involves visualizing the orbits of 55 asteroids within a heliocentric (Sun-centered) reference frame. I hope you have gathered the data required for you to plot the same. I have take some points to plot as-
  • Semi-major axis (a): The average distance from the asteroid to the Sun, defining the size of the orbit.
  • Eccentricity (e): A measure of how much the orbit deviates from a perfect circle.
  • Inclination (i): The tilt of the asteroid's orbit relative to the plane of the Solar System (ecliptic plane).
  • Longitude of the ascending node (Ω): The angle from a reference direction (usually the direction of the vernal equinox) to the point where the asteroid crosses the ecliptic plane going northward.
  • Argument of periapsis (ω): The angle from the ascending node to the orbit's closest point to the Sun (periapsis) within the orbital plane.
  • True anomaly (ν): The angle between the direction of periapsis and the current position of the asteroid on its orbit, measured at the Sun.
To plot the orbits, I have converted the orbital parameters into Cartesian coordinates (x, y, z) for each point along the orbit.
  • In a 3D space with the plot3 function. This will allow you to visualize the inclinations and positions of the orbits relative to one another.
  • 2D projection of these orbits, typically onto the ecliptic plane (x-y plane), using the plot function. To see the relative paths from a top-down perspective
Here's an example MATLAB script that demonstrates how to plot the orbits of asteroids in 3D and 2D, using the orbital parameters. This script assumes that the orbital parameters are stored in a matrix where each row represents an asteroid and the columns correspond to the six orbital parameters.
% Define constants
AU = 149597870.7; % Astronomical unit in km
mu = 1.32712440018e11; % Standard gravitational parameter of the Sun in km^3/s^2
% Example asteroid orbital parameters: [a, e, i, Ω, ω, ν]
% Each row corresponds to an asteroid
asteroid_orbital_params = [
2.767 * AU, 0.0756, 10.593, 80.399, 73.773, 0;
% ... Add rows for each of the 55 asteroids
];
% Number of asteroids
num_asteroids = size(asteroid_orbital_params, 1);
% Prepare figure for plotting
figure;
hold on;
grid on;
xlabel('X (km)');
ylabel('Y (km)');
zlabel('Z (km)');
title('Heliocentric Reference Frame with Asteroid Orbits');
% Plot Sun
plot3(0, 0, 0, 'yo', 'MarkerSize', 10, 'MarkerFaceColor', 'y');
% Plot each asteroid's orbit
for idx = 1:num_asteroids
% Extract orbital parameters for the asteroid
a = asteroid_orbital_params(idx, 1);
e = asteroid_orbital_params(idx, 2);
i = deg2rad(asteroid_orbital_params(idx, 3)); % Convert to radians
Omega = deg2rad(asteroid_orbital_params(idx, 4)); % Convert to radians
omega = deg2rad(asteroid_orbital_params(idx, 5)); % Convert to radians
nu = deg2rad(asteroid_orbital_params(idx, 6)); % Convert to radians
% Number of points to plot the orbit
num_points = 360;
% True anomaly range
nu_range = linspace(0, 2*pi, num_points);
% Preallocate arrays for Cartesian coordinates
x = zeros(1, num_points);
y = zeros(1, num_points);
z = zeros(1, num_points);
% Calculate position for each point on the orbit
for j = 1:num_points
% Current true anomaly
nu_current = nu_range(j);
% Radius at current true anomaly
r = (a * (1 - e^2)) / (1 + e * cos(nu_current));
% Position in the orbital plane
x_orbital = r * cos(nu_current);
y_orbital = r * sin(nu_current);
% Convert to 3D Cartesian coordinates
x(j) = (cos(omega) * cos(Omega) - sin(omega) * sin(Omega) * cos(i)) * x_orbital + (-sin(omega) * cos(Omega) - cos(omega) * sin(Omega) * cos(i)) * y_orbital;
y(j) = (cos(omega) * sin(Omega) + sin(omega) * cos(Omega) * cos(i)) * x_orbital + (-sin(omega) * sin(Omega) + cos(omega) * cos(Omega) * cos(i)) * y_orbital;
z(j) = (sin(omega) * sin(i)) * x_orbital + (cos(omega) * sin(i)) * y_orbital;
end
% Plot the orbit in 3D
plot3(x, y, z);
end
% Adjust plot
axis equal;
view(3); % Set to 3D view
% Plot 2D projection (XY plane)
figure;
hold on;
grid on;
xlabel('X (km)');
ylabel('Y (km)');
title('Heliocentric Reference Frame with Asteroid Orbits (2D Projection)');

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Earth and Planetary Science 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by