Programming for Physics Simulations – How is it Done ?

A lot of people even those who have experience in programming have very wrong ideas about how programs work. Moreover the common notion among the general public is that software development only requires knowledge and experience in programming, however its actually quite the contrary. Software development can have very high level of specialization unlike many other fields including medical, for instance one experienced cardiologist can do the same another cardiologist can. How two highly experienced web developers may not be able to do what the other can. What people don’t understand is programming on its own is just a tool and nearly worthless on its own. Its like maths, whats the use of maths? No use. But apply math in physics, chem, astronomy …… then suddenly it becomes very useful. That’s how programming is too , People often specialize in different fields like maths, physics and  architecture. Like for example if a developer needs to develop a CAD program , he or she would need to have a good knowledge of architecture to effectively develop a suitable product, therefore they would need to specialize in architecture too.

So here in this post I intend to instill some ideas into the minds of new developers. To do this I will explain how to simulate a basic physics concept used in cut the rope from scratch through C++ code. I will keep it as simple as possible. However one needs basic knowledge of trigonometry to follow. Please note that all the programs embedded on this post is either HTML5 or flash, however I will be discussing from C++ point of view.

This is what we are gonna make using C++, Its pretty basic right ? Yeah it is. But even this is not very easy as people might feel, Many people will be thinking oh whats the big deal, you pull it, it follows , whats so hard? Indeed you pull it, it follows but nothing on a computer ever happens on its own, Someone should take care of the physics through code.

Ok so how do I think the right way ?

The right way to think is not to think in terms of code but in terms of logic, Imagine this as a physics situation, there are 2 balls, separated by a distance and an elastic rope joined between them,  larger the distance, faster the ball moves. However even if the rope becomes slack, due to inertia the ball continues to be in motion. So inertia has to be also considered. Gravity needs to be considered to make the thing realistic. Finally the fact that we are working in a 2D space has to be kept in mind.  

On further analysis, You might feel that there is a need for friction / resistance between the balls, or they would never stop moving

Physics solution :

This is the solution to the problem in physics point of view

Here 'T' is the tension force in the elastic rope. 'x' is the extension (elongation) of  the rope, 'A' is the rope’s cross sectional area, 'L' is the natural length of the rope while 'Y' is the Young's modulus.  This is the formula derived from Hooke's law.

We need to resolve these forces into vectors since we are working in 2D space. Then we need to add the gravitational force on the vertical component so it will be T sinθ + mg.

From here its easy, all one needs to do is calculate acceleration in both components using 'F = ma', then find change in velocity using 'v = u + at', And from velocity you can find new position of the balls. But you might have already realized both acceleration and velocities keep on changing , that means to correctly calculate the final velocity, you will need to use integral calculus which is beyond the scope of this post.

However this analysis is correct and can be implemented in programming

Programming Implementation :

Remember I said we need integration to calculate the velocities and position at a given time? But in programming we have a work around i.e  Running several iterations to approximately calculate the velocity and position . This means we split a time interval, say 1 sec into many smaller time intervals Say 0.1 s, Since 0.1s is a very small interval, we can assume the acceleration and velocity doesn't change in that interval . Accuracy can be increased to any desired extent by decreasing the small time interval.

Ok now we are finally ready to get started

Main variables

float x1=200,y1=200,x2=200,y2=200, s=100;
float vx2=0,vy2=0;

(x1,y1) are the initial position of the first ball , the ball which is controlled by the mouse. (x2,y2) is the position of the ball which is connected to the first ball. 'vx2' and 'vy2' are the components of velocity of the second ball in the X and Y directions respectively. 's' is the natural length of the rope.

The following bits of code has to be put in a loop that runs the iterations

1)  Getting the coordinates of the mouse and setting that as the coordinates of the first ball
 

x1=mousex(); 
y1=mousey();

2) Finding distance and angle between the 2 centres.

float d1=distance(x1,y1,x2,y2);
float ang1=angle(x1,y1,x2,y2);

Here distance() is a function which calculates distance between 2 points at an instant of time and angle() is a function which calculates angle between 2 points i.e. angle made by the line joining the 2 points and the horizontal.

3) Checking if the rope is stretched

We calculate distance between the 2 balls to know if the rope is stretched and angle to resolve force into components

//If the rope is stretched
if(d1-s>0) {   //Resolving tension into vectors
    vy2+=(d1-s)*sin(ang1)*.007;
    vx2+=(d1-s)*cos(ang1)*.007;
}

/*  here all the constants clubbed together as 0.007 , you can change it
We are adding the force directly to velocities by assuming balls to be of unit
mass and time interval to also be unit time
*/

4) Adding gravitational force

 vy2+=.23;
 //Gravitational force added on y coordinate
// 0.23 is an experimentally calculated value ( you can experiment too ! )

5) Moving the body and friction

We move the body as follows as we have assumed time interval to be unit time interval

y2+=vy2;
 x2+=vx2;
//Friction reduces velocity by 1%
vy2*=.99;
vx2*=.99;

      

6) Summing up

//draws 2 circles on the screen
circle(x1,y1,25);
circle(x2,y2,25);
//draws a line between them
line(x1,y1,x2,y2);

And we are done !!

Several improvements can be made like having parameters such as length, area, masses of spheres etc. Further more this concept can be expanded on more than 2 balls

Full code :

#include<graphics.h>

#include<math.h>

//Distance is a function which calculates Distance between 2 points

float distance(int x1,int y1,int x2,int y2)

{

    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));

}

//Angle is a function which calculates the angle between 2 points

// ( tan inverse of the Slope of the line joining the 2 points )

float angle(int x1 , int y1, int x2 , int y2)

{

    float ang;

    //In case the 2 points are perpendicular

    if(x1-x2==0)

    {

        ang=-1.57;

    }

    else

    {

        ang = atan((float)(y1 - y2) / (x1 - x2));

    }

    //Quadtrants

    if (ang < 0 && y1> y2)

    {

        ang += 3.14;

    }

    else if (ang > 0 && x1 < x2)

    {

        ang += 3.14;

    }

    return ang;

}

int main()

{

    //declaration of variables

    float x1=200,y1=200,x2=200,y2=200, s=100;

    float vx2=0,vy2=0;

    //Creating a new window

    initwindow (550, 400,"Windows BGI",50,50, false,true);

    while(1)

    {

        //Getting the coordinates of the mouse

        x1=mousex();

        y1=mousey();

        //Finding distance and angle between 2 points

        float d1=distance(x1,y1,x2,y2);

        float ang1=angle(x1,y1,x2,y2);

        //If the rope is stretched

        if(d1-s>0)

        {

            //Resolving tension into vectors

            vy2+=(d1-s)*sin(ang1)*.007;

            vx2+=(d1-s)*cos(ang1)*.007;

        }

        //Gravitational force on y coordinate

        vy2+=.23;

        //Now moving the body

        y2+=vy2;

        x2+=vx2;

        //Friction reduces velocity by 1%

        vy2*=.99;

        vx2*=.99;

        //draws 2 circles on the screen

        circle(x1,y1,25);

        circle(x2,y2,25);

        //draws a line between them

        line (x1,y1,x2,y2);

        //waits for 40 milliseconds

        delay(40);

        //something to do with memory

        swapbuffers();

        //clear the screen so that a new frame can be redrawn

        cleardevice();

    }

    return 0;

}