One , Operation results

Two ,MATLAB program
% Simulation of the orbit of the four planets %Author: Lei Zhen %Date: 2019-04-26 % Universal gravitation formula F=GMm/r^2
% Object motion is the change of position caused by velocity and acceleration clc clear all close all set(0,'defaultfigurecolor','w'); %T
quality x coordinate y coordinate z coordinate x speed y speed z speed x acceleration y acceleration z acceleration Planet diameter % （ Closing speed , Not component velocity ） T=[ 1.9891e30 0 0
0 0 0 0 0 0 0 140e7 % sunlight 1 3.30e23 56.6720e9 0 0 0 47.88e3 0 0 0 0 4880000 % Mercury 2
4.869e24 108.2065e9 0 0 0 35.02e3 0 0 0 0 12103600 % Venus 3 5.98e24 149.5791e9 0 0
0 29.79e3 0 0 0 0 12756300 % earth 4 6.4219e23 226.9024e9 0 0 0 24.13e3 0 0 0 0
6794000 % Mars 5 ]; %%%%%%%%%%%%%%%%%%%2000 year 1 month 1 day 12 Time star heliocentric longitude
Position velocity conversion of celestial bodies %%%%%%%%%%%%%%%%%%%%%%% JD=[253.78494444 182.60408333 100.36890633
359.44775000].*pi/180; %2000 year 1 month 1 day 12 The Yellow meridian of the sun's heart JL=T(2:5,2); T(2:5,2)=JL.*cos(JD');
T(2:5,3)=JL.*sin(JD'); V=T(2:5,6); T(2:5,5)=V.*-sin(JD'); T(2:5,6)=V.*cos(JD');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
si=size(T); temp=10;% Drawing speed and precision control , The bigger, the faster, the more accurate time=365*24*temp*3;%3 circle t=3600/temp;
for k=1:time for i=1:si(1) % day=fix((k/60+12)/24)+rem(k/60+12,24)/24;
fx=0;fy=0;fz=0; for j=1:si(1) if i~=j
r=sqrt((T(i,2)-T(j,2))^2+(T(i,3)-T(j,3))^2+(T(i,4)-T(j,4))^2);%% Linear distance between two planets if
r>T(i,8)/2+T(j,8)/2 x=T(i,2)-T(j,2);% Two planets in x Distance in direction y=T(i,3)-T(j,3);% Two planets in y Distance in direction
z=T(i,4)-T(j,4);% Two planets in z Distance in direction f=6.67e-11*T(i,1)*T(j,1)/r^2;% The gravity of two planets
fx=fx-x/r*f;%x Direction component fy=fy-y/r*f;%y Direction component fz=fz-z/r*f;%z Direction component else
text(T(i,2),T(i,3),' Collision ') axis equal hold on end end end T(i,8)=fx/T(i,1);
%x Acceleration in direction T(i,9)=fy/T(i,1); %y Acceleration in direction T(i,10)=fz/T(i,1);%z Acceleration in direction end for
i=1:si(1) T(i,2)=T(i,2)+T(i,5)*t+0.5*T(i,8)*t^2; % Planetary x Axis coordinates
T(i,3)=T(i,3)+T(i,6)*t+0.5*T(i,9)*t^2;% Planetary y Axis coordinates
T(i,4)=T(i,4)+T(i,7)*t+0.5*T(i,10)*t^2;% Planetary z Axis coordinates T(i,5)=T(i,5)+T(i,8)*t;
% Planetary x Axial speed T(i,6)=T(i,6)+T(i,9)*t; % Planetary y Axial speed
T(i,7)=T(i,7)+T(i,10)*t;% Planetary z Axial speed end if rem(k,24*temp*5)==0 % each 5 Days drawing track , Single loop point control %
h_axes=findobj(gcf,'type','axes'); %% Get the handle of all coordinates in the current graph %
h_children_axes=allchild(h_axes); %% Handle to the child object that gets the coordinates % delete(h_children_axes); %%
Delete all coordinate sub object handles , To achieve the purpose of emptying figure(1),axis equal, axis([-2.5e11,2.5e11,-2.5e11,2.5e11]);
% Sun centric drawing , Demonstration of Rao sun revolution set(gcf,'position',[250,250,700,700]);
set(gca,'position',[0,0,1,1]);
plot(T(1,2),T(1,3),'or','MarkerFaceColor','r','MarkerSize',20); hold on,
plot(T(2,2),T(2,3),'ok','MarkerFaceColor','k','MarkerSize',5); hold on,
plot(T(3,2),T(3,3),'oy','MarkerFaceColor','y','MarkerSize',5); hold on,
plot(T(4,2),T(4,3),'ob','MarkerFaceColor','b','MarkerSize',5); hold on,
plot(T(5,2),T(5,3),'oc','MarkerFaceColor','c','MarkerSize',5); drawnow end end

Technology
Daily Recommendation
views 5