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


#define MAX 10000

#define NUMX 90

#define NUMX1 91

// #define NFREQ 11000

#define NFREQ 30

#define LAST_F_N 27


void record_check();

void normalized_relative_sound_power_per_unit_core_length();

void normalized_relative_sound_power_spectrum_level();


void inputdata1();
void files();
void printspl();
void directivity();

void material();

void step1();
void step2();
void step3();
void step4();
void step5();
void step6();
void step7();
void step8();	

void blocked();  

void franken();

void method2();


void c1();
void c2();
void c3();
void c4();
void c5();

void source_allocate();

void source_allocated_Eldred_Wilby();

void length_diameter();

void SPL_check();


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

const double ref  = 1.0e-12;
const double Aref = 20.e-06;


const double mm_per_inch = 25.4;

const double inch_per_mm = 1./mm_per_inch;


const double m_per_ft = 0.3048;

const double ft_per_m = 1./m_per_ft;


const double m_per_inch = 0.3048/12.;

const double inch_per_m=1./m_per_inch;


const double N_per_lbf = 4.448;

const double lbf_per_N = 1./N_per_lbf;


const double kgpm3_per_lbmpin3 = 27675;

const double kg_per_lbm = 0.45351;

double xtr;

double iring=0.1;

int iflag;

int irad;
int igeo;
int i,ifm,ilast,imat,ind,isi,isn,isys,iw;
int ijk;
int j,num;
int N;

int ikv;

long ie;

long jk;

int ncoord;

int ncoord_core;

double SPL_sbp[200][200];

double radius[200];

double term1,term2,term3,term4;

double xx[NUMX1];
double yy[NUMX1];

double xxr[100];
double yyr[100];

double xxr_core[100];
double yyr_core[100];

double xx_core[100];
double yy_core[100];


double fc[100];
double fcc[100];

double ssum[100];

double power_spectrum[200][200];

double LWSB[200][200];

double sss[200][200];

double ppp[200][200];

double dx;

double xt;

double sum;

long ie_max;

double normalized_aps_per_length[MAX];

double power_scale;

double slope_limit;
double db_per_octave;

double aceff;
double a,cc;
double az;

double cmat,cref;

double de,dei;


double slope,strouhal,E,F,U;
double WOA,LW;
double amp[MAX];
double samp[MAX];

double sn[MAX];
double nrspl[MAX];
double f[MAX],freq[MAX];
double Lwb[MAX];

double Lwb_reference[MAX];

double sns[MAX];
double asap[MAX];
double x,xs,x1,x2;
double alpha;
double r,rx,ry;
double spl[MAX];

double LWS[200];

double SPL_reference[MAX];

double d1,d2,d3,d4,d5;
double an[12],dn[12];

double beta_deg,beta_rad,rho,theta;
double delta;
double cspeed;
double stationdiam;
double scale;

double DI;
double ZZ;
double W;

double Wp[MAX];

double xm2[MAX];



double fr[20],ar[20];
double al,thickness,vdb; 
double k1,k2,len;
double sfr,va,df;
double s[50],psd[50];

double Lwb_record;
double Lwb_recordmax=1.0e+90;


double SPL_record;
double SPL_recordmax=1.0e+90;

int jjk;

int motor;

int i_end_slope;

const double max=32767.;


char sgeo[10][70];

char motor_type[40];

char prefix[25];

char program_name[50];
char Irvine[50];

char filename[8][FILENAME_MAX];

FILE *pFile[8];



void main()
{


	strcpy(program_name," liftoff_hyrbid.cpp, ver 1.0,	March 20, 2007 ");
	strcpy(Irvine,"  By Tom Irvine   Email: tomirvine@aol.com ");


	printf("\n");

	printf("\n\n %s \n\n",program_name);
	printf("\n  %s \n",Irvine);


    printf("\n\n References:  \n");
	printf("\n 1. Acoustic Loads Generated by the Propulsion System. ");
	printf("\n    NASA SP-8072, Monograph N71-33195, 1971. \n");
	printf("\n 2. Wilby & Wilby, Prediction of External Acoustic Field ");
	printf("\n    on Taurus During Liftoff, AARC No. 121, 1991. \n");


    printf("\n  This is a hybrid approach combining methods 1 and 2 ");
	printf("\n  from Reference 1. \n");

//	printf("\n Eldred Method, page 28. \n\n");


	printf("\n Select system");
	printf("\n 1=English  2=metric \n");
	scanf("%d",&isys);

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

	inputdata1();

//	printf("\n\n Enter number of stations of interest:\n");
//	scanf("%d",&isi);

	isi=1;

	for(ikv=1;ikv<=isi;ikv++)
    {
        printf("\n\n ************** \n");

		printf("\n\n Enter output filename prefix for station %d \n",ikv);
		scanf("%s",&prefix);

		files();

        length_diameter();

		fprintf(pFile[1],"\n Station Diameter = %7.3g in", stationdiam*inch_per_m);
		fprintf(pFile[1],"\n                  = %7.3g m", stationdiam);

		fprintf(pFile[1],"\n\n Station = %12.3g inches ",12.*xs/m_per_ft );
	      fprintf(pFile[1],"\n         = %12.3g meters \n\n",xs);

		step1();
		step2();
		step3();
		step4();
		step5();

		fclose(pFile[0]);

		step6();
	
/*
		printf("\n\n Calculate resulting vibration via the Franken method? "); 

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

		scanf("%d",&ifm);

		if(ifm==1)
		{
			printf("\n Franken upper limit curve is used per NASA CR-1302 \n");

			material();

			stationdiam/=m_per_ft;

			franken();
		}
*/
        fclose(pFile[0]);
		fclose(pFile[1]);
		fclose(pFile[2]);
		fclose(pFile[3]);
		fclose(pFile[4]);
		fclose(pFile[5]);


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

		printf("\n\n  Acoustic text file:  %s \n",filename[1]);	
		printf(  "\n  Acoustic plot file:  %s \n",filename[2]);
		printf(  "\n  Acoustic plot file +3 dB:  %s \n\n",filename[3]);


		if(ifm==1)
		{
			printf("\n  Skin vibration PSD file:  %s \n",filename[4]);
			printf("\n  Skin vibration PSD file +3 dB:  %s \n\n",filename[5]);
		}

	}

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


}
void inputdata1()
{


		printf("\n Select first stage motor:	 \n");
		printf("\n 1=SR-19						   ");
		printf("\n 2=SR-19 with two MLRS strap-ons ");
		printf("\n 3=Castor 4B					   ");
		printf("\n 4=Peacekeeper Stage 1	       ");	
		printf("\n 5=other                       \n");
		
		scanf("%d",&motor);

		strcpy(motor_type,"Other motor");

		if(motor==1)  
		{
			strcpy(motor_type,"SR-19 motor");
		}
		if(motor==2)  
		{
			strcpy(motor_type,"SR-19 with two MLRS strap-ons");
		}
		if(motor==3)  
		{
			strcpy(motor_type,"Castor 4B");
		}
		if(motor==4)  
		{
			strcpy(motor_type,"Peacekeeper Stage 1");
		}


		if(motor==1 || motor==2)  // SR-19
		{
			F=50000.;

			if(isys==2)
			{
				F*=N_per_lbf;
			}
		}
		if(motor==3)  // Castor 4B
		{
			F=1.12e+005; 

			if(isys==2)
			{
				F*=N_per_lbf;			
			}
		}
		if(motor==4)  // Peacekeeper Stage 1
		{
			F=5.70e+005; 

			if(isys==2)
			{
				F*=N_per_lbf;			
			}
		}
		if(motor==5)
		{
			if(isys==1)
			{
				printf("\n Enter thrust (lbf)\n");
			}
			else
			{
				printf("\n Enter thrust (N)\n");
			}

			char sthr[16];
			scanf("%s",&sthr);
			sscanf(sthr,"%lf",&F);
		}

		if(isys==1)
		{
			F*=N_per_lbf; 
		}

		if(motor==1 || motor==2)  // SR-19
		{
			U=9254.;

			if(isys==2)
			{
				U*=m_per_ft;
			}
		}
		if(motor==3)  // Castor 4B
		{
			U=8475.; 

			if(isys==2)
			{
				U*=m_per_ft;			
			}
		}
		if(motor==4)  // Castor 4B
		{
			U=9060.; 

			if(isys==2)
			{
				U*=m_per_ft;			
			}
		}
		if(motor==5)
		{

			if(isys==1)
			{
				printf("\n Enter exhaust exit velocity (ft/sec)\n");
			}
			else
			{
				printf("\n Enter exhaust exit velocity (m/sec)\n");
			}
			char svel[16];
			scanf("%s",&svel);
			sscanf(svel,"%lf",&U);
		}

		if(isys==1)
		{
			U*=m_per_ft; 
		}


		if(isys==1)
		{
			printf("\n Enter speed of sound (ft/sec)\n Typical speed = 1120 ft/sec \n");
		}
		else
		{
			printf("\n Enter speed of sound (m/sec)\n  Typical speed = 343 m/sec \n");
		}
		char ssp[16];
		scanf("%s",&ssp);
		sscanf(ssp,"%lf",&cspeed);


		if(isys==1)
		{
			cspeed*=m_per_ft; 
		}


		printf("\n Enter acoustic efficiency in decimal format.  \n (Conservative value = 0.01) \n");
		scanf("%lf",&aceff);

		printf("\n\n Acoustic efficiency = %5.3g percent \n\n",aceff*100.);


		if(motor==1 || motor==2)  // SR-19
		{
			dei=28.5;

			if(isys==2)
			{
				dei*=m_per_inch;
			}
		}
		if(motor==3)  // Castor 4B
		{
			dei=37.5; 

			if(isys==2)
			{
				dei*=m_per_inch;			
			}
		}
		if(motor==4)  // Peacekeeper Stage 1
		{
			dei=60.9; 

			if(isys==2)
			{
				dei*=m_per_inch;			
			}
		}
		if(motor==5)
		{
			if(isys==1)
			{
				printf("\n Enter nozzle diameter (inches)\n");
			}
			else
			{
				printf("\n Enter nozzle diameter (m)\n");
			}

			char snoz[12];
			scanf("%s",&snoz);
			sscanf(snoz,"%lf",&dei);
		}


//	scanf("%lf",&dei);


		if(isys==1)
		{
			dei*=m_per_ft/12.; 
		}
    

		N=1;

		if(motor==4)
		{
			printf("\n\n Enter number of nozzles \n");
			scanf("%d",&N);
		}

		printf("\n Select volume for sound radiation \n");
		printf("\n 1=hemisphere (flat ground plane) ");
		printf("\n 2=spherical  (free space) \n");
		scanf("%d",&irad);



		printf("\n\n Apply Wilby correction factor for pressure reflections "); 
		  printf("\n at the surface of the vehicle? \n");
		printf("1=yes  2=no \n");
		scanf("%d",&iw);
	

		printf("\n Choose geometry for source allocation\n");

		strcpy(sgeo[1],"single nozzle, undeflected (Smith & Brown)");
		strcpy(sgeo[2],"single nozzle, undeflected (Morgan & Young)");
		strcpy(sgeo[3],"multiple nozzles, undeflected");
		strcpy(sgeo[4],"deflected, open scoop");
		strcpy(sgeo[5],"deflected, closed bucket");
		strcpy(sgeo[6],"deflected, single 45 deg plate");
		strcpy(sgeo[7],"deflected, 90 deg flat plate, conical diffuser, or wedge");


		printf("\n 1=  %s",sgeo[1]);
		printf("\n 2=  %s",sgeo[2]);
		printf("\n 3=  %s",sgeo[3]);
		printf("\n 4=  %s",sgeo[4]);
		printf("\n 5=  %s",sgeo[5]); 
		printf("\n 6=  %s",sgeo[6]);
		printf("\n 7=  %s \n",sgeo[7]);


		scanf("%d",&igeo);

		if(igeo !=1 && igeo !=2 && igeo !=3 && igeo !=4 && igeo !=5 && igeo !=6 && igeo !=7)
		{
			igeo =1;
		}



		if(igeo==4 || igeo==5 || igeo==6 || igeo==7)
		{
			
			if(isys==1)
			{
				printf("\n Enter distance (inches) from nozzle exit to ground. \n");
            
				char sx1[16];
				scanf("%s",&sx1);
				sscanf(sx1,"%lf",&x1);

 
				x1*=m_per_ft/12.; 
			
//			fprintf(pFile[1],"\n x1 = %10.2f meters \n",x1);

			}
			else
			{
				printf("\n Enter distance (meters) from nozzle exit to ground.\n");
	        
				char sx1[16];
				scanf("%s",&sx1);
				sscanf(sx1,"%lf",&x1);

			}

			printf("\n Enter angle (deg) rho from deflected exhaust "); 
			printf("\n flow axis to ground (typical value 10 degrees). \n");
	
			char srho[16];
			scanf("%s",&srho);
			sscanf(srho,"%lf",&rho);

        }


		i_end_slope=1;

		db_per_octave = 6.;   // conservative slope limit

		slope_limit=db_per_octave/3.;

//	fprintf(pFile[1],"\n\n  Conservative smoothing is applied at high frequencies. \n");	

		printf("\n\n Enter number of iteration (suggest at least 10000) \n");
		scanf("%ld",&ie_max);
		printf("\n");

}
void step1()
{
}
void step2()
{
	WOA=0.5*aceff*N*F*U;

//    printf("\n 0.5*aceff = %8.4lf \n",0.5*aceff); 	
	
//	WOA=0.005*N*F*U;

	if(motor==2) // SR-19 with two MLRS strap-ons
	{ 
        double MLRS_vel=7600.*m_per_ft;

		double MLRS_thrust=37000.*N_per_lbf;

        power_scale=((2*MLRS_vel*MLRS_thrust)+F*U)/(F*U); 

		WOA*=power_scale;

		printf("\n power scale = %8.3g \n",power_scale);

	}

    printf("\n\n overall acoustic power WOA = %12.3g Watts\n",WOA); 

}
void step3()
{
	LW=120. + 10.*log10(WOA);

	if(igeo ==7)
	{
		LW-=3.;
	}

    printf("\n\n overall acoustic power level LW = %12.3g dB ",LW);
	printf("\n Ref = 1.0e-12 Watts \n\n");


	if(motor==2) // SR-19 with two MLRS strap-ons
	{ 
		fprintf(pFile[1],"\n The power is increased by %3.1lf dB to account for the strap-on motors. ",10.*log10(power_scale));
         printf("\n The power is increased by %3.1lf dB to account for the strap-on motors. ",10.*log10(power_scale));
	}

    fprintf(pFile[1],"\n\n overall acoustic power level LW = %12.3g dB \n",LW);
	fprintf(pFile[1],"\n Ref = 1.0e-12 Watts \n\n");
	
}
void step4()
{
	de=sqrt(N)*dei;

//	printf("\n de= %8.3g \n",de);
}
void step5()
{   
	// power spectrum


	   sn[0]=ref;
	nrspl[0]=-1.0e+12;


	   sn[1]=0.002;
	nrspl[1]=-1.;


	   sn[2]=0.005;
	nrspl[2]=8.;


	   sn[3]=0.01;
	nrspl[3]=10.;	   


	   sn[4]=0.02;
	nrspl[4]=11.;


	   sn[5]=0.03;
	nrspl[5]=10.5;


	   sn[6]=0.05;
	nrspl[6]=9.;


	   sn[7]=0.1;
	nrspl[7]=6.;


	   sn[8]=0.2;
	nrspl[8]=1.;


	   sn[9]=0.5;
	nrspl[9]=-7.5;


	   sn[10]=1.;
	nrspl[10]=-13.5;


	   sn[11]=2.;
	nrspl[11]=-20.;


	   sn[12]=5.0;
	nrspl[12]=-27.;


	   sn[13]=10000.;
	nrspl[13]=-500.;


	num=13;


	for(i=0; i<=num; i++)
	{
		nrspl[i]=(ref)*pow( 10.,((nrspl[i])/10.));
	}


	freq[0]=20.;
	freq[1]=25.;
	freq[2]=31.5;
	freq[3]=40.;
	freq[4]=50.;
	freq[5]=63.;
	freq[6]=80.;
	freq[7]=100.;
	freq[8]=125.;
	freq[9]=160.;
	freq[10]=200.;
	freq[11]=250.;
	freq[12]=315.;
	freq[13]=400.;
	freq[14]=500.;
	freq[15]=630.;
	freq[16]=800.;
	freq[17]=1000.;
	freq[18]=1250.;
	freq[19]=1600.;
	freq[20]=2000.;
	freq[21]=2500.;
	freq[22]=3150.;
	freq[23]=4000.;
	freq[24]=5000.;
	freq[25]=6300.;
	freq[26]=8000.;
	freq[27]=10000.;
	
	ilast = 27;


	for(i=0; i<= ilast; i++)
	{
		fc[i]=freq[i];

        f[i]=freq[i];
		strouhal = freq[i]*de/U;

		for(j=0; j<=(num-1); j++)
		{
	
			if(strouhal==sn[j])
			{
                amp[i]=nrspl[i];
              
				break;
			}
			if( strouhal>sn[j] && strouhal<sn[j+1] )
			{

				slope=log10(nrspl[j+1]/nrspl[j])/log10(sn[j+1]/sn[j]);

				az=log10(nrspl[j]);			
				az+=slope*(log10(strouhal)-log10(sn[j]));
				
				amp[i]=pow(10.,az);

				break;
			}
		}

	}


	for(i=0; i<= ilast; i++)   // Lwb = Sound Power Level 
	{

        double df = (pow(2.,(1./6.))-1./pow(2.,(1./6.)))*freq[i];

		Lwb[i]= 10.*log10(amp[i]/ref) +LW -10.*log10(U/de) +10.*log10( df);

		if( i>=1 && (Lwb[i-1]-Lwb[i] ) > 2. )
		{
            Lwb[i]=Lwb[i-1]-2;
		}

		Lwb_reference[i]=Lwb[i];
	
		printspl();
	}

//	printf("\n\n Step 5 acoustic power spectrum written to step5.spl \n");

}
void step6()
{


	if(igeo==1)  //single nozzle, undeflected (Smith & Brown)  Fig 14 from SP-8072
	{

		 sns[0]=ref;
		 sns[1]=1.0e-03;
		 sns[2]=2.0;
		 sns[3]=1.0e+12;

		 asap[0]=90.;
		 asap[1]=90.;
		 asap[2]=17.;
		 asap[3]=17.*0.999;

		 isn=3;
		
	}
	if(igeo==2)  //single nozzle, undeflected (Morgan & Young)  Fig 14 from SP-8072
	{
		 sns[0]=ref;
		 sns[1]=1.0e-03;
		 sns[2]=2.0e-02;
		 sns[3]=5.0e-02;
		 sns[4]=2.0e-00;
		 sns[5]=1.0e+12;

		 asap[0]=30.;
		 asap[1]=30.;
		 asap[2]=26.;
		 asap[3]=18.;
		 asap[4]=10.;
		 asap[5]=10.*0.999;

		 isn=5;
	}
	if(igeo==3) //multiple nozzles, undeflected   Fig 14 from SP-8072
	{
		 sns[0]=ref;
		 sns[1]=1.0e-03;
		 sns[2]=7.0e-03;
		 sns[3]=1.2e-02;
		 sns[4]=1.0e-01;
		 sns[5]=2.0e-00;
		 sns[6]=1.0e+12;

		 asap[0]=90.;
		 asap[1]=90.;
		 asap[2]=45.;
		 asap[3]=35.;
		 asap[4]=17.;
		 asap[5]=10.;
		 asap[6]=10.*0.999;

		 isn=6;
	}
	if(igeo==4)  //deflected, open scoop    Fig 14 from SP-8072
	{
		 sns[0]=ref;
		 sns[1]=1.0e-03;
		 sns[2]=1.0e-02;
	     sns[3]=4.0e-02;
		 sns[4]=1.0e-01;
	     sns[5]=5.0e-01;
		 sns[6]=1.0e-00;
		 sns[7]=2.0e-00;
		 sns[8]=1.0e+12;

		 asap[0]=35.;
		 asap[1]=35.;
		 asap[2]=32.;
		 asap[3]=25.;
		 asap[4]=18.;
		 asap[5]=5.;
		 asap[6]=2.1;
		 asap[7]=0.2;
		 asap[8]=0.2*0.999;

		 isn=8;
	}
	if(igeo==5) //deflected, closed bucket   Fig 14 from SP-8072
	{
		 sns[0]=ref;
		 sns[1]=1.0e-03;
		 sns[2]=5.0e-02;
		 sns[3]=1.0e-01;
		 sns[4]=1.0e-00;
		 sns[5]=1.0e+12;

		 asap[0]=17.;
		 asap[1]=17.;
		 asap[2]=8.;
		 asap[3]=6.;
		 asap[4]=0.3;
		 asap[5]=0.3*0.999;

		 isn=5;
	}
	
	if(igeo !=6 && igeo !=7)
	{
		source_allocate();
    }
	else
	{
		source_allocated_Eldred_Wilby();
	}

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

	step7();

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

	step8();

	method2();


}
void step7()
{


	iflag=0;

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

	    theta=180.;
		
		strouhal = freq[i]*de/U;

		if(igeo==1 || igeo==2 || igeo==3)  // undeflected
		{
			r=xs+samp[i];    // samp[i] is apparent source allocation position
		}
		else   // deflected
		{


			x2= samp[i] -x1;


			if(x2 > 0. )
			{

				ry= (xs+x1)-samp[i]*sin(PI*rho/180.);

				rx= x2*cos(PI*rho/180.);

				r=sqrt( pow(rx,2.) + pow(ry,2.) );

				beta_rad = atan2(ry,rx);
				beta_deg = 180.*beta_rad/PI;

				theta = 180. - rho - beta_deg;   
			}
			else   // special case where source occurs before deflection
			{

//                  xs  =  station length
//              samp[i] =  apparent source allocation position

//				r = xs + samp[i];
//				r += 2.*fabs((x1-fabs(samp[i])));

				
				r = xs + x1;

                theta = 90.-rho;   
			}

		}

		if(theta>180.){theta=180.;}
		if(theta< 20.){theta= 20.;}

		directivity();


        double ab = 8.; // hemisphere (flat ground plane)

		if(irad==2)     // sphere (free space)
		{
			ab = 11;
		}


		spl[i]=Lwb[i]-10.*log10(pow(r,2.)) -ab +DI;   // SPL equation


		delta=0.;
		if(iw==1)
		{
			blocked();
		}

		spl[i]+=delta;

		SPL_reference[i]=spl[i];



//  format

		if(i==0)
		{
			fprintf(pFile[1],"\n\n Zero dB reference: 20 micro Pa \n");

			fprintf(pFile[1],"\n   f(Hz)    Strouhal   Lwb[i]   source[m]   r[m]  theta(deg)  DI(dB)  delta(dB) spl(dB) \n");
		}


/*
		if(i_end_slope==1 && i>=2)   // conservative smoothing of peak;
		{
			if( (spl[i-1]-spl[i-2]) > 0.1 &&  (spl[i]-spl[i-1]) < 0.   )
			{
				spl[i]=spl[i-1];

				iflag=1;
			}
		}


		if(i_end_slope==1 && i>=1)   // conservative smoothing based on slope limit
		{
			if(spl[i] < spl[i-1]-slope_limit)
			{
				spl[i]=spl[i-1]-slope_limit;
			}
		}

		if(i>=2 && iflag==1 && (-spl[i]+spl[i-1]) > slope_limit/2. )
        {

			spl[i]=spl[i-1]-slope_limit/2;

			iflag=0;
		}
*/

/*
		fprintf(pFile[1],"%9.1f  %8.4f  %9.2f  %7.2f  %7.3f  %7.3f  %7.2f  %7.2f   %7.2f\n",freq[i],strouhal,Lwb[i],samp[i],r,theta,DI,delta,spl[i]);
	    
		fprintf(pFile[2],"%9.4e \t %9.4e \n",freq[i],spl[i]);

		fprintf(pFile[3],"%9.4e \t %9.4e \n",freq[i],spl[i]+3.);
*/	
	}




	
}
void step8()
{
	double sum=0.;

	double ap;

	for(i=0; i<=ilast; i++)
	{
		ap=(Aref)*pow( 10.,((spl[i])/20.));

		sum+=pow(ap,2.);
	}
    sum=sqrt(sum);

	sum=20.*log10(sum/Aref);

	printf("\n\n overall SPL = %12.3g dB   ref 20 micro Pa\n\n",sum);
	fprintf(pFile[1],"\n\n overall SPL = %12.3g dB   ref 20 micro Pa\n\n",sum);
}	

void files()
{

    int iw;

	for(iw=0;iw<=5;iw++)
	{
		strcpy(filename[iw],prefix);
	}


	strcat(filename[0],"_step5.spl");

	strcat(filename[1],"_spl.out");

	strcat(filename[2],"_spl.plt");

	strcat(filename[3],"_spl_3dB.plt");

	strcat(filename[4],"_psd.plt");

	strcat(filename[5],"_psd_3db.plt");



	for(iw=0;iw<=5;iw++)
	{
		pFile[iw]=fopen(filename[iw], "w");
	}	


	fprintf(pFile[1],"\n\n %s \n",program_name);
 	fprintf(pFile[1],"\n %s \n",Irvine);
	
	
	fprintf(pFile[1],"\n\n References:  \n");
	fprintf(pFile[1],"\n 1. Rocket Vehicle Liftoff Acoustics and Skin Vibration");
	fprintf(pFile[1],"\n    Acoustic Loads Generated by the Propulsion System. ");
	fprintf(pFile[1],"\n    NASA SP-8072, Monograph N71-33195, 1971. \n");
	fprintf(pFile[1],"\n 2. Wilby & Wilby, Prediction of External Acoustic Field ");
	fprintf(pFile[1],"\n    on Taurus During Liftoff, AARC No. 121, 1991. \n");


    fprintf(pFile[1],"\n %s \n",motor_type);


	fprintf(pFile[1],"\n Thrust = %8.3g lbf", F*lbf_per_N);
	fprintf(pFile[1],"\n        = %8.3g N \n", F);
	
	fprintf(pFile[1],"\n Exhaust exit velocity = %8.4g ft/sec ", U*ft_per_m);
	fprintf(pFile[1],"\n                       = %8.4g m/sec \n", U);
	
	fprintf(pFile[1],"\n Speed of sound = %8.4g ft/sec ", cspeed*ft_per_m);
	fprintf(pFile[1],"\n                = %8.4g m/sec \n", cspeed);

	fprintf(pFile[1],"\n Acoustic Efficiency = %8.3lf decimal    ", aceff);
	fprintf(pFile[1],"\n                     = %8.3lf percent  \n", aceff*100.);

	fprintf(pFile[1],"\n Nozzle Diameter = %7.3lf in", dei*inch_per_m);
	fprintf(pFile[1],"\n                 = %7.3lf m", dei);
	
	fprintf(pFile[1],"\n\n Number of nozzles = %d ", N);
	

    if(irad==1)
	{
		fprintf(pFile[1],"\n\n  Sound radiation into hemisphere (flat ground plane). \n");
	}
	else
	{
		fprintf(pFile[1],"\n\n  Sound radiation into sphere (free space). \n");
	}


	if(iw==1)
	{
		fprintf(pFile[1],"\n\n  Apply Wilby correction factor for surface reflections.	\n");	
	}


	for(i=1;i<=7;i++)
	{
		if(igeo==i){fprintf(pFile[1],"\n\n %s \n\n",sgeo[i]);}
    }

	
	if(igeo==4 || igeo==5  || igeo==6 || igeo==7)
	{
	
		fprintf(pFile[1],"\n x1 = %10.2f inches  (distance from nozzle exit plane to ground) ",x1*inch_per_m);
		fprintf(pFile[1],"\n    = %10.2f meters  ",x1);
		fprintf(pFile[1],"\n    = %10.2f nozzle exit diameters  \n",x1/dei);
			
		fprintf(pFile[1],"\n angle from deflected flow to ground = %12.1f deg \n",rho);

    }


}
void printspl()
{
//    amp[i]=10.*log10( amp[i]*1.0e+12 );

	fprintf(pFile[0],"%12.2f %12.4g \n",freq[i],Lwb[i]);
}

void source_allocated_Eldred_Wilby()
{
    double sigma;
    double yn;

	printf("\n\n source allocation (Eldred & Wilby)\n");

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

		strouhal = freq[i]*de/U;

		if( strouhal <= 1.87)
		{
			sigma=log10(strouhal) - 0.5645;

            yn=1.61273 + 1.550865/( exp(sigma)-exp(-sigma) ); 

			samp[i] = pow(10.,yn);
		}
		else
		{
			samp[i]=0.1;
		}

	}

	for(i=0; i<= ilast; i++)
	{
		samp[i]*=de;
	}
}

void source_allocate()
{
	printf("\n\n source allocation \n");

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

		strouhal = freq[i]*de/U;

		for(j=0; j<=(num-1); j++)
		{
	
			if(strouhal==sns[j])
			{
                samp[i]=asap[i];   // asap - apparent source allocation position
              
				break;
			}
			if( strouhal>sns[j] && strouhal<sns[j+1] )
			{
				slope=log10(asap[j+1]/asap[j])/log10(sns[j+1]/sns[j]);


				samp[i] = asap[j]*pow( (strouhal/sns[j] ),slope);

				if( i> 0 && samp[i] > samp[i-1] )
				{
//	                    printf(" samp[%d]= %8.3g  asap= %8.3g str= %8.3g  sns[%d]= %8.3g  slope= %8.3g \n",i,samp[i],asap[j],strouhal,j,sns[j],slope);

				
				   printf(" samp[%d]= %8.4g  samp[%d]= %8.4g \n",i,samp[i],i-1,samp[i-1] );
				
				}


				break;
			}
		}

	}
	for(i=0; i<= ilast; i++)
	{
		samp[i]*=de;
	}
}
void directivity()
{

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

	double k1,k2;

	double m1,m2,m3,m4,m5;

	double len;

	int jk;

	c1();
	for(jk=0; jk<10; jk++)
	{
		if(theta >= an[jk] && theta <= an[jk+1] )
		{

			len=an[jk+1]-an[jk];
			k2= (theta   -an[jk])/len;

			k1=1.-k2;
			d1 = k1*dn[jk] + k2*dn[jk+1];

			break;
		}
	}


	c2();
	for(jk=0; jk<10; jk++)
	{
		if(theta >= an[jk] && theta <= an[jk+1] )
		{

			len=an[jk+1]-an[jk];
			k2= (theta   -an[jk])/len;


			k1=1.-k2;
			d2 = k1*dn[jk] + k2*dn[jk+1];

			break;
		}
	}

	c3();
	for(jk=0; jk<10; jk++)
	{
		if(theta >= an[jk] && theta <= an[jk+1] )
		{

			len=an[jk+1]-an[jk];
			k2= (theta   -an[jk])/len;


			k1=1.-k2;
			d3 = k1*dn[jk] + k2*dn[jk+1];

			break;
		}
	}


	c4();
	for(jk=0; jk<10; jk++)
	{
		if(theta >= an[jk] && theta <= an[jk+1] )
		{

			len=an[jk+1]-an[jk];
			k2= (theta   -an[jk])/len;


			k1=1.-k2;
			d4 = k1*dn[jk] + k2*dn[jk+1];

			break;
		}
	}


	c5();
	for(jk=0; jk<10; jk++)
	{
		if(theta >= an[jk] && theta <= an[jk+1] )
		{

			len=an[jk+1]-an[jk];
			k2= (theta   -an[jk])/len;


			k1=1.-k2;
			d5 = k1*dn[jk] + k2*dn[jk+1];

			break;
		}
	}




	m1= 0.004;
	m2= 0.0125;
	m3= 0.04;
	m4= 0.125;
	m5= 0.4;


//	fprintf(pFile[1],"\n %4.1f %4.1f %4.1f %4.1f %4.1f\n",d1,d2,d3,d4,d5);


    if(strouhal <= m1)
	{
		DI =d1;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}
    if(strouhal > m1  && strouhal <= m2)
	{
		len=m2-m1;
		k2=(strouhal-m1)/len;

		k1=1.-k2;
		DI = k1*d1 + k2*d2;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}
    if(strouhal > m2 && strouhal <= m3)
	{
		len=m3-m2;
		k2=(strouhal-m2)/len;

		k1=1.-k2;
		DI = k1*d2 + k2*d3;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}
    if(strouhal > m3   && strouhal <= m4)
	{
		len=m4-m3;
		k2=(strouhal-m3)/len;

		k1=1.-k2;
		DI = k1*d3 + k2*d4;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}
    if(strouhal > m4  && strouhal <= m5)
	{
		len=m5-m4;
		k2=(strouhal-m4)/len;

		k1=1.-k2;
		DI = k1*d4 + k2*d5;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}
	if(strouhal > m5)
	{
		DI =d5;

//		fprintf(pFile[1]," %12.3e %12.3e %12.3e\n",strouhal,theta,DI);
	}

//    printf(" %12.3e %12.3e %12.3e\n",strouhal,theta,DI);

//	printf("\n end directivity \n");

}

void c1()
{


	an[0]=  20.;	
	an[1]=	28.;	
	an[2]=	32.;	
	an[3]=	40.;	
	an[4]=	60.;
	an[5]=	80.;	
	an[6]=	100.;	
	an[7]=	140.;	
	an[8]=	150.;	
	an[9]=	160.;	
	an[10]=  180.; 

	dn[0]=	3.5;
	dn[1]=	7.5;
	dn[2]=	7.5;
	dn[3]=	6.;
	dn[4]=	0.65;
	dn[5]=	-4.5;
	dn[6]=	-9.;
	dn[7]=	-14.;
	dn[8]=	-15.;
	dn[9]=	-16.;
	dn[10]=	-17.2;

}
void c2()
{

	an[0]=  20.;	
	an[1]=	30.;	
	an[2]=	32.;	
	an[3]=	38.;	
	an[4]=	40.;
	an[5]=	43.;	
	an[6]=	60.;	
	an[7]=	100.;	
	an[8]=	140.;	
	an[9]=	160.;	
	an[10]=  180.; 

	dn[0]=	0.;
	dn[1]=	3.5;
	dn[2]=	6.0;
	dn[3]=	6.3;
	dn[4]=	6.0;
	dn[5]=	2.5;
	dn[6]=	-7.5;
	dn[7]=	-13.5;
	dn[8]=	-15.5;
	dn[9]=	-16.;
	dn[10]=	-17.;

}
void c3()
{

	an[0]=  20.;	
	an[1]=	30.;	
	an[2]=	40.;	
	an[3]=	48.;	
	an[4]=	52.;
	an[5]=	60.;	
	an[6]=	100.;	
	an[7]=	120.;	
	an[8]=	140.;	
	an[9]=	160.;	
	an[10]=  180.; 

	dn[0]=	-1.;
	dn[1]=	2.1;
	dn[2]=	5.0;
	dn[3]=	5.5;
	dn[4]=	5.5;
	dn[5]=	3.5;
	dn[6]=	-7.;
	dn[7]=	-9.9;
	dn[8]=	-12.5;
	dn[9]=	-14.5;
	dn[10]=	-16.;
}
void c4()
{
	an[0]=  20.;	
	an[1]=	30.;	
	an[2]=	40.;	
	an[3]=	50.;	
	an[4]=	60.;
	an[5]=	100.;	
	an[6]=	120.;	
	an[7]=	140.;	
	an[8]=	160.;	
	an[9]=	160.;	
	an[10]=  180.; 

	dn[0]=	-2.;
	dn[1]=	1.5;
	dn[2]=	3.5;
	dn[3]=	5.0;
	dn[4]=	4.0;
	dn[5]=	-5.5;
	dn[6]=	-9.;
	dn[7]=	-11.5;
	dn[8]=	-12.8;
	dn[9]=	-13.4;
	dn[10]=	-14.;
}
void c5()
{
	an[0]=  20.;	
	an[1]=	30.;	
	an[2]=	40.;	
	an[3]=	57.;	
	an[4]=	60.;
	an[5]=	80.;	
	an[6]=	100.;	
	an[7]=	120.;	
	an[8]=	140.;	
	an[9]=	160.;	
	an[10]=  180.; 

	dn[0]=	-2.5;
	dn[1]=	0.9;
	dn[2]=	2.8;
	dn[3]=	4.;
	dn[4]=	4.;
	dn[5]=	1.;
	dn[6]=	-4.;
	dn[7]=	-8.;
	dn[8]=	-10.;
	dn[9]=	-12.;
	dn[10]=	-13.;
}
void blocked()    
{

// Wilby correction factor for pressure reflections at the surface of the vehicle.
// Wilby document, equation 15, page 30.

	double sinB;

// strouhal, stationdiam, cspeed,beta

	sinB=sin(beta_rad);

    ZZ=sinB*PI*freq[i]*stationdiam/cspeed;

//	printf("\n ZZ=%12.4e  sinB=%12.4e  beta=%12.4e  freq[i]=%12.4e  stationdiam=%12.4e  cspeed=%12.4e \n",ZZ,sinB,beta,freq[i],PI,stationdiam,cspeed);
//   exit(1);

	if(ZZ <= 0.1){delta=3.;}

	if(ZZ > 0.1 && ZZ < 12.8)
	{

		ZZ=2.*log10(1.25*ZZ)/log10(10.);
		
		delta=1.5*(3.+tanh(ZZ));
	}

	if(ZZ >= 12.8){delta=6.;}

}

void material(void)
{

    cref=200000./12.;


	printf("\n Enter material: \n 1=aluminum  2=steel  3=other \n");
    scanf("%d",&imat);

    if(imat !=1 && imat !=2 && imat != 3)
	{
		printf("\n Input Error \n");

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

	if(imat==1)
	{
		W=0.1;

		cmat=cref;
	}
	if(imat==2)
	{
		W=0.29;

		cmat=cref;
	}
    if(imat==3)
	{

		if(isys==1)
		{

			printf("\n Enter cylindrical skin weight density (lbm/in^3) \n");
   			char sW[16];
			scanf("%s",&sW);
			sscanf(sW,"%lf",&W);

			printf("\n Enter elastic modulus (lbf/in^2) \n");
   			char sE[16];
			scanf("%s",&sE);
			sscanf(sE,"%lf",&E);

		    fprintf(pFile[1],"\n\n Elastic Modulus = %10.3g lbf/in^2 \n",E);

			cmat = sqrt( (E/W)*(32.2/12.) );

		}
		else
		{

			printf("\n Enter cylindrical skin weight density (kg/m^3) \n");
   			char sW[16];
			scanf("%s",&sW);
			sscanf(sW,"%lf",&W);

			printf("\n Enter elastic modulus (Pa) \n");
   			char sE[16];
			scanf("%s",&sE);
			sscanf(sE,"%lf",&E);

		    fprintf(pFile[1],"\n\n Elastic Modulus = %10.3g Pa \n",E);

			cmat = sqrt( (E/W) );

			cmat*=ft_per_m;

	        W/=kgpm3_per_lbmpin3;
		}

	}


	fprintf(pFile[1],"\n\n Speed of Sound in Material = %10.4g ft/sec ",cmat);	
	  fprintf(pFile[1],"\n                            = %10.4g m/sec \n",cmat*m_per_ft);
	
	printf("\n\n Speed of Sound in Material = %10.4g ft/sec ",cmat);
	  printf("\n                            = %10.4g m/sec \n\n",cmat*m_per_ft);

	scale = cmat/cref;

}
void franken(void)
{
	

	if(isys==1)
	{

		printf("\n Enter skin thickness (in) \n");
   		char sthick[16];
		scanf("%s",&sthick);
		sscanf(sthick,"%lf",&thickness);

    }
	else
	{

		printf("\n Enter skin thickness (mm) \n");
   		char sthick[16];
		scanf("%s",&sthick);
		sscanf(sthick,"%lf",&thickness);


		thickness*=inch_per_mm;

	}

	
	fprintf(pFile[1],"\n Skin Thickness = %12.3f inches ",thickness);
	fprintf(pFile[1],"\n                = %12.3f mm \n",thickness*mm_per_inch);

	fprintf(pFile[1],"\n\n Weight Density = %10.4g lbm/in^3 ",W);
	  fprintf(pFile[1],"\n                = %10.4g kg/m^3 \n",W*kgpm3_per_lbmpin3);


	W*=pow(12.,3.);
 
    thickness/=12.;

    W*=thickness;

	fprintf(pFile[1],"\n Weight Density (per area) = %12.3f lbm/ft^2 ",W);
	fprintf(pFile[1],"\n                           = %12.3f kg/m^2 \n",W*kg_per_lbm*pow(ft_per_m,2));


	fprintf(pFile[1],"\n\n  f(Hz) \t PSD(G^2/Hz)  \n");

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

	fr[0]=1.0;
	fr[1]=2.0;
	fr[2]=2.7;
	fr[3]=3.0;
	fr[4]=3.3;
	fr[5]=3.6;
	fr[6]=3.7;
	fr[7]=3.78;
	fr[8]=3.9;
	fr[9]=4.0;
	fr[10]=4.1;
	fr[11]=4.3;
	fr[12]=4.7;

	ar[0]=-158;	
	ar[1]=-146;
	ar[2]=-140.5;
	ar[3]=-136.;
	ar[4]=-130.;
	ar[5]=-121.;
	ar[6]=-119.;
	ar[7]=-118.8;
	ar[8]=-119.;
	ar[9]=-120.;
	ar[10]=-121.7;
	ar[11]=-123.;
	ar[12]=-124.5;


	for(i=0; i<=12; i++)
	{
		fr[i]= log10( (pow(10.,fr[i])*scale) );
	}



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

	int jlast = 12;


	for(i=0; i<= ilast; i++)
	{
		sfr = log10( f[i]*stationdiam );


//		printf("\n %12.4e  %12.4e\n",sfr,f[i]);


        if( sfr <= fr[0])
		{
			va= ar[0];
		}
		if(sfr >= fr[jlast] )
		{
			va= ar[jlast];
		}

		if( sfr > fr[0] && sfr < fr[jlast] )
		{
			
			for(j=0; j<=11; j++)
			{
				if(sfr >= fr[j] && sfr <= fr[j+1] )
				{
					len = fr[j+1] -fr[j];
					k2  = (sfr    -fr[j])/len;

					k1=1.-k2;
					va = k1*ar[j] + k2*ar[j+1];

					break;
				}
			}

		}
	

		vdb = va - 20.*log10(W) + spl[i];

		al = 1.*pow(10.,(vdb/20.));

		df = (pow(2.,(1./6.))-1./pow(2.,(1./6.)))*f[i];

//		psd[i] = 0.5*pow(al,2.)/df;

// assume vdb is in terms of RMS

		psd[i] = pow(al,2.)/df;

/*
		if(f[i] > 98. && f[i] < 102.)
		{
			printf(" f=%8.3lf  psd=%9.3e  df=%9.3lf  va=%9.3lf  20logW=%8.4lf spl=%8.3lf\n",f[i],psd[i],df,va,20.*log10(W),spl[i]);
		}
*/

/*
			fprintf(pFile[1]," %8.3lf \t %9.3e\n",f[i],psd[i]);
			fprintf(pFile[4]," %12.4e \t %12.4e\n",f[i],psd[i]);
			fprintf(pFile[5]," %12.4e \t %12.4e\n",f[i],2.*psd[i]);
*/
		}

	    double grms,grmsb,ra,rb;

		ra=0.;
		rb=0.;

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

			  s[i]=log( psd[i+1]/psd[i] )/log( f[i+1]/f[i] );

              if(s[i] < -1.0001 ||  s[i] > -0.9999 )
			  {	   
                 rb+= ( psd[i+1] * f[i+1]- psd[i]*f[i])/( s[i]+1.);
			  }
		      else
			  {
                 rb+= psd[i]*f[i]*log( f[i+1]/f[i]);
			  }
		 
              if(f[i+1] < 2050.){ra=rb;}

		}
        grms=sqrt(ra);
		grmsb=sqrt(rb);

		printf("\n Ring Frequency = %12.4g Hz \n",cmat/(PI*stationdiam));
		fprintf(pFile[1],"\n Ring Frequency = %12.4g Hz \n",cmat/(PI*stationdiam));


		fprintf(pFile[1],"\n\n Overall Skin Vibration Level (up to 2000 Hz) = %12.4g GRMS\n",grms );
		printf("\n\n Overall Skin Vibration Level (up to 2000 Hz)= %12.4g GRMS\n",grms );

		fprintf(pFile[1],"\n Overall Skin Vibration Level (up to 10,000 Hz) = %12.4g GRMS\n",grmsb );
		printf("\n Overall Skin Vibration Level (up to 10,000 Hz)= %12.4g GRMS\n",grmsb );

}
void length_diameter()
{
	
	printf("\n\n Note:  station length is referenced to nozzle exit plane \n\n");


	if(isys==1)
	{
			printf("\n Enter station length (inches) \n");
	}
	else
	{
			printf("\n Enter station length (meters) \n",xs);
	}

	    char slen[12];
        scanf("%s",&slen);
        sscanf(slen,"%lf",&xs);

		if(isys==1)
		{
			xs*=m_per_ft/12.; 
		}

		if(isys==1)
		{
			printf("\n Enter station diameter (in)\n");
		}
		else
		{
			printf("\n Enter station diameter (m)\n");
		}
	    char sdia[12];
        scanf("%s",&sdia);
        sscanf(sdia,"%lf",&stationdiam);
	    
		if(isys==1)
		{
			stationdiam*=m_per_ft/12.; 
		}

}
void method2()
{

	long i;

	
	dx=de;

	for(i=0;i<=NUMX;i++)
	{
		xm2[i]=(i+1)*de;
	}



	Lwb_recordmax=1.0e+90;

    iflag=0;

	jjk=1;

	for(ie=1;ie<ie_max;ie++)
	{
//		iring = (rand()/max);

		if(iflag==1)
		{
			jjk++;

			if(jjk==long(ie_max/40.))
			{
				iflag=0;
				jjk=0;
			}
		}
		else
		{
			jjk=0;
		}




		xt=(53.*(rand()/max)+7.)*de;     // core length

		if(iflag==1 || ((rand()/max)<0.5 && ie >= 100))
		{
			xt=xtr*(0.975+0.05*rand()/max);
		}

		normalized_relative_sound_power_per_unit_core_length();

		normalized_relative_sound_power_spectrum_level();

		SPL_check();

		record_check();

//		printf(" %ld  %8.4g  %8.4g  %8.4g  %8.4g\n",ie,Lwb_record,Lwb_recordmax,SPL_record,SPL_recordmax);

	}

	double sum1;
	double sum2;

	printf("\n");
	fprintf(pFile[1],"\n");

	for(jk=0;jk<=LAST_F_N;jk++)
	{	
		sum1=0.;
		sum2=0.;

		for(i=0;i<=NUMX;i++)
		{
			sum1+=pow(10,(sss[i][jk]/10.));
			sum2+=pow(10,(ppp[i][jk]/10.));
		}
		sum1=10.*log10(sum1);
		sum2=10.*log10(sum2);

		printf(" %8.4g \t %8.4g \t %8.4g \t %8.4g \t %8.4g \n",fc[jk],sum1,Lwb_reference[jk],sum2,SPL_reference[jk]);

		fprintf(pFile[1]," %8.4g \t %8.4g \n",fc[jk],sum2);
		fprintf(pFile[2]," %8.4g \t %8.4g \n",fc[jk],sum2);
		fprintf(pFile[3]," %8.4g \t %8.4g \n",fc[jk],sum2+3.);
	}

	fprintf(pFile[1],"\n\n xt=%8.4g \n\n",xt);

	for(i=0;i<=ncoord_core;i++)
	{
		fprintf(pFile[1]," %8.4g \t %8.4g \n",xxr_core[i],yyr_core[i]);
	}
	
	fprintf(pFile[1],"\n\n");

	for(i=0;i<=ncoord;i++)  // elk
	{
		fprintf(pFile[1]," %8.4g \t %8.4g \n",xxr[i],yyr[i]);
	}

}
void record_check()
{

//	    printf(" ie=%ld  jjk=%d  iflag=%d  xt=%8.4g  xtr=%8.4g  %8.4g  %8.4g \n",ie,jjk,iflag,xt,xtr,Lwb_record,SPL_record); 

//		getch();

	    double phrase=Lwb_record + SPL_record;

		if( phrase < Lwb_recordmax  )
		{
			iflag=1;

			jjk=0;


			Lwb_recordmax=phrase;

			xtr=xt;


			printf(" %ld  xt=%8.4g  error=%8.4g  \n",ie,xt,Lwb_recordmax);

		
			for(jk=0;jk<=LAST_F_N;jk++)
			{	
				for(i=0;i<=NUMX;i++)
				{
					sss[i][jk]=LWSB[i][jk];
					ppp[i][jk]=SPL_sbp[i][jk];
				}
			}

//            printf("\n");
			for(i=0;i<=ncoord;i++)
			{
				xxr[i]=xx[i];
				yyr[i]=yy[i];

//				printf(" %8.4g \t %8.4g \n",xxr[i],yyr[i]);

			}
//            printf("\n");

			for(i=0;i<=ncoord_core;i++)
			{
				xxr_core[i]=xx_core[i];
				yyr_core[i]=yy_core[i];

				if(xx_core[i]<=0)
				{
					printf("\n error**  xx_core[%ld]=%8.4g \n",i,xx_core[i]);
					exit(1);
				}
			}

		}


}
void SPL_check()
{
	
//	iflag=0;

	long jk;

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

	    theta=180.;
		

		if(igeo==1 || igeo==2 || igeo==3)  // undeflected
		{
			radius[i]=xs+xm2[i];    // samp[i] is apparent source allocation position

//		fprintf(pFile[1]," r=%8.4g  xs=%8.4g  x=%8.4g \n",radius[i],xs,xm2[i]);
		
		
		}
		else   // deflected
		{

			x2= xm2[i] -x1;


			if(x2 > 0. )
			{

				ry= (xs+x1)-xm2[i]*sin(PI*rho/180.);

				rx= x2*cos(PI*rho/180.);

				radius[i]=sqrt( pow(rx,2.) + pow(ry,2.) );

				beta_rad = atan2(ry,rx);
				beta_deg = 180.*beta_rad/PI;

				theta = 180. - rho - beta_deg;   

//				fprintf(pFile[1]," r=%8.4g  rx=%8.4g  ry=%8.4g \n",radius[i],rx,ry);

			}
			else   // special case where source occurs before deflection
			{

//                  xs  =  station length
//              samp[i] =  apparent source allocation position

//				r = xs + samp[i];
//				r += 2.*fabs((x1-fabs(samp[i])));

				
				radius[i] = xs + x1;

//				fprintf(pFile[1]," r=%8.4g  xs=%8.4g  x1=%8.4g \n",radius[i],xs,x1);


                theta = 90.-rho;   
			}

		
		}

		if(theta>180.){theta=180.;}
		if(theta< 20.){theta= 20.;}


//		fprintf(pFile[1],"\n");

		for(jk=0;jk<=LAST_F_N;jk++)
		{

			strouhal = fcc[jk]*de/U;
			directivity();

			double ab = 8.; // hemisphere (flat ground plane)

			if(irad==2)     // sphere (free space)
			{
				ab = 11;
			}

			delta=0.;

			if(iw==1)
			{
				blocked();
			}


			term1=LWSB[i][jk];
			term2=-10.*log10(pow(radius[i],2.));
			term3=-11;
			term4= DI;

			SPL_sbp[i][jk] = term1 + term2 + term3 + term4 + ab + delta;

//			fprintf(pFile[1]," i=%ld fcc=%8.4g  r=%8.4g LWSB=%8.4g t2=%8.4g t3=%8.4g DI=%8.4g %8.4g\n",i,fcc[jk],radius[i],term1,term2,term3,term4,SPL_sbp[i][jk]);

		}

	}

	double sum=0.;

	SPL_record=0.;

	for(jk=0;jk<=LAST_F_N;jk++)
	{		
		for(i=0;i<=NUMX;i++)
		{
			sum+=pow(10,(SPL_sbp[i][jk]/10.));
		}
	
		sum=10.*log10(sum);

		if(fabs(SPL_reference[jk]-sum)>SPL_record)
		{
			SPL_record=fabs(SPL_reference[jk]-sum);
		}
	}
//	sum=10.*log10(sum);

//	SPL_record=10.*log10(SPL_record);


//	printf("\n\n overall SPL = %8.4g dB\n",sum);

}
void normalized_relative_sound_power_per_unit_core_length()
{


//	fprintf(pFile[1],"\n\n  normalized acoustic power per unit core length  \n\n");

//  obtain the normalized acoustic power per unit core length
//  from figure 12


	double xxa[NUMX1];
	double apos,aposa;
	double xa;
	double aaa;
	double length;
	double c1,c2;

	double temp;

	long ia,ib;

	ncoord_core=7;

	if(ncoord_core<2)
	{
		printf("\n ncoord_core error \n");
		exit(1);
	}


//	xx_core[0]=  0.1*(0.8+0.4*rand()/max);         
//	yy_core[0]=-18.7*(0.8+0.4*rand()/max);


		
		xx_core[0]=0.05;
		xx_core[ncoord_core]=10.;

        xx_core[1]=  0.2*rand()/max+xx_core[0];
        xx_core[2]=  1.0*rand()/max+xx_core[0];
        xx_core[3]=  2.0*rand()/max+xx_core[0];
        xx_core[4]=  5.0*rand()/max+xx_core[0];
        xx_core[5]=  8.0*rand()/max+xx_core[0];
        xx_core[6]= 10.0*rand()/max+xx_core[0];


		for(ia=0;ia<=ncoord_core-1;ia++)
		{
			for(ib=ia+1;ib<=ncoord_core;ib++)
			{
				if(xx_core[ia]>xx_core[ib])
				{
					temp=xx_core[ia];
					xx_core[ia]=xx_core[ib];
					xx_core[ib]=temp;
				}
			}
		}

		if( rand()/max < 0.5 )
		{

			yy_core[4]=-3.;

			yy_core[3]=yy_core[4]-10.*(0.8*rand()/max+0.2);
			yy_core[2]=yy_core[3]-10.*(0.8*rand()/max+0.2);
			yy_core[1]=yy_core[2]-10.*(0.8*rand()/max+0.2);
			yy_core[0]=yy_core[1]-10.*(0.8*rand()/max+0.2);

			yy_core[5]=yy_core[4]-10.*(0.8*rand()/max+0.2);
			yy_core[6]=yy_core[5]-10.*(0.8*rand()/max+0.2);
			yy_core[7]=yy_core[6]-10.*(0.8*rand()/max+0.2);

		}
		else
		{
			yy_core[3]=-3.;


			yy_core[2]=yy_core[3]-10.*(0.8*rand()/max+0.2);
			yy_core[1]=yy_core[2]-10.*(0.8*rand()/max+0.2);
			yy_core[0]=yy_core[1]-10.*(0.8*rand()/max+0.2);

			yy_core[4]=yy_core[3]-10.*(0.8*rand()/max+0.2);
			yy_core[5]=yy_core[4]-10.*(0.8*rand()/max+0.2);
			yy_core[6]=yy_core[5]-10.*(0.8*rand()/max+0.2);
			yy_core[7]=yy_core[6]-10.*(0.8*rand()/max+0.2);
		}

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

	for(i=0;i<=ncoord_core;i++)
	{
        if(xx_core[i]<=0)
		{
			printf(" error1a:  xx_core[%ld]=%8.4g \n",i,xx_core[i]);
			exit(1);
		}
        if(yy_core[i]>=0)
		{
			printf(" error2a:  yy_core[%ld]=%8.4g \n",i,yy_core[i]);
			exit(1);
		}
	}
//*********************************************


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

	if(iflag==1 || ((rand()/max)<0.5 && ie >= 100))
	{

		for(i=0;i<=ncoord_core;i++)
		{
			yy_core[i]=yyr_core[i]*(0.975+0.05*(rand()/max)); 
			xx_core[i]=xxr_core[i]*(0.975+0.05*(rand()/max)); 
			

        if(xx_core[i]<=0)
		{
			printf(" error1aa:  xx_core[%ld]=%8.4g   xxr_core[%ld]=%8.4g\n",i,xx_core[i],i,xxr_core[i]);
			exit(1);
		}



		}
	}

	for(i=0;i<=ncoord_core;i++)
	{
        if(xx_core[i]<=0)
		{
			printf(" error1b:  xx_core[%ld]=%8.4g \n",i,xx_core[i]);
			exit(1);
		}
        if(yy_core[i]>=0)
		{
			printf(" error2b:  yy_core[%ld]=%8.4g \n",i,yy_core[i]);
			exit(1);
		}
	}




	for(ia=1;ia<=ncoord_core-1;ia++)
	{
		for(ib=ia+1;ib<=ncoord_core;ib++)
		{
			if(xx_core[ia]>xx_core[ib])
			{
					temp=xx_core[ia];
					xx_core[ia]=xx_core[ib];
					xx_core[ib]=temp;

					temp=yy_core[ia];
					yy_core[ia]=yy_core[ib];
					yy_core[ib]=temp;
			}
		}
	}

	for(i=0;i<=ncoord_core;i++)
	{
        if(xx_core[i]<=0)
		{
			printf(" error1c:  xx_core[%ld]=%8.4g \n",i,xx_core[i]);
			exit(1);
		}
        if(yy_core[i]>=0)
		{
			printf(" error2c:  yy_core[%ld]=%8.4g \n",i,yy_core[i]);
			exit(1);
		}
	}


	for(i=0;i<=ncoord_core;i++)
	{
        if(xx_core[i]<=0)
		{
			printf(" error1d:  xx_core[%ld]=%8.4g \n",i,xx_core[i]);
			exit(1);
		}
        if(yy_core[i]>=0)
		{
			printf(" error2d:  yy_core[%ld]=%8.4g \n",i,yy_core[i]);
			exit(1);
		}

		xxa[i]=pow(10,xx_core[i]);
	}


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

		normalized_aps_per_length[i]=-18.703;
		aaa=normalized_aps_per_length[i]/10.;
		aaa=pow(10.,aaa);
		Wp[i]=aaa*WOA/xt;

		apos=xm2[i]/xt;

		aposa=pow(10,apos);

		if(apos<xxa[0]){apos=xxa[0];}
		if(apos>xxa[ncoord_core]){apos=xxa[ncoord_core];}

// interpolate

		for(j=0;j<=(ncoord_core-1);j++)
		{
			if(aposa>=xxa[j] && aposa<=xxa[j+1])
			{
				xa=aposa-xxa[j];
				length=xxa[j+1]-xxa[j];
				c2=xa/length;
				c1=1.-c2;

				normalized_aps_per_length[i]=c1*yy_core[j] + c2*yy_core[j+1];

				aaa=normalized_aps_per_length[i]/10.;

				aaa=pow(10.,aaa);

				Wp[i]=aaa*WOA/xt;

				break;
			}

		}

	}


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

	for(i=0;i<=NUMX;i++)
	{
		term1=  normalized_aps_per_length[i];
		term2=  LW;
		term3=  10.*log10(dx/xt);

		LWS[i]=term1+term2+term3;       // re  10^(-12) W         
	}

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


/*
	          printf("\n\n W & normalized aps per length \n");
	fprintf(pFile[1],"\n\n W & normalized aps per length \n");


	for(i=0;i<=NUMX;i++) 
	{
		          printf(" %ld \t %8.4g \t %8.4g \t %8.4g \n",i,xm2[i]/xt,Wp[i],normalized_aps_per_length[i]);
		fprintf(pFile[1]," %ld \t %8.4g \t %8.4g \t %8.4g \n",i,xm2[i]/xt,Wp[i],normalized_aps_per_length[i]);
	}

	          printf("\n");
	fprintf(pFile[1],"\n");
*/
}
void normalized_relative_sound_power_spectrum_level()
{

	double xxa[NUMX1];
	double apos,aposa;
	double xa;

	double length;
	double c1,c2;

	double delta_fb;

	double a_sub_o = cspeed;
	double a_sub_e = 1066;   // m/sec  (3500 ft/sec)   

	double oct=pow(2,(1./6.));

	long jk;


//	double temp;

	double ratio = a_sub_e/(U*a_sub_o);

//	printf("\n ratio = %8.4g \n",ratio);

	ncoord=7;


//	iflag=0;

//	for(ie=1;ie<ie_max;ie++)   // moose
//	{

//		****go back to old way;

/*
		xx[0]=0.05*(0.95*rand()/max+0.1);   
		xx[1]=0.5*(0.95*rand()/max+0.1);
		xx[2]=2.0*(0.95*rand()/max+0.1);
		xx[3]=10.*(0.95*rand()/max+0.1);
		xx[4]=50.*(0.95*rand()/max+0.1);
		xx[5]=100.*(0.95*rand()/max+0.1);
		xx[6]=200.*(0.95*rand()/max+0.1);
		xx[ncoord]=500.;

		yy[0]=50*rand()/max;
		yy[1]=50*rand()/max;
		yy[2]=50*rand()/max;
		yy[3]=50*rand()/max;
		yy[4]=50*rand()/max;
		yy[5]=50*rand()/max;
		yy[6]=50*rand()/max;
		yy[ncoord]=50*rand()/max;
*/

		
		xx[0]=0.05;
		xx[ncoord]=500.;

        xx[1]=  1.*rand()/max;
        xx[2]= 10.*rand()/max;
        xx[3]= 30.*rand()/max;
        xx[4]= 60.*rand()/max;
        xx[5]=150.*rand()/max;
        xx[6]=250.*rand()/max;

		int ia,ib;

		double temp;

		for(ia=0;ia<=ncoord-1;ia++)
		{
			for(ib=ia+1;ib<=ncoord;ib++)
			{
				if(xx[ia]>xx[ib])
				{
					temp=xx[ia];
					xx[ia]=xx[ib];
					xx[ib]=temp;
				}
			}
		}


		if( rand()/max < 0.5)
		{

			yy[3]=1.;

			yy[2]=yy[3]-10.*(0.8*rand()/max+0.2);
			yy[1]=yy[2]-10.*(0.8*rand()/max+0.2);
			yy[0]=yy[1]-10.*(0.8*rand()/max+0.2);

			yy[4]=yy[3]-10.*(0.8*rand()/max+0.2);
			yy[5]=yy[4]-10.*(0.8*rand()/max+0.2);
			yy[6]=yy[5]-10.*(0.8*rand()/max+0.2);
			yy[7]=yy[6]-10.*(0.8*rand()/max+0.2);

		}
		else
		{


			yy[4]=1.;

			yy[3]=yy[4]-10.*(0.8*rand()/max+0.2);
			yy[2]=yy[3]-10.*(0.8*rand()/max+0.2);
			yy[1]=yy[2]-10.*(0.8*rand()/max+0.2);
			yy[0]=yy[1]-10.*(0.8*rand()/max+0.2);

			yy[5]=yy[4]-10.*(0.8*rand()/max+0.2);
			yy[6]=yy[5]-10.*(0.8*rand()/max+0.2);
			yy[7]=yy[6]-10.*(0.8*rand()/max+0.2);
		}


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


//		for(ia=1;ia<=ncoord;ia++)
//		{
//			printf(" %8.4g\n",xx[ia]);
//		}
//		exit(1);

		
	    if(iflag==1 || ((rand()/max)<0.5 && ie >= 100))
		{
			int jj;

			for(jj=1;jj<ncoord;jj++)
			{
				xx[jj]=xxr[jj]*(0.975+0.05*rand()/max);
				yy[jj]=yyr[jj]*(0.975+0.05*rand()/max);
			}
			yy[ncoord]=yyr[ncoord]*(0.975+0.05*rand()/max);
		}
   

		for(i=0;i<=ncoord;i++)
		{
			xxa[i]=pow(10,xx[i]);
		}


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

			for(jk=0;jk<NFREQ;jk++)
			{

				apos=(fc[jk]*xm2[i]*a_sub_e)/(U*a_sub_o);

				fcc[jk]=fc[jk];

				apos=(fcc[jk]*xm2[i])*ratio;

				aposa=pow(10,apos);

// interpolate

				for(j=0;j<=(ncoord-1);j++)
				{
					if(aposa>=xxa[j] && aposa<=xxa[j+1])
					{
						xa=aposa-xxa[j];
						length=xxa[j+1]-xxa[j];
						c2=xa/length;
						c1=1.-c2;

						power_spectrum[i][jk]=c1*yy[j] + c2*yy[j+1];
						


						break;
					}
				}

				if(aposa<=xxa[0])
				{
					power_spectrum[i][jk]=yy[0];
					

				}
				if(aposa>=xxa[ncoord])
				{
					power_spectrum[i][jk]=yy[ncoord]; 

				}

//			fprintf(pFile[1],"i=%ld jk=%ld  fc=%8.4g  apos=%8.4g  ps=%8.4g \n",i,jk,fcc,apos,power_spectrum[i][jk]);

			}
		}

//
//        U = fully expanded exit velox   
//  a_sub_o = speed of sound in atmosphere
//  a_sub_e = speed of sound at nozzle exit
//
//  assume a_sub_e = U


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

			for(jk=0;jk<=LAST_F_N;jk++)
			{
				delta_fb=( -(1./oct) +oct )*fc[jk];

//				delta_fb=1;

				if(delta_fb <=0)
				{

					printf("\n      oct = %8.4g   \n",oct);
					printf("\n  fc[%ld] = %8.4g   \n",jk,fc[jk]);
					printf("\n delta_fb = %8.4g   \n",delta_fb);

					exit(1);
				}

				term1=power_spectrum[i][jk];
				
				term2=LWS[i];


				if(term2<=0.)
				{
					printf("warning: LWS[%ld]=%8.4g \n",i,LWS[i]);
				}


				if(xm2[i]<=0.)
				{
					printf("error: xm2[%ld]=%8.4g \n",i,xm2[i]);
					exit(1);
				}


				term3=-10.*log10((U*a_sub_o)/(xm2[i]*a_sub_e));
				
				term4=10.*log10(delta_fb);   

//			The coefficient should probably be 10.  NASA SP-8072 appears to have an error.
//			term4 is zero if delta_fb=1, however, thus making the coefficient a moot point.

				LWSB[i][jk]=term1+term2+term3+term4;

				if(LWSB[i][jk]<=0.)
				{
	

				printf("error1: LWSB[%ld][%ld]=%8.4g \n",i,jk,LWSB[i][jk]);

				printf(" i=%ld jk=%ld t1=%8.4g t2=%8.4g t3=%8.4g t4=%8.4g LWSB=%8.4g \n",i,jk,term1,term2,term3,term4,LWSB[i][jk]);

					exit(1);
				}
//
//			fprintf(pFile[1]," i=%ld jk=%ld t1=%8.4g t2=%8.4g t3=%8.4g t4=%8.4g LWSB=%8.4g \n",i,jk,term1,term2,term3,term4,LWSB[i][jk]);

//          getch();

			}
		}
//**************** scale

		sum=0.;

		for(jk=0;jk<=LAST_F_N;jk++)
		{	

			for(i=0;i<=NUMX;i++)
			{
				sum+=pow(10,(LWSB[i][jk])/10);
				
                
				if(LWSB[i][jk]<=0.)
				{
					printf("error2: LWSB[%ld][%ld]=%8.4g \n",i,jk,LWSB[i][jk]);
					exit(1);
				}

			}

			if(sum<=0){sum=1.0e+90;}
		}
		sum=10.*log10(sum);

		if(sum>=-1000. && sum<=1000.)
		{

		}
		else
		{
			sum=1.0e+90;
		}


		double scale = LW-sum;




		for(jk=0;jk<=LAST_F_N;jk++)
		{	
			for(i=0;i<=NUMX;i++)
			{
				LWSB[i][jk]+=scale;
				
				if(LWSB[i][jk]<=0)
				{
		            printf(" LW=%8.4g  sum=%8.4g  scale=%8.4g \n",LW,sum,scale);

					printf("warning3: LWSB[%ld][%ld]=%8.4g  scale=%8.4g \n",i,jk,LWSB[i][jk],scale);
					
					LWSB[i][jk]=1.0e-90;
				}
			}
		}

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

		Lwb_record=0;

		for(jk=0;jk<=LAST_F_N;jk++)
		{	
			sum=0.;

			for(i=0;i<=NUMX;i++)
			{
				sum+=pow(10,(LWSB[i][jk])/10);	
			}
			sum=10.*log10(sum);


		if(fabs(Lwb_reference[jk]-sum)>Lwb_record)
		{
			Lwb_record=fabs(Lwb_reference[jk]-sum);
		}

//            printf(" fc[%ld]=%8.4g  Lwb=%8.4g   Lwbr=%8.4g \n",jk,fc[jk],sum,Lwb_reference[jk]);

		}

//        getch(); 	


//	}

//**

/*
              printf("\n acoustic power for each frequency \n");
    fprintf(pFile[1],"\n acoustic power for each frequency \n");

	double maxx;
	double maxp;

	double ps;

	for(jk=0;jk<LAST_F_N;jk++)
	{
		sum=0.;

		maxx=0.;
        maxp=-1.0e+90;
		
		for(i=0;i<=NUMX;i++)
		{
            if(LWSB[i][jk]>maxp)
			{
				maxp=LWSB[i][jk];

				maxx=xm2[i];

				ps=power_spectrum[i][jk];
			}

			sum+=pow(10,(LWSB[i][jk])/10);	
		}
		ssum[jk]=10.*log10(sum);

		apos=(fcc[jk]*maxx)*ratio;

	     printf(" fcc=%8.4g \t maxx=%8.4g \t maxp=%8.4g \t apos=%8.4g \t LWSB=%8.4g \n",fcc[jk],maxx,maxp,apos,ssum[jk]);
		fprintf(pFile[1]," fcc=%8.4g \t maxx=%8.4g \t maxp=%8.4g \t apos=%8.4g \t LWSB=%8.4g \n",fcc[jk],maxx,maxp,apos,ssum[jk]);

	}

*/
}