% 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("Rutherforduv pokus\n");
% -----------------------------------------------
   % definice pocatecnich podminek a konstant
u=1.67e-27;     % hmotnostni atomova jednotka
q=1.602e-19;    % elementarni naboj
eps0=8.854e-12; % premitivita vakua
ma=4*u;         % hmotnost alpha-castice
Qa=2*q;         % naboj alpha-castice
Qj=79*q;        % naboj jadra Au
kQQ=1/4/pi/eps0*Qa*Qj; % koeficient v Coulombove zakone
ymax=2e-11;     % maximalni hodnota y pro vykreslovani   
t0=0;
x0=-1e-11;
y0=[-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 9 10]*1e-12;         
vx0=1e6;
vy0=0;
% -----------------------------------------------
dt=1e-20;       % integracni krok
   % -----------------------------------------------
   % parametry obrazku
clearplot
axis([-11 11 -21 21],'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(pm)"
__gnuplot_set__ ylabel "y(pm)"
__gnuplot_set__ title "Rutherforduv pokus"
hold on % pro vice kreivek v jednom obrazku
% -----------------------------------------------
   % cyklus pro nekolik ruznych trajektorii
cas_spusteni = clock();
for j=1:columns(y0) 
% -----------------------------------------------
% vlastni numericka integrace - Eulerova metoda
% -----------------------------------------------
clear x
clear y
t(1)=t0;
x(1)=x0;
y(1)=y0(j);
vx=vx0;
vy=vy0;
i=1;
% -----------------------------------------------
 while ((abs(y(i))<=ymax) & (x(i)>=x0)) 
   i=i+1;
   x(i)=x(i-1)+vx*dt;
   y(i)=y(i-1)+vy*dt;
   r=sqrt(x(i)^2+y(i)^2); 
   ax=kQQ/r^3*x(i-1)/ma;
   ay=kQQ/r^3*y(i-1)/ma;
   vx=vx+ax*dt;
   vy=vy+ay*dt;
   t(i)=t(i-1)+dt;
 endwhile
% -----------------------------------------------
plot(x/1e-12,y/1e-12,'1');
endfor
hold off
% -----------------------------------------------
celkovy_cas = etime(clock(), cas_spusteni);
disp(['Program bezel ',num2str(celkovy_cas),' sekund.']);
% -----------------------------------------------

