disp(' ')
disp(' psd_syn_sine.m  version 1.4  August 16, 2013 ')
disp(' By Tom Irvine   Email:  tomirvine@aol.com ')
disp(' ')
%
close all;
%
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;
%
fig_num=1;
%
tpi=2*pi;
%
disp(' This program synthesizes a time history to satisfy a PSD ');
disp(' ')
disp(' The PSD specification must be in a two-column matrix format: ')
disp(' Freq(Hz) & Accel(G^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(' ')
tmax = input('  Enter the duration (sec) ');
dur=tmax;
disp(' ');
df   = input('  Enter the delta frequency (Hz) ');
%
[dt,freq,amp,slope,grms] = psd_syn_sine_data_entry(freq,amp,tmax);
%
out4 = sprintf('\n  Input PSD overall level = %g ',grms);
disp(out4)
%
%  time vector
%
   nL=round(dur/dt);
%   
   TT=linspace(0,nL*dt,nL); 
%
%  frequency vector
%   
   f1=min(freq);
   f2=max(freq);
   nF=round((f2-f1)/df);
%   
   ff=linspace(f1,f1+nF*df,nF+1);
%
   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
%
%    out1=sprintf(' i=%d  ff=%8.4g  spec=%8.4g  ',i,ff(i),spec(i));
%    disp(out1);
%
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);
%%    out1=sprintf(' j=%d  ff=%8.4g  mag=%8.4g  ph=%8.4g ',j,ff(j),mag(j),ph(j));
%%    disp(out1);
%
    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);
    grmsout = sqrt(mu^2+stddev^2);
%
    scale=grms/grmsout;
    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(' ')
disp(' Select output options ');
disp(' ' );
disp(' Output file to ASCII text? ');
choice=input(' 1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
%%
      [writefname, writepname] = uiputfile('*','Save data as');
	  writepfname = fullfile(writepname, writefname);
	  writedata = [TT' psd_th' ];
	  fid = fopen(writepfname,'w');
	  fprintf(fid,'  %15.7e  %15.7g \n',writedata');
	  fclose(fid);
%%
end
%
disp(' Plot time history in Matlab? ');
mchoice=input(' 1=yes  2=no  ' );
if mchoice == 1 
    figure(fig_num);
    fig_num=fig_num+1;
    plot(TT,psd_th);
%    
    out5 = sprintf(' Synthesized Time History   std dev=%8.3g ',std(psd_th));
    title(out5);   
    ylabel(' Accel (G)');
    xlabel(' Time (sec)');
    disp(' ')
    disp(' Plot histogram? ')
    hchoice = input(' 1=yes 2=no ');
%
    if(hchoice == 1)
        [fig_num]=plot_histogram(psd_th,'Accel(G)',fig_num);
    end   
end    
%
clear time_history;
time_history=[TT;psd_th]'; 
disp(' ')
disp(' The output time history is called "time_history" ');