我有一个绕着恒星运行的行星的模拟。行星上的力方程为:
G⋅M1⋅M2⋅(q1-q2)/ || q1-q2 || 3
其中G是重力常数，M1和M2是恒星和行星的质量，而q1-q2是两个物体的位置 vector 之间的差。
我已经将这个方程式转换成分量形式，其中
Fx =G⋅M1⋅M2⋅(x1-x2)/ || q1-q2 || 3
Fy =G⋅M1⋅M2⋅(y1-y2)/ || q1-q2 || 3
其中x1-x2是 vector 的x分量之差，y1-y2是 vector 的y分量之差。
|| q1-q2 || 3可以重写为sqrt((x1-x2)2+(y1-y2)2)3，我们可以重写为((x1-x2)2+(y1-y2)2) 3/2
这是我的代码
xdif = (planets[0].x - stars[0].x);
ydif = (planets[0].y - stars[0].y);
planets[0].ax = -10*xdif / pow((pow(xdif, 2) + pow(ydif, 2)), 1.5);
planets[0].ay = -10*ydif / pow((pow(xdif, 2) + pow(ydif, 2)), 1.5);
由于某种原因，行星没有按照应有的轨道运动。有人知道我错了吗？
编辑::
忘了说了，所有这些变量都是 double 的，语言是c++
编辑2::这是完整的代码
#define WIN32_LEAN_AND_MEAN
#include <math.h>
#include <windows.h>
#include <stdlib.h>
#include <string>
#include <iostream>
#include <fstream>
#include <chrono>
using namespace std;
#define SCREEN_WIDTH 1920
#define SCREEN_HEIGHT 1080
#define LGREY RGB(200,200,200)
#define BLACK RGB(0, 0, 0)
#define WHITE RGB(255,255,255)
WPARAM w;
HINSTANCE hInst;
LRESULT CALLBACK WindowProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam);
POINT mousePos;
int t;
struct star { double x, y; };
struct Body { double x, y, vx, vy, ax, ay; };
star stars[2];
Body planets[1];
long newtime= std::chrono::system_clock::now().time_since_epoch().count();
long oldtime=newtime;
int signx;
int signy;
double force = 0;
double xdif;
double ydif;
void drawCircle(int x, int y, int w, int h, COLORREF insidecolor, COLORREF bordercolor, HDC hdc) {
HGDIOBJ B1 = CreateSolidBrush(insidecolor);
HPEN P1 = CreatePen(PS_SOLID, 1, bordercolor);
SelectObject(hdc, B1);
SelectObject(hdc, P1);
Arc(hdc, x, y, x + w, y + h, x - 1, y - 1, x + 1, y + 1);
DeleteObject(B1);
DeleteObject(P1);
}
int WINAPI WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow) {
mousePos = { 1200, 600 };
WNDCLASSEX wc;
ZeroMemory(&wc, sizeof(WNDCLASSEX));
wc.cbSize = sizeof(WNDCLASSEX);
wc.style = CS_HREDRAW | CS_VREDRAW;
wc.lpfnWndProc = WindowProc;
wc.hInstance = hInstance;
wc.hCursor = LoadCursor(NULL, IDC_ARROW);
wc.lpszClassName = L"WindowClass";
wc.hbrBackground = CreateSolidBrush(WHITE);
RegisterClassEx(&wc);
HWND hWnd = CreateWindowEx(NULL, L"WindowClass", L"Template", WS_OVERLAPPEDWINDOW, 10, 10, SCREEN_WIDTH, SCREEN_HEIGHT, NULL, NULL, hInstance, NULL);
hInst = hInstance;
//separated by 10 au
stars[0].x = 0;
stars[0].y = 0;
stars[1].x = 10;
stars[1].y = 0;
planets[0].x = -5;
planets[0].y = 0;
planets[0].vx = 0;
planets[0].vy = -1;
planets[0].ax = 0;
planets[0].ay = 0;
ShowWindow(hWnd, nCmdShow);
UpdateWindow(hWnd);
MSG msg;
bool running = TRUE;
while (running) {
if (PeekMessage(&msg, NULL, 0, 0, PM_REMOVE)) {
TranslateMessage(&msg);
DispatchMessage(&msg);
if (msg.message == WM_QUIT) {
running = FALSE;
}
}
}
return msg.wParam;
}
LRESULT CALLBACK WindowProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam) {
PAINTSTRUCT ps;
HDC hdc;
switch (message)
{
case WM_PAINT:
hdc = BeginPaint(hWnd, &ps);
planets[0].x += planets[0].vx;
planets[0].y += planets[0].vy;
planets[0].vx += planets[0].ax;
planets[0].vy += planets[0].ay;
xdif = (planets[0].x - stars[0].x);
ydif = (planets[0].y - stars[0].y);
planets[0].ax = -10.0*xdif / pow((pow(xdif, 2) + pow(ydif, 2)), 1.5);
planets[0].ay = -10.0*ydif / pow((pow(xdif, 2) + pow(ydif, 2)), 1.5);
for (int i = 0; i < 1; i++) {
drawCircle(planets[i].x * 50 + 1920 / 2 - 10, planets[i].y + 1080 / 2 - 10, 20, 20, BLACK, BLACK, hdc);
}
for (int i = 0; i < 2; i++) {//1920x1080
drawCircle(stars[i].x*50+1920/2 - 25, stars[i].y+1080/2 - 25, 50, 50, BLACK, BLACK, hdc);
}
while (newtime < oldtime + 10000000/60) {//60 frames per second
newtime = std::chrono::system_clock::now().time_since_epoch().count();
}
oldtime = newtime;
InvalidateRect(hWnd, NULL, TRUE);
EndPaint(hWnd, &ps);
break;
case WM_DESTROY:
PostQuitMessage(0);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
return 0;
}
使用Visual C++ 2015，可以按以下方式对其进行编译(使用控制台子系统，适用于测试故障应用程序):
cl orbit.cpp -D UNICODE gdi32.lib user32.lib kernel32.lib / link / subsystem:控制台/ entry:WinMainCRTStartup
最佳答案
我使用的公式有问题，我相信分离 vector 的正确方法是使用公式C *(q1-q2)/ || q1-q2 || ^ 2
这是代码:
xdif = (planets[0].x - stars[0].x);
ydif = (planets[0].y - stars[0].y);
planets[0].ax = -.01*xdif / pow(hypot(xdif, ydif), 2);
planets[0].ay = -.01*ydif / pow(hypot(xdif, ydif), 2);
这似乎暂时可行，我认为我可以解决一些不稳定问题。非常感谢!
关于c++ - 在有吸引力的 body 上寻找加速力，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38134783/