I have a working Keplerian orbit simulator, as long as I don’t change any of the orbital values at runtime it all works fine. However, I am trying to add thrust while the simulation is running. Right now I can convert between position vectors (unity world space coordinates) and Keplerian elements except for the value called “True Anomaly.” This value is basically the angle from the lowest point in the orbit (Periapsis) to the current position of the object around the ellipse of the orbit. Right now I calculate that value based on the total time elapsed and the total time it takes to go around the orbit (As Kepler’s law states that an object will move the same amount of degrees in the same amount of time around its orbit). But that calculation is useless once I start to change to a new orbit.
SO! I’m sure there is a way to do this, but I don’t quite have the math background to figure it out.
I have the positional vectors of the orbiting object (ox, oy, oz,) the lowest point of the new orbit (px, py, pz), the highest point of the orbit (ax, ay, az) and the position of the object that it is orbiting (fx, fy, fz) and also every Keplerian element except true anomaly and epoch. I’m very uneducated on ellipse geometry especially in 3 dimensions and I’m not sure that would be the way to go about this. Any insight would be appreciated!
I don’t have specific answer, but you may find some guidance by searching for lectures by the Kerbal Space Program devs, who solved this exact problem. I’m certain they have talked about this when discussing their “rails” physics system.
It’s been a while since I did orbital mechanics. I may be off base here.
Keplers laws are a special case of Newton’s laws. To the point that Newton derived his laws using Keplers laws and observational data as a basis.
However Newtons laws can take into account outside forces acting on the object. Which is when an object is changing orbit such as through firing rockets, ect.
Not quite. Keplers 2nd Law states that a line through the orbiting body to the sun (or other mass the object is orbiting) will sweep the same area in the same unit of time. In layman terms this means the further away the orbiting object is from the sun, the slower it moves.
I looked into it a bit more (and deleted my previous post - sorry I completely misunderstood you). So presumably you are looking at :
true anomaly = arccos( e•r / |e| |r| )
where:
e is the eccentricity vector,
r is the orbital position vector (segment from Focus to Orbit) of the orbiting body.
You’ll probably want to find some Kepler Solver code to find e, and I found this:
pi = Math.PI;
function EccAnom(ec,am,dp,_maxiter) {
// ec=eccentricity, am=mean anomaly,
// dp=number of decimal places
pi=Math.PI;
i=0;
delta=Math.pow(10,-dp);
var E, F;
// some attempt to optimize prediction
if (ec<0.8) {
E=am;
} else {
E= am + Math.sin(am);
}
F = E - ec*Math.sin(E) - am;
while ((Math.abs(F)>delta) && (i<_maxiter)) {
E = E - F/(1.0-(ec* Math.cos(E) ));
F = E - ec * Math.sin(E) - am;
i = i + 1;
}
return Math.round(E*Math.pow(10,dp))/Math.pow(10,dp);
}
function TrueAnom(ec,E,dp) {
S=Math.sin(E);
C=Math.cos(E);
fak=Math.sqrt(1.0-ec^2);
phi = 2.0 * Math.atan(Math.sqrt((1.0+ec)/(1.0-ec))*Math.tan(E/2.0));
return Math.round(phi*Math.pow(10,dp))/Math.pow(10,dp);
}
function MeanAnom(time,_period) {
curr_frame = timeToFrames(time);
if (curr_frame <= _period) {
frames_done = curr_frame;
if (frames_done < 1) frames_done = 1;
} else {
frames_done = curr_frame % _period;
}
_fractime = (frames_done * 1.0 ) / _period;
mean_temp = (2.0*Math.PI) * (-1.0 * _fractime);
return mean_temp;
}
//==============================
// a=semimajor axis, ec=eccentricity, E=eccentric anomaly
// delta = delta digits to exit, period = per., in frames
//----------------------------------------------------------
_eccen = 0.9;
_delta = 14;
_maxiter = 1000;
_period = 300;
_semi_a = 70.0;
_semi_b = _semi_a * Math.sqrt(1.0-_eccen^2);
_meananom = MeanAnom(time,_period);
_eccentricanomaly = EccAnom(_eccen,_meananom,_delta,_maxiter);
_trueanomaly = TrueAnom(_eccen,_eccentricanomaly,_delta);
r = _semi_a * (1.0 - _eccen^2) / (1.0 + (_eccen*Math.cos(_trueanomaly)));
x = r * Math.cos(_trueanomaly);
y = r * Math.sin(_trueanomaly);
_foc=_semi_a*_eccen;
[1460+x+_foc,540+y];
This is exactly what I was looking for! I had read up on the eccentric anomaly but didn’t quite get how to do the “guess and check” style prediction that it required, Ill try and get this system working with my current setup of existing values! thank you soo much!
EDIT:
Unfortunately, after getting a good look into the code an issue still exists
Here it requires the number of “Frames” since start, the problem I have specifically is the orbit itself has changed so the total time the program has been running (“curr_frame = timeToFrames(time);”) doesn’t work quite right, I would have to calculate somehow the amount of time that would have passed to get the object to this position on the new orbit, and that seems rather involved and difficult. My gut reaction is this will be better solved with geometry but Im not quite sure.
Geometry could work. In the picture you’ve drawn there is a fairly simple way to solve for theta in two dimesnsions provided you know the relative positions of Focus, Center, and Object in orbit (which it seems you do). I’m not sure if that’s exactly what you’re looking for, but I’ll write up a solution later and if that seems correct I can probably expand it to 3-dimensions.
Ok, so the method to find theta as shown in your graphic should be as follows:
Find the center of the ellipse by finding the midpoint of the Apoapsis and Periapsis (just add the vectors and divide the resulting x,y,z by 2)
You’ll then need to solve the triangle created by connecting Center-Focus-Object (C-F-O). For this you’ll calculate the norm of each resulting vector, |CF| ; |FO| ; |OC|. Each vector is calculated by subtracting the position vectors i.e. CF = F.position - C.position. Then you should be able to just use Vector3.magnitude.
The important angle phi (see graphic below) is derived using the law of cosines:
Cos(Φ) = |CF|² + |OC|² - |FO|² / 2 |CF| |OC|
take the arccos to find Φ and add it to pi and you should have the theta of your original graphic!
