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
| function [u2]=force_g(t,u,masse)
module=-G*masse*((u(1)^2+u(2)^2)^(-3/2))
u2=[module*u(1); module*u(2)]
endfunction
function [du]=force(t,u,masse0,masse1)
u1=[u(1);u(2)]
du1=[u(3);u(4)]
u2=[u(5);u(6)]
du2=[u(7);u(8)]
ddu1=force_g(t,u1,masse0)
ddu2=force_g(t,u2,masse0)+force_g(t,u2-u1,masse1)
du=[du1;ddu1;du2;ddu2]
endfunction
//constantes
G=0.04;//constante gravitationnelle en unités astronomique*
m0=1000;//masse soleil
m1=1;//masse Jupiter
dt=0.05;// intervalle entre les différentes dates
T=500;//durée de la simulation
dx=0.5;dy=0.5;// taille du dessin des planètes
//coordonnées de la planète au début
x1=5;y1=0;vx1=0;vy1=2.5;
//coordonnées de la comète au début
x2=6;y2=6;vx2=-2;vy2=-0.5;
// valeur initiale de U
u0=[x1; y1; vx1; vy1; x2; y2; vx2; vy2];
//dates de calcul de U
t=[0:dt:T];
//calcul de la solution
u=ode(u0,0,t,list(force,m0,m1));
//récupération des coordonnées correspondant à la position (X,Y) des deux objets
X=[u(1,:)',u(5,:)'];Y=[u(2,:)',u(6,:)'];
//animation graphique
plot2d(X,Y,[5 2],rect=[-10,-10,10,10],frameflag=4)
rect=[min(X),min(Y),max(X),max(Y)]
N=length(t);
plot2d(0,0,-9,rect=rect,frameflag=3)
f=gcf()
winnum=winsid() ;
toolbar(winnum(1),'off');
f.pixmap="on"
for k=2:N
drawlater()
clf()
plot2d(0,0,-9,rect=rect,frameflag=3)
plot2d(X(1:k,:),Y(1:k,:),[5 2],rect=rect,frameflag=0)
P1=[X(k,1)-dx/2;Y(k,1)+dy/2;dx;dy;0;360*64];//la planète en rouge
P2=[X(k,2)-dx/2;Y(k,2)+dy/2;dx;dy;0;360*64];//la comète en bleu
xfarcs([P1,P2],[5,2])
xstring(-4,-6,'t='+string(t(k)))
drawnow()
end |
