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

    // Parametri simulacije
tMax = 200;            // Čas simuliranja [yr]
dt = 0.01;                 // Časovni korak [yr]
nMax = 1000;            // Maksimalno število iteracij

    // 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]

    // 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]

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



// Izpis
n = tMax / dt + 1;
i = 1                   // Š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);



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

    // 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 ;

    // 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;


    end
    R1x = r1x(i);
    R1y = r1y(i);

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

    //plot(t, R1x)
end

plot(r1x, r1y);
plot(r2x, r2y);