forked from MichaelMurry/interplanetary
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gravityflyby.m
36 lines (30 loc) · 1.26 KB
/
gravityflyby.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
function [ psi, hp, rp ] = gravityflyby( target_planet, V_inf_in, ...
V_inf_out )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description: This function computes the turn angle and altitude of
% closest approach to target planet given excess velocity in and out of
% the target planet.
% Inputs:
% target_planet = string variable of the target planet
% V_inf_in = excess velocity into the target planet [km/s^2]
% V_inf_out = excess velocity out of the target planet [km/s^2]
% Planet
% Outputs
% rp = radius of peraispse [km]
% hp = altitude of periaspse [km]
% psi = turn angle [rad]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define the gravitational parameter of target planet [km^3/s^2]
mu = gparam(target_planet);
% Define target planet' equitorial radius [km]
r_0 = eqradius(target_planet);
% Compute turn angle [rad]
psi = acos(dot(V_inf_in{1,1}, V_inf_out{1,1})/...
(norm(V_inf_in{1,1})*norm(V_inf_out{1,1})));
% Compute radius of periapse [km]
rp = mu/norm(V_inf_in{1,1})^2*(1/(cos((pi - psi)/2)) - 1);
% Compute alittude of periapse [km]
hp = rp - r_0;
end