C++ SolarSystem::numberOfBodies方法代码示例

void RK4::integrate(std::valarray<double> &X, std::valarray<double> &V, double dt, SolarSystem mysystem, double G, double eps)
{
    std::valarray<double> k1(1,6*mysystem.numberOfBodies());
    std::valarray<double> k2(1,6*mysystem.numberOfBodies());
    std::valarray<double> k3(1,6*mysystem.numberOfBodies());
    std::valarray<double> k4(1,6*mysystem.numberOfBodies());

    // RK4 integration using vector X from solarysystem class.
    k1 = mysystem.calculateRK4(X, V, G, eps) * dt;
    k2 = mysystem.calculateRK4(X + 0.5 * k1, V, G, eps) * dt;
    k3 = mysystem.calculateRK4(X + 0.5 * k2, V, G, eps) * dt;
    k4 = mysystem.calculateRK4(X + k3, V, G, eps) * dt;
    X += (1.0/6) * (k1 + 2 * (k2 + k3) + k4);

}

void RK4::integrateSolarSystem(SolarSystem &system, double dt, int sunCheck)
{
    // Creating vectors
    int n_bodies = system.numberOfBodies();

    std::vector<double> A  = std::vector<double>(4*n_bodies);
    std::vector<double> K1 = std::vector<double>(4*n_bodies);
    std::vector<double> K2 = std::vector<double>(4*n_bodies);
    std::vector<double> K3 = std::vector<double>(4*n_bodies);
    std::vector<double> K4 = std::vector<double>(4*n_bodies);

    // Setting up A
    /* vector A = [x1, y1, vx1, vy1,    (first body)
                   x2, y2, vx2, vy2,    (second body)
                   ... for all celestial bodies
    */
    for(int i=0;i<n_bodies;i++)
    {
        CelestialBody *body = system.bodies[i];
        A[4*i]  = body->position.x();
        A[4*i+1]= body->position.y();
        A[4*i+2]= body->velocity.x();
        A[4*i+3]= body->velocity.y();
    }

    // Runge-Kutta 4th order integration, start
    K1= dAdt(system,A);

    // Returning acceleration values to CelestialBody objects
    for(int i=0;i<n_bodies;i++)
    {
        CelestialBody *body = system.bodies[i];
        body->acceleration.set(K1[4*i+2], K1[4*i+3], 0);
    }

    // Runge-Kutta 4th order integration, continued
    K1 = mult(K1,dt);
    K2 = dAdt(system,add(A,mult(K1,1/2.0))); K2 = mult(K2,dt);
    K3 = dAdt(system,add(A,mult(K2,1/2.0))); K3 = mult(K3,dt);
    K4 = dAdt(system,add(A,K3))            ; K4 = mult(K4,dt);

    // Combining (K1 + 2*K2 + 2*K3 + K4)/6 into K1
    K1 = add(K1,K4);
    K1 = add(K1,mult(K2,2));
    K1 = add(K1,mult(K3,2));
    K1 = mult(K1,1/6.0);
    A = add(A,K1);

    // Returning new position and velocity values to CelestialBody objects
    for(int i=sunCheck;i<n_bodies;i++)  // For stationary Sun, set startpoint i=1
    {
        CelestialBody *body = system.bodies[i];
        body->position.set(A[4*i],   A[4*i+1], 0);
        body->velocity.set(A[4*i+2], A[4*i+3], 0);
    }
}

std::vector<double> RK4::dAdt(SolarSystem &system, std::vector<double> A)
{
    double pi = 4*std::atan(1.0);
    double G = 4*pi*pi;
    int n_bodies = system.numberOfBodies();
    std::vector<double> dAdt = std::vector<double>(4*n_bodies);

    // Zeroing the vector dAdt
    for(int i=0;i<4*n_bodies;i++)
    {
        dAdt[i]=0;
    }

    // Derivating vector A
    for(int i=0;i<n_bodies;i++)
    {
        // Derivative of position is velocity
        dAdt[4*i] = A[4*i+2];
        dAdt[4*i+1] = A[4*i+3];

        double mass_i = system.bodies.at(i)->mass;
        for(int j=i+1; j<n_bodies; j++)
        {
            // Calculating gravitational force between objects
            double mass_j = system.bodies.at(j)->mass;
            double dX = A[4*j]-A[4*i];    // dX and dY "pointing" from body1 to body2
            double dY = A[4*j+1]-A[4*i+1];
            double dr = std::sqrt(dX*dX+dY*dY);

            double forcefactor = G*mass_i*mass_j/(dr*dr*dr);

            dAdt[4*i+2] += dX*forcefactor/mass_i;
            dAdt[4*i+3] += dY*forcefactor/mass_i;

            dAdt[4*j+2] -= dX*forcefactor/mass_j;
            dAdt[4*j+3] -= dY*forcefactor/mass_j;
        }
    }
return dAdt;
}

int main()
{
    // SET METHOD: 0 for Verlet and 1 for Runge-Kutta
    int method = 1;

    // SET SUN MOVEMENT ABILITY: 0 for moving sun, 1 for stationary sun
    int stationarySun = 1;

    // Set endpoint of time calculations and timestep (dt)
    float number_of_years = 20;
    double timestep = 1e-2;

    // Setting initial solar system and celestial bodies
    const double pi = 4*std::atan(1.0); // Pi with double-precision
    const double G = 4*pi*pi;

    // CelestialBody name(x,     y,vx, vy,               mass,    radius,  name);
    CelestialBody sun    (0,     0, 0, 0,                1,       4.64e-3, "Sun");
//    CelestialBody mercury(0.39,  0, 0, 2*pi/sqrt(0.39),  1.65e-7, 1.63e-5, "Mercury");
//    CelestialBody venus  (0.72,  0, 0, 2*pi/sqrt(0.72),  2.45e-6, 4.03e-5, "Venus");
    CelestialBody earth  (1,     0, 0, 2*pi,             3e-6,    4.25e-5, "Earth");
//    CelestialBody mars   (1.52,  0, 0, 2*pi/sqrt(1.52),  3.2e-7,  2.26e-5, "Mars");
//    CelestialBody jupiter(5.2,   0, 0, 2*pi/sqrt(5.2),   9.5e-4,  4.66e-4, "Jupiter");
//    CelestialBody saturn (9.54,  0, 0, 2*pi/sqrt(9.54),  2.85e-4, 3.88e-4, "Saturn");
//    CelestialBody uranus (19.19, 0, 0, 2*pi/sqrt(19.19), 4.35e-5, 1.69e-4, "Uranus");
//    CelestialBody neptune(30.07, 0, 0, 2*pi/sqrt(30.07), 5.1e-5,  1.64e-4, "Neptune");
//    CelestialBody pluto  (39.48, 0, 0, 2*pi/sqrt(39.48), 6.55e-9, 7.67e-6, "Pluto");

    SolarSystem solSyst;
    solSyst.addCelestialBody(sun);
//    solSyst.addCelestialBody(mercury);
//    solSyst.addCelestialBody(venus);
    solSyst.addCelestialBody(earth);
//    solSyst.addCelestialBody(mars);
//    solSyst.addCelestialBody(jupiter);
//    solSyst.addCelestialBody(saturn);
//    solSyst.addCelestialBody(uranus);
//    solSyst.addCelestialBody(neptune);
//    solSyst.addCelestialBody(pluto);



    RK4 solSystRK;


    int n_steps = number_of_years/timestep;              // Number of calculation points
    int n_bodies = solSyst.numberOfBodies();


    if(method==0)   // Verlet method
    {
        // Defining vectors for saving acceleration values
        vector<double> newAcceleration = vector<double>(2*n_bodies);
        vector<double> avgAcceleration = vector<double>(2*n_bodies);

        // avgAcceleration must be zero at start
        for(int i=0;i<2*n_bodies;i++)
        {
            avgAcceleration[i]=0;
        }

        for(int step=0;step<=n_steps;step++)
        {
            printf("Progress: %4.1f %% \r", 100.0*((double)step)/((double)n_steps));
            solSyst.dumpToFile(timestep, step);

            // Calculating forces for bodies in solar system
            solSyst.calculateForces();

            for (int i=stationarySun;i<n_bodies;i++)
            {
                CelestialBody *body = solSyst.bodies[i];

                // Fetching forces acting on body
                newAcceleration[2*i]   = G*body->force.x()/body->mass;
                newAcceleration[2*i+1] = G*body->force.y()/body->mass;

                // Setting acceleration as mean value of current step and previous step
                avgAcceleration[2*i]   = (avgAcceleration[2*i]   + newAcceleration[2*i])  /2.0;
                avgAcceleration[2*i+1] = (avgAcceleration[2*i+1] + newAcceleration[2*i+1])/2.0;

                // Returning acceleration to class parameter
                body->acceleration.set(avgAcceleration[2*i], avgAcceleration[2*i+1], 0);

                // Calculating new velocity and returning value to class parameter
                double velocity_x = body->velocity.x()+ avgAcceleration[2*i]  * timestep;
                double velocity_y = body->velocity.y()+ avgAcceleration[2*i+1]* timestep;
                body->velocity.set(velocity_x, velocity_y, 0);

                // Position = previous position + v*dt + 1/2 a*dt^2, returned to calss parameter
                double position_x = body->position.x()+ velocity_x * timestep + ( 0.5 * avgAcceleration[2*i]   * timestep * timestep);
                double position_y = body->position.y()+ velocity_y * timestep + ( 0.5 * avgAcceleration[2*i+1] * timestep * timestep);
                body->position.set(position_x, position_y, 0);

            }

//            cout << "Sum kinetic and potential energy: " << solSyst.totalEnergy()<<endl;

            // Collision test
            int colltest = collisionTest(solSyst);
//.........这里部分代码省略.........