Bonjour,
J'ai réaliser un programme permettant la visualisation d'orbite de satellites géostationnaire autour de la Terre, mais malheureusement, je n'arrive pas à obtenir de bonne conditions initiales afin que celui-ci m'offre des résultats cohérent. Pour un système terre/lune, mes CI marchent très bien, mais dès que je passe aux satellites plus rien.
J'utilise la méthode de Verlet afin d'obtenir de nouvelles positions en fonction du temps, je me disais aussi que c'était peut être le pas de temps qui était mauvais.
Je vous envoie en photo le code, et si jamais quelqu'un aurait une idée sur quoi changer (CI ou pas de temps ou ???) pour obtenir des trajectoires exploitables avec des satellites.
En vous remerciant,
BosonXYZT ;)
Voici mon code (désolé je n'ai pas réussir à mettre de photo)
#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
using namespace std;
void Satellites(int& choix, double& m, double& X, double& Y, double& Z, double& VX, double& VY, double& VZ) //fonction définissant la Terre et les satelittes
{
if (choix == 1) // TERRE
{
m =5.9736E+24; //masse
X=-1.315839435E+11; Y=6.763707294E+10; Z = 4.136174231E+05; //positions initiales
VX=0; VY=0; VZ=0; //vitesses initiales
}
else if (choix == 3) // Lune
{
m =7.3477E+22;
X=-1.319042007E+11; Y=6.786886988E+10; Z =-1.649538549E+07;
VX=-614.453; VY=-783.202; VZ=79.439;
}
else if (choix == 2) // SATELLITES
{
m =1000;
X=-36000000; Y=0; Z =0;
VX=0; VY=3074; VZ=0;
}
else
{
cout<<"Veuillez réessayer"<<endl;
}
}
void CalculForces(double FX[],double FY[], double FZ[],double X[],double Y[],double Z[],double m[],int N) //fonction définissant les forces qui s'éxercent sur les différents objets célestes
{
double G=6.6743e-11; //constante gravitationnelle
for(int i=0;i<N;i++)
{
FX[i]=0; // on initialise les forces à 0
FY[i]=0;
FZ[i]=0;
for(int j=0;j<N;j++)
{
if (i!=j)
{
double d = pow(X[i]-X[j],2)+pow(Y[i]-Y[j],2)+pow(Z[i]-Z[j],2);
d = pow(d,1.5);
FX[i]+=G*m[i]*m[j]/d*(X[j]-X[i]); //formules des forces sur chaque directions
FY[i]+=G*m[i]*m[j]/d*(Y[j]-Y[i]);
FZ[i]+=G*m[i]*m[j]/d*(Z[j]-Z[i]);
}
}
}
}
void Verlet(double X[], double Y[], double Z[], double VX[], double VY[], double VZ[], double FX[], double FY[], double FZ[], double m[], int N, double h) //algorithme de propagation avec la methode de verlet
{
CalculForces(FX,FY,FZ,X,Y,Z,m,N); //on appelle la fonction qui calcule les forces
for (int i = 0; i < N; i++)
{
X[i] = X[i] + VX[i]*h + 0.5*FX[i]/m[i]*pow(h,2); //fonction qui calcul les nouvelles positions en fonction du temps
Y[i] = Y[i] + VY[i]*h + 0.5*FY[i]/m[i]*pow(h,2);
Z[i] = Z[i] + VZ[i]*h + 0.5*FZ[i]/m[i]*pow(h,2);
VX[i] = VX[i] + 0.5*h*FX[i]/m[i]; //fonction qui calcul les nouvelles vitesses une première fois en fonction du temps
VY[i] = VY[i] + 0.5*h*FY[i]/m[i];
VZ[i] = VZ[i] + 0.5*h*FZ[i]/m[i];
}
CalculForces(FX,FY,FZ,X,Y,Z,m,N); // on rappelle la fonction qui calcule les forces
for (int i = 0; i < N; i++)
{
VX[i] = VX[i] + 0.5*h*FX[i]/m[i]; //fonction qui calcul les nouvelles vitesses une deuxième fois en fonction du temps
VY[i] = VY[i] + 0.5*h*FY[i]/m[i];
VZ[i] = VZ[i] + 0.5*h*FZ[i]/m[i];
}
enregistrer(N,X,Y,Z); //on enregistre les valeurs des positions dans un fichier
}
int main()
{
int N;
cout << "Combien de satellites voulez-vous visualiser ?" << endl;
cin >> N; //on choisit combien de satelittes on veut étudier
int choix[N];
double X[N], Y[N], Z[N], VX[N], VY[N], VZ[N], FX[N], FY[N], FZ[N], m[N] ;
for (int i = 0; i < N; i++)
{
cout << "Satellites " << i+1 << endl;
do {
cout << "choisir :" << endl;
cout << "1 : Terre" << endl;
cout << "2 : 1 satellites" << endl;
cout << "3 : 2 satellites" << endl;
cout << "4 : 3 satellites" << endl;
cout << "5 : 4 satellites" << endl;
cout << "6 : 5 satellites" << endl;
cout << "7 : 6 satellites" << endl;
cin >> choix[i]; //choix de l'utilisateur
} while (choix[i] <1 || choix[i] > 7);
Satellites(choix[i], m[i], X[i],Y[i],Z[i],VX[i],VY[i],VZ[i]);
}
int dt = 100; //intervalle de temps entre chaque point
fstream f; //ouverture du ficher où l'on ecrira l'énergie totale ainsi que les aires calculée
f.open("etot.txt", ios::out); //on ouvre un fichier pour rentrer les valeurs de l'énergie totale ainsi que celles obtenues avec la loi des aires
double C=0;
for (int t = 0; t <=3600*24*10; t = t + dt)
{
Verlet(X,Y,Z,VX,VY,VZ,FX,FY,FZ,m,N,dt); //on appelle l'algorithme de Verlet pour calculer positions et vitesses
if (t%10==0) //intervalle de temps pour calculer l'énergie totale ainsi que la loi des aires, on prend tout les points ayant un rest de divisiob euclidienne par 10 égal à 0
{
double EC = Ec(VX, VY, VZ, m, N); //on calcul l'énergie cinétique
double EP = Ep(X, Y, Z, m, N);//on calcul l'énergie potentielle
double Etot = EC + EP; //calcul énergie totale
f<<Etot<<'\t'; //ecriture de l'énergie totale dans le fichier
C=loi_aire(m,VX,VY,VZ,X,Y,Z); //on appelle la fonction qui calcul la loi des aires
f<<C<<endl; //ecriture des valeurs calculée par la loi des aires
}
}
f.close(); //fermeture du fichier
return 0;
}
# pow ?
Posté par djibb (site web personnel) . Évalué à 2.
Coucou,
la racine carré, c'est puissance 0,5 non ?
pas 1,5 ?
[^] # Re: pow ?
Posté par lann . Évalué à 1.
Ou sqrt
[^] # Re: pow ?
Posté par BosonXYZT . Évalué à 1.
C'est racine de d3
donc d3/2 = d1.5
# Henri Poincarré
Posté par Luc-Skywalker . Évalué à 4.
génial mathématicien a démontré en 1888 que le système à trois corps était sensible aux conditions initiales.
"Si tous les cons volaient, il ferait nuit" F. Dard
# Quelques remarques
Posté par ml . Évalué à 8. Dernière modification le 12 avril 2023 à 23:45.
En vrac :
Normalement on utilise le barycentre des masses comme origine du repère. On peut facilement prouver que le barycentre a une vitesse nulle (au moins pour 2 corps, pour >2 je me rappelle plus faudrait faire le calcul). Il faut absolument que tu corriges cela, car là tu bouffes de la précision numérique avec un système pas centré.
Fait moins de snapshot (1 sur 10 pas, comme ton calcul de conservation, c’est encore trop fréquent àmha) sinon ça va te bouffer trop de temps. Probablement que tes pas sont trop grand si tu fais un snapshot par pas.
Pour t’aider à trouver la bonne taille de pas, tu peux faire un jeu de simu en réduisant le pas, et voir à quel moment le résultat se stabilise (ie. réduire le pas ne change pas fondamentalement le résultat de la simu).
La loi_aire n’est pas donnée, mais fondamentalement il s’agit de la conservation du moment cinétique (c’est un vecteur). Tu gagnerais à le formaliser comme ça car la conservation du moment cinétique permet de s’assurer que ton orbite est bien dans un plan (par conservation de la direction du moment) et que tu obéis à la loi des aires (par conservation de l’amplitude du moment). De plus as-tu bien vérifié que les deux grandeurs (énergie et moment) sont bien conservées ?
Normalement on normalise un code physique : G vaut 1, toutes les grandeurs (variables, masses, etc.) sont autour de 1. Par exemple la masse du corps principal est fixée à 1 et celle des corps qui orbitent à une fraction de la masse du corps principal. Le corps qui orbite est à une distance 1. C’est l’idée grosso modo. La conversion en unités physique se fait à part… éventuellement.
Pour quelque chose de plus pro. les conditions initiales (CI) devraient être dans un fichier — identique au format de tes snapshots, car des CI sont un snapshot de ta simu (tu pourras ainsi arrêter et reprendre ta simu sur le dernier snapshot si tu veux). Il te faudra éventuellemnt un petit utilitaire pour construire un snapshot ad-hoc pour démarrer une simu (mais pour commencer je te recommande de faire des snapshots ascii tabulé directement exploitable à la main).
Fais des CI simple, dans le plan XY, avec une orbite circulaire, etc. Dans le cas à deux corps, l’orbite se calcule de manière analytique (tu retrouveras facilement sur le web les calculs je pense). Ça te permettra de valider ton programme par comparaison au résultat analytique, avec un jeu de fichiers standards (d’où l’importance du point 6). Avec des CI simples, ça facilitera la visualisation graphique des snapshot aussi : suffira de tracer les points (X,Y), vérifier que Z reste nul, etc.
Il y aura une dérive numérique, même dans le cas à deux corps, surtout avec une méthode “naïve” (si je me souviens bien il y en a qui sont conçues pour conserver les deux grandeurs énergie & moment, et donc plus robustes, mais ça devient beaucoup plus difficile à implémenter). Faut pas s’affoler si l’orbite n’est pas très stable, même après à peine 10 révolutions tu pourrais voir une décalage de l’orbite, signe qu’il faudra réduire la taille du pas (au prix d’une simu plus longue).
Optimisation. La force du corps A sur le corps B est l’opposée de la force du corps B sur le corps A. (À toi de retrouver comme exploiter cela dans ton code.)
Optimisation. La trajectoire d’un corps ne dépend pas de sa masse. (À toi de trouver comment exploiter cela dans ton code.)
[^] # Re: Quelques remarques
Posté par ǝpɐןƃu∀ nǝıɥʇʇɐW-ǝɹɹǝıԀ (site web personnel) . Évalué à 4.
Comme vous faites le tour de pas mal de points, j’en profite pour ajouter mon grain de sel.
Bon amusement.
« IRAFURORBREVISESTANIMUMREGEQUINISIPARETIMPERAT » — Odes — Horace
# Confusion altitude et distance au centre de la Terre
Posté par ǝpɐןƃu∀ nǝıɥʇʇɐW-ǝɹɹǝıԀ (site web personnel) . Évalué à 4. Dernière modification le 14 avril 2023 à 19:22.
Tout est dans le titre. De mémoire, l’altitude des satellites géostationnaires serait de l’ordre de 36000 km, donc leur distance au centre de la Terre, de l’ordre de 42400 km. Ça devrait mieux fonctionner ainsi. Le même code fonctionnerait aussi avec le Soleil, les autres planètes et leurs satellites… Et même en ajoutant des échantillons des nuages d’astéroïdes et petit corps, ça passe assez bien su PC. Le C C++ rapide :-).
« IRAFURORBREVISESTANIMUMREGEQUINISIPARETIMPERAT » — Odes — Horace
Suivre le flux des commentaires
Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.