N체 시뮬레이션

#include <iostream>
#include <vector>
#include <list>
#include <fstream>
#include <string>
#include <boost/numeric/odeint.hpp>
#include "point.h"

typedef point<double, 3> Vector;
typedef std::vector<Vector> container;
typedef std::vector<double> mass;

constexpr double G_const = 6.67430e-11;

struct position
{
    const mass &masses;

    position(const mass &masses) : masses(masses)
    {
    }

    void operator()(const container &p, container &dqdt) const
    {
        for (int l = 0; l < p.size(); ++l)
        {
            dqdt[l] = p[l] / masses[l];
        }
    }
};

struct momentum
{
    const mass &masses;

    momentum(const mass &masses) : masses(masses)
    {
    }

    void operator()(const container &q, container &dpdt) const
    {
        for (int l = 0; l < q.size(); ++l)
        {
            dpdt[l] = 0.0;

            for (int m = 0; m < l; ++m)
            {
                Vector dist = q[m] - q[l];
                double r = abs(dist);
                dist *= (G_const * masses[l] * masses[m] / std::pow(r, 3));
                dpdt[l] += dist;
                dpdt[m] -= dist;
            }
        }
    }
};

Vector COM(const container &x, const mass &m)
{
    double m_all = 0.0;
    Vector temp;

    for (int l = 0; l < x.size(); ++l)
    {
        m_all += m[l];
        temp += m[l] * x[l];
    }

    temp /= m_all;
    return temp;
}

int main()
{
    mass masses = {1.98850e30,      // sun
                   3.302e23,     //mercury
                   4.8685e24,   //venus
                   5.97e24,     //earth
                   6.4171e23,
                   1.898e27,  // jupiter
                   5.683e26,  // saturn
                   8.681e25, // uranus
                   1.024e26, // neptune
                   1.30900e22      // pluto
    };

    container q = {Vector(0.0, 0.0, 0.0),                        // sun
                   Vector(-4.429040230571739e10, -5.122487388712578e10, -1.231269594356008e08), //mercury
                   Vector(9.433420726785934e10, -5.406597866002910e10, -6.185568261650380e9),   //venus
                   Vector(-1.421009056466009e11,  4.252331579033457e10, -1.779123675208539e6),
                   Vector(-2.902570650357215e10,  2.354782966418272e11,  5.646696531297803e9),
                   Vector(5.083257248705834e11, -5.638806876803797e11, -9.030947426587164e9),   // jupiter
                   Vector(8.612728194772711e11, -1.218887557008212e12, -1.308599067382348e10),    // saturn
                   Vector(2.272952271035277e12, 1.890758866689061e12, -2.242502329388785e10),   // uranus
                   Vector(4.412093131015528e12, -7.537562554024563e11, -8.616978441277650e10), // neptune
                   Vector(2.130690844336437e12, -4.655371851305727e12, -1.179723791735542e11)   // pluto
    };

    container p = {Vector(0.0, 0.0, 0.0),                        // sun
                   Vector(2.696326924444061e4, -2.961236470789841e4, 4.893179724867140e3),  //mercury
                   Vector(1.720638694010401e4, 3.023866050286396e4, -5.779240550025033e2),    //venus
                   Vector(-9.034004116118341e3, -2.864142202244999e4,  2.360353661973846e0),
                   Vector(-2.313108968913745e4, -9.080084179894222e2,  5.483935852274529e2),
                   Vector(9.557244459356859e3, 9.374665522774988e3, -2.527570523608271e2), // jupiter
                   Vector(7.363416859553132e3, 5.560644251357978e3, -3.899190600102365e2),   // saturn
                   Vector(-4.396393078151773e3, 4.927180734879181e3, 7.504538606204814e1),   // uranus
                   Vector(8.900553319463884e2, 5.400725982881951e3, -1.323001903350283e2),   // neptune
                   Vector(5.075950613888891e3, 1.112771293399162e3, -1.561804726683429e3)   // pluto
    };//EPHEMERIS 2021/03/04

    Vector Qm = COM(q, masses);
    Vector Pm = COM(p, masses);

    for (int l = 0; l < p.size(); ++l)
    {
        q[l] -= Qm;
        p[l] -= Pm;
    }

    for (int l = 0; l < p.size(); ++l)
        p[l] *= masses[l];

    boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan<container> rkn;

    double max_t = 1000000;
    double dt = 100;
    double t = 0.0;

    std::vector<std::list<Vector>> data(q.size());
    for (int i = 0; i < (max_t / dt); i++)
    {
        rkn.do_step(std::make_pair(position(masses), momentum(masses)), std::make_pair(boost::ref(q), boost::ref(p)), t,
                    dt);
        t += dt;

        for (int l = 0; l < q.size(); ++l)
            data.at(l).push_back(q.at(l));
    }

    std::ofstream out0("sun.txt");
    std::ofstream out1("mercury.txt");
    std::ofstream out2("venus.txt");
    std::ofstream out3("earth.txt");
    std::ofstream out4("mars.txt");
    std::ofstream out5("jupiter.txt");
    std::ofstream out6("saturn.txt");
    std::ofstream out7("uranus.txt");
    std::ofstream out8("neptune.txt");
    std::ofstream out9("pluto.txt");

    for (auto &x : data.at(0))
        out0 << x << std::endl;

    for (auto &x : data.at(1))
        out1 << x << std::endl;

    for (auto &x : data.at(2))
        out2 << x << std::endl;

    for (auto &x : data.at(3))
        out3 << x << std::endl;

    for (auto &x : data.at(4))
        out4 << x << std::endl;

    for (auto &x : data.at(5))
        out5 << x << std::endl;

    for (auto &x : data.at(6))
        out6 << x << std::endl;

    for (auto &x : data.at(7))
        out7 << x << std::endl;

    for (auto &x : data.at(8))
        out8 << x << std::endl;

    for (auto &x : data.at(9))
        out9 << x << std::endl;
}

마지막에 좀 노가다가 있긴 하지만 귀찮아서 그런거고

Solar system - 1.66.0

여기를 참조하였다.

 

천체 정보는

HORIZONS Web-Interface

여기서 볼 수 있다.

 

심플렉틱 적분기를 직접 만들어보려 했으나 부동소수점 오류가 발생하는 것으로 보여서 boost::odeint를 사용했다.

 

SFML로 시각화한 버전도 만들어 볼 예정이다.

OpenGL 버전도 있지만 입력기 문제 때문에 처음부터 다시 뜯어고쳐야 해서 더 진행할지는 잘 모르겠다.

 

이걸 어떻게 해야 멀티코어로 계산할 수 있을까?