%
% Oilerio metodas
%
function main
clc, clear all,
close all
spalva='r';
spalva2='g';
x0=0; % pradinis laiko momentas (s)
%dx=2; % integravimo zingsnis
dx=0.1; % integravimo zingsnis
% dx=0.5; % integravimo zingsnis

global m1 m2 v0 tmax ks ts k1 k2 g;
m1 = 0.3; % pirmo kuno mase
m2 = 0.3; % antro kuno mase
v0 = -60; % pradinis abieju kunu greitis
tmax = 10; % galinis laiko momentas (s)
ks = 0.005; % oro pasipriesinimas abiems
ts = 2; % laikas po kurio issiskiria kunams oro pasipriesinimas
k1 = 0.05; % oro pasipriesinimas pirmam kunui
k2 = 0.01; % oro pasipriesinimas antram kunui
g = 9.8; % laisvo kritimo pagreitis

nsteps=tmax/dx;
figure(1), hold on, grid on;
xlabel('laikas (s)'); ylabel('greitis (v)');
legend('pirmas kūnas', 'antras kūnas');
title('Oilerio metodas');
x=x0;y1=v0;y2=v0;
pntx=x;pnty1=y1;pnty2=y2;
for i=1:nsteps
    dy1=DY1(x,y1);
    y1=y1+dx*dy1;
    dy2=DY2(x,y2);
    y2=y2+dx*dy2;
    x=x+dx;
    plot(x,-y1,[spalva,'.'],'MarkerSize',8)
    plot([pntx,x],[-pnty1,-y1],[spalva,'-']);
    plot(x,-y2,[spalva,'.'],'MarkerSize',8)
    plot([pntx,x],[-pnty2,-y2],[spalva2,'-']);
    pntx=x;pnty1=y1;pnty2=y2;
end
return
function dy1=DY1(x,y1)
    if x > ts
        dy1 = sign(-y1) * k1 * y1.^2 / m1 + 9.8;
    else
        dy1 = sign(-y1) * ks * y1.^2 / (m1 + m2) + 9.8;
    end
return,end

function dy2=DY2(x,y2)
    if x > ts
        dy2 = sign(-y2) * k2 * y2.^2 / m2 + 9.8;
    else
        dy2 = sign(-y2) * ks * y2.^2 / (m1 + m2) + 9.8;
    end
return,end
end