%solar_system_orbits
% Various representations of the solar system.
% LAST UPDATED by Andy French September 2022

function solar_system_orbits

animation_option = 'animation';  %Options are 'animation + avi', 'animation', 'none'
planets_option = 'inner';  %Options are 'inner','outer', 'all'
dimensions_option = '2D';  %Options are '2D' or '3D'
background_option = 'white';  %Background colour is either black or white
object_at_origin = 'earth';  %Plot orbits relative to this object

%Run solar system plotter
solar_sytem_orbit_plotter( animation_option, planets_option, dimensions_option, background_option, object_at_origin );

%%

%Solar system plotter
function solar_sytem_orbit_plotter( animation_option, planets_option, dimensions_option, background_option, object_at_origin  )

%Define planets to plot, and compute orbit time in years to encompass
%largest orbit
if strcmp( planets_option, 'inner' )
    planets = {'mercury','venus','moon','earth','mars','sun'};
    t = linspace(0,30,1e4);
elseif strcmp( planets_option, 'outer' )
    planets = {'sun','jupiter','ganymede','uranus','neptune','pluto'};
    t = linspace(0,1000,5e4);
else
    planets = {'sun','mercury','venus','earth','moon','mars','jupiter','uranus','neptune','pluto'};
    t = linspace(0,300,5e4);
end

if strcmp( background_option, 'black' )==1
    %Create figure with a black background, with white axes
    fig = figure('name','solar system','color',[0 0 0],'InvertHardCopy', 'off');
    axes('nextplot','add','fontsize',14,...
        'color',[0 0 0],'xcolor',[1 1 1],'ycolor',[1 1 1],'zcolor',[1 1 1]);
else
    %Create figure with a white background, with black axes
    fig = figure('name','solar system','color',[1 1 1],'InvertHardCopy', 'off');
    axes('nextplot','add','fontsize',14,...
        'color',[1 1 1],'xcolor',[0 0 0],'ycolor',[0 0 0],'zcolor',[0 0 0]);
end

%

%% Plot elliptical orbit of planets %%

%Initialize arrays for x,y,z coordinates of planets
x = zeros( length(t), length(planets) ); y = zeros(size(x)); z = zeros(size(y));

%Planet marker colour
planet_colour = zeros( length(planets), 3 );

%Compute orbit of object at the origin
[xo,yo,zo,thetao,po] = calc_orbit( object_at_origin, t, '3D' );

%Compute orbits of planets
for n=1:length(planets)
    
    %Calculate orbit, relative to object at the origin
    [xp,yp,zp,theta,p] = calc_orbit( planets{n}, t, '3D' );
    xp = xp - xo; yp = yp - yo; zp = zp - zo;
    
    %Update plot
    if strcmp( dimensions_option, '3D' )
        plot3(xp,yp,zp,'color',p.rgb,'linewidth',2);
    else
        plot(xp,yp,'color',p.rgb,'linewidth',2);
    end
    
    %Update planet x,y,z coordinate array
    x(:,n) = xp; y(:,n) = yp; z(:,n) = zp; planet_colour(n,:) =  p.rgb;
end

%Add plot annotation
axis equal; grid on; box on;
if strcmp( background_option, 'black' )==1; RGB = [1 1 1]; else RGB = [0 0 0]; end
legend(planets,'fontsize',14,'textcolor',RGB,'edgecolor',RGB);
xlabel('x /AU'); ylabel('y /AU');
tit = title([planets_option,' solar system relative to ',object_at_origin],'color',RGB);
if strcmp( dimensions_option, '3D' )
    view(3); axis vis3d;  zlabel('z /AU');
end

%Plot object at the origin
if strcmp( planets_option, 'all' ) == 0
    if strcmp( object_at_origin , 'sun' )
        plot3(0,0,0,'oy','markersize',10,'markerfacecolor','y');
        plot3(0,0,0,'y*','markersize',30);
    else
        if strcmp( dimensions_option, '3D' )
            plot3( 0,0,0,'o','markerfacecolor',po.rgb,'color',po.rgb );
        else
            plot( 0,0,'o','markerfacecolor',po.rgb,'color',po.rgb );
        end
    end
end

%Print PNG file
print( gcf,[planets_option,' solar system ',...
    dimensions_option,' ',object_at_origin,'.png'],'-dpng','-r300');

%

%Animate orbits if option set
if strcmp( animation_option, 'animation' ) || strcmp( animation_option, 'animation + avi' )
    
    %Plot initial planet positions
    for n=1:length(planets)
        if strcmp( dimensions_option, '3D' )
            plt(n) = plot3( x(1,1),y(1,1),z(1,1),'o','markerfacecolor',planet_colour(n,:),'color',planet_colour(n,:) );
        else
            plt(n) = plot( x(1,1),y(1,1),'o','markerfacecolor',planet_colour(n,:),'color',planet_colour(n,:) );
        end
    end
    
    %Open video file if avi option set
    if strcmp( animation_option , 'animation + avi' )
        video= VideoWriter( [planets_option,' solar system.avi'] );
        video.FrameRate = 30;
        open(video);
    end
    
    %Initialize time frame counter
    i = 1; 
    
    %Loop through time t and update planet positions
    while ishandle(fig)     
        for n=1:length(planets)
            if strcmp( dimensions_option, '3D' )
                set( plt(n), 'xdata',x(i,n),'ydata',y(i,n),'zdata',z(i,n) );
            else
                set( plt(n), 'xdata',x(i,n),'ydata',y(i,n) );
            end
        end
        set( tit, 'string',['Solar system: t=',num2str(t(i),3),' years'] );
        drawnow;
        
        %Write video frame if option set
        if strcmp( animation_option, 'animation + avi' )
            print( gcf, 'videoframe.png','-dpng','-r200');
            [image,map] = imread('videoframe.png','png');
            writeVideo(video,image);
        end
        i = i+1;
        if i>length(t)
            i = 1;
        end
    end
end

%Close video file if option set
if strcmp( animation_option , 'animation + avi' )
    close(video);
    delete( 'videoframe.png' );
end

%%

%Compute orbit in x,y,z coordinates
function [x,y,z,theta,p] = calc_orbit( planet, t, dimensions_option )
p = get_planet_data( planet );
theta = angle_vs_time( t, p.P, p.ecc, p.theta0 );
r = p.a*(1-p.ecc^2)./( 1 - p.ecc*cos( theta ) );
x = r.*cos(theta); y = r.*sin(theta); z = zeros(size(x));
if strcmp( dimensions_option, '3D' )
    beta = p.beta*pi/180;
    x = x*cos(beta); z=x*sin(beta);
end

%Special cases for Moons
if strcmp( planet, 'moon' )==1
    [xx,yy,zz,thetaa,pp] = calc_orbit( 'earth', t, dimensions_option );
    x = x + xx; y = y + yy; z = z + zz;
elseif strcmp( planet, 'ganymede' )==1
    [xx,yy,zz,thetaa,pp] = calc_orbit( 'jupiter', t, dimensions_option );
    x = x + xx; y = y + yy; z = z + zz;
end

%%

%Numeric method to compute polar angle vs orbit time
function theta = angle_vs_time( t, P, ecc, theta0 )

%Number of orbits
N = ceil( t(end)/P );
if P==0; theta = zeros(size(t)); return; end

%Angle step for Simpson's rule
dtheta = N*2*pi/length(t);

%

%Define array of polar angles for orbits
theta = theta0 : dtheta : ( 2*pi*N + theta0 );

%Evaluate integrand of time integral
f = (1 - ecc*cos(theta) ).^(-2);

%Define Simpson rule coefficients c = [ 1, 4, 2, 4, 2, 4, ....1 ]
L = length(theta);
isodd = rem( 1:(L-2),2 ); isodd( isodd==1 ) = 4; isodd( isodd==0 ) = 2;
c = [1, isodd, 1];

%Calculate array of times
tt = P*( (1-ecc^2)^(3/2) )*(1/(2*pi))*dtheta*(1/3).*cumsum( c.*f );

%Interpolate the polar angles for the eccentric orbit at the circular orbit
%times
theta = interp1( tt, theta, t, 'spline' );

%%

%Get planet data from data structure
function planet = get_planet_data( planet_name )
[sun,mercury,venus,earth,moon,mars,...
    jupiter,ganymede,saturn,uranus,neptune,pluto] =...
    solar_system_parameters;
eval([ 'planet = ',planet_name,';'] );

%%

%Solar system parameters
function [sun,mercury,venus,earth,moon,mars,...
    jupiter,ganymede,saturn,uranus,neptune,pluto] =...
    solar_system_parameters

% NOTE: don't use rand for theta0 as this function is called numerous times
% during the code. theta0 needs to be fixed for a given simulation run.

%Solar system physical parameters
ME = 5.9742e14; %Mass of the Earth /kg
MS = 1.9891e30; %Solar mass /kg
RS = 6.960e8; %Solar radius /m
RE = 6.37814e6; %Earth radius /m
AU = 1.495979e11; %1AU /m
Yr = 365*24*3600;  %1Yr /s
G = 6.67e-11; %Gravitational constant G /Nm^2 kg^-2

%Sun
sun.m = 0.055;   %Mass in Earth masses
sun.a = 0;  %Semi-major axis in AU
sun.ecc = 0; %Orbital eccentricity
sun.beta = 0; %Orbital inclination /deg
sun.theta0 = 0;  %Initial polar angle at t=0
sun.R = 109.123; %Radius /Earth radii
sun.trot = 0; %Rotational period /days
sun.P = 0; %Orbital period /years
sun.rgb = (1/255)*[255,255,0];

%Mercury
mercury.m = 0.055;   %Mass in Earth masses
mercury.a = 0.387;  %Semi-major axis in AU
mercury.ecc = 0.21; %Orbital eccentricity
mercury.beta = 7.00; %Orbital inclination /deg
mercury.theta0 = 2*pi*0;  %Initial polar angle at t=0
mercury.R = 0.383; %Radius /Earth radii
mercury.trot = 58.646; %Rotational period /days
mercury.P = 0.241; %Orbital period /years
mercury.rgb = (1/255)*[128,128,128];

%Venus
venus.m = 0.815;   %Mass in Earth masses
venus.a = 0.723;  %Semi-major axis in AU
venus.ecc = 0.01; %Orbital eccentricity
venus.beta = 3.39; %Orbital inclination /deg
venus.theta0 = 2*pi*0;  %Initial polar angle at t=0
venus.R = 0.949; %Radius /Earth radii
venus.trot = 243.018; %Rotational period /days
venus.P = 0.615; %Orbital period /years
venus.rgb = (1/255)*[212,204,44];

%Earth
earth.m = 1;   %Mass in Earth masses
earth.a = 1;  %Semi-major axis in AU
earth.ecc = 0.02; %Orbital eccentricity
earth.beta = 0; %Orbital inclination /deg
earth.theta0 = 2*pi*0;  %Initial polar angle at t=0
earth.R = 1; %Radius /Earth radii
earth.trot = 1; %Rotational period /days
earth.P = 1; %Orbital period /years
earth.rgb = (1/255)*[0,0,255];

%Moon
moon.m = 0.0123;   %Mass in Earth masses
moon.a = 0.00257;  %Semi-major axis in AU
moon.ecc = 0.0549; %Orbital eccentricity
moon.beta = 5.145; %Orbital inclination /deg
moon.theta0 = 2*pi*0;  %Initial polar angle at t=0
moon.R = 0.2727; %Radius /Earth radii
moon.trot = 27.322; %Rotational period /days
moon.P = 27.322/365.25; %Orbital period /years
moon.rgb = (1/255)*[169,169,169];

%Mars
mars.m = 0.107;   %Mass in Earth masses
mars.a = 1.523;  %Semi-major axis in AU
mars.ecc = 0.09; %Orbital eccentricity
mars.beta = 1.85; %Orbital inclination /deg
mars.theta0 = 2*pi*0;  %Initial polar angle at t=0
mars.R = 0.533; %Radius /Earth radii
mars.trot = 1.026; %Rotational period /days
mars.P = 1.881; %Orbital period /years
mars.rgb = (1/255)*[255,0,0];

%Jupiter
jupiter.m = 317.85;   %Mass in Earth masses
jupiter.a = 5.202;  %Semi-major axis in AU
jupiter.ecc = 0.05; %Orbital eccentricity
jupiter.beta = 1.31; %Orbital inclination /deg
jupiter.theta0 = 2*pi*0;  %Initial polar angle at t=0
jupiter.R = 11.209; %Radius /Earth radii
jupiter.trot = 0.413; %Rotational period /days
jupiter.P = 11.861; %Orbital period /years
jupiter.rgb = (1/255)*[188,68,68];

%Ganymede
ganymede.m = 0.025;   %Mass in Earth masses
ganymede.a = 1070400e3/AU;  %Semi-major axis in AU (about Jupiter)
ganymede.ecc = 0.0013; %Orbital eccentricity
ganymede.beta = 2.214; %Orbital inclination /deg
ganymede.theta0 = 2*pi*0;  %Initial polar angle at t=0
ganymede.R = 0.4213; %Radius /Earth radii
ganymede.trot = 7.154; %Rotational period /days
ganymede.P = 7.154/365.25; %Orbital period /years
ganymede.rgb = (1/255)*[169,169,169];

%Saturn
saturn.m = 95.159;   %Mass in Earth masses
saturn.a = 9.576;  %Semi-major axis in AU
saturn.ecc = 0.06; %Orbital eccentricity
saturn.beta = 2.49; %Orbital inclination /deg
saturn.theta0 = 2*pi*0;  %Initial polar angle at t=0
saturn.R = 9.449; %Radius /Earth radii
saturn.trot = 0.444; %Rotational period /days
saturn.P = 29.628; %Orbital period /years
saturn.rgb = (1/255)*[244,194,12];

%Uranus
uranus.m = 14.5;   %Mass in Earth masses
uranus.a = 19.293;  %Semi-major axis in AU
uranus.ecc = 0.05; %Orbital eccentricity
uranus.beta = 0.77; %Orbital inclination /deg
uranus.theta0 = 2*pi*0;  %Initial polar angle at t=0
uranus.R = 4.007; %Radius /Earth radii
uranus.trot = 0.718; %Rotational period /days
uranus.P = 84.747; %Orbital period /years
uranus.rgb = (1/255)*[14,242,231];

%Neptune
neptune.m = 17.204;   %Mass in Earth masses
neptune.a = 30.246;  %Semi-major axis in AU
neptune.ecc = 0.01; %Orbital eccentricity
neptune.beta = 1.77; %Orbital inclination /deg
neptune.theta0 = 2*pi*0;  %Initial polar angle at t=0
neptune.R = 3.883; %Radius /Earth radii
neptune.trot = 0.671; %Rotational period /days
neptune.P = 166.344; %Orbital period /years
neptune.rgb = (1/255)*[1,98,255];

%Pluto
pluto.m = 0.003;   %Mass in Earth masses
pluto.a = 39.509;  %Semi-major axis in AU
pluto.ecc = 0.25; %Orbital eccentricity
pluto.beta = 17.5; %Orbital inclination /deg
pluto.theta0 = 2*pi*0;  %Initial polar angle at t=0
pluto.R = 0.187; %Radius /Earth radii
pluto.trot = 6.387; %Rotational period /days
pluto.P = 248.348; %Orbital period /years
pluto.rgb = (1/255)*[149,107,107];

%End of code