Optimization of clubhead loft and swing elevation angles for maximum distance of a golf drive

Pergamon

Compurerx & Srruc~rs Vol. 53. No. I. pp. 19-25. I994 Copyright g; 1994 Eloevier Science Ltd Printed in Great Britain. All rights reserved 0045-7949194 $7.00 + 0.00

00457949(94)EO179-6

OPTIMIZATION OF CLUBHEAD LOFT AND SWING ELEVATION ANGLES FOR MAXIMUM DISTANCE OF A GOLF DRIVE D. C. Winfield? and T. E. Tanf tResearch Engineering, 935 Hiawatha, Memphis, TN 38117, U.S.A. SDepartment of Mechanical Engineering, Memphis State University, TN 38152, U.S.A. (Received 7 July 1993)

Abstract-This paper presents a study on the optimization of loft and swing elevation angles on the impact between the clubhead of a driver and a golfball in order to maximize the distance of a drive. Computer programs were written to simulate the collision between the golfball and clubhead as well as the golfball in flight. A general, three-dimensional impact model using principles of momentum conservation on rigid bodies was used to simulate the impact between the golfball and clubhead to extract the spin and velocity vectors of the ball after impact. An aerodynamic model was then used to simulate ball flight in order to obtain the landing position of the ball. The results of the golfball landing positions generated by the computer simulations were compared to experimental data of golfball landing positions of shots hit by the golfing robot ‘Iron Byron’. The computer models were then used to calculate the optimal loft and swing elevation angles for a particular swing speed, clubhead mass, and golfball aerodynamic properties by making use of a nonlinear optimization routine. Also, the relationship between the maximizing distance for various driver loft angles and swing elevation angles is discussed.

INTRODUCTION

MATHEMATICAL

There are many interesting problems that can be solved using methods of optimization. One involves the optimization of impact conditions between the clubhead of a driver and a golfball. Generally speaking, one of the goals for every golfer on par-four and par-five holes is to hit the longest, straightest drive humanly possible. A statistical study by Cochran and Stobbs stated that, in tournament play, the distance from the tee was more important than accuracy [l]. Extensive studies have been done that modeled the golfball during flight by taking into account the aerodynamic forces on the ball. Some studies also give recommendations on the initial trajectory angle to produce the longest possible drive for a given ball spin rate and velocity [14]. However, none have obtained the initial trajectory angle, ball spin and ball speed, nor determined the loft and swing path of the clubhead through a computer simulation of the impact with a clubhead. One aspect of this study is to analyze the clubhead impact with the golfball in order to obtain the ball spin and velocity and then to model the golfball during flight taking into account the aerodynamic forces. The maximum possible distance of a drive for a given clubhead and swing speed will be determined from the calculated ball spin and velocity resulting from the impact. An optimization routine will be used to find the optimum combination of the chtbhead loft angle and the swing elevation angle in order to maximize the driving distance.

Impact

MODELING

model

The impact between the golf clubhead and golfball will be modeled using rigid body mechanics. It has been shown by Milne that the shaft plays a negligible role during impact [5]. The shaft being a flexible member does not add any appreciable impulse and does not greatly affect the dynamic response of the clubhead during impact. The main purpose of the shaft is to get the clubhead in position at impact with the maximum possible kinetic energy. Figure 1 shows two views of the clubhead. The reference frame attached to the clubhead H at the center of mass C is a set of mutually perpendicular axes, labelled X, Y and 2. Three unit vectors, I, J and K, are directed along the X, Y and Z axes, respectively. The X axis indicates the vertical direction that is perpendicular to the Z axis which is parallel to a straight line along the middle of the fairway. Due to the curvature of the clubhead surface, a set of mutually perpendicular unit vectors, i, j and k, directed along the X, y and z axes, respectively, are introduced. The k unit vector is along a line that is normal to the clubface at the point of impact, and i and j are parallel to two mutually perpendicular tangents at the point of impact. The i, j, k unit vectors can be related to the I, J, K unit vectors by utilizing the loft/elevation angle 0 and the bulge angle 4 as shown in Fig. 1. The loft/elevation angle 0 is the combination of the inherent loft angle of the clubhead I and 19

20

D. C. Winfield and T. E. Tan colliding bodies have the same material properties and that the collision is totally elastic are made. However, the impact momentum model does provide a versatile and elegant way to analyze the dynamics of an inelastic collision and has been used for many applications. The translational velocity V’ at the clubhead center of mass and the angular velocity W” are given by V”=u,I+u,J+u,K

(1)

~“=u~I+u~J$u~K,

(2)

and

where u,, us.. . . , u6are called the generalized speeds of the clubhead [lo]. The translational velocity of the center of mass of the golfball VB and the angular velocity aG are given by VB=~,I+~gJ+~9K Fig. 1. Front and top views of clubhead and the location reference frames X-Y-Z, x’-y’-z’ and x-y-r.

the swing elevation angle t. The elevation angle L is the angle the bottom of the face has with respect to the horizontal Z direction just prior to impact, and is produced and influenced by swing dynamics, shaft properties, and the mass and mass moment of inertias of the clubhead. The transfo~ation matrix between the two sets of unit vectors can be derived through a rotation of the loft/elevation angle B about the Y axis to form an intermediate set of axes, x’, y’ and z’, with i’, j’ and k’ unit vectors, respectively. A second rotation of the bulge angle 4 about the x’ axis results in the final orientation of the i, j and k unit vectors. The dynamic model will be derived from the principles of generalized momentum and generalized impulse for three-dimensional rigid bodies. Some simplifying assumptions of the model include: (1) the contact forces are impulsive in nature, (2) the collision is assumed to be instantaneous, and (3) no defo~ations occur during the collision. It has been shown that golfball contact occurs in a period of about 470 psec [6] and the peak force can be as high as 5000 lb [7]. The golfball, which has a diameter of 1.68 in, can undergo a deformation of around 0.25 in during impact at a rather low impact velocity of around 65 mph [6]. It should be understood that there are probably no reai cases of collisions occurring in nature that do not violate the simplifying assumptions of the generalized momentum and impulse model of colliding bodies. An alternative method to model collision in the normal direction and include the effects of deformation is to use a refinement by Hunt and Crossley [S] of Hertz’s theory of collision. The tangential effects of the collision including deformation have been analyzed by Maw et al. [9]. The process to model the tangential effects of a collision is complex and simplifying assumptions such as the

(3)

and w”=u,,I+u,,J+zc

I2

K3

(4)

where u,, usI . . , u,* are generalized speeds oi the golfball. The kinetic energy of the clubhead K” and the kinetic energy of the golfball KG are given by KH=:+,H.~“.&

+m,V’~V’]

(5)

and KG = f[oC. IG. wc; + m,VB. VB],

(6)

where IH and IG are the inertia dyadics of the clubhead and the golfball, respectively, and pnH and mG are the mass of the clubhead and golfball, respectively. The generalized momentum pi is given by pi = ~KH~~u~+ i3KG/du,,

(7)

where i is an index from I to 12. Hence, there are 12 equations of generalized momentum. The velocity VP at the impact point P on the clubface and the velocity VP of the impact point P’ coincident with P on the golfball is given, respectively, as 63)

and VP’= VB+ oo X $“:B,

(9)

where flc is the position vector from the clubhead center of mass to the impact point P and fiB is the position vector from the golfball center of mass to the impact point P’.

Optimization of clubhead loft for m~imum distance golf drive The position vector 8”” represents the radius of the golfball which is 0.84 in and is always parallel to the k unit vector. The clubhead impact surface will be modeled as an ellipsoid in order to attain flB for different impact locations. The ellipsoid will be defined with respect to the x’, y’, z’ reference frame and the position vector PiB expressed in the i’, j’, k unit vectors will be rpjB = r, i’ + r2jl + r,k’.

(10)

21

speeds ~~(00).The generalized speeds after impact U!(t) are unknown and must be found. It should be noted that the initial value of the generalized speeds u~(0)-uu(O) for the golfball are all zero. Equation (15) yields 12 scalar equations but 15 unknowns +(t), R,, 4 and %I. L(t)... The remaining equations come from using the velocity of separation V, which is the relative velocity of the colliding bodies just after impact and the velocity of approach V, which is the relative velocity just prior to impact. V, and V, are given by

The equations describing the ellipsoid in the x’, y’, z’ reference frame with its origin at the center of mass of the clubhead is

v, = V’(t) -V(t),

(16)

v, = V’(0) - V’(0).

(17)

and r;/Pf. + r;/P;+ + (r, -d

+ Pk,)2/P;. = 1,

(11)

where d is the distance from clubhead center of mass to the face in the k’ direction, P,, is the radius of curvature in the i’ direction, PJ is the radius of curvature in the j’ direction, and Pkpis the radius of curvature in the k’ direction. An impact that occurs at the ‘sweetspot’ of the clubhead is where r, and r, are equal to zero. Off-center impacts can be specified by nonzero values of r, for shots high or low on the clubface or r2 for shots toward the toe or heel of the clubhead. With r, and r2 specified, r, can be found by solving eqn (11) to give

The thirteenth equation comes from the assumption that eqns (16) and (17) have opposite directions parallel to the k unit vector and their magnitudes are proportional to each other. The constant of proportionality, eN, is called the normal coefficient of restitution which represents a loss of kinetic energy due to complex factors such as wave propagation through the colliding bodies and hysteresis losses due to local deformation. The relationship can he expressed as k&V,=

r3 = [(I - ri/Pf. -ri/Pi2,)P$j’i2

+ d - Pp.

The bulge angle Cpcan be found by taking the inverse tangent of the derivative of eqn (11) with respect to rz. The total loft angle 0 is given by the nominal loft of the clubhead E.at the ‘sweetspot’ plus the swing elevation angle 6 and also plus the angle of the face curvature or roll angle of the ellipsoidal shape. The roll angle, which adds loft for impacts high on the face and substracts loft for impacts low on the face, can be found by taking the inverse tangent of the derivative of eqn (11) with respective to r,. During impact the contact forces are assumed to be impulsive in nature. The time integral of the contact forces produces the resultant impulse R in which R=R,I+R,J+R,K.

(13)

The generalized impulse & is then given as Z,= @V;,‘&,) *R - (aVf’j&)

. R,

(14)

where i is the index that ranges from 1 to 12. The impact momentum equations can be written as 1; = Api,

-e,k.V,.

(18)

(12)

(1%

where Appiis the change in momentum due to impact. The value of the momentum prior to impact is found by knowing the initial values of the generalized

Cochran and Stobbs simply state that the normal coefficient of restitution is 0.8 for a putt and 0.67 for a hard drive [l]. Experimental work by Gobush [6] has shown that the golfball is not slipping when impact ceases for clubhead lofts ranging from 20” to 40”. By assuming that the ball is not slipping, the remaining two equations are given as i. V, = 0

(19)

j.V,=O.

f20)

and

If the ball were slipping at the end of impact, friction factors and tangential coefficients of restitution [l l] could be used to more accurately model the tangential aspects of the collision. Equations (18), (19) and (20) provide the last three equations required to form a system of 15 equations and 15 unknowns which can be easily solved by many numerical schemes such as Gaussian elimination. The solution will include the post-impact velocity and spin of the ball as well as that for the clubhead. Aerod~~~mic model

Once contact is broken between the club and the golfball, the ball experiences gravitational forces as well as drag and lift forces. Drag is caused by air resistance during flight while lift is caused primarily

22

D. C. Winfield and T. E. Tan

by the spinning of the golfball. Drag acts in the opposite direction to the velocity of the golfball. The lift force acts perpendicular to the direction of the velocity of the golfball and the direction of golfball spin which is given by n, = n,* x nVBI

(21)

nI, n,, and nvB are the direction cosines of the liftforce, angular velocity of the golfball, and translational velocity of the golfball, respectively. The magnitudes of the drag force Fd and the lift force F, are given by F,, = C,p(VB

Va)A

(22)

F, = C,p(VB VB)A,

(23)

and

where C,, is the drag coefficient, C, is the lift coefficient, p is the density of air, and A is the projected cross-sectional area of the golfball. C,, and C, have been plotted as a function of the magnitude of spin and velocity of the ball by Aoyama [12]. The applied gravitational and aerodynamic forces FG exerted on the golfball during flight are FG = F,n, - FdnvB - m,gI,

(24)

where g is the acceleration due to gravity. The inertia forces of the golfball F*’ are given by F*G = -m,(&+/dfI

+ &,/dtJ

+ &q/atK).

(25)

This represents the mass of the golfball multiplied by the acceleration of the golfball or the time derivative of the velocity of the golfball. By using Kane’s method [IO], the three equations of motions can be derived as (F’ + F*G) I = 0

(26)

(FG + F*‘).

J = 0

(27)

(Fo + F*‘).

K = 0.

(28)

and

During flight, a skin friction force applies a torque that slows down the rotation of the golfball. Some studies by Lieberman [13] have shown that the spin decays exponentially from its initial spin rate wG(t = 0) and is given by uG = ~~(2 = O)exp( -&), where S is the decay coefficient 0.6 to 0.7 for most golfballs.

(29)

that can range from

Equations (26), (27) and (28) can be integrated by using a number of techniques such as the Runge-Kutta method. At each time step, the values for C,, C,, and mG will have to be updated. The integration procedure will be carried out until the golfball lands which is when the I component of the golfball displacement is zero. The distance from the center of the fairway is given by the displacement along the J direction. OPTIMlZATlON

STUDY

Given the number of variables present in the impact and trajectory models of this study, there are many possible optimization problems that can be explored. In order to conduct the optimization process, certain parameters are fixed and others are allowed to the vary so as to minimize a specific objective function. One study is to maximize the distance of a drive down the center of the fairway. An assumption in this study is that the impact will always occur at the the ‘sweetspot’ of the clubhead. The ‘sweetspot’ is the position on the clubface where the normal vector passes through the center of gravity. Hence, the ball will only have backspin and the ellipsoidal clubface shape will have no effect on the ball spins or velocities after the impact. The general three-dimensional impact model could be reduced to a two-dimensional one in order to solve the proposed problem. However, the three-dimensional impact program could be used to solve a wide range of problems involving golf clubhead and ball collision not considered here. The initial velocity vector of the clubhead will be assumed to be in the Z direction rotated slightly through the elevation angle t. It should be noted that the spin induced on ball is due to the direction of the clubhead velocity vector being different from the direction normal to the face at the impact point. Since the velocity vector is assumed to be in the direction of the elevation angle 6, the spin induced on the ball results only from the clubhead loft angle i. For a driver or another wood, this is probably a good assumption. However, for an iron this would be a poor assumption because most golfers frequently attempt to induce more spin by hitting ‘down’ on the ball, especially with short irons. Since a golfball is almost always placed on a tee when using a driver, the elevation angle t accounts for the fact that a golfer can angle the clubhead prior to impact through slight manipulation of the swing dynamics. The parameters to be optimized in this study will be the clubhead loft angle i and the swing elevation angle t. All the other parameters, such as swing speed, clubhead mass and aerodynamic properties, will be specified and remain constant. The optimization process is to minimize an objective function which is defined as the negative of the distance the ball travels down the fairway in the Z direction.

Optimization of clubhead loft for maximum distance golf drive RESULTS

The optimization was performed using a quasiNewton algorithm EO4JAF contained in a subroutine library by the Numerical Algorithms Group (NAG) [14]. The subroutine minimizes a nonlinear objective function with boundary constraints on the variables. The impact equations were solved using Grout’s factorization method. The ball flight equations were numerically integrated using a Runge-Kutta-Mearson subroutine DOZBHF also found in a subroutine library by NAG [14]. The RungeKutta-Mearson method is useful since the integration continues to advance at a variable time step until a stopping criterion is met. In the case of golfball flight, the integration proceeded until the X displacement of the golfball equalled zero indicating that the ball golfball had landed on the ground. In order to validate the computer simulation of impact and ball flight, the golfball landing position struck by a 203-g driver with a loft angle 1 of 10.0” and a swing elevation angle t of 0.0” at a velocity of 103 mph was calculated. The coefficient of restitution was assumed to be 0.7. The calculated landing position of the ball was 197.31 yards in the Z direction. Just after impact the ball was calculated to have a velocity in the vertical or X direction of 19.84 mph and in the Z direction going down the middle of the fairway of 139.16 mph with a rotational velocity or backspin of 38.57 revisec. This meant that initial trajectory angle of the ball was 8.29”. The clubhead slowed down to 71.39 mph in the Z direction after impact. The ball had a total flight time of 4.75 set and had a peak position of 14.18 yards in the vertical X direction at 2.48 set into the flight. The golfing robot ‘Iron Byron’ was then used to impact a golfball at 103 mph with the same driver modeled in the computer at the same swing velocity and swing elevation angle. Twenty trials were performed using the robot in which the distance the ball traveled was 200.3 yards with a standard deviation of 3.61 yards. Hence, the computer simulation predicted a landing position of the golfball that fell within the standard deviation of the landing positions of balls hit by ‘Iron Byron’. However, the velocity of the golfball just after impact and its rotational velocity were difficult to measure and cannot be verified at present. For optimization, a 203-g clubhead and an initial clubhead speed of 110 mph were used. Most drivers are around 197-210 g in mass while a 110 mph swing would be typical on the professional level. The optimization routine determined that the optimum clubhead loft angle and swing elevation angle was 2.74” and 27.44”, respectively. It is interesting that the optimum clubhead loft angle does not fall within the manufactured range of lofts for a driver, which is 7.c12.0”. Also, the swing elevation angle is unrealistic. To achieve a 27.44” elevation angle one would probably need a IO-in tee! The distance achieved by the optimum drive was 250.09 yards.

23

Table I. Optimum elevation angles for various loft angles Loft angle i (“)

Optimum elevation angle t (“)

Distance of drive (yards)

5 I 8 9 10 11 12 13 15

22.09 17.53 15.47 13.54 11.71 9.96 8.29 6.68 3.68

248.48 245.1 I 242.14 239.93 236.68 233.04 229.04 224.15 215.74

The optimization program was adjusted to a single parameter problem in order to find the optimum swing elevation angle for a given clubhead loft to maximize the distance of the drive. The results are shown in Table 1. The optimum elevation angle decreases as the clubhead loft angle increases. Also, as to be expected, the golfball distance decreases as the loft angle increases. By placing the golfball on a tall tee, one could probably attain a swing elevation angle of 10”. For instance, the distance from the shoulder or pivot point of an average golfer to the golfball is about 70 in. In order to achieve an elevation angle of about 10” for the clubhead, the bottom of the clubhead would be a little over an inch off the ground. Hence, it is quite possible to maximize distance with an 11” or 12” driver. The optimization program was also modified to find the optimum loft angle for a given swing elevation angle. The results are shown in Table 2. The optimum loft angle decreases as the elevation angle increases. Also, as the elevation angle increases, the maximum possible driving distance increases. It is interesting to note that for all positive swing elevation angles, the computed loft angles are within the range of the loft angles manufactured for drivers. However, it should be noted that it is no more difficult to produce a swing with a relatively high elevation angle than that with a lower elevation angle. The distance travelled by a golfball was calculated for elevation angles ranging from 0” to 10” at 1” increments and loft angles ranging from 7’ to 12” also

Table 2. Optimum loft angles for various elevation Elevation angle t Optimum loft angle (“1 1 (“1 -10 -5 -4 -3 -2 -1 0 I 2 3 4 5 10

17.14 14.51 14.03 13.65 13.11 12.64 12.21 11.79 11.38 10.99 10.60 10.22 8.39

Distance

angles

of drive

(yards) 188.93 206.08 209.02 211.84 214.52 217.09 219.55 221.92 224. I7 226.33 228.38 230.32 238.55

24

D. C. Winfield and T. E. Tan

MO

7

7.5

8

9 95 10 10.5 85 Clubhead Loft Angle (Lamda)

11

115

12

Fig. 2. Average golfball distance for different clubhead loft angles for a 203 g clubhead swung at 110 mph, with swing elevation angles ranging from 0” to JO”.

at 1” increments. The range of swing elevation angles should encompass the angles that typical golfers produce when driving the ball. The distance of all shots for a given clubhead loft angle was averaged. The results are shown in Fig. 2. The data points in Fig. 2 were curve-fitted to a fifth order polynomial and then the derivative of the polynomial was taken. The optimum loft angle of the clubhead which produced the maximum distance for drives over the given range of elevation angles was determined to be 10.68” corresponding to a 228.01 yard driving distance. It should be noted that the optimum loft angle assumes an even distribution of swing elevation angles ranging from 0” to 10”. CONCLUSIONS

The results showed that to achieve maximum drive distance, the optimum loft angle is beyond the range of manufactured drivers and the elevation angle is too large for a normal swing. The importance of striking the ball with a high clubhead elevation angle was also shown. The results showed that for an even distribution of swing elevation angles ranging from 0” to lo”, the optimum loft of the driver was 10.68” to produce, on the average, the maximum driving distance. The typical golfer certainly does not precisely swing the club at a given elevation angle. Experimentation needs to be performed to determine the range and distribution of swing elevation angles that typical golfers swing the driver. Also, the individual golfer could use the results of the optimization to adjust the loft of his driver to best fit his distribution of swing elevation angles. It is interesting that the optimum loft for the range of swing elevation angles considered is within the typical range of driver loft angles. It should be noted that the results are particular to the aerodynamic characteristics of the specific golfball used in the simulation. As shown by MacDonald and Hanzley (41, different golfballs have different aerodynamic properties. Hence, the results could

easily differ if a golfball with totally different aerodynamic properties were used. The results of the ball landing position for the impact and trajectory programs compared well with the experimental results using ‘Iron Byron’. However, more validation on the models needs to be made. For instance, high speed photography could be used to determine the ball spins and velocities after impact. Also, the coefficient of restitution could also be determined using high speed photography. This study has neglected the rolling distance of the golfball after landing on the ground. That quantity is quite complex to predict because the rolling distance will be affected by factors such as the impact location on the ground, type of grass, condition of the grass, fairway terrain, etc. For instance, on a dry hard surface, a golfball with a higher flight path will certainly have less rolling distance than one with a lower flight path. But, the results presented are valid for a wet and soft fairway when the balls are embedded into the ground upon landing. The impact analysis is quite extensive in that it takes into account the complete inertia dyadic of the clubhead, the three-dimensional rotation of the clubhead, the three-dimensional rotation of golfball, and the clubface bulge and roll. Other types of optimization problems can be explored using this model. One that is currently being studied, is to optimize the face bulge and roll for a given inertia dyadic so that a golfball will land close to the middle of the fairway regardless of the location the golfball impacts the clubface. Also, the interplay between the inertia dyadic of the clubhead and the clubface shape could be investigated. The results of this study will be published at a later date. REFERENCES

1. L. Cochran and J. Stobbs, The Search for the Perjkct Swing, pp. 112-178. The Booklegger, Grassvalley, California (1968). 2. H. Erlichson, Maximum projectile range with drag and lift, with particular application to golf. Am. J. Phys. 51, 357-362 (1983). 3. J. J. McPhee and G. C. Andrews, Effect of sidespin and wind on projectile trajectory, with particular application to golf. Am. J. Phys. 56, 933-939 (1988). 4. W. M. MacDonald and S. Hanzley, The physics of the drive in golf. Am. J. Phys. 59, 213-218 (1991). 5. R. D. Milne, What is the role of the shaft in the golf swing? Science and Golf: Proc. 1st World Scientific Congress of Golf, pp. 252-257. E. & F. N. Spon, New York (1990). 6. W. Gobush, Impact force measurements on golf balls. Science and Golf: Proc. 1st World ScieniiJic Con gress of Golf, pp. 219-224. E. & F. N. Spon, New York (1990). 7. C. E. Sheie, The golf club-ball collision-50,OOOg’s. Science and Golf: Proc. 1st World Scientific Congress of Golf, pp. 1999204. E. & F. N. Spon, New York-(1990j. 8. K. H. Hunt and F. R. E. Crosslev. Coefficient of restitution interpreted as damping -in vibroimpact. ASME J. appl. Mech. 42, 440445 (1975). 9. N. Maw, J. R. Barber and J. N. Fawcett, The oblique impact of elastic spheres. Wear 38, 101-I 14 (1976).

Optimization of clubhead loft for maximum distance golf drive 10. T. R. Kane and D. A. Levison, Dynamics: Theory and Applications, pp. 37-159. McGraw-Hill, New York (1985). Il. R. M. Brach, Tangential restitution in collisions. Computational techni&es for impact. In Computational Techniques for Contact, Impact, Penetration and Perforation of Solids, pp. l-8. ASME, New York (1989).

12. S. Aoyama, A modern method for the measurement of aerodynamic lift and drag on golf balls.

25

Science and Golf: Proc. 1s1 World Scientific Congress

of Go/f, pp. 199-204. E. C F. N. Span, New York (1990). 13. B. Lieberman, Estimating lift and drag coefficients from golf ball trajectories. S&nce and Golf: Proc. 1st World Scientific Congress of Golf, pp. 187-198. E. &. F. N. Spon, New York (1990). 14. Numerical Algorithms Group Workstation Library Manula. Grove Park, Illinois (1991).