How the Function "NewmarkNon" is working and meaning of Keyi and Keyp in this coding?

%----------------------- Mario Paz Example 7.1 ----------------------------
clear;
%--------------------------------------------------------------------------
inpf=input('File name of input file data   ','s');
M=importdata([inpf,'.xlsx']);
f=M.data.Force1;        % Force
t = f(:,1); %Time
F = f(:,2); %Force
Dt = t(2)-t(1); %Time interval %Time interval
m=0.2;   % Mass
k =12.35; %Stiffness
c =0.274; %Damping coefficient
xi =c/(2*sqrt(m*k)); %Damping ratio
omega = sqrt(k/m); %Natural frequency (rad/sec)
Rt = 15; %Forces yielding in tension
Rc = -15; %Forces yielding in compression
%%%%Linear acceleration method (Newmark beta method[Ch.6])
gamma =1/2;
beta = 1/6;
tt= length(t);
%%%Initial calculation
u(1)=0; %Initial condition; Displ.
v(1)=0; %Initial condition; Velocity
a(1)=(F(1)-c*v(1)-k*u(1))/m;
A = m/(beta*Dt)+gamma*c/beta; %A in DFbar = DF + A*v0+B*a0 (Eq. 6.46)
B = m/(2*beta)+Dt*c*((0.5*gamma/beta)-1); %B in DFbar = DF + A*v0+B*a0 (Eq. 6.46)
kbar = k +gamma*c/(beta*Dt)+m/(beta*Dt*Dt); %Eq.6.45
key = 0; %Initial key=0 before Loop over
k_p = k; %Initial stiffness before iteration
%%%Setting up for initial value of Loop over
u0=u;
v0=v;
a0=a;
t =t(1);
%%%% Calculate initial yield points %%%%
ut = Rt/k;
uc = Rc/k;
R(1)=0;
%%%%Iteration for each time step using Newmark beta method
ua =[]; 
va =[]; 
aa=[]; 
ta=[];
for i = 1:(tt-1)
DF=F(i+1)-F(i);
F1 = F(i+1);
[t,u,v,a, kbar, R, keyp, key, Du, k_p] = NewmarkNon( t, DF, Dt, u0, v0, ut, uc, a0, F1, k,c, m, Rt, Rc, R, gamma, beta, key, k_p);
keyi(:,i)=key(:,1)';
keypi(:,i)=keyp(:,1)';
F1i(:,i)=F1(:,1)'; 
Ri(:,i)= R(:,1)';
k_pi(:,i)=k_p(1,:)';
% %%%Creating column of time,displ,velocity and acceleration
ta = [ta; t];
ua = [ua; u];
va = [va; v];
aa = [aa; a];
%%%New yielding point-Increasing displacement (Eq.6.49)
if v < 0 && v0 >0
ut = max(ua);
uc = ut-(Rt-Rc)/k;
else
ut = ut;
uc = uc;
end
%%%New yielding point-Decreasing displacement (Eq.6.50)
if v > 0 && v0 <0
uc = min(ua);
ut = uc+(Rt-Rc)/k;
else
ut = ut;
uc = uc;
end
%%%Setting up parameters for next iteration
u0 = u;
v0 = v;
a0 = a;
t = t;
R = R;
end
result = [ta,F1i',ua,keypi',keyi',va,Ri',aa,k_pi'];
%%%Plot response
figure (1)
plot (ta, ua);
grid on
xlabel ('t(sec)');
ylabel ('Displacement(in.)');

figure (2)
plot (ua,Ri);
grid on
xlabel ('Displacement(in.)');
ylabel ('Restoring Force (kip)');

function [t,u,v,a, kbar, R, keyp, key, Du, k_p] = NewmarkNon( t, DF, Dt, u0, v0, ut, uc,a0, F1, k, c, m, Rt, Rc, R, gamma, beta, key, k_p)
%%%Algorithm-For each time step:(3),(4),(5),(6),(7)
kbar = k_p +gamma*c/(beta*Dt)+m/(beta*Dt*Dt); %Eq.6.57
A = m/(beta*Dt)+gamma*c/beta; %A in DFbar (Eq. 6.46)
B = m/(2*beta)+Dt*c*((0.5*gamma/beta)-1); %B in DFbar (Eq. 6.46)
DFbar = DF + A*v0+B*a0; % I n cremental effective force (Eq.6.59)
Du = DFbar/kbar; % I n cremental displacement (Eq.6.60)
Dudot = gamma*Du/(beta*Dt)-gamma*v0/beta+ Dt*a0*(1-0.5*gamma/beta); %Incremental velocity(Eq.6.61)
u=u0+Du ; %Displacement at the end of time interval (Eq.6.62)
v=v0+Dudot; %Velocity at the end of time interval (Eq.6.63)
%%%Algorithm-For each time step:(2)-(a)
if u<ut && u>uc
keyp = 0;
elseif u > ut
keyp = 1;
else
keyp = -1;
end
%%%Algorithm-For each time step:(2)-(b) and (2)-(c)
if keyp == 1 && v > 0
key = 1;
elseif keyp == -1 && v < 0
key = -1;
else
key = 0;
end
%%%Algorithm-For each time step:(8) (Eq. 6.65)
if key ==0
if (R+(Du)*k)>=0
R = min(R+(Du)*k, Rt);
else
R = max(R+(Du)*k, Rc);
end
elseif key ==1
R = Rt;
else
R = Rc;
end
%%%Algorithm-For each time step:(3) (Eq.6.58)
if key == 0
k_p = k;
else
k_p = 0;
end
%%%Algorithm-For each time step:(8) (Eq.6.64)
a=1/m*(F1-c*v-R);
t=t+Dt;
end