function Result=Central_Difference(Elements,Material,Support,Free,M,C,f,fs)
%Input
% Elements: a structure containing Elements{i}.DOFs and Elements{i}.Material
% where Elements{i}.DOFs=[j k] means element i connect DOF j with k
% and Elements{i}.Material=m assign material m to element i
% Material: a structure containing material properties for bilinear springs
% where Material{m}.k1 is Spring stiffness
% Material{m}.x1 is Spring deformation beyond which the stiffness decreases
% Material{m}.k2 is Reduced stiffness
% Support: a vector of support (Fixed) DOFs of size (nSupport,1)
% Free: a vector of free DOFs of size (nFree,1)
% M:mass matrix (nFree*nFree)
% C:damping matrix (nFree*nFree)
% f:external force matrix(nFree,N)
% fs: sampling frequency
% where N is the length of data points of dynamic force
% Output:
% Result: is a structure consist of
% Result.Displacement: Displacement (nFree*N)
% Result.Velocity: Velocity (nFree*N)
% Result.Acceleration: Acceleration (nFree*N)
%-----------------------------------------
[n, N]=size(f);
dt=1/fs; %Sampling rate
x(:,1)=zeros(n,1); %initial conditions
v(:,1)=zeros(n,1); %initial conditions
[F_i, ~]=Element_Analysis(Elements,Material,Support,Free,x(:,1));
a(:,1)=M\(f(:,1)-C*v(:,1)-F_i); %initial conditions
%Initial Calculations
x_i_1=x(:,1)-v(:,1)*dt+a(:,1)*dt^2/2;
K_hat=M/dt^2+C/(2*dt);
aa=M/dt^2-C/(2*dt);
bb=-2*M/dt^2;
for i=1:1:N-1
F_hat=f(:,i)-aa*x_i_1-bb*x(:,i)-F_i;
x(:,i+1)=K_hat\F_hat;
[F_i, ~]=Element_Analysis(Elements,Material,Support,Free,x(:,i+1));
v(:,i+1)=(x(:,i+1)-x_i_1)/(2*dt);
a(:,i+1)=(x(:,i+1)-2*x(:,i)+x_i_1)/dt^2;
x_i_1=x(:,i);
end
Result.Displacement=x;
Result.Velocity=v;
Result.Acceleration=a;
end
