-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathM3_Interplanetary.asv
122 lines (80 loc) · 3.4 KB
/
M3_Interplanetary.asv
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
116
117
118
clear; clc;
load SOI_OUT.mat X_SOI T_SOI V_SOI
jdt = juliandate(datetime([2043, 12, 28, 16, 44, 0]) + seconds(T_SOI));
%% Object Initilization
uf = UtilityFunctions();
earth = CelestialObject("Earth", 5.9722e24, 6371.0, 149.598e6, 23.44, jdt); % name, mass, planetary radius, heliocentric radius, jdt
mars = CelestialObject("Mars", 0.64169e24, 3389.5, 227.956e6, 0, jdt); % name, mass, planetary radius, heliocentric radius, jdt
% Spacecraft SOI-Exit conditions
V_SOI(3) = 0;
X_SOI(3) = 0;
X_i = X_SOI' + earth.heliocentric_pos;
V_i = V_SOI' + earth.heliocentric_vel;
V_i = 0.996 * V_i;
%% Setup Geometry and Plots
figure(1)
earth_plot = plot(earth.heliocentric_pos(1), earth.heliocentric_pos(2), "b.", "MarkerSize", 10);
hold on
grid on
axis equal
xlim([-3e8 3e8])
ylim([-3e8 3e8])
quiver(earth.heliocentric_pos(1), earth.heliocentric_pos(2), earth.heliocentric_vel(1), earth.heliocentric_vel(2), 1e6, "blue");
mars_plot = plot(mars.heliocentric_pos(1), mars.heliocentric_pos(2), "r.", "MarkerSize", 10);
quiver(mars.heliocentric_pos(1), mars.heliocentric_pos(2), mars.heliocentric_vel(1), mars.heliocentric_vel(2), 1e6, "red")
sun = plot(0, 0, "y.", "MarkerSize", 30);
SC_plot = plot(X_i(1), X_i(2), "r.", "MarkerSize", 10);
SC_quiver = quiver(X_i(1), X_i(2), V_i(1), V_i(2), 1e6, "red");
%% Analytical Orbital Elements
V_ri = dot(uf.hat(X_i), V_i);
h_vector = cross(X_i, V_i);
h = norm(h_vector);
e_vector = (1/earth.mu_sun) * ((norm(V_i)^2 - earth.mu_sun / norm(X_i)) * X_i - norm(X_i) * V_ri * V_i);
eccentricity = norm(e_vector);
r_p = (h^2 / earth.mu_sun) ./ (1 + eccentricity * cosd(0));
r_a = (h^2 / earth.mu_sun) ./ (1 + eccentricity * cosd(180));
a = 0.5 * (r_p + r_a);
T_orbit = sqrt(a^3 * 4 * pi^2 / earth.mu_sun);
t_anomaly_now = uf.angle_between(e_vector, X_i);
E = 2 * atan(sqrt((1-eccentricity)/(+1+eccentricity)) * tand(t_anomaly_now/2));
Me = E - eccentricity * sin(E);
t_periapsis = Me * T_orbit / (2*pi);
arg_periapsis = uf.angle_between(e_vector, [1, 0, 0]);
peri_rot = [cosd(arg_periapsis), -sind(arg_periapsis); sind(arg_periapsis) cosd(arg_periapsis)];
t_anomaly = 0:0.5:360;
r = (h^2 / earth.mu_sun) ./ (1 + eccentricity * cosd(t_anomaly));
X_transfer = r' .* [cosd(t_anomaly'), ....
sind(t_anomaly')];
X_Mars = mars.r_orbit .* [cosd(t_anomaly'), ....
sind(t_anomaly')];
X_Earth = earth.r_orbit .* [cosd(t_anomaly'), ....
sind(t_anomaly')];
for i = 1:length(X_transfer)
X_transfer(i, :) = peri_rot * X_transfer(i, :)';
end
plot(X_transfer(:, 1), X_transfer(:, 2))
plot(X_Mars(:, 1), X_Mars(:, 2), "k")
plot(X_Earth(:, 1), X_Earth(:, 2), "k")
N = 1000;
T = zeros(1, N);
X_SC = zeros(N, 2);
V_SC = X_SC;
dt = 86400/2; % seconds
for i = 1:N
t = T(i);
Me = 2 * pi * (t + t_periapsis) / T_orbit;
t_fun = @(E) E - eccentricity * sin(E) - Me;
E = fzero(t_fun, Me);
t_anomaly_now = rad2deg(2 * atan(tan(E/2) * sqrt((1+eccentricity)/(1-eccentricity))));
r_now = (h^2 / earth.mu_sun) ./ (1 + eccentricity * cosd(t_anomaly_now));
X_SC(i, :) = peri_rot * r_now * [cosd(t_anomaly_now); sind(t_anomaly_now)];
[SC_plot.XData, SC_plot.YData] = [X_SC(i, 1), X_SC(i, 2)];
earth = earth.refresh(dt);
mars = mars.refresh(dt);
earth_plot.XData = earth.heliocentric_pos(1);
earth_plot.YData = earth.heliocentric_pos(2);
mars_plot.XData = mars.heliocentric_pos(1);
mars_plot.YData = mars.heliocentric_pos(2);
pause(0.01)
T(i + 1) = t + dt;
end