Astro-satpass
view release on metacpan or search on metacpan
lib/Astro/Coord/ECI/Utils.pm view on Meta::CPAN
=head1 NAME
Astro::Coord::ECI::Utils - Utility routines for astronomical calculations
=head1 SYNOPSIS
use Astro::Coord::ECI::Utils qw{ :all };
my $now = time ();
print "The current Julian day is ", julianday ($now);
=head2 UN-DEPRECATION NOTICE
In version 0.131, C<date2epoch()> and C<epoch2datetime()> were
deprecated in favor of C<Time::Local::timegm_posix()> and
C<CORE::gmtime()>.
I later realized that this was wrong. The problem with it was that
astronomical dates are often given in the Julian calendar before October
15 1582, which the previously-recognized subroutines do not handle.
As of version 0.132 the deprecation of these subroutines is retracted,
and all reference to it (other than this notice) will be removed.
=head1 DESCRIPTION
This module was written to provide a home for all the constants and
utility subroutines used by B<Astro::Coord::ECI> and its descendants.
What ended up here was anything that was essentially a subroutine, not
a method.
Because figuring out how to convert to and from Perl time bids fair to
become complicated, this module is also responsible for figuring out how
that is done, and exporting whatever is needful to export. See C<:time>
below for the gory details.
This module exports nothing by default. But all the constants,
variables, and subroutines documented below are exportable, and the
following export tags may be used:
=over
=item :all
This imports everything exportable into your name space.
=item :greg_time
This imports all time routines except the discouraged routines
C<time_gm()> and C<time_local()>.
=item :mainstream
This imports everything not deprecated into your name space.
=item :params
This imports the parameter validation routines C<__classisa> and
C<__instance>.
=item :ref
This imports all the C<*_REF> constants.
=item :time
This imports the time routines into your name space. If
L<Time::y2038|Time::y2038> is available, it will be loaded, and both
this tag and C<:all> will import C<gmtime>, C<localtime>, C<time_gm>,
C<time_local>, C<greg_time_gm>, and C<greg_time_local> into your name
space. Otherwise, C<Time::Local|Time::Local> will be loaded, and both
this tag and C<:all> will import C<time_gm>, C<time_local>,
C<greg_time_gm>, and C<greg_time_local> into your name space.
=item :vector
This imports the vector arithmetic routines. It includes anything whose
name begins with C<'vector_'>.
=back
Under Perl 5.6 you may find that, if you use any of the above tags, you
need to specify them first in your import list.
=head2 The following constants are exportable:
AU = number of kilometers in an astronomical unit
JD_OF_EPOCH = the Julian Day of Perl epoch 0
LIGHTYEAR = number of kilometers in a light year
PARSEC = number of kilometers in a parsec
PERL2000 = January 1 2000, 12 noon universal, in Perl time
PI = the circle ratio, computed as atan2 (0, -1)
PIOVER2 = half the circle ratio
SECSPERDAY = the number of seconds in a day
SECS_PER_SIDERIAL_DAY = seconds in a siderial day
SECS_PER_TROPICAL_YEAR = seconds in a tropical year
SPEED_OF_LIGHT = speed of light in kilometers per second
TWOPI = twice the circle ratio
ARRAY_REF = 'ARRAY'
CODE_REF = 'CODE'
HASH_REF = 'HASH'
SCALAR_REF = 'SCALAR'
=head2 The following global variables are exportable:
=head3 $DATETIMEFORMAT
This variable represents the POSIX::strftime format used to convert
times to strings. The default value is '%a %b %d %Y %H:%M:%S' to be
consistent with the behavior of gmtime (or, to be precise, the
behavior of ctime as documented on my system).
=head3 $JD_GREGORIAN
This variable represents the Julian Day of the switch from Julian to
Gregorian calendars. This is used by date2jd(), jd2date(), and the
routines which depend on them, for deciding whether the date is to be
interpreted as in the Julian or Gregorian calendar. Its initial setting
is 2299160.5, which represents midnight October 15 1582 in the Gregorian
calendar, which is the date that calendar was first adopted. This is
slightly different than the value of 2299161 (noon of the same day) used
by Jean Meeus.
If you are interested in historical calculations, you may wish to reset
this appropriately. If you use date2jd to calculate the new value, be
aware of the effect the current setting of $JD_GREGORIAN has on the
interpretation of the date you give.
=head2 In addition, the following subroutines are exportable:
=over 4
=cut
package Astro::Coord::ECI::Utils;
use strict;
use warnings;
our $VERSION = '0.134';
our @ISA = qw{Exporter};
use Carp;
## use Config;
## use Data::Dumper;
use POSIX qw{ floor modf strftime };
use Scalar::Util qw{ blessed };
my @greg_time_routines;
BEGIN {
# NOTE WELL
#
# The logic here should be consistent with the optional-module text
# emitted by inc/Astro/Coord/ECI/Recommend.pm.
#
eval {
require Time::y2038;
Time::y2038->import( qw{ gmtime localtime } );
# sub time_gm
*time_gm = sub {
my @date = @_;
$date[5] = _year_adjust_y2038( $date[5] );
return Time::y2038::timegm( @date );
};
# greg_time_local
*greg_time_gm = sub {
my @date = @_;
$date[5] -= 1900;
return Time::y2038::timegm( @date );
};
# sub time_local
*time_local = sub {
my @date = @_;
$date[5] = _year_adjust_y2038( $date[5] );
return Time::y2038::timelocal( @date );
lib/Astro/Coord/ECI/Utils.pm view on Meta::CPAN
foreach my $mag ( @arg ) {
$sum += 10 ** ( -0.4 * $mag );
}
return -2.5 * log( $sum ) / log( 10 );
}
=item $angle = asin ($value)
This subroutine calculates the arc in radians whose sine is the given
value.
=cut
sub asin {return atan2 ($_[0], sqrt (1 - $_[0] * $_[0]))}
=item $magnitude = atmospheric_extinction ($elevation, $height);
This subroutine calculates the typical atmospheric extinction in
magnitudes at the given elevation above the horizon in radians and the
given height above sea level in kilometers.
The algorithm comes from Daniel W. E. Green's article "Magnitude
Corrections for Atmospheric Extinction", which was published in
the July 1992 issue of "International Comet Quarterly", and is
available online at
L<http://www.icq.eps.harvard.edu/ICQExtinct.html>. The text of
this article makes it clear that the actual value of the
atmospheric extinction can vary greatly from the typical
values given even in the absence of cloud cover.
=cut
# Note that the "constant" 0.120 in Aaer (aerosol scattering) is
# based on a compromise value A0 = 0.050 in Green's equation 3
# (not exhibited here), which can vary from 0.035 in the winter to
# 0.065 in the summer. This makes a difference of a couple tenths
# at 20 degrees elevation, but a couple magnitudes at the
# horizon. Green also remarks that the 1.5 denominator in the
# same equation (a.k.a. the scale height) can be up to twice
# that.
sub atmospheric_extinction {
my ($elevation, $height) = @_;
my $cosZ = cos (PIOVER2 - $elevation);
my $X = 1/($cosZ + 0.025 * exp (-11 * $cosZ)); # Green 1
my $Aray = 0.1451 * exp (-$height / 7.996); # Green 2
my $Aaer = 0.120 * exp (-$height / 1.5); # Green 4
return ($Aray + $Aaer + 0.016) * $X; # Green 5, 6
}
=item $jd = date2jd ($sec, $min, $hr, $day, $mon, $yr)
This subroutine converts the given date to the corresponding Julian day.
The inputs are a Perl date and time; $mon is in the range 0 -
11, and $yr is from 1900, with earlier years being negative. The year 1
BC is represented as -1900.
If less than 6 arguments are provided, zeroes will be prepended to the
argument list as needed.
The date is presumed to be in the Gregorian calendar. If the resultant
Julian Day is before $JD_GREGORIAN, the date is reinterpreted as being
from the Julian calendar.
The only validation is that the month be between 0 and 11 inclusive, and
that the year be not less than -6612 (4713 BC). Fractional days are
accepted.
The algorithm is from Jean Meeus' "Astronomical Algorithms", second
edition, chapter 7 ("Julian Day"), pages 60ff, but the month is
zero-based, not 1-based, and years are 1900-based.
=cut
our $DATETIMEFORMAT;
our $JD_GREGORIAN;
BEGIN {
$DATETIMEFORMAT = '%a %b %d %Y %H:%M:%S';
$JD_GREGORIAN = 2299160.5;
}
sub date2jd {
my @args = @_;
unshift @args, 0 while @args < 6;
my ($sec, $min, $hr, $day, $mon, $yr) = @args;
$mon++; # Algorithm expects month 1-12.
$yr += 1900; # Algorithm expects zero-based year.
($yr < -4712) and croak "Error - Invalid year $yr";
($mon < 1 || $mon > 12) and croak "Error - Invalid month $mon";
if ($mon < 3) {
--$yr;
$mon += 12;
}
my $A = floor ($yr / 100);
my $B = 2 - $A + floor ($A / 4);
my $jd = floor (365.25 * ($yr + 4716)) +
floor (30.6001 * ($mon + 1)) + $day + $B - 1524.5 +
((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24;
$jd < $JD_GREGORIAN and
$jd = floor (365.25 * ($yr + 4716)) +
floor (30.6001 * ($mon + 1)) + $day - 1524.5 +
((($sec || 0) / 60 + ($min || 0)) / 60 + ($hr || 0)) / 24;
return $jd;
}
use constant JD_OF_EPOCH => date2jd (gmtime (0));
=item $epoch = date2epoch ($sec, $min, $hr, $day, $mon, $yr)
This is a convenience routine that converts the given date to seconds
since the epoch, going through date2jd() to do so. The arguments are the
same as those of date2jd().
If less than 6 arguments are provided, zeroes will be prepended to the
argument list as needed.
The functionality is similar to C<Time::Local::timegm()>, but the
arguments will be interpreted according to the Julian calendar if the
date is before L<$JD_GREGORIAN|/$JD_GREGORIAN>.
=cut
sub date2epoch {
my @args = @_;
__subroutine_deprecation();
unshift @args, 0 while @args < 6;
my ($sec, $min, $hr, $day, $mon, $yr) = @args;
return (date2jd ($day, $mon, $yr) - JD_OF_EPOCH) * SECSPERDAY +
(($hr || 0) * 60 + ($min || 0)) * 60 + ($sec || 0);
}
=item $time = decode_space_track_json_time( $string )
This subroutine decodes a time in the format Space Track uses in their
JSON code. This is ISO-8601-ish, but with a possible fractional part and
no zone.
=cut
sub decode_space_track_json_time {
my ( $string ) = @_;
$string =~ m{ \A \s*
( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+
( [0-9]+ ) [^0-9]+ ( [0-9]+ ) [^0-9]+ ( [0-9]+ ) (?: ( [.] [0-9]* ) )? \s* \z }smx
or return;
my @time = ( $1, $2, $3, $4, $5, $6 );
my $frac = $7;
$time[0] = __tle_year_to_Gregorian_year( $time[0] );
$time[1] -= 1;
my $rslt = greg_time_gm( reverse @time );
defined $frac
and $frac ne '.'
and $rslt += $frac;
return $rslt;
}
# my ( $self, $station, @args ) = __default_station( @_ )
#
# This exportable subroutine checks whether the second argument embodies
# an Astro::Coord::ECI object. If so, the arguments are returned
# verbatim. If not, the 'station' attribute of the invocant is inserted
# after the first argument and the result returned. If the 'station'
# attribute of the invocant has not been set, an exception is thrown.
sub __default_station {
my ( $self, @args ) = @_;
if ( ! embodies( $args[0], 'Astro::Coord::ECI' ) ) {
my $sta = $self->get( 'station' )
or croak 'Station attribute not set';
unshift @args, $sta;
}
return ( $self, @args );
}
=item $rad = deg2rad ($degr)
This subroutine converts degrees to radians. If the argument is
C<undef>, C<undef> will be returned.
lib/Astro/Coord/ECI/Utils.pm view on Meta::CPAN
}
=item $seconds = dynamical_delta ($time);
This method returns the difference between dynamical and universal time
at the given universal time. That is,
$dynamical = $time + dynamical_delta ($time)
if $time is universal time.
The algorithm is from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 10, page 78. Meeus notes that this is actually an
observed quantity, and the algorithm is an approximation.
=cut
sub dynamical_delta {
my ($time) = @_;
my $year = (gmtime $time)[5] + 1900;
my $t = ($year - 2000) / 100;
my $correction = .37 * ($year - 2100); # Meeus' correction to (10.2)
return (25.3 * $t + 102) * $t + 102 # Meeus (10.2)
+ $correction; # Meeus' correction.
}
=item $boolean = embodies ($thingy, $class)
This subroutine represents a safe way to call the 'represents' method on
$thingy. You get back true if and only if $thingy->can('represents')
does not throw an exception and returns true, and
$thingy->represents($class) returns true. Otherwise it returns false.
Any exception is trapped and dismissed.
This subroutine is called 'embodies' because it was too confusing to
call it 'represents', both for the author and for the Perl interpreter.
=cut
sub embodies {
my ($thingy, $class) = @_;
my $code = eval {$thingy->can('represents')};
return $code ? $code->($thingy, $class) : undef;
}
=item ($sec, $min, $hr, $day, $mon, $yr, $wday, $yday, 0) = epoch2datetime ($epoch)
This convenience subroutine converts the given time in seconds from the
system epoch to the corresponding date and time. It is implemented in
terms of jd2date (), with the year and month returned from that
subroutine. The day is a whole number, with the fractional part
converted to hours, minutes, and seconds. The $wday is the day of the
week, with Sunday being 0. The $yday is the day of the year, with
January 1 being 0. The trailing 0 is the summer time (or daylight saving
time) indicator which is always 0 to be consistent with gmtime.
If called in scalar context, it returns the date formatted by
POSIX::strftime, using the format string in $DATETIMEFORMAT.
The functionality is similar to C<CORE::gmtime()>, but the result will
be in the Julian calendar if the date is before
L<$JD_GREGORIAN|/$JD_GREGORIAN>.
The input must convert to a non-negative Julian date. The exact lower
limit depends on the system, but is computed by -(JD_OF_EPOCH * 86400).
For Unix systems with an epoch of January 1 1970, this is -210866760000.
Additional algorithms for day of week and day of year come from Jean
Meeus' "Astronomical Algorithms", 2nd Edition, Chapter 7 (Julian Day),
page 65.
=cut
sub epoch2datetime {
my ($time) = @_;
__subroutine_deprecation();
my $day = floor ($time / SECSPERDAY);
my $sec = $time - $day * SECSPERDAY;
($day, my $mon, my $yr, undef, my $leap) = jd2date (
my $jd = $day + JD_OF_EPOCH);
$day = floor ($day + .5);
my $min = floor ($sec / 60);
$sec = $sec - $min * 60;
my $hr = floor ($min / 60);
$min = $min - $hr * 60;
my $wday = ($jd + 1.5) % 7;
my $yd = floor (275 * ($mon + 1) / 9) - (2 - $leap) * floor (($mon +
10) / 12) + $day - 31;
wantarray and return ($sec, $min, $hr, $day, $mon, $yr, $wday, $yd,
0);
# FIXME this seems to be needed under Perl 5.39.8, but not under
# earlier Perls.
local $ENV{TZ} = 'UTC';
return POSIX::strftime( $DATETIMEFORMAT, $sec, $min, $hr, $day,
$mon, $yr, $wday, $yd );
}
=item $time = find_first_true ($start, $end, \&test, $limit);
This function finds the first time between $start and $end for which
test ($time) is true. The resolution is $limit, which defaults to 1 if
not specified. If the times are reversed (i.e. the start time is after
the end time) the time returned is the last time test ($time) is true.
The C<test()> function is called with the Perl time as its only
argument. It is assumed to be false for the first part of the interval,
and true for the rest. If this assumption is violated, the result of
this subroutine should be considered meaningless.
The calculation is done by, essentially, a binary search; the interval
is repeatedly split, the function is evaluated at the midpoint, and a
new interval selected based on whether the result is true or false.
Actually, nothing in this function says the independent variable has to
be time.
=cut
sub find_first_true {
my ($begin, $end, $test, $limit) = @_;
$limit ||= 1;
lib/Astro/Coord/ECI/Utils.pm view on Meta::CPAN
This wrapper is needed because the routines have subtly different
signatures. L<Time::y2038|Time::y2038> C<timegm()> interprets years
strictly as Perl years. L<Time::Local|Time::Local> C<timegm_modern()>
interprets them strictly as Gregorian years. L<Time::Local|Time::Local>
C<timegm()> interprets them differently depending on the value of the
year. Years greater than or equal to 1000 are Gregorian years, but all
others are Perl years, except for the range 0 to 99 inclusive, which are
within 50 years of the current year.
If you are doing historical calculations, see
L<Historical Calculations|Astro::Coord::ECI::Sun/Historical Calculations>
in the L<Astro::Coord::ECI::Sun|Astro::Coord::ECI::Sun> documentation
for a discussion of input and output time conversion.
=item $epoch = greg_time_local( $sec, $min, $hr, $day, $mon, $yr );
This exportable subroutine is a wrapper for either
C<Time::y2038::timelocal()> (if that module is installed),
C<Time::Local::timelocal_modern()> (if that is available), or
C<Time::Local::timelocal()> (if not.)
This subroutine interprets years as Gregorian years.
The difference between this and c<time_local()> is that C<time_local()>
interprets the year the way C<Time::Local::timelocal()> does. For that
reason, this subroutine is preferred over c<time_local()>.
This wrapper is needed for the same reason C<greg_time_gm()> is
needed.
If you are doing historical calculations, see
L<Historical Calculations|Astro::Coord::ECI::Sun/Historical Calculations>
in the L<Astro::Coord::ECI::Sun|Astro::Coord::ECI::Sun> documentation
for a discussion of input and output time conversion.
=item $difference = intensity_to_magnitude ($ratio)
This function converts a ratio of light intensities to a difference in
stellar magnitudes. The algorithm comes from Jean Meeus' "Astronomical
Algorithms", Second Edition, Chapter 56, Page 395.
Note that, because of the way magnitudes work (a more negative number
represents a brighter star) you get back a positive result for an
intensity ratio less than 1, and a negative result for an intensity
ratio greater than 1.
=cut
{ # Begin local symbol block
my $intensity_to_mag_factor; # Calculate only if needed.
sub intensity_to_magnitude {
return - ($intensity_to_mag_factor ||= 2.5 / log (10)) * log ($_[0]);
}
}
=item ($day, $mon, $yr, $greg, $leap) = jd2date ($jd)
This subroutine converts the given Julian day to the corresponding date.
The returns are year - 1900, month (0 to 11), day (which may have a
fractional part), a Gregorian calendar indicator which is true if the
date is in the Gregorian calendar and false if it is in the Julian
calendar, and a leap (or bissextile) year indicator which is true if the
year is a leap year and false otherwise. The year 1 BC is returned as
-1900 (i.e. as year 0), and so on. The date will probably have a
fractional part (e.g. 2006 1 1.5 for noon January first 2006).
If the $jd is before $JD_GREGORIAN, the date will be in the Julian
calendar; otherwise it will be in the Gregorian calendar.
The input may not be less than 0.
The algorithm is from Jean Meeus' "Astronomical Algorithms", second
edition, chapter 7 ("Julian Day"), pages 63ff, but the month is
zero-based, not 1-based, and the year is 1900-based.
=cut
sub jd2date {
my ($time) = @_;
my $mod_jd = $time + 0.5;
my $Z = floor ($mod_jd);
my $F = $mod_jd - $Z;
my $A = (my $julian = $Z < $JD_GREGORIAN) ? $Z : do {
my $alpha = floor (($Z - 1867216.25)/36524.25);
$Z + 1 + $alpha - floor ($alpha / 4);
};
my $B = $A + 1524;
my $C = floor (($B - 122.1) / 365.25);
my $D = floor (365.25 * $C);
my $E = floor (($B - $D) / 30.6001);
my $day = $B - $D - floor (30.6001 * $E) + $F;
my $mon = $E < 14 ? $E - 1 : $E - 13;
my $yr = $mon > 2 ? $C - 4716 : $C - 4715;
return ($day, $mon - 1, $yr - 1900, !$julian, ($julian ? !($yr % 4) : (
$yr % 400 ? $yr % 100 ? !($yr % 4) : 0 : 1)));
## % 400 ? 1 : $yr % 100 ? 0 : !($yr % 4))));
}
=item ($sec, $min, $hr, $day, $mon, $yr, $wday, $yday, 0) = jd2datetime ($jd)
This convenience subroutine converts the given Julian day to the
corresponding date and time. All this really does is converts its
argument to seconds since the system epoch, and pass off to
epoch2datetime().
The input may not be less than 0.
=cut
sub jd2datetime {
my ($time) = @_;
return epoch2datetime(($time - JD_OF_EPOCH) * SECSPERDAY);
}
=item $century = jcent2000 ($time);
Several of the algorithms in Jean Meeus' "Astronomical Algorithms"
are expressed in terms of the number of Julian centuries from epoch
J2000.0 (e.g equations 12.1, 22.1). This subroutine encapsulates
that calculation.
=cut
sub jcent2000 {return jday2000 ($_[0]) / 36525}
=item $jd = jday2000 ($time);
This subroutine converts a Perl date to the number of Julian days
( run in 1.018 second using v1.01-cache-2.11-cpan-99c4e6809bf )