用C++、Easy X模拟三体运动教程

御风而行_652

于 2024-02-05 23:57:20 发布

阅读量1.8k

点赞数 45

文章标签： c++ 开发语言

本文链接：https://blog.csdn.net/2301_82057487/article/details/136049744

版权

为了回馈和支持互联网学习氛围，决定将之前写的一些有趣的程序在过年的时候集中开源。

使用C/C++和Easy X窗口绘制插件模拟三维三体运动。完整代码在最后。

想学的同志们可以先来看一看程序效果，看一下是不是自己想要的东西。

三体运动-五分钟纯享

Easy X是C语言窗口绘制免费插件，安装非常简单。

首先，只要有数学上的算法，那么程序编写都是非常简单的事情，所以要按照“计算机数学”的思维将牛顿引力公式模拟出来。

牛顿引力加运动学的公式在数学上的表达式为a=(G*m2)/d^2，a为星球的加速度，m1和m2为两颗星球的质量，d为距离，G为引力常量，到时候用C语言写出来即可。

接下来是程序整体设计框架和总体流程。

一般一个这样的程序比较大，对于没有类似算法经验的开发人来说，开发时间会比较长，所以为了不要自己用自己的程序把自己看晕，尽量大量用自定义的函数来代替主函数里面的细节算法，以下是从“打开窗口”到“依次画出每个天体”的代码，

    initgraph(wx, wy);			//打开窗口
    setorigin(wx / 2, wy / 2);		//重设原点
    setaspectratio(d, -d);			//放缩，方便展示和调参
    while (true)                  //实时操作开始
    {
        for (short i = 0; i < vn; i++)
        {
            for (short i = 0; i < n; i++)//依次计算每个天体的速度，要先把所有v计算完之后才能移动，若计算完v紧接着就移动，就会对其它天体的计算造成误差
            {
            v(i);         //此时，i为主角天体，带入函数后，i等于函数中的j
            }
            for (short i = 0; i < n; i++)         //依次移动每个天体
            {
                move(i);        //此时，i为主角天体，带入函数后，i等于函数中的j
            }
        }
        
        for (short i = 0; i < n; i++)         //依次画出每个天体
        {
            perspective(i);        //此时，i为主角天体，带入函数后，i等于函数中的j
        }

上面的v(i)和move(i)函数都写到另外的地方，不要在主函数里面添乱。

因为是实时计算，所以主要的程序模块都是在无限循环里面的。当然，如果不想无限循环，也可以用Easy X里面与键盘有关的功能去控制，这一点在程序调试的时候有时候会比较有用。

然后就是天体的数学模型构建。一个天体在数学上只需要vx（X轴速度）, vy（Y轴速度）, vz（Z轴速度）, m（质量）, r（半径）, x（X坐标）, y（Y坐标）, z（坐标）这八个参数来描述，使用结构来将这些参数装在一起，

struct Star
{
    double vx, vy, vz, m, r, x, y, z; 
};

并实例化，此处直接用结构形式的数组，

Star star[20]; //和int [1145];是一个道理。改一下那个20就可以模拟1000个星体的运动，如银河系

那为什么天体的数学参数里没有加速度呢，非常简单，速度是累计的，迭代的，或者说v1=v0+a*△t，此时的速度需要用上一时刻，上一次循环或者说上一帧，的速度才能算出来，意思就是说要储存上一时刻的速度v0，但是加速度a是必须根据此时这一帧的参数计算的，或者说上一时刻的加速度跟这一时刻也没有用，就不用专门开个内存空间来存起来。

这里一定要专门说一下，有一个现象在做天体运动的时候会非常常见，那就是两个天体会转着转着就相互靠近，但是无论牛顿的引力还是相对论的模型里面都肯定不该是这样的，这个疑难杂症的问题非常简单，但发现起来非常不容易。

看图一中的“依次计算每个天体的速度”和“依次画出每个天体”是在两个循环里面的，意思就是说，我的程序是先把所有人的vx, vy, vz速度算完，此时还没有移动天体，没有改变x, y, z，然后再用v1=v0+a*△t来进行位移。

只要这么做了之后，就不会相互慢慢靠近。原理说起来比较麻烦，但是不用大物和高数，就是一个小学问题。

图二是每个天体计算速度的时候紧接着就移动该天体，

天体1是正在计算的天体，也叫主角天体，天体1计算好速度之后就马上移动天体1，从上面的图到下面的图，然后再计算天体2的加速度、速度，但是天体2在计算加速度的时候，a=(G*m2)/d^2，d已经变成了d2，但是本来应该是d1的，因为这个天体运动是没有先后顺序的嘛，d1才是这一帧原本的距离，因为天体1提前移动了，所以天体2计算距离的时候就有误差，这是算法的问题，不是计算机小数保留不够多的问题，所以计算机再好，小数点保留再多，也还是会有这个问题的。

图三是改进后的算法。

计算速度的算法与物理学上的计算流程可以不太一样，这个主要就是数学计算流程上的优化了，但是大致比较简单，比较简单就不用流程图了，不必在此过多纠结，但是可以注意一下细节，比如说下面第三行的加速度，不要放到后面任何一个循环里面，因为本来加速度就是一个算完就不用了，但是下一次算还需要这种临时的变量，这一行代码就不必放到循环里面反复创建双精度浮点变量，在主函数中的循环里面调用这个函数的时候也不用每次都初始化加速度，这样优化可以减少计算量，加快程序速度。

void v(short j)                  //j是主角天体

{

    double a;

    for (short i = 0; i < n; i++)        //叠加投影后的加速度，计算瞬时投影速度

    {

        if (i != j)

        {

            a = G * star[i].m / pow((pow((star[i].x - star[j].x), 2.0) + pow((star[i].y - star[j].y), 2.0) + pow((star[i].z - star[j].z), 2.0)), 1.5);     //这个a不是完整的加速度

            if (a > amax)

            {

                a = amax;                                 //限制它的加速度

            }

            star[j].vx = star[j].vx + (star[i].x - star[j].x) * a * t * (1/vn);

            star[j].vy = star[j].vy + (star[i].y - star[j].y) * a * t * (1 / vn);

            star[j].vz = star[j].vz + (star[i].z - star[j].z) * a * t* (1 / vn);

        }

    }

}

看起来比较复杂，当然，多了依托限制天体最大速度的代码，这是我在调参的时候加的，因为我在调参的时候发现天体的参数一旦设计得不好，天体就容易撞在一起，然后根据a=(G*m2)/d^2，d趋近于零，然后加速度和速度暴增，一下就飞出视野了。

然后是位移，即x=x0+v*t，用代码写出来就是下面的这一段，

void move(short j)

{

    star[j].x = star[j].x + star[j].vx * t * (1 / vn) * 0.001;

    star[j].y = star[j].y + star[j].vy * t * (1 / vn) * 0.001;

    star[j].z = star[j].z + star[j].vz * t * (1 / vn) * 0.001;

}

至此，程序的主体大体构建完成了，之后，再调亿下参数就可以运行了，我大概调了六十多组数据，每组数据都包括三个天体的那几个参数，还有其它像最小时间间隔这些杂七杂八的参数，而且更裂开的是，我不知道要调这么多的数据，没这个经验，但是已经把手动初始化输入写好了，所以这六十多组数据全部是用手敲进去的，而且后期调参的时候，不管怎么去调，老是会出现那个相互靠近的现象，我以为这个时间细分得足够，或者把所有都换成双精度浮点就可以了，然后就一直调。

以下是完整的程序。

#include <iostream>
#include <graphics.h>
#include <windows.h>
#include <ctime>
#include <cstdlib>
using namespace std;

struct Star
{
    double vx, vy, vz, m, r, x, y, z;         //v是瞬时速度，位移不能用用x=x+v0*t+0.5*a*t*t计算，因为无法计算a的平均值
};
Star star[20];
short wx,wy;                                                       //窗口大小
short n;                                                         //天体数目
short t;                                                  //模拟的最小时间间隔
float G;                                                  //引力常量
clock_t clockback;                                                  //记录上一次计算的时间
short clean;                                                         //clear用于间隔清屏的时间
short amax;                                                          //天体最大速度
short clear;
short random;
float d;                                                       //放大倍数
short e;
float vn;                                                    //精度，将时间间隔t分为vn份,如果用short数据类型，（1/vn）计算结果就是0
void v(short j)                                                   //j是主角天体
{
    double a;
    for (short i = 0; i < n; i++)                                         //叠加投影后的加速度，计算瞬时投影速度
    {
        if (i != j)
        {
            a = G * star[i].m / pow((pow((star[i].x - star[j].x), 2.0) + pow((star[i].y - star[j].y), 2.0) + pow((star[i].z - star[j].z), 2.0)), 1.5);//这个a不是完整的加速度
            /*if (a > amax)
            {
                a = amax;
                star[j].vx = star[j].vx + (star[i].x - star[j].x) * a * t;
                star[j].vy = star[j].vy + (star[i].y - star[j].y) * a * t;
                star[j].vz = star[j].vz + (star[i].z - star[j].z) * a * t;
            }
            else 
            {
                star[j].vx = star[j].vx + (star[i].x - star[j].x) * a * t;
                star[j].vy = star[j].vy + (star[i].y - star[j].y) * a * t;
                star[j].vz = star[j].vz + (star[i].z - star[j].z) * a * t;          //速度计算自带正负
            }
            //star[j].vx = star[j].vx + (star[i].x - star[j].x) * a * t;
            //star[j].vy = star[j].vy + (star[i].y - star[j].y) * a * t;
            //star[j].vz = star[j].vz + (star[i].z - star[j].z) * a * t;          //速度计算自带正负*/
            if (a > amax)
            {
                a = amax;                                 //限制它的加速度
            }
            star[j].vx = star[j].vx + (star[i].x - star[j].x) * a * t * (1/vn);
            star[j].vy = star[j].vy + (star[i].y - star[j].y) * a * t * (1 / vn);
            star[j].vz = star[j].vz + (star[i].z - star[j].z) * a * t* (1 / vn);
        }

    }
}
void move(short j)
{
    star[j].x = star[j].x + star[j].vx * t * (1 / vn) * 0.001;
    star[j].y = star[j].y + star[j].vy * t * (1 / vn) * 0.001;
    star[j].z = star[j].z + star[j].vz * t * (1 / vn) * 0.001;
}
void perspective(short i)                                        //透视
{
    solidcircle(star[i].x, star[i].y, star[i].r);
    line(star[i].x/2, star[i].y/2,0,0);
}
void start()                                                   //初始化输入
{
    cout << "输入窗口的横向长：";                                     //初始化开始
    cin >> wx;
    cout << "输入窗口的纵向宽：";
    cin >> wy;
    cout << "输入天体数目：";
    cin >> n;
    cout << "输入模拟时间间隔（毫秒）：";
    cin >> t;
    cout << "输入天体运行的最大加速度：";
    cin >> amax;
    cout << "输入清屏参数：";
    cin >> clear;
    cout << "输入放缩参数：";
    cin >> d;
    cout << "输入引力常量：";
    cin >> G;
    for (short i = 0; i < n; i++)
    {
    cout << "输入天体" << i + 1 << "的x坐标：";
    cin >> star[i].x;
        cout << "输入天体" << i + 1 << "的y坐标：";
        cin >> star[i].y;
        cout << "输入天体" << i + 1 << "的z坐标：";
        cin >> star[i].z;
        cout << "输入天体" << i + 1 << "沿X轴方向的速度：";
        cin >> star[i].vx;
        cout << "输入天体" << i + 1 << "沿Y轴方向的速度：";
        cin >> star[i].vy;
        cout << "输入天体" << i + 1 << "沿Z轴方向的速度：";
        cin >> star[i].vz;
        cout << "输入天体" << i + 1 << "的质量：";
        cin >> star[i].m;
        cout << "输入天体" << i + 1 << "的半径：";
        cin >> star[i].r;
    }                                                     //初始化结束
}
void randomstart()
{
    cout << "输入天体数目：";
    cin >> n;
    amax = 1000;
    srand((unsigned)time(NULL));
    star[n - 1].m = rand()/100;
    star[n].r = star[n].m / 10;
    for (short i = 0; i < n-1; i++)
    {
        /*srand((unsigned)time(NULL) + i);
        star[i].x = rand()/10000000;
        srand((unsigned)time(NULL)+i*2);
        star[i].y = rand() / 10000000;
        srand((unsigned)time(NULL)+i*3);
        star[i].z = rand() / 10000000;
        srand((unsigned)time(NULL)+i*4);
        star[i].vx = rand() / 10000000;
        srand((unsigned)time(NULL)+i*5);
        star[i].vy = rand() / 10000000;
        srand((unsigned)time(NULL)+i*6);
        star[i].vz = rand() / 10000000;
        srand((unsigned)time(NULL)+i*7);
        star[i].m = rand() / 10000000;
        star[i].r = star[i].m/10;*/


        star[i].x = rand()/100;
        star[i].y = rand() / 100;
        star[i].z = rand() / 100;
        star[i].vx = rand() / 10;
        star[i].vy = rand() / 10;
        star[i].vz = rand() / 10;
        star[i].m = rand() / 100;
        star[i].r = star[i].m / 10;
        star[n-1].vx = -star[i].vx* star[i].m/ star[n - 1].m;
        star[n - 1].vy = -star[i].vy * star[i].m/ star[n - 1].m;
        star[n - 1].vz = -star[i].vz * star[i].m/ star[n - 1].m;
        //cout << star[n-1].vx;
        //cin >> n;
    }
}
int main()                      //主函数开始
{
    wx = 1000;
    wy = 600;
    n = 4;
    star[0].x = -100;
    star[0].y = 0;
    star[0].z = 0;
    star[0].vx = 1000;
    star[0].vy = -500;
    star[0].vz = 200;
    star[0].m = 100;
    star[0].r = 10;

    star[1].x = 400;
    star[1].y = -80;
    star[1].z = 50;
    star[1].vx = -2000;
    star[1].vy = -500;
    star[1].vz = -150;
    star[1].m = 100;
    star[1].r = 10;

    star[2].x = -120;
    star[2].y = 25;
    star[2].z = 10;
    star[2].vx = 1000;
    star[2].vy = 1000;
    star[2].vz = -50;
    star[2].m = 100;
    star[2].r = 10;

    star[3].x = 20;
    star[3].y = -150;
    star[3].z = 100;
    star[3].vx = 0;
    star[3].vy = 0;
    star[3].vz = 10;
    star[3].m = 0;
    star[3].r = 5;
    t = 1;
    amax = 20;
    G = 20000;
    clear = 3;
    d = 0.5;
    e = 1;
    vn = 10; cout << 1 / vn;
    cout << "是否使用默认数据（默认则输入1，随机则输入2，手动输入则输入3）：";
    cin >> random;
    if (random == 2)
    {
        randomstart();
    }
    else if(random == 3)
    {
        start();
    }
    
    initgraph(wx, wy);
    setorigin(wx / 2, wy / 2);
    setaspectratio(d, -d);
    while (true)                //实时操作开始
    {
        for (short i = 0; i < vn; i++)
        {
            for (short i = 0; i < n; i++)       //依次计算每个天体的速度，要先把所有v计算完之后才能移动，若计算完v紧接着就移动，就会对其它天体的计算造成误差
            {
            v(i);                                 //此时，i为主角天体，带入函数后，i等于函数中的j
            }
            for (short i = 0; i < n; i++)                        //依次移动每个天体
            {
                move(i);                                 //此时，i为主角天体，带入函数后，i等于函数中的j
            }
        }
        
        for (short i = 0; i < n; i++)                        //依次画出每个天体
        {
            perspective(i);                                 //此时，i为主角天体，带入函数后，i等于函数中的j
        }
        if (random == 1)
        {
        if (star[n - 1].x < wx)
        {
            if (star[n - 1].x > -wx)
            {
                if (star[n - 1].y < wy)
                {
                    if (star[n - 1].y > -wy)
                    {
                        e=1;                                      //要在画完和清屏之间暂停
                    }
                }
            }
        }
        if (e != 1)
        {
            star[n - 1].x = -star[n - 1].x*0.9;
            star[n - 1].y = -star[n - 1].y*0.9;
            star[n - 1].z = -star[n - 1].x*0.9;
        }
        e = 0;
        }
        
        /*int f;
        if (star[3].x < 1000)
        {
            f = f + 1;
        }
        if (star[3].x >- 1000)
        {
            f = f + 2;
        }
        
        
        if (f = 3)
        {
            Sleep(t);
        }
        if (f = 2)
        {
            Sleep(t+1);
        }
        if (f = 1)
        {
            Sleep(t + 1);
        }*/
        
        Sleep(t);                                      //要在画完和清屏之间暂停
        if (clean== clear)                                             //用于减少清屏的频率
        {
            cleardevice();
            clean = 0;
        }
        else
        {
            clean = clean++;
        }
    }
    return 0;
}

以上是全部的代码，包括烂尾的透视函数和调试的废代码（doge）。

以下是牢骚。

这个是我写过的第三个C++程序，当时是高三寒假，从一开始的设想美好，装逼必备，到之后的灰头土脸，到中间的寻找bug，到最后和同学说起这个的时候发现那个相互靠近的问题的原理，到最后做好发B站视频，还有程序优化，中午拿到教室里兄弟们面前大装特装，一路走下来，感觉真的以后如果能只干体力活就能挣钱，那就太好了。

C++因为是自学的，所以不是那么熟练，但是有之前加密器和熵增模拟的程序的铺垫，还是能有些管理上和程序上的经验，比如说一开始是要做三维的，但是真的难度较大而意义不明，做了一晚上还没摸到门道，因为这意味着要手搓一个渲染器诶，虽然其实是可以的，但是毕竟是高三，赶快写完程序比较好，当时是三体动画化带起了一些热潮，高二又在B站大佬的刺激下自学的C++，这不得自己写个程序啊，所以在强烈的好奇心和装逼欲的催促之下，就在寒假开始慢慢长途，而且当时是有我历史上最大的一段学业压力的，很难想象我当时是如何兼顾的，寒假在老家一边调参，一边修改程序，一边写和找bug，最后做了个不三不四的东西就放了，之后高三最后一个学期，和我的好同学，Three Solid Cubes,一个计算机和数学大佬，在风萧之时凭栏的时候聊到这里，突然就悟了，然后就在高三把程序拆开来改，终于没有那个相互靠近的问题了，但是我真的觉得太得不偿失了，一点都不开心，只有偶尔在跟别人聊起这段经历的时候，在我侃侃而谈的时候会有一点骄傲和无奈。