End closure shapes for filament wound pressure vessels

C o l l e g e of A e r o n a u t i c s Memo 7605 May 1976

TECHNISCHE HOGESCHOOL DELFT LUCHTVAART- EN RUIMTEVAARHECHNIEK

BIBLIOTHEEK

Kluyverweg 1 - DELFT

7 r

CRANFIELD INSTITUTE OF TECHNOLOGY

COLLEGE OF AERONAUTICS

END CLOSURE SHAPES FOR FILAMENT WOUND PRESSURE VESSELS

R.C. FEWS R. TETLOW

CRANFIELD INSTITUTE OF TECHNOLOGY

END CLOSURE SHAPES FOR FILAMENT WOUND PRESSURE VESSELS

by

R.C. FEWS MSc, CEng, MIMechE R. TETLOW MSc, CEng, MRAeS

SUMMARY

A computer programme is developed to define efficient filament wound pressure vessel end closure contours, for both geodesic and planar methods of winding,

Some general results are derived and a full description of theory, together with worked examples and programme list, are given in the report.

CONTENTS Page NOTATION 1. INTRODUCTION 1 1.1 Background 1 1.2 Design Assumptions 1 1.3 General Equation 1 l.h Geodesic Winding 2

1.5 Zero Hoop Stress Head Contour 2

1.6 Planar Winding 2 1.7 Other Head Contours 3

2. COMPUTER PROGRAMME 3

2.1 General 3 2.2 Description of Input Data 4

2.3 Programme Output Data 4 2.4 Some General Results 4

3. DISCUSSION 5 3.1 Programme Verification 5 3.2 Programme Limitations 5 3.3 Programme Uses 5 h, REFERENCES 6 APPENDICES A THEORETICAL FUNDAMENTALS

1. Geodesic Wind Pattern 2. Planar Wind Pattern

3. Helix Angle Relationships 4. End Closure Contours

5. Fibre Helix Angle 6. Zero Hoop Stress Head

NOTATION

D Internal diameter of pressure vessel L Overall length of pressure vessel N Stress resultant

P Internal pressure

a Internal radius of cylinder I'-i Boss radius at end 1

Tp Boss radius at end 2 t Total shell th:i cknoss

^1

Normalised coordinates z

-x^ Normalised boss radius at end 1, ie —

o a a Fibre helix angle

(3 1 Fibre and meridian angles of inclination 5 _J (see Fig.5)

a Fibre stress

SUBSCRIPTS

^ Parameter referred to meridional direction 6 Parameter referred to hoop direction

1. INTRODUCTION 1.1 Background

The current programme of work was initiated to establish stress analysis techniques for application to filament wound pressure vessels.

The Ministry of Defence (P.E.) are currently interested in the development of a particular filament wound rocket

motor case. The motor is to be designed using the best theoretical knowledge available and specimens are to be

manufactured for structural test and eventual direct comparison with components fabricated in more conventional materials.

The role of Cranfield Institute of Technology in this procedure is to monitor structural design and investigate general stress analysis methods, as applied to filament wound pressure vessels.

The basic philosophy adopted in this task is to package the structure into selective design areas, identify suitable stress analysis procedure and develop computer programmes to manipulate the various theoretical relationships involved, Pressure vessel end closures form an obvious specialised

package in this work programme and relationships pertinent to the evaluation of efficient filament wound end closure shapes are discussed in this report. A computer programme is defined to identify desirable end closure contours, for both geodesic and planar winding procedures,

1.2 Design Assumptions

In the defination of end closure contours it is assumed that the end closure is to be integrally wound with the

cylindrical portion of the pressure vessel, as indicated in Figure 1. End closure plies are assumed to be helical in nature, with as many plies orientated at -a° to a meridian as at +a. Emphasis is to be placed on the optimum end closure contour and not on its precise state of stress. Netting

analysis assumptions are therefore employed and hence, all contours are evaluated to give zero stress at right angles to the fibres,

1.3 General Equation

Using the relationships for the radii of curvature of a surface of revolution and the condition of equilibrium under internal pressure loading, the following equation can be derived

P ^ ö=i = 2 - tan^a

y'[i + (y')^] ...(1)

See Appendix A

2

-Equation (1) can only be solved if the helix angle (a) can be expressed in terms of the normalised coordinates x and y. This demands a knowledge of the winding method to be employed and in this connection the following are the most commonly encountered in industry.

1.4 Geodesic Winding

With this method fibres are laid on geodesic paths, which on a surface of revolution are defined

by:-X sin a = const ,,,(2) A geodesic path is the shortest distance between two

points on a surface and is the only path along which fibres can be laid, without the tendency to slip under the winding tension (see Figure l ) .

There is often a requirement for polar bosses in the manufacture of filament wound pressure vessels and since the filaments must be tangent to these bosses the constant in equation 2 is defined as X Q

ie:-sin a = ^o/x ...(3) Substituting equation (3) into (l) and solving the

differential equation, results in head contours; which in view of the wind pattern, (equation 3) are normally termed Geodesic Isotensoid Head Contours,

1.5 Zero Hoop Stress Head Contour

If the polar boss normalised radius X Q is zero, the fibre paths follow meridians and the helix angle (a) is

zero over the whole extent of the end closure. When equation (1) is satifled, there can be no stress at right angles to the fibres and for the zero boss radius case, this direction is coincident with the hoop direction. Consequently geodesic contours derived for zero boss radii, are often termed Zero Hoop Stress Head Contours.

1.6 Planar Winding

Geodesic winding machines require fairly complex gearing and as a consequence? even though they produce the most

efficient end closures, they are rarely used in practice, A more popular winding method employs rotating arm type machinery, which lays the fibre in a longitudinal plane tangent to the polar bosses (See Figure 2 ) ,

For this winding procedure, head contour definition is not as straight forward as in the case of the geodesic wind patterns discussed above. The helix angle relationship in terms of the normalised coordinates (x) and (y), can now only be established if the head contour is known and there

is therefore no explicit solution to equation (1). Acceptable solutions can however be derived, by assuming an initial head contour (ie geodesic), evaluating the helix angle relationship on this contour and substituting back into equation (l), to give a modified solution. The helix angle relationship can now be re-evaluated and further iterations carried out,

until convergence to the desired planar contour is achieved, 1.7 Other Head Contours

The head contours discussed above are derived from an internal pressuration loading environment only and any

alternative, or additional loading will result in different optimum shapes.

Most practical pressure vessel designs utilise a low helix angle planar wrap to form the end closures, making good the strength in the parallel region by additional hoop windings. There is usually therefore, a difference in

shell thickness at the junction of end closure and cylinder, and the resulting discontinuity gives rise to both bending and transverse shear forces. These discontinuity forces are dependent on the local shell curvatures and using continuum analysis, it is possible to vary the local head contour to minimise these effects.

2, COMPUTER PROGRAMME 2,1 General

The computer programme derives both geodesic and planar head contours, (where applicable) from a step by step

numerical solution of equation (l). The datum location is the equator of the end closure and normalised radial

coordinates are defined for 0.02 increments in axial

normalised coordinate, working from an initial axial value of 0.0 (equator). For all but the zero hoop stress head, a stationary point in the contour defining, differential

equation (equation l) exists and output is terminated at or around this location. The normalised axial coordinate at the stationary point is of order 0,6 and hence the number of coordinate points generated is approximately 30, For axial coordinates in excess of 0,58, the step size is reduced from 0.02 to 0,005, as variations of radial coordinate are much more sensitive to axial increments in this region and a step of 0.02 may cause data generation to progress beyond the stationary point.

The programme always generates the geodesic head

contour first and then, according to an input marker, iterates on this contour using the planar helix angle relationship,

to thereby generate the planar head contour. When planar head contours are requested, both geodesic and planar results are output by the programme.

4

-2.2 Description of Input Data

Input data to the programme is punched on a single card, containing six real variations

ie:-Rl R2 R AL TO RLIM

Rl is the polar boss radius at end 1 of the pressure vessel.

R2 is the polar boss radius at end 2 of the pressure vessel.

(N.B, contours output always refer to end l)

R is the internal radius of the parallel portion of the pressure vessel.

AL is the parallel length of the pressure vessel. If AL is punched is 0,0, then only the geodesic contour will be evaluated,

TO is the end closure thickness at its equator (end 1) RLIM is the limiting minimum radial coordinate. If RLIM

is greater than 1.2 x Rl, then data generation will cease at RLIM. If RLIM is less than 1.2 x Rl data generation will cease at 1.2 x Rl.

2.3 Programme Output Data

The programme outputs radial and axial coordinates in both normalised and natural form, together with respective helix angles and shell thicknesses. If the parallel length of the vessel is input as 0,0, then output is limited to the geodesic contourj otherwise tables referring to the planar, as well as the geodesic contour are output,

(see the worked examples of Appendix B for actual programme output).

2.4 Some General Results

In Figure 3^ Geodesic and Planar helix angles are compared from equator to polar boss, for a pressure vessel having an X Q value of 0.5 and an L/D ratio of 2, It can be seen that significant differences in helix angle occur, especially in the equator regions, where the planar helix angle is approximately half its geodesic counterpart. These helix angle relationships feature strongly in the head

contour defining equations and are responsible for the

variations in the optimum head shapes achieved from the two winding methods.

In Figure 4 programme derived geodesic and planar head contours are plotted for boss to cylinder rations of 0,2, 0,5 and 0.7. Prom this figure it can be seen that the planar contours are generally less domed than their geodesic

equivalents, although the former will be affected by variations in the L/D ratio of the vessel to be wound,

3. DISCUSSION

3.1 Programme Verification

Geodesic and Planar contours for the three boss to cylinder ratios plotted in Figure 4 have been published in Ref,2. Comparison of these results with the Figure 4

information, revealed no significant differences, 3.2 Programme Limitations

Results for zero hoop stress heads cannot be derived in the immediate vicinity of the polar boss. This is a consequence of the numerical solution procedure used and results from the fact that a discontinuity appears in the

'Runge Kutta' relationship at x = 0, (see Appendix A ) , For the derivation of these contours it is suggested that RLIM (see 2,2) should not be less than 0.2 x R.

3.3 Programme Uses

Integrally wound end closures must always conform to the helical ply requirements of any parallel portion that may

exist in the completed pressure vessel. It is therefore often necessary to design end closures conforming to given equatorial helix angles, as well as given polar boss diameters. It is likely therefore that end closure contours will be different for each new pressure vessel design and they may well be

different at each end of a pressure vessel. Manual definition of these contours can be extremely time consuming and future application of the programme to this task, should therefore result in very real savings in both time and effort,

In addition to its independent use the programme forms an important link in the chain of computer packages currently being developed at Cranfield. In this respect the logic is directly applicable to automatic data generation procedures for the proposed finite element stress analysis programme,

6

-4. REFERENCES

1. ZICKEL, J, Isotensoid Pressure Vessels. A.R.S. Journal Vol,32 1963.

2, DARMS, F,J, Optimum Design for Filament Wound LITVAK, J, Rocket Motor Cases,

S,P,I, Reinforce Plastics Division 6-D Annual Conference Vol.19 1964.

c\

Y

Meridions

Path Of Fibre

FIG 2 PLANAR WIND PATTERN

h

Geodesie Planar 90 80 70 60 I-50 AO Ci < 3 0 dl X 2 0 10 -_L _L X 01 0-2 0-3 OA 0 5 0 6 07 0 8 0-9 1-0 Normalised Radial Co-Ordinate—w~

10 -GEODESIC CONTO U RS 0-7

I 06

<u 7d £ 0 5

"2

o

S 0 A

<

0-3

I 0-2

E 1— o ^ z o i 0-1 0-2 0 3 OA 0 5 0-6 0 7 0-8 Normalised Radial Co-Ordinate PLANAR CONTOURS X • ±_s 2 0 7 1 1 0-6 0» <-• | 0 5 Ó Ö o G-A < 0-3 Norma l ise d p p 1. k > — _ -~ -J . 01 J -02 OA 05 0 6 0-7 0 8 0 9

Normalised Radial Co-Ordinate — ^ 10

12 -AO | 3 0 (b Q .5 X < 20 10

1^^

1 20 30 Radial Distance (mm)—*~ AO 50

APPENDIX A

TPÏEORETICAL FUNDAMENTALS

1 . GENERAL DIFFERENTIAL EQUATION

1.1 Assumptions

In the following theoretical relationships no strength is attributed to the resin (netting analysis) and the fibres are assumed to be loaded to identical stress levels.

1.2 Curvature Equations

The radii of curvature of a surface of revolution in terms of normalised axial and radial coordinates can be defined by

a y*^

...(1)

^9 _. - x L l + (yM^J^

a y'

where the primes indicate differentiation with respect to x. 1.3 Equilibrium Equations "«; % "«( a) Internal pressure

2

2

b) Fibre loads = a^ t cos a 2 ...(2) 2 N Q = a^ t sin a From 2 and 5 r

...(3)

2 -

-i

= tan^ a

^/rf ...(4)

A2

-Substituting for VQ and r, from (1) into (4) gives

r ^ y " 5-=r- = 2 - tan^a

y' [ 1 + (y')^] ...(5)

(see Ref.l for original definition of this equation)

2. HELIX ANGLE EVALUATION 2.1 Geodesic Winding

Before equation (5) can be solved, the helix angle (a) must be expressed in terms of the normalised coordinates x and y. When the fibres follow geodesic paths this task is trivial, as a geodesic path is by

definition:-X sin a = CONST .. .(6) Where finite polar bosses are used in the winding procedure, the constant in (6) can be evaluated from the fact that the fibre must be tangent to the boss

iei-X

sin a = rr^ (See Figure l) 2.2 Planar Winding

The planar winding technique consists of winding fibres in a longitudinal plane tangent to polar bosses (see Figure 2 ) . In defining the helix angle relationship for this wind pattern, use is made of the diagram shown in Figure (5). This diagram

shov;s the tangents to a fibre and its intersecting meridian, at a general point on the end closure surface (x,y,z). The angles of inclination of these two tangents are (5-|) and O2.) for the meridian and (62) and (P2) for the fibre (see Figure 5 ) .

Assuming unit projections of cp and ep in the zx plane.

^^ " cos Sg '' ^^ "" cos 6j^ From triangle abp

P. + P p \

ab = 2 sin < g —- > = cd

From triangles bep and acp

ce =

From triangle ode

2

(H

^ ^2^

+ (tan 5-, - tan 5p)

1

2

and therefore using the Cosine Rule on triangle cep the helix angle is defined from

cos a ₂ cos a-, 2 cos a, 2 ^1 + ^2 4 sin -^ '^ cos Op cos a,

^ —- + (tan a - tan Op)

...(7)

Equation (7) gives the helix angle (a) in terms of the angles of inclination of the tangents to the paths of fibre and meridian, at a point of intersection on the end closure surface. The angles 5 and p must now be expressed in terms of the cartesian coordinate system (x,y,z) and this can be achieved as follows.

A meridian is defined as

^ = Const (see Figure 2)

therefore X = r cos /

y = r sin ^

...(8)

Before further progress can be made, it must be assumed that the head contour is known, so that a third equation can be added to (8)

iel-r = f(z)

The path of the fibre is known, since it is wound in a longitudinal plane tangent to the polar bosses and for the fibre therefore

X = z tan a,

...(9)

Where the origin of (z) is at the intersection of the fibre path with the axis of the vessel and the angle O Q is the helix angle of the fibre over the parallel portion of the vessel, (see Figure 2 ) . The angle ( O Q ) can be calculated from the overall vessel length and the diameter of the polar bosses.

A4

-At any point of intersection between fibre and meridian z tan a = r cos i6 o '^ therefore = cos T ] z tan a ~i y o tan

From Figure 5 it can be seen that

meridian ...(10)

•• * - P i = l f ^ - = - ^ l l

. *. p^ = tan Similarly -1 cos /$ ^ ...(11) fibre = tan a. •'• ^2 " % _...(12)

Also from Figure 5 it can be seen that

tan 5-,

meridian

where z is defined from triangle fbp

K Z le ZT = '1 cos 3,

öy or ÖZ Therefore tan 6. =( T ^ § ^

^z' _meridian

5, = tan

Similarly

-1 _sin

f> cos 3

_{1 ^J}

ör~|

...(13)

5p = tan sin / cos 3p ^

...(14) Substitution of equations (10 to 14) into (7) gives an expression for the helix angle (a), at a general hoop circle (r), in terras of r, z and the local slope of the end closure surface at (r)j ie èr

-^

.*. a = f(r,z)

or in terms of the normalised coordinates used to define the geodesic contour a = g(x,y) where x = -• a z - z (see Figure 2) ...(15)

3. SOLUTION OF DIFFERENTIAL EQUATION 3.1 Geodesic Solution

Prom equations 5 and 6

y" = y'[i + (y')^]

2 2~ 2x2 - 3 x / I 2 2N X(X - X Q ) ie y " = f(x,v) where v = y' Initial Values

When X = i; y = 0 and y' — 'i^

...(16)

Due to the initial value y' = -^^ it is not possible to solve equation (l6) in its original form.

7^

dx Therefore rewriting (l6) 2' dx 1 +

(ST

f(x) ...(l6a)

A6

-...(17) then

ifz = d _ l , _ l _ d v ^ l d ^ = 1 dv

^ dx V • ^2 dx " 7? dy " ^ dy

.*. substituting into l6a

, dx where v = -^

Initial Values

dx When X = 1; y = O and J^ = O

Equations (17) can now be.solved using a Runge Kutta

method and the programme employs a fourth order routine, giving accuracy of the order of the step size to the power of five. 3.2 Planar Solution

From equations (5) and (15)

y " - y ' [ l + Cy')^]i[2- tan2{g(x,y)}] • ___(,g)

As has been mentioned previously the function g(x,y), can only be evaluated for a known head contour. The geodesic solution is therefore first evaluated (see 3.1) and the

function g(x,y) is derived from this solution (see 2.2). Equation (l8) is then solved in a similar manner to the

geodesic solution discussed above, to give a new head shape. The function g(x,y) is then defined from this contour and the

process repeated until convergence on the true planar contour is achieved.

4. SHELL THICKNESS VARIATIONS

Due to the variation of fibre helix angle and the fact that fibres cross hoop circles of decreasing radius as they approach the boss, there will be a tendency for the shell to thicken towards the poles.

This thickness variation is given in Ref.l as the following

X cos a ...(19) Where to is the shell thickness at the equator of the

end closure.

Equation (19) is applicable to both geodesic and planar wound end closures, providing the appropriate helix angle

BI

-APPENDIX B

WORKED EXAMPLES AND PROGRAMME LIST

1. WORKED EXAMPLE 1 1.1 Problem Description

It is desired to identify a Zero Hoop Stress head contour, to suit a pressure vessel of 90 mm diameter, with a minimum wall thickness of 1 mm.

1.2 Input Data

O B O D f S I C CONTO'JR N O l ? M A L I s t O C < ' « O R f i I N A T E S RADIAL OtSTANCE 0. 0. 0, 0. 0. 0, 0. 0. 0. 0. 0, 0. 0. 0. 0, 0. 0. 0. 0. 0. 0. 0, 0. 0. 0. 0. 0. 0. 0. 0. 0. ,1.OOOE .99960E .99840E ,9«fJ9E ,99357E .98992E ,98542E ,98007E ,97383t ,96668E ,95858E ,94949E ,93936e ,92814E ,91576E ,9(>213E ,88716E ,87073E ,85270E ,8S290E ,81109E .78700E ,7«0?SE ,73033E ,«9653E ,657816 ,61249E ,557«9E ,48743É ,38515E ,13710E f.1 OO 00 00 00 00 00 00 0 0 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 OO no 00 0 0 A X I A L n i S T A N C E 0, 0, 0. 0, 0. 0. 0, 0. 0. 0. 0. 0. 0. 0, 0. 0. 0. 0. 0, 0. 0. 0. 0, 0, 0. 0, 0. 0. 0. 0, 0, . l O O O O E 0" ..^OOOOE-01 ,40000?' • 0"^ ,600006-01 , ^iOOOOE-0' ,10000? ,12000E ,14000E ,16000E ,18000E ,200001 ,^2000E ,?4000E ,26000E ,i?8000E ,30000E ,32000E ,34000E , 3 6 0 0 0 E ,38000E ,40000E .42000E ,44000E ,46000E ,'•80005 .SOOOOE ,52000E ,340001 ,''>600'ÖE .580001 ,S8500E 0»' 0' Oo 0" OO 0»' 00 0" 00 0') 00 OO Ou 0.1 00 00 on Ou 00 00 00 OC 00 00 00 00

13 Continued - Ü ; . -NATURAL Cr RACIAL OISTAWCf 0.45000E 0.44982K 0.44928E 0.44838E 0.44710Ê 0.44546Ë 0.44344E 0.44103E 0.43822E 0.435niE 0.431S6E 0.4?727E 0.4?271E 0.4i7*6F 0.4i209E 0.4 506f; 0.3'^922E 0.391R3E 0»3R372E 0.374«0t^ 0.SA499E 0.354l5fc 0.34211E 0.3?8A5E 0.51344E 0.296niK 0.275A2E 0.25096E 0,2i934E 0.17332t 0.61694E 1-0,, 0 2 02 02 0 2 02 0 2 02 02 02 02 02 02 02 02 02 02 02 0 2 0 2 02 02 0 2 02 02 0 2 02 02 02 02 02 01 O I N A T E S AXIAL Dl'JTAMCE O.OOOOOE O.UOOOOE 0.18000E o.^7oooe 0.36000E 0.45000E 0.54000E 0,63000E 0.72000E O.olOOOE 0.90000B 0.99000E 0.108006 0.';1700E 0.12A00B 0.13500E Ö.14400E O.-iSSoOE 0. IA20ÜE 0. I7100E o.i8oooe 0.T8900E 0.19800Ê 0.20700E 0.i:l600E 0.*;2500E o.;;34ooE O.;H3OOE 0.25200E O.?6100E 0.26325E 00

or.

0-i 01 0" 0-, 01 0' 01 0' 0'' 0-02 02 02 07. 02 0^ 0.'. 0^; Od 02 02 02

oa

o:;

0?.

o^

02 02 02 H E L I X ANGLE (DEC) Ó.OOO0OE 0.000OOÊ O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE

o . o o o o a E

O.OOOOOE O.OOOOOB O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOB O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE O.OOOOOE '>0 00 00 00 00 00 00 00 00 00 00

no

00 00 00 0 0 00 00 00 00 00 00 00 0 0 00 00 00 00 00 00 00 S H E L L T H I C K N E S S O . I O O O O E 0.100041 0.100161 0,10036E 0.100651 0.101021 0.101*81 o.to2nSE 0.102691 0,103451 0.104321 0.10532B 0.106461 0.107741 0.109201 0.11085E 0.112721 0.114851 0.117271 0.120061 Ö.tt329f 0.1270AI 0.1S154t 0.1569 31 0.143571 0.18:>02e 0,16527f 0.179318 0.20516f 0.259641 0.259641 » 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01 01

1.3 Continued

The contour defined in the above output is plotted in Figure 6.

2. WORKED EXAMPLE 2 2.1 Problem Description

It is required to establish the end closure shapes for a planar wound pressure vessel, having a maximum disuneter of 90 mm with 45 mm diameter polar bosses. The minimum shell thickness is to be 1 mm and the vessel is to have a parallel length of 126 mm. In addition to the definition of end

closure contours, the variations of both shell thickness and fibre helix angle are to be established from equator to pole, 2.2 Input Data

2.3 Programme Output GEODESIC COi^iOiiR W0pM4LIsED C'i-ORJ.INATES RADIAL DISTANCE O.luOOÖÊ 0 . 9 O 9 6 7 E 0 . 9 O 8 * 7 F 0 . 9 O 6 9 9 E 0 . 9 9 4 6 5 E 0 . 9 0 1 6 2 1 . 0 . 9 « 7 9 i r : 0 . 9 « 3 4 9 F 0 . 9 7 8 3 7 E 0 . 9 7 2 5 3 f : 0 . 9 6 5 9 4 E 0 . 9 5 8 5 9 r 0 . 9 5 0 4 7 F ; 0 . 9 4 l 5 4 r : 0 . 9 3 1 7 7 E 0 . 9 2 1 1 5 f : 0 . 9 ; 9 6 2 f c 0 . 8 0 7 1 6 K 0 . 8 r < 3 7 2 f 0.669Z%i: 0 . 8 5 3 7 0 E 0 . 8 3 7 0 0 ^ : 0 . 8 l 9 1 0 i : 0 . 7 0 9 9 2 . -0 . 7 7 9 3 9 E 0 . 7 ? 7 4 $ E 0 . 7 3 4 0 2 E 0 . 7 ; . 9 0 8 f 0 . 6 « 2 6 3 P 0 . 6 S 4 7 9 i 0 . 6 2 5 8 / F 0 . 6 ) 8 5 4 ! 0 . 6 > 1 1 9 r 0.6».3851-0 . 5 9 6 5 3 f c 01 00 OO on «0 i'O 00 t»0 0 0 ttO (10 OO 0 0 0 0 0 0 00 0 0 0 0 ' 0 0 0 00 0 0 UO | i O OO 00 0 0 0 0 00 UO 00 0 0 ('0 00 0 0 ^XIAL n i S T A N C E 0 0 ^ (\, ( • ; ^ o! (' , 0 , 1 < ^ ( l . <)', i ; , (, _ ' 1 . 0 , 0 . 1*1 f ; , ••( . (> f; , (• . 0 . f) ^ I l ^ ( / 0 , 1 ; , 1' , (\ (' , I I ft , (1 ^ (! [ .'OOOOE . 2 0 0 0 0 E ' . 4 0 0 0 0 E ' .ÓOOOOE' ,:iooooE' . OOOOE . ,2O00E . . 4 0 0 0 E . ; 6 0 0 0 E . . 8 0 0 0 F , 2 0 0 0 0 E 2 2 0 0 0 F , 2 4 0 0 0 E ,,^6000F , r^SOOOF , ÏOOOOF ,:52000E .;4000E 3 6 0 0 0 F , - ' 8 0 0 0 F *.0O00F . ' • 2 0 0 0 F ^ 4 0 0 0 E .*6000F ' • 8 0 0 0 6 •OOOOE S 2 0 0 0 F ;>4000F •6O00F 8O00E S8500E ^ooo*^ > 9 5 0 0 E '•OOOOE , 0 5 0 0 E 0 ; j -O-- 0 1 • 0 ' - 0 ' Ow 0.; 0 0 ' 0' O v 0 Oi 0 • 0'. 0'. 0 0' 0 0 : 0 0 0 0. 0.; 0 • 0 ' 0 0 • 0 0 -0' 0 0 • Ö i '

N A T U R A L f O - O K D I N ^ T E S R A D I A L O I S T A N C K 0.45goor ui 0.4498S!. 02 0.4A940^- 02 0.44865^ 02 0.4475«?t i,p, 0 . 4 ^ 6 2 3? (.2 0.4/*4'56f i.2 0 . 4 / . 2 5 7 ! 02 0 . 4 4 0 2 7 r 02 0 . 4 T 7 6 4 I 02 0 . 4 3 4 6 ? i 02 0 . 4 ^ 1 3 7 * (,2 0 . 4 / 7 7 1 : 02 0 . 4 ; ; 369 f u? 0 . 4 | 9 5 0 ' - 02 0 . 4 i 4 5 2 , 02 0 . 4 0 3 3 ; 02 0 . 4 372^ (12 0 . 3 O 7 6 « i ; 02 0.39111SF: 02 0 . 3 8 4 1 6 ' ;i2 0 . 3 7 6 A i f 02 0 . 3 6 8 5 9 ) 02 0 . 3 5 9 9 6 f 02 0 . 3 5 0 7 3 r o2 0 . 3 4 o a 5 r ..2 0 . 3 3 0 3 1 . 0? 0 . 3 1 9 0 9 ' 02 0 . 3 71«» /i2 0 . 2 9 4 6 6 ) . M2 0 . ? f l 6 4 r 02 0 . 2 7 8 3 4 » 02 0 . ? 7 5 0 4 f 02 0 . 2 7 1 7 3 ! ; 02 0.?'^844< 02 01 0. 0. 0. '5. 0. 0. 0. 0. 0. 0. f . 0. 0. 0. ". 0. 0, 0. 0 . 0. n ^ • 0. < > . , 0 . 'i',' 'K> 0., 0.. 0.. 0.. f ' . . '•'.' ' > . < ^*XjAL STANCE •»OOOOF 0 ; •'OOOOE 0 ^«'lOOF 0 .^'7000F 0 .')6000F 0 450Ö0F 0. :-4000F a ('3000F 0 "'2000F 0 «'lOOOF 0 'OOOOF 0 ^•9000F 0-. oaooF 0, M ^ O O f 0. .2600F 0 •^SOOF Oi. 440ÖP Or 5300F 0/ •6200F 0,: 71 OOF Or^ «OOOF 0 «0 00F 0, 9«00P 0. i!0 70OF 0., 1600F 0. .250OF ;),• .3400F 0 : '.4300P 0/ 520ÖF iV-.61 OOF 0 6325F 0 :6fi50F 0, ; 6 7 7 5 F 0 '7000F 0,: 722SF 0, H l u l X ANQLF (OEG) 0 . 3 0 0 0 0 F 2 0 . 3 0 0 1 1 E v»2 0 , 3 0 0 4 4 F ^'2 0 3 0 1 0 0 F J 2 0'.3017ÖF 2 0 . 3 0 2 8 0 F ,'2 0 . 3 0 4 0 6 F -.2 0 . 3 O 5 5 7 R v2 0 . 3 0 7 3 4 F - 2 0 . 3 0 9 3 9 F 12 0 . 3 1 1 7 4 F . 2 0 . 3 1 4 4 0 E n2 0 . 3 1 7 3 9 F .'? 0 . 3 2 0 7 6 E M 2 0 . 3 2 4 5 3 F )? 0.3287bE ;2 0.3334»SF '2 0.3387UF ..2 0,3445i'F '2 Ö.35114F i2 0 . 3 5 8 5 2 E -12 0.36682F 2 0.37620F '.2 0.386«7F 2 0.3O906F .^2 0.41308F ;2 0 42936F t? 0.44841F )2 O.A7ö'^'3E J? 0.49783F 02 0.53024F 1.2 0.53936E o2 0.548O3F .,2 C . 5 5 8 0 6 E )2 0 . 5 6 9 4 « F 12 S„FLL THICKNFSS O.IOOOOE 01 0 . 1 0 0 0 4 e 01 0 . 1 0 0 1 8 E 01 u . i o o 4 o e 01 0 . 1 0 0 7 2 E 01 " . 1 0 1 1 3 1 01 0 . 1 0 1 6 4 E 01 0.10;?26E 01 0.10?98E 01 0,10^82B 01 0 . 1 0 4 7 0 E 01 0.10.^«9E 01 O.10714E 01 o . l o ^ S S e 01 o , 1 i o i 5 E 01 0 . 1 1 1 9 4 E 01 " . 1 1 3 9 7 E 01 0 . 1 1 6 2 6 E 01 0.11fi85E 01 0 . 1 2 1 7 0 E 01 0 . 1 2 5 1 6 E 01 0.12'"'»02E 01 0 . 1 3 3 4 8 E 01 0 , 1 S « 7 0 e 01 0 . 1 4 4 8 5 E nl 0.15?21E 01 0 . 1 6 1 1 5 E 01 0 . 1 7 ? 2 5 e 01 0 . 1 8 6 3 5 E 01 0 . 2 0 4 8 4 E 01 0 . 2 5 0 0 5 E 01 0.2S''84E 01 0 . 2 4 6 3 8 E 01 0.25f.79E 01 0 . 2 6 2 9 7 E 01

B7 -2.3 Cont B A L A W C E D IN PLANS H Ê A O C 'NTO'IR N O R M A L I S E D Ct'-OH - I W A T E S RADIAL OISTA^iCP O.li/OOOK 0 . 9 9 9 6 1 f ; 0 . 9 9 8 4 5 K 0 . 9 O 6 5 1 f ! O.9038OC 0 . 9 O 0 3 0 J : 0 . 9 A 6 0 0 F 0 . 9 8 0 9 0 f : 0 . 9 7 4 9 8 F 0 . 9 6 8 2 1 F : 0 . 9 / i 0 5 9 ( . 0 . 9 S 2 0 7 f ; 0 . 9 A 2 6 4 J 0 . 9 3 2 2 5 E 0 . 9 2 0 8 7 r 0 . 9 . 846F. 0 . 8 0 4 0 5 r 0 . 8 « 0 3 0 t 0 . 8 6 4 4 4 C 0 . 8 4 7 2 8 e 0 . 8 ? 8 7 6 l 0 . 8 . 8 7 7 r 0 . 7 « 7 2 2 r o.764on: 0 . 7 3 9 ü i r 0 . 7 1 21 7r: 0 . 6 « 3 4 4 F 0.652'<»0F 0 . 6 2 0 F 9 F 0 . 5 « 7 9 5 i : 01 00 0 0 " 0 Of) 0 0 •)0 00 oO 00 0 0 00 00 ('C dO 0 0 00 0 0 00 0 0 t i O 0 0 ('0 0 0 0 0 00 00 00 0 0 0 0 0] 0 , 0 . 0 . 0 . 1) ^ 0 0 h 0 , 0 ' 0 ('i 0 0 0 it f l 11 0 0 (> r 0 f 1 M (1 (1 M X U L STANCF 'OOOOF c'OOOOE4 0 0 0 0 E - «.OOOOE-Ö •0 •0 •0 .OOOOF-O 'OOOOE 0 2 0 0 0 F 0' ;400ÖE 16OOOF 8 0 0 0 F ^OOOOF , : 2 0 0 0 F ^ 4 0 0 0 F ;^6oooE ;.'8oooF . 3 0 0 0 0 F 3 2 0 0 0 F 3 4 0 0 0 E -;6000F 38Ü00E > 0 0 0 0 F , ' • 2 0 0 0 F .'•J^OOOE 4 6 0 0 0 F . 8 0 0 0 F .SOOOOE '•.2000 F .'^OOOE . : . 6 0 0 0 F . •.>8000E 0 0 0' 0 0 0 0 0 0 0' 0 0 0 0' 0 0 0' 0 0 0 0 0 0'

KjAT'JRAL 1 ' < j - 0 , , n i N ' , T f S a A ^ I A L O I S T * K ' C i : 0 , 4 ^ 0 0 1) r 0 . 4 4 9 « i r 0 . 4 4 9 3 0 1 0 . 4 4 8 4 3 i . 0 . 4 4 7 ? 1 r 0 . 4 4 5 6 i K 0 . 4 4 3 7 0 : 0 . 4 4 1 4 1 . 0 . 4 3 8 7 4 0 . 4 3 5 / ' O . : 0 . 4 3 2 2 ^ ^ 0 . 4 , ' 8 4 3 ! : 0 . 4 ? 4 l O e 0 . 4 i 9 5 l ' : 0 . 4 i 4 3 9 ' 0 . 4 , 8 X 0 - . 0 . 4 273.^ 0 . 3 0 6 1 3 ^ 0 . 3 8 9 0 0 i 0 . 3 « l 2 ö i : 0 . 3 7 2 V 4 : 0 . 3 A 3 9 ' > ! 0 . 3 ' ? 4 ? 5 0 . 3 4 3 f t 0 i 0 . 3 ^ 2 S 6 f 0 . 3 < ^ n 4 ö r 0 . 3 7 5 5 r 0 . 2 0 3 ^ 1 / 0 . 2 7 9 t o ^ 0 . 2 6 4 S A 0 2 o 2 ü 2 0 2 f; 2 0 2 i»2 0 2 ••2 <i2 c 2 o 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 .)2 0 2 " 2 0 2 0 2 •i2 01 f 0 , 0.' 0 . 0 . 0 . e . 0 . ' ) . 0 . • 0 . n ^ 0.' 0 , 0 . 0 . 0 . ' \ 0 . ' • • • » 0 . 0 . 0 . fl • 0 . ^ . 0 . " . A X I A L 'sTANCF ' 01) 0 0 P ^'OOOOE i 8 0 0 0 F 2 7 0 0 0 E . ^ 6 0 0 0 F 4 5 0 0 0 F • . 4 0 0 0 F o 3 0 0 0 E . ' 2 0 0 0 F ''.lOOOF '«'OOOOF ' • - 9 0 0 0 F 0 8 0 0 F . 1 7 0 OF 2 6 0 0 F 3 5 0 0 E 4 4 0 Ü F • 5 3 0 0 F <S?0OF 7 1 0 0 E «OOOE • « 9 0 0 F vaooF ; : o 7 o o F : 1 6 0 0 F .•^2500F ; : 3 4 0 0 F ; 4 Ï 0 0 F ; 5 2 0 0 F 6 1 0 0 F A , 0 n . 0 0 0 : 0 O'. n-i O'l 0 0 ' 0 • 0 . , 0 0 . : 0. 0 ' Oi:. 0--: o;-0 / 0 , ' 0 : n.-O.-: 0 : 0. 0 . : 0 . . H E L I X A « 6 L H ( D E G ) 0 . 1 4 0 3 6 F 0 . 1 4 8 1 5 E O . I S S O ^ E 0 . 1 6 3 8 3 F 0 . 1 7 1 7 5 F 0 . 1 7 9 7 6 F 0 . 1 8 7 - S 7 h 0 . 1 9 6 0 9 F 0 . 2 0 4 4 e i F 0 . 2 1 3 0 Ü E 0 . 2 2 1 7 4 F 0 . 2 3 Q 7 1 g 0 . 2 3 9 9 5 E 0 . 2 4 9 5 1 E 0 , ? 5 9 4 3 F 0 . 2 6 9 7 9 E 0 . 2 Ö 0 Ó 4 F 0 . 2 9 2 1 Ü F 0 . 3 O 4 2 5 F 0 , S 1 7 2 4 F 0 . 3 3 1 2 4 F 0 . 3 4 6 4 6 E 0 . 3 6 3 1 ('E 0 , 3 a i 7 2 E 0 . 4 0 2 S 7 F 0 . 4 2 6 2 9 F 0 . 4 5 3 * ^ 2 ^ 0 . 4 « 3 4 O F 0 5 2 ? f ; 9 p 0 . 5 2 5 1 6 F •7 '.'2 0 2 2 '••2 0 2 •••^l 2 2 2 2 •''2 • 2 2 '•:2 '^2 •'2 0 2 '•'2 :->2 •)2 " 2 •;2 >2 = 2 . ' 2 't2 0 2 •!2 ^•2 S H I ^ L L T H I f K N F S S O.IOOOOE 01 0 . 1 O 0 3 9 E 01 0 , 1 0 0 8 a E 01 0 . 1 0 1 4 7 E 01 0 , 1 0 « M 8 E 01 0 . 1 0 2 9 0 8 01 0 . 1 0 3 9 3 B 01 0 . 1 0 4 9 9 B 01 0 . 1 0 6 1 9 6 01 0 . 1 0 7 5 5 B 01 0 . 1 0 ^ 0 6 6 01 Ó . 1 1 0 7 6 i 01 o ; 1 1 2 6 5 8 01 0 . 1 1 4 7 8 E 01 0 . 1 1 7 1 6 e 01 0.11<i-83E 01 0.12,-^85B 01 0 . 1 2 f t 2 6 e 01 ' • . 1 3 0 1 5 6 01 0 . 1 5 ^ 6 1 6 01 0 . 1 3 C ' 7 7 E 01 0 , 1 4 5 8 1 6 01 0 . 1 5 ? 9 5 E 01 0 . 1 6 1 S 2 E 01 0 . 1 ? r o 2 E 01 Ü . 1 8 ' ' . 1 5 6 01 0 . 2 0 ' ' 0 3 E 01 0 . 2 2 ^ - 4 2 6 01 0 . 2 5 ^ 2 8 6 01 0 . 2 7 1 1 5 6 01

- .i'J '

^•0 PROGRAMME OUTPUT

SEND TO(ED,PROGRAM FILE.STORE) PROQRAM(PTTI)

lNPUTl=CRn OUTPUT2=LP0

COMPRESS INTEGER AND LOGICAL END C C MASTER SEGHE^'T C MASTER DOMESHAPE DIMENSION X(in),Y(10),P(10) COMMON X D ( 1 0 n ) , Y D C i n 0 ) , A N G i ( 1 0 0 ) , T ( 1 0 0 ) 10 FORMAT ( 6 F 0 . 0 ) ?0 FORMAT ( 1 X , / 1 ( 1 X , E 1 2 . 5 ) )

30 FORMAT (//2X,16HQE0DES1C CONTOUR,/)

40 FORMAT (2X,23HN0RMALISED CO-ORDINATES,/)

50 FORMAT ( 4 X , 6 H R A D I A L , 8 X , 5 H A X I A L , / , 3 X , 8 H D I S T A N C E , 5 X , 8 H D I S T A M C E , / ) 60 FORMAT (/2X,20HNATUR'\L CO-ORDINATES,/)

70 FORMAT ( 4 X , 6 H R A P I A L , 8 X , 5 H A X I A L . 5 X , 1 1 H H E L I X A N Q L E , 5 X , 5 H S M E L L , / , 13X,8HD1ST»NCE,5X,8HDISTANCE,7X,5H(DEQ),6X,9HTHICKNESS,/)

80 FORMAT ( / / 2 X , 3 0 H B A L A r : r E D IN PLANE HEAD CONTOUR,/)

P I = 3 . 1 4 1 5 9 READ ( 1 , 1 0 ) R 1 , R 2 , R , A L , T 0 , R L I M RLlM = RLlrVR D I S = 1 . 2 * ^ 1 / R •• I F ( R L I M . L T . D I S ) RLIM=D1S Xn=R1/R NL = 0 1 NL=NL+1 YB1T=0.0? Y ( 2 ) = 0 . 0 X ( 2 ) = 1 . 0 P ( 2 ) = 0 . n -1 = 0 C

C NUMERICAL SOLUTION TO D . E . USINQ RUNGE KUTTA C ? 1=1+1 X ( 1 ) = X ( 2 ) P ( 1 ) = P ( 2 ) Y ( 1 ) = Y ( 2 ) X n < I ) = X ( 1 ) Y D ( I ) = Y ( 1 ) I F ( N L . G T . I ) GO TO 3 ANGd ) = A S 1 N ( X 0 / X ( 1 ) ) T ( l ) = ( T 0 » n o S ( A N G ( l ) ) ) / ( X ( l ) » C O S ( A N G ( l ) ) )

3 0 Cont

C

C SOLUTION WITH GEODESIC HELIX A N a E - R A D l U S RELATIONSHIP C AK0=P(1) A M O = F F 1 ( X ( 1 ) , P ( 1 ) , X O ) X X = X ( 1 ) + Y B I T / 2 . 0 » A K 0 PP=P(1)+YB1T/2,0»AM0 AK1=PP AM1=FF1(XX,PP,X0) X X = X ( 1 ) + Y B I T / 2 . 0 * A K 1 P P = P ( 1 ) + Y B l T * A M 1 / 2 . 0 AK2=PP AM2=FF1(XX,PP,X0) XX=X(1)-t-YBlT*AK2 PP=P(1)+YBIT»AM2 AK3=PP AM3=FF1(XX,PP,X0) GO TO 5 C

C SOLUTION WITH PLANAR HELIX A N a E - R A D l U S RELATIONSHIP C 3 AK0=P(1) A M 0 = F r 2 ( X ( 1 ) , P ( 1 ) , N P T S ) X X = X ( 1 ) + Y B I T / 2 . 0 « A K 0 p p = p ( l ) - f Y B l T / 2 . 0 » A M 0 AK1=PP AM1=FF2(XX,FP,NPTS) X X = X ( 1 ) + Y B I T / 2 . 0 » A K 1 F P = P ( 1 ) + Y B l T « A M 1 / 2 . 0 AK2=PP AM2=FF2(XX,PP,NPTS) XX=X(1)+YBIT»AK2 PP=P(1)+YBIT*AM2 AK3=PP AM3=FF?(XX,PP,NPTS) 5 X ( 2 ) = X ( l ) + Y P l T » ( A K 0 + 2 . 0 * A K 1 + 2 . 0 » A K 2 + A K 3 ) / 6 . 0 P ( 2 ) = P ( l ) + Y B l T * ( A M 0 + 2 . 0 « A M 1 + 2 . 0 * A M 2 + A M 3 ) / 6 . 0 I F ( Y ( 1 ) . Q T . 0 . 5 7 ) Y B I T = 0 . 0 0 5 Y ( 2 ) = Y ( 1 ) + Y B I T IF ( X ( 2 ) . L T . R L I M ) GO TO 6 GO TO 2 6 1=1+1 X D ( 1 ) = X ( 2 ) Y D ( I ) = Y ( 2 ) NPTS=1 IF ( N L . G T . 1 ) GO TO 9 l F ( X ( 2 ) . I . T . X O + 0 . 0 0 5 ) GO TO 2 7 A N G ( I ) = A S I N ( X 0 / X ( 2 ) ) T ( I ) = ( T O * C O S ( A N G ( 1 ) ) ) / ( X ( 1 ) « C O S ( A N G ( 1 ) ) )

Bi.l -«,0 TO 28 2 7 ANGd )=0.0 T(l)=0.0 2b WRITE (2,30) GO TO 12 15 WRITE (2,80) AL=0.0 12 WRITE (2,40) WRITE (2,50) DO 14 1=1,NPTS 14 WRITE (2,20) Xr:(l ), YOd ) WRITE (2,60) WRITE (2,70) DO 13 1=1,NPTS XNA=XD(I)«R YNA = Yn(I )*R HA=ANG(1)*180.0/PI 13 WRITE (2,20) XNA,YNA,HA,T(!) 9 IF(AL.LT.0.0001) GO TO 4 CALL HELIX ( R 1 , R 2 , R , A L , N P T S , T O ) ! F ( N L . G T . 6 ) GO TO 15 GO TO 1 4 STOP OK END

LINEAR INTERPOLATION ROUTINE

SUBROUTINE L I N I N T (NPTS,RP,AANG) COMMON R ( 1 0 0 ) , Z ( 1 0 0 ) , A N G ( 1 0 0 ) , T ( 1 0 0 ) DO 1 1=1,NPTS IF ( R ( I ) - R P ) 2, 3 , 1 2 A = R ( I - 1 ) - R ( I ) B = R ( 1 - 1 ) - R P C = A N G ( l ) - A r ' G ( l - l ) AANG=ANG(1-1)+C»B/A GO TO 4 3 AANG=ANG(i) GO TO 4 1 CONTINUE 4 RETURN END

PLANAR HELIX ANGLE-RADIUS FUNCTION

FUNCTION FF2 (X,P,NPTS) CALL L I N I t i T ( r i P T r ; , X , A L F } T 2 = T A N ( A I . F ) T2=T2»T? FX=2.n-T2 P P = ( 1 . 0 + P * P ) / X FF2=-PP*FX RETURN END

? 0 Cont

C GFOOESIC HELIX ANGLE-RADIUS FUNCTION C FUNCTIE" F F K X . P . A O ) F X = ( 2 . 0 * X » X - 3 . 0 » X O » X O ) / ( X * X - X O * X O ) P P = ( 1 . ' ^ - i - P » P ) / X FF1=-PP»FX RETURN END C

C EVALUATION OF THE PLANAR HELIX ANGLE AS A FUNCTION OF RADIUS C GIVEN AN K ; I TI AL CONTOUR C SUBROUTINE HELIX ( R S , R L , R A D , A L , N P T S , T O ) COMMON R ( 1 0 0 ) , Z ( 1 0 0 ) . A N G ( 1 0 0 ) , T ( 1 0 0 ) T A N T = ( R S + R L ) / ( A L + 1 . 2 » R A D ) E2=ATAN(TANT) ZK=RS/TAr:T-0.6»RAD 1 = 1 2 1 = 1 + 1 ZP=Z(1)*PAD+ZK CPHY=ZP*TANT/(R(I )»RAD) PHY=ACns(CPllY) MEND=1 1-^ S L O P = p n L Y ( Z ( I - l ) , Z ( I ) , Z ( I + 1 ) , R ( l - 1 ) , R ( l ) , R ( l + l ) , I . i E N D ) SP!IY=SIM(PHY) P1-ATAN(CPHY*SL0P) T.fi1 = S P ! l Y » C 0 S ( R 1 ) » S L 0 P T A 2 = S P i i Y * C 0 S ( B 2 ) « S L 0 P P3=TA1-TA2 P 3 = P 3 * P 3 A1 = .A.TAN(TA1) A2=-ATAN(TA2) P 1 = 1 . 0 / C O S ( A 2 ) P 1 = P 1 * P 1 P 2 = 1 . 0 / C O S ( A 1 ) P 2 = P 2 * P 2 P 4 = S l N ( ( P 2 + B 1 ) / 2 . 0 ) P 4 = 4 . 0 * P 4 * P 4 R H S = ( P 1 + P ? - P 3 - P 4 ) » C O b ( A 2 ) * C O s ( A 1 ) / 2 . 0 IF ( M E N n . E O . 2 ) 1=1+1 ANG(l)=ACOS(RHS) T ( I ) = T 0 » C 0 S ( B 2 ) / ( R ( I ) » R H S ) IF ( l , E n . M P T £ - 1 ) GO TO 11 IF (MEND.EO.2) GO TO 1 2 GO TO 2 11 MEND=2 GO TO 13 12 ANG(1)=R2 T ( 1 ) = T 0 RETURN END

Cont

B13

-POINT EVALUATION OF CONTOUR SLOPE DP./DZ

FUNCTION P C L Y ( X 1 , X 2 , X 3 , R 1 , R 2 , R 3 , M E N D ) X12=X1«X1 X22=X2»X2 X32=X3»X3 P1=X12-X22 P2=X22-X32 T 0 P = ( R 1 - R 2 ) * P 2 - ( R 2 - R 3 ) » P 1 POT=(X1-X2)»P2-(X2-X3)*P1 B=TOP/BOT A = ( ( R ' I - R 2 ) - B * ( X 1 - X 2 ) ) / P 1 C=P1-B«X1-A*X12 POLY=2.0»A«X2+B ' IF (MEND.En.2) POLY=2.0«A»X3+R POLY=-POLY RETURN END FINISH