pro solsystem
   readmatrix,aa,'sol_system.dat',o=5
   ss=size(aa)
   plot,[0],[0],psym=1,xr=[-40,40],yr=[-40,40],xst=1,yst=1,xtit='x [au]',ytit='y [au]'
;   plot,[0],[0],psym=1,xr=[-2,2],yr=[-2,2],xst=1,yst=1,xtit='x [au]',ytit='y [au]'
 
   FOR i=(ss(2)-1)/8-1,0,-1 DO BEGIN
     print,i
     oplot,aa(*,i*8+1),aa(*,i*8+2),psym=3,col=mycolor(i)
   END
   stop
end