Time integration: Verlet ODE integration method is known to have certain stability but not a big deal of accuracy. However, is fairly easy to implement. An exhaustive explanation on how it works can be found here: http://en.wikipedia.org/wiki/Verlet_integration. In our implementation, it is divided into functions: AccumulateAccelerations and Verlet.
In AccumulateAccelerations we iterate through each mass and add up earthquake readings from a file (per millisecond), gravity acceleration (980 cm/s2), and the internal forces caused by deformation divided by the weight of the mass.
The Verlet procedure, simply updates the position of each particle according to the formula:
x(t)=x(t-1) + v(t)·dt + a(t)·dt^2
Constraint integration: The only constraints included in the current implementation are those of the distance between two particles. By means of the concept of impulse, once two particles move close or apart by the effect of the accelerations, a corrective force is applied on each one of them. This force is applied in a very short lapse of time, hence can be regarded as an impulse. A more detailed explanation can be found in the work by Jan Bender (www.impulse-based.de). The algebraic system of equations that arises so as to satisfy all of the constraints simultaneously is solved by means of a Gauss-Seidel iterations.
Matter integration: For the consideration of matter properties the program uses the Finite Element Method. A quadratic shape function is used according to the implementation taught in http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/. Chapters 20 and 21 explain clearly how to infer the values for the case of an elastic beam.
This provides our system with a rigidity matrix that is easily attached to each of the distance constraints. As the displacements have already been obtained previously, the solution of the force vector is straightforward by means of a simple rigidity matrix-displacements vector operation:
[f] = [K] · [u]
This force vector stores the effect of the deformation into the constraint and is later used to obtain the particle acceleration in the next timestep.