for Games and Simulations
John Schultz
Updated 11/26/2006 12:17:23 AM -0800
Locally Deformable Rigid Bodies with Kinetic Energy Correction
The
evolution of the physics system went something like this:
1. Pure impulse based with bisection to find time of first contact (Baraff,
etc.).
2. (1) with (momentum) weighted averaging of contacts. Started to get slow
(bisection).
3. Tried Jacobsen's rigid body (RB) Verlet suggestions. Reasonably fast, easy to
work with, but springy stacking.
4. Started adding impulses into a hybrid Verlet system.
5. Added simple system to correct kinetic energy so that Verlet RB's (VRB) can
bounce on the ground without losing any energy (coeffrest 1.0). Similar idea to
allow two VRB's to collide without losing energy.
6. Created a simple method to match up a Verlet particle tetrahedron to a
standard inertia tensor so that switching between integrators in real-time was
possible. This was also very useful in validating a rotating RB with various
standard methods such as Euler, Runge-Kutta, etc. The VRB (tetrahedral particle
system) exhibits more or less correct behavior (precession, etc.), while some
standard RB integration methods do not show proper precession behavior, especial
for higher angular velocities. VRB systems start to become less accurate with
higher angular velocities and larger time steps, so switching to a more accurate
method in this case is useful for realism.
7. Tried the Guendelman, Bridson, Fedkiw "GBF" (Nonconvex Rigid Bodies with
Stacking) concepts.
8. I found that if I allowed the RB's to locally deform using the Verlet
tetrahedron, then at the end of the time step set all the tetrahedrons back to
non-deformed, a reasonable solution was possible without iterating. In earlier
methods, the VRB would be relaxed using the standard Verlet constraint particle
method, with net velocity effects changing over time.
9. I then stopped switching between VRB and standard-style RB integrators (I
still switch between two integrators depending on angular velocity).
10. The system works fine with or without GBF style prediction. For each overlap
contact pair projection resolution (deforms tetrahedron and affects RB's
rotation and position), kinetic energy is corrected based on the effective
change in all the new particle positions.
I have not had a chance to optimize this system much, but it's very simple and
requires no complex matrix solvers, etc. Pretty much everything is based on
concepts from physics and particle systems as opposed to non-physics-based
optimization/matrix concepts. Constraints could use more work (where
optimization methods are popular and pretty solid), but it's good enough for my
current driving game. Constraints also cause a small amount of local RB
distortion (improves stability).
The local rigid body distortion concept can probably be taken further and
modeled better than my system; it matches fairly well to what happens in
reality, as everything in the final analysis is rubber. For example, there may
be a better way to model/handle the local distortions, energy correction, etc. I
currently always draw the rigid bodies with orthogonal matrices, though in the
past it was interesting to use a distorted matrix to visualize the RB distortion
(could be useful in a more generalized system for blending between completely
rigid and softer objects). A more advanced system (would require some iteration)
could model semi-soft and hard RB's (which require no iteration).
Kinetic Energy Correction for a Single Rigid Body against Terrain:
flt getTotalKineticEnergy(RigidBody & rb1) {
flt totalKE = 0.f;
for (int i=0; i < 4; i++) {
PhysicalParticle & p1 = physicalParticles[rb1.physicalParticles[i]];
p1.vel = p1.cur - p1.old;
totalKE += halfDeltaTime2Inv*p1.mass*p1.vel.lengthSqr();
} // for
return totalKE;
} // getTotalKineticEnergy
void setTotalKineticEnergy(RigidBody & rb1,flt totalKE,flt newTotalKE) {
flt scale;
if (totalKE < .0000001f) scale = 1.f;
else scale = rSqrt(newTotalKE/totalKE);
for (int i=0; i < 4; i++) {
PhysicalParticle & p1 = physicalParticles[rb1.physicalParticles[i]];
p1.vel = p1.cur - p1.old;
p1.vel *= scale;
p1.old = p1.cur - p1.vel;
} // for
} // setTotalKineticEnergy
Kinetic Energy Correction for Two Colliding Rigid Bodies after projection resolution (particles form one tetrahedron per RB):
// Get total KE before projection resolution
flt getTotalKineticEnergy(RigidBody & rb1,RigidBody & rb2) {
flt totalKE = 0.f;
for (int i=0; i < 4; i++) {
PhysicalParticle & p1 = physicalParticles[rb1.physicalParticles[i]];
PhysicalParticle & p2 = physicalParticles[rb2.physicalParticles[i]];
p1.vel = p1.cur - p1.old;
p2.vel = p2.cur - p2.old;
totalKE += halfDeltaTime2Inv*p1.mass*p1.vel.lengthSqr();
totalKE += halfDeltaTime2Inv*p2.mass*p2.vel.lengthSqr();
} // for
return totalKE;
} // getTotalKineticEnergy
// Correct total KE after projection resolution
void setTotalKineticEnergy(RigidBody & rb1,RigidBody & rb2,flt totalKE,flt newTotalKE) {
flt scale;
if (totalKE < .0000001f) scale = 1.f;
else scale = rSqrt(newTotalKE/totalKE);
for (int i=0; i < 4; i++) {
PhysicalParticle & p1 = physicalParticles[rb1.physicalParticles[i]];
PhysicalParticle & p2 = physicalParticles[rb2.physicalParticles[i]];
p1.vel = p1.cur - p1.old;
p2.vel = p2.cur - p2.old;
p1.vel *= scale;
p2.vel *= scale;
p1.old = p1.cur - p1.vel;
p2.old = p2.cur - p2.vel;
} // for
} // setTotalKineticEnergy
void ContactInfo::applyContactDisplacement(void) {
Vec3 curRB1,currentContactPoint,currentContactNormal;
getCurrentContact(curRB1,currentContactPoint,currentContactNormal);
flt penetrationDist = currentContactNormal ^ (curRB1 - currentContactPoint);
if (penetrationDist >= 0.f) {
return;
} // if
Vec3 delta = currentContactNormal * (-penetrationDist);
if (!rb2 || (rb2->activationType == RTPhysicalProps::ACTIVATION_TYPE_STATIC)) {
rb1->displaceAdjust(curRB1,delta); // <TODO> When friction is off, this can gain energy. <NOTE> Energy comment applies to older version only.
// rb1->updateRigidTransform(); // Must update: not using just bary transforms. <NOTE> Currently performed in displaceAdjust().
} else { // Two Dynamic Objects
// Pretty much has to be based on mass for physics to work correctly.
flt totalMassInv = 1.f/(rb1->linMass + rb2->linMass);
Vec3 delta1 = delta * rb2->linMass * totalMassInv; // Weighted adjustment based on total mass only.
Vec3 delta2 = delta * rb1->linMass * (-totalMassInv); // - = negate sign on scalar component for speed.
// <ENERGY CORRECTOR>
flt totalKEBefore = rb1->physicsSystem->getTotalKineticEnergy(*rb1,*rb2);
rb1->impulseAdjust(curRB1,delta1); // This system works OK: not a cause of stacking jitter (sleep system issue).
rb2->impulseAdjust(curRB1,delta2); // impulseAdjust() displaces and distorts the Verlet tetrahedron using a weighted barycentric method (see Jacobsen's paper).
rb1->updateLinearConstraints(); // Restore tetrahedron (using Newton-Raphson sqrt approx. per Jacobsen).
rb2->updateLinearConstraints(); // ""
// Keep top-most objects from losing too much energy: don't apply KE adjust unless multi-object-contact.
if ((rb1->multiObjectContact || rb2->multiObjectContact) && !(rb1->getCar() || rb2->getCar())) {
rb1->multiObjectContact = rb2->multiObjectContact = true;
flt totalKEAfter = rb1->physicsSystem->getTotalKineticEnergy(*rb1,*rb2);
if (totalKEAfter > totalKEBefore*1.1f) { // 1.1: this always appears to work best.
rb1->physicsSystem->setTotalKineticEnergy(*rb1,*rb2,totalKEAfter,totalKEBefore);
} // if
} // if
rb1->updateRigidTransform(); // Must update: not using just bary transforms.
rb2->updateRigidTransform(); // ""
} // if
} // ContactInfo::applyContactDisplacement
inline void setParticleVelocity(PhysicalParticle & p,const Vec3 & centerOfMass,const Vec3 & linVel,const Vec3 & rotVel) {
p.vel = (rotVel % (p.cur - centerOfMass)) + linVel; // % = cross product
p.old = p.cur - p.vel;
} // setParticleVelocity
void RigidBody::setVelocity(const Vec3 & _linVel,const Vec3 & _rotVel) {
for (int i=0; i < 4; i++) {
PhysicalParticle & p = physicsSystem->physicalParticles[physicalParticles[i]];
setParticleVelocity(p,pos.t,_linVel,_rotVel);
} // if
linVel = _linVel;
rotVel = _rotVel;
linM = linVel*linMass;
rotM = scaleLocal(rotVel,rotMassIdeal,pos);
} // RigidBody::setVelocity
// After obtaining momentum, you can convert to velocity as needed.
void RigidBody::particleMomentum(Vec3 & currentLinM,Vec3 & currentRotM) {
PhysicalParticle & p0 = physicsSystem->physicalParticles[physicalParticles[0]];
PhysicalParticle & p1 = physicsSystem->physicalParticles[physicalParticles[1]];
PhysicalParticle & p2 = physicsSystem->physicalParticles[physicalParticles[2]];
PhysicalParticle & p3 = physicsSystem->physicalParticles[physicalParticles[3]];
p0.vel = p0.cur - p0.old;
p1.vel = p1.cur - p1.old;
p2.vel = p2.cur - p2.old;
p3.vel = p3.cur - p3.old;
Vec3 p0massVel = p0.vel*p0.mass;
Vec3 p1massVel = p1.vel*p1.mass;
Vec3 p2massVel = p2.vel*p2.mass;
Vec3 p3massVel = p3.vel*p3.mass;
// Linear Momentum
currentLinM = p0massVel;
currentLinM += p1massVel;
currentLinM += p2massVel;
currentLinM += p3massVel;
Vec3 cm = (p0.cur+p1.cur+p2.cur+p3.cur)*.25f;
Vec3 p0r = p0.cur - cm; // Subtract center of mass.
Vec3 p1r = p1.cur - cm; // Subtract center of mass.
Vec3 p2r = p2.cur - cm; // Subtract center of mass.
Vec3 p3r = p3.cur - cm; // Subtract center of mass.
// Angular Momentum
currentRotM = p0r % p0massVel; // % = cross product
currentRotM += p1r % p1massVel;
currentRotM += p2r % p2massVel;
currentRotM += p3r % p3massVel;
} // RigidBody::particleMomentum
Spatial Field-CFD (SF-CFD): blended spatial functions effecting particle and rigid body motion
If fast and simple particle motion is desired for fluid
simulation, it should be possible to get away with a simple distance function
from the disturbance, where chaotic, rotational motion increases with the
distance from the disturbance center. For example, as the sub moves along near
the ocean floor, your particles look themselves up in your (borrowing a term
from 3DSMax) SpaceWarp function, and their motion is
adjusted. Particles near the sub hull simply displace (add velocity) away from
the hull, particles further behind the sub are rotated (given velocity
increments supporting the rotation of the field) in vortices left/right based on
their positions.
Thus, as the sub moves through space, particles behind the sub are disturbed in
a reasonably realistic dual-vortex fashion (due to say a twin-screw prop system;
only one main vortex for a single screw). Add another noise lookup function to
perturb flow even more, farther away (to simulate near completely turbulent,
chaotic flow (Reynolds
number > 4000)). Effect magnitude is scaled by sub velocity. The particles
could also be given normals and rotation to effect their flow as if they were
little flat plates (or rigid-body tetrahedrons) (another idea: 'macro particles
(sphere influence)' stored in a 3D spatial hash which store linear and angular
velocity; each particle would look itself up in the 3D spatial hash to see if it
is effected by the 'macro particle's' sphere of influence (picking up some of
the 'macro particle's influence based on simple sphere distance)).
In this way you should be able to very quickly (and hopefully easily) produce
good-enough-for-a-real-time-game fluid motion particle effects without using a
complex Navier-Stokes/grid-system (just a few blended distance functions
(perhaps see 3DSMax SDK SpaceWarp code examples for more info (WSM:
World Space Modifiers))).
If you use SF-CFD along with an implicit volume for the
sub (ellipsoid or rounded cylinder), it will be intuitively easy to implement,
visualize, and debug. If using macro particles (small rigid body spheres) you
can increase apparent pixel density by rendering alpha-blended rotating
billboards of hundreds or thousands of particles (a texture mapped sphere might
also work OK, though the curvature may cause issues with the edges. In the case
of the billboard: the texture map would be created with many particles of
various sizes, alpha-densities, and have a roughly circular perimeter). Texture
coordinate tricks can also be used to swirl and distort the rendered billboard
(instead of a simple quad, use a grid with sufficient density to perturb the
texture in a pleasing manner).
(The following
method can be implemented very easily and quickly given an existing rigid body
simulator)
1. Create a FluidSample class/struct:
struct FluidSample {
Vec3 defPos; // Local space, relative to center of mass (and center of
rotation).
Vec3 curPos; // Transformed to world space.
Vec3 defNormal; // Surface normal at fluid sample point, local space.
Vec3 curNormal; // Surface normal at fluid sample point, world space.
float scale; // This can be anything: area, volume, hand tuned.
};
2. Create FluidSamples for each vertex of the boat hull, or create new
ones at the centroid of the triangles, or create by custom means (by hand,
algorithmically, etc.).
3. Add an operational rigid body integrator (Symplectic Euler (compute
velocity first, use to update position) is fine).
4. Put the boat in the water: first tests use a plane for the water
surface.
5. Transform the FluidSamples to world space. The center of mass/rotation
should be closer to the bottom of the boat for stability. Initial tests can use
a cube or tetrahedron with center of mass at object center.
6. Loop through all of the FluidSamples:
If the sample is underwater, apply a force (point of application: curPos):
a. If the FluidSamples represent volumes, the force is constant (multiplied by
scale).
b. If the FluidSamples represent areas, the force is scale*depth. This
implements a spring.
c. Compute the velocity of curPos (RigidBody member function):
Vec3 RigidBody::pointVel(const Vec3 & p) const {
return (rotVel.cross(p - centerOfMass)) + linVel;
} // RigidBody::pointVel
d. Given the curPos velocity and curNormal, you can now compute fluid
interaction: e.g. flat plate force => 1/2*rho*V^2*angle applied in the direction
of curNormal, surface/parasitic drag (in the direction of velocity projected to
the plane defined by curPos and curNormal). This step will introduce
drag/damping so that more complicated methods (such as dealing with rotational
fluid effects and associated energy loss) aren't necessary for realistic looking
behavior. Forces are added for each FluidSample (point of application: curPos):
void RigidBody::addWorldForce(const Vec3 &
worldPos,const Vec3 & worldForce) {
linForce += worldForce;
rotForce += (worldPos - centerOfMass).cross(worldForce);
} // RigidBody::addWorldForce
To keep it simple, don't process fluid interaction if the surface points away
from the velocity.
7. Integrate the summed forces and goto step 5.
This will create highly realistic boat simulation (including "getting up on
step" for a speed boat). Rudders, skegs, keels, hydrofoil wings, are all modeled
the same way. I used this method in 1989 when starting development of a flight
simulator for the Amiga: this method is also suitable for aircraft and computing
aero effects on cars and other objects (hence the generalized FluidSample). When
using Symplectic Euler, it can be made extremely stable, even with a dt of 1/30,
and 100% fixed point math (would easily run on today's mobile hardware).
Boats are relatively easy to simulate. Ranked in order of progressive difficulty
for realistic physical simulation for games:
1. Spacecraft.
2. Boats.
3. Aircraft.
4. Cars.
5. People/Creatures (autonomous).
You can build a simple "paper" airplane with fluid samples
for:
1. Left wing
2. Right wing
3. Fuselage in front of CG ("cylindrical" fluid sample)
4. Fuselage behind CG ("cylindrical" fluid sample)
5. Vertical stabilizer
6. Left horizontal stabilizer
7. Right horizontal stabilizer
To add control surfaces:
8. Left aileron
9. Right aileron
10. Rudder
11. Left elevator
12. Right elevator
(6&11, and 7&12 can be combined into a stabilator or taileron, the same concept
can be applied to 5&10).
The dot products scaling the normals will do the right thing for the "flat plate
/ airfoil" samples: they'll work correctly for all attitudes: inverted flight,
turning, stalling, (attempting to) fly backwards (proper tail slide behavior),
etc.
Once you have the math worked out for simple flat plate aerodynamics (wrapped up
in a FluidSample), you can then populate your rigid body with FluidSamples to
get generalized behavior for airplanes, boats, cars, flying animals, etc. You
won't have to worry about anything more complicated (unless your application
requires it): IJW.
Re: Amiga: 32/64-bit fixed point integer (1<<14 == 1.0) was used (no floating
point). The FluidSample method is made stable by clamping net forces to a stable
range.
For a simple flat plate "wing", you can just apply the
resulting force in the direction of the surface normal (for the FuildSample).
The resulting vector contains both lift and ("induced") drag. Lift is typically
referred to relative to flow velocity (perpendicular), drag (induced+parasitic,
etc.) in the direction of flow velocity, weight opposite lift, etc.
A typical airplane has its wing set with a small
angle of incidence
so that enough lift is generated while cruising the pilot does not need to
constantly pull back on the stick/yoke to keep the plane flying level.
If you draw the velocity and normal vectors on paper, you can see that you can
totally adjust wing efficiency/behavior by modifing net lift+drag. Wing shape
simply dictates how you should modify efficiency (+/- lift, +/- drag). You can
study lift+drag plots for various wing shapes to help build lookup tables (or
model in other ways).
// % = cross product, ^ = dot product, ~ = normalized copy.
inline Vec3 Vec3::perp(const Vec3 & V) const {
return *this % V % *this;
} // Vec3::perp
// Project V onto normalized *this.
inline Vec3 Vec3::proj(const Vec3 & V) const {
return *this * (*this ^ V);
} // Vec3::proj
// rb.pos is a Mat3t: rotation + position (t = Vec3 position).
// FluidSample.cpp
// Created 12/25/03 John Schultz
// Copyright 2006 Brightland Corporation
#include "physicsSystem.h"
void FluidSample::update(void) {
RigidBody & rb = *physicsSystem->motionObjects[motionObjectIndex]->getRigidBody();
// Calc velocity from linVel, rotVel, and relative pos (pointVel).
pos = rb.pos ^ offset; // offset is a Mat3t in rb local space, where the y column represents the surface normal up direction.
Vec3 relPos = pos.t - rb.pos.t;
vel = rb.linVel + (rb.rotVel % relPos);
vel.clampLength(MPH_TO_MPS(350.f)); // Velocity clamping is important for absolute stability.
Vec3 normal;
if (isCylinder) {
normal = pos.zCol().perp(vel);
flt normLen = normal.length();
if (normLen < .001f) {
normal = 0.f;
} else {
normal *= 1.f/(-normLen);
} // if
} else {
normal = pos.yCol();
} // if
Vec3 vel2 = vel*vel.length();
flt alpha = (vel2 ^ normal)*(-.5f);
force = normal * (alpha * area * physicsSystem->getDeltaTime2Inv());
Vec3 inducedDragDir = ~(-vel);
Vec3 dragOffset = inducedDragDir.proj(force)*efficiency;
force -= dragOffset;
flt forceLen = force.length();
force.clampLength(1e9f); // Not really needed with velocity clamping: after tuning, only one clamp is needed.
rb.addWorldForce(pos.t,force);
} // FluidSample::update
// FluidSample.cpp
Use a Symplectic Integrator, such as a Newton-Störmer-Verlet ( NSV) variant:
v += a*dt
x += v*dt
See Hairer et al's book, Geometric Numerical Integration Structure-Preserving Algorithms for Ordinary Differential Equations.
Use the Velocity Verlet Variant to reduce error if necessary:
flt oldAccel = a; //0.f; // May be self-starting issues... x = v = 0.f; for (i=0; i < count; i++) { x += v*dt + .5f*oldAccel*dt*dt; v += .5f*(oldAccel+a)*dt; oldAccel = a; lprintf("V: %f X: %f\n",v,x); } // for
Example test code for various simple integrators:
struct
Test {
Test() {
int count =
100;
flt x = 0.f;
flt v = 0.f;
flt a = 9.81f;
flt dt = 1.f/flt(count);
int i;
lprintf("Simple Euler:\n");
for (i=0;
i < count; i++) {
x += v*dt;
v += a*dt;
lprintf("V: %f X: %f\n",v,x);
} // for
x = v = 0.f;
lprintf("More Accurate Euler:\n");
for (i=0;
i < count; i++) {
x += v*dt + .5f*a*dt*dt;
v += a*dt;
lprintf("V: %f X: %f\n",v,x);
} // for
x = v = 0.f;
lprintf("NSV: Newton-Stormer-Verlet
('Semi-implicit' Euler):\n");
for (i=0;
i < count; i++) {
v += a*dt;
x += v*dt;
lprintf("V: %f X: %f\n",v,x);
} // for
lprintf("Velocity-less Verlet:\n");
flt xc = 0.f;
flt xo = xc;
for (i=0;
i < count; i++) {
v = xc - xo + a*dt*dt;
xo = xc;
xc += v;
lprintf("V: %f X: %f\n",v/dt,xc);
} // for
lprintf("My NSV variant\n");
x = v = 0.f;
flt dt2 = dt*dt;
for (i=0;
i < count; i++) {
v += a*dt2;
x += v; // v is
prescaled: really a displacement.
lprintf("V: %f X: %f\n",v/dt,x);
} // for
lprintf("Velocity Verlet variant\n");
flt oldAccel = a; //0.f; // May be
self-starting issues...
x = v = 0.f;
for (i=0;
i < count; i++) {
x += v*dt + .5f*oldAccel*dt*dt;
v +=
.5f*(oldAccel+a)*dt;
oldAccel = a;
lprintf("V: %f X: %f\n",v,x);
} // for
}
} test;
Output:
Simple Euler: V: 9.810000 X: 4.855951 More Accurate Euler: V: 9.810000 X: 4.905001 NSV: Newton-Stormer-Verlet ('Semi-implicit' Euler): V: 9.810000 X: 4.954051 Velocity-less Verlet: V: 9.810033 X: 4.954060 My NSV variant V: 9.810008 X: 4.954050 Velocity Verlet variant V: 9.810000 X: 4.905001
Use the following equation to work through your examples
by hand to verify your work:
x' = x + v*dt + .5*a*dt*dt
For example, with dt = 1sec, initial velocity 0, an object will drop .5*9.81m or
4.905m.
v' = v + a*dt, after 1s, v = 9.81m/s.
Exact answers for the above Test program
from direct computation:
I don't see a reason to ever use Explicit Euler. The
simplified (non-staggered-grid) Newton-Störmer-Verlet (NSV) form has the same
computational cost with much greater stability and accuracy. The slightly
more complex Velocity Verlet variant is even more accurate. See also my short
test code example that shows the accuracy of Explicit Euler, Verlet, Velocity
Verlet, and another simple variant.
Standard Runge-Kutta (2,4) is not symplectic, loses energy, and is less
accurate than symplectic methods. 2,4 derivs per step also make them more
expensive and complicated.
See the discussion, graphics, and proofs here. From Hairer et al's book,
Geometric Numerical Integration, page 4, fig. 1.2, excellent graphics show
that:
1. Explicit Euler gains energy: is unstable (effectively underdamped).
2. Implicit Euler loses energy: is overstable (effectively overdamped).
3. Symplectic Euler preserves energy (technically, "geometric structure/area of
the map/graph"): it is the most accurate of the three methods
(Velocity-less Verlet, NSV, etc., are slightly different, but have similar
symplectic behavior. See the analysis, proofs, and graphs in the book).
See also page 9: Symplectic Euler provides the most accurate solution compared
to the exact solution (Kepler problem). Explicit Euler drifts away from the
solution, and the Implicit Midpoint drifts much more than Symplectic Euler. See
page 10 for error (accuracy) comparisions:
method,error in H,error in L,global error
explict Euler,O(th),O(th),O(th^2)
symplectic Euler,O(h),0,O(th)
implicit midpoint,O(h^2),0,O(th^2)
For game simulations, I have implemented, tested, and studied Explicit Euler,
Velocity-less Verlet, NSV, Velocity Verlet, Trapezoid NSV variants, RK2, and
RK4. For my game simulations, NSV is the clear winner for stability, efficiency,
and ease of implementation (Trapezoid/RK-like tricks must be used to improve
high-angular velocity accuracy/behavior for rigid bodies).
The elegance of the Runge-Kutta Taylor series improving the state
curve-approximation did not provide any significant value (the only area it
helped was to provide more stability for spring-like forces and periodic motion
(IIRC, running a single deriv-per-step method at 1/2,1/4 the time step had about
the same effect as RK2,RK4)). Symplectic RK variants exist, but for games, I am
not aware of any method that (significantly) improves on (simplified) NSV:
v += a*dt
x += v*dt
Once you have integrated position, you can map/scale the
position to the correct coordinates. Thus, scale x to convert to cm/pixels, etc.
for rendering.
For pixel applications, you'll need to compute the motion in floating point or
fixed point; probably easiest to compute as float, then convert to int for pixel
rendering (if lots of pixels, convert using fistp in asm as opposed to int px =
int(fx)):
inline int fastFtoI(float f) { int i; __asm { fld f frndint fistp i } return i; } // Method from: http://ldesoras.free.fr/doc/articles/rounding_en.pdf
// fastRoundToInt(f + .5) will return f+1 via fadd st,st(0) and sar i,1 trick. inline int fastRoundToInt(float f) { int i; static const float round_to_nearest = 0.5f; __asm { fld f fadd st,st(0) fadd round_to_nearest fistp i sar i,1 } return i; }