clear, clf, hold off

%*************************************************************
% INPUT PARAMATERS FOR THE MODEL
% Car design stuff
% frontal area of the car
A_car=0.002 %meters squared (a big car)
%weight = 5 oz, an oz is 28.35 grams -- the mass to put into all of the
%equations
M=(5 * 28.35)/1000
% Drag coeffieicnt for the car
Cd=0.2;
% Frictional coefficent axle to wheel ID
% http://www.wps.on.ca/technical/tables/uhmw-polyethylene_tables.htm
mu=0.15
% Wheel stuff, ri=inside diameter, ro=outside diameter, I=mass moment of
% inertia
ri=0.003
ro=0.015
I=.000239
% Angle of track 20 degrees
theta=20*(pi/180);

%********************************************************
% The we live on the surface of the earth parameters
rho_air=1.22;
g=9.81;
% h matters it is the step size
h=0.01;
%*********************************************************

% Initialize loop paramaters
t=0;
n=0;
v=0;
% Starting loop values
t_rec(1)=t;
v_rec(1)=v;

% 'da loop
while t<=20
n=n+1;
fd=(-Cd*.5*rho_air*A_car*v*v)
ff=mu*M*g*(ri/ro)
v=v+h*((fd+ff)/M + g*sin(theta));
t=t+h;
v_rec(n+1)=v;
t_rec(n+1)=t;
end

%*****************************************************************
% Plot the velocity versus time
subplot(2,1,1); plot (t_rec,v_rec);
xlabel('time (s)')
ylabel('velocity (m/s)')

%*****************************************************************
% Plot the position versus time
xzero =0
pos = 0.2*(cumsum(v_rec) -(v_rec + v_rec(1))/2 ) + xzero;
subplot(2,1,2); plot (t_rec,pos);
xlabel('time (s)')
ylabel('Position (m)')