Lực hấp dẫn

  • Video minh hoạ

  • Code chương trình

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
function chuyendongvetinh_2D
clc
close all
clear all

%% CONSTANTS
G = 6.67e-11;
M = 5.97219e24;
R = 6400e3;

%% INPUT DATA
x = 2*R;
y = 0;
v = 6.5e3; % m/s
alpha = 90;

dt = 100; % s
tmax = 50*3600;

%% Xu ly
t = 0;
alpha = alpha/180*pi;
vx = v*cos(alpha);
vy = v*sin(alpha);
x_array = x/R;
y_array = y/R;

L = 1;

%% CALCULATION
while t<tmax
t = t+dt;
r = sqrt(x*x+y*y);
a = G*M/r/r;
ax = -a*x/r;
ay = -a*y/r;
vx = vx+ax*dt;
vy = vy+ay*dt;
x = x+vx*dt;
y = y+vy*dt;

x_array = [x_array x/R];
y_array = [y_array y/R];

if L<r/R
L = r/R;
end
end

%% FIGURE
fig1 = figure('name','Luchapdan','color','white','numbertitle','off');
hold on
viscircles([0 0],1);
xlabel('x [Earth radius]');
ylabel('y [Earth radius]');

axis equal
axis([-L L -L L]);
comet(x_array,y_array);

end