disp(' ')
disp(' disp_psd_syn_sine.m  version 1.0    January 11, 2012 ')
disp(' By Tom Irvine   Email:  tomirvine@aol.com ')
disp(' ')
%
clear a;
clear b;
clear aa;
clear bb;
clear sum;
clear amp;
clear ff;
clear freq;
clear fft_freq;
clear ph;
clear Y;
clear z;
clear psd_th;
clear slope;
clear spec;
clear THM;
clear TT;
clear zz;
%
rand('state',sum(100*clock));
tpi=2.*pi;
%
disp(' This program synthesizes a time history to satisfy a displacement PSD ');
disp(' ')
disp(' The Displacement PSD specification must be in a two-column matrix format: ')
disp(' Freq(Hz) & Disp(m^2/Hz) ')
disp(' ')
%
disp(' Select file input method ');
disp('   1=external ASCII file ');
disp('   2=file preloaded into Matlab ');
disp('   3=Excel file ');
file_choice = input('');
%
if(file_choice==1)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THM = fscanf(fid,'%g %g',[2 inf]);
        THM=THM';
end
if(file_choice==2)
        THM = input(' Enter the matrix name:  ');
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
amp=THM(:,2);
freq=THM(:,1);
%
disp(' Select output time history ');
disp('   1=displacement ');
disp('   2=acceleration ');
ith=input(' ');
%
disp(' ')
tmax = input('  Enter the duration (sec) ');
dur=tmax;
disp(' ');
%% df   = input('  Enter the delta frequency (Hz) ');
%
df=freq(1)/4;
%
damp=amp;
%
if(ith==2)
%   accel
    [dt,freq,amp,slope,arms,drms] = dpsd_syn_data_entry(freq,amp,tmax);
else
%   disp    
    [dt,freq,amp,slope,drms] = psd_syn_data_entry(freq,amp,tmax);    
end    
%
out4 = sprintf('\n  Input PSD overall level = %g meters RMS',drms);
disp(out4)
out5 = sprintf('\n  Input PSD overall level = %g (m/sec)^ RMS',arms);
disp(out5)
%
%  time vector
%
   nL=round(dur/dt);
%   
   TT=linspace(0,nL*dt,nL); 
%
   clear fss;
   clear samp;
   clear fslope;
%
   mmm=max(size(freq));
%
   samp=zeros(1,1+mmm);
   fss=zeros(1,1+mmm);  
   fslope=zeros(1,mmm);   
   samp(1)=amp(1)*0.7943;   % 3 db/oct
   fss(1)=freq(1)*(1/(2^(1/3)));  % 1/3 octave
   for(j=1:mmm)
       samp(j+1)=amp(j);
       fss(j+1)=freq(j);
   end
 %
   for(j=1:mmm-1)
       fslope(j+1)=slope(j);
   end
 %
   nnn=log(samp(2)/samp(1));
   ddd=log(fss(2)/fss(1));
   ratio=nnn/ddd;
   fslope(1)=ratio;
 %
   clear freq;
   clear amp;
   clear slope;
   freq=fss;
   amp=samp;
   slope=fslope;
%
%  frequency vector
%   
   f1=min(freq);
   f2=max(freq);
   nF=round((f2-f1)/df);
   ff=linspace(f1,nF*df,nF);
%
   ph=rand(1,nF)*tpi;
%
%
LS=max(size(freq))-1
spec=zeros(1,nF);
%
for(i=1:nF)
%
    for(j=1:LS)
 %            
	    if( ff(i) >= freq(j) && ff(i) < freq(j+1) )
%				
            fr = (ff(i)/freq(j) );
	        spec(i)=amp(j)*( fr^slope(j) );         
%
		    break;
        end
    end
%
    if( ff(i) == freq(j+1))
       spec(i)=amp(LS+1);    
    end
%
end
%
mag=sqrt(spec);
%
psd_th=zeros(1,nL);
%
ff=ff*tpi;
%
progressbar
for(j=1:nF)
    progressbar(j/nF)
    arg=ff(j)*TT+ph(j);
%
    psd_th=psd_th+mag(j)*sin(arg);
%
end
progressbar(1)
%
out5 = sprintf('\n scale time history ');
disp(out5)   
%
%  scale th
%
mu=mean(psd_th);
stddev=std(psd_th);
armsout = sqrt(mu^2+stddev^2);
%
if(ith==2)
    scale=arms/armsout;
else
    scale=drms/armsout;    
end
psd_th=psd_th*scale;
%
%  Output
%
mu=mean(psd_th);
sd=std(psd_th);
mx=max(psd_th);
mi=min(psd_th);
rms=sqrt(sd^2+mu^2);
%
disp(' ')
n=length(psd_th);
out0 = sprintf(' number of points = %d ',n);
out1 = sprintf(' mean = %8.4g    std = %8.4g    rms = %8.4g \n',mu,sd,rms);
out2 = sprintf(' max  = %9.4g  ',mx);
out3 = sprintf(' min  = %9.4g  \n',mi);
%
disp(out0);
disp(out1);
disp(out2);
disp(out3);
%
disp(' Plot time history in Matlab? ');
mchoice=input(' 1=yes  2=no  ' );
if mchoice == 1 
    figure(1);
    plot(TT,psd_th);
% 
    if(ith==2)
        out5 = sprintf(' Synthesized Acceleration Time History   std dev=%8.3g (m/sec^2) RMS',std(psd_th));   
        ylabel(' Accel (m/sec^2)');
    else
        out5 = sprintf(' Synthesized Displacement Time History   std dev=%8.3g meters RMS',std(psd_th));   
        ylabel(' Disp (m)');        
    end
    title(out5);
%   
    xlabel(' Time (sec)');
    disp(' ')
    disp(' Plot histogram? ')
    hchoice = input(' 1=yes 2=no ');
%
    if(hchoice == 1)
        clear x;
        xx=max(abs(psd_th));
        x=linspace(-xx,xx,15);       
        figure(2);    
        hist(psd_th,x)
        ylabel(' Counts');
        if(ith==2)
            xlabel(' Accel (m/sec^2)');  
        else
            xlabel(' Disp (m)');             
        end
        title(' Histogram ');  
    end   
end    
%
clear accel_time_history;
clear disp_time_history;
%
disp(' ')
if(ith==2)
    accel_time_history=[TT;psd_th]'; 
    disp(' The output acceleration time history is called "accel_time_history" ');    
else
    disp_time_history=[TT;psd_th]'; 
    disp(' The output disp time history is called "disp_time_history" ');
end