function dy = bikeacc(t, y)
% acc: Computes accelerations including drag (Cd), rolling (Cr) and 
% and gravity
% y contains, [x z xdot zdot energy] and dy is the time derivative
%
global Grav Rho Mass Cd Cr Area Slope As Bs Lambda Pr Fr_Max TLength % Define these globals so that they will be accessible in the acceleration routine

dy = zeros(5,1);    % a column vector
%
% get the road slope for our current x position.
theta = atan(Slope + As*cos(2*pi*y(1)/Lambda)*2*pi/Lambda - Bs*sin(2*pi*y(1)/Lambda)*2*pi/Lambda);
%
% Compute forces: Assuming bike remains on road surface
% Drag:
vmag = sqrt(y(3)^2+y(4)^2); thetad = atan2(y(4),y(3));
dxd = -Rho*Cd*vmag^2*Area*cos(thetad)/(2*Mass);
dzd = -Rho*Cd*vmag^2*Area*sin(thetad)/(2*Mass);
dxr = +Grav*Cr*cos(theta) ;
dzr = +Grav*Cr*sin(theta) ;
%
% Get rider force
Fr_mag = Fr_Max ;
if ( vmag > 0 ) 
    Fr_mag = Pr/vmag ;
end
if ( Fr_mag > Fr_Max ) Fr_mag = Fr_Max ; end
height = Slope*y(1) + As*sin(2*pi*y(1)/Lambda) + Bs*cos(2*pi*y(1)/Lambda);
%if ( y(2) > height + 0.1 ) Fr_mag = 0; end
dxp = Fr_mag*cos(theta)/Mass ;
dzp = Fr_mag*sin(theta)/Mass ;
%
% Normal force from surface (Grav is a negative value)
dxn =  Grav*cos(theta)*sin(theta);
dzn = -Grav*cos(theta)^2;
% Get the curvature
dy2dx2 = -As*sin(2*pi*y(1)/Lambda)*(2*pi/Lambda)^2 - Bs*cos(2*pi*y(1)/Lambda)*(2*pi/Lambda)^2;
IRc = (dy2dx2)/(1+theta^2)^(3/2);
dxc = vmag^2*sin(theta)*IRc;
dzc = vmag^2*cos(theta)*IRc;
%fprintf(' %5.2f %9.5f %9.5f %5.2f %6.4f %10.1f \n',vmag, dxc, dzc, dxd, dzd, Rc)
%
% Update the dy velocity acceleration vector
dy(1) = y(3);
dy(2) = y(4);
dy(3) = dxd + dxr + dxp + dxn + dxc ;
dy(4) = dzd + dzr + dzp + dzn + dzc + Grav ;
dy(5) = Fr_mag * vmag;