PHP script to calculate the times of solstices and equinoxes

equinox.php

<?php
/* Algorithms taken from Meeus Astronomical Algorithms, 2nd edition */

/* Run the php script on the command line:
 *   php equinox.php <year> <what>
 * with <what> being MAR, JUN, SEP, or DEC
 *
 * For the the Summer Solstice of 2020:
 *   php equinox.php 2020 JUN
 */

/* Constants */
define('MAR', 0);
define('JUN', 1);
define('SEP', 2);
define('DEC', 3);

/* Read arguments from the command line */
$year = $argv[1];
$which = constant($argv[2]);

/* Tables from Meeus, page 178 */
/* Table 27.A, for the years -1000 to 1000 */
$yearTable0 = [
	MAR => [ 1721139.29189, 365242.13740, 0.06134, 0.00111, 0.00071 ],
	JUN => [ 1721233.25401, 365241.72562, 0.05323, 0.00907, 0.00025 ],
	SEP => [ 1721325.70455, 365242.49558, 0.11677, 0.00297, 0.00074 ],
	DEC => [ 1721414.39987, 365242.88257, 0.00769, 0.00933, 0.00006 ],
];

/* Table 27.B, for the years 1000 to 3000 */
$yearTable2000 = [
	MAR => [ 2451623.80984, 365242.37404, 0.05169, 0.00411, 0.00057 ],
	JUN => [ 2451716.56767, 365241.62603, 0.00325, 0.00888, 0.00030 ],
	SEP => [ 2451810.21715, 365242.01767, 0.11575, 0.00337, 0.00078 ],
	DEC => [ 2451900.05952, 365242.74049, 0.06223, 0.00823, 0.00032 ],
];

/* Meeus, page 177 */
function calculateJDE0($year, $which)
{
	global $yearTable0, $yearTable2000;

	$table = $year < 1000 ? $yearTable0 : $yearTable2000;
	$Y     = $year < 1000 ? ($year / 1000) : (($year - 2000) / 1000);
	$terms = $table[$which];

	$JDE0 = $terms[0] +
		($terms[1] * $Y) +
		($terms[2] * $Y * $Y) +
		($terms[3] * $Y * $Y * $Y) +
		($terms[4] * $Y * $Y * $Y * $Y);

	return $JDE0;
}

/* Meeus, Table 27.C, page 179 */
function calculateS($T)
{
	$table = [
		[ 485, 324.96,   1934.136 ],
		[ 203, 337.23,  32964.467 ],
		[ 199, 342.08,     20.186 ],
		[ 182,  27.85, 445267.112 ],
		[ 156,  73.14,  45036.886 ],
		[ 136, 171.52,  22518.443 ],
		[  77, 222.54,  65928.934 ],
		[  74, 296.72,   3034.906 ],
		[  70, 243.58,   9037.513 ],
		[  58, 119.81,  33718.147 ],
		[  52, 297.17,    150.678 ],
		[  50,  21.02,   2281.226 ],
		[  45, 247.54,  29929.562 ],
		[  44, 325.15,  31555.956 ],
		[  29,  60.93,   4443.417 ],
		[  18, 155.12,  67555.328 ],
		[  17, 288.79,   4562.452 ],
		[  16, 198.04,  62894.029 ],
		[  14, 199.76,  31436.921 ],
		[  12,  95.39,  14577.848 ],
		[  12, 287.11,  31931.756 ],
		[  12, 320.81,  34777.259 ],
		[   9, 227.73,   1222.114 ],
		[   8,  15.45,  16859.074 ],
	];
	
	$sum = 0;
	foreach( $table as $term ) {
		$c = $term[0] * cos(deg2rad( $term[1] + ($term[2] * $T) ));

		$sum += $c;
	}

	var_dump($term);
	return $sum;
}

function deltaDTtoUTFromTable($year)
{
	/* Meeus, page 79 */
	$deltaTtab = [
		 121,  112,  103,   95,   88,   82,   77,   72,   68,   63,   60,   56,   53,   51,   48,   46,   44,   42,   40,   38,
		  35,   33,   31,   29,   26,   24,   22,   20,   18,   16,   14,   12,   11,   10,    9,    8,    7,    7,    7,    7,
		   7,    7,    8,    8,    9,    9,    9,    9,    9,   10,   10,   10,   10,   10,   10,   10,   10,   11,   11,   11, 
		  11,   11,   12,   12,   12,   12,   13,   13,   13,   14,   14,   14,   14,   15,   15,   15,   15,   15,   16,   16,
		  16,   16,   16,   16,   16,   16,   15,   15,   14,   13, 13.1, 12.5, 12.2,   12,   12,   12,   12,   12,   12, 11.9,
		11.6,   11, 10.2,  9.2,  8.2,  7.1,  6.2,  5.6,  5.4,  5.3,  5.4,  5.6,  5.9,  6.2,  6.5,  6.8,  7.1,  7.3,  7.5,  7.6,
		 7.7,  7.3,  6.2,  5.2,  2.7,  1.4, -1.2, -2.8, -3.8, -4.8, -5.5, -5.3, -5.6, -5.7, -5.9,   -6, -6.3, -6.5, -6.2, -4.7,
		-2.8, -0.1,  2.6,  5.3,  7.7, 10.4, 13.3,   16, 18.2, 20.2, 21.1, 22.4, 23.5, 23.8, 24.3,   24, 23.9, 23.9, 23.7, 24.0,
		24.3, 25.3, 26.2, 27.3, 28.2, 29.1,   30, 30.7, 31.4, 32.2, 33.1,   34,   35, 36.5, 38.3, 40.2, 42.2, 44.5, 46.5, 48.5,
		50.5, 52.2, 53.8, 54.9, 55.8, 56.9, 58.3,   60, 61.6,   63,   65, 66.6,
	];

	$i = floor(($year - 1620) / 2);
	$f = (($year - 1620) / 2) - $i;  /* Fractional part of year */
	$dT = $deltaTtab[$i] + (($deltaTtab[$i + 1] - $deltaTtab[$i]) * $f);

	return $dT;
}

/* Meeus, chapter 10 */
function deltaDTtoUT($year)
{
	if (($year >= 1620) && ($year <= 2000)) {
		return deltaDTtoUTFromTable($year);
	}

	$t = ($year - 2000) / 100;

	if ($year < 948) {
		return 2177 + (497 * $t) + (44.1 * $t * $t);
	}
	
	if (($year >= 948 && $year < 1620) || $year >= 2000) {
		$dT = 102 + (102 * $t) + (25.3 * $t * $t);
	}

	if ($year >= 2000 && $year < 2100) {
		$dT += 0.37 * ($year - 2100);
	}

	return $dT;
}

function JDEtoTimestamp($JDE)
{
	$tmp = $JDE;
	$tmp -= 2440587.5;
	$tmp *= 86400;

	return $tmp;
}

/* Meeus, page 177 */
$JDE0 = calculateJDE0($year, $which);
$T   = ($JDE0 - 2451545.0) / 36525;
$W = (35999.373 * $T) - 2.47;
$L = 1 + 0.0334 * cos(deg2rad($W)) + 0.0007 * cos(2 * deg2rad($W));
$S = calculateS($T);

/* Meeus, page 178 */
$JDE = $JDE0 + ((0.00001 * $S) / $L);

/* Convert TD to PHP Date */
$date = JDEtoTimestamp($JDE);
$TD = new \DateTimeImmutable('@' . round($date));

/* Correct DT to UT, and convert DT to PHP Date */
$date += deltaDTtoUT($year);
$UT = new \DateTimeImmutable('@' . round($date));

/* Output */
echo "Dynamical Time: ", $TD->format(DateTimeImmutable::ISO8601), "\n";
echo "UT:             ", $UT->format(DateTimeImmutable::ISO8601), "\n";

1. line 145 should be ( on bottom of page 78 UT = TD - deltaT ):
   $date -= deltaDTtoUT($year);
2. enable additional correction factor ( needed for years 2000-2100 ) in lines 113-115

Bless you so much sir.
The errata remedy is highly appreciated.

Notice that this Script with so much potential can only be run on a Command Line! Please, please, sir, roll it all into a function or Class, so that weak folk like me can use it.

Please Sir,

Adolf

Seems you have forgotten the - signs in some of the numbers from tables 27.A and .B... the tables should be:

	MAR => [ 1721139.29189, +365242.13740, +0.06134, +0.00111, -0.00071 ],
	JUN => [ 1721233.25401, +365241.72562, -0.05323, +0.00907, +0.00025 ],
	SEP => [ 1721325.70455, +365242.49558, -0.11677, -0.00297, +0.00074 ],
	DEC => [ 1721414.39987, +365242.88257, -0.00769, -0.00933, -0.00006 ],

and

	MAR => [ 2451623.80984, +365242.37404, +0.05169, -0.00411, -0.00057 ],
	JUN => [ 2451716.56767, +365241.62603, +0.00325, +0.00888, -0.00030 ],
	SEP => [ 2451810.21715, +365242.01767, -0.11575, +0.00337, +0.00078 ],
	DEC => [ 2451900.05952, +365242.74049, -0.06223, -0.00823, +0.00032 ],

I've made the updates by @pstray and @davidefa from their comments above.

Sorry to dig this up, but this yelds an error if the year is between 1600 and 2000: no $dT in the deltaDTtoUT.
When you say you didn’t bother include that one table page 79, maybe you should have for those years.

You can do something like this:

function deltaDTtoUT($year) {

	/* Meeus, page 79 */
	$deltaTtab = [
		121,    112,   103,    95,    88,    82,    77,    72,    68,   63,    60,  56,  53,  51,  48, 46,
		44,      42,    40,    38,    35,    33,    31,    29,    26,   24,    22,  20,  18,  16,  14, 12,
		11,      10,     9,     8,     7,     7,     7,     7,     7,    7,     8,   8,   9,   9,   9,  9, 9, 10, 10, 10,
		10,      10,    10,    10,    10,    11,    11,    11,    11,   11,    12,  12,  12,  12,  13, 13,
		13,      14,    14,    14,    14,    15,    15,    15,    15,   15,    16,  16,  16,  16,  16, 16,
		16,      16,    15,    15,    14,    13,  13.1,  12.5,  12.2,   12,    12,  12,  12,  12,  12,
		11.9,  11.6,    11,  10.2,   9.2,   8.2,   7.1,   6.2,   5.6,  5.4,   5.3, 5.4, 5.6,
		5.9,    6.2,   6.5,   6.8,   7.1,   7.3,   7.5,   7.6,   7.7,  7.3,   6.2, 5.2, 2.7,
		1.4,   -1.2,  -2.8,  -3.8,  -4.8,  -5.5,  -5.3,  -5.6,  -5.7, -5.9,    -6,
		-6.3,  -6.5,  -6.2,  -4.7,  -2.8,  -0.1,   2.6,   5.3,   7.7, 10.4,  13.3,  16,
		18.2,  20.2,  21.1,  22.4,  23.5,  23.8,  24.3,    24,  23.9, 23.9,  23.7,
		24,    24.3,  25.3,  26.2,  27.3,  28.2,  29.1,    30,  30.7, 31.4,  32.2,
		33.1,    34,    35,  36.5,  38.3,  40.2,  42.2,  44.5,  46.5, 48.5,  50.5,
		52.2,  53.8,  54.9,  55.8,  56.9,  58.3,    60,  61.6,    63,   65,  66.6
	];

	if (($year >= 1620) && ($year <= 2000)) {
		$i = floor(($year - 1620) / 2);
		$f = (($year - 1620) / 2) - $i;  /* Fractional part of year */
		$dT = $deltaTtab[$i] + (($deltaTtab[$i + 1] - $deltaTtab[$i]) * $f);
	} else {
		$t = ($year - 2000) / 100;
		if ($year < 948) {
			$dT = 2177 + (497 * $t) + (44.1 * $t * $t);
		} else {
			$dT = 102 + (102 * $t) + (25.3 * $t * $t);
			if (($year > 2000) && ($year < 2100)) {
				$dT += 0.37 * ($year - 2100);
			}
		}
	}

	return $dT;
}

This way there won’t be any undefined variable dT error, and it will cover all the years.

Anyway, thanks for this script! It will be useful for me :)

I think I've found a better function, but it's in JavaScript: https://github.com/commenthol/astronomia/blob/master/src/deltat.js

In any case, I've updated the function with something similar as @timovn suggested.