%
disp(' ')
disp(' P9550.m   ver 1.0  June 27, 2011  ');
disp(' by Tom Irvine   ');
disp(' Email: tomirvine@aol.com  ');
%
disp(' ');
disp(' This program calculates the P95/50 level of a number of PSD functions. ');
disp(' P95/50 level is the 95th percentile value with a 50 percent confidence. ');
%
clear f;
clear amp;
clear level;
clear THM;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(' The input matrix must have a frequency column followed by ');
disp(' two or more PSD columns ');
disp(' ')
THM = input(' Enter the matrix name:  ');
%
sz=size(THM);
nf=sz(1);
num_col=sz(2)-1;
f=THM(:,1);
amp=THM(:,2:sz(2));
%
disp(' ');
out1=sprintf(' Number of amplitude columns = %d ',num_col);
disp(out1);
%
if(num_col<2)
    disp(' Input error.  Press Ctrl+C to exit ');
    uuu=input(' ');
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' ')
disp(' Select distribution ');
disp(' 1=normal 2=lognormal ');
idist=input(' ');
%
disp(' ')
disp(' Apply Barrett correction for skewness? ');
disp(' 1=yes  2=no ');
ibar=input(' ');
%
[ktol]=P9550_values(num_col);
%
out1=sprintf(' Number of amplitude columns = %d ',num_col);
disp(out1);
out1=sprintf(' ktol = %6.3g ',ktol);
disp(out1);
%
if(idist==2)
	amp=log10(amp);
end
%
level=zeros(nf,1);
%
for(i=1:nf)
%
    mu=mean(amp(i,:));
    sd= std(amp(i,:));
%   
    kt=0;
    sk=0;
    for j=1:num_col
        kt=kt+(amp(i,j)-mu)^4;
        sk=sk+(amp(i,j)-mu)^3;
    end      
    kt=kt/((num_col-1)*sd^4);
    sk=sk/((num_col-1)*sd^3);
%
    if(ibar==1)
       kbar=ktol*exp(0.255*sd-0.0155*(sd^2));
       level(i) = mu + kbar*sd;
    else
       level(i)=mu+ktol*sd;       
    end
%
	if(idist==2)
		level(i)=10^level(i);
    end
%
end
%
sz=size(f);
if(sz(2)>sz(1))
    f=f';
end
%
sz=size(f);
if(sz(2)>sz(1))
    level=level';
end
%
disp(' ')
disp(' Output file: p9550_psd ');
p9550_psd=[f level];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear grms;
clear n;
clear a;
clear s;
clear THF;
clear ytickn;
%
a=level;
n=nf;
%
disp(' ')
%
[s,grms]=calculate_PSD_slopes(f,a);
out5 = sprintf(' Overall Acceleration = %10.3g GRMS',grms);
disp(out5);
%
plot(f,a)
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log','XminorTick','off','YminorTick','off');
%
fmin=min(f);
fmax=max(f);
%
for(i=1:20)
    for(j=1:9)
        level = j*10^(i-10);
        fmax=level;
        if(level>=max(f))
            break;
        end
    end
    if(level>=max(f))
        break;
    end
end
%
for(i=1:20)
    for(j=1:9)
        level = j*10^(i-10);
        fmin=level;
        if(min(f)<=level)
            break;
        end
    end
    if(min(f)<=level)
            break;
    end  
end
%
if(fmin<10)
    fmin=10;
end
%
ymax=max(a);
ymin=min(a);
%
for(i=1:20)
    level = 10^(i-10);
    if(min(a) > level)
        ymin=level;
    end
end
%
for(i=1:20)
    level = 10^(10-i);
    if(max(a) < level)
        ymax=level;
    end
end
%
axis([fmin,fmax,ymin,(ymax*1.01)]);
%
ylabel('Accel (G^2/Hz)');   
xlabel('Frequency (Hz)');
out=sprintf(' Power Spectral Density %7.3g GRMS ',grms);
title(out);
%
set(gca,'XGrid','on','GridLineStyle',':');
set(gca,'YGrid','on','GridLineStyle',':');
%
if(fmin==20.)
    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';})
else
    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';})   
end
%
ya=ymin;
k=1;
for(i=1:10)
    for(j=1:9)
        ytickn(k)=j*ya; 
        k=k+1;
    end
    ya=ya*10.;
    ytickn(k)=ya; 
%    
    if(ya>ymax)
        break;
    end
end
%
set(gca,'ytick',ytickn);
%
nd=fix(log10(ymax/ymin));
%
string1=num2str(ymin);
string2=num2str(10*ymin);
string3=num2str(100*ymin);
string4=num2str(1000*ymin);
string5=num2str(10000*ymin);
string6=num2str(100000*ymin);
string7=num2str(1000000*ymin);
%
if(nd==1)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2})
end
if(nd==2)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2;'';'';'';'';'';'';'';'';string3})
end
if(nd==3)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2;'';'';'';'';'';'';'';'';string3;'';'';'';'';'';'';'';'';string4})
end
if(nd==4)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2;'';'';'';'';'';'';'';'';string3;'';'';'';'';'';'';'';'';string4;'';'';'';'';'';'';'';'';string5})
end
if(nd==5)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2;'';'';'';'';'';'';'';'';string3;'';'';'';'';'';'';'';'';string4;'';'';'';'';'';'';'';'';string5;'';'';'';'';'';'';'';'';string6})
end
if(nd==6)
    set(gca,'YTickLabel',{string1;'';'';'';'';'';'';'';'';string2;'';'';'';'';'';'';'';'';string3;'';'';'';'';'';'';'';'';string4;'';'';'';'';'';'';'';'';string5;'';'';'';'';'';'';'';'';string6;'';'';'';'';'';'';'';'';string7})
end