Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Useful calculation routines and landing phase angle calculator for Kerbal Space Program
function [phi_out] = landingcalc( orbit_alt, landing_pe )
% Returns the landing phase angle, where the ship starts in a circular
% orbit and drops its PE down to landing_pe.
% Set to true to log a bunch of details.
details = false;
% Kerbin planetary data from Wiki
dt = 1;
m = 1.0;
A = 1;
mu = 3.5316e12;
Ratm = 600000+69077.5;
Rmin = 600000;
P0 = 1;
H0 = 5000;
d = 0.2;
SOI = 84159286;
omega_rot = 2*pi/21600; % Rotational frequency (rad/s)
% Measure from the planet's centre
r0 = orbit_alt + Rmin;
rpe = landing_pe + Rmin;
if (rpe > Ratm)
phi_out = 0;
fprintf('rpe > Ratm, no atmosphere encounter.\n');
return;
end
% circ. orbit
v0 = sqrt(mu/r0);
% Step 1: Find the delta-v that when applied retrograde drops you into the
% orbit with the target landing PE.
% I worked this out on paper, too long for this margin to contain, etcetera
% . . .
v1_mag = sqrt(2*mu*(rpe/r0)*(r0-rpe)/(r0^2-rpe^2));
delta_v = abs(v1_mag-v0);
% We now have our orbital parameters for the actual landing orbit.
% (v,r,rpe). I convert these into more useful parameters here.
r = [r0 0];
v = [0 v1_mag];
[ ep ec a hmag rpe rap ] = get_orbit_params( r, v, mu );
% Calculate vcontact and rcontact where the orbit intersects the
% atmosphere.
cos_theta_contact = (1/ec)*(a*(1-ec^2)/Ratm - 1);
theta_contact = acos(cos_theta_contact);
vcontact_mag = sqrt(2*(ep + mu/Ratm));
theta_1 = asin(hmag/(Ratm*vcontact_mag));
vcontact = vcontact_mag * [-cos(theta_1+theta_contact) -sin(theta_1+theta_contact)];
rcontact = Ratm * [cos(theta_contact) sin(theta_contact)];
% Find eccentric anomaly
E = atan2(sqrt(1-ec^2)*sin(theta_contact), ec+cos(theta_contact));
% Find transit time until atmospheric contact.
% Man, "eccentric anomaly" is about archaic as it gets.
% Note that M0 = pi, and the expression is effectively negated to give a
% positive time.
Ttransit = sqrt(a^3/mu)*(pi-E+ec*sin(E));
% Simulate craft in the atmosphere.
F = in_atmo_force( Rmin, P0, H0, d, mu, m, A );
[rfinal vfinal Tatm] = integrate_path( F, m, rcontact, vcontact, dt, Rmin, Ratm );
% Angle through which the planet rotates and angle through which the craft
% travels through the atmosphere.
theta_atm = 180/pi*acos((rcontact*rfinal')/(norm(rcontact)*norm(rfinal)));
theta_rot = 180/pi * omega_rot * (Tatm + Ttransit);
% NOTE: All values measured from planet centre
% Log a bunch of useful orbital data.
if (details)
fprintf('Landing Maneuver Data\n------------------------------\n');
fprintf('Landing maneuver delta-v: %f m/s\n', delta_v);
fprintf('\n');
fprintf('Initial Orbit Data\n------------------------------\n');
fprintf('=> Sp. orbital energy: %f\n', ep);
fprintf('=> Eccentricity: %f\n', ec);
fprintf('=> Sp. angular momentum: %f\n', hmag);
fprintf('=> Semi-major axis: %f\n', a);
fprintf('=> Eccentric anomaly at atmosphere contact: %f\n', E);
fprintf('=> PE: %f m\n', rpe);
fprintf('=> AP: %f m\n', (1+ec)*a);
fprintf('=> Transit time: %f\n', Ttransit);
fprintf('=> Atmos. contact angle: %f degrees\n', theta_contact*180/pi);
fprintf('=> Atmospheric entry velocity: %f m/s\n', norm(vcontact));
fprintf('\n');
fprintf('Atmospheric Maneuver Data\n------------------------------\n');
fprintf('=> Time in atmosphere: %f s\n', Tatm);
fprintf('=> Phase angle of manoeuver: %f degrees\n', ...
theta_atm);
if (norm(rfinal) > Rmin)
fprintf('=> Escaped atmosphere.\n');
fprintf('=> Exit velocity: %f m/s\n', norm(vfinal));
fprintf('=> Final Orbit Data\n------------------------------\n');
fprintf('=> => Sp. orbital energy: %f\n', ep);
fprintf('=> => Eccentricity: %f\n', ec);
fprintf('=> => Sp. angular momentum: %f\n', hmag);
fprintf('=> => Semi-major axis: %f\n', a);
fprintf('=> => PE: %f\n', rpe);
fprintf('=> => AP: %f\n', (1+ec)*a);
else
fprintf('=> Surface impact!\n');
fprintf('=> Impact velocity: %f m/s\n', norm(vfinal));
end
fprintf('\n');
fprintf('Planet rotation angle: %f degrees\n', theta_rot);
end
if (norm(rfinal) > Rmin) % Escaped atmosphere
phi_out = 0; % No landing!
return;
end
% Phase angle to put PE ahead of target for landing.
phi_out = theta_rot+(theta_contact*180/pi-theta_atm);
end
function [ F ] = in_atmo_force( R0, P0, H0, d, mu, m, A )
% Returns a function which gives the total vector force in atmosphere
% as a function of position and velocity.
Kp = 1.2230948554874*0.008;
% Total force = drag + gravity
F = @(r,v) -0.5*Kp*P0*exp((R0-norm(r))/H0)*norm(v)*d*m*A*v - m*(mu/norm(r).^3)*r;
end
function [r v t] = integrate_path( F, m, r0, v0, dt, Rmin, Ratm )
% Integrate the equation of motion to find the path in atmosphere.
% Terminate on collision with surface (sea level) or atmosphere escape.
t = 0;
r = r0;
v = v0;
a = @(rin, vin) F(rin, vin)./m;
firstrun = true;
while firstrun || (norm(r) <= Ratm && norm(r) >= Rmin)
r_old = r;
v_old = v;
a_t = a(r_old, v_old);
r = r_old + v_old*dt + 0.5*a_t*dt.^2;
v_est = v_old + 0.5*dt*(a_t + a(r,v_old+dt*a_t));
v = v_old + 0.5*(a_t + a(r,v_est))*dt;
t = t + dt;
firstrun = false;
end
end
function [ ep ec a hmag rpe rap ] = get_orbit_params( r, v, mu )
% Get useful orbit parameters from position vector, velocity vector, and
% gravitational parameter.
% sp. orbital energy
ep = dot(v, v)/2 - mu/norm(r);
%sp angular momentum
r3d = [r(1) r(2) 0];
v3d = [v(1) v(2) 0];
hmag = cross(r3d, v3d);
hmag = hmag(3);
%eccentricity
ec = sqrt(1+2*ep*hmag.^2./(mu.^2));
%semi-major axis
a = -mu/(2*ep);
%periapse distance
rpe = -a*(ec-1);
% Apoapse distance (valid iff ec < 0)
rap = (1+ec)*a;
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment