-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpset1_dim.m
97 lines (84 loc) · 2.32 KB
/
pset1_dim.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
clear all
close all
set(0,'DefaultLineLineWidth',1.5)
set(0,'DefaultLineMarkerSize',15)
set(0,'DefaultAxesFontSize',22)
set(0,'DefaultTextFontSize',26)
%% Problem Set 1 Problem 9
constants
%% Initial Conditions
mu1 = 328900.56^-1; % From NASA
omega = sqrt((mu.earth + mu.sun)/dis.sun^3);
%% solving L1
syms x
eqn = x - (1-mu1)/(x+mu1)^2 + mu1/(x-1+mu1)^2 == 0;
digitsOld = digits(100);
L1 = vpasolve(eqn,x,[0,1]);
highPrecisionL1 = double(L1) ;
digits(digitsOld)
%% Initial r and v
r = [.98900883730044109; 0; 0.000802140914099732]*dis.sun;
v = [0; 0.010697471710460349-2.896403e-3; 0]*dis.sun*omega; % be careful changing. Very sensitive.
state_sc_earth = [r; v; omega; mu.sun; mu.earth; mu1];
%% solving the orbit
tstart = 0;
tint = 25;
tend = 3600*24*177.15; % convert to nondimensional time
% tend = 3600*24*365*omega; % convert to nondimensional time
options = odeset('RelTol', 1e-3, 'AbsTol', 1e-6);
[t_out, y_out] = ode113(@orbit_prop_dim, [tstart:tint:tend]', state_sc_earth, options);
%% Plotting
figure()
plot(y_out(:,1), y_out(:,2),'k','LineWidth',3)
hold on
% plot(y_out(:,7), y_out(:,8))
scatter(double(L1)*dis.sun,0,100,'filled','r')
% scatter(mu1*dis.sun,0,25,'filled')
scatter((1-mu1)*dis.sun,0,100,'filled','b')
title('Halo Orbit around L1 Point in the Rotating Barycenter Frame')
xlabel('x position, km')
ylabel('y position, km')
% axis equal
grid on
legend('Halo Orbit', 'L1 Point', 'Earth')
hold off
(y_out(1,1) - y_out(end,1)) %check on final error in the distance
figure()
subplot(3,1,1)
plot3(y_out(:,1),y_out(:,2),y_out(:,3),'k','LineWidth',3)
view([0 -90 0])
hold on
scatter3(double(L1)*dis.sun,0,0,100,'filled','r')
xlabel('x, km')
zlabel('z, km')
grid on
hold off
subplot(3,1,2)
plot3(y_out(:,1),y_out(:,2),y_out(:,3),'k','LineWidth',3)
hold on
scatter3(double(L1)*dis.sun,0,0,100,'filled','r')
view([90 0 0])
ylabel('y, km')
zlabel('z, km')
grid on
hold off
subplot(3,1,3)
plot3(y_out(:,1),y_out(:,2),y_out(:,3),'k','LineWidth',3)
view([0 0 90])
hold on
scatter3(double(L1)*dis.sun,0,0,100,'filled','r')
xlabel('x, km')
ylabel('y, km')
grid on
hold off
sgtitle(['Different views of the Halo Orbit around the {\color{red}L1} point'])
%
% hold on
% scatter3(double(L1)*dis.sun,0,0,100,'filled','r')
% title('Halo Orbit')
% xlabel('x position, km')
% ylabel('y position, km')
% zlabel('z position, km')
% axis equal
% grid on
% hold off