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";

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.

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.