星体运动模拟
为了使用C++14和EasyX实现多星体运行系统的图形化轨迹动画，可以按照以下步骤进行：
如果您复制文件的时候提示权限不足，那是因为您没有管理员权限对C盘Program Files文件夹进行写入操作，但是别忘了您能读取，所以请将Dev-cpp整个复制到D盘后进行操作
终版2025.3.9
#include <graphics.h>
#include <bits/stdc++.h>
using namespace std;
// 引力常数，已缩放以适应数值稳定性
double G = 6.674184 * 1e-1; // * 1e10
double step = 200; // 时间步长
struct Object
{
double x = 0, y = 0;
double v_x = 0, v_y = 0;
double a_x = 0, a_y = 0;
double mass = 0;
list<pair<double, double>> trajectory; // 存储轨迹点
bool central = false;
double radius = 0;
};
vector<Object> lists;
// 计算所有天体间的引力加速度
void compute_all_accelerations()
{
// 重置加速度
for (auto &obj : lists)
{
obj.a_x = 0;
obj.a_y = 0;
}
// 双重循环计算每对天体间的相互作用
for (int i = 0; i < lists.size(); ++i)
{
for (int j = i + 1; j < lists.size(); ++j)
{
Object &a = lists[i];
Object &b = lists[j];
double dx = b.x - a.x;
double dy = b.y - a.y;
double r_sq = dx * dx + dy * dy + 1e-11; // 防止除以零
double r = sqrt(r_sq);
double F = G * a.mass * b.mass / r_sq;
double Fx = F * dx / r;
double Fy = F * dy / r;
// 更新加速度（牛顿第三定律）
a.a_x += Fx / a.mass;
a.a_y += Fy / a.mass;
b.a_x -= Fx / b.mass;
b.a_y -= Fy / b.mass;
}
}
for (int i = 0 ; i < lists.size() ; ++ i)
{
if (lists[i].central == true)
{
lists[i].a_x = lists[i].a_y = 0;
}
}
}
// 创建新天体
void create_object(double x, double y, double mass, double vx = 0, double vy = 0, bool central = false) // 位置 * 1e5 质量 * 1e10 速度 * 1e5
{
Object obj;
obj.x = x * 1e5;
obj.y = y * 1e5;
obj.mass = mass;
obj.v_x = vx * 1e5;
obj.v_y = vy * 1e5;
obj.central = central;
obj.radius = min(10, max(3, static_cast<int>(sqrt(obj.mass) * 0.1)));
lists.push_back(obj);
}
void union_object(int x, int y) // x和y是lists里的位置
{
if (lists[x].mass < lists[y].mass)
{
swap(x, y);
}
Object obj;
obj.mass = lists[x].mass + lists[y].mass;
obj.radius = pow(pow(lists[y].radius, 3) + pow(lists[x].radius, 3), 1 / 3.0);
obj.x = lists[x].x;
obj.y = lists[x].y;
obj.central = (lists[x].central or lists[y].central);
if (!obj.central)
{
// 遗传速度 （假设无能量损耗）
// 只能假设大的静止了，不然碰撞完会以相反的方向运动
lists[x].v_x = lists[x].v_y = 0;
obj.v_x = (lists[x].v_x * lists[x].mass + lists[y].v_x + lists[y].mass) / (lists[x].mass + lists[y].mass);
obj.v_y = (lists[x].v_y * lists[x].mass + lists[y].v_y + lists[y].mass) / (lists[x].mass + lists[y].mass);
cerr << lists[y].v_x << ' ' << lists[y].v_y << '\n';
cerr << obj.v_x << ' ' << obj.v_y << '\n';
}
if (x > y)
{
swap(x, y);
}
lists.erase(lists.begin() + y);
lists.erase(lists.begin() + x);
lists.push_back(obj);
}
void Input()
{
cout << "质量将会 * 10^10, 单位m将会 * 10^5\n输入0启用默认演示\n";
int T;
cin >> T;
if (T == 0)
{
create_object(0, 0, 1e10, 0, 0, true);
create_object(- 30, 0, 1, 0, 6e-4, false);
return;
}
cout << "输入格式: x, y, mass, vx, vy, 是否为中心天体(bool)\n";
while (T --)
{
double x, y, mass, vx, vy;
bool central;
cin >> x >> y >> mass >> vx >> vy;
cin >> central;
create_object(x, y, mass, vx, vy, central);
}
}
signed main()
{
initgraph(800, 600, EW_SHOWCONSOLE);
setorigin(400, 300); // 坐标系中心
setaspectratio(1, -1); // Y轴向上为正
const double scale = 1e-5; // 像素缩放比例
const int max_trail = 5000; // 最大轨迹长度
const int clock_pord = 10; // 计算多少次保留一次痕迹（绘制一次屏幕）
// 初始化天体系统
Input();
ExMessage msg;
bool running = true;
long long int clocks = 0;
while (running)
{
// 处理退出事件
if (peekmessage(&msg, EX_KEY) && msg.message == WM_KEYDOWN)
{
if (msg.vkcode == VK_ESCAPE)
{
running = false;
}
}
compute_all_accelerations();
if (clocks % clock_pord == 0)
{
cleardevice();
}
// 更新速度和位置
for (auto &obj : lists)
{
obj.v_x += obj.a_x * step;
obj.v_y += obj.a_y * step;
obj.x += obj.v_x * step;
obj.y += obj.v_y * step;
// cerr << obj.v_x << '\n' << obj.v_y << '\n';
// 更新轨迹
if (clocks % clock_pord == 0)
{
obj.trajectory.emplace_back(obj.x, obj.y);
if (obj.trajectory.size() > max_trail)
{
obj.trajectory.pop_front();
}
}
}
// 判断是否发生碰撞
for (int i = 0 ; i < lists.size() ; i ++)
{
for (int j = i + 1 ; j < lists.size() ; j ++)
{
if (sqrt(pow((lists[i].x - lists[j].x), 2) + pow((lists[i].y - lists[j].y), 2)) < 1e5 * lists[j].radius / 1.05)
{
union_object(i, j);
}
}
}
// 绘制场景
if (clocks % clock_pord == 0)
{
cleardevice();
// 绘制轨迹
for (const auto& obj : lists)
{
COLORREF trail_color = (obj.mass > 2500) ? 0xAAAAAA : 0x6666FF;
setlinecolor(trail_color);
// 使用 polyline 批量绘制轨迹（效率提升）
if (!obj.trajectory.empty())
{
POINT* pts = new POINT[obj.trajectory.size()];
int i = 0;
for (auto& p : obj.trajectory)
{
pts[i].x = static_cast<int>(p.first * scale);
pts[i].y = static_cast<int>(p.second * scale);
i++;
}
setlinecolor(trail_color);
polyline(pts, i);
delete[] pts;
}
}
// 绘制天体
for (const auto &obj : lists)
{
int radius = obj.radius;
COLORREF color = (obj.mass > 1e9) ? RED : CYAN;
setfillcolor(color);
setlinecolor(color);
fillcircle(static_cast<int>(obj.x * scale),
static_cast<int>(obj.y * scale),
radius);
}
FlushBatchDraw();
Sleep(5); // 控制动画速度
}
clocks ++;
}
closegraph();
return 0;
}
/*
-static-libstdc++ -static-libgcc -leasyx -lgdi32 -lole32 -lwinmm
*/
关键说明
-
物理计算优化：
- 使用双重循环计算天体间相互作用力，确保 \(O(\frac{n²}{2})\) 时间复杂度
- 添加软化长度（1e-11）防止除零错误
- 采用牛顿第三定律减少计算量
-
图形渲染：
- 坐标系中心设为窗口中心
- Y轴反转以实现天文学常用坐标系
- 实时轨迹绘制（保留最近500个位置点）
-
参数调校：
G值经过缩放处理以适应数值稳定性step值平衡精度和性能- 示例初始化双星系统参数经过验证
- 剩余的注释里有
Beta （添加帧率控制）:
#include <graphics.h>
#include <bits/stdc++.h>
#include <conio.h>
using namespace std;
// 引力常数，已缩放以适应数值稳定性
double G = 6.674184 * 1e-1; // * 1e10
double step = 200; // 时间步长
struct Object
{
double x = 0, y = 0;
double v_x = 0, v_y = 0;
double a_x = 0, a_y = 0;
double mass = 0;
list<pair<double, double>> trajectory; // 存储轨迹点
bool central = false;
double radius = 0;
};
vector<Object> lists;
// 计算所有天体间的引力加速度
void compute_all_accelerations()
{
// 重置加速度
for (auto &obj : lists)
{
obj.a_x = 0;
obj.a_y = 0;
}
// 双重循环计算每对天体间的相互作用
for (int i = 0; i < lists.size(); ++i)
{
for (int j = i + 1; j < lists.size(); ++j)
{
Object &a = lists[i];
Object &b = lists[j];
double dx = b.x - a.x;
double dy = b.y - a.y;
double r_sq = dx * dx + dy * dy + 1e-11; // 防止除以零
double r = sqrt(r_sq);
double F = G * a.mass * b.mass / r_sq;
double Fx = F * dx / r;
double Fy = F * dy / r;
// 更新加速度（牛顿第三定律）
a.a_x += Fx / a.mass;
a.a_y += Fy / a.mass;
b.a_x -= Fx / b.mass;
b.a_y -= Fy / b.mass;
}
}
for (int i = 0 ; i < lists.size() ; ++ i)
{
if (lists[i].central == true)
{
lists[i].a_x = lists[i].a_y = 0;
}
}
}
// 创建新天体
void create_object(double x, double y, double mass, double vx = 0, double vy = 0, bool central = false) // 位置 * 1e5 质量 * 1e10 速度 * 1e5
{
Object obj;
obj.x = x * 1e5;
obj.y = y * 1e5;
obj.mass = mass;
obj.v_x = vx * 1e5;
obj.v_y = vy * 1e5;
obj.central = central;
obj.radius = min(10, max(3, static_cast<int>(sqrt(obj.mass) * 0.1)));
lists.push_back(obj);
}
void union_object(int x, int y) // x和y是lists里的位置
{
if (lists[x].mass < lists[y].mass)
{
swap(x, y);
}
Object obj;
obj.mass = lists[x].mass + lists[y].mass;
obj.radius = pow(pow(lists[y].radius, 3) + pow(lists[x].radius, 3), 1 / 3.0);
obj.x = lists[x].x;
obj.y = lists[x].y;
obj.central = (lists[x].central or lists[y].central);
if (!obj.central)
{
// 遗传速度 （假设无能量损耗）
// 只能假设大的静止了，不然碰撞完会以相反的方向运动
lists[x].v_x = lists[x].v_y = 0;
obj.v_x = (lists[x].v_x * lists[x].mass + lists[y].v_x + lists[y].mass) / (lists[x].mass + lists[y].mass);
obj.v_y = (lists[x].v_y * lists[x].mass + lists[y].v_y + lists[y].mass) / (lists[x].mass + lists[y].mass);
cerr << lists[y].v_x << ' ' << lists[y].v_y << '\n';
cerr << obj.v_x << ' ' << obj.v_y << '\n';
}
if (x > y)
{
swap(x, y);
}
lists.erase(lists.begin() + y);
lists.erase(lists.begin() + x);
lists.push_back(obj);
}
void Input()
{
cout << "质量将会 * 10^10, 单位m将会 * 10^5\n输入0启用默认演示\n";
int T;
cin >> T;
if (T == 0)
{
create_object(0, 0, 1e10, 0, 0, true);
create_object(- 30, 0, 1, 0, 6e-4, false);
return;
}
cout << "输入格式: x, y, mass, vx, vy, 是否为中心天体(bool)\n";
while (T --)
{
double x, y, mass, vx, vy;
bool central;
cin >> x >> y >> mass >> vx >> vy;
cin >> central;
create_object(x, y, mass, vx, vy, central);
}
}
signed main()
{
timeBeginPeriod(1); // 提高系统定时器时钟分辨率为1毫秒
initgraph(800, 600, EW_SHOWCONSOLE);
setorigin(400, 300); // 坐标系中心
setaspectratio(1, -1); // Y轴向上为正
const double scale = 1e-5; // 像素缩放比例
const int max_trail = 50000; // 最大轨迹长度
const int clock_pord = 10; // 计算多少次保留一次痕迹（绘制一次屏幕）
settextstyle(16, 8, _T("consolas"));
settextcolor(WHITE);
// 初始化天体系统
Input();
ExMessage msg;
bool running = true;
long long int clocks = 0;
int frame_count = 0;
int last_FPS = 0;
// 开始时间，结束时间，频率F
LARGE_INTEGER startCount, endCount, F;
// 获取频率F
QueryPerformanceFrequency(&F);
while (running)
{
// 获取起始计数
QueryPerformanceCounter(&startCount);
// 处理退出事件
if (peekmessage(&msg, EX_KEY) && msg.message == WM_KEYDOWN)
{
if (msg.vkcode == VK_ESCAPE)
{
running = false;
}
}
compute_all_accelerations();
if (clocks % clock_pord == 0)
{
cleardevice();
}
// 更新速度和位置
for (auto &obj : lists)
{
obj.v_x += obj.a_x * step;
obj.v_y += obj.a_y * step;
obj.x += obj.v_x * step;
obj.y += obj.v_y * step;
// cerr << obj.v_x << '\n' << obj.v_y << '\n';
// 更新轨迹
if (clocks % clock_pord == 0)
{
obj.trajectory.emplace_back(obj.x, obj.y);
if (obj.trajectory.size() > max_trail)
{
obj.trajectory.pop_front();
}
}
}
// 判断是否发生碰撞
for (int i = 0 ; i < lists.size() ; i ++)
{
for (int j = i + 1 ; j < lists.size() ; j ++)
{
if (sqrt(pow((lists[i].x - lists[j].x), 2) + pow((lists[i].y - lists[j].y), 2)) < 1e5 * lists[j].radius / 1.05)
{
union_object(i, j);
}
}
}
// 绘制场景
if (clocks % clock_pord == 0)
{
frame_count ++;
cleardevice();
static DWORD last_time = GetTickCount();
if (GetTickCount() - last_time > 1000)
{
// cerr << "FPS: " << frame_count << '\n';
last_FPS = frame_count;
frame_count = 0;
last_time = GetTickCount();
}
// 绘制信息表
// FPS
setaspectratio(1, 1); // Y轴向下为正
outtextxy(- 300, - 200, _T(("FPS:" + to_string(last_FPS)).c_str()));
setaspectratio(1, - 1); // Y轴向上为正
// 绘制轨迹
for (const auto& obj : lists)
{
COLORREF trail_color = (obj.mass > 2500) ? 0xAAAAAA : 0x6666FF;
setlinecolor(trail_color);
// 使用 polyline 批量绘制轨迹（效率提升）
if (!obj.trajectory.empty())
{
POINT* pts = new POINT[obj.trajectory.size()];
int i = 0;
for (auto& p : obj.trajectory)
{
pts[i].x = static_cast<int>(p.first * scale);
pts[i].y = static_cast<int>(p.second * scale);
i++;
}
setlinecolor(trail_color);
polyline(pts, i);
delete[] pts;
}
}
// 绘制天体
for (const auto &obj : lists)
{
int radius = obj.radius;
COLORREF color = (obj.mass > 1e9) ? RED : CYAN;
setfillcolor(color);
setlinecolor(color);
fillcircle(static_cast<int>(obj.x * scale),
static_cast<int>(obj.y * scale),
radius);
}
FlushBatchDraw();
// Sleep(5); // 控制动画速度
// 获取结束计数
QueryPerformanceCounter(&endCount);
// 计算时差
long long elapse = (endCount.QuadPart - startCount.QuadPart) / F.QuadPart * 1000000;
while(elapse < 8000)
{
Sleep(1);
// 重新获取结束时间
QueryPerformanceCounter(&endCount);
// 更新时差
elapse = (endCount.QuadPart - startCount.QuadPart) * 1000000 / F.QuadPart;
}
}
clocks ++;
}
closegraph();
//timeBeginPeriod与timeEndPeriod必须成对存在，并且传一样的值
timeEndPeriod(1);
return 0;
}
/*
-static-libstdc++ -static-libgcc -leasyx -lgdi32 -lole32 -lwinmm
*/
浙公网安备 33010602011771号