disp('  ');
disp(' mdof_arbit_force_newmark.m   ver 1.4  December 21, 2011');
disp('  ');
disp(' by Tom Irvine   Email: tomirvine@aol.com');
disp('  ');
disp(' Reference:  Rao V. Dukkipati, Vehicle Dynamics');
disp('  ');
disp(' This script uses the Newmark-Beta method to solve the following '); 
disp(' equation of motion:   M (d^2x/dt^2) + C dx/dt + K x = F ');
%
disp(' ');
%
clear DI;
clear VI;
clear w;
clear t;
clear M;
clear C;
clear K;
clear accel;
clear length;
clear tt;
clear ff;
clear FI;
clear tint;
clear fint;
clear x;
clear y;
clear accel;
%
disp(' ');
disp(' Enter the units system ');
disp(' 1=English  2=metric ');
iu=input(' ');
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
if(iu==1)
     disp(' Select input mass unit ');
     disp('  1=lbm  2=lbf sec^2/in  ');
     imu=input(' ');
else
    disp(' mass unit = kg ');
end
%
disp(' ');
if(iu==1)
    disp(' damping unit = lbf sec/in ');
else
    disp(' damping unit = N sec/m ');
end
%
disp(' ');
if(iu==1)
    disp(' stiffness unit = lbf/in ');
else
    disp(' stiffness unit = N/m ');
end
%
disp(' ');
disp(' Select file input method ');
disp('   1=file preloaded into Matlab ');
disp('   2=Excel file ');
file_choice = input('');
%
disp(' ');
disp(' Mass Matrix ');
%
if(file_choice==1)
        M = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        M = xlsread(xfile);
%         
end
%
if(iu==1 && imu==1)
  M=M/386;
end
%
disp(' ');
disp(' Damping Matrix ');
%
if(file_choice==1)
        C = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        C = xlsread(xfile);
%         
end
%
disp(' ');
disp(' Stiffness Matrix ');
%
if(file_choice==1)
        K = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        K = xlsread(xfile);
%         
end
% 
disp(' ');
disp(' Enter the initial velocity vector ');
VI=input(' ');
disp(' ');
disp(' Enter the initial displacement vector ');
DI=input(' ');
disp(' ');
% 
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(K,M,1);
%
Tmin=1/(max(fn));
Tmax=1/(min(fn));
%
disp(' ');
disp(' Enter duration(sec)');
dur=input(' ');
%
srr=20*max(fn);
disp(' ');
disp(' Enter sample rate (samples/sec)');
out1=sprintf(' (recommend %8.4g)',srr);
disp(out1);
sr=input(' ');
%
dt=1/sr;
%
NT=round(dur/dt);
%
t=linspace(0,dur,NT);
%
ndof=max(size(M));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[FFI,force_dof,nff]=ODE_force_input(iu,ndof,NT,t,file_choice);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[U,Ud,Udd]=Newmark_force(DI,VI,dt,NT,ndof,M,C,K,FFI,force_dof);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
x=U';
y=Ud';
accel=Udd';
num=ndof;
%
mdof_plot(t,x,y,accel,num,iu)