UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Large deflection analysis of shallow framed structures. Radomske, Brian Arthur 1972

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1972_A7 R33.pdf [ 5.4MB ]
Metadata
JSON: 831-1.0050543.json
JSON-LD: 831-1.0050543-ld.json
RDF/XML (Pretty): 831-1.0050543-rdf.xml
RDF/JSON: 831-1.0050543-rdf.json
Turtle: 831-1.0050543-turtle.txt
N-Triples: 831-1.0050543-rdf-ntriples.txt
Original Record: 831-1.0050543-source.json
Full Text
831-1.0050543-fulltext.txt
Citation
831-1.0050543.ris

Full Text

cl LARGE DEFLECTION ANALYSIS OF SHALLOW FRAMED STRUCTURES by 4* BRIAN A. RADOMSKE B.A.Sc. Univers ity of B r i t i s h Columbia, 1968 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in the Department of C i v i l Engineering We accept this thesis as conforming to the required standard: THE UNIVERSITY OF BRITISH COLUMBIA August 1972 In presenting th i s thes i s in pa r t i a l f u l f i lment o f the requirements for an advanced degree at the Un ivers i ty of B r i t i s h Columbia, I agree that the L ibrary sha l l make i t f r e e l y ava i l ab le for reference and study. I fur ther agree that permission for extensive copying of th i s thes i s for scho lar ly purposes may be granted by the Head of my Department or by his representat ives. It is understood that copying or pub l i ca t i on of th i s thes i s fo r f i nanc i a l gain sha l l not be allowed without my wr i t ten permiss ion. Department of C i v i l Engineering The Univers i ty of B r i t i s h Columbia Vancouver 8, Canada Date September .11, 1972. i ABSTRACT E l a s t i c structures exh ib i t i n s t a b i l i t i e s which ar i se through the occurrence of f i n i t e displacements even when const i tut i ve properties re-main l i nea r . A non-l inear analysis which recognizes rotat ions i n the s t r a i n displacement re lat ionsh ip is formulated fo r analyzing three-dimensional framed structures. A f i n i t e element method i s used whereby the rotations with in each element are r e s t r i c t e d in s i ze by use of a loca l element reference frame attached to the element. Two such coordinate systems are developed. Then an incremental so lut ion technique based on an instantaneous l i n e a r i z a t i o n of a Taylor ser ies expansion of the forces about the displacement con f i g -uration at the beginning of each increment i s developed. The snap-through buckling of shallow frames, arches and domes i s studied with a view to documenting the e f f ec t on the equi l ibr ium paths of the type of moving coordinate frame, the number of elements, and the s ize of the increment step. i i TABLE OF CONTENTS Chapter Page I INTRODUCTION 1 II COORDINATE SYSTEMS 4 1. Configuration Space 4 2. The Global Coordinate Systems 9 3. Element Coordinate Systems 12 III DISPLACEMENT RELATIONSHIPS 17 1. Global Degrees of Freedom for an Element 17 2. Transformation Matrices i f - * and i f - K 19 3. The I n i t i a l Coordinate Frame ( y 0 ) 24 4. The Incremental Structure Deformations 5r 30 IV TANGENTIAL REFERENCE FRAME 33 1. De f in i t i on of the Element Degrees of Freedom 33 2. Derivation of the Trans!ational s 33 3. Derivation of the Rotational Components of s 36 V SECANT REFERENCE FRAME 39 1. The Element Degrees of Freedom § 39 2. Derivation of s as a Function of £ 41 VI FORCE-DISPLACEMENT RELATIONSHIPS 45 1. Element S t i f fness 45 2. Global System Equi l ibr ium 48 3. Derivation of a* and v * Matrices 51 4. Incremental Solution Technique 52 i i i Chapter Page VII BIFURCATION BUCKLING 54 VIII NUMERICAL EXAMPLES 57 1. Wi l l iams ' Toggle 57 2 . A rgyr i s ' Arch 64 3. Wright 's Ret iculated Shell Segment 68 4. Three Dimensional Elbow 74 5. Ring Dome " 77 IXN CONCLUSION 81 BIBLIOGRAPHY 81 APPENDIX 84 LIST OF FIGURES Figure Page 2.1 Global Coordinate Systems 10 2.2 Tangential Element Coordinate System 14 2.3 Secant Element Coordinate System 16 3.1 Global Degrees of Freedom at node j 18 3.2 H - X Transformation 22 3.3 4*^ " Transformation 22 3.4 ^ ^ T r a n s f o r m a t i o n 22 3.5 I n i t i a l Coordinate Frame 26 3.6 S ingu la r i t y of I n i t i a l Coordinate Frame 29 3.7 Incremental Rotations at node j 32 4.1 Deformations S for Tangential Frame 34 5.1 Deformations s fo r Secant Frame 40 5.2 Rotational Degrees of Freedom s 2 and s 3 44 8.1 Wi l l iams ' Toggle with Central Load 58 8.2 Incremental Deflection Solutions for Wi l l iams ' Toggle 63 V Figure Page 8.3 A rgy r i s ' Arch 65 8.4 Load Deflection Curve fo r A rgy r i s ' Arch 67 8.5 Wright 's Shell Segment 70 8.6 Load Deflection Curve of Wright 's Segment with Points 6 f i xed against t rans la t ion 72 8.7 Load Def lect ion Curve of Wright 's Segment with Points B free to t rans late 73 8.8 Three-Dimensional Elbow 75 8.9 Ring Dome 79 8.10 Load Def lect ion Curve for Ring Dome 80 8.11 Deflected conf igurat ion of a Dome Rib 80 vi LIST OF TABLES Table Page 1 Solutions for Wi l l i ams ' Toggle 60 2 Wi l l iams ' Toggle - 0.01 inch displacement increments 61 3 Wi l l iams ' Toggle - 10 elements per ha l f 61 4 A rgy r i s 1 Arch 66 5 Wright 's Shel l Segment 71 6 Three Dimensional Elbow 76 ACKNOWLEDGEMENTS I am grateful fo r a l l the help and guidance given to me by my supervisor, Dr. N.D. Nathan. I also would l i k e to thank the National Research Council of Canada for f i nanc ia l assistance. CHAPTER I INTRODUCTION C las s i ca l s t ructura l analys is implies a unique so lut ion to every s t ructura l problem since the theory i s based on i n f i n i t e s ima l l i n ea r d i s -placements. In f a c t , real structures exh ib i t i n s t a b i l i t i e s associated with non-unique solut ions and, in order to detect these, i t i s necessary to introduce non-l inear analys i s . I n s t a b i l i t i e s may ar i se through non-l inear material properties or through the occurence of f i n i t e displacements even when the cons t i tu t i ve properties remain e l a s t i c . We confine our attent ion to the l a t t e r , the e l a s t i c i n s t a b i l i t i e s associated with f i n i t e displacements. E l a s t i c i n s t a b i l i t i e s can be of three kinds: 1. B i f u r ca t i on , 2. Snap-through, 3. F i n i t e Disturbance. The f i r s t may be detected by the so lut ion of the c l a s s i ca l eigenvalue problem which has been formulated to include rotations of i n f i n i t e s ima l elements in the equi l ibr ium equations of e l a s t i c i t y . Recognition of the second and t h i r d requires a so lut ion technique capable of t rac ing the equi l ibr ium path a f te r the advent of f i n i t e displacements. This necess i-tates the recognit ion of rotations in the s t ra i n displacement r e l a t i o n -ships. Frame structures exh ib i t e i ther b i fu rcat ion or snap-through buckl ing, and we l i m i t ourselves here to the study of such st ructures. In the f i n i t e element analysis of frame structures, large rotations may be dealt with in one of two methods: 1. The member properties may be deduced, with 2 respect to a reference frame f ixed in space, on the basis of the f u l l non- l inear s t r a i n -displacement re la t ions . 2. The member properties may be deduced, with respect to a moving reference frame attached to the members in question. Then, by sub-d iv id ing the structure into a s u f f i c i e n t num-ber of members or "elements", we may r e s t r i c t the rotat ions with in each element r e l a t i ve to i t s own reference frame to any desired extent. The l a t t e r approach i s used here. Member properties are based on the assumption of small rotations and s t ra ins which are small compared to the rotat ions . The assumptions are adequate for the detection of b i fu rcat ion points i n the equi l ibr ium path derived from c l a s s i c a l e l a s t i c i t y , but when f i n i t e displacements are studied, the moving coordinate systems as des-cr ibed above are required. Non-linear problems are then solved by an i n -cremental procedure based on an instantaneous l i nea r i z a t i on with in each increment of a Taylor ser ies expansion of the force vector about the d i s -placements at the beginning of the increment. Early work on the theoret ica l analysis of e l a s t i c post-buckling was performed by Koiter (2). His work, using a continuum mechanics approach, centered on the invest igat ion of imperfection sens i t i ve s t ructures. Br i tvec and Chi 1 vers (5) developed matrix methods based on a potent ia l energy formulation to analyze the i n i t i a l post-buckling curves of r i g i d l y -jo inted plane frames. Martin (6), Supple (7), and Roorda (8), (9), have also done research into i n i t i a l slopes of post-buckling paths with and 3 without imperfections. Snap-through buckling of shallow arches and plane frames have been studied by Argyris (1), Jennings (4) and Williams (10). Ebner and Ucci ferro (12) compared several f i n i t e element methods and t h e i r a p p l i -cations to geometrically nonlinear s t ructura l problems. Their work was confined to planar structures only. A complete summary of the f i n i t e element analys is of nonlinear structures i s given by Ma l l e t t and Marcal (3). The present work takes element s t i f f ne s s matrices from Nathan (11). Two-dimensional structures are studied with a view to documenting computational experience. Snap-through buckling of frames and arches i s invest igated with respect to factors such as element s i ze and number, increment s i z e , and choice of element moving reference system. The necessary re lat ionsh ips are developed for extending the work of previous invest igators into three-dimensions. A simple three-dimensional space dome element i s studied and results compared to theoret ica l work by Wright (13). The snap-through buckling of a large r ing dome i s then studied. No attempt fs made to study b i fu rcat ion buckling of the structures presented. Such modes of i n s t a b i l i t i e s are prevented from occurr ing. 4 CHAPTER II COORDINATE SYSTEMS 1. Configuration Space We r e s t r i c t the class of structures being studied to plane and space frames, i . e . where the length of each member predominates over i t s width and thickness. The structure can then be subdivided into any number of l i n e elements connected to each other at "nodes". Displacements of a l l points with in each element are completely determined by the displacements of the nodes. Thus, i f there are a t o ta l of n nodal displacements or "degrees of freedom" defined by the vector r , the conf igurat ion of the structure w i l l be completely determined by the pos i t ion of the point with coordinates r in an h-space, the so -ca l led "conf igurat ion space". We are at l i b e r t y to choose the locations of our nodes and thus of the number and s izes of our elements. Such choice as th is represents the essent ia l step in the mathematical i d ea l i z a t i on of the s t ructure. Obviously there must be s u f f i c i e n t constraints on the degrees of freedom to provide for global equi l ibr ium of the structure as a whole. We begin by def ining the necessary coordinate systems i n which we measure displacements and forces and then we re late these reference frames to each other. We require a global cartesian reference frame (X) W 1 ' * - n base vectors X. whose o r i g i n and axes or ientat ions are quite a r b i t r a r i l y f i xed in space. We measure a l l nodal coordinates and global forces in th i s f i xed frame. frame we envision an arb i t rary member TK , 5 connecting nodes j to k ; at each of these two nodes are f i xed loca l global reference frames. These two frames, termed the ( X J ) and (x K ) frames with base vectors X and X respect ive ly , have or ig ins f i x e d , with reference to the frame, by pos i t ion vectors j ? J and P k respec-t i v e l y . The pos i t ion vectors describe the distance from the (^X. ) ^ r a m e or i g i n to the i n i t i a l pos i t ion of the nodes j and k respect ive ly. The U J ) and ( X k ) frame axes are f ixed pa r a l l e l to the global reference frame Q C ) axes. Within these loca l global frames we measure the nodal deformations of the structure T . The rotations r1"" are taken about the o r i g i na l ( X ) a x e s i n a defined order and the trans lat ions are in the global axes d i rect ions where: <~3 & and r 1 to 6 At node j of member J K we define an i n i t i a l loca l reference frame . A 0 named with base vectors *j . The coordinate axes of th i s frame are defined by the or ientat ion of the undeformed member i t s e l f ; the base A i / .Jv A „ vectors X of the IX ) frame are related to the y base vectors of ) by an orthonormal transformation ^- : ^ ° a A j where *4- i s a matrix function of the constant angular rotat ions ( r^.° r^ s ° *~fc° ) T def in ing the i n i t i a l o r ientat ion of the member. The der ivat ion of the transformation and the determination of the f ° o \ T i n i t i a l angular rotat ions \ *H *\, J w i l l be covered in Chapter I I I . Now we must re late the pos i t ion of the deformed member 3" < to the loca l global reference frames (xM and ( X K ) by two moving coordinate systems (y ) and (y ) with base vectors \j and y * a f f i xed to I A \ "~ nodes J and K. respect ive ly. The base vectors X are re lated to the A i j base vectors y by an orthogonal transformation V~ whose der ivat ion follows in Chapter I I I : y = t -where ^ i s a matrix function of the large global rotat ions fc* rs* r f o * ) T w h i c h correspond to ( r* s° r i 0 ) T and which define the current or ientat ion of the member tangent at node j . S im i l a r l y , at node K we have 5 = f x K / ^* ^* " * \T where A f is a matrix function of large rotations ( r * l o j . In order to evaluate the F we require to add the e f fec t of deform-ations to the angles def ining the i n i t i a l o r ientat ion of the member. How-ever, £ ° and r*~ are not d i r e c t l y addi t ive and therefore we introduce intermediate measures of deformations r . These are pecul iar to a given - V o member and are given by r ~jr - r , where r 2 lo If. The £ are re lated to the structure nodal deformations £ by an i nc re -mental re lat ionsh ip 8 where Sr. The global t rans lat ions r f o r the element are equivalent to the structure nodal t rans lat ions of r but the element rotat ions r*" can be re lated to the structure nodal rotat ions of r only by the incremental equation given above, where - t r = G 9 and r -3 1 2 Each member has associated with i t a loca l moving coordinate system in which element displacements are measured and in which the s t i f f ne s s matrix i s formulated. The loca l moving frame base vectors are functions of the global deformations r The element s t i f f ne s s matrix can be transformed to t h e Q ^ system and thus, u l t imate ly , the response of the member in global coordinates can be re lated to the response of the structure in loca l element coordinates. 2. The Global Coordinate Systems A right-handed cartesian reference frame ca l l ed the structure global coordinate system (X^ i s shown in Figure 2.1. This a rb i t ra ry reference frame remains f i xed in space and a l l geometry u lt imately relates to th i s coordinate system. Let P J and p be vectors defining the pos i t ion of nodes J and k of member T < respect ively where: 10 A vc FIGURE 2.1 GLOBAL COORDINATE SYSTEMS 11 and (2.1) i P K i f F* and p * are i n i t i a l pos i t ion vectors, then the s t ra i gh t l i ne element J k L has vector d i rec t ion P^- p* . Node coordinate systems (x*) and ( X * ) defined at nodes j and k respect ively have base vectors A j A X J and x where: J A. X - I X A X J - 2 - 3 and X x K - 3 (2.2) A j A < The components of the X and X are themselves vectors, they define the vector d i rec t ion of each of the coordinate axes. We assert that the base vectors J<_ , x , X are the same since the X and X base vectors remained f i xed in o r ientat ion throughout a l l and any load d i sp lace-ment history of the s t ructure. The i n i t i a l pos i t ion vectors and f i x the or ig ins of the ( x ^ and C x k ) systems. We measure the nodal degrees of freedom f* or structure (global) generalized displacement co-ordinates i n these node coordinate systems. 3. Element Coordinate Systems Structura l problems invo lv ing f i n i t e displacements can be solved, as previously s tated, by one of the fo l lowing procedures. F i r s t l y , we can include the e f fect s of f i n i t e rotations wi th in the element boundaries when der iv ing the element s t i f f n e s s . Secondly, we can subdivide the structure into many elements so as to reduce the rotat ions and trans lat ions of the element to with in small acceptable bounds. A less ref ined element s t i f f -ness i s used and the problem of the large rotations i s handled by the use of moving member coordinate systems f ixed to the elements in question. The l a t t e r procedure is employed in th i s work. We use an incremental load and displacement technique to fol low the load-displacement behaviour of the st ructure. At each increment step we recalculate the s t i f f ne s s of each element and reassemble the global s t i f f n e s s , an instantaneous s t i f f -ness tangent to the real load-displacement surface. We define a moving member coordinate system which i s f i xed to each element, and which moves through the global displacements ^ . In f a c t , two such element coordinate systems - the f i r s t ca l l ed the tangential reference frame, and the second ca l l ed the secant reference frame are studied herein. The tangential reference frame . 6j*") has unit base vectors given by ^ where: A t S - 2 _ 3 (2.3) i / -fc\ ' " A . The ^cj ) frame has i t s o r i g i n at node j of the element J K . The y base vector i s tangent to the centroidal axis of the element at nodej . A t The \j base vector i s coincident with the major p r i nc ipa l axis of the - 2 element cross sect ion at node j while the base vector i s coincident with the minor p r inc ipa l axes. S t i f fness matrices are presented fo r rec-tangular cross sections only. Figure 2.2 i l l u s t r a t e s th i s reference frame. A l l element deformations in th i s system w i l l be defined by the dis-placements r e l a t i v e to these reference axes. The secant reference frame ^ t j S j h a s unit base vectors given by tj where: ^ s A s y - 3 (2.4) A 3 In th is frame the M base vector i s defined by the vector jo in ing node ^ to node K. . The pr inc ipa l axes of the cross section at node j are not in general normal to the ^ ax is . A rb i t ra ry def in i t ions are therefore FIGURE 2.2 TANGENTIAL ELEMENT COORDINATE SYSTEM given below for the d i rect ions of the coordinate axes L $ 2 and A f te r given displacements at nodes j and K we can ca lcu late the base vector as: where q = f c i , c i 2 5 3 ^ T a r e t n e projections on the global axes of the i n i t i a l member J K . Figure 2.3 shows the secant coordinate frame. We define the other base vectors by Equations 2.6 and 2.7. A - 2 / / A^. A S \ 7 7~fc A S \ t j 5 = - 3 " (2.6) A S A s A 5 (2.7) ^ 3 = VJ, x y In th i s coordinate system the element deformations are described by displacements measured at both nodes of the element. FIGURE 2.3 SECANT ELEMENT COORDINATE SYSTEM CHAPTER III DISPLACEMENT RELATIONSHIPS 1. Global Degrees of Freedom for an Element Each element has twelve global degrees of freedom, s i x at each node, which are defined by the corresponding structure degrees of freedom £" at that node. At each node there are three trans!at ional degrees of freedom act ing along the (X ) coordinate system axes and three rotat iona l degrees of freedom whose axes of rotat ion depend upon the or ientat ion of the mem-ber in space. Since the intent ion i s to deal with f i n i t e rotat ions , the de f i n i t i on s of these axes of rotat ion i s rather complex. Figure 3.1 shows the s i x degrees of freedom for node J of an element J K where s ing le arrows denote trans lat ions and double arrows denote rotat ions. — t We c a l l the s i x element trans!at ional degrees of freedom r where - -t P" i s defined i n Chapter II as r f c * r 2 r6 These trans lat ions can be re lated to the element or loca l displacements of the element by a vector transformation. However, the f i n i t e global ro ta -tions where r i s defined in Chapter II as i * ^zy — r r = FIGURE 3.1 GLOBAL DEGREES OF FREEDOM AT NODE j do not transform as vectors. Problems of th i s type can be handled by spec i fy ing a pa r t i cu l a r order fo r the rotations or by the use of Euler angles. However, since there i s no s ing le system of Euler angles covering a whole sphere without ambiguity, the former method w i l l be used here. I t i s to be noted that three f i n i t e rotations taken in an order d i f f e ren t from that spec i f i ed w i l l lead to a d i f f e ren t or ientat ion i n space. Therefore, the rotations at each node of our element w i l l be defined as having occurred i n a s p e c i f i c order and this order w i l l be adhered to throughout the development of th is thes i s . I t i s asserted that th i s w i l l lead to no insurmountable d i f f i c u l t i e s provided care i s used in problem input and interpretat ion of re su l t s . Ambiguities can thus be avoided. 2. Transformation Matrices ¥ and ± The H-*^  matrix was introduced i n Chapter II as a matrix function of ~ ^ j rotations r e l a t i ng the base vectors of two coordinate systems, the X. of the (xJ) frame and an a rb i t ra ry frame with base vectors . The U J w i n be the transformation of the X vectors in to the vectors at node j of element T K . 5J = V X J (3.1, We have defined ^ previously as a matrix function of f i n i t e rotations ( n / rs* C ) T where: o 0 a ' The ( ITs define the i n i t i a l member or ientat ion i n space. We propose to order the rotations as fo l lows: f i r s t X* » then t% , and, f i n a l l y , Yl . *- A J F i r s t of a l l , the t*. rotat ion transforms from the base vectors X to a set of instantaneous base vectors X of an intermediate coordinate frame by Equation 3.2: A I , 1 A i X = V- X (3.2) I I V- i s a s ingle rotat ion as shown by Figure 3.2 and ^ i s given by: / o O 0 CosCrf) S n (rH*) o - S » h ( V ) CosCrf) and A r X : - 2 X T - 3 Secondly, the rotat ion transforms from a known set of base vectors A j ATJ X ' to the instantaneous coordinate frame base vectors X mediate frame ( X ^ ) by Equation 3.3: of an i n t e r -(3.3) where C o 5 C r * ) o - S , n ( V ) O / o O C o S ( r 5 * ) 21 and A H X -' A H N - 2 A I X - 3 A I The rotat ion was about an instantaneous axis defined by the X * base vector. F i n a l l y , the r o ta t i on , which i s about the instantaneous Al A j j . X 3 base vector, transforms to the base vectors y of the [y ) frame A T T / n \ "~ from the known base vectors X of the I X ) frame by Equation 3.4: A j m. A , TT_ y = V- x (3.4) where CesCrJ*) S i n ( V ) O - Sin ( V ) O o o 1 The rotat ions r s and are shown in Figures 3.3 and 3.4. tions 3.2, 3.3 and 3.4 we derive Equation 3.5: 9 » = v - 1 V A ; j A : y J = v / J x where From Equa-(3.5) v f J = ^ ^ ^ and Cos{tSi CosCrtf) -r5ir,(r4*)s;n(r/)Ge/< +Sn (V) Sin CC) r J = - Q>s(rSf SifxCci) Cos(<)s.'n(r/)slrtc:r*) + S.nCrv*) a 5tV) Sin (Y s*) - Sn<V) c7>sO/) CosCrA*) G>s(r£) where the order of rotations i s , again, , rg , r * I t can be ea s i l y shown that i s an orthonormal transformation since i t has the fo l lowing properties and [<tf=[+i]T where superscr ipt T denotes the transpose of the matrix and [ i j i s the i den t i t y matrix. Thus An a rb i t ra ry vector v* in the frame with components f -vr =• W l 11 have components in the ) frame given by UJ i n Equation 3.7: to = H" v- (3.7) where to = to, CO: are the components of u) , which are the components of AT measured i n the ( ^ ) frame. At node k of element J K we have a transformation matrix H~ which was introduced in Chapter II as a matrix of rotat ions re la t ing the base vectors of the frame and the base vectors of an a rb i t ra ry ( y ^ ) frame. The V- w i l l be the transformation of the X vectors into the * j vectors by Equation 3.8: A k K A K *d = f- 5 (3.8) The 4~ transformation i s composed of large angle rotations ( ^ o * ^ V ^ T 1 1 1 t n a t s p e c i f i c order and i s s im i l a r in form to CosCr*) C o s f r , 2 * ) -CosCn*)srnCr,|) -r S»o<-n<>*) Q>sCnz) Sin (r*) The I n i t i a l Coordinate Frame In Equations 3.6 and 3.8 we defined the relat ionships between the base vectors of the node global coordinate systems and the a rb i t ra ry coordinate 25 systems ( * j J ) and ( ^ k ) at the j and k nodes of member JK respect-i ve l y . I f we give the base vectors x. and X the fol lowing values: (3.9) k , . . ^ \ ' X where € = Co o i )T — 6 A l j then the base vectors ^ w i l l be the rows of ^ and the base vectors A K tc -^ w i l l - b e the rows of M- , i . e . : l (3.10) C\K * ( v ^ G , . ) i f ' ( 1 , 2 ) 4 - H c J B ) ) T where and 2 I n i t i a l l y , before any deformations have been applied to the s t ructure, the elements w i l l be s t ra ight and the (<j J ) and ( y ^ ) frames w i l l be al igned in space. The transformation matrices and 4.^  reduce to an i n i t i a l transformation matrix tf° which relates the i n i t i a l nodal coordin-ate frame ^ 1 base vectors and the X or X base vectors. We can ca lcu late the i n i t i a l angles, defined as ( r^" P s ° r j ) ' , knowing the i n i t i a l geometry of each member. This requires that we know the coordinates of a l l the nodes as well as a point which defines the i n i t i a l o r ientat ion of the member cross section in space. Figure 3.5 i l l u s t r a t e s the i n i t i a l coordinate system base vectors in re la t ion to the node global coordinate 26 FIGURE 3.5 INITIAL COORDINATE FRAME 27 system base vectors. The ( t j 0 ) frame has base vectors ^ where: A o , A © A c A N T v = ( ^> X ^ ) (3-11) The *j axis i s defined as the minor p r inc ipa l axis for f lexure of the undeformed element cross sect ion. The axis i s defined as the major p r inc ipa l axis and c^ i s the centroidal axis of the element. Now, i f we wr i te : X J = I f J ^ (3.12) with and s „ y 3 3 ) (3.13) then the components of *-J° and cj , t j * , , e t c . , are eas i l y deter-_ X . A O A o mined from the coordinates of the nodes. ( y i s dependent upon y and A o \j and y ie lds no new information) But, from Equation 3.10, we know the ~~ A O A & o components of y and *j in terms of , rs° , and , v i z . Equa-t ion 3.14: A 0 / E S l n r C ) S n ( r 5 0 ) C o s ( r f c ° ) -Cos r r + * ) s ' . o f r s 8 ) ro S Cr 6 ' ) \ T ^ , = ^ C o s C O f c s G - j ) + Cos C ^ ) S i n < C ) + S , n ( V ) S nCr f c") / (3.14) ( S , nC r s " ) - 5 i n C r / ) C o s t s ' ) C o s C r A ° ) G > s r r s ^ ) T Thus we wr ite the fol lowing equations which w i l l be solved for the three angles, r;° , r5- , and r° : 28 <jj = Cos Crs°) Cos (rk°) 1\2 - Co^CrH°) Sin Cr<S) + S m Cr * J ) 5 i " n C V s * ) tbs C ^ ) ai3° = Sin (rH°) Sir, Cr*>') - Cos Cn,') SV n Crs") G>-5>Cr*r) ( 3 . 1 5 ) j* 3 i - S. 'n C r 5 ° ) 3 3 / = - 5 i r ) ( ^ ' ) C o s (' r 6 ' ) C i s C r * * ) c b s C r s " ) Solving Equation 3 .15 y ie ld s 3 . 1 6 : Sin(rsD)= y 3 »° Si'n Cr*°) " ~r~r; provided Co 5 ( r e . " ) 7* O Cos C^s / (Tos Cn,0) = provided C * o s C V s 0 ) £ o o ( 3 . 1 6 ) Ca) S i o ( r b > S i l C o s ^ r 5 c ) 4 - ^ 2 3 » C o s C V ) provided " i 3 3Vo and Cb s ( r s * ) C o s C v t ) = ( b ) o r .\ - ^ 3 ° ro&Crg') + ^ 3 - 3 ^ ^ S i n C rfe, ) — - 1 - , v ~ Cos <- rS / provided ^ ^ | ^ o and C o s G ^ J / o I f ^ 3 3 = o choose Equation (b), but i f u ) 3 2 also = 0} then use the fol lowing procedure. There are two s ingular points fo r which the above equations are not va l i d and that i s when 4 ^ =1" ' ; i t fol lows then that ^*>z° ^ 0 and 4 3 3 ° - o and the C o s O . The Oj") frame would have i t s «-}3° A. A " axis along the * i d i rec t ion as shown in Figure 3 .6. The ^ { base vector 29 FIGURE 3.6 ' SINGULARITY OF INITIAL COORDINATE FRAME w i l l l i e in the X ? ^ 3 plane and thus the component w i l l be zero. Therefore: Cos ( r</) = -S i n ( < V ) = 0 . ^ C O - O ( 3 . i 7 ) S in ( r s ' ) - - X I s»'n ^ n») = y 12. 4. The Incremental Structure Deformations S r The re la t ionsh ip between the intermediate incremental deformations S f and the structure incremental nodal degrees of freedom Sr follows using equations 3.2 and 3.3. I t w i l l be reca l led that £ ~ r r are the i n -crements in the rotations r* a r i s i ng from increments <Sr r i n the structure degrees of freedom P . Thus , S~H , and <§r H are a l l rotations about the x| ax i s ; but and § r 5 are rotations about I j v- c _ the X 2 axis while Sfs i s about the Xz ax i s ; and and a>r f c TT" J are rotations about the axis while S r t i s about the axis. At node k s im i l a r observations can be made. Therefore the transformation from Sr to S i - i s given by where: T o o o o o o o o I p o o o — 2 » o o o 1 O O O Cos C *"icT) S o Cn 6*) I o o I - o I o o o I O o O o = o o o o o o Figure 3.7 shows the incremental rotations for node j of element T K . 32 FIGURE 3.7 INCREMENTAL ROTATIONS AT NODE j CHAPTER IV TANGENTIAL REFERENCE FRAME 1. Def in i t ion of the Element Degrees of Freedom S The displacements of the nodes of an element with reference to the loca l coordinate system, the ' l o c a l degrees of freedom' 5 are functions of the structure degrees of freedom r defined previously i n Chapter I I I . For the tangential reference frame there are s i x degrees of freedom per element, a l l defined at one node as shown in Figure 4.1. The S are kept small with respect to the length of the element by su i tab ly increasing the number of elements required to model a given structure under a large rotat ion or t rans la t i on . The trans lat ions I S, S 2 S 3) act along the tangential reference A-fc v T frame base vectors y respect ively and the rotat ions C 5 S x . ) are rotations about these same base vectors. Actual ly the ( S^. S s S t ) ^ are only approximately about the y base vectors, as w i l l be shown in the fo l lowing sect ions, but the approximations w i l l be found acceptable. The j frame coincides with the ) frame defined in Equation 3 . 5 , i . e . q*= ^  V ( 4 . D 2. Derivation of the Translat ional S Figure 4.1 shows the element in the tangential reference frame before and a f te r the f i n i t e displacements r , in which the element deformations are exaggerated for c l a r i t y . Consider the element J X with j at the o r i g i n of the ( * J ) frame FIGURE 4.1 DEFORMATIONS S FOR TANGENTIAL FRAME 35 i n i t i a l l y and having length L . Before the displacements, the nodal co-ordinates of J in ( X J ) are C o o o " ) and of K. are given by the components of A where: A - ( 3 t a2 a , ) T and ( 4 . 2 ) Af te r the displacements where r f c i s defined in Chapter I I , node j has (x^ ) frame coordinates of and node k has coordinates ments are In the (.y^j frame, the coordinates of j and k before d i sp lace-respect ive ly. I f the displacement were a r i g i d and o o body motion, these coordinates would remain unchanged; node K would not move in the ) frame. The dif ference between the actual motion of node K and that corresponding to a r i g i d body motion i s the deformation vector ( S i S 2 • Therefore, the deformations (S{ S 2 S 3 ) T , seen as components of displacement of k with respect to J in the C-*J) frame, are given by Equation 4.3 and ca l l ed S : whe re / L a , a, r 7 - L (4.3) ( »,o <4-i ( 1 , 2 ) 4 j 0 , 3 ) are the x^ -components of the length of the ele-A t ment which i s or iented along the *j d i r ec t i on . In the Cd ) frame, the loca l coordinate system, the vector 5 r ( 5 | S z S , J i s the S transformed by Equation 4.4. 5 R = S (4.4) 7-^  has been derived in Chapter III and i s composed of rotations given as in Chapter II: r •* r* 3. Derivation of the Rotational Components of £ We have the, fol lowing transformations from Equations 3.5 and 3.8: ^ * I -We can re late base vectors at nodes j and K. from Equations 3.5 and 3.8 A j A k by noting that the base vectors X and x are i d e n t i c a l . Thus: • iT and 3 = *1 i (4.5) | [ * f - K ] which re lates ba se We define the transformation vectors ' j and . This transformation matrix, composed of three rotat ions ( e , ^ 3 ) J w i l l n a v e the fol lowing form: 37 Cos C©0 S.n ( e 3 ) + 5>n(A) s r n ( o 3 ^ - Cos Co2) Sin(O-i) Cos C&t) CzsC&a) - Si"nf&,)SinCefe)Sn(5%) + S r o ^ C o s C O j ) S i n C © 2 ) Although this transformation matrix implies that &\ i s about the (4.6) ax i s , e>2 i s about some intermediate axis between y 2 and «j , and i s about the M K ax i s , the rotations are small and the axes of rotat ion 3 -t can be considered to be the ^ axes. I f we assume that the sine of the angle i s the angle i t s e l f and the cosine of the angle i s un i ty, 3f reduces to: l . o V = l.o © « &2 l . o (4.7) We make the assumption that ( s 5 s t o ) are the rotat ions ( &x e>z e>3 ) respect ive ly. I t i s f e l t that these approximations w i l l not a f fec t the required der ivat ives of s with respect to jr . From Equations 4.5 and 4.7 we obtain Equations 4.8: S 5 = (4.8) In summary, fo r the tangential coordinate system, we have: 5 ; > . r ( a . + r , - ^ q - ^ E , ! ) ^ ^ 2 - t r 8 - ^ ) c f J ( 2 ) 2 ) + ^ 3+rq-r3)v/.Y2J3) (4.9) 39 CHAPTER V SECANT REFERENCE FRAME 1. The Element Degrees of Freedom § In the secant reference frame, l i k e the tangential reference frame, the element w i l l have s i x loca l degrees of freedom S . This number res-u l t s a f t e r subtracting the s i x r i g i d body motions from the twelve global degrees of freedom fo r the element. The s i x retained degrees of freedom consist of an elongation, four bending rotat ions , and a tors ional r o ta t i on ; two bending rotat ions are at node j and the other four are at nodek . Figure 5.1 shows the member before and a f te r f i n i t e displacements r , where the member displacements have been exaggerated fo r c l a r i t y . The f igure shows the base vectors of the ( X ^ ) system transformed through o4* n*yto the W system at node j and through ( Vi©* n * ^ 2 * ) T t 0 t h e ^J^j system at node k by Equations 3.5 and 3.8 respect ive ly . The secant reference frame base vectors y have been defined by Equations 2.5, 2.6 and 2.7. Equation 2.5 eliminates t rans la t iona l r i g i d A 5 ^ 5 body motions. Equation 2.6 sets y orthogonal to y and the tangent A -fc _ A 5 base vector M ( i . e . orthogonal to the plane containing y and the minor p r inc ipa l axis of the cross section i n i t s deformed pos i t ion at j ). Equation 2.7 completes the right-handed orthogonal t r i a d . We w i l l also need a loca l coordinate system at node k al igned with A S the ^ ( base vector in order to define the element degrees of freedom adequately. We c a l l th i s coordinate frame theC^ S ) f r ame with base vectors ^ which are given in Equations 5.1, 5.2 and 5.3. FIGURE 5.1 DEFORMATIONS S FOR SECANT FRAME A s A s 1. A s 2 , A k A - 5 j 3 x g» A „ . A S ^ 3 X ^ S 2 . a A „ A s 2 , s X z | (5.1) (5.2) 41 (5.3) 2. Derivation of s as a Function of r We begin by noting that since the element deformations are required to be smal l , then the ) coordinate system and the ( . H J coo rd i na te systems w i l l be almost coincident. Likewise, the (cj K Jcoord inate frame and ( 2 s ) c o o r d i n a t e frame are almost coincident and therefore the rotat ion transformations re la t ing t he i r base vectors w i l l be composed of small angles. We can say, then, that the fol lowing s imp l i f i ca t i on s to the sca lar products hold: ^ . Cj S = C O S Cdt)= l.o s ince <X, i s considered small QJ . 0 s = cos ( ? o ° r ^ ) = £ « Z . * p = cos (9o° ±«>) = + * 3 c *^ , are also small angles 42 The elongation of the element 3 t i s determined by the extension of the l i ne jo in ing node J to node K : where i_ defines the length of the element, L = J a* -+ a.* -t a* and ( q , a z are the i n i t i a l components of the length vector. From Figure 5.2 i t can be seen that: s2 = - A 5 y • 4 (5.5) s 3 = i * ^ 2 (5.6) At node K we obtain s im i l a r expressions fo r the and Sfc ro-tat ions : ss = - A. K /V S 2 , - 3 (5.7) A K ^ 5 ^ • 5* (5.8) The f i n a l degree of freedom, the angle of twist of the element, can be A S A 5 approximated by the sca lar product of vectors * j and . A S A 5 S = H • 2 (5.9) In summary, for the secant element coordinate system, we have the fol lowing degrees of freedom : A S s 2 ' ^ (5.10) A J A S 43 A s _ 3 A. s W 2 A s (5.10 cont'd) A- S ' ^ 2 FIGURE 5.2 ROTATIONAL DEGREES OF FREEDOM S 2 AND CHAPTER VI FORCE-DISPLACEMENT RELATIONSHIPS 1. Element S t i f fness The element s t i f f ne s s matrices were taken from Nathan (11). They were developed fo r a doubly-symmetric th in pr ismatic element with eleven degrees of freedom: three t r a n s l a t i o n a l , four r o t a t i o n a l , two tors ional and two rate of change of twi s t angle. Nathan's s t i f f ne s s matrices were reduced from eleven degrees of freedom to s i x by e l iminat ing those which are restrained in the present case. The matrices were also transformed to conform to the present coordinate system de f i n i t i on s . We present here only the s t i f f ne s s matrices fo r the cant i lever element since the secant element s t i f f ne s s can be derived from the given matrices. In reference (11) the element forces £>c are f i r s t determined as non-equations are then l i nea r i zed by an expansion of the element forces i n a Taylor ser ies about some instantaneous element displacements S° . By making the displacement increments A S su i tab ly sma l l , terms past the f i r s t can be discarded from Equation 6.1. l i n ea r functions of the displacements These (6.1) or (6.2) where s ; The matrix k can be wr i t ten j? = k „ !5 « » where the s t i f f ne s s matrix k Q i s independent of S ° and the matrix k\ i s l i n e a r i n s" . Use of the matrices given below requires that elongations, shear displace-ments, and rotat ions remain smal l . Displacements must be small compared to the element dimensions and the rotat ions must be much less than un i ty . Nathan's transformed element s t i f fnesses fo l low: A S L O I E E I , L * O O 5 y m O o o G J L-O o Co &I2 O O L 2 -O O O ^-6X3 (6.3) where I 3 i s the moment of i n e r t i a about the cant i lever coordinate system .t »2 ^ axis and I 2 is about the <^ axis. J z is greater than I , T = tors ional moment of i n e r t i a A = element area L = element length E = Young's modulus of e l a s t i c i t y G = Torsional modulus 47 k,(s26)= o o O m o O ST* o O a o o lOL. O . 2 . ^ £ o tOL. o o o 2 , o S i - * O o o o o o _ Co G J ^ o o o o - 2 / £r 2 5 I T 5 o IO L O o o (6.4) (6.5) o o O o o O o o &£-i o u o o o o O o o 5 Z_z o o (6.6) o o O A S O O o -21 5 l_2 o o 2 f)£ • 5 O o o O O O o 5 |_ O O (6.7) It, (sj)* o -I D u O o O O o O zt e i i 5 o D O o 5 L . O Z A £ '5 O o O O O (6.8) 2. Global System Equilibrium We turn now to the problem of deducing the system response in global coordinates from the element equilibrium equations. We suppose that there are h independent global degrees of freedom P . The element or local degrees of freedom S have been expressed in previous chapters as functions of the coordinates ir . Si= ( ( * ) I - 1 ,2 , nri ( 6 9 ) The element incremental displacements follow as: = <3ij Cr) b j k Cr) S r K (6.10) or expressed i n matrix notat ion: (6.11) The transformation matrix <X i s a function of the global coordinates V . By the p r i n c i p l e of v i r t ua l work, the external work of the R forces moving through the corresponding v i r t ua l displacements S r must be equal to the interna l work of the element forces $ moving through the compat-i b l e v i r t ua l displacements S~s . Thus, (6.12) Using the compat ib i l i ty Equation 6.11, we have (6.13) I f 6.13 i s true for an a rb i t ra ry v i r t ua l displacement § V > w e c o n ~ elude that : 5 = « T £ (6.14) We expand the R c in a Taylor ser ies about an i n i t i a l global displace-ment £ ° , on the assumption that there are no s i n g u l a r i t i e s . The domain of r i s su i tab ly r e s t r i c t ed to ensure a one-to-one correspondence between load and displacement. 50 (6.15) I f we r e s t r i c t the increments >±r such that the second and subsequent terms of the expansion can be neglected, we obtain: or in matrix form: (6.16) A R = K A r-where K = —1 r c The j th column of K i s given by (6.17) = zf ^j? a j 5 5 Tjrj + 2>f b r 2 r j b f ^ r j (6.18) where j^ j is a column vector of zeros with a one in the j th place, and, y * -C c = b b d ri Thus the global stiffness K is given by K = + K 2 (6.19) where ~ ^ and the columns of K2 = [ g , ( ^ +" « ^ ) b : i J ? The element stiffnesses te0 and k» are found in Equations 6.3 to 6.8. 3. Derivation of the g and Y Matrices Having outlined the equilibrium equations required to solve the large displacement structural problem, we next consider the a , V and C matrices necessary for the computation of the instantaneous stiffness K of Equation 6.17. Then, given the global displacements r and f and the element forces^ at any stage, we shall be able to calculate-the global instantaneous or tangent stiffness of the structure . This s t i f f -ness defines the tangent plane to the load-displacement surface at the instantaneous displacement configuration. The loads R are assumed to be a continuous twice differentiable single-valued function of the displacements r . This assumption i s , in fact, violated in the case of bifurcation buckling, which will be discussed separately below. The d i f f e r en t i a t i o n of Equations 4.9 and 5.10 for the secant or tang-ent element s t i f f ne s s systems is a mechanical process; however, the resu lts are extremely lengthy and have been included d i r e c t l y i n the com-puter program and there seems l i t t l e point i n reproducing these results here. 4. Incremental Solut ion Technique When the loads are s ingle-valued functions of displacements, an i n -cremental displacement method can be used. This method i s most useful i n fo l lowing the complete load displacement h istory of snap-through buckling problems where the same tota l load occurs at d i f f e ren t displacements. Otherwise an incremental load method may be used. No attempt i s made here to analyze structures whose load-displacement curves involve b i furcat ions of the buckling paths beyond the b i fu rcat ion point. The fo l lowing i l l u s t r a t e s the procedures involved during one increment of the so lut ion technique. At the beginning of an increment, the tota l global displacements £" at the end of the previous increment are used to ca lcu late the a* ' , fc> and V matrices defined by Equations 6.10, 3.18 and 6.18. The element d i s -placements § at the end of the previous increment are used in the s t i f f -ness matrix and then the matrix products of Equation 6.19 are c a l -culated to produce the instantaneous tangent s t i f fnes s matrix K of Equa-t ion 6.17. The s t i f fnes s matrix < i s to be inverted by a band inversion rout ine. For maximum so lut ion e f f i c i ency the band width of the matrix should be a 53 minimum. Both po s i t i v e -de f i n i t e and pos i t i ve semi-def in i te matrices are allowed. A one parameter (X) load system i s applied to the structure and i s assumed to vary l i n e a r l y with displacement during the increment. For the incremental load so lut ion method, the load increment vector which i s supplied by the analyst i s used in equation 6.17 to produce the incremental def lect ion vector. However, for the displacement increment method, a un i t load vector /3 i s applied to produce the def lect ion vector $ (6.20) The load increment X i s ca lculated by l i n e a r l y proportioning the def lect ions J ° and f . A = (6.21) where j> i s the displacement increment spec i f ied by the analyst and ^ i s the displacement at the same degree of freedom as under the un i t load vector /3 Thus we have: * 5 = X @ (6.22) A N D A . i r - X ^ (6.23) These increments of global force and displacement are added to the to ta l force and displacement vectors K and r respect ive ly . Ss i s then calculated by Equation 6.11 and S/ST by Equation 6.2 before proceeding to the next increment of the so lu t ion . CHAPTER VII BIFURCATION BUCKLING Many structures exh ib i t a branching of t he i r equi l ibr ium paths that i s associated with some mode of s t ructura l i n s t a b i l i t y . For most such s t ructures , the determination of th i s ' c r i t i c a l value ' involves using the l i nea r s t i f f ne s s based on the i n i t i a l geometry plus a non l i n ea r c o n t r i -bution based on some small perturbation from the i n i t i a l pos i t ion . Thus the s t i f f ne s s matrices Ko , K| , and should be useful fo r deter-mination of the b i fu rcat ion point. The nature of the equi l ibr ium paths i n the immediate post buckling range i s of in teres t to the analyst, but an invest igat ion of th i s problem w i l l not be attempted here. Unfortunately, the so lut ion of the eigenvalue problem associated with the determination of the b i fu rcat ion point involves the invers ion of a f u l l matrix and this puts a constra int on the s ize of the problem which can be handled. But i f only a few of the smaller eigenvalues are re-quired and the s t i f f ne s s matrices are wel l banded, then an i t e r a t i v e a l -gorithm i s poss ible. The method i s e f f i c i e n t for f inding a few of the lowest eigenvalues of a pos i t i ve de f in i te matrix Ko . No attempt should be made to determine a l l of the eigenvalues by this i t e r a t i v e scheme. A problem in convergence does ar ise when two eigenvalues are the same or almost equal but we assume that this case does not occur often and that, when i t does, i t can be solved by an alternate procedure. The procedure for determining the lowest c r i t i c a l value only w i l l be out l ined here, since the in teres t of the analyst i s usually confined to this one. An i t e r a t i v e algorithm has the advantage that the so lut ion can be obtained to any desired accuracy consistent with the inherent roundoff 55 errors of the computer. We r e s t r i c t the load configuration on the structure to a s ing le para-meter system where the loads increase uniformly according to a factor A . The def lect ion vector S_r corresponding to a normalized force vector i s ca lculated using the l i nea r s t i f f ne s s Ko . ^\ and ] f z a r e then b u i l t using th i s def lect ion vector to y i e l d : k P (5, f S r ' ) 4- (Sr1) ( 7 J ) We form the tota l s t i f f ne s s of Equation 7.2 assuming that i s l i n e a r up to some c r i t i c a l value of , - \ r r • ( Ko + ^ P ) ^ - - S R ( 7 . 2 ) The c l a s s i c a l eigenvalue problem resu l t s when the right-hand side of Equation 7.2 vanishes. On rearrangement we obtain: I K 0 S r = _ K-p S P T ' ~~ " ~ ~ (7.3) where S r - corresponds to the eigenvector and, again, X<rr is the eigen-value. In the i t e r a t i v e scheme, a t r i a l vector a r i s inserted in the right-hand side of Equation 7.3, and the left-hand side is used to solve for an improved eigenvector S r " 1 . The length of the new vector i s a f i r s t approximation of the inverse of the eigenvalue \r- • We normalize S r ' " and repeat th is procedure un t i l the c r i t i c a l value X c r s t a b i l -izes according to some convergence c r i t e r i o n or unt i l the number of i t e r a t i v e cycles exceeds some given number. The rate of convergence to \ r r and S r , the eigenvector, and the time required depends on the s i ze of the problem as wel l as on the r e l a t i ve difference between the f i r s t and second eigenvalues. Derivation of th is procedure i s given in Wilkinson (16). CHAPTER VIII NUMERICAL EXAMPLES 1. Wi l l iams ' Toggle In 1964 Williams (10) extensively studied the snap-through buckling phenomena of a shallow planar frame i l l u s t r a t e d in Figure 8.1. The geo-metric and e l a s t i c properties of th is frame, ca l l ed a ' t ogg le ' by Wi l l iams, are the fo l lowing: £ 1 = 9.27 x 10 3 l b . i n c h 2 A G = 1.885 x 10 6 l b . t- = 26 inch h = 0.32 inch The large def lect ion behaviour of th i s structure under a v e r t i c a l applied load on the center l ine i s characterized by a softening region followed by a hardening region. Ebner and Ucci ferro (12) in a recent paper reviewed and compared f i ve pa r t i c u l a r methods and t he i r appl icat ions to geometric non-l inear planar problems. These f i ve formulations are: 1. Argyr is (1964), (13) 2. Martin (1965), (6) 3. Jennings (1968), (4) 4. Ma l le t and Marcal (1968), (3) 5. Powell (1969), (14). In a l l of these methods the non-l inear force displacement re lat ionsh ips II 13 . E l = q . 2 7 x t o 3 it> i n c h 2 A E = i. ess x to f c l b W s e ot 4. . 3 2 inch FIGURE 8.1 WILLIAMS' TOGGLE WITH CENTRAL LOAD 59 have been formulated in terms of s t i f f ne s s matrices. The so lut ion tech-niques incorporated were: 1. d i rec t so lu t i on , 2. load incrementation, 3. displacement incrementation. Wi l l iams ' toggle was used by Ebner and Ucc i ferro for test ing the f i v e methods. The solut ions for the toggle under a constant center l ine v e r t i -cal load of eighty pounds are given in Table 1. In view of symmetry, only hal f of the structure was modelled and only eight elements were used. The r e s u l t s , as can be expected, varied but Jennings' procedures give exce l lent agreement with Wi l l iams ' so lut ion f o r a l l three so lut ion tech-niques. The ' secant ' and ' t a ngen t i a l ' element solutions for various d i s -placement increment s izes and numbers of elements (modelling hal f of the toggle) are shown in Table 2 and Table 3. For th i s example, the ' secant ' element solut ions converged monotonically to a lower center l ine def lect ion than the ' t a ngen t i a l ' element so lut ion with the number of elements cons-tant at ten per ha l f span. However, the dif ference of 0.002 inches in 0.6159 or 0.3% i s neg l i g i b l e . When the displacement increment s ize was kept constant at 0.01 inch and the number of elements used to model the structure were va r ied , the ' secant ' elements performed extremely wel l even when using only one element. The ' t angen t i a l ' element solut ions showed a more rapid convergence although no solutions are ava i lable for using one or two elements per hal f span. This element behaves poorly when the e l e -ment displacements become too large. It i s probable that the ' t angen t i a l ' element gives a poor r e f l e c t i o n of the ax ia l deformation when element 60 Jennings Powell Mallet-Marcal TABLE 1 SOLUTIONS FOR WILLIAMS' TOGGLE (Under 80 l b . load and 8 elements per ha l f ) Formulation Direct Di rect Direct Center!ine Deflect ion (inches) 0.611 0.611 0.600 Williams Jennings Powell Martin Exact 1 l b . load increments 1 l b . load increments 1 l b . load increments 0.611 0.639 0.639 0.640 Jennings 0.007 inch, d i s p l . increments 0.616 Martin 0.007 inch, d i s p l . increments 0.621 'Tangential 1 element 'Secant ' element 0.007 inch, d i s p l . increments 0.007 inch, d i s p l . increments .6203 .6161 TABLE 2 WILLIAMS' TOGGLE O.Ol inch displacement increments and 80 l b . load 'Tangent ia l ' Elements Centerl ine Deflect ion (inches) 'Secant ' Elements Centerl ine Deflect ion (inches) 1 element per ha l f 2 elements per ha l f 5 elements per ha l f 10 elements per ha l f 20 elements per ha l f 0.6550 0.6200 0.6180 0.6312 0.6308 0.6191 0.6179 0.6179 TABLE 3 WILLIAM'S TOGGLE 10 elements per ha l f and 80 l b . load 'Tangent ia l ' Elements 'Secant ' Elements Centerl ine Deflect ion Centerl ine Deflection (inches) (inches) 0.007 inches per increment 0.010 inches per increment 0.020 inches per increment 0.040 inches per increment 0.6178 0.6200 0.6262 0.6406 0.6159 0.6179 0.6246 0.6387 62 def lect ions are large. This would be important i n the structure under discuss ion. Jennings (4) derived two element s t i f fnesses including varying degrees of geometric non - l i nea r i t y , and was able to obtain good agreement with Wi l l iams ' re su l t s . He noted that where geometric changes were s i g -n i f i c a n t , many more of his less sophist icated elements were required than of his sophist icated non-l inear elements to accurately model the load deformation behaviour of the st ructure. Although Jennings does not spec i -fy what displacement increment s izes he used, the ' t angen t i a l ' element gave s im i l a r resu l t s to Jennings' less sophist icated element so lut ions. The ' secant ' element solut ions appear to agree with the resu l t s by Jennings' better element. In general the ' secant ' elements performed better than the ' t a ngen t i a l ' elements. However, the time required for the ' secant ' element so lut ion was greater than for the corresponding ' t angen t i a l ' element so lut ion. A l so , the d i f ference in the resu lts between using twenty elements per hal f span and ten elements per half span i s smal l , ind icat ing that the u l t ra- f ineness of so lut ion i s not warranted. The reduced increment s ize seems more impor-tant than increasing the number of elements in the st ructure. Some of the load def lect ion paths for th i s structure are shown in Figure 8.2. • 3 .s .Co - 7 center/she deflection (inches) FIGURE 8.2 INCREMENTAL DEFLECTION SOLUTIONS FOR WILLIAMS' TOGGLE as 2. A rgyr i s 1 Arch Ebner and Ucciferro (12) also used a plane arch, f i r s t tested by Argyris (1), to compare the f i ve various formulations c i ted previously. The results of these formulations are given in Table 4 along with the 'tangential 1 and ' secant ' element solut ions fo r the same number of e l e -ments per ha l f span and s ize of displacement increment. The e l a s t i c and geometric properties of the arch shown in Figure 8.3 are as fo l lows: E A = 10 7 l b . E r = 10 7 l b . inches 2 r i s e = 3.14 inches radius of curvature = 400 inches length = 100 inches The arch i s pinned at the ends and the load i s applied v e r t i c a l l y at the center l ine. Since s t ra ight l i ne elements are used, at least ten elements are required to model the ent ire arch as a polygonal arc. The fewer e l e -ments used, the poorer the mathematical model of the intended s t ructure. For th i s problem, there was l i t t l e d i f ference in the results between the ' secant ' element and ' t angen t i a l ' element solut ions. Agreement with the displacement incrementation solutions of Martin, Jennings and Argyris i s good. De f i n i t e l y , smaller def lect ion increment sizes and more i nc re -ments produced ' b e t t e r ' results for the arch of a reasonable number of e le ments than for the same structure composed of many elements but a larger increment s i ze . In the snap-through buckling of this arch, the true so lu -t ion curve can be expected to f a l l below the upper snap value that the L E A = l o 7 l b ET - l O 7 l b i n c b 2 r i s e . = 3.\+ inch. f r fd ius o-f curvat ore. — i n c h . L_ = too i n c h , FIGURE 8.3 ARGYRIS' ARCH TABLE 4 ARGYRIS' ARCH , Number Formulation of elements Increment Upper Snap Lower Snap l b . l b . Martin Load Incr. 5/half 25 lb s . 2300 -Powell Load Incr. 5/half 25 lb s . 2300 -Jennings Load Incr. 5/half 25 lb s . 2300 -Martin D i sp l . Incr. 5/half 0.07 inch 2344 814 10/half 0.07 inch 2360 828 Jennings D i sp l . Incr. 5/half 0.07 inch 2360 836 10/half 0.07 inch 2358 837 Argyri s 10/half 0.157 inch 2450 800 'Secant ' Elements 5/half 0.07 inch 2328 849 . 10/half 0.07 inch 2356 838 10/half 0.157 inch 2460 846 'Tangent ia l ' Elements 5/half 0.07 inch 2325 824 10/half 0.07 inch 2355 835 10/half 0.157 inch 2460 843 FIGURE 8.4 LOAD DEFLECTION CURVE FOR ARGYRIS' ARCH ' secant ' and ' t angen t i a l ' element solut ions pred ict . Nevertheless the solut ions for the arch modelled by f i ve elements per ha l f span predict lower upper snap values than the corresponding so lut ion value fo r the ten elements per ha l f span arch. Quite l i k e l y the polygonal structures are s u f f i c i e n t l y d i f f e ren t for these two cases to account fo r th i s anomaly. One advantage possessed by the displacement incrementation so lut ion technique, compared with the load incrementation so lut ion method, i s the a b i l i t y to fol low the complete load def lect ion curve of a snap-through buckling problem. Several such curves fo r the arch for various i nc re -ment s izes and numbers of elements are shown in Figure 8.4. 3. Wright 's Ret iculated Shel l Segment D.T. Wright (13) in 1965, i n a paper concerning the design and s t a -b i l i t y of r e t i c u l a r domed st ructures , gave some theoret ica l ca lcu lat ions fo r snap-through buckling of a pa r t i cu l a r spherical she l l segment shown in Figure 8.5. He considered the behaviour of th i s segment under a ve r t i c a l load P at node A . In the der ivat ion of the load def lect ion curve he assumed that the she l l was shallow and that the fol lowing rat ios held: j _ ~ S i n — ^ f a n J z L t- J-r L r where L»- = the radius of curvature of the dome L = the length of the element shown in Figure 8.5. When he considered the nodes B to be pinned and ignored members BB which do not contribute any moment res i s tance, he arr ived at an upper bound so lut ion fo r snap-through buckling of the s h e l l . His theoret ica l so lut ion was: P = 3fi£ h' Ch'- h'z) f EI (h-h') • where K and h ' are defined i n Figure 8.5. D i f fe ren t i a t i ng Equation 8.1 with respect to h ' he found that snap through buckling would only occur when the fol lowing condit ion held: r 6 ^ .zc*3 h (8.2) where = , the radius of gyration of the members h = the r i s e of the element. When he considered the case where the j o i n t s B are allowed to t rans late and the members B e> are extens ib le , he arr ived at the lower bound so l u -t ion for the load-def lect ion curve of the she l l given by Equation 8.3. p _ 3 A _ e K ' Ch 2 - h ' 2 ) + E J f b - h ' ) ( 8 3 ) The c r i t i c a l values fo r th i s equation upon d i f f e r en t i a t i n g with re s -pect to h ' occurs when: ^ . 1 8 5 b (8.4) To ve r i f y Equations 8.1 and 8.3 the material properties of the toggle of Figure 8.1 were subst ituted into these equations. The radius of gyrat ion, rg , for th i s structure i s .072 inches. Thus both equations 8.2 and 8.4 are s a t i s f i e d since ^ . 1 0 5 2 fo r equation 8.2 ^ . o 7 - f f ° r equation 8.4 Therefore snap through buckling w i l l occur fo r the pinned structure and fo r the structure which allows members e>e> to extend. Buckling out of plane has been prevented as we are so le ly interested in ve r i f y ing Wright 's ca lcu la t ions . oil members have equal profx?rt*'es A - -Ofc.2 6 i n c h e s 2 I = . 0 0 0 3 0 ^ i nches ' * £ =• 3.o x i o 7 ps i P A j iranslct ion of node 0 FIGURE 8.5 WRIGHT'S SHELL SEGMENT 71 TABLE 5 WRIGHT'S SHELL SEGMENT Upper Lower Snap Snap B B Inextensible l b . l b . Wright 's formulation 99.9 48.5 12 'Secant ' elements 0.01" def lect ion increments 96.8 59.4 12 'Tangential ' elements 0.01" def lect ion increments 96.1 49.0 B B Extensible Wright 's formulation 75.2 73.4 12 'Secant ' elements 0.01" def lect ion increments * * 12 'Tangential 1 elements 0.01" def lect ion increments 71.8 71.6 * Solut ion does not exh ib i t maximum or minimum • • • vVn 'ght ' 5 formula 8.1 12 toocjent i o l e l e m e n t s per s t r u c t u r e , .oi" incr. .1 .2 .3 .<"» .5 .6 -7 V e r t i c a l d e f l e c t i o n ot node A ('inches) FIGURE 8.6 LOAD DEFLECTION CURVE OF WRIGHT'S SEGMENT WITH POINTS B FIXED AGAINST TRANSLATION «0 \HO. 0» o C +J o loo. 0 "40. 2o. W r i g h t s f o r m o l Q S-3 12 t a n g e n t i a l e l e m e n t s per '/<b s t r u c t u r e .01" incr. IZ s e c a n t e l e m e n t s pe r {/<o s t r u c t u r e .01 " incr. / / o. ,1 . 2 . 3 . 5 V e r t i c a l d e f l e c t i o n ct node A f i n c h e s ) .7 FIGURE 8.7 LOAD DEFLECTION CURVE OF WRIGHT'S SEGMENT WITH POINTS FREE TO TRANSLATE CO 74 By symmetry, we need consider only one s i x th of the she l l segment. Twelve elements were used to model th i s structure and def lec t ion i nc re -ments of 0.01 inches were used. Agreement with the formulas 8.1 and 8.3 was f a i r as shown in Table 5. The theoret ica l load-def lect ion curves and the experimental resu l t s are shown in Figures 8.6 and 8.7. 4. Three Dimensional Elbow In order to compare the accuracy of the ' t angen t i a l ' and ' secant ' elements, the structure shown in Figure 8.8, an elbow f i xed at both ends was subjected to an applied force increment. The actual applied forces were compared with those required to equ i l i b ra te calculated interna l forces (Equation 6.14: R = Q t £f ). F i r s t l y , the elbow was modelled by s ixteen elements and twenty v e r t i -cal force increments of -30 were applied at node 9. Secondly, twenty increments of moment about the X^ , axis of 60 were applied at node 9. The resu l t s given by Table 6 compare the known to ta l applied force with the forces ca lculated by equation 6.14 for various increment steps. Also shown are the def lect ions under the applied force and the vector of calculated forces at node 9 of the structure fo r increment step 20. A l l but the applied force component of th i s vector should be zero but the magnitude of the error when compared with the applied force i s smal l . Considering the s ize of the de f l e c t i on , 0.9 in a length of 4, and of the r o ta t i on , 0.265 radians in only twenty increment steps, the resu lts are considered reasonable. 75 x. (0 8 8 -£-3 X , s e c t lon A - A d> ' — © 6> <£> <Z> \<& i no de^ number d> element number <2> V'/ it. ^ 7 7 CQ CO i C o I s -A -E <5 l.o in* iorsion i . o in* s t r o n g b e n d i n g p J a n e /. o in u» e a k t e n d i n g p l a n e . I. O in z = l o o o o . k s J = 3 o o o . k s C FIGURE 8.8 THREE-DIMENSIONAL ELBOW TABLE 6 THREE DIMENSIONAL ELBOW Applied Ver t i ca l Force Increments Increment Step Total Applied Force Force Calculated by 'Tangent ia l 1 Elements 'Secant' Elements Force Deflection Force Def lect ion 5 - 150 - 150.03 - 0.2481 - 149.95 - 0.2492 10 - 300 - 300.18 - 0.4849 - 299.92 - 0.4953 15 - 450 - 450.32 - 0.7023 - 449.90 - 0.7343 20 - 600 - 600.69 - 0.8982 - 599.93 - 0.9633 5 =" § f°r node nine and increment step 20 'Tangential ' (- 4.44, - 9.53, - 600.69, 7.48, 3.79, - 1.18) 'Secant ' (- 4.51, - 10.20, - 599.93, 0.21, - 0.21, 0.08) Node forces ( P x, , Fx 2 , F* 3 , «\ x t } M x n i*lxj ) Applied Moment Increments Increment Step Total Applied Moment Moment Calculated by 'Tangent ia l ' Elements 'Secant ' Elements Moment Rotation Moment Rotation 5 300 299.95 . 0.0655 299.97 0.0655 10 600 599.85 0.1314 599.87 0.1312 15 900 899.74 0.1974 899.64 0.1977 20 1200 1199.60 0.2646 1199.15 0.2654 - ~ - _ for node nine 'Tangent ia l ' (- 0.79, - 0.12, - 0.02, 0.97, 1199.60, - 83.04) 'Secant ' (- 0.23, 0.32, - 0.08, 0.32, 1199.15, - 83.39) 5. Ring Dome The large r ing dome shown in Figure 8.9 represents a p rac t i ca l s t ruc -ture. By symmetry only ha l f of the dome was modelled. This helped reduce the band width of the problem immensely and meant that subdividing the structure into more elements would not ser ious ly increase th i s band width. The e l a s t i c and geometric properties of the r ing dome are as fo l lows: A l l members: I j = 144 i n . 4 ( tors ional moment of i ne r t i a ) T z = 720 i n . 4 (strong bending plane) I 3 = 48.0 in . 1 * (weak bending plane) Area = 4.0 i n . 2 Diameter of dome = 80 f t . Rise at Centre = 8 f t . Radius of Curvature = 104 f t . E = 3000 ksi C-r = 1200 ksi The ve r t i ca l plane was made the strong bending plane and the exter io r nodes were f i x ed . The structure was divided into 108 ' c a n t i l e v e r ' elements resu l t ing in 504 degrees of freedom and a band width of only 42. The ' t angen t i a l ' e l e -ments were chosen pr imar i ly because of cost. The p r inc ipa l intent was to show the load def lect ion behaviour of the dome. For t h i s , the upper r ing was subjected to 16 v e r t i c a l def lect ion increments of 0.25 feet under a un i -form ve r t i c a l load on the upper r ing. The load def lect ion curve and the configuration of a radius of the structure a f te r 16 increments are shown in Figure 8.10 and f igure 8.11. However no attempt has been made to 78 thoroughly analyze a l l the i n s t a b i l i t i e s of th i s s t ructure. Had another mode of f a i l u r e occured other than that shown, then the computer program would have at least indicated th i s when the s t i f f ne s s matrix went negative d e f i n i t e . Although there i s no so lut ion with which to compare these r e su l t s , they are se l f - cons i s tent and there i s every i nd ica t ion that the method was able to handle th i s comparatively large problem successfu l ly . FIGURE 8.9 RING DOME 80 40. r Uppcr ring vertical deflection ( f e e t ) FIGURE 8.10 LOAD DEFLECTION CURVE FOR RING DOME FIGURE 8.11 DEFLECTED CONFIGURATION OF A DOME RIB CHAPTER IX CONCLUSIONS A procedure was developed, using non-l inear s t i f fnes s matrices from Nathan (11), to fo l low the load def lect ion paths of shallow framed s t ruc -tures. The necessary transformation matrices and geometrical r e l a t i o n -ships were formulated for extending previous two-dimensional work to three dimensions. An incremental so lut ion technique was used and shown to be p r a c t i c a l . Two moving element coordinate systems and elements, the ' secant ' and ' t angen t i a l ' systems, were developed and re lated to a f i xed global system. The snap-through buckling paths of plane frames, arches and space frames were studied. The ' secant ' element was found to be more e f f i c i e n t for some large def lect ion problems in that fewer elements were required to reach a su i tab le so lu t i on . However, the ' t angen t i a l ' so lut ion required about one ha l f the computer time of the corresponding. 'secant ' s o lu t i on . Reducing the increment step s i ze was found more e f fec t i ve than increasing the number of elements used to model the s t ructure. 82 BIBLIOGRAPHY 1. Argyr i s , J .H . , "Continua and Discontinua", Proceedings of the Confer-ence on Matrix Methods in Structural Mechanics, Wright Patterson A i r Force Base, October 1965, A i r Force F l i gh t Dynamics Laboratory, Dayton, Ohio, AFFDL-TR-66-20. 2. Ko i ter , W.T., " E l a s t i c S t a b i l i t y and Post-Buckling Behaviour", Non-Linear Problems, Edited by R.E. Langer, Univers i ty of Wisconsin Press, Madison, Wisconsin, 1963. 3. Ma l l e t t , R.H. and Marcal, P.V., " F i n i t e Element Analysis of Nonlinear St ructures " , Journal of the Structural D i v i s i o n , Proceedings of the A.S.C.E., Vo l . 94, No. ST9, September 1968. 4. Jennings, A., "Frame Analysis inc luding Change of Geometry", Journal of the Structural Division, Proceedings of the A.S.C.E., Vol . 94, No. ST3, March 1968. 5. B r i t vec , S.J. and Ch i l ve r , A.H., " E l a s t i c Buckling of R ig id ly -Jo inted Braced Frames", Journal of the Engineering Mechanics D i v i s i o n , Proceedings of the A.S.C.E., Vo l . 89, No. EM6, December 1963. 6. Martin, H.C., "On the Derivation of S t i f fness Matrices for the Analysis of Large Deflect ion and S t a b i l i t y Problems", Proceedings of Conference on Matrix Methods in Structural Mechanics, Wright Patterson A i r Force Base, October 1965, A i r Force F l i gh t Dynamics Laboratory, Dayton, Ohio, pp. 697-715. 7. Supple, W.J., " I n i t i a l Post-Buckling Behaviour of a Class of E l a s t i c Structura l Systems", International Journal of Non-Linear Mechanics, Vol . 4, pp. 23-36, 1969. 8. Roorda, J . , " S t a b i l i t y of Structures with Small Imperfections", Journal of the Engineering Mechanics D i v i s i o n , Proceedings of the A.S.C.E., Vo l . 91, No. EMI, February 1965. 9. Roorda, J . , "The Buckling Behaviour of Imperfect Structural Systems", Journal of the Mechanics of Physics and Solids, Vo l . 13, pp. 267-280, 1965. 10. Wi l l iams, F.W., "An Approach to the Nonlinear Behaviour of the Members of a R ig id ly Jointed Plane Framework with F in i te Def lect ions " , Quarterly Journal of Mechanics and Applied Mathematics, Vo l . 17, 1964, pp. 451-469. 11. Nathan, N.D., " F i n i t e Element Formulation of Geometrically Non-Linear Problems," Ph.D. Dissertation, Univers i ty of Washington, Sea t t l e , Washington, March 1969. 83 12. Ebner, A.M. and Ucc i fer ro, J . J . , "A Theoretical and Numerical Compari-son of E l a s t i c Nonlinear F in i te Element Methods", National Symposium on Computerized Structural Analysis and Design, George Washington Univers i ty , Washington D.C., March 1972. 13. Wright, D.T., "Membrane Forces and Buckling in Ret iculated S h e l l s " , Journal of the Structural D i v i s i o n , Proceedings of the A.S.C.E., Vol. 91, No. ST1, February 1965. 14. Powel l , G.H., "Theory of Nonlinear E l a s t i c S t ructures " , Journal of the Structural Division, Proceedings of the A.S.C.E., Vo l . 95, No. ST12, pp. 2687-2701, 1969. 15. A rgy r i s , J .H . , "Recent Advances in Matrix Methods of Structura l Ana ly s i s " , Progress in Aeronautical Sciences, Vol. 4, MacMillan Co., New York, 1964. 16. Wilkinson, J .H . , The Algebraic Eigenvalue Problem, Oxford (1965). 84 T SNAPS 1! C PROGRAM TO FOLLOW THE LOAD DEFLECTION HISTORY CF FRAMED STRUCTURES 2 ; C ECTH CANTILEVER ANC SECANT TYPE ELEMENTS 3' I M P L I C I T R E A L * 0 ( A - H , 0 - 2 ) 4 i REAL *4 TITLE 5 C CIMENSICN ELCCK CNE AND TWO 6 COMMON /eL K1/ CN,F l ,P2,SMKG,NP 7 COMMON /B LK 2/ C4 , C5 , S4 , S5,C I 0,C11,SI 0 , SI 1 , P R E V , A 8 C CI MENS ION JOINT AMOUNTS CI MENS ION X ( 1 5C ) , Y ( 1 50 ) , 2 (150 ), J N ! 1 50 ) ,NC (15C ,6 ) 10 : C CIMENSICN MEMBER AMOUNTS 11 * , JNL (150 ),JNG( 15C) , N T Y P E t 1 5 0 ) , M L N i 1 5 0 ) 12 * , R R ( 1 5 0 , 1 2 ) , S F F ( 1 5 0 , 6 ) , S S ( 1 5 C , 6 ) , J L N ( 1 5 0 ) , S 4 C . ( 1 5 0 ) , C 4 0 ( 1 5 C ) •13" * , S 5 C i l 5 0 ) , C 5 C ! 1 5 C ) , S 6 0 ( 1 5 0 ) , C 6 C ( 1 5 C ) tXOR IG( 15 C ) , YORIG!150) 14 4 , Z G R I G ( 1 5 0 ) , T R { 150, 12) 15 C CIMENSICN L C £ C AMOUNTS 16 *,EV<600) , C R l e C G ) , C C L ! 1 5 0 ) , P C R I T { 1 5 0 ) 17 C CI MENS ION OTHERS 18 * , T I T L E ( 2 0 ) , R ( 6 ) , A A A 1 1 5 ) , E X ( 1 5 ) , E Y ( 1 5 ) , E Z ( 1 5 ) , E E E ( 1 5 ) , G G G ( 1 5 ) f E S ( 6 ) 19 * , F C C N S ( 6 ) , T C T A L ( 2 5 ) , L D E G < 2 5 ) , N P ! 1 2 ) , S M K 0 ( 6 , 6 ) , C E F L ! 3 ) 20 * , ECR( 12) ,ER{ 1 2 ) ,SMPREV(6,fc),EDS< 6) , EDSF{6 ) , A ( 6 ,12) , NDL(3) 21 * , V ( 1 2 , 6 , 1 2 ) , S V ( 1 2 , 1 2 ) , S M K 1 ( 6 , 6 ) , S S K 0 ( 1 2 , 1 2 ) , S S K 1 ( 1 2 , 1 2 ) , A T C ( 1 2 , 6 ) 22 * , A T l ( 1 2 , f c ) , T F V < 2 5 ) , E L K A { 4 5 ) , B L K B ( 1 1 6 ) , N l N C R ( 2 0 ) f S I N C 0 S ( 6 ) , X Y Z C H 3 ) 23 * , S K 0 ( 2 5 0 C C ) , S K I (25CCC) 24 ECUIVALENCE ( ELKA ! 1 ) , DM ), ( ELKC ( 1 ) , C 4 ) 25 CAT A NI NCR/1,2,3,6-, 1 1 , 2 1 , 3 1 , 4 1 , 5 1 , 6 1 , 7 1 , 8 1 , 9 1 , l O l t l l l , 1 2 1 , 151, 26 * 176 , 2 0 1 , 1CC1/ 27 C * * 28 1 FCRMAT(20A4) 29 2 FCR.-AT ( 2 C I 4 ) '30 3 FORMAT(7I4.3F1C.3) 31 4 FORMAT!6F1C .3 ) 32 23 F C R M A T ( 2 F 1 C . 3 , 2 I 4 ) 3 3 5 FCRMIAT(lHl) 34 6 FORMAT(/20X,• STRUCTURE NUMBER',14,» HAS',14,' JOINTS AND',14,' ME 35 4 M E E R S ' / 2 C X , • THERE ARE',14,' INCREMENTS A N D 1 , 1 4 , ' ELEMENT TYPES'/ 36 *2CX,' NBIF I S ' , 1 4 , ' N7 I S ' , 1 4 , ' ANC N'8 I S ' , 1 4 , ' AN D N E L MT IS ' , 14 37 *,/?0X,« (NELMT IS C MEANS CANTILEVER ELEMENTS, NON ZERO MEANS SECA 38 * NT ELEMENTS) • ) 39 7 FORMAT!/' JOINT KG. TEE 6 CEGREES CF FREEDOM THE I N I T I A L JOINT 40 * C O O R C I N A T E S ' / l l X , ' X Y Z MX MY MZ• , 1 OX,• X« , 9 X , • Y • , 9 X , • Z • 4 1 * ) .42 8 FCRNAT16X, 14, I X , 6 I 4 , 4X ,3 F 1 C . 3 ) ,.43 9 FORMAT!/' MEMBER J N L JNG JOINT IN X1-X3 PLANE MEMBER TYPE 44 * XP YP ZP') 4 5 1C F C R M A T ( 1 X , I 4 , 2 X , I 4 , 2 X , I 4 , 5 X , I 4 , 2 0 X , I 4 , 9 X , 3 D 1 2 . 4 ) 46 11 FORMAT!/' THE CI FFERENT ELEMENT TYFES: NUMBER AREA I ( X 1 ) 47 * I (X 2 ) I!X3 ) E G« ) 48 12 F C R M A K 3 1 X , I 4 . 6 C 1 2 . 4 ) 49 13 FORMAT!/ ' CE F L E C T I C N PARAMETER C I V I C E C BY LENGTH 0 F ' , F 1 0 . 3 , / ' THE 5C * PORTION OF THE EIGENVECTOR LSEO AS A DEFLECTION INCREMENT I S ' , E l 5 1 *0.4,' (THIS IS USED ONLY IF NBIF .NE . C)«) 52 14 FORM A T I / 5 C X , 'JOINT C E G R E E OF FREECfjM NUMBERS') 53 15 FORMAT(49X, I 4,4X,614) 54 16 FORMAT ( / ' T H E TOTAL NUMBER OF CEGREES O F FREEDOM =«,14) 55 17 F C 0 MA T( 14 , 6 F £ . C ) 5 6 18 FCRMAT(214,3FE.C) 57 19 FORMAT! ' THE MAXIMUM HALF BAND WIDTH INCLUDING DIAGONAL ',14/ 85 58 *• T t T A L NUMBER CF LOCATIONS IN S T I F F N E S S NUXNB',16) 5 9 20 FCRMAT(/54X,'I NCREMENT NUMEER ',15) 60 21 FORM AT ( / 'J 1 X, « MEMBER UATA«/50X,' MN JNL JNG IXX I YY , 6 1 * IZZ E G ARE A *) 6 2 22 F C P M A T ( 5 0 X , I 3 , 2 I 4 , 3 E 1 2 . 4 , 3 E 1 1 . 3 ) 63 24 F U R M A T( • ELEMENT TOTAL UI S P L A C E w c N T S R', 1 C X , t F 1 4 . 5 / , 4 0 X ,6F14 .5 ) 64 25 FORMAT(1CX,« ELEMENT TOTAL FORCES S ' , 6X, 6F 1 4 . 4 ) 65 2t FCRMATC R M 1 C CET = ',C14.5) 66 27 FCRMAT(41X, 'CEGREE OF FREEDOM TOTAL FORCE ' / 5 0 ( 4 5 X , 14 7X,D14 .5/) ) , 67 2 8 FORMAT('*',/2 I1C/20A4) 68 29 FCRMAT(6F12.7 ) 69 30 FCRMAT{• (THESE FCRCES ARE NORMALIZED TO A UNIT FORCE VECTOR IN ' , 70 * ' T H E APPLIED DEFLECTION INCREMENT M ETHOD) •) 7 1 31 FCBMAT(2CX,« APPLIED FORCE INCREMENT METHOD USED 1) 72 32 FCRMAT(2CX, * APPLIED CEFLECT ICN INCREMENT MET HOC USED' ) , 7 3 70 C FORMA T( • JOINT A P P L I E D FORCES*/2X , 1 4 , 2 X , 6 F 1 2 . 3 J 74 701 FCRMATC JOINT CCF NDEFL DE FL 1 , 2 I 4, 3 ( I 4 , G 1 2. 3) , G12 . 3 / ) 75 7C2 FCRMAT(2CX,' CESERVEC MEMBERS * ,4X ,2014) 76 703 FORMAT!* PCSITN NUMBERS FCR MEMBER',14,/IX, 1214) 77 801 FCRMAT(20X , * E G AR ST1 ST2 ST3 DM • , 1X,7G12.3,/2OX,• XM YM ZM P I 78 *P2 F3 P4' ,2X,7G12.3) . 79 802 F C R M A T ! 3 0 X ' ELEMENT SMKO ' , / 6 ( I X , 6 G 1 2 . 3 / ) ) eo 704 FCRMAT!2CX,' OBSERVED J O I N T S ' , 5 X , 2 0 I 4 ) 8 1 70 5 FCRMAT(30X,•STRUCTURE DATA'//' JOINT DEGREES OF FREEDOM X 82 * Y Z' ) 83 706 FOPMAT(7I4,3X,3D12.4) 84 707 F0RMAT(30X,» MEMBER NUMBER',15) , 85 7C6 FCRMAT( • ELEMENT DELTA R C ISPL . ', 12F9.5 ) 8 6' 7C9 FCRMAT!10 X t ' ELEMENT DELTA S D I S P L . ' , 7 X , 6 F 1 4 . 5 , / 1 0 X , ' TOTAL S e 87 * ' C I S P L . ' , 1 5 X,6F14.5) 88 804 FCRMAT(/50X, «A MATRIX',/ £ ( 11G 11. 3 , G 1 0 . 3 / ) ) 8 9 £05 FORMAT!//5CX»*TWELVE V MATRICES*,// 1 2 ( 6 ( L L G 1 1 .3,G10.3/)// ) ) - 9C 71 C FOPMAT(1CX,* ELEMENT INCREMENTAL FCRCES*,2X,6F14.4) 9 1 7 1 1 FORMAT(40X,•V*SFF MATR IX ' /, I 2 ! 1 I G 1 1 . 3 , G 1 0 . 3 / ) ) 92 712 FCRMAT(5CX,' K l STIFFNESS'/6(2CX,6C-13.5/ )) ' 93 806 FORMAK 5CX, • ELEMENT STIFFNESS KC+K1'/ 6 ( 20X , 6G1 3 . 5/) ) 94 71 3 FCRMAT!40X, ' AT*K0 *A ' , /12( 1X, 1 2( IPE I C . 3 ) / ) ) 9 5 714 FCRMAT ( 4 0 X , « £T*K1*A' , /1 2 (1 X , 1 2 ( 1? E 10 . 3 ) / ) ) „ 96 809 FORMAT! • SKC« I X , 1 IE 1 1.3) 97 81C FORMAT(• S K I ' , 1X,11E11.3) c 8 715 FCRMAT!' D R I NCR. ' , IX , 10(1PD12 . 3 ) / ( I C X , 1 C ( 1 PC 1 2 . 3 ) ) ) • 9 9 1CC1 READ(5 ,1 ,ENC = 1CCC) T I T L E ICO WR I T E ( 6 , 5 ) 1C 1 WRITE(6,1) T I T L E 1C2 READ(5 ,2)NRS,NJ,NM,NINCRT,NELMT,ITYPE,NBIF , N7,N8 ,NWRT1 t NWRT2, LOCI 1C3 WRITE!6,6) M R S , NJ ,NM ,N INCRT,ITYPE,NBIF,N7,N8,NFLMT 104 I F U O C I . N E . C ) K R I T E ( 6 , 3 1 ) V1C5 I F I L C D I . E C . C )V>RITE(6,32 ) 1 C6 c * * NRS IS THE STRLCTURE NUMEER 1G7 NJ IS THE NUMBER CF JOINTS 1C8 C * * NM IS THE NUMBER CF MEMBERS 1C9 c * * NINC R T IS THE NUMBER OF INCREMENTS 110 ITYPE IS THE NUMBER OF DIFFERENT ELEMENT TYPES 11 1 NEIF = 0 DO NCT FIND LOWEST BUCKLING LOAD 1 12 N 01 F NE C CALCULATE LOWEST 3UCKLINC LOAD 113 r »». c N7 .NE.O OUTPUT CN UNIT 7 N8.NE.0 OUTPUT ON UNIT 8 114 c * * NELMT = ZERC FCR CANTILEVER ELEMENTS 115 = NGN 7 E R C FCR SECANT ELEMENTS 116 O * NV<RT1 IS THE I N C R E M E N T N L ^ B E R AT fcHICH O L T P L T CN U N I T S 7 AND 8 117 C * * S T A R T S 86 118 C * * N V R T 2 IS T h e I N C R E M E N T NUMBER AT WHICH T H I S G U T P U T S T O P S 119 C * * L C D I = 0 F C R A F P L I E C D E F L E C T ION I N C R E M E N T M E T H C C -120 C * * L C D I = NON .0 FOR A P P L I E D F O R C E I N C R E M E N T METHOD 121 C * * R E A C IN T E E J O I N T I N F O R M A T I O N 122 DO 100 1 = 1, NJ 123 R F A C ( 5 , 3 ) J N ( I ) , (ND( I , K ) , K = 1 , 6 ) , X { I) , Y ( I ) ,2{ I ) 124 I C C C O N T I N U E 125 fcRITE(6,7) 126 WRIT E ( 6 » 6 ) ( J N ( l ) , ( N D ( I , K ) , K = l , 6 ) , > ( I ) , Y ( I ) , Z ( I ) , I = l , N J ) 127 C * * R E A C IN I N I T I A L MEMBER D A T A : J N L J NG J N P N T Y P E 126 C * * J N P IS A J O I N T E X I S T I N G IN T H E X 1 - X 3 P L A N E N E C E S S A R Y TO D E F I N E X3 A X I S 129 W R I T E ( 6 , 9 > 1 3 0 C * * 131 CC 101 1=1,NM 132 R E A C ( 5 , 2 ) J N L ( I ) , J N G ( I ) , J N P , N T Y P E ( I ) .1*3 3 R l l )=X( JNP ) - M JNL< I) ) -I 34 R ( 2 ) = Y ( J N P ) -Y { JNL ( I ) ) 1 3 5 R ( 3 ) = Z ( J N P ) - Z ( J N L ( I ) ) 1 3 6 R ( 4 ) = X ( J N G ( I ) ) -X ( J N L ( I ) ) 137 R ( 5 ) = Y ( J N G ( I ) ) -Y ( J N L ( I ) ) 138 R ( 6 ) = Z ( J N G ( I ) ) - Z ( J N L ( I ) ) 139 C A L L C R G S S ( R ( 4 ) , R ( 1 ) , E R ) 1 4 0 C A L L C R O S S ( E R , R ( 4 } , R ( 1 ) ) 141 XF = R ( 1 ) + X ( J N L U ) ) 142 YP = R ( 2 ) + Y ( J N L ( I ) ) 1 4 3 Z P = R ( 3 ) + Z ( J N L ( I ) ) H 4 4 W R I T E ( 6 , 1 C ) I , J N L ( I ) , J N G ( I ) , J N P ( N T Y P E ( I ) , X P , Y P , ZP 14 5 X C R I G ( I ) = X ( J N G ( 1 ) ) - X ( J N L { I ) ) . 1 4 6 Y Q R I G ( I ) = Y ( J N G ( I ) ) - Y ( J N L ( I ) ) 1 4 7 ZCR IG ( I ) = Z ( J N G ( I ) ) - Z ( J N L ( I ) ) 148 X F C = X P - X ( J N L ( I ) ) 149 Y F C = Y P - Y I J N L (I ) ) t l 50 Z P O = Z P - Z ( J N L ( I ) ) 151 C L Z = D S Q R T ( X P C * * 2 + Y P O * * 2 + Z P G * * 2 ) ""152 X1 = C S C R T ( X C R I G ( I ) * * 2 + Y C P IG ( I ) * * 2 + ZCP I G l I ) * * 2 ) i ' 53 XC1 = X C R I G (I ) / X 1 154 YO1 = YCR IG ( I ) / X I 155 Z C 1 = Z C P 1 G ( I ) / X l •*l5b XC2 = X P 0 / D L Z ,157 Y C 2 = Y P t ) / D L Z 158 Z C 2 = Z P 0 / D L Z 1-59 T 1= DSGRT (I . - X C 2 * * 2 ) 160 IF ( C A B S U 1 ) . L T . 1. C - 14) GOTO 1 4 5 161 C 4 C ( I )= Z 0 2 / T 1 -162 S4C ( I ) = - Y C 2 / T 1 163 C 5 C( I ) = T1 ^16 4 S 5 C ( I ) = X0 2 16 5 06 0 ( I ) = X 0 1 / T 1 166 I F ( C A B S ( Z 0 2 ) . L T . 1 . C - 1 2 ) G C T C 2 2 6 167 S 6 C I I ) = ( Y C I * T I * T 1 + Y C 2 * X 0 2 * X 0 1 ) / { Z 0 2 * T 1 ) 168 G C T C i l O l 169 226 I F ( C A B S ( Y C 2 ) . LT . 1 . C - 1 2 ) G 0 T C 2 2 7 170 S 6 C I I ) = - ( Z C 1 * T 1 * T l + Z G 2 * X G l * X G 2 ) / ( T 1 * Y C 2 ) 171 G C T U 1 0 1 17 2 22 7 S 6 C ( 1 ) = Y G 1 ; 173 G C T G 1 0 1 174 145 04 C I I ) = 1 . D C 175 S 4 C ( I )=O.CC 116 C 5 C ( I ) = 0.0C 177 S5C I I ) =xi: 2 178 C6C( I )=-X02 *ZC1 179 * S6C ( I ) = YG1 2 8 0 1C1 CCNT INUE 18 1 WR I T E ( 6 , 1 1 ) 182 CC 102 1=1, ITYPE 163 R E A C ( 5 ,4) * A M I ) , EX( I ) , E Y ( I ) , E Z ( I ) , E E E ( I ) ,GGG(I ) 184 WR I T E ( 6 , 12 ) I , A A A ( I ) ,EX( 1 ) , E Y ( I ) , E Z ( I ) , E E E ( I ) ,GGG(I) •1 8 5 102 CCNT INUE 1,8 6 VRITE OUT I N I T I A L CAT A 1 8 7 READ ( 5 , 2 3 ) PLCTL ,EIGFAC,JCINT,NDCF 188 WRITE(6,13) PLCTL,EIGFAC 189 REAC IN THE DEFLECTICN PARAMETER AND LCAC CONFIGURATION 19C CARE SHOULD BE EXERCISED IN CHOOSING THE DEFLECTICN PARAMETER 19 1 c * * NLPTS IS THE NUMBER OF JOINTS WHERE LOADS ARE APPLIED 19 2 1,9 3 R E A C ( 5 , 2 ) NLPTS 194 WR I T E ( 6 , 1 4 ) 195 I F I N B I F . N E . G ) N B I F = 1 196 NU=1 197 DC 104 J=1,NJ 198 DO 105 1=1,6 199 IF (N D ( J »I)-1 ) 2C0,201,202 2CC , 20C NC ( J , I)=0 2C1 GO TO 105 202 20 1 NCI J , I)=NU •'2X3 NU=NU+1 204 GO TO 10 5 2C5 202 NN=ND(J,I) 2Cb i N C ( J , I ) = N 0 ( N N , I ) 2 07 1C5 CCNTINUE 208 WRITE(6,15) J N ( J ) , ( N D ( J , K ) , K = 1 , 6 ) 2C9 104 CCNTINUE ,2,10 NU=NU-1 "211 WPITE(6,1£) NU '212 INCEX=0 3.13 I F(JCINT.NE.C.ANC.NCCF.NE .0 ) INCEX=NC(JC INT,NCCF ) 2 14 N=l 215 NINCRT=NINCRT + 1 TF=C. CO 2 17 N G C N E = 0 ••2 18 NN=C 249 DO 103 J = l , N L P T S 220 R E A C ( 5 , 1 7 ) J C INT, (FCONS (K ),K=1,6) 22 1 DC 137 K = l , 6 2 2 2 M = NO(JO INT,K) 223 IF(M.EQ.0 JGCTC137 2 24 I F ( F C C N S ( K ) . E G . C . C O ) G C T 0 1 3 7 225 N N = N N + 1 226 LCEG(NN)=M 227 I F ( I N D E X .EC .C ) INCEX = LDEG(1 ) 228 TFV(NN)=FCCNS(K) 229 T CT AL(NN ) = C.CG 230 TF = TF + TFV(NN )**2 2.31 137 CCNT INUE 222 WPIT E ( 6 , 7 0 C ) J CIN T , (FCCN S(K ) ,K= 1 , 6 ) 223 103 CCNTINUE 234 NLFTS=NN 235 I F ( L C D I . N E . C ) G C T C 1 4 4 2 36 TF=OSCRT(TF) 237 DC 138 J=1,NLPTS 88 23E T F V ( J ) = T F V ( J ) / T F 239' 136 CCNTINUE 2 4 0 W R I T E ( 6 , 2 C ) 241i 144 REAC(5,18> J C I NT , N CD F , ( DE F L ( K ) , K = 1, 3 ) 24 2' C** NCCF = 1 I F THE DEFLECT IONS APPLIEC APE TRANSLATIONS 243 C** NCGF = 4 IF THE DEFLECTIONS APPILED ARE ROTATICNS 244 K=C 245 224 JNT=ND(JCINT.NDCF+K ) "46 IF ( JNT.EC.CIGCTC223 2471-, NCL ( K+1)=JNT 2 48 K=K+1 249' I F ( K . E C . 3 ) G C T C 2 2 5 2 5C GCTC224 251 223 NCL(K+1)=1 2 52 D E F L ( K + l )=C.C 25 3 ! K=K+1 254' I F I K . G E . 3 ) G 0 T 0 2 2 5 2 5 5 ; GCTG224 256' 225 C E F L D = D E F L ( 1 ) * * 2 + C E F L ( 2 ) * * 2 + D E F L < 3 ) * * 2 2 57 K R I T E ( 6 , 7 C 1 ) J C I N T , N C G F , ( N D L ( K ) , D E F L ( K ) , K = 1 , 3 ) 258' 2C4 N19=N+19 259' READ(5,2) ( M L N ( J ) ,J=N,N19) 2'fcC C REAC IN THOSE MEMBERS THAT YCU fel SH TC SEE A CC M PL E TE LOAD-DEFLECT I ON 261 C HISTORY 26 2 I F ( M L N ( N ) .EG.C> GO TO 20 3 263 I F ( M L N ( N ) . N E . C . A N C . N 7 . N E . 0 ) W R I T E ( 7 , 7 0 2 ) ( M L N ( J ) , J = N,N19) .264 N=N+20 265 . GC TO 204 266 . 203 CCNTINUE 267 C REAC IN THOSE J C I N T S THAT YCL WISH TC FOLLOW THROUGH THE LOAD CEFLECTI 268 C F ISTORY 269 N=l 27C 2C6 M 9=N + 19 271 REAC(5,2) ( J L N ( J ) , J = K«Kl9 ) '272 I F ( J L N ( N ) . E Q .C ) GC TO 205 27 3 I F ( J L N ( N ) . N E . 0 .ANC.N7.NE .0 )WRITE(7 , 7 0 4 ) { J L N ( J ) , J = N,N19) 2.7 4 N=N+20 275 GC TO 206 ,276 205 CONTINUE 27 7 REWIND 1 278 REWIND 2 279 REWIND 3 2'3 0 NG0NE=O 281 NB=C 282 PCR = DABS(TFV( 1 ) J/TFV ( 1 ) 28 3 IF(N8.EQ .C )WRITE(6 , 2 1 ) 284 C CALCULATE POSITION NUMBERS NP 285 DC 106 K=1,NN 2-36 I = J M L ( K ) 287 J=JNG ( K ) 28 8 NF( 1)=ND( I , 1 ) 289 NP(2 ) =ND(I,2 ) 29 0 N F ( 3 )=ND(1,3) 2S1 N P ( 4 ) = N D ( 1 , 4 ) 29 2 N P ( 5 ) =ND(1,5) 2 9 3 N F ( 6 ) =ND ( 1,6 ) 29 4 NF(7 ) =ND(J, I ) 29 5 NF(B) = ND(J,2 ) 296 N P ( 9 ) = N D ( J , 3 > 2 97 NF ( 10 ) =NC(J,4 ) 89 2c.k -KF( 11)=NC ( J , 5 ) 299 NP ( 1 2 )=!Sli: ( J , (. ) 3C0 L = NTYPE(K ) 2C1 IF(N7.NE.C ) U«ITE(7,7C3) K,NP 2C2 E = E E E ( L ) 3G3 G=GGG(L) 2C4 A R = A AA(L ) 2C5 S T l ^ E X ( L ) .3C6 ST2=EY(L) 307 S T 3 = E Z ( L ) 3C8 I F J G . E C . C )G = E/2.6 2C9 IF ( S T l . E ' J . C . ) ST1 = ST2 + ST3 3 10 XM = X ( J ) -X( I ) 311 Y f - Y ( J ) - Y ( I ) 212 Z M = Z ( J ) - Z ( I ) 2 13 CM=DSORT(XM**2+YK**2+ZM**2) '314 P1=AR*E/CM 315 P2=12.*E*ST3/(CN**3) 316 P3=12.*E*ST2/(CK**3) -317 P4=G*ST1/CM ,318 CALL ELMSTC(P 1,P2 ,F3 ,P4 ,SMK0,Cf,NELKT) 3 19 W R I T E ( 1 ) E L K A '32C C CALCULATE C L A S S I C A L STIFFNESS FOR THE ELEMENT 321 C CALCULATE MEMBER BANC WIDTH 222 11=0 -2 2 3 J 1=1000 ,324 DC 107 L=1,12 325 N=NP(L) "226 I F (N.EQ.O ) GC TC 107 3*27 I F ( I l . L T . N ) l l = N 328 I F ( J l . G T . K ) J1=N 329 101 CONTINUE .220 NB1=I1-J1 + 1 ' 331 I F ( N B . L T .NE1 )NB=NB1 -2 22 I F ( N 8 . N E . C ) W P I T E (8,8C1 ) E,G,AR,ST1,ST2,ST 3,CM,XM,YM,ZM,PI,P2,P3,P4 ,,2 33 IF (N8.NE.0 )WR ITE (6,8C2) SM.KO :334 IFJN8.EQ.0 1WRITE (6,22) K , I , J , ST 1, S T 2 , S T3 , E , G , AR 22 5 l C d CCNTINUE .3 36 NTCTAL=Nl*NB ^237 WRITE(6,19) NB,NT CTA L 338 CC 108 K=1,NL -.3 39 DR(K)=0. 3 40 106 CCNTINUE 341 IUNIT=2 ,34 2 J L M T = 3 24 3 NSEVEN=N7 344 N EIG HT =N 8 .34 5 IF(NWRT2 .EC .C )NWRT 2=N INCRT 3 46 D C 1 4 6 J = 1 , 6 347 E R ( J ) = 0 . 348 E R ( J + 6)=C . 3 4 9 E S ( J ) = 0 . D C 250 146 CCNTINUE .3 51 DO 110 J = 1,12 3 52 CC 110 K=1,NM 353 T R(K,J)=C.CC 354 RR(K,J)=C.CC 355 IF (J.GT .6 ) GC TC 110 2 56 S F F ( K , J J = C . r C 357 S S ( K , J ) = C . C C 90 ! 358 11C ^ CCNTINUE i 359 CC 109 JK = 1,MNCRT : 1 6 0 JKM=JK-1 ?6 1 N8-=0 36 2 N7=0 26 3 IF{I\WRT1.LE.JKM.ANC.NWPT2.GE.JKM,)N7 = NSEVEN 3,64 IF (NWRT1 .LE.JKM.AND.NWPT2.Gt'.JKM)N8 = NE IGHT 3 6 5 I F(JKM.E0.C . dND.N7.E C.0)G0TD229 '•3-66 WRITE(6,2C) JKM 2 6 7 ' 229 I T £ K P = I U M T .r36 8 IUN I T = JUN I T 369 JUNIT=ITEMP ;2 7C M N C = 0 371 CC 111 L=1,2C •27 2 I FIN INCR (L ) .EO.JK )NINC = 1 2.73 I F l N I N C . E C . l ) GC TC 207 .274 111 CCNTINUE 375 207 I F ( N 7 . N E . C ) W R I T E ( 6 , 7 C 5 ) .376 IFCN7.NE.0 1GCTC221 3 77 I F ( J K K . N E . C . / N D . M h C . E C . l ) V > R I T E ( 6 , 7 0 5 ) 378 IF{ JKM.NE.C.£ND.MNC.EC.O.AND.JLN( 1) .NE . 0 ) V.R ITE ( 6 ,70 5 ) 3r79 231 DC 112 JL=1,NJ 3c.eC DC 113 J = l , 3 381 NN=ND(JL,J) 382 IF(NN.EQ.C) GO TO 113 383 GC TO ( 2 0 8 , 2 C 9 , 2 1 0 ) , J -384 2Cc X( J L ) = X ( J D + C R l N N ) 38 5 GC TO 112 386 209 Y ( J L ) = Y ( J L ) + C R ( K N ) 387 GC TO 113 288 21C Z ( J L ) = Z ( J L ) + CR(NN ) 289 113 CONTINUE 390 IJT=0 ,3 91 I F ( J K M . E C . C . * N D . N7 . E C . 0 ) CC TO 112 '292 DC 114 L=1,5CC 593 I F ( J L N ( L ) . E O . C ) G C TO 211 2'94 I F ( J L N ( L ) . E C . J L ) 1 J T = 1 295 IF ( I J T . E Q . l ) G C TC 21 1 •2 9 6 114 CONTINUE ,397 211 I F ( I J T . E C . l .CR.N7.NE.0 .OR .NINC . EQ. 1 )WR I T E ( 6 , 7 C 6 ) J N ( J L ) , 2<8 4 ( N C < J L , L ) , L = 1 , 6 ) , X ( J L ) , Y ( J L ) , Z ( J L ) 399 112 CCNTINUE 4 t 0 DC 115 I=1,NTCTAL 4C1 S K C ( I ) = 0 . J4C2 S K 1 ( I ) = 0 -4C3 115 CCNTINUE 4 04 REWIND 1 4 C 5 REWIND 2 4C6 REWIND 3 4 07 DO 116 M L = 1, N M '4 08 R E A C ( l ) ELK A 4 C 9 I f N = 0 410 DC 117 L=1,5 C C 4 11 I F I N L N I D . E C . O G O TO 212 412 I F ( M L N ( L ) . F C . N L ) I M N = 1 4 13 IF(IMN.EC . I )GC TC 212 4 14 117 CCNTINUE 415 212 IF{JKM.EC.C.*NC.N7.EC .0 )GCTQ228 4 16 I F { IMN .EC . L . CR. N7. NE .0 .CR .M NC.NE.O ) WR ITT ( 6 , 707 ) ML • 4 17 226 CO 120 J = l , 6 9 1 4 18 E C R ( J ) = 0 . C 419 ECP(J + 6)=C.C 42 0 E C S ( J ) = 0 . 421 E C S F ( J ) = 0 . 422 12C CCNTINUE H23 DC 118 J = l , 1 2 424 NS=NP(J) 425 I F ( N S . L E . O JGCTCl 18 ^26 E C R ( J ) = D P ( N S ) 4 2 7 , T R ( M L , J ) = T R ( M L , J ) + E C R < J ) 428 118 CCNTINUE 429 I F ( J K M . E C .C1GCTC2 15 430 IF (N7 .NE .C .CR. I l*N.NE. C) WR H E ( 6,7C8> ( EDR( J ) ,J=1, 12) 431 R E A C U U N I T ) ELKB 432 CC 119 K = l , 6 4 33 DC 119 J=1,12 434 E C S ( K ) = E C S { K ) + A ( K , J ) * E D R ( J ) 435 119 CCNTINUE 436 DC 121 J='l , 6 437 SS(ML,J)=SS(ML»J)+EDS(J) 433 ES(J)=SS(ML»J) 439 121 CCNTINUE AfiO T5=EDR(5) 441 E C R ( 5 ) = T 5 * C 4 + E C R ( 6 ) * S 4 442 ECR(6)=ECR(4 )*S5-T5=>S4*C5 + EDR ( 6 ) * C 4 * C 5 443 T5=EDR(11) 444 E C R ( 1 1 ) = T 5 * C 1 C + E 0 R ( 1 2 ) * S 1 0 445 ECR ( 1 2 ) = EDR(1C)*S11-T5*S10*C11+EDR{ 12)*C10*C 1 1 446 ' C C 1 2 2 J = 1,12 447 P R ( M L , J ) = R R ( f , L , J ) + E D P ( J ) 448 E R ( J ) = R R ( N L , J ) 449 122 CCNTINUE .4 50 IF(N7.NE.C.CR.IMN.NE.0.0R.MNC.NE.C)WRITE(6,7C9) EDS,ES 45 1 IF (N 7.NE.C .OR . IMN .EC . 1 .OR.NI NC. NE. C) WRITE ( 6 , 24) ( T R ( M L , J I ) , J I = 1 , 12) 4 5 2 CC 123 J = l , 6 453 CC 124 K=l,6 454 124 E C S F t J ) = S N P R E V ( J , K ) * E D S ( K ) + E D S F ( J ) 455 SFF(ML , J ) = S F F { M L , J ) + E O S F ( J ) ,4.5 6 122 CCNTINUE 4 57 I F { I M N . E G . 1 . C R . N 7 . N E . C . 0 R . M N C . N E . C ) W R I T E ( 6 , 2 5 ) ( S F F { M L , J ) , J = 1 , 6 ) 458 IF(N7.NE.O) W R I T t ( 6 , 7 10) EDSF 459 215 X Y Z C ( 1 ) = X C R I G ( M L ) 4-60 X Y Z C ( 2) = YDR IG (ML) 461 XYZ0(3)=ZORIG(ML ) 462 SI NCOS(1 )=S4C (ML ) 463 S I N C C S ( 2 ) = C 4 C ( M L ) '464 S INCOS ( 3 ) = S5C (ML ) 465 SIN C 0 S ( 4 ) = C 5 C(ML) 466 SI NCOS(5)=S6C(ML) 4 67 S I N C 0 S ( 6 ) = C 6 C ( M L ) .468 IFINELMT.EQ.C JCALL A VMATC (\», ER , XYZC , S INCOS ) 469 IF(NELMT.NE .CJCALL AVMATS(V,ER,XYZC,SI NCOS) 47C CALL TRANS ( X ) 4 7 1 I F ( N 7 . N E . 0 ) WRITE ( 7 , 8 0 4 ) ( ( A ( J,K) , K = I , 12 ) ,J=1,6) 47 2 I F ( N 8 . N E . 0 ) W R I T E ( 8 , 8 0 5 ) ( ( ( V ( K , J , I ) , K = 1 , 1 2 ) , J = 1 , 6 ) , I = 1 , 1 2 ) 473 DC 125 K=l,12 474 DC 126 J = l , 1 2 475 S V ( J , K ) = 0 . C 4 7 6 DC 1 2 7 1 = 1 , 6 4 77 . S V - ( J , K ) = S V ( J , K ) +V( K , I , J ) * S F F ( K L , I ) 47 8 127 C C N T I N U E 4 7 9 126 CCNT INUE 4 P C 125 C C N T I N U E 48 1 I F ( N 7 . N E . 0 ) W R I T E ( 7 , 7 l l ) ( ( S V ( J , K ) , K = l , 1 2 ) , J = l , 1 2 ) ^82 C A L L E L M S T 1 ( P 1 , P 2,DM , E S , S V K 1 , N E L M T ) 48 3 DC 128 1 = 1 , 6 ^84 DC 128 J = l , 6 ' i B 5 S M P R E V ( I , J ) = S P K C ( 1 , J ) + SMK1 (I , J ) 1 8 6 12 8 CCNT INUE 4E7 I F ( N 7 . M E . C ) V » R I T E ( 7 , 7 1 2 ) ( ( S f K l ( I t J ) , I = l , 6 ) , J = l f 6 J '-:£8 I F ( N 8 . N E . C ) V s R I T E ( 8 , 8 C 6 ) { { S N P R E V ( I , J ) , I = l , 6 ) , J = l , 6 ) 4 8 9 W R I T E ( J U N I T ) 8 L K B 49 C CC 1 2 9 L = l , 1 2 49 1 DC 129 K = l , 1 2 ft 9 2 S S K C ( K , L ) = C . 49 3 S S K 1 ( K , L ) = C . 4 94 I F I L - 6 ) 1 3 0 , 1 3 C , 1 2 9 4 9 5 13C AT C ( K , L ) = C . 4 9 6 A T 1 ( K , L ) = 0 . 4 97 129 CCNT INUE 49 8 DC 131 K = l , 1 2 4 9 9 DQ 131 J = l , 6 5 C C DC 132 1 = 1 , 6 5C 1 A T C ( K , J ) = A T C ( K , J ) + A ( I , K ) * S f K O ( I , J > 5G2 A T l ( K , J ) = A T l ( K , J ) + A ( I , K ) * S M K l ( I , J > 5 C 3 132 CCNT INUE 5C4 131 C C N T I N U E 5 C 5 DC 1 3 3 K = l , 1 2 5 C 6 DQ 133 J = l , 1 2 5 C 7 DC 1 3 4 1 = 1 , 6 5C8 S S K C ( K , J ) = S S K C ( K , J ) + A T O ( K , I ) * / S ( I , J ) 5 C 9 S S K 1 ( K , J ) = S S K 1 ( K , J ) + A T 1 ( K , I ) * A ( I , J ) .5-10 1 3 4 CCNT INUE 5 1 1 1 3 3 C C N T I N U E 5 12 I F I N 7 . N E . O ) W R I T E ( 7 , 7 1 3 ) ( ( S S K 0 ( K , J ) , K = 1 , 1 2 ) , J = 1 , 12 ) 5-13 I F ( N 7 . N E . 0 ) W R I T E I 7 , 7 1 4 ) ( ( S S K l l K, J ) ,K = 1 , 1 2 ) , J = 1 , 12) 5:14 DC 135 L = 1 , 1 2 5 15 DC 136 M = 1 , 1 2 5,1 6 M 1 = NP (L ) 5 17 N 2 = N P ( M ) 5 18 I F ( N 1 . E Q . C ) G C TC 135 5 19 IF ( M 2 . E Q . C ) G C T C 136 5 2 0 I F ( M 1 . C T . M 2 ) G C T C 136 52 1 K= ( M l - 1 ) * ( N E - 1 ) + r*2 5 2 2 S K C ( K ) = S K C ( K ) + S S K O ( L , M ) 5 2 3 S K 1 ( K ) = S K 1 ( K ) - S S K 1 ( L , M ) - S V ( L , M ) 5 24 136 C C N T I N U E 5 2 5 1 3 5 C C N T I N U E 5 2 6 116 CCNT INUE 5 2 7 IF ( N S . ' J E . O ) W R I T E ( 8 , 8 0 9 ) ( SK 0 ( J I ) , J I = 1, NTOT AL ) •5 2S I F ( N 8 . N E . C ) W R I T E ( 8 , 3 1 0 ) ( S K 1 ( J I ) , J I = 1 , N T C T A L ) 5 2 9 I F ( N B I F . N E . C . A N C . J K . L E . 3 ) C £ L L C R I T L D ( S K 0 , S K 1 , N U , NB , J K , D R , 5 3 0 * E V , E I G F A C , P C R ) 531 I F ( J K . G T . l . A N C . N C I F . N E . O . AND . J K . L E .3 )GC TC 2 1 6 5 3 2 DC 141 J = 1 , N T C T A L 5 3 3 S K C ( J ) = S K O ( J ) - S K 1 ( J ) 5 3 4 S K 1 ( J ) = S K C ( J ) 5 3 5 141 C C N T I N U E 92 536 219 CG 139 J = 1 , N U 5 37 C R I J)=0 . 5 3 8 139 C C N r INUE 539 DC 140 J=I,NLPTS -'54 0 D R(L 0 E G ( J ) )= T F V ( J ) 54 1 14C CCNT INUE 5.4 2 CET=1.D-14 543 NSCALE=1 544 IF(NGCNE.EQ.C)CALL DFBAND(SKO,CR,NU,NB,1,DET,AET,JET,NSCALE) 545 IF ( N G C i \ E . E C . C . A N C . C E T . L T . l .D-14 )G0 TC 217 J.S4 6 IF(KGCNE) ?2C,22C,218 5 47 217 NGQNE=l .548 G C TO 219 549 218 CALL DBANDI (SKI,CR,NU ,NB , i,DET ) .550 22C IF(JKM.EQ.C.AND.N7.EC.C JG0T023C 55 1 WRITE(6,26 )CET '552 23C C I V = ( D R ( N D L ( 1 ) ) * C E F L ( 1 ) + D R ( N D L ( 2 ) ) 4 C E F L ( 2 ) + C R { N C L ( 3 ) ) 553 * * C E F L ( 3 ) J/CEFLD 5 54 CIV = 1 ./D IV 555 I F( N B I F . N E . C . A N C . J K . E C . l ) C IV=1. 556 1 F ( L 0 D I . N E . 0 )CIV=1.0 5 57 DO 142 J=1,NU -558 C R ( J ) = D R ( J ) * C I V 559 142 CCNTINUE 5^60 222 DO 143 J = l ,NLPTS 561 TCTAL(J)=TCT£L(J ) + T F V ( J ) * C I V 562 143 CCNTINUE 56 3 I F ( JKM.Eli .0 JGCTC232 ^564 JKL = JK-NB IF t565 JKR=JKL-1 5.66 I F ( J K R . E C . C ) J K R = 1 567 P C R I T ( J K L ) = TCTAL( D/PCR 56 8 D C L ( J K L ) = 0 A B S ( D R ( INDEX)/PLOTL ) + D C L ( J K R ) 569 GCTC233 -570 232 P C R I T ( 1 ) = T C T / ? L ( 1 ) / P C P 571 DOL(1)=DABS(CR(INDEX )/PLOTL) : 5 7 2 233 WRITE(6,27) (LDEC- ( J ) ,TOTAL ( J ) »J=1,NLPTS) ^57 3 221 I F(N7.NE.C)WRITE ( 7 , 7 1 5 ) ( C R ( J I ) , J I = 1 , N U ) 5-74 GC TO 109 575 216 C IV = PCR ?76 I F ( J K . E C . 3 ) C I V = C . 577 GC TO 222 " 578 109 CCNT INUE 579 WR1TE(0,28) N R S , J K L , T I T L E 58 C W R I T E ( C , 2 9 ) ( P C R I T ( l ) , C C L ( I ) , I = l , J K L ) 58 1 GC TO 1001 J 5 8 2 10C0 STOP 5 8 3 ENC 584 SUBROUTINE TRANS(V) 58 5 I M P L I C I T R E A L * 8 ( A - H , 0 - Z ) 586 CCMMON /ELK2/C4,C5,S4,S5,C10,C11,S10,S11,SM,A 567 CI MENS ION SM ( t,6 ) ,A(6 ,12 ) , V( 12,6 , 12 ) ,G(12,12 ) >588 * , B 1 ( 3 , 3 ) , C ( 3 , 3 ) , T ( 6 , 1 2 ) , U ( 6 , 1 2 ) , B ( 6 , 1 2 ) , T M ( 3 , 3 ) 589 CC 1 J= 1 , 12 59C CC 2 K = l ,12 591 G(K,J)=O.DC 592 I F ( K . E Q . J ) G ( K , J ) = 1 . D 0 593 I F ( K . C T . 6 JGCTC2 594 B ( K , J ) = A ( K , J ) 595 A (K , J ) =0.DC 93 5 9 6 .7 C C N T I N U E 597 1 C C N T I N U E 5 , 8 , G ( 5 , 5 ) = C 4 5 9 9 G ( 6 , 6 ) = C 4 * C 5 6 C C G ( 5 , 6 ) = S4 6 0 1 G ( 6 , 5 ) = - S 4 * C 5 6C2 G ( 6 , 4 ) = S 5 6 C 3 G ( l l , l i ) = C l C 6 C 4 G ( 1 2 , 1 2 ) = C 1 1 * C 1 C 6 C 5 G ( 1 1 , 1 2 ) = S 1 C •6-C6 G ( 1 2 , 1 1 ) = - S 1 C * C 1 1 6 C 7 , G ( 1 2 , 1 0 ) = SI 1 6C8 DC 106 N = l , 7 , 6 6 0 9 DO 200 J = l , 3 6 1 0 CC 201 K = l , 3 6 1 1 B 1 ( K , J ) = G ( K + 2 + N , J + 2 + N ) '6 12 201 C C N T I N U E ,613 2 C C C C N T I N U E : 6 \ A CC 10 5 M = l , 4 , 3 6 1 5 DC 102 1 = 1 , 1 2 6 1 6 DC 100 J = l , 3 6 1 7 DC 101 K = l , 3 £ 1 8 T M K , J ) = V ( I , K + M - 1 , J+N + 2 ) 6.19 1C1 C C N T INUE 6t2 0 • I C C C O N T I N U E 6 2 1 DC 150 J = l , 3 6 2 2 DC 151 K = l , 3 6 2 3 C ( K , J ) = O . D C 6 2 4 DC 152 1 = 1 , 3 6 25 C ( K , J ) = C ( K , J ) - » T P ( K , L ) « e i ( L , 6 2 6 152 C C N T INUE 6 27 151 C O N T I N U E 6 2 8 15C C C N T I N U E 6 2 9 DC 1 0 3 J = l , 3 -630 DC 104 K = l , 3 6 3 1 104 V ( I , K + M - 1 , J + N + 2 ) = C ( K , J ) 6 2 2 1 C 3 C C N T I N U E '6'3 3 102 C C N T I N U E 6<34 1C5 C C N T I N U E 6 3 5 1 0 6 C C N T I N U E -4,3 6 N N = 1 6 2 7 N l = 4 "638 N!2 = 6 6 3 9 9 9 CC 40 J = l , 3 6 4 C DC 41 K = l , 3 6 4 1 41 T M ( K , J ) = C . C C 6 4 2 4 0 C O N T I N U E 6 4 3 G C T 0 ( 5 0 , 6 0 , 7 C , 8 C ) , N N 6 4 4 5C T M 2 , 2 ) = - S 4 6 4 5 T M ( 2 , 3 ) = C 4 6 4 6 T M ( 3 , 2 ) = - C 4 ^ C 5 6 4 7 T M ( 3 , 3 ) = - S 4 * C 5 6 4 8 M=4 6 4 9 G 0 T 0 8 1 6 5 0 6C T M ( 3 , 1 ) = C 5 £ 5 I T M 3 , 2 ) = S 4 * S 5 6 52 T M ( 3 , 3 ) = - C 4 * S 5 6 5 3 M = 5 6 5 4 G C T C 8 1 6 5 5 7C T M ( 2 , 2 ) = - S L C 6 5 6 T M ( 2 , 3 ) = C 1 0 6 5 7 TN ( 3 , 2 ) = - C l C * C l 1 9 5 6 5 8 TK ( 3 , 3 ) = - S l C * C l l 6 5 9 M=1C 6 6 0 N 1 = 1 0 6 6 1 N2 =12 66 2 G C V 0 6 1 6 6 3 80 T M ( 3 , 1 ) = C 1 1 r 6 4 T M 3 , 2 ) = S 1 C * S 1 1 c'6 5 TM ( 3 , 3 ) = - C l C * S - l 1 6.66 M = l l 6 6 7 81 CC 9 0 J = l ,6 6 6 8 DC 9 2 K = N 1 , N 2 6 6 9 DO 9 1 L =N 1,N 2 6 7 0 91 A( J , K ) = A ( J , K ) + C ( J , L ) *TK. ( L - N 1 + 1 , K - M l + l ) 671 92 C C N T I N U E (672 9C C C N T I N U E 67 3 DO 9 3 K =N 1, N 2 6 7 A DC 9 4 J = l , 6 6 7 5 V ( M , J , K ) = V ( N , J , K ) + A ( J , K ) 6 7 6 A ( J , K ) = 0 . C C G 7 7 94 C C N T I N U E 6 7 8 93 C C N T I N U E 6 7 9 NN=NN-U 6*80 IF (NN . L E . 4 ) G C T 0 9 9 681 DC 25U K = l , 1 2 6 8 2 DC 251 J = l , 6 ff8 3 T ( J , K ) = V ( 5 , J , K ) 6 8 4 U ( J , K ) = V ( 1 1 , J , K ) 6 8 5 251 C C N T I N U E 6 8 6 ' 2 5 C C C N T INUE 6 8 7 DO 3 J = l , 6 6 8 8 DC 4 K = l , 1 2 6 8 9 V ( 4 , J , K ) = \ M 4 , J , K ) + V ( 6 , J , K ) * S 5 6 9 0 V ( 5 , J , K ) = T ( J , K ) * C 4 - V ( 6 , J , K ) * S 4 * C 5 69 1 V ( 6 , J , K ) = T ( J , K ) * S 4 + V ( 6 , J , K ) * C 4 « C 5 6 92 V ( 1 C , J , K ) = V ( 1 C , J , K ) + V ( 1 2 , J , K ) * S 1 1 6,93 V ( l l , J , K ) = U J , K ) * C 1 0 - V ( 1 2 , J , K ) * S 1 0 * C i l 6 9 4 V ( 1 2 , J , K ) = U J , K ) * S 1 0 + V ( 1 2 , J , K ) * C 1 0 * C 1 1 69 5 CC 5 1 = 1 , 1 2 ,696 A ( J , K ) = A ( J , K ) + B ( J , L ) * G ( L , K ) 6 97 5 C O N T I N U E 6 9 8 4 C C N T I N U E 6.S9 3 C C N T I N U E 7 C 0 R E T U R N 701 ENC 7 0 2 S U B R O U T I N E E L M S T C ( A , G , C , C , K , L , N N ) 7 C 3 R E A L M S A , B , C , C , L , K ( 6 , 6 ) 7C4 DC 1 1 = 1 , 6 ?,C5 CO 2 J = l , 6 7C6 K ( J , I ) = 0 . 7 C 7 2 C C N I I N U E 7 C 8 1 C C N T I N U E 7 C 9 I F ( N N . N E . G J G C T C 1 C 7 1C K ( 1 , 1 ) = A 7 11 K ( 2 , 2 ) = C 7 1 2 K ( ? ,t» ) = - . 5 * C * L 7 13 K ( 3 , * ) = B 7 14 M ? , 5 ) = . 5 * B * L 7 15 K ( 4 , s ) = 0 7 16 K ( 5 , 5 ) = e * L * L / 3 . 7 17 K (6 ,6 ) = C * L * L / 3 . 7 18 GOTO 11 -7 19 . 1C K ( 1 , 1 ) = A 7 20 K( 2 , 2 ) = B * L * L / 3 . 7 2 1 K ( 2 , 5 ) = K ( 2 , 2 ) / 2 . 7 2 2 K ( 3 » 3 ) = C * L * L / 3 . 7 2 3 K ( 3 , 6 ) = K ( 3 , 3 ) / 2 . 7 24 K ( 4 , 4 ) = D 7 2 5 K ( 5 , 5 ) = 8 * L * L / 3 . 72 6 K ( 6 , 6 ) = C * L * l / 3 . 7 2 7 11 CC 3 1 = 1 , 6 7 28 CO 3 J = I , 6 7 29 K ( J , I ) = K ( I , J ) 7 3 C 3 C C N T I N U E 731 R E T U R N .7 2 2 EN C 7 3 3 S U B R O U T I N E E L MST 1 (P 1, P 2 , DM, S , S MK 1, NN ) 7 2 4 R E A L * 8 P 1 , P 2 , C f , S C , S { 6 ) , SMK1 (6 ,6 ) 7 3 5 DC 2 1 = 1 , 6 7 3 6 DC 2 J = l , 6 '737 2 S P K K I , J ) = C . 7 3 8 IF ( N N . N E . 0 1 G C T C 1 C 7 3 9 ; S M K 1 ( 2 , 6 ) = - . 1 * P 1 * S ( 1) 7 4 0 j S M K 1 ( 3 , 3 ) = 1 . 2 * P 1 / D M * S ( 1 ) 7 4 1 ' S F . K K 3 , 5 ) = . 1 * P 1 * S ( 1 ) 7 4 2 S M K 1 ( 4 , 4 ) = P 2 * 0 M / 1 2 . * S ( 1 ) 7 4 3 S M K 1 ( 5 , 5 ) = 2 . / 1 5 . * P 1 <DM*S (1 ) 7 4 4 S M K 1 ( 6 , 6 ) = 2 . / 1 5 . * P 1 * D f * S (1 ) , 7 4 5 S M K 1 ( 2 , 2 ) = 1 . 2 * P 1 / D M * S ( 1 > 7 4 6 S M K H 1 , 3 ) = + . 1 * P 1 * S ( 5 ) +1 . 2 * P 1/DM*S ( 3 ) 7 4 7 S M K 1 ( 1 , 5 ) = + 2 . / 1 5 . * P 1 * C K * S (5 ) + . 1 * P 1 * S ( 3 ) 7 4 8 S M K K 2 , 4 ) = - 4 . 2 / 1 2 . * P 2 * D M * S ( 5 ) - . 5 * P 2 * S ( 3 ) 7 4 9 S M K 1 ( 4 , 6 ) = + 3 . 2 / 1 2 . * P 2 * D M * D M * S ( 5 ) + 4 . 2 / 1 2 . * P 2 * D M * S ( 3 ) 7 5 C S f K K l , 2 ) = - . l * P l * S ( 6 ) + 1 . 2 * P 1 / C M * S { 2 ) 7 5 1 S M K 1 ( 1 , 6 ) = - . 1 * P 1 * S { 2 ) + 2 . / 1 5 . * P 1 * C M * S ( 6 ) 7 5 2 S M K 1 ( 3 , 4 ) = 4 . 2 / 1 2 . * P 2 * D M * S ( 6 ) - , 5 * P 2 * S ( 2 ) •353 S M K 1 ( 4 , 5 ) = 3 . 2 / 1 2 . * F 2 * C M * D F ' * S ( 6 ) - 4 . 2 / 12 . * P 2 * CM*S ( 2 ) 7,54 G C T 0 2 0 0 7 5 5 10 S M K 1 ( 2 , 2 ) = 2 . / 1 5 . * P 1 * D M * S ( 1 ) ,756 SMK1 ( 3 , 3 ) = S I » K 1 ( 2 , 2 ) 7 5 7 S K K 1 ( 2 , 5 ) = - P 1 * C f * S ( 1 ) / 3 0 . 7 5 8 S M K 1 ( 3 , 6 ) = S t K 1 ( 2 * 5 ) 7 5 9 S M K 1 ( 4 , 4 ) = P 2 * C M / 1 2 . * S ( 1 ) 7,6'C S M K 1 ( 5 , 5 ) = 2 . / 1 5 . * F 1 * D f * S (1 ) 7 6 1 S M K 1 ( 6 , 6 ) = 2 . / 1 5 . * P 1 * D M * S ( 1 ) 7 6 2 S M K U 3 , 4 ) = . 8 / 1 2 . * P 2 =*D M*D M*S ( 2) + P 2 * D M * D M / 1 2 . * S I 5) 7 6 3 S M K I ( 4 , 6 ) = 3 . 2 / 1 2 . * F 2 * C M * D M * S ( 5 ) + P 2 * C M * C M / 1 2 . * S ( 2 ) 7 6 4 S M K 1 ( 1 , 3 ) = 2 . / 1 5 . * F l * D P * S ( 3 ) - P l * C M / 3 0 . * S ( 6 ) 7 6 5 S M K 1 ( 1 , 6 ) = 2 . / 1 5 . * P 1 * D M * S ( 6 ) - P 1 * 0 M / 3 C . * S I 3 ) 7 6 6 S M K 1 ( 2 , 4 ) = . 8 / 12 . *P 2 =*DM *DM *S ( 3 ) + P 2 * P M * D M / 1 2 . * S ( 6 ) 7 6 7 S M K 1 ( 4 , 5 ) = 3 . 2 / 1 2 . * F 2 * C M * D M * S ( 6 ) + P 2 < C M * D M / 1 2 . *S ( 3 ) 7 6 8 S M K 1 ( 1 , 2 ) = 2 . / 1 5 . * P l * D M * S ( 2 ) - P 1 * U * / 3 0 . * S ( 5 ) 7 6 9 S M K 1 ( 1 , 5 ) = 2 . / 15 , * P 1 * D M * S ( 5 ) - P 1 * D M / 3 C . * S { 2 ) 77C 2CC DC 3 1 = 1 , 6 77 1 DC 3 J = l , 6 7 7 2 S M K U J , I ) =SMK 1 ( I , J ) 773 3 C C N T I N U E 7 74 R E T U R N 77 5 EN C 7 76 SU E ROUTINE CRITLC(SK0,SK1,N , M,JK,DR, EV,EIGF AC,PCR) 777 I M P L I C I T REAL*8(A-H , C-Z ) 778 DIMENSION C R ( 1 ) , SKO( I ) ,SK I ( 1) ,EV( I) ,ST(25CC) 779 IF ( J K . E G . l )GC TC 1 -' 76 C IF < J K . E C . 3 ) G C TC 10 781 EFS=l.D -e -78 2 EPSV=l.D-4 7 63 IT=100 . 7 64 CALL DVPCWR(S KC , SSI , N , M ,M ,EV,N,EVALI ,1 ,E)»S, EPSV , I T , ST, COND ) ' 785 PCR=EVAL I ^'786 WRITE(6,4) EVALI 7 87 •> 4 F C R MAT(' C R I T I C A L LOAD I S ' , E 1 4 . 5 ) 788 CIV=0.0 789 DO 16 J=1,N t 79C EVR=DABS(EV(J ) ) 791 1 i DIV=DMAX1 (DIN. ,EVR) J 7 9 2 DO 2 J=1,N .793 2 E V ( J ) = E V ( J ) / C I V . 5794 WRITE(6,3) ( E V ( J I ) , J I = 1 , N ) 795 3 FORMiAT(//« EIGENVECTOR*// ICO ( IX , 1 2 G 1 0 . 3 / ) ) 796 DC 5 J = 1,N ' 797 c C R ( J ) = 0 R ( J ) * ( P C R - 1 . ) 798 1 RETURN ,799 10 DO 15 J=1,N t8C0 15 C R ( J ) = E V ( J ) * E I G F A C e c i RETURN 8C2 END "803 SUBROUTINE CRQSS(A,B,C) ,* 804 REAL #8 A,E,C > 8C5 DIMENSION A ( l ),E( 1) ,C(1 ) >8C6 C ( 1 )=A( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) 8 07 C ( 2 ) = A(3 ) * f i (1 )-A( 1 ) * B ( 3 ) 808 C ( 3 ) = A ( 1 ) * 0 ( 2 ) - A ( 2 ) * 3 ( 1) EC9 RETURN "'810 END , 811 SUBROUTINE DEANCI( A , B , N , K, LT , CET ) 812 I M P L I C I T REALMS(A-H , 0 - 2 ) "8 13 CCMMON /ZDE T/ CE,NCN '8 14 COMMON /ZCCN/ CCND 815 C I KENS ION A ( 1 ) , B ( 1 ) 1-8 16 DOUBLE PRECISION CSCRT,DAGS,DSIGN 8 17 I F (M .EQ . 1 ) GC TC ICC 818 MM=M.-1 6 19 NM=N*M 8 2 0 N M 1 = N M - M M 82 1 IF ( L T . N E . l ) GC TO 55 822 CCNC=1.0 823 NCN = C 824 DE = C . 8 2 5; MP=M+1 626 KK=2 827 FAC=DET 828 NRGW=1 829 I F ( A ( 1 J.EQ.C. ) GC TO 60 £ 3 C SLM=O.D0 '83 1 DC 77 I = 1,M 822 77 SUM=SUM + A( 1 )*A( I ) 8 33 S = 1 . O / O S C R T ( C A P S ( A ( 1 ) ) ) 8 34 CE=A( 1 ) 835 A ( 1 ) = S 97 8 3 6 IF ( A ( 1 ) . L T .0 . ) A ( 1 ) = - S £ 3 7 C C N D = C C N D / ( A ( 1 ) * A ( 1 ) * C S C R T ( S U M ) ) 8 3 8 SUM = A ( 2 ) * A ( 2 ) + . A ( M P ) * A ( M P ) 3"W M2=2*M £ 4 C B I G L = D A B S ( A ( 1 ) ) 8 *1 S M L =B IGL 8 4 2 MPP=MP+1 8 4 3 CC 88 I =MPP ,M2 8 4 4 68 SUM = SUM + A( I ) *A ( I ) '•5 A ( 2 ) = A ( 2 ) * D A E S ( A ( 1) ) ,;4b S = A ( M P ) - A ( 2 ) * A ( 2 ) * D S I G N ( 1 . D O , A ( 1 ) ) P<<7 I F ( S . N E . C ) GO TC 16 4 8 NRCW=2 8 4 9 GO TO 60 .850 16 A ( M P ) = 1 . C / C S G R T ( C A E S (S ) ) g-5l C E = 0 E * S >P5 2 CCND=COND / (A (MP) * A ( M P ) * D SGR T ( SUM ) ) 8 5 3 I F ( S . L T . O . ) A ( M P ) = - A ( M P ) 8 5 4 AAA = D A B S ( A ( N P ) ) 8 5 5 I F ( A A A . G T . B I G L ) fiIGL = A A A 8 5 6 I F ( A A A . L T . S M L ) SML=AAA 5 7 IF ( N . E Q . 2 ) GC T C 53 e 5 8 MP=MP+K P 5 9 DC 6 2 J = MP,NM 1, M 8,60! J P = J - M M 8 6 1 MZC=0 8 6 2 I F ( K K . G E . M ) GC T C 1 8 6 3 K K = K K+1 -364 1 1=1 8 5 5 J C = 1 816 6 GC TO 2 8 6 7 1 KK = KK + Mt 8 6 8 L I =KK —MM 8 6 9 J C = K K - N M S 7 C 2 CO 6 5 I = K K , J P , " M v . / 1 IF ( A ( I ) . EQ .C . ) GC TO 6 4 8 7 2 GC TO 66 •Kl 3 64 J C = JC+M &74 65 M,ZC = MZC + 1 8 7 5 A S U M 1 = 0 . 4E76 GC TO 61 8r>7 i t MMZC=MM*MZC 8 7 8 I 1=1 l + H Z C 8J9 KM=KK+MMZC 8"8C S U M = A ( K M ) * A ( K M ) 8 8 1 A (K M ) = A ( K M ) rC A 8 S ( A ( J C ) ) *8 8 2 ASUM1 = A (KM ) * A ( K M ) * D S I G N ( 1 . D O , A ( J C ) ) 8 8 3 I F ( K M . G E . J P > G C TC 6 8 8 4 K J = K M + V M " 3 5 J J = J C ^ 8 6 CO 5 I = K J , J P , M M 8 8 7 S U M = S L M + A ( I ) * A (I ) • 8 8 8 J J = J J + M 8 8 9 A SUM 2 = 0 . 8 9 0 IM=I -MM 89 1 I 1 = 1 1 + 1 £ 9 2 K I = I I + M M Z C 8 9 3 K Z = J C 894 CC 7 K =KM , 1M, f> M 8 9 5 ASUM2 = A S I N 2 + A < K I ) * A ( K ) * D S I G N ( 1 . C 0 , M K Z ) ) 8 9 6 KZ=KZ+M 8 9 7 7 K I = K I + K K 8 98 A ( I ) = ( A ( I ) - A S U f 2 ) * C A B S ( A ( K l ) ) 8^ 3 9 A S I V 1 = A S U H 1 + A( I ) * A ( I ) * D S I G N ( 1 . 0 0 , A ( J J ) ) sec 5 C C N T INUE CG I 6 C C N T I N U E 90 2 JMN=J+KV 9 0 3 CO 3 I = J , J M V 9C4 3 SUK=SUM+A( I ) * A ( I ) £•^5 61 S = A ( J ) - A S I M , r . C 6 I F ( S . N E . C . ) GC TC 63 r l NROW=(J + f N ) / " CQ GC TO 60 9 C 9 63 A ( J ) = 1 . 0 / D S Q R T ( C A B S ( S ) ) 9 10 CCN D = C O N D / ( A ( J ) * A ( J ) * D S O R T ( S U N ) ) 9 4 1 C E = C E * S 9 12 I F ( C A B S ( D E ) . G T . l . E - 1 5 ) GO TO 144 9.13 C E = D E * 1 . E + 1 5 9 1 4 N C N = N C N - 1 5 9 1 5 GC TO 1 4 5 9 16 144 I F ( D A 3 S I D E ) . L T . l . E + 1 5 ) GC TO 145 9 1 7 D E = C E * 1 . E - 1 5 9 18 NCN=NCN+15 9 19 145 C C N T I N U E r ^0 I F ( S . L T . C . C ) A ( J ) = - A ( J ) 21 AAA=DABS ( A ( J ) ) 922 I F ( A A A . G T . E I G L ) E I GL = A A A 9 2 3 I F ( A A A . L I . S N L ) SML = A A A 9 2 4 62 C C N T I N U E 9 2 5 GC TO 53 9,2.6 6C W R I T E ( 6 , 9 9 ) NRCV. 9 2 7 54 D E T = 0 . 9 2 8 99 FO RM A T ( 3 5 HC ERROR C O N D I T I O N E NCCUNT ERED IN R O W , 1 6 ) 9 2 9 R E T U R N 9 3 0 53 DET = S M L / 8 IGL 9 3 1 55 B ( 1 ) = B ( 1 ) * C A E S ( A ( 1 ) ) 9 3 2 KK=1 ?;3 3 K l = l 9-3 4 J = l 9 3 5 L 4 = l <9.36 CC 8 L = 2 , N 9^7 B S U N 1 = 0 . 9 3 8 - L M = L - 1 9 3 9 J = J+M 9 4 0 I F U K . G E . N )GC T C 12 9 4 1 K K = K K + 1 9 4 2 GC TO 13 9 * 3 12 KK=KK+M 944 K1=K1+L 4 5 L 4 = L 4 + M ->46 13 J K = KK 9 4 7 L 5 = L 4 -9 4 8 CC 9 K = K1 , LN 9 4 9 B S U N 1 = 3 S I N 1 + A ( J K ) * B ( K ) * D S I G N ( 1 . C 0 , M L 5 ) ) 9 5 0 JK=JK+MM '•95 1 L 5 = L 5 + M 9 52 9 C C N T INUE 9 5 3 8 B ( L ) = ( B ( L ) - B S U K l ) * C A B S ( A { J ) ) 9 5 4 B (N ) = B (N ) * A ( Nu 1 ) 9 5 5 N N N = N K i 99 9 5 6 N N = N - 1 9 5 7 NC=N 100 9 5 0 DC 10 L = l , N K • 9.5 9.' B SDN 2 = 0 . 9 6 0 N L = N - L 961 N L l = N - L + l 9.6 2 NNN=NMM-M 9 6 3 N J I = N M M 9 6 4 I F ( L . G E . N ) N 0 = N C - 1 •65 DU 11 K = N L 1 , N C - .9&6 NJ 1=NJ1+1 ' 9,7,6 7 B S L M 2 = B S L N 2+A ( N J 1 ) * B ( K ) ' . P , : 6 8 11 C C N T I N U E ' 9 6 9 10 B ( N L ) = ( B ( N L J - E S U N 2 ) * A ( N M M ) 9 7 C R E T U R N 97 1 I CC DO 101 1 = 1 t N •97 2 I F ( A ( I ) . F Q . C . C ) CO TO 60 9-,73 101 E ( I ) = B ( I ) /A ( I ) 9:74 R E T U R N 9 7 5 END 9 7 6 S U B R O U T I N E A V M A T C ( V , R , X X , S N ) 9 7 7 I M P L I C I T R E A L * 8 ( A - H . C - 2 ) 9'78 CCMMON / B L K 2/ C4 , C 5 , S4 , S 5 , C 1 0 , C 11 , S i C , SI 1 , SMPR E V , A 9,79 C I MENS ION V I ( 3 , 3 ) , W J ( 3 , 3 ) , A A { 3 ) , T ( 5) ,P ( 5 ) , Q (6 ) , DW I { 3 ,3 ,3 ) , DW J ( 3 ,3 , <TE Ci * 3 ) , C T ( 5 , 6 ) , B P ( 5 , 6 ) , C C ( 6 , 6 ) , D D W I ( 3 , 2 , 2 , 2 ) , D D f c J ( 3 , 3 , 3 , 3 ) 9 8 1 * , A ( 6 , 1 2 ) , v ( 1 2 f 6 , 1 2 ) , X X ( 3 ) f S N ( 6 ) * R ( 1 2 ) , S M P R E V ( 6 , 6 ) 9 8 2 DO 1 J = l , 1 2 •S8 3 DC 1 1 = 1 ,6 ^ 8 4 DC 1 K = l , 1 2 ,9.8 5 A ( I , J ) = C . C C 9\86 V ( K , I , J ) = O . C C 9 E 7 1 C C N T I N U E 9 8 8 6 C 4 = U C 0 S ( R ( 4 ) ) * S N ( 2 ) - D S I N ( R ( 4 ) ) * S N ( 1 ) 9 89 S4 = C C 0 S ( R ( 4 ) ) * SN( 1 ) + D S I N ( R ( 4 ) ) * S N ( 2 ) -9 9C C 5 = CCO S ( R ( 5 ) ) * S N (4 ) -DS IN(R (5 ) ) * S N ( 2) $ 9 1 S5 = C C O S ( R ( 5 ) ) * S N ( 3 ) + D S I N ( R ( 5 ) ) * S N ( 4 ) " 9 9 2 C6 = D C 0 S ( R ( 6 > ) * S M 6 J - D S I N ( R ( 6 ) ) * S N ( 5 ) •1,9 3 S6 = C C 0 S ( R ( 6 ) ) * SN ( 5 ) + D S IN (R ( 6 ) ) * S N ( 6 ) t9 9 4 C11 = D C 0 S ( R (11 ) ) * S N ( 4 ) - C S 1 N ( R ( 1 1 ) ) * S N ( 3 ) 99 5 S I 1 = DC0 S (R ( 1 1 > ) * S N ( 2 )+D S I N (fl ( 11) ) * SN ( 4 ) i . 9 9 6 C 1 0 = D C O S ( R (1C ) ) * S N ( 2 ) -DS IN (R ( 1C ) ) * S N ( 1 ) % 9 7 S 1 C = D C O S ( R ( 1 C ) ) * S N ( 1 ) + C S I N ( R ( 1 0 ) ) * S N ( 2 ) 9 9 8 S 1 2 = D C 0 S ( R ( 1 2 ) ) * SN ( 5 ) + C S IN ( P, ( 12 ) ) * S N ( 6 ) - § 9 9 C12 = D C O S ( R ( 1 2 ) ) * S N ( 6 ) - D S I N ( R ( 1 2 ) ) * S N ( 5 ) U C C Ml ( 1 ,1 ) = C 5 * C 6 ICC 1 WI ( 1 , 2 ) = - C 5 * S 6 1 C C 2 W I ( 1 , 3 ) = S 5 U C 3 ' w l ( 2 , l ) = C 4 * S 6 + S 4 * S 5 * C 6 ~T.C04 W I ( 2 , 2 ) = C 4 * C 6 - S 4 * S 5 * S t 1-CC5 W I ( 2 , 3 ) = - S 4 * C 5 1CC6 W I ( 3 , 1 ) = S 4 * S 6 - C 4 * S 5 * C f c 1C07 W I ( 3 , 2 ) = S 4 * C 6 + C 4 * S 5 * S 6 •1C08 K l ( 2 , 3 ) = C 4 * C 5 1 C C 9 W J ( 1 , 1 ) = C 1 1 * C 1 2 1C1C W J ( 1 , 2 ) = - C l l * S 1 2 1C11 W J ( 1 , 3 ) = SI I 1 C 1 2 W J ( 2 , 1 ) = C 1 C * S 1 2 + S 1 0 * S 1 1 * C 1 2 1 C 1 3 K J ( 2 , 2 ) = C 1 C * C 1 2 - S 1 0 * S 1 1 * S 1 2 1C14 W J ( 2 , 3 ) = - S l C * C l l 1 C 1 5 W J ( 3 , 1 ) = S 1 C * S 1 2 - C 1 G * S 1 1 * C 1 ? K 16 1C 17 IC 18 .019 JC20 C21 C2 2 C23 C24 •'2 5 i. C 2 6 1C27 ICT28 1C29. l,C3Ci 10:31 K 3 2 1 02 3 1 0 4 IC35 1C36 LC27 Cr3 8| LC>39 I '•1 £4 c 1C41 IC42 l\<3 1C44 LC&4 5 1^46 1C47 1C48 1C49 K 5 C 1^5 1 1G52 H"5 3 l'c54 , 1C55 1£56 1657 1C58 1-C5 9 '1C60 1C61 v?,C6 2 1&>3 1C64 1-C 6 5 . 1C66 ! 1C67 1C6 8 IG6 9 1C70 •1C71 1C7 2 1 C73 1C74 1 C75 WJ(3,2)= 51C*C12+C10*S11*S12 WJ ( 3 , 3 > = C1C-C1 1 AA( I )=XX ( 1 ) + * M 2 ) = X X ( 2 M A A ( 3 ) = XX( 3 ) + R( 7 )-R( 1) R(8 )-R{2 ) R ( 9 ) - R ( 3 ) 101 DO 20 K = l , 3 CC 20 1=1,3 DC 20 J = l , 3 GC TO( 21 , 20* 23 ) ,K EQ . 1 ) C W I ( K , I , J ) EQ.2)CW I (K, I , J ) EQ.3JCWI ( K , I ,J) EQ . l )CWJ(K, I , J ) .EQ.2 )CWJ (K, I , J ) EC.3)DWJ(K, I ,J) 21 IF ( I IF ( I I F ( I. I F ( I IF ( I IF ( I , GCTC20 23 I F ( J . E Q . l ) G W I (K, I , J ) I F ( J . E C . 1 ) D W J ( K , I , J ) IF ( J.EQ.2 )DW I (K , I , J) I F ( J . E Q . 2 ) C W J ( K , I , J ) IF ( J.EQ.3 >CWI(K, I , J ) I F ( J . E Q . 3 ) C W J ( K , I ,J) 2 0 CCNTINUE O.DO -W I ( 3 , J ) ' W 1 ( 2 , J ) = G.DO - W J ( 3 , J ) W J ( 2 , J ) WK 1,2) = WJ ( 1,2 ) -WK 1,1) - W J ( I , 1 ) =0 .DO = 0 .DO CW 1 ( 2 , 1 , 1 ) DWI(2,1,2) DWI ( 2, 1,3) DWK 2,2, 1 ) DWI ( 2 , 2 , 2 ) CVi I ( 2 , 2 , 3 ) DW I ( 2, 3, 1 ) CWI ( 2 , 3 , 2 ) DWI(2,3,3) DWJ(2,1,1) DWJ(2,1,2 ) = 0W J ( 2 , 1 , 3 ) = DWJ(2,2,1 ) = DWJ(2,2,2) DWJ(2 ,2 ,3 ) DWJ(2,3,1) 2) = 3 )=--S5*C6 S5*S6 C5 S4*C5*C6 -S4*C5*S6 S4*S5 -C4*C5*C6 C4*C5*S6 -C4*S5 -S 11*C 12 S11*S12 C l l S 1C*C11*C12 -S 10*C 11*512 S I G * S 1 1 -C1C*C11*C12 C10*C11*S12 -C 1C*S 1 1 W J ( 2 , 3 ) * W I ( 3 , L ) + W J ( 3 , 3 ) * W I ( 2 , 1 ) WJ( 1,3 ) * S 5 *C 6 + W J ( 2 , 3 ) * S 4 * C 5 * C 6 - W J ( 3 , 3 ) * C 4 * C 5 * C 6 W J ( 1 ,3 )*H I ( 1, 2 ) + WJ ( 2 , 3 )*WI ( 2 , 2 )+WJ ( 3, 3 )*WI ( 3,2) T( 1 ) C l 1*W I ( 1, 1 ) + S 1 0 * S 1 1 * W I ( 2 , 1 ) - C 1 0 * S 1 1 * W I ( 3 , 1 ) -W J ( 2 , 3 ) *W I ( 3 , 2 ) +WJ(3,3 ) * W I ( 2 , 2 ) W J t l ,3 )*S'5*S6-W J (2,3 ) *S4 *C5 * S 6 +H J ( 3 , 3 )*C4*C5*S6 WJ( 1 , 3)*WI ( 1 , 1 ) - W J ( 2 , 3 ) * W I ( 2 , 1 ) - W J ( 3 , 3 ) * W I ( 3 , 1 ) -P ( I ) C 1 1 * W I ( 1 , 2 ) + S10*S11*V»I(2,2)-C10*S11*KI(3,2) W J ( 2 , 2 ) * W I ( 3 , 1 ) + W J ( 3 , 2 ) * W I ( 2 , 1 ) W J ( 1 , 2 ) * S 5 * C 6 + W J ( 2 , 2 ) * S 4 * C 5 * C 6 - W J ( 3 , 2 ) * C 4 * C 5 * G 6 U J ( l , 2 ) * W I ( l , 2 ) + f c J ( 2 , 2 ) * W I ( 2 , 2 ) + W J ( 3 , 2 ) * W l ( 3 , 2 ) -G ( 1 ) G(5 ) = S11*S12*WI ( 1, 1 ) - S 1 0 * C l l * S l 2 * V i I (2 ,1 ) + C1C*C11*SI 2*WI C ( 6 ) = - W J ( l , l ) * W I ( l , l ) - W J ( 2 , l ) * W I ( 2 , l ) - K J ( 3 , l ) * W I ( 3 , l ) DC l t» J = l , 3 DC 16 K= 1,2 A ( J , K ) = - W I ( K , J ) DWJ(2,3, DWJ (2 , 3, T( 1) T ( 2 ) T (3) T( 4) T (5 ) P ( 1 ) P( 2) P ( 3 ) P ( 4 ) ? (5 ) Q ( 1 ) C(2 ) 0 ( 3 ) C ( 4 ) (3,1) C76 A ( J , K + 6 ) = - M J , K ) C77 A ( J , 5 ) = A A ( K )*CWI ( 2 , K , J ) + A ( J , 5 J 102 C78 16 C C N T I N U E Cs7<5 A ( J , 4 ) = - A A ( 2 ) * W I ( 3 , J ) + A A ( 3 ) * w l ( 2 » J ) C « 0 A ( 1 , 6 ) = A A ( J ) * W l ( J , 2 ) + A ( l , 6 ) C.B1 A ( 2 , 6 ) = - £ A ( J ) * W I ( J , 1 ) + A ( 2 , 6 ) C82 15 C C M I N U E C 8 3 DC 12 K = l , 3 C84 A ( 5 , K + 3 ) = T ( K ) U35 A ( 4 , K + 3 ) = - P (K ) , Lfc'fe A ( 6 , K + 3 ) = - C (K ) 10.87 1 IF ( K . E G . 2 ) G C T C 1 2 4. £ 8 8 A ( 5 , K + 9 ) = T t K + 3 ) •C8S A ( 4 , K + 9 ) = - P ( K + 3 ) C 9 0 A ( 6 , K + 9 ) = - C ( K + 3) C*91 12 C O N T I N U E G92 A ( 6 , 1 2 ) = - C ( 6 ) :CS'3 DC 84 N = l , 3 .C94 DC 8 4 K = l , 3 C 9 5 DO 8 4 1 = 1 , 3 C96 CC 8 4 J = l , 3 C 9 7 G C T O ( 8 5 , 6 4 , 8 6 ) ,K C § 8 85 IF ( I . E Q . l )CCW K N , K , I , J ) = C . D O xC99 I F ( I . E Q . 2 ) C C W I ( N , K , I, J ) =-OW I (N , 3 , J ) U P O ; I F ( I . E C - . 3 ) C C V I ( N , K , I , J ) = C W I ( N , 2 , J ) 101 I F ( l . E Q . l ) C D W J ( N , K , I , J ) = 0 . D 0 102 I F < I . E Q . 2 ) C C W J ( N , K , I , J ) = - D W J ( N , 3 , J > l'cC3 1 I F U . E C . 3 ) C C W J ( N , K , I , J ) = D W J ( N f 2 , J ) ,104 G C T C 8 4 1105 86 I F ( J . E C . l ) C D W I ( N , K , I , J ) = C W I ( N , I , 2 ) 11X6 I F ( J . E C . l ) D D U J ( N , K , I , J ) = C f c J ( N , I , 2 ) 1 107 I F ( J . E Q . 2 ) D D U ( N , K , I , J ) = - C W I ( N , I , 1 ) 1C8 I F ( J . E Q . 2 ) C 0 W J ( N , K , I , J ) = - D W J ( N , I , 1 ) 1C9 IF ( J . E G . 3 ) C C V I ( N , K , I , J ) = 0 . C O 1 1 0 I F ( J . E Q . 3 ) C C W J ( N , K , I , J ) = 0 . CO .ia 1 84 C C N T I N U E i l 2 CC 87 J= 1 , 3 3 CCWI ( 1 , 2 , 1 , J ) = O . C C 1 M 4 CCWJ ( 1 , 2 , 1, J ) = C . C C 1 1 5 CCWI ( 3 , 2 , J , 3 ) = 0 . C 0 4,16 C C W J ( 3 ,2 , J , 3 ) = C . C 0 117 C C W I ( 3 , 2 , J , I ) = C W I ( 2 , J , 2 ) 118 C C W J ( 3 , 2 , J , 1 ) = D W J ( 2 , J , 2 ) U 1 1 9 CCW I (3 , 2 , J , 2 ) = -DW I ( 2 , J , 1 ) hAZO CCWJ( 3 , 2 , J , 2 ) = - C W J ( 2 , J , l ) 1121 C C W I ( 2 , 2 , 1 , J )= -WI ( 1 , J ) {-122 CCWJ (2 ,2 , 1 , J ) = - W J ( 1, J ) 1123 CCWI ( 1 ,2 , 2 , J ) = - C W I ( 2 ,3 , J ) C C W J ( 1 , 2 , 2 , J ) = - C W J ( 2 , 3 , J ) C C W I ( 1 , 2 , 3 , J )= C W I ( 2 , 2 , J ) C C W J ( 1 , 2 , 3 , J ) = C W J ( 2 , 2 , J ) 87 C C N T I N U E CCWI ( 2 , 2 , 3 , 1 )= C 4 * S 5 * C 6 C C W J ( 2 , 2 , 3 , 1 )= C 1 C * S 1 1 * C 1 2 CCWI (2 ,2 , 3 , 2 ) = - C 4 * S 5 * S 6 CCWJ ( 2 , 2 , 3 , 2 ) = - C 1 0 * S l l * S 1 2 CCWI ( 2 , 2 , 3 , 3 ) = - C 4 * C 5 C C W J { 2 , 2 , 3 , 3 ) = - C l C * C l l CCWI ( 2 , 2 , 2 , 1 ) = - S 4 * S 5 * C 6 C C W J ( 2 , 2 , 2 , 1 ) = - S 1 0 * S 1 1 * C 1 2 1136 C C W I ( 2 , 2 , 2 , 2 )= S 4 * S 5 * S 6 1137: C C W J ( 2 , 2 , 2 , 2 ) = S 1 - C * S U * S 1 2 103 11 381 CCW I ( 2 , 2 , 2 , 3 - ) = S 4 * C 5 V 3 9] C C W J ( 2 , 2 , 2 , 3 ) = C 1 1 * S 1 0 •140| V (5 , 1, 1 ) = S 5 * C 6 1 4 1! V ( 6 , 1 ,1 ) =- W I ( 1 , 2 ) i;42 V ( 4 , 1 , 2 ) = W I ( 3 , 1 ) 143! V ( 5 , 1 , 2 ) = - S 4 * C 5 * C 6 144, V ( 6 , 1 , 2 ) = - W I ( 2 , 2 ) 4 5 | V ( 4 , 1 ,3 )=-W I ( 2 , 1 ) ^ 4 ' 6 , V ( 5 , l , 3 ) = C 4 * C 5 * C 6 I 147 1 V (6 , 1 , 3 ) = - W I ( 3 , 2 ) 1148 V ( 4 , 1 , 4 ) = - A A ( 2 ) * W I ( 2 » 1 ) - A A ( 3 ) * W I ( 3 , 1 ) I'14 9 I V ( 5 , l , 4 ) = A A ( 2 ) * C 4 * C 5* C 6 4 £ A ( 3 ) * S 4 * C 5 * C 6 1150 , V (6 , 1 , 4 ) = - A A ( 2) * W l ( 3 ,2 ) + /5 A ( 3 ) *W I ( 2 ,2 ) 11-5 1 V ( 8 , 1 , 4 ) = - W I ( 3 , 1 ) !il'52 V < 9 , 1 , 4 ) = W I ( 2 , 1 ) 115 3 V ( 5 , 1 , 5 ) = - £ A ( 1 ) * C 5 * C 6 - * A { 2 ) * S 4 * S 5 * C 6 + A A ( 3 ) * C 4 * S 5 * C 6 I £ 5 4 V ( 6 , l , 5 ) = A M 1 ) * S 5 * S 6 - A A ( 2 ) * S 4 * C 5 * S 6 + A A ( 3 ) * C 4 * C 5 * S 6 155 V ( 7 , l , 5 ) = - S 5 * C 6 .156 V ( 8 , l , 5 ) = S 4 * C 5 * C 6 137 V ( 9 , 1 , 5 ) = - C 4 * C 5 * C 6 .158 V ( 6 , 1 , 6 ) = - £ A ( 1 )*WI ( 1 , 1 ) - A A ( 2 ) * W I ( 2 , 1 ) — A A ( 3 ) * W I ( 3 , 1 ) 11,59 V ( 7 , l , 6 ) = W I ( 1 , 2 ) U 6 0 '• V ( 8 , 1 , 6 ) = W I ( 2 , 2 ) 116 1 V ( 9 ,1 , 6 ) = W I ( 3 , 2 ) 162 V ( 5 ,2 , 1 ) = - S 5 * S 6 L6 3 V (6 , 2 , 1 ) = W I ( 1 , 1 ) . l r 64 V ( 4 , 2 , 2 ) = W 1 ( 3 , 2 ) I 1^65 V ( 5 , 2 ,2 )= S 4 * C 5 * S 6 l l r £ 6 V( 6 , 2 , 2 ) = WI ( 2 , 1 ) .167 V ( 4 , 2 , 3 ) = - W I ( 2 , 2 ) 168 V ( 5 , 2 , 3 ) = - C 4 * C 5 * S 6 169 V ( 6 , 2 , 3 ) = W I ( 3 , 1 ) 170 - V ( 4 , 2 , 4 ) = - A A ( 2 ) * W I ( 2 , 2 ) - A A ( 3 ) * W I ( 3 , 2 ) 1<71 V ( 5 , 2 , 4 ) = - « A ( 2 ) * C 4 * C 5 * S 6 - A A ( 3 ) *S 4 * C 5* S t t '172 V I 6 , 2 , 4) = A A ( 2 ) * W 1 ( 3 , 1 ) - * M 3 ) * W 1 ( 2 , 1 ) U<»7 3 V ( 8 , 2 , 4 ) = - W I ( 3 , 2 ) 4>74 V ( 9 , 2 , 4 ) = W I ( 2 , 2 ) 1175 V ( 5 , 2 , 5 ) = 2 A ( 1 ) * C 5 * S 6 + A A (2 ) * S 4 * S 5 * S 6 - A A ( 3 ) * C 4 * S 5 * S 6 . j l 76 V ( 6 , 2 , 5 ) = A A ( 1 ) * S 5 * C 6 - A A( 2 ) * S 4 * C 5 * C 6 + .AA ( 3 ) * C 4 * C 5 * C 6 177 V ( 7 , 2 , 5 ) = S 5 * S 6 178 V ( 8 , 2 , 5 ) = - S 4 * C 5 * S 6 1 ?9 V ( 9 , 2 , 5 ) = C 4 * C 5 * S 6 „ l'SO V ( 6 , 2 , 6 ) = - A A { 1) *W1 ( 1,2 ) -AA{ 2) *Wl ( 2 , 2 ) - A A ( 3 )*WI ( 3 ,2 ) 1181 V ( 7 , 2 , 6 ) = - W I ( 1 , 1 ) U 8 2 V ( 8 , 2 , 6 ) = - W l ( 2 , 1 ) 116 3 V I 9 , 2 , 6 ) = - W I ( 3 , I ) -.-184 V ( 5 , 3 , 1 )= - C 5 I I € 5 . V ( 4 , 3 ,2 ) = W I ( 3 , 3 ) 11-8 6 V ( 5 , 3 , 2 ) = - S 4 * S 5 1187 V ( 4 , 3 , 3 ) = -W I ( 2 , 3 ) '.188 V ( 5 , 3 , 3 )= C 4 * S 5 U « 9 V ( 4 , 3 , 4 ) = - A / ( 2 ) * W I ( 2 , 3 ) + A A ( 3 ) * C w I ( l , 2 , 3 ) 190 V ( 8 , 3 , 4 ) = - W I ( 3 , 3 ) 1:! :> 1 V ( 9 , 3 , 4 ) = W I { 2 , 3 ) 1192 V ( 5 , 3 , 5 ) = - A A ( 1 ) * S5 + A A ( 2 ) * S 4 * C 5 - A A ( ? ) * C 4 * C 5 1193 V ( 7 , 3 , 5 ) = C5 .194 V ( 8 , 3 , 5 ) = S 4 * S 5 .19 5 V ( 9 , 3 , 5 ) = - C 4 * S 5 I ] 196 PC 120 J = l , 3 1 1 9 7 D T ( I , J ) = - W J ( 2 , 3 ) *DW I ( J , 3 , l ) + k J ( 3 , 3 ) * 0 K l ( J , 2 , l ) 104 1 1 9 8 DT (1 , J+3 ) = - C W J ( J , 2, 3 )*V* I ( 3, I ) - C W J ( J , 2, 3 )*KI ( 2 , 1) 1 9 9 D T ( 2 , J ) = + V « J ( 1 , 3 ) * C C W I ( J , 2 , 1 , 1 ) + W J ( 2 , 3 )*DDW I ( J , 2 , 2 , 1 ) + W J ( 3 , 3 ) * C C K 1 ( .iCO * J , 2 , 3 , I ) 1 ? 0 1 D T ( 2 , J + 3 ) = + C W J ( J , l , 3 ) *DW I ( 2 , 1 , 1 ) + D W J ( J , 2 , 3 ) * C W I ( 2 , 2 , 1 ) +DW J ( J , 3 , 3) * 12X2 * C M ( 2 , 3 , 1 ) L 2 C 3 CT (3 , J ) = V t J ( I , 3 ) * D W l ( J , l , 2 ) + k J ( 2 , 3 ) * C U I ( J v 2 , 2 ) + W J ( 3 , 3 ) ± D W l ( J , 3 , 2 ) 12C4 D T ( 3 , J + 3 ) = D W J ( J , 1 , 3 ) * W I ( 1 , 2 ) + C W J ( J , 2 , 3 ) * W I ( 2 , 2 ) + D W J ( J , 3 , 3 ) * W I ( 3 , 2 ) - H 5 D T ( 4 , J ) = - D T ( 1 , J ) 1 • ?C6, CT (4 , J + 3 ) = - C T ( 1, J + 3) 1 ?C7i D T ( 5 , J ) = DV. J ( 2 , 1 , 3) *DW I ( J , 1, 1 ) + DWJ ( 2 , 2 , 3 ) * 0 M ( J , 2 , 1 ) + D W J ( 2 , 3 , 3 ) * .? C 8, * C W I ( J , 3 , 1 ) 1 2 C 9 D T ( 5 , J + 3 ) = C D f c J ( J , 2 , l , 3 ) * w I ( l , l ) + C 0 V i J ( J , 2 , 2 , 3 ) * W I ( 2 , l ) + 0 0 W J ( J , 2 , 3 , 3 1 2 1 0 * ) * W I ( 3 , 1 ) t ? l l C F ( 1 , J ) = - k J ( 2 , 3 )*DW I ( J , 3 , 2 ) + W J ( 3 , 3 ) * D W I ( J , 2 , 2 ) 1 i 12 C P ( l , J + 3 ) = - D V J ( J , 2 , 3 ) * V » I ( 3 , 2 ) + C W J ( J , 3 , 3 ) * W I ( 2 , 2 ) 12, 13 DP ( 2 , J ) = W J ( 1 , 3 ) * D D W I ( J , 2 , 1 ,2 ) + V>J( 2 , 3) * D D WI ( J , 2 , 2 , 2 ) + W J (3 , 3 ) * CCH i ( 121-*, * J , 2 , 3 , 2 ) 12 15, C P ( 2 , J + 3 ) = D f c J ( J , 1 , 3 ) * C K I (2 , 1 , 2 ) + n w j ( J , 2 , 3 ) * C M ( 2 , 2 , 2 ) + D W J ( J , 3 , 3 ) * 121 ; *DK 1 ( 2 , 2 , 2 ) I 2 1 V D P ( 3 , J ) = - W J ( 1 , 3 ) * D W I ( J , 1 , 1 ) - W J ( 2 , 3 ) * D W I ( J , 2 , 1 ) - W J ( 3 , 3 ) * D W I ( J , 3 , 1 ) 12 K C P ( 3 , J + 3 ) = - B W J ( J , l , 3 ) * M ( l , l ) - C W J ( J , 2 , 3 ) * f c I ( 2 , l ) - D W J ( J , 3 , 3 ) * W I ( 3 , 1 ' 1 2 1 9 * ) |l?20'i D P ( 4 , J ) = - D P ( 1 , J ) 1 2 2 1 D P ( 4 , J + 3 ) = - C P ( l . J + 3 ) 12 2 2 D P ( 5 , J ) = D W J ( 2 , l , 2 ) * D W I { J , l , 2 ) + C K J ( 2 , 2 t 3 ) * D W I ( J , 2 , 2 ) + D W J ( 2 , 3 , 3 ) * H 2 3! * D W I ( J , 3 , 2 ) 4Z2-'.. C P ( 5 , J + 3 ) = C C k J ( J , 2 , l , 3 ) * W I ( 1 , 2 ) + C 0 W J ( J , 2 , 2 , 3 )*Vi I ( 2 , 2 ) + D D W J ( J , 2 , 3 , 3 12 2 5 * ) * W I (2 ,2 ) 1 2 / DC ( 1 , J ) =- Vs J ( 2 , 2 ) *DVs I ( J , 3 , 1 ) + W J ( 3 , 2 ) * D W I ( J , 2 , 1 ) 1 2 2 ; C C ( l , J + 3 ) = - C k J ( J , 2 , 2 ) * W I ( 3 , l ) + C W J ( J , 2 , 2 ) * K 1 ( 2 , 1 ) 12 26 DC( 2 , J ) =V* J (1 , 2 ) *CCWI ( J ,2 ,1 ,1 ) + V v J ( 2 , 2 ) # B D W l ( J » 2 , 2 , l ) +W J ( 3t 2 )*CCW I ( J 1 2 2 9 * , 2 , 3 , 1 ) 1-2 3G CC (2 , J+3 ) = CWJ ( J , 1 , 2 ) 4 0 W K 2 , 1, 1) + DW J ( J , 2, 2 ) *D V»I ( 2 , 2 ,1 ) +D W J ( J , 3 , 2 ) * 1 2 3 1 * D U ( 2 , 3 , 1 ) 12 3 2 D C ( 3 , J ) = W J { 1 , 2 ) * C W I { J , 1 , 2 ) ^ W J ( 2 , 2 ) * U W I ( J , 2 , 2 ) + K J ( 3 , 2 ) * D K I ( J , 3 , 2 ) I-C? 3 3 D C ( 3 , J + 3 ) = C t t J ( J , l , 2 ) * W I ( l , 2 ) + C W J ( J , 2 , 2 ) * W I ( 2 , 2 ) + D W J ( J , 3 , 2 ) * W I ( 3 , 2 ) 1-2 3 4, D C ( 4 , J ) = - D Q ( 1 , J ) '12 3 5, DC (4 , J + 3 ) = - 0 C ( 1 , J + 3 ) 1A,2 36 CQ(5,J)=DWJ(2,1,2)*DW I ( J , 1 , 1 ) + C W J ( 2 , 2 , 2 ) *DWI ( J , 2 , 1 ) + D W J ( 2 , 3 , 2 ) * 12,37 * D k I ( J , 3 , 1 ) 12 38| D C ( 5 , J + 3 ) = C D W J ( J , 2 ,1 , 2 ) * W I ( 1 , 1 ) + C D W J ( J , 2 , 2 , 2 ) * W I ( 2 , 1 ) + D D K J ( J , 2 , 3 , 2 1 1 " 3 9 ; * ) * W I ( 3 , 1 ) i 2 4 0 ; CC ( 6 , J ) = - V « J ( I , 1 )*DWI ( J , 1, 1 ) -WJ ( 2 , 1 ) * D W I ( J , 2 , 1 ) - W J ( 3 , 1 ) * D W I ( J , 3 , 1 ) 1 2 4 1 C C ( 6 , J + 3 ) = - D f c J ( J , 1 , 1 ) * f c I ( 1 , 1 ) - C W J ( J , 2 , 1 ) * W I ( 2 , 1 ) - D W J ( J , 3 , 1 ) * W I ( 3 , 1 ' 1 2 4 2 ! * ) ' 1 2 4 3 i 1 2 C C C N T I N U E 1 2 4 4 V ( 4 , 5 , 4 ) = D T ( 1 , 1 ) 1-2 4 5 J \ M 5 , 5 , 4 ) = 0 T ( 1 , 2 ) 124 6 , V ( 6 , 5 , 4 ) = C T ( 1 , 3 ) 1 2 4 7 V ( 1 C , 5 , 4 ) = C T ( 1,4 ) 12 4 8 V( 1 1 , 5 , 4 ) = C T ( 1 , 5 ) 1 2 4 9 , V ( 5 , 5 , 5 ) = C T ( 2 , 2 ) 1 2 5 0 , V ( 6 , 5 , 5 ) = C T ( 2 , 3 ) ' I 2 5 1 V( 1 C , 5 , 5 ) = C T ( 2 , 4 ) 12 52 V ( 1 1 , 5 , 5 ) = 0 T ( 2 , 5 ) 1 2 5 3 V ( 6 , 5 , 6 ) = C T ( 3 , 3 ) 1 2 5 4 V ( 1 C , 5 , 6 ) = D T < 3 , 4 ) 1 2 5 5 V ( 1 1 , 5 , 6 ) = C T ( 3 , 5 ) 12 5 6 V ( 1 0 , 5 , 1 0 ) = C T ( 4 , 4 ) 1 2 57 V( 1 1 , 5 , 1 C ) = C T ( 4 , 5 ) l n c 1258 V ( U , 5 , 1 1 ) = C T ( 5 , 5 ) 2 5 9 DO 31 J = l , 3 126C DC 32 K = l , 3 1--2 61 V ( K + 3 , 4 , J + 3 ) = - C P ( J , K ) v 1'262 I F ( K . G E . 3 >GOTO 32 1 2 6 3 V (K + 9 , 4 , J+3 ) = - C P ( J , K - » 3 ) 1 2 6 4 I F ( J . G E . 3 1 G C T C 3 2 ' 2 6 5 V ( K + 9 , 4 , J + 9 ) = - D P ( J + 3 , K + 3 ) i ' 266 32 C C N T I N U E ' 2 6 7 31 C C N T I N U E . 2 6 6 DC 41 J = l , 3 1 2 6 9 DO 4 2 K = 1 , 3 1 2 7 C V ( K + 3 , 6 , J * 3 ) = - C C ( J , K ) \ \27 1 V( K + 9 ,6 , J 4 3 ) =-DC { J ,K + 3 ) 1 2 7 2 V ( K + 9 , 6 , J + 9 ) = - D G ( J + 3 . K + 3 ) 1*27 3 42 C C N T I N U E l t274 41 C C N T I N U E 1 2 7 5 DC 51 K = l , 1 2 1 2 7 6 DO 52 J = l , 6 1V277 DC 53 1 = 1 , 1 2 12 7 8 V ( I , J , K ) = V ( K , J , I ) ) ; 2 7 9 53 C O N T I N U E );280 , 52 C C N T I N U E 1 2 8 1 51 C C N T I N U E 1 2 8 2 R E T U R N V 2 8 3 ENC 1*284 S U B R O U T I N E AV M ATS (V , R , XX , SN ) j>285 I M P L I C I T R E A L * 8 ( A - H , C - 2 ) it? 86 CCMMON /B L K 2/ C4 , C 5 , S4 , S 5 , C 10 , C 1 L , SI C , SI 1 , SM P R E V , A 1 2 8 7 D I M E N S I O N A ( 6 , 1 2 ) , V ( 1 2 , 6 , 1 2 ) , S M P R E V ( 6 , 6 ) , X X ( 3 > , S N ( 6 ) , R ( 1 2 ) 1 2 8 8 „ D I M E N S I O N W I I ( 3 , 3 ) , W J J { 3 , 3 ), A * ( 3 ) , > S I < 3 ) , BS ( 3 ) , Y S J ( 3 ) , Y SK ( 3 ) , 1 2 8 9 * Z S I ( 3 ) , C S ( 3 ) , Z S J ( 3 ) , Z S K { 3 ) , D W I ( 3 , 3 , 3 ) , D W J ( 3 , 3 , 3 ) , D E L ( 9 ) , D Y S 1 ( 3 , 9 ) , 1 2 9 C * D B S ( 3 , 9 ) , D C S ( 3 , 9 ) , D Y E S J ( 9 ) , D Z E S J ( 9 ) , D Y S J ( 3 , 9 ) , D Z S J ( 3 , 9 ) , D Y S K ( 3 , 9 ) , It25 1 * D Z S K ( 3 , 9 ) , C Y IS I (9 ) , C Z IS I ( 9 ), D Y I S J ( 9 ) , D Z I S J (9 ) , C Y I SK( 9 ) ,DZ I SK( 9) 129 2 D I M E N S I O N C D E L ( 9 , 9 ) , C D Y S I ( 9 , 3 , 9 ) , D C W I ( 3 , 3 , 3 , 3 ) , C D W J ( 3 , 3 , 3 , 3 ) , v . '293 * C C B S ( 9 , 3 , 9 ) , C D C S ( 9 , 3 , 9 ) , D D Y E S J ( 9 , 9 ) ,DD ZE S J ( 9 , 9) , DD YS J ( 9 , 3 , 9 ) , 12 94 * C C Z S J ( 9 , 3 , 9 ) , C C Y S K ( 9 , 3 , 9 ) , D D Z S K ( 9 , 3 , 9 ) , D D Y I S I ( 9 , 9 ) ,DDZ I S I ( 9 , 9 ) , 1 2 9 5 * C C Y I S J ( 9 , 9 ) , C C Z I S J ( 9 , 9 ) , D D Y I S K ( 9 , 9 ) , C C Z I S K ( 9 , 9 ) ,1296 DO 1 J = l , 1 2 1:297 DC 1 1 = 1 , 6 1 2 9 8 DC 1 K=l,12 -2299 A ( I , J ) = C . D C I^3CC V ( K , I , J ) = O . C C 1 3 C 1 1 C C N T I N U E 1 3 0 2 6 C 4 = D C C S ( P ( 4 ) ) * S N ( 2 ) - D S I N ( R ( 4 ) ) * S N ( 1 ) 1*3 0 3 S4 = D C 0 S ( R ( 4 ) ) * SN( 1 H D S I N ( R ( 4 ) ) * S N ( 2) 1304 C5 = CCO S ( R ( 5 ) ) * S N ( 4 ) - D S I N ( R ( 5 ) ) * S N ( 2 ) C 5 S5 = C C G S ( R ( 5 ) ) * S M 3 ) + C S I N ( R ( 5 ) ) * SN<4 ) i 3 C 6 C6= DCO S ( R ( 6 ) ) * S N ( 6 ) - D S I N ( R ( 6 ) ) * S N ( 5) 13 07 S6 = C C 0 S ( R ( 6 ) ) * S K ( 5 ) + D S I N ( R ( 6 ) ) * S N ( 6 ) <13C8 C 1 1 = D C C S ( R ( 1 1 ) ) * S N ( 4 ) - C S I N ( R ( 1 1 ) ) * S N ( 3 ) i 3 C 9 S I 1 = DCO S ( R ( 1 1 ) ) * S N ( 2 ) + C S I N ( R ( I 1 ) ) * S N ( 4 ) 1 3 1 0 C 1 C = D C 0 S ( R ( 1 C ) ) * S N ( 2 ) - C S I N ( R ( 1 C ) ) * S N ( 1 ) 1 3 1 1 S1C = D C 0 S ( R ( 1 C ) ) * S N ( 1 ) + C S IN ( R ( 1 0 ) ) * S N ( 2 ) 1 3 1 2 S 1 2 =DC 0 S ( R ( 1 2 ) ) * S N ( 5 ) +D S I N ( R ( 1 2 ) )* S N ( 6 ) 1 3 1 3 C 1 2 = D C 0 S ( R ( 1 2 ) ) * S N { 6 ) - D S I N ( R ( 1 2 ) ) * S N ( 5 ) 1 3 1 4 W I I ( 1 ,1 ) = C 5 * C 6 1 3 1 5 W I I ( 1 , 2 ) = - C 5 * S 6 I 1 3 1 6 W I I ( 1 ,3 ) = S 5 1 3 1 7 W I I ( 2 , 1 ) = C 4 * S 6 + S 4 * S 5 * C 6 106 1 3 1 8 W I I ( 2 , 2 ) = C 4 * C 6 - S 4 * S 5 * S 6 ' 3 19 W I I ( 2 , 3 ) = - S 4 * C 5 1320 WI I ( 3 , 1 ) = S 4 * S 6 - C 4 * S 5 * C 6 h?.2l W l l ( 3 , 2 ) = S 4 * C 6 + C 4 * S 5 * S 6 , 2 22 W I I ( 3 , 2 ) = C 4 * C 5 1 3 2 3 W J J ( 1 , 1 ) = C 1 1 * C 1 2 1 2 2 4 WJJ ( 1 , 2 ) = - C 1 1 * S 1 2 ' 2 2 5 W J J ( 1 , 3 ) = S l l 1 3 2 6 W J J ( 2 , 1 ) = C 1 C * S 1 2 + S 1 G * S U * C 1 2 ' 2 2 7 W J J ( 2 , 2 ) = C 1 C * C 1 2 - S 1 0 * S U * S 1 2 2 2 8 W J J ( 2 , 3 ) = - S 1 C * C 1 1 1 3 2 9 W J J ( 3 , 1 ) = S I C * S 1 2 - C 1 0 * S 1 1 * C 1 2 1 2 3 C W J J ( 3 , 2 ) = S 1 C * C 1 2 + C 1 0 * S 1 1 * S 1 2 \ 3 2 l W J J ( 3 , 3 ) = C 1 C * C 1 1 1*3 3 2 A A ( 1 ) = X X < 1 ) + R ( 7 ) - R ( l ) ^ 3 2 3 A A ( 2 ) = X X ( 2 M R ( 8 ) - R ( 2 ) f 3 3 4 A A ( 3 ) = X X ( 3 ) + R ( 9 ) - R ( 3 ) 13 3 5 EL =D SC R T ( A A ( 1 ) * * 2 + A A ( 2 ) * * 2 + AA ( 3 ) * * 2 ) 1 3 3 6 Y S I (1 ) = A A ( 1 J / E L 1 2 3 7 Y S I ( 2 ) = A A ( 2 ) / E L 1 2 3 8 - Y S I ( 3 ) = / > A ( 3 ) / E L V 3 3 9 B S ( 1 ) = WI 1 ( 2 , 3 ) * Y S I ( 3 ) - V. I I { 3., 3 ) * Y S I ( 2) ; 3 4 b B S ( 2 ) = W I I ( 3 , 3 ) *YS I (1 )-WI I ( 1 , 3 ) * Y S I ( 2 ) 1 3 4 1 B S ( 3 ) = WII ( 1 , 3 ) * Y S I ( 2 ) - W I I ( 2 , 3 ) * Y S I ( 1 ) 1 3 4 2 Y E S J = D S Q R T ( B S ( 1 ) * * 2 + B S ( 2 >**2 + B S ( 3 ) * * 2 ) 11343 Y S J ( 1 ) = E S ( 1 ) / Y E S J l!34A Y S J ( 2 ) = E S ( 2 ) / Y E S J ^ 3 4 5 Y S J ( 3 ) = B S ( 3 ) / Y E S J £ 3 4 6 C A L L C R O S S ( Y S I , Y S J , Y S K ) 1 2 4 7 Z S I ( I ) = Y S I U ) 1 2 4 8 - Z S I ( 2 ) = Y S I ( 2 ) 1 3 4 9 ZS I ( 3 ) = Y S I ( 3 ) <125C C S ( 1 ) = W J J ( 2 , 3 ) * Y S I ( 2 ) - W J J ( 3 , 3 ) * Y S I ( 2 ) 1 2 5 1 C S ( 2 ) = W v J ( 3 , 3 ) * Y S I ( l ) - W J J ( l , 3 ) * Y S I ( 3 ) 1 3 5 2 C S ( 3 ) = W J J ( 1 , 3 ) *YS I( 2 ) - W J J ( 2 , 3 ) * Y S I ( 1 ) 13 53 Z E S J = C S G R T ( C S ( 1 ) * * 2 + C S ( 2 ) * * 2 + C S ( 3 ) * * 2 ) <125f Z S J ( 1 ) = C S ( 1 ) / Z E S J 1 3 5 5 Z S J ( 2 ) = C S ( 2 ) / Z E S J ,13 5 6 Z S J ( 3 ) = C S (3 ) / Z E S J r,2 57 C A L L C R O S S (ZS I , Z S J , Z S K ) 1 3 5 8 Y I S I = W I I ( l , l ) * Y S I ( l ) + W I I ( 2 , l ) * Y S I ( 2 ) + W I I ( 3 f l ) * Y S I ( 3 ) r 3 5 9 Z I S I = W J J ( 1, 1 ) * Y S I ( 1 ) + W J J ( 2 , 1 ) * Y S I ( 2 ) + W J J ( 3 , l ) * Y S I ( 3 ) 1 3 6 C Y I S K = W I I ( 1 , 1 ) * Y S K ( 1 ) + W I I (2 , i ) * Y S K ( 2 )+W I I ( 3 , 1 ) * Y S K ( 3 ) 136 1 Y I S J = W 1 1 ( 1 , 1 ) * Y S J ( 1)+WI I ( 2 , 1 ) * Y S J ( 2 ) + W I I ( 3 , 1 ) * Y S J ( 3 ) 136 2 Z I S J = W J J ( 1 , 1 ) * Z S J ( 1 ) + W J J ( 2 , 1 ) * Z S J ( 2 ) + W J J ( 3 , 1 ) * Z S J ( 3 ) 113 6 3 Z I S K = W J J ( 1 ,1 ) * Z S K { 1 ) + W J J ( 2 , 1 ) * Z S K ( 2 ) + W J J ( 3 , 1 ) * Z S K ( 3 ) 1 2 6 4 DC 20 K = l , 3 - 1 3 6 5 DO 20 1 = 1,3 1 3 6 6 DC 20 J = l , 3 1 2 6 7 GC T C ( 2 1 , 2 0 , 2 3 ) ,K , 1 2 6 8 21 I F ( I . E Q . 1 ) C W I ( K , I , J ) = 0 . 0 0 1;369 I F ( I . E Q . 2 ) C W I ( K , I , J ) = - W I I ( 3 , J ) 1 3 7 C IF ( I . EQ .3 )CW I ( K , I , J ) = W I I ( 2 , J ) 1 3 7 1 I F ( I . E Q . l ) D W J ( K , I , J ) = 0 . D 0 1 3 7 2 I F ( I . E Q . 2 ) C W J ( K , I , J ) = - W J J ( 3 , J ) 1 3 7 3 IF ( I . E Q . 3 ) C W J ( K , I , J )= W J J ( 2 , J ) 12 74 G C T C 2 0 1 3 7 5 23 I F ( J . E Q . l )CW I ( K , I , J)=W I I ( I , 2 ) 13 76 I F ( J . E Q . 1 ) C W J ( K , I . J ) = W J J ( I , 2 ) 1377 ' I F ( J . E C . 2 )DW I ( K , I , J ) = -W I I ( I, 1 ) 107 378 ! I F ( J . E Q . 2 ) C W J ( K , I , J ) = - W J J ( I , i ) 7 9 ! I F ( J . E Q . 3 ) C W I ( K , I , J ) = O . D 0 ; t 2 8 C I F ( J . E C . 3 ) D W J ( K , I , J ) = 0 .DO 12)81 2 0 C C N T I N U E 13*8 2' DW I ( 2 , 1,1 )= -S 5 * C 6 1383 D W I ( 2 , 1 , 2 ) = S 5 * S 6 1384 C W K 2 , 1 , 3 ) = C5 ' 8 5 DWI ( 2 , 2 , 1 ) = S 4 * C 5 * C 6 iBb CWI ( 2 , 2 , 2 ) = - S 4 * C 5 * S 6 "•£7 D W I ( 2 , 2 , 3 ) = S 4 * S 5 . 8 8 D W I ( 2 , 3 , 1 ) = - C 4 * C 5 * C 6 3 8 9 D W I ( 2 , 3 , 2 ) = C 4 * C 5 * S 6 3 9 C CWI (2 , 3 , 3 ) = - C 4 * S 5 2^91 D W J ( 2 , 1 , 1 ) = - S 1 1 * C 1 2 3 9 2 U W J ( 2 , 1 , 2 ) = S 1 1 * S 1 2 343 3 DWJ ( 2 , 1 , 3 ) = C l l 3-94 D W J ( 2 , 2 , 1 ) = S 1 C * C 1 1 * C 1 2 3 9 5 D W J ( 2 , 2 , 2 ) = - S l C * C l l * S 1 2 3 9 6 D W ' J ( 2 , 2 , 3 ) = S 1 C * S 1 1 3\97 DWJ( 2 , 3 , 1 ) = - C 1 C * C U * C 1 2 3 9 8 D W J ( 2 , 3 , 2 ) = C 1 C * C 1 1 * S 1 2 3^99; D W J ( 2 , 3 , 3 ) = - C 1 C * S 1 1 £ C C CC 30 J = l , 3 4 C 1 D E L ( J ) = - Y S I ( J ) 4 C 2 D E L ( J + 3 ) = 0 . C 4\C3 D E L ( J + 6 ) = - C E L ( J ) 4'C4 DC 31 K = l , 3 4>C5: DEL (K ) = - Y S I (K ) 4"C6 DEL (K + 6 ) = - C E L (K ) 4 C 7 D Y S K J , K ) = + D E L ( J ) * C E L ( K ) / E L 4 C 8 I F U . E Q . K J D Y S I ( J , K ) = D Y S I ( J , K ) - 1 . / E L » 4 C 9 D Y S K J , K + 6 ) = - D Y S I ( J , K ) 4 1 C D Y S I ( J , K + 3 ) = C . 0 4^11 31 C C N T I N U E 4 1 2 3 0 C O N T I N U E M 3 CC 4 0 J = l , 3 4 1 4 DC 4 0 K = l , 3 4 1 5 N = J + 1 4 1 6 I F ( K . E Q . 4 ) f = l 4^17 N = rV + 1 4 1 8 I F ( N . E Q . 4 ) N=1 -V19 D B S ( J , K ) = W I I ( M , 3 ) * 0 Y S I ( N , K ) - W I I (N , 3 ) * D Y S I ( K , K ) 4 2 0 C C S ( J , K ) = W J J { N , 3 ) * D Y S I ( N , K ) - W J J ( N , 2 ) * D Y S I ( . V , K ) 4 2 1 D B S ( J , K + 6 ) = - C E S ( J , K ) 4 2 2 D C S ( J . K + 6 ) = - C C S ( J , K ) . 423 D E S ( J , K + 3 ) = C W I ( K , M , 3 ) * Y S I ( N ) - C W l ( K , N , 3 ) * Y S I ( M ) ;424 D C S ( J , K + 3 ) = C W J ( K , V , 3 ) * Y S I ( N ) - CWJ ( K , N , 3 ) *YS I ( * ) :4 2 5 4C C C N T I N U E .4 26 DO 50 K = l , 9 4 2 7 D Y E S J ( K ) = ( D B S ( 1 , K ) * B S ( 1 )+DSS ( 2 , K ) * E S ( 2 ) + D B S ( 3,K ) * B S ( 3 ) ) / Y E S J 1428 C Z E S J ( K ) = ( n C S ( l , K ) * C S ( l ) + C C S ( 2 , K ) * C S ( 2 ) + C C S ( 3 , K ) * C S ( 3 ) ) / Z C S J 14 29 DO 51 J = l , 3 1 4 3 0 D Y S J ( J , K ) = ( C E S ( J , K ) - D Y E S J ( K ) * Y S J ( J ) ) / Y E S J T 4 3 1 51 C Z S J ( J , K ) = ( C C S ( J , K ) - C Z E S J ( K ) * Z S J ( J ) ) / Z E S J 143 2 5C C C N T I N U E 1 4 3 3 CC 1 5 0 K = l , 9 1 4 3 4 CC 52 J = l , 3 14 3 5 N = J + 1 1 4 3 6 1 4 3 7 1438 •4 39^ ! 4 4 0 ' . 4 4 i 1 4 4 2 ' 1 4 4 3| 114 4 4 4 4 5^  -t 4 6 4 4 7^ *48; 4 4 9; 4 5C| ^ 5 1, ,4521 14 54 4 5 5 [456, L s4 57 14 58 1459; l 4 6 C i .46 1 46 2 •$6 3 . 4 6 4 4 6 5 •Abb' 4 6 7 .468 4 6 9 i,4 7 C 1471 147 2 i,4 7 3 A 7 4 : 4 7 5 1476 1^77 i4*7 8 W 7 9 1480 1481 1482 148 3 1484 K',8 5 I4 86 1487 1488 1 4 8 9 149G OA 9 1 149 2 1 4 9 3 |1494 1 4 9 5 52 1 5C 6C 7C 71 73 72 74 75 76 I F ( M . F 0 . 4 )M=1 N=P + 1 IF ( N . E 0 . 4 JN=1 D Y S K ( J , K ) = Y S I ( K ) * O Y S J (N , K ) + DYS I (M , K )* YS J < N ) - Y S I ( N ) *D Y S J ( K , K ) 1 - C Y S I ( N , K ) * Y S J ( M D Z S K ( J , K ) = Y S I ( N ) * D Z S J ( N » K ) + C Y S I ( M , K ) * Z S J ( N ) - Y S I ( N ) * D Z S J ( . V , K ) 1 - C Y S I ( N , K ) * Z S J ( N ' ) C C N T I N U E C C N T INUE DO 60 K = l , 3 108 C Y I S I ( K ) = V » I I ( 1 , 1 ) * C Y S I ( K ) + h I I ( 2 , 1 ) * D Y S I ( 2 , K ) + W I I ( 3 » 1 ) * D Y S I ( 3 , K ) ) + W J J ( 3 , 1 ) * D Y S I ( 3 , K ) D Z I S I ( K ) = W J J ( 1 , 1 ) * C Y S I ( 1 , K ) + V U J ( 2 , 1 ) * C Y S I ( 2 , K DY IS I (K + 6 ) = - D Y I S I ( K ) CZ IS I (K + 6 ) = - C Z I S I ( K } D Y I S I ( K + 3 ) =D W ( K , 1 , 1 ) * Y S I (1 ) + C W I ( K , 2 , l ) * Y S l ( 2 D Z I S I ( K + 3 ) = C V > J ( K , 1 , 1 ) * Y S I ( 1 )+DV>J ( K , 2 , 1 ) * Y S I ( 2 D Y I S J ( K ) = K I I ( 1 , 1 ) * C Y S J ( 1 , K ) + W I I ( 2 , 1 ) * D Y S J ( 2 , K D Y I S J ( K + 6 ) = - C Y I S J ( K ) D Z I S J ( K ) = U J J ( 1 , 1 ) * C Z S J ( 1 , DZ I S J ( K + 6 ) = - C Z I S J ( K ) . D Y I S J ( K + 3 ) = W I 1 ( 1 , 1 ) * D Y S J ( * ( 3 , K + 3 ) + D M ( K , l , l ) * Y S J ( l ) )+ D W I ( K , 3 , 1 ) * Y S 1 ( 3 ) )+ D W J ( K , 3 , l ) * Y S I ( 3 ) ) + W l I ( 3 , L ) * D Y S J ( 3 , K ) K ) + W J J ( 2 , l ) * 0 Z S J ( 2 , K ) + W J J ( 3 , l ) * D Z S J ( 3 , K ) l , K + 3 ) + W l I ( 2 , l ) * C Y S J ( 2 , K + 3 ) + W I I ( 3 , 1 ) * D Y S J + OWI ( K , 2 , 1 ) * Y S J ( 2 ) + CWI ( K , 3 , 1 ) * Y S J ( 3 ) D Z I S J ( K + 2 ) = V i J J ( l , l ) * D Z S J ( l , K + 3) + toJJ(2,l)*DZSJ(2,K+3)+WJJ(3,I ) * C Z S J 4 ( 3 , K + 3 ) + C W J ( K , l , l ) * Z S J ( l ) C Y I S K ( K ) = M I ( 1 , 1 ) * C Y S K ( 1 , + DWJ(K , 2 » 1 ) * Z S J ( 2 ) + D W J ( K , 3 , 1 ) * Z S J ( 3 ) K ) + W I I ( 2 , 1 ) * D Y S K ( 2 , K ) + W I I ( 3 , 1 ) * D Y S K ( 3 , K ) DZ I SK (K ) =V\ J J ( 1 , 1 ) * D Z S K ( 1, K )+ W J J ( 2 , 1) * D Z SK (2 , K ) + W J J ( 3 ,1 ) *D ZSK (3 , K ) D Y I S M K + 6 ) = - C Y ISK (K ) D Z I S K ( K + 6 ) = - C Z I S K ( K ) DY I S K ( K + 3 ) = H I ( 1 , 1) * C Y S K ( l , K + 3 )+W I I ( 2 , 1 ) * D Y S K ( 2 , K + 3 ) + W I I ( 3 , 1 ) * D Y S * K ( 3 , K + 3)+DVi I ( K , 1 , 1 ) * Y S K ( 1 )+DWI ( K , 2 , 1 ) * Y S K ( 2 ) +DW I ( K , 3 , 1 ) * Y S K ( 3 ) D Z I S K ( K + 3 ) = V> J J ( 1 ,1 ) * C Z S K ( 1 , K + 3 )+W J J ( 2 , 1 ) * D Z S K ( 2 , K + 3)+Vi J J ( 3 , 1) * D Z S * K ( 3 , K + 3 ) + D * . J ( K , l , l ) * Z S K ( l ) + D W J ( K , 2 , l ) * Z S K ( 2 ) + C W J ( K , 3 , 1 ) * Z S K (3 ) C C N T I N U E > DC 7 0 K = 1 , 9 A( 1 ,K ) = D E L ( K ) DC 71 K = l , 9 A ( 2 , K ) = ( - Y I S l * D Y I S K ( K ) + Y I S K * D Y I S I ( K ) ) / Y I S I * * 2 A ( 3 , K ) = ( Y I S I * O Y I S J ( K ) - Y I S J * 0 Y IS I ( K ) ) /Y IS 1**2 M= 1 N=3 CC 7 2 K=M,N A ( 5 , K ) = ( - Z I S I * C Z I S K ( K ) + Z I S K * C Z I S I ( K ) ) / Z I S I * * 2 A ( 6 , K ) = ( + Z I S I * D Z I S J ( K ) - Z I S J * 0 Z I S I ( K ) ) / Z I S I * * 2 A ( 4 , K ) = D Y S K ( 1 ,K ) * Z S J ( 1 ) + C Y S K ( 2 , K ) * Z S J ( 2 ) + D Y S K ( 3 , K ) * Z S J ( 3 ) + C Z S J ( 1 1 , K ) * Y S K ( 1 ) + C Z S J ( 2 , K ) * Y S K ( 2 ) + C Z S J ( 3 , K ) * Y S K ( 3 ) C O N T I N U E IF ( K . E C . l ) G C T C 7 4 G C T C 7 5 N = 7 N = 9 G C T C 7 3 DC 7 6 K = 4 , 6 A ( 5 , K + 6 ) = ( - Z I S I * C Z I S K ( K ) + Z I S K * C Z I S I ( K ) J / Z I S I * * 2 A ( 6 , K + 6 ) = ( Z I S I * C Z I S J ( K ) - Z I S J * C Z I S I ( K ) ) / Z I S l * * 2 A ( 4 , K ) = DYSK ( 1 ,K ) * Z S J ( 1 ) + C Y S K ( 2 , K ) * Z S J ( 2 ) +CYSK ( 3 , K ) * Z S J ( 3 ) A ( 4 , K + 6 ) = D Z S J ( 1 , K ) * Y S K ( 1 )+ D Z S J ( 2 , K ) * Y S K ( 2 ) + C Z S J ( 3 , K ) * Y S K ( 3 ) C C N T I N U E CC 80 J = l , 3 DO 80 K = l , 9 C C E L ( J , K > = - C Y S I ( J , K ) 1 4 9 6 C C E L ( J + 3 , K ) = C . O 1 4 9 7 C C E L ( J + 6 . K ) = - C D E L ( J , K ) 1 4 9 8 8 0 C C N T I N U E .4 9 9 DC 81 1 = 1 , 9 -15C0 CC 81 J = l , 3 1 5 0 1 DO 8 1 K = l , 3 15 C 2 C C Y S I ( I , J , K ) = + C C E L ( J , I ) * O E L ( K ) / E L + C C E L ( K , I ) * C EL ( J ) / E L 15C3 * - C E L I J ) * D E L ( K ) * C E L ( I ) / E l * * 2 1 5 C 4 I F ( J . E Q . K ) C C Y S I ( I , J ,K )=DOYS I (I , J , K ) + D E L ( I ) / ( E L * E L ) C C Y S I ( l , J , K + 6 ) = - C D Y S l ( I , J , K ) .5C6 C C Y S I ( I , J , K + 3 ) = C . C 5C7| 81 C O N T I N U E i 0 8 CC 8 4 N = l , 3 15C9 DC 8 4 K = l , 3 L 5 1 0 DC 84 1 = 1 , 3 l\511 DO 8 4 J = 1, 3 1 5 1 2 G C T 0 ( 8 5 , 8 4 , 86 ),K 1 5 1 3 85 I F ( I . E C . 1 ) C C V > I ( N , K , I , J ) = 0 . C 0 1/514 I F ( I . E G . 2 ) C D W l ( N , K , I , J ) = - C W l ( N , 3 , J ) 1515 I F ( I . E Q . 3 ) C C W I ( N , K , I , J ) = D W I ( N , 2 , J ) 1516 I F ( I . E C . 1 ) C C W J ( N , K , I , J ) = 0 . C G 1 5 1 7 I F ( I . E Q . 2 ) D D W J ( N , K , I » J ) = - D W J ( N ' , 3 , J ) 1518 I F ( I . E Q . 3 ) C C W J ( N , K , I , J ) = D W J ( N , 2 , J ) 1519; G C T C 8 4 l";52C 86 I F ( J . E Q . l )CCWI ( N , K , I , J ) = C W I ( N , I , 2 ) 1521 I F I J . E Q . l ) C D V J ( N , K » I , J ) = 0 W J ( N , I , 2 ) 15 22 I F ( J . E C . 2 )CCW I ( N , K , I , J ) = -CW I ( N , I , I ) l>52 3 I F ( J . E C . 2 ) C D V J ( N , K , I , J ) =-DW J ( N, I , 1 ) 1524 I F ( J . E G . 3 ) C 0 W I ( N , K , I , J ) = C . D 0 15 25, I F ( J . E Q . 3 1 C C W J ( N , K , I, J ) = 0 .CO 15.26 84 C O N T I N U E 1527 DC 87 J = l , 3 * 1528 C C W I ( 1 , 2 , 1 , J ) = 0 . D O 1 5 2 9 C C W J ( 1 , 2 , 1 , J ) = C . D O 1 5 3 0 C C W I ( 3 , 2 , J , 3 ) = C . C C 1531 C C W J ( 3 , 2 , J , 3 ) = C . D C 153 2 CCWI ( 3 , 2 , J , l ) = CWI ( 2 , J , 2 ) ^ 5 3 3 C C W J ( 3 , 2 , J , i ) = D W J ( 2 , J , 2 ) 1*5 3 4 CCWI ( 3 , 2, J , 2 ) = - C W I ( 2 , J , 1 ) 15 35 C C W J ( 3 , 2 , J , 2 ) = - C W J ( 2 , J , l ) 1,536 CCWI (2 , 2 , 1 , J ) = - W I I (1 , J ) V%27 C C W J ( 2 , 2 , 1 , J ) = - W J J ( l , J ) 1538 CCWI ( 1 , 2 , 2 , J ) = - D W I ( 2 , 3 , J ) l p 3 9 C C W J ( 1 , 2 , 2 , J ) = - C W J ( 2 ,3 , J ) 1'540 CCWI ( 1 , 2 , 3 , J ) = C W I ( 2 , 2 , J ) 1 5 4 1 C C W J ( 1 , 2 , 3 , J )= D W J ( 2 , 2 , J ) 154 2' 8 7 ' C C N T I N U E 1 5 4 3 CCWI ( 2 , 2 , 3 , 1 )= C 4 * S 5 * C 6 1 5 4 4 C C W J ( 2 , 2 , 3 , 1 ) = C 1 C * S U * C 1 2 , r 3 4 5 CCWI ( 2 , 2 , 3 , 2 ) = - C 4 * S 5 : * S 6 r5 4 6 C C V » J ( 2 , 2 , 3 , 2 J = - C 1 C * S 1 1 * S 1 2 1547 C C W I ( 2 , 2 , 2 , 3 ) = - C 4 * C 5 1 5 4 8 C C W J ( 2 , 2 , 3 , 3 ) = - C 1 C * C 1 1 1 5 4 9 C C W I ( 2 , 2 , 2 , 1 ) = - S 4 * S 5 * C 6 15 50 C C W J ( 2 , 2 , 2 , 1 ) =-S 1 C * S 1 1 * C 1 2 15 51 C C W I ( 2 , 2 , 2 , 2 )= S 4 * S 5 * S 6 1'552 C C W J ( 2 , 2 , 2 , 2 )= S 1 G * S 1 1 * S 1 2 1 5 5 3 C C W I ( 2 , 2 , 2 , 3 )= S 4 * C 5 1 5 5 4 C C U J ( 2 , 2 , 2 , 3 ) = C 1 1 * S 1 C 1 5 5 5 CC 8 2 1 = 1 , 3 1 5 5 6 DC 82 J = l , 3 1 5 5 7 DC 8 2 K = l , 3 5 5 8 M = J+1 1 1 0 11559 I F ( M . E G . 4 )M=1 , 1 5 6 0 N = M + 1 ' 1 5 6 1 I F ( N . E Q . 4 )N=1 156 2 C C E S ( I , J , K ) = W I I (M , 3 ) *DDYS I ( I , N , K ) -W I I ( N, 3 ) * C D YS I ( I , M, K ) ' 1 5 6 3 C C C S ( I , J , K ) = W J J ( M , 3 ) * C C Y S I ( I , N , K ) - W J J { N , 3 ) * C C Y S I ( I » M, K ) 1 5 6 4 D C P S ( I , J , K + 6 ) = - D D B S ( t , J , K ) 5 6 5 C C C S ( I , J , K + 6 ) = - C C C S I I , J , K ) ; 5 6 6 DCB S ( I , J , K+3 ) = C W I ( K , M f 3 ) * D Y S I ( N . I ) -DW I (K,N ,3 ) * D Y S I (M, I) 5 6 7 ' C D C S l I , J , K + 3 ) = C W J ( K , M , 3 ) * C Y S I ( N , I ) - D W J ( K , N , 3 ) * D Y S I ( M , I ) 568 C C B S ( 1+6 , J,!<) = WI I (M , 3 ) * D D Y S I (I + 6 , N , K )-W 11 ( N , 3 ) * D D Y S I ( 1 + 6 , N , K ) 1569 C C C S ( I + 6 , J , K ) = W J J ( M , 3 ) * D D Y S I ( I + 6 , N , K ) - w J J ( N , 3 ) * C D Y S I ( 1+6, M , K ) 1 5 7 0 C C B S ( 1 + 6 , J . K + 6 ) = - D O B S (1 + 6 , J , K ) l \571 C C C S ( 1+6 , J , K + 6 ) = - D D C S ( I +6, J , K ) 1 5 7 2 C C E S ( I + 6 , J , K 4 3 ) = CWI (K , M , 3 ) * D Y S I ( N , 1+6) - D W I ( K , N , 3 ) * D Y S I ( P , I+6) U573 C D C S ( 1 + 6 , J ,K + 3 ) = CW J { K , f*,3 ) * CYS I {N , I + 6 ) -CW J ( K , N , 3 ) * D Y S I ( M, I +6 ) L 5 7 4 C C B S ( 1 + 3 , J , K ) = D W K I , M , 3 ) * D Y S I ( N , K ) - O W I ( I , N , 3 ) * D Y S I ( M , K) 157 5 C C C S ( 1 + 3 , J , K ) = D W J ( I , M , 3 ) * C Y S I ( N . K ) - D W J I I , N , 3 ) * C Y S I ( M , K) 1 5 7 6 C C B S ( I + 3 , J , K + 3 ) = C C W I ( I , K , f , 3 ) * Y S l ( N J - D D W I ( I , K , N , 3 ) * Y S I ( t f ) 1V577 C C C S ( I + 3 , J , K + 3 ) = DDWJ( I , K , K , 3 ) * Y SI ( N ) - D D W J ( I , K , N , 3 ) * Y S I ( M ) 1 5 7 8 C C E S l 1 + 3 , J , K 4 6 )= - D D B S ( I + 3 , J , K ) 1,579 C C C S ( 1 + 3 , J , K + 6 ) = - CDC S < I +3 , J , K ) i ; 580 82 C C N T I N U E 1 5 8 1 DC 9 0 1 = 1 , 9 1 5 8 2 DC 91 K = l , 9 D-58 3 C C Y E S J ( K , I ) = ( D D B S ( I , 1 , K ) * B S ( 1 ) + C B S ( 1 ,K ) *DBS ( 1 , I ) +DDBS ( I , 2 , K ) * B S 1*5 8 4 * (2 ) + DBS ( 2 , K ) * C B S ( 2 , I ) + D D 8 S ( I , 3 , K ) * B S ( 3 ) + D B S ( 3 , K ) * D B S < 3 , I ) ) / Y E S J 1*585 * - C V E S J ( K ) * C Y E S J ( I ) / Y E S J 15 86 C C Z E S J ( K , I ) = ( C D C S ( I , 1 , K ) * C S ( 1 ) + CCS f 1 , K ) * 0 C S ( 1 , 1 ) + D D C S ( I , 2 , K ) * C S 15 87 * ( 2 ) + D C S ( 2 , K ) * C C S ( 2 , 1 ) + D D C S ( I , 3 , K ) * C S ( 3 ) + D C S ( 3 , K ) * D C S ( 3 , I ) ) / Z E S J 1588 * - C Z E S J ( K ) * C Z E S J ( I ) / Z E S J 1 5 8 9 DC 9 2 J = l , 3 1 5 9 0 C C Y S J ( I , J , K ) = ( C D B S ( I , J , K ) - D D Y E S J ( K , I ) * Y S J ( J ) - C Y E S J ( K ) * D Y S J ( J , I ) )/ £ 5 9 1 * Y E S J - D Y S J ( J , K ) * C Y E S J ( I ) / Y E S J 1 5 9 2 D D Z S J ( I , J , K ) = ( C C C S ( I , J , K ) - C O Z E S J ( K , I ) * Z S J ( J ) - C Z E S J ( K ) * D Z S J ( J , I ) ) / 5 9 3 * Z E S J - D Z S J ( J , K ) * D Z E S J ( I ) / Z E S J i , 5 9 4 92 C C N T I N U E 1 5 9 5 91 C C N T I N U E 1 5 9 6 90 CON T I N U E 1 5 9 7 DC 191 1 = 1 , 9 1 5 9 8 DC 1 9 0 K = l , 9 1^599 DC 93 J = l , 3 1 46C0 M = J + 1 1 6 C 1 I F ( M . E Q . 4 )P=1 1 6 C 2 N=N + 1 1 6 0 3 I F ( N . E Q . 4 )N= 1 1 6 C 4 C C Y S K { I , J , K ) = D Y S l ( M , I ) * D Y S J ( N , K ) + YS I (K ) * D D Y S J ( I ,N , K ) + D D Y S I ( I , M , K ) . . I 6 C 5 * * Y S J ( N ) + D Y S I ( N , K ) * C Y S J ( N , I ) - D Y S l ( N f I ) * C Y S J ( N , K ) - Y S I ( N ) * D D Y S J ( I , M , K 1 6 C 6 * ) - C D Y S I ( I , N , K ) * Y S J ( M ) - D Y S I ( N , K ) * r Y S J ( M , I ) 16 07 C C Z S K ( I , J , K ) = D Y S I ( M , I ) * D Z S J ( N , K ) + Y S I ( M ) * C D Z S J ( I ,N , K ) + D D Y S I ( I , M , K ) * 16C8 * Z S J ( N ) + C Y S I ( N , K ) * C Z S J ( N , I ) - D Y S I ( N , I ) * D Z S J ( y , K ) - Y S I ( N ) * D D Z S J ( I , N , K ) • 1-6C9 * - C D Y S I ( I , N , K ) * Z S J ( M ) - D Y S I ( N , K ) *-CZS J ( N , I) 16 10 93 C O N T I N U E •i:.6 11 19C C C N T I N U E 1 6 1 2 191 C C N T I N U E 16 13 DC 9 4 K = l , 3 1 6 1 4 CO 9 5 1 = 1 , 9 16 15 CCY IS I ( K , I ) = W I I ( 1 , 1 ) * C C Y S I ( I , 1,K )+W I I ( 2, I ) * C C Y S H I , 2 ,K )+WI I ( 3 , 1 )* 1 £ 1 6 * C C Y S I ( I , 3 , K ) 111 16 17 C C Z I S I ( K , I ) = W J J ( 1 , 1) * D D Y S I ( I , 1 , K ) + W J J ( 2 , 1 ) * D C Y S I ( I , 2 , K ) + W J J ( 3 , 1 ) * '.618 * C C Y S I ( I , 3 , K ) 16 19 C C Y I S J ( K , I ) = W I I ( 1 , 1 ) * C C Y S J ( I , 1 , K ) + W 1 1 ( 2 , 1 ) * C C Y S J ( 1, 2, K ) +W I 1 ( 3 , 1) * 1 6 2 0 * C C Y S J ( I , 3 , K ) 1621 C C Z I S J ( K , I ) = W J J ( 1 , 1 ) =*DD ZS J ( I , 1 , K ) + W J J ( 2 , 1 ) * D C Z S J ( I , 2 , K ) + W J J ( 3 , 1 ) * 1^2 2 * C C Z S J ( I , 3 , K ) 16 23 C C Y I S K ( K , I ) = W I I ( 1 , 1 ) * D D Y S K ( I , I ,K ) + W I I ( 2 , 1 ) * D C Y SK( I ,2 » K ) + WI I ( 3 , 1 ) * |l.624 * C C Y S K ( I , 3 , K ) 6 2 5 C C Z I S K ( K , I ) = W J J ( l , l ) * C C Z S K ( I , l , K ) + W ' J J ( 2 , l ) * C C Z S K ( I , 2 , K ) + W J J ( 3 , 1 ) * 6 2 6 * C C Z S K ( 1, 3 , K ) 6 2 7 95 C C N T I N U E 128 94 C C N T I N U E . 6 2 9 DC 1 9 5 K = l , 3 1630 DO 194 1 = 4 , 6 D6 3 1 N = I - 3 16 32 C C Y I S I ( K , I ) = C W I ( N , 1 , 1 ) * D Y S I ( 1 , K ) + D W U N , 2 , 1 ) * C Y S I ( 2 , K J + D W I ( N , 1 6 3 3 * 3 , 1 ) * D Y S I ( 3 , K ) 1,6 3 4 C C Z I S I ( K , I ) = CWJ ( N , 1, 1 ) * D Y S I ( 1 ,K ) +DW J ( N , 2 , I ) * C Y SI ( 2 , K ) +D W J ( N , L 6 3 5 * 3 , 1 ) * U Y S I ( 3 , K ) 1 6 3 6 C C Y I S J ( K , I ) = C C Y I S J ( K , I ) + D W I ( N , l , l ) * D Y S J ( l , K ) + 0 W I ( N , 2 , 1 ) * D Y S J ( 2 t'6 37 * ,K ) + DWl ( N , 3 , 1 ) * C Y S J ( 2 , K ) I t 3 8 C C Z I S J ( K , I ) = C C Z I S J ( K , I ) + C W J ( N , 1 , 1 ) * C Z S J ( 1 , K ) 4 C W J ( N , 2 , 1 ) * 0 Z S J ( 2 , K ) + 1 6 3 9 * D W J ( N , 3 , 1 ) * O Z S J ( 3 , K ) 1;&4C C C Y I S K ( K , I ) = C C Y I S K ( K , I ) + D W I ( N , 1 , 1 ) * D Y S K ( 1 , K ) + D W I (N * 2 «1 ) * D Y S K ( 2 , K ) + 1641 *DWI ( N , 3 , 1 ) * C Y S K { 3 , K ) 16 4 2 C C Z 1 S K ( K , I ) = C C Z I S K ( K , I ) + 0 W J ( N , 1, 1) * D Z S K ( 1 , K ) + D W J ( N , 2 , 1 ) * D Z S K ( 2 , K ) + '.',643 =*DWJ(N , 3 , 1 ) * C Z S K ( 3 , K ) 1644 194 C C N T I N U E 1645; 1 9 5 C C N T I N U E i|646 DO 196 K = l , 3 1647 DC 1 9 7 1 = 1 ,9 1648 96 C C Y I S I (K+6 , I ) = - C C Y l S I ( K , I ) L 6 4 9 C C Z I S K K + fc, I ) = - C D Z I S I ( K , I ) »,65C C C Y I S J ( K + 6 , I ) = - C C Y I S J ( K , I ) 1;651 C C Z I S J ( K + 6 , I ) = - C C Z I S J ( K , I ) 1652 C C Y I S K ( K - < 6 , I ) = - C C Y I S K ( K , I ) 115 3 C C Z I S K ( K + 6 , I ) = - C C Z I S K ( K , I ) i '654 M = K + 3 1 6 5 5 C C Y I S I ( M , I ) = C W I ( K , l , l ) * C Y S I l l , I ) + D W I ( K , 2 , l ) * C Y S I ( 2 , I ) + D W I ( K , 3 , l ) * ' 6 5 6 * D Y S I ( 3 , I ) 16 57 C C Z I S I ( M , I ) = C W J ( K , 1 , 1 ) * D Y S I ( 1 , I ) + D W J ( K , 2 , 1 ) * C Y S I ( 2 , I ) + 0 W J ( K , 3 , 1 ) * 1658 * C Y S I ( 3 , I ) ! £ 5 9 C C Y I S J ( M , I ) = WII ( 1 , 1 ) * C C Y S J ( I , 1 , M ) + WI I ( 2 , I ) * C CYS J ( 1 ,2 ,M)+WI I ( 3 , 1 ) * 1 ;66C * C C Y S J ( I , 3 , M ) + C W I ( K , 1 , 1 ) * D Y S J ( 1 , I ) + C W I ( K , 2 , 1 ) * D Y S J ( 2 , I ) +DWI ( K , 3 , 1) * 1 6 6 1 * C Y S J ( 3 , I ) 1 6 6 2 C C Z I S J ( M , I ) = W J J ( 1 , 1 ) * D D Z S J ( I , 1 , M ) + W J J ( 2 , 1 ) * C C Z S J ( I , 2 , K ) + W J J ( 3 , 1 ) * 166 3 * C C Z S J ( I , 3 , N ) + C W J ( K , 1 , 1 ) * D Z S J ( 1,I ) + D H J ( K , 2 , 1 ) * C Z S J ( 2 , I ) + D W J ( K , 3 , 1 ) * 1 6 6 4 ' * C Z S J ( 3 , I ) J ( : 6 5 CC YI SK( M, I ) = WI I ( 1 ,1 ) * C C Y S K ( I , 1 ,Ki) + W I I (2 , 1 ) * C C Y S K ( I, 2 ,M ) +W I I ( 3 , 1) * IS 66 * C C Y S K ( I , 3 , N ) + DWI ( K , 1 ,1 ) * D Y S K ( 1,1 )+DWl ( K , 2 , 1 ) * C Y S K (2 , I )+DWI (K. , 3 , 1) * 1667 * C Y S K ( 3 , I ) 1 6 6 8 C C Z I S K ( M , I ) = W J J ( 1 , 1 ) * C C Z S K ( I , 1 , M ) + K J J ( 2 , 1 ) * C C Z S K ( I , 2 , M ) + W J J ( 3 , 1 ) * 1 A 6 9 * C C Z S K ( I , 3 , V) + CWJ( K , 1 , 1 ) * D Z S K ( 1, I ) + D W J ( K , 2 , 1 ) * C Z S K ( 2 , I ) + D W J ( K , 3 , 1 ) * 1 6 7 C * C Z S K ( 3 , I ) 1,671 1 9 7 C C N T I N U E 1 6 7 2 196 C O N T I N U E 1 6 7 3 CC 198 K = l , 3 1 6 7 4 DC 199 1 = 4 , 6 1 6 7 5 M=K+3 1676 677 !l 6 7 8 6 7 9' .680 ' 168 1 "6-8 2! 1683 6 8 4 •85 a 8*6, ^87 i 6T6 a j 6891 69 C, Sv9 1 b92 &S3 6 9 4 6 9 5 6 9 6 6 9 7 6 V98 ; fr>99; T C C i 7 0 1 7 0 2 AC 3 7 0 4 7 C 5 > 06 7 C 7 7 0 8 7 C 9 '7 10 .71 1 7 12 f l 3 714 7 1 5 7 1 6 t i l 7 7 1 8 tl9 i20 7 2 1 7 2 2 7;2 3 7 2 4 .> 25 7 26 7 2 7 . 7 2 8 n 29 . 7 3 0 J ? 3 1 1/3 2 1 7 2 3 1 7 3 4 11735 N = I - 3 H 2 C C Y I S I ( K , I ) = C C W I ( N , K , l , l ) * Y S I ( l ) + C C W I ( N , K , 2 , i ) * Y S l ( 2 J+DDWI ( N » K » 3» 1 * ) * Y S I < 2 ) D C 2 I S I (M, I ) = CCWJ( N,K , 1,1 ) *YS I ( 1 ) + D C W J ( N , K , 2 , 1 ) * Y S I ( 2 ) +DD W J ( N , K , 3 , 1 • ) * Y S I ( 3 ) C C Y I S J ( M , I ) = C C Y I S J ( M , I ) + C W l ( N , l , l M C Y S J ( l , r M + C W I ( N , 2 , l ) * D Y S J ( 2 , M ) + • DWI (M , 3 , l . ) * D Y S J ( 3 ,M)+DDWI ( N , K , 1, 1) * Y S J (1 ) + DDWI ( N , K , 2 , 1 ) * Y S J ( 2 ) + • CCWI ( N , K , 3 , 1 ) * Y S J ( 3 ) C C Z I S J ( M , I ) = C C Z I S J ( M , I ) + D W J ( N , l , l ) * D Z S j ( l , n + C W J ( N , 2 , l ) * D Z S J ( 2 , M ) + • DWJ ( N , 3 ,1 1 * D Z S J ( 2 ,M) + D D W J ( N , K , 1 , 1 1 * Z S J (1 ) + 00 V» J ( N » K , 2 ,1 ) * Z S J ( 2 ) + • C C w J ( N , K , 3 , 1 ) * Z S J ( 3 ) C C Y I S K ( M , I ) = C C Y I S K ( K , D + C W I ( N , 1 , 1 ) * C Y S K ( I , K ) + C W I ( N , 2 , 1 ) * D Y S K ( 2 , M ) + • D W I ( N , 3 , 1 ) * C Y S K ( 2 » M ) + D D W I ( N , K , 1 , I ) * Y SKI 1}+ D0 W I ( N , K , 2 , 1 ) * Y S K ( 2 ) + * C C W I ( N , K , 3 , 1 ) * Y S K ( 3 ) C C Z I S K ( t f , I ) = C C Z I S K ( M , I ) + C W J ( N , 1 , 1 ) * C Z S K 1 1 , M ) + C W J ( N , 2 , 1 ) * D Z S K ( 2 , M ) + * D W J ( N , 3 , 1 ) * D Z S K ( 2 , M ) + D 0 W J ( N , K , 1 , 1 ) * Z S K ( 1 ) + D D W J ( N , K , 2 , 1 ) * Z S K ( 2 ) + * C C W J ( N , K , 3 , 1 ) * Z S K ( 3 ) 1 9 9 C C M I N U E 1 9 8 C C N T I N U E DO 110 1 = 1 , 9 DC 111 K = l , 9 V ( I ,1 , K ) = C C E L l I , K ) 111 C C N T I N U E 1 1 C C C N T I N U E DC 1 1 3 1 = 1 , 9 DC 114 K = l , 9 V I I , 2 , K ) = ( - Y I S I * D D Y I S K I K , I ) - D Y I S H I ) * D Y I S K ( K ) + Y I S K * D D Y I S I ( K , I ) + • D Y I S K ( IJ * C Y I S I ( K ) ) / ( Y I S I * * 2 ) - 2 . * A ( 2 , K ) * D Y I S I ( I ) / Y I S I V ( I , 3 , K ) = ( Y I S I * C C Y I S J ( K , I J + O Y I S I ( I ) * 0 Y I S J ( K ) - Y I S J * D D Y IS I ( K , I ) -• D Y I S J ( I ) * D Y I S I ( K ) ) / Y I S 1 * * 2 - 2 . * A { 3 , K ) * D Y I S I ( I J / Y I S I M= I N = K I F ( I . G E . 4 . A N C . I . L E . 6 ) K = I + 6 I F ( K . G E . 4 . A N C . K . L E . 6 ) N = K + 6 V ( M , 5 , N ) = ( - Z I S I * C C Z I S K ( K , I ) - D Z I S I ( 1 ) * D Z I S K I K ) + ZI SK*DDZ. IS I ( K, I ) + * D Z I SK I I ) * D Z IS I (K ) ) / Z l S I * * 2 - 2 . * M 5 , M * 0 Z I S I ( I J / Z I S I V ( M , 6 , N ) = ( Z I SI * C D Z I S J ( K , I )+DZ I SI ( I ) * Q Z I S J ( K ) - Z I S J * D D Z I S I ( K , I ) -* D Z I S J ( I ) * C Z IS I (K ) ) /Z IS I < * 2 - 2 . V A ( 6 , N ) * D Z I S I ( I ) / Z I S I 114 C C N T I N U E 113 C C N T I N U E M=l N=3 M l = l N 1 = 3 1 1 9 DC 1 1 5 I=f ,N DC 116 K = M , M V ( I , 4 , K ) = C C Y S K ( I , 1 , K ) * Z S J ( 1 ) + C Y S K ( 1 ,K ) * D Z S J ( 1 , 1 ) +DDY S K ( I , 2 , K ) * Z S * J ( 2 ) + D Y S K ( 2 , K ) * C Z S J ( 2 , I ) + D D Y S K ( I , 3 , K ) * Z S J ( 3 ) + D Y S K ( 3 , K ) * D Z S J ( 3 , I ) • + C C Z S J ( I , l , K ) * Y S K ( l ) + C Z S J ( l , K ) * C Y S K ( l , I) + D D Z S J ( I , 2 , K ) * Y S K ( 2 ) + D Z S J * ( 2 , K ) * D Y S K ( 2 , I ) + C C Z S J ( I , 3 , K ) * Y S K ( 3 ) + D Z S J ( 3 , K ) * C Y S K ( 3 , I ) 1 1 6 C O N T I N U E 1 1 5 C C N T I N U E I F ( . E G . 1 ) G C T C 1 1 7 GOTO 118 1 1 7 K=7 N = 9 G C T 0 1 1 9 118 M= 1 N = 3 1 F ( H 1 . E G . I 1 G C T C 1 3 0 1743 174 4! ' 745 / 4 6 7 4 7 '748 7 4 9 .750 751 7 5 21 i,2 53i St 5 4 L7 5 5 1756 7 57 .1*58 1759, i"?6o: .761 176 2 7 6 3 .764 L76 5;' [7,6 6 . 7 6 7 1768 1769 1-7 7 0 1?71 1772 U7 7 3 77 4 1775 '.77 6 1^77 1778 ' 7 7 9 17 8 C 178 1 1782 17 8 3 78 4 85 V7 6 6 1787 .7 8 8 1 " 8 9 2 CF 13C 124 121 12C 122 123 126 131 132 134 1 3 3 1 3 5 136 138 1 3 7 M 1 = 1 N l = 3 G C T C 1 2 4 M l = 7 N 1 = 9 G C T 0 1 1 9 DC 120 1 = 4 , 6 DC 121 K = M , M V < I , 4 , K ) = C C Y S K ( 1 , 1 , * ( 3 ) + C Z S J ( 1 , K ) * C Y S K ( C C N T I N U E C C N T I N U E I F ( M l . E G . 1 1 G C T C 1 2 2 GC TC 1 2 3 M l = 7 N l = 9 G C T C 1 2 4 DC 126 1 = 4 , 6 DO 126 K = N 1 , N 1 V ( I + 6 , 4 , K ) = C Y S K ( 1 , K * ( 3 , 1 ) + D D 2 S J ( I , 1 , K ) * Y CONT INUE I F ( H 1 . E 0 . 7 1 G C T C 1 3 1 G C T C 1 3 2 M l = l N l = 3 G C T C 1 2 3 DC 133 I = P 1 , M DO 134 K = 4 , 6 V ( 1,4 ,K )= C C Y S K ( 1 , 1 , * ( 2 ) + D Y S K ( 2 , K ) * C Z S j ( 2 V ( I , 4 , K + 6 ) = C C Z S J ( 1 , 1 * K ( 2 ) + D Z S J ( 2 , K ) * C Y S K ( C C N T I N U E C C N T INUE I F ( M 1 . E Q . 1 1 G C T 0 1 3 5 G C T 0 1 3 6 M l = 7 N 1 = 9 G C T 0 1 3 2 CC 137 1 = 4 , 6 DO 138 K = 4 , 6 V ( I , 4 , K ) = C C Y S K I 1 , 1 , * ( 3 ) V ( I + 6 , 4 , K ) = C Y 5 K ( 1 , K ) * 3 , I ) V ( I , 4 , K + 6 ) = C Z S J ( 1 , K * ( 3 , I ) V ( I + 6 , 4 , K + 6 ) = C C Z S J ( I * Y S K ( 3 ) C C N T I N U E CCNT INUE RETURN EN C 113 K ) * Z S J ( 1 ) + C D Y S K { I , 2 , K ) * Z S J ( 2 ) + D D Y S K ( I , 3 , K ) * Z S J 1 , I ) + D Z S J ( 2 , K ) * C Y S K ( 2 , I ) + C Z S J ( 3 , K ) * D Y S K ( 3 , 1 ) ) * D Z S J ( l , I ) + O Y £ K ( 2 , K ) * D Z S J ( 2 , I ) + O Y S K ( 3 , K ) * D Z S J S K ( 1 J + C D Z S J l I , 2 ,K 5 * Y S K ( 2 ) + C C Z S J ( I , 3 , K ) *YSK ( 3 ) K ) * Z S J ( 1 ) + C Y S K ( 1 , K ) * D Z S J ( 1 , 1 ) + D D Y S K ( I , 2 , K ) * Z S J , I ) + D D Y S K < I , 3 , K ) * Z S J ( 3 ) + C Y S K ( 3 , K ) * D Z S J ( 3 , 1 ) , K ) * Y S K ( 1) + D Z S J ( l , K ) * D Y S K ( 1 , I ) + D D Z S J ( I , 2 , K ) * Y S 2 , I ) + D D Z S J ( I , 3 , K ) * Y S K ( 3 ) + D Z S J ( 3 , K ) * D Y S K ( 3 , I ) K ) * Z S J ( 1 ) + C D Y S K { I , 2 , K ) * Z S J ( 2 ) + D D Y S K ( I , 3 , K ) * Z S J * C Z S J ( 1 , I J + C Y S K ( 2 , K ) * D Z S J ( 2 , I ) + D Y S K ( 3 , K ) * D Z S J ( ) * C Y S K ( 1 , I ) + D Z S J ( 2 , K ) * D Y S K ( 2 , I )+D ZS J ( 3 , K ) *D YSK , 1 , K ) * Y S K ( 1 ) + 0 C Z S J ( I , 2 , K ) * Y S K ( 2 ) + D D Z S J ( I , 3 , K ) * F I LE I C 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0050543/manifest

Comment

Related Items