就尺寸和质量而言，是否可以进行逼真的n体太阳系仿真? [英] Is it possible to make realistic n-body solar system simulation in matter of size and mass?

问题描述

重要说明:该问题与"PhysX"有 完全没有 关系，PhysX是一种计算机游戏物理系统(可用于街机游戏中的物理学，例如球类游戏等)； PhysX是Unity3D和其他游戏引擎的内置系统. PhysX与此处完全无关.

/////////////////////UPDATE(从下往上阅读)//////////////////////

我一直在记录值并搜索确切的问题在哪里，我想我找到了它. 我的代码中有类似的内容

Velocity += Acceleration * Time.deltaTime;
position += Velocity * Time.deltaTime;

现在的加速度约为0,0000000000000009 ..随着代码的流动，速度会按照预期的速度增加，浮子没有问题. 但是在开始时，地球的初始位置是(0,0,23500f)您可以在我最后给出的图表中看到这一点.

现在，当我将速度* timedelta(此时约为0,00000000000000005)添加到23500的位置时，基本上不会添加它.位置仍然是(0，0，23500)而不是(0,0，23500.00000000000005)，因此地球不会移动，因此加速度不会改变.

如果我将地球的初始位置设置为0,0,0并仍然将加速度设置为0.0000000000000000009以假设其位置为(0,0,23500) 然后，它添加"速度* timedelta. 它变为(0,0,000000000000000000005)并保持增加.当float为0时，添加如此小的值没有问题.但是，如果浮点数约为23500，则不会将小数值相加.

我不知道这到底是统一的问题还是c#的浮动.

这就是为什么我不能使它以较小的值工作.如果我能克服这个问题，我的问题将得到解决.

///////////////////////////////////////////////////////////////////////////////

我一直在开发n体物理学来模拟我们的太阳系，所以我一直在收集数据以使其尽可能逼真.但是数据大小有问题.我搜寻了互联网的每一个微小部分，却找不到一个能解释人们如何克服这一点的单一解释. (如果他们这样)，所以我在这里尝试射击.

因此，为了保持行星之间的距离，半径和质量"之比固定，我创建了一个Excel文件来计算所有数据. (因为为什么有人会把如果地球上有那个"半径图"放到地球的质量上，这是为什么呢?) 我将把ss作为附件.它从根本上标准化"或换句话说缩放"行星的每个属性到给定的参考.在这种情况下，我将其称为地球半径".

我团结一致地工作，而且您知道，不能统一使用太大"或太小的"值.所以我不得不将太阳系缩小，很多！"

因此，我使用牛顿万有引力定律，即F = GMm/r ^ 2， 为了简单起见，我直接针对所有其他物体的给定物体计算a = GM/r ^ 2.

因此，地球朝着太阳"的重力加速度的实际值大约为0,000006 km/s ^ 2，这即使是一个很小的值，也可以统一使用，但它可以起作用.但是，要获得此值，1我需要将地球的半径(比例)设置为6371单位，将太阳的比例设置为696342 !！

所以我说，让地球半径为1，以统一单位为单位. 所以，当半径改变时，一切都改变了，质量，距离... 我保持了行星的密度，并从具有新半径的新体积中计算了质量. 所有计算都在附件中.

问题是，当我将地球半径设为1时，对太阳的重力加速度变为0,0000000000009 真是小得可笑.当然，Unity不能使用该值.

因此，如果我改为增加地球的半径，那么太阳的质量和半径会变得非常大，然后又无法使用.

我不知道其他人如何解决这个问题，他们如何解决这个问题，但是正如我从这里看到的那样，似乎无法对太阳系进行逼真的n体模拟. (至少一致)

所以我需要10个代表才能发布图像-_-，我将改为提供链接. http://berkaydursun.com/solar_system_simulator/data.PNG 还有一个目录是具有n体计算但具有UNREALISTIC值的工作中的实验性太阳系模拟. 它工作得很好，甚至看起来也接近真实，但是不，它没有正确的比率^^ 如果您希望 http://berkaydursun.com/solar_system_simulator/

，可以在此处进行测试

哇，我几乎每个段落都以"So" ^^

开头

解决方案

我也对sol系统进行了编程，所以这里是我的见解:

1. 渲染

   我使用 OpenGL 和 1:1 缩放.所有单位均在 SI 中，因此 [m，s，kg，...] .该问题开始于 Z缓冲区.通常的 Z缓冲区位宽为16/24/32 bit，与您所需的位置相去甚远.我正在从 0.1m到1000 AU 进行渲染，那么该如何克服呢?

   我确实通过同时结合 Z-sorting 和 Z-buffering 的3个平截头进行渲染来管理它(由于透明的环，Z排序是必需的...并且其他效果).所以首先我将最远的部分渲染到zfar=1000AU.天穹投影到z=750AU距离处，然后清除 Z缓冲区，并渲染到zfar=0.1AU的对象.然后再次清除 Z缓冲区，并渲染关闭的对象，直到zfar=100000 m.

   要完成这项工作，您必须拥有尽可能精确的投影矩阵. gluPerspective的cotangens不精确，因此它需要修复相关的元素(花了我很长时间才能发现). Z near值取决于 Z缓冲区的位宽.如果编码正确，则即使使用变焦10000x，也能很好地工作.我从矿井主视图实时地将该程序用作矿井望远镜:)的对象导航/搜索器.我结合了3D星星，天体，船，真实地面(通过DTM和卫星纹理).它甚至可以输出红蓝浮雕:).可以从表面，大气，空间...(不仅限于地球)进行渲染.没有其他第三方库，则使用OpenGL.这是它的样子:

   如您所见，它可以在任何海拔高度或放大情况下正常工作，就像大气散射着色器

2. 模拟

   我没有使用 n体重力模拟，因为为此，您需要很多非常难以获取的数据(以所需的精度几乎是不可能的).必须非常精确地进行计算.

   我改用 Kepler方程，因此请参见以下内容:

   • 求解开普勒方程
   • C ++实现.

   如果您仍要使用重力模型，请使用 NASA 中的 JPL地域.我认为他们那里也有C/C ++的源代码，但是他们使用的地图与我的地图不兼容，因此对我来说是不可用的.

   通常，开普勒方程具有较大的误差，但不会随时间增加太多.重力模型更精确，但其误差会随着时间而增加，您需要不断更新天体数据以使其起作用...

[edit1]集成精度

您当前的实现方式如下:

 // object variables
double  acc[3],vel[3],pos[3];
// timer iteration
double dt=timer.interval;
for (int i=0;i<3;i++)
 {
 vel[i]+=acc[i]*dt;
 pos[i]+=vel[i]*dt;
 }
 

问题是，当您添加非常小的价值和非常大的价值时，它们就是 在添加之前移至相同的指数，这将舍入大量数据 为了避免这种情况，只需将其更改为:

 // object variables
double          vel0[3],pos0[3]; // low
double          vel1[3],pos1[3]; // high
double  acc [3],vel [3],pos [3]; // full
// timer iteration
double dt =timer.interval;
double max=10.0; // precision range constant
for (int i=0;i<3;i++)
 {
 vel0[i]+=acc[i]*dt; if (fabs(vel0[i]>=max)) { vel1[i]+=vel0[i]; vel0[i]=0.0; } vel[i]=vel0[i]+vel1[i];
 pos0[i]+=vel[i]*dt; if (fabs(pos0[i]>=max)) { pos1[i]+=pos0[i]; pos0[i]=0.0; } pos[i]=pos0[i]+pos1[i];
 }
 

现在xxx0已集成到max，整个内容已添加到xxx1

四舍五入仍然存在，但不再累积.您必须选择max值，以确保积分本身是安全的，并且添加的xxx0+xxx1也必须是安全的.因此，如果一次拆分的数字相差太大，则拆分两次或更多次...

• 喜欢:xxx0+=yyy*dt; if (fabs(xxx0>max0))... if (fabs(xxx1>max1))...

[Edit2]星标

• 星星的颜色我见过的最好的星星可视化
• 星标B-V颜色索引到明显的RGB颜色所有星标目录都使用B-V索引
• 使用Stelar目录，还有星标交叉引用链接
• Skybox:组合不同的恒星数据

[Edit3]进一步提高Newton D'ALembert集成精度

迭代积分的基本问题是，由于在积分步骤dt期间位置发生变化，因此基于当前身体位置施加基于重力的加速度将导致更大的轨道，这在天真的积分中没有考虑.为了解决这个问题，请看这张图片:

假设我们的身体处于圆形轨道并且处于0度位置.我没有使用基于当前位置的加速度方向，而是使用了0.5*dt之后的位置.这样就增加了一点点加速度，从而导致更高的精度(与开普勒轨道相对应).通过此调整，我得以成功地从开普勒轨道转换为牛顿D'Alembert用于2体系统. (下一步是对n体执行此操作).只有来自不受潮汐影响和/或卫星影响的2个身体系统才可能与来自太阳系的真实数据进行粗略的关联.要构造自己的虚构数据，可以使用开普勒圆轨道和均衡重力的等向力:

G = 6.67384e-11;
v = sqrt(G*M/a);                           // orbital speed
T = sqrt((4.0*M_PI*M_PI*a*a*a)/(G*(m+M))); // orbital period

其中，a是圆形轨道半径，m是体重，M是焦点体重(太阳).为了使精度保持在可接受的公差范围内(对我来说)，积分步骤dt应为:

dt = 0.000001*T

因此要放置新的测试主体，只需将其放在:

pos = (a,0,0)
vel = (0,sqrt(G*M/a),0)

主焦点(太阳)在:

pos = (0,0,0)
vel = (0,0,0)

这会将您的身体置于圆形轨道上，以便您可以比较开普勒与牛顿D'Alembert来评估模拟的精度.

Important note: this question has utterly no relation to "PhysX", which is a computer-game-physics system (useful for the physics in arcade games such as ball games, etc); PhysX is a system built-in to Unity3D and other game engines; PhysX is totally irrelevant here.

//////////////////// UPDATE (read bottom first) /////////////////////

I have been logging the values and searching where the exact problem is, and i think i found it. I have something like this in my code

Velocity += Acceleration * Time.deltaTime;
position += Velocity * Time.deltaTime;

The acceleration is something like 0,0000000000000009.. right now. As the code flows, velocity increases as it should be, no problem with the float. But in the beggining, the initial position of earth is (0,0,23500f) You can see this in the chart i gave at the end.

Well now when i add speed * timedelta (which is something like 0,00000000000000005 at this point) to position which is 23500, it basicly doesn't adds it. position is still (0, 0, 23500) not something like (0,0, 23500.00000000000005), thus the earth doesn't move, thus the acceleration doesn't change.

If i set the initial position of earth to 0,0,0 and still, set the acceleration to 0.0000000000000000009 to assuume it's position is (0,0,23500) It then "ADDS" the velocity * timedelta. It becomes something like (0,0,000000000000000000005) and keep increases. When the float is 0, there is no problem with adding such small value. But if the float is something like 23500, then it doesn't adds up the small values.

I don't know if it's exactly unity's problem or c#'s float.

And that's why i can't make it work with small values. If i can overcome this, my problem will be solved.

///////////////////////////////////////////////////////////////////////////////

i have been develeping n-body phyics to simulate our solar system, so i have been gathering data around to make it as realistic as possible. But there is a problem with the data size. I searched every tiny bit of internet and i couldn't find a single explanation how people overcomes this. (If they so) So i trying my shot here.

So, to keep the ratio of distance, radius and "mass" between planets fixed, i created an excel file to calculate all the datas. (Because why the heck would someone put "what would be the earth's mass if it had "that" radius chart" on the internet?) I will give the ss as attachment. It basicly "normalizes" or in other words "scales" every property of a planet to a given reference. In this case, i took the reference as "earth's radius."

I work in unity, and you know, you can't work with "too big" or "too small" values in unity. So i had to scale the solar system down, "a lot!"

So i use the Newton's law of universal gravitation, which is F = GMm/r^2, to make it simple, i am directly calculating the a = GM/r^2, for a given body from all other bodies.

So, the real value of earth's gravitational acceleration "towards sun" is roughly 0,000006 km/s^2, which is even incredibly small value to work with in unity, but it could work. Yet, to get this value,1 i need to set earth's radius (scale) to 6371 unit, and sun to scale of 696,342!, which is TOO big to render it in unity.

So i said, let the earth's radius be 1, in unity units. So, when the radius changes, everything changes, the mass, the distance... I kept the density of the planet and calculate the mass from the new volume with the new radius. All the calculations are in the attachment.

So the thing is, when i take the earth's radius as 1, the gravitational accelaration towards sun becomes is something like 0,0000000000009 which is ridiculously small. And of course Unity doesn't work with that value.

So, if i increase the earth's radius instead, then the mass and radius of the Sun gets ridiculously big and then again, i can't work with it.

I don't know how other people fixed this, what they did to overcome this problem but as i see from here, it looks impossible to make realistic n-body simulation of solar system. (in unity atleast)

So i need to have 10 rep to post images -_-, i will give link instead. http://berkaydursun.com/solar_system_simulator/data.PNG Also one directory up is the working experimental solar system simulation with n-body calculations but with UNREALISTIC values. It works quite well, and it even looks somehow close to real, but no, it doesn't have the right ratios ^^ You can test it here if you wish http://berkaydursun.com/solar_system_simulator/

Edit: WoW I almost started every paragraph with "So" ^^

解决方案

I did program sol system simulation too so here are mine insights:

1. rendering

   I use OpenGL with 1:1 scaling. All units are in SI so [m,s,kg,...]. The problem gets started with Z-buffer. The usual Z-buffer bit wide is 16/24/32 bit which is nowhere near what you need. I am rendering from 0.1m up to 1000 AU so how to overcome this?

   I did manage it by rendering with 3 frustrums at once combining Z-sorting and Z-buffering (Z-sort is necessary because of transparent rings ...and other effects). So first I render most distant parts up to zfar=1000AU. Sky dome is projected at z=750AU distance then clear the Z-buffer and render objects up to zfar=0.1AU. Then clear the Z-buffer again and render close objects up to zfar=100000 m.

   To get this work you have to have as precise projection matrix as you can. The gluPerspective has unprecise cotangens so it needs to repair concerned elements (get me a long time to spot that). Z near value is dependent on the Z-buffer bit width. When coded properly then this works good even with zoom 10000x. I use this program as navigation/searcher of objects for mine telescope :) in real time from mine home view. I combine 3D stars, astro bodies, ships, real ground (via DTM and satellite texture). Its capable even of red-cyan anaglyph output :). Can render from surface,atmosphere,space ... (not just locked to Earth). No other 3th party lib then OpenGL is used. Here is how it looks like:

   As you can see it works fine on any altitude or zoom the atmosphere is done like this atmosphere scattering shader

2. simulation

   I am not using n-body gravity simulation because for that you need a lot of data that is very very hard to get (and almost impossible in desired precision). The computations must be done very precisely.

   I use Kepler's equation instead so see these:

   • Solving Kepler's equation
   • C++ implementation.

   If you still want to use gravity model then use JPL horizons from NASA. I think they have also source codes in C/C++ there but they are using incompatible reference frame with mine maps so it is unusable for me.

   In general Kepler's equation has bigger error but it is not increasing with time too much. The gravity model is more precise but its error is rising with time and you need to update the astro body data continuously to make it work ...

[edit1] integration precision

your current implementation is like this:

// object variables
double  acc[3],vel[3],pos[3];
// timer iteration
double dt=timer.interval;
for (int i=0;i<3;i++)
 {
 vel[i]+=acc[i]*dt;
 pos[i]+=vel[i]*dt;
 }

The problem is when you are adding very small and very big value then they are shifted to the same exponent before addition which will round off significant data To avoid this just change it to this:

// object variables
double          vel0[3],pos0[3]; // low
double          vel1[3],pos1[3]; // high
double  acc [3],vel [3],pos [3]; // full
// timer iteration
double dt =timer.interval;
double max=10.0; // precision range constant
for (int i=0;i<3;i++)
 {
 vel0[i]+=acc[i]*dt; if (fabs(vel0[i]>=max)) { vel1[i]+=vel0[i]; vel0[i]=0.0; } vel[i]=vel0[i]+vel1[i];
 pos0[i]+=vel[i]*dt; if (fabs(pos0[i]>=max)) { pos1[i]+=pos0[i]; pos0[i]=0.0; } pos[i]=pos0[i]+pos1[i];
 }

Now xxx0 is integrated up to max and the whole thing is added to xxx1

The rounding is still there but it is not cumulative anymore. You have to select max value that the integration itself is safe and also the addition xxx0+xxx1 have to be safe. So if the numbers are too different for one spliting then split twice or more ...

• like: xxx0+=yyy*dt; if (fabs(xxx0>max0))... if (fabs(xxx1>max1))...

[Edit2] Stars

• The Colors of the Stars The best star visualization I ever saw
• Star B-V color index to apparent RGB color all star catalogues uses B-V index
• Using Stelar catalogs also star names cross-reference link is there
• Skybox: combine different star data

[Edit3] Improving Newton D'ALembert integration precision even more

The basic problem with iterative integration is that applaying gravity based acceleration based on current body position will result in bigger orbits because durring integration step dt the position changes a bit which is not accounted in the naive integration. To remedy this take a look at this picture:

Let assume our body is at circular orbit and in the 0 deg position. Instead of using acceleration direction based on current position I used position after 0.5*dt. This augment the acceleration small bit resulting in much much higher precision (correspondence to Kepler orbits). With this tweak I was able to successfully convert from Kepler orbit into Newton D'Alembert for 2 body system. (doing this for n-body is the next step). Of coarse correlating with real data from our solar system is only possible for 2 body system not affected by tidal effects and or moons. To construct own fictional data you can use Kepler circular orbit and contripedal force equalizing gravity:

G = 6.67384e-11;
v = sqrt(G*M/a);                           // orbital speed
T = sqrt((4.0*M_PI*M_PI*a*a*a)/(G*(m+M))); // orbital period

where a is the circular orbit radius m is body mass , M is focal body mass (sun). To maintain precision in acceptable tolerance (for me) the integration step dt should be:

dt = 0.000001*T

So to put new body for testing just put it at:

pos = (a,0,0)
vel = (0,sqrt(G*M/a),0)

While main focal body (Sun) is at:

pos = (0,0,0)
vel = (0,0,0)

This will place your body in circular orbit so you can compare Kepler versus Newton D'Alembert to evaluate precision of your simulation.

这篇关于就尺寸和质量而言，是否可以进行逼真的n体太阳系仿真?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！