модель — Solar System Project, C ++. Нельзя получить отрицательные ускорения (планета не будет вращаться вокруг)

Во-первых, я использую что-то очень близкое к методу Эйлера, чтобы вычислить положение моей планеты. Я знаю, что это не самый точный метод, но я уже почти неделю работаю над использованием скоростного верлета и не могу заставить его работать. Моя проблема в том, что я не могу заставить свою планету вращаться вокруг Солнца, позиции х или у постоянно увеличиваются. Любая помощь приветствуется! Спасибо!
Вот мой код:

void updatePosition(CelestialObject object1, CelestialObject object2 )

{ // -----------------------------X calculations-----------------------------------
//calc force
float forceX = forceFuncX(object1.getX(), object2.getX(),object1.getY(), object2.getY(), object1.getMass(), object2.getMass());`

//accel calc
float AX = accelerationFuncX(forceX,object1.getX(), object2.getX(),object1.getMass());
agk::PrintC("Accel X: ");
float AXprint = AX*dt;
agk::Print(AXprint);

//velocity
float VX  = object1.getVX();
VX = VX + AX*dt;
agk::PrintC("Velocity X: ");
agk::Print(VX);

//positionCalc
float X = object1.getX();
X = X + VX*dt;
agk::PrintC("Position X: ");
agk::Print(X);

//-------------------------Y calculations------------------------------------
//force
float forceY = forceFuncY(object1.getX(), object2.getX(),object1.getY(), object2.getY(), object1.getMass(), object2.getMass());

//accel
float AY = accelerationFuncY(forceY,object1.getY(),object2.getY(), object1.getMass()); //y
agk::PrintC("Accel Y: ");
float AYprint = AY*dt;
agk::Print(AYprint);

//velocity
float VY = object1.getVY();
VY = VY + AY*dt;
agk::PrintC("Velocity Y: ");
agk::Print(VY);

//position
float Y = object1.getY();
Y = Y + VY*dt;
agk::PrintC("Position Y: ");
agk::Print(Y);

object1.setPosition(X, Y);
agk::CreateParticles(X,Y);
}

Вот функции, которые он вызывает:

double forceFuncX(float object1x,float object2x,float object1y, float object2y, double mass1, double mass2)
{
float d = object1x - object2x;
float r = sqrt(pow(object2x - object1x,2) + pow(object2y-object1y,2));;
//float r = sqrt(pow(object1x-object2x,2)+pow(object1y-object2y,2));
//double F = (G*(mass1*mass2))/pow(d,2);
float F = (G*(mass1*mass2))/(r*r);

return F;
}double forceFuncY(float object1x,float object2x,float object1y, float object2y, double mass1, double mass2)
{
float d = object1y - object2y;
float r = sqrt(pow(object2x - object1x,2) + pow(object2y-object1y,2));;
//float r = sqrt(pow(object1x-object2x,2)+pow(object1y-object2y,2));
//double F = (G*(mass1*mass2))/pow(d,2);
float F = (G*(mass1*mass2))/(r*r);

return F;
}

а также

float accelerationFuncX(float force, float object1x, float object2x, double mass) //gives the acceleration of an object
{
float accel = (force*(object2x-object1x))/mass;
return accel;
}

float accelerationFuncY(float force, float object1y, float object2y, double mass)
{
float accel =(force*(object2y-object1y))/mass;
return accel;
}

0

Решение

Поскольку вы используете декартову систему координат, вы должны обнаружить, что ваша сила вдоль x в некоторых конфигурациях положительна, а в других — отрицательна и одинакова вдоль y. В частности, сила F — это вектор, который указывает от центра масс планеты к центру масс Солнца, вокруг которого он вращается. Итак, есть две проблемы с кодом:

  1. Вы возвращаете Fy и Fx, оба всегда положительны и равны гравитационной силе между двумя телами: F = G * m1 * m2 / r ^ 2. Но F — длина вектора силы. Тогда компоненты Fx и Fy вектора силы будут иметь вид Fx = — F cos alpha и Fy = — F sin alpha, где alpha — угол планеты относительно оси x. Например, если альфа = 0, а солнце находится в центре вашей системы координат (x = y = 0), то Fx = -F (поскольку планета находится в положительном x, а солнце в x = 0, следовательно, сила указывает на отрицательное значение Икс).

  2. Ускорение a равно F / m, тогда как у вас есть vec F * (delta x) / m, где vec F — сила вектор и дельта х — это расстояние х между планетой и солнцем. А находится в том же направлении, что и вектор силы. Таким образом, у вас должно быть что-то вроде ax = Fx / m1 и ay = Fy / m1, если вы хотите вычислить ускорение массы m1 (скажем, планеты).

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

1

Другие решения

Чтобы избежать использования угла и преобразования его вперед и назад с помощью тригонометрических функций, используйте векторную силу

F = (Fx, Fy) = — (x, y) * G * m1 * m2 / r ^ 3.

Кроме того, используйте учебник «Движущиеся звезды вокруг» для справочной реализации различных симплектических и не очень симплектических методов интеграции и Хайрер-идр. для теории и истории позади этого.

0

По вопросам рекламы [email protected]