disp(' ')
disp(' psd_integ.m   ver 1.9   October 3, 2012')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This program calculates the overall level of ')
disp(' a power spectral density function. ')
disp(' ')
%
close all;
clear f;
clear a;
clear v;
clear d;
clear input_amp;
clear s;
clear omega;
%
disp(' Select units ');
disp('  1=English: accel(G),       vel(in/sec), disp(in)  ');
disp('  2=metric : accel(G),       vel(m/sec),  disp(mm)  ');
disp('  3=metric : accel(m/sec^2), vel(m/sec),  disp(mm)  ');
iunit=input(' ');
%
if(iunit==1)
   accel_RMS_label='GRMS';
    accel_3s_label='G 3-sigma';
     vel_RMS_label='in/sec RMS';
      vel_3s_label='in/sec 3-sigma';
    disp_RMS_label='in RMS';
     disp_3s_label='in 3-sigma';  
end
%
if(iunit==2)
   accel_RMS_label='GRMS';
    accel_3s_label='G 3-sigma';
     vel_RMS_label='m/sec RMS';
      vel_3s_label='m/sec 3-sigma';
    disp_RMS_label='mm RMS';
     disp_3s_label='mm 3-sigma';  
end
%
if(iunit==3)
   accel_RMS_label='m/sec^2 RMS';
    accel_3s_label='m/sec^2 3-sigma';
     vel_RMS_label='m/sec RMS';
      vel_3s_label='m/sec 3-sigma';
    disp_RMS_label='mm RMS';
     disp_3s_label='mm 3-sigma';    
end
%
disp(' Select input PSD type ');
disp('  1=acceleration ');
disp('  2=velocity     ');
disp('  3=displacement ');
itype=input(' ');
%
    if(itype==1)
        if(iunit==1)
            PSD_label='Accel(G^2/Hz)';
        end
        if(iunit==2)
            PSD_label='Accel(G^2/Hz)';            
        end  
        if(iunit==3)
            PSD_label='Accel(m/sec^2)^2/Hz';
        end        
    end
    if(itype==2)
        if(iunit==1)
            PSD_label='Vel(in/sec)';
        end
        if(iunit==2)
            PSD_label='Vel(m/sec)';            
        end  
        if(iunit==3)
            PSD_label='Vel(m/sec)';
        end             
    end
    if(itype==3)
        if(iunit==1)
            PSD_label='Disp(in)';
        end
        if(iunit==2)
            PSD_label='Disp(mm)';            
        end  
        if(iunit==3)
            PSD_label='Disp(mm)';
        end        
    end
%
disp(' Select data input method: ')
disp('  1=Open external ASCII file    ');
disp('  2=Preloaded input file  ');  
disp('  3=Keyboard entry        ');
disp('  4=Excel file ');
m = input('');
disp(' ')
%
if(m~=1 && m~=2 && m~=4)
    m=3;
end
%
if(m==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename); 
    fid = fopen(filename,'r');
    THF = fscanf(fid,'%g %g',[2 inf]);
    THF=THF';   
    f=THF(:,1);
    input_amp=THF(:,2);
    n=max(size(f)); 
%
        if(itype==1)        
            a=input_amp;
        end 
        if(itype==2)
            v=input_amp;            
        end
        if(itype==3)
            d=input_amp;         
        end
%
end
if(m==2)
    disp(' The input file must be pre-loaded into Matlab.')
    disp(' ')
%    
    out1=sprintf(' The input file must have two columns: frequency(Hz) & %s',PSD_label);
    disp(out1);
%    
    disp(' ')
    THF = input(' Enter the input PSD filename:  ');
    f=THF(:,1);
    input_amp=THF(:,2);
    n=max(size(f));
%
        if(itype==1)        
            a=input_amp;
        end 
        if(itype==2)
            v=input_amp;            
        end
        if(itype==3)
            d=input_amp;         
        end
%
end
if(m==3)
    n= input(' Enter the number of input points: ');
    disp(' ')
    out7=sprintf(' The frequency unit is (Hz).  The amplitude unit is: %s ',PSD_label);
    sprintf(out7);
    disp(' ')
    for i=1:n
        out1=sprintf(' Enter frequency %d:  ',i);
        f(i)=input(out1);
        out2=sprintf(' Enter amplitude %d:  ',i);
        input_amp(i)=input(out2);  
        disp(' ')
    end
%
        if(itype==1)        
            a=input_amp;
        end 
        if(itype==2)
            v=input_amp;            
        end
        if(itype==3)
            d=input_amp;         
        end
%
end
if(m==4)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THF = xlsread(xfile);
        f=THF(:,1);
        n=max(size(f));        
        if(itype==1)        
            a=THF(:,2);
        end 
        if(itype==2)
            v=THF(:,2);            
        end
        if(itype==3)
            d=THF(:,2);         
        end        
%         
end
%
omega=2*pi*f;
%
%
%% disp(' Select units ');
%% disp(' 1=English: accel(G),       vel(in/sec), disp(in)  ');
%% disp(' 2=metric : accel(G),       vel(m/sec),  disp(mm)  ');
%% disp(' 3=metric : accel(m/sec^2), vel(m/sec),  disp(mm)  ');
%% iunit=input(' ');
%
%
%   Convert acceleration to G;
%
for i=1:n
%
        if(itype==1)  %% accel PSD input
%
            if(iunit==3)
                a(i)=a(i)/9.81^2;
            end
 %           
        end
%
        if(itype==2)  %% vel PSD input
            a(i)=v(i)*omega(i)^2;
%            
            if(iunit==1)
                a(i)=a(i)/386^2;
            end
            if(iunit==2)
                a(i)=a(i)/9.81^2;
            end
 %
            if(iunit==3)
                a(i)=a(i)/9.81^2;
            end
 %           
        end
        if(itype==3)   %% disp PSD input
            a(i)=d(i)*omega(i)^4;  
%            
            if(iunit==1)
                a(i)=a(i)/386^2;
            end
            if(iunit==2)
                a(i)=a(i)/(9.81*1000)^2;
            end
            if(iunit==3)
                a(i)=a(i)/(9.81*1000)^2;
            end
 %                  
        end    
%        
end
%
%   a is now in G^2;
%
disp(' ')
%
[s,grms]=calculate_PSD_slopes(f,a);
%
if(iunit==3)
    grms=grms*9.81;
end
out5 = sprintf('\n Overall Acceleration = %10.3g %s',grms,accel_RMS_label);
disp(out5);
out8 = sprintf('                      = %10.3g %s',3*grms,accel_3s_label);
disp(out8);
%
vrms = APSD_VPSD(f,a);
%
meters_per_inch=0.0254;
mm_per_inch=25.4;
%
if(iunit==2 || iunit==3)
    vrms=vrms*meters_per_inch;
end
%
out6 = sprintf('\n Overall Velocity     = %10.3g %s',vrms,vel_RMS_label);
disp(out6);
%
drms = APSD_DPSD(f,a);
%    
if(iunit==2 || iunit==3)
    drms=drms*mm_per_inch;
end
%
out7 = sprintf('\n Overall Displacement = %10.3g %s',drms,disp_RMS_label);
disp(out7);
out7 = sprintf('                      = %10.3g %s\n',3*drms,disp_3s_label);
disp(out7);
%  
disp(' ')
disp(' Plot input PSD ?');
choice = input(' 1=yes  2= no  ');
%
if(itype==1)
    ymin=min(a);
    ymax=max(a);
end
%
if(itype==2)
    ymin=min(v);
    ymax=max(v);
end
%
if(itype==3)
    ymin=min(d);
    ymax=max(d);
end
%
pymax=ymax;
pymin=ymax/100.;
%
for i=20:-1:-20
    if(ymax <= 10^i )
        pymax = 10^i;
    end
end
for i=20:-1:-20
    if(ymin <= 10^i )
        pymin = 10^(i-1);
    end
end
%
sz=size(f);
if(sz(2)>sz(1))
    f=f';
end
%
if(choice==1)
   figure(1);
%   
   if(itype==1)
        sz=size(a);
        if(sz(2)>sz(1))
            a=a';
        end
        plot(f,a);
   end
 %
    if(itype==2)
        sz=size(v);
        if(sz(2)>sz(1))
            v=v';
        end
        plot(f,v);
   end
 %
    if(itype==3)
        sz=size(d);
        if(sz(2)>sz(1))
            d=d';
        end
        plot(f,d);
   end
 %
   out1=sprintf('%s',PSD_label);
   ylabel(out1);
   xlabel('Frequency (Hz)');
   if(itype==1)
        out=sprintf(' Power Spectral Density %7.3g %s ',grms,accel_RMS_label);
   end
   if(itype==2)
        out=sprintf(' Power Spectral Density %7.3g %s ',vrms,vel_RMS_label);
   end   
   if(itype==3)
        out=sprintf(' Power Spectral Density %7.3g %s ',drms,disp_RMS_label);
   end
   title(out);
   grid;
   set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
   axis([f(1),f(n),pymin,pymax]);
   if(f(n)==2000 && f(1)==10)
        set(gca,'MinorGridLineStyle','none', 'XTickMode', 'manual');
        set(gca,'xtick',[10 100 1000 2000]);
   end
   if(f(n)==2000 && f(1)==20)
        set(gca,'MinorGridLineStyle','none', 'XTickMode', 'manual');
        set(gca,'xtick',[20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 2000]);
        set(gca,'XTicklabel',{'20';'';'';'';'';'';'';'';'100';'';'';'';'';'';'';'';'';'1000';'2000'});
        h=gca;
        set(h,'XMinorTick','off');
   end
end
%
%  Rename matlab matrix with descriptive name
%
    clear aa;
    aa=size(f);
    clear power_spectral_density;
% 
    if(itype==1)
        power_spectral_density=[f a];
    end
    if(itype==2)
        power_spectral_density=[f v];
    end
    if(itype==3)
        power_spectral_density=[f d];
    end  
%    
    out5 = sprintf('\n\n The PSD matrix is renamed as "power_spectral_density" ');
    disp(out5);   
%  