Как рассчитать продолжительность светового дня в любой день года

Недавно мне в голову пришла одна идея: всем известно что в течение года время светового дня меняется - зимой день короче, летом длиннее. При этом есть еще и точки минимальной и максимальной продолжительности дня - дни солнцестояния зимнего и летнего, а так же есть дни равноденствия осеннего и весеннего - время когда день равен ночи. Можно предположить что в дни равноденствия, например 21 марта - день весеннего равноденствия, когда день равен ночи, продолжительность дня составит 12 часов. На самом деле это почти так, но не так, на это время влияет куча факторов, и если посмотреть на сайте https://www.timeanddate.com/sun/kazakstan/astana то можно увидеть что реально продолжительность дня отличается от 12 часов. Ну да ладно. Просто мне пришла идея - можно ли рассчитать время продолжительности дня в любое время года. На самом деле можно, именно так и работает приведенный выше сайт, но мне было интересно самому разобраться и вот что я нарыл по этому поводу:

Для меня было совершенно очевидно, что изменение времени происходят циклически, из года в года, почти как синусоида. И потому решил сделать такое уравнение на основе синусоиды. Дальше пришлось вспоминать курс математики, открыл википедию, и нашел там формулу синусоиды

  Y=a+b*sin(c*x+d)

где a - смещение по оси OY, b - амплитуда (по оси OY),c - растяжение по оси OX, d - cмещение по оси OX

Осталось только рассчитать все эти коэффициенты. На сайте взял данные по максимальной продолжительности дня в Астане, и минимальной, это позволяет определить коэффициенты смещения и амплитуды. Потом нужно было весь год уместить в один период синусоиды 2*пи. Ну и в конце-концов, максимум нужно было сместить по оси ОХ чтобы летнее солнцестояние пришлось на 21 июня.

Все данные накидал в Excel и вывел такую формулу

Y=12.231+4.338sin(0.01721x-1.3728)

Но после того как заполнил всю таблицу на полный год, увидел большие расхождения относительно данных, взятых на сайте, и поэтому решил добавить данные с сайта, и построить дополнительный график. Искажения синусоиды происходят в следствие наклона земной оси, и эллиптической орбиты Земли . 

Что касается расчетов, то амплитуда рассчитывается как максимальное значение минус минимальное и деленное на 2, смещение по вертикальной оси, как максимум минус амплитуда, либо минимум + амплитуда. Коэффициент растяжения как я и говорил - нужно уместить в период 2*пи весь год - 365 дней, решается пропорцией. Смещение по оси Х определяется так: весь год 365 дней - это полный период, но точка максимума на синусоиде, находится в 1/4 периода, значит нужно разделить 365 на 4 получилось 91,25, но максимум в нашей системе - это день 171, и потому нужно отнять от 171-91,25=79,75, ну и поскольку все это перекладывается опять же на координатную сетку с масштабированную под 2пи, то через пропорцию находим коэффициент смещения.