Sign in to follow this  
Concentrate

Simulating pendulum movement

Recommended Posts

I know there is something wrong here. I am trying to simulate a simple pendulum by am failing. Here is the relevancy.
struct Pendulum
{
	RigidBall ball;
	Line rope;
	const int LEN;	
	
	float changeInAngle;
	const float restAngle;

	Pendulum() : LEN( 6 ), restAngle( Constants::DEGTORAD * 270 )
	{
		rope = Line( Vector3D( 0.0f,-LEN,0.0f) );	
		rope.setWhiteColor();
		ball = RigidBall(1,1,0.0f,-LEN,0.0f);
		ball.setYellowColor();
		ball.setMass( 1.0f );		

		changeInAngle = 75 * Constants::DEGTORAD;
		
		ball.posVec = Vector3D( LEN * cos( ( restAngle + changeInAngle ) ),
								LEN * sin( ( restAngle + changeInAngle ) ),
								0.0f);			
	}
	void calculate(float deltaTime)
	{
		static const float EXTRA_FORCE = 5;

		static float g = Constants::GRAVITY / LEN * EXTRA_FORCE;		
		//F = mg*sin(theta) = ma;
		//a = F/ m
		//a = mgsin(theta) / m	
		//a = g*sin(theta)
		ball.accelVec = Vector3D( 0.0f, g * sin( changeInAngle ) , 0.0f );

		//V = v_0 + a*t
		ball.velVec += ball.accelVec * deltaTime;

		// X = v*dt + 1/2 a*dt^2
		ball.posVec += ball.velVec * deltaTime;// + 1.0f/2.0f * ball.accelVec * deltaTime * deltaTime;
		
		changeInAngle = atan2( ball.posVec.getY(), ball.posVec.getX() ) - restAngle;
		
		rope.setEndPoint( ball.posVec );
	}
	void draw(float deltaTime)
	{
		calculate(deltaTime);
		rope.draw();	
		ball.draw();
		//draw axis
		glBegin(GL_LINES);
			glColor3f(1.0f,0.0f,0.0f);
			glVertex3f(-7.0f,0.0f,0.0f);
			glVertex3f(7.0f,0.0f,0.0f);
			glVertex3f(0.0f,-7.0f,0.0f);
			glVertex3f(0.0f,7.0f,0.0f);
		glEnd();
		
	}
}pendulum;



I know some logic is wrong. Can you help ? When I run it the pendulum does not oscillate towards the x = 0 line. It just decays down the y-axis. after a few millisecond. [Edited by - Concentrate on December 25, 2009 10:23:42 PM]

Share this post


Link to post
Share on other sites
I though my calculation was correct but apparently not.

Can someone check my calculation.

void calculate(float deltaTime)
{
static const float EXTRA_FORCE = 5;

static float g = Constants::GRAVITY / LEN * EXTRA_FORCE;
//F = mg*sin(theta) = ma;
//a = F/ m
//a = mgsin(theta) / m
//a = g*sin(theta)
ball.accelVec = Vector3D( g * sin( changeInAngle ) ,g * cos(changeInAngle) , 0.0f );

//V = v_0 + a*t
ball.velVec += ball.accelVec * deltaTime;


// X = v*dt + 1/2 a*dt^2
ball.posVec += ball.velVec * deltaTime;// + 1.0f/2.0f * ball.accelVec * deltaTime * deltaTime;

changeInAngle = atan2( ball.posVec.getY(), ball.posVec.getX() ) - restAngle;

rope.setEndPoint( ball.posVec );
}




restAngle is 270 degrees in radians.

Share this post


Link to post
Share on other sites
It appears that you're calculating the force on the ball based on a change in angle.

The force on the ball due to gravity will be a constant m*g in the vertical direction, whatever the position of the ball. The acceleration vector will be (0,-m*g,0), and the change in velocity accelVec*deltaSecs, if the Y axis is pointing up.

Share this post


Link to post
Share on other sites
>>It appears that you're calculating the force on the ball based on a change in angle.

Yes.

>>The force on the ball due to gravity will be a constant m*g in the vertical direction, whatever the position of the ball. The acceleration vector will be (0,-m*g,0), and the change in velocity accelVec*deltaSecs, if the Y axis is pointing up.

Thanks. How will I account, for the horizontal force then ?

Share this post


Link to post
Share on other sites
The horizontal force comes from the rope. You're probably simulating it as a straight rod (which is simplest).

Divide the vertical force (m*g*deltaSecs) into a force along the direction of the rope and a force normal to the rope, according to the angle the rope makes with the vertical. The force along the direction of the rope will be cancelled by the force of the rope on the ball. The normal force (perpendicular to the rope) will cause acceleration normal to the rope.

Share this post


Link to post
Share on other sites
Quote:
Original post by Buckeye
The horizontal force comes from the rope. You're probably simulating it as a straight rod (which is simplest).

Divide the vertical force (m*g*deltaSecs) into a force along the direction of the rope and a force normal to the rope, according to the angle the rope makes with the vertical. The force along the direction of the rope will be cancelled by the force of the rope on the ball. The normal force (perpendicular to the rope) will cause acceleration normal to the rope.



Is that tangential force
 mg*sin(theta)?

Share this post


Link to post
Share on other sites
If, by "tangential," you mean perpendicular to the direction of the rope, then yes. And if your "theta" is the angle of the rope with respect to the Y axis, Y pointing down, then that force would be m*g*sin(theta).

Make sure you apply the sign of the angle correctly - the force will always tend to swing the ball back toward that Y axis.

Share this post


Link to post
Share on other sites
Sorry for the late response.

I think my calculations are correct but the result shows otherwise. I am trying
what you said, but its acting weird.

This is what I have now :

void calculate(float deltaTime)
{
static const float EXTRA_FORCE = 10;

static float g = Constants::GRAVITY / LEN * EXTRA_FORCE;

//F = mg*sin(theta) = ma;
//a = F/ m
//a = mgsin(theta) / m
//a = g*sin(theta)
ball.forceVec = Vector3D( ball.getMass() * g * sin(changeInAngle) , ball.getMass() * g , 0.0f);

ball.accelVec = ball.forceVec/ ball.getMass();

//V = v_0 + a*t
ball.velVec += ball.accelVec * deltaTime;


// X = v*dt + 1/2 a*dt^2
ball.posVec += ball.velVec * deltaTime + 1.0f/2.0f * ball.accelVec * deltaTime * deltaTime;


changeInAngle = atan2( ball.posVec.getY(), ball.posVec.getX() ) - restAngle;

rope.setEndPoint( ball.posVec );
}



Here is whats happening , http://www.youtube.com/watch?v=uzPGJD39zPI

Share this post


Link to post
Share on other sites
pendulum

Maybe this picture will help. In this picture, theta is the angle between the Y-axis and the rope.

mg points down along the Y-axis (always). The tangential force [mg sin(theta)] is not opposed and is what causes the acceleration. The force on the ball from the rope can only be along the direction of the rope and is always equal to mg*cos(theta).

Divide the rope force into x- and y-components (shown on the right):

x-force = -mg cos(theta) sin(theta)
y-force = -mg cos(theta) cos(theta)

At an angle theta, the force on the ball is:

Vector3D( -m*g*cos(theta)*sin(theta), m*g - (-m*g*cos(theta)*cos(theta)), 0);

The Y component = mg + (1-cos^2(theta)). 1-cos^2 = sin^2, so:

force = Vector3D( -m*g*cos(theta)*sin(theta), m*g*sin(theta)*sin(theta), 0);

Testing some values-

When the ball is hanging straight down, theta = 0 and force = (0,0,0). That is, there is no tangential force and the rope supports the ball. The velocity is not affected.

When the ball is straight out to the right, theta = pi/2 and force = (0, mg, 0). Gravity is pulling straight down on the ball and acceleration is maximum.


Share this post


Link to post
Share on other sites
Thanks for the explanation. Shouldn't the second triangle (the one on the right of your picture) be the same as the first one ( the one on the left of your picture) because I thought they should be similar triangles?

I must be doing something wrong because I adjusted the force from your derivation
and its similar to what I had in the youtube video (plus/minus).

This is what I have for the current definition,

struct Pendulum
{
RigidBall ball;
Line rope;
const int LEN;

float changeInAngle;
const float restAngle;

Pendulum() : LEN( 6 ), restAngle( Constants::DEGTORAD * 270 )
{
rope = Line( Vector3D( 0.0f,-LEN,0.0f) );
rope.setWhiteColor();
ball = RigidBall(1,1,0.0f,-LEN,0.0f);
ball.setYellowColor();
ball.setMass( 1.0f );

changeInAngle = 75 * Constants::DEGTORAD;

ball.posVec = Vector3D( LEN * cos( ( restAngle + changeInAngle ) ),
LEN * sin( ( restAngle + changeInAngle ) ),
0.0f);
}
void calculate(float deltaTime)
{
static const float EXTRA_FORCE = 10;

static float g = Constants::GRAVITY / LEN * EXTRA_FORCE;
static float m = ball.getMass();
float theta = changeInAngle;
//F = mg*sin(theta) = ma;
//a = F/ m
//a = mgsin(theta) / m
//a = g*sin(theta)
ball.forceVec = Vector3D( -m*g*cos(theta)*sin(theta), m*g*sin(theta)*sin(theta), 0);


ball.accelVec = ball.forceVec/ ball.getMass();

//V = v_0 + a*t
ball.velVec += ball.accelVec * deltaTime;


// X = v*dt + 1/2 a*dt^2
ball.posVec += ball.velVec * deltaTime + 1.0f/2.0f * ball.accelVec * deltaTime * deltaTime;


changeInAngle = atan2( ball.posVec.getY(), ball.posVec.getX() ) - restAngle;

rope.setEndPoint( ball.posVec );
}
void draw(float deltaTime)
{
calculate(deltaTime);
rope.draw();
ball.draw();

}
}pendulum;




I don't think its any problem with my display function because I do not
manipluate anything except the time in there.

Share this post


Link to post
Share on other sites
The triangle on the left shows the decomposition of the mg vector into a tangent vector and a vector in the direction of the rope. The triangle on the right shows the decomposition of the vector in the direction of the rope.

Your calculation of the angle theta and ball.posVec are both incorrect. theta = atan2( ball.x, ball.y ). Don't add " -restAngle." theta is the angle of the rope with the y-axis as shown in the picture.

I think you're incorrectly adding the acceleration to both the velocity vector and the position vector.

I'm pretty sure it should be:
pos += vel*time + 1/2 a*time*time

then
vel += a*time

Share this post


Link to post
Share on other sites
Quote:
Original post by Buckeye

Your calculation of the angle theta and ball.posVec are both incorrect. theta = atan2( ball.x, ball.y ). Don't add " -restAngle." theta is the angle of the rope with the y-axis as shown in the picture.


I was doing this just like the spring-mass, but instead of the change in
position, its change in angle. So I thought that the rest angle is 270 degrees,
so if there is a movement then there is a change in angle, which is
currentAngle - restAngle? I am not sure why I need to remove the restAngle, but
i'll give it a try.

Quote:

I think you're incorrectly adding the acceleration to both the velocity vector and the position vector.
I'm pretty sure it should be:
pos += vel*time + 1/2 a*time*time

then
vel += a*time


Isn't that what I have? Maybe I should use a different approach? And by the way, thanks for your responses. You seem to be the only helpful one around here, usually other people are more friendly and tend to help out as well. Thank you.

Share this post


Link to post
Share on other sites
You can include the restAngle if you want, but it's just a constant offset, so why complicate the calculation? As long as the angle you use for the calc is the angle between the rope position and the Y-axis, it's all the same.

You have:
vel_1 = vel_0 + a*time
pos += vel_1*time + 1/2a*time*time

which is
pos += (vel_0 +a*time)*time + 1/2a*time*time

so you've accounted for the acceleration's effect on the position twice in one time interval. I'm suggesting you reverse the order of those two equations. I probably wasn't clear on that.

That is:
you have-
vel += a*time
pos += vel*time + 1/2 a*time*time

I think-
pos += vel*time + 1/2 a*time*time
vel += a*time

is correct.

Share this post


Link to post
Share on other sites
Thanks for the explanation. I tried it but it won't budge. You know I am a little
disappointed. I think we did the calculations correctly, we implemented it
correctly, but the physics simulation says otherwise. If you wan't and have time, I can zip it to you so you can play around with it. Although I will admit, that its not pretty.

Share this post


Link to post
Share on other sites
I didn't properly take into account the length of the rope.

The derivation of the equation is based on using cylindrical coordinates as cartesian coordinates are very difficult to work with in this case.

This is the update function. "fSecs" is your deltaTime or time-per-frame. theta0 is the angle of the rope with the y-axis at time zero.

You can, of course, make gravity and/or mass and/or theta0 whatever you like. Start out with float time = 0.

float theta;
float theta0 = 1.5708f; // pi/2 - horizontal
float grav = 32.0f;
float mass = 1.0f;
float ropeLen = 10.0f;
float cosFac = sqrt(grav/mass)*time;
theta = theta0 * cos(cosFac);
pendLoc.x = ropeLen*sin(theta);
pendLoc.y = -ropeLen*cos(theta);
time += fSecs;

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this