% 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("Pruzinove kyvadlo\n");
% -----------------------------------------------
% definice pocatecnich podminek a konstant
l=0.3;     % delka kyvadla
k=11.5;    % tuhost gumy
g=9.81;    % tihove zrychleni
m=0.12;    % hmotnost zavazi
t0=0;
%x0=-.1;
%y0=-.3;         
x0=.3;
y0=-0.2;         
vx0=0;
vy0=0;
% -----------------------------------------------
dt=1e-2;       % integracni krok
mez=20;        % mez integrace
% -----------------------------------------------
% vlastni numericka integrace - Eulerova metoda
% -----------------------------------------------
cas_spusteni = clock();
t(1)=t0;
x(1)=x0;
y(1)=y0;
vx=vx0;
vy=vy0;
i=1;
% -----------------------------------------------
 while (t(i)<=mez) 
   i=i+1;
   r=sqrt(x(i-1)^2+y(i-1)^2);
   ax=-k*(r-l)/r*x(i-1)/m;
   ay=-g-k*(r-l)/r*y(i-1)/m;
   vx=vx+ax*dt;
   vy=vy+ay*dt; 
   x(i)=x(i-1)+vx*dt;
   y(i)=y(i-1)+vy*dt;
   t(i)=t(i-1)+dt;
 endwhile
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% -----------------------------------------------
   % parametry obrazku
clearplot
axis([min(x) max(x) min(y) 0], '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
__gnuplot_set__ ytics
% popisky obrazku
__gnuplot_set__ xlabel "x(m)"
__gnuplot_set__ ylabel "y(m)"
__gnuplot_set__ title "Pruzinove kyvadlo"
plot(x,y);
% %%%%%%%%%%%%%%%%%%%%%%%% End of file %%%%%%%%%%%%%%%%%%%%%%%%
