Time.pm

package AstroScript::Time;

use warnings;
use strict;

our $VERSION = '1.00';

use Exporter qw/import/;

our $SEC_PER_DAY = 86400; # Seconds per day
our $SEC_PER_CEN = 3155760000;
our $J2000 = 2451545; # Standard Julian Date for 1.1.2000 12:00
our $J1900 = 2415020; # Standard Julian Date for  31.12.1899 12:00 (astronomical epoch 1900.0)
our $SOLAR_TO_SIDEREAL = 1.002737909350795; # Difference in between Sidereal and Solar hour (the former is shorter)
our $GREGORIAN_START = 15821004; # Start of Gregorian calendar (YYYYMMDD)
our $JD_UNIX_EPOCH = _gmtime2jd(gmtime(0));  # Standard Julian date for the beginning of Unix epoch, Jan 1 1970 on most Unix systems

our %EXPORT_TAGS = (
    all => [
        qw/jd_cent after_gregorian cal2jd jd2cal jd0 unix2jd jd2mjd mjd2jd
          jd2unix jdnow t1900 jd2dt jd2te jd2lst t2dt
          $SEC_PER_DAY $SEC_PER_CEN $J2000 $J1900 $GREGORIAN_START $JD_UNIX_EPOCH/
    ],
);

our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );

use POSIX;
use AstroScript::MathUtils qw/polynome ddd frac to_range/;

my @DT_TERMS = ( # terms for Delta-T
    [ 10.4, -80.8, 413.9,  -572.3 ],           # 1825-
    [ 6.6,  46.3,  -358.4, 18.8 ],             # 1850-
    [ -3.9, -10.8, -166.2, 867.4 ],            # 1875-
    [ -2.6, 114.1, 327.5,  -1467.4 ],          # 1900-
    [ 24.2, -6.3,  -8.2,   483.4 ],            # 1925-
    [ 29.3, 32.5,  -3.8,   550.7 ],            # 1950-
    [ 45.3, 130.5, -570.5, 1516.7 ],           # 1975-
    [ 0.25, 45.3,  130.5,  -570.5, 1516.7 ]    # 2000-2005
);

sub after_gregorian {
    my $y = shift;
    my $m = shift;
    my $d = shift;
    my %arg = (gregorian_start => $GREGORIAN_START, @_);
    return 0 unless defined $arg{gregorian_start};
    polynome( 100, $d, $m, $y ) >= $arg{gregorian_start};
}

sub cal2jd {
    my $ye = shift;
    my $mo = shift;
    my $da = shift;
    my %arg = (gregorian_start => $GREGORIAN_START, @_);

    my $j = $da + 1720996.5;
    my ( $m, $y ) = ( $mo > 2 ) ? ( $mo, $ye ) : ( $mo + 12, $ye - 1 );
    if ( after_gregorian( $ye, $mo, $da, %arg ) ) {
        $j += int( $y / 400 ) - int( $y / 100 ) + int( $y / 4 );
    }
    else {
        $j += int( ( $y + 4716 ) / 4 ) - 1181;
    }
    $j + 365 * $y + floor( 30.6001 * ( $m + 1 ) )
}

sub jd2cal {
    my $jd = shift;
    my %arg = (gregorian => 1, @_);

    my ( $f, $i ) = modf($jd - $J1900 + 0.5);
    if ($arg{gregorian} && $i > -115860) {
        my $a = floor($i / 36524.25 + 9.9835726e-1) + 14;
        $i += 1 + $a - floor($a / 4);
    }

    my $b = floor($i / 365.25 + 8.02601e-1);
    my $c = $i - floor(365.25 * $b + 7.50001e-1) + 416;
    my $g = floor($c / 30.6001);
    my $da = $c - floor(30.6001 * $g) + $f;
    my $mo = $g - ($g > 13.5 ? 13 : 1);
    my $ye = $b + ($mo < 2.5 ? 1900 : 1899);
    $ye, $mo, $da
}

sub jd0 {
    my $j = shift;
    floor( $j - 0.5 ) + 0.5;
}

sub unix2jd {
    $JD_UNIX_EPOCH + $_[0] / $SEC_PER_DAY
}

sub jd2unix {
    int( ($_[0] - $JD_UNIX_EPOCH) * $SEC_PER_DAY )
}

sub _gmtime2jd {
    cal2jd(
        $_[5] + 1900,
        $_[4] + 1,
        $_[3] + ddd(@_[2, 1, 0]) / 24
    )
}


sub jdnow {
    _gmtime2jd(gmtime())
}

sub jd2mjd {
    $_[0] - $J2000;
}


sub mjd2jd {
    $_[0] + $J2000;
}

# converts Julian date to period in centuries since epoch
# Arguments:
# julian date
# julian date corresponding to the epoch start
sub _t {
    my ( $jd, $epoch ) = @_;
    ( $jd - $epoch ) / 36525;
}


sub jd_cent {
    _t( $_[0], $J2000 );
}


sub t1900 {
    _t( $_[0], $J1900 );
}


sub t2dt {
    my $t2000 = shift;
    my ( $dt, $valid ) = ( 0, 0 );
    unless ( ( $t2000 < -1.75 ) or ( $t2000 > 0.05 ) ) {
        $valid = 1;
        my $i = int( floor( $t2000 / 0.25 ) );
        my $t = $t2000 - $i * 0.25;
        my @terms = @{ $DT_TERMS[ $i + 7 ] };
        $dt = polynome( $t, @terms );
    }
    $dt, $valid
}

sub jd2dt {
    my $j = shift;
    t2dt( jd_cent($j) );
}


sub jd2te {
    my $j = shift;
    my ( $dt, $valid ) = jd2dt($j);
    my $t = jd_cent($j);
    if ($valid) {
        $t += $dt / $SEC_PER_CEN;
    }
    $t, $dt;
}


sub jd2gst {
    my $jh = shift;
    my $j0 = jd0($jh);
    my $s0 = polynome( t1900($j0), 0.276919398, 100.0021359, 0.000001075 );
    24 * ( frac($s0) + abs( $jh - $j0 ) * $SOLAR_TO_SIDEREAL );
}

sub jd2lst {
    my ( $jd, $lon ) = @_;
    $lon //= 0;
    to_range( jd2gst($jd) - $lon / 15, 24 );
}


1;

__END__

=pod

=encoding UTF-8

=head1 NAME

AstroScript::Time - Time-related routines

=head1 VERSION

Version 0.01

=head1 SYNOPSIS

    use AstroScript::Time qw/:all/;

    # Convert Gregorian (new-style) date to old-style date
    my $j = cal2jd(1799, 6, 6); # Standard Julian date of A.Pushkin's birthday
    my $d = jd2cal($j, gregorian => 0) # (1799, 5, 25) = May 26, 1799.

    # Julian date in centuries since epoch 2000.0
    my $t = jd_cent($j); # -2.0056810403833
    ...

=head1 DESCRIPTION

Library of date/time manipulation routines for practical astronomy. Most of them
are based on so called I<Julian date (JD)>, which is the number of days elapsed
since mean UT noon of B<January 1st 4713 BC>. This system of time measurement is
widely adopted by the astronomers.

=head2 JD and MJD

Many routines use Modified Julian date, which starts at B<2000 January 0>
(2000 January 1.0) as the starting point.

=head2 Civil year vs astronomical year

There is disagreement between astronomers and historians about how to count the
years preceding the year 1. Astronomers generally use zero-based system. The
year before the year +1, is the year zero, and the year preceding the latter is
the year -1. The year which the historians call 585 B.C. is actually the year
-584.

In this module all subroutines accepting year assume that B<there is year zero>.
Thus, the sequence of years is: C<BC -3, -2, -1, 0, 1, 2, 3, AD>.

=head2 Date and Time

Time is represented by fractional part of a day. For example, 7h30m UT
is C<(7 + 30 / 60) / 24 = 0.3125>.

=head3 Gregorian calendar

I<Civil calendar> in most cases means I<proleptic Gregorian calendar>. it is
assumed that Gregorian calendar started at Oct. 4, 1582, when it was first
adopted in several European countries. Many other countries still used the older
Julian calendar. In Soviet Russia, for instance, Gregorian system was accepted
on Jan 26, 1918. See:
L<https://en.wikipedia.org/wiki/Gregorian_calendar#Adoption_of_the_Gregorian_Calendar>


=head1 EXPORTED CONSTANTS

=over

=item * C<$SEC_PER_DAY> seconds per day (86400)

=item * C<$SEC_PER_CEN> seconds per century (3155760000)

=item * C<$J2000> Standard Julian date for start of epoch 2000,0 (2451545)

=item * C<$J1900> Standard Julian date for start of epoch 1900,0 (2415020)

=item * C<$GREGORIAN_START> Start of Gregorian calendar, YYYYMMDD (15821004)

=item * C<$JD_UNIX_EPOCH> Standard Julian date for start of the Unix epoch

=back

=head1 EXPORTED FUNCTIONS

=over

=item * L</jd_cent($jd)>

=item * L</after_gregorian($year, $month, $date)>

=item * L</cal2jd($year, $month, $date)>

=item * L</jd2cal($jd)>

=item * L</jd0($jd)>

=item * L</unix2jd($seconds)>

=item * L</jd2unix($jd)>

=item * L</jdnow()>

=item * L</jd2mjd($jd)>

=item * L</mjd2jd($mjd)>

=item * L</jd_cent($jd)>

=item * L</t1900($jd)>

=item * L</t2dt($t)>

=item * L</jd2dt($jd)>

=item * L</jd2te($jd)>

=item * L</jd2gst($jd)>

=item * L</jd2lst($jd, $lng)>

=back

=head1 FUNCTIONS

=head2 jd_cent($jd)

Convert Standard Julian Date to centuries passed since epoch 2000.0

=head2 after_gregorian($year, $month, $date, gregorian_start => $YYYYMMDD )

Does the given date fall to period after Gregorian calendar?

=head3 Positional Arguments

=over

=item * B<year> (astronomic, zero-based)

=item * B<month> (1-12)

=item * B<date> UTC date (1-31) with hours and minutes as decimal part

=back

=head3 Optional Named Arguments

=over

=item *

B<gregorian_start> — start of Gregorian calendar. Default value is
B<15821004> If the date is Julian ("old style"), use C<undef> value.
To provide non-standard start of Gregorian calendar, provide a number
in format YYYYMMDDD, e.g. C<19180126> for Jan 26, 1918.

=back

=head3 Returns

I<true> or I<false>=.

=head2 cal2jd($year, $month, $date)

Convert civil date/time to Standard Julian date.

If C<gregorian_start> argument is not provided, it is assumed that this is a date
of I<Proleptic Gregorian calendar>, which started at Oct. 4, 1582.

=head3 Positional Arguments:

=over

=item * B<year> (astronomic, zero-based)

=item * B<month> (1-12)

=item * B<date> UTC date (1-31) with hours and minutes as decimal part

=back

=head3 Optional Named Arguments

=over

=item *

gregorian_start — start of Gregorian calendar. Default value is
B<15821004> If the date is Julian ("old style"), use C<undef> value.
To provide non-standard start of Gregorian calendar, provide a number
in format YYYYMMDDD, e.g. C<19180126> for Jan 26, 1918.

=back

=head3 Returns

Standard Julian date

=head2 jd2cal($jd)

Convert Standard Julian date to civil date/time

=head3 Positional Arguments

Standard Julian Date

=head3 Optional Named Arguments

=over

=item * gregorian — if i<true>, the result will be old-style (Julian) date

=back

=head3 Returns

A list corresponding to the input values of L</cal2jd($year, $month, $date)> function.
The date is given in the proleptic Gregorian calendar system unless B<gregorian>
flag is set to I<true>.

=head2 jd0($jd)

Given Standard Julian Date, calculate Standard Julian date for midnight of the same date.

=head2 unix2jd($seconds)

Given Unix time, in seconds, convert it to Standard Julian date.

=head2 jd2unix($jd)

Given a Standard Julian Date, convert it to Unix time, in seconds since start of
Unix epoch.

If JD falls before start of the epoch, result will be negative and thus, unusable
for Unix-specific functions like B<localtime()>.

=head2 jdnow()

Standard Julian date for the current moment.

=head2 jd2mjd($jd)

Standard to Modified Julian date.

=head2 mjd2jd($mjd)

Modified to Standard Julian date.

=head2 jd_cent($jd)

Given aI<Standard Julian date>, calculate time in centuries since epoch 2000.0.

=head2 t1900($jd)

Given a I<Standard Julian date>, calculate time in centuries since epoch 2000.0.

=head2 t2dt($t)

Calculate I<Delta-T> in seconds — difference between I<UTC> and I<Ephemeris (Dynamic)time>.
The approximation spans the years from 1825 to 2005

=head3 Argument

Time in centuries since epoch 2000.0

=head3 Returns

=over

=item * Delta-T, in seconds, or zero if given time is out of range

=item * flag: true if given time is in range and Delta-T was calculated.

=back

=head2 jd2dt($jd)

Same as L<t2dt($t)>, but argument is Standard Julian date.

=head2 jd2te($jd)

Given Standard Julian date, calculate T-period corrected for Delta-T.

=head3 Returns

=over

=item * Time in centuries since epoch 2000.0, in Ephemeris (Dynamic) time.

=item * flag: true if given time is in range and Delta-T was calculated.

=back

=head2 jd2gst($jd)

Given I<Standard Julian date>, calculate I<True Greenwich Sidereal time>.

=head2 jd2lst($jd, $lng)

Givan I<Standard Julian date>, calculate true I<Local Sidereal time>.

=head3 Arguments

=over

=item * $jd — Standard Julian date

=item * $lng — Geographic longitude, negative for Eastern longitude, 0 by default

=back


=head1 AUTHOR

Sergey Krushinsky, C<< <krushinsky at gmail.com> >>

=head1 LICENSE AND COPYRIGHT

Copyright 2010-2019 Sergey Krushinsky.

This program is free software; you can redistribute it and/or modify it
under the terms of either: the GNU General Public License as published
by the Free Software Foundation; or the Artistic License.

See http://dev.perl.org/licenses/ for more information.

=cut