Как сделать траекторию движения планет?

uses crt, graphabc;
type
   info=record
   m,x,dy,y,dx,old:real;
   title:string;
   rad,cl:integer;
   end;
const
 dt=500;
 G=6.67e-11;
var
  xx,yy,cosA,sinA,zoom,r,ax,ay,t : real;
  x0,y0: integer;
  mas: array [1..4] of info;
  i,j:byte;
BEGIN
  SetWindowSize(640,480);
  SetBrushColor(clBlack);
    FillRect(0,0,640,480);
 x0:=WindowWidth div 2;
 y0:=WindowHeight div 2;
   zoom:=2E-10;
  {Сoлнце}
 with mas[1] do
  begin
   title:='Сoлнце';
   m:=1.989e+30;
   x:=0;
   old:=0.1;
   y:=0;
   dx:=0;
   dy:=0;
   rad:=8;
   cl:=(clYellow);
 with mas[2] do
   begin
    title:='Меркурий';
    m:=3.303e+23;
    x:=57.85E+9;
    y:=0;
    old:=0.1;
    dx:=0;
    dy:=4.787E+4;
    cl:=(clPurple);
    rad:=4;
   end;
{Венера}
 with mas[3] do
  begin
   title:='Венера';
   m:=4.869e+24;
   x:=108.11E+9;
   y:=0;
   dx:=0;
   old:=0.1;
   dy:=3.502E+4;
   cl:=(clLime);
   rad:=4;
    end;
{Земля}
 with mas[4] do
  begin
   title:='Земля';
   m:=5.974e+24;
   x:=149.46E+9;
   y:=0;
   dx:=0;
   old:=0.1;
   dy:=2.978E+4;
   cl:=(clBlue);
      rad:=4;
    end;
      Repeat
  t:=t+dt;
  for i:=1 to 4 do
  begin
   ax:=0;
   ay:=0;
   for j:=1 to 4 do
    if i<>j then
     begin
      R:=sqrt(sqr(mas[i].x-mas[j].x)+sqr(mas[i].y-mas[j].y));
      cosA:=-(mas[i].x-mas[j].x)/(R);
      sinA:=-(mas[i].y-mas[j].y)/(R);
      ax:=ax+G*mas[j].m*cosA/sqr(R);
      ay:=ay+G*mas[j].m*sinA/sqr(R);
     end;
   mas[i].dx:=mas[i].dx+dt*ax;
   mas[i].dy:=mas[i].dy+dt*ay;
   mas[i].x:=mas[i].x+mas[i].dx*dt;
   mas[i].y:=mas[i].y+mas[i].dy*dt;
   lockdrawing;
   setbrushcolor(mas[i].cl);
   Circle(round(mas[i].x* zoom)+x0,round(mas[i].y* zoom)+y0,mas[i].rad);
    Redraw;
    end;
  until KeyPressed;
    end;
     end.