จำลองการเคลื่อนที่ของดาวเคราะห์ 1 ดวง โดยมีตย.

[Spoil] คลิกเพื่อดูข้อความที่ซ่อนไว้โปรแกรมตัวอย่างสำหรับการจำลองวงโคจรของโลก
function  
PC417_Kepler_main ;

% Note : G*Ms = 4*pi^2
% v1 = sqrt ([G*Ms*(1-e)]/[a*(1+e)])

mov=avifile('Kepler.avi');
mov.quality=100;
mov.fps=30;
N=10000 ;
dt=1/365 ; % define time

t=zeros(1,N); % store time
x=zeros(1,N); % x position
y=zeros(1,N);
z=zeros(1,N);
vx=zeros(1,N);
vy=zeros(1,N);

%data= zeros (5,N) ;

t(1)= 0;
x(1)= 1;
y(1)= 0;
vx(1)=0;
vy(1)=2*pi;
for i=1:N-1
    [t(i+1) x(i+1) y(i+1) vx(i+1) vy(i+1)]=PC417_Kepler_compute(dt,t(i),x(i),y(i),vx(i),vy(i));

end


[SX SY SZ]=sphere (20);
SX=SX*0.125;
SY=SY*0.125;
SZ=SZ*0.125;

for i=1:400
    fig=figure (33) ; clf ;
    set(fig,'DoubleBuffer','on');
    
    sun=surf(SX*2, SY*2, SZ*2);
    set (sun,'FaceColor','r');
    hold on;

star=surf ( SX+x(i),SY+y(i),SZ );
    set (star,'FaceColor','b');
    plot3 (x(1:i),y(1:i),z(1:i));
    
    axis equal ;
    
    axis ([-1.5 1.5 -1.5 1.5 -0.75 0.75]);
    set (fig,'Position',[200 200 320 240]);
    set(fig,'Color',([1 1 1]));
    ff=getframe(fig);

mov=addframe(mov,ff);
    %set (gca,'Position',[200 200 400 400]);
end

kk=close(mov);

disp ('*');

function  [tj xj yj vxj vyj]=PC417_Kepler_compute (dt,ti,xi,yi,vxi,vyi);
r=sqrt(xi^2+yi^2);
C0=4*pi^2/(r^3);

tj=ti+dt;
vxj=vxi-C0*xi*dt;
vyj=vyi-C0*yi*dt;

xj=xi+vxj*dt;
yj=yi+vyj*dt;

แต่ลองเอาไปใส่ที่โปรแกรมแล้วรันไม่ขึ้นค่ะ
โดยโค้ดข้างบนจะขึ้นการจำลองการเคลื่อนที่และแสดงวงโคจรรอบดวงอาทิตย์