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

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?

[David Hill]

Did you search the file exchange? Seems like there are a few files that might be helpful.

[Giuseppe]

Hi @David Hill yes, I did. But for now I found only Earth-centered functions.

[James Tursa]

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]

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):

Answer 1

[Karan Singh]

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)');