1

Particle Systems Okke Schrijvers

Abstract—In this paper we explore different techniques that are used to implement a 2D particle simulation. We start by given an overview of the system from which the structure of the rest of the report emerges. The problem of a particle system is defined as a ODE so we start by exploring numerical integration schemes to solve the ODE. We then continue by describing unconstrained motion of the particles using forces, including mouse interaction. Then constraints are added to enforce legal positions. We continue by looking at collision detection and finally we handle unconstrained motion of rigid bodies. We conclude the paper by looking at different applications and the effect of parameter settings like the timestep and spring stiffness.

I. I NTRODUCTION In the course 2IV15 S IMULATION IN C OMPUTER G RAPH different ways to simulate natural phenomena are explored. In this paper we shall discuss the first project of this course, which is about 2D particle systems. A particle system is a physics simulation where bodies have a position and velocity, but no shape or orientation. The simulation can be manipulated by applying forces to the particles. Even though the particles have no orientation on their own, one can still get interesting results by connecting particles using springs or constraints. One could for instance create cloth or a soft body by connecting a set of particles in particular ways. In this paper we shall look at some of the technical details of implementing such a particle system. We also take it a step further and add both shape and orientation to the bodies, creating in fact a rigid body simulation. We start in section 2 by examining the internal structure of the simulation, this is somewhat implementation oriented and discusses how one could setup the code so that it is easy to add features. In section 3 we handle the numerical integration of force to velocity and then position based on Newtons formula. Section 4 handles unconstrained motion of the particles. In essence it described how different forces operate on the particles. Section 5 handles constrained motion, i.e. it discusses how we can enforce that (positional) constraints on the particles cannot be violated. Section 6 is concerned with collision detection and response. Section 7 details the conversion from particles to rigid bodies. In particular we focus on the effect of an orientation and angular velocity on the unconstrained motion of rigid bodies. We don’t consider constrained motion. Finally we look at some of the applications of the techniques we discussed in section 8. ICS

II. S YSTEM OVERVIEW The particle simulation can be divided into two parts ([WB97c]). On the one hand we have a number of particles, forces that act upon particles, constraints and collision objects. This part of the simulation is usually referred to as the physics

or dynamics world or system. On the other hand we have a solver. The solver is responsible for taking the current state of the system, and updating it based on Newtons famous equation F = ma. In this section we shall handle both of them briefly and in addition we define an interface between the two parts.

A. System The main responsibility of the system is to maintain the set of particles that are part of the simulation. Each particle has a position, velocity, force accumulator and a mass. The mass is a static property of a particle that is used during the simulation to derive the acceleration of the particle. The force accumulator contains the total force that is applied to a particle at a certain time. This is described as a vector where the direction of the vector is the direction of the force, and the magnitude of the vector is the strength of the force. Since we have a 2D simulation, we have a 2D vector. The velocity of a particle is also a 2D vector and it is the derivative of the position. And finally the position is obviously also a 2D vector. Besides the particles themselves, the system also maintains several other objects that influence the behavior of the particles. Force objects are objects that apply a force to one, two or more particles. Unary forces are forces like gravity and viscous drag. They act on each particle individually. n-ary forces like springs apply forces on a fixed set of particles, based on the properties of those particles. In the case of springs, the force is determined by the position of the two particles. Force objects know how to apply the force to the particles, therefore the system only needs to call the APPLY F ORCE function of each force object. In section 4, we handle the force objects in more detail. Constraints are used to really enforce a position on particles. Although this sounds similar to a force object, there is a subtle difference. A force object calculates a force based on the properties of a particle, which then has to compete with all the other forces in the system. A constraint calculates which force it needs to apply in order to remove all the force that will generate an incorrect position, thereby overriding all the other forces in the system. How this is implemented is described in detail in section 5. Collision objects represent solid objects that cannot be penetrated. Whenever a particle does penetrate a collision object, a counteracting force is applied to the particle such that its velocity is reversed and the collision is resolved. We handle the collision detection and resolution in more detail in section 6.

2

System - Particles

A. Euler’s method

Solver derivEval

Numerical integration schemes aim at solving an ODE by approximation. Euler’s method is based directly on the Taylor expansion of a function. The second degree Taylor expansion of a function f is:

- doStep(system, dt)

- Force objects - Constraints

getParticles

- Collision objects

f (x + h) = f (x) + hf 0 (x) + O(h2 )

Fig. 1: A schematic overview of the separation between system and solver.

(2)

So we can approximate a timestep by simply looking at the first derivative of f and leave an error of order O(h2 ). B. Midpoint method

B. Solver The solver is responsible for solving the equation F = ma. ¨ = F/m, from which it is clear This can be rewritten to x that this is a second order differential equation. We can create a first order ODE by having an additional variable for the velocity, which gives us x˙ = v and v˙ = F/m. When we rewrite this in vector notation we can really see that this is a first order ODE:     v d x = dt v F/m

(1)

The solver itself is now responsible for advancing the simulation. It does so by exposing a function DO S TEP(system, dt) that advanced the given system by time step dt. How this is done is described in section 3.

C. Interface When separating a system into two different parts, it is important to define an interface between them. Since the solver
v needs to solve Equation 1, it needs to be able to evaluate F/a
x as well as have access to v . The evaluation of the derivative is done through the function DERIV E VAL of the system. For access to the position and velocity of the particle, we simply give the solver access to the entire set of particles. A schematic overview can be seen in Figure 1.

III. S OLVER In this section we discuss the solver of the particle simulation. The solver for advancing the simulation
is responsible
v d x by solving dt = . This cannot be done analytically v F/m as the precise forces are not known beforehand, e.g. consider a user that interacts with the particles, so this has to be done numerically. In this section we shall first examine several numerical integration schemes that have been implemented, namely Euler’s method, the Midpoint method and Fourth order Runge-Kutta. We do not give full derivation, for that we refer the interested reader to [WB97b]. Additionally in the final section of this paper we provide some real test data on the performance of each method. This is done because it relies on some of the concepts in the following sections.

The Midpoint method first takes half a step in the direction of the derivative, it then evaluates the derivative at that point and uses that to take a step. The update rule now becomes: h f (x0 )) 2 The methode requires two evaluations of the derivative and the error that you make is O(h3 ) shown in [WB97b]. x(t0 + h) = x(t0 ) + h(f (x0 +

C. Fourth order Runge-Kutta So the Midpoint method improves the O(h2 ) error to O(h3 ), but we can do even better than that. The fourth order RungeKutta method has an error of O(h5 ). The update rule is: k1 = h(f (x0 , t0 )) k1 h k2 = h(f (x0 + , t0 + )) 2 2 k2 h k3 = h(f (x0 + , t0 + )) 2 2 k4 = h(f (x0 + k3, t0 + h)) 1 1 1 1 x(t0 + h) = x(t0 ) + k1 + k2 + k3 + k4 6 3 3 6 So Runge-Kutta 4 requires four evaluations of the derivative. IV. F ORCES In this section we handle unconstrained motion of particles, i.e. we calculate forces based on a particle’s position but we do not create hard constraints. We start of by discussing the structure of how forces are implemented, we then continue by briefly handling each of the implemented constraints. Note that when handling the forces individually we focus on how the forces are calculated, rather than on the specific implementation. A. Force structure Since forces influence the evaluation of the derivative of the system, it is only natural that is part if the DERIV E VAL function, as given in Listing 1. Clearing the force accumulators is simply setting the force accumulators of all the particles to the zero vector. Applying the forces, means that each force object calculates the force on all the particles that it influences. In order to do so in a generic way, we have create a force super

3

class that has a function APPLY F ORCES(particles), where particles is the set of particles in the system. Forces like gravity act on each of the particles, other forces like springs store the particles on which they act. Each of the force objects knows how to apply force on the particles, and adds this to the force accumulator F . The validity of this approach can be found in any elementary physics book, e.g. [FLS63], and we shall not repeat the argument it here. Finally we copy v and F/m to the newState which is used by the solver. DERIV E VAL (newState)

1 2 3

Clear force accumulators Apply forces Copy to newState

When a user interacts with the simulation it is undesirable to directly displace the particles. Forces that act on position and velocity might act in an unexpected way and constraints assume that the current situation is valid, which might not be the case if we act directly on the particles. This is why rather than updating the particle immediately, we apply a spring force on the particle, again using Hooke’s law. We now however do not include a damping factor as the mouse force is not present for a long time. Instead we apply Hooke’s Law immediately by taking the difference between mouse position and particle position as the displacement, effectively taking a rest length of 0. This gives us the following equation for the force: f = −ks (x − m)

(4)

where x is the particle position and m is the mouse position.

Listing 1: Incorporating force objects in the derivEval function.

B. Gravity When near the surface of a large body, e.g. a planet, gravity can be seen as a constant pull in a certain direction: f = mg. Therefore the gravity force object has one field to hold the 2D vector g and for each particle pi we add a force of mi g. C. Viscous drag Viscous drag is an approximation of air resistance on particles. It slows particles down by applying a force, proportional to the velocity of the particle, but in the opposite direction. The viscous drag force object has one field, the drag coefficient kd . The force that is applied on particle pi is then −kd vi . D. Spring A spring is the first binary force that we handle. A spring is defined between two particles and it has a rest length r, a spring constant ks and a damping constant kd . In essence the spring want to keep the two particles a distance r apart. If they are closer, a force is exerted on both particles driving them away. If they are further away, they are pulled together. We implement springs using Hooke’s Law that is given by f = −kx, where x is the displacement vector with regard to the equilibrium position and k is the spring constant. Using this equation in the case of two particles at positions a and b and including a damping factor, we get: ˙l · l fa = − ks (|l| − r) + kd |l|

E. Mouse force

!

l , |l|

fb = −fa

(3)

where l = a − b, as proposed in [WB97c]. The last part of the equation represents the direction from b to a. The part multiplied by ks is the displacement with regard to the equilibrium position. The part that it multiplied by the damping factor can be seen as the projection of the relative velocities of the particles on the line joining them. A small force is applied in the opposite direction.

V. C ONSTRAINTS In the previous section we have seen how we can apply forces to particles to get a nice simulation for unconstrained motion. In many cases however, we would like to impose constraints on particles, for instance we would like to make sure that two particle are always exactly a certain distance apart. If we use a spring to do this, than we always have a certain delay. At first the particles move apart, then the spring will generate a counteracting force, and then the particles move closer together. If we want to decrease this delay we might create a stiffer spring, but this can become unstable easily. One of the problems is that the forces generated by the force objects have to compete with all other forces in the simulation. So what we are going to do instead is calculate the force that is needed to counteract all other forces that are in a direction violating the constraint. In the first subsection we handle the general structure of the constraints functionality. We then handle the three different constraint that have been implemented. A. Constraints structure By now we have mentioned twice that constraints calculate a counteracting force on the particles to ensure that the constraints are not violated, but we have not discussed how exactly we are going to do this. Our approach follows the description of [WB97a]. We start by defining a scalar energy function C for the constraint, that is 0 if and only if the constraint is valid. Furthermore C˙ is 0 if and only the particle has a legal velocity and C¨ = 0 if and only if the particle has a legal acceleration. We enforce the constraints using the method of lagrange multipliers. This is a well known method of implementing constraints. We shall not repeat the complete argument here, for that we refer the interested reader to [WB97a], but we need to solve the following equation, where the last two terms are added to prevent numerical drift: ˙ JWJT λ = −J˙ q˙ − JWQ − ks C − kd C

(5)

4

Where C is a vector containing the scalar values of all ˙ is a vector with all the values of the the constraints, C time derivative of C. W is a diagonal matrix containing the inverse masses of the particles, J is the Jacobian matrix of the constraint functions. This is a sparse matrix that has the dimensions of the number of constraints by twice the number of particles (for a 2D simulation), and it holds the partial derivatives of the constraints with respect to the particles that appear in the constraints. Q is a vector containing the forces on the particles, and q is a vector containing all the particle positions. Finally λ contains the lagrange multipliers. It should be clear that for each constraint, we need to be able ˙ its to evaluate the scalar function C itself, its time derivative C, ˙ partial derivatives J and the time derivative of the Jacobian J. Then Equation 5 simplifies to a linear system Ax = b where A and b are known and we need to find x. Here we can use the sparsity of JWJT and apply the Conjugate Gradient method to efficiently solve this system. We also need to update the DERIV E VAL function to include the constraints. Since we need the current forces on the object in Equation 5, this should be done after the calculation of all other forces. The updated function can be found in Listing 2.

∂C ∂C −2(p1x − p2x ), ∂p1 = 2(p1y − p2y ), ∂p2 = −2(p1y − p2y ). y y For each derivative we have apply the chain rule. Finally taking the time derivative is means simply replacing the position by the velocity in each equation.

D. Horizontal line constraint The horizontal line constraint is the easiest constraint that we have implemented. The constraint function is simply C(p) = py − h where h is the height of the line. It is easy to ∂C ˙ ˙ see that C(p) = vy , ∂C ∂x = 0, ∂y = 1 and J = 0. VI. C OLLISION DETECTION AND RESPONSE So by now we have got a nice simulation where we can influence particles using forces like springs and the mouse, and we can even constrain the particles to a certain position, line or distance using constraints. In this section we add collision objects that particles bounce off of. As always we start by describing the structure of the implementation, and then we handle the two different collision objects separately.

DERIV E VAL (newState)

1 2 3 4

Clear force accumulators Apply forces Enforce constraints Copy to newState

Listing 2: Incorporating constraints in the derivEval function. Constraints are implemented in a base class that has func˙ J and J. ˙ Enforcing the constraints tions that give C, C, now means filling in everything in Equation 5 and calling a Conjugate Gradient solver to solve the system. We then obtain ˆ = JT λ. the constraining force Q B. Circular wire constraint A circular wire constraint keeps a particle on a circle without friction. It has a center c and a radius r. We can then define the energy function C(p) = 12 ((p − c) · (p − c) − r2 ). ˙ Taking the time derivative gives C(p) = 12 v · (p − c) + 1 2 (p − c) · v = v · (p − c). Taking the partial derivatives is easier when rewriting C to C(p) = 12 (p2x + p2y − r2 ). Then ∂C ∂C ∂x = px and ∂y = py . Taking the time derivative of that simply gives vx and vy . C. Rod constraint The rod constraint keeps two particles a certain distance r apart. The energy function is given by C(p1, p2) = (p1x − p2x )2 + (p1y − p2y )2 − r2 . Taking the time d ˙ derivative yields C(p1, p2) = dt (p1 − p2) · (p1 − p2) − 2 r = (v1 − v2) · (p1 − p2) + (p1 − p2) · (v1 − v2) = 2(p1 − p2) · (v1 − v2). We now have four partial derivatives to evaluate. Fortunately ∂C ∂C = 2(p1x − p2x ), ∂p2 = they are quite straightforward. ∂p1 x x

A. Collision detection structure As the title of this section suggests, we are really dealing with two different parts here. One the one hand we need to be able to detect collisions, on the other we need to respond to that collision by pushing the particle away. At first we have tried to implement the collision response by adding a constraint to force the particle outside the collision shape, but this sometimes caused particles to stay attached to the collision object. Additionally when updating the velocity directly it was observed that other constraints would be violated. We are unsure exactly why this happened, but we have opted to use a penalty method instead. Collision objects are implemented in a base class that has a function GENERATE C OLLISION F ORCES(particles) that iterates over the particles, and if a collision between a particle and a collision object is detected, a correcting force is generated. Unlike the other forces and constraints, the calculation of the correcting force is identical for all collision objects: f = n · k · d2

(6)

where n is the normal of the collision, i.e. it is the direction of the shortest path out of the collision object, k is a stiffness constant and d is the penetration depth. We also need to update our DERIV E VAL function to include the penalty force of collisions. Since constraints need all forces currently acting on the particles, it should come as no surprise that the collision detection and response needs to be added before that. The update function is given in Listing 3 In the following two subsections we describe how we calculate the normal and penetration depth for a half-plane and a disk.

5

by looking at how the orientation of an object changes over time and we conclude by looking at how forces act on the orientation of bodies. We do not give justifications for the statements we make, the interested reader is referred to the classical text of Feynman, [FLS63].

DERIV E VAL (newState)

1 2 3 4 5

Clear force accumulators Apply forces Perform collision detection and response Enforce constraints Copy to newState

Listing 3: Incorporating collision detection in the derivEval function.

y

n n

d

p x

O

Fig. 2: A half-plane collision object in 2D. The normal n can be displaced to the origin. We now find the distance of the line to the origin by taking the projection of a point on the line, on the normal: d = n · p.

B. Half-plane collision detection A half-plane collision shape can be seen as a line separating the plane into two parts. Particles are only allowed on one side of the line. A half-plane can thus be described by a normal n, and a distance to the origin d. This can be seen in Figure 2. Checking for a collision is now easy; if p is the position of the particle then we have a collision if n · p > d. The normal to the collision is the same as the normal of the half-plane, and the penetration depth n · p − d. C. Disk collision detection A disk is described by a center c and a radius r. If p is the position of the particle, then we have a collision if and only if |p − c| < r. The collision normal is the normalized vector n ˆ = p − c and the penetration depth is d = r − |p − c|. VII. R IGID BODIES In this section we handle what is needed to extend the unconstrained particle motion demo to get rigid body simulation. The main difference between a particle simulation and a rigid body simulation, is that the dynamic physics objects of the simulation have a shape rather than being a point mass. Because they have shape, they also have an orientation. In this section we shall first look at how to represent the shape and orientation of an object in 2D. We then continue

A. Shape and orientation As a base for our rigid body simulation, we take the particle class of the particle simulation. We augment this with a field for the shape and a field for the orientation. We have implemented only one shape which is a rectangle. The rectangle is defined around a center point called the center of mass. We assume that the distribution of mass is uniform, so the rectangle can be described by its half-extends. Since we are only dealing with two dimensions, we can describe the orientation of a rigid by as a scalar that represents the radians the object has turned counterclockwise with regard to its original orientation. B. Torque, inertia and angular velocity Now that we have an orientation that can change during the simulation, it is useful to look at how this orientation changes. It turns out that there is almost a complete analogue to the position. The time derivative of position x is linear velocity v, whereas the time derivative of orientation R is angular velocity ω. The time derivative of linear velocity v is mass m times force f , whereas the time derivative of angular velocity ω is the moment of inertia I times torque τ . The moment of inertia can be seen as a quantity that describes how difficult it is to rotate an object. It is defined to be the area integral of the distance R to the center of mass squared times the mass at that point: A |r|m. In the case of a rectangle in 2D with uniform mass distribution, it can be derived that this is I = mwh 12 where m is the mass, w is the width and h is the height of the rectangle. Now that we have an orientation and an angular velocity, we need to feed this to the ODE
solver.v For
a particle system, the d x = ODE to be solve was dt v F/m . This is now extended to:     x v    d   v  =  f /m  (7) dt  R   ω  −1 ω I τ C. Applying a force When applying a force at the end of an object has a different influence than when a force is applied at the center. This is because applying a force somewhere besides the center of mass introduces torque. This means that we have to make sure that if we apply a force somewhere, we separate this into a (linear) force on the center of mass and torque. Without giving any justification, we again refer the reader to [FLS63], the torque that is generated is by applying force f at the point with regard to the center of mass r is τ = rx fy − ry fx . Which is the 2D equivalent to taking the cross product between r and f . In the ODE solver, this torque is eventually transferred to a different orientation.

6

We look at three demos specifically, namely a 2D cloth simulation using springs, a 2D soft body simulation again using springs and particles and finally three rigid bodies that are interconnected using springs. A. Suspended cloth

Fig. 3: A series of rigid bodies is connected by springs.

D. Updating the spring In the particle simulation, springs between two particles were always attached at the center of mass. If we have a rigid body, we are free to choose an anchor point on each of the two rigid bodies. The precise behavior of the spring is now determined by more than just the position and velocity of the center of mass. When determining the position of the anchor point, we now also need to take the orientation of the rigid body into account. If r is the offset of the anchor point in default position and orientation, then the current position p0 of that point is:   cos R − sin R p0 = r+p (8) sin R cos R Additionally rather than only using the linear velocity, if the anchor point is not at the center of mass, we also need to take the angular velocity into account. The added linear velocity, the point at updated offset r0 is   −ry ω· (9) rx So to put it all together, when calculating the force generated by the spring, we need an anchor point for each rigid body. We then calculate the position of that anchor point in world coordinates using Equation 8. We calculate the linear velocity of those two anchor point by adding Equation 9 to the linear velocity of the center of mass. And finally we apply the force generated by the spring to the anchor points rather than the center of mass. A screenshot of the rigid body simulation with springs can be seen in Figure 3. VIII. A PPLICATIONS In this section we shall look at how we can use the techniques described in this report to get interesting demos.

We implemented 2D cloth using only particles and springs, the results can be seen in Figure 4. In [Pro95], three different types of springs are identified to simulate cloth in 3D; structural springs, shear springs and flexion springs. Since the simulation is only in 2 dimensions, we found that only using structural springs gives the best cloth-like effect. Using the shear and flexion springs gave rise to a stiff simulation. For the simulation, we determine the stiffness of the spring by the height at which the spring appears. So higher springs are stiffer than the lower springs. Although this is not exactly physical, it is more convincing because it leads to less deformation on the higher springs. The stiffness of the lower springs is set to 1.0 and this is scaled linearly to 25.0 for the highest springs. The downside of this stiffness is that Euler integration very quickly becomes infeasible. If we use a timestep of 0.005 is is still quite easy to let the simulation blow up. Given that the spring stiffness is at most 25, it seems odd since we are nowhere near the theoretical bound of hk > 1, since hk = 0.125. The problem is that there are always multiple springs connected to each particle. Springs can enhance each other by acting on the particles in turn. If the timestep is reduced to 0.003 no explosions have been observed. When using the Midpoint method, we can safely increase the timestep to 0.03 without any problems. Incidentally for the Midpoint method, 0.05 is too large a step. The simulation at this point does not actually blow up, but it start shaking uncontrollably, so we cannot speak of a correct physics simulation. If we look at the error we make using Euler’s Method, it is in the order of h2 , so for a timestep of 0.005, the error is in the order of 2.5 × 10−5 . The error of the Midpoint method is h3 so for a timestep of 0.03, this is about 2.7×10−5 , nearly the same as for Euler’s error! It might seem wrong that the Midpoint method is able to function with a larger error, but this is probably due to a lower constant factor in the error. Finally when using Runge-Kutta fourth order, we can safely increase the timestep to 0.3 with a completely stable simulation. However increasing this to even 0.4 leads to an immediate explosion. So changing from Euler to Midpoint requires an additional evaluation of the derivative, effectively doubling the computation time, but we can make a ten times larger timestep. Then again by changing to Runge-Kutta 4, we go from two evaluations to four, again doubling the computation time, but again we can make a timestep that is 10 times as large. Therefore effectively, if we use Runge-Kutta 4, the simulation can run at 25 times the speed of Eulers Method. B. Man soft body The cloth in the previous used only structural springs to create a soft cloth. If we however also use the flexion and shear springs as defined in [Pro95], we can get a lot stiffer figure that can be shaped at will. This is the basis for the man

7

Fig. 4: A piece of cloth is suspended from a rail. It consists of 25 × 25 and 2 × 25 × 24 structural springs.

Fig. 5: Using stiffer springs and flexion and shear springs can be used to create arbitrary 2D soft bodies.

[WB97c]

figure that can be seen in Figure 5. Especially in combination with collision objects and mouse interaction, one does not only have a nice simulation, but one that is also (perhaps a bit sadistically) fun! The figure consists of 56 particles and 189 springs, each with a stiffness of 25. Although we have less springs than the cloth simulation, all particles are more interconnected so we cannot use the same timestep as for the cloth. Decreasing the timestep to 0.15 and using Runge-Kutta 4 leads to a stable simulation. What is especially nice to see, is that all the properties described in the section on rigid bodies can also be observed in the soft body. If you drag along the man at his center of mass he hardly rotates, if you grab him at the side, you always introduce rotation. In fact we have used the fact that we can approximate the forces on a body by only looking at the external forces when working on the rigid body, but we have not named it as such (it can however be found in [FLS63]). C. Rigid body chain The final application is a chain of rigid bodies, connected by springs, as can be seen in Figure 3. The effect of the anchor point is clear from the R EFERENCES [FLS63]

Richard Feynman, Robert Leighton, and Matthew Sands, The feynmen lectures on physics: Volume 1, 2nd edition ed., The Feynman Lectures on Physics, vol. 1, Addison-Wesley, Boston, 1963. [Pro95] Xavier Provot, Deformation constraints in a mass-spring model to describe rigid cloth behavior, IN GRAPHICS INTERFACE, 1995, pp. 147–154. [WB97a] Witkin, A. and Baraff, D., Constrained Dynamics, Physically Based Modeling: Principles and Practice, 1997. [WB97b] , Differential Equation Basics, Physically Based Modeling: Principles and Practice, 1997.

, Particle System Dynamics, Physically Based Modeling: Principles and Practice, 1997.