191 lines
6.1 KiB
C++
191 lines
6.1 KiB
C++
#include <cstring>
|
||
#include <iostream>
|
||
#include <cmath>
|
||
#include <vector>
|
||
#include <cstdlib>
|
||
#include <iomanip>
|
||
|
||
using namespace std;
|
||
|
||
// 常量定义
|
||
const double G = 6.67430e-11; // 引力常数 (m^3 kg^-1 s^-2)
|
||
const double DT = 0.01; // 时间步长
|
||
const int TMAX = 100; // 总时间步数
|
||
|
||
// 三维向量结构体
|
||
struct Vec3 {
|
||
double x, y, z;
|
||
|
||
Vec3() : x(0), y(0), z(0) {}
|
||
Vec3(double x, double y, double z) : x(x), y(y), z(z) {}
|
||
|
||
Vec3 operator+(const Vec3& other) const {
|
||
return Vec3(x + other.x, y + other.y, z + other.z);
|
||
}
|
||
|
||
Vec3 operator-(const Vec3& other) const {
|
||
return Vec3(x - other.x, y - other.y, z - other.z);
|
||
}
|
||
|
||
Vec3 operator*(double scalar) const {
|
||
return Vec3(x * scalar, y * scalar, z * scalar);
|
||
}
|
||
|
||
Vec3 operator/(double scalar) const {
|
||
return Vec3(x / scalar, y / scalar, z / scalar);
|
||
}
|
||
|
||
double magnitude() const {
|
||
return sqrt(x*x + y*y + z*z);
|
||
}
|
||
};
|
||
|
||
// 天体结构体
|
||
struct Body {
|
||
double mass; // 质量
|
||
Vec3 position; // 位置
|
||
Vec3 velocity; // 速度
|
||
Vec3 force; // 受力
|
||
};
|
||
|
||
// 计算第i个物体所受的引力
|
||
Vec3 compute_force(int i, const vector<Body>& bodies) {
|
||
Vec3 total_force(0, 0, 0);
|
||
|
||
for (size_t j = 0; j < bodies.size(); j++) {
|
||
if (i == j) continue; // 跳过自己
|
||
|
||
// 计算从物体i指向物体j的向量
|
||
Vec3 r_vec = bodies[j].position - bodies[i].position;
|
||
double distance = r_vec.magnitude();
|
||
|
||
// 避免除以零（物体重合的情况）
|
||
if (distance < 1e-10) continue;
|
||
|
||
// 计算引力大小: F = G * m_i * m_j / r^2
|
||
double force_magnitude = G * bodies[i].mass * bodies[j].mass / (distance * distance);
|
||
|
||
// 计算力的方向（单位向量）
|
||
Vec3 force_direction = r_vec / distance;
|
||
|
||
// 累加力（考虑方向）
|
||
total_force = total_force + force_direction * force_magnitude;
|
||
}
|
||
|
||
return total_force;
|
||
}
|
||
|
||
int main(int argc, char** argv) {
|
||
// 可以从命令行参数获取天体数量
|
||
int n = 4; // 默认4个天体
|
||
bool verbose = false;
|
||
if (argc > 1) {
|
||
n = atoi(argv[1]);
|
||
}
|
||
if (argc > 2) {
|
||
verbose = (strcmp(argv[2], "--verbose") == 0 || strcmp(argv[2], "-v") == 0);
|
||
}
|
||
|
||
cout << "N体问题串行模拟" << endl;
|
||
cout << "天体数量: " << n << endl;
|
||
cout << "时间步长: " << DT << " s" << endl;
|
||
cout << "总步数: " << TMAX << endl;
|
||
cout << "----------------------------------------" << endl;
|
||
|
||
// 初始化天体系统
|
||
vector<Body> bodies(n);
|
||
vector<Body> bodies_new(n);
|
||
|
||
// 初始化天体数据（简化版：模拟太阳系内行星）
|
||
// 使用简化的单位系统以便观察效果
|
||
double mass_scale = 1e24; // 质量缩放因子
|
||
double dist_scale = 1e8; // 距离缩放因子
|
||
double vel_scale = 1e3; // 速度缩放因子
|
||
|
||
// 中心天体（类似太阳）
|
||
bodies[0].mass = 1000 * mass_scale;
|
||
bodies[0].position = Vec3(0, 0, 0);
|
||
bodies[0].velocity = Vec3(0, 0, 0);
|
||
|
||
// 其他天体（类似行星）
|
||
for (int i = 1; i < n; i++) {
|
||
bodies[i].mass = (1.0 + i * 0.5) * mass_scale;
|
||
double angle = 2.0 * M_PI * i / n;
|
||
double radius = (1.0 + i * 0.5) * dist_scale;
|
||
|
||
bodies[i].position = Vec3(
|
||
radius * cos(angle),
|
||
radius * sin(angle),
|
||
0.0
|
||
);
|
||
|
||
// 给予切向速度以形成轨道
|
||
double orbital_speed = sqrt(G * bodies[0].mass / radius);
|
||
bodies[i].velocity = Vec3(
|
||
-orbital_speed * sin(angle),
|
||
orbital_speed * cos(angle),
|
||
0.0
|
||
);
|
||
}
|
||
|
||
// 输出初始状态
|
||
cout << fixed << setprecision(6);
|
||
if(verbose){
|
||
cout << "\n初始状态:" << endl;
|
||
for (int i = 0; i < n; i++) {
|
||
cout << "天体 " << i << ": 质量=" << bodies[i].mass/mass_scale << "e24 kg, "
|
||
<< "位置=(" << bodies[i].position.x/dist_scale << ", "
|
||
<< bodies[i].position.y/dist_scale << ", "
|
||
<< bodies[i].position.z/dist_scale << ")e8 m" << endl;
|
||
}
|
||
}
|
||
|
||
// 主循环：N体模拟
|
||
cout << "\n开始模拟..." << endl;
|
||
time_t start_time = clock();
|
||
for (int t = 0; t < TMAX; t++) {
|
||
// 第一步：计算所有物体新的速度和位置
|
||
for (int i = 0; i < n; i++) {
|
||
// 计算第i个物体所受的力
|
||
Vec3 F = compute_force(i, bodies);
|
||
|
||
// 计算新速度: v^(t+1) = v^t + F * dt / m
|
||
Vec3 v_new = bodies[i].velocity + F * DT / bodies[i].mass;
|
||
|
||
// 计算新位置: x^(t+1) = x^t + v^(t+1) * dt
|
||
Vec3 x_new = bodies[i].position + v_new * DT;
|
||
|
||
// 保存到临时数组
|
||
bodies_new[i].mass = bodies[i].mass;
|
||
bodies_new[i].position = x_new;
|
||
bodies_new[i].velocity = v_new;
|
||
}
|
||
|
||
// 第二步：更新所有物体的速度和位置
|
||
for (int i = 0; i < n; i++) {
|
||
bodies[i].position = bodies_new[i].position;
|
||
bodies[i].velocity = bodies_new[i].velocity;
|
||
}
|
||
|
||
// 每10步输出一次状态
|
||
if (verbose && (t + 1) % 10 == 0) {
|
||
cout << "时间步 " << t + 1 << ":" << endl;
|
||
for (int i = 0; i < n; i++) {
|
||
cout << " 天体 " << i << ": "
|
||
<< "位置=(" << bodies[i].position.x/dist_scale << ", "
|
||
<< bodies[i].position.y/dist_scale << ", "
|
||
<< bodies[i].position.z/dist_scale << ")e8 m, "
|
||
<< "速度=(" << bodies[i].velocity.x/vel_scale << ", "
|
||
<< bodies[i].velocity.y/vel_scale << ", "
|
||
<< bodies[i].velocity.z/vel_scale << ")e3 m/s" << endl;
|
||
}
|
||
}
|
||
}
|
||
time_t end_time = clock();
|
||
double elapsed_secs = double(end_time - start_time) / CLOCKS_PER_SEC;
|
||
cout << "\n模拟用时: " << elapsed_secs << " 秒" << endl;
|
||
cout << "\n模拟完成!" << endl;
|
||
|
||
return 0;
|
||
}
|