[ODE] Atmospheric drag

[Oktay Ahiska]

Howdy folks,

I am trying to simulate atmospheric drag with the following code:
*pIter points to an internal structure with some info about the
the current object, including its current velocity and surface
area (these are all spheres) and its dBbody and dMass data.

	// initielze with negative of current normalized velocity vector
	dVector3 dragForce;
	dragForce[0] = -(*pIter)->velNorm[0];
	dragForce[1] = -(*pIter)->velNorm[1];
	dragForce[2] = -(*pIter)->velNorm[2];

         double dragAmt;
         double vm = (*pIter)->velMag;   // pre-computed velocity magnitude

         // surfaceArea: pre-computed surface area of the sphere
	// atmosphericDrag: drag constant, passed to us from caller
	
         dragAmt = vm * vm * (*pIter)->surfaceArea * atmosphericDrag;

         if (dragAmt > 0.) {
	      // clamp force if excessive
	      double dragMax = vm*(*pIter)->mass.mass;
               if (dragAmt > dragMax) {
                     dragAmt = dragMax;
               }

               dBodyAddForce( (*pIter)->body,
			      dragAmt * dragForce[0],
			      dragAmt * dragForce[1],
			      dragAmt * dragForce[2]);
          }

So... I know dragAmt is not really accurate, but I can tweak atmosphericDrag
till it looks "right" but is this otherwise ok?  Am I right to cap off my
drag at  vm*(*pIter)->mass.mass?

Any comments will be very much appreciated, including "you can already do that
by..." or some such.

Thanks in advance,
Oktay