%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% taperull.m:
% Oving 10 FY1001/TFY4145 Mekanisk fysikk 2012:
% Rullende og slurende taperull paa kvartsirkel.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nullstiller alle variable
clear;      
% Leser eksperimentelle data fra filen taperull.txt som inneholder
% 3 kolonner: t, x, y i enheter hhv sekund, cm, cm.
% Det er 2 linjer med overskrift, derfor 'headerlines',2
[tideksp,xeksp,yeksp] = textread('taperull.txt','%f %f %f','headerlines',2);
% Funksjonen length(tabell) returnerer antall elementer i tabell
lengde = length(tideksp);
% Lager tabeller med eksperimentelle hastighetsverdier.
% Starter med aa lage "tomme tabeller".
Vxeksp = zeros(1,lengde);
Vyeksp = zeros(1,lengde);
for i = 2:lengde-1, % For forste og siste element beregnes v_exp her ikke.
	Vxeksp(i) = (xeksp(i+1) - xeksp(i-1))/(tideksp(i+1) - tideksp(i-1));
    Vyeksp(i) = (yeksp(i+1) - yeksp(i-1))/(tideksp(i+1) - tideksp(i-1));
end
% Lager v_exp for siste element for aa unngaa skjemmende nullverdi i plott
Vxeksp(lengde) = (xeksp(lengde) - xeksp(lengde-1))/(tideksp(lengde) - tideksp(lengde-1));
Vyeksp(lengde) = (yeksp(lengde) - yeksp(lengde-1))/(tideksp(lengde) - tideksp(lengde-1));
% Lager tabeller med eksperimentelle verdier for phi og r+R:
phieksp = zeros(1,lengde);
reksp = zeros(1,lengde);
for i = 1:lengde,
    phieksp(i)=atan(xeksp(i)/yeksp(i));
    reksp(i)=sqrt(xeksp(i)*xeksp(i)+yeksp(i)*yeksp(i));
end

% Deretter numeriske beregninger basert paa teori
% Felles parametre (for alle legemer studert eksperimentelt)
g=9.81;      % tyngdens akselerasjon
phimax=pi/2; % maksimal vinkel
dt=0.01;     % tidsintervall, varieres etter behov
R=0.79;      % sirkelradius 79 cm
% Parametre for: taperull
mys=0.4;    % statisk friksjonskoeffisient, kan justeres
myk=0.4;    % kinetisk friksjonskoeffisient, kan justeres
r=0.0375;     % radius taperull
c=5/9;       % legemets treghetsmoment pr masseenhet er i0=c*r*r
% Diverse analytiske uttrykk og sammenhenger. 
% Sammenhold med ovingsteksten.
% Krefter er pr masseenhet.
% N2 langs banen: g*sin(phi) - f = A
% Normalkraft: n = g*cos(phi)-V*V/(r+R)
% Friksjonskraft, ren rulling: f = c*A
% Friksjonskraft, sluring: f = myk*n

% Initialiseringer:
phi(1)   = 0.045;   % Startvinkel for taperull
tid(1)   = 0.0;     % Starttid
V(1)     = 0.02;     % Starthastighet (CM)
OMEGA(1) = V(1)/(r+R);     % Startvinkelhastighet (CM om origo)
n(1)     = g*cos(phi(1));  % Normalkraft ved start
A(1)     = g*sin(phi(1))/(c+1); % Akselerasjon ved start
ALPHA(1) = A(1)/(r+R);  % Vinkelaks. om origo ved start
%Kartesiske komponenter av posisjon og hastighet:
x(1) = (r+R)*sin(phi(1));
y(1) = (r+R)*cos(phi(1));
Vx(1) = V(1)*cos(phi(1));
Vy(1) = -V(1)*sin(phi(1));

j=1;   % Tabellindeks
jmax=20/dt;  % For forsikring mot uendelig loop
jrulling=1;  % Startverdi for jrulling
jkontakt=1;  % Startverdi for jkontakt
while (n(j) > 0 & phi(j) < phimax & j < jmax)
  V(j+1) = V(j) + A(j)*dt;   % Ny V
  OMEGA(j+1) = V(j+1)/(r+R); % Ny OMEGA
  phi(j+1) = phi(j) + OMEGA(j)*dt; % Ny phi. Osv nedover...
  x(j+1) = (r+R)*sin(phi(j+1));
  y(j+1) = (r+R)*cos(phi(j+1));
  Vx(j+1) = V(j+1)*cos(phi(j+1));
  Vy(j+1) = -V(j+1)*sin(phi(j+1));
  n(j+1) = g*cos(phi(j+1)) - V(j+1)*V(j+1)/(r+R);
  fmax(j+1) = mys*n(j+1); % Maksimal friksjonskraft
  A(j+1) = g*sin(phi(j+1))/(c+1); % Ny akselerasjon, rulling antatt
  ALPHA(j+1) = A(j+1)/(r+R);   % Ny vinkelaks. om origo, rulling antatt
  f(j+1) = c*A(j+1);   % Friksjonskraft, rulling antatt
  % Sjekker om beregnet f, med rulling antatt, har blitt storre
  % enn fmax. Hvis ja: Sett f = myk*n og rett opp A og ALPHA.
  % Hvis nei: Da HAR vi ren rulling fortsatt, og vi oker jrulling
  % med 1 og gaar videre til neste tidsverdi
  if f(j+1) > fmax(j+1)
      f(j+1) = myk*n(j+1);  % Friksjonskraft ved sluring
      A(j+1) = g*sin(phi(j+1)) - f(j+1); % Akselerasjon ved sluring
      ALPHA(j+1) = A(j+1)/(r+R);
  else
      jrulling=j+1; % Sparer paa siste j-verdi der vi hadde ren rulling
  end;
  jkontakt=j+1; % Sparer paa siste j-verdi der vi hadde kontakt
  tid(j+1) = dt*j; % Neste tidsverdi
  j = j+1; % Neste tabellindeks-verdi
end

% Naar vi kommer hit, betyr det at vi har naadd en vinkel som gir
% normalkraft n < 0. Sjekken paa om phi < phimax og j < jmax er for
% helt sikkert avslutte loopen selv om noe gaar galt.

% Herfra og ut handler det om aa visualisere resultatene.
% Her er det valgt aa lage 1 figur med 2x2=4 plott.
% Studer syntaksen og laer av det. Og endre etter onske!

figure(1);

subplot(2,2,1);
plot(tid,phi*180/pi,tid,OMEGA*180/pi,tid,n,tideksp,phieksp*180/pi);
title('\phi, \Omega og N som funksjon av tid')
xlabel('t (s)');
%matlab-spesifikt (neste 2 linjer)
handle=legend('\phi (grader)','\Omega (grader/s)','N (N/kg)','location','northwest');
set(handle, 'Box', 'off');
grid on;

subplot(2,2,2);
plot(phi*180/pi,OMEGA*180/pi,phi*180/pi,ALPHA*180/pi,phi*180/pi,n)
title('\Omega, \alpha og N som funksjon av \phi')
xlabel('\phi (grader)');
%matlab-spesifikt (neste 2 linjer)
handle=legend('\Omega (grader/s)','\alpha (grader/s^2)','N (N/kg)','location','northwest');
set(handle, 'Box', 'off');
grid on;
label=strcat('Slurer fra \phi=',num2str(phi(jrulling)*180/pi));
text(10,-20,label);
label=strcat('Mister kontakten ved \phi=',num2str(phi(jkontakt)*180/pi));
text(10,-50,label);
% Fjerne % i de neste linjene hvis utskrift av verdier ved
% rulleslutt og kontaktslutt onskes
% phi(jrulling)*180/pi,OMEGA(jrulling)*180/pi,V(jrulling)
% phi(jkontakt)*180/pi,OMEGA(jkontakt)*180/pi,V(jkontakt)

subplot(2,2,3);
% Plotter hastigheter i cm/s som funksjon av tiden,
% baade teoretisk og eksperimentelt (begge numerisk beregnet).
plot(tid,Vx*100,tid,Vy*100,tideksp,Vxeksp,tideksp,Vyeksp)
title('V_x og V_y som funksjon av t')
xlabel('t (s)');
handle=legend('V_x (cm/s)','V_y (cm/s)','V_x^{exp}','V_y^{exp}','location','southwest');
set(handle, 'Box', 'off');
grid on;


subplot(2,2,4);
% Plotter hastigheter i cm/s som funksjon av vinkel,
% baade teoretisk og eksperimentelt (begge numerisk beregnet).
plot(phi*180/pi,Vx*100,phi*180/pi,Vy*100,phieksp*180/pi,Vxeksp,phieksp*180/pi,Vyeksp)
title('V_x og V_y som funksjon av \phi')
xlabel('\phi (grader)');
handle=legend('V_x (cm/s)','V_y (cm/s)','V_x^{exp}','V_y^{exp}','location','southwest');
set(handle, 'Box', 'off');
grid on;


% Det var det!
