/*
 *  =================================================================== 
 *  Compute buckling loads and mode shapes for slender elastic rod.
 *  =================================================================== 
 *  Copyright (C) 1995-96 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: Mark Austin                                   June 1995
 *  =================================================================== 
 */

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

/* Define section and material properties of elastic rod */

static float fE     =  200;   /* Young's modulus of elasticity   */
static float fI     = 0.10;   /* Cross section moment of inertia */
static float fLength = 5.0;   /* Length of the rod               */

int main( void ) {
enum { NoDof = 5, One = 1 , NoEigen = 2 };
MATRIX *spA, *spB, *spEigenvalue, *spEigenvector;
int   iRow, iCol, iMin, iMax;
float fBuckleLoad;

   /* [a] : Define profile of "A" and "B" matrices */
 
   int *ipProfile = (int *) safeCalloc( NoDof, sizeof(int), __FILE__, __LINE__ );
   ipProfile[ 0 ] = 1;
   for(iRow = 1; iRow < NoDof; iRow = iRow + 1 )
       ipProfile[ iRow ] = 2;

   /* [b] : Allocate memory for skyline matrices */

   spA = matAllocSkyline( "A", DoubleArray, NoDof , NoDof , ipProfile );
   spB = matAllocSkyline( "B", DoubleArray, NoDof , NoDof , ipProfile );

   /* [c] : Fill-in and print contents of "A" and "B" skyline matrices */

   for(iCol =    1; iCol <= NoDof; iCol = iCol + 1 ) {
   for(iRow = iCol; iRow >=     1; iRow = iRow - 1 ) {

       iMin = MIN( iRow, iCol);
       iMax = MAX( iRow, iCol);

       if(iCol == iRow) {
          spA->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] =  2.0;
          spB->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] = 10.0;
       } 
       else if(abs(iCol - iRow) == 1) {
          spA->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] = -1.0;
          spB->uMatrix.dpp[ iMax-1 ][ iMax-iMin+1 ] =  1.0;
       } 
   }
   }

   matPrint( spA );
   matPrint( spB );

   /* [d] : Allocate memory for eigenvalues and eigenvectors */

   spEigenvalue  = matAllocIndirect( "Eigenvalues", DoubleArray, NoEigen,     One );
   spEigenvector = matAllocIndirect("Eigenvectors", DoubleArray,   NoDof, NoEigen );

   /* [e] : Use subspace iteration to compute eigenvalues and eigenvectors */

   matSubspace( spA, spB, spEigenvalue, spEigenvector , NoEigen );

   /* [f] : Compute and print buckling loads and mode shapes */

   printf("\n*** Buckling Loads and Mode Shapes \n");

   for (iRow = 1; iRow <= NoEigen; iRow = iRow + 1 ) {
        fBuckleLoad = spEigenvalue->uMatrix.dpp[ NoEigen-iRow ][ 0 ];
        printf("\nMode %3d : Eigenvalue    = %8.5f\n", iRow, fBuckleLoad );

        fBuckleLoad *= 12*fE*fI/(fLength/6.0)/(fLength/6.0);
        printf("         : Buckling Load = %8.5f\n", fBuckleLoad );

        printf("Mode Shape : %8.5f\n", spEigenvector->uMatrix.dpp[ 0 ][ NoEigen-iRow ]);
        for (iCol = 2; iCol <= NoDof; iCol = iCol + 1 ) {
             printf("           : %8.5f\n", spEigenvector->uMatrix.dpp [iCol-1][ NoEigen-iRow ]);
        }
   }

}