Sorry to the usual audience of this blog, but this article aims to document a very specific topic: Pacejka and ODE. Since I was unable to find any resource on the web, I hope that, thanks to search engines, it may help people trying to mix these two things together.

While working on preliminary tests for ManiaCrash, I had to understand and implement Pacejka Magic Formula, the famous tire model used to simulate quite accurately tires behavior. You may have a look at what is probably one of the best explaination and vulgarization about Pacejka, on Racer website. If it’s quite easy to find some documentation about Pacejka equations, I had hard times finding tips about how implementing it, especially in my particular case, where I need to « integrate » Pacejka into the physics engine I use in Raydium: ODE.

Doing a such thing implies a few points: you must « extract » inputs for Magic Formula from ODE, apply Pacejka on these inputs, and « send » the result to ODE. And each of these steps reveals to be an issue.

While the point of this article is not to explain Pacejka basics, let’s have a quick look at the main two forces, longitudinal and lateral, since the usual last one (aligning moment) is then very easy to understand.

• longitudinal formula: computes Fx force the tire will apply on the road, « along » the car. Inputs: slip ratio (SR) and load (Fz)
• lateral formula: compute Fy force the tire will apply sideways on the road (giving, more or less, the G force the car will support during the corner). Inputs: slip angle (SA) and load (Fz)

Obviously, the best place for all this is the dCollide() result loop, for each tire/road contact. If you’re using Raydium, see raydium_ode_CollideCallback.

Inputs

As a reminder, slip ratio gives how much the tire is « sliding » on the road along its X (It’s related to accelerating and braking situations, in other words) and is often described as the ratio between current angular velocity of the tire and free rolling speed, where I prefer to explain it as the the relative speed of the tire at road contact point. On the other hand, slip angle is the angle between wheel heading (local X axis) and its moving direction. The last input, Fz, is the vertical load of the tire.

The challenge is then to find the best way to get these informations from ODE, and here’s what I’ve come up with. Please note I’ve removed unnecessary variable declarations so the code is easier to read.

First, let’s try to determine the slip angle. The basic idea is to get the velocity of the wheel, turn it into a vector, get the wheel heading and then apply a scalar product on these two.

//// Slip Angle

dVector3 wh, wt; // wheel heading, wheel travel
dReal sa; // slip angle

// 1 - get wheel travel vector
// velocity of the wheel at (0,0,0)
dBodyGetRelPointVel(wheel.body,0,0,0,wt);
// same as wt=dBodyGetLinearVel(wheel.body);

// 2 - get wheel rotation (Y) axis (so we get wheel heading + 90, see step 4)
dBodyVectorToWorld(raydium_ode_element[e1].body,0,1,0,wh);

// we normalize wt (wh is already normalized)
vnorm(wt);

// 3 - determine angle between these two vectors, using scalar product.
// We need a result in the first quadrant only: SA>90 means nothing, since
// a tyre does not have forward/backward, but only a longitudinal axis.
scal=(wt[0]*wh[0]) + (wt[1]*wh[1]) + (wt[2]*wh[2]);
sa=180*(acos(scal))/PI;

// 4 - here, we "switch" from wheel Y to X, and limit SA to [0,90] range
sa=abs(sa-90);

Then comes the slip ratio. Using the « relative speed at contact point » idea and with the magic of ODE, the result is surprisingly simple !

//// Slip Ratio
dReal sr; // slip ratio

// 1 - relative speed at contact point (res is in world coords)
// "n" is a pointer to the dContact structure returned by ODE for the current contact
dBodyGetPosRelPoint (wheel.body,n->geom.pos[0],n->geom.pos[1],n->geom.pos[2],pos);
dBodyGetRelPointVel (wheel.body,pos[0],pos[1],pos[2],res);

// 2 - get velocity along wheel heading
// done using dot product : (fdir1 . relative speed)
sr=raydium_math_abs(
(n->fdir1[0]*res[0]) +
(n->fdir1[1]*res[1]) +
(n->fdir1[2]*res[2]) );

// note: fdir1 is the contact direction. See below.

The last thing is to get the tire load. We use some sort of temporal cohesion and get the last ODE step result as an input. It’s done thru ODE joint feeback, and if here I use Raydium functions to do so, it’s very easy to code your own system using dJointSetFeedback(). Of course, the result should probably be scaled before sending it to Pacejka, since you may not use real word units in your application/game.

//// Tire Load
dJointFeedback *jf;
dVector3 force; // local force produced by the tyre during last physics step

// let's find the generated "tire force" during last step
raydium_ode_contact_feedback_save(wheel); // we ask Raydium to store joint feedback during the next step ...
jf=raydium_ode_contact_feedback_get(wheel); // ... and we fetch the saved feedback data from the previous step.
// note: again, fdir1 is the contact direction. See below.
dRFrom2Axes(wheel_matrix,n->fdir1[0],n->fdir1[1],n->fdir1[2],wh[0],wh[1],wh[2]);
dMULTIPLY1_331(force,wheel_matrix,jf->f1); // "world to wheel"
// force[2] is the vertical load

Outputs

Once you’ve computed Fx and Fy using your Pacejka code (again, you’ll find a lot of documentations and sample implementations for this all around Internet), you need to send it back to ODE. As a simplification, Pacejka formula outputs a « maximum » force that the tire can support, but the situation may require way less force than that. And if it requires more, we must make the tire sliding. How can we explain this to ODE since the usual friction pyramid model does not allow this ?

This part was quite tricky to me, since I discovered that I was not that good at understanding the whole ODE contact model. In facts, the answer was somewhere in the ODE userguide, where it’s said that when dContactApprox « is not specified then the constant-force-limit approximation is used (and mu is a force limit). » Great, that’s just what we need: constant force limit model !

// "n" is a pointer to the dContact ODE structure for the current contact
n->surface.mode = dContactSoftERP | dContactSoftCFM | dContactMu2;
// note that there's no dContactMu1 in ODE (it's implied by mu2)

Since we use mu1 and mu2, another important thing then is to determine contact direction. This is also needed by previous code, where « fdir1 » was used. In facts, this is probably one of the first things to do in your contact callback.

//// contact direction (n->fdir1), based on wheel heading

// cross product of normal and wheelY (fdir2)
n->surface.mode |= dContactFDir1;
n->fdir1[0] = wh[1]*n->geom.normal[2] - wh[2]*n->geom.normal[1];
n->fdir1[1] = wh[2]*n->geom.normal[0] - wh[0]*n->geom.normal[2];
n->fdir1[2] = wh[0]*n->geom.normal[1] - wh[1]*n->geom.normal[0];

// wh is wheel heading dVector3 + 90 degrees. See slip angle, above

Then, you just have to set mu1 and mu2. Again, you may have to scale these values to suit your world units.

mu1=fx*PACEJKA_MC_SCALE;
mu2=fy*PACEJKA_MC_SCALE;

if(isnan(mu1) || mu1<0) mu1=0;
if(isnan(mu2) || mu2<0) mu2=0;

n->surface.mu=mu1; n->surface.mu2=mu2;

That’s it !

This article is a first try, comments are welcome. A few topics are left unexplained, as Pacejka computations and the need for a combined model, for instance. I’ll get on this later. Currently, ManiaCrash source code is not public yet (until the very first alpha), but I’ll add a link to suitable files once done.