Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Gravity simulator in C++

  1. Aug 12, 2012 #1
    Hi everyone, I am new here, nice to meet you.

    For fun I am trying to make a gravity simulator in c++ and using SDL.

    I believe there is a problem in my merge function when I calculate the new velocity and trajectory of two colliding bodies.

    Also, does my main code look accurate. The part where I calculate the gravitational force, determine acceleration, and sum the acceleration from one body to all the other bodies.

    This is my first time attempting a project like this and my professors are gone for the summer so I am kind of in the dark.

    I am using euler integration for now, it is all I know and can implement.

    Thank you,
    Michaela

    if you want to run the code
    the collision code is here
    http://sdl-collide.sourceforge.net/

    Code (Text):

    /*
     * main.cpp
     *
     *  Created on: Jul 28, 2012
     *      Author: michaela
     *      Fnet = (G * m1 * m2) / r^2
     */

    #include "SDL_collide.h"
    #include <SDL/SDL.h>
    #include <SDL/SDL_gfxPrimitives.h>
    #include <list>
    #include <math.h>
    #include <iostream>

    const int SCREEN_HEIGHT = 600;
    const int SCREEN_WIDTH = 800;
    const int SCREEN_BPP = 32;

    const double G = 6.67e-11;
    const double DELTA_TIME = 0.1;

    //Hardcoded number of bodies, used if no argument is specified
    int bodies = 100;

    struct celestial_body
    {
            int count;
            double mass;
            double radius;
            double x, y;
            double vX, vY;
            double aX, aY;
            SDL_Color color;

            celestial_body(int newCount)
            {
                count = newCount;

                //Random mass
                mass = rand() % 100000000;

                radius = 0.5;

                //Random position on screen
                x = fmod((double)rand(), (double)SCREEN_WIDTH);
                y = fmod((double)rand(), (double)SCREEN_HEIGHT);

                vX = vY = 0;

                //Random color and avoid very dark colors
                color.r = (int)mass % (rand() % 10000) % 255;
                color.g = (int)mass % (rand() % 10000) % 255;
                color.b = (int)mass % (rand() % 10000) % 255;

                if(color.r <= 20 && color.g <= 20 && color.b <= 20)
                {
                    color.r += 25;
                    color.g += 25;
                    color.b += 25;
                }
            }
            void draw(SDL_Surface* screen)
            {
                filledCircleRGBA(screen, x, y, radius, color.r, color.g, color.b, 255);
            }
            void update()
            {
                vX = vX + (aX * DELTA_TIME);
                vY = vY + (aY * DELTA_TIME);

                x = x + (vX * DELTA_TIME);
                y = y + (vY * DELTA_TIME);
            }
            void print()
            {
                std::cout << "Mass: " << mass << "\nX: " << x << "\nY: " << y <<
                             "\nVx : " << vX << "   Vy: " << vY <<
                             "\nColor: " << (int)color.r << " " << (int)color.g << " " << (int)color.b << "\n\n";
            }
    };

    void merge(celestial_body &a, celestial_body &b)
    {
        //Radius and color merge is arbitrary (for now)
        a.radius += b.radius;
        a.color.r = (a.color.r + b.color.r) / 2;
        a.color.g = (a.color.g + b.color.g) / 2;
        a.color.b = (a.color.b + b.color.b) / 2;

        //Total velocity for a and b
        double vA = sqrt((a.vX * a.vX) + (a.vY * a.vY));
        double thetaA = atan2(a.vY, a.vX);
        if(thetaA < 0)
            vA *= -1.0;

        double vB = sqrt((b.vX * b.vX) + (b.vY * b.vY));
        double thetaB = atan2(b.vY, b.vX);
        if(thetaB < 0)
            vB *= -1.0;

        //New velocity of collided bodies
        double velocity = ((a.mass * vA) + (b.mass * vB)) / (a.mass + b.mass);

        //Angle of trajectory
        double theta = atan2(((a.mass * a.vY) + (b.mass * b.vY)), ((a.mass * a.vX) + (b.mass * b.vX)));

        std::cout  << "vA: " << vA << "   vB: " << vB << "\n" << velocity << "  :  " << theta << std::endl;

        //Vector components of velocity
        a.vX = velocity * cos(theta);
        a.vY = velocity * sin(theta);

        //Masses are added together (mass final)
        a.mass += b.mass;

        a.print();
    }

    void cleanup(std::list<celestial_body> &cb)
    {
        cb.clear();
        SDL_Quit();
    }

    int main(int argc, char** argv)
    {
        //Process arguments
        if(argc == 2)
        {
            int tempBodies = atoi(argv[1]);
            if(tempBodies > 0)
                bodies = tempBodies;
        }

        SDL_Surface* screen;
        SDL_Event event;

        std::list<celestial_body> cb;
        std::list<celestial_body>::iterator it;
        std::list<celestial_body>::iterator ij;

        double distance, diffX, diffY, tempAx, tempAy;
        double Fnet;    //Net Force on body
        double theta;   //Angle between two points in 2-D space
        double accel;   //Net acceleration of body

        bool quit = false;

        if(SDL_Init(SDL_INIT_EVERYTHING) == -1)
            quit = true;

        SDL_WM_SetCaption("Gravity Simulator", NULL);
        SDL_ShowCursor(1);

        screen = SDL_SetVideoMode(SCREEN_WIDTH, SCREEN_HEIGHT, SCREEN_BPP, SDL_HWSURFACE | SDL_DOUBLEBUF);
        if(screen == NULL)
            quit = true;

        for(int i = 0; i < bodies; ++i)
        {
            celestial_body newBody(i);
            cb.push_back(newBody);
            newBody.print();
        }

        std::cout << "Number of bodies: " << bodies << "\n\n";

        while(!quit)
        {
            while(SDL_PollEvent(&event))
            {
                if(event.type == SDL_QUIT)
                    quit = true;
                if(event.type == SDL_KEYDOWN)
                    if(event.key.keysym.sym == SDLK_ESCAPE)
                        quit = true;
            }

    //      SDL_FillRect(screen, &screen->clip_rect, 0x000000);

            for(it = cb.begin(); it != cb.end(); ++it)
                (*it).draw(screen);

            SDL_Flip(screen);

            for(it = cb.begin(); it != cb.end(); ++it)
            {
                tempAx = tempAy = 0;
                for(ij = cb.begin(); ij != cb.end(); ++ij)
                {
                    //Simply makes sure a body is not calculated against itself
                    if((*it).count != (*ij).count)
                    {
                        diffX = (*it).x - (*ij).x;
                        diffY = (*it).y - (*ij).y;

                        if(!SDL_CollideBoundingCircle((*it).x, (*it).y, (*it).radius, (*ij).x, (*ij).y, (*ij).radius, 0))
                        {
                            //Determine the distance between two bodies
                            distance = sqrt( (diffX * diffX) + (diffY * diffY) );

                            //Determine the net gravitational force between the two
                            Fnet = ((G * (*it).mass * (*ij).mass)) / (distance * distance);

                            //Determine the angle from a to b in 2-Dimensional space
                            theta = atan2(diffY, diffX);

                            //Determine the total acceleration
                            accel = Fnet / (*it).mass;

                            //Find acceleration from vector components
                            //Add them together for each affecting body
                            tempAx += -(accel * cos(theta));
                            tempAy += -(accel * sin(theta));
                        }
                        else
                        {
                            if((*it).mass > (*ij).mass)
                            {
                                merge((*it), (*ij));
                                ij = cb.erase(ij);
                            }
                            else
                            {
                                merge((*ij), (*it));
                                it = cb.erase(it);
                            }
                            std::cout << "Number of bodies: " << cb.size() << "\n";
                        }
                    }
                }
                (*it).aX = tempAx;
                (*it).aY = tempAy;
            }

            for(it = cb.begin(); it != cb.end(); ++it)
                (*it).update();
        }

        cleanup(cb);

        return 0;
    }

     
     
    Last edited: Aug 12, 2012
  2. jcsd
  3. Aug 15, 2012 #2

    chiro

    User Avatar
    Science Advisor

    Hey creativemfs and welcome to the forums.

    I've noticed you are calculating velocity vectors but you are not using forces per se: in other words, instead of introducing forces that are updated at every time interval and introduced upon either a collision event or by some non-contact force (like gravity) you are modifying the velocity vector directly.

    Usually what happens in a physics engine is that the velocity vector doesn't change unless a force is applied (just like Newtons Laws of motion). So in a contact situation you can calculate the force on each body and then use a numerical integrator (like Euler's method) to use the force over some small dt to get the new velocity and then use this new velocity to update the position of said object.

    So in short: you have three things for each object (vectors): position, velocity, acceleration. You only modify acceleration (unless you want to say do debugging, initialization or to hack your variables) and any introduced acceleration gets applied immediately within the time delta. The integrator will take the acceleration components and apply them within the time-delta and then set the acceleration back to zero. If a new event is called that applies more acceleration (like another collision or by having a non-contact force like gravity) then the next time step will incorporate that new force.

    This is the basic idea of how a Newtonian physics engine works: you typically look at forces in one form or another and then you apply these forces to generate velocity changes and finally changes in displacement.

    The question then remains: what kind of collision do you want? It seems like you were just averaging the two velocities together with the corresponding masses to be used as weights.

    If you want to retain this, then you should calculate the impulse and then consider that relative to the time frame (i.e. the time delta in your simulation) to get the force. (So if your impulse is I then calculate I/dt where dt is your delta to get the force).

    Now you take this force and then apply it to your object: I noticed you've merged objects together so you can just delete the other one and apply the force to the one in question.

    In terms of masses, you don't supply these to the numerical integrator: the only thing that the integrator should use is the forces: you can use the mass to create forces, but ultimately it's the job to just look at forces and how these translate into positions at each time-step.
     
  4. Aug 15, 2012 #3

    Mark44

    Staff: Mentor

    Michaela,
    Thanks for putting your code inside HTML code tags! It's rare that new members think to do this.
     
  5. Aug 15, 2012 #4
    thank you for your reply, I did make this modified version. I decided starting with a simpler program will be best.

    This is more like a black hole that sucks everything in.

    So chiro, basically I want two bodies that collide to stick together (mass adds together) and travel at the new velocity in x and y. So I don't know if I need to even calculate a 'theta' because an individual x and y velocity will be traveling at some angle.

    Can you give me an equation, or some code examples, for calculating the velocity of the collided bodies? How do I code impulse?

    I have only taken one physics class in college, but I remember a little about impulse, I will look it up in my book and see if I can figure out where to go.

    oh and your welcome mark. I spend a lot of time on cplusplus.com so code tags are nothing new to me :)
    Code (Text):

    /*
     * main.cpp
     *
     *  Created on: Jul 28, 2012
     *      Author: michaela
     *      Fnet = (G * m1 * m2) / r^2
     *
     *      Optional parameters ./gravity_simulator [--random] [num bodies]
     *
     *      As of now, only the main body affects the other bodies, they do not affect each other
     *      Click on the main body and move the mouse to move it, click again to release
     */

    #include "SDL_collide.h"
    #include <SDL/SDL.h>
    #include <SDL/SDL_gfxPrimitives.h>
    #include <list>
    #include <math.h>
    #include <iostream>

    const int SCREEN_HEIGHT = 600;
    const int SCREEN_WIDTH = 800;
    const int SCREEN_BPP = 32;

    const double G = 6.67e-1;       //Scaled down by e-10
    const double DELTA_TIME = 0.1;  //Arbitrary
    const double RADIUS = 3.0;

    //Hardcoded number of bodies, used if no argument is specified
    int bodies = 100;

    struct celestial_body
    {
            int count;
            double mass;
            double radius;
            double x, y;
            double vX, vY;
            double aX, aY;
            SDL_Color color;

            celestial_body(int newCount = 0, double newRadius = 0.5, double newMass = 0, int newX = 0, int newY = 0)
            {
                count = newCount;

                if(newMass == 0)
                    mass = (rand() % 1000) + 1;
                else
                    mass = newMass;

                radius = newRadius;

                if(newX == 0)
                    x = fmod((double)rand(), (double)SCREEN_WIDTH);
                else
                    x = newX;
                if(newY == 0)
                    y = fmod((double)rand(), (double)SCREEN_HEIGHT);
                else
                    y = newY;

                vX = vY = 0;
                aX = aY = 0;

                //Random color and avoid very dark colors
                color.r = (int)mass % (rand() % 10000) % 255;
                color.g = (int)mass % (rand() % 10000) % 255;
                color.b = (int)mass % (rand() % 10000) % 255;

                if(color.r <= 20 && color.g <= 20 && color.b <= 20)
                {
                    color.r += 25;
                    color.g += 25;
                    color.b += 25;
                }
            }
            void draw(SDL_Surface* screen)
            {
                filledCircleRGBA(screen, x, y, radius, color.r, color.g, color.b, 255);
            }
            void update()
            {
                vX += aX * DELTA_TIME;
                vY += aY * DELTA_TIME;

                x += vX * DELTA_TIME;
                y += vY * DELTA_TIME;
            }
            void print()
            {
                std::cout << "Count: " << count << "\nMass: " << mass << "\nX: " << x << "\nY: " << y <<
                             "\nColor: " << (int)color.r << " " << (int)color.g << " " << (int)color.b << "\n\n";
            }
    };

    void spawn_body(std::list<celestial_body>& cb, int newX = 0, int newY = 0)
    {
        static int count = 1;
        celestial_body new_body(count++, RADIUS, 0, newX, newY);
        cb.push_back(new_body);
        new_body.print();
    }

    void cleanup(std::list<celestial_body>& cb)
    {
        cb.clear();
        SDL_Quit();
    }

    int main(int argc, char** argv)
    {
        SDL_Surface* screen;
        SDL_Event event;

        std::list<celestial_body> cb;
        std::list<celestial_body>::iterator it;
        std::list<celestial_body>::iterator ij;

        double distance, diffX, diffY;
        double Fnet;    //Net Force on body
        double theta;   //Angle between two points in 2-D space
        double accel;   //Net acceleration of body

        bool quit = false;
        bool random = false;
        bool selected = false;

        //Process arguments
        if(argc == 3)
        {
            if(strcmp(argv[1], "--random") == 0)
                random = true;
            int tempBodies = atoi(argv[2]);
            if(tempBodies > 0)
                bodies = tempBodies;
        }

        if(SDL_Init(SDL_INIT_VIDEO) == -1)
            quit = true;

        SDL_WM_SetCaption("Gravity Simulator", NULL);
        SDL_ShowCursor(1);

        screen = SDL_SetVideoMode(SCREEN_WIDTH, SCREEN_HEIGHT, SCREEN_BPP, SDL_HWSURFACE | SDL_DOUBLEBUF);
        if(screen == NULL)
            quit = true;

        //Initialize main static body
        celestial_body main(0, 8.0, 1000, SCREEN_WIDTH / 2, SCREEN_HEIGHT / 2);
        main.print();

        if(random)
            for(int i = 0; i < bodies; ++i)
                spawn_body(cb);

        while(!quit)
        {
            while(SDL_PollEvent(&event))
            {
                if(event.type == SDL_QUIT)
                    quit = true;
                if(event.type == SDL_KEYDOWN)
                    if(event.key.keysym.sym == SDLK_ESCAPE)
                        quit = true;
                if(event.type == SDL_MOUSEBUTTONDOWN)
                {
                    if(SDL_CollideBoundingCircle(main.x, main.y, main.radius, event.button.x, event.button.y, 1, 0))
                        selected = !selected;
                    else
                        spawn_body(cb, event.button.x, event.button.y);
                }
            }

            if(selected)
            {
                main.x = event.motion.x;
                main.y = event.motion.y;
            }

            SDL_FillRect(screen, &screen->clip_rect, 0x000000);

            main.draw(screen);
            for(it = cb.begin(); it != cb.end(); ++it)
                (*it).draw(screen);

            SDL_Flip(screen);

            for(it = cb.begin(); it != cb.end(); ++it)
            {
                diffX = main.x - (*it).x;
                diffY = main.y - (*it).y;

                //Determine the distance between two bodies
                distance = sqrt((diffX * diffX) + (diffY * diffY));

                //Determine the net gravitational force between the two
                Fnet = ((G * (*it).mass * main.mass)) / (distance * distance);

                //Determine the angle from a to b in 2-Dimensional space
                theta = atan2(diffY, diffX);

                //Determine the total acceleration
                accel = Fnet / (*it).mass;

                //Find acceleration from vector components
                (*it).aX = accel * cos(theta);
                (*it).aY = accel * sin(theta);

                //Check for collision between new_body and main
                //main gains the mass of colliding body
                if(SDL_CollideBoundingCircle(main.x, main.y, main.radius, (*it).x, (*it).y, (*it).radius, 0))
                {
                    main.mass += (*it).mass;
                    it = cb.erase(it);
                }
            }

            for(it = cb.begin(); it != cb.end(); ++it)
                (*it).update();

            SDL_Delay(10);
        }

        cleanup(cb);

        return 0;
    }

     
     
    Last edited: Aug 15, 2012
  6. Aug 15, 2012 #5
    I see that (F*dT) / M = dV but how do I get this to apply to two colliding bodies?

    Also was my method of calculating the acceleration due to the gravitational pull of the other bodies and summing it correct?
     
    Last edited: Aug 15, 2012
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Gravity simulator in C++
  1. Quantum Simulation (Replies: 5)

  2. To C or not to C (Replies: 31)

  3. Physics simulation (Replies: 5)

  4. C or C++? (Replies: 8)

  5. Simulating the world (Replies: 14)

Loading...