disp(' ')
disp(' sine_to_narrowband.m   ver 1.1   August 20, 2013')
disp(' by Tom Irvine   Email: tom@vibrationdata.com')
disp(' ')
disp(' This program converts a sine tone into a narrowband PSD. ');
disp(' ');
disp(' The SDOF natural frequency is assumed to be the same as  ');
disp(' the sine tone frequency.');
disp('  ');
disp(' Suggested alpha values for fatigue exponent b ');
disp('  ');
disp('   alpha=1.9  for b>=6.4 ');
disp('   alpha=1.8  for b>=5.5 ');
disp('   alpha=1.7  for b>=4.5 ');
disp('   alpha=1.6  for b>=4.0 ');
disp(' ');
disp(' A smaller alpha yields a higher PSD level ');
%
close all;
%
clear ff;
clear fn;
clear amp;
clear Q;
clear damp;
clear nF;
%
disp(' ');
alpha=input(' Enter alpha ');
%
disp(' ');
fn=input(' Enter sine frequency (Hz) ');
amp=input(' Enter sine amplitude (G) ');
disp(' ')
Q=input(' Enter Q ');
damp=1/(2*Q);
%
disp(' ')
f1=input(' Enter band lower frequency (Hz) ');
f2=input(' Enter band upper frequency (Hz) ');
disp(' ')
%
nF=20;
%
df=(f2-f1)/(nF-1);
%   
ff=linspace(f1,f2,nF);
%
suma=0.; 
%
for j=1:nF 
%
    rho = ff(j)/fn;
	tdr=2.*damp*rho;
%
    c1= tdr^ 2.;
	c2= (1.- (rho^2.))^ 2.;
%
	trans = (1.+ c1 ) / ( c2 + c1 );
%
    opsd=trans*1;
    
    if(j>1 && j<nF)
        suma = suma + opsd*df;
    else
        suma = suma + opsd*df/2;        
    end    
%
end
a_vrs=sqrt(suma);
%       
psd_level=(Q*amp/(alpha*a_vrs))^2;
%
out1=sprintf('\n Narrowband PSD = %8.4g G^2/Hz ',psd_level);
disp(out1);