view release on metacpan or  search on metacpan

lib/Astro/Utils.pm  view on Meta::CPAN

    my $jd = _calc_jd($k, $year);
    return _jd_to_utc($jd);
}

sub _calc_utc_solstice {
    my ($k, $year) = @_;

    die "ERROR: Invalid type for solstice.\n" unless (defined $k && ($k =~ /^[1|3]$/));
    my $jd = _calc_jd($k, $year);
    return _jd_to_utc($jd);
}

sub _calc_tdt_equinox {
    my ($k, $year) = @_;

    die "ERROR: Invalid type for equinox.\n" unless (defined $k && ($k =~ /^[0|2]$/));
    my $jd = _calc_jd($k, $year);
    return _jd_to_tdt($jd);
}

sub _calc_tdt_solstice {
    my ($k, $year) = @_;

    die "ERROR: Invalid type for solstice.\n" unless (defined $k && ($k =~ /^[1|3]$/));
    my $jd = _calc_jd($k, $year);
    return _jd_to_tdt($jd);
}

sub _calc_jd {
    my ($k, $year) = @_;

    # Astronmical Algorithms, Chapter 27.
    my $jde0 = _jde0($k, $year);
    my $t    = ($jde0 - 2451545.0) / 36525;
    my $w    = 35999.373 * $t - 2.47;
    my $dl   = 1 + 0.0334 * _cos($w) + 0.0007 * _cos(2 * $w);
    my $s    = _periodic_term_24($t);

    return ($jde0 + ( (0.00001 * $s) / $dl ));
}

sub _jd_to_utc {
    my ($jd) = @_;

    my ($yr, $mon, $day, $hr, $min, $sec) = _process($jd);

    my $date = DateTime->new(year => $yr, month => $mon, day => $day, hour => $hr, minute => $min, second => $sec);

    # Astronmical Algorithms, Chapter 10, page 79 (Table 10.A).
    my $first = 1620;
    my $last  = 2002;
    my @tbl   = (121,112,103,95,88,82,77,72,68,63,60,56,53,51,48,46,44,42,40,38,             # 1620
                 35,33,31,29,26,24,22,20,18,16,14,12,11,10,9,8,7,7,7,7,                      # 1660
                 7,7,8,8,9,9,9,9,9,10,10,10,10,10,10,10,10,11,11,11,                         # 1700
                 11,11,12,12,12,12,13,13,13,14,14,14,14,15,15,15,15,15,16,16,                # 1740
                 16,16,16,16,16,16,15,15,14,13,                                              # 1780
                 13.1,12.5,12.2,12.0,12.0,12.0,12.0,12.0,12.0,11.9,11.6,11.0,10.2,9.2,8.2,   # 1800
                 7.1,6.2,5.6,5.4,5.3,5.4,5.6,5.9,6.2,6.5,6.8,7.1,7.3,7.5,7.6,                # 1830
                 7.7,7.3,6.2,5.2,2.7,1.4,-1.2,-2.8,-3.8,-4.8,-5.5,-5.3,-5.6,-5.7,-5.9,       # 1860
                 -6.0,-6.3,-6.5,-6.2,-4.7,-2.8,-0.1,2.6,5.3,7.7,10.4,13.3,16.0,18.2,20.2,    # 1890
                 21.1,22.4,23.5,23.8,24.3,24.0,23.9,23.9,23.7,24.0,24.3,25.3,26.2,27.3,28.2, # 1920
                 29.1,30.0,30.7,31.4,32.2,33.1,34.0,35.0,36.5,38.3,40.2,42.2,44.5,46.5,48.5, # 1950
                 50.5,52.5,53.8,54.9,55.8,56.9,58.3,60.0,61.6,63.0,63.8,64.3                 # 1980, 2002 last entry
    );

    my $delta_t = 0;
    my $t = ($yr - 2000) / 100;

    if ($yr >= $first && $yr <= $last) {
        if ($yr % 2) {
            $delta_t = ($tbl[($yr - $first - 1) / 2] + $tbl[($yr - $first + 1) / 2] ) / 2;
        }
        else {
            $delta_t = $tbl[($yr - $first) / 2];
        }
    } elsif ($yr < 948) {
        $delta_t = 2177 + 497 * $t + (44.1 * _pow($t, 2));
    } elsif ($yr >= 948) {
        $delta_t =  102 + 102 * $t + (25.3 * _pow($t, 2));
        if ($yr >= 2000 && $yr <= 2100) {
            $delta_t += 0.37 * ($yr - 2100);
        }
    }

    $date->subtract(seconds => $delta_t);

    return sprintf("%04d-%02d-%02d %02d:%02d:%02d",
                   $date->year, $date->month, $date->day, $date->hour, $date->minute, $date->second);
}

sub _jd_to_tdt {
    my ($jd) = @_;

    my ($yr, $mon, $day, $hr, $min, $sec) = _process($jd);
    return sprintf("%04d-%02d-%02d %02d:%02d:%02d", $yr, $mon, $day, $hr, $min, $sec);
}

sub _process {
    my ($jd) = @_;

    # Astronmical Algorithms Chapter 7, page 63.
    my ($a, $alpha);
    my $z = int($jd + 0.5);
    my $f = ($jd + 0.5) - $z;

    if ($z < 2299161) {
        $a = $z;
    }
    else {
        $alpha = int(($z - 1867216.25) / 36524.25);
        $a = $z + 1 + $alpha - int($alpha / 4);
    }

    my $b   = $a + 1524;
    my $c   = int(($b - 122.1) / 365.25);
    my $d   = int(365.25 * $c);
    my $e   = int(($b -$d) / 30.6001);
    my $dt  = $b - $d - int(30.6001 * $e) + $f;
    my $mon = $e - (($e < 13.5)?(1):(13));
    my $yr  = $c - (($mon > 2.5)?(4716):(4715));
    my $day = int($dt);