[C++&Py]三体运动粗略模拟

LuLuWuWei1120 · 发表于 2026-3-29 23:24

本帖最后由 LuLuWuWei1120 于 2026-3-30 01:35 编辑

"傻孩子们，快——跑——啊!”
看了三体，闲来无事，模拟一下三体运动。运动精度明显不高，但是看出物理规律来还是可以的。至于碰撞，在下无力回天，初版想得太简单，后来还是老老实实交给ai来做，还加了一大堆限制条件。推荐就按照我的配置来，效果还算行。可视化部分，因为没有学习，暂时用blender里的脚本编辑器(在blender里打开blenderClient.py)完成。使用socket在c++计算部分和python展示部分互相传数据。数学积分部分交给ai的。Interface界面循环可能还有问题，但是本人暂时没有时间修改，不影响单次使用。

预览图：

以下是代码：

***/include

#Vec3.h

[C++] 纯文本查看 复制代码

#pragma once
#include "Vec3.h"

namespace ThreeBody { 
    class Body {
    public:
        double radius;
        double mass;
        Vec3 position;
        Vec3 velocity;
        Vec3 acceleration;
        Vec3 oldAcceleration;
        Vec3 force;

        Body(double r, double m, const Vec3& position, const Vec3& velocity);
        void applyForce(const Vec3& f);
        void clearForce();
        void updateAcceleration();
        void saveOldAcceleration();
    };
}

#Config.h

[C++] 纯文本查看 复制代码

#pragma once

namespace ThreeBody {

        enum class IntegratorType
        {
                Euler,
                Verlet,
                RK4
        };

        struct Config 
{
        double G = 5.0;
        double dt = 0.1;
        double collisionRate = 1.0;

        //积分类型
        IntegratorType integratorType = IntegratorType::RK4;
        
        //是否开启随机速度
        bool randomVelocity = false;

};
}

#dataserver.h

[C++] 纯文本查看 复制代码

#pragma once
#include "UserInterface.h"
#include "Body.h"
#include <string>


int startTcpServer(unsigned short port, ThreeBody::UserInterface & ui);

void getData(ThreeBody::UserInterface & ui, std::string& data, bool& firstSend);

#Integrator.h

[C++] 纯文本查看 复制代码

#pragma once

#include <vector>
#include <functional>
#include "Body.h"

namespace ThreeBody{
//积分器基类
class Integrator{
    public:
    virtual ~Integrator() = default;
    virtual void integrate(std::vector<Body>& bodies,double dt,std::function<void(std::vector<Body>&)>recomputeAcceleration) = 0;
};

//显式欧拉积分器（一阶）
class EulerIntegrator : public Integrator{
    public:
    void integrate(std::vector<Body>& bodies,double dt,std::function<void(std::vector<Body>&)>recomputeAcceleration) override;
};

//Velocity Verlet积分器（二阶）
class VerletIntegrator : public Integrator{
    public:
    void integrate(std::vector<Body>& bodies,double dt,std::function<void(std::vector<Body>&)>recomputeAcceleration) override;
};

//RK4积分器（四阶）
class RK4Integrator : public Integrator{
    public:
    void integrate(std::vector<Body>& bodies,double dt,std::function<void(std::vector<Body>&)>recomputeAcceleration) override;
};

}

#PhysicsEngine.h

[C++] 纯文本查看 复制代码

#pragma once
#include "Body.h"
#include "Config.h"
#include <vector>

namespace ThreeBody { 
    class PhysicsEngine {
    public:
        explicit PhysicsEngine(const Config& config);
        void computeForce(std::vector<Body>& bodies) const;
        void collision(std::vector<Body>& bodies) const;
        
    private:
        Config myConfig;
    };
}

#Simulation.h

[C++] 纯文本查看 复制代码

#pragma once
#include <vector>
#include <memory>
#include <atomic>
#include <thread>
#include "Body.h"
#include "Config.h"
#include "PhysicsEngine.h"
#include "Integrator.h"

namespace ThreeBody { 
    class Simulation {
    public:
        Simulation(const std::vector<Body>& initialBodies, const Config& config,
                   std::unique_ptr<Integrator> integrator,
                   std::unique_ptr<PhysicsEngine> physicsEngine);
        ~Simulation();
        void start();
        void stop();
        void step();
        bool isRun() const { return running_; }
        const std::vector<Body>& getBodies() const { return bodies_; }
        
    private:
        void runloop();
        double betterDt() const; 
        
        std::vector<Body> bodies_;
        Config config_;
        std::unique_ptr<Integrator> integrator_;
        std::unique_ptr<PhysicsEngine> physicsEngine_;
        std::atomic<bool> running_;
        std::thread simulationThread_;
        double dt_;
    };
}

#UserInterface.h

[C++] 纯文本查看 复制代码

#pragma once

#include <memory>
#include "Simulation.h"
#include "Config.h"

namespace ThreeBody {

class UserInterface {
public:
    explicit UserInterface(const Config& config);
    void run();
    
    bool isRunning() const { return running_; }
    const std::vector<Body>& getBodies() const;

private:
    Config config_;
    bool running_ = false;
    Simulation* currentSimulation_ = nullptr;

    void randomMode();
    void stableMode();
    void customMode();

    std::unique_ptr<Simulation> createSimulation(const std::vector<Body>& bodies);
};

}

#Body.h

[C++] 纯文本查看 复制代码

#pragma once
#include "Vec3.h"

namespace ThreeBody { 
    class Body {
    public:
        double radius;
        double mass;
        Vec3 position;
        Vec3 velocity;
        Vec3 acceleration;
        Vec3 oldAcceleration;
        Vec3 force;

        Body(double r, double m, const Vec3& position, const Vec3& velocity);
        void applyForce(const Vec3& f);
        void clearForce();
        void updateAcceleration();
        void saveOldAcceleration();
    };
}

***/src

#blenderClient.py

[Python] 纯文本查看 复制代码

import socket
import bpy
import mathutils

# 全局变量
data_socket = None
spheres = []
first_receive = True
timer_running = False

IP = '127.0.0.1'
SERVER_PORT = 8888

def clear_scene():
    """清空场景中的所有物体"""
    bpy.ops.object.select_all(action='SELECT')
    bpy.ops.object.delete()

def create_spheres():
    """创建三个球体"""
    global spheres
    spheres = []
    for i in range(3):
        bpy.ops.mesh.primitive_uv_sphere_add(radius=0.5, location=(0, 0, 0))
        sphere = bpy.context.active_object
        sphere.name = f"Body_{i+1}"
        spheres.append(sphere)
    return spheres

def update_sphere_positions(positions):
    """更新球体位置"""
    global spheres
    for i, pos in enumerate(positions[:3]):
        if i < len(spheres):
            if len(pos) == 4:
                x, y, z, radius = pos
                spheres[i].location = mathutils.Vector([x, y, z])
                spheres[i].scale = (radius, radius, radius)
            else:
                spheres[i].location = mathutils.Vector(pos)

def parse_data(data_str):
    """解析数据，只返回最后一帧"""
    global first_receive
    batches = data_str.strip().split('\n')
    last_batch = None
    for batch in batches:
        if not batch.strip():
            continue
        batch_data = []
        has_radius = first_receive
        items = batch.split(';')
        for item in items:
            cleaned_item = item.strip()
            if not cleaned_item:
                continue
            coords = cleaned_item.split(',')
            if has_radius and len(coords) != 4:
                continue
            if not has_radius and len(coords) != 3:
                continue
            try:
                x = float(coords[0].strip())
                y = float(coords[1].strip())
                z = float(coords[2].strip())
                if has_radius:
                    radius = float(coords[3].strip())
                    batch_data.append([x, y, z, radius])
                else:
                    batch_data.append([x, y, z])
            except ValueError:
                continue
        if batch_data:
            last_batch = batch_data
            if first_receive:
                first_receive = False
    return last_batch

def cleanup():
    """清理资源"""
    global data_socket, timer_running
    timer_running = False
    if data_socket:
        try:
            data_socket.close()
        except:
            pass
        data_socket = None

def start_client():
    """启动客户端连接"""
    global data_socket, timer_running
    try:
        data_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
        data_socket.settimeout(5.0)
        data_socket.connect((IP, SERVER_PORT))
        data_socket.setblocking(False)
        timer_running = True
        print(f"已连接到服务器 {IP}:{SERVER_PORT}")
        return True
    except Exception as e:
        print(f"连接失败: {e}")
        return False

def socket_update_timer():
    """定时器回调函数"""
    global data_socket, timer_running
    if not timer_running or not data_socket:
        return 1.0
    try:
        try:
            recved = data_socket.recv(4096)
        except BlockingIOError:
            return 0.016  # 约60fps
        except Exception as e:
            print(f"接收错误: {e}")
            cleanup()
            return 1.0
        
        if not recved:
            print("服务器断开连接")
            cleanup()
            return 1.0
        
        data_str = recved.decode('utf-8')
        batch_data = parse_data(data_str)
        if batch_data and len(batch_data) >= 3:
            update_sphere_positions(batch_data[:3])
            bpy.context.view_layer.update()
        
        return 0.016  # 约60fps
    except Exception as e:
        print(f"定时器错误: {e}")
        cleanup()
        return 1.0

def register_timer():
    """注册Blender定时器"""
    bpy.app.timers.register(socket_update_timer)

def unregister_timer():
    """注销Blender定时器"""
    try:
        bpy.app.timers.unregister(socket_update_timer)
    except:
        pass


print("启动Blender客户端...")
clear_scene()
print("场景清空")
create_spheres()
if start_client():
    register_timer()
    print("客户端启动完成")
    print("脚本已启动，将在后台持续运行")
else:
    print("客户端启动失败")

# 保持脚本运行
def keep_running():
    """保持脚本运行"""
    import time
    while True:
        time.sleep(1.0)

# 注意：在Blender中，脚本会在执行完毕后停止
# 但定时器会继续运行，直到Blender关闭

#Body.cpp

[C++] 纯文本查看 复制代码

#include "Body.h"

namespace ThreeBody { 

Body::Body(double r, double m, const Vec3& p, const Vec3& v) 
    : radius(r), mass(std::max(m, 1.0)), position(p), velocity(v), 
      acceleration(0, 0, 0), oldAcceleration(0, 0, 0), force(0, 0, 0) {}

void Body::applyForce(const Vec3& f) { 
    force += f; 
}

void Body::clearForce() { 
    force = Vec3(0, 0, 0); 
}

void Body::updateAcceleration() { 
    double effectiveMass = std::max(mass, 1.0); 
    acceleration = force / effectiveMass; 
}

void Body::saveOldAcceleration() { 
    oldAcceleration = acceleration; 
}

} // namespace ThreeBody

#dataServer.cpp

[C++] 纯文本查看 复制代码

#include "dataServer.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <winsock2.h>
#include <ws2tcpip.h>
#include <thread>
#include <chrono>
#include <vector>
#include <string>

#pragma comment(lib, "ws2_32.lib")


void getData(ThreeBody::UserInterface& ui, std::string& data, bool& firstSend) {
    data.clear();
    if (!ui.isRunning()) {
        return;
    }

    std::vector<ThreeBody::Body> temp = ui.getBodies();
    for (size_t i = 0; i < temp.size(); ++i) {
        const auto& body = temp[i];
        if (firstSend) {
            // 首次发送时包含半径信息
            data += std::to_string(body.position.x) + "," + std::to_string(body.position.y) + "," + std::to_string(body.position.z) + "," + std::to_string(body.radius) + ";";
        } else {
            // 后续发送只包含位置信息
            data += body.position.toString() + ";";
        }
    }
    // 移除最后一个分号
    if (!data.empty() && data.back() == ';') {
        data.pop_back();
    }
    // 添加换行符
    data += "\n";
    
    // 首次发送后设置为false
    if (firstSend) {
        firstSend = false;
    }
}

int startTcpServer(unsigned short port, ThreeBody::UserInterface & ui) {
    // 1. 初始化Winsock
    WSADATA wsaData;
    int ret = WSAStartup(MAKEWORD(2, 2), &wsaData);
    if (ret != 0) {
        printf("WSAStartup失败: %d\n", ret);
        return -1;
    }

    // 2. 创建监听套接字
    SOCKET listenSocket = socket(AF_INET, SOCK_STREAM, 0);
    if (listenSocket == INVALID_SOCKET) {
        printf("socket失败: %d\n", WSAGetLastError());
        WSACleanup();
        return -1;
    }

    // 3. 绑定本地IP和端口
    struct sockaddr_in saddr;
    memset(&saddr, 0, sizeof(saddr));
    saddr.sin_family = AF_INET;
    saddr.sin_port = htons(port);
    saddr.sin_addr.s_addr = INADDR_ANY;

    ret = bind(listenSocket, (struct sockaddr*)&saddr, sizeof(saddr));
    if (ret == SOCKET_ERROR) {
        printf("bind失败: %d\n", WSAGetLastError());
        closesocket(listenSocket);
        WSACleanup();
        return -1;
    }

    // 4. 设置监听
    ret = listen(listenSocket, 10);
    if (ret == SOCKET_ERROR) {
        printf("listen失败: %d\n", WSAGetLastError());
        closesocket(listenSocket);
        WSACleanup();
        return -1;
    }

    printf("服务器启动，监听端口 %d，等待连接...\n", port);

    // 主循环，支持客户端重连
    while (true) {
        // 5. 接受客户端连接
        struct sockaddr_in caddr;
        int addrlen = sizeof(caddr);
        SOCKET acceptSocket = accept(listenSocket, (struct sockaddr*)&caddr, &addrlen);
        if (acceptSocket == INVALID_SOCKET) {
            printf("accept失败: %d\n", WSAGetLastError());
            continue;
        }

        // 6. 打印客户端信息
        char ip[32];
        inet_ntop(AF_INET, &caddr.sin_addr, ip, sizeof(ip));
        printf("客户端连接 - IP:%s, PORT:%d\n", ip, ntohs(caddr.sin_port));

        // 每次新连接重置firstSend标志
        bool firstSend = true;

        // 7. 通信循环
        while (true) {
            fd_set readfds;
            FD_ZERO(&readfds);
            FD_SET(acceptSocket, &readfds);

            // 设置超时时间0，立刻返回，实现非阻塞检测
            timeval timeout = { 0, 0 };

            // 检测套接字状态
            int selectRet = select(0, &readfds, NULL, NULL, &timeout);
            if (selectRet > 0) {
                // 客户端断开：recv返回0
                char buff[1];
                int recvRet = recv(acceptSocket, buff, 1, MSG_PEEK);
                if (recvRet <= 0) {
                    printf("客户端已断开连接\n");
                    break;
                }
            }

            std::string data;
            getData(ui, data, firstSend);

            if (!data.empty()) {
                int sent = send(acceptSocket, data.c_str(), data.size(), 0);
                if (sent < 0) {
                    printf("发送失败: %d\n", WSAGetLastError());
                    break;
                }
            }

            std::this_thread::sleep_for(std::chrono::milliseconds(15));
        }

        // 8. 关闭客户端套接字
        closesocket(acceptSocket);
        printf("等待新的客户端连接...\n");
    }

    // 9. 关闭监听套接字（理论上不会执行到这里）
    closesocket(listenSocket);
    WSACleanup();

    return 0;
}

#Integrator.cpp

[C++] 纯文本查看 复制代码

#include "Integrator.h"

namespace ThreeBody {

void EulerIntegrator::integrate(std::vector<Body>& bodies, double dt, std::function<void(std::vector<Body>&)> recomputeAcceleration) {
    recomputeAcceleration(bodies);
    for (auto& body : bodies) {
        body.velocity += body.acceleration * dt;
        body.position += body.velocity * dt;
    }
}

void VerletIntegrator::integrate(std::vector<Body>& bodies, double dt, std::function<void(std::vector<Body>&)> recomputeAcceleration) {
    for (auto& body : bodies) {
        body.velocity += body.acceleration * (dt * 0.5);
    }
    for (auto& body : bodies) {
        body.position += body.velocity * dt;
    }
    recomputeAcceleration(bodies);
    for (auto& body : bodies) {
        body.velocity += body.acceleration * (dt * 0.5);
    }
}

void RK4Integrator::integrate(std::vector<Body>& bodies, double dt, std::function<void(std::vector<Body>&)> recomputeAcceleration) {
    size_t n = bodies.size();
    std::vector<Vec3> origPos(n), origVel(n);
    for (size_t i = 0; i < n; ++i) {
        origPos[i] = bodies[i].position;
        origVel[i] = bodies[i].velocity;
    }

    recomputeAcceleration(bodies);
    std::vector<Vec3> a1(n);
    for (size_t i = 0; i < n; ++i) a1[i] = bodies[i].acceleration;

    for (size_t i = 0; i < n; ++i) {
        bodies[i].position = origPos[i] + origVel[i] * (dt * 0.5);
        bodies[i].velocity = origVel[i] + a1[i] * (dt * 0.5);
    }
    recomputeAcceleration(bodies);
    std::vector<Vec3> a2(n);
    for (size_t i = 0; i < n; ++i) a2[i] = bodies[i].acceleration;

    for (size_t i = 0; i < n; ++i) {
        bodies[i].position = origPos[i] + origVel[i] * (dt * 0.5);
        bodies[i].velocity = origVel[i] + a2[i] * (dt * 0.5);
    }
    recomputeAcceleration(bodies);
    std::vector<Vec3> a3(n);
    for (size_t i = 0; i < n; ++i) a3[i] = bodies[i].acceleration;

    for (size_t i = 0; i < n; ++i) {
        bodies[i].position = origPos[i] + origVel[i] * dt;
        bodies[i].velocity = origVel[i] + a3[i] * dt;
    }
    recomputeAcceleration(bodies);
    std::vector<Vec3> a4(n);
    for (size_t i = 0; i < n; ++i) a4[i] = bodies[i].acceleration;

    for (size_t i = 0; i < n; ++i) {
        Vec3 newVel = origVel[i] + (a1[i] + a2[i] * 2.0 + a3[i] * 2.0 + a4[i]) * (dt / 6.0);
        bodies[i].velocity = newVel;
        bodies[i].position = origPos[i] + (origVel[i] + newVel) * (dt * 0.5);
    }
}

}

#main.cpp

[C++] 纯文本查看 复制代码

#include "UserInterface.h"
#include "Config.h"
#include "Simulation.h"
#include <iostream>
#include <thread>
#include <chrono>
#include "dataServer.h"

//终端测试
#include <windows.h>
#include <locale.h>

void testConfig(ThreeBody::Config & c) {
    std::cout << "调试模式，检查功能" << std::endl;
    //打印config
    std::cout << "===== Config 参数 =====\n";
    std::cout << "G: " << c.G << std::endl;
    std::cout << "dt: " << c.dt << std::endl;
    std::cout << "collisionRate: " << c.collisionRate << std::endl;
    std::cout << "integratorType: " << static_cast<int>(c.integratorType) << std::endl;
    std::cout << "randomVelocity: " << (c.randomVelocity ? "开启" : "关闭") << std::endl;
    std::cout << "========================\n";
}

void printBody(ThreeBody::UserInterface& ui){
    while (true) {
        if (ui.isRunning()) {
            const auto& bodies = ui.getBodies();
            system("cls");
            std::cout << "===== 天体信息 =====\n";
            for (size_t i = 0; i < bodies.size(); ++i) {
                const auto& body = bodies[i];
                std::cout << "天体 " << i + 1 << "：\n";
                std::cout << "  位置: (" << body.position.x << ", " << body.position.y << ", " << body.position.z << ")\n";
                std::cout << "  速度: (" << body.velocity.x << ", " << body.velocity.y << ", " << body.velocity.z << ")\n";
                std::cout << "  加速度: (" << body.acceleration.x << ", " << body.acceleration.y << ", " << body.acceleration.z << ")\n";
                std::cout << "\n";
            }
           std::this_thread::sleep_for(std::chrono::milliseconds(100)); // 每100ms更新一次
        } 
    }
};



void testSimulation(ThreeBody::UserInterface & ui) {
    while (true) {
        if (ui.isRunning()) {
            // 模拟开始，显示天体信息
            const auto& bodies = ui.getBodies();
            system("cls");
            std::cout << "===== 天体信息 =====\n";
            for (size_t i = 0; i < bodies.size(); ++i) {
                const auto& body = bodies[i];
                std::cout << "天体 " << i + 1 << "：\n";
                std::cout << "  位置: (" << body.position.x << ", " << body.position.y << ", " << body.position.z << ")\n";
                std::cout << "  速度: (" << body.velocity.x << ", " << body.velocity.y << ", " << body.velocity.z << ")\n";
                std::cout << "  加速度: (" << body.acceleration.x << ", " << body.acceleration.y << ", " << body.acceleration.z << ")\n";
                std::cout << "\n";
            }
            std::this_thread::sleep_for(std::chrono::milliseconds(100)); // 每100ms更新一次
        } else {
            // 模拟未开始或已停止，等待
            std::this_thread::sleep_for(std::chrono::milliseconds(500));
        }
    }
}


int main() {
    //终端测试
    setlocale(CP_ACP, "zh_CN.UTF-8");
    SetConsoleOutputCP(65001);
    SetConsoleCP(65001);

    ThreeBody::Config config;
    ThreeBody::UserInterface ui(config);
    
    // 启动测试函数
    testConfig(config);
    
    // 启动数据服务器（在后台线程中运行）
    std::thread serverThread([&ui]() {
        startTcpServer(8888, ui);
    });
    serverThread.detach();
    
    // 启动测试函数（在后台线程中运行，模拟开始后显示信息）
    std::thread testThread([&ui]() {
        testSimulation(ui);
    });
    testThread.detach();
    
    // 运行用户界面
    ui.run();
    
    return 0;
}

#PhysicsEngine.cpp

[C++] 纯文本查看 复制代码

#include "PhysicsEngine.h"
#include <cmath>

namespace ThreeBody {

PhysicsEngine::PhysicsEngine(const Config& config) : myConfig(config) {}

void PhysicsEngine::computeForce(std::vector<Body>& bodies) const {
    for (auto& body : bodies) {
        body.clearForce();
    }

    const double G = myConfig.G;
    const double softening = 1e-6;

    for (size_t i = 0; i < bodies.size(); ++i) {
        for (size_t j = i + 1; j < bodies.size(); ++j) {
            Vec3 distVec = bodies[j].position - bodies[i].position;
            double distSq = distVec.lengthSquared();
            double dist = std::sqrt(distSq);
            if (dist < 1e-12) continue;

            Vec3 direction = distVec / dist;
            double forceMagnitude = (G * bodies[i].mass * bodies[j].mass) / (distSq + softening);

            Vec3 force = direction * forceMagnitude;

            bodies[i].applyForce(force);
            bodies[j].applyForce(-force);
        }
    }
}

void PhysicsEngine::collision(std::vector<Body>& bodies) const {
    double e = myConfig.collisionRate;
    const double maxVelocity = 100.0;
    const double epsilon = 1e-12;
    
    for (size_t i = 0; i < bodies.size(); ++i) {
        for (size_t j = i + 1; j < bodies.size(); ++j) {
            Vec3 delta = bodies[j].position - bodies[i].position;
            double distSq = delta.lengthSquared();
            double dist = std::sqrt(distSq);
            double minDist = bodies[i].radius + bodies[j].radius;

            if (dist < minDist) {
                // 计算法向量，确保不会除零
                Vec3 n(1, 0, 0); // 默认方向
                if (distSq > epsilon) {
                    n = delta / dist;
                }
                
                double vRel = (bodies[j].velocity - bodies[i].velocity).dot(n);

                if (vRel > 0) continue;

                double m1 = bodies[i].mass, m2 = bodies[j].mass;
                // 确保质量不为零
                if (m1 < epsilon || m2 < epsilon) continue;
                
                double invMass1 = 1.0 / m1;
                double invMass2 = 1.0 / m2;
                double impulseMagnitude = -(1 + e) * vRel / (invMass1 + invMass2);
                Vec3 impulse = n * impulseMagnitude;

                // 应用冲量
                bodies[i].velocity -= impulse * invMass1;
                bodies[j].velocity += impulse * invMass2;

                // 限制速度，避免数值不稳定
                double velLength1 = bodies[i].velocity.length();
                if (velLength1 > maxVelocity && velLength1 > epsilon) {
                    bodies[i].velocity = bodies[i].velocity / velLength1 * maxVelocity;
                }
                double velLength2 = bodies[j].velocity.length();
                if (velLength2 > maxVelocity && velLength2 > epsilon) {
                    bodies[j].velocity = bodies[j].velocity / velLength2 * maxVelocity;
                }

                // 精确的位置校正，确保球体完全分离
                double overlap = minDist - dist;
                if (overlap > 0) {
                    Vec3 correction = n * (overlap * 0.5);
                    bodies[i].position -= correction;
                    bodies[j].position += correction;
                }
            }
        }
    }
}

}

#Simulation.cpp

[C++] 纯文本查看 复制代码

#include "Simulation.h"
#include <iostream>
#include <thread>
#include <chrono>

namespace ThreeBody {

Simulation::Simulation(const std::vector<Body>& initialBodies, const Config& config,
                      std::unique_ptr<Integrator> integrator,
                      std::unique_ptr<PhysicsEngine> physicsEngine)
    : bodies_(initialBodies), config_(config), integrator_(std::move(integrator)),
      physicsEngine_(std::move(physicsEngine)), running_(false), dt_(config_.dt) {
}

Simulation::~Simulation() {
    stop();
}

void Simulation::start() {
    if (running_) return;
    running_ = true;
    simulationThread_ = std::thread(&Simulation::runloop, this);
}

void Simulation::stop() {
    running_ = false;
    if (simulationThread_.joinable()) simulationThread_.join();
}

void Simulation::runloop() {
    while (running_) {
        step();
        std::this_thread::sleep_for(std::chrono::milliseconds(10));
    }
}

void Simulation::step() {
    integrator_->integrate(bodies_, dt_, [this](std::vector<Body>& bodies) {
        physicsEngine_->computeForce(bodies);
        for (auto& body : bodies) {
            body.updateAcceleration();
        }
    });
    physicsEngine_->collision(bodies_);
}

}

#UserInterface.cpp

[C++] 纯文本查看 复制代码

#include "UserInterface.h"
#include "Simulation.h"
#include "Integrator.h"
#include <iostream>
#include <random>
#include <memory>
#include <limits>

namespace ThreeBody {

    UserInterface::UserInterface(const Config& config) : config_(config) {}

    void UserInterface::run() {
        int choice = -1;
        while (choice != 4) {
            std::cout << "\n===== 三体模拟器 =====\n";
            std::cout << "1. 随机模式\n";
            std::cout << "2. 稳定模式（特殊解）\n";
            std::cout << "3. 自定义模式\n";
            std::cout << "4. 退出\n";
            std::cout << "请选择: ";
            
            // 处理输入错误
            if (!(std::cin >> choice)) {
                std::cin.clear(); // 重置输入状态
                std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n'); // 清除输入缓冲区
                std::cout << "无效输入，请重试。\n";
            } else {
                switch (choice) {
                case 1: randomMode(); break;
                case 2: system("cls"); stableMode(); break;
                case 3: customMode(); break;
                case 4: std::cout << "再见！\n"; break;
                default: std::cout << "无效输入，请重试。\n"; break;
                }
            }
            system("pause");
            system("cls");
        }
    }

    std::unique_ptr<Simulation> UserInterface::createSimulation(const std::vector<Body>& bodies) {
        std::unique_ptr<Integrator> integrator;
        switch (config_.integratorType) {
        case IntegratorType::Euler: integrator = std::make_unique<EulerIntegrator>(); break;
        case IntegratorType::Verlet: integrator = std::make_unique<VerletIntegrator>(); break;
        case IntegratorType::RK4: integrator = std::make_unique<RK4Integrator>(); break;
        }
        auto physicsEngine = std::make_unique<PhysicsEngine>(config_);
        auto sim = std::make_unique<Simulation>(bodies, config_, std::move(integrator), std::move(physicsEngine));
        currentSimulation_ = sim.get();
        running_ = true;
        return sim;
    }
    
    const std::vector<Body>& UserInterface::getBodies() const {
        static std::vector<Body> emptyBodies;
        if (currentSimulation_) {
            return currentSimulation_->getBodies();
        }
        return emptyBodies;
    }

    void UserInterface::randomMode() {
        std::random_device rd;
        std::mt19937 gen(rd());
        std::uniform_real_distribution<> massDist(1.0, 10.0);
        std::uniform_real_distribution<> radiusDist(0.5, 3.0);
        std::uniform_real_distribution<> posDist(-40.0, 40.0);
        std::uniform_real_distribution<> velDist(-2.0, 2.0);

        std::vector<Body> bodies;
        for (int i = 0; i < 3; ++i) {
            double mass = massDist(gen);
            double radius = radiusDist(gen);
            Vec3 position(posDist(gen), posDist(gen), posDist(gen));
            Vec3 velocity(0, 0, 0);
            if (config_.randomVelocity) {
                velocity = Vec3(velDist(gen), velDist(gen), velDist(gen));
            }
            bodies.emplace_back(radius, mass, position, velocity);
        }

        auto sim = createSimulation(bodies);
        sim->start();

        std::cout << "随机模式已启动，按 Enter 键停止模拟...\n";
        std::cin.ignore();
        std::cin.get();
        sim->stop();
        running_ = false;
        currentSimulation_ = nullptr;
    }

    void UserInterface::stableMode() {
        std::cout << "\n===== 稳定模式 =====\n";
        std::cout << "1.拉格朗日等边三角形解\n";
        std::cout << "2.8字形轨道解\n";
        std::cout << "3.拉格朗日点L4/L5解(限制性三体）\n";
        std::cout << "请选择: ";
        int choice;
        std::cin >> choice;

        std::vector<Body> bodies;
        double G = config_.G;
        switch (choice) {
        case 1: {
            double m = 1.0;
            double r = 0.2;
            double side = 5.0;
            double v = std::sqrt(G * m / side);
            bodies = {
                Body(r, m, Vec3(0, 0, 0), Vec3(-v / 2, v * std::sqrt(3) / 2, 0)),
                Body(r, m, Vec3(side, 0, 0), Vec3(-v / 2, -v * std::sqrt(3) / 2, 0)),
                Body(r, m, Vec3(side / 2, side * std::sqrt(3) / 2, 0), Vec3(v, 0, 0))
            };
            break;
        }
        case 2: {
            bodies = {
                Body(0.1, 1.0, Vec3(0.347111, 0.532728, 0), Vec3(0.932813, 0.864731, 0)),
                Body(0.1, 1.0, Vec3(-0.347111,-0.532728, 0), Vec3(0.932813, 0.864731, 0)),
                Body(0.1, 1.0, Vec3(0.000000, 0.000000, 0), Vec3(-1.865626,-1.729462, 0))
            };
            break;
        }
        case 3: {
            double M = 10.0;
            double m = 0.1;
            double R = 5.0;
            double omega = std::sqrt(G * M / (R * R * R));
            Vec3 pos1(-R * M / (2 * M), 0, 0);
            Vec3 pos2(R * M / (2 * M), 0, 0);
            double L4_x = 0;
            double L4_y = R * std::sqrt(3) / 2;
            bodies = {
                Body(1.0, M, pos1, Vec3(0, -omega * pos1.x, 0)),
                Body(1.0, M, pos2, Vec3(0, -omega * pos2.x, 0)),
                Body(0.1, m, Vec3(L4_x, L4_y, 0), Vec3(-omega * L4_y, omega * L4_x, 0))
            };
            break;
        }
        default:
            return;
        }

        auto sim = createSimulation(bodies);
        sim->start();
        std::cout << "模拟运行中，按回车停止...\n";
        std::cin.ignore();
        std::cin.get();
        sim->stop();
        running_ = false;
        currentSimulation_ = nullptr;
    }
    void UserInterface::customMode() {
        std::vector<Body> bodies;
        for (int i = 0; i < 3; ++i) {
            std::cout << "\n设置第 " << (i + 1) << " 个天体:\n";
            double radius, mass, posX, posY, posZ, velX, velY, velZ;
            std::cout << "半径: "; std::cin >> radius;
            std::cout << "质量: "; std::cin >> mass;
            std::cout << "位置 (x y z): "; std::cin >> posX >> posY >> posZ;
            std::cout << "速度 (vx vy vz): "; std::cin >> velX >> velY >> velZ;
            bodies.emplace_back(radius, mass, Vec3(posX, posY, posZ), Vec3(velX, velY, velZ));
        }

        auto sim = createSimulation(bodies);
        sim->start();
        std::cout << "模拟运行中，按回车停止...\n";
        std::cin.ignore();
        std::cin.get();
        sim->stop();
        running_ = false;
        currentSimulation_ = nullptr;
    }

}

#Vec3.cpp

[C++] 纯文本查看 复制代码

#include "Vec3.h"
#include <iostream>

namespace ThreeBody {
    std::ostream& operator<<(std::ostream& os, const Vec3& v) {
        os << "(" << v.x << "," << v.y << "," << v.z << ")";
        return os;
    }
}

dickwolf · 发表于 2026-3-30 17:34

高手，可以考虑防三体做游戏了策略游戏

shenqihailuo · 发表于 2026-4-1 16:19

厉害，可以做个动画看看

luckfalconguy · 发表于 2026-4-1 18:52

你的脑洞也可以，能想到这个设计实现也不容易。

帐号		自动登录	找回密码
密码			注册[Register]

[C&C++ 原创] [C++&Py]三体运动粗略模拟

免费评分

浏览过的版块