disp(' ');
disp(' sine_to_psd_spec.m   ver 1.1   October 19, 2011 ');
disp(' ');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' This program calculates an equivalent PSD for a sine input, ');
disp(' where the PSD is intended to be narrowband ');
disp(' ');
%
clear t;
clear f;
clear n;
clear amp;
clear freq;
clear dur;
clear sr;
clear dt;
clear f1;
clear f2;
clear damp;
clear syn;
%
amp  =input(' Enter sine amplitude (G)  ');
freq =input(' Enter sine frequency (Hz) ');
dur  =input(' Enter duration (sec)      ');
disp(' ');
%
Q = input(' Enter the amplification factor for the SDOF response (typically Q=10) ');
nfe = input(' Enter the fatigue exponent (typical n=4) ');
%
f1=1.0e+99;
while(f1>=freq)
    disp(' ');
    out1=sprintf(' The starting frequency for the PSD narrowband must be < %6.4g Hz',freq);
    disp(out1);
    f1=input(' Enter the starting frequency(Hz) ');
end
%
f2=-1;
while(f2<=freq)
    disp(' ');
    out1=sprintf(' The ending frequency for the PSD narrowband must be > %6.4g Hz',freq);
    disp(out1);
    f2=input(' Enter the end frequency(Hz)  ');
end
%
damp=1./(2.*Q);
%
start_level = 1; % GRMS^2/Hz
start_rms=sqrt(1*(f2-f1));
%
out1=sprintf('\n Narrowband overall level = %8.4g GRMS \n',start_rms);
disp(out1);
%
for(kv=1:1)
%
start_rms = sqrt( start_level*(f2-f1));
%
[pps,syn,t,npoints,dt] = sine_to_psd_spec_TH(f1,f2,dur,start_rms);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[acceleration]=sine_to_psd_spec_ARBIT(npoints,freq,t,syn,dt,damp);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[Drandom] = sine_to_psd_spec_RF(acceleration,nfe);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
cycles = dur*freq;
Dsine = cycles*(amp*Q)^nfe;   % apply to one-sided amplitude
%
ratio=Dsine/Drandom;
nr=ratio^(1/nfe);
%%% NewAverageAmp=nr*AverageAmp;
%
intermediate_level = start_level*nr^2;
%
%%% out1=sprintf('\n  Dsine=%8.4g       Drandom=%8.4g ',Dsine,Drandom);
%%% disp(out1);
%%% out1=sprintf('\n  current_level=%8.4g  start_level=%8.4g  ratio=%8.4g',intermediate_level,start_level,ratio);
%%% disp(out1);
%
clear qax;
qax=[f1 intermediate_level;f2 intermediate_level];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
start_level=intermediate_level;
end
%
intermediate_RMS=sqrt(intermediate_level*(f2-f1));
%
sine_RMS=amp/sqrt(2);
%
if(intermediate_RMS<sine_RMS)
    r2=(sine_RMS/intermediate_RMS)^2;
    final_level=intermediate_level*r2;
else
    final_level=intermediate_level;   
end
qax=[f1 final_level;f2 final_level];
disp(' ');
disp(' * * * * * * * * * * * * * * * ');
disp(' ');
out1=sprintf(' Q=%7.3g, fatigue exponent=%6.4g  ',Q,nfe);
disp(out1);
disp(' ');
disp(' Sine Input ');
out1=sprintf(' %7.3g G, %7.4g Hz sine, %8.4g sec  ',amp,freq,dur);
disp(out1);
disp(' ');
disp(' Equivalent PSD ');
disp(' Freq(Hz) Accel(G^2/Hz)');
for(i=1:2)
    out1=sprintf(' %8.4g \t %7.4g',qax(i,1),qax(i,2));
    disp(out1);
end
final_RMS=sqrt(final_level*(f2-f1));
%
out1=sprintf('\n Overall level = %7.3g GRMS \n',final_RMS);
disp(out1);
%