Metaballs and Gravity game; For past two days I made my first OpenGL project, I know, It's fixed pipeline, but my focus was on gravity algorithm, so, tell what do you think, and how can I make it faster?

main.cpp

#pragma comment(lib, "glew32.lib")
#include &lt;stdlib.h&gt;
#include &lt;string&gt;
#include &lt;vector&gt;
#include &lt;iostream&gt;
#include &lt;fstream&gt;
#include "GL/glew.h"
#define GLEW_STATIC
#include &lt;glut.h&gt;
#define WIDTH   800
#define HEIGHT  640
#define AU      (149.6e6 * 1000)
#define SCALE   (250 / AU)
#define G       6.67428e-11

int scale = 10;

struct Vec2 {
    double x;
    double y;
};

using namespace std;

struct Object {
    int id;
    Vec2 pos;
    Vec2 vel;
    double mass;
    double r;
};

Object b = Object{ 0, Vec2{ 2 * AU, 1 * AU, }, Vec2{ 0, 20 * 1000 }, 1.98892 * pow(10, 30), 40 };
Object b2 = Object{ 1, Vec2{ 1 * AU, 1 * AU }, Vec2{ 0, -20 * 1000 }, 1.98892 * pow(10, 30), 40 };
vector&lt;Object&gt; bodies;
int bodyIndex = 0;


void init() {
    glClearColor(0, 0, 0, 1);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0, WIDTH, 0, HEIGHT);
    glewInit();

    glPointSize(10);

    bodies.push_back(b);
    bodies.push_back(b2);
    bodyIndex = bodies.size();
}

Vec2 attraction(Object body, Object other) {
    // Compute the distance of the other body
    double dx = other.pos.x - body.pos.x;
    double dy = other.pos.y - body.pos.y;
    double d = sqrt(pow(other.pos.x - body.pos.x, 2) + pow(other.pos.y - body.pos.y, 2));

    // Report an error if the distance is zero; otherwise we'll
    // get a ZeroDivisionError exception further down
    if (d == 0) std::cout &lt;&lt; "[ ! ] Collision between objects is not possible;" &lt;&lt; std::endl;

    // Compute the force of attraction
    double f = G * body.mass * other.mass / pow(d, 2);

    //Compute the direction of the force
    double theta = atan2(dy, dx);
    return Vec2{ cos(theta) * f, sin(theta) * f };
}

void render() {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);

    int timestep = 24 * 3600;
    glBegin(GL_POINTS);
    for (auto b : bodies) {
        Vec2 totalForce{ 0, 0 };
        for (auto other : bodies) {
            if (bodies.at(b.id).id != bodies.at(other.id).id) {
                Vec2 force = attraction(bodies.at(b.id), bodies.at(other.id));
                totalForce.x += force.x;
                totalForce.y += force.y;
            }
        }

        bodies.at(b.id).vel.x += totalForce.x / bodies.at(b.id).mass * timestep;
        bodies.at(b.id).vel.y += totalForce.y / bodies.at(b.id).mass * timestep;
        bodies.at(b.id).pos.x += bodies.at(b.id).vel.x * timestep;
        bodies.at(b.id).pos.y += bodies.at(b.id).vel.y * timestep;

        for (int x = 0; x &lt; WIDTH; x += scale) {
            for (int y = 0; y &lt; HEIGHT; y += scale) {
                float sum = 0;
                for (auto b : bodies) {
                    sum += (b.r / (sqrt(pow(b.pos.x*SCALE - x, 2) + pow(b.pos.y*SCALE - y, 2))*0.001))*0.0001;
                    glColor3f(sum, sum, sum);
                    glVertex2i(x, y);
                }
            }
        }
    }
    glEnd();

    glutSwapBuffers();
}

void OnMouseClick(int button, int state, int x, int y){
    if (button == GLUT_LEFT_BUTTON &amp;&amp; state == GLUT_DOWN){
        Object X = Object{ bodyIndex, Vec2{ x / SCALE, (HEIGHT - y) / SCALE, }, Vec2{ 0, 0 }, 3.302 * pow(10, 13), 40 };
        bodies.push_back(X);
        bodyIndex++;
    }
}

void timer(int) {
    glutPostRedisplay();
    glutTimerFunc(1000 / 60, timer, 0);
}

int main(int argc, char* argv[]) {
    glutInit(&amp;argc, argv);
    glutInitWindowSize(WIDTH-2, HEIGHT);
    glutInitWindowPosition(500, 200);
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
    glutCreateWindow("Gravity");

    init();
    glutDisplayFunc(render);
    glutMouseFunc(OnMouseClick);
    glutTimerFunc(50, timer, 50);

    glutMainLoop();
    return 0;
}
Refactoring code structure

As of right now, you have:


A struct named Object.
A global std::vector&lt;Object&gt; of Objects.
A function named attraction, which (I assume) calculates the attraction between two separate Objects.
A function named render which (I assume) renders all of the objects in bodies.


While this structure works fine, it's somewhat rigid and will probably cause you issues in the future if you ever try to expand this project by any significant amount. Instead of having this structure, I'd recommend that you design something more akin to the following:


A struct (or class) named Object which represents a gravitational body, with the following member functions:


CalculateAttraction(), for calculating object attractions.
Update(), for updating anything related to the object, like position, etc.
Render(), for rendering the object.

A struct (or class) named ObjectManager which is used to store and manage objects, with the following member functions:


AddBody(), for adding a new gravitational body.
Update(), for updating the object manager and its stored Objects.
Render(), for rendering the object manager and its stored Objects.

A top-level Update() function for GLUT to use.
A top-level Render() function for GLUT to use.


While you don't have to follow this exact structure (you may also need additional behaviors included with the ones I provided), it will make managing and expanding this project much easier.

Using const instead of #define

In C++ (and even later versions of C), it's always better to use const (or constexpr in some special cases) for defining constant values rather than a macro. Using macros can often lead to obscure and confusing issues popping up, especially if you're overzealous in your use of them. The macros in your code (excluding GLEW_STATIC) should be declared as:

const int WIDTH = 800;
const int HEIGHT = 800;
const float AU = 149.6e6 * 1000;
const float SCALE = 250 / AU;
const float G = 6.67428e-11;


Instead of:


#define WIDTH   800
#define HEIGHT  640
#define AU      (149.6e6 * 1000)
#define SCALE   (250 / AU)
#define G       6.67428e-11



Nitpicks


You have a small bug in your code where you #define GLEW_STATIC. In order to use GLEW properly as a static library, GLEW_STATIC must be declared before you include glew.h. You'll also likely have to change "glew32.lib" to "glew32s.lib". Example:

#pragma comment(lib, "glew32s.lib")

...

#define GLEW_STATIC
#include "glew.h"

Please don't ever use using namespace std; in your C++ code. It's a very bad idea. See this Stack Overflow post for the reason why.
The variable bodyIndex is poorly named and unnecessary. std::vector&lt;T&gt; already has a size() function which returns the number of elements in the vector.
While you mentioned that you know you used the fixed function pipeline, I'm going to recommend against it anyways. It's deprecated in a lot of graphics hardware and sometimes requires some special workarounds to get working. It's also much less flexible than the programmable function pipeline.
The various #includes and #defines at the beginning of your code are are grouped together in one big cluster, which isn't very visually appealing. I'd highly recommend separating them into distinct groups, like so:

#pragma comment(lib, "glew32.lib")
...

#include &lt;stdlib.h&gt;
#include &lt;string&gt;
#include &lt;vector&gt;
#include &lt;iostream&gt;
#include &lt;fstream&gt;
...

#define GLEW_STATIC
#include "GL/glew.h"
#include &lt;glut.h&gt;
...

#define WIDTH   800
#define HEIGHT  640
#define AU      (149.6e6 * 1000)
#define SCALE   (250 / AU)
#define G       6.67428e-11
...


As you can see, there are four distinct groups here, #pragmas, standard #includes, 3rd-party #includes and #defines.