disp(' ');
disp(' propeller_psd_grms.m  ver 1.2  May 8, 2013 ');
disp(' ')
%
disp(' This program finds the overall GRMS of the MIL-STD-810G, ');
disp(' Figure 514.6D-2, Category 13, Propeller Aircraft Vibration ');
disp(' Exposure PSD. ');
%
clear f;
clear fl;
clear fu;
clear L;
%
close all;
%
%
P=0.010;
%
while(1)
    disp(' ');
    disp(' Select location ');
    disp('  1=In fuselage or wing forward of propeller ');
    disp('  2=Within one propeller blade radius of propeller passage plane');
    disp('  3=In fuselage or wing aft of propeller');
    disp('  4=In engine compartment, empennage, or pylons');
%
    isl=input(' ');
%
    if(isl==1 || isl==2 || isl==3 || isl==4 )
        break;
    end
%
end
%
fl=zeros(4,1);
fu=zeros(4,1);
L=zeros(4,1);
%
if(isl==1)
    L(1)=0.10;
end
%
if(isl==2)
    L(1)=1.20;
end
%
if(isl==3)
    L(1)=0.30;
end
%
if(isl==4)
    L(1)=0.60;
end
%
disp(' ');
f(1)=input(' Enter the blade passing frequency (Hz): ');
%
f(2)=2*f(1);
f(3)=3*f(1);
f(4)=4*f(1);
%
for i=2:4
   fr=f(i)/f(1);
   delta=-6;
   L(i)=L(1)*fr^((delta/10)/log10(2));
   if(L(i)<P)
       L(i)=P;
   end
end
%
for i=1:4
    fl(i)=f(i)*0.95;
    fu(i)=f(i)*1.05;
end
%
area=0;
for i=1:4
    df=fu(i)-fl(i);
    area=area+df*L(i);
end
%
df=fl(1)-15;
area=area+df*P;
%
for i=1:3
    df=fl(i+1)-fu(i);
    area=area+df*P;
end
%
df=2000-fl(4);
area=area+df*P;
%
grms=sqrt(area);
%
out1=sprintf('\n Overall level = %7.3g GRMS ',grms);
disp(out1);
%
clear ff;
clear aa;
%
ff=[15 fl(1) fl(1) fu(1) fu(1) ...
       fl(2) fl(2) fu(2) fu(2) ...
       fl(3) fl(3) fu(3) fu(3) ...       
       fl(4) fl(4) fu(4) fu(4) 2000];       
%
aa=[P P L(1) L(1) P ...
      P L(2) L(2) P ...
      P L(3) L(3) P ...       
      P L(4) L(4) P P]; 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear length;
n=length(ff);
%
disp(' ');
disp(' PSD ');
disp('  Freq      Accel ');
disp('  (Hz)    (G^2/Hz)');
%
for i=1:n
    out1=sprintf('  %5.1f    %6.4f  ',ff(i),aa(i));
    disp(out1);
end
%
clear propeller_psd;
propeller_psd=[ff' aa'];
%
disp(' ');
disp(' Output array:  propeller_psd ');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure(1);
plot(ff,aa);
out2=sprintf('Power Spectral Density  %7.3g GRMS overall',grms);
title(out2);
xlabel('Frequency (Hz) ');
ylabel('Accel (G^2/Hz)');
%
if(isl==2)
    axis([10,2000,0.001,10]);
else
    axis([10,2000,0.001,1]);  
end
%
set(gca,'XGrid','on','GridLineStyle',':');
set(gca,'YGrid','on','GridLineStyle',':');
set(gca,'xtick',[10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 2000])
set(gca,'XTickLabel',{'10';'';'';'';'';'';'';'';'';'100';'';'';'';'';'';'';'';'';'1000';'2000';})
set(gca,'MinorGridLineStyle',':','GridLineStyle',':','XScale','log','YScale','log','XminorTick','on','YminorTick','off');