%
clear, clc

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers. 1.0
% 10.2014

%% Differentialgleichung loesen
y0 = [0.5*pi;0];                                % Anfangswerte
tmax=20;                                        % Maximale Simulationsdauer
Vorf = 4; % omega0^2 -> g/l
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
tic
[t,y] = ode45(@(t,y) realesPendelDGL(t,y,Vorf),[0,tmax],y0,options);
toc
%% Visualisieren
x1=sin(y(:,1));
x2=-cos(y(:,1));
figure(1),plot(x1,x2),hold on, plot(x1(1),x2(1),'pr'),axis equal, shg
figure(2)
subplot(2,1,1)
plot(t,y(:,1)/pi*180),hold on
plot(t,y0(1)/pi*180*cos(sqrt(Vorf)*t),'m'),shg
%% Bestimmen der Eigenfrequenzen
% die Eventfunktion sorgt fuer die genaue Berechnung einer Periode
tic
winkel = [0.01 0.1 1, 10:10:170, 175, 179];  % Startwinkel
n=0;
for k=winkel
    n=n+1;
    y0 = [k/180*pi;0];
    options = odeset('RelTol',1e-8,'AbsTol',1e-8,'event',@(t,y) realesPendelevent(t,y,y0(1)));
    [t2,y2,te,ye,ie] = ode45(@(t,y) realesPendelDGL(t,y,Vorf),[0,tmax],y0,options);
    omega(n) = 2*pi/(te(2)-te(1));
end
figure(2),
subplot(2,1,2)
plot(winkel,omega),shg
figure(2),subplot(2,1,1)
plot(t2,y2(:,1)/pi*180,'g')
toc