%
%  mdof_arbit_base_core.m  ver 1.7  January 14, 2014
%
function[t,acc,vel,rel_disp,abase]=...
    mdof_arbit_base_core(fn,damp,PF,ModeShapes,T1,T2,ngw,iu)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)
        FS = input(' Enter the matrix name:  ','s');
        THM=evalin('caller',FS);
end
if(file_choice==3)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        THM = xlsread(xfile);
%         
end
%
t=THM(:,1);
f=THM(:,2);
abase=f;
%
tmx=max(t);
tmi=min(t);
n = length(f);
nt=n;
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)
%
if(((dtmax-dtmin)/dt)>0.01)
    disp(' ')
    disp(' Warning:  time step is not constant.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
sz=size(ModeShapes);
%
num_dof=sz(1);
num_modes=sz(2);
%
Q=ModeShapes;
MST=ModeShapes';
%
%
disp(' ');
out1=sprintf(' Enter number of modes to include in analysis (max=%d) ',num_modes);
disp(out1);
ninc=input('  ');
%
if(ninc>num_modes)
   ninc=num_modes;
end
%
num_modes=ninc;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Calculate Filter Coefficients
%
disp(' ');
disp(' Calculate Filter Coefficients ');
% 
omegan=2*pi*fn;
[a1,a2,df1,df2,df3,vf1,vf2,vf3,af1,af2,af3]=...
              ramp_invariant_filter_coefficients(num_modes,omegan,damp,dt);
%
%  Numerical Engine
%
disp(' ')
disp(' Calculating response...');
%
nx=zeros(nt,num_modes);
nv=zeros(nt,num_modes);
na=zeros(nt,num_modes);
%
out1=sprintf('\n nt=%d  num_modes=%d\n',nt,num_modes);
disp(out1);
%
progressbar;
%
for j=1:num_modes
%
   progressbar(j/num_modes);
%
   amodal=-PF(j)*abase;
%
%  displacement
%
    d_forward=[   df1(j),  df2(j), df3(j) ];
    d_back   =[     1, -a1(j), -a2(j) ];
    d_resp=filter(d_forward,d_back,amodal);
%    
%  velocity
%
    v_forward=[   vf1(j),  vf2(j), vf3(j) ];
    v_back   =[     1, -a1(j), -a2(j) ];
    v_resp=filter(v_forward,v_back,amodal);
%    
%  acceleration
%   
    a_forward=[   af1(j),  af2(j), af3(j) ];
    a_back   =[     1, -a1(j), -a2(j) ]; 
    a_resp=filter(a_forward,a_back,amodal);
%
    nx(:,j)=d_resp;  % displacement
    nv(:,j)=v_resp;  % velocity
    na(:,j)=a_resp;  % acceleration  
%
end
clear d_resp;
clear v_resp;
clear a_resp;
pause(0.5);
progressbar(1);
%
Qs=Q(:,1:num_modes);
%
clear Q;
%
%
disp(' Calculate d');
%
d=zeros(nt,num_dof);
%
progressbar;
for i=1:nt
    progressbar(i/nt);
      d(i,1:num_dof)=(Qs*(nx(i,1:num_modes))')';  
end
clear nx;
pause(0.5);
progressbar(1);
%
%
disp(' Calculate v');
%
v=zeros(nt,num_dof);
progressbar;
for i=1:nt
    progressbar(i/nt);
      v(i,1:num_dof)=(Qs*(nv(i,1:num_modes))')';
end
clear nv;
pause(0.5);
progressbar(1);
%
%
disp(' Calculate acc');
%
acc=zeros(nt,num_dof);
progressbar;
for i=1:nt
    progressbar(i/nt);
    acc(i,1:num_dof)=(Qs*(na(i,1:num_modes))')';    
end
clear na;
pause(0.5);
progressbar(1);
%
clear Qs;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Transformation back to xd xf
%
%%    dT=TT*d';
%%    vT=TT*v';
%%    accT=TT*acc';
%
disp(' Transformation back to xd xf ');
%
dvelox=dt*cumtrapz(abase);
ddisp=dt*cumtrapz(dvelox);    
%
sz=size(T1);
base_ddisp_th =(ones(sz(2),1))*ddisp'; 
%
base_dvelox_th=(ones(sz(2),1))*dvelox';
%
base_drive_th =(ones(sz(2),1))*abase';
%
  dT=(T1*base_ddisp_th  + T2*d')';
  vT=(T1*base_dvelox_th + T2*v')';
accT=(T1*base_drive_th  + T2*acc')';
%
clear rd;
clear rv;
clear ra;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Put in original order
%
disp(' Put in original order ');
%
ZdT=[ base_ddisp_th'   dT ];
clear base_ddisp_th;
clear dT;
%
ZvT=[ base_dvelox_th'  vT ];
clear base_dvelox_th;
clear vT;
%
ZaccT=[ base_drive_th'  accT ];
clear base_drive_th;
clear accT;
%
clear d;   
clear v;
clear acc;
%
sz=size(ZdT);
num_dof=sz(2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
d=zeros(nt,num_dof);
%  
progressbar;
for i=1:num_dof
    progressbar(i/num_dof);
%
        d(1:nt,ngw(i))=    (ZdT(1:nt,i));        
%
end
clear zdT;
pause(0.5);
progressbar(1);
%
v=zeros(nt,num_dof);
%  
progressbar;
for i=1:num_dof
    progressbar(i/num_dof);
%       
        v(1:nt,ngw(i))=    (ZvT(1:nt,i));     
%
end
clear zvT;
pause(0.5);
progressbar(1);
%
acc=zeros(nt,num_dof);
%  
progressbar;
for i=1:num_dof
    progressbar(i/num_dof);
%       
        acc(1:nt,ngw(i))=    (ZaccT(1:nt,i));     
%
end
clear zaccT;
pause(0.5);
progressbar(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
rel_disp=zeros(nt,num_dof);
progressbar;
for i=1:num_dof
    rel_disp(:,i)=d(:,i)-ddisp(:);
end
pause(0.5);
progressbar(1);
%
clear ddisp;
clear dvelox;
%
if(iu==1)
    G_conv=386.;
else
    G_conv=9.81;
end
%
rel_disp=rel_disp*G_conv;
vel=v*G_conv;
%