%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rullende metallsylinder paa kvartsirkel.
% Oving 10 Mekanisk fysikk 2011.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;       % nullstiller alle variable
% Leser data fra filen sylinder.txt som inneholder
% 3 kolonner: t, x, y i hhv sekund, cm, cm
% samt to linjer med overskrift foer dataene kommer
[tideksp,xeksp,yeksp] = textread('sylinder.txt','%f %f %f','headerlines',2);
lengde = length(tideksp);  %Bestem hvor mange linjer med t,x,y data
%Slik lager vi eksperimentell hastighets-vektor:
Vxeksp = zeros(1,lengde); % Lager tomme vektorer forst
Vyeksp = zeros(1,lengde);
for i = 2:lengde-1, %Det forste og siste elementet beregner vi ikke hastighet for.
	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
%Slik lager vi eksp. akselerasjons-vektor:
Axeksp = zeros(1,lengde);
Ayeksp = zeros(1,lengde);
for i = 3:lengde-2, %Bruker ikke det forste og siste elementet i v
	Axeksp(i) = (Vxeksp(i+1) - Vxeksp(i-1))/(tideksp(i+1) - tideksp(i-1));
    Ayeksp(i) = (Vyeksp(i+1) - Vyeksp(i-1))/(tideksp(i+1) - tideksp(i-1));
end
%Slik lager vi eksp. vektor for vinkelen phieksp og avstanden reksp:
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 teorien!
mys=0.25;     % statisk friksjonskoeffisient
myk=0.25;     % kinetisk friksjonskoeffisient
g=9.81;      % tyngdens akselerasjon
phimax=pi/2; % maksimal vinkel
dt=1/120;    % tidsintervall (tilsvarer 120 fps)
%%dt=0.01;   % testverdi for sammenligning med 2010
R=0.79;      % sirkelradius 79 cm
%R=49.961; r=0.039; %testverdier
r=0.019;     % rullende legemes radius (metallsylinder)
c=0.9;       % legemets treghetsmoment pr masseenhet er i0=c*r*r
%c=0;        % testverdi
% Analytiske uttrykk og sammenhenger. Alle krefter er pr masseenhet.
% N2 langs banen: g*sin(phi) - f = A
% N2 for rotasjon om CM: f*r = i0*alpha
% Normalkraft: n = g*cos(phi)-V*V/(r+R)
% Friksjonskraft, ren rulling: f = c*A
% Friksjonskraft, sluring: f = myk*n
% Friksjonskraft, fritt fall: f = 0 (Neglisjer luftmotstand og skrueffekt)

% Initialiseringer:
phi(1)   = 0.027;   % Startvinkel for sylinder i Jon9
%phi(1)=0.01;       % testverdi
tid(1)   = 0.0;     % Starttid
V(1)     = 0.0;     % Starthastighet (CM)
omega(1) = 0.0;     % Startvinkelhastighet (rotasjon om CM)
OMEGA(1) = 0.0;     % 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;      % Vinkelaks. om CM ved start
ALPHA(1) = A(1)/(r+R);  % Vinkelaks. om origo ved start
%Kartesisk:
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));
Ax(1) = A(1)*cos(phi(1));
Ay(1) = -A(1)*sin(phi(1));

j=1;   % Vektorindeks
jmax=20/dt;  % For aa unngaa uendelig loop
jrull=1;  % Startverdi for jrull
jkontakt=1; % Startverdi for jkontakt
while (n(j) >= 0 & phi(j) <= phimax & j <= jmax)
  V(j+1) = V(j) + A(j)*dt;   % Ny hastighet V
  OMEGA(j+1) = V(j+1)/(r+R);     % Ny OMEGA
  phi(j+1) = phi(j) + OMEGA(j)*dt; % Ny vinkel phi
  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); % Ny n
  fmax(j+1) = mys*n(j+1); % Maksimal friksjonskraft
  A(j+1) = g*sin(phi(j+1))/(c+1); % Ny akselerasjon
  ALPHA(j+1) = A(j+1)/(r+R);   % Ny vinkelaks. om origo
  f(j+1) = c*A(j+1);    % Faktisk friksjonskraft, ren rulling
  omega(j+1) = V(j+1)/r;   % Ny omega om CM, hvis ren rulling
  alpha(j+1) = A(j+1)/r;   % Ny alpha om CM, hvis ren rulling
  if f(j+1) >= fmax(j+1) % Sluring! Da maa f,A,omega,alpha,ALPHA endres!!
      f(j+1) = myk*n(j+1);  % Friksjonskraft ved sluring
      A(j+1) = g*sin(phi(j+1)) - f(j+1); % Akselerasjon ved sluring
      omega(j+1) = omega(j) + dt*f(j)/(c*r); % Ny omega (generelt)
      alpha(j+1) = f(j+1)/(c*r);             % Ny alpha (generelt)
      ALPHA(j+1) = A(j+1)/(r+R);    % ....siden A ble endret...
  else
      jrull=j+1; % Sparer paa siste j-verdi der vi hadde ren rulling
  end;
  Ax(j+1) = A(j+1)*cos(phi(j+1));
  Ay(j+1) = -A(j+1)*sin(phi(j+1));
  jkontakt=j+1; % Sparer paa siste j-verdi der vi hadde kontakt
  tid(j+1) = dt*j;
  j = j+1;
end
disp('phislur ='); 
phi(jrull)

figure(1);

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

subplot(3,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');
%ylabel('\Omega  og  \alpha og  N')
%matlab-spesifikt (2 linjer)
handle=legend('\Omega (grad/s)','\alpha (grad/s^2)','N (N/kg)','location','north');
set(handle, 'Box', 'off');
grid on;
%%text(phi(jrull)*180/pi,0.12,'\downarrow','FontSize',24)
label=strcat('Slurer fra \phi=',num2str(phi(jrull)*180/pi));
text(10,-20,label);
label=strcat('Tar av ved \phi=',num2str(phi(jkontakt)*180/pi));
text(10,-50,label);
%phi(jrull)*180/pi,OMEGA(jrull)*180/pi,V(jrull)  % Skriv ut verdier ved slurestart
%phi(jkontakt)*180/pi,OMEGA(jkontakt)*180/pi,V(jkontakt) % Verdier ved kontaktslutt

subplot(3,2,3);
%Plott hastigheter i cm/s
plot(phi*180/pi,Vx*100,phi*180/pi,Vy*100,phieksp*180/pi,Vxeksp,phieksp*180/pi,Vyeksp)
title('Vx og Vy som funksjon av \phi')
xlabel('\phi');
handle=legend('Vx (cm/s)','Vy (cm/s)','Vxeksp','Vyeksp','location','north');
set(handle, 'Box', 'off');
grid on;

subplot(3,2,4);
%Plott akselerasjoner i cm/s*s
plot(phi*180/pi,Ax*100,phi*180/pi,Ay*100,phieksp*180/pi,Axeksp,phieksp*180/pi,Ayeksp)
title('Ax og Ay som funksjon av \phi')
xlabel('\phi');
handle=legend('Ax (cm/s^2)','Ay (cm/s^2)','Axeksp','Ayeksp','location','north');
set(handle, 'Box', 'off');
grid on;

subplot(3,2,5);
%Plott hastigheter i cm/s
plot(tid,Vx*100,tid,Vy*100,tideksp,Vxeksp,tideksp,Vyeksp)
title('Vx og Vy som funksjon av t')
xlabel('t (s)');
handle=legend('Vx (cm/s)','Vy (cm/s)','Vxeksp','Vyeksp','location','north');
set(handle, 'Box', 'off');
grid on;

subplot(3,2,6);
%Plott akselerasjoner i cm/s*s
plot(tid,Ax*100,tid,Ay*100,tideksp,Axeksp,tideksp,Ayeksp)
title('Ax og Ay som funksjon av t')
xlabel('t (s)');
handle=legend('Ax (cm/s^2)','Ay (cm/s^2)','Axeksp','Ayeksp','location','north');
set(handle, 'Box', 'off');
grid on;

