disp(' ');
disp(' jet_mixing_noise.m   ver 1.0  by Tom Irvine');
disp(' ');
disp(' This program calculates the OASPL for a subsonic gas jet');
disp(' Reference:  Beranek & Ver, Noise & Vibration Control Engineering ');
disp(' ');
%
clear W;
clear M;
%
disp(' Select length unit ');
disp(' 1 = meters ');
disp(' 2 = inches ');
iul=input(' ');
disp(' ');
%
if(iul==1)
    DN=input(' Enter nozzle exit diameter (m) ');
else
    DN=input(' Enter nozzle exit diameter (inch) ');
    DN=DN/39.37;
end
%
disp(' ');
%
if(iul==1)
    R=input(' Enter distance from center of nozzle exit (m) ');
else
    R=input(' Enter distance from center of nozzle exit (inch) ');
    R=R/39.37;
end
%
disp(' ');
disp(' Select speed unit ');
disp(' 1 = m/sec ');
disp(' 2 = ft/sec ');
disp(' 3 = in/sec ');
disp(' 4 = miles/hr ');
iu=input(' ');
%
if(iu==1)
    unit='m/sec';
    out2=sprintf(' (Typically  343 m/sec) ');
end
%
if(iu==2)
    unit='ft/sec';
    out2=sprintf(' (Typically  1125 ft/sec) ');
end
%
if(iu==3)
    unit='in/sec';
    out2=sprintf(' (Typically  13500 in/sec) '); 
end
if(iu==4)
    unit='miles/hr';
    out2=sprintf(' (Typically  767 miles/hr) ');
end
%
out1=sprintf(' Enter the speed of sound (%s)',unit);
%
disp(out1);
disp(out2);
c=input(' ');
%
disp(' ');
out3=sprintf(' Enter the gas jet speed (%s)',unit);
disp(out3);
U=input(' ');;
%
M=U/c;
%
if(M<1)
%
    Mc=0.62*M;
%
    if(iu==2)  % ft/sec
        U=U/3.28;
    end
    if(iu==3)  % in/sec
        U=U/39.37;
    end
    if(iu==4)  % miles/hour
        U=U/2.24;
    end
%
    disp(' ');
    disp(' Select flow rate unit ');
    disp(' 1 = kg/sec ');
    disp(' 2 = lbm/sec ');
%
    iuf=input(' ');
%
    disp(' ');
%
    if(iuf==1)
        disp(' Enter the mass flow rate (kg/sec)');
    else
        disp(' Enter the mass flow rate (lbm/sec)');
    end
%
    m=input(' ');
%
    if(iuf==2)
        m=m/2.205;
    end
%    
    Wms=0.5*m*U^2;
%
    out0=sprintf('   M = %8.4g Mach   ',M);
    out1=sprintf('   U = %8.4g m/sec  ',U);
    out2=sprintf('   m = %8.4g kg/sec \n',m);
    out3=sprintf(' Wms = %8.4g watts  ',Wms);
    out4=sprintf('     = %8.4g kW  ',Wms/1000);  
%
    disp(out0);
    disp(out1);
    disp(out2);
    disp(out3);
    disp(out4);  
%
    disp(' ');
    disp(' Assume density ratio = 1 ');
    disp(' Assume ambient ground pressure.');
%    
    dens_ratio=1;
%
    w=(3*M^3.5)/(0.6 + M^3.5);
    w=w-1;
%
    Wratio=((4.5e-05)*(dens_ratio)^(w-1))*M^4.5;
    Wratio=Wratio/((1-Mc^2)^2);
%
    out4=sprintf('\n W/Wms = %8.4g  \n',Wratio);
    disp(out4);
%
    W=Wratio*Wms;
    out4=sprintf('   W = %8.4g watts  ',W);
    out5=sprintf('     = %8.4g kW  ',W/1000); 
    disp(out4);
    disp(out5);
    Wref=1.0e-12;
    dB=10*log10(W/Wref);
%    
    out4=sprintf('\n Overall Power Level = %8.4g dB   ref 1.0e-12 Watt \n',dB); 
    disp(out4);
%
    A= (1/4)*pi*DN^2;
%
%
    arg=M^7.5;
    arg=arg/(1 -0.1*M^2.5 +0.015*M^4.5);
    OASPL_90 = 139.5 + 10*log10(A/R^2) + 10*log10(arg);  
%
    disp(' ');
    disp(' *** SPL zero dB = ref 20 micro Pa *** ')
%
    for(j=1:350)
        theta_deg=((j-1)*0.5);
        theta=theta_deg*pi/180;
%
        theta_prime_deg=0.26*(180-theta_deg)*M^0.1;
%        
        arg1=1-(Mc*cos(theta)/(1+Mc^5)^0.2);
        arg2=1- (1/(10^(40.56-theta_prime_deg)+4e-06));
%        
        OASPL(j)=OASPL_90 -30*log10(arg1)-1.67*log10(arg2);
%
        if(j==1)
            OASPL_max=OASPL(j);
            theta_max=theta_deg;
        else
            if(OASPL(j)>OASPL_max)
                OASPL_max=OASPL(j);  
                theta_max=theta_deg;   
            end
        end
%
    end
%
    out4=sprintf('\n\n Maximum overall SPL occurs at %8.4g deg relative to jet axis',theta_max);
    disp(out4);   
    out5=sprintf('\n OASPL = %8.4g dB   at %6.3g deg relative to jet axis',max(OASPL),theta_max);    
    disp(out5);
    out6=sprintf('\n OASPL = %8.4g dB   at 90 deg relative to jet axis\n',OASPL_90);    
    disp(out6);      
%
else
    disp(' ');
    disp(' Error:  flow is supersonic ');
end