Ok i started to do code a physic simulation that i want to implement after in my 3d engine. Right now i have coded all linear part... and i kinda stuck on the rotation part...
heres the formulas im using so far
x(t+dt) = x(t) + v(t) * dt
p(t+dt) = ( M * v(t) ) + ( F(t) * dt )
v(t+dt) = p(t+dt)/M
where M = mass
dt = delta Time
x(t) = position(vector) at time t
x(t+dt) = position(vector) at time t + dt
v(t) = velocity(vector) at time t
v(t+dt) = velocity(vector) at time t + dt
p(t+dt) = linear momentum(vector at time t + dt
F(t) = force(vector) at time t
and heres my object class representing this:
class Object
{
public:
//mass
float mass;
//position
float x;
float y;
float z;
//velocity
float vx;
float vy;
float vz;
//linear momentum
float px;
float py;
float pz;
//force
float fx;
float fy;
float fz;
public:
void ComputeNext( float deltaTime );
};
inline void Object::ComputeNext( float deltaTime )
{
//considering the force and mass are constant...
//compute the position at time + deltaTime
//x(t+dt) = x(t) + v(t)*dt
this->x += this->vx * deltaTime;
this->y += this->vy * deltaTime;
this->z += this->vz * deltaTime;
//these calculations will be used for the next iteration...
//compute the linear momentum at time + deltaTime
//p(t+dt) = M * v(t) + F(t) * dt
this->px = ( this->mass * this->vx ) + ( this->fx * deltaTime );
this->py = ( this->mass * this->vy ) + ( this->fy * deltaTime );
this->pz = ( this->mass * this->vz ) + ( this->fz * deltaTime );
//compute the linear velocity at time + deltaTime
//v(t+dt) = p(t+dt)/M
this->vx = this->px/this->mass;
this->vy = this->py/this->mass;
this->vz = this->pz/this->mass;
}
i have tested my code running a little simulation with this code:
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include "object.h"
int main()
{
Object obj;
//affect random initial state
srand( (unsigned)time( NULL ) );
obj.mass = rand() % 100;
obj.x = rand() % 10;
obj.y = rand() % 10;
obj.z = rand() % 10;
obj.vx = rand() % 10;
obj.vy = rand() % 10;
obj.vz = rand() % 10;
obj.px = obj.vx/obj.mass;
obj.py = obj.vy/obj.mass;
obj.pz = obj.vz/obj.mass;
//constant gravity force
obj.fx = 0;
obj.fy = -9.8f * obj.mass;
obj.fz = 0;
float t = 0.0f;
printf( "Constants:\n==========\nMass:%0.2f\nForce:(%0.2f,%0.2f,%0.2f)\n", obj.mass, obj.fx, obj.fy, obj.fz );
do
{
printf( "\nTime:%0.2f\tPos:(%0.2f;%0.2f;%0.2f)\tVelocity:(%0.2f;%0.2f;%0.2f)\tL. Momentum:(%0.2f,%0.2f,%0.2f)", t, obj.x, obj.y, obj.z, obj.vx, obj.vy, obj.vz, obj.px, obj.py, obj.pz );
float dt = 0.02f;
obj.ComputeNext( dt );
t += dt;
} while( t < 1 );
int c = scanf("\n\nPress any key to quit");
return 0;
}
so this code simulates the animation from 0 to second 1 at 50fps (0.2sec per frame) and it seemed to work fine
so now i want to include the rotation in it and i dont know really what formulas to use to do this...
i have read the papers on
Physic based modelling of SYGGRAPH 95 but im kinda confused about how to calculate the rotation from the angular momentum and the Inertia tensor... well basicly it is because i dont understand his notation(like the dot to represent the derivative) and his array state in his implentation which doesnt seem to take into account the previous position,velocity,etc. when he computes the final values... thats why i have been reformulationg a bit the formulas from the linear part so i could get what i want but the angular part is much more harder to figure out... so anyone could try to help me out to get the formulas or some code that illustrates how to do it... what i need this:
R(t+dt) = R(t) + w(t)*R(t)*dt (i think it should be something like this but not sure.)
w(t+dt) = w(t) + I(t)L(t)*dt (again im not sure about that one...)
L(t+dt)?
EDIT: ok trying to base mysrlf on the equations i had for the linear part i cam up with something like L(t+dt) = I(t)*w(t) + T(t)*dt (is that correct?) and could w(t+dt)= L(t+dt) * I-1(t)?
where dt= deltaTime
R(t) is rotation(matrix) at time t
R(t+dt) is the rotation(matrix) at time t + dt
w(t) is angular velocity(vector) at time t
w(t+dt) is angular velocity(vector) at time t + dt
L(t) is angular momentum(vector) at time t
L(t+dt) is angular momentum(vector) at time t + dt
I(t) is inertia tensor(matrix) at time t
I-1(t) is inertia tensor matrix inverse at time t
T(t) is Torque at time t
actually the rotation formula with a quarternion would be better...
Yann L PWNS Carmack
[edited by - sross on June 28, 2003 6:18:19 PM]
