%TFY4106 Fysikk 2014.
%Øving 7, oppgave 2: Numerisk løsning av bevegelsen
%til en matematisk pendel. Eulermetoden ("Forward Euler").
%
clear all;
format long;
%g = tyngdens akselerasjon (m/s^2)
g = 9.81;
%L = pendelens lengde (m)
L = 1.0;
%Ser her på tilfelle (B):
%Starthastighet v(1)=6.2641 m/s gir max theta nær 180 grader.
%Tidssteg så lite som dt = 1e-6 s viser seg å være nødvendig
%for å gi periodisk oppførsel over (minst) 2 hele perioder.
%For tilfelle (A), bruk f.eks: dt = 0.0001, T = 2*pi/sqrt(g/L), v(1) = 1.0 m/s
%For tilfelle (C), bruk f.eks: dt = 0.0001, T = 1.5 s, v(1) = 6.5 m/s
%dt = 0.0001;
dt = 0.000001;
%dt = 0.0001;
%t = tidstabell fra 0 til 2*perioden (s)
%t = 0:dt:2*2*pi/sqrt(g/L);
t = 0:dt:17.20;
%t = 0:dt:3.0;
%Nullstill tabeller for vinkel theta(t) og hastighet v(t):
theta = 0*t;
v = 0*t;
%Initialbetingelser for theta (rad) og v (m/s)
theta(1) = 0;
v(1) = 6.2641;
E(1) = v(1)^2/2+g*L*(1-cos(theta(1)));
%Løsning med Euler-integrasjon
for n = 1:(length(t)-1)
    theta(n+1) = theta(n) + v(n)*dt/L;
    v(n+1) = v(n) - g*dt*sin(theta(n));
    %Mekanisk energi:
    E(n+1) = v(n+1)^2/2+g*L*(1-cos(theta(n+1)));
end;
%Plotter theta(t)
figure(1);
plot(t,theta*180/pi,'b');
xlabel('t (s)');
ylabel('\theta (grader)');
title('Tilfelle (B)');
figure(2);
%Plotter E(t)/m = (K+U)/m = v^2/2 + g*L*(1-cos(theta)).
%Setter grenser for vertikal akse til 0 og 1.5*maksimal energi.
%(max(E) returnerer største verdi i tabellen E)
plot(t,E);
axis([t(1) t(length(t)) 0 1.5*max(E)]);
xlabel('t (s)');
ylabel('E/m (m^2/s^2)');
title('Tilfelle (B)');
%
%Ekstra (IKKE etterspurt i øvingen):
%Analyse av perioden T som funksjon av v0, evt theta_max:
v0 = [0.01 0.1 1 2 3 4 5 6 6.1 6.2 6.25 6.263 6.2641];
theta_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(theta_max*180/pi,periode);
xlabel('\theta_{max} (grader)');
ylabel ('T (s)');
axis([0 180 0 10]);
%
%Kommentar til slutt: Med, som her, 17.2 millioner tidssteg får man
%vel en fornemmelse av at programmet ikke er spesielt effektivt.
%Mye raskere ville det nok ha gått med bruk av Matlabs innebygde
%difflignings-løsere.