C++模拟2D牛顿力学效果的示例代码

namespace NewtonianPhysics
{
	using Scalar = long double;//标量
	const Scalar G = 6.67259e-11l;//万有引力常数
	namespace _2D
	{
		class Vector//矢量
		{
		public:
			Scalar x;
			Scalar y;
			Vector() :x(0.l), y(0.l) {}
			Vector(Scalar inX, Scalar inY) :x(inX), y(inY) {}
			friend Vector operator+(const Vector& a, const Vector& b)
			{
				return Vector(a.x + b.x, a.y + b.y);
			}
			friend Vector operator-(const Vector& a, const Vector& b)
			{
				return Vector(a.x - b.x, a.y - b.y);
			}
			friend Vector operator*(const Vector& vec, Scalar n)
			{
				return Vector(vec.x * n, vec.y * n);
			}
			friend Vector operator*(Scalar n, const Vector& vec)
			{
				return Vector(vec.x * n, vec.y * n);
			}
			Vector& operator*=(Scalar n)
			{
				x *= n;
				y *= n;
				return *this;
			}
			Vector& operator+=(const Vector& right)
			{
				x += right.x;
				y += right.y;
				return *this;
			}
			Vector& operator-=(const Vector& right)
			{
				x -= right.x;
				y -= right.y;
				return *this;
			}
			Scalar LengthSq() const
			{
				return (x * x + y * y);
			}
			Scalar Length() const
			{
				return (std::sqrt(LengthSq()));
			}
			void Normalize()
			{
				float length = Length();
				x /= length;
				y /= length;
			}
			static Vector Normalize(const Vector& vec)
			{
				Vector temp = vec;
				temp.Normalize();
				return temp;
			}
		};
		class MassPoint
		{
		public:
			Scalar m;//质量/kg
			Vector location;//当前位置
			Vector v;//速度/(m/s)
			Vector f;//受到的合力/(N)
			MassPoint(Scalar mass, Vector _location = { 0,0 }, Vector velocity = { 0,0 })
				:m(mass), location(_location), v(velocity) {}
			void Update(Scalar deltaTime)
			{
				//由F=ma得出加速度并更新速度
				v += f * (1 / m) * deltaTime;
				//由速度得出位移并更新坐标
				location += v * deltaTime;
			}
		};
	}
}

using namespace NewtonianPhysics;
void Print(const _2D::MassPoint& p)
{
	setfillcolor(YELLOW);
	setlinecolor(RED);
	setlinestyle(PS_SOLID, 1);
	auto x = (p.location.x + 2e7) * 800 / 4e7, y = (2e7 - p.location.y) * 800 / 4e7;
	fillcircle(x, y, 5);
}
int main()
{
	initgraph(800, 800); 
	_2D::MassPoint a(5e30l, { 0,0 }), b(1e4, { 1e7,0 }, { 0,5e6 });
	while (1)
	{
		cleardevice();
		Print(a);
		Print(b);
		Sleep(100);
		//F=GMm/r²
		b.f = _2D::Vector::Normalize(a.location - b.location) * (G * a.m * b.m / (a.location - b.location).LengthSq());		a.Update(0.05);
		b.Update(0.1);
	}
	closegraph();
	return 0;
}