c++模拟质点在万有引力下的运动轨迹

$$\mathcal{P}hysics\mathcal{I}\mathcal{S}the\mathcal{M}ost\mathcal{B}eautiful$$

通过c++模拟质点在万有引力下的运动轨迹

0.$\mathcal{F}orward\mathcal{I}troduction$

应某集团首领要求：

通过代码实现天体间万有引力作用的模拟运动轨迹
时限一晚上

通过模拟键盘和鼠标的点击在画图中绘制图像

1.$\mathcal{P}re\mathcal{K}nowledge1\to \mathcal{T}o\mathcal{C}ontrol\mathcal{K}eyboard$

直接上代码

keybd_event(num/*Key_Code*/, 0, 0, 0);//按下按键
keybd_event(num/*Key_Code*/, 0, KEYEVENTF_KEYUP, 0);//放开按键

num是你所要按下的键的键码，字母的键码就是本身的Ascll码，键码表如下（来自网络）：
cnblogs-万一

2.$\mathcal{P}re\mathcal{K}nowledge2\to \mathcal{T}o\mathcal{C}ontrol\mathcal{M}ouse$

仍然上代码

POINT p;
GetCursorPos(&p);//获取鼠标坐标
printf("%d %d\n", p.x, p.y);

SetCursorPos(x, y);//移动到(x,y)
mouse_event(MOUSEEVENTF_LEFTDOWN, x, y, 0, 0);//按下左键
mouse_event(MOUSEEVENTF_LEFTUP, x, y, 0, 0);//放开左键
mouse_event(MOUSEEVENTF_RIGHTDOWN, x, y, 0, 0);//按下右键
mouse_event(MOUSEEVENTF_RIGHTUP, x, y, 0, 0);//放开右键

鼠标的坐标$(p.x,p.y)$以像素为单位，且一定为整数

3.$\mathcal{P}re\mathcal{K}nowledge3\to \mathcal{U}niversal\mathcal{G}ravitatio{n}^{\prime }s\mathcal{L}aw$

根据中学知识知道
$$F_g=\frac{G{m}_1{m}_2}{R^2}$$
$$G\approx 6.67\times 10^{-11}N\cdot m^2\cdot k{g}^{-2}$$
因为$G$的精度要求较大，所以可以适当缩放单位大小，如：

$m\to km$,$kg\to 10000t$

4.$\mathcal{P}re\mathcal{K}nowledge4\to \mathcal{N}ewto{n}^{\prime }s\mathcal{S}econd\mathcal{L}aw$

$F=ma$

5.$\mathcal{S}o\mathcal{D}o\mathcal{I}t\to \mathcal{P}hysics$

我们可以把时间划分成很小很小的单位，如$10ms$，每次以上一次的坐标近似表示这$10ms$内位置
$F=\frac{G{m}_1{m}_2}{R^2}$用$\mathcal{O}(n^2)$的复杂度两两之间的引力
再算出每个质点的向量和$a=\frac{F}{m},\Delta v=at,\Delta x=vt$($t$取$0.01s$)

6.$\mathcal{T}urn\mathcal{I}nto\mathcal{C}ode$

此代码仅兼容部分设备qwq

#include <cstdio>
#include <iostream>
#include <cmath>
#include <windows.h>
#include <winuser.h>
#define N 100
#define KEY_DOWN(VK_NONAME) ((GetAsyncKeyState(VK_NONAME) & 0x8000) ? 1:0)

char st[100];
const int ox = 795, oy = 492, unt = 10;
const long double G = 0.667408;

struct Node{
	long double x, y, vx, vy, fx, fy, m;
}a[N];

inline void cli_key(int id) {
	keybd_event(id, 0, 0, 0);
	keybd_event(id, 0, KEYEVENTF_KEYUP, 0);
}

inline void cli_mouse(int x, int y) {
	SetCursorPos(x, y);
	mouse_event(MOUSEEVENTF_LEFTDOWN, x, y, 0, 0);
	mouse_event(MOUSEEVENTF_LEFTUP, x, y, 0, 0);
}

inline void prt_st(char *st) {
	for (int i = 0; i < strlen(st); ++i) {
		int id;
		if (st[i] >= 'a' && st[i] <= 'z') id = st[i] - 'a' + 'A';
		else id = st[i];
		cli_key(id);
	}
}

inline void prt_line(int x1, int y1, int x2, int y2) {
	SetCursorPos(x1, y1);
	mouse_event(MOUSEEVENTF_LEFTDOWN, x1, y1, 0, 0);
	SetCursorPos(x2, y2);
	mouse_event(MOUSEEVENTF_LEFTUP, x2, y2, 0, 0);
}

inline void cli_pencil() {
	cli_mouse(269, 78);
}

inline void getxy() {
	Sleep(5000);
	POINT p;
	GetCursorPos(&p);
	printf("%d %d\n", p.x, p.y);
}

int main() {
	int n, time;
	puts("time(s)?");
	scanf("%d", &time);
	puts("The scale is 100metres");
	printf("The number of particle(n <= 20):");
	scanf("%d", &n);
	for (int i = 1; i <= n; ++i) {
		printf("The X & Y(*100m) of %d-th particle:", i);
		std :: cin >> a[i].x >> a[i].y;
		printf("The Vx & Vy(m/s) of %d-th particle:", i);
		std :: cin >> a[i].vx >> a[i].vy;
		printf("The M(*10000000kg) of %d-th particle:", i);
		std :: cin >> a[i].m;

	}
	keybd_event(91, 0, 0, 0);
	keybd_event('R', 0, 0, 0);
	keybd_event('R', 0, KEYEVENTF_KEYUP, 0);
	keybd_event(91, 0, KEYEVENTF_KEYUP, 0);
	Sleep(1000);
	prt_st("mspaint");
	cli_key(13);
	Sleep(1000);
	keybd_event(18, 0, 0, 0);
	keybd_event(32, 0, 0, 0);
	keybd_event('X', 0, 0, 0);
	keybd_event('X', 0, KEYEVENTF_KEYUP, 0);
	keybd_event(32, 0, KEYEVENTF_KEYUP, 0);
	keybd_event(18, 0, KEYEVENTF_KEYUP, 0);
	Sleep(1000);
	keybd_event(17, 0, 0, 0);
	keybd_event('E', 0, 0, 0);
	keybd_event('E', 0, KEYEVENTF_KEYUP, 0);
	keybd_event(17, 0, KEYEVENTF_KEYUP, 0);
	Sleep(1000);
	prt_st("1580");
	cli_key(9);
	prt_st("660");
	cli_key(13);
	cli_pencil();
//	getxy();
	prt_line(5, 492, 1585, 492);
	prt_line(795, 162, 795, 822);
	prt_line(ox + unt, oy, ox + unt, oy - 7);
	for (int k = 1; k <= time / 0.01; ++k) {
		for (int i = 1; i <= n; ++i) {
			a[i].fx = a[i].fy = 0;
			for (int j = 1; j <= n; ++j) 
				if (i != j) {
					long double sqrR = (a[i].x - a[j].x) * (a[i].x - a[j].x) + (a[i].y - a[j].y) * (a[i].y - a[j].y);
					long double F = G * (a[i].m * a[j].m) / sqrR;
					if (F > 10000000) F = 10000000;
					a[i].fx += F * (a[j].x - a[i].x) / sqrt(sqrR);
					a[i].fy += F * (a[j].y - a[i].y) / sqrt(sqrR);
				}
		}
		for (int i = 1; i <= n; ++i) {
			printf("%.4Lf %.4Lf %.4Lf %.4Lf %.4Lf %.4Lf %.4Lf\n", a[i].x, a[i].y, a[i].vx, a[i].vy, a[i].m, a[i].fx, a[i].fy);
			a[i].vx += a[i].fx / a[i].m * 0.01;
			a[i].vy += a[i].fy / a[i].m * 0.01;
			long double xx = a[i].x + a[i].vx * 0.01;
			long double yy = a[i].y + a[i].vy * 0.01;
			int xx1 = (int)(a[i].x * unt) + ox, yy1 = (int)(a[i].y * unt) + oy, xx2 = (int)(xx * unt) + ox, yy2 = (int)(yy * unt) + oy;
			int flg = 0;
			if (xx2 <= 10) {
				a[i].x = 78;
				flg = 1;
			}
			if (xx2 >= 1580) {
				a[i].x = -78;
				flg = 1;
			}
			if (yy2 <= 167) {
				a[i].y = 32;
				flg = 1;
			}
			if (yy2 >= 817) {
				a[i].y = -32;
				flg = 1;
			}
			if (flg) continue;
			prt_line(xx1, yy1, xx2, yy2);
			a[i].x = xx;
			a[i].y = yy;
		}
		for (int i = 1; i <= 1 * 100000; ++i) {
			if (KEY_DOWN(27)) {
				Sleep(300);
				while (!KEY_DOWN(27));
				Sleep(300);
				break;
			}
		}
//		Sleep(10);
	}
	return 0;
}

7.$\mathcal{P}ictures$

data1

30
3
0 0
0 0
1000000
20 0
0 30
20000
-20 0
0 -30
20000

data2

240
2
7 0
0 10
10000
-7 0
0 -10
10000

咕咕咕。。