// Projektna naloga - Binarna zvezda
clear
// Vhodni podatki

    // Parametri simulacije
tMax = 200;             // Čas simuliranja [yr]
dt = 0.01;              // Časovni korak [yr]
n = tMax / dt + 1;      // Število korakov
epsilon = 5e-03;

    // Konstante
G0 = 6.67e-11;          // Gravitacijska konstanta [N m^2 / kg^2]
M0 = 1.98892e30;        // Enota za maso sonca [kg]
au = 149597870700;      // Astronomska enota [m]
yr = 365.25 * 24 * 3600 // Leto [s]

    // Pretvorba gravitacijske konstante
G = G0 * M0 * yr ^ 2 / au ^ 3 // [au^3 / M yr^2]

    // Zvezda 1
r1 = zeros(2);          // Položaj zvezde 1 [au]
r1(1) = 0;
r1(2) = 0;
v1 = zeros(2);          // Hitrost zvezde 1 [au/yr]
v1(1) = 1.857;
v1(2) = 0.058;
m1 = 2.063;             // Masa zvezde 1 [M0]
F1 = zeros(2);          // Sila na zvezdo 1 [M0 au / yr^2]
a1 = zeros(2);          // Pospešek zvezde 1 [au / yr^2]
v1abs = zeros(n + 1)    // Absolutna hitrost zvezde 1 [au / yr]

    // Zvezda 2
r2 = zeros(2);          // Položaj zvezde 2 [au]
r2(1) = 7.48;
r2(2) = 5.46;
v2 = zeros(2);          // Hitrost zvezde 2 [au/yr]
v2(1) = -2.216;
v2(2) = -0.1175;
m2 = 1.018;             // Masa zvezde 2 [M0]
F2 = zeros(2);          // Sila na zvezdo 2 [M0 au / yr^2]
a2 = zeros(2);          // Pospešek zvezde 2 [au / yr^2]
v2abs = zeros(n + 1)    // Absolutna hitrost zvezde 2 [au / yr]

u = zeros(2);           // Smerni vektor od zvezde 1 do zvezde 2
e = zeros(2);           // Enotski vektor



// Izpis

i = 0                   // Števec
counter = 0             // Števec
F1x = zeros(n);
F1y = zeros(n);
F2x = zeros(n);
F2y = zeros(n);

a1x = zeros(n);
a1y = zeros(n);
a2x = zeros(n);
a2y = zeros(n);

v1x = zeros(n);
v1y = zeros(n);
v2x = zeros(n);
v2y = zeros(n);
v1x(1) = v1(1);
v1y(1) = v1(2);
v2x(1) = v2(1);
v2y(1) = v2(2);

r1x = zeros(n);
r1y = zeros(n);
r2x = zeros(n);
r2y = zeros(n);
r1x(1) = r1(1);
r1y(1) = r1(2);
r2x(1) = r2(1);
r2y(1) = r2(2);
rRelx = zeros(n);
rRely = zeros(n);
rRelx(1) = r2x(1) - r1x(1);
rRely(1) = r2y(1) - r1y(1);

ts = zeros(n + 1)       // Čas
v1abs(1) = sqrt(v1(1) ^ 2 + v1(2) ^ 2);
v2abs(1) = sqrt(v2(1) ^ 2 + v2(2) ^ 2);

p1 = zeros(1, 3)
p2 = zeros(1, 3)

function odvod = izracunOdvoda(x, y, q)
    odvod = (y(q) - y(q - 1))/(x(q) - x(q - 1));
endfunction

function r = izracunRazdalje (r1, r2)
    x = r2(1) - r1(1);
    y = r2(2) - r1(2);
    r = sqrt(x ^ 2 + y ^ 2);
endfunction

function F1 = izracunSile(e, r)
    F1(1) = G * m1 * m2 * e(1) / r^2;
    F1(2) = G * m1 * m2 * e(2) / r^2;
endfunction


// Program

for t = 0 : dt : tMax

    i = i + 1
    ts(i) = t;
    // Smerni vektor
    u(1) = r2(1) - r1(1);
    u(2) = r2(2) - r1(2);

    // Razdalja med zvezdama
    r = izracunRazdalje(r1, r2);

    // Enotski vektor
    e(1) = u(1) / r;
    e(2) = u(2) / r;

    // Izračun sile
    F1 = izracunSile(e, r);
    F2(1) = -F1(1);
    F2(2) = -F1(2);

    // Izračun pospeška
    a1(1) = F1(1) / m1;
    a1(2) = F1(2) / m1;
    a2(1) = F2(1) / m2;
    a2(2) = F2(2) / m2;

    // Hitrost in lega imata podane začetne vrednosti
    if(i > 1) then
    // Izračun hitrosti
    v1(1) = v1x(i - 1) + a1(1) * dt;
    v1(2) = v1y(i - 1) + a1(2) * dt;
    v2(1) = v2x(i - 1) + a2(1) * dt;
    v2(2) = v2y(i - 1) + a2(2) * dt;
    v1abs(i) = sqrt(v1(1) ^ 2 + v1(2) ^ 2);
    v2abs(i) = sqrt(v2(1) ^ 2 + v2(2) ^ 2);
    // Izračun lege
    r1(1) = r1x(i - 1) + v1(1) * dt;
    r1(2) = r1y(i - 1) + v1(2) * dt;
    r2(1) = r2x(i - 1) + v2(1) * dt;
    r2(2) = r2y(i - 1) + v2(2) * dt;


    end

    // Zapis podatkov za graf
    F1x(i) = F1(1);
    F1y(i) = F1(2);
    F2x(i) = F2(1);
    F2y(i) = F2(2);

    a1x(i) = a1(1);
    a1y(i) = a1(2);
    a2x(i) = a2(1);
    a2y(i) = a2(2);

    if(i > 1) then
    // Izračun hitrosti
    v1x(i) = v1(1);
    v1y(i) = v1(2);
    v2x(i) = v2(1);
    v2y(i) = v2(2);

    // Izračun lege
    r1x(i) = r1x(i - 1) + v1(1) * dt;
    r1y(i) = r1y(i - 1) + v1(2) * dt;
    r2x(i) = r2x(i - 1) + v2(1) * dt;
    r2y(i) = r2y(i - 1) + v2(2) * dt;

    rRelx(i) = r2x(i) - r1x(i);
    rRely(i) = r2y(i) - r1y(i);
    end
    R1x = r1x(i);
    R1y = r1y(i);

    R2x = r2x(i);
    R2y = r2y(i);

    //plot(t, R1x)
end

    // Izračun periode
q = 2
odvod = 100;
while (q < n)
    odvod = izracunOdvoda(r1x, r1y, q);
    if(abs(odvod) < epsilon) then
        counter = counter + 1;
        p1(counter, 1) = r1y(q);
        p1(counter, 2) = q;
        p1(counter, 3) = r1x(q);
        q = q + 20;
     else
        q = q + 1;
     end
end
vsota = 0;
velikost = size(p1);
velikost = velikost(1);
for i = 3 : velikost
    vsota = vsota + p1(i, 3) - p1(i - 2, 3);
end
perioda = vsota / (velikost - 2);
disp(perioda);

// Izris grafov
f1 = figure(1);
plot(r1x, r1y, 'r');
plot(r2x, r2y, 'b');
xlabel("x [au]","fontsize",4,"color","white");
ylabel("y [au]","fontsize",4,"color","white");
title("Lega zvezd","fontsize",4,"color","white");
legend("Zvezda 1", "Zvezda 2");
colordef('black')
f1.background = 0;
/*
f2 = figure(2);
plot(ts, v1abs, 'r');
plot(ts, v2abs, 'b');
xlabel("v [au/yr]","fontsize",4,"color","white");
ylabel("t [yr]","fontsize",4,"color","white");
title("Hitrost zvezd","fontsize",4,"color","white");
legend("Zvezda 1", "Zvezda 2");
colordef('black')
f2.background = 0;

f3 = figure(3);
plot(rRelx, rRely);
xlabel("$\Delta x [yr]$","fontsize",4,"color","white");
ylabel("$\Delta y [yr]$","fontsize",4,"color","white");
title("Relativna lega zvezd","fontsize",4,"color","white");
colordef('black')
f3.background = 0;