disp(' ');
disp(' helicopter.m  ver 1.1  April 8, 2013 ');
disp(' ');
disp(' by Tom Irvine ');
disp('  ');
disp(' MIL-STD-180G,Table 514.6C-X, Helicopter Vibration Exposure ');
disp(' Vibration induced by the main rotor and blade passing frequencies. ');
disp(' ');
%
close all;
%
fig_num=1;
%
while(1)
%
    clear f;
    clear a;
%
    f=zeros(4,1);
    a=zeros(4,1);
%
    while(1)
        disp(' Select Excitation Source');
        disp('  1=Main Rotor ');
        disp('  2=Tail Rotor ');
        disp('  3=Drive Train ');
        is=input(' ');
        disp(' ');    
        if(is==1 || is==2 || is==3)
            break;
        end
    end
%
    if(is==1 || is==2)
%
        if(is==1)
            fr= input(' Enter Main Rotor Hub Frequency (Hz) ');        
        end
        if(is==2)
            fr= input(' Enter Tail Rotor Hub Frequency (Hz) ');          
        end    
%
        disp(' ');
        n= input(' Enter number of blades ');
%
        f(1)=    fr;
        f(2)=  n*fr;
        f(3)=2*n*fr;
        f(4)=3*n*fr;
%
        while(1)
            disp(' ');
            disp(' Selected location ');
            disp('  1=General ');
            disp('  2=Instrument Panel ');
            disp('  3=External Stores ');
            iloc=input(' ');
%
            if(iloc==1 || iloc==2 || iloc==3)
                break;
            end
%
        end
%
        disp(' ');
        disp('  Sine Tones ');
        disp(' ');
        disp('  Freq(Hz)  Accel(G) ');
%
        for i=1:4
            a(i)=0.;
%
            if(iloc==1) % General
%
                if( f(i)>= 3 && f(i) <10)
                    a(i)=0.70/(10.70-f(i));
                end
                if( f(i)>= 10 && f(i) <25)
                    a(i)=0.10*f(i);
                end
                if( f(i)>= 25 && f(i) <40)
                    a(i)=2.50;
                end
                if( f(i)>= 40 && f(i) <50)
                    a(i)=6.50 - 0.10*f(i);
                end
                if( f(i)>= 50 && f(i) <=500)
                    a(i)=1.50;
                end
                W0=0.0010;
                W1=0.010;
                FT=500;
                grms=1.8;
%     
            end
            if(iloc==2) % Instrument Panel
%
                if( f(i)>= 3 && f(i) <10)
                    a(i)=0.70/(10.70-f(i));
                end
                if( f(i)>= 10 && f(i) <25)
                    a(i)=0.070*f(i);
                end
                if( f(i)>= 25 && f(i) <40)
                    a(i)=1.750;
                end
                if( f(i)>= 40 && f(i) <50)
                    a(i)=4.550 - 0.070*f(i);
                end
                if( f(i)>= 50 && f(i) <=500)
                    a(i)=1.050;
                end
                W0=0.0010;
                W1=0.010;
                FT=500;  
                grms=1.8;
%             
            end
            if(iloc==3) % External Stores
                if( f(i)>= 3 && f(i) <10)
                    a(i)=0.70/(10.70-f(i));
                end
                if( f(i)>= 10 && f(i) <25)
                    a(i)=0.150*f(i);
                end
                if( f(i)>= 25 && f(i) <40)
                    a(i)=3.750;
                end
                if( f(i)>= 40 && f(i) <50)
                    a(i)=9.750 - 0.150*f(i);
                end
                if( f(i)>= 50 && f(i) <=500)
                    a(i)=2.250;
                end
                W0=0.0020;
                W1=0.020;
                FT=500;          
                grms=2.5;
%          
            end
%
            if(f(i)< 3)		
                disp(' frequency is too low. ');
            end		
            if(f(i)>500)		
                disp(' frequency is too high. ');
            end  
% 
            if(a(i)>0)
                out1=sprintf(' %8.4g   %6.2g  ',f(i),a(i));
                disp(out1);
            end
%
        end
%
    end
    if(is==3)
%        
        fr=input(' Enter drive frequency (Hz) ');
%
        f(1)=fr;
        f(2)=2*fr;
        f(3)=3*fr;
        f(4)=4*fr;
%
        disp(' ');
        disp('  Sine Tones ');
        disp(' ');
        disp('  Freq(Hz)  Accel(G) ');
%
        for i=1:4
            a(i)=0;
            if(f(i)<5)
                disp(' Frequency too low ');
            end
            if(f(i)>=5 && f(i)<=50)
                a(i)=0.1*f(i);
            end
            if(f(i)>50 && f(i)<=2000)
                a(i)=5.0+0.010*f(i);
            end
            if(f(i)>2000)
                disp(' Frequency too high ');
            end
%
            if(a(i)>0)
                out1=sprintf(' %8.4g   %6.2g  ',f(i),a(i));
                disp(out1);
            end
%
        end
        W0=0.0020;
        W1=0.020;
        FT=2000;         
        grms=3.8;        
%
    end
%
    helicopter_psd=[ 10  W0 ; 100 W1  ; 300 W1 ; FT W0 ];
%
    out1=sprintf('\n Random Vibration Overall Level = %8.4g GRMS\n',grms);
    disp(out1);
%
    figure(fig_num);
    fig_num=fig_num+1;
%    
    ff=[ f(1) f(1) f(1) ...
         f(2) f(2) f(2) ...
         f(3) f(3) f(3) ...
         f(4) f(4) f(4)    ];
%     
    aa=[0.001 a(1) 0.001 ...
        0.001 a(2) 0.001 ...
        0.001 a(3) 0.001 ...
        0.001 a(4) 0.001  ];
%    
    [AX,hline1,hline2] = plotyy(helicopter_psd(:,1),helicopter_psd(:,2),ff,aa,'plot','plot');
%
%
    set(AX(1),'Xscale','log','XminorTick','off','MinorGridLineStyle','none','xtick',[]);
    set(AX(1),'Yscale','log');
%    
    set(AX(2),'XScale','log','XminorTick','off','MinorGridLineStyle','none','xtick',[]);
    set(AX(2),'YScale','log');
%    
    set(AX(2),'YLim',[0.01 10]); 
    set(AX(2),'ytick',[0.01 0.1 1 10]);    
%
    if(is==3)
       set(AX(1),'XLim',[2 2000]) 
       set(AX(2),'XLim',[2 2000])
       set(gca,'xtick',[2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 2000])
       set(gca,'XTickLabel',{'2';'';'';'';'';'';'';'';'10';'';'';'';'';'';'';'';'';'100';'';'';'';'';'';'';'';'';'1000';'2000';})      
%
       set(AX(1),'YLim',[0.001 1]); 
       set(AX(1),'ytick',[0.001 0.01 0.1 1]);
%
    else
       set(AX(1),'XLim',[2 500])   
       set(AX(2),'XLim',[2 500]) 
       set(gca,'xtick',[2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 200 300 400 500])
       set(gca,'XTickLabel',{'2';'';'';'';'';'';'';'';'10';'';'';'';'';'';'';'';'';'100';'';'';'';'500';})             
%
       set(AX(1),'YLim',[0.001 1]); 
       set(AX(1),'ytick',[0.001 0.01 0.1 1 10]);
%
    end
    set(AX(1),'YGrid','on','GridLineStyle',':');
    xlabel('Frequency (Hz)');
    title('Sine-on-Random Helicopter Vibration Exposure ');
    set(get(AX(1),'Ylabel'),'String','PSD Accel(G^2/Hz)'); 
    set(get(AX(2),'Ylabel'),'String','Accel(G)');
%
    disp(' ');
    disp(' Another case?  1=yes 2=no ');
    ic=input(' ');
%    
    if(ic==2)
        break;
    end
end
%