PHP script to calculate the times of solstices and equinoxes
<?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); | |
} | |
/* This is optional | |
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