I hope this is not too off topic for the list.
While profiling
DateTime::Indic with Devel::NYTProf, I discovered that my calls to
DateTime::Event::Lunar ::new_moon_before and ::new_moon_after were major
bottlenecks. So I decided to replace them with code which is ultimately
derived from Meeus's Astronomical Algorithms chapter 32. It makes a huge
difference in speed but there is one problem.
Take a look at the output of the short example program below. In most
cases, the DateTime::Event::Lunar functions and my newmoon function give
almost the same results. But every now again newmoon() is off by an
entire lunar cycle. Why is this? How can I fix it?
#!/usr/bin/perl
use strict;
use warnings;
use DateTime;
use DateTime::Event::Lunar;
use DateTime::Event::Sunrise;
use Math::Trig qw/ deg2rad /;
use POSIX qw/ floor /;
my $date = DateTime->new(year => 2008, month => 1, day => 1);
my $sun = DateTime::Event::Sunrise->new(
latitude => '18.96',
longitude => '72.82',
altitude => 0,
);
# year 2008
for (1 .. 366) {
# time of sunrise on $date adjusted for Indian Standard Time
my $sunrise = $sun->sunrise_datetime($date)->set_time_zone('Asia/Kolkata');
# old method
my $pnm1 = DateTime::Event::Lunar->new_moon_before(
datetime => $sunrise,
)->jd;
my $nnm1 = DateTime::Event::Lunar->new_moon_after(
datetime => $sunrise,
)->jd;
# new method
# 0 = previous new moon before $sunrise
my $pnm2 = newmoon($sunrise->jd, 0);
# 1 = next new moon before $sunrise
my $nnm2 = newmoon($sunrise->jd, 1);
my $pnmdelta = $pnm2 - $pnm1;
my $nnmdelta = $nnm2 - $nnm1;
print "$sunrise $pnmdelta $nnmdelta\n";
$date->add(days => 1);
}
use constant J1900 => 2_415_020.0;
use constant synodic_month => 29.530_588_68;
sub newmoon {
my ($jdate, $arg) = @_;
my $k = floor((($jdate - J1900) / 365.25) * 12.3685) + $arg;
# time in Julian centuries since J1900
my $t = ($jdate - J1900) / 36525.0;
# square for frequent use
my $t2 = $t*$t;
# cube for frequent use
my $t3 = $t2*$t;
# mean time of phase
my $jdnv = (2_415_020.759_33) + synodic_month * $k + 0.000_117_8 * $t2
- 0.000_000_155 * $t3
+ 0.000_33 * sin(deg2rad(166.56 + 132.87 * $t - 0.009_173 * $t2));
# Sun's mean anomaly
my $m = deg2rad(359.224_2 + 29.105_356_08 * $k - 0.000_033_3 * $t2
- 0.000_003_47 * $t3);
# Moon's mean anomaly
my $mprime = deg2rad(306.025_3 + 385.816_918_06 * $k + 0.010_730_6 * $t2
+ 0.000_012_36 * $t3);
# Moon's argument of latitude
my $f = deg2rad(21.296_4 + 390.670_506_46 * $k - 0.001_652_8 * $t2
- 0.000_002_39 * $t3);
# Correction for new moon
my $djd = (0.1734 - 0.000_393 * $t) * sin($m) + 0.002_1 * sin(2 * $m);
$djd = $djd - 0.406_8 * sin($mprime) + 0.016_1 * sin(2 * $mprime);
$djd = $djd - 0.000_4 * sin(3 * $mprime) + 0.010_4 * sin(2 * $f) ;
$djd = $djd - 0.005_1 * sin($m + $mprime) - 0.007_4 * sin($m - $mprime);
$djd = $djd + 0.000_4 * sin(2 * $f + $m) - 0.000_4 * sin(2 * $f - $m);
$djd = $djd - 0.000_6 * sin(2 * $f + $mprime) + 0.001 *
sin(2 * $f - $mprime);
$djd = $djd + 0.000_5 * sin($m + 2 * $mprime);
$jdnv += $djd;
return $jdnv;
}
--
Jaldhar H. Vyas <[email protected]>
