%
clear all; clc;
% parametry Cauchyovy ulohy
% (kmitani zavazi na pruzine s tlumenim)
% m...hmotnost zavazi,
% k...tuhost pruziny
% d...velikost tlumeni
%-----------------------------%
t0  = 0; % pocatecni cas
m  = 0.001;
k  = 4.2;
d  = 0.001; 
%-----------------------------%
T  = 3.0;
h  = 0.001;
n  = int32(T/h);
pi = 3.1415;
%-----------------------------% 
% pocatecni podminka,sloupec
% timto prevedeme rovnici 2 radu na soustavu
% zavedeme nove promene Z = [z1(t),z2(t)] = [x(t), dx(t)/dt] 
% a F = F(t,Z(t))

Z = [ 0.05; 
      0.00 ];

% RK4
for i=1:n 
   % k vektoru Z[2x1] budeme dale pridavat dalsi vektory
   % pro i=1...n, takze nakonec budeme mit pole 2 x n, 
   % kde kazdy sloupec bude reprezentovat vektor Z v case i   
   %-----------------------------% 
   t  = t0 + i * h;

   %Z_n = [Z_n(1), Z_n(2)] -> takto pristupujeme ke slozkam vektoru Z_n
   Z_n = Z(:,i);   

   %-----------------------------%    
   %k1 = F(t_n,Z_n)
   k1 = [Z_n(2);-d/m*Z_n(2)-k/m*Z_n(1) + 1/(10*t+1)*sin(50*t*pi)];
   
   % k2 = F(t_n + h/2, Z_n + h/2 * k1) = F(t_pom, Z_pom)
   t_pom = t + h/2; Z_pom = Z_n + h/2 * k1;   
   k2 = [Z_pom(2);-d/m*Z_pom(2)-k/m*Z_pom(1) + 1/(10*t+1)*sin(50*t*pi)];
   
   % k3 = F(t_n + h/2, Z_n + h/2 * k2) = F(t_pom, Z_pom)
   t_pom = t + h/2; Z_pom = Z_n + h/2 * k2;
   k3 = [Z_pom(2);-d/m*Z_pom(2)-k/m*Z_pom(1) + 1/(10*t+1)*sin(50*t*pi)];
  
   % k4 = F(t_n + h, Z_n + h * k3) = F(t_pom, Z_pom)
   t_pom = t + h; Z_pom = Z_n + h * k3;
   k4 = [Z_pom(2);-d/m*Z_pom(2)-k/m*Z_pom(1) + 1/(10*t+1)*sin(50*t*pi)];
   
   Z(:,i+1) = Z(:,i) + h/6 * (k1 + 2 * k2 + 2 * k3 + k4);    
end

tt = 0:h:T;
figure (1);
plot(tt, Z(1,:)); % graf trajektorie
figure (2);
plot(Z(1,:), Z(2,:)); % fazovy portret
