deltatr=deltart-phi; 
deltafr=deltarf-phi; 
% Computing longitudinal, lateral, and spin creepages at right wheel tread and flange contact
% patches 
etaxrt=a*(x4/V)+(1-rrt/r0);
etayrt=((x3/V)-x2*(rrt/r0)+rrt*dphi/V)/cos(deltatr);
etasprt=((x4/V)+(phi/r0))*cos(deltatr)+(1/r0)*sin(deltatr);
etaxrf=a*(x4/V)+(1-rrf/r0);
etayrf=((x3/V)-x2*(rrf/r0)+rrf*dphi/V)/cos(deltafr);
etasprf=((x4/V)+(phi/r0))*cos(deltafr)+(1/r0)*sin(deltafr);

% Computing longitudinal, lateral, and spin creepages at left wheel tread contact patch
etaxl=-a*(x4/V)+(1-rl/r0);
etayl=((x3/V)-x2*(rl/r0)+rl*dphi/V)/cos(deltal);
etaspl=((x4/V)+(phi/r0))*cos(deltal)-(1/r0)*sin(deltal);

% Obtaining creep forces and moments at left and right wheel contact patches from function
% 'twopt_creep'

[Fcxrt,Fcyrt,Fczrt,Mcxrt,Mcyrt,Mczrt,Fcxrf,Fcyrf,Fczrf,Mcxrf,Mcyrf,Mczrf,Fcxl,Fcyl,Fczl,Mcxl,Mcyl,Mczl]=twopt_creep(x1,etaxrt,etayrt,etasprt,etaxrf,etayrf,etasprf,etaxl,etayl,etaspl,deltatr,deltafr,deltal);

% Correcting signs of the vertical creep force components when the vehicle has two-point
% right wheel contact.

Fczrt=-Fczrt; Fczrf=-Fczrf; Fczl=-Fczl;
Mcyrt=-Mcyrt; Mcyrf=-Mcyrf; Mcyl=-Mcyl;

%Reversing the signs of lateral creep force values to fit in the 'twopt_normal' function under
% "going right" condition

yfc=-yfc;Fcyrt=-Fcyrt;Fcyrf=-Fcyrf;Fcyl=-Fcyl;

% Reversing the signs of x1, x2, x3, and x4 to allow for mirror effect of right and left
% traverse

x1=-x1;x2=-x2;x3=-x3;x4=-x4;

% Obtaining normal forces at left and right wheel contact patches from function
% 'twopt_normal'

[Fnyrt,Fnzrt,Fnyrf,Fnzrf,Fnyl,Fnzl]=twopt_normal(x1,x2,x3,x4,deltatr,deltafr,deltal,rrt,rrf,rl,Fcxrt,Fcyrt,Fczrt,Fcxrf,Fcyrf,Fczrf,Fcxl,Fcyl,Fczl,Mcyrt,Mcyrf,Mcyl);

% Correcting the signs of the lateral normal forces to adjust for the mirror effect and
% restoring the signs of the lateral creep forces

Fnyrt=-Fnyrt; Fnyrf=-Fnyrf; Fnyl=-Fnyl; Fcyrt=-Fcyrt; Fcyrf=-Fcyrf; Fcyl=-Fcyl;
% Restoring the signs of x1, x2, x3, and x4
x1=-x1;x2=-x2;x3=-x3;x4=-x4;
% The two-point equations in state space form for right wheel two-point contact
xdot1=x3;
xdot2=x4; 
xdot3=1/mw*(Fcyrt+Fcyrf+Fcyl+Fnyrt+Fnyrf+Fnyl+Fsuspyw-mw*g*lambda*x1/a);
xdot4=1/Iwz*((-Iwy*V/r0*dphi)+a*(Fcxrt+Fcxrf-Fcxl)+x2*((a-rrt*tan(deltatr))*(Fcyrt+Fnyrt)+(a-rrf*tan(deltafr))*(Fcyrf+Fnyrf)-(a-rl*tan(deltal))*(Fcyl+Fnyl))+Mczrt+Mczrf+Msuspzw);
xdot5=(1/crail)*(-Fnyl-Fcyl-krail*x5);
xdot6=(1/crail)*(-Fnyrt-Fnyrf-Fcyrt-Fcyrf-krail*x6);


% Single-point contact equations -Right wheel tread contact and Left wheel flange contact 
elseif (x1-x5)>(yfc+yfctol) 
% Assigning variable names to rolling radii at left and right wheel contact patches
rl=rlf; 
rr=rrt;
% Computing tangent slopes at left and right wheel contact patches
lambdal=tan(deltalf);
lambdar=tan(deltart);
% Computing wheelset roll angle and rate of change of wheelset roll angle
phi=(rl-rr)/(2*a);
dphi=(lambdal+lambdar)*(x3/(2*a)); 
% Computing effective contact angles at left and right wheel contact patches by
% compensating for wheelset roll angle
deltal=deltalf+phi;
deltar=deltart-phi;
% Computing longitudinal, lateral, and spin creepages at left wheel flange contact patch
etaxl=-a*(x4/V)+(1-rl/r0);
etayl=((x3/V)-x2*(rl/r0)+rl*dphi/V)/cos(deltal);
etaspl=((x4/V)+(phi/r0))*cos(deltal)-(1/r0)*sin(deltal);

% Computing longitudinal, lateral, and spin creepages at right wheel tread contact patch
etaxr=a*(x4/V)+(1-rr/r0);
etayr=((x3/V)-x2*(rr/r0)+rr*dphi/V)/cos(deltar);
etaspr=((x4/V)+(phi/r0))*cos(deltar)+(1/r0)*sin(deltar);

% Obtaining creep forces and moments at left and right wheel contact patches from function
% 'onept_creep'

[Fcxl,Fcyl,Fczl,Mcxl,Mcyl,Mczl,Fcxr,Fcyr,Fczr,Mcxr,Mcyr,Mczr]=onept_creep(etaxl,etayl,etaspl,etaxr,etayr,etaspr,deltal,deltar);

% Obtaining normal forces at left and right wheel contact patches from function
% 'onept_normal'

[Fnyl,Fnzl,Fnyr,Fnzr]=onept_normal(x2,x4,deltal,deltar,rl,rr,Fcxl,Fcxr,Fcyl,Fcyr,Fczl,Fczr,Mcyl,Mcyr);

% The single-point equations in state space form 
xdot1=x3;
xdot2=x4;
xdot3=1/mw*(Fcyl+Fcyr+Fnyl+Fnyr+Fsuspyw-mw*g*lambda*x1/a);
xdot4=1/Iwz*((-Iwy*V/r0*dphi)-a*(Fcxl-Fcxr)-x2*((a-rl*tan(deltal))*(Fcyl+Fnyl)-(a-rr*tan(deltar))*(Fcyr+Fnyr))+Mczl+Mczr+Msuspzw);
xdot5=(1/crail)*(-Fnyl-Fcyl-krail*x5);
xdot6=(1/crail)*(-Fnyr-Fcyr-krail*x6);


% Two point contact equations -Right wheel tread contact and Left wheel two-point 
% contact 
elseif (x1-x5)>yfc 
% Assigning variable name to rolling radius at right wheel contact patch 
rr=rrt; 
% Computing tangent slopes at left and right wheel contact patches 
lambdar=tan(deltart);
lambdalt=tan(deltalt);
lambdalf=tan(deltalf);

% Computing wheelset roll angle and rate of change of wheelset roll angle 
phi=(rlt-rr)/(2*a);
dphi=(lambdalf+lambdar)*(x3/(2*a));

% Computing effective contact angles at left and right wheel contact patches by
% compensating for wheelset roll angle

deltatl=deltalt+phi;
deltafl=deltalf+phi;
deltar=deltart-phi;

% Computing longitudinal, lateral, and spin creepages at right wheel tread contact patch

etaxr=a*(x4/V)+(1-rr/r0);
etayr=((x3/V)-x2*(rr/r0)+rr*dphi/V)/cos(deltar);
etaspr=((x4/V)+(phi/r0))*cos(deltar)+(1/r0)*sin(deltar);

% Computing longitudinal, lateral, and spin creepages at left wheel tread and flange contact
% patches

etaxlt=-a*(x4/V)+(1-rlt/r0);
etaylt=((x3/V)-x2*(rlt/r0)+rlt*dphi/V)/cos(deltatl);
etasplt=((x4/V)+(phi/r0))*cos(deltatl)-(1/r0)*sin(deltatl);
etaxlf=-a*(x4/V)+(1-rlf/r0);
etaylf=((x3/V)-x2*(rlf/r0)+rlf*dphi/V)/cos(deltafl);
etasplf=((x4/V)+(phi/r0))*cos(deltafl)-(1/r0)*sin(deltafl);

% Obtaining creep forces and moments at left and right wheel contact patches from function
% 'twopt_creep'

[Fcxlt,Fcylt,Fczlt,Mcxlt,Mcylt,Mczlt,Fcxlf,Fcylf,Fczlf,Mcxlf,Mcylf,Mczlf,Fcxr,Fcyr,Fczr,Mcxr,Mcyr,Mczr]=twopt_creep(x1,etaxlt,etaylt,etasplt,etaxlf,etaylf,etasplf,etaxr,etayr,etaspr,deltatl,deltafl,deltar);

% Obtaining normal forces at left and right wheel contact patches from function
% 'twopt_normal'

[Fnylt,Fnzlt,Fnylf,Fnzlf,Fnyr,Fnzr]=twopt_normal(x1,x2,x3,x4,deltatl,deltafl,deltar,rlt,rlf,rr,Fcxlt,Fcylt,Fczlt,Fcxlf,Fcylf,Fczlf,Fcxr,Fcyr, Fczr,Mcylt,Mcylf,Mcyr);

% The two-point equations in state space form for left wheel two-point contact
xdot1=x3;
xdot2=x4;
xdot3=1/mw*(Fcylt+Fcylf+Fcyr+Fnylt+Fnylf+Fnyr+Fsuspyw-mw*g*lambda*x1/a);
xdot4=1/Iwz*((-Iwy*V/r0*dphi)-a*(Fcxlt+Fcxlf-Fcxr)-x2*((a-rlt*tan(deltatl))*(Fcylt+Fnylt)+(a-rlf*tan(deltafl))*(Fcylf+Fnylf)-(a-rr*tan(deltar))*(Fcyr+Fnyr))+Mczlt+Mczlf+Msuspzw);
xdot5=(1/crail)*(-Fnylt-Fnylf-Fcylt-Fcylf-krail*x5);
xdot6=(1/crail)*(-Fnyr-Fcyr-krail*x6);
% Single-point contact equations -Right wheel tread contact and Left wheel tread contact 
else 
% Assigning variable names to rolling radii at left and right wheel contact patches
rl=rlt;
rr=rrt;

% Computing tangent slopes at left and right wheel contact patches 
lambdal=tan(deltalt);
lambdar=tan(deltart);
% Computing wheelset roll angle and rate of change of wheelset roll angle
phi=(rl-rr)/(2*a);
dphi=(lambdal+lambdar)*(x3/(2*a)); 
% Computing effective contact angles at left and right wheel contact patches by
% compensating for wheelset roll angle
deltal=deltalt+phi;
deltar=deltart-phi;
% Computing longitudinal, lateral, and spin creepages at left wheel tread contact patch
etaxl=-a*(x4/V)+(1-rl/r0);
etayl=((x3/V)-x2*(rl/r0)+rl*dphi/V)/cos(deltal);
etaspl=((x4/V)+(phi/r0))*cos(deltal)-(1/r0)*sin(deltal);
% Computing longitudinal, lateral, and spin creepages at right wheel tread contact patch
etaxr=a*(x4/V)+(1-rr/r0);
etayr=((x3/V)-x2*(rr/r0)+rr*dphi/V)/cos(deltar);
etaspr=((x4/V)+(phi/r0))*cos(deltar)+(1/r0)*sin(deltar);
% Obtaining creep forces and moments at left and right wheel contact patches from function
% 'onept_creep'
[Fcxl,Fcyl,Fczl,Mcxl,Mcyl,Mczl,Fcxr,Fcyr,Fczr,Mcxr,Mcyr,Mczr]=onept_creep(etaxl,etayl,etaspl,etaxr,etayr,etaspr,deltal,deltar);

% Obtaining normal forces at left and right wheel contact patches from function
% 'onept_normal'

[Fnyl,Fnzl,Fnyr,Fnzr]=onept_normal(x2,x4,deltal,deltar,rl,rr,Fcxl,Fcxr,Fcyl,Fcyr,Fczl,Fczr,Mcyl,Mcyr);

% The single point equations in state space form 
xdot1=x3;
xdot2=x4;
xdot3=1/mw*(Fcyl+Fcyr+Fnyl+Fnyr+Fsuspyw-mw*g*lambda*x1/a);
xdot4=1/Iwz*((-Iwy*V/r0*dphi)-a*(Fcxl-Fcxr)-x2*((a-rl*tan(deltal))*(Fcyl+Fnyl)-(a-rr*tan(deltar))*(Fcyr+Fnyr))+Mczl+Mczr+Msuspzw);
xdot5=(1/crail)*(-Fnyl-Fcyl-krail*x5);
xdot6=(1/crail)*(-Fnyr-Fcyr-krail*x6);

end