Проект моделирования солнечной системы (справка по скорости)

Есть некоторые проблемы с физикой и некоторые проблемы с кодом.

Во-первых, проблема с физикой. Предполагая, что мы не моделируем альтернативную вселенную, где законы физики различны, закон всемирного тяготения Ньютона говорит, что F = G * m1 * m2 / (r * r). Однако сила — это вектор, а не скаляр, поэтому она имеет величину и направление.

Из чего рассчитывается код forceFuncX на самом деле это величина, а не просто составляющая силы, параллельной оси X. forceFuncY имеет тот же недостаток.

Далее идет расчет на ускорение. Физика говорит, что это = F / m. Масса — это скаляр, но и ускорение, и сила — это векторы. Таким образом, чтобы вычислить a_x, мы могли бы использовать F_x / m или F * cos () / М. Потому что () (с быть углом от одного CelesitalObject для другого в 2D-пространстве) = dx / r, мы можем сделать это a_x = F * dx / (m * r), что почти, но не совсем то, что вы получили в своих вычислениях (в делителе отсутствует r).

Альтернативный подход заключается в использовании std::complex но я не буду показывать этот подход с предположением, что вы можете расширить эту модель до трех измерений.

Это приводит нас к проблемам с кодом. Во-первых, поскольку вы используете C ++ и пишете симуляцию физической системы дискретных объектов, имеет смысл, что вы определили CelestialObject учебный класс. Что менее важно, так это то, что ваши функции вызываются путем выделения отдельных частей этих объектов и последующего вызова функций в стиле C. Мы можем улучшить код, используя эти объекты лучше. Для начала, так как вы еще не опубликовали один, вот CelestialObject класс, основанный на интерфейсе, который я вывел из вашего кода:

class CelestialObject
{
public:
CelestialObject(std::string name, float mass, float X, float Y,
float VX, float VY)
: myname(name), m(mass), x(X), y(Y), vx(VX), vy(VY) {}
void setPosition(float X, float Y) { x=X; y=Y; }
void setVelocity(float VX, float VY) { vx=VX; vy=VY; }
float getMass() const { return m; }
float getX() const { return x; }
float getY() const { return y; }
float getVX() const { return vx; }
float getVY() const { return vy; }
friend std::ostream& operator<<(std::ostream& out,
const CelestialObject& obj) {
return out << obj.myname << '\t' << obj.x << '\t' << obj.y
<< '\t' << obj.vx << '\t' << obj.vy << std::endl;
}
private:
std::string myname;
float m, x, y;
float vx, vy;
};

Далее несколько вспомогательных функций:

// returns square of distance between objects
float distance_sq(const CelestialObject &a, const CelestialObject &b)
{
// distance squared is (dy^2 + dx^2)
return pow(a.getY()-b.getY(),2) + pow(a.getX()-b.getX(),2);
}

// returns magnitude of the force between the objects
float force(const CelestialObject &a, const CelestialObject &b)
{
//  F=(G * m1 * m1)/(r^2) in the direction a->b and b->a
return G*a.getMass()*b.getMass()/distance_sq(a, b);
}

// returns the angle from a to b
float angle(const CelestialObject &a, const CelestialObject &b)
{
return atan2f(b.getY()-a.getY(),b.getX()-a.getX());
}

Наконец актуальный верлет:

void updatePosition(CelestialObject &a, CelestialObject &b )
{
float F = force(a,b);
float theta = angle(a,b);
float accela = F/a.getMass();
float accelb = -F/b.getMass();

// now that we have the acceleration of both objects, update positions
// x = x +v *dt + a*dt*dt/2
//   = x + dt * (v + a*dt/2)
a.setPosition(
a.getX() + dt * (a.getVX() + accela*cos(theta)*dt/2),
a.getY() + dt * (a.getVY() + accela*sin(theta)*dt/2)
);
b.setPosition(
b.getX() + dt * (b.getVX() + accelb*cos(theta)*dt/2),
b.getY() + dt * (b.getVY() + accelb*sin(theta)*dt/2)
);
// get new acceleration a'
F = force(a,b);
float thetap = angle(a,b);
float accelap = F/a.getMass();
float accelbp = -F/b.getMass();
// and update velocities
// v = v + (a + a')*dt/2
a.setVelocity(
a.getVX() + (accela*cos(theta) + accelap*cos(thetap))*dt/2,
a.getVY() + (accela*sin(theta) + accelap*sin(thetap))*dt/2
);
b.setVelocity(
b.getVX() + (accelb*cos(theta) + accelbp*cos(thetap))*dt/2,
b.getVY() + (accelb*sin(theta) + accelbp*sin(thetap))*dt/2
);
}

И, наконец, простой тестовый код.

#include <string>
#include <iostream>
#include <vector>
#include <cmath>

const float G(6.67e-11);  // N*(m/kg)^2
const float dt(0.1);      // s
// all of the other code goes here...
int main()
{
CelestialObject anvil("anvil", 70, 370, 0, 0, 0);
CelestialObject earth("earth", 5.97e+24, -6.378e6, 0, 0, 0);
std::cout << "Initial values:\n" << earth << anvil;
std::cout << "Dropping an anvil from the top of a 370m building...\n""It should hit the ground in about 8.7 seconds.\n";
int t;
for (t=0; anvil.getX() > 0; ++t) {
std::cout << dt*t << '\t' << anvil;
updatePosition(anvil, earth);
}
std::cout << "Final values at t = " << dt*t << " seconds:\n"<< earth << anvil;
return 0;
}

В тестовом коде используется временной шаг 0,1 с, который слишком мал для вашей солнечной системы, но подходит для этого быстрого теста, который заключается в том, чтобы увидеть, получим ли мы приемлемый результат для известной системы. В данном случае я выбрал систему из двух тел, состоящую из планеты Земля и наковальни. Этот код имитирует падение наковальни с вершины 370-метрового здания, которое, как мы можем легко рассчитать, ударит о землю примерно через 8,7 с, если пренебречь сопротивлением воздуха. Чтобы не усложнять координаты, я решил поместить начало координат (0,0) на поверхность земли и рассмотреть верхнюю часть здания в точке (370,0). Когда код скомпилирован и запущен, он производит следующее:

Initial values:
earth   -6.378e+06  0   0   0
anvil   370 0   0   0
Dropping an anvil from the top of a 370m building...
It should hit the ground in about 8.7 seconds.
0   anvil   370 0   0   0
0.1 anvil   369.951 -4.27834e-09    -0.97877    -8.55668e-08
0.2 anvil   369.804 -1.71134e-08    -1.95754    -1.71134e-07
0.3 anvil   369.56  -3.85051e-08    -2.93631    -2.567e-07
...
8.3 anvil   32.8567 -2.9474e-05 -81.2408    -7.1023e-06
8.4 anvil   24.6837 -3.01885e-05    -82.2197    -7.18787e-06
8.5 anvil   16.4127 -3.09116e-05    -83.1985    -7.27345e-06
8.6 anvil   8.04394 -3.16432e-05    -84.1774    -7.35902e-06
Final values at t = 8.7 seconds:
earth   -6.378e+06  3.79705e-28 9.98483e-22 8.72901e-29
anvil   -0.422744   -3.23834e-05    -85.1563    -7.4446e-06

Как видите, это похоже на работу, но есть проблемы. Первая проблема состоит в том, что, поскольку объекты должны двигаться только вдоль оси X, все компоненты Y должны быть равны 0. Это не так, потому что этот код не очень хорошо спроектирован с точки зрения численного анализа. Выполнение сложений и вычитаний чисел с плавающей запятой, когда одно число большое, а другое маленькое, является одной из проблем. Другое использование atan2f функция, которая возвращает только float а затем с помощью этого результата в cos() а также sin(), Тригонометрические функции на самом деле лучше избегать, если это возможно.

Наконец, эта программа в настоящее время работает только с двумя объектами. Добавление третьего было бы болезненным с такой схемой, поэтому лучше было бы обработать std::vector<CelestialObject> сначала рассчитав равнодействующая сила на каждом объекте, учитывая позиции и массы всех остальных. Я оставлю это вам, но это должно по крайней мере дать вам старт в правильном направлении.