close all
clc
clear all

set(0,'DefaultTextFontname', 'Times New Roman'); set(0,'DefaultTextFontSize', 12);
set(0,'DefaultAxesFontName', 'Times New Roman'); set(0,'DefaultAxesFontSize', 12);
markersize = 2;
linewidth = 1.2;

%===========================================================================================
%                        Sunce, Jupiter, Zemlja i Halejeva kometa
%===========================================================================================

% Konstante
G = 6.67e-11;       % gravitaciona konstanta
auM = 149597870700;     % astronomska jedinica u metrima
auM_dS = auM / (24*60*60);		% astronomska jedinica u metrima kroz dan u sekundama
% Inicijalizacija
tEnd = 75 * 365 * 24 * 60 * 60;      % vremenski period [sekunde]
dt = 1800;      % vremenski korak  [sekunde]

t = 0:dt:tEnd;
stepsNo = length(t);     % ukupan broj koraka
disp(stepsNo);

m1 = 1.989e30;      % masa Sunca [kg]
m2 = 1.898e27;       % masa Jupitera [kg]
m3 = 2.2e14;        % masa Halejeve komete [kg]
m4 = 5.97219e24;		% masa Zemlje [kg]

% Polozaji i brzine

x1 = zeros(1,stepsNo);		x2 = zeros(1,stepsNo);		x3 = zeros(1,stepsNo);		x4 = zeros(1,stepsNo);              
y1 = zeros(1,stepsNo);		y2 = zeros(1,stepsNo);		y3 = zeros(1,stepsNo);      y4 = zeros(1,stepsNo);
z1 = zeros(1,stepsNo);		z2 = zeros(1,stepsNo);		z3 = zeros(1,stepsNo);      z4 = zeros(1,stepsNo);

vx1 = zeros(1,stepsNo);		vx2 = zeros(1,stepsNo);		vx3 = zeros(1,stepsNo);     vx4 = zeros(1,stepsNo);     
vy1 = zeros(1,stepsNo);		vy2 = zeros(1,stepsNo);		vy3 = zeros(1,stepsNo);     vy4 = zeros(1,stepsNo);
vz1 = zeros(1,stepsNo);		vz2 = zeros(1,stepsNo);		vz3 = zeros(1,stepsNo);     vz4 = zeros(1,stepsNo);

% Sunce
x1(1) = 3.555487376887400e-3 * auM;
y1(1) = 3.443726955832962e-3 * auM;
z1(1) = -1.595926442259376e-4 * auM;
vx1(1) = -2.039381680455059e-6 * auM_dS;
vy1(1) = 6.825851787592621e-6 * auM_dS;
vz1(1) = 4.190413544166179e-8 * auM_dS;

% Jupiter
x2(1) = -5.425729474512933 * auM;
y2(1) = -4.893805883481592e-1 * auM;
z2(1) = -1.233739639493046e-1 * auM;
vx2(1) = 5.899930831654800e-4 * auM_dS;
vy2(1) = -7.158477914441830e-3 * auM_dS;
vz2(1) = 1.655140736264202e-5 * auM_dS;

% Halejeva kometa
x3(1) =  -2.046674700505002E+01 * auM;
y3(1) = 2.578916112487541E+01 * auM;
z3(1) = -9.880486589051440E+00 * auM;
vx3(1) = 4.770882574116445E-05 * auM_dS;
vy3(1) = 7.985846752969086E-04 * auM_dS;
vz3(1) = -1.192740886574741E-04 * auM_dS;

% Zemlja
x4(1) =  8.923201619757489E-01 * auM;
y4(1) = 4.523540753710130E-01 * auM;
z4(1) = -1.781780901718696E-04 * auM;
vx4(1) = -8.031919308950355E-03 * auM_dS;
vy4(1) = 1.529775700370547E-02 * auM_dS;
vz4(1) = -5.667177098913429E-07 * auM_dS;

%% Rastojanja
r12 = sqrt((x1(1)-x2(1))^2 + (y1(1)-y2(1))^2 + (z1(1)-z2(1))^2);
r21 = r12;
r13 = sqrt((x1(1)-x3(1))^2 + (y1(1)-y3(1))^2 + (z1(1)-z3(1))^2);
r31 = r13;
r23 = sqrt((x2(1)-x3(1))^2 + (y2(1)-y3(1))^2 + (z2(1)-z3(1))^2);
r32 = r23;
r14 = sqrt((x1(1)-x4(1))^2 + (y1(1)-y4(1))^2 + (z1(1)-z4(1))^2);
r41 = r14;
r24 = sqrt((x2(1)-x4(1))^2 + (y2(1)-y4(1))^2 + (z2(1)-z4(1))^2);
r42 = r24;
r34 = sqrt((x3(1)-x4(1))^2 + (y3(1)-y4(1))^2 + (z3(1)-z4(1))^2);
r43 = r34;

for i=2:stepsNo

    F21 = G*m1*m2/r21^2;
    F12 = F21;            % F12 = G*m1*m2/r12^2
    
    F31 = G*m1*m3/r31^2;
    F13 = F31;            % F13 = G*m1*m3/r13^2
     
    F32 = G*m2*m3/r32^2;
    F23 = F32;            % F23 = G*m2*m3/r23^2
    
	F41 = G*m1*m4/r41^2;
	F14 = F41;
	
	F42 = G*m2*m4/r42^2;
	F24 = F42; 
	
	F43 = G*m3*m4/r43^2;
	F34 = F43;
	
    %% Racunanje komponenti sila
    
    F12x = F12/r12 * (x1(i-1)-x2(i-1));
    F12y = F12/r12 * (y1(i-1)-y2(i-1));
    F12z = F12/r12 * (z1(i-1)-z2(i-1));
    
    F13x = F13/r13 * (x1(i-1)-x3(i-1));
    F13y = F13/r13 * (y1(i-1)-y3(i-1));
    F13z = F13/r13 * (z1(i-1)-z3(i-1));
    
    F21x = -F12x;           % F21x = F21/r21 * (x2(i-1)-x1(i-1))
    F21y = -F12y;           % F21y = F21/r21 * (y2(i-1)-y1(i-1))
    F21z = -F12z;           % F21z = F21/r21 * (z2(i-1)-z1(i-1))
    
    F23x = F23/r23 * (x2(i-1)-x3(i-1));
    F23y = F23/r23 * (y2(i-1)-y3(i-1));
    F23z = F23/r23 * (z2(i-1)-z3(i-1));
    
    F31x = -F13x;                       
    F31y = -F13y;                       
    F31z = -F13z;           
    
    F32x = -F23x;                     
    F32y = -F23y;                       
    F32z = -F23z;              
    
	F14x = F14/r14 * (x1(i-1)-x4(i-1));
	F14y = F14/r14 * (y1(i-1)-y4(i-1));
	F14z = F14/r14 * (z1(i-1)-z4(i-1));
	
	F41x = -F14x;
	F41y = -F14y;
	F41z = -F14z;
	
	F24x = F24/r24 * (x2(i-1)-x4(i-1));
	F24y = F24/r24 * (y2(i-1)-y4(i-1));
	F24z = F24/r24 * (z2(i-1)-z4(i-1)); 
	
	F42x = -F24x;
	F42y = -F24y;
	F42z = -F24z;
	
	F34x = F34/r34 * (x3(i-1)-x4(i-1));
	F34y = F34/r34 * (y3(i-1)-y4(i-1));
	F34z = F34/r34 * (z3(i-1)-z4(i-1));
	
	F43x = -F34x;
	F43y = -F34y;
	F43z = -F34z;
	
    %% Promena brzina
    vx1(i) = (F21x+F31x+F41x)/m1 * dt + vx1(i-1);
    vy1(i) = (F21y+F31y+F41y)/m1 * dt + vy1(i-1);
    vz1(i) = (F21z+F31z+F41z)/m1 * dt + vz1(i-1);
    
    vx2(i) = (F12x+F32x+F42x)/m2 * dt + vx2(i-1);
    vy2(i) = (F12y+F32y+F42y)/m2 * dt + vy2(i-1);
    vz2(i) = (F12z+F32z+F42z)/m2 * dt + vz2(i-1);
    
    vx3(i) = (F13x+F23x+F43x)/m3 * dt + vx3(i-1);
    vy3(i) = (F13y+F23y+F43y)/m3 * dt + vy3(i-1);
    vz3(i) = (F13z+F23z+F43z)/m3 * dt + vz3(i-1);

    vx4(i) = (F14x+F24x+F34x)/m4 * dt + vx4(i-1);
    vy4(i) = (F14y+F24y+F34y)/m4 * dt + vy4(i-1);
    vz4(i) = (F14z+F24z+F34z)/m4 * dt + vz4(i-1);
        
    %% Promena polozaja
    
    x1(i) = vx1(i) * dt + x1(i-1);
    y1(i) = vy1(i) * dt + y1(i-1);
    z1(i) = vz1(i) * dt + z1(i-1);
    
    x2(i) = vx2(i) * dt + x2(i-1);
    y2(i) = vy2(i) * dt + y2(i-1);
    z2(i) = vz2(i) * dt + z2(i-1);
    
    x3(i) = vx3(i) * dt + x3(i-1);
    y3(i) = vy3(i) * dt + y3(i-1);
    z3(i) = vz3(i) * dt + z3(i-1);

    x4(i) = vx4(i) * dt + x4(i-1);
    y4(i) = vy4(i) * dt + y4(i-1);
    z4(i) = vz4(i) * dt + z4(i-1);	
    
    %% Promena rastojanja
    
    r12 = sqrt((x1(i)-x2(i))^2 + (y1(i)-y2(i))^2 + (z1(i)-z2(i))^2);
    r21 = r12;
    r13 = sqrt((x1(i)-x3(i))^2 + (y1(i)-y3(i))^2 + (z1(i)-z3(i))^2);
    r31 = r13;
    r23 = sqrt((x2(i)-x3(i))^2 + (y2(i)-y3(i))^2 + (z2(i)-z3(i))^2);
    r32 = r23;
	r14 = sqrt((x1(i)-x4(i))^2 + (y1(i)-y4(i))^2 + (z1(i)-z4(i))^2);
	r41 = r14;
	r24 = sqrt((x2(i)-x4(i))^2 + (y2(i)-y4(i))^2 + (z2(i)-z4(i))^2);
	r42 = r24;
	r34 = sqrt((x3(i)-x4(i))^2 + (y3(i)-y4(i))^2 + (z3(i)-z4(i))^2);
	r43 = r34;
	
end;

plot3(x1,y1,z1,'r-',x2,y2,z2,'k-',x3,y3,z3,'m-',x4,y4,z4,'g-', 'Markersize', markersize);
view(2);
axis ('equal');
title('Orbite tela u ravni ekliptike');
xlabel('x [m]');
ylabel('y [m]');
legend('Sunce','Jupiter','Halejeva kometa', 'Zemlja');
