这依然是与《三体》有关的一篇文章。空间中三个星体在万有引力作用下的运动被称之为三体问题，参见我的上一篇文章：三体运动的程序模拟。而这一节，对三体问题进行了扩展，实现了空间中N个星体在万有引力下的运动模拟。
程序中使用了两个物理定律：
（1）万有引力定律
这是牛顿发现的:任意两个质点有通过连心线方向上的力相互吸引。该引力的大小与它们的质量乘积成正比，与它们距离的平方成反比，与两物体的化学本质或物理状态以及中介物质无关。
以数学表示为：
其中：
: 两个物体之间的引力
: 万有引力常数
: 物体1的质量
: 物体2的质量
: 两个物体之间的距离
: 两个物体之间的引力
: 万有引力常数
: 物体1的质量
: 物体2的质量
: 两个物体之间的距离依照国际单位制，F的单位为牛顿(N)，m1和m2的单位为千克(kg)，r 的单位为米(m)，常数G近似地等于6.67 × 10−11 N m2 kg−2（牛顿米的平方每千克的平方）。当然程序中G不可能设置为这么小的一个数,用户可以自己设置万有引力常数,以实现合理的运动.
(2)势能与动能守恒定律
-
引力位能或引力势能是指物体因为大质量物体的万有引力而具有的位能，其大小与其到大质量的距离有关。
其中G为万有引力常数，M、m分别为两物体质量，r为两者距离。
动能，简单的说就是指物体因运动而具有的能量。数值上等于（1/2）Mv^2. 动能是能量的一种，它的国际单位制下单位是焦耳(J)，简称焦。
势能与动能守恒即在整个运动过程中,其总能量是不变的,即N体的动能与势能之和为固定值.
本来我只打算使用万有引力定律,没有考虑势能动能守恒定律.但程序运行后发现,由于没有使用微积分,即使采样间隔时间再小,误差也很大.这让我想起当年学大学物理的时候,大一上学期开高数,然后下学期开了物理.现在对大学物理讲的东西,只记得是用高数中的微积分解物理问题.可那时,我高数学得就是个渣,所以物理也没办法学好.当然就算物理学好了,现在也早忘光了.没有用微积分,这个程序中算法误差不小,模拟一下吧,不要太计较.
核心代码也发一下,希望有兴趣的人可以对其做出优化:
1 /****************************************************************
2
3 File name : NBodyKernel.h
4 Author : 叶峰
5 Version : 1.0a
6 Create Date : 2014/09/19
7 Description :
8
9 *****************************************************************/
10
11 // --------------------------------------------------------------------------------------
12
13 #ifndef _NBodyKernel_H_
14 #define _NBodyKernel_H_
15
16 // INCLUDES -----------------------------------------------------------------------------
17
18 #include "..\..\..\Why\YCommon_h\WhyMath.h"
19 #include "..\..\..\Why\YInterface\Others\YdD3D9Math.h"
20
21 // --------------------------------------------------------------------------------------
22
23 #define YD_ORBIT_VERTICES_COUNT 1024
24
25 // --------------------------------------------------------------------------------------
26
27 class Star
28 {
29 public:
30 Vector3 vPosition; // 位置
31 Vector3 vVelocity; // 速度
32 Vector3 vForce; // 受力
33 Yreal fWeight; // 质量
34 Yreal fRadiusScale; // 缩放
35 Yreal fKineticEnergy; // 动能
36 Ycolor color; // 颜色
37 Ybool bOrbitVisible; // 轨迹可见
38 Yuint id;
39 Vector3 verticesOrbit[YD_ORBIT_VERTICES_COUNT];
40 Star* pNext;
41
42 Star()
43 {
44 fWeight = 0.0f;
45 fRadiusScale = 1.0f;
46 fKineticEnergy = 0.0f;
47 color = 0xffffffff;
48 bOrbitVisible = YD_FALSE;
49 id = -1;
50 memset(verticesOrbit, 0, sizeof(verticesOrbit));
51 pNext = NULL;
52 }
53
54 // 根据速度计算动能
55 void CalculateKineticEnergyByVelocity()
56 {
57 Yreal sqv = D3DXVec3LengthSq(&vVelocity);
58 fKineticEnergy = 0.5f*fWeight*sqv;
59 }
60
61 // 根据动能计算速度
62 void CalculateVelocityByKineticEnergy()
63 {
64 float v = sqrt(2*fKineticEnergy/fWeight);
65 float w = D3DXVec3Length(&vVelocity);
66 vVelocity *= v/w;
67 };
68 };
69
70 // --------------------------------------------------------------------------------------
71
72 class CNBodyKernel
73 {
74 public:
75 CNBodyKernel();
76
77 ~CNBodyKernel();
78
79 // 更新N体
80 void UpdateNBody(Yreal deltaTime);
81
82 // 清空N体
83 void Clear();
84
85 // 设置引力系数
86 void SetGravitationCoefficient(Yreal c);
87
88 // 获得引力系数
89 Yreal GetGravitationCoefficient() const
90 {
91 return m_fGravitationCoefficient;
92 }
93
94 // 返回星体数目
95 Yuint GetStarsCount() const
96 {
97 return m_starsCount;
98 }
99
100 // 返回星体链表
101 const Star* GetStarsList() const
102 {
103 return m_starsList;
104 }
105
106 // 返回星体对象
107 const Star* GetStar(Yuint index) const;
108
109 // 移除星体
110 bool RemoveStar(Yuint index);
111
112 // 添加星体
113 void AddStar(
114 const Vector3& vPosition,
115 const Vector3& vVelocity,
116 const Yreal fWeight,
117 Ycolor color);
118
119
120 // 设置星体的重量
121 void SetStarWeight(Yuint index, Yreal weight);
122
123 // 设置星体的轨迹是否可见
124 void SetStarOrbitVisible(Yuint index, bool visible);
125
126 // 设置星体的颜色
127 void SetStarColor(Yuint index, Ycolor color);
128
129 Yuint GetCurrentOrbit() const
130 {
131 return m_currentOrbit;
132 }
133
134 const Vector3* GetStarPosition(Yuint index) const;
135
136 Yreal GetStarWeight(Yuint index) const;
137
138 bool GetStarOrbitVisible(Yuint index) const;
139
140 Ycolor GetStarColor(Yuint index) const;
141
142 Yuint GetStarID(Yuint index) const;
143
144 void GetSpaceCenter(Vector3& vCenter) const;
145
146 void GetWeightCenter(Vector3& vCenter) const;
147
148 private:
149 // 更新星体的位置,速度,动能
150 void UpdateStar(Star& star, Yreal t);
151
152 // 根据动能重新计算速度
153 void RecalculateStarVelocity(Star& star);
154
155 // 计算每个星体的当前受力
156 void CalculateStarsForce();
157
158 // 计算总势能
159 void CalculateAmountOfPotentialEnergy();
160
161 // 计算总动能
162 void CalculateAmountOfKineticEnergy();
163
164 // 计算总能量
165 void CalculateAmountEnergy()
166 {
167 CalculateAmountOfPotentialEnergy();
168 CalculateAmountOfKineticEnergy();
169 m_fAmountEnergy = m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
170 }
171
172 private:
173 Yreal m_fGravitationCoefficient; // 引力系数
174 Yreal m_fAmountOfPotentialEnergy; // 当前总势能
175 Yreal m_fAmountOfKineticEnergy; // 当前总动能
176 Yreal m_fAmountEnergy; // 当前能量m_fAmountOfPotentialEnergy + m_fAmountOfKineticEnergy;
177
178 Star* m_starsList; // N体链表
179 Yuint m_starsCount; // N体数目
180
181 Yuint m_currentOrbit;
182 };
183
184 // --------------------------------------------------------------------------------------
185
186 #endif
1 /****************************************************************
2
3 File name : NBodyKernel.cpp
4 Author : 叶峰
5 Version : 1.0a
6 Create Date : 2014/09/19
7 Description :
8
9 *****************************************************************/
10
11 // --------------------------------------------------------------------------------------
12
13 #include "..\..\..\Why\YCommon_h\WhyMemoryDefines.h"
14 #include "NBodyKernel.h"
15 #include <assert.h>
16
17 // --------------------------------------------------------------------------------------
18
19 #define YD_ESP 0.00001f
20
21 // --------------------------------------------------------------------------------------
22
23 CNBodyKernel::CNBodyKernel()
24 {
25 m_fGravitationCoefficient = 100.0f;
26 m_fAmountOfPotentialEnergy = 0.0f;
27 m_fAmountOfKineticEnergy = 0.0f;
28 m_fAmountEnergy = 0.0f;
29
30 m_starsList = NULL;
31 m_starsCount = 0;
32
33 m_currentOrbit = 0;
34 }
35
36 CNBodyKernel::~CNBodyKernel()
37 {
38 Clear();
39 }
40
41 // 清空N体
42 void CNBodyKernel::Clear()
43 {
44 Star* pDel = m_starsList;
45 Star* pNext;
46 while (pDel)
47 {
48 pNext = pDel->pNext;
49 YD_DELETE(pDel);
50 pDel = pNext;
51 }
52
53 m_starsList = NULL;
54 m_starsCount = 0;
55 }
56
57 // 获得引力系数
58 void CNBodyKernel::SetGravitationCoefficient(Yreal c)
59 {
60 m_fGravitationCoefficient = c;
61
62 // 计算总能量
63 CalculateAmountEnergy();
64 }
65
66 // 计算每个星体的当前受力
67 void CNBodyKernel::CalculateStarsForce()
68 {
69 // 清空所有引力
70 Star* pStar = m_starsList;
71 while (pStar)
72 {
73 pStar->vForce = Vector3(0.0f, 0.0f, 0.0f);
74 pStar = pStar->pNext;
75 }
76
77 // 计算引力
78 Star* pCurrent = m_starsList;
79 Star* pNext;
80 while (pCurrent)
81 {
82 pNext = pCurrent->pNext;
83 while (pNext)
84 {
85 Vector3 vAB = pNext->vPosition - pCurrent->vPosition;
86 float disSqAB = D3DXVec3LengthSq(&vAB);
87 if (disSqAB < YD_ESP)
88 {
89 disSqAB = YD_ESP;
90 }
91 float disAB = sqrt(disSqAB);
92 Vector3 vForceDir = vAB/disAB;
93 Yreal fForce = m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disSqAB;
94 Vector3 vForce = vForceDir*fForce;
95
96 pCurrent->vForce += vForce;
97 pNext->vForce -= vForce;
98
99 pNext = pNext->pNext;
100 }
101 pCurrent = pCurrent->pNext;
102 }
103 }
104
105 void CNBodyKernel::UpdateStar(Star& star, Yreal t)
106 {
107 // 加速度
108 const Vector3 acc = star.vForce/star.fWeight;
109
110 star.vPosition.x = star.vPosition.x + star.vVelocity.x*t + 0.5f*acc.x*t*t;
111 star.vPosition.y = star.vPosition.y + star.vVelocity.y*t + 0.5f*acc.y*t*t;
112 star.vPosition.z = star.vPosition.z + star.vVelocity.z*t + 0.5f*acc.z*t*t;
113
114 star.vVelocity.x += acc.x*t;
115 star.vVelocity.y += acc.y*t;
116 star.vVelocity.z += acc.z*t;
117
118 // 重新计算动能
119 star.CalculateKineticEnergyByVelocity();
120 }
121
122 void CNBodyKernel::UpdateNBody(Yreal deltaTime)
123 {
124 // 计算每个星体的当前受力
125 CalculateStarsForce();
126
127 // 更新星体的位置与速度
128 Star* pStar = m_starsList;
129 while (pStar)
130 {
131 UpdateStar(*pStar, deltaTime);
132 pStar = pStar->pNext;
133 }
134
135 // 能量守恒
136 CalculateAmountOfKineticEnergy();
137 CalculateAmountOfPotentialEnergy();
138
139 // 理论上的动能
140 Yreal fKineticEnergy = m_fAmountEnergy - m_fAmountOfPotentialEnergy;
141 if (fKineticEnergy < 0.0f)
142 {
143 fKineticEnergy = 0.0f;
144 }
145
146 Yreal r = fKineticEnergy/m_fAmountOfKineticEnergy;
147 pStar = m_starsList;
148 while (pStar)
149 {
150 if (r > YD_ESP)
151 {
152 pStar->fKineticEnergy *= r;
153 }
154 else
155 {
156 pStar->fKineticEnergy = 1.0f;
157 }
158 // 根据动能计算速度
159 pStar->CalculateVelocityByKineticEnergy();
160 pStar = pStar->pNext;
161 }
162
163 // 更新轨迹点
164 pStar = m_starsList;
165 while (pStar)
166 {
167 pStar->verticesOrbit[m_currentOrbit] = pStar->vPosition;
168 pStar = pStar->pNext;
169 }
170 m_currentOrbit++;
171 if (m_currentOrbit >= YD_ORBIT_VERTICES_COUNT)
172 {
173 m_currentOrbit = 0;
174 }
175 }
176
177 // 返回星体对象
178 const Star* CNBodyKernel::GetStar(Yuint index) const
179 {
180 if (index >= m_starsCount)
181 {
182 return NULL;
183 }
184
185 const Star* pStar = m_starsList;
186 while (index && pStar)
187 {
188 index--;
189 pStar = pStar->pNext;
190 }
191
192 return pStar;
193 }
194
195 // 移除星体
196 bool CNBodyKernel::RemoveStar(Yuint index)
197 {
198 if (index >= m_starsCount)
199 {
200 return false;
201 }
202
203 if (index == 0)
204 {
205 Star* pStar = m_starsList->pNext;
206 YD_DELETE(m_starsList);
207 m_starsList = pStar;
208 }
209 else
210 {
211 Star* pStar = m_starsList;
212 Yuint t = index - 1;
213 while (t && pStar)
214 {
215 t--;
216 pStar = pStar->pNext;
217 }
218 assert(pStar);
219
220 Star* pDel = pStar->pNext;
221 pStar->pNext = pDel->pNext;
222 YD_DELETE(pDel);
223 }
224
225 m_starsCount--;
226
227 Yuint id = 0;
228 Star* pStar = m_starsList;
229 while (pStar)
230 {
231 pStar->id = id;
232 id++;
233 pStar = pStar->pNext;
234 }
235
236 // 重算能量
237 CalculateAmountEnergy();
238
239 return true;
240 }
241
242 // 添加星体
243 void CNBodyKernel::AddStar(
244 const Vector3& vPosition,
245 const Vector3& vVelocity,
246 const Yreal fWeight,
247 Ycolor color)
248 {
249 Star* pNewStar = YD_NEW(Star);
250 pNewStar->vPosition = vPosition;
251 pNewStar->vVelocity = vVelocity;
252 pNewStar->fWeight = fWeight;
253 pNewStar->fRadiusScale = yf_pow(fWeight, 1.0f/3.0f);
254 pNewStar->CalculateKineticEnergyByVelocity();
255 pNewStar->color = color;
256 pNewStar->id = m_starsCount;
257 for (Yuint i = 0; i < YD_ORBIT_VERTICES_COUNT; i++)
258 {
259 pNewStar->verticesOrbit[i] = vPosition;
260 }
261
262 if (!m_starsList)
263 {
264 m_starsList = pNewStar;
265 }
266 else
267 {
268 Star* pStar = m_starsList;
269 while (pStar->pNext)
270 {
271 pStar = pStar->pNext;
272 }
273 pStar->pNext = pNewStar;
274 }
275
276 m_starsCount++;
277
278 // 重算能量
279 CalculateAmountEnergy();
280 }
281
282 // 设置星体的重量
283 void CNBodyKernel::SetStarWeight(Yuint index, Yreal weight)
284 {
285 if (index >= m_starsCount)
286 {
287 return;
288 }
289
290 Star* pStar = m_starsList;
291 while (index && pStar)
292 {
293 index--;
294 pStar = pStar->pNext;
295 }
296
297 pStar->fWeight = weight;
298 pStar->fRadiusScale = yf_pow(weight, 1.0f/3.0f);
299 pStar->CalculateKineticEnergyByVelocity();
300
301 // 重算能量
302 CalculateAmountEnergy();
303 }
304
305 // 设置星体的轨迹是否可见
306 void CNBodyKernel::SetStarOrbitVisible(Yuint index, bool visible)
307 {
308 if (index >= m_starsCount)
309 {
310 return;
311 }
312
313 Star* pStar = m_starsList;
314 while (index && pStar)
315 {
316 index--;
317 pStar = pStar->pNext;
318 }
319
320 pStar->bOrbitVisible = visible;
321 }
322
323 // 设置星体的颜色
324 void CNBodyKernel::SetStarColor(Yuint index, Ycolor color)
325 {
326 if (index >= m_starsCount)
327 {
328 return;
329 }
330
331 Star* pStar = m_starsList;
332 while (index && pStar)
333 {
334 index--;
335 pStar = pStar->pNext;
336 }
337
338 pStar->color = color;
339 }
340
341 // 计算总动能
342 void CNBodyKernel::CalculateAmountOfKineticEnergy()
343 {
344 // 动能
345 m_fAmountOfKineticEnergy = 0.0f;
346 Star* pStar = m_starsList;
347 while (pStar)
348 {
349 //pStar->CalculateKineticEnergyByVelocity();
350 m_fAmountOfKineticEnergy += pStar->fKineticEnergy;
351 pStar = pStar->pNext;
352 }
353 }
354
355 // 计算总势能
356 void CNBodyKernel::CalculateAmountOfPotentialEnergy()
357 {
358 m_fAmountOfPotentialEnergy = 0.0f;
359
360 Star* pCurrent = m_starsList;
361 Star* pNext;
362 while (pCurrent)
363 {
364 pNext = pCurrent->pNext;
365 while (pNext)
366 {
367 Vector3 vAB = pCurrent->vPosition - pNext->vPosition;
368 float disAB = D3DXVec3Length(&vAB);
369 float fPotentialEnergyAB = -m_fGravitationCoefficient*pCurrent->fWeight*pNext->fWeight/disAB;
370 m_fAmountOfPotentialEnergy += fPotentialEnergyAB;
371 pNext = pNext->pNext;
372 }
373 pCurrent = pCurrent->pNext;
374 }
375 }
376
377 const Vector3* CNBodyKernel::GetStarPosition(Yuint index) const
378 {
379 if (index >= m_starsCount)
380 {
381 return 0;
382 }
383
384 const Star* pStar = m_starsList;
385 while (index && pStar)
386 {
387 index--;
388 pStar = pStar->pNext;
389 }
390
391 assert(pStar);
392 return &pStar->vPosition;
393 }
394
395 Yreal CNBodyKernel::GetStarWeight(Yuint index) const
396 {
397 if (index >= m_starsCount)
398 {
399 return 0.0f;
400 }
401
402 const Star* pStar = m_starsList;
403 while (index && pStar)
404 {
405 index--;
406 pStar = pStar->pNext;
407 }
408
409 assert(pStar);
410 return pStar->fWeight;
411 }
412
413 bool CNBodyKernel::GetStarOrbitVisible(Yuint index) const
414 {
415 if (index >= m_starsCount)
416 {
417 return false;
418 }
419
420 const Star* pStar = m_starsList;
421 while (index && pStar)
422 {
423 index--;
424 pStar = pStar->pNext;
425 }
426
427 assert(pStar);
428 return pStar->bOrbitVisible != YD_FALSE;
429 }
430
431 Ycolor CNBodyKernel::GetStarColor(Yuint index) const
432 {
433 if (index >= m_starsCount)
434 {
435 return 0;
436 }
437
438 const Star* pStar = m_starsList;
439 while (index && pStar)
440 {
441 index--;
442 pStar = pStar->pNext;
443 }
444
445 assert(pStar);
446 return pStar->color;
447 }
448
449 Yuint CNBodyKernel::GetStarID(Yuint index) const
450 {
451 if (index >= m_starsCount)
452 {
453 return 0;
454 }
455
456 const Star* pStar = m_starsList;
457 while (index && pStar)
458 {
459 index--;
460 pStar = pStar->pNext;
461 }
462
463 assert(pStar);
464 return pStar->id;
465 }
466
467 void CNBodyKernel::GetSpaceCenter(Vector3& vCenter) const
468 {
469 vCenter = Vector3(0.0f, 0.0f, 0.0f);
470 if (!m_starsCount)
471 {
472 return;
473 }
474
475 Star* pStar = m_starsList;
476 while (pStar)
477 {
478 vCenter += pStar->vPosition;
479 pStar = pStar->pNext;
480 }
481
482 vCenter /= (Yreal)m_starsCount;
483 }
484
485 void CNBodyKernel::GetWeightCenter(Vector3& vCenter) const
486 {
487 vCenter = Vector3(0.0f, 0.0f, 0.0f);
488 if (!m_starsCount)
489 {
490 return;
491 }
492
493 Yreal weights = 0.0f;
494 Star* pStar = m_starsList;
495 while (pStar)
496 {
497 vCenter += pStar->vPosition*pStar->fWeight;
498 weights += pStar->fWeight;
499 pStar = pStar->pNext;
500 }
501
502 vCenter /= weights;
503 }
程序启动后,会出现4个球体在运动。不知道大家有没有注意到：Win7的开机动画就是四个球体的运动，所以我这程序的默认也是生成四个球体。
通过UI控件，可以添加删除星体。没有上限，当然星体越多程序越慢。
可以显示每一个星体的运动轨迹
用户可以实时修改任意一个星体的质量
鼠标右键用于控制视角
键盘U用于开关UI用户界面.
通过UI用户界面可以设置三个球体的质量,设置万有引力系数,设置天体运行速度,设置球体的显示大小.
键盘0~9用于开关第N个球体运动轨迹的显示
键盘G,用于开关三维网格的显示
键盘C,用于开关坐标轴的显示
键盘P,用于暂停
键盘R,用于重置.
F11 全屏切换
软件下载地址:http://files.cnblogs.com/WhyEngine/NBody.zip
时间： 2024-10-09 17:06:00