Bonjour,
je développe en C un moteur de jeu entièrement de zéro, je souhaite n'utiliser aucune librairie existante.
Du coup j'en arrive à la physique, j'utilise la méthode Runge-Kutta d'ordre 4 (RK4) et elle me donne de très bon résultat pour beaucoup d'utilisations.
Mais j'ai bien envie de faire une démo spatiale, et j'aimerai bien pouvoir intégrer la gravitation newtonienne à mon RK4 plutôt que de faire de la simple et fausse rotation car ça me permettrait par exemple de gérer des collisions entre objets en orbite.
Le problème : la simulation d'orbite ne fonctionne pas très bien, pour le moment je teste avec un satellite seul, il tourne bien autour de ma planète mais son orbite s'agrandit peu à peu et il fini par partir dans l'espace. Bien entendu j'utilise des types 'double' dans mon code pour éviter des erreurs d'arrondi.
Je ne sais pas vraiment si je dois mettre mon code ici puisque c'est tout langages confondus. L'intégrateur RK est tout ce qu'il y a des plus classique, voilà ce qu'il donne en pseudo-code :
Voilà comment je calcule la vitesse initiale (d'après Wikipédia : http://fr.wikipedia.org/wiki/Vitesse_orbitale) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 deriv Evaluate0(local state) { deriv output output.force = state.force output.velocity = state.velocity return output } deriv Evaluate(local state, t, dt, d) { state.position = state.position + d.velocity * dt state.momentum = state.momentum + d.force * dt state.velocity = state.momentum * (1.0 / state->mass) deriv output output.force = state.force output.velocity = state.velocity return output } void Integrate(ref state) { a = Evaluate0(*state) b = Evaluate(*state, t, dt*0.5f, a) c = Evaluate(*state, t, dt*0.5f, b) d = Evaluate(*state, t, dt, c) state->position = state->position + (a.velocity + (b.velocity + c.velocity) * 2.0 + d.velocity) * 1.0f/6.0 * dt state->momentum = state->momentum + (a.force + (b.force + c.force) * 2.0 + d.force) * 1.0f/6.0 * dt state->velocity = state->momentum * (1.0 / state->mass) } void ApplyForce(ref state, vec) { state->force = state->force + vec } void ApplyGravity(ref s1, ref s2) { const G = 6.67234E-11 d = dist(s1->position, s2->position) P = G * (s1->mass * s2->mass) / (d * d) force = P * (s2->position - s1->position) / d ApplyForce(s1, force) ApplyForce(s2, -force) }
Par exemple pour un rayon orbital de 20 000 km, j'ai une vitesse initiale v0 de 4464.18 m/s
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2 satellite.momentum.y = sqrt(G * earth.mass / R) * sat.mass // sqrt étant la racine carré, G la constante gravitationnelle, et R le rayon orbitale du satellite
Je suppose être dans le vrai si je dis que c'est un problème d'arrondi au niveau des calculs ? Comme faire pour éviter ça sans perdre en performance, y a-t-il des moyens de simplifier les calculs ?
Voilà, si quelqu'un aurait une idée ou un lien magique à propos de ce genre de chose, je n'ai rien trouvé sur google à propos du newtonien avec RK4..
Merci
Partager