MATLAB行星运动仿真程序

clear all;close all;G=6.67*10^-11;% The universal gravitational constantm = [1.989e30,3.5844e23,4.89868e24,5.974e24,6.5714e23,1.89854e27,5.68725e26,8.72204e25,1.02753e26];% An array of massesn = length(m); Number_of_planets = n% The number of massesx_p = [0,58340100000,1.07705e11,1.4959e11,2.27377e11,7.77868e11,1.42709e12,2.87512e12,4.49668e12];% An array of x positionsy_p= [0,0,0,0,0,0,0,0,0];% An array of y positionsx_v = [0,0,0,0,0,0,0,0,0];% An array of x velocitiesy_v = [0,47856.46,34961.72,29780,23883.56,12924.52,9618.94,6789.84,5419.96];% An array of y velocitiesT= 365*24*60*60;%(2*pi*G*m(1))/(x_v(1)^3);dt=360*60*60*10;colordef black;ph = plot(x_p,y_p,'.','MarkerSize',30);xlabel('distance(km)');ylabel('distance(km)');title('Planetary motion');for a=0:dt:1000*T%loop for all masses with respect to timefor k=1:n%loop for individual masses calculating neccessary quantitiesdx = x_p - x_p(k);% difference in x positionsdy = y_p - y_p(k);% difference in y positionsmag = (dx.^2 + dy.^2).^0.5;% magnitude of the distance between the 2 massesf = ((G*m(k)*m)./(mag.^2));% total forcefx_old =(f.*dx./mag);fx_old(k)= 0;fx = sum (fx_old);%Summing the total force in x direction for each planetfy_old =(f.*dy./mag);fy_old(k)= 0;fy = sum (fy_old);%Summing the total force in y direction for each planeta_x = fx/m(k);% calculating acceleration change in x directiona_y = fy/m(k);% calculating acceleration change in y directionx_v(k) = x_v(k) + dt*a_x;% calculating velocity change in x directiony_v(k) = y_v(k) + dt*a_y;% calculating velocity change in y directionx_p_new(k)= x_p(k)+dt*x_v(k);% calculating new x positionsy_p_new(k)= y_p(k)+dt*y_v(k);% calculating new y positionsendx_p=x_p_new;y_p=y_p_new;% overwriting old positions with new positionspause(0.1)plot(x_p,y_p,'.','MarkerSize',30)axis([-(10^15) 10^15 -(10^15) 10^15])xlabel('distance(km)');ylabel('distance(km)');title('Planetary motion');drawnow% Plotting graphend