'***************************************************************************************
'* This program calculates the deformations and stresses in a laminated material made  *
'* of n unidirectional composite layers, knowing the resulting imposed efforts.        *
'*                                                                                     *
'* The 6x6 stiffness matrix to define has the following coefficients:                  *
'*                                                                                     *
'*     A11 A12 A16 B11 B12 B16                                                         *
'*     A12 A22 A26 B12 B22 B26        ! A  B !                                         *
'*     A16 A26 A66 B16 B26 B66    or  !      !                                         *
'*     B11 B12 B16 D11 D12 D16        ! B  D !                                         *
'*     B12 B22 B26 D12 D22 D26                                                         *
'*     B16 B26 B66 D16 D26 D66                                                         *
'*                                                                                     *
'* The imposed efforts vector is (Nx,Ny,Nxy,Mx,My,Mxy)                                 *
'*                                                                                     *
'* The deformation vector is (EPs xx, EPs yy, gam xy, Kx, Ky, Kxy)                     *
'*                                                                                     *
'* ----------------------------------------------------------------------------------- *
'* SAMPLE RUN:                                                                         *
'*   The output file Compos04.lst contains (extract):                                  *
'*                                                                                     *
'*   CALCULATE DEFORMATIONS AND STRESSES                                               *
'*        IN A LAMINATED MATERIAL                                                      *
'*                                                                                     *
'*   Imposed resulting efforts:                                                        *
'*                                                                                     *
'*     NX =0.10000E+07  NY =0.50000E+06  NXY =0.25000E+06                              *
'*     MX =0.00000E+00  MY =0.00000E+00  MXY =0.00000E+00                              *
'*                                                                                     *
'*   Number of layers:            4                                                    *
'*                                                                                     *
'*     Layer  1 : angle  15.0 deg. - thickness   1.5 mm                                *
'*                                                                                     *
'*     Layer  2 : angle -30.0 deg. - thickness   1.0 mm                                *
'*                                                                                     *
'*     Layer  3 : angle -15.0 deg. - thickness   1.5 mm                                *
'*                                                                                     *
'*     Layer  4 : angle  30.0 deg. - thickness   1.0 mm                                *
'*                                                                                     *
'*   Material Parameters:                                                              *
'*                                                                                     *
'*     El =  38.0    Et =   9.0                                                        *
'*   NUlt =  0.32    Glt=   3.6                                                        *
'*                                                                                     *
'*                                                                                     *
'*   BASIC REDUCED STIFFNESS MATRIX IN GPa:                                            *
'*                                                                                     *
'*   0.38945E+02   0.29516E+01   0.00000E+00                                           *
'*   0.29516E+01   0.92237E+01   0.00000E+00                                           *
'*   0.00000E+00   0.00000E+00   0.36000E+01                                           *
'*                                                                                     *
'*   Layer #:           1                                                              *
'*                                                                                     *
'*   REDUCED STIFFNESS MATRIX IN GPa:                                                  *
'*                                                                                     *
'*   0.35212E+02   0.46931E+01   0.67316E+01                                           *
'*   0.46931E+01   0.94731E+01   0.69862E+00                                           *
'*   0.67316E+01   0.69862E+00   0.53416E+01                                           *
'* ---/---                                                                             *
'*   Stresses in each layer in main axes:                                              *
'*                                                                                     *
'*   Layer #           1                                                               *
'*                                                                                     *
'*   For z = -0.25000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.28274E+03  MPa                                                      *
'*    sigma T :  0.78080E+02  MPa                                                      *
'*    sigma LT:  0.45159E+02  MPa                                                      *
'*                                                                                     *
'*   For z = -0.10000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.28858E+03  MPa                                                      *
'*    sigma T :  0.70533E+02  MPa                                                      *
'*    sigma LT:  0.34701E+02  MPa                                                      *
'*                                                                                     *
'*   Layer #           2                                                               *
'*                                                                                     *
'*   For z = -0.10000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.86434E+02  MPa                                                      *
'*    sigma T :  0.10576E+03  MPa                                                      *
'*    sigma LT:  0.57365E+01  MPa                                                      *
'*                                                                                     *
'*   For z =  0.00000E+00 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.11193E+03  MPa                                                      *
'*    sigma T :  0.96965E+02  MPa                                                      *
'*    sigma LT:  0.83890E+01  MPa                                                      *
'*                                                                                     *
'*   Layer #           3                                                               *
'*                                                                                     *
'*   For z =  0.00000E+00 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.15147E+03  MPa                                                      *
'*    sigma T :  0.90075E+02  MPa                                                      *
'*    sigma LT:  0.21129E+02  MPa                                                      *
'*                                                                                     *
'*   For z =  0.15000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.19264E+03  MPa                                                      *
'*    sigma T :  0.76371E+02  MPa                                                      *
'*    sigma LT:  0.19346E+02  MPa                                                      *
'*                                                                                     *
'*   Layer #           4                                                               *
'*                                                                                     *
'*   For z =  0.15000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.33321E+03  MPa                                                      *
'*    sigma T :  0.51876E+02  MPa                                                      *
'*    sigma LT:  0.87729E+01  MPa                                                      *
'*                                                                                     *
'*   For z =  0.25000E+01 mm:                                                          *
'*                                                                                     *
'*    sigma L :  0.31791E+03  MPa                                                      *
'*    sigma T :  0.50191E+02  MPa                                                      *
'*    sigma LT:  0.14086E+01  MPa                                                      *
'* ----------------------------------------------------------------------------------- *
'* REFERENCE: "MATERIAUX COMPOSITES - Comportement mécanique et analyse des structures *
'*             By J.-M. Berthelot - MASSON 1996".                                      *
'*                                                                                     *
'*                                         Basic Release By J-P MOREAU - August 1997.  *
'*                                                      (www.jpmoreau.fr)              *
'***************************************************************************************
'PROGRAM COMPOS04
DEFDBL A-H, O-Z
DEFINT I-N

' MAXCOU = maximum number of layers
      NMAX = 3
      MAXCOU = 10

DIM NULT AS DOUBLE

' angle and thickness of layer
DIM TH(MAXCOU), EP(MAXCOU)

DIM H(MAXCOU)

DIM Q0(NMAX, NMAX), Q(MAXCOU, NMAX, NMAX)
DIM A(NMAX, NMAX), B(NMAX, NMAX), D(NMAX, NMAX)
DIM R(6, 6), X(6, 6), E0(6), EDEB(6), EFIN(6)
DIM N(6) AS DOUBLE

DIM NX AS DOUBLE
DIM NY AS DOUBLE
DIM NXY AS DOUBLE
DIM MX AS DOUBLE
DIM MY AS DOUBLE
DIM MXY AS DOUBLE

DIM T(NMAX, NMAX), AK(MAXCOU, NMAX), BK(MAXCOU, NMAX)
DIM AK1(MAXCOU, NMAX), BK1(MAXCOU, NMAX)
DIM TEMP(NMAX), TEMP1(NMAX)

DIM Title AS STRING

      PI = 4# * ATN(1#)

'Output formats
      F$ = "####.####"
      F1$ = "##########.##"
      F2$ = "####.#### mm"
      F3$ = "####.#### MPa"
    
' Material constants in GPa
      EL = 38#
      ET = 9#
      NULT = .32#
      GLT = 3.6#

' imposed efforts
      NX = 1000000#
      NY = 500000#
      NXY = 250000#
      MX = 0#
      MY = 0#
      MXY = 0#

      N(1) = .000000001# * NX
      N(2) = .000000001# * NY
      N(3) = .000000001# * NXY
      N(4) = .000000001# * MX
      N(5) = .000000001# * MY
      N(6) = .000000001# * MXY

      FOR I = 1 TO 6
	E0(I) = 0#
	EDEB(I) = 0#
	EFIN(I) = 0#
      NEXT I

' print data to output file
      OPEN "compos04.lst" FOR OUTPUT AS #2
      CLS
      PRINT #2,
      PRINT #2, "  CALCULATE DEFORMATIONS AND STRESSES"
      PRINT #2, "       IN A LAMINATED MATERIAL"
      PRINT #2,
      PRINT #2, "  Imposed resulting efforts:"
      PRINT #2,
      PRINT #2, "   NX="; NX; " NY = "; NY; " NXY = "; NXY
      PRINT #2, "   MX="; MX; " MY = "; MY; " MXY = "; MXY

' number of layers in laminate 
      NCOUCHES = 4

' define angles and thicknesses of layers 1-4   
      TH(1) = 15#
      EP(1) = .0015#
      TH(2) = -30#
      EP(2) = .001#
      TH(3) = -15#
      EP(3) = .0015#
      TH(4) = 30#
      EP(4) = .001#
	   
      PRINT #2,
      PRINT #2, "  Number of layers: "; NCOUCHES
      PRINT #2,
      FOR I = 1 TO NCOUCHES
	PRINT #2, "    Layer #"; I; ": angle "; TH(I); " deg. - thickness "; 1000# * EP(I); " mm"
	PRINT #2,
' put angles in radians
	TH(I) = TH(I) * PI / 180#
      NEXT I

      PRINT #2,
      PRINT #2, "  Material Parameters:"
      PRINT #2,
      PRINT #2, "  EL="; EL; "  ET="; ET
      PRINT #2, "  NULT="; NULT; "  GLT="; GLT
      PRINT #2,
     
' calculate the reduced stiffness matrix (angle = 0)
      GOSUB 500 'CALL MATNUL(Q0,NMAX)
      GOSUB 1000 'CALL CAL_Q0(Q0)
      Title = "  BASIC REDUCED STIFFNESS MATRIX IN GPa:"
      GOSUB 600 'CALL MATPR2("  BASIC REDUCED STIFFNESS MATRIX IN GPa:", Q0, NMAX)

' calculate the reduced stiffness matrices for different angles <> 0
      FOR L = 1 TO NCOUCHES
	IF TH(L) = 0# THEN
	  K = L
	  GOSUB 510 'CALL MATEGAL(K,Q,Q0)
	ELSE
	  K = L
	  TH = TH(L)
	  GOSUB 2000 'CALL CAL_Q(K,Q0,Q,TH)          
	END IF
	PRINT
	PRINT #2, "  Layer #"; L; ":"
	Title = "  REDUCED STIFFNESS MATRIX IN GPa:": K = L
	GOSUB 610 'CALL MATPR1("  REDUCED STIFFNESS MATRIX IN GPa:",K,Q,NMAX)
      NEXT L
     
' calculate A matrix
' A = sum of 1e3*EP(k)*Q(K)
      'CALL MATNUL(A,NMAX)
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  A(I, J) = 0#
	NEXT J
      NEXT I
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  FOR K = 1 TO NCOUCHES
	    A(I, J) = A(I, J) + 1000# * (EP(K) * Q(K, I, J))
	  NEXT K
	  IF ABS(A(I, J)) < .00000001# THEN A(I, J) = 0#
	NEXT J
      NEXT I
      Title = "  A MATRIX IN 1E6 N/M:"
      GOSUB 601 'CALL MATPR2('  A MATRIX IN 1E6 N/M:',A,3)
		      
' calculate B matrix
' B = sum of 1e6*(1/2)*(h(k)*h(k)-h(k-1)*h(k-1))*Q(k)
      'CALL MATNUL(B,NMAX)
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  B(I, J) = 0#
	NEXT J
      NEXT I
      HM = 0#
      FOR I = 1 TO NCOUCHES
	HM = HM + EP(I)
      NEXT I
' average height of stack
      HM = HM / 2#
' lower height layer #1
      H(1) = -HM
      FOR I = 2 TO NCOUCHES + 1
	H(I) = H(I - 1) + EP(I - 1)
      NEXT I
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  FOR K = 1 TO NCOUCHES
	    B(I, J) = B(I, J) + 1000000# * .5 * ((H(K + 1) * H(K + 1) - H(K) * H(K)) * Q(K, I, J))
	  NEXT K
	  IF ABS(B(I, J)) < .00000001# THEN B(I, J) = 0#
	NEXT J
      NEXT I
      Title = "  B MATRIX IN 1E3 N:"
      GOSUB 602 'CALL MATPR2("  B MATRIX IN 1E3 N:",B,NMAX)
    
' calculate D matrix
' D = sum of 1e9*(1/3)*(h(k)*h(k)*h(k)-h(k-1)*h(k-1)*h(k-1))*Q(k)
      'CALL MATNUL(D,NMAX)
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  D(I, J) = 0#
	NEXT J
      NEXT I
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  FOR K = 2 TO NCOUCHES + 1
	    D(I, J) = D(I, J) + 1000000000# * (1# / 3#) * ((H(K) * H(K) * H(K) - H(K - 1) * H(K - 1) * H(K - 1)) * Q(K - 1, I, J))
	  NEXT K
	  IF ABS(D(I, J)) < .00000001# THEN D(I, J) = 0#
	NEXT J
      NEXT I
      Title = "  D MATRIX IN N.DM:"
      GOSUB 603 'CALL MATPR2("  D MATRIX IN N.DM:",D,NMAX)
    
' assemble stiffness matrix of the laminate
'          ! A B !
'      R = !     !
'          ! B D !
      'CALL MATNUL(R,6)
      FOR I = 1 TO 6
	FOR J = 1 TO 6
	  R(I, J) = 0#
	NEXT J
      NEXT I
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  R(I, J) = 1000000# * A(I, J)
	NEXT J
      NEXT I
      FOR I = 1 TO NMAX
	FOR J = 4 TO NMAX + 3
	  R(I, J) = 1000# * B(I, J - 3)
	NEXT J
      NEXT I
      FOR I = 4 TO NMAX + 3
	FOR J = 1 TO NMAX
	  R(I, J) = 1000# * B(I - 3, J)
	NEXT J
      NEXT I
      FOR I = 4 TO NMAX + 3
	FOR J = 4 TO NMAX + 3
	  R(I, J) = D(I - 3, J - 3)
	NEXT J
      NEXT I
      Title = "  STIFFNESS R MATRIX:"
      GOSUB 620 'CALL MATPR2("  STIFFNESS R MATRIX:",R,6)
      PRINT #2,

' Direct inversion of R matrix, stored in R
      N = 6: M = 0
      GOSUB 3000 'CALL MATINV(6,0,R,X,DET)
      PRINT #2,
      PRINT #2, "  DETERMINANT OF R: "; DET
      FOR I = 1 TO 6
	FOR J = 1 TO 6
	  R(I, J) = 1000000000# * R(I, J)
	NEXT J
      NEXT I
      Title = "  INVERSE MATRIX OF R:"
      GOSUB 620 'CALL MATPR2("  INVERSE MATRIX OF R:", R, 6)

' Calculate deformations eps0xx, eps0yy, gam0xy,kx,ky,kxy
' in main axes, stored in vector E0(6)

      GOSUB 520 'CALL MATMUL(R,N,E0,6)
      PRINT #2,
      PRINT #2, "  Membrane deformations (in mm) and curbatures:"
      PRINT #2,
      PRINT #2, "    eps0xx: "; : PRINT #2, USING F$; 1000# * E0(1)
      PRINT #2, "    eps0yy: "; : PRINT #2, USING F$; 1000# * E0(2)
      PRINT #2, "    gam0xy: "; : PRINT #2, USING F$; 1000# * E0(3)
      PRINT #2, "    kx    : "; : PRINT #2, USING F$; E0(4)
      PRINT #2, "    ky    : "; : PRINT #2, USING F$; E0(5)
      PRINT #2, "    kxy   : "; : PRINT #2, USING F$; E0(6)

' Deformations in reference axes (x,y)
'     ! epsxx !   ! eps0xx !     ! kx  !
'     ! epsyy ! = ! eps0yy ! + z ! ky  !
'     ! gamxy !   ! gam0xy !     ! kxy !

      PRINT #2,
      PRINT #2, "  Deformations in reference axes (x,y):"
      PRINT #2, "   For z = "; 1000# * H(1); " mm:"
      FOR I = 1 TO NMAX
	EDEB(I) = E0(I) + H(1) * E0(I + 3)
      NEXT I
      PRINT #2, "    epsxx: "; : PRINT #2, USING F2$; 1000# * EDEB(1)
      PRINT #2, "    epsyy: "; : PRINT #2, USING F2$; 1000# * EDEB(2)
      PRINT #2, "    gamxy: "; : PRINT #2, USING F2$; 1000# * EDEB(3)
      PRINT #2, "   For z = "; 1000# * H(NCOUCHES + 1); " mm:"
      FOR I = 1 TO NMAX
	EFIN(I) = E0(I) + H(NCOUCHES + 1) * E0(I + 3)
      NEXT I
      PRINT #2, "    epsxx: "; : PRINT #2, USING F2$; 1000# * EFIN(1)
      PRINT #2, "    epsyy: "; : PRINT #2, USING F2$; 1000# * EFIN(2)
      PRINT #2, "    gamxy: "; : PRINT #2, USING F2$; 1000# * EFIN(3)

' Deformations in each layer in main axes of layer
      PRINT #2,
      PRINT #2, "  Deformations in each layer in main axes:"
      FOR K = 1 TO NCOUCHES
	TH = TH(K)
	GOSUB 4000 'CALL INITT(T,TH(K))
	FOR I = 1 TO NMAX
	  TEMP(I) = E0(I)
	NEXT I
	GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  AK(K, I) = TEMP1(I)
	NEXT I
	FOR I = 1 TO NMAX
	  TEMP(I) = E0(I + 3)
	NEXT I
	GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  BK(K, I) = TEMP1(I)
	NEXT I
	PRINT #2,
	PRINT #2, "  Layer #"; K; ":"
	PRINT #2, "   For z = "; 1000# * H(K); " mm:"
	FOR I = 1 TO NMAX
	  EDEB(I) = AK(K, I) + H(K) * BK(K, I)
	NEXT I
	PRINT #2, "    epsxx: "; : PRINT #2, USING F2$; 1000# * EDEB(1)
	PRINT #2, "    epsyy: "; : PRINT #2, USING F2$; 1000# * EDEB(2)
	PRINT #2, "    gamxy: "; : PRINT #2, USING F2$; 1000# * EDEB(3)

	PRINT #2, "   For z = "; 1000# * H(K + 1); " mm:"
	FOR I = 1 TO NMAX
	  EFIN(I) = AK(K, I) + H(K + 1) * BK(K, I)
	NEXT I
	PRINT #2, "    epsxx: "; : PRINT #2, USING F2$; 1000# * EFIN(1)
	PRINT #2, "    epsyy: "; : PRINT #2, USING F2$; 1000# * EFIN(2)
	PRINT #2, "    gamxy: "; : PRINT #2, USING F2$; 1000# * EFIN(3)
      NEXT K

' Stresses in each layer in reference axes

      PRINT #2,
      PRINT #2, "  Stresses in each layer in reference axes:"
      FOR K = 1 TO NCOUCHES
	FOR I = 1 TO NMAX
	  TEMP(I) = E0(I)
	NEXT I
	FOR I = 1 TO NMAX
	  FOR J = 1 TO NMAX
	    A(I, J) = Q(K, I, J)
	  NEXT J
	NEXT I
	GOSUB 531 'CALL MATMUL(A,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  AK(K, I) = TEMP1(I)
	NEXT I
	FOR I = 1 TO NMAX
	  TEMP(I) = E0(I + 3)
	NEXT I
	GOSUB 531 'CALL MATMUL(A,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  BK(K, I) = TEMP1(I)
	NEXT I
	PRINT #2,
	PRINT #2, "  Layer #"; K; ":"
	PRINT #2, "   For z = "; 1000# * H(K); " mm:"
	FOR I = 1 TO NMAX
	  EDEB(I) = 1000# * (AK(K, I) + H(K) * BK(K, I))
	NEXT I
	PRINT #2, "    sigma L : "; : PRINT #2, USING F3$; EDEB(1)
	PRINT #2, "    sigma T : "; : PRINT #2, USING F3$; EDEB(2)
	PRINT #2, "    sigma LT: "; : PRINT #2, USING F3$; EDEB(3)
	PRINT #2, "   For z = "; 1000# * H(K + 1); " mm:"
	FOR I = 1 TO NMAX
	  EFIN(I) = 1000# * (AK(K, I) + H(K + 1) * BK(K, I))
	NEXT I
	PRINT #2, "    sigma L : "; : PRINT #2, USING F3$; EFIN(1)
	PRINT #2, "    sigma T : "; : PRINT #2, USING F3$; EFIN(2)
	PRINT #2, "    sigma LT: "; : PRINT #2, USING F3$; EFIN(3)
      NEXT K

' Stresses in each layer in main axes

      PRINT #2,
      PRINT #2, "  Stresses in each layer in main axes:"
      FOR K = 1 TO NCOUCHES
	TH = TH(K)
	GOSUB 4100 'CALL INITT1(T,TH(K))
	FOR I = 1 TO NMAX
	  TEMP(I) = AK(K, I)
	NEXT I
	GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  AK1(K, I) = TEMP1(I)
	NEXT I
	FOR I = 1 TO NMAX
	  TEMP(I) = BK(K, I)
	NEXT I
	GOSUB 530 'CALL MATMUL(T,TEMP,TEMP1,3)
	FOR I = 1 TO NMAX
	  BK1(K, I) = TEMP1(I)
	NEXT I
	PRINT #2,
	PRINT #2, "  Layer #"; K; ":"
	PRINT #2, "   For z = "; 1000# * H(K); " mm:"
	FOR I = 1 TO NMAX
	  EDEB(I) = 1000# * (AK1(K, I) + H(K) * BK1(K, I))
	NEXT I
	PRINT #2, "    sigma L : "; : PRINT #2, USING F3$; EDEB(1)
	PRINT #2, "    sigma T : "; : PRINT #2, USING F3$; EDEB(2)
	PRINT #2, "    sigma LT: "; : PRINT #2, USING F3$; EDEB(3)
	PRINT #2, "   For z = "; 1000# * H(K + 1); " mm:"
	FOR I = 1 TO NMAX
	  EFIN(I) = 1000# * (AK1(K, I) + H(K + 1) * BK1(K, I))
	NEXT I
	PRINT #2, "    sigma L : "; : PRINT #2, USING F3$; EFIN(1)
	PRINT #2, "    sigma T : "; : PRINT #2, USING F3$; EFIN(2)
	PRINT #2, "    sigma LT: "; : PRINT #2, USING F3$; EFIN(3)
      NEXT K
     
      PRINT #2,
      CLOSE #2
      PRINT "  END OF PROGRAM."

END  'of main program

'UTILITY ROUTINES

' set matrix Q0(N,N) to zero
500  'SUBROUTINE MATNUL(Q0,NMAX)
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  Q0(I, J) = 0#
	NEXT J
      NEXT I
RETURN

' A(K,3,3) = B(3,3)
510  'SUBROUTINE MATEGAL(K,Q,Q0)
      FOR II = 1 TO NMAX
	FOR JJ = 1 TO NMAX
	  Q(K, II, JJ) = Q0(II, JJ)
	NEXT JJ
      NEXT II
RETURN

' multiply a 6x6 matrix by a 6x1 vector
520  'SUBROUTINE MATMUL(R, N, E0, 6)                                             
      FOR I = 1 TO 6
	SUM = 0#
	FOR K = 1 TO 6
	  SUM = SUM + R(I, K) * N(K)
	NEXT K
	E0(I) = SUM
      NEXT I
RETURN

' multiply a 3x3 matrix by a 3x1 vector
530  'SUBROUTINE MATMUL(T, TEMP, TEMP1, 3)                                             
      FOR I = 1 TO 3
	SUM = 0#
	FOR K1 = 1 TO 3
	  SUM = SUM + T(I, K1) * TEMP(K1)
	NEXT K1
	TEMP1(I) = SUM
      NEXT I
RETURN

531  'SUBROUTINE MATMUL(A, TEMP, TEMP1, 3)
      FOR I = 1 TO 3
	SUM = 0#
	FOR K1 = 1 TO 3
	  SUM = SUM + A(I, K1) * TEMP(K1)
	NEXT K1
	TEMP1(I) = SUM
      NEXT I
RETURN

' print Q0(3,3) matrix
600  'SUBROUTINE MATPR2(TITLE, Q0, NMAX)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  PRINT #2, USING F$; Q0(I, J);
	NEXT J
	PRINT #2,
      NEXT I
RETURN

' print A(3,3) matrix
601  'SUBROUTINE MATPR2(TITLE, A, NMAX)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  PRINT #2, USING F$; A(I, J);
	NEXT J
	PRINT #2,
      NEXT I
RETURN

' print B(3,3) matrix
602  'SUBROUTINE MATPR2(TITLE, B, NMAX)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  PRINT #2, USING F$; B(I, J);
	NEXT J
	PRINT #2,
      NEXT I
RETURN

' print D(3,3) matrix
603  'SUBROUTINE MATPR2(TITLE, B, NMAX)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR I = 1 TO NMAX
	FOR J = 1 TO NMAX
	  PRINT #2, USING F$; D(I, J);
	NEXT J
	PRINT #2,
      NEXT I
RETURN

' print Q(K,3,3) matrix
610  'SUBROUTINE MATPR1(TITLE, K, Q, NMAX)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR II = 1 TO NMAX
	FOR JJ = 1 TO NMAX
	  PRINT #2, USING F$; Q(K, II, JJ);
	NEXT JJ
	PRINT #2,
      NEXT II
RETURN

' print R(6,6) matrix
620  'SUBROUTINE MATPR2(TITLE, R, 6)
      PRINT #2,
      PRINT #2, Title
      PRINT #2,
      FOR I = 1 TO 6
	FOR J = 1 TO 6
	  PRINT #2, USING F1$; R(I, J);
	NEXT J
	PRINT #2,
      NEXT I
RETURN

'  calculate cos(t) power 4
700  'REAL*8 FUNCTION DCOS4(T)       
      A = COS(TH)
      DCOS4 = A * A * A * A
RETURN
      
'  calculate sin(t) power 4
701  'REAL*8 FUNCTION DSIN4(T)       
      A = SIN(TH)
      DSIN4 = A * A * A * A
RETURN
  
'  calculate the reduced stiffness matrix for theta=0
1000 'SUBROUTINE CAL_Q0(Q0)
      Q0(1, 1) = EL / (1# - ((ET / EL) * NULT * NULT))
      Q0(2, 2) = (ET / EL) * Q0(1, 1)
      Q0(1, 2) = NULT * Q0(2, 2)
      Q0(2, 1) = Q0(1, 2)
      Q0(3, 3) = GLT
RETURN
  
'  calculate the reduced stiffness Q matrix for theta <> 0
'  with the condition that the Q matrix for theta=0 is available
'  called here Q2
2000 'SUBROUTINE CAL_Q(K, Q0, Q, TH)
      GOSUB 700: GOSUB 701 'calculate dcos4, dsin4
      TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH)
      Q(K, 1, 1) = Q0(1, 1) * DCOS4 + Q0(2, 2) * DSIN4 + 2# * (Q0(1, 2) + 2# * Q0(3, 3)) * TEMP
      Q(K, 1, 2) = (Q0(1, 1) + Q0(2, 2) - 4# * Q0(3, 3)) * TEMP + Q0(1, 2) * (DCOS4 + DSIN4)
      Q(K, 2, 1) = Q(K, 1, 2)
      TEMP = SIN(TH) * COS(TH) * COS(TH) * COS(TH)
      Q(K, 1, 3) = (Q0(1, 1) - Q0(1, 2) - 2# * Q0(3, 3)) * TEMP
      TEMP = SIN(TH) * SIN(TH) * SIN(TH) * COS(TH)
      Q(K, 1, 3) = Q(K, 1, 3) + (Q0(1, 2) - Q0(2, 2) + 2# * Q0(3, 3)) * TEMP
      Q(K, 3, 1) = Q(K, 1, 3)
      TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH)
      Q(K, 2, 2) = Q0(1, 1) * DSIN4 + 2# * (Q0(1, 2) + 2# * Q0(3, 3)) * TEMP + Q0(2, 2) * DCOS4
      TEMP = SIN(TH) * SIN(TH) * SIN(TH) * COS(TH)
      Q(K, 2, 3) = (Q0(1, 1) - Q0(1, 2) - 2# * Q0(3, 3)) * TEMP
      TEMP = COS(TH) * COS(TH) * COS(TH) * SIN(TH)
      Q(K, 2, 3) = Q(K, 2, 3) + (Q0(1, 2) - Q0(2, 2) + 2# * Q0(3, 3)) * TEMP
      Q(K, 3, 2) = Q(K, 2, 3)
      TEMP = SIN(TH) * SIN(TH) * COS(TH) * COS(TH)
      Q(K, 3, 3) = (Q0(1, 1) + Q0(2, 2) - 2# * (Q0(1, 2) + Q0(3, 3))) * TEMP + Q0(3, 3) * (DSIN4 + DCOS4)
RETURN

3000 'SUB MATINV(N,M,R,X,DET)
'*******************************************
'*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
'*  with Gauss-Jordan method using full    *
'*  pivoting at each step. During the pro- *
'*  cess, original A and B matrices are    *
'*  destroyed to spare storage location.   *
'* --------------------------------------- *
'* INPUTS:    A   MATRIX N*N               *
'*            B   MATRIX N*M               *
'* --------------------------------------- *
'* OUTPUS:    A   INVERSE OF A N*N         *
'*            DET  DETERMINANT OF A        *
'*            B   SOLUTION MATRIX N*M      *
'* --------------------------------------- *
'* NOTA - If M=0 inversion of A matrix     *
'*        only.                            *
'*******************************************                                 
DIM PC(N) AS INTEGER
DIM PL(N) AS INTEGER
DIM CS(N)

      EPSMACH = .0000000001#
						      
' INITIALIZATIONS:                                                          
      DET = 1#
      FOR I = 1 TO N
	PC(I) = 0
	PL(I) = 0
	CS(I) = 0#
      NEXT I
	  
' Main loop                                                                         
      FOR K = 1 TO N
' Searching greatest pivot:                                           
	PV = R(K, K)
	IK = K
	JK = K
	PAV = ABS(PV)
	FOR I = K TO N
	  FOR J = K TO N
	    IF ABS(R(I, J)) > PAV THEN
	      PV = R(I, J)
	      PAV = ABS(PV)
	      IK = I
	      JK = J
	    END IF
	  NEXT J
	NEXT I
									       
' Search terminated, Pivot is PV at location I=IK, J=JK.                 
									       
' Save Pivot location:                                     
	PC(K) = JK
	PL(K) = IK
									       
'Determinant DET is actualised
'If DET=0, ERROR MESSAGE and STOP
'Machine dependant EPSMACH equals here 1D-10                                      
									       
	IF IK <> K THEN DET = -DET
	IF JK <> K THEN DET = -DET
	DET = DET * PV
	IF ABS(DET) < EPSMACH THEN
' ERROR MESSAGE AND STOP                                                  
	  PRINT "  STOP IN MATINV: DETERMINANT TOO SMALL!"
	  END
	END IF
									       
' Locate Pivot in K,K:                                           
	IF IK <> K THEN
	  FOR I = 1 TO N
' EXCHANGE LINES IK and K of R:
	    TT = R(IK, I)
	    R(IK, I) = R(K, I)
	    R(K, I) = TT
	  NEXT I
	END IF
      IF M <> 0 THEN
	FOR I = 1 TO M
	  TT = X(IK, I)
	  X(IK, I) = X(K, I)
	  X(K, I) = TT
	NEXT I
      END IF
' Pivot is located at right line.                                             
	IF JK <> K THEN
	  FOR I = 1 TO N
' EXCHANGE COLUMNS JK and K of R :
	    TT = R(I, JK)
	    R(I, JK) = R(I, K)
	    R(I, K) = TT
	  NEXT I
	END IF
' Pivot is located at right column.                                           
' Pivot is now located at K,K.                                            
									       
' STORE COLUMN K OF R IN VECTOR CS
' ET PUT COLUMN TO ZERO                                          
	FOR I = 1 TO N
	  CS(I) = R(I, K)
	  R(I, K) = 0#
	NEXT I
									       
	CS(K) = 0#
	R(K, K) = 1#
' MODIFY LINE K OF R:
	IF ABS(PV) < EPSMACH THEN
	  PRINT #2, "  STOP MATINV: PIVOT TOO SMALL!"
	  END
	END IF
	FOR I = 1 TO N
	  R(K, I) = R(K, I) / PV
	NEXT I
	IF M <> 0 THEN
	  FOR I = 1 TO M
	    X(K, I) = X(K, I) / PV
	  NEXT I
	END IF
' MODIFY OTHER LINES OF R:
	FOR J = 1 TO N
	  IF J = K THEN GOTO 20
	  FOR I = 1 TO N
' MODIFY LINE J OF R :
	    R(J, I) = R(J, I) - CS(J) * R(K, I)
	  NEXT I
	  IF M <> 0 THEN
	    FOR I = 1 TO M
	      X(J, I) = X(J, I) - CS(J) * X(K, I)
	    NEXT I
	  END IF
20      NEXT J
' LINE K OVER.                                             
      NEXT K 'End K loop                                                           
									       
' MATRIX IS NOW INVERSED - MATRIX MUST BE REARRANGED                      
									       
' EXCHANGE LINES                                                         
      FOR I = N TO 1 STEP -1
	IK = PC(I)
	IF IK = I THEN GOTO 30
' EXCHANGE LINES I and PC(I) of R:
	FOR J = 1 TO N
	  TT = R(I, J)
	  R(I, J) = R(IK, J)
	  R(IK, J) = TT
	NEXT J
	IF M <> 0 THEN
	  FOR J = 1 TO M
	    TT = X(I, J)
	    X(I, J) = X(IK, J)
	    X(IK, J) = TT
	  NEXT J
	END IF
' NO MORE EXCHANGE NECESSARY, GO TO NEXT LINE                                                   
30    NEXT I
									      
' EXCHANGE COLUMNS                                                       
									       
      FOR J = N TO 1 STEP -1
	JK = PL(J)
	IF JK = J THEN GOTO 40
' EXCHANGE COLUMNS J and PL(J) of R:
	FOR I = 1 TO N
	  TT = R(I, J)
	  R(I, J) = R(I, JK)
	  R(I, JK) = TT
	NEXT I
' NO MORE EXCHANGE NECESSARY, GO TO NEXT COLUMN
40    NEXT J
' REARRANGEMENT TERMINATED.                                                     
RETURN 'MATINV                                                                 


' Matrix to change axes (x,y) for deformations

4000 'SUB INITT(T,TH)
      T(1, 1) = COS(TH) * COS(TH)
      T(1, 2) = SIN(TH) * SIN(TH)
      T(1, 3) = SIN(TH) * COS(TH)
      T(2, 1) = SIN(TH) * SIN(TH)
      T(2, 2) = COS(TH) * COS(TH)
      T(2, 3) = -SIN(TH) * COS(TH)
      T(3, 1) = -2# * SIN(TH) * COS(TH)
      T(3, 2) = 2# * SIN(TH) * COS(TH)
      T(3, 3) = COS(TH) * COS(TH) - SIN(TH) * SIN(TH)
RETURN

' Matrix to change axes (x,y) for stresses

4100 'SUB INITT1(T,TH)
      T(1, 1) = COS(TH) * COS(TH)
      T(1, 2) = SIN(TH) * SIN(TH)
      T(1, 3) = 2# * SIN(TH) * COS(TH)
      T(2, 1) = SIN(TH) * SIN(TH)
      T(2, 2) = COS(TH) * COS(TH)
      T(2, 3) = -2# * SIN(TH) * COS(TH)
      T(3, 1) = -SIN(TH) * COS(TH)
      T(3, 2) = SIN(TH) * COS(TH)
      T(3, 3) = COS(TH) * COS(TH) - SIN(TH) * SIN(TH)
RETURN

' End of file Compos04.bas