-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathexample_osculating2mean.m
117 lines (103 loc) · 3.8 KB
/
example_osculating2mean.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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
%% Example osculating2mean
%% Load osculating position timeseries
% Time-series of roughly 10 orbits of a satelite in LEO with atmospheric
% drag; cannon ball solar radiation pressure; third body perturbations from
% the Sun, Moon, Mars, Venus; spherical harmonic gravity of degree and
% order 12
load('./data/output.mat','x');
Ts = 10; %(s)
% Compute osculating orbital elements
OE_osc = zeros(6,size(x,2));
for t = 1:size(x,2)
OE_osc(:,t) = rv2OEOsc(x(:,t));
end
% Compute mean orbital elements
OE_mean_EcksteinUstinov = zeros(6,size(x,2));
OE_mean_EcksteinUstinovKaula = zeros(6,size(x,2));
degree = 10;
for t = 1:size(x,2)
OE_mean_EcksteinUstinov(:,t) = OEOsc2OEMeanEU(OE_osc(:,t));
OE_mean_EcksteinUstinovKaula(:,t) = OEOsc2OEMeanEUK((t-1)*Ts,OE_osc(:,t),degree);
end
%% Plot results
% Constants
mu = 3.986004418e14; %(m^3 s^-2)
RE = 6378.137e3; %(m)
J2 = 1082.6267e-6;
% Compute parameters
semiMajorAxis = mean(OE_osc(1,:));
incl = mean(OE_osc(5,:));
n = sqrt(mu/(semiMajorAxis)^3);
gamma = (J2/2)*(RE/semiMajorAxis)^2;
Omega_dot = -3*gamma*n*cos(incl);
arg_perigee_dot = (3/2)*gamma*n*(5*cos(incl)^2-1);
M_dot = (3/2)*gamma*n*(3*cos(incl)^2-1);
u_dot = n + M_dot + arg_perigee_dot;
T = size(x,2);
figure;
hold on
ylabel('$a$ (m)','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,OE_osc(1,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinov(1,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinovKaula(1,1:T));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
figure;
hold on
ylabel('$a$ (m)','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,OE_osc(1,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinov(1,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinovKaula(1,1:T));
ylim([min(OE_mean_EcksteinUstinov(1,1:T)) max(OE_mean_EcksteinUstinov(1,1:T))])
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
aux = zeros(3,T);
aux(1,:) = OE_osc(2,:)-u_dot*(0:1:T-1)*Ts;
aux(2,:) = OE_mean_EcksteinUstinov(2,:)-u_dot*(0:1:T-1)*Ts;
aux(3,:) = OE_mean_EcksteinUstinovKaula(2,:)-u_dot*(0:1:T-1)*Ts;
aux(aux < -pi) = aux(aux <- pi) - 2*pi*floor((aux(aux <- pi)+pi)/(2*pi));
figure;
hold on
ylabel('$u-\dot{u}t$ (rad)','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,aux(1,:));
plot(0:Ts:(T-1)*Ts,aux(2,:));
plot(0:Ts:(T-1)*Ts,aux(3,:));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
figure;
hold on
ylabel('$e_x$','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,OE_osc(3,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinov(3,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinovKaula(3,1:T));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
figure;
hold on
ylabel('$e_y$','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,OE_osc(4,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinov(4,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinovKaula(4,1:T));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
figure;
hold on
ylabel('$i$ (rad)','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,OE_osc(5,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinov(5,1:T));
plot(0:Ts:(T-1)*Ts,OE_mean_EcksteinUstinovKaula(5,1:T));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');
aux = zeros(3,T);
aux(1,:) = OE_osc(6,1:T)-Omega_dot*(0:1:T-1)*Ts;
aux(2,:) = OE_mean_EcksteinUstinov(6,1:T)-Omega_dot*(0:1:T-1)*Ts;
aux(3,:) = OE_mean_EcksteinUstinovKaula(6,1:T)-Omega_dot*(0:1:T-1)*Ts;
aux(aux > pi) = aux(aux > pi) - 2*pi*floor((aux(aux > pi)+pi)/(2*pi));
figure;
hold on
ylabel('$\Omega-\dot{\Omega}t$ (rad)','Interpreter','latex')
xlabel('$t$ (s)','Interpreter','latex')
plot(0:Ts:(T-1)*Ts,aux(1,:));
plot(0:Ts:(T-1)*Ts,aux(2,:));
plot(0:Ts:(T-1)*Ts,aux(3,:));
legend('Osculating','EcksteinUstinov','EcksteinUstinovKaula');