function sdot = eval_ebike_rhs(t, s, r, p)
% EVAL_EBIKE_RHS - Returns the time derivative of the states, i.e. evaluates
% the right hand side of the explicit ordinary differential equations.
%
% Syntax: sdot = eval_rhs(t, s, r, p)
%
% Inputs:
%   t - Scalar value of time, size 1x1.
%   s - State vector at time t, size mx1 where m is the number of
%       states.
%   r - Input vector at time t, size ox1 were o is the number of
%       inputs.
%   p - Constant parameter vector, size px1 were p is the number of
%       parameters.
% Outputs:
%   sdot - Time derivative of the states at time t, size mx1.

% unpack the states
theta = s(1);  % rotation angle of the wheel [rad]
omega = s(2);  % angular rate of the wheel [rad/s]
i = s(3);  % motor current [A]
x = s(4);  % x distance traveled [m]

% unpack the input
V = r(1);  % voltage applied to motor [V]

% unpack the constant parameters
m = p(1);
R = p(2);
Cr = p(3);
Cd = p(4);
rho = p(5);
A = p(6);
g = p(7);
J = p(8);
bm = p(9);
Kt = p(10);
L = p(11);
Rw = p(12);
X = p(13);
H = p(14);

% calculate the road angle at this x location on the sine curve road
wr = 1/2/pi/X;
slope = H*wr*cos(wr*x);
alpha = atan(slope);

% caculate some intermediate quantities for evaluating the ODEs
v = omega*R;  % v-omega relationship from no-slip tires

% forces acting on e-bicyclist
Fd = 1/2*rho*Cd*A*v^2;  % drag force [N]
Fr = Cr*m*g*cos(alpha);  % rolling resistant [N]
Fg = m*g*sin(alpha);  % gravity force [N]

% evaluate the four 1st order explicit ordinary differential equations
thetadot = omega;
omegadot = (Kt*i - bm*omega - Fd*R - Fr*R - Fg*R) / (2*J + m*R^2);
idot = (-Rw*i - Kt*omega + V)/L;
xdot = v*cos(alpha);

% pack the state derivatives in a vector
sdot = [thetadot; omegadot; idot; xdot];

end

% This m file simulates the linear motion of a simple electric bicycle model
% traveling over a 2D sinusoidal road. The sinusoidal road is defined by its
% amplitude, H, and the horizontal distance corresponding to the "period" of
% the road. The states are wheel angle theta, wheel angular rate omega,
% motor current i, and horizontal distance traveled x. The no-slip
% constraint provides that the forward speed v is equal to the angular rate
% time the wheel speed.

% define numerical constants
m = 100;  % mass of bike and rider [kg]
R = 0.3;  % wheel radius [m]
Cr = 0.005;  % coefficient of rolling resistance [unitless]
Cd = 1.1;  % coefficient of drag [unitless]
rho = 1.2;  % density of air [kg/m^3]
A = 0.5;  % frontal area [m^2]
g = 9.81;  % acceleration due to gravity [m/s^2]
J = 0.24;  % wheel moment of inertia [kg m^2]
bm = 0.2;  % hub friction coefficient [N*m*s]
Kt = 1.5;  % motor torque constant [N*m/A]
L = 0.5;  % motor inductance [H]
Rw = 1;  % motor winding resistance [Ohm]
X = 5;  % x distance for one period of sine curve road [m]
H = 3;  % amplitude of sine curve road [m]

% pack constant parameters into a 1d vector
p = [m; R; Cr; Cd; rho; A; g; J; bm; Kt; L; Rw; X; H];

% apply a constant voltage (throttle applies 50 volts)
V = 50;  % voltage [V]

% create an anoymous function to pass in inputs and parameters
rhs = @(t, s) eval_ebike_rhs(t, s, V, p);

% integrate the dynamic equations over 60 seconds with the initinal
% condition at zero
[t, res] = ode45(rhs, linspace(0, 60), [0, 0, 0, 0]);

% extract the states from the integration results
theta = res(:, 1);
omega = res(:, 2);
i = res(:, 3);
x = res(:, 4);

% plot the variables of interest versus time
subplot(6, 1, 1);
plot(t, 180/pi*theta);
ylabel('\theta [deg]');
xlabel('Time [s]');

subplot(6, 1, 2);
plot(t, 180/pi*omega);
ylabel('\omega [deg/s]');
xlabel('Time [s]');

subplot(6, 1, 3);
plot(t, i);
ylabel('i [A]');
xlabel('Time [s]');

subplot(6, 1, 4);
plot(t, x);
ylabel('x [m]');
xlabel('Time [s]');

subplot(6, 1, 5);
plot(t, omega*R);
ylabel('v [m/s]');
xlabel('Time [s]');

subplot(6, 1, 6)
plot(t, V*i)
xlabel('Time [s]');
ylabel('Power [W]')