-
Notifications
You must be signed in to change notification settings - Fork 3
/
int_orb_eq.m
67 lines (58 loc) · 1.64 KB
/
int_orb_eq.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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
function [rx_vect, ry_vect, rz_vect, vx_vect, vy_vect, vz_vect] = int_orb_eq(kep,mu,time)
% int_orb_eq.m
%
% PROTOTYPE:
% [rx_vect, ry_vect, rz_vect, vx_vect, vy_vect, vz_vect] = int_orb_eq(kep,mu,time)
%
% DESCRIPTION:
% This function implements the integration of the dynamics orbit equations
%
% INPUT:
% a[1] Semimajor axis
% e[1] Eccentricity
% i[1] Inclination
% OMG[1] Right ascension
% omg[1] Pericenter anomaly
% theta[1] True anomaly
% mu
%
% OUTPUT:
% rx_vect[3] Vector of all x components of position vectors
% ry_vect[3] Vector of all y components of position vectors
% rz_vect[3] Vector of all z components of position vectors
% vx_vect[3] Vector of all x components of velocity vectors
% vy_vect[3] Vector of all y components of velocity vectors
% vz_vect[3] Vector of all z components of velocity vectors
%
%
% AUTHOR:
% Alfonso Collogrosso
%
% Using kep2car function to find position and velocity at initial condition
[r,v]=kep2car(kep,mu);
% Intitial position and velocity vector's components
rx0 = r(1);
ry0 = r(2);
rz0 = r(3);
vx0 = v(1);
vy0 = v(2);
vz0 = v(3);
% Initial condition vector
X0 =[rx0, ry0, rz0, vx0, vy0, vz0];
% Time Setting
if nargin <= 7
Period = ((2*pi)/sqrt(mu))*kep(1)^(3/2);
time = (0:86400:Period);
end
% Ode setting
options = odeset('Reltol',1e-13,'AbsTol',1e-14);
% Integration with ode113 function
[~,X] = ode113(@(t,X) dyn_orb_eq(t,X,mu),time,X0,options);
% Definition of the Output after from integration values
rx_vect = X(:,1);
ry_vect = X(:,2);
rz_vect = X(:,3);
vx_vect = X(:,4);
vy_vect = X(:,5);
vz_vect = X(:,6);
end