...смысл в том чтобы используя формулы получить значения склонения, которые в свою очередь используются для ...
поэтому: как найти эти самые исходные формулы?
Смотря какая точность требуется. Если пренебречь эллиптичностью земной орбиты (а она невелика), то формула зависимости склонения от порядкового дня в году будет элементарной: delta = eps*sin((D-D0)/365.25), где eps - угол наклона эклиптики к экватору (~ 23°26'), D - порядковый номер текущего (т.е. для которого ищется склонение Солнца) дня в году, D0 - порядковый номер дня весеннего равноденствия (примерно 81 - это можно уточнить). С порядковыми днями удобно работать через юлианские, но если особой точности не требуется, то есть разные пути радикального упрощения таких вычислений.
А для большей точности уже придётся учитывать неравномерность орбитального движения Земли - в первом приближении движущееся во эклиптике Солнце заставить двигаться с переменной скорость (на среднюю скорость наложить соответствующую синусоиду), а затем эклиптические координаты перевести в экваториальные. И тогда, уже после перевода в горизонтальные, не забыть ещё и рефракцию учесть...