def drawangle(expr A,X,B,r,w,col) =    draw (X + r*unitvector(A-X))         { r*unitvector(A-X) rotated 90 }         ..         { r*unitvector(B-X) rotated 90 }         (X + r*unitvector(B-X)) withpen pencircle scaled w withcolor col;enddef;def drawcircle(expr X,r,w,col,shade,angle) =    draw (X + (r,0))..(X + (0,r))..(X - (r,0))..(X - (0,r))..cycle withpen pencircle scaled w withcolor col;    if shade: fill (X + (r,0))--(X - (r,0))..(X - (0,r))..(X + (r,0))..cycle rotatedaround (X,angle); fi;enddef;def drawcircleS(expr X,r,w,col,angle) =    drawcircle(X,r,w,col,true,angle);enddef;s = 2cm;r = 2.2mm;path p;pair A, B, C, D;pair GEE;  % greatest eastern elongationpair GWE;  % greatest western elongationpair EQ;  % eastern quadraturepair WQ;  % western quadraturepair Sun;  % position of the Sun, defines the origin  Sun = (0,0);pair Earth;  % position of the Earthpair planet;numeric GEEangle, GWEangle, EQangle, WQangle, planetangle;GEEangle = 300;GWEangle = 60;EQangle = 311.8103;WQangle = 48.1897;  % arccos(2/3)pair Conj;  % conjunctionpair Opp;  % oppositionpair SupConj;  % superior conjunctionpair InfConj;  % inferior conjunction% draw the circles and define points of interestfor i=0 upto 2:  A := (s + i*s,0);  B := A rotated 90;  C := B rotated 90;  D := C rotated 90;  p := A..B..C..D..cycle;  draw p withpen pencircle scaled 1pt;    if i=1:    Earth := D;  fi;    if i=0:    GEE := D rotated GEEangle;    GWE := D rotated GWEangle;  fi;    if i=2:    EQ := D rotated EQangle;    WQ := D rotated WQangle;  % arccos(2/3)  fi;    % naming special points for further labeling references  if i=0:    SupConj := B;    InfConj := D;  fi;  if i=2:    Conj := B;    Opp := D;  fi;endfor;% draw small circles at points of interestdrawcircleS(Conj,r,1pt,black,180);drawcircleS(SupConj,r,1pt,black,180);drawcircleS(Opp,r,1pt,black,0);drawcircleS(InfConj,r,1pt,black,0);drawcircle(Earth,2.5mm,1pt,black,false,0);drawcircleS(GEE,r,1pt,black,GEEangle);drawcircleS(GWE,r,1pt,black,GWEangle);drawcircleS(EQ,r,1pt,black,EQangle);drawcircleS(WQ,r,1pt,black,WQangle);drawcircle(Sun,1.6mm,1pt,black,false,0);pair out; out = (0,0.6cm);pair inlft; inlft = (0,2.5mm) rotated 15;pair inrt; inrt = (0,2.5mm) rotated (-15);path q;for i=0 step 30 until 360:  q := (inlft--out--inrt) rotated i;  draw q;endfor;numeric planetAngle; planetAngle = 330;planet := D rotated planetAngle;drawcircleS(planet,r,1pt,black,planetAngle);% draw the dashed linesdraw B--D withpen pencircle scaled 1pt dashed evenly;draw EQ--WQ withpen pencircle scaled 1pt dashed evenly;draw GEE--Earth--GWE withpen pencircle scaled 1pt dashed evenly;draw Sun--planet--Earth withpen pencircle scaled 1pt dashed evenly;% draw the Sun and the Earthdrawangle(Earth,planet,Sun,1.0cm,1pt,black);drawangle(Sun,Earth,planet,0.7cm,1pt,black);% labellinglabel(btex $\alpha$ etex, planet) shifted (1cm,0.8cm);label(btex $\varepsilon$ etex, Earth) shifted (-0.8cm,0.6cm);lbl = 0.2cm;label.lrt("Earth", Earth) shifted (lbl,-lbl);label.urt("Sun", Sun) shifted (lbl,3*lbl);label.bot("Opposition", Opp) shifted (0,-lbl);label.top("Conjunction", Conj) shifted (0,lbl);label.top("Superior conjunction", SupConj) shifted (0,lbl);label.top("Inferior", InfConj) shifted (0,3*lbl);label.top("conjunction", InfConj) shifted (0,lbl);label.lrt("Western quadrature", WQ) shifted (0,-lbl);label.llft("Eastern quadrature", EQ) shifted (0,-lbl);label.top("Greatest western", GWE) shifted (8*lbl,lbl);label.top("elongation", GWE) shifted (6*lbl,-lbl);label.top("Greatest eastern", GEE) shifted (-8*lbl,lbl);label.top("elongation", GEE) shifted (-6*lbl,-lbl);% arrowspair OPO,EO,IPO;OPO := Conj rotated 32;EO := Earth rotated (180+42);IPO := SupConj rotated 80;pair sft; sft := (-lbl,lbl);drawarrow ((OPO + 3*sft)--(OPO + sft));drawarrow ((EO + 10*sft)--(EO + sft));drawarrow ((IPO + 18*sft)--(IPO + sft));label.top("Outer planet's", (OPO + 3*sft)) shifted (3*lbl,1.8*lbl);label.top("orbit", (OPO + 3*sft)) shifted (0,0.5*lbl);label.top("Earth's orbit", (EO + 10*sft)) shifted (0,0.5*lbl);label.top("Inner planet's", (IPO + 18*sft)) shifted (0,0.7*lbl);label.top("orbit", (IPO + 18*sft)) shifted (-3*lbl,-0.5*lbl);