Finding the force of acceleration on an attractive body

Question

I have a simulation with a planet orbiting a star. The equation for force on the planet is:

G⋅M₁⋅M₂⋅(q₁−q₂)/||q₁−q₂||³

Where G is the gravitational constant, M₁ and M₂ are the masses of the star and planet, and q₁−q₂ is the difference between the position vectors of the two bodies.

I have converted this equation into component form where

F_x = G⋅M₁⋅M₂⋅(x₁−x₂)/||q₁−q₂||³

F_y = G⋅M₁⋅M₂⋅(y₁−y₂)/||q₁−q₂||³

with x₁−x₂ being the difference between the x component of the vectors, and y₁−y₂ being the difference between the y components of the vectors.

The ||q₁−q₂||³ can be rewritten as sqrt((x₁−x₂)²+(y₁−y₂)²)³ which can we rewritten as ((x₁−x₂)²+(y₁−y₂)²)^3/2

This is my code

    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);

For some reason, the planet does not follow an orbit like it should. Does anyone know what I got wrong?

EDIT:: Forgot to mention, all of these variables are doubles, and the language is c++

EDIT 2:: Here is the complete code

    #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;
    }

With Visual C++ 2015 it can be compiled as follows (with console subsystem, suitable for testing a malfunctioning application):

cl orbit.cpp -D UNICODE gdi32.lib user32.lib kernel32.lib /link /subsystem:console /entry:WinMainCRTStartup

Answer 1

I had a problem with the formula I was using, I believe the correct way to separate the vectors is using the equation C*(q1-q2)/||q1-q2||^2

This is the code:

    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);

This seems to work for now, there is some instability which I think I'll be able to deal with. Thanks a lot!

Answer 2

Not an answer, but this won't fit in a comment.

You can use <chrono> in a much safer, more accurate, more efficient, and more readable manner if you make the following simple adjustments:

auto newtime = std::chrono::system_clock::now();
auto oldtime = newtime;
// ...
using frames = std::chrono::duration<
                                std::chrono::system_clock::rep, std::ratio<1, 60>>;
while (newtime < oldtime + frames{1}) {//60 frames per second
    newtime = std::chrono::system_clock::now();
}
oldtime = newtime;