计算牛顿重力模拟中的2D轨道路径

发布于 2025-01-21 07:17:32 字数 1061 浏览 4 评论 0原文

我目前正在研究一个项目（娱乐），以模拟牛顿重力，并希望能够可视化围绕单个吸引子绕的物体的未来路径。在任何给定的仿真步骤中，我都知道

对象相对于吸引子的中心的当前位置是2D向量作为
对象的当前速度作为2D矢量，
吸引子的质量
是使用这些参数的对象的质量，

我希望能够能够在模拟中最多可以预测n个步骤，以便我可以在轨道沿每个点之间画一条线。

我的模拟中的位置和速度以m和m/s进行测量，而质量是无单位的。重力常数设置为1，每个物体的质量为1000，对于吸引子，对象为10。

模拟本身产生了我想要的结果，现在我只想能够预测未来的道路。

在自己的研究中，我发现我需要使用 keplerian elements 在某种程度上。各种示例，问题等...我在stackexchange上发现，其他地方并没有提供足够的解释以使我能够在模拟中使用它，或者它们特定于三维几何形状并在我尝试计算时提供不正确的输出他们（nan和无限）。

例如，在某些地方，我读到半肌轴和总轨道能量应该保持恒定，并且从这些轴上，我可以获得偏心率和各种其他属性，但是当使用上述提及的任何给定步骤计算它们时我可以访问的属性，它们变化很大或产生如此荒谬或低的数字，以至于它们本质上是没有用的。我解决了最后一个问题，并具有一个半高轴值，在模拟的上下文中有意义，但是随着步骤的变化又有变化。

例如：

初始条件

吸引子

质量= 1000
位置= 0,0,0

object

质量= 10
位置= 5，0
速度= 0，5

这会产生椭圆形轨道。半肌轴应保持恒定，因为没有其他力在对象上作用，但是，半高轴的值与2-4米不同。

原文

I am currently working on a project (for fun) for simulating Newtonian gravity and want to be able to visualize the future path of an object orbiting around a single attractor. At any given simulation step I know

The object's current position relative to the center of the attractor as a 2D Vector
The object's current velocity as a 2D Vector
The mass of the attractor
The mass of the object

Using these parameters I would like to be able to predict up to N steps ahead in the simulation so that I can draw a line between each point along the orbit.

The Position and Velocity in my simulation are measured in m and m/s while the mass is unitless. The Gravitational Constants is set to 1 and the masses for each object is 1000, for the attractor, and 10 for the object.

The simulation itself produces the result that I want, now I just want to be able to predict the future path.

During my own research, I have found that I need to use the Keplerian Elements in some way. The various examples, questions etc... I have found on stackexchange and elsewhere do not provide sufficient explanation for me to be able to work it into my simulation or they are specific for 3-dimensional geometry and provide incorrect outputs when I attempt to calculate them (NaN and Infinities).

For example, in a few places I have read that the semi-major axis and the total orbital energy should remain constant and from these I can get the eccentricity and various other properties, but when calculating them at any given step using the above-mentioned properties I have access to, they vary wildly or produce numbers that are so ridiculously high or low that they are essentially useless. I resolved the last problem and have a semi-major axis value that makes sense in the context of my simulation, but again it varies from step to step.

As an example:

Initial conditions

Attractor

Mass = 1000
Position = 0,0

Object

Mass = 10
Position = 5, 0
Velocity = 0, 5

This produces an elliptical orbit. The semi-major axis should remain constant as no other forces are acting on the object, however, the value of the semi-major axis varies from 2-4 meters.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

西瑶 2025-01-28 07:17:32

continue

I have figured it out and to any future visitors to this question I am (hopefully) going to make your life much, much easier.

If you are like me, fairly lay in terms of astrodynamics then the sheer number of equations, their very technical explanations, and usage is a little daunting. With a little bit of effort, I now know how to iterate equations to be able to calculate future orbital positions.

What we will need are the velocity and positional coordinates of the object you want to predict as well as the mass of the orbital body and the mass of the object itself.

These need to be in m and m/s for no other reason than for this implementation they are assumed to be m and m/s but the mass can be unitless (at first).

Terms:

Attractor - The nearest massive body that is attracting the Object.

Object - The orbital body that is being attracted to the Attractor.

Step 1: Calculate the True Mass of your objects

True mass is defined as the KG mass of your Attractor and Object. To find this value you need to convert your unitless mass, using the gravitational constant, into kilograms.

If you don't know how to do this, first calculate the ratio of unitless mass to the gravitational constant in your simulation and use the value of that ratio to find true mass using the real gravitational constant.

(I cannot post images or else I would include the equations)

After this point when I refer to the gravitational constant, I will be using our real universe's value.

Step 2: Calculate μ

μ is a ratio of masses and the gravitational constant. If the mass of the attractor is large enough μ = GM, where G is the gravitational constant and M is the mass of the attractor. For my own implementation I found that μ can be calculated as:

μ = G ( M * M * M) / ((M+m)*(M+m))

Where:

M is the mass of the Attractor

m is the mass of the Object

G is the gravitational constant.

It is very very important that this value is correct. If what you are rendering looks weird, it is probably because this value was not calculated correctly and it would be the first thing I check. In my initial attempts, this was off by an order of magnitude, it should have been around 100 but it was outputting as 1000 because the mass unit conversion value I calculated was 10x larger than it should have been.

Step 3: Calculate the Orbital Elements

First up is the Semi-Major Axis. This is the radius from the furthest point in the orbit to the center (Aprox). It is represented with the symbol a and calculated as:

a = -(μ * len(pos)) / (len(pos) * len(vel) * len(vel) - (2 * μ))

Man I wish I could post images so these equations looked better, sorry!

Next up is the period of the orbit. This is used later so that we can step through the orbit incrementally. It is represented by the symbol T and its formula is:

T = 2 * π * sqrt((a * a * a) / μ)

Next, we need angular momentum. Now, I know that the post says 2D but in my own code, I am using Vector3's to store my information in case I decide to make it 3D later. This is important because normally to get the angular momentum you would take the cross product of the velocity and position. Since we are in 2D there is no true cross product for vectors. Instead, you can calculate the momentum as:

p = pos.x * vel.y - pos.y * vel.x

Now onto orbital energy, represented by the symbol ε. We use this to find the eccentricity of the orbit.

ε = (len(vel) * len(vel)) / 2 - (μ / len(pos))

After that we can get eccentricity, represented with the symbol e:

e = sqrt(1 + ((2 * ε * p * p) / (μ * μ)))

The second to last thing we need to calculate is the eccentricity vector, represented with ev:

ev = ((((len(vel) * len(vel)) / μ) - (1 / len(pos))) * pos) - ((dot(pos, vel) / μ) * vel)

I know, there are a lot of brackets but I want to make it clear what goes where. Blame StackOverflow for not letting me post equation images.

The final value is the argument of periapsis. Because we are in 2D we need to use this formula to calculate it:

w = atan2(ev.y, ev.x)

we will use this to orient our orbital path correctly.

Step 4: Keplers Equation

Now is where we get to predict the orbital positions using all the information we calculated above.

First, we need n. n is the radians / second that the Object is traveling around the Attractor. It is fairly straightforward to calculate, it is 2π/T where T is the period. This value is used when we calculate the mean anomaly in the predictive step.

Now we get into the loop. You can set the number of steps you want to integrate over to whatever you want, the higher the number the closer each point along the orbit will be. In my case, I have a step of 100, so in my for loop declaration, I will loop 100 times before exiting. In addition, you will want to define a max tries value because we need to do an integration (don't worry, not as scary as it sounds) to be able to get the correct Eccentric Anomaly.

Now that we are in the loop, we need to calculate the current mean anomaly:

ma = n * (T * (i / steps))

Where:

T is the period of the orbit in seconds

i is the current iterative step

n is speed from earlier (radians/second)

Now we do the integration to solve Kepler's Equation. This can be done using Newton's method. Basically, we take Kepler's Equation and its derivative and iterate over them until the delta value between the two of them is arbitrarily low.

This is why we need max tries, it is possible that this while loop iterates forever and never converges on a value in which case we need to bail out so we don't get stuck.

//Eccentric Anomaly
var ea = ma;
while(tries < maxTries)
{
    tries++;
    //Keplers equation
    var dx = ea - e * sin(ea) - ma;
    //Derivative of Keplers equation
    var dy = 1 - e * cos(ea);
    var dt = dx / dy;
    ea -= dt;
    if(abs(dt) < 10^-8) break;
}

See, the integration wasn't so hard. This loop is an implementation of Newton's method and will oscillate toward to correct value. It should be noted that his method is quick and dirty and produces okay results. But this is the only method I understand so this is why I used it.

Step 5: Calculating Orbital Position

Next up is where we calculate the orbital position. This is simple and I stole it directly from this post.

p = a * (cos(ea) - e)

q = a * sin(ea) * sqrt(1 - (e * e))

You can plot these two values as X and Y and you will get an orbit that looks correct in some cases, but we need to rotate it using the argument of periapsis we calculated earlier. Again, stolen from this post.

x = cos(w) * p - sin(w) * q

y = sin(w) * p - cos(w) * q

And that's it. You're done. You have the X and Y coordinate at the current step and you can plot it on the screen however you like. As the for loop continues to iterate, you will get the full orbital path from start to finish.

I hope this essay helps someone out there!

All the equations I used can be found on Wikipedia just by going to this page and navigating to the corresponding element.

回复收藏 0 原文

~没有更多了~