#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>

#define MAX 300

void enter_data(void);
void files(void);

void bc(void);
void symmetry(void);
void global1(void);
void print_global_before_BC(void);
void local(void);
void axial_stiffness(void);

void inverse_iteration(void);
void iteration(void);

void interactive(void);
void inverse_iteration(void);
void kshift(void);

// ********************************************************************************


void matrix_times_vector(long nd, double aaa[], double w[][MAX],double bbb[] );
void vector_times_vector(long nd, double *csum, double aaa[], double bbb[]);

void copy_matrix(long nd, double aaa[][MAX], double bbb[][MAX]);
void zero_matrix(long nd, double aaa[][MAX]);

void zero_vector(long nd, double aaa[]);

void print_matrix(long nd, double aaa[][MAX]);
void print_vector(long nd, double aaa[]);

// ********************************************************************************


const double TPI = 2.*atan2(0.,-1);
const double  PI =    atan2(0.,-1);

const long numberElements = 24;

long i,ii,ij,ik,j,jj,jk,k,kk;
long m,n,ne,nn;

int knv=0;

int apply;
int iflag;
int nx;
int ipm;

int iload;

int inew;

double elementLength;

double aa,ah,aj,aq;
double ajl;
double ccc;
double det;
double eigen;
double kd;

double mc[30];

double length,inertia[30];
double modulus,mpv,area,mpl,h,load;

double length_mass,cm;

double scale;
double shift=0.;

double x,xx;

double kl[30][5][5],ml[30][5][5];
double kg[MAX][MAX],mg[MAX][MAX],gmg[MAX][MAX];
double a[MAX][MAX];
double xi[MAX];
double bb[MAX];
double ac[MAX][MAX];
double ww[MAX];
double b[MAX];

double ODfixed,IDfixed;
double ODfree,IDfree;

double OD,ID;


double model_mass=0.;


double R,r,volume,volume1,volume2;
double bm;
long mif;



FILE *pFile[20];
char filename[10][25];


void main()
{

	ne = numberElements;

// assume evenly spaced


    printf("\n\n cantilever_mass_axial.cpp \n");
	printf("\n ver 1.2    Jan 26, 2005 ");
	printf("\n by Tom Irvine   Email: tomirvine@aol.com \n");

	printf("\n This program calculates the natural frequency of a");
	printf("\n cantilever beam with an added point mass at the free");
	printf("\n end of the beam and with an applied axial tensile force. \n\n");

 //	printf("\n files \n");

	files();

//	printf("\n enter data \n");

    enter_data();

//	printf("\n local \n");


	while(1)
	{

		local();

//	printf("\n global1 \n");

		global1();

		if(iload==1)
		{

			axial_stiffness();
		}

//	printf("\n symmetry \n");

		symmetry();


//	printf("\n print global before BC \n");

		print_global_before_BC();

//	printf("\n bc \n");

		bc();

//	printf("\n inverse iteration \n");

		inverse_iteration();

		printf("\n\n Run model again with new axial load? \n");
		printf("   1=yes  2= no \n");

		scanf("%d",&inew);

		if(inew != 1 )
		{
			break;
		}
		else
		{
			printf("\n Enter axial tensile load (lbf) \n");
			scanf("%lf",&load); 

			fprintf(pFile[5],"\n\n New load Case, Axial load = %12.4g lbf ",load);
		}
	}


	printf("\n\n The analysis data is written to file: %s \n",filename[5]);


    printf("\n Press any key to exit.\n");
    getch();


}
void print_global_before_BC(void)
{
	for(i=1; i<=nn; i++)
	{
		for(j=1; j<=nn; j++)
		{
			fprintf(pFile[2]," %11.4g",mg[i][j]);
		}
		fprintf(pFile[2],"\n");
	}

	for(i=1; i<=nn; i++)
	{
		for(j=1; j<=nn; j++)
		{
			fprintf(pFile[2]," %11.4g",kg[i][j]);
		}
		fprintf(pFile[2],"\n");
	}
}
void symmetry(void)
{

	//  symmetric

	for(i=1; i<=nn; i++)
	{
		for(j=i+1; j<=nn; j++)
		{
			kg[j][i]=kg[i][j];
			mg[j][i]=mg[i][j];
		}
	}

}
void bc(void)
{

	// cantilever bc


	for(i=3; i<=nn; i++)
	{
		for(j=3; j<=nn; j++)
		{
			fprintf(pFile[3]," %11.4g",mg[i][j]);
		}
		fprintf(pFile[3],"\n");
	}
	for(i=3; i<=nn; i++)
	{
		for(j=3; j<=nn; j++)
		{
			fprintf(pFile[3]," %11.4g",kg[i][j]);
		}
		fprintf(pFile[3],"\n");
	}

	fclose(pFile[3]);


	n=0;

	for(i=3; i<=nn; i++)
	{
		for(j=3; j<=nn; j++)
		{
			mg[i-3][j-3]=mg[i][j];
			kg[i-3][j-3]=kg[i][j];
		}
		n++;
	}


}

void files(void)
{
   
  
   strcpy(filename[1],"km.out");
   pFile[1]=fopen(filename[1], "w");

   strcpy(filename[2],"kmg.out");
   pFile[2]=fopen(filename[2], "w");

   strcpy(filename[3],"kmbc.out");
   pFile[3]=fopen(filename[3], "w");

   strcpy(filename[4],"eigen_data.out");
   pFile[4]=fopen(filename[4], "w");

   strcpy(filename[5],"cantilever.out");
   pFile[5]=fopen(filename[5], "w");


   printf("\n");

   for(i=1; i<=5; i++)
   {
      if(pFile[i] == NULL )
      {
         printf(" Failed to open file: %s \n", filename[i]);
         fclose(pFile[i]);
      }
      else
      {
//         printf(" File: %s opened. \n", filename[i]);
      }
   }

    fprintf(pFile[5],"\n\n cantilever_mass_axial.cpp \n");
	fprintf(pFile[5],"\n ver 1.1  by Tom Irvine   Email: tomirvine@aol.com \n");


}

void global1(void)
{

	for(i=0; i<100; i++)
	{
		for(j=0; j<100; j++)
		{
			kg[i][j]=0.;
			mg[i][j]=0.;
		}
	}

//	printf("ref3b  k11=%lf \n",kl[1][1]);

	for(m=1; m<=numberElements; m++  )
	{

		i=2*m-1;

		kg[i][i]  +=kl[m][1][1];
		kg[i][i+1]+=kl[m][1][2];
		kg[i][i+2]+=kl[m][1][3];
		kg[i][i+3]+=kl[m][1][4];
		
        j=i+1;

		kg[j][j]  +=kl[m][2][2];
		kg[j][j+1]+=kl[m][2][3];
		kg[j][j+2]+=kl[m][2][4];

        j=i+2;

		kg[j][j]  +=kl[m][3][3];
		kg[j][j+1]+=kl[m][3][4];

        j=i+3;

		kg[j][j]  +=kl[m][4][4];


		mg[i][i]  +=ml[m][1][1];
		mg[i][i+1]+=ml[m][1][2];
		mg[i][i+2]+=ml[m][1][3];
		mg[i][i+3]+=ml[m][1][4];
		
        j=i+1;

		mg[j][j]  +=ml[m][2][2];
		mg[j][j+1]+=ml[m][2][3];
		mg[j][j+2]+=ml[m][2][4];

        j=i+2;

		mg[j][j]  +=ml[m][3][3];
		mg[j][j+1]+=ml[m][3][4];

        j=i+3;

		mg[j][j]  +=ml[m][4][4];

    }
	if( ipm==1)
	{
		nx=nn-1;

		mg[nx][nx]+=cm;
	}

			
}
void enter_data(void)
{

	printf("\n Note:  all length dimensions are in inches. \n");

	printf("\n Enter length \n");
	scanf("%lf",&length);



	printf("\n\n Enter O.D. at fixed end. \n");
    scanf("%lf",&ODfixed);



	printf("\n Enter I.D. at fixed end. \n");
    scanf("%lf",&IDfixed);
 

	if(ODfixed <= IDfixed)
	{
		printf("\n Data input error \n");
		printf("\n The OD must be greater than the ID \n");
		exit(1);
	}


	printf("\n\n Enter O.D. at free end. \n");
    scanf("%lf",&ODfree);
   

	printf("\n Enter I.D. at free end. \n");
    scanf("%lf",&IDfree);

  
	if(ODfree <= IDfree)
	{
		printf("\n Data input error \n");
		printf("\n The OD must be greater than the ID \n");
		exit(1);
	}

 
    R = (ODfixed/2.);
	r = (ODfree/2.);
    volume1 = PI*( pow(R,2)+R*r+pow(r,2.))*(length/3.);
    R = (IDfixed/2.);
	r = (IDfree/2.);
    volume2 = PI*( pow(R,2)+R*r+pow(r,2.))*(length/3.);

	volume=volume1-volume2;

	printf("\n\n Total material volume = %12.4g in^3\n",volume);
    fprintf(pFile[5],"\n\n Total material volume = %12.4g in^3\n",volume);


//	printf("\n Enter area moment of inertia (inch^4)\n");
//	scanf("%lf",&inertia);

//	printf("\n Enter cross-section area (in^2) \n");
//	scanf("%lf",&area);


	int imat;

	printf("\n Enter material ");
	printf("\n 1=aluminum ");
	printf("\n 2=steel ");
	printf("\n 3=other \n");

	scanf("%d",&imat);


	if( imat==1)
	{
        modulus = 1.0e+07;
		mpv = 0.1;
    }

	if( imat==2)
	{
        modulus = 3.0e+07;
		mpv = 0.28;
    }

	if( imat != 1 && imat !=2)
	{
		printf("\n Enter elastic modulus (lbf/in^2)\n");
		scanf("%lf",&modulus);

		mif=0;

		while(1)
		{
			printf("\n Enter mass input format \n 1=total mass 2=density \n");
			scanf("%ld",&mif);

			if(mif==1 || mif==2){break;}
		}

		if(mif==2)
		{
			printf("\n\n Enter material mass per volume (lbm/in^3) \n");
			scanf("%lf",&mpv);
		}
		if(mif==1)
		{
			printf("\n\n Enter beam mass not including end mass (lbm) \n");
			scanf("%lf",&bm);
			mpv=bm/volume;
		}
    }

    fprintf(pFile[5],"\n\n length = %8.4g inch\n",length);
	fprintf(pFile[5],"\n O.D. at fixed end = %8.4g inch\n",ODfixed);
    fprintf(pFile[5],"\n I.D. at fixed end = %8.4g inch\n",IDfixed);
    fprintf(pFile[5],"\n O.D. at free end = %8.4g inch\n",ODfree);
    fprintf(pFile[5],"\n I.D. at free end = %8.4g inch\n",IDfree);

    fprintf(pFile[5],"\n\n Elastic Modulus  = %8.4g (lbf/in^2)\n",modulus);
	  fprintf(pFile[5],"\n Material Density = %8.4g (lbm/in^3)\n",mpv);


//    printf("\n Beam mass = %8.3g lbm \n",length*area*mpv);


	mpv/=386.;

        printf("\n\n Add axial load at free end? ");
        printf("\n 1=yes  2=no \n");
  

        scanf("%d",&iload);

		load = 0.;

		if(iload==1)
		{
			printf("\n Enter axial tensile load (lbf) \n");
			scanf("%lf",&load); 
        }

     
   
 //       printf("\n Enter number of elements. (suggest 8 or more) \n");
//		scanf("%ld",&ne);

	

        printf("\n\n Add point mass at free end? ");
        printf("\n 1=yes  2=no \n");
  

        scanf("%d",&ipm);

        cm=0.;

        if(ipm==1)
        {       

//			printf("\n Enter distance from fixed end to concentrated mass (inch)\n");
//			scanf("%lf",&length_mass);

			length_mass=length;

			printf("\n Enter concentrated mass (lbm)\n");
			scanf("%lf",&cm);

         	cm/=386.;

        }

	
//	mpl=mpv*area;

//	printf("\n mpl= %12.4g\n",mpl);

    h=length/double(ne);

//  printf("\n h= %12.4g\n",h);  

//	mc=h*mpl/420.;

//  printf("\n mc= %12.4g\n",mc); 

	nn=2 + (ne*2);

	n=nn;

	double xx;

	elementLength = length/numberElements;

    double ODslope = (ODfixed-ODfree)/length;
    double IDslope = (IDfixed-IDfree)/length;


	xx = elementLength/2.;

	for(i=1;i<=numberElements;i++)
    {

		OD = ODfixed - xx *ODslope;
		ID = IDfixed - xx *IDslope;

		area= (PI/4.)*(pow(OD,2.)-pow(ID,2.));

        inertia[i]=(PI/64.)*(pow(OD,4.)-pow(ID,4.));

		mpl=mpv*area;
		mc[i]=h*mpl/420.;


//		printf(" %ld %12.4g %12.4g %12.4g\n",i,xx,OD,ID);


		xx+=elementLength;

		model_mass+=area*elementLength*mpv*386.;
	}

}
void local(void)
{
	for(i=1;i<=numberElements;i++)
	{

		ml[i][1][1]=156.*mc[i];
		ml[i][1][2]=22.*mc[i];
		ml[i][1][3]=54.*mc[i];
		ml[i][1][4]=-13.*mc[i];

		ml[i][2][2]=4.*mc[i];
		ml[i][2][3]=13.*mc[i];
		ml[i][2][4]=-3.*mc[i];

		ml[i][3][3]=156.*mc[i];
		ml[i][3][4]=-22.*mc[i];

		ml[i][4][4]=4.*mc[i];


		double kc=modulus*inertia[i]/pow(h,3.);


		kl[i][1][1]=12.*kc;

//	printf("ref1  k11=%lf \n",kl[1][1]);


		kl[i][1][2]=6.*kc;
		kl[i][1][3]=-12.*kc;
		kl[i][1][4]=6.*kc;

		kl[i][2][2]=4.*kc;
		kl[i][2][3]=-6.*kc;
		kl[i][2][4]=2.*kc;

		kl[i][3][3]=12.*kc;
		kl[i][3][4]=-6.*kc;

		kl[i][4][4]=4.*kc;

	}
/*
	for(i=1; i<=4; i++)
	{
		for(j=1; j<=4; j++)
		{
			fprintf(pFile[1]," %12.4g ",ml[i][j]);
		}
		fprintf(pFile[1],"\n");
	}
	for(i=1; i<=4; i++)
	{
		for(j=1; j<=4; j++)
		{
			fprintf(pFile[1]," %12.4g ",kl[i][j]);
		}
		fprintf(pFile[1],"\n");
	}
*/
	fclose(pFile[1]);
}

//************************************************************************

void inverse_iteration(void)
{

/*
	printf("\n inverse_iteration.cpp, ver 1.1, February 19, 2004");
	printf("\n\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com ");


    printf("\n\n This program uses inverse iteration to determine the lowest ");
	printf("\n eigenvalue and vector of the generalized eigen problem: \n\n");
    printf("\n     K - lambda M = 0");
	printf("\n\n where K and M are symmetric matrices.");

*/
	interactive();

    if(apply==2)
	{
		kshift();
	}

	iteration();


}

void iteration(void)
{

	double ddd=0.;

// printf("\n initialization" );

       
//       printf("\n n= %ld \n",n);

       for(i=0;i<MAX;++i)
       {
		  xi[i]=pow(double(i),1.2);
       }

	   zero_vector(MAX, bb);
	   zero_vector(MAX, ww);

	   zero_matrix(MAX, ac);

/****  copy data into reserve matrices  **/


//*****************************************

	long ijkl;

	for(ijkl=0; ijkl<40; ijkl++)
	{

	   copy_matrix(n, a, kg);
	   copy_matrix(n, ac, kg);
	   copy_matrix(n, gmg, mg);
     
	   matrix_times_vector(n, b, gmg, xi );
 

/****** begin forward elimination  *********/

//       printf("\n begin forward elimination \n" );

       for(i=0;i<n;++i)
       {
//	       printf(" %d ",i );
	       aq=fabs(a[i][i]);

	       for(jj=i+1;jj<n;++jj)
	       {
//	          printf("\n switching");
	          ah=fabs(a[jj][i]);

	          if(aq<ah)
	          {
		         aq=fabs(a[jj][i]);

		         for(kk=0;kk<n;++kk)
		         {
		             ww[kk]=a[jj][kk];
		              a[jj][kk]=a[i][kk];
		              a[i][kk]=ww[kk];
		         }
		         aj=b[jj];
		         b[jj]=b[i];
		         b[i]=aj;
	          }
	        }


	       for(j=i+1;j<n;++j)
	       {
	           if(  fabs( a[i][i] ) < 1.0e-26  )
	           {
		           printf("\n diagonal term  %d %d equals zero ",i,i);
		           
				   printf("\n  %12.4g\n",a[i][i] );
				   exit(1);
	            }
	        }

	        ik=i+1;
	        for(ii=ik;ii<n;++ii)
	        {
	            aa=a[ii][i]/a[i][i];
	            b[ii]+=-b[i]*a[ii][i]/a[i][i];

	            for(j=i;j<n;++j)
	            {
		           a[ii][j]+=-aa*a[i][j];
		           if( fabs(a[ii][j]) <1.0e-10)
		           {
		               a[ii][j]=0.;
		           }
	            }

			}

       }

/********* back substitution ***********************************/

//       printf("\n begin back substitution ");

       b[n-1]/=a[n-1][n-1];

       for(i=n-2;i>=0;i=i-1)
       {
	        for(k=i+1;k<n;++k)
	        {
//	             printf("\n %d %d",i,k);
	             b[i]+=(-b[k]*a[i][k]);
	        }
	        b[i]/=a[i][i];
       }

/***************************************************************/


// xt M x


	   matrix_times_vector(n, bb, gmg, b );


	   vector_times_vector(n, &ccc, bb, b);

	   ccc =sqrt(ccc);


	   for(i=0;i<n;++i)
       {
		   b[i]/=ccc; 
		   xi[i]=b[i];
	   }

	   matrix_times_vector(n, bb, ac, b );

	   vector_times_vector(n, &ccc, bb, b);

	   matrix_times_vector(n, bb, gmg, b );


	   ddd=0.;

	   vector_times_vector(n, &ddd, bb, b );

//	   printf(" ccc = %lf   ddd = %lf \n",ccc,ddd);

	   ccc/=ddd;

 

   }

   if(knv==0)
   {

	fprintf(pFile[5],"\n\n Beam  mass             = %12.4g lbm \n",model_mass);
	fprintf(pFile[5]," Point mass at Free End = %12.4g lbm \n",cm*386.);
	fprintf(pFile[5]," Total mass             = %12.4g lbm \n",model_mass+cm*386.);


	printf("\n\n Beam  mass             = %12.4g lbm \n",model_mass);
	printf(" Point mass at Free End = %12.4g lbm \n",cm*386.);
	printf(" Total mass             = %12.4g lbm \n",model_mass+cm*386.);

   }


// eigenvector 

 


	fprintf(pFile[4],"\n\n Eigenvalue \n\n");

	   ccc+=shift;

//	   printf("   lambda = %12.5g \n",ccc);
//	   printf("    omega = %12.5g \n",sqrt(ccc));
	  

	fprintf(pFile[4],"     lambda = %12.4g \n",ccc);
	fprintf(pFile[4],"      omega = %12.4g \n",sqrt(ccc));
	fprintf(pFile[4],"         fn = %12.4g Hz \n",sqrt(ccc)/TPI);


	if(knv==0)
	{
		fprintf(pFile[5],"\n\n Axial Tensile Load at Free End = %8.4g lbf ",load);
    }

   knv++;


             printf("\n\n Fundamental frequency fn = %12.3g Hz \n\n",sqrt(ccc)/TPI);
   fprintf(pFile[5],"\n\n Fundamental frequency fn = %12.3g Hz \n\n",sqrt(ccc)/TPI);

	
//   printf("\n\n Eigenvector \n\n");
fprintf(pFile[4],"\n\n Eigenvector \n\n");

       h=length/double(ne);

	   fprintf(pFile[4]," 0. \t  0.  \t 0. \n");

	   k=1;
       for(ij=0;ij<n;ij+=2)
       {
//	         printf("  %d %12.5e \n",ij,b[ij]);
	        fprintf(pFile[4]," %12.5e \t %12.5e \t %12.5e \n",k*h,b[ij],b[ij+1]);
       
			 k++;
	   
	   }
	  

}
void interactive(void)
{

  apply = 1;
  shift = 0.;

}
void kshift(void)
{
		for(i=0; i< n; i++)
		{

			for(j=0; j<n; j++)
			{

               kg[i][j]-= shift*mg[i][j];

			}

		}
}

// *************************************************************************************

void matrix_times_vector(long nd, double aaa[], double bbb[][MAX],double ccc[] )
{
	   long i,j;

	   for(i=0;i<nd;++i)
       {
		    aaa[i]=0.;

	        for(j=0;j<n;++j)
	        {
	             aaa[i]+=bbb[i][j]*ccc[j];
	        }
       }
}
void vector_times_vector(long nd, double *csum, double aaa[], double bbb[])
{
	   long i;

	   *csum=0.;

	   for(i=0;i<nd;++i)
       {
            *csum+=aaa[i]*bbb[i];
	   }


}
void copy_matrix(long nd, double aaa[][MAX], double bbb[][MAX])
{

	   long i,j;

       for(i=0;i<nd;++i)
       {
	      for(j=0;j<nd;++j)
	      {
	           aaa[i][j]=bbb[i][j];
	      }
       }
}
void zero_matrix(long nd, double aaa[][MAX])
{
	   long i,j;

       for(i=0;i<nd;++i)
       {
			for(j=0;j<nd;++j)
			{

	           aaa[i][j]=0.;
			}
       }
}
void zero_vector(long nd, double aaa[])
{
	   long i;

       for(i=0;i<nd;++i)
       {
			aaa[i]=0.;
       }
}
void print_matrix(long nd, double aaa[][MAX])
{
	   long i,j;

	   printf("\n");

	   for(i=0;i<nd;++i)
       {
			for(j=0;j<nd;++j)
			{
	            printf(" %8.3g",aaa[i][j]);
			}
			printf("\n");
       }	  
}
void print_vector(long nd, double aaa[])
{
	   long i;

	   printf("\n");

	   for(i=0;i<nd;++i)
       {
	      printf(" %8.3g \n",aaa[i]);

       }	  
}
void axial_stiffness(void)
{
    double ka = load/(30.*elementLength);

//	printf("\n ka = %12.4g \n");


	if(fabs(ka) > 1.0e-6)
	{

	for(m=1; m<=numberElements; m++  )
	{

		i=2*m-1;

		kg[i][i]  +=36*ka;
		kg[i][i+1]+=3*ka;
		kg[i][i+2]+=-36*ka;
		kg[i][i+3]+=3*ka;
		
        j=i+1;

		kg[j][j]  +=4*ka;
		kg[j][j+1]+=-3*ka;
		kg[j][j+2]+=-1*ka;

        j=i+2;

		kg[j][j]  +=36*ka;
		kg[j][j+1]+=-3*ka;

        j=i+3;

		kg[j][j]  +=4*ka;
	}
	}
}