Hi, dt-ers.
I've been lurking on this list for some time now, and so far I
understand that while a few people have attempted, nobody has come up
with a lunar, solar, or lunisolar calenders (at least I don't remember
seeing it on this list). And so I was just fooling around with some code
to see if I can at least start creating some of the basic building
blocks for it.
# And if somebody already has implemented this, please let me know so
# I don't waste time ;)
My first stab at it was implementing the function to get the solar
longitude at time t, from [1] p.182. I've done it, and I'm rather proud
of the fact that I actually took the time to understand all this with no
prior knowledge of astronomical calculations, but the thing is, the
results are a bit off:
[EMAIL PROTECTED] lib]$ perl test.pl
R.D. | Solar Longitude | Solar Longitude (from table) |
-214193 | 119.964425 | 119.474975
-61387 | 254.767876 | 254.252390
25469 | 181.939866 | 181.435260
49217 | 189.168174 | 188.662093
171307 | 289.601843 | 289.089403
210155 | 59.601590 | 59.119357
253427 | 228.828880 | 228.316498
[EMAIL PROTECTED] lib]$
The numbers are consistently off by 0.5 degrees or so. eek.
This must be something so stupid in my code, but I've been looking at
the code for way too long. If anybody notices anything wrong, I would
appreciate it. A lot.
I have two files attatched. Astro::Sun is the module that implements the
logic, and test.pl generates the partial table from [1] p.400 (errata 285).
To whomever that might take a look at this code, I thank you very much
in advance.
--d
P.S. please don't bother too much with the module naming. That's just
temporary...
[1] Edward M. Reingold, Nachum Dershowitz
"Calendrical Calculations (Millenium Edition)", 2nd ed.
Cambridge University Press, Cambridge, UK 2002
# REFERENCES:
# [1] Edward M. Reingold, Nachum Dershowitz
# "Calendrical Calculations (Millenium Edition)", 2nd ed.
# Cambridge University Press, Cambridge, UK 2002
package Astro::Sun;
use strict;
use DateTime;
use Math::Trig qw(deg2rad);
# I got the following from (DateTime->new(...)->utc_rd_values)[0]
use constant RD_1900_JAN_1 => 693596;
use constant RD_1810_JAN_1 => 660724;
use constant RD_MOMENT_J2000 => 730120.5;
# [1] Table 12.1 p.183 (zero-padded to align veritcally)
use constant SOLAR_LONGITUDE_ARGS => [
# X Y Z
# left side of table 12.1
[ 403406, 270.54861, 0.9287892 ],
[ 119433, 63.91854, 35999.4089666 ],
[ 3891, 317.84300, 71998.2026100 ],
[ 1721, 240.05200, 36000.3572600 ],
[ 350, 247.23000, 32964.4678000 ],
[ 314, 297.82000, 445267.1117000 ],
[ 242, 166.79000, 3.1008000 ],
[ 158, 3.50000, -19.9739000 ],
[ 129, 182.95000, 9038.0293000 ],
[ 99, 29.80000, 33718.1480000 ],
[ 86, 249.20000, -2280.7730000 ],
[ 72, 257.80000, 31556.4930000 ],
[ 64, 69.90000, 9037.7500000 ],
[ 38, 197.10000, -4444.1760000 ],
[ 32, 65.30000, 67555.3160000 ],
[ 28, 341.50000, -4561.5400000 ],
[ 27, 98.50000, 1221.6550000 ],
[ 24, 110.00000, 31437.3690000 ],
[ 21, 342.60000, -31931.7570000 ],
[ 18, 256.10000, 1221.9990000 ],
[ 14, 242.90000, -4442.0390000 ],
[ 13, 151.80000, 119.0660000 ],
[ 12, 53.30000, -4.5780000 ],
[ 10, 205.70000, -39.1270000 ],
[ 10, 146.10000, 90073.7780000 ],
# right side of table 12.1
[ 195207, 340.19128, 35999.1376958 ],
[ 112392, 331.26220, 35998.7287385 ],
[ 2819, 86.63100, 71998.4403000 ],
[ 660, 310.26000, 71997.4812000 ],
[ 334, 260.87000, -19.4410000 ],
[ 268, 343.14000, 45036.8840000 ],
[ 234, 81.53000, 22518.4434000 ],
[ 132, 132.75000, 65928.9345000 ],
[ 114, 162.03000, 3034.7684000 ],
[ 93, 266.40000, 3034.4480000 ],
[ 78, 157.60000, 29929.9920000 ],
[ 68, 185.10000, 149.5880000 ],
[ 46, 8.00000, 107997.4050000 ],
[ 37, 250.40000, 151.7710000 ],
[ 29, 162.70000, 31556.0800000 ],
[ 27, 291.60000, 107996.7060000 ],
[ 25, 146.70000, 62894.1670000 ],
[ 21, 5.20000, 14578.2980000 ],
[ 20, 230.90000, 34777.2430000 ],
[ 17, 45.30000, 62894.5110000 ],
[ 13, 115.20000, 107997.9090000 ],
[ 13, 285.30000, 16859.0710000 ],
[ 10, 126.60000, 26895.2920000 ],
[ 10, 85.90000, 12297.5360000 ]
];
# [1] p.19
sub _dt_to_moment
{
my $dt = shift;
my($rd, $seconds) = $dt->utc_rd_values;
($rd + $seconds / (24 * 3600));
}
# [1] p172
sub _dynamical_moment_from_dt
{
my $dt = shift;
_dt_to_moment($dt) + ephemeris_correction($dt);
}
# [1] p171 + errata 158
sub ephemeris_correction
{
my $dt = shift;
# we need a gregorian calendar, so make sure $dt is just 'DateTime'
$dt = DateTime->from_object(object => $dt);
my $year = $dt->year;
# XXX - possible optimization
my $c = RD_1900_JAN_1 -
(DateTime->new(year => $year, month => 7, day => 1)->utc_rd_values)[0];
my $x = 0.5 + RD_1810_JAN_1 -
(DateTime->new(year => $year, month => 1, day => 1)->utc_rd_values)[0];
my $correction;
if (1988 <= $year && $year <= 2019) {
$correction = ($year - 1933) / (24 * 3600);
} elsif (1900 <= $year && $year <= 1987) {
$correction = -0.00002 + 0.000297 * $c + 0.025184 * ($c ** 2) -
0.181133 * ($c ** 3) + 0.553040 * ($c ** 4) -
0.861938 * ($c ** 5) + 0.677066 * ($c ** 6) -
0.212591 * ($c ** 7);
} elsif (1800 <= $year && $year <= 1899) {
$correction = -0.00009 + 0.003844 * $c + 0.083563 * ($c ** 2) -
0.865736 * ($c ** 3) + 4.867575 * ($c ** 4) -
15.845535 * ($c ** 5) + 31.332267 * ($c ** 6) -
38.291999 * ($c ** 7) + 28.316289 * ($c ** 8) +
11.636204 * ($c ** 9) + 2.043794 * ($c ** 10);
} elsif (1700 <= $year && $year <= 1799) {
$correction = (8.118780842 - 0.005092142 * ($year - 1700) +
0.003336121 * (($year - 1700) ** 2) -
0.0000266484 * (($year = 1700) ** 3)) / ( 24 * 3600 );
} elsif (1620 <= $year && $year <= 1699) {
$correction = (196.58333 - 4.0675 * ($year - 1600) +
0.0219167 * (($year - 1600) ** 2)) / ( 24 * 3600 );
} else {
$correction = (($x ** 2 / 41048480 ) - 15) / ( 24 * 3600 );
}
return $correction;
}
# [1] p.172
sub julian_centuries
{
my $dt = shift;
return (_dynamical_moment_from_dt($dt) - RD_MOMENT_J2000) / 36525;
}
# [1] p.182
sub nutation
{
my $dt = shift;
my $c = julian_centuries($dt);
my $A = 124.90 - 1934.134 * $c + 0.002063 * ($c ** 2);
my $B = 201.11 + 72001.5377 * $c + 0.00057 * ($c ** 2);
return (-0.004778 * sin(deg2rad($A))) + (-0.0003667 * sin(deg2rad($B)));
}
# [1] p.183
sub aberration
{
my $dt = shift;
my $c = julian_centuries($dt);
return 0.0000974 * cos(deg2rad(177.63 + 35999.01848 * $c)) - 0.0005575;
}
# [1] p.182
sub solar_longitude
{
my $dt = shift;
my $c = julian_centuries($dt);
my $big_ugly_number = 0;
foreach my $numbers (@{ SOLAR_LONGITUDE_ARGS() }) {
$big_ugly_number += $numbers->[0] *
sin(deg2rad($numbers->[1] + $numbers->[2] * $c))
}
my $longitude = 282.7771834 +
36000.76953744 * $c +
0.000005729577951308232 * $big_ugly_number;
# The "mod" operator here truncates our results UP, which is bad.
# we first get the floor, take the mod, and then add the fractional
# parts to it
my $pre_processed_number = ($longitude + aberration($dt) + nutation($dt));
my $floor = int($pre_processed_number);
return $floor % 360 + $pre_processed_number - $floor;
}
1;
use strict;
use DateTime;
use Astro::Sun;
my @rd2solar_longitude = (
[ -214193, 119.474975 ],
[ -61387, 254.252390 ],
[ 25469, 181.435260 ],
[ 49217, 188.662093 ],
[ 171307, 289.089403 ],
[ 210155, 59.119357 ],
[ 253427, 228.316498 ],
);
print " R.D. | Solar Longitude | Solar Longitude (from table) |\n";
# 1234567890 123456789012345
foreach my $data (@rd2solar_longitude) {
my $dt = DateTime->new(year => 1, month => 1, day => 1);
$dt += DateTime::Duration->new(days => $data->[0]);
printf( "%9d | %0.6f | %0.6f\n",
$data->[0],
Astro::Sun::solar_longitude($dt),
$data->[1]
);
}
