%TFY4108 Fysikk hausten 2014
%Øving 4, oppgave 4: Numerisk løsning av bevegelsen
%til en matematisk pendel vha Eulermetoden ("Forward Euler").
%
clear all;
format long;
%
dtkroll = 0.0001; % tidssteget 
tkroll = 0:dtkroll:3*pi; % tidsintervallet
dtheta0 = pi/4000; % steget for maksvinkelen
theta0 = 0:dtheta0:pi/4; % intervallet for maksvinkelen
%
theta = 0*tkroll;
Omega = 0*tkroll;
relavvik = 0*theta0
%
for j = 1:length(theta0)
	%Initialbetingelser for theta og Omega
	theta(1) = j*dtheta0;   % vinkel ved tid tkroll = 0
	Omega(1) = 0.0;    % reskalert vinkelhastigheit ved tid tkroll = 0
	%Løsning med Euler-integrasjon
	for n = 1: length(tkroll)
    		Omega(n+1) = Omega(n) - dtkroll*sin(theta(n));
    		theta(n+1) = theta(n) + dtkroll*Omega(n);
    		if theta(n+1)<0
			relavvik(j) = 100*(1-pi/(2*n*dtkroll));
			break 
    		end
	end
end
%
%Plotter relativt avvik som funksjon av maksvinkel
hold on;
plot(theta0*180/pi,relavvik,'b');
xlabel('Maksvinkel (grader)');
ylabel('Relativt avvik (prosent)');