%loop_the_loop
%
% Note to loop the loop (i.e. maintain contact with loop at all times)
% h >= 5*r/2;
%
% LAST UPDATED by Andy French Jan 2019

function loop_the_loop

%Number of trajectory points and fontsize, linewidth
N = 500;
fsize = 18;
linewidth = 2;

%Strength of gravity and radius of circular loop
g = 9.81;
r = 1;

%Initial car height above base of loop, to land at base of loop
%following departure from loop and projectile motion
h = 7*r/4;

%Angle the particle departs from loop
theta = acos( (2/3)*(1-h/r) );
phi = theta - pi/2;

%Velocity at point of departure
v = sqrt( 2*g*(h-r)/3 );

%Time of flight in trajactory phase
T = r*cos(phi)/(v*sin(phi));

%Time array /s for projectile phase
t = linspace(0,T,N);

%Projectile phase x,y coordinates
x = r*cos(phi) -v*t*sin(phi);
y = r + r*sin(phi)+v*t*cos(phi) - 0.5*g*(t.^2);

%Loop x,y coordinates
a = linspace(0,2*pi,N);
xc = r*cos(a);
yc = r + r*sin(a);

%Plot loop and trajectory
plot(xc,yc,'r','linewidth',linewidth);
axis equal;
hold on;
xlim([-r,r]);
ylim([0,2*r]);
plot(x,y,'b--','linewidth',linewidth);
xlabel('x','fontsize',fsize);
ylabel('y','fontsize',fsize);
title('Trajectory where h=7r/4','fontsize',fsize);
set(gca,'fontsize',fsize,'Box','on','Xgrid','on','Ygrid','on',...
    'Xminorgrid','on','Yminorgrid','on','GridLineStyle','-');
print(gcf,'loop_the_loop.png','-dpng','-r300');
close(gcf);

%End of code
