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 elongation
pair GWE;  % greatest western elongation
pair EQ;  % eastern quadrature
pair WQ;  % western quadrature
pair Sun;  % position of the Sun, defines the origin
  Sun = (0,0);
pair Earth;  % position of the Earth
pair planet;
numeric GEEangle, GWEangle, EQangle, WQangle, planetangle;
GEEangle = 300;
GWEangle = 60;
EQangle = 311.8103;
WQangle = 48.1897;  % arccos(2/3)

pair Conj;  % conjunction
pair Opp;  % opposition
pair SupConj;  % superior conjunction
pair InfConj;  % inferior conjunction

% draw the circles and define points of interest
for 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 interest
drawcircleS(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 lines
draw 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 Earth
drawangle(Earth,planet,Sun,1.0cm,1pt,black);
drawangle(Sun,Earth,planet,0.7cm,1pt,black);

% labelling
label(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);

% arrows
pair 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);