% numericke reseni pruzneho kyvadla
%
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.28;     % delka kyvadla
k=12.5;    % tuhost gumy
g=9.81;    % tihove zrychleni
m=0.12;    % hmotnost zavazi
t0=0;
lg=l+m*g/k
omp=sqrt(g/lg)
Omr=sqrt(k/m)
2*omp/Omr
d0=0.0;
%x0=-.1;
%y0=-.3;      
phi0=0.05785;
%phi0=5e-3;   
x0=(lg+d0)*sin(phi0);
y0=-(lg+d0)*cos(phi0);         
vx0=0;
vy0=0;
% -----------------------------------------------
h=5e-3;       % integracni krok
mez=50;        % 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);
   k1x=-k*(r-l)/r*x(i-1)/m;
   k1y=-g-k*(r-l)/r*y(i-1)/m;
   xp=x(i-1)+vx*2*h/3+k1x*2*h^2/9;
   yp=y(i-1)+vy*2*h/3+k1x*2*h^2/9;
   rp=sqrt(xp^2+yp^2);
   k2x=-k*(rp-l)/rp*xp/m;
   k2y=-g-k*(rp-l)/rp*yp/m;
   x(i)=x(i-1)+vx*h+(k1x+k2x)*h^2/4;
   y(i)=y(i-1)+vy*h+(k1y+k2y)*h^2/4;
   vx=vx+(k1x+3*k2x)*h/4;
   vy=vy+(k1y+3*k2y)*h/4;
   t(i)=t(i-1)+h;
 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 5
__gnuplot_set__ mytics 5
__gnuplot_set__ decimalsign ","
__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 %%%%%%%%%%%%%%%%%%%%%%%%
