% Skli på halvkule uten friksjon
% TFY4115 Fysikk 2012 Øving 9 Opg. 1c.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear; % Nullstille alle variabler og deres dimensjoner

% Definerer parametrene (dt og thetamarg kan endres før hver kjøring av programmet):
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thetamarg = 0.01; % Vinkelmargin ved start
dt = 0.01;  % Tidsinkrement i sek. Bør være max 0.02
g = 9.81;   % Tyngdens aks i m/s^2
R = 50;     % Kuleradius i m
pihalf = pi/2; % Maks vinkel, 90 grader. Eg. unødvendig når test FN>0

% Initialiseringer:
theta(1) = thetamarg; % Minstevinkel for å starte å gli
t(1) = 0.0;
omega(1) = 0.0; 
FN(1) = g*cos(theta(1));  % Normalkraft/masse
alfa(1) = g*sin(theta(1))/R; % Vinkelakselerasjon

n = 1; nmax = 100/dt;  % nmax for å stoppe utilsiktet uendelig løkke
while (FN(n)>=0 & theta(n)<pihalf & n<nmax)  
  omega(n+1) = omega(n) + alfa(n)*dt; % Ny vinkelfart, bruker alfa fra forrige loop
  theta(n+1) = theta(n) + omega(n)*dt;  % Ny vinkel.
  FN(n+1) = g*cos(theta(n+1)) - omega(n+1)*omega(n+1)*R; % Normalkraft/masse
  alfa(n+1) = g*sin(theta(n+1))/R;  % Vinkelakselerasjon
  t(n+1) = t(n)+dt;
  n = n+1;
end

figure(1); % Figurnummer
% 2 subplot side om side:
subplot(1,2,1); % theta, omega, FN s.f.a. t
plot( t,theta*180/pi, t,omega*180/pi,':', t,FN,'-.'); % Vise vinkel og vinkelfart i grader
title('\theta og \omega som funksjon av tid')
xlabel('t/s'); 
ylim([0 90]); % Skalaangivelse 
handle=legend('\theta/grad', '\omega/(grad\cdot s^{-1})',  'F_N/(N/kg)', 'location','north');
set(handle, 'Box', 'off'); % Ikke boks rundt legend

subplot(1,2,2); % omega, alfa, FN s.f.a. theta 
plot( theta*180/pi,omega*180/pi,':', theta*180/pi,alfa*180/pi,'-',  theta*180/pi, FN,'-.')
title('\omega, \alpha og F_N som funksjon av \theta')
xlabel('\theta/grad');
handle=legend('\omega/(grad\cdot s^{-1})', '\alpha/(grad\cdot s^{-2})', 'F_N/(N/kg)', 'location', 'northeast');
set(handle, 'Box', 'off');

% Rapportere noen verdier til skjermen: 
fprintf(1,'theta(1)=%5.2f dt=%5.3f nmax=%i\n', theta(1)*180/pi, dt, nmax);
fprintf(1,'FN=0 når n=%i t=%5.2f theta=%5.2f omega=%4.2f v=%5.2f \n\n',n, t(n), theta(n)*180/pi, omega(n)*180/pi, omega(n)*R);
