• Introduction
• Section I: Coordinates Transformation
• Section 2: Orbits in the universal famework
• Section 3: Plot the orbits
• Section 4: Export to VRML
• References
• Licence Agreement

A Virtual 3D Solar System

by Yi Xie,
Astronomy Department of Najing University, China,
E-mail: shea@citiz.net ,
20/01/2004 Yi Xie

NOTE: This application demonstrates how Maple 9 can be used in making a virtual 3D solar system.

Introduction

In solar system, all planetary orbits, except Pluto 's, are nearly in a plane and circular. Here, we use a little knowledge of Celestial Mechanics and power of Maple to make a virtual 3D solar system and illustrate these characteristics.

>	restart; with(linalg): with(plots): with(plottools):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Section I: Coordinates Transformation

 In 2-body problem, such as Sun and a planet, it exists 6 orbital elements that can describe planetary orbital configuration and motion of  the planet. The orbital elements are semi-major axis (a), eccentricity (e), orbital inclination (i), longitude of ascending node ( ), augument of perihelion ( ) and mean anomaly (M). And in 2-body problem, we can easily find that the orbit of planet is in a 2D invariant plane, so  each planet has a different invariant plane. To get an outline picture of solar system, we should convert these invariant planes to a universal framework through coordinates transformation by matrix rotation. See reference for details about 2-body problem and coordinates transformation.

For jth planet, we define 3 matrixes ,  and  for rotation.

>

P1:=matrix([[cos(omega[j]),-sin(omega[j]),0],             [sin(omega[j]), cos(omega[j]),0],             [0,0,1]]); P2:=matrix([[1,0,0],             [0,cos(i[j]),-sin(i[j])],             [0,sin(i[j]), cos(i[j])]]); P3:=matrix([[cos(Omega[j]),-sin(Omega[j]),0],             [sin(Omega[j]), cos(Omega[j]),0],             [0,0,1]]);

For jth planet,  the position vector in invariant plane is A , A =( , , 0), where  is eccentric anomaly and the position vector B in the universal famework is , where .

>	A:=matrix([[a[j]*(cos(E[j])-e[j])],[a[j]*sqrt(1-e[j]^2)*sin(E[j])],[0]]);

>	R:=multiply(P3,P2,P1);

>	B:=multiply(R,A);

Section 2: Orbits in the universal famework

Now, we have got the position vector in the universal famework, which only depends on  for other elements are known. It means we have got the parameter equations of orbits.

Here, we set 9 planetary orbital elements in 5 sequences. Note: the sequence [9] is  longitude of perihelion, instead of augument of perihelion. But we can convert them to augument of perihelion through some spherical trigonometry knowledge.

>

a:=[0.38709893,0.72333199,1.00000011,1.52366231,5.20336301,9.53707032,19.19126393,30.06896348,39.48168677]: e:=[0.20563069,0.00677323,0.01671022,0.09341233,0.04839266,0.05415060,0.04716771,0.00858587,0.24880766]: i:=[7.00487,3.39471,0.00005,1.85061,1.30530,2.48446,0.76986,1.76917,17.14175]: Omega:=[48.33167,76.68069,-11.26064,49.57854,100.55615,113.71504,74.22988,131.72169,110.30347]:omega:=[77.45645,131.53298,102.94719,336.04084,14.75385,92.43194,170.96424,44.97135,224.06676]:

Convert degree to radian.

>	i:=map(x->convert(x,units,deg,rad),i): Omega:=map(x->convert(x,units,deg,rad),Omega): omega:=map(x->convert(x,units,deg,rad),omega):

Calculate augument of perihelion from  longitude of perihelion.

>	for j to 9 do omega[j]:=arcsin(sin(omega[j]-Omega[j])/sin(arccos(sin(i[j])*cos(omega[j]-Omega[j])))): od:

Parameter equations of 9 planetary orbits.

>	x:=array(1..9): y:=array(1..9): z:=array(1..9):

>	for j to 9 do x[j]:=B[1,1]: y[j]:=B[2,1]: z[j]:=B[3,1]: od:

Section 3: Plot the orbits

As we have got the orbit equations, we can plot them directly.

>	pic:=array(1..19):

Plot the orbits.

>	for j to 9 do  pic[j]:=spacecurve([subs(E[j]=E,x[j]),subs(E[j]=E,y[j]),subs(E[j]=E,z[j])],                     E=0..2*Pi,color=COLOR(HUE,j/10.)): od:

Plot the Sun.

>	pic[10]:=plot3d(0.05,0..2*Pi,0..Pi,style=PATCHNOGRID,coords=spherical,color=red):

Texts on each orbit.

>

pic[11]:=textplot3d([subs(E[1]=0,x[1]),subs(E[1]=0,y[1]),subs(E[1]=0,z[1]),"Mercury"]): pic[12]:=textplot3d([subs(E[2]=0,x[2]),subs(E[2]=0,y[2]),subs(E[2]=0,z[2]),"Venus"]): pic[13]:=textplot3d([subs(E[3]=0,x[3]),subs(E[3]=0,y[3]),subs(E[3]=0,z[3]),"Earth"]): pic[14]:=textplot3d([subs(E[4]=0,x[4]),subs(E[4]=0,y[4]),subs(E[4]=0,z[4]),"Mars"]): pic[15]:=textplot3d([subs(E[5]=0,x[5]),subs(E[5]=0,y[5]),subs(E[5]=0,z[5]),"Jupiter"]): pic[16]:=textplot3d([subs(E[6]=0,x[6]),subs(E[6]=0,y[6]),subs(E[6]=0,z[6]),"Saturn"]): pic[17]:=textplot3d([subs(E[7]=0,x[7]),subs(E[7]=0,y[7]),subs(E[7]=0,z[7]),"Uranus"]): pic[18]:=textplot3d([subs(E[8]=0,x[8]),subs(E[8]=0,y[8]),subs(E[8]=0,z[8]),"Neptune"]): pic[19]:=textplot3d([subs(E[9]=0,x[9]),subs(E[9]=0,y[9]),subs(E[9]=0,z[9]),"Pluto"]):

>	pic1:=convert(pic,'list'):

display(pic1,scaling=CONSTRAINED);

>	fig:=display(pic1,scaling=CONSTRAINED):

Section 4: Export to VRML

For semi-major axises vary largely from about 0.4A.U. to nearly 40A.U., so it is difficult to distinguish inner planets in a whole picture of solar system as above. But VRML can solve this problem and it can show some characteristics mentioned at first. In exported VRML file, the color is not good. Of course, you can modify it.
You can visit my VRML about a 3D solar system,   http://www.hjsm.net/astro/sheaweb/files/9planets.wrl .

>	#vrml(fig,"9planets.wrl");

>

References

Murray, C. D. and Dermott, S. F. 1999: "Solar System Dynamics", Cambridge University Press, Cambridge.

Licence Agreement

20/01/2004
Yi Xie