disp(' ')
disp(' integrate_th.m   ver 1.5   August 16, 2011')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
disp(' This script integrates a time history via the trapezoidal rule. ');
disp(' ')
disp(' The input time history must have two columns: time(sec) & amplitude')
disp(' ')
%
clear velocity;
clear displacement;
clear integrated;
clear integrated_sm;
clear THM;
clear t;
clear y;
clear v;
clear vv;
clear c;
clear a;
figure_number=1;
%
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
%
disp(' ');
disp(' Enter input signal type ');
disp('  1=acceleration ');
disp('  2=velocity ');
ist=input(' ');
%
t=THM(:,1);
y=THM(:,2);
%
figure(figure_number);
figure_number=figure_number+1;
plot(t,y);
grid on;
if(ist==1)
    title(' Acceleration, Raw Input Signal');
else
    title(' Velocity, Raw Input Signal');   
end
xlabel('Time (sec)');
%
tmx=max(t);
tmi=min(t);
n = max(size(y));
a=zeros(n,1);
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
disp(' ')
disp(' Time Step ');
dtmin=min(diff(t));
dtmax=max(diff(t));
%
out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
out5 = sprintf(' dt     = %8.4g sec  ',dt);
out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
disp(out4)
disp(out5)
disp(out6)
%
disp(' ')
disp(' Sample Rate ');
out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
disp(out4)
disp(out5)
disp(out6)
%
ncontinue=1;
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.  Continue calculation? 1=yes 2=no ')
    ncontinue=input(' ');
end
%
if(ncontinue==1)
    for(j=1:2)
        disp(' ');
        disp(' Select trend removal prior to integration ')
        disp('   1=mean remove ');
        disp('   2=mean filtered removal ');
        disp('   3=none ');
        irm=input(' ');
        if(irm==1)
            y=y-mean(y);
%
            figure(figure_number);
            figure_number=figure_number+1;
            plot(t,y);
            if(ist==1)
                title(' Acceleration with Mean Removed Prior to Integration');
            else
                title(' Velocity with Mean Removed Prior to Integration');   
            end
            grid on;
            xlabel('Time (sec)');
        end
%
        if(irm==2)
            [y,figure_number]=integrate_th_mf(1,ist,t,y,figure_number);
        end
%
        [y,figure_number]=fade(1,ist,t,y,figure_number);
%
        n=max(size(y));
        v=zeros(n,1);
        v(1)=y(1)*dt/2;
        for(i=2:(n-1))
            v(i)=v(i-1)+y(i)*dt;
        end
        v(n)=v(n-1)+y(n)*dt/2;
%
        disp(' ');
        disp(' multiply by constant?  1=yes  2=no ');
        im=input(' ');
        if(im==1)
            disp(' ');
            c=input(' Enter constant  ');
            v=v*c;
        end
%
        figure(figure_number);
        figure_number=figure_number+1;
%        
        plot(t,v);
        xlabel('Time (sec)');
        y=v';
        if(ist==1)
            title(' Velocity via Integration');
        else
            title(' Displacement via Integration');   
        end     
        grid on;
        vmax=max(v);
        vmin=min(v);
        disp(' ');
        if(ist==1)
            disp(' Velocity via Integration');
        else
            disp(' Displacement via Integration');   
        end   
        out1=sprintf('\n max=%8.4g  min=%8.4g ',vmax,vmin);
        disp(out1);
%
        integrated(:,1)=t';
        integrated(:,2)=v';
%
        disp(' ');
        disp(' The output time history is called: intergated ');
%
        [v,figure_number]=integrate_th_mf(2,ist,t,v,figure_number);
%
        [v,figure_number]=fade(2,ist,t,v,figure_number);
        y=v';  
%%
        sz=size(t);
        if(sz(2)>sz(1))
            t=t';
        end
        sz=size(v);
        if(sz(2)>sz(1))
            v=v';
        end    
%%
        disp(' ');
        if(ist==1)
            velocity=zeros(n,2);
            velocity(:,1)=t;
            velocity(:,2)=v;
        else
            displacement=zeros(n,2);
            displacement(:,1)=t;
            displacement(:,2)=v;           
        end
%
        if(ist==2)
            break;
        end        
%
        disp(' ');
        disp(' integrate again? ');
        disp('  1=yes  2=no ');
        iai=input(' ');
        if(iai==2)
            break;
        end
        ist=2;
        clear v;
        clear vv;
        clear c;
        clear a;
    end
end