% M-file name: single_wheelset.m
% M-file type: Script (main) file

% This program simulates the dynamic behavior of a single wheelset rolling on a straight track.
% The wheelset is modeled using six degrees of freedom, namely, lateral and yaw wheelset
% displacements, lateral and yaw wheelset velocities, and left and right rail displacements.

% Initial conditions

% x10: Initial Lateral Displacement of wheelset (m)
% x20: Initial Yaw Displacement of wheelset (rad)
% x30: Initial Lateral Velocity of wheelset (m/sec)
% x40: Initial Yaw Velocity of wheelset (rad/sec)
% x50: Initial Lateral Displacement of left rail (m)
% x60: Initial Lateral Displacement of right rail (m)

x10=0.00635;x20=0.0010;x30=0.10;x40=0.0;x50=0.0000;x60=0.0000;

% Globalizing all the variables

% Values for global variables need to be specified in the main file alone even though they may
% be used in several function files.

% Global variables

% V: Forward velocity of wheelset (m/sec)
% lambda: Wheel conicity
% a: Half of track gage (m)
% r0: Centered rolling radius of the wheel (m)
% yfc: Flange clearance or flange width (m)
% yfctol: Lateral tolerance added to yfc in order to facilitate numerical simulation (m)
% mw: Mass of wheelset (kg)
% Iwz: Yaw principal mass moment of inertia of wheelset (kg-m2)
% Iwy: Pitch principal mass moment of inertia of wheelset (kg-m2)
% krail: Effective lateral stiffness of rail (N/m)
% crail: Effective lateral damping of rail (N/m)
% g: Acceleration due to gravity (m/s2)
% muN: Product of coefficient of friction between wheel and rail (mu) and the normal load on
% the axle (N)
% f11: Lateral creep coefficient (N)

% f12: Lateral/Spin creep coefficient (N-m)
% f22: Spin creep coefficient (N-m2)
% f33: Longitudinal creep coefficient (N)
% kpx: Primary longitudinal stiffness coefficient (N/m)
% cpx: Primary longitudinal damping coefficient (N-s/m)
% kpy: Primary lateral stiffness coefficient (N/m)
% cpy: Primary lateral damping coefficient (N-s/m)
% dp: Half of lateral distance between primary longitudinal springs (m)
% N: Axle load (N)

global V lambda a r0 yfc yfctol mw Iwz Iwy krail crail g muN f11 f12 f22 f33 kpx cpx kpy cpy dp N;

V=20;lambda=0.125;a=0.716;r0=0.3556;yfc=0.0080;yfctol=0.0010;mw=1751;Iwz=761;Iwy=130;...
krail=14.6e7;crail=14.6e4;g=9.81;muN=12000;f11=9430000;f12=1.2e3;
f22=1e3;f33=10230000;kpx=9.12e5;cpx=8376.9;kpy=5.84e5;cpy=9048.2;dp=0.61;N=100000;

% Solving the system of differential equations

% The system of differential equations is integrated from t=0 sec to t=10 sec with the above
% initial conditions. The state space variables to be solved are specified in the function file
% 'equations.m' . The solver automatically chooses a suitable time step for integration.

tic,[t,x]=ode45('equations',[0,50],[x10;x20;x30;x40;x50;x60]);toc

% Plotting time response and phase portraits 
% Plotting lateral and yaw displacements vs. time 
subplot(2,1,1);plot(t,x(:,1),'r')
xlabel('Time')
ylabel('Lateral Displacement')
title('Lateral Displacement vs Time')
grid on 
subplot(2,1,2);plot(t,x(:,2),'r')
xlabel('Time')
ylabel('Yaw Displacement')
title('Yaw Displacement vs Time')
grid on 
figure 
% Plotting left and right rail displacements vs. time 
subplot(2,1,1);plot(t,x(:,5),'r')
xlabel('Time')
ylabel('Left Rail Displacement')
title('Left Rail Displacement vs Time')
grid on
subplot(2,1,2),plot(t,x(:,6),'r')
xlabel('Time')
ylabel('Right Rail Displacement')
title('Right Rail Displacement vs Time')
grid on

figure

% Plotting lateral displacement vs. lateral velocity

plot(x(:,1),x(:,3),'r')
xlabel('Lateral Displacement')
ylabel('Lateral Velocity')
title('Lateral Displacement -Velocity Phase Portrait')
grid on
figure

% Plotting yaw displacement vs. yaw velocity

plot(x(:,2),x(:,4),'r')
xlabel('Yaw Displacement')
ylabel('Yaw Velocity')
title('Yaw Displacement -Velocity Phase Portrait')
grid on