I have a code for simulation of earth sun and moon i need to incorporate mercury how would i do it please help.

function [] = SunEarthMoon(years,framerate)

% Inputs:
%   years     = Observation Time Period (years)
%   framerate = Video Speed (typically between 1 and 1000)

%%Clean Up
close all
clc

%%Initializaion

x_earth = 147300000000; % [m]

v_earth = 30257; % [m/s]

r_sat = 384748000; % earth surface [m]

r_earth = 6367000; % earth radius [m]

v_sat = 1023; % relative velocity from earth [m/s]

a = 5.145; % Angle to vertical (y) axis 

b = 90; % Angle to horizontal (x) axis in xz plane

x_earth_o = [x_earth; 0; 0];

x_sun_o = [0; 0; 0];

x_sat_o = [x_earth+r_sat; 0; 0];

v_earth_o = [0; v_earth; 0];

v_sun_o = [0; 0; 0];

v_sat_o = v_sat*[cos(pi/180*b)*sin(pi/180*a); cos(pi/180*a); sin(pi/180*b)*sin(pi/180*a)] + v_earth_o;

interval= years*[0 31536000];

%%Error Control

h = [0.01 36000];

tol = 100000;

Options.AbsTol = tol;

Options.MaxStep = h(2);

Options.InitialStep = h(1);

%%Analysis

ao = [x_earth_o; v_earth_o; x_sun_o; v_sun_o; x_sat_o; v_sat_o]; 

[t, x] = ode45(@earthfinal,interval,ao,Options);

for i = 1:length(t)

      R1(i) = (x(i,13)-x(i,1));

      R2(i) = (x(i,14)-x(i,2));

      R3(i) = (x(i,15)-x(i,3));

      R(i) = sqrt(R1(i)^2+R2(i)^2+R3(i)^2);

end

T_index_earth = find([1; x(:,4)].*[x(:,4);1]<=0);

T_index_moon = find([1; R2(:)].*[R2(:); 1]<=0);

for i = 4:length(T_index_earth)

      T_earth_semi(i-3) = (t(T_index_earth(i)-1)-t(T_index_earth(i-2)-1))/24/60/60;

end

T_earth = mean(T_earth_semi);

for i = 4:length(T_index_moon)

      T_moon_semi(i-3) = (t(T_index_moon(i)-1)-t(T_index_moon(i-2)-1))/24/60/60;

end

T_moon = mean(T_moon_semi);

D_earth = 0;

for i = 2:(T_index_earth(4)-1)

      D_earth = D_earth + sqrt((x(i,1)-x(i-1,1))^2+(x(i,2)-x(i-1,2))^2+(x(i,3)-x(i-1,3))^2);

end

D_moon = 0;

for i = 2:(T_index_moon(4)-1)

      D_moon = D_moon + sqrt((R1(i)-R1(i-1))^2+(R2(i)-R2(i-1))^2+(R3(i)-R3(i-1))^2);

end

%%Plots

q = framerate;

scrsz = get(0,'ScreenSize');

figure('position', [0.05*scrsz(3) 0.05*scrsz(4) 0.75*scrsz(3) 0.85*scrsz(4)])

set(gcf,'name','Sun, Earth, and Moon Orbits')

for i = 1:length(t)/q

      subplot(2,2,1)

plot3(x(1:i*q,1),x(1:i*q,2),x(1:i*q,3),'g',x(1:i*q,7),x(1:i*q,8),x(1:i*q,9),'r',x(1:i*q,13),x(1:i*q,14),x(1:
i*q,15),'b')

      axis(1.1*[min(x(:,1)) max(x(:,1)) min(x(:,2)) max(x(:,2)) 2*min(x(:,15)) 2*max(x(:,15))])

      xlabel('Universal X Coordinate (m)')

      ylabel('Universal Y Coordinate (m)')

      zlabel('Universal Z Coordinate (m)')

      title('Relative Orbits')

      legend('Earth','Sun','Moon')

      hold on

      plot3(x(i*q,1),x(i*q,2),x(i*q,3),'g-o',x(i*q,7),x(i*q,8),x(i*q,9),'r-
  o',x(i*q,13),x(i*q,14),x(i*q,15),'b-o')

      hold off

      subplot(2,2,2)

      plot3(R1(1:i*q),R2(1:i*q),R3(1:i*q),'b',zeros(1,i*q),zeros(1,i*q),zeros(1,i*q),'g')

      axis(1.5*[min(R1) max(R1) min(R2) max(R2) min(R3) max(R3)])

      xlabel('Universal X Coordinate (m)')

      ylabel('Universal Y Coordinate (m)')

      zlabel('Universal Z Coordinate (m)')

      title('Relative Moon Orbit About Earth')

      hold on

      plot3(R1(i*q),R2(i*q),R3(i*q),'b-o',0,0,0,'g-o')

      text(0,1.45*max(R2),1.40*max(R3),sprintf('Orbital Period, T = %3.5g days',T_moon))

      text(0,1.45*max(R2),1.15*max(R3),sprintf('Orbital Circumference, D = %3.5g gigameters',D_moon*1e-9))

      hold off

      subplot(2,2,3)

      plot(x(1:i*q,1),x(1:i*q,2),'g',x(1:i*q,7),x(1:i*q,8),'r')

      axis(1.5*[min(x(:,1)) max(x(:,1)) min(x(:,2)) max(x(:,2))])

      xlabel('Universal X Coordinate (m)')

      ylabel('Universal Y Coordinate (m)')

      title('Relative Earth Orbit About Sun')

      hold on

      plot(x(i*q,1),x(i*q,2),'g-o',x(i*q,7),x(i*q,8),'r-o')

      text(1.45*min(x(:,1)),1.40*max(x(:,2)),sprintf('Orbital Period, T = %3.5g days',T_earth))

      text(1.45*min(x(:,1)),1.25*max(x(:,2)),sprintf('Orbital Circumference, D = %3.5g gigameters',D_earth*1e-
  9))

      text(1.45*min(x(:,1)),1.40*min(x(:,2)),sprintf('Time, t = %3.3g days',round(t(q*i)/24/60/60)))

      hold off

      subplot(2,2,4)

      plot(t(1:i*q)/(60*60*24),R(1:i*q)/1000,'b')

      axis([t(1)/24/60/60 t(end)/24/60/60 0.999*min(R)/1000 1.001*max(R)/1000])

      xlabel('Time,t (days)')

      ylabel('Orbit Radius, R (km)')

      title('Moon-Earth Distance')

      hold on

      plot(t(i*q)/(60*60*24),R(i*q)/1000,'b-o')

      hold off

      drawnow

end

end

%%Differential Equation Function

function [udot]= earthfinal(t,u)

m_earth = 5.9742e24; % [kg]

m_sun = 1.98892e30; % [kg]

m_sat = 11110; % [kg]

G = 6.67300e-11; %[(m)^3(kg)^-1(s)^-2];

d_earth_sun = sqrt((u( 7,1)-u(1,1))^2+(u( 8,1)-u(2,1))^2+(u( 9,1)-u(3,1))^2);

d_earth_sat = sqrt((u(13,1)-u(1,1))^2+(u(14,1)-u(2,1))^2+(u(15,1)-u(3,1))^2);

d_sun_sat =   sqrt((u(13,1)-u(7,1))^2+(u(14,1)-u(8,1))^2+(u(15,1)-u(9,1))^2);

% Earth motion

udot( 1,1) = u(4,1);

udot( 2,1) = u(5,1);

udot( 3,1) = u(6,1);

udot( 4,1) = m_sun*G*(u(7,1)-u(1,1))/d_earth_sun^3 + m_sat*G*(u(13,1)-u(1,1))/d_earth_sat^3;

udot( 5,1) = m_sun*G*(u(8,1)-u(2,1))/d_earth_sun^3 + m_sat*G*(u(14,1)-u(2,1))/d_earth_sat^3;

udot( 6,1) = m_sun*G*(u(9,1)-u(3,1))/d_earth_sun^3 + m_sat*G*(u(15,1)-u(3,1))/d_earth_sat^3;

% Sun Motion

udot( 7,1) = u(10,1);

udot( 8,1) = u(11,1);

udot( 9,1) = u(12,1);

udot(10,1) = m_earth*G*(u(1,1)-u(7,1))/d_earth_sun^3 + m_sat*G*(u(13,1)-u(7,1))/d_sun_sat^3;

udot(11,1) = m_earth*G*(u(2,1)-u(8,1))/d_earth_sun^3 + m_sat*G*(u(14,1)-u(8,1))/d_sun_sat^3;

udot(12,1) = m_earth*G*(u(3,1)-u(9,1))/d_earth_sun^3 + m_sat*G*(u(15,1)-u(9,1))/d_sun_sat^3;

% Satellite Motion

udot(13,1) = u(16,1);

udot(14,1) = u(17,1);

udot(15,1) = u(18,1);

udot(16,1) = m_earth*G*(u(1,1)-u(13,1))/d_earth_sat^3 + m_sun*G*(u(7,1)-u(13,1))/d_sun_sat^3;

udot(17,1) = m_earth*G*(u(2,1)-u(14,1))/d_earth_sat^3 + m_sun*G*(u(8,1)-u(14,1))/d_sun_sat^3;

udot(18,1) = m_earth*G*(u(3,1)-u(15,1))/d_earth_sat^3 + m_sun*G*(u(9,1)-u(15,1))/d_sun_sat^3;

end

[John D'Errico]

It is a small planet. Hardly worth the effort. But you have already written code for a 3 body problem. So a 4th body is just a few more interactions. You do it the same way you did the others.

Or, maybe you got this code from someone else, and have no idea how it was written. In that case, you would need to spend some serious time hacking their code, because you will need to figure it out. But that is no different from what you are asking us to do, to spend a few hours figuring out exactly what to do with poorly documented code that we have never seen before. So why not sit down and start doing the work?

[Arslan Ahmad]

you got me bro

[rowrowyourboat]

Were you able to incorporate Mercury's orbit? I have the same problem - thanks

Answer 1

[Conor Kelly]

I have the same code but for some reason it won't compile. What version of matlab is this in? Thank you.

[Walter Roberson]

That code is for MATLAB 5.1 or later.

A small number of the posted lines were wrapped onto the next line and that needed to be fixed.