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;
	}
	return $sum;
}

/* Meeus, chapter 10 */
function deltaDTtoUT($year)
{
	$t = ($year - 2000) / 100;

	if ($year < 948) {
		$dT = 2177 + (497 * $t) + (44.1 * $t * $t);
	}
	
	/* There is a table on page 79 for the years 1620 to 1998, which I didn't
	 * bother to include here */
	
	if (($year >= 948 && $year < 1600) || $year >= 2000) {
		$dT = 102 + (102 * $t) + (25.3 * $t * $t);
	}

	/* to avoid a discontinuity at A.D. 2000, it is advised to add the
	 * correction +0.37 x (year - 2100) for the years 2000 to 2100 */
	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";

I believe we have among us, folks who can take this script to a level where we can run it normally in Calendar Scripts to calculate Equinox for any year on the fly.

Please it will be extremely helpful in dealing with Lunar and Lunisolar Calendars. Java Developers have adapted and now use it in their Scripts.
Please let someone take up the challenge of abstracting / encapsulating the script so we can simply stick a method or function in-script and compute any year's equinox or solstice date and time for calendation purposes.

It could be even part of the DateTime Class. As we have Easter Date function, let's have the equinox / solstice function or method.

Adolf

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.