Chapter 8, Section 8.17, Figure 8.30

We start by inputting the orbital radii of Io, Europa and Ganymede in units of the planetary radius.  The data are taken from

r1=421600./71398.;
r2=670900./71398.;
r3=1070000./71398.;

Set up the two ratios of mean motions.

n12=0.5;
n23=n12/2;

Define functions to calculate the position vector of each satellite as a function of time.

pos1[t_]:=(x=r1*Cos[t+Pi]; y=r1*Sin[t+Pi]; {x,y})
pos2[t_]:=(x=r2*Cos[n12*t]; y=r2*Sin[n12*t]; {x,y})
pos3[t_]:=(x=r3*Cos[n23*t]; y=r3*Sin[n23*t]; {x,y})

Define a Circle for the orbital path of each satellite ...

orb1=Circle[{0,0},r1];
orb2=Circle[{0,0},r2];
orb3=Circle[{0,0},r3];

.. and the planet Jupiter (with unit radius).

jup=Circle[{0,0},1];

Define a function Satellites to do all the plots together for a particular time.

Satellites[t_] :=
Show[Graphics[{jup,orb1, orb2,orb3,
               AbsolutePointSize[7],
               Disk[pos1[t],1],
               Disk[pos2[t],1],
               Disk[pos3[t],1]},
               AspectRatio->Automatic,
               PlotRange->{{-25,25},{-18,18}}]]

Define the specific times we are interested in.  These times give the whole sequence of conjunctions between pairs of satellites.

time = {0, 4*Pi/3, 2Pi, 4Pi, 6Pi, 6Pi+2Pi/3}

Plot the orbital positions at each time, putting each plot in an array element.

Do[plot[i]=Satellites[time[[i]]],{i,6}]

Display all six plots together.

Show[GraphicsArray[
{{plot[1],plot[2],plot[3]},{plot[4],plot[5],plot[6]}}]]