CODE

struct SpringNode {
CVector3<double> position;           // My current position
CVector3<double> velocity;           // My velocity

CVector3<double> force;

double mass;                                       //
std::vector<SpringNode*> neighbors;               // Neighbors
std::vector<double> neighborDistances;             // How far am I away from my neighbors
};

struct AnchorNode {
SpringNode *node;
CVector3<double> position;           // Anchor position
};


class ISpringSystem {

public:

/*!
 */
void ComputeSingleForce( SpringNode const *a, SpringNode const *b );

/*!
 */
void AddAnchorNode( int node, CVector3<double> const &anchorPosition );

/*!
 */
void Animate();


/*!
 */
void Draw();



/*  The derived class will connect the springs for us.
 */
virtual void ConnectSprings( )=0;

protected:
double _springK;     // Elasticity coefficient. The larger K is, the stiffer the spring.
double _inertialDistance;                            // Inertial or resting distance of the springs
CVector3<double> _force;    // Temporary force vector for ComputeSingleForce
CVector3<double> _vel;     // Temporary velocity vector for ComputeSingleForce
std::vector<SpringNode> _springNodes;   // Vector of spring nodes
std::vector<AnchorNode> _anchorNodes;   // Vector of anchor nodes
std::vector<SpringNode*> _movingNodes;   // Vector of moving nodes
double _mass;      // Mass of each spring node
double _energyLoss;     // Energy loss to occur as the simulation progresses
CTimer *_timer;      // Timer used in physics calculations

};


#define SPRING_TOLERANCE           (0.0000000005)

void ISpringSystem::ComputeSingleForce( SpringNode const *a, SpringNode const *b )
{
// Calculate vector between b and a (giving the direction and magnitude of the vector they make up)
_force = b->position - a->position;

// Calculate the relative velocity between a and b
_vel = b->velocity - a->velocity;

//
if ( _force.Length() - _inertialDistance < SPRING_TOLERANCE )
{
 _force.Zero();
 return;
}

// Calculate the intensity of the spring ( -k * deltaDistance )
// This is where the spring constant comes into play. A higher K will yield
// a stiffer spring.
double intensity = -_springK * (_force.Length() - _inertialDistance);

// Normalize the force direction
_force.Normalize();

// Apply the internal force and the friction force
//
_force += ((_force/_force.Length())*intensity) - (40000 * _vel);
}

void ISpringSystem::AddAnchorNode( int node, CVector3<double> const &anchorPosition )
{
// Setup the anchor nodes
_anchorNodes.resize(_anchorNodes.size()+1);
_anchorNodes[_anchorNodes.size()-1].node = &_springNodes[node];
_anchorNodes[_anchorNodes.size()-1].position = anchorPosition;
_movingNodes.erase( _movingNodes.begin() + node );
}



void ISpringSystem::Animate()
{
static bIsInit=false;


// Initialize the timer
if ( !bIsInit ) _timer->Init();


// Get the elapsed seconds
double elapsed = _timer->GetElapsedSeconds();


// Gravitational force
// To compute the gravitational force we use F = ma
// (m = mass, a =  acceleration due to gravity (in units per second which is g=9.81m/s^2)
double fg = -9.81 * _mass;

// Compute all external forces on the springs
for( int i=0; i<_movingNodes.size(); ++i )
{
 // Zero out the force vector
 _movingNodes[i]->force.Zero();

 // Factor in the gravitational force. Since it moves only in the y direction, we
 // only modify the y component of the vector.
 _movingNodes[i]->force[1] = fg;
}


// Compute internal forces on the springs
for( int i=0; i<_movingNodes.size(); ++i )
{
 // Compute the force the neighbor nodes are exerting on each node
 for( int j=0; j<_movingNodes[i]->neighbors.size(); ++j )
 {
  //
  ComputeSingleForce( _movingNodes[i]->neighbors[j], _movingNodes[i] );
  _movingNodes[i]->force += _force;                       // Add the force exerted to this nodes force vector
 }
}


// Increment velocity and position of nodes
for( int i=0; i<_movingNodes.size(); ++i )
{
 _movingNodes[i]->velocity += (_movingNodes[i]->force/_mass)*elapsed;      // Add
 _movingNodes[i]->position += _movingNodes[i]->velocity;       // Add the velocity to pos
 _movingNodes[i]->velocity *= _energyLoss;
}

}



void ISpringSystem::Draw()
{
// Push the current attribute list onto the OpenGL attribute stack
glPushAttrib(GL_CURRENT_BIT);
glLoadIdentity();


// Draw as lines
glBegin(GL_LINES);
glDisable(GL_CULL_FACE);

glColor3f(1.0f, 0.0f, 0.0f);

for( int i=0; i<_springNodes.size()-1; ++i )
{
 glVertex3f( _springNodes[i].position[0], _springNodes[i].position[1], _springNodes[i].position[2] );
 glVertex3f( _springNodes[i+1].position[0], _springNodes[i+1].position[1], _springNodes[i+1].position[2] );
}

glEnd();

// Pop the last element off OpenGL attribute stack
glPopAttrib();
}



// String spring system class
class CStringSpringSystem : public ISpringSystem {

public:
void ConnectSprings();

};


void CStringSpringSystem::ConnectSprings( )
{
_springNodes.resize( 80 );
_movingNodes.resize( 80 );

_timer = new CHiResTimer;
_timer->Init();

_springK = 8000;                    // Rigidity of springs
_inertialDistance = .05;            // Springs are 5 cm apart
_mass = .05;      // Each mass is .05 kg
_energyLoss = 1-.00001;    // Amount of energy the system will lose as time elapses


// Set springNode distances along the positive x-axis
for ( int i=0; i<_springNodes.size(); ++i )
{
 _springNodes[i].position = CVector3<double>( 0.0f + i*_inertialDistance, 1.0f, -5.0f );
 _movingNodes[i] = &_springNodes[i];
 _movingNodes[i]->velocity.Zero();
}

// Set up the first spring node
_springNodes[0].neighbors.resize(1);
_springNodes[0].neighborDistances.resize(1);
_springNodes[0].neighbors[0] = &_springNodes[1];
_springNodes[0].neighborDistances[0] = _inertialDistance;


// Set up the last spring node
_springNodes[_springNodes.size()-1].neighbors.resize(1);
_springNodes[_springNodes.size()-1].neighborDistances.resize(1);
_springNodes[_springNodes.size()-1].neighbors[0] = &_springNodes[_springNodes.size()-2];
_springNodes[_springNodes.size()-1].neighborDistances[0] = _inertialDistance;


// Set up all nodes in between
for( int i=1; i<_springNodes.size()-1; ++i )
{
 _springNodes[i].neighbors.resize(2);
 _springNodes[i].neighborDistances.resize(2);

 //
 _springNodes[i].neighbors[0] = &_springNodes[i-1];
 _springNodes[i].neighborDistances[0] = _inertialDistance;

 //
 _springNodes[i].neighbors[1] = &_springNodes[i+1];
 _springNodes[i].neighborDistances[1] = _inertialDistance;
}

// Make the middle node an anchor node (i.e. it does not move)
AddAnchorNode( (_springNodes.size()-1)/2, CVector3<double>( 0.0f, _springNodes[(_springNodes.size()-1)/2].position[1], -5.0f ) );
}