disp(' ')
disp(' arb_base_ode45.m ')
disp(' ver 1.0  December 2, 2011 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the response of a  ')
disp(' single-degree-of-freedom system to an arbitrary ')
disp(' base input using the ode45 method. ')
disp(' ')
disp(' The ode45 method is a numerical method.')
disp(' ')
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear n;
clear t;
clear x;
clear w;
clear y;
clear Y;
clear amp;
%
disp(' select units')
iunits=input(' 1=English  2=metric ');
%
if(iunits==1)
    disp(' ')
    v0 = input(' input initial relative velocity (in/sec) ');
    disp(' ')
    d0 = input(' input initial relative displacement (in) ');
else
    disp(' ')
    v0 = input(' input initial relative velocity (m/sec) ');
    disp(' ')
    d0 = input(' input initial relative displacement (m) ');    
end
%
disp(' ')
fn=input(' Enter the natural frequency (Hz)  ');
%
disp(' ')
idamp=input(' Enter damping format:  1= damping ratio   2= Q  ');	
%
disp(' ')
if(idamp==1)
    damp=input(' Enter damping ratio (typically 0.05) ');
else
    Q=input(' Enter the amplification factor (typically Q=10) ');
    damp=1./(2.*Q);
end
%
disp(' ')
%
%
omegan=2.*pi*fn;
omegad=omegan*sqrt(1.-damp^2);
%
%   State Space Solution
%
dt=(1/64)/fn;
sr=1/dt;
%
w(1)=d0;
w(2)=v0;
%
clear b;
clear c;
%
b=-omegan^2;
c=-2.*damp*omegan;
domegan=damp*omegan;
%
clear A;
clear B;
clear omega;
%
disp(' The input time history must have two columns: time(sec) & accel(G)')
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
%
t=THM(:,1);
y=THM(:,2);
%
if(iunits==1)
    y=y*386;
else
    y=y*9.81;
end
%
tmx=max(t);
tmi=min(t);
n = length(y);
%
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)
%
    disp(' ')
    disp(' begin FFT ')
    disp(' ')
%
    n = max(size(y));
%
    for(i=1:18)
       if( 2^i > n )
          break;
       end
       N=2^i;
    end
% 
    clear aaa;
    clear amp;
    aaa=y(1:N);
    amp=flipud(aaa);
    Y = fft(amp,N);
%
    clear complex_FFT;
    complex_FFT=zeros(N,3);
    clear FF;
%
    df=1./(N*dt);
    out4 = sprintf(' df  = %8.4g Hz  ',df);
    disp(out4)
%
    FF=linspace(0,df*(N-1),N);
%
    complex_FFT(:,1)=FF';
    complex_FFT(:,2)=real(Y)/N;
    complex_FFT(:,3)=imag(Y)/N;
%
    NH=round(N/2)-1;
%
    clear omega;
    clear A;
    clear B;
%
    omega=2*pi*complex_FFT(1:NH,1);
    A=complex_FFT(1:NH,2);
    B=complex_FFT(1:NH,3);
    A=A*2;
    B=B*2;
%
%
    tstart=0;
    tend=N*dt;
    [t,w] = ode45(@arb_base_function,[tstart tend],w,[],[],b,c,A,B,omega,NH);
%
    figure(1);
    plot(t,w(:,1));
    title('Relative Displacement Response');
    if(iunits==1)
       ylabel('Relative Displacement (in)');
    else
       ylabel('Relative Displacement (m)');
    end
    xlabel('Time (sec)');
    grid on;
%
    clear xdd;
    xdd=c*w(:,2)+b*w(:,1);
% 
    if(iunits==1)
       xdd=xdd/386;
    else
       xdd=xdd/9.81;        
    end
%
    figure(2);
    plot(t,xdd);
    title('Acceleration Response');
    if(iunits==1)
        ylabel('Accel (G)');
    else
        ylabel('Accel (m/sec^2)');
    end
    xlabel('Time (sec)');
    grid on;
%
  disp(' ');
  disp(' Relative Response ');
  if(iunits==1)
  out10=sprintf('\n  maximum = %10.2f in',max(w(:,1)));
    out11=sprintf('  minimum = %10.2f in',min(w(:,1)));
  else
  out10=sprintf('\n  maximum = %10.2f mm',1000*max(w(:,1)));
    out11=sprintf('  minimum = %10.2f mm',1000*min(w(:,1)));      
  end
  disp(out10);
  disp(out11);
%
  disp(' ');
  disp(' Acceleration Response ');
  out10=sprintf('\n  maximum = %10.2f G',max(xdd));
    out11=sprintf('  minimum = %10.2f G',min(xdd));
  disp(out10);
  disp(out11);
%
end