% numericke reseni Foucaultovo kyvadlo
%
clear;
clc;
# vymazeme nastaveni vykreslovani
__gnuplot_raw__ "reset;";
# nastavuje opet vykresleni plnou carou
__gnuplot_set__ data style lines;
__gnuplot_set__ xrange [:]
__gnuplot_set__ yrange [:]
% -----------------------------------------------
printf("kmity Foucaultova kyvadla\n");
% -----------------------------------------------
% definice pocatecnich podminek
t0=0;
x0=0;
y0=0;         
vx0=3;
vy0=0;
global phi=49/180*pi;       % zemepisna sirka
global om=2*pi/(86400);     % uhlova rychlost rotace Zeme
global l=68;                % delka kyvadla
global g=9.81;
printf("doba kmitu pro 68 m je %f\n",2*pi*sqrt(l/g));
% -----------------------------------------------
% vlastni numericka integrace - preddefinovana procedura programu octave
% -----------------------------------------------
% resi jako ctyrrozmerny vektor xp o slozkach
% xp(1)=x, xp(2)=vx, xp(3)=y, xp(4)=vy
% piseme diff. rovnice
% xdot(1)=1. der. x
% xdot(2)=2. der. x
% xdot(3)=1. der. y
% xdot(4)=2. der. y
function [xdot]=fce(xp,t)    
 global phi om g l % prebira globalne definovane konstanty
  xdot=zeros(4,1);
  xdot(1)=xp(2);
  xdot(2)=-g/l*xp(1)+2*om*sin(phi)*xp(4);  
  xdot(3)=xp(4);
  xdot(4)=-g/l*xp(3)-2*om*sin(phi)*xp(2);
endfunction  
% -----------------------------------------------
% meze integrace
dm=0;
hm=60;
% -----------------------------------------------
cas_spusteni = clock();
% cas jako nezavisla promenna
t=linspace(dm,hm,1000);
yp=lsode("fce",[x0,vx0,y0,vy0],t);
x=yp(:,1);
y=yp(:,3);
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% -----------------------------------------------
% vykresleni zavislosti
% vykresleni osy x
%axis('equal');
__gnuplot_set__ zeroaxis lt -1
__gnuplot_set__ nokey
__gnuplot_set__ tics out
__gnuplot_set__ grid
__gnuplot_set__ mxtics
__gnuplot_set__ mytics
__gnuplot_set__ xtics axis
__gnuplot_set__ ytics axis
__gnuplot_set__ noborder
% popisky obrazku
__gnuplot_set__ xlabel "x(m) - smer na jih"
__gnuplot_set__ ylabel "y(m) - smer na vychod"
__gnuplot_set__ title "Kmity Foucaultova kyvadla"
plot(x,y);
