How to plot/visualize (elliptical) orbits of celestial objects in heliocentric reference frame?
66 views (last 30 days)
Show older comments
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 Comments
James Tursa
on 3 Feb 2022
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?
Accepted Answer
Karan Singh
on 30 Jan 2024
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 Comments
More Answers (0)
See Also
Categories
Find more on Earth and Planetary Science in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!