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



#define EXIT \
printf("\n\n Press any key to exit.\n"); \
getch(); \
printf("\n"); \
exit(1); \


#define PRINTPSD \
printf("\n        Freq      PSD         Slope      Slope \n        (Hz)      (G^2/Hz)     (N)       (dB/octave) \n"); \
for(i=0; i<nbreak; i++)\
{\
if( i < (nbreak-1))\
{\
s1=log(apsd_sam[i+1]/apsd_sam[i])/log(f_sam[i+1]/f_sam[i]);\
s2=s1*10.*log10(2.);\
printf("%12.2f %12.4f %12.4f %12.4f\n",f_sam[i],apsd_sam[i],s1,s2 );\
}\
else\
{\
printf("%12.2f %12.4f\n",f_sam[i],apsd_sam[i],s1,s2 );}\
}\




#define ABS_SLOPE \
double sss=0.;\
double fr;\
for(i=0; i< (nbreak-1); i++)\
{\
fr=f_sam[i+1]/f_sam[i];\
sss=0.;\
sss=log(apsd_sam[i+1]/apsd_sam[i])/log(fr);\
if(sss > slopec)\
{\
apsd_sam[i+1]=apsd_sam[i]*pow(fr,slopec);\
}\
if(sss < -slopec)\
{\
apsd_sam[i+1]=apsd_sam[i]*pow(fr,-slopec);\
}\
}\
if(initial==1 && apsd_sam[0] > apsd_sam[1])\
{\
apsd_sam[0]=apsd_sam[1];\
}\





#define FSORT \
for(i=1;i<=(nbreak-2);i++)\
{\
for(j=i+1; j<=(nbreak-1); j++)\
{\
if(f_sam[j] < f_sam[i] )\
{\
temp=f_sam[j];\
f_sam[j]=f_sam[i];\
f_sam[i]=temp;\
}\
}\
}\





#define FSPACE \
double ialarm = 0;\
double fnum = (f_ref[n_ref-1]-f_ref[0])/200.;\
for(i=0; i<= (nbreak-2); i++)\
{\
if (  fabs(f_sam[i+1] - f_sam[i] ) < fnum )\
{\
ialarm = 1;\
break;\
}\
}\
if(ialarm == 1 )\
{\
double df=(f_ref[n_ref-1]-f_ref[0])/(n_ref-1);\
for(i=1; i<=(nbreak-2); i++)\
{\
f_sam[i]=f_sam[i-1]+df;\
}\
}\
f_sam[nbreak-1]=f_ref[n_ref-1];\
FSORT\



#define PRINTRMS \
printf("\n   Trial: drms=%10.4g  vrms=%10.4g  grms=%10.4g ",drms,vrms,grms); \
printf("\n  Record: drms=%10.4g  vrms=%10.4g  grms=%10.4g \n",drmslow,vrmslow,grmslow); \





void read_input_vrs();

void make_input_vrs();



void read_input_parameters();



void generate_sample_psd();

void generate_sample_psd2();



void interp_sam();

void vrs_sample();

void compare();

void scalepsd();

void grms_sam();

void writefile();

void checklow();

void initzero();

void writebest();

void vrms_sam();

void drms_sam();

void vrs_best();

void interp_best();



FILE *pFile[10];

char filename[10][FILENAME_MAX];



const double max=32767.;

const double pi=atan2(0.,-1.);

const double tpi=2.*pi;

double plow=1.0e+099;

double dam;



const double octave=-1.+pow(2.,(1./24.));





const long MAX=10000;

const long MAX2=100000;



long i,ik,j;



long nbreak;  // number of breakpoints

long ntrials; // number of trails

long ntrials2;

long n_ref;   // number of reference coordinates



long ikbest;



double scale;

double grms;

double grmslow=1.0e+50;

double vrms;

double vrmslow=1.0e+50;

double drms;

double drmslow=1.0e+50;



double f_ref[MAX];   // reference natural frequency

double a_ref[MAX];   // reference vrs(GRMS)





double fin[MAX];     // input frequency

double vrsin[MAX];	 // input PSD

double inslope[MAX]; // input slope



double interp_vrsin[MAX];  // interpolated input PSD





double f_sam[MAX];  // frequency of sample breakpoint

double apsd_sam[MAX];  // acceleration PSD amplitude of sample breakpoint

double slope[MAX];



double f_samfine[MAX];   // frequency of sample breakpoint, interpolated

double apsd_samfine[MAX];   // acceleration PSD of sample breakpoint, interpolated

double opsd[MAX];

double vrs_samfine[MAX];



double xf[MAX];

double xapsd[MAX];

double xvrs[MAX];



double xffine[MAX];

double xapsdfine[MAX];

double xslope[MAX];



double vpsd_sam[MAX];


double s;
double ra;

double slopec;

double s1,s2;



double rho;

double tden,tnum;

double sum;



double f1,f2;



int ipeak;



long ic;

long initial;


long nin = 0;

long nnn=0;




int main()
{


	strcpy(filename[1],"big.out");

	pFile[1]=fopen(filename[1],"w");



	strcpy(filename[6],"input.vrs");

	pFile[6]=fopen(filename[6],"w");



	strcpy(filename[7],"interp_p.in");

	pFile[7]=fopen(filename[7],"w");




	initzero();



	read_input_vrs();        // Read reference VRS from file



	read_input_parameters();  // Read input parameters from keyboard



//	make_input_vrs();		 // Convert the input PSD to a VRS




	for(ik=0; ik<ntrials; ik++)
	{


	   generate_sample_psd();  // Generate the sample psd

	   interp_sam();           // Interpolate the sample psd at 1/24 octave

	   vrs_sample();           // Calculate the vrs of the sample psd

	   compare();              // Compare the sample vrs with the reference vrs

	   scalepsd();              // scale the psd



	   grms_sam();             // calculate the grms value



	   vrms_sam();



	   drms_sam();



	   writefile();           // write the scaled psd and its grms value to a file



       checklow();



	   PRINTRMS



	 //  getch();

	}



	nnn=0;



	for(ik=0; ik<ntrials2; ik++)
	{



	   generate_sample_psd2();  // Generate the sample psd



	   interp_sam();           // Interpolate the sample psd at 1/24 octave



	   vrs_sample();           // Calculate the vrs of the sample psd



	   compare();              // Compare the sample vrs with the reference vrs



	   scalepsd();             // scale the psd



	   grms_sam();             // calculate the grms value



	   vrms_sam();



       drms_sam();



	   writefile();           // write the scaled psd and its grms value to a file



       checklow();



	   PRINTRMS



	//   getch();



	}



	printf("\n\n________________________________________________________________________");



	for( i=0; i < nbreak; i++)
	{

		    f_sam[i]=xf[i];

		 apsd_sam[i]=xapsd[i];


    }

	printf("\n\nOptimum Case\n");





	PRINTPSD



	printf("\n drms=%10.4g  vrms=%10.4g  grms=%10.4g \n",drmslow,vrmslow,grmslow);







	printf("\n________________________________________________________________________");



//    printf("\n\n\n         Input VRS written to file:  input.vrs ");



	printf("\n\n       Optimum PSD written to file:  best.psd ");


	printf("\n\n       Optimum VRS written to file:  best.vrs ");

      printf("\n       This file has both GRMS and Three Sigma Levels. ");



	printf("\n\n\n Please call these files into your own graphics program. ");


    fclose(pFile[1]);
	fclose(pFile[6]);
    fclose(pFile[7]);



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

	getch();





}

void read_input_vrs(void)
{

    double aa,ff;


	printf("\n\n\n\n\n\n\n\n vrs_to_psd.cpp   Ver 1.2     August 18, 2011");
	printf("\n\n by Tom Irvine ");
	printf("\n Email:  tomirvine@aol.com ");

	printf("\n\n This program converts a vibration response spectrum to a");
 	  printf("\n power spectral density.");


	printf("\n\n The optimum envelope is the that which yields the lowest GRMS and ");
	printf("\n Velocity RMS values. ");


	printf("\n\n The acceleration is assumed to be the base excitation for an array");
 	  printf("\n of independent, parallel single-degree-of-freedom systems.");


	printf("\n The natural frequency is an independent variable.");




	printf("\n\n The user must input a vibration response spectrum. ");

	printf("\n The file must be ASCII text with two columns: ");

	printf("\n\n Natural Freq(Hz)   Accel (GRMS  or  Three Sigma ) ");

	printf("\n\n The format is free, but no header lines are allowed. \n");



	printf("\n________________________________________________________________________");



	printf("\n\n Input starting analysis frequency (Hz):  ");

	scanf("%lf",&f1);


	printf("\n Input ending analysis frequency (Hz):  ");

	scanf("%lf",&f2);



	printf("\n\n Enter name of input VRS file: \n");



	scanf("%s",&filename[0]);



	printf("\n\n Enter the amplitude type: ");

    printf("\n 1=GRMS 2= Three Sigma \n");

	scanf("%d",&ipeak);




    int idamp;

    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",&dam);
	}
	else
	{
	     printf("\n Enter amplification factor Q (typically 10) ");
         scanf("%lf",&dam);

	     dam = 1./(2.*dam);
	}




	pFile[0]=fopen(filename[0],"rb");



	if( pFile[0] == NULL )
	{

		printf("\n Input File Error. \n");



		EXIT

	}





	i=0;

    while( fscanf(pFile[0],"%lf %lf", &ff, &aa ) > 0 )
	{



		if(ipeak==2){aa/=3.;}



		long ierror = 0;


		if( ff < 0.)
		{

		    printf("\n Input error:  each frequency must be > 0. ");

			ierror = 999;

		}



		if( ff == 0 )
		{

			ff = 0.1;

		}



		if( ff > 0.)
        {

//           a_ref[i]=aa;

//	       f_ref[i]=ff;



		     fin[i]=ff;

		   vrsin[i]=aa;





		   if( fin[i] < fin[i-1] )
		   {

			   printf("\n Input error:  frequencies must ascend. ");


			   EXIT

		   }




		   if(i > 0 )
		   {

		      inslope[i-1]= log(vrsin[i]/vrsin[i-1])/log(fin[i]/fin[i-1]);

		   }

		   if( aa <= 0. )
		   {

			   printf("\n Input error:  each PSD amplitude must be > 0.");

			   ierror = 999;

		   }



		   if( ierror == 999 )
		   {

			   EXIT

		   }

		   i++;

		   if( i > MAX){break;}

		}



	}

		nin = i;



	printf("\n Number of input coordinates = %ld ",nin);





	// Interpolate Input PSD





	if( f1 < fin[0] ){f1 = fin[0]; }

	if( f2 > fin[nin-1] ){f2 = fin[nin-1]; }


	f_ref[0]=f1;



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

		f_ref[i]=f_ref[i-1]*pow(2.,(1./24.));





		if(	 f_ref[i] > f2 ){break;}



	}

	f_ref[i-1] = f2;



    n_ref=i;





	interp_vrsin[0]=vrsin[0];



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

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

		   if( f_ref[i] >= fin[j] && f_ref[i] < fin[j+1] )
		   {

		         interp_vrsin[i]= vrsin[j]*pow((f_ref[i]/fin[j]),inslope[j]);

		   }

		}

	}



	for( i=0; i<n_ref; i++)
	{

		fprintf(pFile[7],"%lf %lf\n",f_ref[i],interp_vrsin[i]);

		a_ref[i]=interp_vrsin[i];

	}



}

void read_input_parameters(void)
{





	printf("\n\n Input number of trials for phase 1 (brute-force phase): ");

    scanf("%ld", &ntrials);



	printf("\n Input number of trials for phase 2 (refinement phase): ");

    scanf("%ld", &ntrials2);





    if(ntrials > MAX2)
	{

		printf("\n The maximum number of trials is 100,000 \n");

		ntrials=MAX2;

	}



   printf("\n\n Enter number of breakpoints for envelope (min=3 max=12): ");

   scanf("%ld",&nbreak);



   if(nbreak < 3  )nbreak=3;

   if(nbreak > 12 )nbreak=12;


   printf("\n\n Constrain slopes?  1=yes 2=no ");

   scanf("%ld",&ic);




   if(ic==1)
   {

	   printf("\n Enter absolute slope limit (dB/octave): ");

	   scanf("%lf",&slopec);

   }
   else
   {

	   slopec = 60.;

   }



   slopec=fabs(slopec);

   slopec=(slopec/10.)/log10(2.);



//   printf("\n\n n limit = %lf \n",slopec);



   printf("\n\n Further constrain initial slope to be positive ?  1=yes 2=no ");

   scanf("%ld",&initial);



}

void generate_sample_psd()
{

//    long index;

	long iflag=99;

	long ijk;



	double temp;



	f_sam[0]=f_ref[0];

	f_sam[nbreak-1]=f_ref[n_ref-1];



       // generate some random numbers for frequency



	   for(ijk=0; ijk < 10000; ijk++)
	   {

	      for(i=1; i<= (nbreak-2); i++)
		  {

//		       index = long( n_ref*rand()/max );

//			   if(index == n_ref){index=n_ref-2;}

			  double ww = log10(f1)+ (log10(f2)-log10(f1))*rand()/max;

			   f_sam[i]=pow(10.,ww);

		  }



	   // sort the frequencies





         FSORT




		  // check frequencies for adequate spacing


//		 FSPACE


	   }


	   // generate some random numbers for amplitude


	   for(i=0; i< nbreak; i++)
	   {

		    apsd_sam[i]=rand()/max;

	   }





	   double amin=1.0e-12;



	   for(i=0; i< nbreak; i++)
	   {

	      if(apsd_sam[i] < amin){apsd_sam[i]=amin;}

	      if(apsd_sam[i] > (1./amin)){apsd_sam[i]=(1./amin);}

	   }







    // absolute slopes are limited to slopec



	ABS_SLOPE



	printf("\n\n Phase 1, Trial %ld, PSD Coordinates \n",ik);



	PRINTPSD







}

void interp_sam(void)
{


	f_samfine[0]=f_sam[0];

	apsd_samfine[0]=apsd_sam[0];





	for(i=0; i< (nbreak-1);i++)
	{

		slope[i]=log(apsd_sam[i+1]/apsd_sam[i])/log(f_sam[i+1]/f_sam[i]);

	}



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

		f_samfine[i]=f_ref[i];



        for(j=0; j < (nbreak-1); j++)
		{

			if( ( f_samfine[i] >= f_sam[j] ) &&  ( f_samfine[i] <= f_sam[j+1] )  )
			{

				apsd_samfine[i]=apsd_sam[j]*pow( ( f_samfine[i] / f_sam[j] ), slope[j] );

				break;

			}

		}



	}



}

void interp_best(void)
{

    xffine[0]=xf[0];

	xapsdfine[0]=xapsd[0];



	for(i=0; i< (nbreak-1);i++)
	{

		xslope[i]=log(xapsd[i+1]/xapsd[i])/log(xf[i+1]/xf[i]);

	}



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

		xffine[i]=f_ref[i];



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

			if( ( xffine[i] >= xf[j] ) &&  ( xffine[i] <= xf[j+1] )  )
			{

				xapsdfine[i]=xapsd[j]*pow( ( xffine[i] / xf[j] ), xslope[j] );

				break;

			}

		}



	}



}

void vrs_sample(void)
{


	double tdr, trans;



    for(i=0; i< n_ref; i++)
	{

		sum=0.;



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

              // f_ref[i] is the natural frequency

			  // f_ref[j] is the forcing frequency



	          rho=f_ref[j]/f_ref[i];



			  tdr=2.*dam*rho;



			  tden=pow((1.-pow(rho,2)),2)+ pow(tdr,2);

			  tnum=1.+pow(tdr,2);



			  trans=tnum/tden;


			  opsd[j]=trans*apsd_samfine[j];

		}

        grms=0.;

        ra=0.;

        s=0.;

        for( j=0; j < n_ref; j++)
        {
			  s=log( opsd[j+1]/opsd[j] )/log( f_ref[j+1]/f_ref[j] );

              if(s < -1.0001 ||  s > -0.9999 )
			  {
                 ra+= ( opsd[j+1] * f_ref[j+1]- opsd[j]*f_ref[j])/( s+1.);
			  }
		      else
			  {
                 ra+= opsd[j]*f_ref[j]*log( f_ref[j+1]/f_ref[j]);
			  }
            grms=sqrt(ra);
        }

		vrs_samfine[i]=grms;

	}

}

void compare(void)
{

	scale=0.;



	for(i=0; i< n_ref; i++)
	{



		    if( vrs_samfine[i] < 1.0e-20)
			{

			    printf("\n error ");

			    exit(1);

			}



            if(  (a_ref[i]/vrs_samfine[i]) > scale )
			{

			    scale=(a_ref[i]/vrs_samfine[i]);

			}


	}


}

void scalepsd()
{

	scale=pow(scale,2.);

	for(i=0; i< nbreak; i++)
	{

		apsd_sam[i]*=scale;

	}


}

void grms_sam()
{

	      double ra=0.;

		  double s[MAX];



          grms=0.;



		  for(i=0; i < MAX; i++)

		  {

			  s[i]=0.;

		  }





	      for( i=0; i < nbreak-1; i++)
		  {



			  s[i]=log( apsd_sam[i+1]/apsd_sam[i] )/log( f_sam[i+1]/f_sam[i] );



              if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			  {

                 ra+= ( apsd_sam[i+1] * f_sam[i+1]- apsd_sam[i]*f_sam[i])/( s[i]+1.);

			  }

		      else
			  {

                 ra+= apsd_sam[i]*f_sam[i]*log( f_sam[i+1]/f_sam[i]);

			  }

		  }

          grms=sqrt(ra);



}





void writefile()
{


    fprintf(pFile[1],"\n%ld %12.5g %12.5g\n",ik,vrms,grms );



	for( i=0; i < nbreak; i++)
	{

         fprintf(pFile[1],"%12.5e %12.5e\n",f_sam[i],apsd_sam[i]);

    }



}

void writebest()
{



	for( i=0; i < nbreak; i++)
	{

         fprintf(pFile[2],"%12.5e %12.5e\n",xf[i],xapsd[i]);



		    f_sam[i]=xf[i];

		 apsd_sam[i]=xapsd[i];



    }

	for( i=0; i < n_ref; i++)
	{

         fprintf(pFile[3],"%12.5e %12.5e %12.5e\n",xffine[i],xvrs[i],3.*xvrs[i]);

    }



}

void checklow()
{

	   if( (drms*vrms*grms) < plow)
	   {

		//   printf("\n\n New Optimum PSD Found.\n");

		   plow=(drms*vrms*grms);


		   drmslow=drms;

		   vrmslow=vrms;

		   grmslow=grms;



		   ikbest=ik;



		   nnn=ntrials;



		   	for( i=0; i < nbreak; i++)
			{

                xf[i]=f_sam[i];

				xapsd[i]=apsd_sam[i];

            }

			for( i=0; i < (nbreak-1); i++)
			{

                xslope[i]=slope[i];

            }



            strcpy(filename[2],"best.psd");

	        pFile[2]=fopen(filename[2],"w");



	        strcpy(filename[3],"best.vrs");

	        pFile[3]=fopen(filename[3],"w");



	        interp_best();        // Interpolate the best psd at 1/24 octave

	        vrs_best();           // Calculate the vrs of the best psd



	        writebest();



            fclose(pFile[2]);

	        fclose(pFile[3]);



	   }

}

void initzero()
{

	for(i=0; i<nbreak; i++)
	{

	  f_ref[i]=0.;

      a_ref[i]=0.;



      f_sam[i]=0.;

      apsd_sam[i]=0.;

      slope[i]=0.;



      f_samfine[i]=0.;

      apsd_samfine[i]=0.;



      vrs_samfine[i]=0.;



      xf[i]=0.;

      xapsd[i]=0.;

    }

}

void vrms_sam()
{

	      double ra=0.;

		  double s[MAX];




          double conv=pow(386.,2.)/(4.*pow(pi,2.));



          vrms=0.;



		  for(i=0; i < MAX; i++)
		  {

			  s[i]=0.;

		  }





		  for(i=0; i < nbreak; i++)
		  {

             vpsd_sam[i]=apsd_sam[i]*conv/pow(f_sam[i],2.);

          }



	      for( i=0; i < nbreak-1; i++)
		  {



			  s[i]=log( vpsd_sam[i+1]/vpsd_sam[i] )/log( f_sam[i+1]/f_sam[i] );



              if(s[i] < -1.0001 ||  s[i] > -0.9999 )

			  {

                 ra+= ( vpsd_sam[i+1] * f_sam[i+1]- vpsd_sam[i]*f_sam[i])/( s[i]+1.);

			  }

		      else

			  {

                 ra+= vpsd_sam[i]*f_sam[i]*log( f_sam[i+1]/f_sam[i]);

			  }

		  }

          vrms=sqrt(ra);



}

void drms_sam()
{

	      double ra=0.;

		  double s[MAX];



		  double dpsd_sam[MAX];

          double conv=1./(4.*pow(pi,2.));



          drms=0.;



		  for(i=0; i < MAX; i++)
		  {

			  s[i]=0.;

		  }





		  for(i=0; i < nbreak; i++)
		  {

             dpsd_sam[i]=vpsd_sam[i]*conv/pow(f_sam[i],2.);

          }


	      for( i=0; i < nbreak-1; i++)
		  {


			  s[i]=log( dpsd_sam[i+1]/dpsd_sam[i] )/log( f_sam[i+1]/f_sam[i] );


              if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			  {

                 ra+= ( dpsd_sam[i+1] * f_sam[i+1]- dpsd_sam[i]*f_sam[i])/( s[i]+1.);

			  }
		      else
			  {

                 ra+= dpsd_sam[i]*f_sam[i]*log( f_sam[i+1]/f_sam[i]);

			  }

		  }

          drms=sqrt(ra);



}

void vrs_best(void)
{


	double tdr, trans;

    for(i=0; i< n_ref; i++)
	{

		sum=0.;



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

              // f_ref[i] is the natural frequency

			  // f_ref[j] is the forcing frequency


	          rho=f_ref[j]/f_ref[i];


			  tdr=2.*dam*rho;



			  tden=pow((1.-pow(rho,2)),2)+ pow(tdr,2);

			  tnum=1.+pow(tdr,2);



			  trans=tnum/tden;

			  opsd[j]=trans*xapsdfine[j];


		}

        grms=0.;

        ra=0.;

        s=0.;

        for( j=0; j < n_ref; j++)
        {
			  s=log( opsd[j+1]/opsd[j] )/log( f_ref[j+1]/f_ref[j] );

              if(s < -1.0001 ||  s > -0.9999 )
			  {
                 ra+= ( opsd[j+1] * f_ref[j+1]- opsd[j]*f_ref[j])/( s+1.);
			  }
		      else
			  {
                 ra+= opsd[j]*f_ref[j]*log( f_ref[j+1]/f_ref[j]);
			  }
            grms=sqrt(ra);
        }

		xvrs[i]=grms;

	}

}

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

void generate_sample_psd2()

{



	double temp;




	f_sam[0]=f_ref[0];

	f_sam[nbreak-1]=f_ref[n_ref-1];



	for(i=1; i<= (nbreak-2); i++)
	{
		f_sam[i]=xf[i]*(0.9+0.2*rand()/max);
	}


	    // sort the frequencies


          FSORT

//		  FSPACE
/*
	for(i=0; i<= (nbreak-1); i++)
	{
		printf(" %8.4g \t %8.4g \n",f_sam[i],xf[i]);
	}
	getch();
*/


	   // generate some random number for amplitude




	   double db;


	   for(i=0; i< nbreak; i++)
	   {


			db=12.*(1.-(double(ik)/double(ntrials2)));

			db*=rand()/max;



			if( rand()/max < 0.5 ){db*=-1.;}



			apsd_sam[i]=xapsd[i]*pow(10.,(db/10.));

	   }



    // add constraint



	ABS_SLOPE



    printf("\n\n Phase 2, Trial %ld, PSD Coordinates \n",ik);


	PRINTPSD



}

void make_input_vrs(void)
{



	double tdr, trans;



    for(i=0; i< n_ref; i++)
	{

		sum=0.;



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

              // f_ref[i] is the natural frequency

			  // f_ref[j] is the forcing frequency



	          rho=f_ref[j]/f_ref[i];



			  tdr=2.*dam*rho;



			  tden=pow((1.-pow(rho,2)),2)+ pow(tdr,2);

			  tnum=1.+pow(tdr,2);



			  trans=tnum/tden;



              double dfi=f_ref[j]*octave;



			  opsd[j]=trans*interp_vrsin[j];

		}
        grms=0.;

        ra=0.;

        s=0.;

        for( j=0; j < n_ref; j++)
        {
			  s=log( opsd[j+1]/opsd[j] )/log( f_ref[j+1]/f_ref[j] );

              if(s < -1.0001 ||  s > -0.9999 )
			  {
                 ra+= ( opsd[j+1] * f_ref[j+1]- opsd[j]*f_ref[j])/( s+1.);
			  }
		      else
			  {
                 ra+= opsd[j]*f_ref[j]*log( f_ref[j+1]/f_ref[j]);
			  }
            grms=sqrt(ra);
        }

		a_ref[i]=grms;

	    fprintf(pFile[6]," %12.4g \t %12.4g\n",f_ref[i],a_ref[i]  );

	}
}