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