%TFY4115 Fysikk hoesten 2013.
%Oving 4, oppgave 1: Numerisk losning av bevegelsen
%til en matematisk pendel. Eulermetoden ("Forward Euler").
%
%Starthastighet v(1) som gir maksimal alpha naer 180 grader,
%derfor redusert tidssteg dt = 0.000001.
clear all;
format long;
%g = tyngdens akselerasjon (m/s^2)
g = 9.81;
%L = pendelens lengde (m)
L = 1.0;
%dt = tidssteget (s)
dt = 0.000001;
%t = tidsintervallet (s)
t = 0:dt:10;
%Definerer vinkel alpha(t) og hastighet v(t)
alpha = 0*t;
v = 0*t;
%Initialbetingelser for alpha (rad) og v (m/s)
alpha(1) = 0;
v(1) = 6.2641;
%Losning med Euler-integrasjon
for n = 1:(length(t)-1)
    alpha(n+1) = alpha(n) + v(n)*dt/L;
    v(n+1) = v(n) - g*dt*sin(alpha(n));
end;
%Plotter alpha(t)
figure(1);
plot(t,alpha*180/pi,'b');
xlabel('t (s)');
ylabel('vinkelutslag (grader)');
figure(2);
%Plotter E(t)/m = K+U = v^2/2 + g*L*(1-cos(alpha)).
%Setter grenser for vertikal akse til 0 og 2*maksimal energi.
%(max(alpha) returnerer stoerste verdi i tabellen alpha)
plot(t,v.^2/2+g*L.*(1-cos(alpha)));
axis([t(1) t(length(t)) 0 v(1)^2/2+g*L*(1-cos(max(alpha)))]);
xlabel('t (s)');
ylabel('E/m (m^2/s^2)');
%Analyse av perioden T som funksjon av v0, evt alpha_max:
v0 = [0.01 0.1 1 2 3 4 5 6 6.1 6.2 6.25 6.263 6.2641];
alpha_max = [0.003 0.032 0.322 0.653 1.003 1.391 1.855 ...
    2.566 2.692 2.869 3.024 3.104 3.132];
periode = [2.01 2.01 2.02 2.06 2.14 2.28 2.54 3.41 3.70 ...
    4.29 5.29 6.88 8.60];
figure(3)
plot(alpha_max*180/pi,periode);
xlabel('max vinkel (grader)');
ylabel ('periode (s)');
axis([0 180 0 10]);
