The first time I implemented a spring model, it fascinated me for hours. The effect itself is amazingly realistic, and its implementation is fairly simple. Fortunately, I found a lot of articles, references, and source code to help with my research. Nevertheless, as I went further down the road, I noticed that in most cases these references limited to the standard applications of spring models -- string, cloth and jelly.

This article reviews the implementation of a spring model from its simplest form to more sophisticated applications, taking the subject a step beyond the material available in most references.

Spring Basics

Before modeling springs with the computer, you need to understand the basic principles of springs from your classic physics book.

As you compress or extend a spring, it creates a force that is opposed to direction fo the force that you are applying. This force is mathematically equated by the formula:

F = -k*Dx
Dx = xc-xi

Where F is the resultant force, k is the spring coefficient, and Dx the distance from its inertial position (xc = current distance, and xi = distance at the inertial position).

Mathematically, you can think of a spring as being two points separated from a distance x. This is its inertial position. In other words, if you don't move these points, no force will be generated because Dx = 0. If you try to compress the distance between these points the Dx will be negative, generating a positive force (expansion). If you try to separate these points the Dx will be positive, generating a negative force (compression). The picture below illustrates these cases.

Pa -/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\- Pb (Dx = 0, F = 0)
Pa -/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\- Pb (Dx < 0, F > 0)
Pa -/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\- Pb (Dx > 0, F < 0)

Springs (using "ascii art")

The k coefficient is called the elasticity coefficient, and weights the spring's final force. In other words, a spring with bigger k is stiffer because it creates a larger force from its inertial state. Conversely, a spring with a smaller k is more flexible because it creates a smaller force from its inertial state.

Spring Force Implementation

Implement the spring force is a direct application of the first two equations. In other words, first compute the distance of the spring's extremities (Pa and Pb). This is done through Pythagoras' theorem. Once you have the distance, you subtract from the inertial distance getting the Dx. Then, multiply this scalar value to the k coefficient and finally use the inverse of this value to compute the force.

As simple as it sounds, the implementation in C code has its details. First, the distance Dx, and the k coefficient are scalars, but the force itself is a vector. Therefore, it requires a direction, which can be from Pa to Pb, or from Pb to Pb. As of right now, the code below always use the direction Pb to Pa, or mathematically Fdir = Pb-Pa.

As the force's direction is normalized, you can use the distance you already computed for the Dx for the normalization. You must be careful is to avoid division by zero when normalizing the vector. This is a nasty condition and it usually does not happen, but you have to do something about it in case it does. You have two options in this particular case: assign a random force to the points, or null out the force. Theoretically, the force should be big enough that as the points get closer they will never occupy the same space. Therefore, it seems to be a good option to randomly assign a massive force in this case. However, this can cause undesirable artifacts, because some points inside the spring model might fly a part with extreme velocity when compared with others. Although nulling out seems wrong, it doesn't cause any major problems and preserves the spring's stability. The code below shows the implementation; it's a little redundant for sake of clarity. I null out the force in case of potential division by zero.

Listing 1: Computing A Spring's Force

#define SPRING_TOLERANCE (0.0000000005f)

void SpringComputeSingleForce (VECTOR *force, VECTOR *pa, VECTOR *pb, float k, float delta_distance)
{
VECTOR force_dir;
float intensity;
float distance;
float delta;

//get the distance
float dx = pb->x-pa->x;
float dy = pb->y-pa->y;
float dz = pb->z-pa->z;

distance = (float) sqrt ((dx*dx)+(dy*dy)+(dz*dz));

if (distance < SPRING_TOLERANCE)
{
force->x =
force->y =
force->z = 0.0f;

return;
}

force_dir.x = dx;
force_dir.y = dy;
force_dir.z = dz;

//normalize
force_dir.x /= distance;
force_dir.y /= distance;
force_dir.z /= distance;

//force
delta = distance-delta_distance;
intensity = k*delta;

//store
force->x = force_dir.x*intensity;
force->y = force_dir.y*intensity;
force->z = force_dir.z*intensity;
}