%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               File Name: Numerical_Evaluation_of_Dynamic_Response_Linear_Interpolation_SDOF               Written by: Koral Eren               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%d = Displacement
%v = Velocity
%a = Acceleration

%m = Mass
%k = Stiffness
%X = Damping ratio (must be smaller than 1) (underdamped) (almost all structures)
%Dt = Time step
%P = External force (can be an equation or -mass x ground acceleration for earthquake)
%u0 = Initial displacement ( can be assumed 0 for earthquake)
%udot0 = Initial velocity ( can be assumed 0 for earthquake)

function [d,v,a] = Solve_SDOF_Linear(m,k,X,Dt,P,u0,udot0)

	Wn = sqrt(k/m);
    	Wd = Wn * sqrt(1-X^2);

    	L=length(P);

    	% Recurrence equations
    	A = exp(-X*Wn*Dt)*(X/sqrt(1-X^2)*sin(Wd*Dt)+cos(Wd*Dt));

    	B = exp(-X*Wn*Dt)*((1/Wd)*sin(Wd*Dt));

    	C = (1/k)*(((2*X)/(Wn*Dt))+exp(-X*Wn*Dt)*(((1-2*X^2)/(Wd*Dt)-X/sqrt(1-X^2))*sin(Wd*Dt)-(1+(2*X)/(Wn*Dt))*cos(Wd*Dt)));

    	D = (1/k)*(1-((2*X)/(Wn*Dt))+exp(-X*Wn*Dt)*((2*X^2-1)/(Wd*Dt)*sin(Wd*Dt)+(2*X)/(Wn*Dt)*cos(Wd*Dt)));

    	Aprime = -exp(-X*Wn*Dt)*(Wn/sqrt(1-X^2)*sin(Wd*Dt));

    	Bprime = exp(-X*Wn*Dt)*(cos(Wd*Dt)-X/sqrt(1-X^2)*sin(Wd*Dt));

    	Cprime = (1/k)*(-(1/Dt)+exp(-X*Wn*Dt)*(((Wn/sqrt(1-X^2))+(X/(Dt*sqrt(1-X^2))))*sin(Wd*Dt)+(1/Dt)*cos(Wd*Dt)));

    	Dprime =(1/(k*Dt))*(1-exp(-X*Wn*Dt)*(X/sqrt(1-X^2)*sin(Wd*Dt)+cos(Wd*Dt)));
 
    	d=zeros(L,1);
    	v=zeros(L,1);
    	a=zeros(L,1);
 
    	d(1) = u0;
    	v(1) = udot0;
    	a(1) = 0;
 
    	for i=2:L
        
        	d(i) = A * d(i-1) + B * v(i-1)+ C * P(i-1) + D * P(i);
        	v(i) = Aprime * d(i-1) + Bprime * v(i-1) + Cprime * P(i-1) + Dprime * P(i);
        	a(i) = (P(i)-2*X*sqrt(m*k)*v(i)-k*d(i))/m;
        
    	end

end