%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% volleyball.m:
% Volleyball paa kvartsirkel.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nullstiller alle variable
clear;      
% Leser eksperimentelle data fra filen volleyball.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('volleyball.txt','%f %f %f','headerlines',2);
% Omregning av posisjon fra cm til m:
xeksp = xeksp/100;
yeksp = yeksp/100;
% Funksjonen length(tabell) returnerer antall elementer i tabell
lengde = length(tideksp);
% Lager tabeller med eksperimentelle hastighetsverdier (m/s).
% 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 phieksp og reksp med eksperimentelle verdier for 
% vinkelen phi (radianer) og avstanden fra origo til ballens massesenter (m):
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
% Lager omegaeksp = Veksp/reksp (1/s):
omegaeksp = sqrt(Vxeksp.*Vxeksp + Vyeksp.*Vyeksp)./reksp;
% Lager alphaeksp = d omegaeksp/dt (1/s^2), tilsvarende som Vxeksp:
alphaeksp = zeros(1,lengde);
for i = 2:lengde-1,
	alphaeksp(i) = (omegaeksp(i+1) - omegaeksp(i-1))/(tideksp(i+1) - tideksp(i-1));
end
alphaeksp(lengde)=(omegaeksp(lengde)-omegaeksp(lengde-1))/(tideksp(lengde)-tideksp(lengde-1));

% Numeriske beregninger basert paa teori
% Diverse parametre:
g=9.81;      % tyngdens akselerasjon (m/s^2)
dt=0.01;     % tidsintervall (s) (juster etter behov)
R=0.79;      % sirkelradius (m)
mys=0.8;     % statisk friksjonskoeffisient (juster etter behov)
myk=0.7;     % kinetisk friksjonskoeffisient (juster etter behov)
r=0.10;      % radius volleyball (m)
c=2/3;       % 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

% Startverdier (dvs ved tidspunkt t = tid(1) = 0.0):
% Vinkel, fart osv gjelder for ballens CM.
tid(1)   = 0.0;     % Starttid
phi(1)   = 0.055;   % Startvinkel (rad)
V(1)     = 0.040;     % Starthastighet (m/s)
OMEGA(1) = V(1)/(r+R); % Vinkelhastighet ved t=0 (1/s)
n(1)     = g*cos(phi(1));  % Normalkraft ved t=0 (N/kg = m/s^2)
A(1)     = g*sin(phi(1))/(c+1); % Akselerasjon ved t=0 (m/s^2)
ALPHA(1) = A(1)/(r+R);  % Vinkelakselerasjon ved t=0 (1/s^2)
%Kartesiske komponenter av posisjon (m) og hastighet (m/s) ved t=0:
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));
% Forbereder oppstart av while-loop. Betingelsene for aa fortsette
% while-loopen er (1) at normalkraften paa ballen er storre enn null,
% (2) at vinkelen phi er mindre enn en fastsatt maksverdi phimax, og
% (3) at antall tidssteg j er mindre enn en fastsatt maksverdi jmax.
j=1;   % Tabellindeks
phimax=pi/2; % Maksimal vinkel (rad)
jmax=20/dt;  % Maksimalt antall tidssteg; forsikring mot uendelig loop
jrulling=1;  % Startverdi for jrulling; sluttverdi for jrulling gir 
             % tidspunkt for overgang til sluring 
jkontakt=1;  % Startverdi for jkontakt; sluttverdi for jkontakt gir
             % tidspunkt for overgang til fritt fall
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 vinkelhastighet OMEGA
  phi(j+1) = phi(j) + OMEGA(j)*dt; % Ny vinkel phi, osv
  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); % Akselerasjon, rulling antatt
  ALPHA(j+1) = A(j+1)/(r+R);      % Vinkelakselerasjon
  f(j+1) = c*A(j+1);              % Friksjonskraft, rulling antatt
  % Sjekker om beregnet f, rulling antatt, har blitt storre enn fmax. 
  % Hvis ja: Sett f = myk*n og endre 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, sluring
      A(j+1) = g*sin(phi(j+1)) - f(j+1); % Akselerasjon, sluring
      ALPHA(j+1) = A(j+1)/(r+R);         % Vinkelakselerasjon
  else
      jrulling=j+1; % Tar vare paa siste j-verdi med ren rulling
  end;
  jkontakt=j+1;     % Tar vare paa siste j-verdi der ballen har kontakt
  tid(j+1) = dt*j;  % Neste tidsverdi
  j = j+1;          % Neste tabellindeks-verdi
end

% Ute av while-loopen.
% Det betyr normalt 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 while-loopen selv om noe "gaar galt".)

% Visualisering av resultatene.

% Figur 1: Eksperimentell og numerisk beregnet vinkel.
figure(1);
plot(tideksp,phieksp*180/pi,tid,phi*180/pi);
title('Figur 1. Eksperimentell og numerisk beregnet vinkel')
xlabel('t (s)');
handle=legend('\phi_{exp} (grader)','\phi_{num} (grader)','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 2: Eksperimentell og numerisk beregnet vinkelhastighet.
figure(2);
plot(phieksp*180/pi,omegaeksp*180/pi,phi*180/pi,OMEGA*180/pi);
title('Figur 2. Eksperimentell og numerisk beregnet vinkelhastighet')
xlabel('\phi (grader)');
handle=legend('\Omega_{exp} (grader/s)','\Omega_{num} (grader/s)','location','northwest');
set(handle, 'Box', 'off');
grid on;

% Figur 3: Eksperimentell og numerisk beregnet vinkelakselerasjon.
figure(3);
plot(phieksp*180/pi,alphaeksp*180/pi,phi*180/pi,ALPHA*180/pi)
title('Figur 3. Eksperimentell og numerisk beregnet vinkelakselerasjon.')
xlabel('\phi (grader)');
handle=legend('\alpha_{exp} (grader/s^2)','\alpha_{num} (grader/s^2)','location','southwest');
set(handle, 'Box', 'off');
grid on;

% Figur 4: Numerisk beregnet normalkraft pr masseenhet.
figure(4);
plot(phi*180/pi,n)
title('Figur 4. Numerisk beregnet normalkraft pr masseenhet.')
xlabel('\phi (grader)');
ylabel('N (N/kg)');
grid on;

% Figur 5: Ballens massesenter, eksperimentelt.
figure(5);
plot(phieksp*180/pi,reksp*100)
title('Figur 5. Ballens massesenter, eksperimentelt.')
xlabel('\phi (grader)');
ylabel('R_{CM} (cm)');
grid on;

% End
