% Code for Problem 2 in HW 4:

global a b track_type

m=10;          % mass of the car  [kg]
g=9.8;         % gravity  [m/s^2]

dt=0.01;       % step size in time  [s]

% -------- Set up the parameters determining the shape of the track  -------------
% ----- (see program  hw4_p2_track.m  for more detail):  -------
a=3;           % a parameter that determines the height and length of the track
               %  (it actually equals the length of the track, which is assumed to start at x=0)
b=0.1;         % the parameter that can control the height of the track independently of  a;
               %  the height is simply proportional to  b
track_type=1;  % 1 => steeper up, gentler sloping down
               % 2 => vice versa
xtrack=[0:a/100:a];
[ytrack junk1 junk2 junk3]=hw4_p2_track(xtrack);
               % set up the shape of the track; for this we don't need yp, ypp, and Q
      
% -------------  Prepare the initial conditions of the car:  -------------------
V0=1;         % initial velocity (along the track)  [m/s]
X=0;            % car's initial location on the track
[y yp ypp Q]=hw4_p2_track(X);  
Y=y;            % (Note: Matlab is case-sensitive.)
VX=V0/Q;        % car's initial x-velocity  [m/s]
auxRHS=(g+ypp*VX^2)/Q;
                % expression appearing on the rhs's of (10b) and (17)
R=m*auxRHS;     % normal reaction force (at X=0)

% -------------  Plot the track and the car on it at t=0:  ----------------------
figure(53); 
plot(xtrack,ytrack,'k')
hold on
% Axis limits along x- and y-axes:
axlim_xmin=-0.1;
axlim_xmax=a*2;
axlim_ymin=-0.02*max(ytrack);
axlim_ymax=2*max(ytrack);
plot(X,Y,'o','markersize',10,'markerfacecolor','g')
axis([axlim_xmin axlim_xmax axlim_ymin axlim_ymax])
title(['R= ' num2str(R) ',   X=' num2str(X) ',   Y=' num2str(Y)], 'fontsize', 16)
text(0.7*axlim_xmax, 0.85*axlim_ymax, ['Track type ' int2str(track_type)],'fontsize',14)
hold off
disp('press "Enter" to continue')
pause

% -------  Prepare to record the normal reaction force as a function of X:  ------
Xrecord=X;
Rrecord=R;


% -----  Compute and plot the position of the car until it reaches the ground  -----
while Y > -1000*eps & R > 0 & X < a   
    
    % Solve the equations of motion on the track by the modified Euler method:
    Xb = X + dt*VX;
    RHS = -auxRHS*yp/Q;      % rhs of (9b)
    VXb = VX+dt*RHS;
    [y yp ypp Q] = hw4_p2_track(Xb); 
    auxRHS = (g+ypp*VXb^2)/Q;
    X = X+(dt/2)*(VX+VXb);
    VX = VX+(dt/2)*(RHS-auxRHS*yp/Q);
    [y yp ypp Q] = hw4_p2_track(X); 
    Y=y;
    auxRHS = (g+ypp*VX^2)/Q;
    R=m*auxRHS;          % normal reaction force [N] (Newtons)
    
    % Record the normal reaction force:
    Xrecord=[Xrecord X];
    Rrecord=[Rrecord R];
    
    % Plot the car and the track:
    %  (REMINDER: When the car is on the track, it's green; off the track - it's red.)
    figure(53); 
	plot(xtrack,ytrack,'k')
	hold on
	plot(X,Y,'o','markersize',10,'markerfacecolor','g')
	axis([axlim_xmin axlim_xmax axlim_ymin axlim_ymax])
	title(['R= ' num2str(R) ',   X=' num2str(X) ',   Y=' num2str(Y)], 'fontsize', 16)
	text(0.7*axlim_xmax, 0.85*axlim_ymax,['Track type ' int2str(track_type)],'fontsize',14)
	hold off
	pause(0.1) 

end

% -----------------  Plot the car when it is in the air: ----------------
if R <= 0
    X0=X;
    Y0=Y;
    VX0=VX;
    VY0=VX*yp;
    t=0;
    while Y > -1000*eps
        t=t+dt;
        X=X0+VX0*t;
        Y=Y0+VY0*t-g*t^2/2;
        [aux_ytrack junk1 junk2 junk3]=hw4_p2_track(X);
        
        figure(53); 
		plot(xtrack,ytrack,'k')
		hold on
		plot(X,Y,'o','markersize',10,'markerfacecolor','r')
		axis([axlim_xmin axlim_xmax axlim_ymin axlim_ymax])
		title(['Y_{CAR} - Y_{TRACK} = ' num2str(Y-aux_ytrack)], 'fontsize', 13)
		text(0.7*axlim_xmax, 0.85*axlim_ymax,['Track type ' int2str(track_type)],'fontsize',14)		
        hold off
		pause(0.1) 
    end
end

% Plot the reaction force:
figure(531); 
plot(xtrack,0.8*ytrack*max(abs(Rrecord))/max(ytrack),':',Xrecord, Rrecord,'r')
axis([-0.3 4 -0.1 1.1*max(Rrecord)])
legend({'track (height not to scale)','reaction force,  Newtons'}, 'fontsize', 12)
xlabel('X-coordinate,  meters','fontsize',14)