%
clear THF;
clear THS;
%
clear a1;
clear a2;
clear accel;
clear acceleration;
clear alpha;
clear amp;
clear b1;
clear b2;
clear beta;
clear damp_ratio;
clear disp;
clear displacement;
clear f;
clear fn;
clear omega;
clear ramp;
clear rrr;
clear rrrold;
clear sa;
clear sb;
clear sc;
clear shock_response;
clear srs_spec;
clear sz;
clear t;
clear td;
clear upper;
clear velocity;
clear velox;
clear wavelet_low;
clear wavelet_up;
clear x;
clear xabs;
clear xmax;
clear xmin;
clear yb; 
%
%
tpi=2*pi;
%
disp(' ');
disp(' th_from_wavelet_table.m  ver 1.3  January 1, 2011 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' This program generates an acceleration time ');
disp(' history from a wavelet table. ');
disp(' ');
%   
disp(' Select wavelet table input format ');
disp(' ');
disp(' 1=  index  freq(Hz)  accel(G)  number of half-sines  delay(sec)');
disp(' 2=  index  accel(G)  freq(Hz)  number of half-sines  delay(sec) ');
iform=input(' ');
%
disp(' ');
disp(' Select data input method ')
disp('  1=data preloaded in matlab ')  
disp('  2=external ASCII file ');
disp('  3=Excel file ');
%
    method=input('');                 
    disp(' ')  
%
    if(method==1)
        THF = input(' Enter the input filename:  ');    
    end
%    
    if(method==2)
        [filename, pathname] = uigetfile('*.*');
        filename = fullfile(pathname, filename);
        fid = fopen(filename,'r');
        THF=load(filename);            
    end 
%
    if(method==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
        THF = xlsread(xfile);
    end   
%
    sz=size(THF);
%
    if(iform==1)
        f=THF(:,2);
    	amp=THF(:,3);
    else
        amp=THF(:,2);
    	f=THF(:,3);
    end
%
    NHS=THF(:,4);
    td=THF(:,5);
%
    disp(' ');
    disp(' Enter time shift (sec) ');
    delay=input(' ');
    td=td+delay;
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if(min(f)<0)
        disp(' Frequency input error ');
        ppp=input(' ');
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
%
	last_wavelet=sz(1);
%
	out1=sprintf(' last_wavelet= %d \n',last_wavelet);
    disp(out1);  
    disp(' ');
    dur=input(' Enter duration (sec) ');
    disp(' ');
    sr=input(' Enter Sample rate '); 
%
dt=1/sr;
nt=round(dur/dt);
%
	out1=sprintf(' nt= %d \n',nt);
    disp(out1);
	out1=sprintf(' dt= %9.5g sec \n',dt);
    disp(out1);  
%
disp(' ');
disp(' multiply velocity and displacement by 386 ? ');
disp(' 1=yes 2=no ');
%
imul=input(' ');
%
isu=1;
if(imul==1)
	isu=386.;
end
disp(' ');
%
beta=tpi*f;
for(i=1:last_wavelet)
    alpha(i)=beta(i)/double(NHS(i));
    upper(i)=td(i)+(NHS(i)/(2.*f(i))); 
end
%
disp(' ref 1 ');
%
for(i=1:last_wavelet)
%    
    wavelet_low(i)=round( 0.5 +   (td(i)/dur)*nt);
     wavelet_up(i)=round(-0.5 +(upper(i)/dur)*nt);   
%    
    if(wavelet_low(i)==0)
        wavelet_low(i)=1;       
    end
    if(wavelet_up(i)>nt)
        wavelet_up(i)=nt;       
    end   
% 
end
%
disp(' ref 2 ');
%
t=linspace(0,nt*dt,nt);  
%
%
out=sprintf('\n\n Output Time History Files: \n');
out1=sprintf('      acceleration ');
out2=sprintf('      velocity ');
out3=sprintf('      displacement ');
disp(out);
disp(out1);
disp(out2);
disp(out3);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    accel=zeros(nt,1);
    velox=zeros(nt,1);
    disp=zeros(nt,1);   
%     
 	for(i=1:last_wavelet)
%%        progressbar(i/last_wavelet) % Update figure       
%
        sa=zeros(nt,1);
        sb=zeros(nt,1);
        sc=zeros(nt,1);
%
        ia=wavelet_low(i);
        ib=wavelet_up(i);
%       
        sa(ia:ib)=sin( alpha(i)*( t(ia:ib)-td(i) ) );
        sb(ia:ib)=sin(  beta(i)*( t(ia:ib)-td(i) ) );
        sc=amp(i)*sa.*sb;
%
		accel(ia:ib)=accel(ia:ib)+sc(ia:ib);
%
    end
%    
 	for(i=1:last_wavelet)    
        APB(i)=alpha(i)+beta(i);
        AMB(i)=alpha(i)-beta(i);
    end
%    
 	for(i=1:last_wavelet)
%
        sa=zeros(nt,1);
        sb=zeros(nt,1);
        sc=zeros(nt,1);  
%
        ia=wavelet_low(i);
        ib=wavelet_up(i);
%
        sa(ia:ib)=sin( APB(i)*( t(ia:ib)-td(i) ) )/APB(i);
        sb(ia:ib)=sin( AMB(i)*( t(ia:ib)-td(i) ) )/AMB(i);
        sc=amp(i)*(-sa+sb)*0.5;
%          
		velox(ia:ib)=velox(ia:ib)+sc(ia:ib);
%				
    end
%
 	for(i=1:last_wavelet)
%
        sa=zeros(nt,1);
        sb=zeros(nt,1);
        sc=zeros(nt,1);  
%
        ia=wavelet_low(i);
        ib=wavelet_up(i);
%
        sa(ia:ib)=(-1+cos(APB(i)*( t(ia:ib)-td(i) ) ))/((APB(i))^2);
        sb(ia:ib)=(-1+cos(AMB(i)*( t(ia:ib)-td(i) ) ))/((AMB(i))^2);
        sc=amp(i)*(sa-sb)*0.5;
%          
		disp(ia:ib)=disp(ia:ib)+sc(ia:ib);
%				
    end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=t';
figure(1);
plot(t,accel);
grid on;
title('Acceleration');
xlabel('Time (sec)');
acceleration=[t,accel];
%
velox=velox*isu;
disp=disp*isu;   
%
figure(2);
plot(t,velox);
xlabel('Time (sec)');
grid on;
title('Velocity');
velocity=[t,velox];
%
figure(3);
plot(t,disp);
xlabel('Time (sec)');
grid on;
title('Displacement');
displacement=[t,disp];
%
clear sa;
clear sb;
clear sc;
clear alpha;
clear beta;
clear t;
clear disp;
clear velox;
clear disp;
clear ia;
clear ib;
clear THF;
clear amp;
clear wavelet_low;
clear wavelet_up;
clear NHS;
clear f;
clear td;
clear AMB;
clear AMP;
clear out1;
clear out2;
clear out3;
clear out4;
clear isu;
clear i;
clear dur;
clear dt;
clear method;