% Main_pendel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parametrar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L           =     5;    % Lengde
g           =  9.81;    % Gravitasjon   
N           =  1000;    % Antall tidsskritt
h           =  0.02;    % Skrittlengde

Theta0      =  pi/3;    % Initialbetingelse
dt_Theta    =     0;    % Initialbetingelse fart

Y_0=[Theta0;dt_Theta];  % Initialverdier i nye variabler.
YY(:,1)=Y_0;            % Förste kollonne i en matrise som vil
                        % brukes for å putte resultatene i.
                        % Antal kolloner vil öke for var tidsskritt.
                        % YY vil få störelsen 2 x (n+1) -matrise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start, iterasjon med Euler metode
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N
    YY(:,i+1)=YY(:,i)+h*Pendel(YY(:,i),g,L);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Löser ligningen med inbygd metode 
% (ode45) med höy nöjaktighet 
% ('RelTol'=1e-13,'AbsTol'=1e-16)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_span=[0:h:h*N];   % Tider [0,h,2h,...,N*h] da 
                    % vi vill aproksimere funksjonen.
options = odeset('RelTol',1e-13,'AbsTol',1e-16);
[T_ode45,Y_ode45] = ode45(@(t,Y) Pendel(Y,g,L),T_span,Y_0,options);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot:
%---------------------------------------------------
% Plotter å lager film der Eulers metode (YY) 
% sammanlignes med Matlabs metode (y_ode45).
% I filmen som lages (Eksempel_film.avi) ser vi
% at amplituden fra Eulers metode öker fra
% numersk feil.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Opner en avi-fil med navnet Eksempel_film.avi
mov = avifile('Eksempel_film.avi','compression','Cinepak','fps',30);

for i=1:(N+1)
    
    % Plotter den numeriske lösningen fra Matlabs 'ode45'
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hold off % Sletter forige figuren
    x_o=L*sin(Y_ode45(i,1)); % x-verdi for tidsskritt i
    y_o=-L*cos(Y_ode45(i,1)); % y-verdi for tidskritt i
    % plotte med röd (r) strek (-) å sirkel (o) '-or'
    plot(([0,x_o]),[0,y_o],'-or') 

    % Plotter den numeriske lösningen med Eulers metode
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hold on % Sparer forige figur
    x=L*sin(YY(1,i)); % x-verdi for tidsskritt i
    y=-L*cos(YY(1,i));
    % plotte med svart (k) strek (-) å sirkel (o) '-or'
    % Tyckelsen for linjen settes til 2 ('LineWidth',2)
    plot(([0,x]),[0,y],'-ok','LineWidth',2)
    
    % Axis setter skalerer akslarna i figuren 
    axis ([-(1.1*L) (1.1*L) -(1.1*L) (0.5*L)]); 
    
    %%%
    F = getframe(gca); 
    mov = addframe(mov,F);  % Legg til figur til film
    %%%
    
    pause(0.02)
    
end
mov = close(mov);  % Steng film

% Slut