%solar_system_spirograph
% Plot lines between planet orbital positions (in x,y plane of ecliptic)
% at equal time intervals.
% Inspired by https://engaging-data.com/planetary-spirograph/.
%
% LAST UPDATED by Andy French August 2022

function solar_system_planet_spirograph(varargin)

%Planet default inputs
if nargin==0
    planet1 = 'earth' ; planet2 = 'mars';
elseif nargin==2
    planet1 = varargin{1}; planet2 = varargin{2};
else
    disp('Incorrect number of inputs. Planets chosen are Earth and Mars');
    planet1 = 'earth' ; planet2 = 'mars';
end
show_axis = 0;

%Get planet data
p1 = get_planet_data( planet1 ); p2 = get_planet_data( planet2 );

%Define planet position timings
tmax = 10*max([p1.P,p2.P]); %Max time /years
dt = tmax/1234;  %Time interval in years

%Define time array /years
t = 0 : dt : tmax;

%Create spirograph figure
figure('name','solar system','color',[1 1 1]);
axes('nextplot','add','fontsize',14);

%Compute orbits
[x1,y1,z1,theta1,p1] = calc_orbit( planet1, t, '3D' );
[x2,y2,z2,theta2,p2] = calc_orbit( planet2, t, '3D' );

%Loop through time t and plot line joining planets
for n=1:length(t)
    plot([x1(n),x2(n)],[y1(n),y2(n)],'color','k','linewidth',0.1);
    plot(x1(n),y1(n),'.','color',p1.rgb);
    plot(x2(n),y2(n),'.','color',p2.rgb);
end
axis equal; grid on; box on;
xlabel('x /AU'); ylabel('y /AU'); title( [planet1,' ',planet2,' spirograph'] );
if show_axis==0; axis off; end
print( gcf,[planet1,' ',planet2,' spirograph.png'],...
    '-dpng','-r300'); close(gcf);

%%

%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

%%

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

%Angle step for Simpson's rule
dtheta = (1e-6)*2*pi*max(t)/P;

%

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

%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
function planet = get_planet_data( planet_name )
[sun,mercury,venus,earth,mars,...
    jupiter,saturn,uranus,neptune,pluto] =...
    solar_system_parameters;
eval([ 'planet = ',planet_name,';'] );

%%

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

%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*rand;  %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*rand;  %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*rand;  %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];

%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*rand;  %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*rand;  %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];

%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*rand;  %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*rand;  %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*rand;  %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*rand;  %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];

%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^-3

%End of code