What is wrong with this function?

[Jaldhar H. Vyas]

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 <jald...@braincells.com>