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


//************ complex number handling ********************
typedef struct 
{
    double r,i;
	
}doublecomplex;

#define  COMPLEXMAG(x,y,z) z=sqrt( pow( x,2) + pow(y,2) );
#define COMPLEXMAG2(x,y,z) z=    ( pow( x,2) + pow(y,2) );


#define COMPLEXRECIP(ar, ai, xr, xi)              \
{                                                 \
           xr=  ar/( pow( ar,2) + pow(ai,2) );    \
           xi= -ai/( pow( ar,2) + pow(ai,2) );    \
}                                                 \

#define COMPLEXMUL(ar, ai, br, bi,cr,ci)          \
{                                                 \
           cr=  ar*br -ai*bi;                     \
           ci=  ar*bi +ai*br;                     \
}                                                 \

#define COMPLEXDIV(nr, ni, dr, di,cr,ci)          \
{                                                 \
           cr =  nr*dr +ni*di;                    \
           cr/= ( pow( dr,2) + pow(di,2) );       \
           ci =  ni*dr -nr*di;                    \
           ci/= ( pow( dr,2) + pow(di,2) );       \
}                                                 \

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

void files(void);
void input_data(void);

doublecomplex H;
doublecomplex d;
doublecomplex n;

double damp,fn,f;

double k,m;

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

double df,rho,Q;

long i,idamp;

int iunits,imethod;

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


void main()
{

	  printf("\n\n force_transfer.cpp   version 1.1   July 27, 2005 \n");
	  printf("\n By Tom Irvine   Email:  tomirvine@aol.com \n");

	  printf("\n This program generates a set of transfer functions ");
	  printf("\n for a single-degree-of-freedom system subjected to ");
	  printf("\n an excitation force.\n");


	  input_data();

      files();


	  double fmax = 8*fn;

      double z;

	  double omega,omegan;
	  double omega2,omegan2;

	  double phase;

      i=1;

	  while(1)
	  {

			f=(i-1)*df;

			rho = f/fn;

			 omega= tpi*f;

			omegan= tpi*fn;

			 omega2=pow(omega,2);
			omegan2=pow(omegan,2);

// nd compliance  kx/F

			d.r = (1.-pow(rho,2.));
			d.i = 2*damp*rho;

			COMPLEXDIV(1., 0., d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[1]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[2]," %12.6e %12.6e %12.6e \n",f,z,phase);

//  compliance  x/F

			n.r = omegan2/k;
			n.i = 0.;
			d.r = omegan2-omega2;
			d.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[3]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[4]," %12.6e %12.6e %12.6e \n",f,z,phase);

//  dynamic stiffness  F/x

	
			d.r = omegan2/k;
			d.i = 0.;
			n.r = omegan2-omega2;
			n.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[5]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[6]," %12.6e %12.6e %12.6e \n",f,z,phase);


//  mobility   v/F

			n.r = 0.;
			n.i = omega*omegan2/k;
			d.r = omegan2-omega2;
			d.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[7]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[8]," %12.6e %12.6e %12.6e \n",f,z,phase);


//  mechanical impedance  F/v

			d.r = 0.;
			d.i = omega*omegan2/k;
			n.r = omegan2-omega2;
			n.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[9]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[10]," %12.6e %12.6e %12.6e \n",f,z,phase);


//  accelerance  a/F

			n.r = -omega2*omegan2/k;
			n.i = 0.;
			d.r = omegan2-omega2;
			d.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[11]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[12]," %12.6e %12.6e %12.6e \n",f,z,phase);


//  apparent mass  F/a


			d.r = -omega2*omegan2/k;
			d.i = 0.;
			n.r = omegan2-omega2;
			n.i = 2*damp*omega*omegan;

			COMPLEXDIV(n.r, n.i, d.r, d.i,H.r,H.i);
			COMPLEXMAG(H.r,H.i,z)

			phase = -atan2(H.i,H.r)*180./pi;

			fprintf(pFile[13]," %12.6e %12.6e %12.6e \n",f,H.r,H.i);
			fprintf(pFile[14]," %12.6e %12.6e %12.6e \n",f,z,phase);

// **

			i++;

			if( f > fmax ){break;}

	  }

      printf("\n\n Output Files: \n\n");

      	 
	  for(i=1;i<=14;i+=2)
	  {
			printf("  %s \n",filename[i]);
			printf("  %s \n\n",filename[i+1]);
	  }	


	  printf("\n (nd = non-dimensional) \n");

	  printf("\n Press any key to exit");
	  getch();
}
void files(void)
{
	  strcpy(filename[1],"nd_compliance_complex.out");
	  strcpy(filename[2],"nd_compliance_magnitude_phase.out");
      
	  strcpy(filename[3],"compliance_complex.out");
	  strcpy(filename[4],"compliance_magnitude_phase.out");
	  strcpy(filename[5],"dynamic_stiffness_complex.out");
	  strcpy(filename[6],"dynamic_stiffness_magnitude_phase.out");


	  strcpy(filename[7],"mobility_complex.out");
	  strcpy(filename[8],"mobility_magnitude_phase.out");
	  strcpy(filename[9],"mechanical_impedance_complex.out");
	  strcpy(filename[10],"mechanical_impedance_magnitude_phase.out");


	  strcpy(filename[11],"accelerance_complex.out");
	  strcpy(filename[12],"accelerance_magnitude_phase.out");
	  strcpy(filename[13],"apparent_mass_complex.out");
	  strcpy(filename[14],"apparent_mass_magnitude_phase.out");

	 
	  for(i=1;i<=14;i++)
	  {
			pFile[i]=fopen(filename[i],"w");
	  }	

}
void input_data(void)
{
	  printf("\n Select units             ");
	  printf("\n 1=English   2=metric    \n");
	  scanf("%d",&iunits);

	  if(iunits !=1 && iunits !=2 ){iunits=1;}

	  printf("\n\n Select input data method:      ");
	  printf("\n 1 = mass & stiffness             ");
	  printf("\n 2 = mass & natural frequency     ");
	  printf("\n 3 = stiffness & natural frequency \n");
	  scanf("%d",&imethod);

	  if(imethod !=1 && imethod !=2 && imethod !=3){imethod=1;}


	  if(iunits==1)
	  {
		if(imethod==1)
		{
			printf("\n Enter mass (lbm) \n");
			scanf("%lf",&m);
			m/=386.;
			printf("\n Enter stiffness (lbf/in) \n");
			scanf("%lf",&k);

			fn=sqrt(k/m)/tpi;
		}
		if(imethod==2)
		{
			printf("\n Enter mass (lbm) \n");
			scanf("%lf",&m);
			m/=386.;
			printf("\n Enter natural frequency (Hz) \n");
			scanf("%lf",&fn);

			k=m*pow((tpi*fn),2);
		}
		if(imethod==3)
		{
			printf("\n Enter stiffness (lbf/in) \n");
			scanf("%lf",&k);
			printf("\n Enter natural frequency (Hz) \n");
			scanf("%lf",&fn);

			m=k/pow((tpi*fn),2);
			m/=386.;			
		}
	  }

	  if(iunits==2)
	  {
		if(imethod==1)
		{
			printf("\n Enter mass (kg) \n");
			scanf("%lf",&m);
			printf("\n Enter stiffness (N/m) \n");
			scanf("%lf",&k);

			fn=sqrt(k/m)/tpi;
		}
		if(imethod==2)
		{
			printf("\n Enter mass (kg) \n");
			scanf("%lf",&m);
			printf("\n Enter natural frequency (Hz) \n");
			scanf("%lf",&fn);

			k=m*pow((tpi*fn),2);
		}
		if(imethod==3)
		{
			printf("\n Enter stiffness (N/m) \n");
			scanf("%lf",&k);
			printf("\n Enter natural frequency (Hz) \n");
			scanf("%lf",&fn);

			m=k/pow((tpi*fn),2);
		}
	  }
 
	  printf("\n Enter damping format. \n 1= damping ratio   2= Q \n");
	  scanf("%ld",&idamp);

	  if(idamp == 1 )
	  {
	     printf("\n Enter damping ratio (typically 0.05) ");
	      scanf("%lf",&damp);
	  }
	  else
	  {
	     printf("\n Enter amplification factor Q (typically 10) ");
	     scanf("%lf",&Q);
	     damp = 1./(2.*Q);
	  }

	  if(damp >= 1.0)
	  {
		  printf("\n Error:  damping value must be < 1 \n");
		
	  }
	
	  printf("\n\n Enter the frequency step (Hz)  ( Typically <= 1 ) \n");
	  scanf("%lf",&df);


	  printf("\n\n\n Input Summary:\n\n");
	  
	  if(iunits==1)
	  {
			printf(" m=%8.4g lbm     \n",m*386.);
			printf(" k=%8.4g lbf/in  \n",k);
	  }
	  if(iunits==2)
	  {
			printf(" m=%8.4g kg     \n",m);
			printf(" k=%8.4g N/m    \n",k);
	  }
	  printf("\n     fn=%8.4g Hz      \n",fn);
	     printf("      Q=%8.4g Hz      \n",Q);
	     printf("   damp=%8.4g Hz      \n",damp);


}