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
| clear
clc
G=6.67*10^-11;
//MASSES (kg)
M_Terre=5.98*10^24;//Masse de la Terre (kg)
M_Soleil=1.989*10^30;
//RAYONS (m)
R_Terre=6.371*10^6;//Rayon de la Terre (en m)
R_Soleil=696340*10^3;
//COORDONNEES DE CHAQUE ASTRE (sous la forme (x,y,z))
//Soleil c'est (0,0,0)
coordTerre=[149597870*10^3 0 0];//Distance terre-soleil puisque l'origine du repère est le soleil
function u_point=f(t,u)//vecteur u_point régissant les equations du mouvement
//Norme des vecteurs r avec pour origine chaque astre
r=sqrt(u(1)^2+u(2)^2+u(3)^2);//Norme r ayant pour origine l'origine du repère donc le premier astre(soleil)
r2=sqrt((u(1)-coordTerre(1))^2+(u(2)-coordTerre(2))^2+(u(3)-coordTerre(3))^2);//r2 en partant de la terre
M=[M_Soleil M_Terre];//Vecteur des masses
r_cube=[r^3 r2^3];//Les rayons puissance 3
vect_somme=M./r_cube;
K=-G*sum(vect_somme);
if r<R_Soleil then
u_point=[0;0;0;0;0;0];
elseif r2<R_Terre then
u_point=[0;0;0;0;0;0];
else
A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;K 0 0 0 0 0;0 K 0 0 0 0;0 0 K 0 0 0];//matrice A de l'équation u_point=A*u+B
B=[0;0;0;G*M_Terre*coordTerre(1)/r2^3;G*M_Terre*coordTerre(2)/r2^3;G*M_Terre*coordTerre(3)/r2^3];//Ici c'est une somme mais elle n'est pas écrite car il n'y a que le soleil et la terre, il faudra faire la somme si il y a un 3eme astre
u_point=A*u+B;//Utiliser une autre fonction à cause de B
end
endfunction
function [x,y,z]=sphere(phi,teta)//calcul de la sphere
x=cos(phi).*cos(teta);
y=cos(phi).*sin(teta);
z=sin(phi);
endfunction
function f=afficherAstre(origine,rayon)
//origine est un vecteur [x y z] representant les coordonnées de l'astre
//rayon est une valeur donnant le rayon de l'astre
[x,y,z]=eval3dp(sphere,linspace(-%pi/2,%pi/2,40),linspace(0,2*%pi,20));//calcul des facettes de la sphere
//on va positionner et dimensionner notre astre
[nbLignes_x,nbColonnes_x]=size(x);
x=x+origine(1)*ones(nbLignes_x,nbColonnes_x);//on positionne l'origine de la sphere
[nbLignes_y,nbColonnes_y]=size(y);
y=y+origine(1)*ones(nbLignes_y,nbColonnes_y);//on positionne l'origine de la sphere
[nbLignes_z,nbColonnes_z]=size(z);
z=z+origine(1)*ones(nbLignes_z,nbColonnes_z);//on positionne l'origine de la sphere
//on dimensionne notre astre pour avoir le bon rayon
x=x.*rayon;
y=y.*rayon;
z=z.*rayon;
clf();
plot3d(x,y,z);
endfunction
function u=rotation(distance,v0,hours)
//distance est la distance entre l'objet et l'astre (en km)
//v0 est la vitesse initiale (m/s)
//hours est le temps de la simulation en heures
u0=[R_Soleil+distance; 0; 0;v0/sqrt(3);v0/sqrt(3);v0/sqrt(3)];//vecteur v0 initial de la forme [x y z x_point y_point z_point]
distance=distance*1000;
t=linspace(0,3600*hours,2000);
u=ode(u0,0,t,f);//resolution du systeme differentiel
afficherAstre([0 0 0],R_Soleil);//affichage de l'astre
//tracer la trajectoire
comet3d(u(1,:),u(2,:),u(3,:),"colors",3);
endfunction
//DEBUT DU PROGRAMME
//Parametres initiaux
altitude=10^10;//en km sans compter le rayon du soleil donc c'est bien l'altitude
vitesse=20000;
duree=1000;
rotation(altitude,vitesse,duree); |
