/*************************************************************************************************************************************/
/*                                                                                                                                   */
/*        T E S T   D U   C A L C U L   D E S   M O M E N T S   D ' I N E R T I E  :                                                 */
/*                                                                                                                                   */
/*                                                                                                                                   */
/*        Author of '$xtc/inertie.01$c' :                                                                                            */
/*                                                                                                                                   */
/*                    Jean-Francois COLONNA (LACTAMME, AAAAMMJJhhmmss).                                                              */
/*                                                                                                                                   */
/*************************************************************************************************************************************/

extern    double    sqrt();

#define   dimX      3
#define   dimY      3

#define   MATRICE(x,y)                                                                                                                  \
                    matrice[y][x]

#define   Xc                                                                                                                            \
                    ((dimX-1)/2)
#define   Yc                                                                                                                            \
                    ((dimY-1)/2)

#define   X                                                                                                                             \
                    (x - Xc)
#define   Y                                                                                                                             \
                    (y - Yc)

#define   NIVEAU_1  1
#define   NIVEAU_2  2
#define   NIVEAU_3  3

#define   inertie(niveau_bas,niveau_central,niveau_haut)                                                                                \
                    {                                                                                                                   \
                    double    matrice[dimY][dimX];                                                                                      \
                    int       x,y;                                                                                                      \
                                                                                                                                        \
                    double    Ixx=0,Ixyx=0,Iyy=0;                                                                                       \
                                                                                                                                        \
                    double    a,b,c;                                                                                                    \
                    double    discriminant;                                                                                             \
                                                                                                                                        \
                    double    s1=0,s2=0;                                                                                                \
                                                                                                                                        \
                    double    v1X,v1Y;                                                                                                  \
                    double    v2X,v2Y;                                                                                                  \
                                                                                                                                        \
                    printf("\n\n");                                                                                                     \
                                                                                                                                        \
                    for       (x=0 ; x<dimX ; x++)                                                                                      \
                              {                                                                                                         \
                              for       (y=0 ; y<dimY ; y++)                                                                            \
                                        {                                                                                               \
                                        if        (y < Yc)                                                                              \
                                                  {                                                                                     \
                                                  MATRICE(x,y) = niveau_bas;                                                            \
                                        /* Initialisation de la partie basse de la matrice...                                        */ \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                                                                                                                        \
                                        if        (y == Yc)                                                                             \
                                                  {                                                                                     \
                                                  MATRICE(x,y) = niveau_central;                                                        \
                                        /* Initialisation de la ligne mediane de la matrice...                                       */ \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                                                                                                                        \
                                        if        (y > Yc)                                                                              \
                                                  {                                                                                     \
                                                  MATRICE(x,y) = niveau_haut;                                                           \
                                        /* Initialisation de la partie haute de la matrice...                                        */ \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  }                                                                                     \
                                        }                                                                                               \
                              }                                                                                                         \
                                                                                                                                        \
                    for       (y=0 ; y<dimY ; y++)                                                                                      \
                              {                                                                                                         \
                              printf("\n");                                                                                             \
                              for       (x=0 ; x<dimX ; x++)                                                                            \
                                        {                                                                                               \
                                        printf("   %.4f",MATRICE(x,y));                                                                 \
                                        }                                                                                               \
                              }                                                                                                         \
                                                                                                                                        \
                    for       (x=0 ; x<dimX ; x++)                                                                                      \
                              {                                                                                                         \
                              for       (y=0 ; y<dimY ; y++)                                                                            \
                                        {                                                                                               \
                                        Ixx  = Ixx + ((Y*Y)*MATRICE(x,y));                                                              \
                                        Ixyx = Ixyx + ((X*Y)*MATRICE(x,y));                                                             \
                                        Iyy  = Iyy + ((X*X)*MATRICE(x,y));                                                              \
                                        /* Calcul des moments d'inertie.                                                             */ \
                                        }                                                                                               \
                              }                                                                                                         \
                                                                                                                                        \
                    printf("\n Ixx=%.4f",Ixx);                                                                                          \
                    printf("\n Ixyx=%.4f",Ixyx);                                                                                        \
                    printf("\n Iyy=%.4f",Iyy);                                                                                          \
                                                                                                                                        \
                    a = 1;                                                                                                              \
                    b = -(Ixx+Iyy);                                                                                                     \
                    c = (Ixx*Iyy) - (Ixyx*Ixyx);                                                                                        \
                    discriminant = (b*b) - (4*a*c);                                                                                     \
                                                                                                                                        \
                    if        (discriminant >= 0)                                                                                       \
                              {                                                                                                         \
                              s1 = ((-b) + sqrt(discriminant)) / (2*a);                                                                 \
                              s2 = ((-b) - sqrt(discriminant)) / (2*a);                                                                 \
                                        /* Calcul des deux valeurs propres...                                                        */ \
                              printf("\n s1=%.4f",s1);                                                                                  \
                              printf("\n s2=%.4f",s2);                                                                                  \
                                                                                                                                        \
                              if        (discriminant >= 0)                                                                             \
                                        {                                                                                               \
                                        double    l1,l2;                                                                                \
                                                                                                                                        \
                                        double    v3X,v3Y;                                                                              \
                                        double    v4X,v4Y;                                                                              \
                                                                                                                                        \
                                        double    l3,l4;                                                                                \
                                                                                                                                        \
                                        if        (Ixyx != 0)                                                                           \
                                                  {                                                                                     \
                                                  v1X = 1;                                                                              \
                                                  v1Y = -((Ixx-s1) / Ixyx);                                                             \
                                                  l1 = sqrt((v1X*v1X) + (v1Y*v1Y));                                                     \
                                                  v1X = v1X / l1;                                                                       \
                                                  v1Y = v1Y / l1;                                                                       \
                                                                                                                                        \
                                                  v2X = 1;                                                                              \
                                                  v2Y = -((Ixx-s2) / Ixyx);                                                             \
                                                  l2 = sqrt((v2X*v2X) + (v2Y*v2Y));                                                     \
                                                  v2X = v2X / l2;                                                                       \
                                                  v2Y = v2Y / l2;                                                                       \
                                        /* Calcul des deux vecteurs propres...                                                       */ \
                                                                                                                                        \
                                                  printf("\n v1=(%.4f,%.4f)",v1X,v1Y);                                                  \
                                                  printf("\n v2=(%.4f,%.4f)",v2X,v2Y);                                                  \
                                                                                                                                        \
                                                  v3X = v1X + v2X;                                                                      \
                                                  v3Y = v1Y + v2Y;                                                                      \
                                                  l3 = sqrt((v3X*v3X) + (v3Y*v3Y));                                                     \
                                                  v3X = v3X / l3;                                                                       \
                                                  v3Y = v3Y / l3;                                                                       \
                                                                                                                                        \
                                                  v4X = v1X - v2X;                                                                      \
                                                  v4Y = v1Y - v2Y;                                                                      \
                                                  l4 = sqrt((v4X*v4X) + (v4Y*v4Y));                                                     \
                                                  v4X = v4X / l4;                                                                       \
                                                  v4Y = v4Y / l4;                                                                       \
                                        /* Calcul de deux vecteurs a 45 degres des vecteurs propres...                               */ \
                                                                                                                                        \
                                                  printf("\n v3=(%.4f,%.4f)",v3X,v3Y);                                                  \
                                                  printf("\n v4=(%.4f,%.4f)",v4X,v4Y);                                                  \
                                                  }                                                                                     \
                                        else                                                                                            \
                                                  {                                                                                     \
                                                  printf("\n Ixyx est nul");                                                            \
                                                  }                                                                                     \
                                        }                                                                                               \
                              else                                                                                                      \
                                        {                                                                                               \
                                        printf("\n le discriminant est nul");                                                           \
                                        }                                                                                               \
                              }                                                                                                         \
                    else                                                                                                                \
                              {                                                                                                         \
                              printf("\n le discriminant est negatif");                                                                 \
                              }                                                                                                         \
                                                                                                                                        \
                    }                                                                                                                   \
                                        /* Procedure de calcul des moments d'inertie d'une certaine matrice...                       */

int       main()
     {

     inertie(NIVEAU_2,NIVEAU_2,NIVEAU_2);
     inertie(NIVEAU_1,NIVEAU_2,NIVEAU_1);
     inertie(NIVEAU_1,NIVEAU_2,NIVEAU_3);
     }