Simulateur perso : bug introuvable (résolu)

Bonjour à tous, j'ai besoin de vous pour m'aider à répondre à une question qui me chiffonne depuis un petit moment.

Voici le contexte. Je développe pendant mon temps libre un "simulateur" informatique de mécanique spatiale. Le simulateur en question applique, sur chaque "corps" mis en scène, les forces gravitationnelles engendrées par la masse des autres corps. Ce calcul est fait en moyenne 300 fois par seconde, et se répète inlassablement.

J'en suis arrivé à un point où j'ai modélisé le système solaire (soleil + 9 planètes), et où j'introduis un petit corps en orbite autour de la Terre : une boule de pétanque :-)

La boule de pétanque est située à 300 km d'altitude et possède une vitesse initiale égale à 8000 m.s-1, orientée à l'horizontale par rapport à la surface terrestre. J'en arrive au problème.

A) Je lance  la simulation en soumettant la boule uniquement à l'attraction de la Terre : elle va suivre une orbite très régulière de 300 km périapse (le point de départ) pour environ ~1500 km apoapse. Elle peut tourner comme ça très longtemps y compris si la Terre est elle-même en mouvement autour du soleil.

B) Je lance la simulation en soumettant la boule à l'attraction de tous les corps célestes  du système solaire : cette fois-ci l'orbite va être tellement perturbée que je crois que j'ai un bug. En une seule orbite, le périapse de la boule aura perdu environ 50 km, et en dix heures environ, la boule se sera écrasée sur Terre. Ça me paraît grotesque comme résultat, avec mon logiciel l'ISS aurait coulé depuis belle lurette.

Je précise que oui, je simule le freinage atmosphérique (infinitésimal à cette altitude sur un si petit objet), que le simulateur fonctionne en 2D (imaginez que tous les corps évoluent dans le plan de l'écliptique) et je précise également que la vitesse initiale de la boule est en direction du soleil, ce qui fait que quand je démarre la simulation B, l'attraction du soleil si petite soit-elle semble attirer la boule suffisamment pour faire monter significativement son apoapse. Quand on arrive à l'apoapse, le soleil freine la boule et fait donc baisser son périapse.

Si j'ai besoin de vous à ce stade de mon debuggage infructueux, c'est déjà pour me rassurer : pouvez-vous me confirmer que le résultat de la simulation B est aberrant par rapport à ce qu'on observerait dans la réalité, et pouvez-vous me dire si la pseudo explication que je me fais quant à l'effet du soleil sur l'orbite B semble plausible ?

Je suis tout disposé à vous donner plus de paramètres où à lancer la simulation différemment si cela peut vous aider à y voir plus clair. Merci d'avance à qui pourra me faire progresser dans ma recherche.

EDIT : je change le sujet, quand je vois les premières réponses il n'y a pas trop de doute et ça confirme mon intuition : c'est bien un bug.
EDIT 01/09/2014 : je change le sujet, le bug est probablement identifié?
EDIT 05/09/2014 : je change le sujet, le bug est résolu, merci à tous !

Je te confirme que c'est un bug: en LEO, l'influence du soleil sur la stabilité court-terme des orbites est vraiment faible.
Cela provient principalement du fait que le delta de gravité aux points extrêmes de l'orbite sont très faibles (parce que le soleil est très loin, indépendamment de sa masse) relativement au potentiel terrestre.
Donc avec ou sans soleil, ton orbite devrait probablement être davantage impactée par la précision numérique de tes calculs différentiels que par la présence du soleil.

Merci pour ton message, ça confirme donc ce que je pensais. Pour l'étape suivante, je vais préparer un scénario de test épuré de tout ce qui est inutile, et je posterai les calculs des deux où trois premières millisecondes de vol, peut-être qu'un détail émergera pour expliquer l'erreur.

De toute façon ce que tu fais là (hors freinage) tient en quelques lignes de code, tu auras vite fait de tout passer en revue...

Si ça tourne correctement lorsqu'il n'y a que la Terre et que ça dérape en ajoutant d'autres corps, ça peut être simplement des paramètres faux, erreurs d'unités par exemple (distance des planètes données en km plutôt qu'en m, masse fausse d'un rapport 1000, etc.)

Space Opera a écrit:De toute façon ce que tu fais là (hors freinage) tient en quelques lignes de code, tu auras vite fait de tout passer en revue...

Puisses-tu avoir raison...

lambda0 a écrit:Si ça tourne correctement lorsqu'il n'y a que la Terre et que ça dérape en ajoutant d'autres corps, ça peut être simplement des paramètres faux, erreurs d'unités par exemple (distance des planètes données en km plutôt qu'en m, masse fausse d'un rapport 1000, etc.)

J'y avais pensé, mais je ne trouve pas. Et si c'est ça, ce sera tellement énorme que je serai tout honteux de ne pas l'avoir vu.

Bon, j'ai généré des traces intermédiaires du programme, en essayant de vous donner le maximum d'info utile. Voici les paramètres :

Repère orthonormé X, Y.
Position des corps exprimée en coordonnées polaires (r, θ)
r = distance au centre
θ = angle par rapport à l'axe X, dans le sens trigonométrique, en radians
Les forces, accélérations et vitesses sont également exprimées en coordonnées polaires.
Incrément de temps entre chaque calcul : 10 ms (0.01 seconde)

Corps célestes :

- SOLEIL

Position centrale : (0, 0)
Déplacement : aucun
Masse : 1.98855E+30 kg

- TERRE

Position : (1.5080513000128702E8 km, -0.0958663954751644 rad soit -5,493°, en gros la Terre est "à droite", légèrement sous l'axe X)
Déplacement : aucun (ni ne suit son orbite naturelle, ni ne tombe vers le soleil)
Masse : 5.9736E+24 kg
Rayon : 6371 km
Atmosphère : aucune (pas de freinage atmosphérique)

Objet en orbite :

- Boule de pétanque

Position initiale par rapport au centre de la Terre (6671 km, PI/2) soit 300 km d'altitude, "en haut" de la Terre
Pour des raisons de commodités, je donnerai la position de la boule par rapport à la Terre
Masse = 700 g (0.7 kg)

Résumé grossier de la situation initiale :

Code:

           Soleil
              |
              v
             ___
            /   \
Axe X -----(     )---------loin-loin-loin-------------------------------------------------------------------
            \___/                                                             .       <-- Boule de pétanque
                                                                              o       <-- Terre

Simulation A : seule la gravité de la Terre est appliquée

Calcul	Temps (s)	Gravité Terre (m.s-2)	Poids (N)		Accélération (m.s-2)		Delta vitesse (m.s-1)		Vitesse résultante (m.s-1)		Delta position (m)		Position résultante (m)
			r	θ	r	θ	r	θ	r	θ	r	θ	r	θ
0	0,00								8000	PI			6671000	PI/2
1	0,01	8,958392	6,270874	4,712389	8,958392	-1,570796	0,089584	-1,570796	8000,000001	-3,141581	80	-3,141581	6671000	1,570808
2	0,02	8,958392	6,270874	4,712401	8,958392	-1,570784	0,089584	-1,570784	8000,000001	-3,14157	80	-3,14157	6671000	1,57082
3	0,03	8,958392	6,270874	4,712413	8,958392	-1,570772	0,089584	-1,570772	8000,000001	-3,141559	80	-3,141559	6671000	1,570832
4	0,04	8,958392	6,270874	4,712425	8,958392	-1,57076	0,089584	-1,57076	8000,000002	-3,141548	80	-3,141548	6671000	1,570844
5	0,05	8,958392	6,270874	4,712437	8,958392	-1,570748	0,089584	-1,570748	8000,000002	-3,141537	80	-3,141537	6671000	1,570856
6	0,06	8,958392	6,270874	4,712449	8,958392	-1,570736	0,089584	-1,570736	8000,000002	-3,141525	80	-3,141525	6671000	1,570868
7	0,07	8,958392	6,270874	4,712461	8,958392	-1,570724	0,089584	-1,570724	8000,000002	-3,141514	80	-3,141514	6671000	1,57088
8	0,08	8,958392	6,270874	4,712473	8,958392	-1,570712	0,089584	-1,570712	8000,000002	-3,141503	80	-3,141503	6671000	1,570892
9	0,09	8,958392	6,270874	4,712485	8,958392	-1,5707	0,089584	-1,5707	8000,000002	-3,141492	80	-3,141492	6671000	1,570904
10	0,10	8,958392	6,270874	4,712497	8,958392	-1,570688	0,089584	-1,570688	8000,000002	-3,141481	80	-3,141481	6671000	1,570916

Simulation B : la gravité de la Terre et du soleil sont appliquées

         

Calcul	Temps (s)	Gravité soleil (m.s-2)	Poids exercé par le soleil (N)		Gravité Terre (m.s-2)	Poids exercé par la Terre (N)		Poids résultant (N)		Accélération (m.s-2)		Delta vitesse (m.s-1)		Vitesse résultante (m.s-1)		Delta position (m)		Position résultante (m)
			r	θ		r	θ	r	θ	r	θ	r	θ	r	θ	r	θ	r	θ
0	0,00													8000,000000	PI			6671000,000000	PI/2
1	0,01	0,005836	0,004085	3,04577	8,958392	6,270874	4,712389	6,270484	-1,571445	8,957835	-1,571445	0,089578	-1,571445	8000,000059	-3,141581	80,000001	-3,141581	6670999,999584	1,570808
2	0,02	0,005836	0,004085	3,04577	8,958392	6,270874	4,712401	6,270484	-1,571433	8,957835	-1,571433	0,089578	-1,571433	8000,000117	-3,14157	80,000001	-3,14157	6670999,999231	1,57082
3	0,03	0,005836	0,004085	3,04577	8,958392	6,270874	4,712413	6,270484	-1,571421	8,957835	-1,571421	0,089578	-1,571421	8000,000176	-3,141559	80,000002	-3,141559	6670999,998942	1,570832
4	0,04	0,005836	0,004085	3,04577	8,958392	6,270874	4,712425	6,270484	-1,571409	8,957835	-1,571409	0,089578	-1,571409	8000,000234	-3,141548	80,000002	-3,141548	6670999,998717	1,570844
5	0,05	0,005836	0,004085	3,04577	8,958392	6,270874	4,712437	6,270484	-1,571397	8,957835	-1,571397	0,089578	-1,571397	8000,000292	-3,141537	80,000003	-3,141537	6670999,998555	1,570856
6	0,06	0,005836	0,004085	3,04577	8,958392	6,270874	4,712449	6,270484	-1,571385	8,957835	-1,571385	0,089578	-1,571385	8000,00035	-3,141525	80,000004	-3,141525	6670999,998457	1,570868
7	0,07	0,005836	0,004085	3,04577	8,958392	6,270874	4,712461	6,270484	-1,571373	8,957835	-1,571373	0,089578	-1,571373	8000,000409	-3,141514	80,000004	-3,141514	6670999,998423	1,57088
8	0,08	0,005836	0,004085	3,04577	8,958392	6,270874	4,712473	6,270484	-1,571361	8,957835	-1,571361	0,089578	-1,571361	8000,000467	-3,141503	80,000005	-3,141503	6670999,998452	1,570892
9	0,09	0,005836	0,004085	3,04577	8,958392	6,270874	4,712485	6,270484	-1,571349	8,957835	-1,571349	0,089578	-1,571349	8000,000525	-3,141492	80,000005	-3,141492	6670999,998544	1,570904
10	0,10	0,005836	0,004085	3,04577	8,958392	6,270874	4,712497	6,270484	-1,571337	8,957835	-1,571337	0,089578	-1,571337	8000,000583	-3,141481	80,000006	-3,141481	6670999,998701	1,570916

(désolé pour la mise en page ci-dessus, mais plus j'y touche, plus c'est déformé, donc j'arrête là)

Dans ta simulation, tu prends la Terre immobile ???
Ca pourrait être la cause de tes soucis.

Non, c'est juste pour essayer de trouver ce bug. Je n'observe pas de différence notable avec la Terre immobile ou en mouvement.

NB : les plus observateurs remarqueront peut-être que la vitesse de la boule augmente quand elle devrait diminuer (elle s'éloigne de son périapse) : ceci est dû à l'intégration qui pose quelques soucis sur les premiers calculs, au bout de 0.15 secondes tout rentre dans l'ordre.

Je viens de me coder ça en quelques lignes (6 petite lignes pour être exact), et ça me donne des résultats normaux. Je n'ai pas tous les résultats dans ton calcul, par exemple dans ton cas avec soleil, j'ai la Terre et la boule qui foncent vers le soleil si la Terre a une vitesse initiale nulle. Et l'orbite de la boule est bien circulaire.
Bon cependant, il faut absolument faire tes calculs en double précision, puisque du 10^30kg flirtent avec du 0.7kg. Peut-être la source de tes soucis ? Pour info chez moi, j'ai moins d'erreur numérique en considérant que la boule a une masse nulle (je n'explicite pas le fait qu'elle attire ni la Terre ni le soleil), ça simplifie les équations et ça diverge moins.

Oui je suis en double précision, j'ai un peu tronqué les chiffres pour que ça reste affichable sur le forum.

Je précise que pour la simulation B, je fais en sorte que la terre ne tombe pas sur le soleil, elle reste strictement immobile (mais comme dit plus haut, si elle est en mouvement ça ne change pas grand chose).

Trouves-tu comme moi la gravité du soleil égale à 0,005836 m.s-2 ? Ca me parait beaucoup.

Autre point : si tu trouves une orbite circulaire c'est étrange, tu devrais être à 300 km péri et ~1500 km apo.

Space Opera a écrit:Je viens de me coder ça en quelques lignes (6 petite lignes pour être exact), et ça me donne des résultats normaux. Je n'ai pas tous les résultats dans ton calcul, par exemple dans ton cas avec soleil, j'ai la Terre et la boule qui foncent vers le soleil si la Terre a une vitesse initiale nulle. Et l'orbite de la boule est bien circulaire.
Bon cependant, il faut absolument faire tes calculs en double précision, puisque du 10^30kg flirtent avec du 0.7kg. Peut-être la source de tes soucis ? Pour info chez moi, j'ai moins d'erreur numérique en considérant que la boule a une masse nulle (je n'explicite pas le fait qu'elle attire ni la Terre ni le soleil), ça simplifie les équations et ça diverge moins.

C'est quoi ton programme si c'est pas indiscret ?

Thierz a écrit:Oui je suis en double précision, j'ai un peu tronqué les chiffres pour que ça reste affichable sur le forum.

Je précise que pour la simulation B, je fais en sorte que la terre ne tombe pas sur le soleil, elle reste strictement immobile (mais comme dit plus haut, si elle est en mouvement ça ne change pas grand chose).

Trouves-tu comme moi la gravité du soleil égale à 0,005836 m.s-2 ? Ca me parait beaucoup.

Autre point : si tu trouves une orbite circulaire c'est étrange, tu devrais être à 300 km péri et ~1500 km apo.

Effectivement, j'avais pris une orbite circulaire comme condition initiale.
Pour la force du soleil, pas besoin d'ordinateur, l'expression analytique te donne le résultat:
F=mA
Donc A = G*masse_soleil / d²
A = 6.6742e-11*1.989e30/(1.496e11)² = 0.0059315857 m/s². Soit à peu près ce que tu trouves. Mais je répète que l'important n'est pas cette valeur, par exemple le trou noir central de notre galaxie t'accélère encore plus, et ça n'a aucun impact sur les orbites de nos satellites, pour la simple raison que la Terre et le satellite sont tous les deux accélérés (au premier ordre) de la même grandeur par notre soleil (ou par le centre de notre galaxie). Il est donc essentiel de laisser "tomber" la Terre sur le soleil (sous forme d'orbite ou autre) dans tes simus, sinon effectivement l'effet du soleil devient important et complètement faux. La simu avec une Terre fixe par rapport au soleil mais où le soleil attire quand même le satellite reviendrait à dire que la Terre a un énorme moteur qui la fait accélérer de 0.0059m/s² en permanence en direction opposée au soleil... pas étonnant que tes orbites en soit perturbées dans ce cas !

lionel a écrit:C'est quoi ton programme si c'est pas indiscret ?

C'est pas du tout indiscret puisque ça m'a pris 2 minutes à faire: 3 lignes de scripts matlab te font résoudre toutes les équa diff que tu veux, et par extension les problèmes à N corps (les joies des fonctions "odeXX"). Et 3 lignes pour expliciter les conditions initiales. Si tu veux vraiment le scripter sans l'aides des "ode", coder un Runge-Kutta n'est pas bien compliqué non plus...

Space Opera a écrit:Mais je répète que l'important n'est pas cette valeur, par exemple le trou noir central de notre galaxie t'accélère encore plus, et ça n'a aucun impact sur les orbites de nos satellites, pour la simple raison que la Terre et le satellite sont tous les deux accélérés (au premier ordre) de la même grandeur par notre soleil (ou par le centre de notre galaxie). Il est donc essentiel de laisser "tomber" la Terre sur le soleil (sous forme d'orbite ou autre) dans tes simus, sinon effectivement l'effet du soleil devient important et complètement faux. La simu avec une Terre fixe par rapport au soleil mais où le soleil attire quand même le satellite reviendrait à dire que la Terre a un énorme moteur qui la fait accélérer de 0.0059m/s² en permanence en direction opposée au soleil... pas étonnant que tes orbites en soit perturbées dans ce cas !

Je commence à me rendre compte que j'ai encore certaines lacunes dans la compréhension de la mécanique céleste (et ça tombe bien, puisque j'avais commencé le développement de ce simulateur pour mieux comprendre, là je sens que je vais apprendre un truc, et que le bug c'est MOI :D).

Tu pourrais STP m'expliquer cette histoire de "premier ordre" (ou me trouver un lien qui explique cette notion, j'ai cherché sur Google mais je ne dois pas avoir les bons termes) ?

Ce que je crois comprendre, c'est que dans le système solaire, tous les objets sont attirés "au premier ordre" par le soleil, et que ma boule est avant tout un satellite du soleil : sa rotation autour de la Terre est secondaire. Mais du coup, comment reconnait-on cet ordre ? Quand je regarde les chiffres au niveau de ma boule, je vois que l'attraction du soleil est très très faible par rapport à celle de la Terre... Question subsidiaire : quand je disais que ma boule a une vitesse initiale de 8000 m.s-1, dans mon esprit c'était "par rapport à la Terre", mais en réalité elle est probablement par rapport au soleil, ce qui fait qu'elle avance en fait bien moins vite que la Terre ? Il est possible que mon erreur vienne de là...

(mode réflexion, je pense tout haut) --> car si je dis que ma boule a une vitesse initiale de 8000 m.s-1, abstraction faite de la Terre, ça voudrait dire qu'en fait elle va assez lentement, et que donc elle est sur une orbite elliptique autour du soleil... Il se trouve que quand je lance la simu, la Terre est à côté et donc elle perturbe cette orbite elliptique, récupère la boule mais sur une orbite terrestre pourrie qui finit dans l'eau rapidement ! Je ne suis sûr de rien, il faut que je bouquine un peu et que je fasse quelques tests.

Tu te fais un peu des nœuds au cerveau là. La mécanique céleste est beaucoup plus simple qu'il n'y parait.

Donc pour résumer:
- Si dans tes simus tous les objets ne bougent pas comme ils sont censés bouger avec la gravitation, c'est qu'il y a une force qui l'en empêche. Ainsi, si la Terre n'est pas attirée par le soleil mais que la boule l'est, tu te rends bien compte qu'il y a un problème dans la formulation du problème :)
- Oublie cette histoire de "premier ordre". Pour dire les choses autrement, le soleil est tellement loin que quelle que soit la position de la boule sur son orbite, elle est attirée avec la même accélération (en norme) par le soleil. Et la Terre aussi est attirée avec la même accélération. D'où le fait que l'impact du soleil est négligeable sur l'orbite de ton objet. Pour dire les choses encore autrement, la valeur de l'accélération induite par le soleil n'a que peu d'intérêt: même avec une accélération plus grosse encore que l'attraction de la Terre ça ne changerait pas grand chose, de toute façon la Terre tomberait à la même vitesse que la boule sur le soleil, donc l'effet de ce dernier est presque invisible sur l'orbite de la boule.
- Une petite image pour comprendre ton bug: le fameux principe de la plume et du marteau sur la Lune ou dans un ascenseur qui tombe: on se fiche pas mal de savoir ce qui attire l'ensemble, les 2 objets vus depuis l'ascenseur ont l'air fixes, et leur interaction se fait en oubliant complètement qu'ils tombent quelque part. Et ton erreur qui consiste à modéliser que la Terre ne bouge pas quand elle devrait tomber sur le soleil, ça revient à dire que tu mets une plume (la boule) et un marteau (la Terre) dans un ascenseur qui tombe, la plume tombe normalement mais tu accroches le marteau au plafond de la cage d'ascenseur ! Donc forcément, vue du marteau (la Terre), la plume tombe attirée par le sol... et ne se comporte pas du tout comme elle est censée se comporter, à savoir rester gentiment à côté. Il faut donc absolument faire tout bouger avec la gravité, sinon GROS problème.
- Et de la même façon qu'étudier la chute d'un objet dans un train qui avance à vitesse constante en ligne droite, ça revient à faire la même chose à l'arrêt (quasiment), étudier une orbite en LEO avec ou sans soleil revient aussi (quasiment) à la même chose. Dans un cas, tout le monde tombe et tourne autour du soleil, dans l'autre cas le soleil n'existe pas. Vu de la Terre, ça change quasiment rien.

Bref, je sais pas comment expliquer ça encore autrement, là j'arrive un peu au bout des façons simples que j'ai trouvé.

Tes explications sont très claires,  pas de souci.

Je vais refaire des tests mais je crois que même avec la Terre en mouvement ça ne marchait pas convenablement. Il y a des raccourcis que j'ai pris en programmant qui sont peut-être fautifs, par exemple mon soleil est immobile cad nullement impacté par la Terre, et ma boule de pétanque n'exerce aucune force sur les corps célestes. Je veux avoir les idées claires sur cette notion de vitesse absolue ou relative aussi. Bref, je fais quelques tests et je vous tiens informés.

Space Opera a écrit:Tu te fais un peu des nœuds au cerveau là. La mécanique céleste est beaucoup plus simple qu'il n'y parait.
Donc pour résumer:...

- Et de la même façon qu'étudier la chute d'un objet dans un train qui avance à vitesse constante en ligne droite, ça revient à faire la même chose à l'arrêt (quasiment), étudier une orbite en LEO avec ou sans soleil revient aussi (quasiment) à la même chose. Dans un cas, tout le monde tombe et tourne autour du soleil, dans l'autre cas le soleil n'existe pas. Vu de la Terre, ça change quasiment rien.

Bref, je sais pas comment expliquer ça encore autrement, là j'arrive un peu au bout des façons simples que j'ai trouvé.

Dans ce dernier exemple tu es presque arrivé à la relativité galiléenne. Tu en es presque arrivé à l'exemple de Galilée avec ses papillons qui volent dans la cale d'un bateau
faisant la route "Venise Alep" (oui, je sais Alep n'est pas en bordure de mer). 

En tous cas, Space Opera merci pour ton effort pédagogique.

Astro-notes a écrit:En tous cas, Space Opera merci pour ton effort pédagogique.

Oui merci Space Opera, je crois que tu m'as mis sur la bonne voie. Je n'ai pas eu le temps de retravailler sur le simulateur, mais tes explications ont fait travailler mon cerveau en tâche de fond toute la journée, et je crois que je sais où se trouve mon bug. Je pense qu'il est dû au fait que je modifie la position de la boule relativement à la Terre, et non en absolu (l'absolu étant le soleil dans ma simulation). En clair, si mon idée est juste, je fais subir deux fois l'attraction du soleil à la boule :
1) une première fois quand j'applique la somme des forces exercées par la Terre + le soleil à la position de la boule par rapport à la Terre
2) une seconde fois quand j'applique la force exercée par le soleil à la position de la Terre par rapport au soleil (la boule étant positionnée par rapport à la Terre)

Un peu comme si, au lieu de faire marcher un homme dans un train, je fais marcher un homme dans un train qui est lui-même dans un train. Vu du train qui est le plus à l'extérieur, l'homme va trop vite, c'est ce qui se passe pour ma boule dont l'apogée augmente énormément et le périgée diminue tout autant.

Je testerai cette hypothèse demain et vous tiendrai informés du résultat

Bon courage Thierz.

Souvent, quand ce n'est pas simple, c'est que le problème a été mal posé.

Il n'est absolument pas besoin de considérer les corps comme bougeant les uns par rapport aux autres. Au premier ordre, ça ne change rien à ton problème (le soleil fait en gros 1deg/jour en longitude).
Par ailleurs, quand on fait un propagateur d'orbite, il faut commencer par se poser la question du repère choisi ainsi que des forces prépondérantes.
Vu ton orbite de style LEO, avant de mettre les 9 planètes du système solaire et le drag, il aurait été préférable de mettre en premier les perturbations liées au potentiel terrestre (et spécialement les termes en J2) car même le degré 3 du potentiel a beaucoup plus d'influence que la Lune... Mais les harmoniques sphériques du potentiel terrestre ne sont effectivement pas les plus aisées à prendre en compte...

Pour en revenir au problème de base, le repère choisi est une grosse partie du problème. C'est bien plus simple de prendre un repère inertiel centré sur la Terre. On ne bascule dans un repère héliocentrique qu'une fois qu'on a quitté la sphère d'influence de la Terre. En tout cas, c'est ce qui est fait traditionnellement...
Après, le choix d'une représentation polaire n'est pas pour moi le plus judicieux.
Prends plutôt un classique repère cartésien avec un vecteur d'état de type (Px,Py,Pz,Vx,Vy,Vz): l'accélération liée au potentiel terrestre est alors du genre -muEarth/norm(P)^3 * P (avec P = vecteur position).
Et pour la force perturbatrice du soleil, il suffit d'avoir le vecteur sat-sun (Psat_sun = Psun-Psat) et de faire la même manip. Facile en cartésien. Mais attention! Ici, l'accélération est du genre a = muSun * (Psat_sun/norm(Psat_sun)^3 - Psun/norm(Psun)^3).
C'est l'accélération "différentielle" qu'il faut prendre et dans ce cas, tu t'en fous des mouvements de la Terre vs le Soleil!

L'accélération globale subie par ton SV est alors: a = aEarth + Asun + Amoon + ...

Ensuite:
- dans ton tableau, ça me fait peur quand on parle de "poids" au lieu de force. Par ailleurs, l'accélération liée aux potentiels gravitationnels ne dépend pas de la masse donc autant utiliser directement l'accélération.
- Quel est l'intégrateur que tu as choisi? RK4?

En résumé, mes conseils:
- Mets toi dans un repère inertiel centré Terre et directement en "3D" (genre J2000 ou ITRF). Ca ne coute pas plus cher de se mettre directement en 3D et ça te permettra d'étendre rapidement ton propagateur par la suite. Dans un premier temps, rien ne t'interdit de mettre ta composante Z à zéro.
- Je te conseille de te mettre en cartésien avec un vecteur d'état à intégrer X=(Px,Py,Pz,Vx,Vy,Vz).
- Ne pas oublier de prendre le différentiel d'accélération dans le cas d'un 3ème corps.

Ensuite il suffit d'intégrer dX/dt = fdyn(X,t)

A+

krys a écrit:Souvent, quand ce n'est pas simple, c'est que le problème a été mal posé.

Il n'est absolument pas besoin de considérer les corps comme bougeant les uns par rapport aux autres. Au premier ordre, ça ne change rien à ton problème (le soleil fait en gros 1deg/jour en longitude).
Par ailleurs, quand on fait un propagateur d'orbite, il faut commencer par se poser la question du repère choisi ainsi que des forces prépondérantes.
Vu ton orbite de style LEO, avant de mettre les 9 planètes du système solaire et le drag, il aurait été préférable de mettre en premier les perturbations liées au potentiel terrestre (et spécialement les termes en J2) car même le degré 3 du potentiel a beaucoup plus d'influence que la Lune... Mais les harmoniques sphériques du potentiel terrestre ne sont effectivement pas les plus aisées à prendre en compte...

Pour en revenir au problème de base, le repère choisi est une grosse partie du problème. C'est bien plus simple de prendre un repère inertiel centré sur la Terre. On ne bascule dans un repère héliocentrique qu'une fois qu'on a quitté la sphère d'influence de la Terre. En tout cas, c'est ce qui est fait traditionnellement...
Après, le choix d'une représentation polaire n'est pas pour moi le plus judicieux.
Prends plutôt un classique repère cartésien avec un vecteur d'état de type (Px,Py,Pz,Vx,Vy,Vz): l'accélération liée au potentiel terrestre est alors du genre -muEarth/norm(P)^3 * P (avec P = vecteur position).
Et pour la force perturbatrice du soleil, il suffit d'avoir le vecteur sat-sun (Psat_sun = Psat-Psun) et de faire la même manip. Facile en cartésien. Mais attention! Ici, l'accélération est du genre a = muSun * (Psat_sun/norm(Psat_sun)^3 - Psun/norm(Psun)^3).
C'est l'accélération "différentielle" qu'il faut prendre et dans ce cas, tu t'en fous des mouvements de la Terre vs le Soleil!

L'accélération globale subie par ton SV est alors: a = aEarth + Asun + Amoon + ...

Ensuite:
- dans ton tableau, ça me fait peur quand on parle de "poids" au lieu de force. Par ailleurs, l'accélération liée aux potentiels gravitationnels ne dépend pas de la masse donc autant utiliser directement l'accélération.
- Quel est l'intégrateur que tu as choisi? RK4?

En résumé, mes conseils:
- Mets toi dans un repère inertiel centré Terre et directement en "3D" (genre J2000 ou ITRF). Ca ne coute pas plus cher de se mettre directement en 3D et ça te permettra d'étendre rapidement ton propagateur par la suite. Dans un premier temps, rien ne t'interdit de mettre ta composante Z à zéro.
- Je te conseille de te mettre en cartésien avec un vecteur d'état à intégrer X=(Px,Py,Pz,Vx,Vy,Vz).
- Ne pas oublier de prendre le différentiel d'accélération dans le cas d'un 3ème corps.

Ensuite il suffit d'intégrer dX/dt = fdyn(X,t)

A+

Bonjour krys, je suis archi désolé mais une bonne moitié de ton message est du chinois pour moi (j'ai mis le chinois en jaune dans la citation ci-dessus, pour que tu te rendes compte de tout ce qu'il me reste à découvrir). J'ai réalisé ce simulateur sur la base de mes connaissances en physique et mathématiques, de niveau scolaire, pour justement comprendre comment tout fonctionne, et je progresse en comparant mes résultats avec la réalité (pas très scientifique comme approche, j'en conviens).

Pour l'instant, je réfléchis à comment je vais concilier mon hypothèse d'hier dans mon code actuel, tout en le conservant "élégant" et maintenable, et surtout sans tout chambouler à ce stade. Je sens que je ne suis plus très loin de quelque chose qui ressemble à la réalité.

Merci en tout cas pour tes conseils, je ne doute pas que je me repencherai dessus quand j'aurai un peu plus de "bouteille" :-)

J'ai profité de ma pause méridienne pour faire rapidement le propagateur en question en Matlab.
Sans Matlab, il faudrait que tu fasses l'intégrateur (ici j'utilise par facilité "ode45" de Matlab) mais en codant un simple RK4 (Runge Kutta d'ordre 4 - cf Wikipedia) et en prenant un pas d'une minute, ça doit être assez précis.
Tu verras que la fonction qui calcule les accélérations fait quelques lignes.
La position du soleil est fixe pour l'instant (fonction getPsun) mais il suffit de la modifier pour décrire le mouvement du soleil dans le repère centré Terre.

Je calcule l'orbite avec juste la force centrale puis en y ajoutant la perturbation liée au soleil et ensuite j'affiche l'écart en X/Y.
Après une journée, tu verras que l'écart en position arrive en gros à 30km (au max) dans mon cas. Les oscillations sur l'écart sont liées à l'impact de la perturbation sur l'excentricité et l'argument du périgée. Le demi-grand axe n'est pas impacté puisque la perturbation du soleil est une force conservative.
Le proxy de ma boite m'interdit de mettre les plots obtenus mais si quelqu'un a Matlab, il pourra les mettre :-)

Bye

Code:

function trajectory( )
muEarth = 398600.4415e+9;
muSun = 1.32712440018e+20;

Xini = [ 6378137+300000, 0, 0, 0, 8000, 0 ]';

options = odeset('RelTol', 1e-10, 'AbsTol', 1e-10);
useSun = false;
[ tout, XoutWithoutSun ] = ode45(@fdyn, 0:10:86400, Xini, options);
useSun = true;
[ ~, Xout ] = ode45(@fdyn, 0:10:86400, Xini, options);


figure(1)
subplot(2,1,1)
plot(XoutWithoutSun(:, 1)/1000, XoutWithoutSun(:, 2)/1000);

xlabel('X [ km ]');
ylabel('Y [ km ]');
grid on
axis equal
subplot(2,1,2)
plot(Xout(:, 1)/1000, Xout(:, 2)/1000);
xlabel('X [ km ]');
ylabel('Y [ km ]');
grid on
axis equal

% ecart en position
figure(2)
plot(tout, (XoutWithoutSun(:, 1:3) - Xout(:, 1:3)));
legend('X', 'Y', 'Z');
xlabel('Time [ s ]');
ylabel('Position [ km ]');
grid on

    function Xdot = fdyn(t, X)
       
        % Acceleration from Earth
        aEarth = -muEarth/norm(X(1:3))^3 * X(1:3);
       
        if (useSun)
            % Accel from Sun
            Psun = getPsun(t);
            PsatSun = Psun - X(1:3);
            aSun = muSun * (PsatSun / norm(PsatSun)^3 - Psun / norm(Psun)^3);
        else
            aSun = [ 0; 0; 0 ];
        end
        % total accel
        aTotal = aEarth + aSun;

        Xdot = [ X(4:6); aTotal ];
    end

    function Psun = getPsun(t)
       Psun = [ 150e9; 0; 0 ];
    end
end

Je n'utilise pas Matlab, c'est Space Opera qui en parle dans ses messages ci-dessus. Pour ma part, j'utilise Java, j'ai codé l'intégrateur.

Thierz a écrit:Pour l'instant, je réfléchis à comment je vais concilier mon hypothèse d'hier dans mon code actuel, tout en le conservant "élégant" et maintenable, et surtout sans tout chambouler à ce stade. Je sens que je ne suis plus très loin de quelque chose qui ressemble à la réalité.

Merci en tout cas pour tes conseils, je ne doute pas que je me repencherai dessus quand j'aurai un peu plus de "bouteille" :-)

Personnellement, je te conseille de remettre à plat tout de suite ton problème.
Ce n'est pas quand tu auras pondu des centaines de lignes de code qu'il faudra changer ton fusil d'épaule car tu te rends compte que c'est une galère de prendre en compte une nouvelle perturbation... (comme on dit: "on ne fait pas d'un âne un cheval de course" ;-)).
Donc je te conseille au minimum de changer et de passer au moins en cartésien. Car tu vas vite te prendre la tête à calculer tes perturbations en r "point" et teta "point" alors que c'est en général simple en cartésien (exemple typique: le drag...). Et ça va être d'autant plus galère si tu restes avec ton repère héliocentrique alors que la trajectoire est une orbite autour de la Terre.

Je pense que tu vas galérer en polaire et en plus, ce ne sera pas facilement transposable en sphérique si un jour, tu veux attaquer le problème en 3D.

Thierz a écrit:Je n'utilise pas Matlab, c'est Space Opera qui en parle dans ses messages ci-dessus. Pour ma part, j'utilise Java, j'ai codé l'intégrateur.

Si tu utilises Java et que tu ne veux pas réinventer la roue, tu peux utiliser la librairie Apache Commons-Math qui te permettra de gérer facilement les vecteurs 3D et elle propose toute une série d'intégrateurs en partant de la méthode d'Euler jusqu'à du Dormand-Prince d'ordre 8, voire plus.