Берем уравнение сферической тригонометрии:
cos(z) = sin(Lat) * sin(Dec) + cos(Lat) * cos(Dec) * cos(t)
z - зенитное расстояние
Lat - широта места наблюдения
Dec - склонение
t - часовой угол
Приплетаем уравнение:
S = Ra + t
S - местное звездное время
Ra - прямое восхождение
t - часовой угол
Меняем в первом уравнении t = s - RA и раскладываем косинус часового угла:
cos(z) = sin(Lat)*sin(Dec) + cos(Lat)*cos(Dec)*cos(S)*cos(Ra) + cos(Lat)*cos(Dec)*sin(S)*sin(Ra)
Вектор звезды в экваториальной системе координат
X(1) = cos(Ra)*cos(Dec)
X(2) = sin(Ra)*cos(Dec)
X(3) = sin(Dec))
Следовательно:
cos(z)=sin(Lat)*
X(3) + cos(Lat)*cos(S) *
X(1) + cos(Lat)*sin(S)*
X(2)Будем решать систему таких уравнений.
Измерений зенитного угла должно быть не меньше трех штук (3 неизвестных потому что), желательно достаточно разнесенных.
Пусть известно N значений времени и зенитного угла. Время переведем в местное звездное, то есть у нас имеется множество значений S(i) и z(i), i=1 до N.
Решаем матричное уравнение : A *
X = y
A(i,1) = cos(Lat)*cos(S(i))
A(i,2) = cos(Lat)*sin(S(i))
A(i,3) = sin(Lat)
y(i) = cos(z(i))
Вектор X даст искомые экваториальные координаты.
tg(Ra) = X(2)/X(1)
tg(Dec) = X(3) / sqrt(X(1)^2 + X(2)^2)
Тест: Задаем Ra = 23.5793 Dec = 55.9123
Дату и координаты места наблюдения не скажу

Сгенерируем якобы измеренные значения:
S(1) = 352.111796315506 z(1) = 23.425130328515
S(2) = 7.15286495434327 z(2) = 16.124207531734
S(3) = 22.1939335931805 z(3) = 12.2870295719788
Измерений три штуки с шагом в 1 час. Можно взять и 100 измерений, но так как ошибки измерений мы не внесли, то ничего минимизировать не надо и хватит 3 штук.
Решили систему уравнений каким-либо численным методом (если измерений всего три, то можно и аналитически вывести решение) и получили вектор x:
x(1) = 0.513666802206147
x(2) = 0.224194541636816
x(3) = 0.828180670996113
Разложили вектор на углы:
Ra = 23.5793000000001
Dec = 55.9123000000005