% script file ex27a.m
% Compute outer solar system
% last changes: 12.01.02
% need to include monitoring of Hamiltonian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear all
close all

global m 
global G

% constants:
G = 2.95912208286e-4;

% masses:
m = [1.00000597686     ;... % Sun
     0.000954786104043 ;... % Jupiter
     0.000285583733151 ;... % Saturn
     0.0000437273164546;... % Uranus
     0.0000517759138449;... % Neptune
     1/(1.3e+8)        ];   % Pluto

% initial positions:
q0 = [  0           0           0        ; ... % Sun
       -3.5023653  -3.8169847  -1.5507963; ... % Jupiter
        9.0755314  -3.0458353  -1.6483708; ... % Saturn
        8.3101420 -16.2901086  -7.2521278; ... % Uranus
       11.4707666 -25.7294829 -10.8169456; ... % Neptune
      -15.5387357 -25.2225594  -3.1902382];    % Pluto

% initial velocities:
v0 = [0           0           0         ; ... % Sun
      0.00565429 -0.00412490  0.00190582; ... % Jupiter
      0.00168318  0.00483525  0.00192462; ... % Saturn
      0.00354178  0.00137102  0.00055029; ... % Uranus
      0.00288930  0.00114527  0.00039677; ... % Neptune
      0.00276725 -0.00170702 -0.00136504];    % Pluto

% momentum = mass*velocity
% initial momentum - each row gets multiplied with the
% corresponding mass:
p0 = diag(m)*v0; 

%p0 = v0; % p is velocity

% time domain:
h = 10;
T_end = 200000;
t = [0:h:T_end];
n = length(t);

Q.E = zeros(n,6,3); 
P.E = zeros(n,6,3);

H.E = zeros(n,1); % Hamiltonian

for i = 1:6  
  Q.E(1,i,:) = q0(i,:);
  P.E(1,i,:) = p0(i,:);
end

% initial Hamiltonian
T = 0;
for i = 1:6
  pe = zeros(1,3);
  pe(:) = P.E(1,i,:);
  T = T + (pe*pe')/m(i);
end
T = T/2;
  
U = 0;
for i = 2:6
  for j = 1:i-1
    difer = zeros(1,3);
    difer(:) = Q.E(1,i,:)-Q.E(1,j,:);  
    U = U + m(i)*m(j)/(norm(difer));
  end
end
U = -G*U;  
H.E(1) = T+U;    


% Euler
for time = 1:n-1
  for k = 1:6
    Q.E(time+1,k,:) = Q.E(time,k,:)+h*P.E(time,k,:)/m(k);
    
    dU = 0;
    for i=1:6
      if i~=k
	difer = zeros(1,3);
	difer(:) = Q.E(time,k,:)-Q.E(time,i,:);
	r = norm(difer);
	dU = dU + m(k)*m(i)*(Q.E(time,k,:)-Q.E(time,i,:))/(r^3);
      end
    end
    dU = G*dU;
    
    P.E(time+1,k,:) = P.E(time,k,:)-h*dU;
  end
  
  % compute Hamiltonian
  
  T = 0;
  for i = 1:6
     pe = zeros(1,3);
     pe(:) = P.E(time+1,i,:);
     T = T + (pe*pe')/m(i);
  end
  T = T/2;
  
  U = 0;
  for i = 2:6
    for j = 1:i-1
      difer = zeros(1,3);
      difer(:) = Q.E(time+1,i,:)-Q.E(time+1,j,:);
      U = U + m(i)*m(j)/(norm(difer));
    end
  end
  U = -G*U;
  
  H.E(time+1) = T+U;    
end

% plot the result
figure('name','Euler')
subplot(2,1,1)
hold on
for k = 1:6
  plot3(Q.E(:,k,1),Q.E(:,k,2),Q.E(:,k,3))
end
view(3)
grid on

subplot(2,1,2)
plot(t,H.E')
title('Hamiltonian')


% Sympletic Euler

Q.SE = zeros(n,6,3); 
P.SE = zeros(n,6,3)
H.SE = zeros(n,1); % Hamiltonian

for i = 1:6  
  Q.SE(1,i,:) = q0(i,:);
  P.SE(1,i,:) = p0(i,:);
end

H.SE(1) = H.E(1);

for time = 1:n-1
  for k = 1:6
    Q.SE(time+1,k,:) = Q.SE(time,k,:)+h*P.SE(time,k,:)/m(k);
  end
  
  for k = 1:6
    dU = 0;
    for i=1:6
      if i~=k
	difer = zeros(1,3);
	difer(:) = Q.SE(time+1,k,:)-Q.SE(time+1,i,:);
	r = norm(difer);
	dU = dU + m(k)*m(i)*(Q.SE(time+1,k,:)-Q.SE(time+1,i,:))/(r^3);
      end
    end
    dU = G*dU;
    
    P.SE(time+1,k,:) = P.SE(time,k,:)-h*dU;
  end
  
  % compute Hamiltonian
  
  T = 0;
  for i = 1:6
    pe = zeros(1,3);
    pe(:) = P.SE(time+1,i,:);
    T = T + (pe*pe')/m(i);
  end
  T = T/2;
  
  U = 0;
  for i = 2:6
    for j = 1:i-1
      difer = zeros(1,3);
      difer(:) = Q.SE(time+1,i,:)-Q.SE(time+1,j,:); 
      U = U + m(i)*m(j)/(norm(difer));
    end
  end
  U = -G*U;
  
  H.SE(time+1) = T+U;      
end

% plot the result
figure('name','Sympletic Euler')
subplot(2,1,1)
hold on
for k = 1:6
  plot3(Q.SE(:,k,1),Q.SE(:,k,2),Q.SE(:,k,3))
end
view(3)
grid on

subplot(2,1,2)
plot(t,H.SE')
title('Hamiltonian')

