Here is how to create this example within Mathematica:

independentVariables = {ax -> 0, ay -> 0, az -> 1, 
   bx -> 0, by -> -1, bz -> 0};

r = Cross[Cross[{px, py, pz}, {qx, qy, qz}], {px, py, pz}];

dependentVariables = {al -> Sqrt[ax*ax + ay*ay + az*az], 
   px -> ax/al, py -> ay/al, pz -> az/al, 
   ax -> px, ay -> py, az -> pz,
   bl -> Sqrt[bx*bx + by*by + bz*bz], 
   qx -> bx/bl, qy -> by/bl, qz -> bz/bl, 
   bx -> qx, by -> qy, bz -> qz,
   rx -> r[[1]], ry -> r[[2]], rz -> r[[3]], 
   rl -> Sqrt[rx*rx + ry*ry + rz*rz], 
   sx -> rx/rl, sy -> ry/rl, sz -> rz/rl, 
   angle -> ArcCos[px*qx + py*qy + pz*qz]};

ntheta = 16; dtheta = Pi/ntheta; nphi = 16; dphi = 2Pi/nphi; ni = 16;

scene = Graphics3D[{{PointSize[0.02], Point[{ax, ay, az}], 
   Point[{bx, by, bz}], Line[{{ax, ay, az}, {px, py, pz}}], 
   Line[{{bx, by, bz}, {qx, qy, qz}}]}, 
   {AbsoluteThickness[0.5], GrayLevel[0.66], 
   Table[Line[{{Sin[theta]Cos[phi], Sin[theta]Sin[phi], 
   Cos[theta]}, {Sin[theta]Cos[phi + dphi], 
   Sin[theta]Sin[phi + dphi], Cos[theta]}}], 
   {theta, dtheta, Pi - dtheta, dtheta}, {phi, 0, 2Pi - dphi, dphi}], 
   Table[Line[{{Sin[theta]Cos[phi], Sin[theta]Sin[phi], 
   Cos[theta]}, {Sin[theta + dtheta]Cos[phi], 
   Sin[theta + dtheta]Sin[phi], Cos[theta + dtheta]}}], 
   {phi, 0, 2Pi - dphi, dphi}, {theta, dtheta, Pi - 2dtheta, dtheta}]}, 
   {Thickness[0.01], 
   Table[Line[{Cos[angle*i/ni]*{px, py, pz} + 
   Sin[angle*i/ni]*{sx, sy, sz}, 
   Cos[angle*(i + 1)/ni]*{px, py, pz} + 
   Sin[angle*(i + 1)/ni]*{sx, sy, sz}}], {i, 0, ni - 1}], 
   Thickness[0.005], GrayLevel[0.25], 
   Table[Line[{Cos[angle + (2 Pi - angle)*i/2/ni]*{px, py, pz} + 
   Sin[angle + (2 Pi - angle)*i/2/ni]*{sx, sy, sz}, 
   Cos[angle + (2 Pi - angle)*(i + 1)/2/ni]*{px, py, pz} + 
   Sin[angle + (2 Pi - angle)*(i + 1)/2/ni]*{sx, sy, sz}}], 
   {i, 0, 2 ni - 1}]}}, 
   Boxed -> False];

Show[N[scene] //. Join[independentVariables, dependentVariables]]