Chapter 2, Section 2.8, Figure 2.14

We start by a function r using Eq.(2.20) to give us the radial distance for given values of semi-major axis, a, eccentricity, e, and true anomaly, f.

r[a_,e_,f_]:=a(1-e^2)/(1+e Cos[f]);

Now define a fixed viewpoint, vp, from which to view the 3D plots.

vp = {1.300, -2.400, 2.000};

Define the first transformation matrix, P1.  This rotates an angle ω about the z-axis. (See Eq.(2.119).)

P1={{Cos[ω],-Sin[ω],0},{Sin[ω],Cos[ω],0},{0,0,1}};

Display P1 as a matrix.

MatrixForm[P1]

Define the second transformation matrix, P2.  This rotates an angle i about the x-axis. (See Eq.(2.119).)

P2={{1,0,0},{0,Cos[i],-Sin[i]},{0,Sin[i],Cos[i]}};

Display P2 as a matrix.

MatrixForm[P2]

Define the third transformation matrix, P3.  This rotates an angle Ω about the z-axis. (See Eq.(2.120).)

P3={{Cos[Ω],-Sin[Ω],0},{Sin[Ω],Cos[Ω],0},{0,0,1}};

Display P3 as a matrix.

MatrixForm[P3]

Now multiply these together to get the complete rotation matrix, P, as required.  (See Eq.(2.122).)

P=P3.P2.P1;

Now define the values of each of the rotation angles.

Ωtry=40°; itry=60°; ωtry=70°;

Define the three original axes ax, ay, az (of length 1.8) and the location of the origin, orig.

ax={1.8,0,0}; ay={0,1.8,0}; az={0,0,1.8}; orig={0,0,0};

Calculate the locations of the transformed axes X, Y, Z by applying the matrix P.

X=P.ax/.{Ω->Ωtry,i->itry,ω->ωtry};
Y=P.ay/.{Ω->Ωtry,i->itry,ω->ωtry};
Z=P.az/.{Ω->Ωtry,i->itry,ω->ωtry};

Now define graphisc objects for the final axes ...

OX=Graphics3D[{Line[{orig,X}]}];
OY=Graphics3D[{Line[{orig,Y}]}];
OZ=Graphics3D[{Line[{orig,Z}]}];

... and the original axes.

aX=Graphics3D[{Line[{orig,ax}]}];
aY=Graphics3D[{Line[{orig,ay}]}];
aZ=Graphics3D[{Line[{orig,az}]}];

Put these together in two graphics objects finalaxes and originalaxes.

finalaxes={aX,aY,aZ}; originalaxes={OX,OY,OZ};

Create a graphics object to denote the ecliptic plane.

ecliptic=
Graphics3D[
{GrayLevel[1],Polygon[{{1.6,1.6,-0.001},{-1.6,1.6,-0.001},
{-1.6,-1.6,-0.001},{1.6,-1.6,-0.001}}]}];

Check what the final, transformed axes are going to look like by creating an ellipse with a=1.0 and e=0.5 ...

ell3D=
(xyz=P.{r[a,e,t] Cos[t], r[a,e,t] Sin[t], 0}
/.{i->60°,ω->70°,Ω->40°,a->1.0,e->0.5};
Graphics3D[
{GrayLevel[0.9],Polygon[Table[N[xyz],{t,0,2Pi,Pi/60}]]}]);

... and plotting this ellipse together with the two sets of axes and the ecliptic plane.

Show[ell3D,originalaxes,finalaxes,ecliptic,
 Boxed->False,SphericalRegion->True,
 Lighting->False,
 ViewPoint->vp
]

Create the ellipse for the starting position with rotations set to zero.

ell3D0=
(xyz=P.{r[a,e,t] Cos[t], r[a,e,t] Sin[t], 0}
/.{i->0°,ω->0°,Ω->0°,a->1.0,e->0.5};
Graphics3D[
{GrayLevel[0.9],Polygon[Table[N[xyz],{t,0,2Pi,Pi/60}]]}]);

Show this ellipse together with both sets of axes (coincident in this case).

p1=
Show[ell3D0,finalaxes,ecliptic,
 Boxed->False,SphericalRegion->True,
 PlotRange->{{-1.8,1.8},{-1.8,1.8},{-1.8,1.8}},
 Lighting->False,
 ViewPoint->vp
]

Create the ellipse for the positions after the first rotation.

ell3D1=
(xyz=P.{r[a,e,t] Cos[t], r[a,e,t] Sin[t], 0}
/.{i->0°,ω->40°,Ω->0°,a->1.0,e->0.5};
Graphics3D[
{GrayLevel[0.9],Polygon[Table[N[xyz],{t,0,2Pi,Pi/60}]]}]);

Work out the locations of the new axes and create graphics objects for them.

ax1=P1.ax/.ω->70°;
ay1=P1.ay/.ω->70°;
az1=P1.az/.ω->70°;
aX1=Graphics3D[{Line[{orig,ax1}]}];
aY1=Graphics3D[{Line[{orig,ay1}]}];
aZ1=Graphics3D[{Line[{orig,az1}]}];
axes1={aX1,aY1,aZ1};

Show the ellipse together with both sets of axes.

p2 = 
Show[ell3D1,axes1,finalaxes,ecliptic,
 Boxed->False,SphericalRegion->True,
 PlotRange->{{-1.8,1.8},{-1.8,1.8},{-1.8,1.8}},
 Lighting->False,
 ViewPoint->vp
]

Create the ellipse for the positions after the second rotation.

ell3D2=
(xyz=P.{r[a,e,t] Cos[t], r[a,e,t] Sin[t], 0}
/.{i->60°,ω->70°,Ω->0°,a->1.0,e->0.5};
Graphics3D[
{GrayLevel[0.9],Polygon[Table[N[xyz],{t,0,2Pi,Pi/60}]]}]);

Work out the locations of the new axes and create graphics objects for them.

ax2=P2.ax1/.i->60°;
ay2=P2.ay1/.i->60°;
az2=P2.az1/.i->60°;
aX2=Graphics3D[{Line[{orig,ax2}]}];
aY2=Graphics3D[{Line[{orig,ay2}]}];
aZ2=Graphics3D[{Line[{orig,az2}]}];
axes2={aX2,aY2,aZ2};

Show the ellipse together with both sets of axes.

p3 = 
Show[ell3D2,axes2,finalaxes,ecliptic,
 Boxed->False,SphericalRegion->True,
 PlotRange->{{-1.8,1.8},{-1.8,1.8},{-1.8,1.8}},
 Lighting->False,
 ViewPoint->vp
]

Create the ellipse for the positions after the third rotation.

ell3D3=
(xyz=P.{r[a,e,t] Cos[t], r[a,e,t] Sin[t], 0}
/.{i->60°,ω->70°,Ω->40°,a->1.0,e->0.5};
Graphics3D[
{GrayLevel[0.9],Polygon[Table[N[xyz],{t,0,2Pi,Pi/60}]]}]);

Work out the locations of the new axes and create graphics objects for them.

ax3=P3.ax2/.Ω->40°;
ay3=P3.ay2/.Ω->40°;
az3=P3.az2/.Ω->40°;
aX3=Graphics3D[{Line[{orig,ax3}]}];
aY3=Graphics3D[{Line[{orig,ay3}]}];
aZ3=Graphics3D[{Line[{orig,az3}]}];
axes3={aX3,aY3,aZ3};

Show the ellipse together with both sets of axes.

p4 = 
Show[ell3D3,axes3,finalaxes,ecliptic,
 Boxed->False,SphericalRegion->True,
 PlotRange->{{-1.8,1.8},{-1.8,1.8},{-1.8,1.8}},
 Lighting->False,
 ViewPoint->vp
]

Display all the four stages as a graphics array.

Show[GraphicsArray[
  {{p1,p2},{p3,p4}}],
  GraphicsSpacing->0.001]