function axle
% kinematic model of the axle (differential drive)

r = .125;   %radius of each actuated wheel (m)
d = .350;   %track width (m)
m = 25;     %mass of the robot (kg)
K = 1000;   %Slip stiffness of the wheels (N/rad)

omega_left = input('angular velocity of left wheel (rad/s)           : ');
omega_right = input('angular velocity of right wheel (rad/s)          : ');

%Note: c.o.m located at the middle point of the rear track 

% Kinematic steering (non-holonomic approach)
v_left=omega_left*r;                %Hypothesis of pure rolling
v_right=omega_right*r;              %Hypothesis of pure rolling
v_G = (v_left + v_right)/2;         %Velocity of the c.o.m. along the direction of rolling
speed = abs(v_G)*3.6;               %Speed of the robot (km/h)
omega_robot = (v_right-v_left)/d;   %Angular velocity of the robot [rad/s]

if abs(omega_robot) < 1e-7
   disp(['The mobile robot travels straight at the speed of ', num2str(speed), ' km/h (no slip)' ]);  
else if speed < 1e-7
   disp(['The mobile robot is spinning on the spot at the angular speed of ', num2str(abs(omega_robot*180/pi)), ' °/s (no slip)']); 
else    
   R = abs(v_G)/omega_robot;                %Radius of the path (m)
   disp(['The mobile robot travels at the speed of ', num2str(speed), ' km/h' ]);  
   disp(['along a circle whose radius is ',num2str(R), ' meters;']);
   disp(['its yaw rate is ',num2str(abs(omega_robot*180/pi)), ' °/s']);

% Slip evaluation
   F=m*v_G^2/R;                        %Centrifugal force [N]
   delta = F/K*180/pi;                 %Slip angle [°]
   disp(['The value of the slip angle is ',num2str(delta),' degrees.']);
   if abs(delta)>2 
      disp('This value is rather high, the accuracy of the non-holonomic model may be unsatisfactory.')
   end

% Plot of the path
   theta=0:pi/100:2*pi;
   hold on;
   plot(R*cos(theta)-R,R*sin(theta),'-b');  % Non-holonomic path
   plot([-d/2,d/2],[0,0],'r');
   plot([-d/2,-d/2],[-r,r],'r');
   plot([d/2,d/2],[-r,r],'r');
   plot(0,.1,'^r');
   axis([-2*abs(R)-.1, 2*abs(R)+.1, -2*abs(R)-.1, 2*abs(R)+.1]);
   axis square;
end   

end