c++ движение планеты в системе двойной звезды или движение трех тел

[anpolol]

программа должна выводить координаты движения звезд и планеты в файл (рисую в gnuplot).
во-первых, не понимаю, почему выходит, что все три объекта движутся прямолинейно (в численном методе Верле я не до конца разобралась, похоже что ошибка там, потому что формула как для прямолинейного движения, хотя препод говорил, что "вроде верно"), во-вторых, не очень понимаю какие начальные условия задавать, в каком отношении друг к другу (ну например скорость и координаты и масса, в каком отношении они должны быть?). Еще препод сказал, что начальные векторы скорости должны быть по касательной к траектории, то есть насколько я понимаю, перпендикулярно к прямой, соединяющей две звезды, я вроде так и сделала, но все равно все неправильно выводит, может неправильно я и сделала. В общем, не понимаю, что делать(

#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
#include <ctime>
#include<time.h>
#include<stdlib.h>

using namespace std;


double const G = 6;
struct tVector {
  double x,y,z, M;
};
class TWOSTARS
{
public:

    double F1x, F2x, F1y, F2y, F1z, F2z, ex1, ex2, ey1, ey2, ez1, ez2, exs, eys, ezs, FSy, FSz, FSx;
    double rast1, rast2, rasts;


    tVector star1;//координаты первой звезды
    tVector star2;
    tVector Vstar1;//координаты скороти планеты звезды
    tVector Vstar2;
    tVector planet;
    tVector Vplanet;
    tVector speedupS1;//ускорение 1 звезды
    tVector speedupS2;
    tVector speedupPlanet;//ускорение планеты



    TWOSTARS ()//начальные координаты
    {
         star1.x=15;
        star1.x=15;
        star1.z= 15;

        star2.x=20;
        star2.y= 20 ;
        star2.z= 20 ;

        Vstar1.x=1 ;
        Vstar1.y=1 ;
        Vstar1.z=-2 ;


        Vstar2.x= -2;
        Vstar2.y= -2;
        Vstar2.z= 4;




        planet.x= 500;
        planet.y= 500;
        planet.z= 500;

        Vplanet.x= 1;
        Vplanet.y=1 ;
        Vplanet.z= 1;
        star1.M= 20;
        star2.M= 20;
        planet.M= 5;
    }


    void changeforce() //сила гравитационного взаимодействия между объектами
    {

        rast1 = sqrt( (planet.x-star1.x)*(planet.x-star1.x)+(planet.y-star1.y)*(planet.y-star1.y)+(planet.z-star1.z)*(planet.z-star1.z));// расстояние между планетой и первой звездой
        rast2 = sqrt( (planet.x-star2.x)*(planet.x-star2.x)+(planet.y-star2.y)*(planet.y-star2.y)+(planet.z-star2.z)*(planet.z-star2.z));// расстояние между планетой и второй звездой
        rasts =  sqrt( (star1.x-star2.x)*(star1.x-star2.x)+(star1.y-star2.y)*(star1.y-star2.y)+(star1.z-star2.z)*(star1.z-star2.z));// расстояние между звездами

        ex1=(planet.x-star1.x)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось x, нужна чтоб у проекций сил гравитации было направление
        ey1=(planet.y-star1.y)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось y
        ez1=(planet.z-star1.z)/rast1;//проекция вектора соединяющего положение планеты и звезды 1 на ось z

        ex2=(planet.x-star2.x)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось x
        ey2=(planet.y-star2.y)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось y
        ez2=(planet.z-star2.z)/rast2;//проекция вектора соединяющего положение планеты и звезды 2 на ось z

        exs=(star1.x-star2.x)/rasts;//проекция вектора соединяющего положение звезд на ось x
        eys=(star1.y-star2.y)/rasts;//проекция вектора соединяющего положение звезд  на ось y
        ezs=(star1.z-star2.z)/rasts;//проекция вектора соединяющего положение звезд на ось z

        F1x=(G*star1.M*planet.M*ex1)/((planet.x-star1.x)*(planet.x-star1.x));//проекция силы гравитации между платнетой и звездой 1 на ось х
        F1y=(G*star1.M*planet.M*ey1)/((planet.y-star1.y)*(planet.y-star1.y));//проекция на ось y
        F1z=(G*star1.M*planet.M*ez1)/((planet.z-star1.z)*(planet.z-star1.z));//проекция на ось z

        F2x=(G*star2.M*planet.M*ex2)/((planet.x-star2.x)*(planet.x-star2.x));//проекция силы гравитации между платнетой и звездой 2 на ось х
        F2y=(G*star2.M*planet.M*ey2)/((planet.y-star2.y)*(planet.y-star2.y));//проекция  на ось y
        F2z=(G*star2.M*planet.M*ez2)/((planet.z-star2.z)*(planet.z-star2.z));//проекция  на ось z

        FSx=(G*star1.M*star2.M*exs)/((star1.x-star2.x)*(star1.x-star2.x));//проекция силы гравитации между звездами  на ось х
        FSy=(G*star1.M*star2.M*eys)/((star1.y-star2.y)*(star1.y-star2.y));//проекция  на ось y
        FSz=(G*star1.M*star2.M*ezs)/((star1.z-star2.z)*(star1.z-star2.z));//проекция  на ось z

    }
    void speedup()//просчитывается ускорение
    {
        speedupS1.x= (F1x-FSx)/star1.M;
        speedupS1.y= (F1y-FSy)/star1.M;
        speedupS1.z=(F1z-FSz)/star1.M;

        speedupS2.x= (F2x+FSx)/star2.M;
        speedupS2.y= (F2x+FSx)/star2.M;
        speedupS2.z=(F2x+FSx)/star2.M;

        speedupPlanet.x= (-F1x-F2x)/planet.M;
        speedupPlanet.y= (-F1x-F2x)/planet.M;
        speedupPlanet.z=(-F1x-F2x)/planet.M;


    }

    void SaveFileStar1(char filename[])//сохранение в файл координат
    {
        ofstream fout(filename, ios::app);
        fout<< setw(20) << star1.x << setw(20) << star1.y << setw(20) << star1.z << "\n";
    }

    void SaveFileStar2(char filename[])
    {
        ofstream fout(filename, ios::app);
        fout<< setw(20) << star2.x << setw(20) << star2.y << setw(20) << star2.z << "\n";
    }

    void SaveFilePlanet(char filename[])
    {
        ofstream fout(filename, ios::app);
        fout<< setw(20) << planet.x << setw(20) << planet.y << setw(20) << planet.z << "\n";
    }


    void coordinates()//изменение координат и скорости тел
    {
        double const step=1;//шаг по времени

        int time = 1000;//все время в течение которого рассматриваем движение
        int amount = (time/step);//сколько шагов нужно сделать в цикле

        for (int t=0; t<amount; t++)
        {
            changeforce();//изменяем силу
            speedup();//находим ускорение

            Vstar1.x+=speedupS1.x*step;
            Vstar1.y+=speedupS1.y*step;
            Vstar1.z+=speedupS1.z*step;

            Vstar2.x+=speedupS2.x*step;
            Vstar2.y+=speedupS2.y*step;
            Vstar2.z+=speedupS2.z*step;

            Vplanet.x+=speedupPlanet.x*step;
            Vplanet.y+=speedupPlanet.y*step;
            Vplanet.z+=speedupPlanet.z*step;


            star1.x+=Vstar1.x*step+0.5*speedupS1.x*step*step;
            star1.y+=Vstar1.y*step+0.5*speedupS1.y*step*step;
            star1.z+=Vstar1.z*step+0.5*speedupS1.z*step*step;

            star2.x+=Vstar2.x*step+0.5*speedupS2.x*step*step;//координаты звезды 2
            star2.y+=Vstar2.y*step+0.5*speedupS2.y*step*step;
            star2.z+=Vstar2.z*step+0.5*speedupS2.z*step*step;

            planet.x+=Vplanet.x*step+0.5*speedupPlanet.x*step*step;//координаты планеты
            planet.y+=Vplanet.y*step+0.5*speedupPlanet.y*step*step;
            planet.z+=Vplanet.z*step+0.5*speedupPlanet.z*step*step;

            SaveFileStar1("Star1"); //сохранение в файл координат
            SaveFileStar2("Star2");
            SaveFilePlanet("Planet");

        }
    }

};





int main()
{
    ofstream Star1("Star1", ios::trunc); //чистим файл если там что-то есть
    Star1.close(); //закрываем его
    ofstream Star2("Star2", ios::trunc);
    Star2.close();
    ofstream Planet("Planet", ios::trunc);
    Planet.close();

    TWOSTARS planet;//используем наш класс
    planet.coordinates();//вызываем в классе функцию которая просчитывает координаты

    return 0;
}

[type_Oleg]

Сообщение от anpolol

.. Еще препод сказал, что начальные векторы скорости должны быть по касательной к траектории..

Это всегда так. Любой вектор скорости , в любой момент времени - по касательной к траектории. На то он и вектор скорости.
А вот это -

Сообщение от anpolol

.. то есть насколько я понимаю, перпендикулярно к прямой, соединяющей две звезды ..

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

У вас типичная задача 3-х тел.
Создайте структуру - " вектор " , состоящую из 3-х переменных double

struct tVector {
  double X,Y,Z;
};

И у вас - координаты ( векторы положения) обеих звезд, координаты планеты, скорости обеих звезд, скорости планеты, ускорения звезд и планеты - все будет типа tVector

Можете создать другую, более полную, с массой тела

struct tCorps {  
  double X,Y,Z, M;
};

Создайте функции для
- расчета ускорения тела от притяжения 2-х остальных
- расчета изменения координат и скорости тела от этого ускорения

И программа будет более компактной, более удобочитаемой и более удобной в отладке.

PS Я делал что-то подобное, только не на Си, а на Паскале, и не по методу Верле, а по более крутому - Рунге-Кутты..

[anpolol]

я все переделала, только вот этого не поняла:

Сообщение от type_Oleg

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

потому что слово "барицентрический" узнала только сегодня и если я понимаю его определение верно, то вроде , мой ответ "да"

Сообщение от type_Oleg

PS Я делал что-то подобное, только не на Си, а на Паскале, и не по методу Верле, а по более крутому - Рунге-Кутты..

я метод Верле с трудом понимаю, а вы про Рунге-Кутта(

[Smitt&Wesson]

Сообщение от anpolol

я все переделала, только вот этого не поняла:

потому что слово "барицентрический" узнала только сегодня и если я понимаю его определение верно, то вроде , мой ответ "да"


я метод Верле с трудом понимаю, а вы про Рунге-Кутта(

Не лезь в барицентрический, это тебе не нужно. Это - высшая математика. Там, дальше - трансцендентые уравнения высшего порядка. А ещё дальше, вообще - пипец.

И вообще, anpolol, думаешь, кто-нибудь будет разбираться в твоём говнокоде? Ха-ха-ха.
Своих, хватает.
Проблемный участок и список ошибок. Больше ничего и не нужно. Если понадобится, скажем.
А ещё лучше, приведи уравнение, по которому считал. Нам чё, его из твоего говнокода, в котором ошибок, как грязи, выводить? Ты, в своём уме?

[type_Oleg]

anpolol, по поводу направления вектора начальной скорости планеты.
У вас он как раз направлен прочь от барицентра двух звезд. Потому что вектор положения r (500;500;500) и вектор скорости v(1;1;1) - одинаково направлены.
См. рисунок.
В таком случае движение будет очень простое - или планета просто улетает по прямой почти прочь ( если скорость больше параболической) , или отлетает по прямой, потом летит обратно и врезается в звезду ( если одна звезда. Если двойная звезда, то сложнее ..)
Условие перпендикулярности - из определения векторного произведения:
r*v = 0
Например, вектор v = (1;-0,5;-0,5) будет перпендикулярным вашему r, потому что 500*1+500*(-0.5)+500*(-0.5) = 0

И еще - откуда у вас исходные данные, из условий, или сами придумали?
Посмотрите вот тут - https://ru.wikipedia.org/wiki/Орбитальная_скорость или тут https://ru.wikipedia.org/wiki/Первая...еская_скорость
Там если не понятно μ = G*M

У вас для планеты, при ваших условиях - начальная скорость получается больше параболической, она улетит по гиперболе при любом направлении начального вектора скорости.
v параб. пл ~ 0.53 ; v нач.пл. = (1^2+1^2+1^2)^0.5 ~ 1.73

А у звезд - наоборот маловата скорость. Они будут вращаться по очень вытянутому эллипсу, и очень сильно сближаться. При сильном сближении резко возрастает погрешность расчета, особенно у такого простого метода. И скорее всего по расчету они тоже разлетятся.

И потом, получается что у всей системы начальный импульс не равен 0. Поэтому барицентр всей системы постоянно куда-то движется в этой абсолютной системе координат. Это неудобно для для отображения.