Compute the trajectory of a particle under the gravity of a massive central body

trajectory.m

% In this script we find the orbital path a body under gravity.
% The mass of this particle is insignifcant compared to the mass
% of the central body.

% Clear Existing Values
clear all;
clf;

%##------------- Start Configuration ----------------##%

%% Initial Conditions
% Initial Position
x = [0, 0.909];
% Initial Velocity
v = [2.5, 0];

%% Environment
M  = 0.826; % Central mass
G  = 4.5;   % Gravitational Constant

%% Simulation and plot
% Number of revolutions
n = 1;
% Number of data points (plot)
N = 1000;

%##------------- End Configuration ----------------##%

% x = initial postion (cartesian)
% v = initial velocity (cartesian)

% We use polar coordinates as it is more convienient for central forces
% r = distance to body from centre
% theta = angle to x axis
% u = 1/r
% the force is of type:
% F = K/r^2;

% Initial r
r_0  = norm(x);
% Initial u
u_0  = 1/r_0;
% Initial theta
theta_0 = acos(dot(x,[1,0])/norm(x));
% Angular momentum per unit mass h = r^2*(dtheta/dt) = norm(cross(v,x)) = constant
% but we take a vector h
h    = cross([x, 0], [v, 0]);
% rdot = dr/dt = norm(projection of v on x)
rdot = dot(v,x)/norm(x);
% The Standard Gravitational Constant
mu = G*M;

% Constant in the equation of motion
A = mu / norm(h)^2;

% The Equation of motion is:
% d^2u/dtheta^2 = A/u^(2+k) - u
% The gereral Solution is:
% u = A + M*cos(theta + phi), M and phi constants

%% Finding the particular solution for x and v.
% the eccentricity vector
e = cross([v, 0], h) / mu - [x, 0] / norm(x);
% Evaluation of the constants
M = A * norm(e);
phi = acos((u_0 - A)/M) - theta_0;
% Use particular solution to compute data points
theta = linspace (theta_0, theta_0 + 2*pi*n, N)';
u = A + M .* cos(theta + phi);
r = 1./u;

%% Set up data polt
subplot(1,2,1);
hold on;
plot (theta, u);
plot (theta, r);
% axis limits
p_max = max(max(u(:,1),max(r)));
p_min = min(min(u(:,1),min(r)));
margin = 0.5;
% finishing touches
axis([min(theta), max(theta), p_min - margin, p_max + margin] );
title("r and u");
xlabel("theta")
legend("u", "r");
grid("on");
hold off;

%% set up trejectory plot
subplot(1,2,2);
polar(theta, r, "--b");
% axis limits
lim    = max(r) + margin;
no     = 2; % number of grid spacings
tticks = 0:30:(360-30);
set (gca, "rtick", lim/no:lim/no:lim, "ttick", tticks);
axis([-lim, lim, -lim, lim], "equal");
title("Trejectory");
grid("on");

please refer to this repo for a detailed explanation