WindStorm.h

00001 #include <stdlib.h>
00002 #include <math.h>
00003 #include <windows.h>
00004 #include "debug.h"
00005 
00006 #define TRUE_WIND_VELOCITY
00007 
00008 extern Debug* dbg;
00009 
00010 void init();
00011 void deInit();
00012 void InitParticles();
00013 void QuitApplication();
00014 float WindStormStrength(float fAltitude);
00015 void WindVelocity(float *pfLocation, float *pfVelocity);
00016 void UpdatePhysics();
00017 void PhysicsStep();
00018 
00019 class Particle
00020 {
00021         public:
00022                 Particle();
00023 
00024         public:
00025                 float m_pfPosition[3];  
00026                 float m_pfVelocity[3];  
00027                 float m_fMass;                  
00028                 float m_fHalfRhoSref;   
00029 };
00030 
00031 float g_fPhysicsTimeStep = 0.02f;                                       // lock physics to 50 updates/second
00032 float g_fAccumulatedTimeStep = 0.0f;                            // accumulated time since last physics update
00033 float g_fSimTime = 0.0f;                                                        // current simulation time
00034 float g_OneOverTwoPi = 1.0f/(8.0f * atanf(1.0f));       // inverse of 2*pi.
00035 float g_pfStormCenter[] = {0.0f, 0.0f};                         // 2D center of storm at ground level.
00036 float g_fFluidDensity = 1.225;                                          // Air density around sea level, kg/meters cubed
00037 float g_fPI = 4.0f * atanf(1.0f);                                       // PI
00038 float g_fGravitationalAccel = -9.8f;                            // rate of gravitational acceleration, in Z direction, meters/second-squared
00039 float g_fRandMax(float(RAND_MAX));                                      // maximum value returned by rand()
00040 float g_fS0 = 0.5f;                                                                     // particle strength at altitude = 0 meters
00041 float g_fS500 = 50.0f;                                                          // particle strength at altitude = 500 meters
00042 LARGE_INTEGER g_liTimerFrequency;                                       // Timer frequency
00043 LARGE_INTEGER g_liPreviousCounter;                                      // Value of performance counter at last update.
00044 
00045 const unsigned int g_uiGoalNumParticles = 400;      // desired number of particles
00046 unsigned int g_uiNumParticles = 0;                                      // Number of particles being simulated
00047 Particle *g_pParticles = 0;                                                     // Array of particles for a *very* simple particle system
00048 
00049 void init()
00050 {
00051         if (FALSE == QueryPerformanceFrequency(&g_liTimerFrequency))
00052         {
00053             dbg->log("FALSE == QueryPerformanceFrequency(&g_liTimerFrequency)");
00054                 return;
00055         }
00056 
00057         InitParticles();
00058 
00059         if (FALSE == QueryPerformanceCounter(&g_liPreviousCounter))
00060         {
00061             dbg->log("FALSE == QueryPerformanceCounter(&g_liPreviousCounter)");
00062                 return;
00063         }
00064 }
00065 
00066 Particle::Particle() : m_pfPosition(), m_pfVelocity(),
00067         m_fMass(1.0f), m_fHalfRhoSref(1.0f)
00068 {
00069         m_pfPosition[0] = 0.0f;
00070         m_pfPosition[1] = 0.0f;
00071         m_pfPosition[2] = 0.0f;
00072 
00073         m_pfVelocity[0] = 0.0f;
00074         m_pfVelocity[1] = 0.0f;
00075         m_pfVelocity[2] = 0.0f;
00076 }
00077 
00078 void InitParticles()
00079 {
00080         g_pParticles = new Particle[g_uiGoalNumParticles];
00081         float fRadius, Sref;
00082         if (0 != g_pParticles)
00083         {
00084                 g_uiNumParticles = g_uiGoalNumParticles;
00085                 for (unsigned int i = 0; i < g_uiGoalNumParticles; i++)
00086                 {
00087                         g_pParticles[i].m_pfPosition[0] = -2.5f + 5.0f * float(rand()) / g_fRandMax;
00088                         g_pParticles[i].m_pfPosition[1] = -2.5f + 5.0f * float(rand()) / g_fRandMax;
00089                         g_pParticles[i].m_pfPosition[2] = 500.0f * float(rand()) / g_fRandMax;
00090 
00091                         g_pParticles[i].m_fMass = 0.05f + .05f * float(rand()) / g_fRandMax;
00092                         fRadius = 0.2f + 0.2f*float(rand()) / g_fRandMax;
00093                         Sref = g_fPI * fRadius * fRadius;
00094                         g_pParticles[i].m_fHalfRhoSref = 0.5 * Sref * g_fFluidDensity;
00095                 }
00096         }
00097 }
00098 
00099 void deInit()
00100 {
00101         delete [] g_pParticles;
00102 }
00103 
00104 // calculate the storm strength at the sent altitude
00105 float WindStormStrength(float fAltitude)
00106 {
00107         return(g_fS0 + ((g_fS500 - g_fS0) * fAltitude * fAltitude / 250000.0f));
00108 }
00109 
00110 // calculate the wind velocity due to the storm, at the sent
00111 // location
00112 void WindVelocity(float *pfLocation, float *pfVelocity)
00113 {
00114         if (0 == pfLocation || 0 == pfVelocity)
00115                 return;
00116 
00117         float fDX = pfLocation[0] - g_pfStormCenter[0];
00118         float fDY = pfLocation[1] - g_pfStormCenter[1];
00119 
00120         float fR = sqrtf(fDX * fDX + fDY * fDY);
00121         if (fR < 1.e-4)
00122                 fR = 1.0f;
00123 
00124 // for a true potential vortex, the velocity approaches
00125 // infinity at the core. We need to limit this to avoid
00126 // undesirable effects. We limit it here by artificially
00127 // forcing particles inside a certain radius to behave
00128 // as though they were at that limit radius. Another
00129 // approach would be to have the vortex strength fade
00130 // towards zero at the core.
00131         float fOneOverR = (1.0f / fR);
00132         if (fOneOverR < 1.0f)
00133         {
00134                 fOneOverR = 1.0f;
00135         }
00136 
00137 // TRUE_WIND_VELOCITY calculates the actual velocity at a point due to
00138 // a potential vortex. If this is NOT defined, a numerical approximation
00139 // to the velocity is calculated instead, based on the physics time step,
00140 // that is designed to keep a particle with zero mass on an exact circular
00141 // path (to floating point precision). The latter case can produce better
00142 // results, since it prevents particle drift due to anything other than
00143 // inertia. The former case, using the true potential vortex velocity,
00144 // will result in particles that drift away from a circular path just
00145 // due to truncation error in the time stepped solution.
00146 //#define TRUE_WIND_VELOCITY
00147 #ifdef TRUE_WIND_VELOCITY
00148         float fCommonFactor = WindStormStrength(pfLocation[2]) * fOneOverR * g_OneOverTwoPi;
00149 
00150         pfVelocity[0] =  pfLocation[1] * fCommonFactor;
00151         pfVelocity[1] = -pfLocation[0] * fCommonFactor;
00152 #else
00153         float fTangentialVelocity = WindStormStrength(pfLocation[2]) * fOneOverR * g_OneOverTwoPi;
00154         float fDeltaAngle = fTangentialVelocity * g_fPhysicsTimeStep;
00155         float fAngle = atan2(fDY, fDX) + fDeltaAngle;
00156         float fNewX = fR * cos(fAngle);
00157         float fNewY = fR * sin(fAngle);
00158         pfVelocity[0] = (fNewX - pfLocation[0]) / g_fPhysicsTimeStep;
00159         pfVelocity[1] = (fNewY - pfLocation[1]) / g_fPhysicsTimeStep;
00160 #endif // TRUE_WIND_VELOCITY
00161 
00162 // there is no vertical wind velocity due to the storm in this example. You
00163 // could add it to get interesting effects, though.
00164         pfVelocity[2] = 0.0f;
00165 }
00166 
00167 // Update the physics as necessary. This uses a simple
00168 // frame-rate independent technique!
00169 void UpdatePhysics()
00170 {
00171         LARGE_INTEGER PC;
00172 
00173         if (FALSE == QueryPerformanceCounter(&PC))
00174         {
00175             dbg->log("FALSE == QueryPerformanceCounter(&PC)");
00176                 return;
00177         }
00178 
00179         LONGLONG DeltaCounter = PC.QuadPart - g_liPreviousCounter.QuadPart;
00180 
00181         g_liPreviousCounter = PC;
00182 
00183         float fDeltaTime = float(DeltaCounter) / float(g_liTimerFrequency.QuadPart);
00184 
00185         g_fAccumulatedTimeStep += fDeltaTime;
00186 
00187         while (g_fAccumulatedTimeStep > g_fPhysicsTimeStep)
00188         {
00189                 PhysicsStep();
00190                 g_fAccumulatedTimeStep -= g_fPhysicsTimeStep;
00191         }
00192 }
00193 
00194 // Take one physics step. This is where the actual physics integration
00195 // step is performed. The time step is always fixed at g_fPhysicsTimeStep
00196 void PhysicsStep()
00197 {
00198         float fVwind[3], fForce[3], fDragMagnitude, fVrelwind_squared, fSpeed;
00199         for (unsigned int i = 0; i < g_uiNumParticles; i++)
00200         {
00201                 WindVelocity(g_pParticles[i].m_pfPosition, fVwind);
00202 
00203 // DO_TEST treats the particles as though they were massless, in which
00204 // case they simply take on the velocity of the windfield instead of
00205 // reacting physically to applied forces. This is useful to check and
00206 // ensure that the simulation does not exhibit significant numerical
00207 // drift (error).
00208 #ifndef MASSLESS_PARTICLES
00209                 // turn fVwind into relative wind and compute its
00210                 // magnitude squared
00211                 fVwind[0] -= g_pParticles[i].m_pfVelocity[0];
00212                 fVwind[1] -= g_pParticles[i].m_pfVelocity[1];
00213                 fVwind[2] -= g_pParticles[i].m_pfVelocity[2];
00214                 fVrelwind_squared = fVwind[0] * fVwind[0] + fVwind[1] * fVwind[1] + fVwind[2] * fVwind[2];
00215                 fSpeed = sqrtf(fVrelwind_squared);
00216 
00217                 // calculate the drag force. Assume Reynold's Number
00218                 // is between 2000 and 200,000, and, therefore,
00219                 // drag coefficient is roughly a constant at 0.4.
00220                 const float fCD = 0.4;
00221 
00222                 fDragMagnitude = g_pParticles[i].m_fHalfRhoSref * fVrelwind_squared * fCD;
00223 
00224                 // fVwind[*] / fSpeed gives the direction of the drag force
00225                 fForce[0] = fDragMagnitude * fVwind[0] / fSpeed;
00226                 fForce[1] = fDragMagnitude * fVwind[1] / fSpeed;
00227                 fForce[2] = fDragMagnitude * fVwind[2] / fSpeed;
00228 
00229                 // increment the force to add in the weight of the particle.
00230 #ifdef INCLUDE_GRAVITY
00231                 fForce[2] += g_fGravitationalAccel * g_pParticles[i].m_fMass;
00232 #endif // INCLUDE_GRAVITY
00233 
00234                 // increment the velocity based on acceleration due to applied
00235                 // forces. This is a simple Euler integration, which is easy
00236                 // but certainly NOT the best choice for a robust game physics
00237                 // implementation!
00238                 g_pParticles[i].m_pfVelocity[0] += g_fPhysicsTimeStep * fForce[0] / g_pParticles[i].m_fMass;
00239                 g_pParticles[i].m_pfVelocity[1] += g_fPhysicsTimeStep * fForce[1] / g_pParticles[i].m_fMass;
00240                 g_pParticles[i].m_pfVelocity[2] += g_fPhysicsTimeStep * fForce[2] / g_pParticles[i].m_fMass;
00241 
00242 
00243 #else // massless particle option
00244                 g_pParticles[i].m_pfVelocity[0] = fVwind[0];
00245                 g_pParticles[i].m_pfVelocity[1] = fVwind[1];
00246                 g_pParticles[i].m_pfVelocity[2] = fVwind[2];
00247 #endif // DO_TEST
00248 
00249                 // Final step is to update position based on current velocity.
00250                 // Since we are using the updated velocity, this is not simple
00251                 // Euler integration, since it uses the updated (t + dt) velocity.\
00252                 // But, still conditionally stable.
00253                 g_pParticles[i].m_pfPosition[0] += g_fPhysicsTimeStep * g_pParticles[i].m_pfVelocity[0];
00254                 g_pParticles[i].m_pfPosition[1] += g_fPhysicsTimeStep * g_pParticles[i].m_pfVelocity[1];
00255                 g_pParticles[i].m_pfPosition[2] += g_fPhysicsTimeStep * g_pParticles[i].m_pfVelocity[2];
00256 
00257                 // respawn the particle if it falls to the ground (due to gravity)
00258                 if (g_pParticles[i].m_pfPosition[2] <= 0.0f)
00259                 {
00260                         g_pParticles[i].m_pfPosition[0] = -2.5f + 5.0f * float(rand()) / g_fRandMax;
00261                         g_pParticles[i].m_pfPosition[1] = -2.5f + 5.0f * float(rand()) / g_fRandMax;
00262                         g_pParticles[i].m_pfPosition[2] = 500.0f * float(rand()) / g_fRandMax;
00263                 }
00264         }
00265 }

Generated on Sun Dec 2 17:09:58 2007 for Swarm by  doxygen 1.4.6-NO