clear, close, clc
format long
g = 9.82;
l = 0.5;
fel = 1;
%% fråga a: beräkna periodtiden T för motsvarande theta0, gör tabell.
n = 10; %antal trapetser, valfritt SÄTT FIXERAT JUST NU.
Ttab = []; %Variabel med alla perioder
runs = 1;
vinklar = [5:5:90];
for theta0 = pi/36:pi/36:pi/2 % intervall för vinkeln som massan släpps från
    k = sin(theta0/2); %k värdet givet i uppgift
    f=@(alpha)1/(sqrt(1-k^2*sin(alpha)^2)); %funktionen i integralen
    a = 0; %lower limit
    b = pi/2; %upper limit
    %trapetsmetoden säger; hitta ett avstånd h mellan punkter för att få
    %regelbundna trapetser under kurvan för att kunna beräkna dess
    %area/integral.
    h = (b-a)/(n); %regelbundna steglängd mellan trapetserna
    s = 0.5*(f(a)+f(b)); %summan av första delen av trapetsmetoden.
    for i = 1:n-1
        s = s + f(a + i*h);
    end %antal trapetser varierar beroende på n som vi väljer
    % större n -> mer exakt svar
    Itheta0 = h * s; %integralen klar här.
    T = 4*sqrt(l/g)*Itheta0;
    Ttab = [Ttab T];
end
T1 = table(Ttab',vinklar');
T1.Properties.VariableNames = {'Periodtid','Vinklar'};
disp(T1)
plot(Ttab(5),s)
figure
%% b) för små vinklar theta kan approximeringen theta = sin(theta)
%som ger periodtiden Ta = 2*pi*sqrt(l/g)
%relativa felet;
Ta = 2*pi*sqrt(l/g);
rel_fel = abs(Ttab-Ta)./Ttab;
rel_fel = rel_fel';
felprocent = rel_fel.*100 %här står procentuellt fel
%% c) Runge-kuttas metod
g = 9.82;
l = 0.5;
h = 0.1;
t = [0:h:2]; %tiden
T = t(end);
vt = [0 2];
theta25 = (25*pi)/180; %begynnelsevinkel 25
theta50 = (50*pi)/180; % 50
v25 = [0]; %begynnelse hastighet
v50 = [0];
fv = @(t, theta, v) -g/l * sin(theta); %dv/dt
ftheta= @(t, theta, v) v; %dtheta/dt
for i =1:T/h
    k1v25 = h*fv(t(i), theta25(i), v25(i));
    k1theta25 = h*ftheta(t(i), theta25(i), v25(i));
    k2v25 = h*fv(t(i)+h/2, theta25(i) + k1theta25/2, v25(i) + k1v25/2);
    k2theta = h*ftheta(t(i)+h/2, theta25(i) + k1theta25/2, v25(i) + k1v25/2);
    k3v25 = h*fv(t(i)+h/2, theta25(i) + k2theta/2, v25(i)+k2v25/2);
    k3theta = h*ftheta(t(i)+h/2, theta25(i) + k2theta/2, v25(i)+k2v25/2);
    k4v25 = h*fv(t(i)+h,theta25(i)+k3theta, v25(i)+k3v25);
    k4theta = h*ftheta(t(i)+h,theta25(i)+k3theta, v25(i)+k3v25);
    tot_theta =theta25(i) + (1/6)*(k1theta25 + 2*k2theta + 2*k3theta + k4theta);
    theta25= [theta25 tot_theta];
    tot_v25 = v25(i) + (1/6)*(k1v25 + 2*k2v25 + 2*k3v25 + k4v25);
    v25 = [v25  tot_v25];

     k1v50 = h*fv(t(i), theta50(i), v50(i));
    k1theta50 = h*ftheta(t(i), theta50(i), v50(i));
    k2v50 = h*fv(t(i)+h/2, theta50(i) + k1theta50/2, v50(i) + k1v50/2);
    k2theta50 = h*ftheta(t(i)+h/2, theta50(i) + k1theta50/2, v50(i) + k1v50/2);
    k3v50 = h*fv(t(i)+h/2, theta50(i) + k2theta50/2, v50(i)+k2v50/2);
    k3theta50 = h*ftheta(t(i)+h/2, theta50(i) + k2theta50/2, v50(i)+k2v50/2);
    k4v50 = h*fv(t(i)+h,theta50(i)+k3theta50, v50(i)+k3v50);
    k4theta50 = h*ftheta(t(i)+h,theta50(i)+k3theta50, v50(i)+k3v50);
    tot_theta =theta50(i) + (1/6)*(k1theta50 + 2*k2theta50 + 2*k3theta50 + k4theta50);
    theta50= [theta50 tot_theta];
    tot_v50 = v50(i) + (1/6)*(k1v50 + 2*k2v50 + 2*k3v50 + k4v50);
    v50 = [v50  tot_v50];

end
subplot(2,1,1);
plot(t,theta25,'linewidth',2)
xlabel('Tid (s)','fontsize',15)
ylabel('Vinkel (rad)','fontsize',15)
hold on
plot(t,theta50,'linewidth',2)
title('Theta(t)','fontsize',15)
legend('25','50')

% periodtiden

%% ode45 ---------  uppgift e
subplot(2,1,2);
theta25_0 = (25*pi)/180; %begynnelsevinkel 25
theta50_0 = (50*pi)/180;
[tt,yy] = ode45(@(tt,yy) diff(tt,yy),[0 2],[theta25_0 0] );
[ttt,yyy] = ode45(@(tt,yy) diff(tt,yy),[0 2],[theta50_0 0] );
plot(tt, yy(:,1),'linewidth',2)
title('Ode45 Check', 'fontsize',15)
xlabel('Tid (s)','fontsize',15)
ylabel('Vinkel (rad)','fontsize',15)
hold on
plot(ttt,yyy(:,1),'linewidth',2)
legend('25','50')


%% funktion för ode45--------
function dydt = diff(tt,yy)
dydt= zeros(2, 1);
dydt(1)=yy(2);
dydt(2)=(-9.82/0.5)*sin(yy(1));
end