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

#define MAX 100

#define FINE 400

void files(void);
void finish(void);
void freefree(void);
void fixedfree(void);
void fixedfixed(void);
void zero(void);
void torque_core(void);

void freefree_engine(void);
void fixedfree_engine(void);
void fixedfixed_engine(void);

void choice1(void);
void choice2(void);

void middle_entry(void);

void printvector(void);

	long ij;
	double dd;

int i,ibc,idisks;
int jk;
int kv;
int n,nk;

int ishafts;

int iflag =999;

double k[MAX],j[MAX];

double top=0.;

double delta;
double omega,omegab,omegabb;
double om,omb,ombb;

double sum;

double T,Tb,temp;

double x[MAX];

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

double f1,oct;

double omega_before;

char filename[5][FILENAME_MAX];
FILE *pFile[5];

char date_string[30];
char ver_string[30];

double maxfreq;

double lbmft2_per_kgm2 = 23.734;
double lbmin2_per_kgm2 = 3417.8;

double inlbfrad_per_Nmrad = 8.8504;
double ftlbfrad_per_Nmrad = 0.73753;


long NT=32000;

// 23.734 lbm ft^2 = kg m^2

long i_inertia_unit;
long i_stiffness_unit;

double inertia_scale;
double stiffness_scale;

char inertia_label[30];
char stiffness_label[30];

void main()
{
	files();


	strcpy(date_string,"December 21, 2005");
     strcpy(ver_string,"Holzer ver 2.1");

	printf("\n\n\n %s  %s \n",ver_string,date_string);
	printf("\n By Tom Irvine \n Email:  tomirvine@aol.com");
	printf("\n\n\n Holzer's Method for Torsional Systems \n");

	fprintf(pFile[1],"\n\n\n %s  %s \n",ver_string,date_string);
	fprintf(pFile[1],"\n By Tom Irvine \n Email:  tomirvine@aol.com");
	fprintf(pFile[1],"\n\n\n Holzer's Method for Torsional Systems \n");

  
	i_inertia_unit=0;


	while(1)
	{
		if( i_inertia_unit==1 || i_inertia_unit==2 || i_inertia_unit==3){break;}

		printf("\n\n Select inertia unit:  ");
		printf("\n 1=lbm ft^2   2=lbm in^2   3=kg m^2 \n");
		scanf("%d",&i_inertia_unit);
    }


	if(i_inertia_unit==1)
	{
		inertia_scale=1./lbmft2_per_kgm2;
		strcpy(inertia_label,"(lbm ft^2)");

	}
	if(i_inertia_unit==2)
	{
		inertia_scale=1./lbmin2_per_kgm2;
		strcpy(inertia_label,"(lbm in^2)");
	}
	if(i_inertia_unit==3)
	{
		inertia_scale=1.;
		strcpy(inertia_label,"(kg m^2)");
	}



	while(1)
	{
		if( i_stiffness_unit==1 || i_stiffness_unit==2 || i_stiffness_unit==3 ){break;}

		printf("\n\n Select stiffness unit:  ");
		printf("\n 1=ft lbf/rad   2=in lbf/rad   3=N m/rad \n");
		scanf("%d",&i_stiffness_unit);
    }
	if(i_stiffness_unit==1)
	{
		stiffness_scale=1./ftlbfrad_per_Nmrad;
		strcpy(stiffness_label,"(ft lbf/rad)");

	}
	if(i_stiffness_unit==2)
	{
		stiffness_scale=1./inlbfrad_per_Nmrad;
		strcpy(stiffness_label,"(in lbf/rad)");
	}
	if(i_stiffness_unit==3)
	{
		stiffness_scale=1.;
		strcpy(stiffness_label,"(N m/rad)");
	}

	

//    while(1)
//	{


		f1=0.001;
		oct=pow(2.,(1./48));

		omega_before=f1;

		printf("\n\n\n Choose Boundary Conditions: \n"); 
		

		printf("\n  1=Free-Free   ");  
		printf("\n  2=Fixed-Free  ");  
		printf("\n  3=Fixed-Fixed \n\n"); 

		scanf("%d",&ibc);

		if(ibc != 1 && ibc !=2 && ibc != 3)
		{
			printf("\n\n Input Error \n");
			finish();
		}
		else
		{

			printf("\n\n Enter number of disks \n");
			scanf("%d",&idisks);

			if(idisks >= MAX)
			{
				printf("\n Input Error \n");
				finish();
			}

		}

		omega=0.;

	
		choice1();
	
//	}


	printf("\n\n Output files: \n");
	printf("\n Holzer.txt - text file");
	printf("\n Holzer.plt - plot file \n");
	
	finish();

}
void finish(void)
{
	printf("\n\n Press any key to exit.\n");
	getch();
	exit(1);
}
void freefree(void)
{
	fprintf(pFile[1],"\n BC:  Free-Free \n");

	for(i=1; i<idisks; i++)
	{
         middle_entry();
    } 
	printf("\n\n Enter inertia of disk  %d (kg m^2)\n",idisks);
	scanf("%lf",&j[idisks]);
	j[idisks]*=inertia_scale;



	for(i=1; i<=idisks; i++)
	{
		x[i]=1.;
    }

    omega = 0.;
   
	kv=1;
	
	printvector();
    
	kv=2;

	omega=f1;

//    printf(" n=%ld  NT=%ld  om%12.4g\n",n,NT,omega);

	for(n=1; n<=NT; n++)
	{


	    omega=omega_before*oct;

		temp=omega;
	 
		freefree_engine();
       
		torque_core();

		omega_before=temp;

		if(omega > 1.0e+05){break;}


		if( kv == idisks+1 && omega > 2.*top )
		{
			break;
		}
	}
}
void freefree_engine(void)
{
		zero();

		x[1] = 1.;

		for(i=2; i<=idisks; i++)
		{
	
	 		sum=0.;
			for( nk = 1; nk <=(i-1) ; nk++)
			{
				sum+=j[nk]*x[nk];
			}

			x[i] = x[i-1] - (pow(omega,2.)/k[i-1])*sum;
//            printf("omega=%12.4g \t x=%12.4g \n",omega,x[i]);
		}

		T=0.;

		for(i=1; i<= idisks; i++)
		{
			T+= j[i]*pow(omega,2.)*x[i];

//             printf("omega=%12.4g \t T=%12.4g \n",omega,T);
		}
//	    printf("%12.4g \t %12.4g \n",omega,T);
//		fprintf(pFile[2],"%12.4g \t %12.4g \n",omega,T);		

//        printf("\n");
//		getch();
}

void fixedfree(void)
{

	fprintf(pFile[1],"\n BC:  Fixed-Free \n");

	printf("\n\n Let disk 1 be the disk at the free end \n");
	fprintf(pFile[1],"\n\n Let disk 1 be the disk at the free end \n");
	

	for(i=1; i<=idisks; i++)
	{    
		middle_entry();         		
    } 

	kv=1;
	for(n=1; n<=NT; n++)
	{
	    omega=omega_before*oct;

		temp=omega;
	 
		fixedfree_engine();
       
		torque_core();

		omega_before=temp;

		if(omega > 1.0e+05){break;}


		if( kv == idisks+1 && omega > 2.*top )
		{
			break;
		}
	}
}
void fixedfree_engine(void)
{
		zero();

	  
		x[1] = 1.;

		for(i=2; i<=(idisks+1); i++)
		{
	
	 		sum=0.;
			for( nk = 1; nk <=(i-1) ; nk++)
			{
				sum+=j[nk]*x[nk];
			}

			x[i] = x[i-1] - (pow(omega,2.)/k[i-1])*sum;
			T=x[i];

		}
 
	   
}
void fixedfixed(void)
{
	fprintf(pFile[1],"\n BC:  Fixed-Fixed \n");


	printf("\n\n Enter stiffness of shaft 1 %s\n",stiffness_label);
	scanf("%lf",&k[1]);
	k[1]*=stiffness_scale;


	for(i=1; i<=idisks; i++)
	{
		printf("\n\n Enter inertia of disk  %d %s \n",i,inertia_label);
		scanf("%lf",&j[i]);    
		j[i]*=inertia_scale;		

		printf("\n\n Enter stiffness of shaft %d %s\n",i+1,stiffness_label);
		scanf("%lf",&k[i+1]);
		k[i+1]*=stiffness_scale;	
	} 
	

    kv=1;
	for(n=1; n<=NT; n++)
	{
	    omega=omega_before*oct;

		temp=omega;
	 
		fixedfixed_engine();
       
		torque_core();

		omega_before=temp;

		if(omega > 1.0e+05){break;}


		if( kv == idisks+1 && omega > 2.*top )
		{
			break;
		}
	}


}
void fixedfixed_engine(void)
{
		zero();

		x[1] = 1.;

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

			sum=0.;
			for( nk = 1; nk <=i ; nk++)
			{
				sum+=j[nk]*x[nk];
			}

	
			x[i+1] = x[i] + (1./k[i+1]) * ( k[1]*x[1] - (pow(omega,2.)*sum ));
            
        }

		T = x[idisks+1];


}

void files(void)
{
	strcpy(filename[1],"Holzer.txt");
	pFile[1]=fopen(filename[1], "w");
		
	strcpy(filename[2],"Holzer.plt");
	pFile[2]=fopen(filename[2], "w");

}
void zero(void)
{
	int ijk;

	for(ijk = 0; ijk < MAX; ijk++)
	{
		x[ijk]=0.;
	}
}
void printvector(void)
{

   if(kv==1)
   {
	   if(ibc==1)
	{
		ishafts=idisks-1;
	}
	if(ibc==2)
	{
		ishafts=idisks;
	}
	if(ibc==3)
	{
		ishafts=idisks+1;
	}


	for(i=1; i<=idisks; i++)
	{
		fprintf(pFile[1],"\n Disk %d inertia\n",i);
        fprintf(pFile[1],"\n = %9.4g kg m^2",j[i]);
        fprintf(pFile[1],"\n = %9.4g lbm in^2",j[i]*lbmin2_per_kgm2 );
        fprintf(pFile[1],"\n = %9.4g lbm ft^2 \n",j[i]*lbmft2_per_kgm2 );
    } 
  
	for(i=1; i<=ishafts; i++)
	{
		fprintf(pFile[1],"\n Shaft %d stiffness\n",i);
        fprintf(pFile[1],"\n = %9.4g N m/rad",k[i]);
        fprintf(pFile[1],"\n = %9.4g in lbf/rad",k[i]*inlbfrad_per_Nmrad );
        fprintf(pFile[1],"\n = %9.4g ft lbf/rad \n",k[i]*ftlbfrad_per_Nmrad );
    } 

    printf("\n\n");

   }

   if(kv <= idisks)
   {

		int iv;

		printf("\n omega %d = %9.4g rad/sec ",kv, omega);
		 printf("\n         = %9.4g Hz \n", omega/TPI);


		fprintf(pFile[1],"\n omega %d = %9.4g rad/sec ",kv, omega);
		 fprintf(pFile[1],"\n         = %9.4g Hz \n", omega/TPI);



		fprintf(pFile[1],"\n\n Displacement Vector %d\n",kv);

		for(iv=1; iv<=idisks; iv++)
		{
			fprintf(pFile[1],"\n %12.4g",x[iv]);
		}
		fprintf(pFile[1],"\n\n\n");
   }
   top=omega;
}
void torque_core(void)
{

	    double temp2=T;

	    fprintf(pFile[2],"%12.4g \t %12.4g \n\n",omega,T);

		if(n>=2 && kv<=idisks)
		{
			if(T*Tb < 0. )
			{

                double Trecord =9.0e+100;
				double omrecord = 0.;

                delta=omega-omega_before; 

//				printf(" o=%12.4g  ob=%12.4g  del=%12.4g\n",omega,omega_before,delta);

				omega=omega_before;

				dd=delta/FINE;

				for(ij=0; ij<=FINE; ij++)
				{
					choice2();


					if(fabs(T)<Trecord)
					{
						Trecord=fabs(T);
					    omrecord = omega;
			 
					}

					omega+=dd;
				}

					
				omega=omrecord;

                choice2();

				printvector();

				kv++;

			}
			if(T*Tb == 0. && Tb !=0.)
			{

				printf("z  om=%12.4g  T=%12.4g\n",omega,T);

		     	printvector();

				kv++;

			}
	
		}

	
		
		Tb=temp2;


}
void choice1(void)
{

	if(ibc == 1)
	{
		freefree();
	}
	if(ibc == 2)
	{
		fixedfree();
	}
	if(ibc == 3)
	{
		fixedfixed();
	}
}
void choice2(void)
{

	if(ibc == 1)
	{
		freefree_engine();
	}
	if(ibc == 2)
	{
		fixedfree_engine();
	}
	if(ibc == 3)
	{
		fixedfixed_engine();
	}
}
void middle_entry()
{
		printf("\n\n Enter inertia of disk  %d %s\n",i,inertia_label);
		scanf("%lf",&j[i]);
		j[i]*=inertia_scale;

		
		printf("\n\n Enter stiffness of shaft %d %s\n",i,stiffness_label);
		scanf("%lf",&k[i]);
		k[i]*=stiffness_scale;

}