Jump to content

A bit of An Orbital Issue


LordVe

Recommended Posts

I seem to be having quite a bit of issue with a Non-Gravitational Orbital Simulation. Well, saying it is a bit is understating it. I can't seem to get the orbit to actually work.

With the Planet having a Semi-Major of 57.909100, a period of 4.092346 degrees/sec, and an Ecc of 0, so that it is a Circle, it goes way out in an elliptical orbit. I am unsure what I am doing wrong.

Here are the Equations I am using:


Finding the Coords:
t = Degrees * Math.Degrees2Radians //Convert to Radians
x = a cos(t) (Radius along the Semi-Major Axis)
z = b sin(t) (Radius along the Semi-Minor Axis)

Calculating the Semi-Minor as:
b = a sqrt(1 - Ecc^2)

Link to comment
Share on other sites

You do have the code for computing semi-minor before you are using it to compute the z coordinate, right? Otherwise, it'd be useful to see a bit more of your actual code.

Also, keep in mind that x = 0, z = 0 is going to be center of ellipse, not focus. In other words, barycenter of your system is not going to be at (0, 0). I don't know which one you are going for.

Link to comment
Share on other sites

Ok, I might be wrong, but it looks to me like OrbitalMechanics.FindSemiMinor() is called when object is created, and only then do you set the SemiMajor axis. So SemiMinor would be computed wrong. You can check that by displaying or printing to console the values of SemiMajor and SemiMnior axis on every tick.

If that's really the problem, the simplest way to fix it would be to create a Body.Initialize() function and call it after you set all the other parameters. So it'd be something like this.

Planet = (GameObject)UnityEngine.Object.Instantiate(sphere);
Planet.name = "Mercury";
Planet.SetActive(true);
transform = Planet.transform;
body = Planet.AddComponent<Body>(); //I'm guessing Body.Start() gets called here.
body.Root = GameObject.Find("Sol").transform;
body.SemiMajor_Axis = 57.909100f;
body.Size = 2439.7f;
body.degrees = 4.092346f;
body.Initialize(); //This call should set up body.SemiMinor_Axis
transform.parent = SolSystem.transform;
Planet.renderer.material.color = Color.red;

I might be wrong, though. I haven't worked with Unity too much. But in either case, taking a look at what the values of SemiMinor and SemiMajor actually are would help you diagnose the problem.

Link to comment
Share on other sites

Are you looking at the mass of the sun? The position of a planet and its speed are not enough to generate an orbit, you also need the mass of the central body to be just right to make a circular orbit. Unless the orbit is hardcoded regardless of the sun's mass, in which case this is not relevant.

Link to comment
Share on other sites

You seem to be conflating the Mean Anomaly, the Eccentric Anomaly, and the True Anomaly. The three are equivalent if your orbit is circular, if your orbit is elliptical, they are only equivalent at periapse (where all three are 0), and apoapse (when all three are equal to pi)

Link to comment
Share on other sites

@K^2: Doing that doesn't make a difference, but you are right, I should be sure to do it after I set the Major axis.

@metaphor: I am attempting to do it on rails like KSP does. Gravity is not a factor.

@maltesh: I may be using the wrong terminology in defining my variables, but I am copying KSP for simplicity, and they call it Mean Anomaly.

For reference, this is where I got my Coordinate equations... I just don't understand why I can't get it to work. I can set both the Major and Minor axis to 3, and it still goes far out of the orbit I intended, and for some reason, never negative on the z axis... I would look into how Unity does the RotateAround method, but since it calls and external dll to do it, I can't.

So, for those following along at home, I am trying to create a function to give me the Coordinates of a body based only, for now, on the Semi-Major axis, the degrees per second, and the Eccentricity of the orbit. I will of course need to add in the Inclination, but that is more equations, and if I can't get it to function with just 2, I can imagine how it will do with more.

Link to comment
Share on other sites

Hm. What about this bit of code in Body.Update()

transform.Translate(Coords);

Does this call translate object to Coords or by Coords? Because if it's the later, it would explain a lot, including why z doesn't go negative.

Link to comment
Share on other sites

I'd start with making sure that your orbit has the star at the focus rather than center. As I mentioned earlier, your equations set the star in the geometrical center of the ellipse. That's not a problem for a circular orbit, since the foci coincide with center for circle, but the moment you add a bit of eccentricity, it will look wrong.

You have equation for ellipse with origin on center:

x = a*cos(t);
z = b*sin(t);

Instead, you want equation where origin is at the focus. This is achieved by simply shifting the ellipse by focal distance.

x = a*e + a*cos(t);
z = b*sin(t);

Inclination (i) is simply a rotation around the Z axis (since your major axis runs along X). In general, a rotation around Z axis looks like this:

x_after_rotation = cos(angle) * x_before_rotation - sin(angle) * y_before_rotation;
y_after_rotation = sin(angle) * x_before_rotation + cos(angle) * y_before_rotation;

But prior to inclination y is zero. So if you apply this formula, the equation for an inclined ellipse looks like this:

x = cos(i) * (a*e + a*cos(t));
y = sin(i) * (a*e + a*cos(t));
z = b*sin(t);

The z coordinate doesn't change since Z is the axis we are rotating around.

In principle, there are two more angles you need to set the orbit. That's the argument of periapsis and longitude of ascending node. But if you are modeling a solar system, where all inclinations are small, you can just add a rotation around the Y axis. If you do so, the full equation takes the following form.

x = cos(f) * cos(i) * (a*e + a*cos(t)) + sin(f) * b*sin(t);;
y = sin(i) * (a*e + a*cos(t));
z = cos(f) * b*sin(t) - sin(f) * cos(i) * (a*e + a*cos(t));

Here f is the angle by which you rotate position of the major axis. For most notable objects in the Solar system, this is more than adequate.

Oh, and maltesh is absolutely correct about the fact that your equations are for eccentric anomaly (variable t in all of the above), not mean anomaly. What that means is that for highly eccentric orbits your objects won't traverse the orbit at quite the right rate. Mean anomaly and eccentric anomaly are closely related, however, so it's not very difficult to recompute one from the other. I can show you how to do that if you'll want to correct for this.

Link to comment
Share on other sites

First off, Thank you.

All I know is that t is called the Parameter... Beyond that... So I assumed it was the same as the Mean Anomaly. So, Calculating it as Degrees/Second * Degrees2Radians won't produce the proper 'orbital speed'? IE a Body with an e of .2 and one of 0 won't reach the intercepts at the same time.

Link to comment
Share on other sites

Parameter is a general term for the independent variable in the parametric equation. The simplest parametric equation for ellipse is written with eccentric anomaly as parameter. That's the equation you used. But it's possible to use a different parameter, including the mean anomaly. You'll just need a different equation. (Or a map from mean anomaly to eccentric anomaly.)

Mean anomaly is defined in such a way that it advances by equal amounts over equal time intervals. In other words, if you want to define the rate at which toe object traverses orbit in terms of degrees/second, you are talking about mean anomaly. So if you still want to use parametric equations based on eccentric anomaly, you need to first convert mean anomaly to eccentric anomaly.

But like I said, the correction is small if the eccentricity is small, so I'd focus on getting everything else to work first.

Link to comment
Share on other sites

I am guessing, but the fact that my venus is going clockwise, that the inclination is in radians. Am I correct?

So, what is the conversion from what I am using to the 'proper' one? If I am going to do this, I would like to do it right. There is a Comet with an Ecc of .7 so, Ecc will be an issue. ^_^

Edited by LordVe
Link to comment
Share on other sites

Yes, all angles are in radians. So if your input is in degrees, don't forget to convert them.

It's easy to go from eccentric anomaly to mean anomaly. m = t - e * sin(t). But you need to go backwards and get t from m. I'm not sure there is an exact solution for that, but there are some iterative methods you can use. Off the top of my head, this should work.

//Given m and e, compute t.
t = 0;
for(j = 0; j < 100; j++)t = m + e * sin(t);

The higher the number in "j < 100" the better the approximation is going to be. There is probably a better way.

Link to comment
Share on other sites

Yes, all angles are in radians. So if your input is in degrees, don't forget to convert them.

It's easy to go from eccentric anomaly to mean anomaly. m = t - e * sin(t). But you need to go backwards and get t from m. I'm not sure there is an exact solution for that, but there are some iterative methods you can use. Off the top of my head, this should work.

//Given m and e, compute t.
t = 0;
for(j = 0; j < 100; j++)t = m + e * sin(t);

The higher the number in "j < 100" the better the approximation is going to be. There is probably a better way.

There is no analytical solution to the equation*, so iterative methods are typical for solving it. The method you've given, which is fixed-point iteration, will converge reasonably quickly (likely 100 iterations will put you far past the precision of most floating point types) if it does converge (which it will so long as e < 1; it can also converge if e >=1, but not all the time). You can immediately knock one iteration off of your loop by starting with t = m rather than 0.

Other common methods for solving this equation is to use Newton's method (which is technically still fixed-point iteration, but we change the form of the function to hopefully get better convergence rates) on f(t) = t - m - e * sin(t), which yields the following scheme:

f'(t) = 1 - e * cos(t)

t[n+1] = t[n] - f(t[n]) / f'(t[n]) = t[n] - (t[n] - M - e * sin(t[n])) / (1 - e * cos(t[n]))

Or, in a slightly more readable code form:

t = m;
for( int j = 0; j < 100; ++j ) t -= ( t - m - e * sin(t) ) / ( 1 - e * cos( t ) );

This will actually converge very quickly, again 100 iterations is likely not necessary. In fact 7 iterations should be sufficient for e < 0.5. If e is near 1 things will get worse, but I never saw it take more than about 35 iterations to get within machine precision. Ideally your loop would look more like this:


t0 = 0;
t1 = m;
epsilon = 1e-15;
j = 0;
while( abs( t1 - t0 ) > epsilon && j < 100 )
{
t0 = t1;
t1 -= ( t1 - m - e * sin(t1) ) / ( 1 - e * cos(t1) );
++j;
}

*NOTE: Strictly speaking there kind of is an analytical solution for it, but you'll get it in terms of a Taylor series. I'm not sure that this gives any computational advantages; the standard iterative approaches may still be quicker to compute to the desired accuracy.

Link to comment
Share on other sites

While I know I have already asked a lot, I wonder if you can also tell me how to find the the rotation on the Y, IE to find f in K^2's equation.

Also, since the equations are for bodies orbiting origin, How do I do a Moon?

Edited by LordVe
Link to comment
Share on other sites

This thread is quite old. Please consider starting a new thread rather than reviving this one.

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

×
×
  • Create New...