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

#define MAX 10000

#define NUMX 200

#define NUMX1 201

// #define NFREQ 11000

#define NFREQ 30


void input_data(void);

void step1(void);
void step2(void);
void step3(void);
void step4(void);
void step5(void);
void step6(void);
void step7(void);
void step8(void);
void step9(void);
void step10(void);
void step11(void);
void step12(void);



void files(void);

void onethird(void);

void directivity(void);

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

void blocked(void);

void length_diameter(void);

void step5_reference(void);

int ncoord;

double ssum[NFREQ];

double an[100],dn[100];

double beta_rad;

double c,de,Me;
double xt;
double xs;

double xfinal;
double dx;
double x[NUMX1];

double rx,ry;

double r;

double x2;

double beta_deg;

double theta;

double delta;

double strouhal;

double sum;

double DI;

double ZZ;

double normalized_aps_per_length[NUMX1];

double power_spectrum[NUMX1][NFREQ];

double LWSB[NUMX1][NFREQ];

double sss[NUMX1][NFREQ];

double W[NUMX1];

double LWS[NUMX1];

double fc[100];

double radius[NUMX1];

double SPL_sbp[NUMX1][NFREQ];

double SPL_bp[NFREQ];

double SPL[NFREQ];

double Lwb_ref[NFREQ];

double xxr[100],yyr[100];

long i,j,jk;

long ilast;

long LAST_F_N;

double WOA;

int motor;

int iw;

int i_end_slope;

long ie_max=10000;

double fcc[NFREQ];

double db_per_octave;
double slope_limit;

double term1,term2,term3,term4;


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 aceff;
int N;
double F;
double U;

double MLRS_vel;
double MLRS_thrust;

double power_scale;

double LW;

double d1,d2,d3,d4,d5;

int isys;

int irad;
int igeo;

double cspeed;

double dei;

double stationdiam;

double x1;
double rho;

const double max=32767.;

char motor_type[40];

char sgeo[10][70];

char prefix[25];

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




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

void main()
{


	strcpy(program_name," liftoff_h_undeflected.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 Select system");
	printf("\n 1=English  2=metric \n");
	scanf("%d",&isys);

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


	input_data();

	long ikv=1;

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

	printf("\n\n Enter output filename prefix \n");

	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);

	onethird();

	fprintf(pFile[1],"\n\n step 1 \n");
	step1();
	fprintf(pFile[1],"\n\n step 2 \n");
	step2();
	fprintf(pFile[1],"\n\n step 3 \n");
	step3();
	fprintf(pFile[1],"\n\n step 4 \n");
	step4();
	fprintf(pFile[1],"\n\n step 5 \n");
	step5();
	fprintf(pFile[1],"\n\n step 6 \n");
	step6();
	fprintf(pFile[1],"\n\n step 7 \n");
	step7();
	fprintf(pFile[1],"\n\n step 8 \n");
	step8();
	fprintf(pFile[1],"\n\n step 9 \n");
	step9();
	fprintf(pFile[1],"\n\n step 10 \n");
	step10();
	fprintf(pFile[1],"\n\n step 11 \n");
	step11();

}
void step1(void)
{
// step 1 - determine the flow axis relative to the vehicle & stand

}
void step2(void)
{

// step 2 - estimate the overall acoustic power

	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);
    fprintf(pFile[1],"\n\n overall acoustic power WOA = %12.3g Watts\n",WOA); 	

}
void step3(void)
{

// step 3 - calculate the overall sound power level
	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(void)
{

// step 4  - computer an exit nozzle diamter if the vehicle has more than one nozzle
	
	de=sqrt(N)*dei;

	fprintf(pFile[1],"\n de = %8.4g m    ",de);
	fprintf(pFile[1],"\n    = %8.4g ft \n",de*ft_per_m);
}
void step5(void)
{
// step 5  - determine the length of the supersonic core

// de = nozzle exit diameter
// Me = fully expanded nozzle exit mach number


	Me=U/cspeed;

	xt= 3.45*pow((1.+0.38*Me),2)*de;


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


	          printf("\n supersonic core length  xt=%8.4g m    ",xt);
	          printf("\n                         xt=%8.4g ft \n",xt*ft_per_m);

	fprintf(pFile[1],"\n supersonic core length  xt=%8.4g m    ",xt);
	fprintf(pFile[1],"\n                         xt=%8.4g ft \n",xt*ft_per_m);	


}
void step6(void)
{

// step 6  - divide the flow into a number of slices
	
	xfinal = 4*xt;

	dx=xfinal/NUMX;

	fprintf(pFile[1],"\n dx = %8.4g m    ",dx);
	fprintf(pFile[1],"\n    = %8.4g ft \n",dx*ft_per_m);

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

}
void step7(void)   
{

	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 xx[NUMX1];
	double yy[NUMX1];
	double xxa[NUMX1];
	double apos,aposa;
	double xa;
	double aaa;
	double length;
	double c1,c2;


	xx[0]=0.1;         
	yy[0]=-18.703;
 
	xx[1]=0.12198;     
	yy[1]=-17.573;
 
	xx[2]=0.15953;     
	yy[2]=-15.999;
 
	xx[3]=0.2188;     
	yy[3]=-14.014;
 
	xx[4]=0.31887;     
	yy[4]=-11.891;
 
	xx[5]=0.41701;      
	yy[5]=-10.18;
 
	xx[6]=0.51776;     
	yy[6]=-8.7767;
	 
	xx[7]=0.65693;     
	yy[7]=-7.3732;
 
	xx[8]=0.81919;     
	yy[8]=-6.0383;
 
	xx[9]=1.0394;     
	yy[9]=-4.6348;
 
	xx[10]=1.2042;      
	yy[10]=-3.779;
 
	xx[11]=1.4073;     
	yy[11]=-3.0596;
 
	xx[12]=1.638;     
	yy[12]=-2.9545; 

	xx[13]=1.9233;      
	yy[13]=-3.293; 

	xx[14]=2.1722;     
	yy[14]=-4.0076; 

	xx[15]=2.38;     
	yy[15]=-4.9275; 

	xx[16]=2.6313;     
	yy[16]=-6.7688; 

	xx[17]=2.8592;     
	yy[17]=-8.7469; 

	xx[18]=3.0793;     
	yy[18]=-10.008; 

	xx[19]=3.3606;     
	yy[19]=-12.089; 

	xx[20]=3.6837;     
	yy[20]=-14.237; 

	xx[21]=3.9334;      
	yy[21]=-15.84; 

	xx[22]=4.2553;     
	yy[22]=-17.614; 

	xx[23]=4.5436;     
	yy[23]=-19.114; 

	xx[24]=4.7259;     
	yy[24]=-20.001; 

	int ncoord=24;


	for(i=0;i<=ncoord;i++)
	{
		xxa[i]=pow(10,xx[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);
		W[i]=aaa*WOA/xt;

		apos=x[i]/xt;

		aposa=pow(10,apos);

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

// 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;

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

				aaa=normalized_aps_per_length[i]/10.;

				aaa=pow(10.,aaa);

				W[i]=aaa*WOA/xt;

				break;
			}

		}

	}


	          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,x[i]/xt,W[i],normalized_aps_per_length[i]);
		fprintf(pFile[1]," %ld \t %8.4g \t %8.4g \t %8.4g \n",i,x[i]/xt,W[i],normalized_aps_per_length[i]);
	}

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

}
void step8(void)
{

// LWS  acoustic power for each slice

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


	double term1,term2,term3;

	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(" x[%ld]=%8.4g \t LWS=%8.4g \n",i,x[i],LWS[i]);
        fprintf(pFile[1]," %ld \t %8.4g \t %8.4g \t %8.4g \t %8.4g \n",i,term1,term2,term3,LWS[i]);
	}


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

	fprintf(pFile[1],"\n Total acoustic power = %8.4g dB\n",sum);


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

}
void step9(void)
{
	step5_reference();

	double record;
	double recordmax=1.0e+90;

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

	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;

//	Lwb_ref[i]

	long ie;

//	double temp;

	double ratio = a_sub_e/(U*a_sub_o);

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

	ncoord=5;

	long iflag=0;

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

		xx[0]=0.05;      
		yy[0]=0.; 

//		int jp;

//		double mmm;

//		double dec;

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

		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[ncoord]=500.;

		yy[1]=30*rand()/max;
		yy[2]=50*rand()/max;
		yy[3]=50*rand()/max;
		yy[4]=50*rand()/max-10.;
		yy[ncoord]=60*rand()/max-30.;
/*
		for(jp=1;jp<=ncoord;jp++)
		{
            double dec_start=-1.;
			double dec_end=2.;

			dec=dec_end-dec_start;
            
			mmm=dec_start+(dec/double(ncoord-1))*(jp-1);

			xx[jp]=5.*pow(10.,mmm);


//			printf(" jp=%d  mmm=%8.4g xx=%8.4g\n",jp,mmm,xx[jp]);
			
			yy[jp]=(rand()/max)*40.;

			if(rand()/max<0.5)
			{
				yy[jp]=-yy[jp];
			}
		}
*/

		int ia,ib;

		double temp;

		for(ia=1;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;
				}
			}
		}
//		for(ia=1;ia<=ncoord;ia++)
//		{
//			printf(" %8.4g\n",xx[ia]);
//		}
//		exit(1);

		
		if(ie>long(ie_max/10.) && iflag==1 && (rand()/max)<0.6)
		{
			int jj;

			for(jj=1;jj<ncoord;jj++)
			{
				xx[jj]=xxr[jj]*(0.99+0.02*rand()/max);
				yy[jj]=yyr[jj]*(0.99+0.02*rand()/max);
			}
			yy[ncoord]=yyr[ncoord]*(0.99+0.02*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]*x[i]*a_sub_e)/(U*a_sub_o);

				fcc[jk]=fc[jk];

				apos=(fcc[jk]*x[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];

				term3=-10.*log10((U*a_sub_o)/(x[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;
//
//			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);	
			}
		}
		sum=10.*log10(sum);

		double scale = LW-sum;


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

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

		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);

			record+=fabs(sum-Lwb_ref[jk]);

		}

		if(record < recordmax)
		{
			iflag=1;

			recordmax=record;

			printf(" %ld  %8.4g  %8.4g \n",ie,record,recordmax);

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

			for(i=0;i<=ncoord;i++)
			{
				xxr[i]=xx[i];
				yyr[i]=yy[i];
			}
		}


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

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

//**


              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=x[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]);

	}


}
void step10(void)
{

	int iflag=0;

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

	    theta=180.;
		

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

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

			x2= x[i] -x1;


			if(x2 > 0. )
			{

				ry= (xs+x1)-x[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  x=%8.4g \n",radius[i],rx,ry,x[i]);

			}
			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  x[i]=%8.4g\n",radius[i],xs,x1,x[i]);


                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;
			}

			ab-=3; // correction factor

			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.;

	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);

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

	                              
}
void step11(void)  
{

// compute the sound-pressure level in each frequency band at each point
// by logarithmic summation of contributions from each of the slices


	double sum;

	printf("\n\n Final SPL \n");

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

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

			sum+=pow(10,(SPL_sbp[i][jk])/10.);

		}
		SPL_bp[jk]=10*log10(sum);

//		printf(" %ld  %8.4g\n",jk,SPL_bp[jk]);
	}


	printf("\n write output data...\n");

	double rms=0.;

	sum=0.;

	for(i=1; i<LAST_F_N; i++)
	{
		sum+=pow(10,(SPL_bp[i])/10.);

		fprintf(pFile[1],"\n \t %10.2lf \t %12.4e",fc[i],SPL_bp[i]);
		fprintf(pFile[2],"\n \t %10.2lf \t %12.4e",fc[i],SPL_bp[i]);
		fprintf(pFile[3],"\n \t %10.2lf \t %12.4e",fc[i],SPL_bp[i]+3.);  // elk

		printf("\n \t %10.2lf \t %8.4g",fc[i],SPL_bp[i]);
	}
	rms=10*log10(sum);


              printf("\n\n Overall Level = %12.4g dB \n",rms);
    fprintf(pFile[1],"\n\n Overall Level = %12.4g dB \n",rms);

	printf("\n");

	for(i=0;i<=ncoord;i++)
	{
		printf(" %ld  %8.4g  %8.4g \n",i,xxr[i],yyr[i]);
	}



}

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 input_data()
{

		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);
*/
		irad=2;


		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]);

		printf("\n");

		scanf("%d",&igeo);

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

		if(igeo !=1 && igeo !=2 && igeo !=3)
		{
			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");	

}
void onethird(void)
{

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

	LAST_F_N=28;
}
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*fc[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 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 step5_reference()
{   
	// power spectrum

	double f[50];

	double sn[50];
	double nrspl[50];

	double freq[50];

	double amp[50];

	double slope;
	double az;

	const double ref  = 1.0e-12;


	   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.;


	int num=13;


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


	freq[0]=20.;	
	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++)
	{
        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_ref[i]= 10.*log10(amp[i]/ref) +LW -10.*log10(U/de) +10.*log10( df);

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

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

	

}