#!/usr/bin/octave

function v01=v01(A)
    global omega0;
    global v0a;
    omega2=-A^2*omega0/16;
    omega=omega0+omega2;
    v01=omega*A + A^3/192*3*omega - v0a;
end


function Ampl=Ampl(v)
    global v0a;
    v0a=v;
    Ampl=fsolve(@v01,1.);
end

function omega2=omega2(A)
    global g;
    global l;
    global omega0;
    omega2=-g*A^2/16/l/omega0;
end

global omega0;
global g;
global l;

g=9.81;
l=1;

omega0=sqrt(g/l);


# Megoldjuk az ingamozgás egyenletét különböző amplitúdókra
# 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

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

t0=0;
t1=2.5;

# 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:.6:3
    z0=[0,v0];
    [t z] = ode45(@xdot,[t0,t1],z0,vopt);
    x=z(:,1);
    Am=Ampl(v0);
    omega2=-Am^2*omega0/16;
    omega=omega0+omega2;
    xa0=v0/omega0*sin(omega0*t);
    xa=Am*sin(omega*t) - Am^3/192*sin(3*omega*t);
# Kirajzoljuk, az első esetben feliratozva (a színek magyarázata), utána anélkül
    if (v0==0)
	plot(t,x,"-@x1;egzakt;");
	hold on;
	plot(t,xa,"-@+2;elso korr.;");
	hold on;
	plot(t,xa0,"-@o3;harmonikus;");
	hold on;
    else
	plot(t,x,"-@x1;;");
	hold on;
	plot(t,xa,"-@+2;;");
	hold on;
	plot(t,xa0,"-@o3;;");
	hold on;
    endif
end

print("inga.jpg","-djpg","-color","-landscape")
#print("ingamozgas.ps","-dps","-landscape","-solid")
