A general expression for the great circle of radius R, passing through points A= R ( ax, ay, az) and B= R (bx, by, bz) (where (ax, ay, az) is the unit vector a along OA and (bx, by, bz) is the unit vector b along OB) is given by
x/R = [Cosψ]. ax + [Sinψ].{(az2+ay2)bx - (azbz+ayby)ax}/sinδ
y/R= [Cosψ].ay + [Sinψ]{(az2+ax2)by - (azbz+axbx)ay}/sinδ
z/R= [Cosψ].az + [Sinψ]{(ay2+ax2)bz - (ayby+axbx)az}/sinδ
Cos δ = axbx + ayby + azbz (δ is the angle between OA and OB vectors)
These are derived by setting up a new, right handed, orthonormal set of basis vectors e1, e2, e3 oriented such that e1 is along OA, and e3 is perpendicular to the great circle. In this basis the equation for any point Q on the gt circle is simply:
Q/R = Cosψ e1 + Sinψ e2
as ψ increases from 0 to δ we sweep along the great circle from A to B. For the complete circle
0≤ψ≤2π radians. To get the equations for the x, y, z components of the circle we just relate the new basis e1 etc to the old basis i,j,k. Explicitly:
e1= a = axi + ayj + azk
e2= e3 X e1
e3= a X b /Sinδ =
(1/sinδ).
\begin{vmatrix}
i&j&k\\
a_x&a_y&a_z\\
b_x&b_y&b_z\end{vmatrix}
(nb remember the sinδ factor: since the unit vectorsa and b are not in general orthogonal a X b= sinδ ≠ 1 so to get a unit vector e3 we must divide by sinδ)
then x/R = i.Q/R etc