% plotSolarSystem.m
% Plots the trajectories returned by solarSystem
% J.M.Soler, Jan.2011

% Set times for trajectory
clear, clf
t1 = 0;     % initial time for trajectories, from Jan.1,2022, in days
t2 = 700;   % final time, in days
dt = 1;     % integration time interval, in days

% Find trajectories
[r,body] = solarSystem( t1, t2, dt );

% Plot trajectories
[nd,nb,nt] = size(r);   % number of dimensions, bodies, and times
x=zeros(nt,1); y=x; z=x;
for ib = 1:nb
    x(:)=r(1,ib,:);
    y(:)=r(2,ib,:);
    z(:)=r(3,ib,:);
    plot3(x,y,z), hold on
%    t = (1:30:nt);                 % once every month
%    plot3(x(t),y(t),z(t),'o')      % plot the body position as an 'o'
end
axis equal, grid on
xlabel('x (au)'), ylabel('y (au)'), zlabel('z (au)')

% Find eccentricities
for ib = 1:nb
    d = sqrt(r(1,ib,:).^2+r(2,ib,:).^2+r(3,ib,:).^2);
    dmin(ib) = min(d);
    dmax(ib) = max(d);
end
eccentricities = [body,num2str((dmax-dmin)'./(dmax+dmin)')]

% Find orbital periods (times of closest approach to initial point)
% Find distances to initial positions
for it = 1:nt
    delta_r(:,:,it) = r(:,:,it)-r(:,:,1);
end
% Find velocities
v(:,:,2:nt-1) = (r(:,:,3:nt)-r(:,:,1:nt-2)) / (2*dt);
v(:,:,1) = v(:,:,2);
v(:,:,nt) = v(:,:,nt-1);
% Find the times delta_t to/from closest approach, such that 
% delta_x=v*delta_t => delta_t=delta_x/v=delta_r*v/v^2
delta_t = zeros(nb,nt);
delta_t(:,:) = sum(delta_r.*v,1) ./ sum(v.^2,1);
% Find periods
period = zeros(nb,1);
for ib = 1:nb
    % Find the time interval containing the first zero of delta_t
    it = find(delta_t(ib,1:nt-1)<0 & delta_t(ib,2:nt)>0,1);
    if (numel(it)==0)   % no zero of delta_t found
        fprintf('%s%s\n','Trajectory shorter than period of ',body(ib,:))
    else
        % Linear interpolation between it and it+1
        % d=d1+(d2-d1)*(t-t1)/dt=0 => t=t1-dt*d1/(d2-d1)
        period(ib) = (it-1)*dt - ...
                     dt*delta_t(ib,it)/(delta_t(ib,it+1)-delta_t(ib,it));
    end
end
orbital_periods_in_days = [body,num2str(period)]
