Моделирование N тела в С++

Я пытаюсь реализовать версию OpenMP двумерного моделирования n-body.

Но есть проблема: я предполагаю, что начальная скорость и ускорение каждой частицы равны нулю. Когда частицы сначала собираются вместе, они рассеиваются на высокой скорости и не собираются снова.

Это не похоже на Закон Ньютона, верно?

Может кто-нибудь объяснить, почему это происходит и как я могу исправить ошибку?

Вот часть моего кода:

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {

    # pragma omp for schedule(static)
        for ( i = 0; i < g_N; ++i ) {
            g_pv[i].a_x = 0.0;
            g_pv[i].a_y = 0.0;
            for ( j = 0; j < g_N; ++j ) {
                if (i == j)
                    continue;
                double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
            } 
            g_pv[i].v_x += period * g_pv[i].a_x;
            g_pv[i].v_y += period * g_pv[i].a_y;
        }
    # pragma omp for schedule(static)
        for ( int i = 0; i < g_N; ++i ) {
            g_pv[i].pos_x += g_pv[i].v_x * period;
            g_pv[i].pos_y += g_pv[i].v_y * period;
        }
    }
}

Не беспокойтесь о OpenMP, просто рассматривайте его как последовательную версию. OpenMP не сильно повлияет на результат.

Изменить: чтобы уточнить, вот весь код (в этой части могут быть некоторые ошибки, но описанная выше проблема должна произойти в приведенном выше разделе кода)

# include <iostream>
# include <fstream>
# include <iomanip>
# include <cmath>
# include <vector>
# include <cstdlib>
# include <omp.h>

# include <GL/glew.h>
# include <GL/freeglut.h>
# include <GL/gl.h>

using namespace std;

/* the size of the opengl window */
# define WIDTH 2000
# define HEIGHT 2000

/* define the global constants */
const double G = 6.67 * pow(10, -11);
// const double G = 6.67;
const double e = 0.00001;
const double period = 1;

/* define the structure of particle */
struct particle
{
    double m;
    double pos_x;
    double pos_y;
    double v_x;
    double v_y;
    double a_x;
    double a_y;

    particle(double m = 0, double pos_x = 0, double pos_y = 0, 
            double v_x = 0, double v_y = 0, double a_x = 0, double a_y = 0)
    {
        this->m         = m;
        this->pos_x     = pos_x;
        this->pos_y     = pos_y;
        this->v_x       = v_x;
        this->v_y       = v_y;
        this->a_x       = a_x;
        this->a_y       = a_y;
    }
};

/* define the global data */
int g_N;                        // number of particles
vector<particle> g_pv;          // particle vector

void setUp();
void update();
void display();

int main(int argc, char ** argv) {

    /* set up the window */
    glutInit(&argc, argv);
    glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
    glutInitWindowSize (WIDTH, HEIGHT);
    glutInitWindowPosition (100, 100);
    glutCreateWindow ("openmp");

    /* initialize */
    setUp();

    glutDisplayFunc(display);
    glutMainLoop();

    return 0;
}

/* read the input data */
void setUp() {
    glMatrixMode (GL_PROJECTION);
    glLoadIdentity ();
    /* Sets a 2d projection matrix
     * (0,0) is the lower left corner (WIDTH, HEIGHT) is the upper right */
    glOrtho (0, WIDTH, 0, HEIGHT, 0, 1);
    glDisable(GL_DEPTH_TEST);

    ifstream inFile;
    inFile.open("input_25.txt");

    inFile >> g_N;
    g_pv.resize(g_N);
    for ( int i = 0; i < g_N; ++i )
    {
        inFile >> g_pv[i].m >> g_pv[i].pos_x >> g_pv[i].pos_y
               >> g_pv[i].v_x >> g_pv[i].v_y >> g_pv[i].a_x >> g_pv[i].a_y; 
    }
    inFile.close();
}

/* display in openGL */
void display(void) {
    glClear(GL_COLOR_BUFFER_BIT);
    for(int i = 0; i < g_pv.size(); ++i) {
        /* Get the ith particle */
        particle p = g_pv[i];

        /* Draw the particle as a little square. */
        glBegin(GL_QUADS);
        glColor3f (1.0, 1.0, 1.0);
        glVertex2f(p.pos_x + 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y - 2);
        glVertex2f(p.pos_x + 2, p.pos_y - 2);
        glEnd();
    }   
    update();
    glutPostRedisplay();
    glFlush ();
}

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {
        /* compute the force */
        # pragma omp for schedule(static)
            for ( i = 0; i < g_N; ++i ) {
                g_pv[i].a_x = 0.0;
                g_pv[i].a_y = 0.0;
                for ( j = 0; j < g_N; ++j ) {
                    if (i == j)
                        continue;
                    double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                    g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                    g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
                } 
                g_pv[i].v_x += period * g_pv[i].a_x;
                g_pv[i].v_y += period * g_pv[i].a_y;
            }

        /* compute the velocity */
        # pragma omp for schedule(static)
            for ( int i = 0; i < g_N; ++i ) {
                g_pv[i].pos_x += g_pv[i].v_x * period;
                g_pv[i].pos_y += g_pv[i].v_y * period;
            }
    }
}