Допустим, у нас есть эклиптическая система координат, в ней и работаем.
Одна из точек большого круга (БК) лежит на эклиптике. Координаты этой точки: долгота L1, широта b1(b1=0, потому что на эклиптике).
Возьмем в качестве второй точки БК точку кульминации этого БК в эклиптической системе координат. Ее координаты: L2 = L1 + 90, b2 = u (u - это угол, под которым БК пересекает эклиптику, например 30 градусов).
Вектор v1 на сфере: (cos(L1), sin(L1), 0)
Вектор v2 на сфере: (cos(L1+90)*cos(u), sin(L1+90)*cos(u), sin(u))
Нормаль плоскости (она же является вектором полюса БК): n = v1 x v2 (векторное произведение векторов, принадлежащих БК).
n = (A, B, C) = (sin(L1)*sin(u), -cos(L1)*sin(u), cos(u))
Уравнение плоскости, проходящей через центр сферы: A*x+B*y+C*z=0
В сферической эклиптической системе координат:
sin(L1)*sin(u)*cos(L)*cos(b) - cos(L1)*sin(u)*sin(L)*cos(b) + cos(u)*sin(b)=0
где L, b - эклиптические координаты произвольной точки, лежащей на искомом БК.
Делим это уравнение на cos(b) и cos(u) и получаем:
tg(b) = tg(u) * sin(L - L1)
Изменяя L от 0 до 360 градусов, можно вычислить все соответствующие широты как b = arctg(tg(u) * sin(L - L1))