% function dy = foucault_dgl(t,y,g,l,omega,lam)
% Simulation und Visualisierung des Foucaultschen Pendels

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers. 2.1
% 07.2015


omega = 2*pi/(24*3600);   % Winkelgeschwindigkeit der Erde
g = 9.81;                 % Erdbeschleunigung
l = 60;                   % m Pendellaenge
%
lam = 70.0;               % geographische Breite
lam = lam*pi/180;
%
% vorf = 2.*omega.*sin(lam);

y0 = [0.01;0;0;0];       % Anfangswert: Auslenkung des Pendels
tmax = 12*3600;          % Laufzeit

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

tic
[t,y] = ode113(@(t,y) foucault_dgl(t,y,g,l,omega,lam),[0, tmax], y0,options);
toc

plot(y(:,1),y(:,2)),axis equal,shg
%% Visualisierung der Drehung der Schwingungsebene
figure,hold on,axis equal,xlim([-y0(1),y0(1)]),ylim([-y0(1),y0(1)]),shg
tic
comet(y(:,1),y(:,2),0.005),shg
toc
%%
% figure,xlim([-y0(1),y0(1)]),ylim([-y0(1),y0(1)]),shg
% hold on
% for k=1:length(y)
%     plot(y(k,1),y(k,2),'k.')
% end
% hold off
