From 4b4c77c4204a87d96ed653903afe937d490e878b Mon Sep 17 00:00:00 2001
From: Robert Osfield <robert@openscenegraph.com>
Date: Fri, 24 Feb 2006 13:56:07 +0000
Subject: [PATCH] Added clamping of the wind force contribution to acceletion
 of particle so that the de-acceleration never exceed the actual wind vector
 itself.

---
 src/osgParticle/FluidProgram.cpp | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/src/osgParticle/FluidProgram.cpp b/src/osgParticle/FluidProgram.cpp
index b3767f7ca..5465360ab 100644
--- a/src/osgParticle/FluidProgram.cpp
+++ b/src/osgParticle/FluidProgram.cpp
@@ -32,16 +32,25 @@ void osgParticle::FluidProgram::execute(double dt)
             float Volume = Area*radius*four_over_three;
         
             // compute force due to gravity + boyancy of displacing the fluid that the particle is emersed in.
-            osg::Vec3 force = _acceleration * (particle->getMass() - _density*Volume);
-
+            osg::Vec3 accel_gravity = _acceleration * ((particle->getMass() - _density*Volume) * particle->getMassInv());
+            
             // compute force due to friction
+            osg::Vec3 velBefore = particle->getVelocity();
             osg::Vec3 relative_wind = particle->getVelocity()-_wind;            
-            force -= relative_wind * Area * (_viscosityCoefficient + _densityCoefficeint*relative_wind.length());            
+            osg::Vec3 wind_force = - relative_wind * Area * (_viscosityCoefficient + _densityCoefficeint*relative_wind.length());
+            osg::Vec3 wind_accel = wind_force * particle->getMassInv();
+
+            double critical_dt2 = relative_wind.length2()/wind_accel.length2();
+            double compenstated_dt = dt;
+            if (critical_dt2 < dt*dt)
+            {
+                // osg::notify(osg::NOTICE)<<"** Could be critical: dt="<<dt<<" critical_dt="<<sqrtf(critical_dt2)<<std::endl;
+                compenstated_dt = sqrtf(critical_dt2)*0.8f;
+            }
+
+            particle->addVelocity(accel_gravity*dt + wind_accel*compenstated_dt);
             
-            // divide force by mass to get acceleration.
-            force *= particle->getMassInv()*dt;
-            
-            particle->addVelocity(force);
+
         }
     }
 }