/**************************************************************************************
* 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 =   1.0000e+06  NY =   5.0000e+05  NXY =   2.5000e+05                        *
*     MX =   0.00000000  MY =   0.00000000  MXY =   0.00000000                        *
*                                                                                     *
*  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:                                             *
*    38.9445     2.9516     0.0000                                                    *
*     2.9516     9.2237     0.0000                                                    *
*     0.0000     0.0000     3.6000                                                    *
*                                                                                     *
*  Layer #:1                                                                          *
*                                                                                     *
*  REDUCED STIFFNESS MATRIX IN GPa:                                                   *
*    35.2120     4.6931     6.7316                                                    *
*     4.6931     9.4731     0.6986                                                    *
*     6.7316     0.6986     5.3416                                                    *
*                                                                                     *
*  Layer #:2                                                                          *
*                                                                                     *
*  REDUCED STIFFNESS MATRIX IN GPa:                                                   *
*    26.2896     8.1763    -9.4512                                                    *
*     8.1763    11.4292    -3.4183                                                    *
*    -9.4512    -3.4183     8.8247                                                    *
*                                                                                     *
*  Layer #:3                                                                          *
*                                                                                     *
*  REDUCED STIFFNESS MATRIX IN GPa:                                                   *
*    35.2120     4.6931    -6.7316                                                    *
*     4.6931     9.4731    -0.6986                                                    *
*    -6.7316    -0.6986     5.3416                                                    *
*  --/--                                                                              *
*  STIFFNESS R MATRIX:                                                                * 
*  1.5822e+008  3.0432e+007  0.0000e+000 -1.3384e+004  5.2247e+003 -1.6154e+003       *
*  3.0432e+007  5.1278e+007  0.0000e+000  5.2247e+003  2.9342e+003  5.9258e+003       * 
*  0.0000e+000  0.0000e+000  3.3674e+007 -1.6154e+003  5.9258e+003  5.2247e+003       *
* -1.3384e+004  5.2247e+003 -1.6154e+003  3.2738e+002  6.4271e+001  6.0686e+001       *
*  5.2247e+003  2.9342e+003  5.9258e+003  6.4271e+001  1.0732e+002  1.5438e+001       *
* -1.6154e+003  5.9258e+003  5.2247e+003  6.0686e+001  1.5438e+001  7.1025e+001       *
*                                                                                     *
*  DETERMINANT OF R:  4.2688e+029                                                     *
*                                                                                     *
*  INVERSE MATRIX OF R:                                                               *
*  7.2071e+000 -4.3217e+000  6.9155e-002  4.1544e+002 -5.2541e+002  2.7865e+002       *
* -4.3217e+000  2.2297e+001  2.7923e-001 -1.8738e+002 -4.2073e+001 -1.8099e+003       *
*  6.9155e-002  2.7923e-001  3.0508e+001  1.0321e+003 -1.9209e+003 -2.7302e+003       *
*  4.1544e+002 -1.8738e+002  1.0321e+003  4.0518e+006 -2.0576e+006 -3.0655e+006       *
* -5.2541e+002 -4.2073e+001 -1.9209e+003 -2.0576e+006  1.0747e+007 -4.4512e+005       *
*  2.7865e+002 -1.8099e+003 -2.7302e+003 -3.0655e+006 -4.4512e+005  1.7154e+007       *
*  --/--                                                                              *
*  Stresses in each layer in main axes:                                               *
*                                                                                     *
*  Layer #1                                                                           *
*   For z =     -2.50000 mm:                                                          *
*    sigma L :    282.73562  MPa                                                      *
*    sigma T :     78.07995  MPa                                                      *
*    sigma LT:     45.15936  MPa                                                      *
*   For z =     -1.00000 mm:                                                          *
*    sigma L :    288.58230  MPa                                                      *
*    sigma T :     70.53304  MPa                                                      *
*    sigma LT:     34.70107  MPa                                                      *
*                                                                                     *
*  Layer #2                                                                           *
*   For z =     -1.00000 mm:                                                          *
*    sigma L :     86.43422  MPa                                                      *
*    sigma T :    105.75930  MPa                                                      *
*    sigma LT:      5.73649  MPa                                                      *
*   For z =      0.00000 mm:                                                          *
*    sigma L :    111.92641  MPa                                                      *
*    sigma T :     96.96499  MPa                                                      *
*    sigma LT:      8.38896  MPa                                                      *
*                                                                                     *
*  Layer #3                                                                           *
*   For z =      0.00000 mm:                                                          *
*    sigma L :    151.46586  MPa                                                      *
*    sigma T :     90.07486  MPa                                                      *
*    sigma LT:     21.12949  MPa                                                      *
*   For z =      1.50000 mm:                                                          *
*    sigma L :    192.64459  MPa                                                      *
*    sigma T :     76.37100  MPa                                                      *
*    sigma LT:     19.34600  MPa                                                      *
*                                                                                     *
*  Layer #4                                                                           *
*   For z =      1.50000 mm:                                                          *
*    sigma L :    333.21160  MPa                                                      *
*    sigma T :     51.87584  MPa                                                      *
*    sigma LT:      8.77294  MPa                                                      *
*   For z =      2.50000 mm:                                                          *
*    sigma L :    317.90586  MPa                                                      *
*    sigma T :     50.19097  MPa                                                      *
*    sigma LT:      1.40860  MPa                                                      *
*                                                                                     *
* ----------------------------------------------------------------------------------- *
* REFERENCE: "MATERIAUX COMPOSITES - Comportement mécanique et analyse des structures *
*             By J.-M. Berthelot - MASSON 1996".                                      *
*                                                                                     *
*                                           C++ Release By J-P MOREAU - August 1997.  *
*                                                      (www.jpmoreau.fr)              *
**************************************************************************************/
#include <stdio.h>
#include <math.h>

// MAXCOU = maximum number of layers
#define  NMAX      3
#define  MAXCOU   10
#define  MACH_EPS  2e-16

//index 0 not used except for inversion of R
typedef double MAT[NMAX+1][NMAX+1];
typedef double MAT1[MAXCOU+1][NMAX+1][NMAX+1];
typedef double MAT2[MAXCOU+1][NMAX+1];
typedef double MAT6[7][7];

double  EL,ET,HM,NULT,GLT,PI;

// TH, EP = angle and thickness of layers 
double TH[MAXCOU+1], EP[MAXCOU+1], H[MAXCOU+1];

MAT    A,B,D, Q0, T;  
MAT1   Q;
MAT6   R, X; 
double N[7], E0[7], EDEB[7], EFIN[7];
double DET,NX,NY,NXY,MX,MY,MXY;
MAT2   AK, AK1, BK, BK1;

double TEMP[NMAX], TEMP1[NMAX];
int    I,J,K, NCOUCHES;

FILE   *fp;
     

// set a matrix (N,N) to zero
    void MATNUL(MAT A,int N) {   
      int I,J;
      for (I=1; I<=N; I++)
        for (J=1; J<=N; J++)
          A[I][J]=0.0;
    }     
     
// A(K,3,3) = B(3,3)
    void MATEGAL(int K, int N, MAT1 A, MAT B) {
      int I,J;
      for (I=1; I<=N; I++)
        for (J=1; J<=N; J++)
          A[K][I][J]=B[I][J];
    }

// multiply a 3x3 matrix by a 3x1 vector  
    void MATMUL(MAT A, double *B, double *C,int N) {
      double SUM; int I,K;
      for (I=1; I<=N; I++) {
        SUM= 0.0;
        for (K=1; K<=N; K++)
          SUM += A[I][K]*B[K];
        C[I]=SUM;
      }
    }

// multiply a 6x6 matrix by a 6x1 vector  
    void MATMUL6(MAT6 A, double *B, double *C,int N) {
      double SUM; int I,K;
      for (I=1; I<=N; I++) {
        SUM= 0.0;
        for (K=1; K<=N; K++)
          SUM += A[I][K]*B[K];
        C[I]=SUM;
      }
    }

// print a NxN matrix           
    void MATPR2(char *title, int N, MAT A) {
      int I,J; 
      fprintf(fp,"\n%s\n", title);
      for (I=1; I<=N; I++) { 
        for (J=1; J<=N; J++) fprintf(fp," %10.4f", A[I][J]);
        fprintf(fp,"\n");
      }
    }

    void MATPR1(char *title, int K, int N, MAT1 A) {
      int I,J; 
      fprintf(fp,"\n%s\n", title);
      for (I=1; I<=N; I++) { 
        for (J=1; J<=N; J++) fprintf(fp," %10.4f", A[K][I][J]);
        fprintf(fp,"\n");
      }
    }

    void MATPR6(char *title, int N, MAT6 A) {
      int I,J; 
      fprintf(fp,"\n%s\n", title);
      for (I=1; I<=N; I++) { 
        for (J=1; J<=N; J++) fprintf(fp," %12.4e", A[I][J]);
        fprintf(fp,"\n");
      }
    }

    double dcos4(double t) {
      double a;
      a=cos(t);
      return (a*a*a*a);
    }

    double dsin4(double t) {
      double a;
      a=sin(t);
      return (a*a*a*a);
    }
  
//  calculate the reduced stiffness matrix for theta=0
    void CAL_Q0(MAT Q) {
      Q[1][1] = EL/(1.0-((ET/EL)*NULT*NULT));
      Q[2][2] = (ET/EL)*Q[1][1];
      Q[1][2] = NULT*Q[2][2];
      Q[2][1] = Q[1][2];
      Q[3][3] = GLT;
    }

// calculate the reduced stiffness Q matrix for theta <> 0 
// with the condition that the Q matrix for theta=0 is available
// called here Q2
    void CAL_Q(int K, MAT Q2, MAT1 Q, double TH) {  
      double TEMP;
      TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH);
      Q[K][1][1]= Q2[1][1]*dcos4(TH)+Q2[2][2]*dsin4(TH)+2.0*(Q2[1][2]+2.0*Q2[3][3])*TEMP;
      Q[K][1][2]= (Q2[1][1]+Q2[2][2]-4.0*Q2[3][3])*TEMP+Q2[1][2]*(dcos4(TH)+dsin4(TH));
      Q[K][2][1]=Q[K][1][2];
      TEMP = sin(TH)*cos(TH)*cos(TH)*cos(TH);
      Q[K][1][3]=(Q2[1][1]-Q2[1][2]-2.0*Q2[3][3])*TEMP;
      TEMP = sin(TH)*sin(TH)*sin(TH)*cos(TH);
      Q[K][1][3]=Q[K][1][3]+(Q2[1][2]-Q2[2][2]+2.0*Q2[3][3])*TEMP;
      Q[K][3][1]=Q[K][1][3];
      TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH);
      Q[K][2][2]= Q2[1][1]*dsin4(TH)+2.0*(Q2[1][2]+2.0*Q2[3][3])*TEMP+Q2[2][2]*dcos4(TH);
      TEMP = sin(TH)*sin(TH)*sin(TH)*cos(TH);
      Q[K][2][3]=(Q2[1][1]-Q2[1][2]-2.0*Q2[3][3])*TEMP;
      TEMP = cos(TH)*cos(TH)*cos(TH)*sin(TH);
      Q[K][2][3]=Q[K][2][3]+(Q2[1][2]-Q2[2][2]+2.0*Q2[3][3])*TEMP;
      Q[K][3][2]=Q[K][2][3];
      TEMP = sin(TH)*sin(TH)*cos(TH)*cos(TH);
      Q[K][3][3]= (Q2[1][1]+Q2[2][2]-2.0*(Q2[1][2]+Q2[3][3]))*TEMP+Q2[3][3]*(dsin4(TH)+dcos4(TH));
    }

/******************************************                                     
*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
*  with Gauss-Jordan method using full    *
*  pivoting at each step. During the pro- *
*  cess, original AA and BB matrices are  *
*  destroyed to spare storage location.   *
* --------------------------------------- *
* INPUTS:    AA   MATRIX N*N              *                                     
*            BB   MATRIX N*M              *                                     
* --------------------------------------- *                                     
* OUTPUS:    AA   INVERSE OF AA N*N       *                                     
*            DET  DETERMINANT OF AA       *                                     
*            BB   SOLUTION MATRIX N*M     *                                     
* --------------------------------------- *                                     
* NOTA - If M=0 inversion of AA matrix    *
*        only.                            *                                     
******************************************/
void MATINV(int N,int M,MAT6 AA,MAT6 BB,double *DET)   { 	
      double PC[7],PL[7],CS[7];                                                  
      double PV,PAV,temp,TT;   
      int I,IK,J,JK,K;
                                                         
// Initializations :                                                             
      *DET= 1.0; 
      for (I=0; I<N; I++)  {                                                                
        PC[I]= 0.0;                                                                
        PL[I]= 0.0;                                                                
        CS[I]= 0.0;                                                                
      }              
// Main loop                                                                         
      for (K=0; K<N; K++) {                                                                  
// Searching greatest pivot :                                               
        PV=AA[K][K];                                                              
        IK=K;                                                                    
        JK=K;                                                                    
        PAV= fabs(PV);                                                            
        for (I=K; I<N; I++)                                                                
          for (J=K; J<N; J++)  {     
            temp = fabs(AA[I][J]);                                                        
            if (temp > PAV) {                                      
              PV=AA[I][J];                                                        
              PAV= fabs(PV);                                                      
              IK=I;                                                              
              JK=J;                                                              
            }                                                               
          }                                                                 
                                                                               
// Search terminated, the pivot is in location I=IK, J=JK.                     
// Memorizing pivot location:  		  

	PC[K]=JK;                                                                
    PL[K]=IK;                                                                
                                                                               
// DETERMINANT DET is actualised                                              
// If DET=0, ERROR MESSAGE and STOP                                           
// Machine dependant EPSMACH equals here 1e-20
                                                                               
        if (IK!=K) *DET=-*DET;                                                   
        if (JK!=K) *DET=-*DET;                                                   
        *DET=*DET*PV;  
        temp= fabs(*DET);                                                            
        if (temp < MACH_EPS) {                                          
          fprintf(fp,"\n  The determinant equals ZERO !!!\n");                                                              
          return;                                                                  
        }                                                                   
                                                                               
// POSITIONNING PIVOT IN K,K:		
        if(IK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
// EXCHANGE LINES IK and K of matrix AA:			
            TT=AA[IK][I];                                                         
            AA[IK][I]=AA[K][I];                                                    
            AA[K][I]=TT;                                                          
          }                                                                 
                                                                           
        if(M!=0)                                                           
          for (I=0; I<M; I++) {                                                               
            TT=BB[IK][I];                                                           
            BB[IK][I]=BB[K][I];                                                      
            BB[K][I]=TT;                                                            
          }                                                                   
                                                                           
// PIVOT is at correct line		  
        if(JK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
// EXCHANGE COLUMNS JK and K of matrix AA:
            TT=AA[I][JK];                                                         
            AA[I][JK]=AA[I][K];                                                    
            AA[I][K]=TT;                                                          
          }                                                                 
                                                                           
// The PIVOT is at correct column.                                              
// and is located in K,K.                                                   
                                                                               
// Column K of matrix AA is stored in CS vector                             
// then column K is set to zero. 		  
        for (I=0; I<N; I++) {                                                                
          CS[I]=AA[I][K];                                                         
          AA[I][K]= 0.0;                                                          
        }                                                                   
                                                                               
        CS[K]= 0.0;                                                                
        AA[K][K]= 1.0;                                                              
// Line K of matrix AA is modified:		
        temp= fabs(PV);                                          
        if(temp < MACH_EPS) {                                            
          fprintf(fp,"\n  PIVOT TOO SMALL - STOP\n");                               
          return;                                                                  
        }                                                                   
        for (I=0; I<N; I++)                                                                
          AA[K][I]=AA[K][I]/PV;                                                    
                                                                           
        if (M!=0)                                                         
          for (I=0; I<M; I++)                                                             
            BB[K][I]=BB[K][I]/PV;                                                  
                                                                           
// Other lines of matrix AA are modified:		  
        for (J=0; J<N; J++) {                                                                
          if (J==K) J++;                                                  
          for (I=0; I<N; I++)                                                              
// Line J of matrix AA is modified:
            AA[J][I]=AA[J][I]-CS[J]*AA[K][I];
          if (M!=0)                                                       
            for (I=0; I<M; I++)                                                          
              BB[J][I]=BB[J][I]-CS[J]*BB[K][I];                                     
        }                                                                   
// Line K is ready
      } //End of K loop

// MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA
                                                                               
// EXCHANGE LINES                                                                        
      for (I=N-1; I>=0; I--) {                                                               
        IK=(int) PC[I];                                                                
        if (IK==I) goto fin1;                                                   
// EXCHANGE LINES I and PC(I) of matrix AA:
        for (J=0; J<N; J++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[IK][J];                                                      
          AA[IK][J]=TT;                                                           
        }
        if (M!=0)                                                         
          for (J=0; J<M; J++) {                                                             
            TT=BB[I][J];                                                          
            BB[I][J]=BB[IK][J];                                                    
            BB[IK][J]=TT;                                                         
          }                                                                 
// NO MORE EXCHANGE is NECESSARY                                                      
// Go to next line                                                 
        fin1: ;     			  
      } // for i                                                                     
                                                                               
// EXCHANGE COLUMNS:  	  
      for (J=N-1; J>=0; J--) {                                                                         
        JK=(int) PL[J];                                                                
        if (JK==J) goto fin2;                                                   
// EXCHANGE COLUMNS J ET PL(J) of matrix AA:
        for (I=0; I<N; I++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[I][JK];                                                      
          AA[I][JK]=TT;                                                           
        } 
        fin2: ;                                                                  
// NO MORE EXCHANGE is NECESSARY                                                      
// Go to next column.
      }                                                                     
// REARRANGEMENT TERMINATED.                                                        
      return;  
  } //Matinv()             


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

    void INITT(MAT T, double 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.0*sin(TH)*cos(TH);
      T[3][2]=2.0*sin(TH)*cos(TH);
      T[3][3]=cos(TH)*cos(TH)-sin(TH)*sin(TH);
    }

// Matrix to change axes (x][y) for stresses

    void INITT1(MAT T, double TH) {
      T[1][1]=cos(TH)*cos(TH);
      T[1][2]=sin(TH)*sin(TH);
      T[1][3]=2.0*sin(TH)*cos(TH);
      T[2][1]=sin(TH)*sin(TH);
      T[2][2]=cos(TH)*cos(TH);
      T[2][3]=-2.0*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);
    }


void main()  {

      PI = 4.0*atan(1.0);

// Material constants in GPa

      EL = 38.0;
      ET =  9.0;
      NULT =  0.32;
      GLT = 3.6;

// imposed efforts
      NX=1e6;
      NY=0.5e6;
      NXY=0.25e6;
      MX=0.0;
      MY=0.0;
      MXY=0.0;

      N[1]=1e-9*NX;
      N[2]=1e-9*NY;
      N[3]=1e-9*NXY;
      N[4]=1e-9*MX;
      N[5]=1e-9*MY;
      N[6]=1e-9*MXY;

      for (I=1; I<7; I++) {
        E0[I]=0.0;
        EDEB[I]=0.0;
        EFIN[I]=0.0;
      }

// print data to output file
      fp=fopen("compos04.lst","w");
      fprintf(fp,"\n  CALCULATE DEFORMATIONS AND STRESSES\n");
      fprintf(fp,"       IN A LAMINATED MATERIAL\n\n");    

      fprintf(fp,"  Imposed resulting efforts:\n");
      fprintf(fp,"     NX = %12.4e", NX);
      fprintf(fp,"  NY = %12.4e", NY);
      fprintf(fp,"  NXY = %12.4e\n", NXY);
      fprintf(fp,"     MX = %12.4e", MX);
      fprintf(fp,"  MY = %12.4e", MY);
      fprintf(fp,"  MXY = %12.4e\n", MXY);

// number of layers in laminate    
      NCOUCHES=4;

// define angles and thicknesses of layers 1-4      
      TH[1]=15.0;
      EP[1]=1.5e-3;
      TH[2]= -30.0;
      EP[2]=1e-3;
      TH[3]=-15.0;
      EP[3]=1.5e-3;
      TH[4]= 30.0;
      EP[4]=1e-3;
           
      fprintf(fp,"\n  Number of layers: %d\n\n", NCOUCHES);

      for (I=1; I<=NCOUCHES; I++) {
		fprintf(fp,"   Layer #%d: angle %6.1f deg. - thickness %6.1f mm\n\n", I, TH[I], 1000*EP[I]);
        //put angles in radians
        TH[I] = TH[I]*PI/180.0;
      }

      fprintf(fp,"\n  Material Parameters:\n");
      fprintf(fp,"     El = %6.1f    Et = %6.1f\n", EL, ET);
      fprintf(fp,"   NUlt = %6.2f    Glt= %6.1f\n\n", NULT, GLT);
           
// reduced stiffness matrix (angle = 0)
      MATNUL(Q0,NMAX); 
      CAL_Q0(Q0);  
      MATPR2 ("  BASIC REDUCED STIFFNESS MATRIX IN GPa:", NMAX, Q0);
     
// calculate reduced stiffness matrices for different angles <> 0
      for (I=1; I<=NCOUCHES; I++) {
        if (TH[I] == 0.0)  MATEGAL(I,NMAX,Q,Q0);
        else  CAL_Q(I,Q0,Q,TH[I]);
        fprintf(fp,"\n  Layer #%d\n", I);
        MATPR1("  REDUCED STIFFNESS MATRIX IN GPa:", I, NMAX, Q);
      }  

// calculate A matrix
// A = sum of 1e3*EP(k)*Q(K)
      MATNUL(A, NMAX);
      for (I=1; I<=NMAX; I++)
	    for (J=1; J<=NMAX; J++) {
          for (K=1; K<=NCOUCHES; K++)
            A[I][J]=A[I][J]+1000*(EP[K]*Q[K][I][J]);
          if (fabs(A[I][J]) < 1e-8) A[I][J]=0.0;
        }  
      MATPR2("  A MATRIX IN 1E6 N/M:", NMAX, A);

// calculate B matrix
// B = sum of 1e6*(1/2)*(h(k)*h(k)-h(k-1)*h(k-1))*Q(k)
      MATNUL(B,NMAX);
      HM=0.0;
      for (I=1; I<=NCOUCHES; I++) HM += EP[I];
// average height of stack
      HM=HM/2.0;
// lower height layer #1
      H[1] = -HM;
      for (I=2; I<=NCOUCHES+1; I++)  H[I]=H[I-1]+EP[I-1];
      for (I=1; I<=NMAX; I++)
	    for (J=1; J<=NMAX; J++) {
          for (K=1; K<=NCOUCHES; K++)
            B[I][J]=B[I][J]+1e6*0.5*((H[K+1]*H[K+1]-H[K]*H[K])*Q[K][I][J]);
          if (fabs(B[I][J]) < 1e-8)  B[I][J]=0.0;
		}
      MATPR2("  B MATRIX IN 1E3 N:", NMAX, B);                           
    
// 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)
      MATNUL(D,NMAX);
      for (I=1; I<=NMAX; I++)
	    for (J=1; J<=NMAX; J++) {
          for (K=2; K<=NCOUCHES+1; K++)
            D[I][J]=D[I][J]+1e9*(1.0/3.0)*((H[K]*H[K]*H[K]-H[K-1]*H[K-1]*H[K-1])*Q[K-1][I][J]);
          if (fabs(D[I][J]) < 1e-8)  D[I][J]=0.0;
        }  
      MATPR2("  D MATRIX IN N.DM:", NMAX, D);
    
/* assemble the stiffness matrix of the laminate 
           ! A B !
       R = !     !
           ! B D !           */
      for (I=1; I<7; I++)
        for (J=1; J<7; J++)
          R[I][J]=0.0;
      for (I=1; I<=NMAX; I++)
	    for (J=1; J<=NMAX; J++)
          R[I][J] = 1e6*A[I][J];
      for (I=1; I<=NMAX; I++)
        for (J=4; J<=NMAX+3; J++)
          R[I][J] = 1e3*B[I][J-3];
      for (I=4; I<=NMAX+3; I++)
        for (J=1; J<=NMAX; J++)
          R[I][J] = 1e3*B[I-3][J];
      for (I=4; I<=NMAX+3; I++)
        for (J=4; J<=NMAX+3; J++)
          R[I][J] = D[I-3][J-3];
      MATPR6("  STIFFNESS R MATRIX:", 6, R);

// start index in R from zero before using MATINV
	  for (I=0; I<6; I++)
        for (J=0; J<6; J++)
			R[I][J] = R[I+1][J+1];

// Direct inversion of R matrix, stored in R
      MATINV(6,0,R,X,&DET);

      fprintf(fp,"\n  DETERMINANT OF R: %12.4e\n\n", DET);
      
// restart index in R from one to continue
	  for (I=6; I>0; I--)
        for (J=6; J>0; J--)
			R[I][J] = R[I-1][J-1];

	  for (I=1; I<7; I++)
        for (J=1; J<7; J++)
          R[I][J]=1e9*R[I][J];
      MATPR6("  INVERSE MATRIX OF R:", 6, R);

// Calculate déformations eps0xx, eps0yy, gam0xy,kx,ky,kxy
// in main axes, stored in matrix E0(6)
      MATMUL6(R,N,E0,6);
      fprintf(fp,"\n  Membrane deformations (in mm) and curbatures:\n\n");

      fprintf(fp,"    eps0xx: %12.5f\n", 1e3*E0[1]);
      fprintf(fp,"    eps0yy: %12.5f\n", 1e3*E0[2]);
      fprintf(fp,"    eps0xy: %12.5f\n", 1e3*E0[3]);
      fprintf(fp,"    kx    : %12.5f\n", E0[4]);
      fprintf(fp,"    ky    : %12.5f\n", E0[5]);
      fprintf(fp,"    kxy   : %12.5f\n", E0[6]);

/* Deformations in reference axes (x,y)
      ! epsxx !   ! eps0xx !    ! kx  !
      ! epsyy ! = ! eps0yy ! +z ! ky  !
      ! gamxy !   ! gam0xy !    ! kxy !    */
      fprintf(fp,"\n  Deformations in reference axes (x,y):\n");
      fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[1]);
      for (I=1; I<=NMAX; I++) EDEB[I]=E0[I]+H[1]*E0[I+3];
      fprintf(fp,"    epsxx: %12.5f mm\n", 1e3*EDEB[1]);
      fprintf(fp,"    epsyy: %12.5f mm\n", 1e3*EDEB[2]);
      fprintf(fp,"    epsxy: %12.5f mm\n", 1e3*EDEB[3]);
      fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[NCOUCHES+1]);
      for (I=1; I<=NMAX; I++) EFIN[I]=E0[I]+H[NCOUCHES+1]*E0[I+3];
      fprintf(fp,"    epsxx: %12.5f mm\n", 1e3*EFIN[1]);
      fprintf(fp,"    epsyy: %12.5f mm\n", 1e3*EFIN[2]);
      fprintf(fp,"    epsxy: %12.5f mm\n", 1e3*EFIN[3]);

// Deformations in each layer in main axes of layer
      fprintf(fp,"\n  Deformations in each layer in main axes:\n");
      for (K=1; K<=NCOUCHES; K++) {
        INITT(T,TH[K]);
        for (I=1; I<=NMAX; I++)  TEMP[I]=E0[I];
        MATMUL(T,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) AK[K][I]=TEMP1[I];
        for (I=1; I<=NMAX; I++) TEMP[I]=E0[I+3];
        MATMUL(T,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) BK[K][I]=TEMP1[I];
        fprintf(fp,"\n  Layer #%d\n", K);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K]);
        for (I=1; I<=NMAX; I++)  EDEB[I]=AK[K][I]+H[K]*BK[K][I];
        fprintf(fp,"    epsxx: %12.5f mm\n", 1e3*EDEB[1]);
        fprintf(fp,"    epsyy: %12.5f mm\n", 1e3*EDEB[2]);
        fprintf(fp,"    epsxy: %12.5f mm\n", 1e3*EDEB[3]);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K+1]);
        for (I=1; I<=NMAX; I++)  EFIN[I]=AK[K][I]+H[K+1]*BK[K][I];
        fprintf(fp,"    epsxx: %12.5f mm\n", 1e3*EFIN[1]);
        fprintf(fp,"    epsyy: %12.5f mm\n", 1e3*EFIN[2]);
        fprintf(fp,"    epsxy: %12.5f mm\n", 1e3*EFIN[3]);
      }

// Stresses in each layer in reference axes

      fprintf(fp,"\n  Stresses in each layer in reference axes:\n");
      for (K=1; K<=NCOUCHES; K++) {
        for (I=1; I<=NMAX; I++) TEMP[I]=E0[I];
		for (I=1; I<=NMAX; I++)
          for (J=1; J<=NMAX; J++)
            A[I][J]=Q[K][I][J];
        MATMUL(A,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) AK[K][I]=TEMP1[I];
        for (I=1; I<=NMAX; I++) TEMP[I]=E0[I+3];
        MATMUL(A,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) BK[K][I]=TEMP1[I];
        fprintf(fp,"\n  Layer #%d\n", K);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K]);
        for (I=1; I<=NMAX; I++) EDEB[I]=1e3*(AK[K][I]+H[K]*BK[K][I]);
        fprintf(fp,"    sigma L : %12.5f MPa\n", EDEB[1]);
        fprintf(fp,"    sigma T : %12.5f MPa\n", EDEB[2]);
        fprintf(fp,"    sigma LT: %12.5f MPa\n", EDEB[3]);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K+1]);
        for (I=1; I<=NMAX; I++) EFIN[I]=1e3*(AK[K][I]+H[K+1]*BK[K][I]);
        fprintf(fp,"    sigma L : %12.5f MPa\n", EFIN[1]);
        fprintf(fp,"    sigma T : %12.5f MPa\n", EFIN[2]);
        fprintf(fp,"    sigma LT: %12.5f MPa\n", EFIN[3]);
      }

// Stresses in each layer in main axes

      fprintf(fp,"\n  Stresses in each layer in main axes:\n");
      for (K=1; K<=NCOUCHES; K++) {
        INITT1(T,TH[K]);
        for (I=1; I<=NMAX; I++) TEMP[I]=AK[K][I];
        MATMUL(T,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) AK1[K][I]=TEMP1[I];
        for (I=1; I<=NMAX; I++) TEMP[I]=BK[K][I];
        MATMUL(T,TEMP,TEMP1,NMAX);
        for (I=1; I<=NMAX; I++) BK1[K][I]=TEMP1[I];
        fprintf(fp,"\n  Layer #%d\n", K);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K]);
        for (I=1; I<=NMAX; I++)  EDEB[I]=1e3*(AK1[K][I]+H[K]*BK1[K][I]);
        fprintf(fp,"    sigma L : %12.5f MPa\n", EDEB[1]);
        fprintf(fp,"    sigma T : %12.5f MPa\n", EDEB[2]);
        fprintf(fp,"    sigma LT: %12.5f MPa\n", EDEB[3]);
        fprintf(fp,"   For z = %12.5f mm:\n", 1e3*H[K+1]);
        for (I=1; I<=NMAX; I++)  EFIN[I]=1E3*(AK1[K][I]+H[K+1]*BK1[K][I]);
        fprintf(fp,"    sigma L : %12.5f MPa\n", EFIN[1]);
        fprintf(fp,"    sigma T : %12.5f MPa\n", EFIN[2]);
        fprintf(fp,"    sigma LT: %12.5f MPa\n", EFIN[3]);
      }
     
      fprintf(fp,"\n");
      fclose(fp);
      printf("\n  END OF PROGRAM.\n\n");

}

// End of file Compos04.cpp