/*
 *  =========================================================================
 *  Compute Engineering Properties in Channel Cross Section.
 *                                                                     
 *  Note : The cross section shape is assembled from three simple rectangular
 *         shapes.                                                                     
 *  =========================================================================
 *  Copyright (C) 1996 by Mark Austin and David Mazzoni.
 *                                                                     
 *  This software is provided "as is" without express or implied warranty.
 *  Permission is granted to use this software on any computer system,
 *  and to redistribute it freely, subject to the following restrictions:
 * 
 *  1. The authors are not responsible for the consequences of use of
 *     this software, even if they arise from defects in the software.
 *  2. The origin of this software must not be misrepresented, either
 *     by explicit claim or by omission.
 *  3. Altered versions must be plainly marked as such, and must not
 *     be misrepresented as being the original software.
 *  4. This notice is to remain intact.
 *                                                                     
 *  Written by: M.A. Austin                                      August, 1996
 *  =========================================================================
 */

#include <stdio.h>
#include <math.h>
#include "matrix.h"
#include "vector.h"
#include "miscellaneous.h"
#include "polygon.h"

/* 
 *   Define [x,y] coordinates of four corners in polygon rectangle.
 */

static struct {
       float fXcoord1, fYcoord1;    /* [x,y] coordinates of node 1 */
       float fXcoord2, fYcoord2;    /* [x,y] coordinates of node 2 */
       float fXcoord3, fYcoord3;    /* [x,y] coordinates of node 3 */
       float fXcoord4, fYcoord4;    /* [x,y] coordinates of node 4 */
       } saPolygon[] = {
              3.0,   3.0,  13.0,  3.0,  13.0,   2.0,   3.0,   2.0,
              3.0,   6.0,   4.0,  6.0,   4.0,   3.0,   3.0,   3.0,
              3.0,   7.0,   7.0,  7.0,   7.0,   6.0,   3.0,   6.0,
       };

#define NoSubPolygons (sizeof(saPolygon))/sizeof(saPolygon[0])

int main( void ) {
MATRIX  *spInertia, *spEigenvalue, *spEigenvector;
POLYGON *spPolygon1 = (POLYGON *) NULL; 
POLYGON *spPolygon2; 
EDGE *spEdge;
float Xcoord, Ycoord;
float fArea, fIxx, fIyy, fIxy;
float  fXcentroid, fYcentroid;
int ii, ij;

   /* [a] : Build polygon assembly by walking along saPolygonCoords and incrementally */
   /*       adding sub-polygons to polygon data structure                             */

      for( ii = 1 ; ii <= NoSubPolygons; ii = ii + 1) {
      if( ii == 1 ) {
          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord1,
                                    saPolygon[ii-1].fYcoord1, 0.0 );
          spPolygon1 = polyCreate("Polygon 1", spEdge );

          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord2,
                                    saPolygon[ii-1].fYcoord2, 0.0 );
          polyAddEdge( spPolygon1, spEdge );
          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord3,
                                    saPolygon[ii-1].fYcoord3, 0.0 );
          polyAddEdge( spPolygon1, spEdge );
          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord4,
                                    saPolygon[ii-1].fYcoord4, 0.0 );
          polyAddEdge( spPolygon1, spEdge );
      } else {

          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord1,
                                    saPolygon[ii-1].fYcoord1, 0.0 );
          spPolygon2 = polyCreate("Polygon 1", spEdge );

          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord2,
                                    saPolygon[ii-1].fYcoord2, 0.0 );
          polyAddEdge( spPolygon2, spEdge );
          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord3,
                                    saPolygon[ii-1].fYcoord3, 0.0 );
          polyAddEdge( spPolygon2, spEdge );
          spEdge  = polyCreateEdge( saPolygon[ii-1].fXcoord4,
                                    saPolygon[ii-1].fYcoord4, 0.0 );
          polyAddEdge( spPolygon2, spEdge );

          polyAppend( spPolygon1 , spPolygon2 );
      }
      }

      polyPrint( spPolygon1 );

   /* [b] : Compute and print engineering properties */

      fArea = polyProperty( polyPropertyArea, spPolygon1 );
      fIxx  = polyProperty( polyPropertyIxx,  spPolygon1 );
      fIyy  = polyProperty( polyPropertyIyy,  spPolygon1 );
      fIxy  = polyProperty( polyPropertyIxy,  spPolygon1 );

      fXcentroid = polyProperty( polyPropertyXcentroid, spPolygon1 )/fArea;
      fYcentroid = polyProperty( polyPropertyYcentroid, spPolygon1 )/fArea;

      printf("\n\n");
      printf("*** Engineering Properties\n");
      printf("*** ======================\n\n");
   
      printf("*** Area                  = %6.2f\n", fArea );
      printf("*** Centroid [x,y]        = [ %6.2f, %6.2f ] \n\n", fXcentroid, fYcentroid );

      printf("*** I_xx (about origin)   = %6.2f\n", fIxx );
      printf("*** I_yy (about origin)   = %6.2f\n", fIyy );
      printf("*** I_xy (about origin)   = %6.2f\n", fIxy );

   /* [c] : Compute engineering properties about polygon centroid */

      printf("\n");
      printf("*** I_xx (about centroid) = %6.2f\n", (fIxx - fYcentroid*fYcentroid*fArea));
      printf("*** I_yy (about centroid) = %6.2f\n", (fIyy - fXcentroid*fXcentroid*fArea));
      printf("*** I_xy (about centroid) = %6.2f\n", (fIxy - fXcentroid*fYcentroid*fArea));

   /* [d] : Setup inertia matrix for engineering properties */

      spInertia = matAllocIndirect("Inertia",  DoubleArray, 2, 2 );

      spInertia->uMatrix.dpp[ 0 ][ 0 ] =   fIxx - fYcentroid*fYcentroid*fArea;
      spInertia->uMatrix.dpp[ 1 ][ 1 ] =   fIyy - fXcentroid*fXcentroid*fArea;
      spInertia->uMatrix.dpp[ 0 ][ 1 ] = -(fIxy - fYcentroid*fXcentroid*fArea);
      spInertia->uMatrix.dpp[ 1 ][ 0 ] = -(fIxy - fYcentroid*fXcentroid*fArea);

      matPrint( spInertia );

   /* [e] : Use householder algorithm to compute eigenvalues/eigenvectors */

      spEigenvalue  = matAllocIndirect( "Eigenvalues", DoubleArray, 2, 1);
      spEigenvector = matAllocIndirect("Eigenvectors", DoubleArray, 2, 2);

      matHouseholder( spInertia, spEigenvalue, spEigenvector );

      printf("\n*** Principal Component Analysis\n");

      matPrint( spEigenvalue  );
      matPrint( spEigenvector );
}