I would like to write a program able to simulate a keplerian orbite given an initial position and velocity (r0, v0) and able to manage accelerations in different directions.
To simulate the orbit, I use the expressions given here and manage particular orbits (without inclination, circular ...). The algorithm is quite simple :
Given r0 and v0 I calculate the orbital elements (a, e, i, Ω, 𝜔, 𝜈) and the coordinates of the satellite in its orbital plan following those formulas : x_orb = a * (cos(E) - e) | y_orb = a*sqrt(1-e²)*sin(E)
I calculate the mean anomaly M, then the eccentric anomaly E and finally the true anomaly 𝜈. Given 𝜈, I obtain the coordinates in the orbital plan given the equations above. Finally, I obtain the coordinate r and v by a multiplication with 3 rotation matrices given here (formulas (12), (13), (14))
All orbits are perfectly simulated.
Problems come when I try to apply a ΔV to my satellite : during the main loop, I just modify the velocity and recalculate all orbital elements as if the actual r and v were initial position and velocity. Then, the upper algorithm is followed again.
Maybe I'm missing something silly, but it doesn't work at all : when I recalculate the orbital elements and relaunch the algorithm, my satellite abruptly changes its position and follow a new orbit.
Is my way to "simulate" the acceleration bad or am I just missing something ?