flockanimatortornado.cpp

#include "debug.h"
#include "flock.h"
#include "boid.h"
#include "flockanimatortornado.h"

using namespace irr;
using namespace core;
using namespace scene;

extern Debug* dbg;

FlockAnimatorTornado::FlockAnimatorTornado(IrrlichtDevice* d)
{
    // ctor
    device = d;
    pi = 3.1415926535897932384626433832795f;
    oneOverTwoPi = 1.0f / (2.0f * pi);
    fluidDensity = 1.225f;   
    gravitationalAccel = -9.8f;
    strength_0 = 0.5f;
    strength_500 = 250.0f;
    stormCenter = vector3df(0.0f, 0.0f, 0.0f);
        mass = 0.01f;                                       
        halfRhoSref = pi * 3 * 3 * 0.5f * fluidDensity;     

    dbg->log("create FlockAnimatorTornado");
}

FlockAnimatorTornado::~FlockAnimatorTornado()
{
    // dtor
    dbg->log("FlockAnimatorTornado deleted");
}

void FlockAnimatorTornado::animateNode(ISceneNode* node, u32 timeMs)
{
    // for each Boid in Flock
    Flock* f = ((Flock*) node);
    stormCenter = f->getTarget();
    list<Boid*>::Iterator it = f->getBoidList().begin();
    for (; it != f->getBoidList().end(); ++it)
    {
        vector3df force = getForce((*it), timeMs);
        force *= f->getTornadoWeight();
        (*it)->velocity += force;

        // debug
        if (it == f->getBoidList().begin())
        {
            dbg->show("first boid pos", (*it)->position);
            dbg->show("first boid force", force);
        }
    }

}

// calculate the storm strength at the sent altitude
f32 FlockAnimatorTornado::getWindStrength(f32 altitude)
{
        return (strength_0 + ((strength_500 - strength_0) * altitude * altitude / 250000.0f));
}

// calculate the wind velocity due to the storm, at the sent location
vector3df FlockAnimatorTornado::getWindVelocity(vector3df location, u32 timeMs)
{
    // return the resulting velosity;
    vector3df result;

        f32 dX = location.X - stormCenter.X;
        f32 dZ = location.Z - stormCenter.Z;

        f32 r = sqrtf(dX * dX + dZ * dZ);
        if (r < 1.e-4)
        {
                r = 1.0f;
        }

    // for a true potential vortex, the velocity approaches
    // infinity at the core. We need to limit this to avoid
    // undesirable effects. We limit it here by artificially
    // forcing particles inside a certain radius to behave
    // as though they were at that limit radius. Another
    // approach would be to have the vortex strength fade
    // towards zero at the core.
        f32 oneOverR = 1.0f / r;
        if (oneOverR > 1.0f)
        {
                oneOverR = 1.0f;
        }

    // TRUE_WIND_VELOCITY calculates the actual velocity at a point due to
    // a potential vortex. If this is NOT defined, a numerical approximation
    // to the velocity is calculated instead, based on the physics time step,
    // that is designed to keep a particle with zero mass on an exact circular
    // path (to floating point precision). The latter case can produce better
    // results, since it prevents particle drift due to anything other than
    // inertia. The former case, using the true potential vortex velocity,
    // will result in particles that drift away from a circular path just
    // due to truncation error in the time stepped solution.
    //#define TRUE_WIND_VELOCITY
#define TRUE_WIND_VELOCITY
#ifdef TRUE_WIND_VELOCITY
        f32 commonFactor = getWindStrength(location.Z) * oneOverR * oneOverTwoPi;

        result.X = dZ * commonFactor;          // original: location.Z * commonFactor
        result.Z = -dX * commonFactor;         // original: location.X * commonFactor
        vector3df dbgCheck=result;              

#else
    f32 physicsTimeStep = timeMs / 1000.f;
        f32 tangentialVelocity = getWindStrength(location.Y) * oneOverR * oneOverTwoPi;
        f32 deltaAngle = tangentialVelocity * physicsTimeStep;
        f32 angle = atan2(dZ, dX) + deltaAngle;
        f32 newX = r * cos(angle);
        f32 newZ = r * sin(angle);
        result.X = (newX - location.X) / physicsTimeStep;
        result.Z = (newZ - location.Z) / physicsTimeStep;
#endif // TRUE_WIND_VELOCITY

    // there is no vertical wind velocity due to the storm. We could add it to
    // get interesting effects, though.
        result.Y = 0.0f;

        return result;
}

vector3df FlockAnimatorTornado::getForce(Boid* b, u32 timeMs)
{
        vector3df result;
        vector3df wind;
        vector3df force;
        f32 dragMagnitude;
        f32 relwind_squared;
        f32 speed;

        wind = getWindVelocity(b->position, timeMs);

        // turn wind into relative wind and compute its magnitude squared
        wind -= b->velocity;

        relwind_squared = wind.X * wind.X + wind.Y * wind.Y + wind.Z * wind.Z;
    speed = sqrtf(relwind_squared);

    // calculate the drag force. Assume Reynold's Number
        // is between 2000 and 200,000, and, therefore,
        // drag coefficient is roughly a constant at 0.4.
        const float cd = 0.4f;

    dragMagnitude = halfRhoSref * relwind_squared * cd;

    // wind / speed gives the direction of the drag force
    force = wind * dragMagnitude / speed;

        // increment the force to add in the weight of the particle.
        force.Y += gravitationalAccel * mass;

    // increment the velocity based on acceleration due to applied
        // forces. This is a simple Euler integration, which is easy
        // but certainly NOT the best choice for a robust game physics
        // implementation!
        result = force / mass;

    return result;
}

---