Во-первых, я использую что-то очень близкое к методу Эйлера, чтобы вычислить положение моей планеты. Я знаю, что это не самый точный метод, но я уже почти неделю работаю над использованием скоростного верлета и не могу заставить его работать. Моя проблема в том, что я не могу заставить свою планету вращаться вокруг Солнца, позиции х или у постоянно увеличиваются. Любая помощь приветствуется! Спасибо!
Вот мой код:
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;
}
Поскольку вы используете декартову систему координат, вы должны обнаружить, что ваша сила вдоль x в некоторых конфигурациях положительна, а в других — отрицательна и одинакова вдоль y. В частности, сила F — это вектор, который указывает от центра масс планеты к центру масс Солнца, вокруг которого он вращается. Итак, есть две проблемы с кодом:
Вы возвращаете 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, следовательно, сила указывает на отрицательное значение Икс).
Ускорение a равно F / m, тогда как у вас есть vec F * (delta x) / m, где vec F — сила вектор и дельта х — это расстояние х между планетой и солнцем. А находится в том же направлении, что и вектор силы. Таким образом, у вас должно быть что-то вроде ax = Fx / m1 и ay = Fy / m1, если вы хотите вычислить ускорение массы m1 (скажем, планеты).
Исправление кода на основе вышеизложенного должно быть простым, поэтому я оставлю это на ваше усмотрение, но если не ясно, просто оставьте комментарий.
Чтобы избежать использования угла и преобразования его вперед и назад с помощью тригонометрических функций, используйте векторную силу
F = (Fx, Fy) = — (x, y) * G * m1 * m2 / r ^ 3.
Кроме того, используйте учебник «Движущиеся звезды вокруг» для справочной реализации различных симплектических и не очень симплектических методов интеграции и Хайрер-идр. для теории и истории позади этого.