Last active
December 16, 2023 10:57
-
-
Save derickr/f32dd7a05d5c0a099db4e449111f5ccd to your computer and use it in GitHub Desktop.
PHP script to calculate the times of solstices and equinoxes
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
<?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"; |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I've made the updates by @pstray and @davidefa from their comments above.