Как запретить программе возвращать «inf»?

Я работаю над знаменитой задачей классической физики: - задачей трёх тел. Я составил следующую программу: -

#include <math.h>
#include <iostream>
using namespace std;

class System_variables
{
    private:

        // These are the system variables:- position, velocity and accleration
        long double body_1_pos[3];
        long double body_2_pos[3];
        long double body_3_pos[3];

        long double body_1_v[3];
        long double body_2_v[3];
        long double body_3_v[3];

        long double acc_body_1[3];
        long double acc_body_2[3];
        long double acc_body_3[3];

        // Needed constants
        const long double body_mass = 6 * pow(10, 24);
        double time_step;
        const long double G = 6.67 * pow(10, -11);

    public:


        System_variables(long double body_1_pos_i[3], long double body_2_pos_i[3], long double body_3_pos_i[3], long double body_1_v_i[3], long double body_2_v_i[3], long double body_3_v_i[3])
        {
            // Initiates the system variables, position and velocity.
            for (int m = 0; m<=2; m++)
            {
                body_1_pos[m] = body_1_pos_i[m];
                body_2_pos[m] = body_2_pos_i[m];
                body_3_pos[m] = body_3_pos_i[m];

                body_1_v[m] = body_1_v_i[m];
                body_2_v[m] = body_2_v_i[m];
                body_3_v[m] = body_3_v_i[m];
            }
        }
        void init_time_step(double step)
        {
            // Initializes time step.
            time_step = step;
        }
        void calculate_acc()
        {   
            // Calculates the accleration on the bodies using newtons equations.
            for (int m=0; m<=2; m++)
            {
                acc_body_1[m] = G * body_mass/pow((body_2_pos[m] - body_1_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_1_pos[m]), 2);
                acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);
                acc_body_3[m] = G * body_mass/pow((body_1_pos[m] - body_3_pos[m]), 2) + G * body_mass/pow((body_2_pos[m] - body_3_pos[m]), 2);
            };
            return;
        }

        void apply_acc_to_v()
        {   
            // Applys the accleration to the velocities of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_v[m] += acc_body_1[m] * time_step;
                body_2_v[m] += acc_body_2[m] * time_step;
                body_3_v[m] += acc_body_3[m] * time_step;
            }
            return;
        }

        void apply_v_to_pos()
        {   
            // Applys the velocities to the positions of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_pos[m] += body_1_v[m] * time_step;
                body_2_pos[m] += body_2_v[m] * time_step;
                body_3_pos[m] += body_3_v[m] * time_step;
            }
            return;
        }

        long double get_body_1_pos(int index)
        {
            // Returns the position of body 1 based on index.
            return body_1_pos[index];
        }

        long double get_body_2_pos(int index)
        {   
            // Returns the position of body 2 based on index.
            return body_2_pos[index];
        }

        long double get_body_3_pos(int index)
        {
            // Returns the position of body 3 based on index.
            return body_3_pos[index];
        }

        long double get_body_1_v(int index)
        {
            // Returns the velocity of body 1 based on index.
            return body_1_v[index];
        }

        long double get_body_2_v(int index)
        {
            // Returns the velocity of body 2 based on index.
            return body_2_v[index];
        }

        long double get_body_3_v(int index)
        {
            // Returns the velocity of body 3 based on index.
            return body_3_v[index];
        }


};

int main()
{
    long double array_1[3] = {100,100,100}; 
    long double array_2[3]  = {10000000000, 10000000000, 10000000000};
    long double array_3[3] = {10000000000, 10000000000, -10000000000};
    long double array_0[3] =  {0,0,0};

    System_variables system(array_1, array_2, array_3, array_0, array_0, array_0);
    system.init_time_step(100);

    for (int time = 0; time <= 10; time++)
    {
        cout << "x position of body 1: " << system.get_body_1_pos(0) << endl 
        << "y position of body 1: " << system.get_body_1_pos(1) << endl 
        << "z position of body 1: " << system.get_body_1_pos(2) << endl << endl;

        cout << "x position of body 2: " << system.get_body_2_pos(0) << endl 
        << "y position of body 2: " << system.get_body_2_pos(1) << endl 
        << "z position of body 2: " << system.get_body_2_pos(2) << endl << endl;

        cout << "x position of body 3: " << system.get_body_3_pos(0) << endl 
        << "y position of body 3: " << system.get_body_3_pos(1) << endl 
        << "z position of body 3: " << system.get_body_3_pos(2) << endl << endl;

        system.calculate_acc();
        system.apply_acc_to_v();
        system.apply_v_to_pos();
    }
}

Следует отметить, что вместо использования классического типа данных double я использовал long double, потому что, когда я использовал тип данных double, через некоторое время выводился inf, поэтому я использовал long double, думая, что он будет хранить десятичные дроби в больших размерах. точность. Но это не сработало. Я тоже пробовал float (хотя точность меньше), но результат всегда один и тот же: -

x position of body 1: 100
y position of body 1: 100
z position of body 1: 100

x position of body 2: 1e+10
y position of body 2: 1e+10
z position of body 2: 1e+10

x position of body 3: 1e+10
y position of body 3: 1e+10
z position of body 3: -1e+10

x position of body 1: 100.08
y position of body 1: 100.08
z position of body 1: 100.08

x position of body 2: inf
y position of body 2: inf
z position of body 2: 1e+10

x position of body 3: inf
y position of body 3: inf
z position of body 3: -1e+10

x position of body 1: 100.16
y position of body 1: 100.16
z position of body 1: 100.24

x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10

x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10

x position of body 1: -nan(ind)
y position of body 1: -nan(ind)
z position of body 1: 100.48

x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10

x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10

Я подумываю об использовании указателей, но сомневаюсь, что это правильная идея. Что вы думаете?