#!/usr/bin/octave

# Megadjuk a gyorsulást, az erőtörvényből
# A másodrendű diffegyenletet egy két elsőrendűből álló rendszerré alakítjuk
function xdot=xdot(t,z)
    global g;
    global l;
    global tomeg;
    x=z(1);
    v=z(2);
    gyorsulas = -1/l * g *sin(x);
    xdot=[v,gyorsulas];
endfunction

# Értéket adunk az állandóknak
global g;
global l;
global tomeg;
g=9.81;
l=1/2;
tomeg=1;

# Felcímkézzük a tengelyeket
xlabel("t");
ylabel("x");
hold on;

t0=0;
t1=3;

# A differenciálegyenlet-megoldó RK-módszer paraméterei
vopt = odeset ("RelTol", 1e-8, "AbsTol", 1e-8, "NormControl", "on","InitialStep",1e-2,"MaxStep",0.1);

# Kiszámoljuk az inga mozgását különböző kezdősebességek esetén, mindig a nulla kitéréstől indítva az ingát
for v0=0:2:12
    z0=[0,v0];
    [t z] = ode45(@xdot,[t0,t1],z0,vopt);
    x=z(:,1);
    plot(t,x);
    hold on;
end
# Kirajzoljuk a határesetet is, amikor az energia éppen elég a maximum eléréséhez
z0=[0,2*sqrt(g/l)];
[t z] = ode45(@xdot,[t0,t1],z0,vopt);
x=z(:,1);
plot(t,x);
    
print("ingamozgas.ps","-dps","-color","-landscape")
