%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% losning6.m:
% Oving 6 FY1001/TFY4145 Mekanisk fysikk 2014:
% Rullende og slurende bordtennisball paa kvartsirkel.
% Plott tilsvarende losning6.pdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nullstiller alle variable
clear;      
% Leser eksperimentelle data fra filen bordtennisball.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('bordtennisball.txt','%f %f %f','headerlines',2);
% Funksjonen length(tabell) returnerer antall elementer i tabell
lengde = length(tideksp);
dteksp = (tideksp(lengde)-tideksp(1))/(lengde-1);
% Lager tabeller med eksperimentelle hastighetsverdier og energiverdier.
% Starter med aa lage "tomme tabeller".
Vxeksp = zeros(1,lengde);
Vyeksp = zeros(1,lengde);
Veksp = zeros(1,lengde);
Keksp = zeros(1,lengde);
Ueksp = zeros(1,lengde);
Eeksp = zeros(1,lengde);
% Trenger verdiene for c og g allerede her:
c = 2/3;
g = 9.81;
for i = 2:lengde-1, % For forste og siste element beregnes det 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));
    Veksp(i) = sqrt(Vxeksp(i)*Vxeksp(i)+Vyeksp(i)*Vyeksp(i));
    % Keksp=kinetisk energi gyldig ved ren rulling
    % Faktorer 100: omregning fra cm til m
    Keksp(i) = (c+1)*Veksp(i)*Veksp(i)/(2*100*100);
    Ueksp(i) = g*yeksp(i)/100;
    Eeksp(i) = Keksp(i)+Ueksp(i);
end
% Lager eksperimentelle verdier ogsaa for forste og siste element for 
% aa unngaa "skjemmende" nullverdier i grafene
Vxeksp(lengde)=(xeksp(lengde)-xeksp(lengde-1))/(tideksp(lengde)-tideksp(lengde-1));
Vyeksp(lengde)=(yeksp(lengde)-yeksp(lengde-1))/(tideksp(lengde)-tideksp(lengde-1));
Veksp(lengde)=sqrt(Vxeksp(lengde)*Vxeksp(lengde)+Vyeksp(lengde)*Vyeksp(lengde));
Keksp(lengde)=(c+1)*Veksp(lengde)*Veksp(lengde)/(2*100*100);
Ueksp(lengde)=g*yeksp(lengde)/100;
Eeksp(lengde)=Keksp(lengde)+Ueksp(lengde);
% I forste element legger vi for enkelhets skyld inn samme verdi som
% i element nr 2:
Vxeksp(1)=Vxeksp(2);
Vyeksp(1)=Vyeksp(2);
Veksp(1)=Veksp(2);
Keksp(1)=Keksp(2);
Ueksp(1)=Ueksp(2);
Eeksp(1)=Eeksp(2);
% Lager ogsaa 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 (allerede definert ovenfor)
phimax=pi/2; % maksimal vinkel
% Neste linjer sorger for at numerisk tidssteg dt blir lik eksperimentelt
% tidssteg dteksp dividert med et heltall storre enn eller lik 1, og paa
% en slik maate at numerisk tidssteg dt ikke blir storre enn 0.01 sek.
% Kommandoen floor(tall) avrunder tall til naermeste heltall som er
% mindre enn eller lik tall.
dtfaktor = 1+floor(dteksp/0.01);
dt=dteksp/dtfaktor;
R=0.785;     % sirkelradius 78.5 cm
% Parametre for: bordtennisball
mys=0.25;    % statisk friksjonskoeffisient, kan justeres
myk=0.25;    % kinetisk friksjonskoeffisient, kan justeres
r=0.019;     % radius bordtennisball
c=2/3;       % legemets treghetsmoment pr masseenhet er i0=c*r*r
% Diverse analytiske uttrykk og sammenhenger, kun for oversiktens skyld 
% 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.047;   % Startvinkel for bordtennisball
tid(1)   = 0.0;     % Starttid
V(1)     = 0.015;     % Starthastighet (for CM)
OMEGA(1) = V(1)/(r+R);     % Startvinkelhastighet (for CM om origo)
omega(1) = V(1)/r;         % Startvinkelhastighet (om CM)
% vrelativ = V - omega*r = hastigheten til legemets kontaktpunkt mot
% underlaget; vrelativ = 0 saa lenge vi har ren rulling, deretter blir
% vrelativ positiv naar legemet begynner aa slure
vrelativ(1) = V(1) - omega(1)*r;
n(1)     = g*cos(phi(1));  % Normalkraft ved start
A(1)     = g*sin(phi(1))/(c+1); % Akselerasjon ved start
f(1)     = c*A(1); % Friksjonskraft ved start (rulling antatt)
ALPHA(1) = A(1)/(r+R);  % Vinkelakselerasjon 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));
% Initialisering av energier, numerisk:
kintrans(1) = V(1)*V(1)/2;
kinrot(1) = c*r*r*omega(1)*omega(1)/2;
potensiell(1) = g*y(1);
% pwork(j) = friksjonsarbeid pr tidsenhet utfort i kontaktpunktet 
% i tidssteg nr j; dette blir her et slags effekttap pr masseenhet,
% med enheten J/(s*kg) = W/kg
pwork(1) = 0.0;
% fwork(j) = sum_(i=1)^(i=j) pwork(i)*dt = totalt utfort friksjonsarbeid
% fram til og med tidssteg nr j
fwork(1) = 0.0;
% totalmekanisk = reell mekanisk energi, siden kinrot er uttrykt ved
% r*omega (som ved sluring er mindre enn V)
totalmekanisk(1) = kintrans(1)+kinrot(1)+potensiell(1);
% totalenergi = total mekanisk energi + totalt utfort friksjonsarbeid;
% totalenergi skal vaere konstant for alle tidssteg inklusive sluring
totalenergi(1) = totalmekanisk(1)+fwork(1);
% rulleenergi = total mekanisk energi HVIS vi antar ren rulling;
% rulleenergi vil vaere lik totalenergi inntil legemet begynner aa slure;
% ved sluring vil rulleenergi bli storre enn totalenergi siden uttrykket
% for rulleenergi antar at omega*r=V; man kan sammenligne Eeksp og
% rulleenergi, siden begge storrelser antar ren rulling.
rulleenergi(1) = (c+1)*V(1)*V(1)/2 + potensiell(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
  % Hvis ren rulling: omega = V/r:
  omega(j+1) = V(j+1)/r;
  vrelativ(j+1) = V(j+1) - omega(j+1)*r;
  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));
  % Energiforhold (se forklaringer ovenfor):
  kintrans(j+1) = V(j+1)*V(j+1)/2;
  kinrot(j+1) = c*r*r*omega(j+1)*omega(j+1)/2;
  potensiell(j+1) = g*y(j+1);
  rulleenergi(j+1) = (c+1)*V(j+1)*V(j+1)/2 + potensiell(j+1);
  totalmekanisk(j+1) = kintrans(j+1)+kinrot(j+1)+potensiell(j+1);
  n(j+1) = g*cos(phi(j+1)) - V(j+1)*V(j+1)/(r+R);
  % Siste n-verdi kommer til aa bli negativ slik programmet er skrevet.
  % Vi endrer derfor n til null hvis den har blitt negativ.
  if n(j+1) < 0
      n(j+1)=0;
  end
  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)
      % Vi maa rette opp diverse storrelser hvis sluring har begynt
      % Vi maa ogsaa regne ut friksjonsarbeidet for gjeldende tidssteg.
      omega(j+1) = omega(j) + f(j)*dt/(c*r);
      vrelativ(j+1) = V(j+1)-omega(j+1)*r;
      f(j+1) = myk*n(j+1);  % Friksjonskraft ved sluring
      pwork(j+1) = f(j)*vrelativ(j); % Effekttap pga sluring
      fwork(j+1) = fwork(j) + pwork(j+1)*dt; % Akkumulert friks.arb.
      kinrot(j+1) = c*r*r*omega(j+1)*omega(j+1)/2;
      totalmekanisk(j+1) = kintrans(j+1)+kinrot(j+1)+potensiell(j+1);
      totalenergi(j+1) = totalmekanisk(j+1)+fwork(j+1);
      A(j+1) = g*sin(phi(j+1)) - f(j+1); % Akselerasjon ved sluring
      ALPHA(j+1) = A(j+1)/(r+R);
  else
      pwork(j+1) = 0.0;  % ren rulling, dermed intet effekttap
      fwork(j+1) = 0.0; % se forrige linje
      totalenergi(j+1) = totalmekanisk(j+1)+ fwork(j+1);
      jrulling=j+1; % Sparer paa siste j-verdi der vi hadde ren rulling
  end;
  jkontakt=j; % 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 aa avslutte loopen selv om noe gaar galt.

% Herfra og ut handler det om aa visualisere resultatene.
% Her lages forst en figur med 2x2=4 plott.

figure(1);

subplot(2,2,1);
plot(tid,phi*180/pi,tid,OMEGA*180/pi,tid,n,tideksp,phieksp*180/pi);
title('\phi, \Omega, N og \phi^{exp} som funksjon av tid')
xlabel('Tid (s)');
%matlab-spesifikt (neste 2 linjer)
handle=legend('\phi (grader)','\Omega (grader/s)','N (N/kg)','\phi^{exp} (grader)','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('\phi_r=',num2str(phi(jrulling)*180/pi),' t_r=',num2str(tid(jrulling)));
text(10,300,label);
label=strcat('\phi_k=',num2str(phi(jkontakt)*180/pi),' t_k=',num2str(tid(jkontakt)));
text(10,200,label);
% Fjerne % i de neste linjene hvis utskrift av ulike 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('Hastighetskomponenter som funksjon av tid')
xlabel('Tid (s)');
ylabel('Hastighet (cm/s)');
handle=legend('V_x^{num}','V_y^{num}','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('Hastighetskomponenter som funksjon av vinkel')
xlabel('\phi (grader)');
ylabel('Hastighet (cm/s)');
handle=legend('V_x^{num}','V_y^{num}','V_x^{exp}','V_y^{exp}','location','southwest');
set(handle, 'Box', 'off');
grid on;

% Figur 2:
% Animasjon av bordtennisballens bevegelse nedover kvartsirkelen,
% numerisk (blaa sirkler, 'ob') og eksperimentelt (roede stjerner, '*r'):
figure(2);
axis([0 125*R 0 125*R])
antnum = jkontakt;
anteksp = lengde;
antanim = floor(antnum/dtfaktor);
hold on;
xlabel('x (cm)');
ylabel('y (cm)');
title('Rulling av bordtennisball nedover kvartsirkel');
grid on;
for j=1:antanim
    numtid = tid((j-1)*dtfaktor+1);
    eksptid = tideksp(j);
    xekspanim = xeksp(j);
    yekspanim = yeksp(j);
    xanim = x((j-1)*dtfaktor+1)*100;
    yanim = y((j-1)*dtfaktor+1)*100;
    % Kontrollerer at tidspunktet er likt numerisk og eksperimentelt:
    if abs(eksptid-numtid) >= dt
        print('Ulike tidspunkter!');
    end
    plot(xanim,yanim,'ob',xekspanim,yekspanim,'*r');  
    pause(0.1);  % Forsinker framvisningen i 0.1 s (juster etter behov)
end
handle=legend('Numerisk','Eksperimentelt','location','northeast');
set(handle, 'Box', 'off');

% Figur 3:
% Numerisk og eksperimentell hastighet som funksjon av vinkel
figure(3);
plot(phi*180/pi,V*100,'ob-',phieksp*180/pi,Veksp,'*r-');
title('Numerisk og eksperimentell hastighet som funksjon av vinkel')
xlabel('Vinkel (grader)');
ylabel('Hastighet (cm/s)');
handle=legend('V^{num}','V^{exp}','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 4:
% Eksperimentell avstand fra kvartsirkelens origo til legemets CM,
% som funksjon av vinkel
figure(4);
plot(phieksp*180/pi,reksp,'ob-');
title('Eksperimentell avstand fra origo til CM som funksjon av vinkel')
xlabel('Vinkel (grader)');
ylabel('Avstand fra origo til CM (cm)');
grid on;

% Figur 5:
% Energiforholdene som funksjon av tiden:
figure(5);
plot(tid,totalmekanisk,tid,totalenergi,tideksp,Eeksp,tid,rulleenergi);
title('Numerisk og eksperimentell energi som funksjon av tid')
xlabel('Tid (s)');
ylabel('Energi pr masseenhet (J/kg)');
handle=legend('E^{mek}_{tot}','E_{tot}','E_{exp} (rulling antatt!)','E_{num} (rulling antatt!)','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 6:
% Energiforholdene som funksjon av vinkelen:
figure(6);
plot(phi*180/pi,totalmekanisk,phi*180/pi,totalenergi,phieksp*180/pi,Eeksp,phi*180/pi,rulleenergi);
title('Numerisk og eksperimentell energi som funksjon av vinkel')
xlabel('Vinkel (grader)');
ylabel('Energi pr masseenhet (J/kg)');
handle=legend('E^{mek}_{tot}','E_{tot}','E_{exp} (rulling antatt!)','E_{num} (rulling antatt!)','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 7:
% Akkumulert friksjonsarbeid som funksjon av vinkelen:
figure(7);
plot(phi*180/pi,fwork);
title('Totalt friksjonsarbeid som funksjon av vinkel');
xlabel('Vinkel (grader)');
ylabel('Arbeid pr masseenhet (J/kg)');
handle=legend('Totalt friksjonsarbeid','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 8:
% Effekttap pga sluring som funksjon av vinkelen:
figure(8);
plot(phi*180/pi,pwork);
title('Effekttap pga sluring som funksjon av vinkel');
xlabel('Vinkel (grader)');
ylabel('Effekttap pr masseenhet (W/kg)');
handle=legend('Effekttap','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 9:
% Relativhastighet i kontaktpunktet som funksjon av vinkelen:
figure(9);
plot(phi*180/pi,vrelativ*100);
title('Kontaktpunktets hastighet som funksjon av vinkel');
xlabel('Vinkel (grader)');
ylabel('Hastighet (cm/s)');
handle=legend('V_{rel}','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 10:
% Beregnet friksjonskraft f som funksjon av vinkelen:
figure(10);
plot(phi*180/pi,f);
title('Friksjonskraft som funksjon av vinkel');
xlabel('Vinkel (grader)');
ylabel('Kraft pr masseenhet (N/kg)');
handle=legend('f','location','northwest');
set(handle, 'Box', 'off');
grid on;
