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. U n i v e r s i t y o f 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 i n the Department of Civil  Engineering  We accept t h i s thesis as conforming to the required standard:  THE UNIVERSITY OF BRITISH COLUMBIA August 1972  In p r e s e n t i n g t h i s  thesis  in p a r t i a l  f u l f i l m e n t o f the requirements f o r  an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, the L i b r a r y s h a l l I  f u r t h e r agree t h a t p e r m i s s i o n  for  for  r e f e r e n c e and  f o r e x t e n s i v e copying o f t h i s  that  study. thesis  s c h o l a r l y purposes may be granted by the Head of my Department o r  by h i s of  make i t f r e e l y a v a i l a b l e  I agree  this  written  representatives. thesis  i s understood that copying o r p u b l i c a t i o n  f o r f i n a n c i a l gain s h a l l  permission.  Department o f  Civil  Engineering  The U n i v e r s i t y o f B r i t i s h Vancouver 8, Canada  Date  It  Columbia  September .11, 1972.  not be allowed without my  i  ABSTRACT E l a s t i c s t r u c t u r e s e x h i b i t i n s t a b i l i t i e s which a r i s e through the occurrence of f i n i t e displacements even when c o n s t i t u t i v e properties r e main l i n e a r .  A n o n - l i n e a r analysis which recognizes r o t a t i o n s i n the  s t r a i n displacement r e l a t i o n s h i p i s formulated f o r analyzing threedimensional framed s t r u c t u r e s . A f i n i t e element method i s used whereby the r o t a t i o n s w i t h i n each element are r e s t r i c t e d in s i z e by use of a l o c a l element reference frame attached to the element.  Two such coordinate systems are developed.  Then  an incremental s o l u t i o n technique based on an instantaneous l i n e a r i z a t i o n of a Taylor s e r i e s expansion of the forces about the displacement c o n 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 e c t on the e q u i l i b r i u m paths of the type of moving coordinate frame, the number of elements, and the s i z e of the increment step.  ii  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  III  IV  V  VI  12  DISPLACEMENT RELATIONSHIPS  17  1.  Global Degrees of Freedom f o r an Element  2.  Transformation Matrices  3.  The I n i t i a l Coordinate Frame ( y )  24  4.  The Incremental Structure Deformations 5r  30  and  if-*  if-  17 19  K  0  TANGENTIAL REFERENCE FRAME  33  1.  D e f i n i t i o n of the Element Degrees o f Freedom  33  2.  Derivation of the T r a n s ! a t i o n a l  33  3.  Derivation of the Rotational Components of s  s  36  SECANT REFERENCE FRAME  39  1.  The Element Degrees of Freedom §  39  2.  Derivation of s  as a Function of  £  41  FORCE-DISPLACEMENT RELATIONSHIPS  45  1.  Element S t i f f n e s s  45  2.  Global System E q u i l i b r i u m  48  3.  Derivation of  4.  Incremental S o l u t i o n Technique  a*  and  v*  Matrices  51 52  iii  Chapter  Page  VII  BIFURCATION BUCKLING  54  VIII  NUMERICAL EXAMPLES  57  1.  W i l l i a m s ' Toggle  57  2.  A r g y r i s ' Arch  64  3.  Wright's R e t i c u l a t e d 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-  3.3  4*^"  3.4  Transformation  22  Transformation  22  ^^Transformation  22  X  3.5  I n i t i a l Coordinate Frame  26  3.6  S i n g u l a 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  f o r Tangential Frame  34  5.1  Deformations  f o r Secant Frame  40  5.2  Rotational Degrees o f Freedom  8.1  W i l l i a m s ' Toggle with Central Load  8.2  Incremental D e f l e c t i o n Solutions f o r W i l l i a m s ' Toggle  s  s  2  and  s  3  44 58  63  V  Figure  Page  8.3  A r g y r i s ' Arch  65  8.4  Load D e f l e c t i o n Curve f o r A r g y r i s ' Arch  67  8.5  Wright's Shell Segment  70  8.6  Load D e f l e c t i o n Curve of W r i g h t ' s Segment with Points 6 f i x e d against t r a n s l a t i o n  72  8.7  Load D e f l e c t i o n Curve of Wright's Segment w i t h Points B f r e e to t r a n s l a t e  73  8.8  Three-Dimensional Elbow  75  8.9  Ring Dome  79  8.10  Load D e f l e c t i o n Curve f o r Ring Dome  80  8.11  Deflected c o n f i g u r a t i o n of a Dome Rib  80  vi  LIST OF TABLES  Table  Page  1  Solutions f o r W i l l i a m s ' Toggle  2  W i l l i a m s ' Toggle - 0.01 inch displacement  60  increments  61  3  W i l l i a m s ' Toggle - 10 elements per h a l f  61  4  Argyris  Arch  66  5  Wright's S h e l l Segment  71  6  Three Dimensional Elbow  76  1  ACKNOWLEDGEMENTS  I am g r a t e f u l f o 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 n a n c i a l assistance.  CHAPTER I INTRODUCTION C l a s s i c a l s t r u c t u r a l a n a l y s i s implies a unique s o l u t i o n to every s t r u c t u r a l problem since the theory i s based on i n f i n i t e s i m a l l i n e a r d i s placements.  In f a c t , real structures e x h i b i t i n s t a b i l i t i e s associated  with non-unique s o l u t i o n s and, i n order to detect these, i t i s to introduce n o n - l i n e a r  necessary  analysis.  I n s t a b i l i t i e s may a r i s e through n o n - l i n e a r material properties or through the occurence of f i n i t e displacements even when the c o n s t i t u t i v e properties remain e l a s t i c .  We confine our a t t e n t i o n 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.  Bifurcation,  2.  Snap-through,  3.  F i n i t e Disturbance.  The f i r s t may be detected by the s o l u t i o n of the c l a s s i c a l  eigenvalue  problem which has been formulated to include r o t a t i o n s of i n f i n i t e s i m a l elements i n the e q u i l i b r i u m equations of e l a s t i c i t y .  Recognition of the  second and t h i r d requires a s o l u t i o n technique capable of t r a c i n g the e q u i l i b r i u m path a f t e r the advent of f i n i t e displacements.  This n e c e s s i -  tates the r e c o g n i t i o n of r o t a t i o n s i n the s t r a i n displacement r e l a t i o n ships.  Frame s t r u c t u r e s e x h i b i t e i t h e r b i f u r c a t i o n or snap-through  b u c k l i n g , and we l i m i t ourselves here to the study of such s t r u c t u r e s . In the f i n i t e element analysis of frame s t r u c t u r e s , large rotations may be dealt with i n one of two methods: 1.  The member properties may be deduced, with  2  respect to a reference frame f i x e d in space, on the basis of the f u l l n o n - l i n e a r s t r a i n displacement r e l a t i o n s . 2.  The member properties may be deduced, with respect to a moving reference frame attached to the members i n question.  Then, by sub-  d i v i d i n g the s t r u c t u r e i n t o a s u f f i c i e n t number of members or "elements", we may r e s t r i c t the r o t a t i o n s w i t h i n each element r e l a t i v e 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 r o t a t i o n s and s t r a i n s which are small compared to the rotations.  The assumptions  are adequate f o r the detection of b i f u r c a t i o n  points i n the e q u i l i b r i u m 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 s t u d i e d , the moving coordinate systems as desc r i b e d above are required.  Non-linear problems are then solved by an i n -  cremental procedure based on an instantaneous  l i n e a r i z a t i o n w i t h i n each  increment of a Taylor s e r i e s expansion of the force vector about the d i s placements at the beginning of the increment. Early work on the t h e o r e t i c a l a n a l y s i s of e l a s t i c post-buckling was performed by K o i t e r (2).  His work, using a continuum mechanics approach,  centered on the i n v e s t i g a t i o n of imperfection s e n s i t i v e s t r u c t u r e s . B r i t v e c and Chi 1 vers (5) developed matrix methods based on a p o t e n t i a l energy formulation to analyze the i n i t i a l post-buckling curves of r i g i d l y j o i n t e d plane frames.  Martin (6), Supple ( 7 ) , and Roorda ( 8 ) , ( 9 ) , have  a l s o done research i n t o 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  U c c i f e r r o (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 r u c t u r a l problems. confined to planar structures o n l y .  T h e i r work was  A complete summary of the f i n i t e  element a n a l y s i s of nonlinear s t r u c t u r e s i s given by M a l l e t t and Marcal (3).  The present work takes element s t i f f n e s s matrices from Nathan  (11).  Two-dimensional s t r u c t u r e s are studied with a view to documenting computational experience.  Snap-through buckling of frames and arches  is  i n v e s t i g a t e d with respect to f a c t o r s such as element s i z e and number, increment s i z e , and choice of element moving reference system. The necessary r e l a t i o n s h i p s are developed f o r extending the work of previous i n v e s t i g a t o r s i n t o three-dimensions.  A simple three-dimensional  space dome element i s studied and r e s u l t s compared to t h e o r e t i c a l work by Wright (13). studied.  The snap-through buckling of a large r i n g dome i s then  No attempt fs made to study b i f u r c a t i o n buckling of the  s t r u c t u r e s presented. occurring.  Such modes of i n s t a b i l i t i e s are prevented from  4  CHAPTER II COORDINATE SYSTEMS 1.  Configuration Space We r e s t r i c t the class o f 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 s t r u c t u r e can then be subdivided i n t o any number of  l i n e elements connected to each other at "nodes".  Displacements of a l l  points w i t h i n each element are completely determined by the displacements of the nodes.  Thus, i f there are a t o t a l of n  nodal displacements or  "degrees of freedom" defined by the vector r , the c o n f i g u r a t i o n of the s t r u c t u r e w i l l be completely determined by the p o s i t i o n of the point with coordinates r  i n an h-space, the s o - c a l l e d " c o n f i g u r a t i o n space".  We are at l i b e r t y to choose the l o c a t i o n s of our nodes and thus of the number and s i z e s of our elements.  Such choice as t h i s represents the  e s s e n t i a l step i n the mathematical i d e a l i z a t i o n of the s t r u c t u r e . Obviously there must be s u f f i c i e n t c o n s t r a i n t s on the degrees of freedom to provide f o r global e q u i l i b r i u m of the s t r u c t u r e as a whole. We begin by d e f i n i n g the necessary coordinate systems i n which we measure displacements and forces and then we r e l a t e these reference frames to each other. We require a global c a r t e s i a n reference frame ( X ) X.  W1  '*-  n  base vectors  whose o r i g i n and axes o r i e n t a t i o n s are q u i t e a r b i t r a r i l y f i x e d i n  space.  We measure a l l nodal coordinates and global forces i n t h i s f i x e d  frame. frame we envision an a r b i t r a r y member TK  ,  5  connecting nodes j  to k ;  at each of these two nodes are f i x e d l o c a l  global reference frames.  These two frames, termed the  frames with base vectors  X  with reference to the tively.  and  X  (X ) J  and(x ) K  r e s p e c t i v e l y , have o r i g i n s  frame, by p o s i t i o n vectors  j?  J  and P  fixed, k  The p o s i t i o n vectors describe the distance from the (^X. ) ^  o r i g i n to the i n i t i a l p o s i t i o n of the nodes j  respectively.  r a m e  The  (X )  frame axes are f i x e d p a r a l l e l to the global reference  frame Q C ) axes.  Within these l o c a l global frames we measure the nodal  U ) J  and  and k  respec-  k  deformations of the s t r u c t u r e T . original ( X )  i  a x e s  n  a  3  &  and  r 1  to  r"" 1  are taken about the  defined order and the t r a n s l a t i o n s  the global axes d i r e c t i o n s where:  <~  The rotations  are i n  6  At node j of member J K .  we define an i n i t i a l l o c a l reference frame A  named  with base vectors  0  *j .  The coordinate axes o f t h i s frame  are defined by the o r i e n t a t i o n o f the undeformed member i t s e l f ;  / .Jv  Ai  vectors )  X  A„  o f the I X ) frame are r e l a t e d to the y  by an orthonormal transformation ^°  where  the base  base vectors o f  ^- :  A j  a  *4- i s a matrix f u n c t i o n of the constant angular r o t a t i o n s  ( r^.°  r ^ ° *~ ° ) s  d e f i n i n g the i n i t i a l o r i e n t a t i o n o f the member.  T  fc  The d e r i v a t i o n of the  transformation and the determination of the  f i n i t i a l angular r o t a t i o n s  °  \ *H  o\ T *\, J  w i l l be covered i n  Chapter I I I . Now we must r e l a t e the p o s i t i o n of the deformed member 3" < to the l o c a l global reference frames systems  (y )  (xM  and ( X ) by two moving coordinate  and (y ) with base vectors  K  \j  nodes J  and K. r e s p e c t i v e l y . A i  base vectors  y  where fc*  ^ r  =  The base vectors  "~  X  by an orthogonal transformation  follows i n Chapter  y  and y * a f f i x e d to A \  I  are r e l a t e d to the  j  V~  whose d e r i v a t i o n  III:  t  -  i s a matrix function of the large global  s*  r  fo*)  T  w h i c h  correspond t o  (  rotations  r* °  r i ) 0  s  T  and  which define the current o r i e n t a t i o n of the member tangent at node j . S i m i l a r l y , at node K we have  5 =f K  where A f  x  i s a matrix function of large rotations  /  (  r*  ^* ^* l o  "*\T j .  In order to evaluate the  F  we require to add the e f f e c t of deform-  ations to the angles d e f i n i n g the i n i t i a l o r i e n t a t i o n of the member. ever, £ ° and  How-  r*~ are not d i r e c t l y a d d i t i v e and therefore we introduce  intermediate measures of deformations member and are given by  -  r ~jr  V  r  .  These are p e c u l i a r to a given  o  - r  , where  r  2  lo  The £  If.  are r e l a t e d to the s t r u c t u r e nodal deformations  mental r e l a t i o n s h i p  £ by an i n c r e -  8  where Sr.  The global t r a n s l a t i o n s  r  s t r u c t u r e nodal t r a n s l a t i o n s of r  f o r the element are equivalent to the but the element r o t a t i o n s  r e l a t e d to the s t r u c t u r e nodal r o t a t i o n s of r equation given above, where  -t  r  =  G  r*" can be  only by the incremental  9  and  3  r12  Each member has associated with i t a l o c a l moving coordinate system i n which element displacements are measured and i n which the s t i f f n e s s matrix i s formulated.  The l o c a l moving frame base vectors are functions of the  global deformations  r  The element s t i f f n e s s matrix can be transformed to t h e Q ^ system and thus, u l t i m a t e l y , the response of the member i n global coordinates can be r e l a t e d to the response of the s t r u c t u r e i n l o c a l element coordinates. 2.  The Global Coordinate Systems A right-handed c a r t e s i a n reference frame c a l l e d the s t r u c t u r e global  coordinate s y s t e m ( X ^  i s shown in Figure 2.1.  This a r b i t r a r y reference  frame remains f i x e d i n space and a l l geometry u l t i m a t e l y r e l a t e s to t h i s coordinate system. nodes J  and k  Let  P  J  of member T <  and  p  be vectors d e f i n i n g the p o s i t i o n of  r e s p e c t i v e l y where:  10  A vc  FIGURE 2.1  GLOBAL COORDINATE SYSTEMS  11  and  (2.1)  i  P if  F*  and  K  p*  are i n i t i a l p o s i t i o n v e c t o r s , then the s t r a i g h t  element J k L has vector d i r e c t i o n  P^-  (x*) and ( X * ) defined at nodes j Aj A X J and x where: A. J X  p* and k  .  line  Node coordinate systems  r e s p e c t i v e l y have base vectors  - I  X  A X  J  -  2  -  3  and  (2.2) x  K  X  - 3  A j  A <  X  The components of the  X  and  are themselves v e c t o r s , they define  the vector d i r e c t i o n of each of the coordinate axes. base vectors J<_ , x  ,  X  We assert that the  are the same since the X  and  X  base  vectors remained f i x e d i n o r i e n t a t i o n throughout a l l and any load d i s p l a c e ment h i s t o r y of the s t r u c t u r e . f i x the o r i g i n s of the  (x^  The i n i t i a l p o s i t i o n vectors and C x k )  systems.  and  We measure the nodal  degrees of freedom f* or s t r u c t u r e (global) generalized displacement coordinates i n these node coordinate systems. 3.  Element Coordinate Systems S t r u c t u r a l problems i n v o l v i n g f i n i t e displacements can be solved, as  p r e v i o u s l y s t a t e d , by one of the f o l l o w i n g procedures.  F i r s t l y , we can  include the e f f e c t s of f i n i t e r o t a t i o n s w i t h i n the element boundaries when d e r i v i n g the element s t i f f n e s s .  Secondly, we can subdivide the s t r u c t u r e  i n t o many elements so as to reduce the r o t a t i o n s and t r a n s l a t i o n s of the element to w i t h i n small acceptable bounds.  A l e s s r e f i n e d element s t i f f -  ness i s used and the problem of the large r o t a t i o n s i s handled by the use o f moving member coordinate systems f i x e d to the elements i n question. The l a t t e r procedure i s employed in t h i s work.  We use an incremental  load and displacement technique to follow the load-displacement behaviour of the s t r u c t u r e .  At each increment step we r e c a l c u l a t e the s t i f f n e s s of  each element and reassemble the global s t i f f n e s s , an instantaneous ness tangent to the real load-displacement surface.  stiff-  We define a moving  member coordinate system which i s f i x e d 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 c a l l e d the tangential reference frame, and the second  c a l l e d the secant reference frame are studied h e r e i n . The tangential reference frame . 6j*") has u n i t base vectors given by ^  where: A  t  S  (2.3)  - 2  _3  /  -fc\  The  ^cj ) frame has i t s o r i g i n at node  The  \j - 2  '  j of the element J K .  "  A  .  The y i base vector i s tangent to the c e n t r o i d a l axis of the element a t n o d e j . A t  base vector i s coincident with the j while the  element cross s e c t i o n at node w i t h the minor p r i n c i p a l axes.  base vector i s coincident  S t i f f n e s s matrices are presented f o r r e c -  tangular cross sections only. frame.  major p r i n c i p a l axis of the  Figure 2.2 i l l u s t r a t e s t h i s reference  A l l element deformations i n t h 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 j h a s unit base vectors given by S  tj  where:  ^  s  (2.4) A  y -  In t h i s frame the ^ to node K. .  s 3  A 3  M  base vector i s defined by the vector j o i n i n g node  The p r i n c i p a l axes of the cross section at node j  i n general normal to the  ^  axis.  are not  A r b i t r a r y d e f i n i t i o n s are therefore  FIGURE 2.2  TANGENTIAL ELEMENT COORDINATE SYSTEM  given below f o r the d i r e c t i o n s of the coordinate axes A f t e r given displacements at nodes j  and K  L  $  2  and  we can c a l c u l a t e the  base vector as:  where  q  f ci, ci  =  2  5 ^  T  a r e  t n e  3  the i n i t i a l member J K .  projections on the global axes of  Figure 2.3 shows the secant coordinate frame.  We  define the other base vectors by Equations 2.6 and 2.7. A  tj  5  -  2  A  ^  =  S  "  3  A  =  3  / / A^.  s  VJ, x  A  A \ S  7  7~fc  A \  (2.6)  S  5  y  (2.7)  In t h 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 f o r an Element Each element has twelve global degrees of freedom, s i x a t each node,  which are defined by the corresponding s t r u c t u r e degrees of freedom £" at that node.  At each node there are three t r a n s ! a t i o n a l degrees of freedom  a c t i n g along the (X ) coordinate system axes and three r o t a t i o n a l degrees of freedom whose axes of r o t a t i o n depend upon the o r i e n t a t i o n of the member i n space.  Since the i n t e n t i o n i s to deal with f i n i t e r o t a t i o n s , the  d e f i n i t i o n s of these axes of r o t a t i o n i s rather complex.  Figure 3.1 shows  the s i x degrees of freedom f o r node J  where s i n g l e  of an element J K  arrows denote t r a n s l a t i o n s and double arrows denote r o t a t i o n s . —t We c a l l the s i x element t r a n s ! a t i o n a l degrees of freedom  r  where  - -t  P"  i s defined i n Chapter II r  as  2  r * f c  r  6  These t r a n s l a t i o n s can be r e l a t e d to the element or l o c a l displacements of the element by a vector transformation. tions  where r  i s defined in Chapter II  — r  r = i*  However, the f i n i t e global r o t a -  ^zy  as  FIGURE 3.1  GLOBAL DEGREES OF FREEDOM AT NODE  j  do not transform as vectors.  Problems of t h i s type can be handled by  s p e c i f y i n g a p a r t i c u l a r order f o r the r o t a t i o n s or by the use of E u l e r angles.  However, since there i s no s i n g l e system of Euler angles covering  a whole sphere without ambiguity, the former method w i l l be used here.  It  i s to be noted that three f i n i t e r o t a t i o n s taken i n an order d i f f e r e n t from that s p e c i f i e d w i l l  lead to a d i f f e r e n t o r i e n t a t i o n i n space.  Therefore, the r o t a t i o n s 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 t h i s order w i l l be adhered to throughout the development of t h i s t h e s i s .  I t i s asserted that t h 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 i n problem input and i n t e r p r e t a t i o n of r e s u l t s . 2.  Transformation Matrices ¥ The  H* -^  and  Ambiguities can thus be avoided. ±  matrix was introduced i n Chapter II  as a matrix function of  ~  ^j  r o t a t i o n s r e l a t i n g the base vectors of two coordinate systems, the of the The  (x ) J  UJ  w  in  frame and an a r b i t r a r y  frame with base vectors  be the transformation of the  X  vectors a t node j  of element T K  5 =V  X. .  vectors i n t o the  .  J  We have defined ( n/  r* s  C)  ^  X  J  previously as a matrix function of f i n i t e r o t a t i o n s where:  T  o  0  a'  (3.1,  The  (  ITs  define the i n i t i a l member o r i e n t a t i o n i n space.  propose to order the rotations as f o l l o w s : finally,  Yl  first  X*  » then t%  We  , and,  .  F i r s t of a l l , the  t*.  *-  A  r o t a t i o n transforms from the base vectors  to a s e t of instantaneous base vectors  X  X  J  of an intermediate coordinate  frame by Equation 3.2: A  X  I  =  ,1  V-  A  X  i  (3.2)  I  V-  I  i s a s i n g l e r o t a t i o n as shown by Figure 3.2 and ^ /  o  0  CosCrf)  o  -S»h(V)  i s given by:  O  Sn  (r *) H  CosCrf)  and  A  X  r  :  -  2  X  T  - 3  Secondly, the r o t a t i o n transforms from a known set of base vectors Aj ATJ X ' to the instantaneous coordinate frame base vectors X of an i n t e r mediate frame  ( X ^ ) by Equation 3.3: (3.3)  where  Co C 5  O  r*)  o /  O  -S,n(V) o  Co (r *) S  5  21  and  ' A H  N  A H  X  -  -2  AI  X  - 3  The  r o t a t i o n was about an instantaneous axis defined by the  base vector.  *  F i n a l l y , the  A  3  base v e c t o r , transforms to the base vectors A  from the known base vectors A A , TT_  j  m. y = V- x  X  I  r o t a t i o n , which i s about the instantaneous  Al  X  A  y  j.  j  o f the [y  )  frame  / n \ "~  TT  X  of the  IX  ) frame by Equation 3.4: (3.4)  where Sin(V)  CesCrJ*)  O  (V)  - Sin  r  s  1  o  o  The r o t a t i o n s  and  O  are shown i n Figures 3.3 and 3.4.  From Equa-  tions 3.2, 3.3 and 3.4 we derive Equation 3.5:  9» A  v-  =  ;  yJ  =  v/J  j  A  V  1  :  x  where  vf  J =  ^  ^  ^  (3.5)  and Cos{tSi  CosCrtf) -r5ir,(r *)s;n(r/)Ge/< +Sn 4  r  J  =  - Q>s(r  S i n (Y *) s  Sin  CC)  Cos(<)s.'n(r/)s c:r*) lrt  SifxCci)  Sf  (V)  -  Sn<V)  where the order o f rotations i s , again, I t can be e a s i l y shown that  c7>sO/)  + S.nCrv*)  a tV)  CosCr *)  G>s(r£)  A  5  , rg , r *  i s an orthonormal transformation  since i t has the f o l l o w i n g properties  [<tf=[+i]  T  and  where s u p e r s c r i p t T denotes the transpose of the matrix and [ i j i s the i d e n t i t y matrix.  Thus  An a r b i t r a r y vector v* i n the  frame with components  f  -vr =•  Wl  11 have components i n the  to = H" v-  ) frame given by UJ i n Equation 3.7: (3.7)  where  to =  to, CO:  are the components of u) , which are the components o f AT measured i n the (^)  frame. At node k  of element J K we have a transformation matrix H~  which  was introduced i n Chapter II as a matrix of r o t a t i o n s r e l a t i n g the base vectors of the frame. *j  frame and the base vectors of an a r b i t r a r y  The V-  w i l l be the transformation of the  X  (y^)  vectors i n t o the  vectors by Equation 3.8: k  A  *d The ( ^o* ^  f-  = 4~  V  ^ T  K  A  K  5  (3.8)  transformation i s composed o f large angle r o t a t i o n s 1  1  CosCr*)  1  t  n  a  s p e c i f i c order and i s s i m i l a r i n form to  t  Cosfr, *) 2  -CosCn*)srnCr,|)  Sin  -r S»o<-n<>*) Q>sCnz)  (r*)  The I n i t i a l Coordinate Frame In Equations 3.6 and 3.8 we defined the r e l a t i o n s h i p s between the base vectors of the node global coordinate systems and the a r b i t r a r y coordinate  25  systems  (*j )  and  J  ively.  (^ ) k  at the j  and k  I f we give the base vectors  X  k  ,  .  .  ^  nodes of member JK  and X  x.  respect-  the f o l l o w i n g values:  (3.9)  \ '  where  i)  T  € = Co —  o  6  Al  then the base vectors A K  ^  ^  w i l l - b e the rows o f  j  w i l l be the rows of  ^  tc  and the base vectors  -  M- , i . e . :  l C\  (3.10)  K  *  ( v^G,.)  if'(1,2)  4-Hc B))  T  J  where and 2  I n i t i a l l y , before any deformations have been a p p l i e d to the s t r u c t u r e , the elements w i l l be s t r a i g h t and the (<j J ) aligned i n space. initial  The transformation matrices  frames w i l l be  and 4.^  transformation matrix tf° which r e l a t e s the i n i t i a l  ate frame ^  1 base vectors and the  c a l c u l a t e the i n i t i a l initial  and ( y ^ )  X  or  angles, defined as  geometry of each member.  X  ( r^"  reduce to an nodal c o o r d i n -  base vectors.  We can  P ° r j ) ' , knowing the s  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 i e n t a t i o n of the member cross s e c t i o n i n 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 i n r e l a t i o n to the node global coordinate  26  FIGURE 3.5  INITIAL COORDINATE FRAME  27  system base vectors. The  The  *j  ( t j ) frame has base vectors  ^  0  A o  ,  A ©  v =  ( ^>  A  A  c  X  T  N  ^  where:  )  (- ) 3  11  axis i s defined as the minor p r i n c i p a l axis f o r f l e x u r e of the  undeformed element cross s e c t i o n . p r i n c i p a l axis and  c^  The  axis i s defined as the major  i s the c e n t r o i d a l axis of the element.  Now, i f we w r i t e : XJ  =  I f  ^  J  (3.12)  with (3.13)  y  s „  and  then the components of  3 3  *-J° and cj _  )  , tj* ,  , e t c . , are e a s i l y d e t e r . A O  X  mined from the coordinates of the nodes.  Ao  ( y  i s dependent upon  y  and  A o  \j  and y i e l d s no new information)  ~~  A O  components of tion A  ^ ,  =  and  *j  But, from Equation 3.10, we know the  &  o  i n terms of  , r° s  , and  , v i z . Equa-  3.14: /  0  y  A  SlnrC)Sn(r )Cos(r °) 5  E  ^CosCOfcsG-j)  +  Cos  C^)  0  f c  Sin < C )  -Cosrr *)s'.ofrs )ro Cr ')\ 8  +  +S,n  S  ( V ) S Cr ") n  T  6  /  fc  (3.14) (  S,nCr ") s  -5inCr/)  Costs')  CosCr °) A  G>srr ^)  T  s  Thus we w r i t e the f o l l o w i n g equations which w i l l be solved f o r the three angles,  r;° , r 5  , and r°  :  28 <jj  =  Cos  1\2  -  Co^Cr °)  Cr °)  C o s (r °)  s  k  H  Sin  ai3° =  (r °) H  Sin  Cr<S)  +  Sm  Cr*J)  Sir,  Cr*>') -  Cos  Cn,')  5i"nCVs*)  tbs  C^)  SV n Crs") G>-5>Cr*r) (3.15)  j*3i -  S.'n  Cr °)  33/ =  - 5 i r ) ( ^ ' )  5  Cis  Cr**)  C o s ('  cbs  r ') 6  Cr ") s  Solving Equation 3 . 1 5 y i e l d s 3 . 1 6 : y »°  Sin(r )= D  3  s  Si'n Cr*°)  ~r~r;  "  provided  Co 5 ( r e . " )  provided  C*osCV  7*  O  Cos C^s /  (Tos Cn, ) = 0  0 s  )£o (3.16)  o  Ca)  Sio(r  >  b  Cos^r  Sil  c 5  )  4- ^ 2  Cos  provided Cos ( b )  o  "i  3 3  and  Vo  3»  C V )  Cb s ( r * ) s  Cvt)= .\ C rfe, ) —  r  Sin  ^ 3 ° ro&Crg') + ^ 3 - 3 ^ 1 , v ~ Cos <- S / and CosG^J/o  ^  r  provided  If  ^^|^o  ^ 3 3 = o choose Equation ( b ) , but i f  u  ) 3 2  also = 0  }  then use the f o l l o w i n g procedure. There are two s i n g u l a r points f o r which the above equations are not v a l i d and that i s when 433°- o  4 ^ =1" '  and the  Cos  ; i t f o l l o w s then that O  .  *i  0  and  The Oj") frame would have i t s  «-} ° 3  A "  A.  axis along the  ^*>z° ^  d i r e c t i o n as shown i n Figure 3 . 6 .  The ^  {  base vector  29  FIGURE 3.6 '  SINGULARITY OF INITIAL COORDINATE FRAME  w i l l l i e i n the  X  ^ 3 plane and thus the  ?  component w i l l be zero.  Therefore:  .  Cos  ( </)  Sin  ( <V)  r  =  0  ^ C O - O Sin  ( 3  (r ')--  s»'n ^ 4.  =  .  i 7 )  XI  s  n») y =  12.  The Incremental Structure Deformations  Sr  The r e l a t i o n s h i p between the intermediate incremental deformations Sf  and the s t r u c t u r e incremental nodal degrees of freedom Sr  using equations 3.2 and 3.3. crements i n the r o t a t i o n s  r*  I t w i l l be r e c a l l e d that £ ~ a r i s i n g from increments  s t r u c t u r e degrees of freedom P . r o t a t i o n s about the I the X axis while 2  axis;  x| Sf  s  are r o t a t i o n s about the At node k  Thus  , S~  , and  H  but  and § r j i s about the X axis; 5  z  TT"  axis while  Sr  t  are the i n -  r r  <Sr  r  <§r  i n the H  c  fc  i s about the  to S i -  are a l l  are r o t a t i o n s about v_ and and a>r  s i m i l a r observations can be made.  Therefore the transformation from Sr  follows  i s given by  J  axis.  where:  o  T o  o  o  o  o p  o  o  I  o  o  o  o  o  »  — 2  o  1  O  C o s C *"icT)  O  I -  o =  O  I  o  o  o  I  o  o  o  I  O  o  O  o  o  o  o  o  o  S o Cn *) 6  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.  D e f i n i t i o n of the Element Degrees of Freedom S The displacements of the nodes of an element with reference to the  l o c a l coordinate system, the ' l o c a l degrees of freedom' 5 of the s t r u c t u r e degrees of freedom r  are functions  defined previously i n Chapter  III.  For the t a n g e n t i a l reference frame there are s i x degrees of freedom per element, a l l defined at one node as shown i n Figure 4.1.  The S  are kept  small with respect to the length of the element by s u i t a b l y i n c r e a s i n g the number of elements required to model a given s t r u c t u r e under a large r o t a t i o n or t r a n s l a t i o n . The t r a n s l a t i o n s  I S,  S  S )  2  3  act along the tangential reference  A-fc  frame base vectors  y  vT  r e s p e c t i v e l y and the r o t a t i o n s  are r o t a t i o n s about these same base vectors. are only approximately about the  y  C  5  x.)  S  A c t u a l l y the ( S^. S  s  2.  frame coincides with the  t  base vectors, as w i l l be shown i n the  f o l l o w i n g s e c t i o n s , but the approximations w i l l be found acceptable. j  S )^  The  ) frame defined i n Equation 3 . 5 , i . e .  q*= ^ V  (4.D  Derivation of the T r a n s l a t i o n a l S Figure 4.1 shows the element in the tangential reference frame before  and a f t e r the f i n i t e displacements r , i n which the element deformations are exaggerated f o r 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 . ordinates of J  in  (X )  Before the displacements, the nodal co-  are C o  J  o  o")  and of K. are given by the  components of A where:  (  A -  3  t  a2  T  a,  and  )  A f t e r the displacements j  has  (x^ )  In the  (4.2)  where  r  i s defined i n Chapter I I ,  and node k  frame coordinates of  (.y^j  f c  frame, the coordinates of j  node  has coordinates  and k  before d i s p l a c e -  o r e s p e c t i v e l y . I f the displacement were a r i g i d o body motion, these coordinates would remain unchanged; node K would not ments are  and  move i n the  ) frame.  K (  The d i f f e r e n c e between the actual motion of node  and t h a t corresponding to a r i g i d body motion i s the deformation vector Si S  • Therefore, the deformations  2  components of displacement of k  with respect to J  are given by Equation 4.3 and c a l l e d S a,  r  S  2  S  i n the  ( »,o  7  <4-i (  L  a, /  {  4 are the  )  3  C-* ) J  T  , seen as frame,  :  -  whe re  (S  j  (4.3)  1,2)  0,3)  x^ -components of the length of the e l e -  L At  ment which i s o r i e n t e d along the *j  l o c a l coordinate system, the vector transformed by Equation 4.4.  direction. 5  r  (5|  In the Cd S  z  S,J  ) frame, the i s the  S  5  =  R  7-^  (4.4)  S  has been derived i n Chapter III  given as i n Chapter  and i s composed o f r o t a t i o n s  II:  r •* r* 3.  Derivation of the Rotational Components of  £  We have the, f o l l o w i n g transformations from Equations 3.5 and 3.8:  ^  *  I  -  We can r e l a t e base vectors at nodes j  and K. from Equations 3.5 and 3.8  A j  by noting t h a t the base vectors X  *1  •iT  and  A k  and  x  are i d e n t i c a l .  Thus:  i  (4.5)  3 =  We define the transformation vectors rotations  'j  and ( e,  . ^3)  |[*f-K]  which r e l a t e s ba se  This transformation m a t r i x , composed of three J will  n a v e  the f o l l o w i n g form:  37  Cos C©0 S.n ( e ) 3  + C o s C& ) Cos Co )  Sin(O-i)  2  Sin  3  CzsC&a)  t  -  5>n(A) s r n ( o ^  - Si"nf&,)SinCefe)Sn(5%) +  (4.6)  Sro^CosCOj)  C©2)  Although t h i s transformation matrix implies t h a t &\ i s about the a x i s , e>  i  s  2  about some intermediate axis between y  i s about the M  K  2  and  «j  , and  a x i s , the r o t a t i o n s are small and the axes of r o t a t i o n -t  3  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 u n i t y , 3f to:  reduces  l.o  V  l.o  =  ©«  (4.7)  l.o  &2  We make the assumption that (  s  ( &  I t i s f e l t t h a t these approximations  x  e>  e>  z  3  )  respectively.  5  s  to  )  w i l l not a f f e c t the required d e r i v a t i v e s of s Equations 4.5 and 4.7 we obtain Equations  are the r o t a t i o n s  with respect to jr  .  From  4.8:  (4.8) S  5  =  In summary, f o r the t a n g e n t i a l coordinate system, we have:  5 ;  >.r  (a.+r,-^ q - ^ E , ! ) ^  ^ -tr -^)c J( 2  8  f  2 )  2 ) + ^ +r -r )v/.Y2 3) 3  q  3  J  (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 l o c a l degrees of freedom S .  This number r e s -  u l t s a f t e r s u b t r a c t i n g the s i x r i g i d body motions from the twelve global degrees of freedom f o r the element.  The s i x r e t a i n e d degrees of freedom  c o n s i s t of an e l o n g a t i o n , four bending r o t a t i o n s , and a t o r s i o n a l r o t a t i o n ; two bending r o t a t i o n s are at node j  and the other four are at nodek  Figure 5.1 shows the member before and a f t e r f i n i t e displacements r where the member displacements have been exaggerated f o r c l a r i t y . f i g u r e shows the base vectors of the o*  the W  n*yto  4  ( Vi©* n *  ^2*)  T  t 0  t h e  (X^)  body motions. base vector M  -fc  j and through  have been defined by  y  Equation 2.5 eliminates t r a n s l a t i o n a l r i g i d  Equation 2.6 sets A  The  system at node k by Equations 3.5 and  3.8 r e s p e c t i v e l y . The secant reference frame base vectors Equations 2.5, 2.6 and 2.7.  ,  system transformed through  system at node  ^J^j  .  A 5  y  orthogonal to  ^  y  5  and the tangent A  _  ( i . e . orthogonal to the plane containing y  5  and the  minor p r i n c i p a l axis of the cross s e c t i o n i n i t s deformed p o s i t i o n at j Equation 2.7 completes the right-handed orthogonal  triad.  We w i l l also need a l o c a l coordinate system at node k A  the  ^  aligned with  S  (  base vector i n order to define the element degrees of freedom  adequately. ^  ).  We c a l l t h i s coordinate frame t h e C ^ ) f r a m e S  which are given i n Equations 5.1, 5.2 and 5.3.  with base vectors  FIGURE 5.1  DEFORMATIONS S  FOR SECANT FRAME  41  s  A  s  A  (5.1)  1. A  j  As 2 ,  A-5  k  x g»  3  A  A  „.  ^3  ^ 2.  2.  A  S  D e r i v a t i o n of s  X  „  2 ,  a  (5.2)  S  s  A  s  X z|  (5.3)  as a Function of  r  We begin by noting that since the element deformations are required to be s m a l l , then the  ) coordinate system and  systems w i l l be almost c o i n c i d e n t .  the(.HJcoordinate  Likewise, the (cj J c o o r d i n a t e frame K  and ( 2 ) c o o r d i n a t e frame are almost c o i n c i d e n t and t h e r e f o r e the r o t a t i o n s  transformations r e l a t i n g t h e i r base vectors w i l l be composed of small angles.  We can say, then, that the f o l l o w i n g s i m p l i f i c a t i o n s to the  s c a l a r products hold: ^  Q  . Cj S  J  . 0 .  =  s  *p  =  =  C O S Cd )= t  l.o since  cos ( ? o ° r ^ ) cos (9o°  ±«>) c*^  ,  <X, i s considered small  = £« Z =  +  * 3  are also small angles  42  The elongation of the element the l i n e j o i n i n g node J  3  t  i s determined by the extension of  to node K :  where i_ defines the length of the element, L = J (q, a  a*  -+ a.* -t a*  and  are the i n i t i a l components of the length vector.  z  From Figure 5.2 i t can be seen t h a t : A 5  s = -  y  2  s =  i  3  4  •  *  (5.5)  (5.6)  ^2  At node K we obtain s i m i l a r expressions  and Sfc  f o r the  ro-  tations: A. K  s = -  /V S  (5.7)  2 , - 3  s  ^ 5  A K  ^  • 5*  (5.8)  The f i n a l degree of freedom, the angle of t w i s t of the element, can be AS A approximated by the s c a l a r product of vectors * j and . 5  AS  S  =  A  H  •  2  5  (5.9)  In summary, f o r the secant element coordinate system, we have the f o l l o w i n g degrees of freedom  :  A S  s  '  2  A  J  ^ A S  (5.10)  43  A. s  s _3  A  W2 A  A'  ^  s  S  2  (5.10 cont'd)  FIGURE 5.2  ROTATIONAL DEGREES OF FREEDOM  S  2  AND  CHAPTER VI FORCE-DISPLACEMENT RELATIONSHIPS 1.  Element S t i f f n e s s The element s t i f f n e s s matrices were taken from Nathan (11).  They  were developed f o r a doubly-symmetric t h i n p r i s m a t i c element w i t h eleven degrees o f 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 t o r s i o n a l  and two r a t e o f change of t w i s t angle.  Nathan's s t i f f n e s s matrices were  reduced from eleven degrees o f freedom to s i x by e l i m i n a t i n g those which are r e s t r a i n e d i n the present case.  The matrices were a l s o transformed  to conform to the present coordinate system d e f i n i t i o n s .  We present here  only the s t i f f n e s s matrices f o r the c a n t i l e v e r element since the secant element s t i f f n e s s can be derived from the given m a t r i c e s . In reference (11) the element forces £>c are f i r s t determined as nonl i n e a r functions of the displacements  These  equations are then l i n e a r i z e d by an expansion o f the element forces i n a Taylor s e r i e s about some instantaneous element displacements S ° . making the displacement increments  AS  By  s u i t a b l y s m a l l , terms past the  f i r s t can be discarded from Equation 6.1.  (6.1)  (6.2)  or where s;  The matrix k matrix  k  can be w r i t t e n  j? =  k „  !5 «  » where the s t i f f n e s s  i s independent of S ° and the matrix k\  Q  is linear in  s"  .  Use o f the matrices given below requires that elongations, shear d i s p l a c e ments, and r o t a t i o n s remain s m a l l .  Displacements must be small compared  to the element dimensions and the r o t a t i o n s must be much l e s s than u n i t y . Nathan's transformed element s t i f f n e s s e s f o l l o w : AS L IEEI, L*  O O O  o  O  o  O  where ^  I  3  5y  O  (6.3)  GJ  o  L-  m  L-  Co &I2  O  O  O  2  ^-6X3  O  i s the moment of i n e r t i a about the c a n t i l e v e r coordinate system  axis and  I  2  .t  i s about the <^  axis.  J  z  »2  T  =  t o r s i o n a l 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  is greater than  I,  47  o  o  O  o  O  o  O  o  o  o  tOL.  m  ST*  (6.4)  a o  ^ £  .2.  O  lOL.  o  o  2,  o S i - *  k,(s )= 6  2  o  O  o  o  o  o  o  o  (6.5)  _ Co G J ^  o  o IT  5  o  O  IO L  £r  - 2 /  2  5  o o  o o  O o  O  o o o  &£-  o  o  O  o  o  io u  (6.6)  o  5  o  Z_  z  o  o  o O  o A S  o  f)£  2 •5  O  O  (6.7)  -21  O  o  o  O  o  o  O  o  l_2  5  5  O  |_  O  O  o  -  O I D u  o  It, ( s j ) *  O  o D Z  A£ '5  2.  O  (6.8)  zt e i i  O  o  5  O  o o  O  5  L.  O  O O  O  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  The element incremental displacements follow as:  ( 6  9 )  =  <3ij Cr)  b  j  k  Cr)  S r  (6.10)  K  or expressed i n matrix n o t a t i o n : (6.11) The transformation matrix <X i s a f u n c t i o n of the global coordinates V .  By the p r i n c i p l e of v i r t u a l work, the external work of the R  moving through the corresponding v i r t u a l displacements S r to the i n t e r n a l work of the element forces $ i b l e v i r t u a l displacements S~s .  forces  must be equal  moving through the compat-  Thus, (6.12)  Using the c o m p a t i b i l i t y Equation 6.11, we have (6.13) I f 6.13 i s true f o r an a r b i t r a r y v i r t u a l displacement § V  >  w e  c o n  ~  elude t h a t : 5  = «  T  £  We expand the R c  (6.14) i n a Taylor s e r i e s 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 . of r  The domain  i s s u i t a b l y r e s t r i c t e d 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 o b t a i n :  (6.16)  or i n matrix form:  A  where  R  =  K  K  A  (6.17)  r-  =  — r 1  The j  =  th column of K  zf ^j? 55  c  i s given by  aj Tjrj  +  2>f  b r  2rj  bf  ^rj  (6.18)  where j ^ j is a column vector of zeros with a one in the j and,  y *  C  c  -  b d  =  b  ri  Thus the global stiffness K where  th place,  =  K  +  K  ~  and the columns of  is given by (6.19)  2  ^  K  =  2  ,  [ g  The element stiffnesses  te  0  ( and  ^ k»  +" «  ^  )  b :  iJ  ?  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 of Equation 6.17.  Then, given the global displacements  the element f o r c e s ^  r and f  K  and  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 e n t i a t i o n of Equations 4.9 and 5.10 f o r the secant or tangent element s t i f f n e s s systems is a mechanical process;  however, the  r e s u l t s are extremely lengthy and have been included d i r e c t l y i n the computer program and there seems l i t t l e point i n reproducing these r e s u l t s here. 4.  Incremental S o l u t i o n Technique When the loads are s i n g l e - v a l u e d functions of displacements, an i n -  cremental displacement method can be used.  This method i s most useful i n  f o l l o w i n g the complete load displacement h i s t o r y of snap-through  buckling  problems where the same t o t a l load occurs at d i f f e r e n t displacements. Otherwise an incremental load method may be used.  No attempt i s made here  to analyze s t r u c t u r e s whose load-displacement curves involve b i f u r c a t i o n s of the buckling paths beyond the b i f u r c a t i o n point. The f o l l o w i n g i l l u s t r a t e s the procedures involved during one increment of the s o l u t i o n technique. At the beginning of an increment, the t o t a l global displacements  £"  at the end of the previous increment are used to c a l c u l a t e the a * ' , fc> and V  matrices defined by Equations 6.10, 3.18 and 6.18.  placements § ness matrix  The element d i s -  at the end of the previous increment are used i n the s t i f f and then the matrix products of Equation 6.19 are c a l -  culated to produce the instantaneous  tangent s t i f f n e s s matrix K of Equa-  t i o n 6.17. The s t i f f n e s s matrix <  i s to be inverted by a band i n v e r s i o n r o u t i n e .  For maximum s o l u t i o n e f f i c i e n c y the band width of the matrix should be a  53  minimum.  Both p o s i t i v e - d e f i n i t e and p o s i t i v e s e m i - d e f i n i t e matrices are  allowed. A one parameter (X) load system i s applied to the s t r u c t u r e and i s assumed to vary l i n e a r l y with displacement during the increment.  For the  incremental load s o l u t i o n method, the load increment vector which i s supplied by the analyst i s used i n equation 6.17 to produce the incremental d e f l e c t i o n vector. load vector /3  However, f o r the displacement increment method, a u n i t  i s applied to produce the d e f l e c t i o n vector  $ (6.20)  The load increment deflections J °  X  and  i s c a l c u l a t e d by l i n e a r l y proportioning the f .  A = where j>  (6.21)  i s the displacement increment s p e c i f i e d by the analyst and ^  i s the displacement at the same degree of freedom as load vector  under the u n i t  /3  Thus we have: * 5 = A N D  A. i r  -  X  @  (6.22)  X ^  (6.23)  These increments of global force and displacement are added to the t o t a l force and displacement vectors  K  then c a l c u l a t e d by Equation 6.11 and S/ST  and r  r e s p e c t i v e l y . Ss  by Equation 6.2 before  proceeding to the next increment of the s o l u t i o n .  is  CHAPTER VII BIFURCATION BUCKLING Many structures e x h i b i t a branching of t h e i r e q u i l i b r i u m paths that i s associated with some mode of s t r u c t u r a l i n s t a b i l i t y .  For most such  s t r u c t u r e s , the determination of t h i s ' c r i t i c a l v a l u e ' involves using the l i n e a r s t i f f n e s s based on the i n i t i a l geometry plus a non l i n e a r c o n t r i bution based on some small perturbation from the i n i t i a l p o s i t i o n . the s t i f f n e s s matrices  Ko  , K|  , and  mination of the b i f u r c a t i o n p o i n t .  Thus  should be useful f o r deter-  The nature of the e q u i l i b r i u m paths  i n the immediate post buckling range i s of i n t e r e s t to the a n a l y s t , but an i n v e s t i g a t i o n of t h i s problem w i l l not be attempted here. Unfortunately, the s o l u t i o n of the eigenvalue problem associated with the determination of the b i f u r c a t i o n point involves the i n v e r s i o n of a f u l l matrix and t h i s puts a c o n s t r a i n t on the s i z e of the problem which can be handled.  But i f only a few of the smaller eigenvalues are r e -  quired and the s t i f f n e s s matrices are w e l l banded, then an i t e r a t i v e a l gorithm i s p o s s i b l e .  The method i s e f f i c i e n t f o r f i n d i n g a few of the  lowest eigenvalues of a p o s i t i v e d e f i n i t e matrix  Ko .  No attempt should  be made to determine a l l of the eigenvalues by t h i s i t e r a t i v e scheme.  A  problem i n convergence does a r i s e when two eigenvalues are the same or almost equal but we assume that this case does not occur often and t h a t , when i t does, i t can be solved by an a l t e r n a t e procedure. The procedure f o r determining the lowest c r i t i c a l value only w i l l be o u t l i n e d here, since the i n t e r e s t of the analyst i s usually confined to t h i s one.  An i t e r a t i v e algorithm has the advantage that the s o l u t i o n 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 c o n f i g u r a t i o n on the s t r u c t u r e to a s i n g l e parameter system where the loads increase uniformly according to a f a c t o r A . The d e f l e c t i o n vector  S_r  corresponding to a normalized force vector  i s c a l c u l a t e d using the l i n e a r s t i f f n e s s Ko  .  ^\  and ] f z  a r e  then b u i l t using t h i s d e f l e c t i o n vector to y i e l d : k  (5,  P  f S r ' )  (Sr )  4-  1  (  We form the t o t a l s t i f f n e s s of Equation 7.2 assuming that l i n e a r up to some c r i t i c a l value of (  Ko  +  ^  P  )  ^  -  , -\rr  J  )  is  •  S R  -  7  ( 7  .2)  The c l a s s i c a l eigenvalue problem r e s u l t s when the right-hand side of Equation 7.2 vanishes. I  T where  K ' ~~ 0  Sr "  On rearrangement we o b t a i n : =  _ K-p  S P  ~  ~  (7.3)  S r - corresponds to the eigenvector and, again,  X<rr  i s 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 i n s e r t e d i n the  right-hand side of Equation 7.3, and the l e f t - h a n d s i d e i s used to solve f o r 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 Sr'"  \r-  and repeat t h i s procedure u n t i l the c r i t i c a l value  • We normalize Xc  r  stabil-  izes according to some convergence c r i t e r i o n or u n t i l the number of i t e r a t i v e cycles exceeds some given number. \ r r  The rate of convergence to  and S r , the eigenvector, and the time required depends on the s i z e  of the problem as w e l l as on the r e l a t i v e d i f f e r e n c e between the f i r s t and second eigenvalues. Derivation of t h i s procedure i s given i n Wilkinson  (16).  CHAPTER VIII NUMERICAL EXAMPLES 1.  Williams'  Toggle  In 1964 Williams (10) e x t e n s i v e l y studied the snap-through phenomena of a shallow planar frame i l l u s t r a t e d in Figure 8.1.  buckling The geo-  metric and e l a s t i c properties of t h i s frame, c a l l e d a ' t o g g l e ' by W i l l i a m s , are the f o l l o w i n g : £1  =  9.27 x 10  AG  =  1.885 x 10  t-  =  26 inch  h  =  0.32 inch  lb. inch  3  6  2  lb.  The large d e f l e c t i o n behaviour of t h i s s t r u c t u r e under a v e r t i c a l applied load on the c e n t e r l i n e i s c h a r a c t e r i z e d by a softening region followed by a hardening region. Ebner and U c c i f e r r o (12) i n a recent paper reviewed and compared f i v e p a r t i c u l a r methods and t h e i r a p p l i c a t i o n s to geometric n o n - l i n e a r planar problems.  These f i v e formulations are: 1.  Argyris  (1964), (13)  2.  Martin (1965), (6)  3.  Jennings  4.  M a l l e t and Marcal (1968), (3)  5.  Powell (1969), (14).  (1968), (4)  In a l l of these methods the n o n - l i n e a r force displacement r e l a t i o n s h i p s  13  .  II  E l  =  q.27  A E  =  i. ess  W s e  FIGURE 8.1  ot  4.  x  t o  it>  3  x to  fc  .32  i n c h  2  lb inch  WILLIAMS' TOGGLE WITH CENTRAL LOAD  59  have been formulated i n terms of s t i f f n e s s matrices.  The s o l u t i o n t e c h -  niques incorporated were: 1.  direct solution,  2.  load incrementation,  3.  displacement incrementation.  W i l l i a m s ' toggle was used by Ebner and U c c i f e r r o f o r t e s t i n g the f i v e methods.  The s o l u t i o n s f o r the toggle under a constant c e n t e r l i n e v e r t i -  cal load of eighty pounds are given i n Table 1.  In view of symmetry, only  h a l f of the s t r u c t u r e was modelled and only eight elements were used. The r e s u l t s , as can be expected, v a r i e d but Jennings'  procedures  give  e x c e l l e n t agreement with W i l l i a m s ' s o l u t i o n f o r a l l three s o l u t i o n t e c h niques.  The ' s e c a n t ' and ' t a n g e n t i a l ' element s o l u t i o n s f o r various  dis-  placement increment s i z e s and numbers of elements (modelling h a l f of the toggle) are shown i n Table 2 and Table 3.  For t h i s example, the ' s e c a n t '  element s o l u t i o n s converged monotonically to a lower c e n t e r l i n e d e f l e c t i o n than the ' t a n g e n t i a l ' element s o l u t i o n with the number of elements constant at ten per h a l f span.  However, the d i f f e r e n c e of 0.002 inches i n  0.6159 or 0.3% i s n e g l i g i b l e .  When the displacement increment s i z e was  kept constant at 0.01 inch and the number of elements used to model the s t r u c t u r e were v a r i e d , the ' s e c a n t ' elements performed extremely well even when using only one element.  The ' t a n g e n t i a l ' element s o l u t i o n s showed a  more rapid convergence although no s o l u t i o n s are a v a i l a b l e f o r using one or two elements per h a l f span.  This element behaves poorly when the e l e -  ment displacements become too l a r g e .  It i s probable that the ' t a n g e n t i a l '  element gives a poor r e f l e c t i o n of the a x i a l deformation when element  60  TABLE 1 SOLUTIONS FOR WILLIAMS' TOGGLE (Under 80 l b . load and 8 elements per h a l f )  Formulation  Center!ine Deflection (inches)  Jennings  Direct  0.611  Powell  Di r e c t  0.611  Mallet-Marcal  Direct  0.600  Williams  Exact  0.611  Jennings  1 l b . load increments  0.639  Powell  1 l b . load increments  0.639  Martin  1 l b . load increments  0.640  Jennings  0.007 inch, d i s p l . increments  0.616  Martin  0.007 i n c h , d i s p l . increments  0.621  0.007 i n c h , d i s p l . increments  .6203  0.007 i n c h , d i s p l . increments  .6161  'Tangential 'Secant'  1  element  element  TABLE 2 WILLIAMS' TOGGLE O.Ol inch displacement increments and 80 l b . load  ' T a n g e n t i a l ' Elements Centerline D e f l e c t i o n (inches)  'Secant' Elements Centerline Deflection (inches)  per h a l f  0.6312  2 elements per h a l f  0.6308  1 element  5 elements per h a l f  0.6550  0.6191  10 elements per h a l f  0.6200  0.6179  20 elements per h a l f  0.6180  0.6179  TABLE 3 WILLIAM'S TOGGLE 10 elements per h a l f and 80 l b . load  ' T a n g e n t i a l ' Elements Centerline Deflection (inches)  'Secant' Elements Centerline Deflection (inches)  0.007 inches per increment  0.6178  0.6159  0.010 inches per increment  0.6200  0.6179  0.020 inches per increment  0.6262  0.6246  0.040 inches per increment  0.6406  0.6387  62  d e f l e c t i o n s are l a r g e .  This would be important i n the s t r u c t u r e under  discussion. Jennings  (4) derived two element s t i f f n e s s e s i n c l u d i n g varying  degrees of geometric n o n - l i n e a r i t y , and was able to obtain good agreement with W i l l i a m s ' r e s u 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 s o p h i s t i c a t e d elements were required than of his s o p h i s t i c a t e d n o n - l i n e a r elements to accurately model the load deformation behaviour of the s t r u c t u r e .  Although Jennings does not s p e c i -  f y what displacement increment s i z e s he used, the ' t a n g e n t i a l ' element gave s i m i l a r r e s u l t s to Jennings'  less s o p h i s t i c a t e d element s o l u t i o n s .  The ' s e c a n t ' element s o l u t i o n s appear to agree with the r e s u l t s by Jennings'  b e t t e r element.  In general the ' s e c a n t ' elements performed b e t t e r than the ' t a n g e n t i a l ' elements.  However, the time required f o r the ' s e c a n t ' element s o l u t i o n  was greater than f o r the corresponding  ' t a n g e n t i a l ' element s o l u t i o n .  Also,  the d i f f e r e n c e i n the r e s u l t s between using twenty elements per h a l f span and ten elements per h a l f span i s s m a l l , i n d i c a t i n g that the u l t r a - f i n e n e s s of s o l u t i o n i s not warranted.  The reduced increment s i z e seems more impor-  tant than increasing the number of elements i n the s t r u c t u r e . Some of the load d e f l e c t i o n paths f o r t h i s s t r u c t u r e are shown i n Figure 8.2.  •3  center/she  FIGURE 8.2  .s deflection  .Co  - 7  (inches)  INCREMENTAL DEFLECTION SOLUTIONS FOR WILLIAMS' TOGGLE as  2.  Argyris  Arch  1  Ebner and U c c i f e r r o (12) also used a plane arch, f i r s t tested by A r g y r i s ( 1 ) , to compare the f i v e various formulations c i t e d p r e v i o u s l y . The r e s u l t s o f these formulations are given i n Table 4 along with the 'tangential  1  and ' s e c a n t ' element s o l u t i o n s f o r the same number of e l e -  ments per h a l f span and s i z e of displacement increment. The e l a s t i c and geometric properties of the arch shown i n Figure 8.3 are as f o l l o w s : EA  =  10  7  lb.  E r  =  10  7  l b . inches  rise  =  3.14 inches  radius of curvature  =  400 inches  length  =  100 inches  2  The arch i s pinned a t the ends and the load i s applied v e r t i c a l l y at the centerline.  Since s t r a i g h t l i n e elements are used, at l e a s t ten elements  are required to model the e n t i r e arch as a polygonal arc.  The fewer e l e -  ments used, the poorer the mathematical model of the intended s t r u c t u r e . For t h i s problem, there was l i t t l e d i f f e r e n c e i n the r e s u l t s between the ' s e c a n t ' element and ' t a n g e n t i a l ' element s o l u t i o n s .  Agreement with  the displacement incrementation s o l u t i o n s of M a r t i n , Jennings and Argyris i s good.  D e f i n i t e l y , smaller d e f l e c t i o n increment sizes and more i n c r e -  ments produced ' b e t t e r ' r e s u l t s f o r the arch of a reasonable number of e l e ments than f o r the same s t r u c t u r e composed of many elements but a l a r g e r increment s i z e .  In the snap-through buckling of this a r c h , the true s o l u -  t i o n curve can be expected to f a l l below the upper snap value that the  L  E  A  ET  =  lo  7  -  lO  7  rise. = frfdius L_ =  lb lb  3.\+ o-f  too  incb  2  inch.  curvat  ore.  —  inch,  FIGURE 8.3  ARGYRIS' ARCH  inch.  TABLE 4 ARGYRIS' ARCH ,  Formulation  Number of elements  Increment  Upper Snap  Lower Snap  lb.  lb.  Martin Load Incr.  5/half  25 l b s .  2300  Powell Load Incr.  5/half  25 l b s .  2300  Jennings Load Incr.  5/half  25 l b s .  2300  -  Martin D i s p l . Incr.  5/half  0.07 inch  2344  814  10/half  0.07 inch  2360  828  5/half  0.07 inch  2360  836  10/half  0.07 inch  2358  837  10/half  0.157 inch  2450  800  5/half  0.07 inch  2328  849  . 10/half  0.07 inch  2356  838  10/half  0.157 inch  2460  846  5/half  0.07 inch  2325  824  10/half  0.07 inch  2355  835  10/half  0.157 inch  2460  843  Jennings D i s p l . Incr.  Argyri s  'Secant'  Elements  'Tangential'  Elements  FIGURE 8.4  LOAD DEFLECTION CURVE FOR ARGYRIS' ARCH  'secant'  and ' t a n g e n t i a l ' element s o l u t i o n s p r e d i c t .  Nevertheless the  s o l u t i o n s f o r the arch modelled by f i v e elements per h a l f span p r e d i c t lower upper snap values than the corresponding s o l u t i o n value f o r the ten elements per h a l f span arch.  Quite l i k e l y the polygonal s t r u c t u r e s are  s u f f i c i e n t l y d i f f e r e n t f o r these two cases t o account f o r t h i s anomaly. One advantage possessed by the displacement incrementation s o l u t i o n technique, compared with the load incrementation s o l u t i o n method, i s the a b i l i t y to f o l l o w the complete load d e f l e c t i o n curve of a snap-through buckling problem.  Several such curves f o r the arch f o r various  ment s i z e s and numbers of elements are shown in Figure 3.  Wright's  incre-  8.4.  R e t i c u l a t e d S h e l l Segment  D.T. Wright (13) i n 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 s t r u c t u r e s , gave some t h e o r e t i c a l c a l c u l a t i o n s f o r snap-through Figure 8.5. load P  buckling of a p a r t i c u l a r spherical s h e l l segment shown i n  He considered the behaviour of t h i s segment under a v e r t i c a l  a t node A .  In the d e r i v a t i o n of the load d e f l e c t i o n curve he  assumed that the s h e l l was shallow and that the f o l l o w i n g r a t i o s h e l d : j_  ~  Sin —  L twhere  f a n Jz  J-r  L»- = L  ^  =  L r  the radius of curvature of the dome the length of the element shown i n Figure 8.5.  When he considered the nodes B  to be pinned and ignored members BB which  do not contribute any moment r e s i s t a n c e , s o l u t i o n f o r snap-through was:  he a r r i v e d at an upper bound  buckling of the s h e l l .  His t h e o r e t i c a l s o l u t i o n  P  3fi£  =  where  K  Ch'-  h'  and  h'  f  h' ) z  EI  (h-h')  •  are defined i n Figure 8.5.  D i f f e r e n t i a t i n g Equation 8.1 with respect to  h'  he found that snap  through buckling would only occur when the f o l l o w i n g c o n d i t i o n h e l d : r  6  .zc*3 h  ^  where  (8.2)  = h  , the radius of gyration of the members = the r i s e of the element.  When he considered the case where the j o i n t s B and the members B e>  are allowed to t r a n s l a t e  are e x t e n s i b l e , he a r r i v e d at the lower bound s o l u -  t i o n f o r the l o a d - d e f l e c t i o n curve of the s h e l l given by Equation 8.3.  p_  A _ e K ' Ch - h ' ) 2  3  2  E Jfb-h')  +  ( ) 8 3  The c r i t i c a l values f o r t h i s equation upon d i f f e r e n t i a t i n g with r e s pect to  h ' occurs when: ^  . 185  (8.4)  b  To v e r i f y Equations 8.1 and 8.3 the material properties of the toggle of Figure 8.1 were s u b s t i t u t e d i n t o these equations. The radius of g y r a t i o n , r g , f o r t h i s s t r u c t u r e i s .072 inches.  Thus  both equations 8.2 and 8.4 are s a t i s f i e d since ^  .1052  f o r equation  8.2  ^  .o7-f  f ° equation  8.4  r  Therefore snap through buckling w i l l occur f o r the pinned s t r u c t u r e and f o r the s t r u c t u r e which allows members e>e> to extend.  Buckling out of  plane has been prevented as we are s o l e l y i n t e r e s t e d i n v e r i f y i n g Wright's calculations.  oil members A -  have  -Ofc.2 6  equal  inches  profx?rt*'es  2  I = . 0 0 0 3 0 ^ inches'* £ =• 3.o x i o ps i 7  P A  j  iranslct ion of node 0  FIGURE 8.5  WRIGHT'S SHELL SEGMENT  71  TABLE 5 WRIGHT'S SHELL SEGMENT  Upper Snap  Lower Snap  lb.  lb.  99.9  48.5  12 'Secant' elements 0.01" d e f l e c t i o n increments  96.8  59.4  12 'Tangential ' elements 0.01" d e f l e c t i o n increments  96.1  49.0  75.2  73.4  *  *  71.8  71.6  BB Wright's  Inextensible  formulation  BB Wright's  Extensible  formulation  12 'Secant' elements 0.01" d e f l e c t i o n increments 12 'Tangential  *  1  elements 0.01" d e f l e c t i o n increments  S o l u t i o n does not e x h i b i t maximum or minimum  • • • vVn'ght'5 12  .1  toocjent  formula iol  .2 Vertical  FIGURE 8 . 6  8.1  e l e m e n t s  .3  per  structure,  .<"»  d e f l e c t i o n ot n o d e A  .oi"  .5  incr.  .6  -7  ('inches)  LOAD DEFLECTION CURVE OF WRIGHT'S SEGMENT WITH POINTS B FIXED AGAINST TRANSLATION  Wrights «0  formolQ  12  tangential  IZ  secant  S-3  elements  per  '/<b  structure  .01"  incr.  \HO. elements  per  /<o  {  structure  .01 " incr.  / /  0»  o C +J  loo.  o 0  "40.  2o.  o.  ,1  .2  Vertical FIGURE 8.7  .3  deflection  .5  ct  node  A  .7  finches)  LOAD DEFLECTION CURVE OF WRIGHT'S SEGMENT WITH POINTS  FREE TO TRANSLATE CO  74  By symmetry, we need consider only one s i x t h of the s h e l l  segment.  Twelve elements were used to model t h i s s t r u c t u r e and d e f l e c t i o n i n c r e ments of 0.01 inches were used. was f a i r as shown i n Table 5.  Agreement with the formulas 8.1 and 8.3 The t h e o r e t i c a l l o a d - d e f l e c t i o n curves and  the experimental r e s u l t s are shown i n Figures 8.6 and 8.7.  4.  Three Dimensional  Elbow  In order to compare the accuracy of the ' t a n g e n t i a l ' and ' s e c a n t ' elements, the s t r u c t u r e shown i n Figure 8.8, an elbow f i x e d at both ends was subjected to an applied force increment.  The actual applied forces  were compared with those required to e q u i l i b r a t e c a l c u l a t e d i n t e r n a l forces (Equation 6.14:  R = Q  t  £f  ).  F i r s t l y , the elbow was modelled by sixteen elements and twenty v e r t i cal f o r c e 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  r e s u l t s given by Table 6 compare the known t o t a l applied force with the forces c a l c u l a t e d by equation 6.14 f o r various increment steps.  Also  shown are the d e f l e c t i o n s under the applied force and the vector of c a l c u l a t e d forces at node 9 of the s t r u c t u r e f o r increment step 20.  All  but the applied f o r c e component of t h i s vector should be zero but the magnitude of the e r r o r when compared with the applied force i s s m a l l . Considering the s i z e of the d e f l e c t i o n , 0.9 i n a length of 4, and of the r o t a t i o n , 0.265 radians in only twenty increment steps, the r e s u l t s are considered reasonable.  75  x. (0 8 8  d> ' — ©  6>  <£>  <Z> \<&  <2>  CO i  i  no de^  number  d> element  C o  V'/  number  X, ^7  7  CQ  it.  -£-3  s e c t lon l.o  in*  i.o  in*  s t r o n g  in  u» e a k  I s - /. o A - I. O E <5  = =  A - A  in  iorsion  tending  k sJ  3 ooo.  ksC  THREE-DIMENSIONAL  pJane  plane.  z  lo o o o .  FIGURE 8.8  b e n d i n g  ELBOW  TABLE 6 THREE DIMENSIONAL ELBOW Applied V e r t i c a l Force Increments  Total Applied Force  Increment Step  Force Calculated by 'Tangential Force  1  Elements Deflection  'Secant' Force  Deflection  Elements  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 , , Fx , F* x  2  3  , «\  x t }  Mx  n  i*lxj )  Applied Moment Increments Moment Calculated by  Increment Step  Total Applied Moment  5  300  299.95  10  600  15 20  - ~ -  _  'Tangential' Moment  Elements Rotation  'Secant' Moment  0.0655  299.97  0.0655  599.85  0.1314  599.87  0.1312  900  899.74  0.1974  899.64  0.1977  1200  1199.60  0.2646  1199.15  0.2654  .  f o r node nine  'Tangential'  (- 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)  Elements Rotation  5.  Ring Dome The large r i n g dome shown i n Figure 8.9 represents a p r a c t i c a l s t r u c -  ture.  By symmetry only h a l f of the dome was modelled.  This helped reduce  the band width of the problem immensely and meant that subdividing the s t r u c t u r e i n t o more elements would not s e r i o u s l y increase t h i s band width. The e l a s t i c and geometric properties of the r i n g dome are as f o l l o w s : All  members:  Ij T I  = 144 z  in.  = 720  in.  3  =  48.0  Area  =  4.0 i n .  (strong bending plane)  4  i n . * (weak bending plane) 1  2  Diameter of dome Rise a t Centre  ( t o r s i o n a l moment of i n e r t i a )  4  = =  Radius of Curvature E  =  3000 ksi  C-r  =  1200 ksi  80 f t . 8 ft. =  104 f t .  The v e r t i c a l plane was made the strong bending plane and the e x t e r i o r nodes were f i x e d . The s t r u c t u r e was divided i n t o 108 ' c a n t i l e v e r ' elements r e s u l t i n g in 504 degrees of freedom and a band width of only 42. ments were chosen p r i m a r i l y because of c o s t .  The ' t a n g e n t i a l ' e l e -  The p r i n c i p a l i n t e n t was to  show the load d e f l e c t i o n behaviour of the dome.  For t h i s , the upper r i n g  was subjected to 16 v e r t i c a l d e f l e c t i o n increments of 0.25 f e e t under a u n i form v e r t i c a l load on the upper r i n g .  The load d e f l e c t i o n curve and the  c o n f i g u r a t i o n of a radius of the s t r u c t u r e a f t e r 16 increments are shown in Figure 8.10  and f i g u r e 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 t h i s s t r u c t u r e .  Had another  mode of f a i l u r e occured other than that shown, then the computer program would have at l e a s t i n d i c a t e d t h i s when the s t i f f n e s s matrix went negative d e f i n i t e .  Although there i s no s o l u t i o n with which to compare  these r e s u l t s , they are s e l f - c o n s i s t e n t and there i s every i n d i c a t i o n that the method was able to handle t h i s comparatively large problem successfully.  FIGURE 8.9  RING DOME  80  40. r  Uppcr  FIGURE 8.10  FIGURE 8.11  ring  vertical  deflection  (  f  e  LOAD DEFLECTION CURVE FOR RING DOME  DEFLECTED CONFIGURATION OF A DOME RIB  e  t  )  CHAPTER IX CONCLUSIONS A procedure was developed, using n o n - l i n e a r s t i f f n e s s matrices from Nathan (11), to f o l l o w the load d e f l e c t i o n paths of shallow framed s t r u c tures.  The necessary transformation matrices and geometrical  relation-  ships were formulated f o r extending previous two-dimensional work to three dimensions. An incremental s o l u t i o n technique was used and shown to be p r a c t i c a l . Two moving element coordinate systems and elements, the ' s e c a n t '  and  ' t a n g e n t i a l ' systems, were developed and r e l a t e d to a f i x e d global system. The snap-through buckling paths of plane frames, arches and space frames were s t u d i e d .  The ' s e c a n t ' element was found to be more e f f i c i e n t  f o r some large d e f l e c t i o n problems i n that fewer elements were required to reach a s u i t a b l e s o l u t i o n .  However, the ' t a n g e n t i a l ' s o l u t i o n required  about one h a l f the computer time of the corresponding.'secant' Reducing the increment step s i z e was found more  solution.  e f f e c t i v e than increasing  the number of elements used to model the s t r u c t u r e .  82  BIBLIOGRAPHY 1.  A r g y r i s , J . H . , "Continua and D i s c o n t i n u a " , Proceedings ence on Matrix  Methods in Structural  Mechanics,  of the Confer-  Wright Patterson A i r  Force Base, October 1965, A i r Force F l i g h t Dynamics Laboratory, Dayton, Ohio, AFFDL-TR-66-20.  2.  3.  K o i t e r , W.T., " E l a s t i c S t a b i l i t y and Post-Buckling Behaviour", NonLinear Problems, Edited by R.E. Langer, U n i v e r s i t y of Wisconsin Press, Madison, Wisconsin, 1963. M a l l e t t , R.H. and Marcal, P.V., " F i n i t e Element Analysis of Nonlinear  S t r u c t u r e s " , Journal  of the Structural  D i v i s i o n , Proceedings of  A.S.C.E., V o l . 94, No. ST9, September 1968.  the  4.  Jennings, A., "Frame Analysis i n c l u d i n g Change of Geometry", Journal of the Structural Division, Proceedings of the A.S.C.E., V o l . 94, No. ST3, March 1968.  5.  B r i t v e c , S . J . and C h i l v e r , A . H . , " E l a s t i c Buckling of R i g i d l y - J o i n t e d Braced Frames", Journal  of the Engineering  Mechanics  Division,  Proceedings of the A.S.C.E., V o l . 89, No. EM6, December 1963. 6.  M a r t i n , H.C., "On the Derivation of S t i f f n e s s Matrices f o r the Analysis of Large D e f l e c t i o n and S t a b i l i t y Problems", Proceedings Conference  on Matrix  Methods in Structural  Mechanics,  of  Wright Patterson  A i r Force Base, October 1965, A i r Force F l i g h 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  S t r u c t u r a l Systems", International  V o l . 4, pp. 23-36, 1969. 8.  Journal  of Non-Linear  Roorda, J . , " S t a b i l i t y of Structures with Small Journal  of the Engineering  Mechanics  Mechanics,  Imperfections",  D i v i s i o n , Proceedings of  A.S.C.E., V o l . 91, No. EMI, February 1965.  9.  Roorda, J . , "The Buckling Behaviour of Imperfect S t r u c t u r a l Systems", Journal  1965.  10.  of the Mechanics  of Physics  and Solids,  Vol.  13, pp. 267-280,  W i l l i a m s , F.W., "An Approach to the Nonlinear Behaviour of the Members of a R i g i d l y Jointed Plane Framework with F i n i t e D e f l e c t i o n s " , Quarterly  Journal  1964, pp. 451-469.  11.  the  of Mechanics  and Applied  Mathematics,  Vol.  17,  Nathan, N.D., " F i n i t e Element Formulation of Geometrically Non-Linear Problems," Ph.D. Dissertation, U n i v e r s i t y of Washington, S e a t t l e , Washington, March 1969.  83 12.  Ebner, A.M. and U c c i f e r r o , J . J . , "A Theoretical and Numerical Comparison of E l a s t i c Nonlinear F i n i t e Element Methods", National Symposium on Computerized  Structural  Analysis  and Design,  U n i v e r s i t y , Washington D.C., March 1972.  13.  George Washington  Wright, D.T., "Membrane Forces and Buckling i n R e t i c u l a t e d 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 . ,  V o l . 91, No. ST1, February 1965.  14.  P o w e l l , G.H., "Theory of Nonlinear E l a s t i c S t r u c t u r e s " , Journal of the Structural Division, Proceedings of the A.S.C.E., V o l . 95, No. ST12, pp. 2687-2701, 1969.  15.  A r g y r i s , J . H . , "Recent Advances i n Matrix Methods o f S t r u c t u r a l A n a l y s i s " , Progress  New York, 1964. 16.  in Aeronautical  W i l k i n s o n , J . H . , The Algebraic  Sciences,  Eigenvalue  V o l . 4, MacMillan Co.,  Problem,  Oxford  (1965).  84 T SNAPS 1! C 2 C 3' 4 i 5 C 6 7 8 C ;  10 : 11 12 •13" 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 '30 31 32 33 34 35 36 37 38 39 40 41 .42 ,.43 44 45 46 47 48 49 5C 51 52 53 54 55  C  PROGRAM TO FOLLOW THE LOAD D E F L E C T I O N H I S T O R Y CF FRAMED S T R U C T U R E S ECTH C A N T I L E V E R ANC SECANT T Y P E E L E M E N T S IMPLICIT REAL*0(A-H,0-2) REAL *4 T I T L E C I M E N S I C N ELCCK CNE AND TWO COMMON / e L K 1 / C N , F l , P 2 , S M K G , N P COMMON /B LK 2/ C4 , C5 , S4 , S 5 , C I 0 , C 1 1 , S I 0 , SI 1 , P R E V , A C I MENS ION J O I N T AMOUNTS CI MENS ION X ( 1 5C ) , Y ( 1 50 ) , 2 ( 1 5 0 ), J N ! 1 50 ) ,NC (15C , 6 ) C I M E N S I C N MEMBER AMOUNTS * , JNL (150 ),JNG( 15C) , N T Y P E t 1 5 0 ) , M L N i 1 5 0 ) *,RR(150,12),SFF(150,6),SS(15C,6),JLN(150),S4C.(150),C40(15C) * , 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 ) X O R IG( 15 C ) , Y O R I G ! 1 5 0 ) 4 , Z G R I G ( 1 5 0 ) , T R { 150, 12) C I M E N S I C N L C £ C AMOUNTS *,EV<600) , C R l e C G ) , C C L ! 1 5 0 ) , P C R I T { 1 5 0 ) CI MENS ION OTHERS *,TITLE(20),R(6),AAA115),EX(15),EY(15),EZ(15),EEE(15),GGG(15) ES(6) *,FCCNS(6),TCTAL(25),LDEG<25),NP!12),SMK0(6,6),CEFL!3) * , ECR( 12) ,ER{ 1 2 ) , S M P R E V ( 6 , f c ) , E D S < 6) , E D S F { 6 ) , A ( 6 ,12) , N D L ( 3 ) *,V(12,6,12),SV(12,12),SMK1(6,6),SSK0(12,12),SSK1(12,12),ATC(12,6) *,ATl(12,fc),TFV<25),ELKA{45),BLKB(116),NlNCR(20) SINC0S(6),XYZCH3) * , S K 0 ( 2 5 0 C C ) , S K I (25CCC) E C U I V A L E N C E ( ELKA ! 1 ) , DM ), ( E L K C ( 1 ) , C 4 ) CAT A NI N C R / 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 , 1 5 1 , * 176 , 2 0 1 , 1 C C 1 / t  C C  f  f  C**  1 2 3 4 23 5 6  FCRMAT(20A4) FCR.-AT ( 2 C I 4 ) FORMAT(7I4.3F1C.3) F O R M A T ! 6 F 1 C .3 ) FCRMAT(2F1C.3,2I4) FCRMIAT(lHl) F O R M A T ( / 2 0 X , • S T R U C T U R E N U M B E R ' , 1 4 , » H A S ' , 1 4 , ' J O I N T S A N D ' , 1 4 , ' ME 4 M E E R S ' / 2 C X , • THERE A R E ' , 1 4 , ' INCREMENTS A N D , 1 4 , ' ELEMENT T Y P E S ' / * 2 C X , ' N B I F I S ' , 1 4 , ' N7 I S ' , 1 4 , ' ANC N'8 I S ' , 1 4 , ' AN D N E L MT IS ' , 14 *,/?0X,« (NELMT I S C MEANS C A N T I L E V E R E L E M E N T S , NON ZERO MEANS S E C A * NT E L E M E N T S ) • ) FORMAT!/' J O I N T KG. T E E 6 C E G R E E S C F FREEDOM THE I N I T I A L J O I N T * COORCINATES'/llX,' X Y Z MX MY MZ• , 1 OX,• X« , 9 X , • Y • , 9 X , • Z • * ) F C R N A T 1 6 X , 1 4 , I X , 6 I 4 , 4X ,3 F 1 C . 3 ) FORMAT!/' MEMBER JN L JNG J O I N T IN X 1 - X 3 PLANE MEMBER T Y P E * XP YP ZP') FCRMAT(1X,I4,2X,I4,2X,I4,5X,I4,20X,I4,9X,3D12.4) FORMAT!/' THE CI F F E R E N T ELEMENT T Y F E S : NUMBER AREA I (X1) * I (X 2 ) I!X3 ) E G« ) FCRMAK31X, I4.6C12.4) FORMAT!/ ' C E F L E C T I C N PARAMETER C I V I C E C BY L E N G T H 0 F ' , F 1 0 . 3 , / ' THE * P O R T I O N OF THE E I G E N V E C T O R L S E O AS A D E F L E C T I O N INCREMENT IS',El * 0 . 4 , ' ( T H I S I S USED ONLY IF N B I F .NE . C ) « ) FORM A T I / 5 C X , ' J O I N T C E G R E E OF FREECfjM N U M B E R S ' ) FORMAT(49X, I 4,4X,614) FORMAT(/' T H E TOTAL NUMBER OF C E G R E E S O F FREEDOM = « , 1 4 ) F C MA T( 14 , 6 F £ . C ) 1  7  8 9 1C 11 12 13  14 15 16 17  0  56 57 58 59 60 , 6 1 62 63 64 65 66 , 67 68 69 70 71 72 ,73 74 75 76 77 78 . 79 eo 81 82 83 84 , 85 8 6' 87 88 89 - 9C 91 92 ' 93 94 95 „ 96 97 8 • 99 ICO 1C 1 1C2 1C3 104 1C5 1 C6 1G7 1C8 1C9 110 11 1 1 12 113 114 115 c  18 19 20 21 22 24 25 2t 27 28 29 30 31 32 70 C 701 7C2 703 801 802 704 70 5 706 707 7C6 7C9 804 £05 71 C 71 1 712 806 71 3 714 809 81C 715 1CC1  V  c** C* * c**  r  »».  c c**  FCRMAT(214,3FE.C) FORMAT! ' THE MAXIMUM HALF BAND WIDTH I N C L U D I N G DIAGONAL ' , 1 4 / 85 *• T t T A L NUMBER CF L O C A T I O N S IN S T I F F N E S S N U X N B ' , 1 6 ) F C R M A T ( / 5 4 X , ' I NCREMENT NUMEER ' , 1 5 ) I YY FORM AT ( / ' J 1 X, « MEMBER U A T A « / 5 0 X , ' MN J N L JNG IXX * IZZ E G ARE A *) FCPMAT(50X,I3,2I4,3E12.4,3E11.3) F U R M A T( • ELEMENT TOTAL UI S P L A C E c N T S R', 1 C X , t F 1 4 . 5 / , 4 0 X , 6 F 1 4 .5 ) F O R M A T ( 1 C X , « ELEMENT TOTAL F O R C E S S ' , 6 X , 6F 1 4 . 4 ) F C R M A T C R M 1 C CET = ' , C 1 4 . 5 ) F C R M A T ( 4 1 X , 'CEGREE OF FREEDOM TOTAL FORCE ' / 5 0 ( 4 5 X , 14 7 X , D 1 4 . 5 / ) ) FORMAT('*',/2 I 1 C / 2 0 A 4 ) FCRMAT(6F12.7 ) F C R M A T { • ( T H E S E F C R C E S ARE N O R M A L I Z E D TO A U N I T FORCE VECTOR IN ' , * ' T H E A P P L I E D D E F L E C T I O N INCREMENT M ETHOD) •) F C B M A T ( 2 C X , « A P P L I E D FORCE INCREMENT METHOD U S E D ) F C R M A T ( 2 C X , * A P P L I E D C E F L E C T ICN INCREMENT MET HOC U S E D ' ) FORMA T( • JOINT A P P L I E D FORCES*/2X , 14,2X,6F12.3J FCRMATC JOINT CCF NDEFL DE FL , 2 I 4, 3 ( I 4 , G 1 2. 3)G12 , .3/ ) F C R M A T ( 2 C X , ' C E S E R V E C MEMBERS * ,4X , 2 0 1 4 ) FORMAT!* P C S I T N NUMBERS FCR M E M B E R ' , 1 4 , / I X , 1 2 1 4 ) F C R M A T ( 2 0 X , * E G AR ST1 S T 2 S T 3 DM • , 1 X , 7 G 1 2 . 3 , / 2 O X , • XM YM ZM P I *P2 F3 P 4 ' , 2 X , 7 G 1 2 . 3 ) F C R M A T ! 3 0 X ' ELEMENT SMKO ' , / 6 ( I X , 6 G 1 2 . 3 / ) ) FCRMAT!2CX,' OBSERVED J O I N T S ' , 5 X , 2 0 I 4 ) X F C R M A T ( 3 0 X , • S T R U C T U R E DATA'//' J O I N T D E G R E E S OF FREEDOM * Y Z' ) FOPMAT(7I4,3X,3D12.4) F0RMAT(30X,» MEMBER N U M B E R ' , 1 5 ) F C R M A T ( • ELEMENT DELTA R C I S P L . ', 1 2 F 9 . 5 ) F C R M A T ! 1 0 X t ' ELEMENT DELTA S D I S P L . ' , 7 X , 6 F 1 4 . 5 , / 1 0 X , ' TOTAL S e * 'CISPL.', 15X,6F14.5) F C R M A T ( / 5 0 X , «A M A T R I X ' , / £ ( 11G 1 1 . 3 , G 1 0 . 3 / ) ) F O R M A T ! / / 5 C X » * T W E L V E V M A T R I C E S * , / / 1 2 ( 6 ( L L G 1 1 .3,G10.3/)// ) ) F O P M A T ( 1 C X , * ELEMENT I N C R E M E N T A L F C R C E S * , 2 X , 6 F 1 4 . 4 ) F O R M A T ( 4 0 X , • V * S F F MATR IX ' /, I 2 ! 1 I G 1 1 . 3 , G 1 0 . 3 / ) ) FCRMAT(5CX,' K l STIFFNESS'/6(2CX,6C-13.5/ )) F O R M A K 5 C X , • ELEMENT S T I F F N E S S K C + K 1 ' / 6 ( 2 0 X , 6G1 3 . 5/) ) F C R M A T ! 4 0 X , ' A T * K 0 *A ' , / 1 2 ( 1X, 1 2( IPE I C . 3 ) / ) ) FCRMAT ( 4 0 X , « £ T * K 1 * A ' , /1 2 (1 X , 1 2 ( 1? E 10 . 3 ) / ) ) FORMAT! • SKC« I X , 1 IE 1 1.3) FORMAT(• S K I ' , 1 X , 1 1 E 1 1 . 3 ) FCRMAT!' D R I NCR. ' , I X , 1 0 ( 1 P D 1 2 . 3 ) / ( I C X , 1 C ( 1 PC 1 2 . 3 ) ) ) R E A D ( 5 ,1 ,ENC = 1CCC) T I T L E WR I T E ( 6 , 5 ) WRITE(6,1) TITLE R E A D ( 5 , 2 ) N R S , N J , N M , N I N C R T , N E L M T , I T Y P E , N B I F , N7,N8 ,NWRT1t NWRT2, L O C I W R I T E ! 6 , 6 ) M R S , N J ,NM ,N I N C R T , I T Y P E , N B I F , N 7 , N 8 , N F L M T IFUOCI.NE.C )KRITE(6,31 ) I F I L C D I . E C . C )V>RITE(6,32 ) NRS IS THE S T R L C T U R E NUMEER NJ I S THE NUMBER CF J O I N T S NM I S THE NUMBER CF MEMBERS NINC R T I S THE NUMBER OF INCREMENTS I T Y P E I S THE NUMBER OF D I F F E R E N T ELEMENT T Y P E S N E I F = 0 DO NCT F I N D LOWEST B U C K L I N G LOAD N 01 F NE C C A L C U L A T E LOWEST 3 U C K L I N C LOAD N7 .NE.O OUTPUT CN U N I T 7 N 8 . N E . 0 OUTPUT ON U N I T 8 NELMT = ZERC FCR C A N T I L E V E R ELEMENTS = NGN 7 E C FCR SECANT ELEMENTS w  1  1  R  116 117 118  O * C* * C**  NV<RT1 IS STARTS N V R T 2 IS  119 -120 121 122 123 124 125 126  C** C** 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 L C D I = N O N .0 FOR A P P L I E D F O R C E I N C R E M E N T R E A C IN T E E J O I N T INFORMATION DO 1 0 0 1 = 1, N J R F A C ( 5 , 3 ) J N ( I ) , ( N D ( I , K ) , K = 1 , 6 ) , X { I) , Y ( I CONTINUE fcRITE(6,7) WRIT E ( 6 » 6 ) (JN(l),(ND(I,K),K=l,6),>(I),  127 126 129 130 131 132 .1*3 3  C** C**  ICC  INCREMENT  NL^BER  AT  The  INCREMENT  NUMBER  AT  fcHICH  OLTPLT  CN  UNITS  R E A C IN INITIAL M E M B E R D A T A : J N L J NG J N P IS A JOINT EXISTING IN T H E X 1 - X 3 WRITE(6,9>  THIS  GUTPUT  R(3)=Z(JNP)-Z  136 137 138 139 140 141  R ( 4 ) = X ( J N G ( I ) )-X ( J N L ( I ) ) R ( 5 ) = Y ( J N G ( I ) )-Y (JNL ( I ) ) R ( 6 ) = Z ( J N G ( I ) )-Z ( J N L ( I ) ) CALL C R G S S ( R ( 4 ) , R ( 1 ) ,ER) CALL C R O S S ( E R , R ( 4 } ,R(1 ) ) XF = R ( 1 ) + X ( J N L U ) )  142  YP = R ( 2 ) + Y ( J N L ( I ) )  Y(I),Z(I),I=l,NJ)  (JNL(I ) )  ZP=R(3)+Z(JNL( I ) ) WRITE(6,1C) I , J N L (I ) , J N G ( I ) , J N P N T Y P E ( I ) , X P , Y P , ZP XCRIG(I)= X (JNG ( 1 ))-X(JNL{ I ) )  .146 147 148 149 tl 50 151  YQRIG(I)=Y(JNG( I ) )-Y(JNL(I ) ) ZCR IG( I ) = Z ( J N G ( I ) )-Z ( J N L ( I ) ) XFC=XP-X(JNL( I ) ) Y F C = Y P - Y I J N L (I ) ) ZPO=ZP-Z(JNL( I ) ) CLZ=DSQRT(X P C * * 2 + YPO**2+ZPG**2)  ""152 i'53 154 155  X1 = C S C R T ( X C R I G ( I ) * * 2 + Y C P IG ( I ) * * 2 + Z C P I G l I ) * * 2 X C 1 = X C R I G (I ) / X 1 YO1 = YCR IG( I ) / X I ZC1 = ZCP1G(I ) / X l XC2 = X P 0 / D L Z  (  ,157 158 1-59 160  YC2 = YPt)/DLZ ZC2=ZP0/DLZ T 1= D S G R T ( I . - X C 2 * * 2 ) I F ( C A B S U 1 ) . L T . 1. C - 1 4 ) G O T O 1 4 5  161 -162 163 ^16 4 16 5 166 167  C 4 C ( I )= Z02/T1 S4C ( I ) = - Y C 2 / T 1 C 5 C( I ) = T1 S5C ( I ) = X0 2 06 0 ( I ) = X 0 1 / T 1 I F ( C A B S ( Z 0 2 ) . L T . 1. C - 1 2 ) G C T C 2 2 6 S6CII)=(YCI*TI*T1+YC2*X0 2*X01)/{Z02*T1)  170 171 17 2 173 174  175  GCTCilOl I F ( C A B S ( Y C 2 ) . LT . 1 . C - 1 2 ) G 0 T C 2 2 7 S6CI I ) = - ( Z C 1 * T 1 * T l + Z G 2 * X G l * X G 2 ) / ( T 1 * Y C 2 ) GCTU101  22 7 145  S6C(1)=YG1 GCTG101 04 C I I ) = 1 . D C  S 4 C ( I )=O.CC  STOPS  JNP NTYPE PLANE NECESSARY  143 H44 14 5  226  8  ) ,2{ I )  R ( 2 ) = Y ( J N P )-Y {JNL ( I ) )  135  168 169  AND  METHCC METHOD  CC 1 0 1 1=1,NM REAC(5,2) JNL ( I ) , J N G ( I ) , J N P , N T Y P E ( I ) R l l )=X( J N P ) - M JNL< I) )  •*l5b  7  86 WHICH  C**  -I 3 4  ;  THE  )  TO  DEFINE  X3  AXIS  116 177  178 179 * 280 18 1 182 163 184 •1 8 5 1,8 6 187 188 189 19C 19 1 19 2 1,9 3 194 195 196 197 198 199 2CC , 2C1 202 •'2X3 204 2C5 2Cb i 2 07 208 2C9 ,2,10 "211 '212 3.13 2 14 215 2 17 ••2 18 249 220 22 1 2 22 223 2 24 225 226 227 228 229 230 2.31 222 223 234 235  C 5 C ( I ) = 0.0C S5C I I ) =xi: 2 C6C( I )=-X02 *ZC1 S6C ( I ) = YG1 CCNT INUE 1C1 WR I T E ( 6 , 1 1 ) CC 102 1=1, I T Y P E R E A C ( 5 ,4) * A M I ) , E X ( I ) , E Y ( I ) , E Z ( I ) , E E E ( I ) , G G G ( I ) WR I T E ( 6 , 12 ) I , A A A ( I ) ,EX( 1 ) , E Y ( I ) , E Z ( I ) , E E E ( I ) , G G G ( I ) CCNT INUE 102 V R I T E OUT I N I T I A L CAT A READ ( 5 , 2 3 ) P L C T L , E I G F A C , J C I N T , N D C F WRITE(6,13) PLCTL,EIGFAC REAC I N THE D E F L E C T I C N PARAMETER AND LCAC C O N F I G U R A T I O N CARE SHOULD BE E X E R C I S E D IN CHOOSING THE D E F L E C T I C N PARAMETER c * * N L P T S I S THE NUMBER OF J O I N T S WHERE L O A D S ARE A P P L I E D  20C 20 1  202 1C5 104  137 103  R E A C ( 5 , 2 ) NLPTS WR I T E ( 6 , 1 4 ) IFINBIF.NE.G)NBIF=1 NU=1 DC 104 J = 1 , N J DO 105 1=1,6 IF (N D ( J »I)-1 ) 2 C 0 , 2 0 1 , 2 0 2 NC ( J , I ) = 0 GO TO 1 0 5 NCI J , I ) = N U NU=NU+1 GO TO 10 5 NN=ND(J,I) NC(J,I)= N0(NN,I) CCNTINUE WRITE(6,15) J N ( J ) , ( N D ( J , K ) , K = 1 , 6 ) CCNTINUE NU=NU-1 W P I T E ( 6 , 1 £ ) NU INCEX=0 I F ( J C I N T . N E . C . A N C . N C C F . N E .0 ) I N C E X = N C ( J C I N T , N C C F ) N=l NINCRT=NINCRT + 1 T F = C . CO NGCN E =0 NN=C DO 103 J = l , N L P T S R E A C ( 5 , 1 7 ) J C INT, (FCONS(K ),K=1,6) DC 137 K = l , 6 M = NO(JO INT,K) IF(M.EQ.0 JGCTC137 IF(FCCNS(K).EG.C.CO)GCT0137 N N =NN + 1 LCEG(NN)=M I F ( I N D E X .EC .C ) INCEX = L D E G ( 1 ) TFV(NN)=FCCNS(K) T CT A L ( N N ) = C.CG TF = TF + T F V ( N N )**2 CCNT INUE WPIT E ( 6 , 7 0 C ) J C I N T , ( F C C N S ( K ) ,K= 1 , 6 ) CCNTINUE NLFTS=NN IF(LCDI.NE.C)GCTC144  2 36 237 23E 239' 240 241i 24 2' 243 244 245 "46 2471-,  TF=OSCRT(TF) DC 138 J = 1 , N L P T S TFV(J)=TFV(J)/TF 136 CCNTINUE WR I TE ( 6, 2C ) 144 R E A C ( 5 , 1 8 > J C I NT , N CD F , ( DE F L ( K ) , K = 1, 3 ) C** NCCF = 1 I F THE D E F L E C T IONS A P P L I E C APE T R A N S L A T I O N S C** NCGF = 4 IF THE D E F L E C T I O N S A P P I L E D ARE R O T A T I C N S K=C 224 JNT=ND(JCINT.NDCF+K) IF ( JNT.EC.CIGCTC223 NCL(K+1)=JNT  2 48  K=K+1  249' 2 5C 251 2 52  IF(K.EC.3)GCTC225 GCTC224 NCL(K+1)=1 D E F L ( K + l )=C.C  25 3  !  254' 255 256' 2 57 258' 259' 2'fcC 261 26 2 263 .264 265 . 266 . 267 268 269 27C 271 '272 27 3 2.7 4 275 ,276 27 7 278 279 2'3 0 281 282 28 3 284 285 2-36 287 28 8 289 29 0 2S1 29 2 29 3 29 4 29 5 ;  223  88  K=K+1  IFIK.GE.3)G0T0225 GCTG224 225 CEFLD=DEFL(1)**2+CEFL(2)**2+DEFL<3)**2 KRITE(6,7C1)JCINT,NCGF,(NDL(K),DEFL(K),K=1,3) 2C4 N19=N+19 R E A D ( 5 , 2 ) ( M L N ( J ) ,J=N,N19) C REAC IN THOSE MEMBERS THAT YCU fel SH TC SEE A CC M PL E TE L O A D - D E F L E C T I ON C HISTORY I F ( M L N ( N ) .EG.C> GO TO 20 3 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) N=N+20 GC TO 204 203 CCNTINUE C R E A C IN THOSE J C I N T S THAT YCL WISH TC FOLLOW THROUGH THE LOAD C E F L E C T I C F ISTORY N=l 2C6 M 9 = N + 19 R E A C ( 5 , 2 ) ( J L N ( J ) , J = K«Kl9 ) I F ( J L N ( N ) . E Q .C ) GC TO 2 0 5 I F ( J L N ( N ) . N E . 0 .ANC.N7.NE .0 ) W R I T E ( 7 , 7 0 4 ) { J L N ( J ) , J = N , N 1 9 ) N=N+20 GC TO 206 205 CONTINUE REWIND 1 REWIND 2 REWIND 3 NG0NE=O NB=C PCR = D A B S ( T F V ( 1 ) J/TFV ( 1 ) I F ( N 8 . E Q .C ) W R I T E ( 6 , 2 1 ) C C A L C U L A T E P O S I T I O N NUMBERS NP DC 106 K=1,NN I = JML(K) J=JNG(K) NF( 1)=ND( I , 1 ) NP(2 ) =ND(I,2 ) N F ( 3 )=ND(1,3) NP(4)=ND(1,4) N P ( 5 ) =ND(1,5) N F ( 6 ) =ND ( 1,6 ) NF(7 ) =ND(J, I ) NF(B) = ND(J,2 )  296 2 97 2 .k  NP(9)=ND(J,3> NF ( 10 ) = N C ( J , 4 ) -KF( 11)=NC ( J , 5 )  299  NP ( 1 2 )=!Sli: ( J , (. )  c  3C0 2C1 2C2 3G3 2C4 2C5 .3C6 307 3C8 2C9 3 10 311 212 2 13 '314 315 316 -317 ,318 3 19 '32C 321 222 -2 2 3 ,324 325 "226 3*27 328 329 .220 ' 331 -2 22 ,,2 33 334 22 5 .3 36 ^237 338 -.3 39 3 40 341 ,34 2 24 3 344 .34 5 3 46 347 348 349 250 .3 51 3 52 353 354 355 :  89  L = NTYPE(K ) IF(N7.NE.C)U«ITE(7,7C3) K,NP E = EEE(L ) G=GGG(L) A R = A AA(L ) STl^EX(L) ST2=EY(L) ST3=EZ(L) I F J G . E C . C )G = E / 2 . 6 I F ( S T l . E ' J . C . ) ST1 = ST2 + ST3 XM = X ( J ) - X ( I ) Yf-Y(J)-Y( I ) ZM=Z(J)-Z(I) CM=DSORT(XM**2+YK**2+ZM**2) P1=AR*E/CM P2=12.*E*ST3/(CN**3) P3=12.*E*ST2/(CK**3) P4=G*ST1/CM C A L L E L M S T C ( P 1,P2 ,F3 ,P4 , S M K 0 , C f , N E L K T ) WR ITE ( 1 ) E LKA C C A L C U L A T E C L A S S I C A L S T I F F N E S S FOR THE ELEMENT C C A L C U L A T E MEMBER BANC WIDTH 11=0 J 1=1000 DC 107 L = 1 , 1 2 N=NP(L) I F (N.EQ.O ) GC TC 107 IF(Il.LT.N)ll=N I F ( J l . G T . K ) J1=N 101 CONTINUE NB1=I1-J1 + 1 I F ( N B . L T .NE1 )NB=NB1 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 I F ( N 8 . N E . 0 )WR I T E ( 6 , 8 C 2 ) SM.KO I F J N 8 . E Q . 0 1WRITE ( 6 , 2 2 ) K , I , J , ST 1, S T 2 , S T3 , E , G , AR lCd CCNTINUE NTCTAL=Nl*NB W R I T E ( 6 , 1 9 ) NB,NT CTA L CC 108 K=1,NL DR(K)=0. 106 CCNTINUE IUNIT=2 JLMT =3 NSEVEN=N7 N E I G HT =N 8 I F ( N W R T 2 .EC .C )NWRT 2=N INCRT DC146J=1,6 ER(J)=0. E R ( J + 6)=C . ES(J)=0.DC 146 CCNTINUE DO 110 J = 1 , 1 2 CC 110 K=1,NM TR(K,J)=C.CC RR(K,J)=C.CC IF ( J . G T .6 ) GC TC 110  2 56 357 ! 358 i 359 160 ?6 1 36 2 26 3 3,64 36 5 '•3-66 267' . 36 8 369 2 7C 371 •27 2 2.73 .274 375 .376 3 77 378 3 79 3c.eC 381 382 383 -384 38 5 386 387 288 289 390 ,3 91 '292 593 2'94 295 •2 9 6 ,397 2<8 399 4t0 4C1 4C2 4C3 4 04 4C5 4C6 4 07 '4 08 4 C9 410 4 11 412 4 13 4 14 415  11C  :  229  r  ;  r  111 207  231  2Cc 209 21C 113  114 211 112  J  115  117 212  SFF(K,JJ=C.rC SS(K,J)=C.CC ^ CCNTINUE CC 1 0 9 JK = 1 , M N C R T JKM=JK-1 N8-=0 N7=0 IF{I\WRT1.LE.JKM.ANC.NWPT2.GE.JKM,)N7 = N S E V E N I F (NWRT1 . L E . J K M . A N D . N W P T 2 . G t ' . J K M ) N 8 = NE IGHT I F ( J K M . E 0 . C . dND.N7.E C . 0 ) G 0 T D 2 2 9 W R I T E ( 6 , 2 C ) JKM IT£KP=IUMT IUN I T = JUN I T JUNIT=ITEMP MNC = 0 CC 111 L = 1 , 2 C I F I N INCR (L ) . E O . J K ) N I N C = 1 I F l N I N C . E C . l ) GC TC 207 CCNTINUE IF(N7.NE.C)WRITE(6,7C5) IFCN7.NE.0 1GCTC221 IF(JKK.NE.C./ND.MhC.EC.l)V>RITE(6,705) I F { J K M . N E . C . £ N D . M N C . E C . O . A N D . J L N ( 1) .NE . 0 ) V.R I T E ( 6 ,70 5 ) DC 112 J L = 1 , N J DC 1 1 3 J = l , 3 NN=ND(JL,J) I F ( N N . E Q . C ) GO TO 1 1 3 GC TO ( 2 0 8 , 2 C 9 , 2 1 0 ) , J X( J L ) = X ( J D + C R l N N ) GC TO 112 Y(JL)=Y(JL)+CR(KN) GC TO 1 1 3 Z ( J L ) = Z ( J L ) + CR(NN ) CONTINUE IJT=0 I F ( J K M . E C . C . * N D . N7 . E C . 0 ) CC TO 112 DC 114 L = 1 , 5 C C I F ( J L N ( L ) . E O . C ) G C TO 2 1 1 IF(JLN(L).EC.JL)1JT=1 I F ( I J T . E Q . l ) G C TC 21 1 CONTINUE I F ( I J T . E C . l . C R . N 7 . N E . 0 .OR .NINC . EQ. 1 )WR I T E ( 6 , 7 C 6 ) J N ( J L ) , 4(NC<JL,L),L=1,6),X(JL),Y(JL),Z(JL) CCNTINUE DC 115 I = 1 , N T C T A L SKC(I)=0. SK1(I)=0CCNTINUE REWIND 1 REWIND 2 REWIND 3 DO 116 M L = 1, N M R E A C ( l ) ELK A If N= 0 DC 117 L = 1 , 5 C C I F I N L N I D . E C . O G O TO 2 1 2 IF(MLN(L).FC.NL)IMN=1 I F ( I M N . E C . I )GC TC 2 1 2 CCNTINUE I F { J K M . E C . C . * N C . N 7 . E C .0 ) G C T Q 2 2 8  90  4 16 • 4 17 4 18 419 42 0 421 422 H23 424 425 ^26 427, 428 429 430 431 432 4 33 434 435 436 437 433 439 AfiO 441 442 443 444 445 446 ' 447 448 449 .4 50 45 1 4 52 453 454 455 ,4.5 6 4 57 458 459 4-60 461 462 463 '464 465 466 4 67 .468 469 47C 471 47 2 473 474 475  I F { IMN .EC . L . CR. N7. NE .0 .CR . M NC.NE.O ) WR ITT ( 6 , 7 0 7 ) ML CO 120 J = l , 6 ECR(J ) =0.C ECP(J+6)=C.C ECS(J)=0. ECSF(J)=0. 12C CCNTINUE DC 118 J = l , 1 2 NS=NP(J) I F ( N S . L E . O J G C T C l 18 ECR(J)=DP(NS ) TR(ML,J)=TR(ML,J)+ECR<J) 118 CCNTINUE I F ( J K M . E C . C 1 G C T C 2 15 I F (N7 .NE .C .CR. I l*N.NE. C) WR H E ( 6,7C8> ( EDR( J ) , J = 1 , 12) R E A C U U N I T ) ELKB CC 119 K = l , 6 DC 119 J = 1 , 1 2 ECS(K)=ECS{K)+A(K,J)*EDR(J ) 119 CCNTINUE DC 121 J = ' l , 6 SS(ML,J)=SS(ML»J)+EDS(J) ES(J)=SS(ML»J) 121 CCNTINUE T5=EDR(5) ECR(5)=T5*C4+ECR(6)*S4 E C R ( 6 ) = E C R ( 4 )*S5-T5=>S4*C5 + EDR ( 6 ) * C 4 * C 5 T5=EDR(11) ECR(11)=T5*C1C+E0R(12)*S10 ECR ( 1 2 ) = E D R ( 1 C ) * S 1 1 - T 5 * S 1 0 * C 1 1 + E D R { 1 2 ) * C 1 0 * C 1 1 C C 1 2 2 J = 1,12 PR(ML,J)=RR(f L,J)+EDP(J) ER(J)=RR(NL,J) 122 CCNTINUE IF(N7.NE.C.CR.IMN.NE.0.0R.MNC.NE.C)WRITE(6,7C9) EDS,ES I F (N 7.NE.C .OR . IMN .EC . 1 .OR.NI NC. NE. C) WRITE ( 6 , 2 4 ) ( T R ( M L , J I ) , J I = 1 , 1 2 ) CC 1 2 3 J = l , 6 CC 124 K = l , 6 124 ECSFtJ)=SNPREV(J,K)*EDS(K)+EDSF(J) SFF(ML , J)=SFF{ML,J) + EOSF(J) 122 CCNTINUE IF{IMN.EG.1.CR.N7.NE.C.0R.MNC.NE.C)WRITE(6,25)(SFF{ML,J),J=1,6) I F ( N 7 . N E . O ) W R I T t ( 6 , 7 1 0 ) EDSF 215 XYZC(1)=XCRIG(ML) X Y Z C ( 2) = YDR IG ( M L ) XYZ0(3)=ZORIG(ML ) S I N C O S ( 1 )=S4C (ML ) SINCCS(2)=C4C(ML) S INCOS ( 3 ) = S5C (ML ) SIN C 0 S(4)= C 5 C(ML) SI N C O S ( 5 ) = S 6 C ( M L ) SINC0S(6)=C6C(ML) I F I N E L M T . E Q . C J C A L L A VMATC (\», ER , XYZC , S INCOS ) I F ( N E L M T . N E . C J C A L L A V M A T S ( V , E R , X Y Z C , S I NCOS) CALL T R A N S ( X ) I F ( N 7 . N E . 0 ) WRITE ( 7 , 8 0 4 ) ( ( A ( J , K ) , K = I , 12 ) , J = 1 , 6 ) IF(N8.NE.0)WRITE(8,805)(((V(K,J,I),K=1,12),J=1,6),I=1,12) DC 125 K = l , 1 2 DC 126 J = l , 1 2 SV(J,K)=0.C 226  9 1  ,  4 76  DC 1 2 7 1=1,6 .SV-(J,K) = SV(J,K)  4 77  +V( K , I , J ) * S F F ( K L , I )  92  127 126  CCNTINUE CCNT INUE  125  48 1  CCNTINUE IF(N7.NE.0)WRITE(7,7ll)((SV(J,K),K=l,12),J=l,12)  ^82  CALL  48 3 ^84  DC  128  1=1,6  DC  128  J=l,6  47 8 4  79 C  4 P  'iB5 186 4E7  E L M S T 1 ( P 1 , P 2,DM , E S,S V K 1 , N E L M T )  SMPREV(I,J)=SPKC( 12 8  1 , J ) + SMK1 (I  , J)  CCNT INUE I F ( N 7 . M E . C ) V » R I T E ( 7 , 7 1 2 ) ( ( S f K l ( I J ) , I = l , 6 ) , J = l 6 J 6) IF(N8.NE.C)VsRITE(8,8C6){{SNPREV(I,J),I=l,6),J=l, t  '-:£8 48 9  WRITE(JUNIT) 8LKB CC 129 L=l,12  49 C  DC  49 1 ft 9 2 49 3  129  K=l,12  S S K C ( K , L )= C. SSK1(K,L)=C. I F I L - 6 ) 1 3 0 , 1 3 C , 129  4 94 495 496  13C  97 49 8  129  AT C ( K , L ) = C . AT1(K,L)=0.  4 99  CCNT INUE DC 1 3 1 K=l,12 DQ 131 J=l,6  5CC 5C 1  DC 1 3 2 1=1,6 ATC(K,J)=ATC(K,J)+A(I,K)*SfKO(I,J>  4  f  ATl(K,J)=ATl(K,J)+A(I,K)*SMKl(I,J>  5G2 5C3  132  CCNT INUE  5C4  131  CCNTINUE DC 1 3 3 K=l,12 DQ 133 J=l,12 DC 1 3 4 1=1,6 SSKC(K,J)=SSKC(K,J)+ATO(K,I)*/S(I,J) SSK1(K,J)=SSK1(K,J)+AT1(K,I)*A(I,J)  134  CCNT INUE  133  CCNTINUE IFIN7.NE.O)WRITE(7,713)  ( ( S S K 0 ( K , J ) , K = 1 , 1 2 ) , J = 1 , 12)  IF(N7.NE .0)WRITEI7,714)  ( ( S S K l l K, J ) ,K = 1 , 1 2 ) , J = 1 , 1 2 )  5C5 5C6 5C7 5C8 5C9 .5-10 511 5 12 5-13 5:14  DC  135  L=1,12  DC 1 3 6 M=1,12 M 1 = NP ( L ) N2=NP(M )  5 15 5,1 6 5 17  IF(N1.EQ.C)GC IF ( M 2 . E Q . C ) G C  5 18 5 19  TC TC  135 136  520  IF(M1.CT.M2)GC  52 1 522  K= ( M l - 1 ) * ( N E - 1 ) + r*2 S K C ( K ) = SKC(K )+ S S K O ( L , M ) SK1(K)= SK1(K)-SSK1(L,M)-SV(L,M)  TC  136  523 5 24  136  CCNTINUE  525 526  135 116  CCNTINUE CCNT INUE IF ( N S . ' J E . O )WR I T E ( 8 , 8 0 9 )  ( SK 0 ( J I ) , J I = 1, N T O T A L )  IF(N8.NE.C)WRITE  ( S K 1 ( J I ) , J I = 1,NTCTAL)  527 •5 2 S 529  (8,310)  IF(NBIF.NE.C.ANC.JK.LE.3)C£LL  530 531 532  DC  533 534  SKC(J)=SKO(J)-SK1(J) SK1(J ) =SKC(J )  5 35  C R I T L D ( S K 0 , S K 1 , N U , NB , J K , D R ,  * E V , E IGFAC,PCR ) I F ( J K . G T . l . A N C . N C I F . N E . O . AND . J K . L E . 3 ) G C  141  141  J=1,NTCTAL  CCNTINUE  TC  216  536 5 37 538 539 -'54 0 54 1 5.4 2 543 544 545 .S4 6 5 47 .548 549 .550 55 1 '552 553 5 54 555 556 5 57 -558 559 5^60 561 562 56 3 ^564 565 5.66 567 56 8 569 -570 571 572 ^57 3 5-74 575 ?76 577 " 578 579 58 C 58 1 582 583 584 58 5 586 567 >588 589 59C 591 592 593 594 595  219 139  14C  J  217 218 22C 23C  142 222 143  t  :  J  232 233 221 216  109  10C0  CG  139  J=1,NU  C R I J)=0 . C C N r INUE DC 140 J = I , N L P T S D R ( L 0 E G ( J ) )= T F V ( J ) CCNT INUE CET=1.D-14 NSCALE=1 IF(NGCNE.EQ.C)CALL DFBAND(SKO,CR,NU,NB,1,DET,AET,JET,NSCALE) I F ( N G C i \ E . E C . C . A N C . C E T . L T . l .D-14 )G0 TC 2 1 7 IF(KGCNE) ?2C,22C,218 NGQNE=l G C TO 2 1 9 C A L L DBANDI ( S K I , C R , N U , N B , i , D E T ) IF(JKM.EQ.C.AND.N7.EC.C JG0T023C W R I T E ( 6 , 2 6 )CET CIV=(DR(NDL(1))*CEFL(1)+DR(NDL(2))4CEFL(2)+CR{NCL(3)) * * C E F L ( 3 ) J/CEFLD C I V = 1 ./D IV I F( N B I F . N E . C . A N C . J K . E C . l ) C I V = 1 . 1F(L0DI.NE.0 )CIV=1.0 DO 1 4 2 J=1,NU CR(J)=DR( J)*CI V CCNTINUE  DO  143  J=l,NLPTS  TCTAL(J)=TCT£L(J ) + T F V ( J ) * C I V CCNTINUE I F ( J K M . E l i .0 J G C T C 2 3 2 J K L = JK-NB IF JKR=JKL-1 IF(JKR.EC.C)JKR=1 P C R I T ( J K L ) = TCTAL( D/PCR DCL(JKL)=0ABS(DR( INDEX)/PLOTL )+ DCL(JKR) GCTC233 PCRIT(1)=TCT/?L(1)/PCP DOL(1)=DABS(CR(INDEX )/PLOTL) W R I T E ( 6 , 2 7 ) (LDEC- ( J ) ,TOTAL ( J ) » J = 1 , N L P T S ) IF(N7.NE.C)WRITE (7,715) ( C R ( J I ) , J I = 1 , N U ) GC TO 109 C IV = PCR IF(JK.EC.3)CIV=C. GC TO 222 CCNT INUE WR1TE(0,28) NRS,JKL,TITLE WRITE(C,29)(PCRIT(l),CCL(I),I=l,JKL) GC TO 1 0 0 1 STOP ENC SUBROUTINE T R A N S ( V ) IMPLICIT R EA L*8( A- H , 0- Z ) CCMMON /ELK2/C4,C5,S4,S5,C10,C11,S10,S11,SM,A CI MENS ION SM ( t , 6 ) , A ( 6 ,12 ) , V( 12,6 , 12 ) , G ( 1 2 , 1 2 ) *,B1(3,3),C(3,3),T(6,12),U(6,12),B(6,12),TM(3,3) CC 1 J= 1 , 12 CC 2 K = l ,12 G(K,J)=O.DC IF(K.EQ.J)G(K,J)=1.D0 I F ( K . C T . 6 JGCTC2 B(K,J)=A(K,J) A (K , J ) =0.DC  93  596 597 5,8, 599 6CC  .7 1  601 6C2 6C3 6C4 6C5 •6-C6 6C7  CCNTINUE CCNTINUE G(5,5)=C4 G(6,6)=C4*C5 G ( 5 , 6 ) = S4 G(6,5)=-S4*C5 G(6,4)=S5 G ( l l , l i ) = C l C G(12,12)=C11*C1C G(11,12)=S1C G(12,11)=-S1C*C1 1  ,  6C8 609 610 611 '6 1 2 ,613 :6\A  G ( 1 2 , 1 0 ) = SI 1  201 2CC  DC 1 0 6 N=l,7,6 DO 2 0 0 J=l,3 CC 2 0 1 K=l,3 B1(K,J)=G(K+2+N,J+2+N) CCNTINUE CCNTINUE CC 10 5 M=l,4,3  615  DC  616 617 £18 6.19 62 0 • 621 622 623 624  DC 1 0 0 J=l,3 DC 1 0 1 K=l,3 T M K , J ) = V ( I , K + M - 1 , J+N + 2) C C N T INUE CONTINUE DC 1 5 0 J=l,3 DC 1 5 1 K=l,3 C(K,J)=O.DC DC 1 5 2 1=1,3  t  1C1 ICC  6 25  102  1=1,12  C(K,J)=C(K,J)-»TP(K,L)«ei(L,  626 6 27 628 629 -630 631  152 151 15C  622  1C3  '6'3 3 6<34 635 -4,3 6 627 "638 639  102 1C5 106  104  99  64C  C C N T INUE CONTINUE CCNTINUE DC 1 0 3 J=l,3 DC 1 0 4 K=l,3 V(I,K+M-1,J+N+2)=C(K,J) CCNTINUE CCNTINUE CCNTINUE CCNTINUE NN=1 Nl=4 N!2 = 6 CC 4 0 J=l,3 DC  641  41  642 643 644 645 646  40 5C  647 6 48  41  CONTINUE GCT0(50,60,7C,8C),NN TM2,2)=-S4 TM(2,3)=C4 TM(3,2)=-C4^C5 TM(3 M=4  649  K=l,3  TM(K,J)=C.CC  ,3)=-S4*C5  G0T081  650 £ 5I 6 52 653 654  6C  655  7C  TM(3,1)=C5 TM3,2)=S4*S5 TM(3,3)=-C4*S5 M=5 GCTC81 TM(2,2)=-SLC  656  TM(2,3)=C10  657  TN ( 3 , 2 ) = - C l C * C l  658  TK  659 660  M=1C N1=10  661  N2 = 1 2 GCV06 1 80  TM(3,1)=C11  r64  T M 3,2)  c'6 5  TM ( 3 , 3 ) = - C l C * S - l 1  6.66  =S1C*S11  M=ll  667  CC  90  668  DC  92  669 670  91  DO 9 1 L =N 1,N 2 A ( J , K ) = A ( J , K ) + C ( J , L ) *TK. ( L - N 1 + 1 , K - M  671  92  CCNTINUE  (672  81  J = l ,6  9C  K=N1,N2  DO  93  67A  DC  94  675  V(M,J,K)=  676  A(J,K)=0.CC  G77  94  678  93  679  l+l)  CCNTINUE  67 3  K =N 1, N 2 J=l,6 V(N,J,K)+A(J,K)  CCNTINUE CCNTINUE NN=NN-U  6*80  I F (NN . L E . 4 ) G C T 0 9 9  681  DC  25U  682  DC  251  ff8 3  T(J,K)=V(5,J,K)  684  K=l,12 J= l,6  U(J,K)=V(11,J,K)  685 686  9  (3,3)=-SlC*Cll  66 2 663  1  '  251  CCNTINUE  25C  C C N T INUE  687  DO  3  688  DC  4  689  J= l,6 K=l,12  V(4,J,K)=\M4,J,K)+V(6,J,K)*S5  690  V(5,J,K)=T(J,K)*C4-V(6,J,K)*S4*C5  69 1  V ( 6 , J , K ) = T ( J , K )*S4 + V ( 6 , J , K ) * C 4 « C 5  6 92  V(1C,J,K)=V(1C,J,K)+V(12,J,K)*S11  6,93  V ( l l , J , K ) = U J , K ) * C 1 0 - V ( 1 2 , J , K ) * S 1 0 * C i l  694  V( 1 2 , J , K ) = U  69 5 ,696  CC 5 1 = 1 , 1 2 A(J,K)=A(J,K)+B(J,L)*G(L,K)  6 97  5  698  4  CCNTINUE  6.S9  3  CCNTINUE  7C0  J , K )*S10 + V ( 1 2 , J , K ) * C 1 0 * C 1 1  CONTINUE  RETURN  701  ENC  702  SUBROUTINE  7C3  REALMS  7C4  DC  1  1=1,6  ?,C5  CO  2  J= l,6  7C6  K(J , I  7C7  2  7C8  1  ELMSTC ( A , G, C , C , K , L , N N )  A,B,C,C,L,K(6,6)  )=0.  CCNIINUE CCNTINUE  7C9  IF(NN.NE.GJGCTC1C  7 1C  K(1,1)=A  7 11 712 7 13 7 14 7 15  K(2,2)=C K ( ? ,t» ) = - . 5 * C * L K(3,*)=B M ? , 5 )  = .5*B*L  K ( 4 , s ) = 0  5  7 16  K(5,5)=e*L*L/3.  7 17 7 18 -7 1 9 . 7 20 721 722 723 7 24 725  K (6 , 6 ) = C * L * L / 3 . G O T O 11 K(1,1)=A K( 2 , 2 ) = B * L * L / 3 . K(2,5)=K(2,2)/2. K(3»3)=C*L*L/3. K(3,6)=K(3,3 )/2. K(4,4)=D K(5,5)=8*L*L/3.  72 6 727 7 28 7 29 73C 731 .7 2 2 733 724 735 736 '737  1C  11  3  2  738 739 740 j 741 ' 742 743  IF ( N N . N E . 0 1 G C T C 1 C SMK1(2,6)= - . 1 * P 1 * S ( 1) SMK1(3,3)= 1 .2 * P 1 / D M *S( 1 ) SF.KK 3,5)= .1*P1*S(1) SMK1(4,4)= P2*0M/12.*S(1) SMK1(5,5)= 2 . / 1 5 . *P 1<DM*S(1 )  ;  744 ,745 746 747 748 749 75C 751 752 •353 7,54 755 ,756 757 758 759 7,6'C 761 762 763 764 765 766 767 768 769 77C 77 1 772 773 7 74 77 5  K(6,6)=C*L*l/3 . CC 3 1=1,6 CO 3 J = I , 6 K(J,I)=K(I,J) CCNTINUE RETURN EN C S U B R O U T I N E E L MST 1 (P 1 , P 2 , D M , S , S MK 1, NN ) R E A L * 8 P 1 , P 2 , C f , S C , S { 6 ) , SMK1 (6 , 6 ) DC 2 1=1,6 DC 2 J = l , 6 SPKK I ,J)=C.  10  2CC  3  SMK1(6,6)= 2 . / 1 5 . * P 1 * D f * S (1 ) SMK1(2,2)= 1.2*P 1/DM*S(1> SMKH1,3)= +.1*P1*S(5) +1 . 2 * P 1 / D M * S ( 3 ) SMK1(1,5)= + 2 . / 1 5 . * P 1 * C K * S (5 ) + .1*P1*S(3) S M K K 2 , 4 ) =- 4 . 2 / 1 2 . * P 2 * D M * S ( 5 ) .5*P2*S(3) SMK1(4,6)= +3 . 2 / 1 2 . * P 2 * D M * D M * S (5) + 4.2/12.*P2*DM*S(3) S f K K l ,2)=-.l*Pl*S(6) +1.2*P1/CM*S{2) SMK1(1,6 )=- . 1 * P 1 * S { 2 ) + 2 ./15.*P1*CM*S(6) SMK1(3,4)= 4 . 2 / 1 2 . * P 2 * D M * S ( 6 ),5*P2*S(2) SMK1(4,5)= 3 . 2 / 1 2 . * F2*CM*DF'*S ( 6 )4 . 2 / 12 . * P 2 * C M * S ( 2 ) GCT0200 SMK1(2,2)= 2./15.*P1*DM*S(1) SMK1 ( 3 , 3 ) = SI»K1(2,2) SKK1(2,5)= - P 1* C f * S ( 1 ) / 3 0 . SMK1(3,6)= S t K 1 ( 2* 5) SMK1(4,4)= P2*CM/12.*S(1 ) SMK1(5,5)= 2 . / 1 5 . * F 1 * D f * S (1 ) SMK1(6,6)= 2./15.*P1*DM*S(1 ) SMKU3,4)= . 8 / 1 2 . * P 2 =*D M*D M*S ( 2 ) + P 2 * D M * D M / 1 2 . * S I 5) SMKI(4,6)=3.2/12.*F2*CM*DM*S(5)+P2*CM*CM/12.*S(2) SMK1(1,3)= 2./15.*Fl*DP*S(3)-Pl*CM/30.*S(6) SMK1(1,6)= 2 . / 1 5 . *P 1 * D M * S ( 6 ) - P 1 * 0 M / 3 C . * S I 3) SMK1(2,4)= . 8 / 12 . * P 2 =*DM *DM *S ( 3 ) + P 2 * P M * D M / 1 2 . * S ( 6 ) 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 ) SMK1(1,2)= 2./15.*Pl*DM*S(2)P1*U*/30.*S(5 ) SMK1(1,5)= 2 . / 15 , * P 1 * D M * S ( 5 ) P 1*DM/3C .*S{2) DC 3 1=1,6 DC 3 J = l , 6 S M K U J , I ) =SMK 1 ( I , J ) CCNTINUE RETURN EN C  7 76 777 778 779 -' 76 C 781 -78 2 7 63 . 7 64 ' 785 ^'786 7 87 •> 788 789 t 79C 791 792 .793 . 5794 795 796 ' 797 798 ,799 8C0 eci 8C2 "803 ,* 8 0 4 > 8C5 >8C6 8 07 808 EC9 "'810 , 811 812 "8 13 '8 14 815 1-8 16 8 17 818 6 19 820 82 1 822 823 824 8 2 5; 626 827 828 829 £3 C '83 1 822 8 33 8 34 835  4  1i  J  t  2 3 c  1 10 15  SU E R O U T I N E C R I T L C ( S K 0 , S K 1 , N , M , J K , D R , E V , E I G F A C , P C R ) I M P L I C I T R E A L * 8 ( A - H , C-Z ) T(25CC) D I M E N S I O N C R ( 1 ) , SKO( I ) ,SK I ( 1 ) ,EV( ,IS) IF ( J K . E G . l )GC TC 1 I F < J K . E C . 3 ) G C TC 10 EFS=l.D-e EPSV=l.D-4 IT=100 C A L L D V P C W R ( S KC , S S I , N , M ,M , E V , N , E V A L ,1 I ,E)»S, EPSV , I T , ST, COND ) PCR=EVAL I WRITE(6,4) EVALI F C R M A T ( ' C R I T I C A L LOAD I S ' , E 1 4 . 5 ) CIV=0.0 DO 16 J = 1 , N EVR=DABS(EV(J ) ) DIV=DMAX1 (DIN. ,EVR) DO 2 J = 1 , N EV(J)=EV(J)/CIV WRITE(6,3) ( E V ( J I ) , J I = 1 , N ) FORMiAT(//« E I G E N V E C T O R * / / I C O ( IX , 1 2 G 1 0 . 3 / ) ) DC 5 J = 1,N CR(J)=0 R(J)*(PCR-1.) RETURN DO 15 J = 1 , N CR(J)=EV(J)*EIGFAC RETURN END SUBROUTINE C R Q S S ( A , B , C ) REAL #8 A , E , C D I M E N S I O N A ( l ) , E ( 1) , C ( 1 ) C ( 1 )=A( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) C ( 2 ) = A ( 3 ) * f i ( 1 )-A( 1 ) * B ( 3 ) C ( 3 ) = A ( 1 ) * 0 ( 2 ) - A ( 2 ) * 3 ( 1) RETURN END S U B R O U T I N E D E A N C I ( A , B , N , K, LT , CET ) IMPLICIT REALMS(A-H , 0 - 2 ) CCMMON /ZDE T/ CE,NCN COMMON / Z C C N / CCND C I KENS ION A ( 1 ) , B ( 1 ) DOUBLE P R E C I S I O N C S C R T , D A G S , D S I G N I F (M .EQ . 1 ) GC TC I C C MM=M.-1 NM=N*M NM 1 =NM-MM I F ( L T . N E . l ) GC TO 5 5 CCNC=1.0 NCN = C DE = C . MP=M+1 KK=2 FAC=DET NRGW=1 I F ( A ( 1 J . E Q . C . ) GC TO 60 SLM=O.D0 DC 7 7 I = 1,M 77 SUM=SUM + A( 1 ) * A ( I ) S=1.O/OSCRT(CAPS(A(1))) CE=A( 1 ) A ( 1) =S  97  836 £37 838  IF ( A ( 1 ) . L T .0 . ) A ( 1 ) = - S CCND=CCND/(A(1)*A(1)*CSCRT(SUM)) SUM = A ( 2 ) * A ( 2 ) + . A ( M P ) * A ( M P )  3"W £4C 8*1 842 843 844 '•5  M2=2*M BIGL=DABS(A(1)) S M L =B I G L MPP=MP+1 CC 8 8 I=MPP,M2 SUM = SUM + A ( I ) * A ( I ) A ( 2 ) = A ( 2 ) * D A E S ( A ( 1) )  68  ,;4b P<<7 48 849 .850  16  g-5l >P5 2 853 854 855 856 57  S=A(MP)-A(2)*A(2)*DSIGN(1 . D O , A ( 1)) IF(S.NE.C) GO T C 16 NRCW=2 GO TO 6 0 A ( M P ) = 1 . C / C S G R T ( C A E S (S ) ) CE=0E*S C C N D = C O N D / (A ( M P ) * A ( M P ) * D SGR T ( SUM ) ) IF(S.LT.O.) A(MP)=-A(MP) AAA = D A B S ( A ( N P ) ) I F ( A A A . G T . B I G L ) fiIGL = A A A I F ( A A A . L T . S M L ) SML=AAA I F ( N . E Q . 2 ) GC T C 5 3  e58 P59  MP=MP+K DC 6 2 J = M P , N M 1 , M  8,60! 861 862  J P = J-MM MZC=0 IF(KK.GE.M)  863 -364 855 816 6 867 868 869 S7C / 1 872 •Kl 3 &74  2  GC T O 2 KK = KK + Mt L I =KK —MM JC=KK-NM CO 6 5 I=KK,JP,"M  64 65  I F ( A ( I ) . EQ . C . ) GC GC T O 6 6 J C = JC+M M,ZC = MZC + 1  it  ASUM1=0. G C T O 61 MMZC=MM*MZC  1  v.  875 4E76 8r>7 878 8J9 8"8C 881 *8 8 2 883 884 "35 ^86 887 •888  GC T C  1  K K = K K+1 1 1=1 JC=1  TO 6 4  I 1=1 l + H Z C KM=KK+MMZC SUM=A(KM)*A(KM) A (K M ) = A ( K M ) r C A 8 S ( A ( J C ) ) ASUM1 = A(KM ) * A ( K M ) * D S I G N ( 1 . D O , A ( J C ) ) I F ( K M . G E . J P > G C TC 6 KJ =KM+V M JJ=JC CO 5 I=KJ,JP,MM S U M = S L M + A ( I ) * A (I ) JJ=JJ+M  8 89 890 89 1  A SUM 2 = 0 . IM=I-MM I 1 =1 1 + 1  £92 893 894  KI=II+MMZC KZ=JC CC 7 K =KM , 1 M , f> M  895  ASUM2 =  ASIN2+A<KI)*A(K)*DSIGN(1.C0,MKZ))  896 897  KZ=KZ+M KI=KI+KK  7  A(I)  8 98  = (A(I)-ASUf2)*CABS  A S I V 1 = A S U H 1 + A( I ) * A (  8^3 9  sec  5  CCNT INUE  CGI  6  CCNTINUE JMN=J+KV  90 2  CO  903 9C4  3  , r  I =J , J M V  SUK=SUM+A( I ) * A ( I IF(S.NE.C. )  .C6 r  l  NROW=(J + f N  CQ 63  9C9  I)* DSIGN(1.00 ,A ( J J ) )  )  S =A ( J) - A S I M  61  £•^5  3  99  ( A (Kl ) )  GC  TC  63  )/"  GC TO 60 A(J)=1.0/DSQRT(CABS(S) )  9 10  CCN D = C O N D / ( A ( J ) * A ( J ) * D S O R T ( S U N ) )  941 9 12  CE=CE*S IF(CABS(DE).GT.l.E-15)  9.13 914  CE=DE*1.E+15  915  GC 144  TO  145  TO  145  CCNTINUE  ^0  IF(S.LT.C.C)  21 922  AAA=DABS ( A (J ) ) IF(AAA.GT.EIGL)  923  I F ( A A A . L I .SNL )  924  62  CCNTINUE  925 9,2.6  6C  GC T O 5 3 WRITE(6,99)  927  54  928  99  929  A ( J ) = - A ( J ) E I GL = A A A SML = A A A  NRCV.  DET=0. FO RM A T ( 3 5 HC E R R O R  CONDITION  E NCCUNT ERED  IN  RETURN  930  53  931  55  DET = S M L / 8 IGL B(1)=B(1)*CAES(A(  932  KK=1  ?;3 3 9-3 4  Kl=l  1))  J=l L4=l  935 <9.36  CC  9^7  8  L = 2,N  BSUN1=0. -  LM=L-1 J = J+M  939 940  I F U K . G E . N )GC  941  KK=KK+ 1  94 2  GC 12  944 45 ->46  GC  145  IF(DA3SIDE).LT.l.E+15) NCN=NCN+15  9 18 9 19  9*3  144  DE=CE*1.E-15  917  938  TO  NCN=NCN-15  9 16  r  GO  TO  TC  12  13  KK=KK+M K1=K1+L L4=L4+M  13  J K = KK  947  L5 = L4  -9 4 8  CC 9 K = K1 , LN BSUN1=3SIN1 + A ( J K ) * B ( K ) * D S I G N ( 1 . C 0 , M L 5 ) )  949 950 '•95 1  JK=JK+MM L5=L5+M  9 52  9  953 954  8  955  CCNT INUE B(L) = (B(L ) - B S U K l ) * C A B S ( A { J ) B (N ) = B ( N ) * A ( Nu 1 ) NN N=NK i  )  ROW,16)  956  NN=N-1  957  NC=N  950  DC  100  9.5 9.'  B SDN 2 = 0 .  960  NL=N-L  10  L=l,NK •  961  NLl=N-L+l  9.6 2  NNN=NMM-M  963  NJ I =NMM  964  I F ( L . G E . N )N0=NC- 1  •65  DU  - .9&6  11  K=NL1,NC  NJ 1=NJ1+1  ' 9,7,6 7  B S L M 2 = B S L N 2+A  (NJ1)*B(K)  '.P, 68  11  CCNTINUE  '969  10  B(NL) = (B(NLJ-ESUN2)*A(NMM )  ICC  DO  :  97C 97 1  RETURN  •97 2 9-,73  101  1=1tN  IF(A(I).FQ.C.C ) 101  9:74  CO  TO  60  E ( I ) = B ( I )/A ( I ) RETURN  975  END  976  SUBROUTINE  977  IMPLICIT  9'78  CCMMON  9,79  C I M E N S ION  AVMATC(V,R,XX,SN )  REAL*8(A-H.C-2 )  / B L K 2 / C 4 , C 5 , S 4 , S 5 , C 1 0 , C 11 , S i C , S I 1 , SMPR E V , A 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  98 1  *,A(6,12),v(12f6,12),XX(3)fSN(6)*R(  982 •S8 3 ^ 8 4 ,9.8 5  DO  1  DC  1  1 = 1,6  DC  1  K= l,12 C.CC  V(K,I,J)=O.CC  9E7  1  CCNTINUE  988  6  C4=UC0S(R(4  9 89 -9 9 C $9 1  12),SMPREV  J=l,12  A(I,J)=  9\86  C(6,6),DDWI(3,2,2,2),DDfcJ(3,3,3,3)  ) ) *SN(2)-DS IN(R(4 ) )*SN(1)  S4 = C C 0 S ( R ( 4 )  )*SN(  1)+DSIN(R(4))*SN(2)  C 5 = C C O S ( R ( 5 ) ) * S N (4 ) - D S S5 = C C O S ( R ( 5 )  IN(R  (5 ) ) * S N ( 2)  )*SN(3)+DSIN(R(5)  C6 = D C 0 S ( R ( 6  •1,9 3  S 6 = C C 0 S ( R ( 6 ) ) * SN ( 5 ) + D S IN (R ( 6 ) ) * S N ( 6 )  t9 9 4 99 5 i.996 %97 998 -§99  C11 = D C 0 S ( R SI  >)*SM  (11 ) ) * S N ( 4 ) - C S 1 N ( R ( 1 1 )  C10 = DCOS(R (1C ) ) *SN ( 2 )-DS S12 =D C 0 S ( R ( 1 2 )  IN (R ( 1C ) ) * S N ( 1 )  ) * SN ( 5 ) + C S IN ( P, ( 1 2 )  C12 = D C O S ( R ( 1 2 ) ) * S N ( 6 ) - D S Ml ( 1 ,1 ) = C 5 * C 6 WI  'wl(2,l)=  ~T.C04  WI(2,2)=  1-CC5  WI(2,3)=  1CC6  WI(3,1)=  1C07  WI(3,2)=  •1C08  Kl (2,3)=  S5 C4*S6 + S4*S5*C6 C4*C6-S4*S5*St - S4 * C 5 S4*S6-C4*S5*Cfc S4*C6+C4*S5*S6 C4*C5  WJ ( 1, 1)=C1 1 * C 1 2 WJ(1,2)=-Cll*S12  1C11  WJ(1,3)=  1C12  WJ(2,1)=  C1C*S12+S10*S11*C12  1C13  KJ(2,2)=  C1C*C12-S10*S11*S12  1C14  WJ(2,3)=-SlC*Cll  1C15  WJ(3,1)=  SI  )*SN(6)  IN(R(12))*SN(5)  (1,2)=-C5*S6  W I ( 1,3)=  U C 3  1C1C  )*SN(3)  S1C=DCOS(R(1C))*SN(1)+CSIN(R(10))*SN(2)  ICC 1  1CC 9  IN(R(6))*SN(5)  1 = D C 0 S (R ( 1 1 > ) * S N ( 2 )+D S I N (fl ( 1 1 ) ) * SN ( 4 )  U C C 1 C C2  6J-DS  )*SN(4)  " 992  I  S1C*S12-C1G*S11*C1?  (6,6)  K 16 1C 17 IC 18 .019 C20 C21 C2 2 C23 C24 •'2 5 i. C 2 6 1C27 ICT28 1C29. l,C3Ci 10 31 K 32 1 02 3 10 4 IC35 1C36 LC27 C 3 8| LC>39 I J  21  :  r  '•1 £4 c 1C41 IC42 l\<3  1C44 LC&4 5 1^46 1C47 1C48 1C49 K 5C 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  23  20  WJ(3,2)= 51C*C12+C10*S11*S12 WJ ( 3 , 3 > = C 1 C - C 1 1 AA( I )=XX ( 1 )R( + 7 )-R( 1) 101 * M 2 ) = X X ( 2 M R(8 )-R{2 ) A A ( 3 ) = XX( 3 ) + R ( 9 ) - R ( 3 ) DO 20 K = l , 3 CC 20 1=1,3 DC 20 J = l , 3 GC TO( 21 , 2 0 * 2 3 ) ,K I F ( I EQ . 1 ) C W I ( K , I , J ) O.DO I F ( I EQ.2)CW I ( K , I , J ) -W I ( 3 , J ) I F ( I. E Q . 3 J C W I ( K , I , J ) ' W 1 ( 2 , J ) I F ( I E Q . l ) C W J ( K , I , J ) = G.DO I F ( I .EQ.2 )CWJ ( K , I , J ) - W J ( 3 , J ) I F ( I E,C . 3 ) D W J ( K , I , J ) W J ( 2 , J ) GCTC20 I F ( J . E Q . l ) G W I ( K , I , J ) W K 1,2) I F ( J . E C . 1 ) D W J ( K , I , J ) = WJ ( 1,2 ) I F ( J . E Q . 2 )DW I (K , I , J )-WK 1,1) IF(J.EQ.2)CWJ(K,I,J) -WJ(I,1) I F ( J . E Q . 3 > C W I ( K , I , J =0 ) .DO I F ( J . E Q . 3 ) C W J ( K , I , J ) = 0 .DO CCNTINUE CW 1 ( 2 , 1 , 1 ) - S 5 * C 6 S5*S6 DWI(2,1,2) DWI ( 2, 1,3) C5 S4*C5*C6 DWK 2 , 2 , 1 ) DWI ( 2 , 2 , 2 ) - S 4 * C 5 * S 6 CVi I ( 2 , 2 , 3 ) S4*S5 DW I ( 2, 3, 1 )- C 4 * C 5 * C 6 CWI ( 2 , 3 , 2 ) C4*C5*S6 DWI(2,3,3) -C4*S5 D W J ( 2 , 1 , 1 ) -S 11*C 12 DWJ(2,1,2 ) =S 1 1 * S 1 2 0WJ(2,1,3)= C l l DWJ(2,2,1 ) =S 1C*C11*C12 D W J ( 2 , 2 , 2 ) -S 10*C 1 1 * 5 1 2 D W J ( 2 ,2 ,3 ) SIG* S 1 1 DWJ(2,3,1) -C1C*C11*C12 DWJ(2,3, 2) = C10*C11*S12 DWJ (2 , 3,3 )=--C 1C*S 1 1 WJ(2,3)*WI(3,L)+WJ(3,3)*WI(2,1) T( 1) T ( 2 ) 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 T ( 3 ) W J ( 1 ,3 )*H I ( 1, 2 ) + WJ ( 2 , 3 )*WI ( 2 , 2 )+WJ ( 3, 3 )*WI ( 3,2) T( 1 ) T( 4) T (5 ) 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 ) P ( 1 ) -W J ( 2 , 3 ) *W I ( 3 , 2 ) + W J ( 3 , 3 ) * W I ( 2 , 2 ) P( 2) W J t l ,3 )*S'5*S6-W J ( 2 , 3 ) *S4 *C5 * S 6 +H J ( 3 , 3 ) * C 4 * C 5 * S 6 P(3) WJ( 1 , 3 ) * W I ( 1 , 1 ) - W J ( 2 , 3 ) * W I ( 2 , 1 ) - W J ( 3 , 3 ) * W I ( 3 , 1 ) P ( 4 ) -P ( I ) ? (5 ) C 1 1 * W I ( 1 , 2 ) + S10*S11*V»I(2,2)-C10*S11*KI(3,2) WJ(2,2)*WI(3,1)+WJ(3,2)*WI(2,1) Q ( 1 ) C(2 ) WJ(1,2)*S5*C6 +WJ(2,2)*S4*C5*C6-WJ(3,2)*C4*C5*G6 0(3) UJ(l,2)*WI(l,2)+fcJ(2,2)*WI(2,2)+WJ(3,2)*Wl(3,2) C ( 4 ) -G ( 1 ) G(5)= S 1 1 * S 1 2 * W I ( 1, 1 ) - S 1 0 * C l l * S l 2 * V i I (2 ,1 ) + C 1 C * C 1 1 * S I 2*WI ( 3 , 1 ) C(6)=-WJ(l,l)*WI(l,l)-WJ(2,l)*WI(2,l)-KJ(3,l)*WI(3,l) DC l t » J = l , 3 DC 16 K = 1,2 A ( J , K ) = - W I( K ,J)  C76 C77  A( J , K + 6 ) = - M J , K ) A ( J , 5 ) = A A ( K )*CWI ( 2 , K , J ) + A ( J , 5 J  102  C78 C7<5 C«0 C.B1  16  CCNTINUE A ( J , 4 ) = - A A ( 2 ) * W I ( 3 , J ) + A A ( 3 ) * w l ( 2 » J ) A (1,6 )= AA(J)*Wl(J,2)+A(l,6) A ( 2 , 6 ) = - £ A ( J ) * W I ( J , 1)+A(2,6)  C82 C83 C84  15  CCMINUE DC 1 2 K = l , 3 A(5,K+3)= T(K) A ( 4 , K + 3 ) = - P (K ) A ( 6 , K + 3 ) = - C (K ) IF ( K . E G . 2 ) G C T C 1 2 A(5,K+9)= TtK+3) A(4,K +9)=-P(K +3 ) A ( 6 , K +9 ) = - C ( K + 3) CONTINUE  s  U35 , Lfc'fe 10.87 1 4. £ 8 8 •C8S C90 C*91  12  G92 :CS'3 .C94  A(6,12)= -C(6 ) DC 8 4 N = l , 3 DC 8 4 K = l , 3  C95 C96 C97 C§8 xC99 UPO 101 102 l'cC3 ,104 1105 11X6 1 107 1C8 1C9 110 .ia 1  DO  85 ;  1  84  1=1,3  CC 8 4 J=l,3 G C T O ( 8 5 , 6 4 , 8 6 ),K I F ( I . E Q . l )CCW K N , K , I , J ) = C . D O IF(I.EQ.2)CCWI(N,K, I, J ) =-OW I (N , 3 , J ) IF( I.EC-.3)CCV I (N,K,I ,J)= CWI(N,2,J) IF( l . E Q . l ) C D W J ( N , K , I , J ) = 0.D0 IF<I.EQ.2)CCWJ(N,K,I,J )=-DWJ(N,3,J> I F U . E C . 3 ) C C W J ( N , K , I ,J)= DWJ(N 2,J) f  86  84  i l 2  GCTC84 IF(J.EC.l)CDWI (N,K,I,J)= CWI(N,I,2) IF(J.EC.l)DDUJ(N,K,I ,J) = CfcJ(N,I,2) IF(J.EQ.2)DDU(N,K,I,J)=-CWI(N,I,1) IF(J.EQ.2)C0WJ(N,K,I,J)=-DWJ(N,I,1) IF ( J . E G . 3 )CCV I ( N , K , I , J ) = 0 . C O IF(J.EQ.3)CCWJ(N,K,I CCNTINUE  , J ) = 0 . CO  CC 8 7 J = 1 , 3 CCWI ( 1 , 2 , 1 , J ) = O . C C C C W J ( 1 , 2 , 1, J ) = C . C C CCWI ( 3 , 2 , J , 3 ) = 0 . C 0 C C W J ( 3 ,2 , J , 3 ) = C . C 0 CCWI(3 , 2 , J , I ) = CWI(2 , J , 2 ) C C W J ( 3 , 2 , J , 1 ) = DWJ(2, J,2 ) CCW I ( 3 , 2 , J , 2 ) = - D W I ( 2 , J , 1 ) CCWJ( 3 , 2 , J , 2 ) = - C W J ( 2 , J , l ) C C W I ( 2 , 2 , 1 , J )=-WI( 1,J) C C W J ( 2 , 2 , 1 , J ) = - W J ( 1, J ) C C W I ( 1 , 2 , 2 , J ) = - C W I ( 2 ,3 , J )  3 1M4 115 4,16 117 118 U119  h ZO A  1121 {-122 1123  87  CCWJ( 1 , 2 , 2 , J )= - C W J ( 2 , 3 , J ) C C W I ( 1 , 2 , 3 , J )= C W I ( 2 , 2 , J ) CCWJ(1,2,3,J)= CWJ(2,2,J) CCNTINUE CCWI ( 2 , 2 , 3 , 1 CCWJ(2,2,3,1 CCWI ( 2 , 2 , 3 , 2 CCWJ ( 2 , 2 , 3 , 2 CCWI ( 2 , 2 , 3 , 3 CCWJ{2,2,3,3 CCWI ( 2 , 2 , 2 , 1 CCWJ(2,2,2,1  )= C 4 * S 5 * C 6 )= C 1 C * S 1 1 * C 1 2 )=-C4*S5*S6 )=-C10*Sll*S12 )= - C 4 * C 5 )=-ClC*Cll )= - S 4 * S 5 * C 6 )=-S10*S11*C12  1136 1137: 11 3 8 1 V 3 9]  C C W I ( 2 , 2 , 2 , 2 )= S4*S5*S6 C C W J ( 2 , 2 , 2 , 2 )= S1-C*SU*S12 CCW I ( 2 , 2 , 2 , 3 - ) = S4*C5 CCWJ(2,2,2,3)=C11*S10  •140| 1 4 1! i;42 143! 144,  103  V (5 , 1, 1 ) = S 5 * C 6 V ( 6 , 1 ,1 ) =- W I ( 1 , 2 ) V(4,1,2)=WI(3,1) V(5 , 1 ,2)=-S4*C5*C6 V(6,1,2)=-WI( 2,2)  45| ^4'6, I 147 1 1148  V ( 4 , 1 , 3 )=-W I ( 2 , 1 ) V(5,l,3)= C4*C5*C6 V (6 , 1 , 3 ) = - W I ( 3 , 2 ) V ( 4 , 1 , 4 ) = - A A ( 2 ) * W I ( 2 » 1 ) - A A ( 3 ) * W I ( 3 , 1 )  I'14 9 I 1150, 11-5 1 !il'52 115 3 I £54  V(5,l,4)= A A ( 2 ) * C 4 * C 5* C 6 4 £ A ( 3 ) * S 4 * C 5 * C 6 V ( 6 , 1 , 4 ) = - A A ( 2) * W l ( 3 , 2 ) + /5 A ( 3 ) *W I ( 2 , 2 ) V ( 8 , 1,4)=-WI ( 3 , 1 ) V<9,1,4)= WI(2,1) V ( 5 , 1 , 5 ) = - £ A ( 1)*C5*C6-*A{2 )*S4*S5*C6+AA(3 )*C4*S5*C6 V(6,l,5)= AM 1)*S5*S6-AA(2 )*S4*C5*S6+AA(3)*C4*C5*S6  155  V(7,l,5)=  -S5*C6  .156  V(8,l,5)=  S4*C5*C6  137 .158 11,59 U60 • '  V (9, 1,5)=-C4*C5*C6 V ( 6 , 1 , 6 ) = - £ A ( 1 )*WI ( 1 , 1 ) - A A ( 2 ) * W I ( 2 , 1 ) V(7,l,6)= WI(1,2) V (8,1,6)= WI(2,2)  116 1 162 L6 3  V ( 9 ,1 , 6 ) = W I ( 3 , 2 ) V ( 5 ,2 , 1 ) = - S 5 * S 6 V (6 , 2 , 1 ) = W I ( 1 , 1 )  .l 64 I 1^65 ll £6 .167 168  V(4,2,2)= W1(3,2) V ( 5 , 2 , 2 )= S 4 * C 5 * S 6 V ( 6 , 2 , 2 ) = WI ( 2 , 1 ) V(4,2,3)=-WI( 2,2) V (5,2,3) =-C4*C5*S6  r  r  169 170 1<71  -  —AA ( 3 ) * W I (3,1)  V(6,2,3)= WI(3,1) V(4,2,4)=-AA(2)*WI(2,2)-AA(3)*WI(3,2) V ( 5 , 2 , 4 ) = - « A ( 2 ) * C 4 * C 5 * S 6 - A A ( 3 ) *S 4 * C 5 * S t  t'172 U<»7 3 4>74 1175  V I 6 , 2 , 4) = A A ( 2 ) * W 1 ( 3 , 1 ) - * M 3 ) * W 1 ( 2 , 1 ) V(8,2,4)=-WI(3,2) V(9,2,4)= WI(2,2) V(5,2,5)= 2 A ( 1 ) * C 5 * S 6 + A A (2 ) * S 4 * S 5 * S 6 - A A ( 3  .jl76 177 178 1?9  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 V(7,2,5)= S5*S6 V(8,2,5)=-S4*C5*S6 V(9,2,5)= C4*C5*S6  „ l'SO 1181 U82 116 3 -.-184 II € 5 11-8 6  V ( 6 , 2 , 6 ) = - A A { 1) *W1 ( 1 , 2 ) - A A { 2 ) * W l ( 2 , 2 ) - A A ( 3 )*WI V(7,2,6)=-WI (1, 1) V(8,2,6)= -Wl(2,1) V I 9 , 2 , 6 ) = - W I ( 3 , I) V ( 5 , 3 , 1 )= - C 5 V ( 4 , 3 ,2 ) = W I ( 3 , 3 ) V(5,3,2)=-S4*S5  .  )*C4*S5*S6  1187 '.188  V(4,3,3)= V ( 5 , 3 , 3 )=  U«9 190 1:! > 1 1192  V(4,3,4)= -A/(2)*WI(2,3)+AA(3)*CwI(l,2,3) V(8,3,4)=-WI(3,3) V (9 , 3, 4 )= W I { 2,3) V ( 5 , 3 , 5 ) = - A A ( 1 )*S5 +AA(2)*S4*C5 -AA(?)*C4*C5  1193 .194 .19 5  V(7,3,5)= C5 V(8,3,5)=S4*S5 V(9,3,5)=-C4*S5  :  -W I ( 2 , 3 ) C4*S5  ( 3,2 )  I ] 196 1197 1198 199  PC 120 J = l , 3 D T ( I , J ) = - W J ( 2 , 3 ) *DW I ( J , 3 , l ) + k J ( 3 , 3 ) * 0 K l ( J , 2 , l ) DT (1 , J + 3 ) = - C W J ( J , 2 , 3 )*V* I ( 3 , I ) - C W J ( J , 2, 3 ) * K I ( 2 , 1) 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  v  -H5  DT(4, J)=-DT(1,J) CT (4 , J + 3 ) = - C T ( 1, J + 3 ) 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 ) *  1 • ?C6, 1 ?C7i .? C 8, 12C9 1210 t ? l l 1 i 12 12, 13 121-*, 12 15, 121; I21V 12 K '1219 |l?20'i 1221 12 2 2 H 2 3!  *CWI(J,3,1) DT(5,J+3)=CDfcJ(J,2,l,3)*wI(l,l) + C0ViJ(J,2,2,3)*WI(2,l)+0 0WJ(J,2,3,3 *)*WI(3,1) C F ( 1 , J ) = - k J ( 2 , 3 )*DW I ( J , 3 , 2 ) + W J ( 3 , 3 ) * D W I ( J , 2 , 2 ) CP(l,J+3)=-DVJ(J,2,3)*V»I(3,2) + CWJ(J,3,3)*WI(2,2) * JDP, 2( 2, ,3 J, )2=)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 ) * C C H i ( 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)+DWJ(J,3,3)* *DK 1 ( 2 , 2 , 2 ) DP(3,J)=-WJ(1,3)*DW I ( J , 1 , 1 )-WJ(2,3)*DWI(J,2, 1)-WJ(3 , 3)*DWI(J,3,1) 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)*WI ( 3 , 1 * ) DP(4 , J ) =-DP(1,J ) DP(4,J+3)=-CP(l.J+3) DP(5,J)=DWJ(2,l,2)*DWI{J,l,2)+CKJ(2,2 3)*DWI(J,2,2)+DWJ(2,3,3)* *DWI(J,3,2) t  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 )*ViI( 2 , 2 ) + D D W J ( J , 2 , 3 , 3 * ) * W I (2 ,2 ) D C ( 1 , J ) =- Vs J ( 2 , 2 ) *DVs I ( J , 3 , 1 ) + W J ( 3 , 2 ) * D W I ( J , 2 , 1 ) CC(l,J+3)=-CkJ(J,2,2)*WI(3,l )+CWJ(J,2,2)*K1(2,1) D C ( 2 , J ) =V* J (1 , 2 ) * C C W I ( J , 2 ,1 ,1 ) + V v J ( 2 , 2 ) # B D W l ( J » 2 , 2 , l ) +W J ( 3t 2 )*CCW I ( J  12 2 5 12/ 122; 12 2 6 1229 1-2 3 G 1231 12 3 2  *,2,3,1) 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 ) * *DU ( 2,3,1 ) 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) DC(3,J+3)=CttJ(J,l,2)*WI(l,2 )+ CWJ(J,2,2)*WI(2,2)+DWJ(J,3,2)*WI(3,2) DC(4,J)= -DQ(1,J) DC (4 , J + 3 ) = - 0 C ( 1 , J + 3 ) 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 ) *  I-C? 3 3 1-2 3 4, ' 1 2 3 5,  1A,2 3 6 12,37 12 3 8 | "39 ;  *DkI ( J , 3 , 1 ) 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  11  i240; 1241  1255  )*CCK1(  * J , 2, 3, I ) 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) * * C M ( 2,3,1 ) CT ( 3 , J ) = VtJ(I,3)*DWl(J,l,2)+kJ(2,3)*CUI(J 2,2)+WJ(3,3)±DWl(J,3,2) 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 )  .iCO 1?01 12X2 L2C3 12C4  '1242! '124 3 i 1244 1-2 4 5 J 124 6 , 1247 12 4 8 1249 , 1250 , ' I 2 5 1 12 52 1253 1254  104  12C  *)*WI(3, 1) CC ( 6 , J ) = - V « J ( I , 1 ) * D W I ( J , 1 , 1 ) - W J ( 2 , 1 ) * D W I ( J , 2 , 1 ) - W J ( 3 , 1 ) * D W I ( J , 3 , 1 ) C C ( 6 , J +3)= - D f c J ( J , 1 , 1 ) * f c I(1,1)-CWJ(J,2,1)*WI(2,1)-DWJ(J,3,1)*WI(3,1 * ) CCNTINUE V(4,5,4)=DT(1,1) \M5,5,4)=0T(1,2) V(6,5,4)=CT(1,3) V ( 1 C , 5 , 4 ) = C T ( 1,4 ) V( 1 1 , 5 , 4 ) = C T ( 1 , 5 ) V(5,5,5)= CT(2,2 ) V(6,5,5)= CT(2,3) V( 1 C , 5 , 5 ) = C T ( 2 , 4 ) V(11,5,5)=0T(2,5) V(6,5,6)= CT(3,3) V(1C,5,6)=DT<3,4) V(11,5,6)=CT(3,5)  12 5 6  V(10,5,10)=CT(4,4)  1 2 57  V( 1 1 , 5 , 1 C ) = C T ( 4 , 5 )  1258  V ( U , 5 , 1 1 ) =CT(5,5)  259  DO  31  J= l,3  126C  DC  32  K=l,3  1--2 6 1  l  I F ( K . G E . 3 >GOTO  32  1263  V (K + 9 , 4 , J + 3 ) = - C P ( J , K - » 3 )  1264  IF(J.GE.31GCTC32  '26 5  V(K+9,4,J+9)=-DP(J+3,K+3  i'266  32  '267  31  )  CCNTINUE CCNTINUE  .266  DC  41  J= l,3  1269  DO  42  K=1,3  127C  V(K+3,6,J*3)=-CC(J,K)  \\27 1  V( K + 9 , 6 , J 4 3 ) =-DC  1272  V ( K + 9 , 6 , J + 9) = -DG( J + 3.K + 3)  1*27 3  42  CCNTINUE  lt274  41  CCNTINUE  1275  DC  51  K=l,12  1276  DO  52  J= l,6  1V277  DC  53  12 7 8  V(I,J,K)=  ) 279  53  ;  1281  c  V(K+3,4,J+3)=-CP(J,K)  v 1'262  );280  n  ,  { J ,K + 3 )  1=1,12 V(K,J,I)  CONTINUE  52  CCNTINUE  51  CCNTINUE  1282  RETURN  V28 3  ENC  1*284  SUBROUTINE  j>285  IMPLICIT  it? 8 6  CCMMON  1287  DIMENSION  A ( 6 , 1 2 ) , V ( 1 2 , 6 , 1 2 ) , SMPREV(6,6) ,XX(3>,SN(6)  DIMENSION  W I I ( 3 , 3 ) , W J J { 3 , 3 ), A * ( 3 ) , > S I < 3 ) , BS ( 3 ) , Y S J ( 3 ) , Y SK ( 3 ) ,  1288  „  AV M A T S (V , R , XX , SN )  RE A L*8( A-H , C- 2 )  /B L K 2/  C 4 , C 5 , S 4 , S 5 , C 1 0 , C 1 L , SI C , SI 1 , SM P R E V , A ,R(12)  1289  *ZSI(3),CS(3),ZSJ(3),ZSK{3),DWI(3,3,3),DWJ(3,3,3),DEL(9),DYS1(3,9),  129C  *DBS(3,9),DCS(3,9),DYESJ(9 ) ,DZESJ(9) ,DYSJ(3,9),DZSJ(3,9),DYSK(3,9) ,  It25 1  * D Z S K ( 3 , 9 ) , C Y IS I ( 9 ) , C Z I S I ( 9 ), D Y I S J ( 9 ) , D Z I S J ( 9 ) , C Y I S K ( 9 ) , D Z I S K ( 9 )  129 2 v  DIMENSION  .'293  CDEL(9,9),CDYSI(9,3,9),DCWI(3,3,3,3),CDWJ(3,3,3,3),  * C C B S ( 9 , 3 , 9 ) , C D C S ( 9 , 3 , 9 ) , D D Y E S J ( 9 , 9 ) , D D ZE S J ( 9 , 9) , DD Y S J ( 9 , 3 , 9 ) ,  12 9 4  *CCZSJ(9,3,9),CCYSK(9,3,9),DDZSK(9,3,9),DDYISI  1295  *CCYISJ(9,9),CCZISJ(9,9)  ,1296  DO  1 J=  1:297  DC  1  ,DDYISK(9,9)  l,12  1=1,6  DC 1 K = l , 1 2  1298 -2299  A(I,J)=  I^3CC  V ( K , I , J ) = O.CC  C.DC  13C1  1  CCNTINUE  1302  6  C4=DCCS(P(4))*SN(2)-DSIN(R(4  1*3 0 3  S4 = D C 0 S ( R ( 4 )  1304  C5 = CCO  C5  )*SN(  1HDSIN(R(4)  ) )*SN(1) )*SN(  2)  S(R(5))*SN(4)-DSIN(R(5))*SN(2)  S5 = C C G S ( R ( 5 ) ) * S M 3 ) + C S I N ( R ( 5 ) ) * S N < 4 )  i3C6  C 6 = DCO S ( R ( 6 ) ) * S N ( 6 ) - D S  13 0 7  S6 = C C 0 S ( R ( 6 ) ) * S K ( 5 ) + D S I N ( R ( 6 ) ) * S N ( 6 )  <13C8  I N ( R ( 6 ) )*SN(  5)  C11=DCCS(R(11))*SN(4)-CSIN(R(11))*SN(3)  i3C9  SI  1310  C1C =D C 0 S ( R ( 1 C ) ) * S N ( 2 ) - C S I N ( R ( 1 C ) )*SN(1 )  1 = DCO S ( R ( 1 1 ) ) * S N ( 2 ) + C S I N ( R ( I 1 ) ) * S N ( 4 )  1311  S 1 C = D C 0 S ( R ( 1 C ) ) * S N ( 1 ) + C S IN ( R ( 1 0 ) ) * S N ( 2 )  1312  S 1 2 =DC 0 S ( R ( 1 2 ) ) * S N ( 5 ) +D S I N ( R ( 1 2 ) ) * S N ( 6 )  1313  C12 =D C 0 S ( R ( 1 2 ) ) * S N { 6 ) - D S I N ( R ( 1 2 ) ) * S N ( 5 )  1314  WII(1  1315  WII(1,2)=-C5*S6  ,1  )=C5*C6  ( 9 , 9 ) ,DDZI S I ( 9 , 9 )  ,CCZISK(9,9)  ,  I  1316 1317 1318 ' 3 19 1320  WII( 1,3 ) = S 5 WII(2,1)= C4*S6+S4*S5*C6 WII(2,2)= C4*C6-S4*S5*S6 W I I ( 2 , 3 )= - S 4 * C 5 WI I ( 3 , 1 ) = S4*S6-C4*S5*C6 Wll(3,2)= S4*C6+C4*S5*S6 WII(3,2)=C4*C5 WJJ( 1,1)=C11*C12 WJJ( 1 , 2 ) = - C 1 1 * S 1 2 WJJ(1,3)= S l l WJJ(2,1)= C1C*S12 + S1G*SU*C12 WJJ(2,2)= C1C*C12-S10*SU*S12  h?.2l , 2 22 1323 1224 '225 1326 ' 227 228 1329 123C \3 2 l 1*3 3 2 ^3 2 3 f334 13 3 5 1336 123 7 1238 V339 ;34b 1341 134 2 11343 l!34A ^345 £346 1247 1248 1349 <125C 1251 1352 13 5 3 <125f 1355 ,13 5 6 r,2 5 7 1358 359 136C 136 1 136 2 113 6 3 1264 -1365 1366 1267 ,1268  -  -  WJJ(2,3)=-S1C*C11 WJJ(3,1)= SIC*S12-C 10*S11*C12 WJJ(3,2)= S1C*C12+C10*S11*S12 WJJ(3,3)=C1C*C11 AA(1)=XX<1)+ R(7)-R(l) AA(2)=XX(2M R(8)-R(2) AA(3)=XX(3)+ R(9)-R(3) E L =D S C R T ( A A ( 1 ) * * 2 + A A ( 2 ) * * 2 + AA ( 3 ) * * 2 ) Y S I (1 ) = A A ( 1 J / E L YSI(2)=AA(2)/EL YSI(3)= />A(3)/EL B S ( 1 ) = WI 1 ( 2 , 3 ) * Y S I ( 3 ) - V. I I { 3., 3 ) * Y S I ( 2 ) B S ( 2 ) = W I I ( 3 , 3 ) * Y S I (1 )-WI I ( 1 , 3 ) * Y S I ( 2 ) B S ( 3 ) = WII ( 1 , 3 ) * Y S I ( 2 ) - W I I ( 2 , 3 ) * Y S I ( 1 ) Y E S J = D S Q R T ( B S ( 1 ) * * 2 + B S ( 2 >**2 + B S ( 3 ) * * 2 ) YSJ(1)= ES(1)/YESJ YSJ(2)= ES(2)/YESJ YSJ(3)= BS(3)/YESJ CALL CROSS(YS I , Y S J , Y S K ) ZSI(I)=YSIU) ZSI(2)=YSI(2) ZS I ( 3 ) = Y S I ( 3 ) CS(1)= WJJ(2,3)*YSI(2)-WJJ(3,3)*YSI(2) CS(2)= WvJ(3,3)*YSI(l)-WJJ(l,3)*YSI(3) C S ( 3 ) = W J J ( 1 , 3 ) * Y S I( 2 ) - W J J ( 2 , 3 ) * Y S I ( 1 ) Z E S J = C S G R T ( C S ( 1 )* * 2 +C S ( 2 ) * * 2 + C S ( 3 ) * * 2 ) ZSJ(1)=CS(1)/ZESJ ZSJ(2)=CS(2)/ZESJ Z S J( 3 ) =C S (3 ) / Z E S J C A L L CROSS (ZS I , Z S J , Z S K ) YISI=WII(l,l)*YSI(l)+WII(2,l)*YSI(2)+WII(3 l)*YSI(3) 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 ) 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 ) 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 ) Z ISJ=WJJ(1,1 )*ZSJ( 1)+WJJ(2,1)*ZSJ( 2)+WJJ(3,1)*ZSJ(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 ) DC 2 0 K = l , 3 DO 2 0 1 = 1 , 3 DC 2 0 J = l , 3 GC T C ( 2 1 , 2 0 , 2 3 ) , K IF(I.EQ.1)CWI(K,I,J)=0.00 f  r  1;369 137C 1371 1372 1373 12 7 4 1375  106  21  23  IF( I . E Q . 2 ) C W I ( K , I , J ) = - W I I ( 3 , J ) I F ( I . EQ . 3 )CW I ( K , I , J ) = W I I ( 2 , J ) IF(I.EQ.l)DWJ(K,I,J)=0.D0 IF(I.EQ.2)CWJ(K,I,J)=-WJJ(3,J) I F ( I . E Q . 3 ) C W J ( K , I , J )= W J J ( 2 , J ) GCTC20 I F ( J . E Q . l )CW I ( K , I , J)=W I I ( I , 2 )  13 7 6 1377' 378! 79 t28C 12)81 13*8 2' 1383 !  ;  20  1384 '85 iBb "•£7 .88 389 39C 2^91 392 343 3 3-94 395 396 3\97 398 3^99; £CC 4C1 4C2 4\C3  IF(J.EQ.1)CWJ(K,I.J)=WJJ(I,2) I F ( J . E C . 2 )DW I ( K , I , J ) = -W I I ( I, 1 ) IF(J.EQ.2)CWJ(K,I,J)=-WJJ(I,i) IF(J.EQ.3)CWI(K,I,J)=O.D0 I F ( J . E C . 3 ) D W J ( K , I , J ) = 0 .DO CCNTINUE DW I ( 2 , 1,1 ) = - S 5 * C 6 DWI(2,1,2)= S5*S6 CWK 2,1,3)= C5 DWI ( 2 , 2 , 1 ) = S 4 * C 5 * C 6 CWI ( 2 , 2 , 2 ) = - S 4 * C 5 * S 6 DWI(2,2,3)= S4*S5 DWI(2,3,1)=-C4*C5*C6 DWI(2,3,2)= C4*C5*S6 CWI ( 2 , 3 , 3 ) = - C 4 * S 5 DWJ(2,1,1)=-S11*C12 UWJ(2,1,2)= S11*S12 DWJ ( 2 , 1 , 3 ) = C l l DWJ(2,2,1)= S1C*C11*C12 DWJ(2,2,2)=-SlC*Cll*S12 DW'J(2,2,3)= S1C*S11 DWJ( 2 , 3 , 1 ) = - C 1 C * C U * C 1 2 DWJ(2,3,2)= C1C*C11*S12 DWJ(2,3,3)=-C1C*S11 CC 30 J=l,3 DEL(J)= -YSI(J) DEL(J+3)=0.C DEL(J+6)=-CEL (J )  4'C4 4>C5: 4"C6  DC 3 1 K = l , 3 D E L (K ) = - Y S I (K ) D E L (K + 6 ) = - C E L (K )  4C7 4C8 4C9 41C  DYSKJ,K)=+DEL(J)*CEL(K)/EL IFU.EQ.KJDYSI(J,K)=DYSI(J,K)-1./EL D Y S K J,K +6 )=-DYSI ( J , K ) DYSI(J,K +3 )=C .0  4^11 412 M 3 414 415 416 4^17 418 -V19 420 421 422 .423 ;424 :4 2 5 .4 2 6 427 1428 14 2 9 1430 T431 143 2 1433 1434 14 3 5  107  31 30  »  CCNTINUE CONTINUE CC 4 0 J=l,3 DC 4 0 K=l,3 N = J+1 IF(K.EQ.4)f=l N = rV + 1 I F ( N . E Q . 4 ) N=1 DBS(J,K)= W I I ( M , 3 ) * 0 Y S I ( N , K ) - W I I (N , 3 ) * D Y S I ( K , K ) CCS(J,K)= WJJ{N,3)*DYSI(N,K)-WJJ(N,2)*DYSI(. ,K) DBS(J,K+6)=-CES(J,K) DCS(J.K+6)=-CCS(J,K) DES(J,K+3)= CWI(K,M,3)*YSI(N)-CWl(K,N,3)*YSI(M) DCS(J,K+3)= C W J ( K , V , 3 ) * Y S I ( N )-CWJ ( K , N ,3 ) *YS I ( * ) CCNTINUE DO 5 0 K = l , 9 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 V  4C  51 5C  CZESJ(K)=(nCS(l,K)*CS(l)+CCS(2,K)*CS(2)+CCS(3,K)*CS(3))/ZCSJ DO 5 1 J=l,3 DYSJ(J,K)=(CES(J,K)-DYESJ(K)*YSJ(J))/YESJ CZSJ(J,K)=(CCS(J,K)-CZESJ(K)*ZSJ(J))/ZESJ CCNTINUE CC 1 5 0 K=l,9 CC 5 2 J=l,3 N= J+ 1  1436  I F ( M . F 0 . 4 )M=1 N=P + 1 IF ( N . E 0 . 4 JN=1  1437 1438 •4 39^  DYSK(J,K )=  ! 440 '.44  1  i  Y S I ( K ) * O Y S J (N , K ) + D Y S I (M , K )* Y S J < N ) - Y S I ( N ) * D Y S J ( K , K )  -CYSI(N,K)*YSJ(M DZSK(J,K)=  1442  1  108  YSI  (N)*DZSJ(N»K)+CYSI  (M,K)*ZSJ(N)-YSI  (N)*DZSJ(.V,K)  -CYSI (N,K)*ZSJ(N' )  ' 1 4 4 3|  52  CCNTINUE  114 4 4  1 5C  CCNT INUE  4 4 5^  DO  -t 4 6  K ) + h I I ( 2 , 1 ) * D Y S I ( 2 , K ) + W I I ( 3 » 1 ) * D Y S I ( 3 , K ) C Y I S I ( K ) = V » I I (1,1 )*CYSI( DZISI(K)=WJJ(1,1)*CYSI(1,K)+VUJ(2,1)*CYSI(2,K )+WJJ(3,1)*DYSI (3,K)  4 4 7^  60  K=l,3  *48;  DY IS I (K + 6 ) = - D Y I S I ( K )  4 4 9;  C Z IS I (K + 6 ) = - C Z I S I ( K }  4 5C|  D Y I S I ( K + 3 ) =D W  ^ 5 1,  (K,1,1)*YSI  (1  DZISI(K +3)=CV>J(K,1,1)*YSI(1  ,4521  ) + C W I ( K , 2 , l ) * Y S l ( 2 )+ ( K , 2 , 1 ) * Y S I ( 2 )+  )+DV>J  DYISJ(K)=KII(1,1)*CYSJ(1,K)+WII(2,1)*DYSJ(2,K  DWI(K,3,1)*YS1 DWJ(K,3,l  (3)  )*YSI(3)  )+WlI(3,L)*DYSJ(3,K)  DYISJ(K+6)=-CYISJ(K) 14  DZISJ(K)=UJJ(1,1)*CZSJ(1, K)+WJJ(2,l)*0ZSJ(2,K)+WJJ(3,l)*DZSJ(3,K)  54  455  DZ I S J ( K + 6 ) = - C Z I S J ( K )  [456, L 4 57 s  14  58  D Z I S J ( K + 2) = V i J J ( l , l ) * D Z S J ( l , K + 3) +  1459; .46 1  D Y I S M K + 6 ) = - C Y I S K (K )  •$6 3  DZISK(K+6)=-CZISK(K)  .464  DY I S K ( K + 3 ) =  465  DZISK(K +3 ) =  4 67  1471 73  70 71  DEL ( K )  K=l,9  A(3,K)=(YISI*OYISJ(K)-YISJ*0Y  IS I ( K  ) ) / Y IS  1**2  M= 1 N=3 73  CC  72  K=M,N  1^77  A(5,K)=(-ZISI*CZISK(K)+ZISK*CZISI(K))/ZISI**2  i4*7 8  A(6,K)=(+ZISI*DZISJ(K)-ZISJ*0ZISI(K))/ZISI**2  W79  A(4,K)=  1480  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 )  1,K)*YSK(1)+ 72  IF ( K . E C . l  148 3  GCTC75 74  N= 9  I4 8 6  GCTC7 3  1487  75  )GCTC74  DC  76  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  1489  A(6,K+6)=(ZISI*CZISJ(K)-ZISJ*CZISI(K))/ZISl**2  149G  A(4,K)=  OA 9 1 149 2  DYSK ( 1 ,K ) * Z S J ( 1 ) + C Y S K ( 2 , K ) * Z S J ( 2 )  A(4,K+6)= 76  D Z S J ( 1 , K ) * Y S K ( 1 )+  CCNTINUE  1493  CC  80  J=l,3  |1494  DO  80  K=l,9  1495  +CZSJ(1  +CZSJ(3 , K)*YSK(3)  N=7  K',8 5  1488  CZSJ(2,K)*YSK(2)  CONTINUE  1482 1484  (3)  A(2,K)=(-YISl*DYISK(K)+YISK*DYISI(K ) )/YISI**2 71  :475  1481  + C W J ( K , 3 , 1 )*ZSK  >  K=1,9  A( 1 ,K) = DC  +DW I ( K , 3 , 1 ) * Y S 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 ) + V i J J ( 3 , 1) * D Z S  CCNTINUE DC  7C  A 74 1476  ( 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) +D * . J ( K , l , l ) * Z S K ( l ) + D W J ( K , 2 , l ) * Z S K ( 2 ) 6C  147 2 i,4  H I  * K ( 3 , K + 3 ) + D V i I ( K , 1 , 1 ) * Y S K ( 1 )+DWI ( K , 2 , 1 ) * Y S K ( 2 )  •Abb'  469  )*CZSJ  D Z 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 Z S K ( 3 , K )  46 2  i,4 7 C  toJJ(2,l)*DZSJ(2,K+3)+WJJ(3,I  + DWJ(K , 2 » 1 ) * Z S J ( 2 ) + D W J ( K , 3 , 1 ) * Z S J (3) 4(3,K+3)+CWJ(K,l,l)*ZSJ(l) C Y I S K ( K ) = M I ( 1 , 1 ) * C Y S K ( 1 , K)+WII(2,1)*DYSK(2,K)+WII(3,1)*DYSK (3,K)  l46Ci  .468  .  DYISJ(K+3)=WI 1(1,1) *DYSJ ( 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 ) *(3,K+3)+DM(K,l,l)*YSJ(l)  CCEL(J  ,K>=-CYSI(J,K)  +CYSK ( 3 , K ) * Z S J ( 3 )  DZSJ(2,K)*YSK(2)+  CZSJ(3 , K)*YSK(3)  1496  CCEL(J+3,K)=C.O  1497  CCEL(J+6.K)=-CDEL(J,K)  1498  80  .4 9 9  CCNTINUE DC  81  1=1,9  -15C0  CC  81  J=l,3  1501  DO  8 1  K=l,3  15 C 2  CCYSI(I,J,K)=+CCEL(J,I)*OEL(K)/EL+CCEL(K,  15C3  * - C E L I J ) * D E L ( K )*CEL( I  15C4  IF(J.EQ.K)  I ) * C EL ( J ) / E L  )/El**2  C C Y S I ( I , J ,K ) = D O Y S I (I , J , K ) + D E L ( I ) / ( E L * E L )  CCYSI(l,J,K+6)=-CDYSl(I,J,K) .5C6 5C7|  CCYSI ( I , J,K+3 81  i08  CC  84  N=l,3  15C9  DC  84  K=l,3  L 510  DC  l\511  DO  1512 1513  )=C.C  CONTINUE  84 84  1=1,3 J = 1, 3  G C T 0 ( 8 5 , 8 4 , 86 85  1/514  ),K  IF(I.EC.1)CCV>I(N,K,I,J)=0.C0 IF(I.EG.2)CDWl(N,K,I,J)=-CWl(N,3,J)  1515  IF(I.EQ.3)CCWI(N,K,I ,J) =  1516 1517  IF(I.EQ.2)DDWJ(N,K,I  1518  IF( I . E Q . 3 ) C C W J ( N , K , I , J )=  1519;  GCTC84  l";52C  DWI(N,2,J)  IF(I.EC.1)CCWJ(N,K,I,J)=0.CG  86  DWJ(N,2,J)  )CCWI ( N , K , I , J ) =  CWI(N,I,2)  1521  I F I J . E Q . l ) C D V J ( N , K » I , J)=  0WJ(N,I,2)  15 2 2  I F ( J . E C . 2 )CCW I ( N , K , I , J  l>52 3  I F ( J . E C . 2 ) C D V J ( N , K , I , J ) =-DW J ( N , I , 1 )  1524  ) = -CW I ( N , I ,  I F ( J . E Q . 3 1 C C W J ( N , K , I, J ) = 0 . C O 84  1527  CONTINUE DC  1528  87  J=l,3  *  CCWI(1,2,1,J)=0.DO  1529  CCWJ(1,2,1,J  1530  CCWI(3,2,J, 3)=C.CC  1531  )=C.DO  C C W J ( 3 , 2 , J , 3 )= C.DC  153 2  CCWI ( 3 , 2 , J , l  ^533  C C W J ( 3 , 2 , J , i ) =DWJ(2, J,2 )  ) = CWI ( 2 , J , 2 )  1*5 3 4  CCWI ( 3 , 2, J , 2 ) = - C W I ( 2 , J , 1 )  15 3 5  C C W J ( 3 , 2 , J , 2 )= - C W J ( 2 , J , l )  1,536  C C W I ( 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 (  lp39  C C W J ( 1 , 2 , 2 , J )= - C W J ( 2 ,3 , J )  1'540  CCWI ( 1 , 2 , 3 , J ) =  CWI(2,2,J)  1541  C C W J ( 1 , 2 , 3 , J )=  DWJ(2,2,J)  154 2'  87'  2,3,J)  CCNTINUE  1543  CCWI ( 2 , 2 , 3 , 1 )=  C4*S5*C6  1544  CCWJ(2,2,3,1)=  C 1 C * S U * C 1 2  , r  345  CCWI ( 2 , 2 , 3 , 2  )=-C4*S5 *S6 :  r5 4 6  CCV»J(2,2,3,2J=-C1C*S11*S12  1547  CCWI(2,2,2,3)=-C4*C5  1548 1549  I )  IF(J.EG.3)C0WI(N,K,I,J)=C.D0  15 25, 15.26  IF( J . E Q . l  »J)=-DWJ(N',3,J)  CCWJ(2,2,3,3)=-C1C*C11 CCWI(2,2,2,1)=-S4*S5*C6  15 5 0  C C W J ( 2 , 2 , 2 , 1 ) =-S  15 5 1  C C W I ( 2 , 2 , 2 , 2 )=  S4*S5*S6  1'552  C C W J ( 2 , 2 , 2 , 2 )=  S1G*S11*S12  1553  C C W I ( 2 , 2 , 2 , 3 )=  S4*C5  1554  C C U J ( 2 , 2 , 2 , 3 )=  C11*S1C  1555  CC  82  1=1,3  1C * S 1 1 * C 1 2  1556 1557 558 11559 ,1560 '1561 156 2 '1563 1564  DC 8 2 J = l , 3 DC 8 2 K = l , 3 M = J+1 I F ( M . E G . 4 )M=1 N=M +1 I F ( N . E Q . 4 )N=1 C C E S ( I , J , K ) = W I I (M , 3 ) * D D Y S I ( I , N , K ) - W I I ( N , 3 ) * C D Y S I ( I , M , K ) 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 ) D C P S ( I , J , K + 6 )=-DDBS( t , J , K )  565 566 567' 568 1569 1570 l\571 1572 U573 L57 4 157 5 1576 1V577 1578  1592 59 3 i 594 1595 1596 1597 1598 1^599 1 6C0 ,  82  92 91 90  16C1 16C2 1603 16C4 ..I6C5 16C6 16 0 7 16C8  1614 16 15  0  f  4  • 1-6C9 16 1 0 •i:.6 11 1612 16 13  1  CCCS(I,J,K+6 )=-CCCS I I,J,K ) D C B S ( I , J , K+3 ) = C W I ( K , M 3 ) * D Y S I ( N . I ) - D W I ( K , N , 3 ) * D Y S I ( M , I) CDCSl I , J , K +3 )= C W J ( K , M , 3 ) * C Y S I ( N , I ) -DWJ(K,N,3)*DYSI(M, I) 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 ) CCCS(I+6,J,K)=WJJ(M,3)*DDYSI(I+6,N,K)-wJJ(N,3)*CDYSI( 1+6, M , K ) C C B S ( 1 + 6 , J . K + 6 )= -DOB S(1 + 6 , J , K ) CCCS ( 1+6, J , K + 6 ) = - D D C S ( I +6, J , K ) 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) C D C S ( 1 + 6 , J ,K + 3 ) = CW J { K , f*,3 ) * C Y S I {N , I + 6 ) - C W J ( K , N , 3 ) * D Y S I ( M , I +6 ) 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) 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) CCBS(I+3,J,K+3)= CCWI(I,K,f,3)*YSl(N J-DDWI(I,K,N,3)*YSI(tf) C C C S ( I + 3 , J , K + 3 ) = D D W J ( I , K , K , 3 ) * Y SI ( N ) - D D W J ( I , K , N , 3 ) * Y S I ( M ) C C E S l 1 + 3 , J , K 4 6 )= - D D B S ( I + 3 , J , K )  ;  1,579 i;580 1581 1582 D-58 3 1*5 8 4 1*585 15 86 15 8 7 1588 1589 1590 £591  1  C C C S ( 1 + 3 , J , K + 6 ) = - C D C S < I +3 , J , K ) CCNTINUE DC 9 0 1 = 1 , 9 DC 9 1 K = l , 9 C C Y E S J ( K , I ) = ( D D B S ( I , 1 , K ) * B S ( 1 ) + C B S ( 1 ,K ) * D B S ( 1 , I ) + D D B S ( I , 2 , K ) * B S * ( 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 *-CVESJ(K)*CYESJ( I )/YESJ 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 *(2)+DCS(2,K)*CCS(2,1)+DDCS(I,3,K)*CS(3)+DCS(3,K)*DCS(3,I))/ZESJ *-CZESJ(K)*CZESJ( I )/ZESJ DC 9 2 J = l , 3 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) )/ *YESJ -DYSJ(J,K)*CYESJ( I )/YESJ DDZSJ(I,J,K)=(CCCS(I,J,K)-COZESJ(K,I)*ZSJ(J)-CZESJ(K)*DZSJ(J,I))/ * Z E S J - DZSJ(J,K)*DZESJ( I )/ZESJ CCNTINUE CCNTINUE CON T I N U E DC 1 9 1 1 = 1 , 9 DC 1 9 0 K = l , 9 DC 9 3 J = l , 3 M= J+1 I F ( M . E Q . 4 )P=1 N=N + 1 I F ( N . E Q . 4 )N= 1 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 ) * * Y S J ( N ) + D Y S I ( N , K ) * C Y S J ( N , I ) - D Y S l ( N I ) * C Y S J ( N , K ) - Y S I ( N ) *DDYSJ ( I , M ,K *)-CDYSI(I,N,K)*YSJ(M)-DYSI(N,K)*rYSJ(M,I) 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 ) * *ZSJ(N)+CYSI(N,K)*CZSJ(N,I)-DYSI(N,I)*DZSJ(y,K)-YSI(N)*DDZSJ(I,N,K) f  93 19C 191  * - C D Y S I ( I , N , K ) * Z S J ( M ) - D Y S I ( N , K ) *-CZS J ( N , I ) CONTINUE CCNTINUE CCNTINUE DC 9 4 K = l , 3 CO 9 5 CCY ISI  1=1,9 ( 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£16 16 17 '.618 1 6 19 1620 1621 1^2 2 16 2 3 |l.624 625 626 627 128 .629 1630 D6 3 1 16 3 2  *CCYSI(I,3,K) 111 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 ) * *CCYSI(I,3,K) 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) * *CCYSJ(I,3,K ) 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 ) * *CCZSJ(I,3,K) 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 S K ( I , 2 » K ) + WI I ( 3 , 1 ) * *CCYSK( I,3,K ) CCZISK(K,I)=WJJ(l,l)*CCZSK(I,l,K)+ W'JJ(2,l)*CCZSK(I,2,K)+WJJ(3,1)* * C C Z S K ( 1, 3 , K ) 95 94  CCNTINUE CCNTINUE DC 1 9 5 K = l , 3 DO 1 9 4 1 = 4 , 6 N=I-3 CCYISI(K,I)=CWI(N,1,1)*DYSI(1,K)  + D W U N , 2 , 1 ) * C Y S I ( 2 , K J+DWI(N,  1633 1,6 3 4  * 3 , 1 ) *DY S I ( 3 , K ) 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 ,  L635 1636 t'6 3 7 It 3 8 1639 1;&4C 1641 16 4 2  * 3 , 1)*UYSI(3,K) CCYISJ(K,I)=CCYISJ(K,I)+DWI(N,l,l)*DYSJ(l,K)+0WI(N,2 ,1)*DYSJ(2 * ,K ) + DWl ( N , 3 , 1 ) * C Y S J ( 2 , K ) CCZISJ(K,I)=CCZISJ(K,I)+CWJ(N,1,1)*CZSJ(1,K)4CWJ(N,2,1)*0ZSJ (2,K) + *DWJ(N,3,1)*OZSJ(3,K) CCYISK(K,I)=CCYISK(K,I)+DWI(N,1,1)*DYSK(1,K)+DWI (N * 2 «1 ) * D Y S K ( 2 , K ) + *DWI ( N , 3 , 1 ) * C Y S K { 3 , K ) 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 1644 1645; i|646 1647 1648 L649 »,65C 1;651 1652 115 3 i'654  194 195  96  1655 '656 16 5 7 1658 !£59 1 66C  =*DWJ(N , 3 , 1 ) * C Z S K ( CCNTINUE CCNTINUE DO 1 9 6 K = l , 3  3,K)  DC 1 9 7 1 = 1 , 9 C C Y I S I (K+6 , I ) = - C C Y l S I ( K , I ) C C Z I S K K + fc, I ) = - C D Z I S I ( K , I ) C C Y I S J ( K + 6,I ) =- C C Y I S J ( K , I ) C C Z I S J ( K +6,I ) = - C C Z I S J (K, I ) CCYISK(K-<6,I)=-CCYISK(K,I) C C Z I SK(K + 6,I ) = - C C Z I S K ( K , I ) M=K+ 3 CCYISI(M,I)=CWI(K,l,l)*CYSIll,I )+DWI(K,2,l)*CYSI(2,I)+DWI(K,3,l)* *DYSI(3,I) CCZISI(M,I)=CWJ(K,1,1)*DYSI(1,I)+DWJ(K,2,1)*CYSI(2,I)+0WJ(K,3,1)* *CYSI(3,I) C C Y I S J ( M , I ) = WII ( 1 , 1 ) * C C Y S J ( I , 1 , M ) + WI I ( 2 , I ) * C C Y S J ( 1 , 2 , M ) + W I I ( 3 , 1 ) * * 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) *  ;  1661 1662 166 3 1664' J :6 5 IS 66  *CYSJ(3,I) CCZISJ(M,I)=WJJ(1,1) *CCZSJ( I,3,N)+CWJ(K, *CZSJ(3,I) C C Y I S K ( M , I ) = WI I ( 1 ,1 * C C Y S K ( I , 3 , N ) + DWI ( K ,  1667 1668 1A69 167C  *CYSK(3,I) CCZISK(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) +WJJ(3,1)* * C C Z S K ( I , 3 , V) + C W J ( 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 ) * *CZSK(3,I)  (  1,671 1672 1673 1674 1675  197 196  CCNTINUE CONTINUE CC 1 9 8 K=l,3 DC 1 9 9 1 = 4 , 6 M=K+3  *DDZSJ(I,1,M)+WJJ(2,1)*CCZSJ(I,2,K)+WJJ(3,1)* 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 ) * ) * 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) * 1 ,1 ) * D Y S K ( 1,1 ) + D W l ( K , 2 , 1 ) * C Y S K ( 2 , I )+DWI (K. , 3 , 1) *  N=I-3 CCYISI(K,I)=CCWI(N,K,l,l)*YSI(l)+CCWI(N,K,2,i)*YSl(  1676 677 !l 6 7 8 6 7 9' .680' 168 1 "6-8 2! 1683 684 •85  *)*YSI<2) D C 2 I S I ( M , I ) = C C W J ( N , K , 1 , 1 ) * Y S I ( 1 ) + D C W J ( N , K , 2 , 1 ) * Y S I ( 2 ) +DD W J ( N , K , 3 , 1 • )*YSI (3) CCYISJ(M,I)=CCYISJ(M,I)+CWl(N,l,lMCYSJ(l,rM+CWI(N,2,l)  • DWJ ( N , 3 , 1 1 * D Z S J ( 2 , M ) + D D W J ( N , K , 1 , 1 1 * Z S J ( 1 ) + 0 0 V» J ( N » K , 2 ,1 ) * Z S J ( 2 ) + • C C w J ( N , K , 3 , 1 )*ZSJ (3) C C Y I S K ( M , I ) = C C Y I S K ( K , D + C W I (N , 1,1)*CYSK(I,K)+CWI (N,2,1)*DYSK (2,M) + • D W I ( N , 3 , 1 ) * C Y S K ( 2 » M ) + D D W I ( N , K , 1 , I ) * Y S K I 1}+ D0 W I ( N , K , 2 , 1 ) * Y S K ( 2 ) + *CCWI(N,K,3,1 )*YSK(3)  ^87i 6T6 a j 6891 69 C, Sv9 1  694 695 696 69 7 6 98 ; fr>99; TCCi V  701 702 AC 3 7 04 7C5 > 06 7C7 708 7C9 '7 10 .71 1 7 12 fl3 714 715 716 til 7 718  199 198  111 11C  721 722  N=K IF(I.GE.4.ANC.I.LE.6)K=I+6 IF(K.GE.4.ANC.K.LE.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 S K 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  114 113  119  V ( M , 6 , N ) = ( Z I SI * C D Z I S J ( K , I ) + D Z I SI ( I ) * Q Z I S J ( K ) - Z I S J * D D Z I S * 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 CCNTINUE CCNTINUE M=l N=3 Ml = l N 1= 3  I ( K, I )-  DC 1 1 5 I = f , N DC 1 1 6 K = M , M V ( I , 4 , K ) = C C Y S K ( I , 1 , K ) * Z S J ( 1 ) +CYSK(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) +DDYSK(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)*DYSK(2,I)+CCZSJ(I,3,K)*YSK(3)+DZSJ(3,K)*CYSK(3,I)  116 115  CONTINUE CCNTINUE  117  I F ( . E G.1)GCTC117 GOTO 1 1 8 K=7 N=9 GCT0119  n 29 .730 J?31 1/3 2 1723 1734 11735  DO 110 1 = 1 , 9 DC 1 1 1 K = l , 9 V ( I ,1 , K ) = C C E L l I , K ) CCNTINUE CCNTINUE  • 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 ) • DYISJ(I)*DYISI(K ))/YIS 1**2-2.*A{3,K)*DYISI(I J/YISI M= I  7;2 3 724 .> 2 5 7 26 727 .728  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) CCMINUE CCNTINUE  DC 1 1 3 1 = 1 , 9 DC 114 K = l , 9 VII,2,K)=(-YISI*DDYISKIK,I )-DYISH I)*DYISK(K) + YISK*DDYISI(K,I)+  tl9  i20  *DYSJ (2 , M ) +  • DWI (M , 3 , l . ) * D Y S J ( 3 , M ) + D D W I ( 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 ) CCZISJ(M,I)=CCZISJ(M,I)+DWJ(N,l,l)*DZSj(l,n+CWJ(N,2,l)*DZSJ(2,M)+  a 8*6,  b92 &S3  H2 2 J+DDWI ( N » K » 3 » 1  118  M= 1 N=3 1F(H1.EG.I1GCTC130  M 1 = 1 Nl=3 GCTC 13C  113  124  Ml=7 N 1=9 GCT0119 DC  120  1743  DC  121  174 4!  V<I,4,K)=  124  *(3)  '745  1=4,6 K=  M , M CCYSK(  +CZSJ(1,K  /46  121  CCNTINUE  747  12C  CCNTINUE  1,1, K ) * Z SJ (1) + C D Y S K { I , 2 , K ) * Z S J ( 2 ) + D D Y S K ( I ,3,K)*ZSJ )*CYSK( 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 )  '748  IF(Ml.EG.11GCTC122  749  GC T C 1 2 3  .750  122  751  Ml = 7 Nl=9  7 5 21  GCTC124  i,2 53i  123  St 5 4  L7 5 5  DC  126  1=4,6  DO  126  K=N1,N1  V(I+6,4,K)=  1756  CYSK(1,K ) * D Z S J ( l , I ) +O Y £ K ( 2 , K ) * D Z S J ( 2 , I ) + O Y S K ( 3 , K )*DZSJ ) + D D 2 S J ( I , 1 , K ) * Y 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 ) * Y S K ( 3 ) CONT INUE  *(3,1  7 57 .1*58 1759,  126  IF(H1.E0.71GCTC131 GCTC132  i"?6o:  131  .761  Ml=l Nl =3  176 2  GCTC123 132  76 3  DC  133  I = P 1 , M  .764  DO  134  K=4,6  L76 5;'  V ( 1 , 4 , K )=  CCYSK ( 1,1, 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 *(2)+DYSK(2,K)*CZSj(2 ,I)+DDYSK<I,3,K)*ZSJ(3)+CYSK(3,K)*DZSJ(3,1) V ( I , 4 , K + 6 ) = C C Z S J ( 1 , 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 * K ( 2 ) + D Z S J ( 2 , K )*CYSK( 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 )  [7,6 6 .767 1768 1769  134  CCNTINUE  1-7 7 0  133  CCNT INUE  1?71  IF(M 1 . E Q . 1 1 G C T 0 1 3 5  1772  GCT0136  U7 7 3  135  77 4  Ml = 7 N 1= 9  1775  GCT0132 136  '.77 6 1^77  CC  137  1=4,6  DO  138  K=4,6  V( I , 4 , K ) =  1778  178 1  V(I,4,K+6)=  1782  J +C Y S K ( 2 , K ) * D Z S J ( 2 , I  )+DYSK(3,K)*DZSJ(  C Z S J ( 1 , K ) * C Y S K ( 1 , I ) + D Z S J ( 2 , K ) * D Y S K ( 2 , I )+D Z S J ( 3 , K ) * D Y S K  * ( 3, I )  78 4  V(I+6,4,K+6)=CCZSJ(  85  *YSK(3)  V7 6 6  138  CCNTINUE  1787  137  CCNT INUE  .7 8 8  RETURN  1"89  EN C  IC  ,3,K)*ZSJ  * 3, I )  17 8 3  CF  1,1,K)*ZSJ(1)+CDYSK{ I ,2,K)*ZSJ(2)+DDYSK(I  V ( I + 6 , 4 , K ) = C Y 5 K ( 1 , K ) * C Z S J (1,I  17 8 C  2  CCYSKI  * (3 )  '779  F I LE  I , 1, K ) * Y S K ( 1)+0CZSJ( I,2,K  )*YSK(2)+DDZSJ(I  ,3,K)*  

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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050543/manifest

Comment

Related Items