Skip to content

Instantly share code, notes, and snippets.

# alterbaron/landingcalc.m

Created May 24, 2013 14:05
Show Gist options
• Save alterbaron/5643760 to your computer and use it in GitHub Desktop.
Useful calculation routines and landing phase angle calculator for Kerbal Space Program
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
to join this conversation on GitHub. Already have an account? Sign in to comment