=head1 NAME
Astro::Coord::ECI - Manipulate geocentric coordinates
=head1 SYNOPSIS
use Astro::Coord::ECI;
use Astro::Coord::ECI::Sun;
use Astro::Coord::ECI::TLE;
use Astro::Coord::ECI::Utils qw{rad2deg};
# 1600 Pennsylvania Avenue, in radians, radians, and KM
my ($lat, $lon, $elev) = (0.678911227503559,
-1.34456123391096, 0.01668);
# Record the time
my $time = time ();
# Set up observer's location
my $loc = Astro::Coord::ECI->geodetic ($lat, $lon, $elev);
# Uncomment to turn off atmospheric refraction if desired.
# $loc->set( refraction => 0 );
# Instantiate the Sun.
my $sun = Astro::Coord::ECI::Sun->universal ($time);
# Figure out if the Sun is up at the observer's location.
my ($azimuth, $elevation, $range) = $loc->azel ($sun);
print "The Sun is ", rad2deg ($elevation),
" degrees above the horizon.\n";
See the L<Astro::Coord::ECI::TLE|Astro::Coord::ECI::TLE> documentation
for an example involving satellite pass prediction.
=head1 NOTICE
In version [% next_version %] I very belatedly took cognizance of the
fact that in the Southern hemisphere the seasons are reversed. The
result is an incompatible change in the text description of the equinox
and solstice events generated by the
L<Astro::Coord::ECI::Sun|Astro::Coord::ECI::Sun> object if its
C<station> attribute is set to a location south of the Equator. If this
attribute is not set at all you get the Northern-hemisphere text.
=head1 DESCRIPTION
This module was written to provide a base class for a system to predict
satellite visibility. Its main task is to convert the
Earth-Centered Inertial (ECI) coordinates generated by the NORAD models
into other coordinate systems (e.g. latitude, longitude, and altitude
above mean sea level), but a other functions have accreted to it.
In addition, a few support routines have been exposed for testing, or
whatever.
All distances are in kilometers, and all angles are in radians
(including right ascension, which is usually measured in hours).
Times are normal Perl times, whether used as universal or dynamical
time. Universal time is what is usually meant, unless otherwise stated.
Known subclasses include B<Astro::Coord::ECI::Moon> to predict the
position of the Moon, B<Astro::Coord::ECI::Star> to predict the
position of a star, or anything else that can be considered fixed on
the celestial sphere, B<Astro::Coord::ECI::Sun> to predict the position
of the Sun, B<Astro::Coord::ECI::TLE> to predict the position of a
satellite given the NORAD orbital parameters, and
B<Astro::Coord::ECI::TLE::Iridium> (a subclass of
Astro::Coord::ECI::TLE) to predict Iridium flares.
B<Note> that in version 0.022_01 the velocity code got a substantial
rework, which is still in progress. I am attempting give useful
velocities where they are available, and no velocities at all where they
are not. Unfortunately I have yet to locate any worked problems, so the
velocity code is, in the most literal sense, untested.
In general, locations specified in Earth-fixed coordinates are
considered to share the rotational velocity of the Earth, and locations
specified in inertial coordinates are considered to have undefined
velocities unless a velocity was specified explicitly or produced by the
object's model. This involved a change in the behavior of the eci()
method when velocity was not specified but location was. See the next
paragraph for my excuse.
B<Caveat user:> This class and its subclasses should probably be
considered alpha code, meaning that the public interface may not be
completely stable. I will try not to change it, but given sufficient
reason I will do so. If I do so, I will draw attention to the change
in the documentation.
=head2 Methods
The following methods should be considered public.
=over 4
=cut
package Astro::Coord::ECI;
use strict;
use warnings;
our $VERSION = '0.093';
use Astro::Coord::ECI::Utils qw{:all};
use Carp;
use Data::Dumper;
use POSIX qw{floor strftime};
use Storable ();
use constant NO_CASCADING_STATIONS =>
q{Cascading 'station' attributes are not supported};
use constant TWO_DEGREES => deg2rad( 2 );
my %mutator; # Attribute mutators. We define these just after the
# set() method, for convenience.
my %known_ellipsoid; # Known reference ellipsoids. We define these
# just before the reference_ellipsoid() method for
# convenience.
my %static = ( # The geoid, etc. Geoid gets set later.
# almanac_horizon => set when instantiated, so that the following
# _almanac_horizon => gets set automatically by the mutator.
angularvelocity => 7.292114992e-5, # Of surface of Earth, 1998. Meeus, p.83
debug => 0,
diameter => 0,
edge_of_earths_shadow => 1,
horizon => deg2rad (20),
refraction => 1,
twilight => deg2rad (-6),
);
my %savatr; # Attribs saved across "normal" operations. Set at end.
my @kilatr = # Attributes to purge when setting coordinates.
qw{_need_purge _ECI_cache inertial local_mean_time specified}; #?
=item $coord = Astro::Coord::ECI->new ();
This method instantiates a coordinate object. Any arguments are passed
to the set() method once the object has been instantiated.
=cut
sub new {
my ( $class, @args ) = @_;
my $self = bless { %static }, ref $class || $class;
$self->{inertial} = $self->__initial_inertial();
@args and $self->set( @args );
exists $self->{almanac_horizon}
or $self->set( almanac_horizon => 0 );
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug - Instantiated ", Dumper ($self);
};
return $self;
}
=item $angle = $coord->angle ($coord2, $coord3);
This method calculates the angle between the $coord2 and $coord3
objects, as seen from $coord. The calculation uses the law of
haversines, and does not take atmospheric refraction into account. The
return is a number of radians between 0 and pi.
The algorithm comes from "Ask Dr. Math" on Drexel's Math Forum,
L<http://mathforum.org/library/drmath/view/51879.html>, which attributes
it to the Geographic Information Systems FAQ,
L<http://www.faqs.org/faqs/geography/infosystems-faq/>, which in turn
attributes it to R. W. Sinnott, "Virtues of the Haversine," Sky and
Telescope, volume 68, number 2, 1984, page 159.
Prior to version 0.011_03 the law of cosines was used, but this produced
errors on extremely small angles. The haversine law was selected based
on Jean Meeus, "Astronomical Algorithms", 2nd edition, chapter 17
"Angular Separation".
This method ignores the C<station> attribute.
=cut
sub angle {
my $self = shift;
my $B = shift;
my $C = shift;
(ref $B && $B->represents (__PACKAGE__)
&& ref $C && $C->represents (__PACKAGE__))
or croak <<eod;
Error - Both arguments must represent @{[__PACKAGE__]} objects.
eod
my ($ra1, $dec1) = $self->{inertial} ?
$B->equatorial ($self) : $self->_angle_non_inertial ($B);
my ($ra2, $dec2) = $self->{inertial} ?
$C->equatorial ($self) : $self->_angle_non_inertial ($C);
my $sindec = sin (($dec2 - $dec1)/2);
my $sinra = sin (($ra2 - $ra1)/2);
my $a = $sindec * $sindec +
cos ($dec1) * cos ($dec2) * $sinra * $sinra;
return 2 * atan2 (sqrt ($a), sqrt (1 - $a));
}
sub _angle_non_inertial {
my $self = shift;
my $other = shift;
my ($x1, $y1, $z1) = $self->ecef ();
my ($x2, $y2, $z2) = $other->ecef ();
my $X = $x2 - $x1;
my $Y = $y2 - $y1;
my $Z = $z2 - $z1;
my $lon = atan2 ($Y, $X);
my $lat = mod2pi (atan2 ($Z, sqrt ($X * $X + $Y * $Y)));
return ($lon, $lat);
}
=item $which = $coord->attribute ($name);
This method returns the name of the class that implements the named
attribute, or undef if the attribute name is not valid.
=cut
sub attribute {return exists $mutator{$_[1]} ? __PACKAGE__ : undef}
=item ($azimuth, $elevation, $range) = $coord->azel( $coord2 );
This method takes another coordinate object, and computes its azimuth,
elevation, and range in reference to the object doing the computing.
The return is azimuth in radians measured clockwise from North (always
positive), elevation above the horizon in radians (negative if
below), and range in kilometers.
As a side effect, the time of the $coord object may be set from the
$coord2 object.
Atmospheric refraction is taken into account using the
C<correct_for_refraction()> method if the C<$coord> object's
C<refraction> attribute is true. The invocant's
C<correct_for_refraction()> method is the one used; that is, if
C<$coord> and C<$coord2> have different C<correct_for_refraction()>
methods, the C<$coord> object's method is used.
This method is implemented in terms of azel_offset(). See that method's
documentation for further details.
=item ( $azimuth, $elevation, $range ) = $coord->azel();
This method computes the azimuth, elevation, and range if the C<$coord>
object as seen from the position stored in the C<$coord> object's
C<station> attribute. An exception will be thrown if the C<station>
attribute is not set.
=cut
sub azel { ## no critic (RequireArgUnpacking)
@_ > 2
and croak q{Too many arguments};
@_ = _expand_args_default_station( @_ );
$_[2] = $_[2] ? 1 : 0;
goto &azel_offset;
}
=item ($azimuth, $elevation, $range) = $coord->azel_offset($coord2, $offset);
This method takes another coordinate object, and computes its azimuth,
elevation, and range in reference to the object doing the computing.
The return is azimuth in radians measured clockwise from North (always
positive), elevation above the horizon in radians (negative if
below), and range in kilometers.
If the optional C<$offset> argument is provided, the elevation is offset
upward by the given fraction of the radius of the C<$coord2> object.
Thus, an C<$offset> of C<1> specifies the upper limb of the object, C<0>
specifies the center of the object, and C<-1> specifies the lower limb.
This offset is applied before atmospheric refraction.
As a side effect, the time of the $coord object may be set from the
$coord2 object.
By default, atmospheric refraction is taken into account in the
calculation of elevation, using the C<correct_for_refraction()> method.
This better represents the observed position in the sky when the object
is above the horizon, but produces a gap in the data when the object
passes below the horizon, since I have no refraction equations for rock.
See the C<correct_for_refraction()> documentation for details.
The invocant's C<correct_for_refraction()> method is the one used; that
is, if C<$coord> and C<$coord2> have different
C<correct_for_refraction()> methods, the C<$coord> object's method is
used.
If you want to ignore atmospheric refraction (and not have a gap in your
data), set the L<refraction|/refraction> attribute of the $coord object
to any value Perl considers false (e.g. 0).
The algorithm for position is the author's, but he confesses to having
to refer to T. S. Kelso's "Computers and Satellites"
column in "Satellite Times", November/December 1995, titled "Orbital
Coordinate Systems, Part II" and available at
F<http://celestrak.com/columns/v02n02/> to get the signs right.
If velocities are available for both bodies involved in the calculation,
they will be returned after the position (i.e. you get a six-element
array back instead of a three-element array). The velocity of a point
specified in Earth-fixed coordinates (e.g. geodetic latitude, longitude,
and altitude) is assumed to be the rotational velocity of the Earth at
that point. The returned velocities are azimuthal velocity in radians
per second (B<not> radians of azimuth, which get smaller as you go
farther away from the Equator), elevational velocity in radians per
second, and velocity in recession in kilometers per second, in that
order.
If velocities are available for both bodies B<and> the C<$coord2> object
has its C<frequency> attribute set, the returned array will contain
seven elements, with the seventh being the Doppler shift.
The algorithm for recessional velocity comes from John A. Magliacane's
C<Predict> program, available at
L<http://www.qsl.net/kd2bd/predict.html>. The transverse velocity
computations are the author's, but use the same basic approach: vector
dot product of the velocity vector with a unit vector in the appropriate
direction, the latter generated by the appropriate (I hope!) vector
cross products.
Velocity conversions to C<azel()> appear to me to be mostly sane. See
L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below, for details.
If velocities are available I<and> you have provided a non-zero value
for the C<frequency> attribute, you will get the Doppler shift as the
seventh element of the returned array. The I<caveats> about velocity in
recession apply to the Doppler shift as well. The frequency can come
from either the C<$coord> or C<$coord2> object, but the C<$coord2> is
preferred, and getting frequency from the C<$coord> object is being put
through a deprecation cycle and removed.
=item ( $azimuth, $elevation, $range ) = $coord->azel_offset( $offset );
This method computes the azimuth, elevation, and range if the C<$coord>
object as seen from the location stored in the C<$coord> object's
C<station> attribute. The functionality is as above, except for the fact
that in the above version the station is the invocant, whereas in this
version the orbiting body is the invocant.
This version also returns velocities if they are available for both
bodies, and Doppler shift if in addition the C<frequency> attribute of
the invocant is set.
=cut
sub azel_offset {
my ( $self, $trn2, $offset ) = _expand_args_default_station( @_ );
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug azel_offset - ", Dumper ($self, $trn2, $offset);
};
# _local_cartesian() returns NEU coordinates. Converting these to
# spherical gets Azimuth (clockwise from North), Elevation, and
# Range.
my ( $azimuth, $elevation, $range, @velocity ) =
_convert_cartesian_to_spherical(
$self->_local_cartesian( $trn2 )
);
# If the velocity and frequency are defined, we provide the Doppler
# shift as well.
if ( defined $velocity[2] ) {
my $freq = $trn2->get( 'frequency' );
if ( not defined $freq ) {
$freq = $self->get( 'frequency' );
defined $freq
and croak 'Calculation of Doppler shift based ',
'on the frequency attribute of the observing ',
'station is not allowed';
}
if ( defined $freq ) {
$velocity[3] = - $freq * $velocity[2] / SPEED_OF_LIGHT;
}
}
# Adjust for upper limb and refraction if needed.
$offset
and $elevation += atan2(
$trn2->get( 'diameter' ) * $offset / 2,
$range,
);
$self->{refraction} and
$elevation = $self->correct_for_refraction( $elevation );
return ( $azimuth, $elevation, $range, @velocity );
}
=item $coord2 = $coord->clone ();
This method does a deep clone of an object, producing a different
but identical object.
At the moment, it's really just a wrapper for Storable::dclone.
=cut
sub clone {
my ( $self ) = @_;
return Storable::dclone( $self );
}
=item $elevation = $coord->correct_for_refraction ($elevation);
This method corrects the given angular elevation for atmospheric
refraction. This is done only if the elevation has a reasonable chance
of being visible; that is, if the elevation before correction is not
more than two degrees below either the 'horizon' attribute or zero,
whichever is lesser. Sorry for the discontinuity thus introduced, but I
did not wish to carry the refraction calculation too low because I am
uncertain about the behavior of the algorithm, and I do not have a
corresponding algorithm for refraction through magma.
This method can also be called as a class method, in which case the
correction is applied only if the uncorrected elevation is greater than
minus two degrees. It is really only exposed for testing purposes (hence
the cumbersome name). The azel() method calls it for you if the
L<refraction|/refraction> attribute is true.
The algorithm for atmospheric refraction comes from Thorfinn
Saemundsson's article in "Sky and Telescope", volume 72, page 70
(July 1986) as reported Jean Meeus' "Astronomical Algorithms",
2nd Edition, chapter 16, page 106, and includes the adjustment
suggested by Meeus.
=cut
sub correct_for_refraction {
my $self = shift;
my $elevation = shift;
# If called as a static method, our effective horizon is zero. If
# called as a normal method, our effective horizon is the lesser of
# the 'horizon' setting or zero.
my $horizon = ref $self ? min( 0, $self->get( 'horizon' ) ) : 0;
# We exclude anything with an elevation <= 2 degrees below our
# effective horizon; this is presumed to be not visible, since the
# maximum deflection is about 35 minutes of arc. This is not
# portable to (e.g.) Venus.
if ( $elevation > $horizon - TWO_DEGREES ) {
# Thorsteinn Saemundsson's algorithm for refraction, as reported
# in Meeus, page 106, equation 16.4, and adjusted per the
# suggestion in Meeus' following paragraph. Thorsteinn's
# formula is in terms of angles in degrees and produces
# a correction in minutes of arc. Meeus reports the original
# publication as Sky and Telescope, volume 72 page 70, July 1986.
# In deference to Thorsteinn I will point out:
# * The Icelanders do not use family names. The "Saemundsson"
# simply means his father's name was Saemund.
# * I have transcribed the names into 7-bit characters.
# "Thorsteinn" actually does not begin with "Th" but with
# the letter "Thorn." Similarly, the "ae" in "Saemund" is
# supposed to be a ligature (i.e. squished into one letter).
my $deg = rad2deg ($elevation);
my $correction = 1.02 / tan (deg2rad ($deg + 10.3/($deg + 5.11))) +
.0019279;
$self->get ('debug') and print <<eod;
Debug correct_for_refraction
input: $deg degrees of arc
correction: $correction minutes of arc
eod
# End of Thorsteinn's algorithm.
$correction = deg2rad ($correction / 60);
$elevation += $correction;
}
return $elevation;
}
=item $angle = $coord->dip ();
This method calculates the dip angle of the horizon due to the
altitude of the body, in radians. It will be negative for a location
above the surface of the reference ellipsoid, and positive for a
location below the surface.
The algorithm is simple enough to be the author's.
=cut
sub dip {
my $self = shift;
my (undef, undef, $h) = $self->geodetic ();
my (undef, undef, $rho) = $self->geocentric ();
return $h >= 0 ?
- acos (($rho - $h) / $rho) :
acos ($rho / ($rho - $h));
}
=item $coord = $coord->dynamical ($time);
This method sets the dynamical time represented by the object.
This method can also be called as a class method, in which case it
instantiates the desired object.
The algorithm for converting this to universal time comes from Jean
Meeus' "Astronomical Algorithms", 2nd Edition, Chapter 10, pages 78ff.
=item $time = $coord->dynamical ();
This method returns the dynamical time previously set, or the
universal time previously set, converted to dynamical.
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 10, pages 78ff.
=cut
sub dynamical {
my ($self, @args) = @_;
unless (@args) {
ref $self or croak <<eod;
Error - The dynamical() method may not be called as a class method
unless you specify arguments.
eod
return ($self->{dynamical} ||= $self->{universal} +
dynamical_delta ($self->{universal} || croak <<eod));
Error - Universal time of object has not been set.
eod
}
if (@args == 1) {
my $time = shift @args;
$self = $self->new () unless ref $self;
$self->{_no_set}++; # Supress running the model if any
$self->universal ($time - dynamical_delta ($time));
$self->{dynamical} = $time;
--$self->{_no_set}; # Undo supression of model
$self->_call_time_set (); # Run the model if any
} else {
croak <<eod;
Error - The dynamical() method must be called with either zero
arguments (to retrieve the time) or one argument (to set the
time).
eod
}
return $self;
}
=item $coord = $coord->ecef($x, $y, $z, $xdot, $ydot, $zdot)
This method sets the coordinates represented by the object in terms of
L</Earth-Centered, Earth-fixed (ECEF) coordinates> in kilometers, with
the x axis being latitude 0 longitude 0, the y axis being latitude 0
longitude 90 degrees east, and the z axis being latitude 90 degrees
north. The velocities in kilometers per second are optional, and will
default to zero. ECEF velocities are considered to be relative to the
surface of the Earth; when converting to ECI, the rotational velocity of
the Earth will be added in.
B<Note> that prior to version 0.022_01, the documentation said that the
default velocity was the rotational velocity of the Earth. This was not
correct B<in ECEF coordinates.> The functionality of the code itself in
this respect did not change.
The object itself is returned.
This method can also be called as a class method, in which case it
instantiates the desired object.
=item ($x, $y, $z, $xdot, $ydot, $zdot) = $coord->ecef();
This method returns the object's L</Earth-Centered, Earth-fixed (ECEF)
coordinates>.
If the original coordinate setting was in an inertial system (e.g. eci,
equatorial, or ecliptic) B<and> the absolute difference between the
current value of 'equinox_dynamical' and the current dynamical() setting
is greater than the value of $Astro::Coord::ECI::EQUINOX_TOLERANCE, the
coordinates will be precessed to the current dynamical time before
conversion. Yes, this should be done any time the equinox is not the
current time, but for satellite prediction precession by a year or
so does not seem to make much difference in practice. The default value
of $Astro::Coord:ECI::EQUINOX_TOLERANCE is 365 days. B<Note> that if
this behavior or the default value of $EQUINOX_TOLERANCE begins to look
like a bug, it will be changed, and noted in the documentation.
Some velocity conversions involving C<ecef()> appear to me to be mostly
sane. See L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below, for
details.
=cut
sub ecef {
my ($self, @args) = @_;
$self = $self->_check_coord (ecef => \@args);
unless (@args) {
my $cache = ($self->{_ECI_cache} ||= {});
return @{$cache->{fixed}{ecef}} if $cache->{fixed}{ecef};
return $self->_convert_eci_to_ecef () if $self->{inertial};
croak <<eod;
Error - Object has not been initialized.
eod
}
@args == 3 and push @args, 0, 0, 0;
if (@args == 6) {
foreach my $key (@kilatr) {delete $self->{$key}}
## $self->{_ECI_cache}{fixed}{ecef} = \@args;
$self->{_ECI_cache} = {fixed => {ecef => \@args}};
$self->{specified} = 'ecef';
$self->{inertial} = 0;
} else {
croak <<eod;
Error - The ecef() method must be called with either zero arguments (to
retrieve coordinates), three arguments (to set coordinates,
with velocity defaulting to zero) or six arguments.
eod
}
return $self;
}
=item $coord = $coord->eci ($x, $y, $z, $xdot, $ydot, $zdot, $time)
This method sets the coordinates represented by the object in terms of
L</Earth-Centered Inertial (ECI) coordinates> in kilometers, time being
universal time, the x axis being 0 hours L</Right Ascension> 0 degrees
L</Declination>, y being 6 hours L</Right Ascension> 0 degrees
L</Declination>, and z being 90 degrees north L</Declination>. The
velocities in kilometers per second are optional. If omitted, the object
will be considered not to have a velocity. B<This is a change introduced
with version 0.022_01.> Before that, the velocity defaulted to 0.
The time argument is optional if the time represented by the object
has already been set (e.g. by the universal() or dynamical() methods).
The object itself is returned.
This method can also be called as a class method, in which case it
instantiates the desired object. In this case the time is not optional.
The algorithm for converting from ECI to geocentric coordinates and
back is based on the description of ECI coordinates in T. S. Kelso's
"Computers and Satellites" column in "Satellite Times",
September/October 1995, titled "Orbital Coordinate Systems, Part I"
and available at F<http://celestrak.com/columns/v02n01/>.
=item ($x, $y, $z, $xdot, $ydot, $zdot) = $coord->eci($time);
This method returns the L</Earth-Centered Inertial (ECI) coordinates>
of the object at the given time. The time argument is actually
optional if the time represented by the object has already been set.
If you specify a time, the time represented by the object will be set
to that time. The net effect of specifying a time is equivalent to
($x, $y, $z, $xdot, $ydot, $zdot) = $coord->universal($time)->eci()
If the original coordinate setting was in a non-inertial system (e.g.
ECEF or geodetic), the equinox_dynamical attribute will be set to the
object's dynamical time.
Velocities will be returned if they were originally provided. B<This is
a change introduced in version 0.022_01.> Prior to that version,
velocities were always returned.
Some velocity conversions involving C<eci()> appear to me to be mostly
sane. See L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below, for
details.
=cut
sub eci {
my ($self, @args) = @_;
$self = $self->_check_coord (eci => \@args);
unless (@args) {
my $cache = ($self->{_ECI_cache} ||= {});
return @{$cache->{inertial}{eci}} if $cache->{inertial}{eci};
return $self->_convert_ecef_to_eci () if $self->{specified};
croak <<eod;
Error - Object has not been initialized.
eod
}
## @args == 3 and push @args, 0, 0, 0;
if (@args == 3 || @args == 6) {
foreach my $key (@kilatr) {delete $self->{$key}}
$self->{_ECI_cache} = {inertial => {eci => \@args}};
$self->{specified} = 'eci';
$self->{inertial} = 1;
} else {
croak <<eod;
Error - The eci() method must be called with either zero or one
arguments (to retrieve coordinates), three or four arguments
(to set coordinates, with velocity defaulting to zero), or
six or seven arguments.
eod
}
return $self;
}
=item $coord = $coord->ecliptic_cartesian( $X, $Y, $Z, $time );
This method sets the Cartesian L</Ecliptic> coordinates represented by
the object in terms of C<X> (kilometers toward the vernal equinox), C<Y>
(90 degrees along the ecliptic from C<X>), and C<Z> (toward ecliptic
North). The time is universal time. The object itself is returned.
The time argument is optional if the time represented by the object has
already been set (e.g. by the universal() or dynamical() methods).
This method can also be called as a class method, in which case it
instantiates the desired object. In this case the time is not optional.
B<Caveat:> you can pass in velocities (before the C<$time>) but they are
unsupported, meaning that I can not at this point say whether they will
be transformed sanely, much less correctly. B<Caveat programmer>.
=item ( $X, $Y, $Z) = $coord->ecliptic_cartesian( $time );
This method returns the Cartesian ecliptic coordinates of the object at
the given time. The time is optional if the time represented by the
object has already been set (e.g. by the universal() or dynamical()
methods).
B<Caveat:> velocities are unsupported by this method. That means you may
(or may not, depending on the coordinates originally set) get them back,
but I have not looked into whether they are actually correct.
=cut
sub ecliptic_cartesian {
my ( $self, @args ) = @_;
$self = $self->_check_coord( ecliptic_cartesian => \@args );
if ( @args ) {
@args % 3
and croak <<'EOD';
Error - The ecliptic_cartesian() method must be called with either zero
or one arguments (to retrieve coordinates), or three or four
arguments (to set coordinates). There is currently no six or
seven argument version.
EOD
my $epsilon = obliquity( $self->dynamical() );
my $sin_epsilon = sin $epsilon;
my $cos_epsilon = cos $epsilon;
my @eci;
for ( my $inx = 0; $inx < @args; $inx += 3 ) {
push @eci, $args[ $inx ];
push @eci, - $args[ $inx + 2 ] * $sin_epsilon +
$args[ $inx + 1 ] * $cos_epsilon;
push @eci, $args[ $inx + 2 ] * $cos_epsilon +
$args[ $inx + 1 ] * $sin_epsilon;
}
$self->eci( @eci );
$self->{_ECI_cache}{inertial}{ecliptic_cartesian} = \@args;
$self->{inertial} = 1;
$self->{specified} = 'ecliptic_cartesian';
return $self;
} else {
$self->{_ECI_cache}{inertial}{ecliptic_cartesian}
and return @{
$self->{_ECI_cache}{inertial}{ecliptic_cartesian}
};
my @eci = $self->eci();
my $epsilon = obliquity( $self->dynamical() );
my $sin_epsilon = sin $epsilon;
my $cos_epsilon = cos $epsilon;
my @ecliptic_cartesian;
for ( my $inx = 0; $inx < @eci; $inx += 3 ) {
push @ecliptic_cartesian, $eci[ $inx ];
push @ecliptic_cartesian, $eci[ $inx + 2 ] * $sin_epsilon +
$eci[ $inx + 1 ] * $cos_epsilon;
push @ecliptic_cartesian, $eci[ $inx + 2 ] * $cos_epsilon -
$eci[ $inx + 1 ] * $sin_epsilon;
}
return @{
$self->{_ECI_cache}{inertial}{ecliptic_cartesian} =
\@ecliptic_cartesian
};
}
}
=item $coord = $coord->ecliptic ($latitude, $longitude, $range, $time);
This method sets the L</Ecliptic> coordinates represented by the object
in terms of L</Ecliptic latitude> and L</Ecliptic longitude> in radians,
and the range to the object in kilometers, time being universal time.
The object itself is returned.
The time argument is optional if the time represented by the object has
already been set (e.g. by the universal() or dynamical() methods).
The latitude should be a number between -PI/2 and PI/2 radians
inclusive. The longitude should be a number between -2*PI and 2*PI
radians inclusive. The increased range (one would expect -PI to PI) is
because in some astronomical usages latitudes outside the range + or -
180 degrees are employed. A warning will be generated if either of these
range checks fails.
This method can also be called as a class method, in which case it
instantiates the desired object. In this case the time is not optional.
The algorithm for converting from ecliptic latitude and longitude to
right ascension and declination comes from Jean Meeus'
"Astronomical Algorithms", 2nd Edition, Chapter 13, page 93.
=item ($latitude, $longitude, $range) = $coord->ecliptic ($time);
This method returns the ecliptic latitude and longitude of the
object at the given time. The time is optional if the time represented
by the object has already been set (e.g. by the universal() or
dynamical() methods).
B<Caveat:> velocities are not returned by this method.
=cut
sub ecliptic {
my ($self, @args) = @_;
$self = $self->_check_coord( ecliptic => \@args );
if ( @args ) {
@args % 3
and croak 'Arguments are in threes, plus an optional time';
$self->{_ECI_cache}{inertial}{ecliptic} = \@args;
$self->ecliptic_cartesian(
_convert_spherical_x_to_cartesian( @args ) );
$self->{specified} = 'ecliptic';
$self->{inertial} = 1;
return $self;
} else {
return @{ $self->{_ECI_cache}{inertial}{ecliptic} ||= [
_convert_cartesian_to_spherical_x(
$self->ecliptic_cartesian() )
]
};
}
}
=item $longitude = $coord->ecliptic_longitude();
This method returns the ecliptic longitude of the body at its current
time setting, in radians. It is really just a convenience method, since
it is equivalent to C<< ( $coord->ecliptic() )[1] >>, and in fact that
is how it is implemented.
=cut
sub ecliptic_longitude {
my ( $self ) = @_;
return ( $self->ecliptic() )[1];
}
=item $coord->equatorial ($rightasc, $declin, $range, $time);
This method sets the L</Equatorial> coordinates represented by the
object (relative to the center of the Earth) in terms of
L</Right Ascension> and L</Declination> in radians, and the range to the
object in kilometers, time being universal time. The object itself is
returned.
If C<equatorial()> is called in this way, the C<station> attribute will
be ignored, for historical reasons.
The right ascension should be a number between 0 and 2*PI radians
inclusive. The declination should be a number between -PI/2 and PI/2
radians inclusive. A warning will be generated if either of these range
checks fails.
The time argument is optional if the time represented by the object
has already been set (e.g. by the universal() or dynamical() methods).
This method can also be called as a class method, in which case it
instantiates the desired object. In this case the time is not optional.
You may optionally pass velocity information in addition to position
information. If you do this, the velocity in right ascension (in radians
per second), declination (ditto), and range (in kilometers per second in
recession) are passed after the position arguments, and before the $time
argument if any.
=item ( $rightasc, $declin, $range ) = $coord->equatorial( $time );
This method returns the L</Equatorial> coordinates of the object at the
relative to the center of the Earth. The C<station> attribute is
ignored. The time argument is optional if the time represented by the
object has already been set (e.g. by the universal() or dynamical()
methods).
If velocities are available from the object (i.e. if it is an instance
of Astro::Coord::ECI::TLE) the return will contain velocity in right
ascension, declination, and range, the first two being in radians per
second and the third in kilometers per second in recession.
B<Caveat:> these velocities are believed by me to be sane B<if> they are
derived from ECI coordinates. This method does not support setting
velocities. See L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below,
for details.
=item ($rightasc, $declin, $range) = $coord->equatorial( $coord2 );
This method is retained (for the moment) for historical reasons, but
C<equatorial_apparent()> is preferred.
This method returns the apparent equatorial coordinates of the object
represented by $coord2, as seen from the location represented by
$coord.
As a side effect, the time of the $coord object may be set from the
$coord2 object.
If the L<refraction|/refraction> attribute of the $coord object is
true, the coordinates will be corrected for atmospheric refraction using
the correct_for_refraction() method.
If velocities are available from both objects (i.e. if both objects are
Astro::Coord::ECI::TLE objects) the return will contain velocity in
right ascension, declination, and range, the first two being in radians
per second and the third in kilometers per second in recession.
B<Caveat:> these velocities are believed by me B<not> to be correct.
=cut
sub equatorial {
my ( $self, @args ) = @_;
my $body;
@args
and embodies( $args[0], __PACKAGE__ )
and $body = shift @args;
$self = $self->_check_coord( equatorial => \@args );
my $time;
$body or $time = $self->universal;
if ( @args ) {
@args % 3
and croak 'Arguments must be in threes, with an optional time';
$body
and croak 'You may not set the equatorial coordinates ',
'relative to an observer';
## my ($ra, $dec, $range, @eqvel) = @args;
$args[0] = _check_right_ascension( 'right ascension' => $args[0] );
$args[1] = _check_latitude( declination => $args[1] );
foreach my $key (@kilatr) {delete $self->{$key}}
$self->{_ECI_cache}{inertial}{equatorial} = \@args;
$self->eci(
_convert_spherical_to_cartesian( @args ) );
$self->{specified} = 'equatorial';
$self->{inertial} = 1;
return $self;
} elsif ( $body ) {
return $self->_equatorial_reduced( $body );
} else {
return @{ $self->{_ECI_cache}{inertial}{equatorial} ||= [
_convert_cartesian_to_spherical( $self->eci() )
] };
}
}
=item ($rightasc, $declin, $range) = $coord->equatorial_apparent();
This method returns the apparent equatorial coordinates of the C<$coord>
object as seen from the object in the C<station> attribute of the
C<$coord> object.
This method will return velocities if available, but I have no idea
whether they are correct.
=cut
sub equatorial_apparent {
my ( $self ) = @_;
my $station = $self->get( 'station' )
or croak 'Station attribute is required';
return $station->_equatorial_reduced( $self );
}
=item my ($rasc, $decl, $range, $v_rasc, $v_decl, $v_r) = $coord->equatorial_unreduced($body);
This method computes the unreduced equatorial position of the second ECI
object as seen from the first. If the second argument is undefined,
computes the equatorial position of the first object as seen from the
center of the Earth. Unlike the equatorial() method itself, this method
is an accessor only. This method would probably not be exposed except
for the anticipation of the usefulness of $range and $v_r in satellite
conjunction computations, and the desirability of not taking the time to
make the two atan2() calls that are unneeded in this application.
The 'unreduced' in the name of the method is intended to refer to the
fact that the $rasc and $decl are not the right ascension and
declination themselves, but the arguments to atan2() needed to compute
them, and $v_rasc and $v_decl are in km/sec, rather than being divided
by the range to get radians per second.
This method ignores the C<'station'> attribute.
The returned data are:
$rasc is an array reference, which can be converted to the right
ascension in radians by mod2pi(atan2($rasc->[0], $rasc->[1])).
$decl is an array reference, which can be converted to the declination
in radians by atan2($decl->[0], $decl->[1]).
$range is the range in kilometers.
$v_rasc is the velocity in the right ascensional direction in kilometers
per second. It can be converted to radians per second by dividing by
$range.
$v_decl is the velocity in the declinational direction in kilometers per
second. It can be converted to radians per second by dividing by $range.
$v_r is the velocity in recession in kilometers per second. Negative
values indicate that the objects are approaching.
The velocities are only returned if they are available from the input
objects.
B<Caveat:> these velocities are believed by me B<not> to be correct.
=cut
sub equatorial_unreduced {
# Unpack the two objects.
my ($self, $body) = @_;
# Compute the relative positions if there are in fact two objects;
# otherwise just get the absolute position.
my @pos;
if ($body) {
@pos = $body->eci();
my @base = $self->eci($body->universal());
my $limit = @pos < @base ? @pos : @base;
foreach my $inx (0 .. $limit - 1) {
$pos[$inx] -= $base[$inx];
}
splice @pos, $limit;
} else {
$self->{_ECI_cache}{inertial}{equatorial_unreduced} and return
@{$self->{_ECI_cache}{inertial}{equatorial_unreduced}};
@pos = $self->eci();
}
# Rotate the coordinate system so that the second body lies in the
# X-Z plane. The matrix is
# +- -+
# | cos rasc -sin rasc 0 |
# | sin rasc cos rasc 0 |
# | 0 0 1 |
# +- -+
# where rasc is -atan2(y,x). This means sin rasc = -y/h and
# cos rasc = x/h where h = sqrt(x*x + y*y). You postmultiply
# the matrix by the vector to get the new vector.
my $hypot_rasc = sqrt($pos[0] * $pos[0] + $pos[1] * $pos[1]);
# Now we need to rotate the coordinates in the new X-Z plane until
# the second body lies along the X axis. The matrix is
# +- -+
# | cos decl 0 -sin decl |
# | 0 1 0 |
# | sin decl 0 cos decl |
# +- -+
# where decl is -atan2(z,x') (in the new coordinate system), or
# -atan2(z,h) in the old coordinate system. This means that sin decl
# = z/h' and cos decl = x/h' where h' = sqrt(h*h + z*z). Again you
# postmultiply the matrix by the vector to get the result.
my $hypot_decl = sqrt($hypot_rasc * $hypot_rasc + $pos[2] * $pos[2]);
# But at this point we have the equatorial coordinates themselves in
# terms of the original coordinates and the various intermediate
# quantities needed to compute the matrix. So we only need the
# matrix if we need to deal with the velocities. For this, the
# velocity rotation matrix is
# +- -+ +- -+
# | cos decl 0 -sin decl | | cos rasc -sin rasc 0 |
# | 0 1 0 | x | sin rasc cos rasc 0 | =
# | sin decl 0 cos decl | | 0 0 1 |
# +- -+ +- -+
#
# +- -+
# | cos decl cos rasc -cos decl sin rasc -sin decl |
# | sin rasc cos rasc 0 |
# | sin decl cos rasc -sin decl sin rasc cos decl |
# +- +
my @rslt = ([$pos[1], $pos[0]], [$pos[2], $hypot_rasc], $hypot_decl);
if (@pos >= 6) {
my $cos_rasc = $pos[0]/$hypot_rasc;
my $sin_rasc = -$pos[1]/$hypot_rasc;
my $cos_decl = $hypot_rasc/$hypot_decl;
my $sin_decl = -$pos[2]/$hypot_decl;
push @rslt,
$sin_rasc * $pos[3] + $cos_rasc * $pos[4],
$sin_decl * $cos_rasc * $pos[3]
- $sin_decl * $sin_rasc * $pos[4] + $cos_decl * $pos[5],
# Computationally the following is the top row of the
# matrix, but it is swapped to the bottom in the output for
# consistency of returned data sequence.
$cos_decl * $cos_rasc * $pos[3]
- $cos_decl * $sin_rasc * $pos[4] - $sin_decl * $pos[5];
}
$body
or $self->{_ECI_cache}{inertial}{equatorial_unreduced} = \@rslt;
return @rslt;
}
sub _equatorial_reduced {
my ( $self, $body ) = @_;
my @rslt = $self->equatorial_unreduced( $body );
$rslt[0] = mod2pi( atan2( $rslt[0][0], $rslt[0][1] ) );
$rslt[1] = atan2( $rslt[1][0], $rslt[1][1] );
if (@rslt >= 6) {
$rslt[3] /= $rslt[2];
$rslt[4] /= $rslt[2];
}
return @rslt;
}
=item $coord = $coord->equinox_dynamical ($value);
This method sets the value of the equinox_dynamical attribute, and
returns the modified object. If called without an argument, it returns
the current value of the equinox_dynamical attribute.
Yes, it is equivalent to $coord->set (equinox_dynamical => $value) and
$coord->get ('equinox_dynamical'). But there seems to be a significant
performance penalty in the $self->SUPER::set () needed to get this
attribute set from a subclass. It is possible that more methods like
this will be added, but I do not plan to eliminate the set() interface.
=cut
sub equinox_dynamical {
if (@_ > 1) {
$_[0]{equinox_dynamical} = $_[1];
return $_[0];
} else {
return $_[0]{equinox_dynamical};
}
}
=item $coord = $coord->geocentric($psiprime, $lambda, $rho);
This method sets the L</Geocentric> coordinates represented by the
object in terms of L</Geocentric latitude> psiprime and L</Longitude>
lambda in radians, and distance from the center of the Earth rho in
kilometers.
The latitude should be a number between -PI/2 and PI/2 radians
inclusive. The longitude should be a number between -2*PI and 2*PI
radians inclusive. The increased range (one would expect -PI to PI) is
because in some astronomical usages latitudes outside the range + or -
180 degrees are employed. A warning will be generated if either of these
range checks fails.
This method can also be called as a class method, in which case it
instantiates the desired object.
B<This method should not be used with map coordinates> because map
latitude is L</Geodetic latitude>, measured in terms of the tangent of
the reference ellipsoid, whereas geocentric coordinates are,
essentially, spherical coordinates.
The algorithm for conversion between geocentric and ECEF is the
author's.
=item ($psiprime, $lambda, $rho) = $coord->geocentric();
This method returns the L</Geocentric latitude>, L</Longitude>, and
distance to the center of the Earth.
=cut
sub geocentric {
my ($self, @args) = @_;
$self = $self->_check_coord (geocentric => \@args);
unless (@args) {
if ( ! $self->{_ECI_cache}{fixed}{geocentric} ) {
## my ($x, $y, $z, $xdot, $ydot, $zdot) = $self->ecef;
my ($x, $y, $z) = $self->ecef;
my $rsq = $x * $x + $y * $y;
my $rho = sqrt ($z * $z + $rsq);
my $lambda = atan2 ($y, $x);
my $psiprime = atan2 ($z, sqrt ($rsq));
$self->get ('debug') and print <<eod;
Debug geocentric () - ecef -> geocentric
inputs:
x = $x
y = $y
z = $z
outputs:
psiprime = $psiprime
lambda = $lambda
rho = $rho
eod
$self->{_ECI_cache}{fixed}{geocentric} =
[ $psiprime, $lambda, $rho ];
}
return @{ $self->{_ECI_cache}{fixed}{geocentric} };
}
if (@args == 3) {
my ($psiprime, $lambda, $rho) = @args;
$psiprime = _check_latitude(latitude => $psiprime);
$lambda = _check_longitude(longitude => $lambda);
my $z = $rho * sin ($psiprime);
my $r = $rho * cos ($psiprime);
my $x = $r * cos ($lambda);
my $y = $r * sin ($lambda);
$self->get ('debug') and print <<eod;
Debug geocentric () - geocentric -> ecef
inputs:
psiprime = $psiprime
lambda = $lambda
rho = $rho
outputs:
x = $x
y = $y
z = $z
eod
$self->ecef ($x, $y, $z);
$self->{_ECI_cache}{fixed}{geocentric} = \@args;
$self->{specified} = 'geocentric';
$self->{inertial} = 0;
} else {
croak <<eod;
Error - Method geocentric() must be called with either zero arguments
(to retrieve coordinates) or three arguments (to set
coordinates). There is currently no six argument version.
eod
}
return $self;
}
=item $coord = $coord->geodetic($psi, $lambda, $h, $ellipsoid);
This method sets the L</Geodetic> coordinates represented by the object
in terms of its L</Geodetic latitude> psi and L</Longitude> lambda in
radians, and its height h above mean sea level in kilometers.
The latitude should be a number between -PI/2 and PI/2 radians
inclusive. The longitude should be a number between -2*PI and 2*PI
radians inclusive. The increased range (one would expect -PI to PI) is
because in some astronomical usages latitudes outside the range + or -
180 degrees are employed. A warning will be generated if either of these
range checks fails.
The ellipsoid argument is the name of a L</Reference Ellipsoid> known
to the class, and is optional. If passed, it will set the ellipsoid
to be used for calculations with this object. If not passed, the
default ellipsoid is used.
This method can also be called as a class method, in which case it
instantiates the desired object.
The conversion from geodetic to geocentric comes from Jean Meeus'
"Astronomical Algorithms", 2nd Edition, Chapter 11, page 82.
B<This is the method that should be used with map coordinates.>
=item ($psi, $lambda, $h) = $coord->geodetic($ellipsoid);
This method returns the geodetic latitude, longitude, and height
above mean sea level.
The ellipsoid argument is the name of a L</Reference ellipsoid> known
to the class, and is optional. If not specified, the most-recently-set
ellipsoid will be used.
The conversion from geocentric to geodetic comes from Kazimierz
Borkowski's "Accurate Algorithms to Transform Geocentric to Geodetic
Coordinates", at F<http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm>.
This is best viewed with Internet Explorer because of its use of Microsoft's
Symbol font.
=cut
sub geodetic {
my ($self, @args) = @_;
# Detect and acquire the optional ellipsoid name argument. We do
# this before the check, since the check expects the extra
# argument to be a time.
my $elps = (@args == 1 || @args == 4) ? pop @args : undef;
$self = $self->_check_coord (geodetic => \@args, $elps ? $elps : ());
# The following is just a sleazy way to get a consistent
# error message if the ellipsoid name is unknown.
$elps && $self->reference_ellipsoid ($elps);
# If we're fetching the geodetic coordinates
unless (@args) {
# Return cached coordinates if they exist and we did not
# override the default ellipsoid.
return @{$self->{_ECI_cache}{fixed}{geodetic}}
if $self->{_ECI_cache}{fixed}{geodetic} && !$elps;
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug geodetic - explicit ellipsoid ", Dumper ($elps);
};
# Get a reference to the ellipsoid data to use.
$elps = $elps ? $known_ellipsoid{$elps} : $self;
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug geodetic - ellipsoid ", Dumper ($elps);
};
# Calculate geodetic coordinates.
my ($phiprime, $lambda, $rho) = $self->geocentric;
my $r = $rho * cos ($phiprime);
my $b = $elps->{semimajor} * (1- $elps->{flattening});
my $a = $elps->{semimajor};
my $z = $rho * sin ($phiprime);
# The $b is _not_ a magic variable, since we are not inside
# any of the specialized blocks that make $b magic. Perl::Critic
# is simply confused.
$b = - $b ## no critic (RequireLocalizedPunctuationVars)
if $z < 0; # Per Borkowski, for southern hemisphere.
# The following algorithm is due to Kazimierz Borkowski's
# paper "Accurate Algorithms to Transform Geocentric to Geodetic
# Coordinates", from
# http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm
my $bz = $b * $z;
my $asq_bsq = $a * $a - $b * $b;
my $ar = $a * $r;
my $E = ($bz - $asq_bsq) / $ar; # Borkowski (10)
my $F = ($bz + $asq_bsq) / $ar; # Borkowski (11)
my $Q = ($E * $E - $F * $F) * 2; # Borkowski (17)
my $P = ($E * $F + 1) * 4 / 3; # Borkowski (16)
my $D = $P * $P * $P + $Q * $Q; # Borkowski (15)
my $v = $D >= 0 ? do {
my $d = sqrt $D;
my $onethird = 1 / 3;
my $vp = ($d - $Q) ** $onethird - # Borkowski (14a)
($d + $Q) ** $onethird;
$vp * $vp >= abs ($P) ? $vp :
- ($vp * $vp * $vp + 2 * $Q) / # Borkowski (20)
(3 * $P);
} : do {
my $p = - $P;
sqrt (cos (acos ($Q / # Borkowski (14b)
sqrt ($p * $p * $p)) / 3) * $p) * 2;
};
my $G = (sqrt ($E * $E + $v) # Borkowski (13)
+ $E) / 2;
my $t = sqrt ($G * $G + ($F - $v * $G) # Borkowski (12)
/ (2 * $G - $E)) - $G;
my $phi = atan2 (($a * (1 - $t * $t)) / # Borkowski (18)
(2 * $b * $t), 1); # equivalent to atan (arg1)
my $h = ($r - $a * $t) * cos ($phi) + # Borkowski (19)
($z - $b) * sin ($phi);
# End of Borkowski's algorthm.
# Cache the results of the calculation if they were done using
# the default ellipsoid.
$self->{_ECI_cache}{fixed}{geodetic} = [$phi, $lambda, $h] unless $elps;
# Return the results in any event.
$self->get ('debug') and print <<eod;
Debug geodetic: geocentric to geodetic
inputs:
phiprime = $phiprime
lambda = $lambda
rho = $rho
intermediates:
z = $z
r = $r
E = $E
F = $F
P = $P
Q = $Q
D = $D
v = $v
G = $G
t = $t
outputs:
phi = atan2 (a * (1 - t * t), 2 * b * t)
= atan2 (@{[$a * (1 - $t * $t)]}, @{[2 * $b * $t]})
= $phi (radians)
h = (r - a * t) * cos (phi) + (z - b) * sin (phi)
= @{[$r - $a * $t]} * cos (phi) + @{[$z - $b]} * sin (phi)
= $h (kilometers)
eod
return ($phi, $lambda, $h);
}
# If we're setting the geodetic coordinates.
if (@args == 3) {
# Set the ellipsoid for the object if one was specified.
$self->set (ellipsoid => $elps) if $elps;
# Calculate the geocentric data.
my ($phi, $lambda, $h) = @args;
$phi = _check_latitude(latitude => $phi);
$lambda = _check_longitude(longitude => $lambda);
my $bovera = 1 - $self->{flattening};
# The following algorithm appears on page 82 of the second
# edition of Jean Meeus' "Astronomical Algorithms."
my $u = atan2 ($bovera * tan ($phi), 1);
my $rhosinlatprime = $bovera * sin ($u) +
$h / $self->{semimajor} * sin ($phi);
my $rhocoslatprime = cos ($u) +
$h / $self->{semimajor} * cos ($phi);
my $phiprime = atan2 ($rhosinlatprime, $rhocoslatprime);
my $rho = $self->{semimajor} * ($rhocoslatprime ?
$rhocoslatprime / cos ($phiprime) :
$rhosinlatprime / sin ($phiprime));
# End of Meeus' algorithm.
# Set the geocentric data as the coordinates.
$self->geocentric ($phiprime, $lambda, $rho);
# Cache the geodetic coordinates.
$self->{_ECI_cache}{fixed}{geodetic} = [$phi, $lambda, $h];
$self->{specified} = 'geodetic';
$self->{inertial} = 0;
# Else if the number of coordinates is bogus, croak.
} else {
croak <<eod;
Error - Method geodetic() must be called with either zero arguments
(to retrieve coordinates) or three arguments (to set
coordinates). There is currently no six argument version.
eod
}
# Return the object, wherever it came from.
return $self;
}
=item $value = $coord->get ($attrib);
This method returns the named attributes of the object. If called in
list context, you can give more than one attribute name, and it will
return all their values.
If called as a class method, it returns the current default values.
See L</Attributes> for a list of the attributes you can get.
=cut
{ # Begin local symbol block.
my %accessor = (
sun => sub {
## my ( $self, $name ) = @_;
my ( $self ) = @_; # Name unused
return ( $self->{sun} ||= $self->SUN_CLASS()->new() );
},
);
sub get {
my ($self, @args) = @_;
ref $self or $self = \%static;
my @rslt;
foreach my $name (@args) {
exists $mutator{$name} or croak <<eod;
Error - Attribute '$name' does not exist.
eod
if ($accessor{$name}) {
push @rslt, $accessor{$name}->($self, $name);
} else {
push @rslt, $self->{$name};
}
}
return wantarray ? @rslt : $rslt[0];
}
} # End local symbol block
=item @coords = $coord->enu();
This method reports the coordinates of C<$coord> in the East-North-Up
coordinate system, as seen from the position stored in the C<station>
attribute of C<$coord>.
Velocity conversions to C<enu()> appear to me to be mostly sane. See
L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below, for details.
=cut
sub enu { # East, North, Up
my ( $self ) = @_;
my @vector = $self->neu();
@vector[ 0, 1 ] = @vector[ 1, 0 ]; # Swap North and East
@vector > 3 # If we have velocity,
and @vector[ 3, 4 ] = @vector[ 4, 3 ]; # Ditto
return @vector;
}
=item @coords = $coord->neu();
This method reports the coordinates of C<$coord2> in the North-East-Up
coordinate system, as seen from the position stored in the C<station>
attribute of C<$coord>. This is a B<left-handed> coordinate system.
Velocity conversions to C<neu()> appear to me to be mostly sane. See
L<A NOTE ON VELOCITIES|/A NOTE ON VELOCITIES>, below, for details.
=cut
sub neu { # North, East, Up
my ( $self ) = @_;
my $station = $self->get( 'station' )
or croak 'Station attribute required';
return $station->_local_cartesian( $self );
}
=item @coords = $coord->heliocentric_ecliptic_cartesian()
This method returns the Heliocentric ecliptic Cartesian position of the
object. You can optionally pass time as an argument; this is equivalent
to
@coords = $coord->universal( $time )
->heliocentric_ecliptic_cartesian();
At this time this object is not able to derive Heliocentric ecliptic
Cartesian coordinates from other coordinates (say, ECI).
=item $coord->heliocentric_ecliptic_cartesian( $X, $Y, $Z )
This method sets the object's position in Heliocentric ecliptic
Cartesian coordinates. Velocities may not be set at the moment, though
this is planned. You may also pass an optional time as the last
argument.
The invocant is returned.
=cut
sub heliocentric_ecliptic_cartesian {
my ( $self, @args ) = @_;
$self = $self->_check_coord( heliocentric_ecliptic_cartesian => \@args );
if ( @args == 3 ) {
$self->{_ECI_cache} = {
inertial => {
heliocentric_ecliptic_cartesian => \@args,
},
};
$self->{specified} = 'heliocentric_ecliptic_cartesian';
$self->{inertial} = 1;
$self->_convert_heliocentric_ecliptic_cartesian_to_eci();
} elsif ( @args ) {
croak 'heliocentric_ecliptic_cartesian() wants 0, 1, 3 or 4 arguments';
} else {
$self->{_ECI_cache}{inertial}{heliocentric_ecliptic_cartesian}
and return @{
$self->{_ECI_cache}{inertial}{heliocentric_ecliptic_cartesian}
};
return $self->_convert_eci_to_heliocentric_ecliptic_cartesian();
}
return $self;
}
=item @coords = $coord->heliocentric_ecliptic();
This method returns the Heliocentric ecliptic coordinates of the object.
=cut
sub heliocentric_ecliptic {
my ( $self, @args ) = @_;
$self = $self->_check_coord( heliocentric_ecliptic => \@args );
if ( @args ) {
@args % 3
and croak 'Arguments must be in threes, plus an optional time';
$self->{_ECI_cache}{inertial}{heliocentric_ecliptic} = \@args;
$self->heliocentric_ecliptic_cartesian(
_convert_spherical_x_to_cartesian( @args ) );
$self->{specified} = 'heliocentric_ecliptic';
return $self;
} else {
return @{
$self->{_ECI_cache}{inertial}{heliocentric_ecliptic} ||= [
_convert_cartesian_to_spherical_x(
$self->heliocentric_ecliptic_cartesian() ) ] }
}
return $self;
}
# ( $temp, $X, $Y, $Z, $Xprime, $Yprime, $Zprime ) =
# $self->_local_cartesian( $trn2 );
# This method computes local Cartesian coordinates of $trn2 as seen from
# $self. The first return is intermediate results useful for the azimuth
# and elevation. The subsequent results are X, Y, and Z coordinates (and
# velocities if available), in the North, East, Up coordinate system.
# This is a left-handed coordinate system, but so is the azel() system,
# which it serves.
sub _local_cartesian {
my ( $self, $trn2 ) = @_;
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug local_cartesian - ", Dumper( $self, $trn2 );
};
my $time = $trn2->universal();
$self->universal( $time );
# Pick up the ECEF of the body and the station
my @tgt = $trn2->ecef();
my @obs = $self->ecef();
# Translate the ECEF coordinates to the station.
my $limit = int( min( scalar @tgt, scalar @obs ) / 3 ) * 3 - 1;
foreach my $inx ( 0 .. $limit ) {
$tgt[$inx] -= $obs[$inx];
}
splice @tgt, $limit + 1;
# Pick up the latitude and longitude.
my ( $lat, $lon ) = $self->geodetic();
my $sinlat = sin $lat;
my $coslat = cos $lat;
my $sinlon = sin $lon;
my $coslon = cos $lon;
# Rotate X->Y to longitude of station,
# then Z->X to latitude of station
@tgt[ 0, 1 ] = (
$tgt[0] * $coslon + $tgt[1] * $sinlon,
- $tgt[0] * $sinlon + $tgt[1] * $coslon,
);
@tgt[ 0, 2 ] = (
$tgt[0] * $sinlat - $tgt[2] * $coslat,
$tgt[0] * $coslat + $tgt[2] * $sinlat,
);
$tgt[0] = - $tgt[0]; # Convert Southing to Northing
if ( @tgt > 5 ) {
@tgt[ 3, 4 ] = (
$tgt[3] * $coslon + $tgt[4] * $sinlon,
- $tgt[3] * $sinlon + $tgt[4] * $coslon,
);
@tgt[ 3, 5 ] = (
$tgt[3] * $sinlat - $tgt[5] * $coslat,
$tgt[3] * $coslat + $tgt[5] * $sinlat,
);
$tgt[3] = - $tgt[3]; # Convert South velocity to North
}
return @tgt;
}
=item $coord = $coord->local_mean_time ($time);
This method sets the local mean time of the object. B<This is not
local standard time,> but the universal time plus the longitude
of the object expressed in seconds. Another definition is mean
solar time plus 12 hours (since the solar day begins at noon).
You will get an exception of some sort if the position of the
object has not been set, or if the object represents inertial
coordinates, or on any subclass whose position depends on the time.
=item $time = $coord->local_mean_time ()
This method returns the local mean time of the object. It will raise
an exception if the time has not been set.
Note that this returns the actual local time, not the GMT equivalent.
This means that in formatting for output, you call
strftime $format, gmtime $coord->local_mean_time ();
=cut
sub local_mean_time {
my ($self, @args) = @_;
ref $self or croak <<eod;
Error - The local_mean_time() method may not be called as a class method.
eod
unless (@args) {
$self->{universal} || croak <<eod;
Error - Object's time has not been set.
eod
$self->{local_mean_time} = $self->universal () +
_local_mean_delta ($self)
unless defined $self->{local_mean_time};
return $self->{local_mean_time};
}
if (@args == 1) {
my $time = shift @args;
$self->{specified} or croak <<eod;
Error - Object's coordinates have not been set.
eod
$self->{inertial} and croak <<eod;
Error - You can not set the local time of an object that represents
a position in an inertial coordinate system, because this
causes the earth-fixed position to change, invalidating the
local time.
$self->can ('time_set') and croak <<eod;
Error - You can not set the local time on an @{[ref $self]}
object, because when you do the time_set() method will just
move the object, invalidating the local time.
eod
$self->universal($time - _local_mean_delta($self));
$self->{local_mean_time} = $time;
} else {
croak <<eod;
Error - The local_mean_time() method must be called with either zero
arguments (to retrieve the time) or one argument (to set
the time).
eod
}
return $self;
}
=item $time = $coord->local_time ();
This method returns the local time (defined as solar time plus 12 hours)
of the given object. There is no way to set the local time.
This is really just a convenience routine - the calculation is simply
the local mean time plus the equation of time at the given time.
Note that this returns the actual local time, not the GMT equivalent.
This means that in formatting for output, you call
strftime $format, gmtime $coord->local_time ();
=cut
sub local_time {
my $self = shift;
my $dyntim = $self->dynamical();
return $self->local_mean_time() + equation_of_time($dyntim);
}
=item ( $maidenhead_loc, $height ) = $coord->maidenhead( $precision );
This method returns the location of the object in the Maidenhead Locator
System. Height above the reference ellipsoid is not part of the system,
but is returned anyway, in kilometers. The $precision is optional, and
is an integer greater than 0.
The default precision is 4, but this can be modified by setting
C<$Astro::Coord::ECI::DEFAULT_MAIDENHEAD_PRECISION> to the desired
value. For example, if you want the default precision to be 3 (which it
probably should have been in the first place), you can use
no warnings qw{ once };
local $Astro::Coord::ECI::DEFAULT_MAIDENHEAD_PRECISION = 3;
Note that precisions greater than 4 are not defined by the standard.
This method extends the system by alternating letters (base 24) with
digits (base 10), but this is unsupported since the results will change,
possibly without notice, if the standard is extended in a manner
incompatible with this implementation.
Conversion of latitudes and longitudes to Maidenhead Grid is subject to
truncation error, perhaps more so since latitude and longitude are
specified in radians. An attempt has been made to minimize this by using
Perl's stringification of numbers to ensure that something that looks
like C<42> is not handled as C<41.999999999385>. This probably amounts
to shifting some grid squares very slightly to the north-west, but in
practice it seems to give better results for points on the boundaries of
the grid squares.
=item $coord->maidenhead( $maidenhead_loc, $height );
This method sets the geodetic location in the Maidenhead Locator System.
Height above the reference ellipsoid is not part of the system, but is
accepted anyway, in kilometers, defaulting to 0.
The actual location set is the center of the specified grid square.
Locations longer than 8 characters are not defined by the standard. This
method extends precision by assuming alternate letters (base 24) and
digits (base 10), but this will change, possibly without notice, if the
standard is extended in a manner incompatible with this implementation.
The object itself is returned, to allow call chaining.
=cut
{
our $DEFAULT_MAIDENHEAD_PRECISION = 4;
my $alpha = 'abcdefghijklmnopqrstuvwxyz';
sub maidenhead {
my ( $self, @args ) = @_;
if ( @args > 1 || defined $args[0] && $args[0] =~ m/ \D /smx ) {
my ( $loc, $alt ) = @args;
defined $alt or $alt = 0;
my $precision = length $loc;
$precision % 2
and croak "Invalid Maidenhead locator; length must be even";
$precision /= 2;
my $lat = 0.5;
my $lon = 0.5;
my @chars = split qr{}smx, $loc;
foreach my $base ( reverse _maidenhead_bases( $precision ) ) {
foreach ( $lat, $lon ) {
my $chr = pop @chars;
if ( $base > 10 ) {
( my $inx = index $alpha, lc $chr ) < 0
and croak 'Invalid Maidenhead locator; ',
"'$chr' is not a letter";
my $limit = @chars > 1 ? $base - 1 : $base;
$inx > $limit
and croak 'Invalid Maidenhead locator; ',
"'$chr' is greater than '",
substr( $alpha, $limit, 1 ), "'";
$_ += $inx;
$_ /= $base;
} else {
$chr =~ m/ \D /
and croak 'Invalid Maidenhead locator; ',
"'$chr' is not a digit";
$_ += $chr;
$_ /= $base;
}
}
}
$self->geodetic(
deg2rad( $lat * 180 - 90 ),
deg2rad( $lon * 360 - 180 ),
$alt
);
return $self;
} else {
my ( $precision ) = @args;
defined $precision
and $precision > 0
or $precision = $DEFAULT_MAIDENHEAD_PRECISION;
my @bases = reverse _maidenhead_bases( $precision );
# In order to minimize truncation errors we multiply
# everything out, and then work right-to-left using the mod
# operator. We stringify to take advantage of Perl's smarts
# about when something is 42 versus 41.999999999583.
my $multiplier = 1;
foreach ( @bases ) {
$multiplier *= $_;
}
my ( $lat, $lon, $alt ) = $self->geodetic();
$lat = ( $lat + PIOVER2 ) * $multiplier / PI;
$lon = ( $lon + PI ) * $multiplier / TWOPI;
$lat = floor( "$lat" );
$lon = floor( "$lon" );
my @rslt;
foreach my $base ( @bases ) {
foreach ( $lat, $lon ) {
my $inx = $_ % $base;
push @rslt, $base > 10 ?
substr( $alpha, $inx, 1 ) : $inx;
$_ = floor( $_ / $base );
}
}
@rslt
and $rslt[-1] = uc $rslt[-1];
@rslt > 1
and $rslt[-2] = uc $rslt[-2];
return (
join( '', reverse @rslt ),
$alt,
);
}
}
sub _maidenhead_bases {
my ( $precision ) = @_;
my @bases;
foreach my $inx ( 1 .. $precision ) {
push @bases, $inx % 2 ? 24 : 10;
}
$bases[0] = 18;
return @bases;
}
}
=item $value = $coord->mean_angular_velocity();
This method returns the mean angular velocity of the body in radians
per second. If the $coord object has a period() method, this method
just returns two pi divided by the period. Otherwise it returns the
contents of the angularvelocity attribute.
=cut
sub mean_angular_velocity {
my $self = shift;
return $self->can ('period') ?
TWOPI / $self->period :
$self->{angularvelocity};
}
=item $time = $coord->next_azimuth( $body, $azimuth );
This method returns the next time the given C<$body> passes the given
C<$azimuth> as seen from the given C<$coord>, calculated to the nearest
second. The start time is the current time setting of the C<$body>
object.
=item $time = $coord->next_azimuth( $azimuth );
This method returns the next time the C<$coord> object passes the given
C<$azimuth> as seen from the location in the C<$coord> object's
C<station> attribute, calculated to the nearest second. The start time
is the current time setting of the C<$coord> object.
=cut
sub next_azimuth {
my ( $self, $body, $azimuth ) = _expand_args_default_station( @_ );
ref $self or croak <<'EOD';
Error - The next_azimuth() method may not be called as a class method.
EOD
$body->represents(__PACKAGE__) or croak <<"EOD";
Error - The argument to next_azimuth() must be a subclass of
@{[__PACKAGE__]}.
EOD
my $want = shift;
defined $want and $want = $want ? 1 : 0;
my $denom = $body->mean_angular_velocity -
$self->mean_angular_velocity;
## my $retro = $denom >= 0 ? 0 : 1;
($denom = abs ($denom)) < 1e-11 and croak <<eod;
Error - The next_azimuth() method will not work for geosynchronous
bodies.
eod
my $apparent = TWOPI / $denom;
my $begin = $self->universal;
my $delta = floor( $apparent / 8 );
my $end = $begin + $delta;
my $begin_angle = mod2pi(
( $self->azel( $body->universal( $begin ) ) )[0] - $azimuth );
my $end_angle = mod2pi(
( $self->azel( $body->universal( $end ) ) )[0] - $azimuth );
while ( $begin_angle < PI || $end_angle >= PI ) {
$begin_angle = $end_angle;
$begin = $end;
$end = $end + $delta;
$end_angle = mod2pi(
( $self->azel( $body->universal( $end ) ) )[0] - $azimuth );
}
while ($end - $begin > 1) {
my $mid = floor (($begin + $end) / 2);
my $mid_angle = mod2pi(
( $self->azel( $body->universal( $mid ) ) )[0] - $azimuth );
( $begin, $end ) = ( $mid_angle >= PI ) ?
( $mid, $end ) : ( $begin, $mid );
}
$body->universal ($end);
$self->universal ($end);
return $end;
}
=item ($time, $rise) = $coord->next_elevation ($body, $elev, $upper)
This method calculates the next time the given body passes above or
below the given elevation (in radians) as seen from C<$coord>. The
C<$elev> argument may be omitted (or passed as undef), and will default
to 0. If the C<$upper> argument is true, the calculation will be based
on the upper limb of the body (as determined from its C<angulardiameter>
attribute); if false, the calculation will be based on the center of the
body. The C<$upper> argument defaults to true if the C<$elev> argument
is zero or positive, and false if the C<$elev> argument is negative.
The algorithm is successive approximation, and assumes that the
body will be at its highest at meridian passage. It also assumes
that if the body hasn't passed the given elevation in 183 days it
never will. In this case it returns undef in scalar context, or
an empty list in list context.
=item ($time, $rise) = $coord->next_elevation( $elev, $upper );
This method calculates the next time the C<$coord> object passes above
or below the given elevation (in radians) as seen from the position
found in the C<$coord> object's C<station> attribute. The C<$elev>
argument may be omitted (or passed as undef), and will default to 0. If
the C<$upper> argument is true, the calculation will be based on the
upper limb of the body (as determined from its C<angulardiameter>
attribute); if false, the calculation will be based on the center of the
body. The C<$upper> argument defaults to true if the C<$elev> argument
is zero or positive, and false if the C<$elev> argument is negative.
The algorithm is successive approximation, and assumes that the
body will be at its highest at meridian passage. It also assumes
that if the body hasn't passed the given elevation in 183 days it
never will. In this case it returns undef in scalar context, or
an empty list in list context.
=cut
use constant NEVER_PASS_ELEV => 183 * SECSPERDAY;
sub next_elevation {
my ( $self, $body, $angle, $upper ) = _expand_args_default_station( @_ );
ref $self or croak <<'EOD';
Error - The next_elevation() method may not be called as a class method.
EOD
$body->represents(__PACKAGE__) or croak <<"EOD";
Error - The first argument to next_elevation() must be a subclass of
@{[__PACKAGE__]}.
EOD
$angle ||= 0;
defined $upper or $upper = $angle >= 0;
$upper = $upper ? 1 : 0;
my $begin = $self->universal;
my $original = $begin;
my ( undef, $elev ) = $self->azel_offset(
$body->universal( $begin ), $upper );
my $rise = ( $elev < $angle ) || 0;
my ($end, $above) = $self->next_meridian ($body, $rise);
my $give_up = $body->NEVER_PASS_ELEV ();
while ( 1 ) {
my ( undef, $elev ) = $self->azel_offset( $body, $upper );
( ( $elev < $angle ) || 0 ) == $rise
or last;
return if $end - $original > $give_up;
$begin = $end;
($end, $above) = $self->next_meridian ($body, $rise);
}
while ($end - $begin > 1) {
my $mid = floor (($begin + $end) / 2);
my ( undef, $elev ) = $self->universal( $mid )->
azel_offset( $body->universal( $mid ), $upper );
( $begin, $end ) =
($elev < $angle || 0) == $rise ? ($mid, $end) : ($begin, $mid);
}
$self->universal ($end); # Ensure consistent time.
$body->universal ($end);
return wantarray ? ($end, $rise) : $end;
}
=item ( $time, $above ) = $coord->next_meridian( $body, $want )
This method calculates the next meridian passage of the given C<$body>
over (or under) the location specified by the C<$coord> object. The
C<$body> object must be a subclass of Astro::Coord::ECI.
The optional C<$want> argument should be specified as true (e.g. 1) if
you want the next passage above the observer, or as false (e.g. 0) if
you want the next passage below the observer. If this argument is
omitted or undefined, you get whichever passage is next.
The start time of the search is the current time setting of the
C<$coord> object.
The returns are the time of the meridian passage, and an indicator which
is true if the passage is above the observer (i.e. local noon if the
C<$body> represents the Sun), or false if below (i.e. local midnight if
the C<$body> represents the Sun). If called in scalar context, you get
the time only.
The current time of both C<$coord> and C<$body> objects are left at the
returned time.
The algorithm is by successive approximation. It will croak if the
period of the C<$body> is close to synchronous, and will probably not
work well for bodies in highly eccentric orbits. The calculation is to
the nearest second, and the time returned is the first even second after
the body crosses the meridian.
=item ( $time, $above ) = $coord->next_meridian( $want )
This method calculates the next meridian passage of the C<$coord> object
over (or under) the location specified by the C<$coord> object's
C<station> attribute.
The optional C<$want> argument should be specified as true (e.g. 1) if
you want the next passage above the observer, or as false (e.g. 0) if
you want the next passage below the observer. If this argument is
omitted or undefined, you get whichever passage is next.
The start time of the search is the current time setting of the
C<$coord> object.
The returns are the time of the meridian passage, and an indicator which
is true if the passage is above the observer (i.e. local noon if the
C<$coord> object represents the Sun), or false if below (i.e. local
midnight if the C<$coord> object represents the Sun). If called in
scalar context, you get the time only.
The current time of both C<$coord> and its C<station> are left at the
returned time.
The algorithm is by successive approximation. It will croak if the
period of the C<$body> is close to synchronous, and will probably not
work well for bodies in highly eccentric orbits. The calculation is to
the nearest second, and the time returned is the first even second after
the body crosses the meridian.
=cut
sub next_meridian {
my ( $self, $body, $want ) = _expand_args_default_station( @_ );
ref $self or croak <<'EOD';
Error - The next_meridian() method may not be called as a class method.
EOD
$body->represents(__PACKAGE__) or croak <<"EOD";
Error - The argument to next_meridian() must be a subclass of
@{[__PACKAGE__]}.
EOD
defined $want and $want = $want ? 1 : 0;
my $denom = $body->mean_angular_velocity -
$self->mean_angular_velocity;
my $retro = $denom >= 0 ? 0 : 1;
($denom = abs ($denom)) < 1e-11 and croak <<'EOD';
Error - The next_meridian() method will not work for geosynchronous
bodies.
EOD
my $apparent = TWOPI / $denom;
my $begin = $self->universal;
my $delta = floor ($apparent / 16);
my $end = $begin + $delta;
my ($above, $opposite) =
mod2pi (($body->universal($begin)->geocentric)[1]
- ($self->universal($begin)->geocentric)[1]) >= PI ?
(1 - $retro, PI) : ($retro, 0);
($begin, $end) = ($end, $end + $delta)
while mod2pi (($body->universal($end)->geocentric)[1] -
($self->universal($end)->geocentric)[1] + $opposite) < PI;
if (defined $want && $want != $above) {
$above = $want;
$opposite = $opposite ? 0 : PI;
($begin, $end) = ($end, $end + $delta)
while mod2pi (($body->universal($end)->geocentric)[1] -
($self->universal($end)->geocentric)[1] + $opposite) < PI;
}
while ($end - $begin > 1) {
my $mid = floor (($begin + $end) / 2);
my $long = ($body->universal($mid)->geocentric)[1];
my $merid = ($self->universal($mid)->geocentric)[1];
($begin, $end) =
mod2pi ($long - $merid + $opposite) < PI ?
($mid, $end) : ($begin, $mid);
}
$body->universal ($end);
$self->universal ($end);
return wantarray ? ($end, $above) : $end;
}
=item $coord = $coord->precess ($time);
This method is a convenience wrapper for precess_dynamical(). The
functionality is the same except that B<the time is specified in
universal time.>
=cut
sub precess {
$_[1]
and $_[1] += dynamical_delta( $_[1] );
goto &precess_dynamical;
}
=item $coord = $coord->precess_dynamical ($time);
This method precesses the coordinates of the object to the given
equinox, B<specified in dynamical time.> The starting equinox is the
value of the 'equinox_dynamical' attribute, or the current time setting
if the 'equinox_dynamical' attribute is any false value (i.e. undef, 0,
or ''). A warning will be issued if the value of 'equinox_dynamical' is
undef (which is the default setting) and the object represents inertial
coordinates. As of version 0.013_02, B<the time setting of the object is
unaffected by this operation.>
As a side effect, the value of the 'equinox_dynamical' attribute will be
set to the dynamical time corresponding to the argument.
As of version 0.061_01, this does nothing to non-inertial
objects -- that is, those whose position was set in Earth-fixed
coordinates.
If the object's 'station' attribute is set, the station is also
precessed.
The object itself is returned.
The algorithm comes from Jean Meeus, "Astronomical Algorithms", 2nd
Edition, Chapter 21, pages 134ff (a.k.a. "the rigorous method").
=cut
sub precess_dynamical {
my ( $self, $end ) = @_;
$end
or croak "No equinox time specified";
# Non-inertial coordinate systems are not referred to the equinox,
# and so do not get precessed.
$self->get( 'inertial' )
or return $self;
defined ( my $start = $self->get( 'equinox_dynamical' ) )
or carp "Warning - Precess called with equinox_dynamical ",
"attribute undefined";
$start ||= $self->dynamical ();
my $sta;
if ( $sta = $self->get( 'station' ) and $sta->get( 'inertial' )
) {
$sta->get( 'station' )
and croak NO_CASCADING_STATIONS;
$sta->universal( $self->universal() );
$sta->precess_dynamical( $end );
}
$start == $end
and return $self;
my ($alpha0, $delta0, $rho0) = $self->equatorial ();
my $T = jcent2000($start);
my $t = ($end - $start) / (36525 * SECSPERDAY);
# The following four assignments correspond to Meeus' (21.2).
my $zterm = (- 0.000139 * $T + 1.39656) * $T + 2306.2181;
my $zeta = deg2rad((((0.017998 * $t + (- 0.000344 * $T + 0.30188)) *
$t + $zterm) * $t) / 3600);
my $z = deg2rad((((0.018203 * $t + (0.000066 * $T + 1.09468)) * $t +
$zterm) * $t) / 3600);
my $theta = deg2rad(((( - 0.041833 * $t - (0.000217 * $T + 0.42665))
* $t + (- 0.000217 * $T - 0.85330) * $T + 2004.3109) * $t) /
3600);
# The following assignments correspond to Meeus' (21.4).
my $sindelta0 = sin ($delta0);
my $cosdelta0 = cos ($delta0);
my $sintheta = sin ($theta);
my $costheta = cos ($theta);
my $cosdelta0cosalpha0 = cos ($alpha0 + $zeta) * $cosdelta0;
my $A = $cosdelta0 * sin ($alpha0 + $zeta);
my $B = $costheta * $cosdelta0cosalpha0 - $sintheta * $sindelta0;
my $C = $sintheta * $cosdelta0cosalpha0 + $costheta * $sindelta0;
my $alpha = mod2pi(atan2 ($A , $B) + $z);
my $delta = asin($C);
$self->equatorial ($alpha, $delta, $rho0);
$self->set (equinox_dynamical => $end);
return $self;
}
=item Astro::Coord::ECI->reference_ellipsoid($semi, $flat, $name);
This class method can be used to define or redefine reference
ellipsoids.
Nothing bad will happen if you call this as an object method, but
it still just creates a reference ellipsoid definition -- the object
is unaffected.
It is not an error to redefine an existing ellipsoid.
=item ($semi, $flat, $name) = Astro::Coord::ECI->reference_ellipsoid($name)
This class method returns the definition of the named reference
ellipsoid. It croaks if there is no such ellipsoid.
You can also call this as an object method, but the functionality is
the same.
The following reference ellipsoids are known to the class initially:
CLARKE-1866 - 1866.
semimajor => 6378.2064 km, flattening => 1/294.98.
Source: http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
GRS67 - Geodetic Reference System 1967.
semimajor => 6378.160 km, flattening => 1/298.247.
Source: http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
GRS80 - Geodetic Reference System 1980.
semimajor => 6378.137 km, flattening => 1/298.25722210088
(flattening per U.S. Department of Commerce 1989).
Source: http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
IAU68 - International Astronomical Union, 1968.
semimajor => 6378.160 km, flattening => 1/298.25.
Source: http://maic.jmu.edu/sic/standards/ellipsoid.htm
IAU76 - International Astronomical Union, 1976.
semimajor => 6378.14 km, flattening => 1 / 298.257.
Source: Jean Meeus, "Astronomical Algorithms", 2nd Edition
NAD83 - North American Datum, 1983.
semimajor => 6378.137 km, flattening => 1/298.257.
Source: http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
(NAD83 uses the GRS80 ellipsoid)
sphere - Just in case you were wondering how much difference it
makes (a max of 11 minutes 32.73 seconds of arc, per Jean
Meeus).
semimajor => 6378.137 km (from GRS80), flattening => 0.
WGS72 - World Geodetic System 1972.
semimajor => 6378.135 km, flattening=> 1/298.26.
Source: http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
WGS84 - World Geodetic System 1984.
semimajor => 6378.137 km, flattening => 1/1/298.257223563.
Source: http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
Reference ellipsoid names are case-sensitive.
The default model is WGS84.
=cut
# Wish I had:
# Maling, D.H., 1989, Measurements from Maps: Principles and methods of cartometry, Pergamon Press, Oxford, England.
# Maling, D.H., 1993, Coordinate Systems and Map Projections, Pergamon Press, Oxford, England.
# http://www.gsi.go.jp/PCGIAP/95wg/wg3/geodinf.htm has a partial list of who uses
# what in the Asia/Pacific.
%known_ellipsoid = ( # Reference Ellipsoids
'CLARKE-1866' => { # http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
semimajor => 6378.2064,
flattening => 1/294.9786982,
},
GRS67 => { # http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
semimajor => 6378.160,
flattening => 1/298.247167427,
},
GRS80 => { # http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
semimajor => 6378.137, # km
flattening => 1/298.25722210088, # U.S. Dept of Commerce 1989 (else 1/298.26)
},
IAU68 => { # http://maic.jmu.edu/sic/standards/ellipsoid.htm
semimajor => 6378.160,
flattening => 1/298.25,
},
IAU76 => { # Meeus, p. 82.
semimajor => 6378.14,
flattening => 1/298.257,
},
NAD83 => { # http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
semimajor => 6378.137, # km
flattening => 1/298.257,
},
sphere => { # Defined by me for grins, with semimajor from GRS80.
semimajor => 6378.137, # km, from GRS80
flattening => 0, # It's a sphere!
},
WGS72 => { # http://biology.usgs.gov/fgdc.metadata/version2/spref/horiz/geodet/faq.htm
semimajor => 6378.135, # km
flattening => 1/298.26,
},
WGS84 => { # http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
semimajor => 6378.137,
flattening => 1/298.257223563,
},
);
foreach my $name (keys %known_ellipsoid) {
$known_ellipsoid{$name}{name} = $name;
}
sub reference_ellipsoid {
my ( undef, @args ) = @_; # Invocant unused
my $name = pop @args or croak <<eod;
Error - You must specify the name of a reference ellipsoid.
eod
if (@args == 0) {
$known_ellipsoid{$name} or croak <<eod;
Error - Reference ellipsoid $name is unknown to this software. Known
ellipsoids are:
@{[join ', ', sort keys %known_ellipsoid]}.
eod
} elsif (@args == 2) {
$known_ellipsoid{$name} = {
semimajor => $args[0],
flattening => $args[1],
name => $name,
};
} else {
croak <<eod;
Error - You must call the reference_ellipsoid class method with either one
argument (to fetch the definition of a known ellipsoid) or three
arguments (to define a new ellipsoid or redefine an old one).
eod
}
return @{$known_ellipsoid{$name}}{qw{semimajor flattening name}}
}
=item $coord->represents ($class);
This method returns true if the $coord object represents the given class.
It is pretty much like isa (), but if called on a container class (i.e.
Astro::Coord::ECI::TLE::Set), it returns true based on the class of
the members of the set, and dies if the set has no members.
The $class argument is optional. If not specified (or undef), it is
pretty much like ref $coord || $coord (i.e. it returns the class
name), but with the delegation behavior described in the previous
paragraph if the $coord object is a container.
There. This took many more words to explain than it did to implement.
=cut
sub represents {
return defined ($_[1]) ?
## $_[0]->represents()->isa($_[1]) :
__classisa($_[0]->represents(), $_[1]) ? 1 : 0 :
(ref $_[0] || $_[0]);
}
=item $coord->set (name => value ...);
This method sets various attributes of the object. If called as a class
method, it changes the defaults.
For reasons that seemed good at the time, this method returns the
object it was called on (i.e. $coord in the above example).
See L</Attributes> for a list of the attributes you can set.
=cut
use constant SET_ACTION_NONE => 0; # Do nothing.
use constant SET_ACTION_RESET => 1; # Reset the coordinates based on the initial setting.
sub set {
my ($self, @args) = @_;
my $original = $self;
ref $self
or $self = \%static;
@args % 2
and croak 'Error - The set() method requires an even ',
'number of arguments.';
my $action = 0;
while (@args) {
my $name = shift @args;
exists $mutator{$name}
or croak "Error - Attribute '$name' does not exist.";
CODE_REF eq ref $mutator{$name}
or croak "Error - Attribute '$name' is read-only";
$action |= $mutator{$name}->($self, $name, shift @args);
}
$self->{_need_purge} = 1
if ref $self && $self->{specified} && $action & SET_ACTION_RESET;
return $original;
}
# The following are the mutators for the attributes. All are
# passed three arguments: a reference to the hash to be set,
# the hash key to be set, and the value. They must return the
# bitwise 'or' of the desired action masks, defined above.
%mutator = (
almanac_horizon => \&_set_almanac_horizon,
angularvelocity => \&_set_value,
debug => \&_set_value,
diameter => \&_set_value,
edge_of_earths_shadow => \&_set_value,
ellipsoid => \&_set_reference_ellipsoid,
equinox_dynamical => \&_set_value, # CAVEAT: _convert_eci_to_ecef
# accesses this directly for
# speed.
flattening => \&_set_custom_ellipsoid,
frequency => \&_set_value,
horizon => \&_set_value,
id => \&_set_id,
inertial => undef,
name => \&_set_value,
refraction => \&_set_value,
specified => undef,
semimajor => \&_set_custom_ellipsoid,
station => \&_set_station,
sun => \&_set_sun,
twilight => \&_set_value,
);
{
# TODO this will all be nicer if I had state variables.
my %special = (
## horizon => sub { return $_[0]->get( 'horizon' ); },
height => sub { return $_[0]->dip(); },
);
sub _set_almanac_horizon {
my ( $hash, $attr, $value ) = @_;
defined $value
or $value = 0; # Default
if ( $special{$value} ) {
$hash->{"_$attr"} = $special{$value};
} elsif ( looks_like_number( $value ) ) {
$hash->{"_$attr"} = sub { return $_[0]->get( $attr ) };
} else {
croak "'$value' is an invalid value for '$attr'";
}
$hash->{$attr} = $value;
return SET_ACTION_NONE;
}
}
# Get the actual value of the almanac horizon, from whatever
# source.
sub __get_almanac_horizon {
goto $_[0]{_almanac_horizon};
}
# If you set semimajor or flattening, the ellipsoid name becomes
# undefined. Also clear any cached geodetic coordinates.
sub _set_custom_ellipsoid {
$_[0]->{ellipsoid} = undef;
$_[0]->{$_[1]} = $_[2];
return SET_ACTION_RESET;
}
sub _set_station {
my ( $hash, $attr, $value ) = @_;
if ( defined $value ) {
embodies( $value, 'Astro::Coord::ECI' )
or croak "Attribute $attr must be undef or an ",
'Astro::Coord::ECI object';
$value->get( 'station' )
and croak NO_CASCADING_STATIONS;
}
$hash->{$attr} = $value;
return SET_ACTION_RESET;
}
use constant SUN_CLASS => 'Astro::Coord::ECI::Sun';
sub _set_sun {
## my ( $self, $name, $value ) = @_;
my ( $self, undef, $value ) = @_; # Name unused
embodies( $value, $self->SUN_CLASS() )
or croak 'The value of the sun attribute must represent the Sun';
ref $value
or $value = $value->new();
return SET_ACTION_NONE;
}
# Unfortunately, the TLE subclass may need objects reblessed if
# the ID changes. So much for factoring. Sigh.
sub _set_id {
$_[0]{$_[1]} = $_[2];
$_[0]->rebless () if $_[0]->can ('rebless');
return SET_ACTION_NONE;
}
# If this is a reference ellipsoid name, check it, and if it's
# OK set semimajor and flattening also. Also clear any cached
# geodetic coordinates.
sub _set_reference_ellipsoid {
my ($self, $name, $value) = @_;
defined $value or croak <<eod;
Error - Can not set reference ellipsoid to undefined.
eod
exists $known_ellipsoid{$value} or croak <<eod;
Error - Reference ellipsoid '$value' is unknown.
eod
$self->{semimajor} = $known_ellipsoid{$value}{semimajor};
$self->{flattening} = $known_ellipsoid{$value}{flattening};
$self->{$name} = $value;
return SET_ACTION_RESET;
}
# If this is a vanilla setting, just do it.
sub _set_value {
$_[0]->{$_[1]} = $_[2];
return SET_ACTION_NONE;
}
=item $coord->universal ($time)
This method sets the time represented by the object, in universal time
(a.k.a. CUT, a.k.a. Zulu, a.k.a. Greenwich).
This method can also be called as a class method, in which case it
instantiates the desired object.
=item $time = $coord->universal ();
This method returns the universal time previously set.
=cut
sub universal {
my ( $self, $time, @args ) = @_;
@args
and croak <<'EOD';
Error - The universal() method must be called with either zero
arguments (to retrieve the time) or one argument (to set the
time).
EOD
if ( defined $time ) {
return $self->__universal( $time, 1 );
} else {
ref $self or croak <<'EOD';
Error - The universal() method may not be called as a class method
unless you specify arguments.
EOD
defined $self->{universal}
or croak <<'EOD';
Error - Object's time has not been set.
EOD
return $self->{universal};
}
}
# Set universal time without running model
sub __universal {
my ( $self, $time, $run_model ) = @_;
defined $time
or confess 'Programming error - __universal() requires a defined time';
$self = $self->new () unless ref $self;
defined $self->{universal}
and $time == $self->{universal}
and return $self;
$self->__clear_time();
$self->{universal} = $time;
$run_model
and $self->_call_time_set();
return $self;
}
sub __clear_time {
my ( $self ) = @_;
delete $self->{universal};
delete $self->{local_mean_time};
delete $self->{dynamical};
if ($self->{specified}) {
if ($self->{inertial}) {
delete $self->{_ECI_cache}{fixed};
} else {
delete $self->{_ECI_cache}{inertial};
}
}
return;
}
#######################################################################
#
# Internal
#
# $coord->_call_time_set ()
# This method calls the time_set method if it exists and if we are
# not already in it. It is a way to avoid endless recursion if the
# time_set method should happen to set the time.
sub _call_time_set {
my $self = shift;
$self->can ('time_set') or return;
my $exception;
unless ($self->{_no_set}++) {
eval {$self->time_set (); 1;}
or $exception = $@;
}
--$self->{_no_set} or delete $self->{_no_set};
if ($exception) {
$self->__clear_time();
# Re-raise the exception now that we have cleaned up.
die $exception; ## no critic (RequireCarping)
}
return;
}
# $coord->_check_coord (method => \@_)
# This is designed to be called "up front" for any of the methods
# that retrieve or set coordinates, to be sure the object is in
# a consistent state.
# * If $self is not a reference, it creates a new object if there
# are arguments, or croaks if there are none.
# * If the number arguments passed (after removing self and the
# method name) is one more than a multiple of three, the last
# argument is removed, and used to set the universal time of
# the object.
# * If there are arguments, all coordinate caches are cleared;
# otherwise the coordinates are reset if needed.
# The object itself is returned.
sub _check_coord {
my ($self, $method, $args, @extra) = @_;
unless (ref $self) {
@$args or croak <<eod;
Error - The $method() method may not be called as a class method
unless you specify arguments.
eod
$self = $self->new ();
}
$self->{debug} and do {
local $Data::Dumper::Terse = 1;
print "Debug $method (", Dumper (@$args, @extra), ")\n";
};
{
my $inx = 0;
local $Data::Dumper::Terse = 1;
foreach (@$args) {
croak <<eod unless defined $_;
Error - @{[ (caller (1))[3] ]} argument $inx is undefined.
Arguments are (@{[ join ', ',
map {my $x = Dumper $_; chomp $x; $x} @$args ]})
eod
$inx++;
}
}
$self->universal (pop @$args) if @$args % 3 == 1;
if ($self->{specified}) {
if (@$args) {
# Cached coordinate deletion moved to ecef(), so it's only done once.
} elsif ($self->{_need_purge}) {
delete $self->{_need_purge};
my $spec = $self->{specified};
my $data = [$self->$spec ()];
foreach my $key (@kilatr) {delete $self->{$key}}
$self->$spec (@$data);
}
}
return $self;
}
# _check_latitude($arg_name => $value)
#
# This subroutine range-checks the given latitude value,
# generating a warning if it is outside the range -PIOVER2 <=
# $value <= PIOVER2. The $arg_name is used in the exception, if
# any. The value is normalized to the range -PI to PI, and
# returned. Not that a value outside the validation range makes
# any sense.
sub _check_latitude {
(-&PIOVER2 <= $_[1] && $_[1] <= &PIOVER2)
or carp (ucfirst $_[0],
' must be in radians, between -PI/2 and PI/2');
return mod2pi($_[1] + PI) - PI;
}
# $value = _check_longitude($arg_name => $value)
#
# This subroutine range-checks the given longitude value,
# generating a warning if it is outside the range -TWOPI <= $value
# <= TWOPI. The $arg_name is used in the exception, if any. The
# value is normalized to the range -PI to PI, and returned.
sub _check_longitude {
(-&TWOPI <= $_[1] && $_[1] <= &TWOPI)
or carp (ucfirst $_[0],
' must be in radians, between -2*PI and 2*PI');
return mod2pi($_[1] + PI) - PI;
}
# _check_right_ascension($arg_name => $value)
#
# This subroutine range-checks the given right ascension value,
# generating a warning if it is outside the range 0 <= $value <=
# TWOPI. The $arg_name is used in the exception, if any. The value
# is normalized to the range 0 to TWOPI, and returned.
sub _check_right_ascension {
(0 <= $_[1] && $_[1] <= &TWOPI)
or carp (ucfirst $_[0],
' must be in radians, between 0 and 2*PI');
return mod2pi($_[1]);
}
# @cartesian = _convert_spherical_to_cartesian( @spherical )
#
# This subroutine converts spherical coordinates to Cartesian
# coordinates of the same handedness.
#
# The first three arguments are Phi (in the X-Y plane, e.g. azimuth or
# longitude, in radians), Theta (perpendicular to the X-Y plane, e.g.
# elevation or latitude, in radians) and Rho (range). Subsequent
# triplets of arguments are optional, and can represent derivitaves of
# position (velocity, acceleration ... ) at that point.
#
# The return is the corresponding X, Y, and Z coordinates. If more than
# three triplets of arguments are specified, they will be converted to
# spherical coordinates as though measured at that point, and returned
# in the same order. That is, if you supplied Phi, Theta, and Rho
# velocities, you will get back X, Y, and Z velocities. velocities, you
# will get back Phi, Theta, and Rho velocities, in that order.
sub _convert_spherical_to_cartesian {
my @sph_data = @_;
@sph_data
and not @sph_data % 3
or confess( 'Programming error - Want 3 or 6 arguments' );
# The first triplet is position. We extract it into its own array,
# then compute the corresponding Cartesian coordinates using the
# definition of the coordinate system.
my @sph_pos = splice @sph_data, 0, 3;
my $range = $sph_pos[2];
my $sin_theta = sin $sph_pos[1];
my $cos_theta = cos $sph_pos[1];
my $sin_phi = sin $sph_pos[0];
my $cos_phi = cos $sph_pos[0];
my $diag = $range * $cos_theta;
my $X = $diag * $cos_phi;
my $Y = $diag * $sin_phi;
my $Z = $range * $sin_theta;
# Accumulate results.
my @rslt = ( $X, $Y, $Z );
# If we have velocity (accelelation, etc) components
if ( @sph_data ) {
# We compute unit vectors in the Cartesian coordinate system.
# x = cos Theta cos Phi r - sin Theta cos Phi theta - sin Phi phi
# y = cos Theta sin Phi r - sin Theta sin Phi theta - cos Phi phi
# z = sin Theta r + cos Theta theta
my $X_hat = [ - $sin_phi, - $sin_theta * $cos_phi,
$cos_theta * $cos_phi ];
my $Y_hat = [ $cos_phi, - $sin_theta * $sin_phi,
$cos_theta * $sin_phi ];
my $Z_hat = [ 0, $cos_theta, $sin_theta ];
while ( @sph_data ) {
# Each triplet is then converted by projecting the Cartesian
# vector onto the appropriate unit vector. Azimuth and
# elevation are also converted to length by multiplying by
# the range. NOTE that this is the small-angle
# approximation, but should be OK since we assume we're
# dealing with derivatives of position.
my @sph_info = splice @sph_data, 0, 3;
$sph_info[0] *= $range;
$sph_info[1] *= $range;
push @rslt, vector_dot_product( $X_hat, \@sph_info );
push @rslt, vector_dot_product( $Y_hat, \@sph_info );
push @rslt, vector_dot_product( $Z_hat, \@sph_info );
}
}
return @rslt;
}
# @cartesian = _convert_spherical_x_to_cartesian( @spherical )
#
# This subroutine is a convenience wrapper for
# _convert_spherical_to_cartesian(). The only difference is that the
# arguments are Theta (perpendicular to the X-Y plane), Phi (in the X-Y
# plane) and Rho, rather than Phi, Theta, Rho. This subroutine is my
# penance for not having all the methods that involve spherical
# coordinates return their arguments in the same order.
sub _convert_spherical_x_to_cartesian {
my @args = @_;
for ( my $inx = 0; $inx < @args; $inx += 3 ) {
@args[ $inx, $inx + 1 ] = @args[ $inx + 1, $inx ];
}
return _convert_spherical_to_cartesian( @args );
}
# @spherical = _convert_cartesian_to_spherical( @cartesian )
#
# This subroutine converts three-dimensional Cartesian coordinates to
# spherical coordinates of the same handedness.
#
# The first three arguments are the X, Y, and Z coordinates, and are
# required. Subsequent triplets af arguments are optional, and can
# represent anything (velocity, acceleration ... ) at that point.
#
# The return is the spherical coordinates Phi (in the X-Y plane, e.g.
# azimuth or longitude, in radians), Theta (perpendicular to the X-Y
# plane, e.g. elevation or latitude, in radians), and Rho (range). If
# more than three triplets of arguments are specified, they will be
# converted to spherical coordinates as though measured at that point,
# and returned in the same order. That is, if you supplied X, Y, and Z
# velocities, you will get back Phi, Theta, and Rho velocities, in that
# order.
sub _convert_cartesian_to_spherical {
my @cart_data = @_;
@cart_data
and not @cart_data % 3
or confess( 'Programming error - Want 3 or 6 arguments' );
# The first triplet is position. We extract it into its own array,
# then compute the corresponding spherical coordinates using the
# definition of the coordinate system.
my @cart_pos = splice @cart_data, 0, 3;
my $range = vector_magnitude( \@cart_pos );
my $azimuth = ( $cart_pos[0] || $cart_pos[1] ) ?
mod2pi( atan2 $cart_pos[1], $cart_pos[0] ) :
0;
my $elevation = $range ? asin( $cart_pos[2] / $range ) : 0;
# Accumulate results.
my @rslt = ( $azimuth, $elevation, $range );
# If we have velocity (accelelation, etc) components
if ( @cart_data ) {
# We compute unit vectors in the spherical coordinate system.
#
# The "Relationships among Unit Vectors" at
# http://plaza.obu.edu/corneliusk/mp/rauv.pdf (and retained in
# the ref/ directory) gives the transformation both ways. With
# x, y, and z being the Cartesian unit vectors, Theta and Phi
# being the elevation (in the range 0 to pi, 0 being along the +
# Z axis) and azimuth (X toward Y, i.e. right-handed), and r,
# theta, and phi being the corresponding unit vectors:
#
# r = sin Theta cos Phi x + sin Theta sin Phi y + cos Theta z
# theta = cos Theta cos Phi x + cos Theta sin Phi y - sin Theta z
# phi = - sin Phi x + cos phi y
#
# and
#
# x = sin Theta cos Phi r + cos Theta cos Phi theta - sin Phi phi
# y = sin Theta sin Phi r + cos Theta sin Phi theta - cos Phi phi
# z = cos Theta r - sin Theta theta
#
# It looks to me like I get the Theta convention I'm using by
# replacing sin Theta with cos Theta and cos Theta by sin Theta
# (because Dr. Cornelius takes 0 as the positive Z axis whereas
# I take zero as the X-Y plane) and changing the sign of theta
# (since Dr. Cornelius' Theta increases in the negative Z
# direction, whereas mine increases in the positive Z
# direction).
#
# The document was found at http://plaza.obu.edu/corneliusk/
# which is the page for Dr. Kevin Cornelius' Mathematical
# Physics (PHYS 4053) course at Ouachita Baptist University in
# Arkadelphia AR.
my $diag = sqrt( $cart_pos[0] * $cart_pos[0] +
$cart_pos[1] * $cart_pos[1] );
my ( $sin_theta, $cos_theta, $sin_phi, $cos_phi );
if ( $range > 0 ) {
$sin_theta = $cart_pos[2] / $range;
$cos_theta = $diag / $range;
if ( $diag > 0 ) {
$sin_phi = $cart_pos[1] / $diag;
$cos_phi = $cart_pos[0] / $diag;
} else {
$sin_phi = 0;
$cos_phi = 1;
}
} else {
$sin_theta = $sin_phi = 0;
$cos_theta = $cos_phi = 1;
}
# phi = - sin Phi x + cos phi y
# theta = - sin Theta cos Phi x - sin Theta sin Phi y + cos Theta z
# r = cos Theta cos Phi x + cos Theta sin Phi y + sin Theta z
my $az_hat = [ - $sin_phi, $cos_phi, 0 ];
my $el_hat = [ - $sin_theta * $cos_phi, - $sin_theta * $sin_phi,
$cos_theta ];
my $rng_hat = [ $cos_theta * $cos_phi, $cos_theta * $sin_phi,
$sin_theta ];
while ( @cart_data ) {
# Each triplet is then converted by projecting the Cartesian
# vector onto the appropriate unit vector. Azimuth and
# elevation are also converted to radians by dividing by the
# range. NOTE that this is the small-angle approximation,
# but since we assume we're dealing with derivitaves, it's
# OK.
my @cart_info = splice @cart_data, 0, 3;
push @rslt, vector_dot_product( $az_hat, \@cart_info ) / $range;
push @rslt, vector_dot_product( $el_hat, \@cart_info ) / $range;
push @rslt, vector_dot_product( $rng_hat, \@cart_info );
}
}
return @rslt;
}
# @spherical = _convert_cartesian_to_spherical_x( @cartesian )
#
# This subroutine is a convenience wrapper for
# _convert_cartesian_to_spherical(). The only difference is that the
# return is Theta (perpendicular to the X-Y plane), Phi (in the X-Y
# plane) and Rho, rather than Phi, Theta, Rho. This subroutine is my
# penance for not having all the methods that involve spherical
# coordinates return their arguments in the same order.
sub _convert_cartesian_to_spherical_x {
my @sph_data = _convert_cartesian_to_spherical( @_ );
for ( my $inx = 0; $inx < @sph_data; $inx += 3 ) {
@sph_data[ $inx, $inx + 1 ] = @sph_data[ $inx + 1, $inx ];
}
return @sph_data;
}
# @eci = $self->_convert_ecef_to_eci ()
# This subroutine converts the object's ECEF setting to ECI, and
# both caches and returns the result.
sub _convert_ecef_to_eci {
my $self = shift;
my $thetag = thetag ($self->universal);
my @data = $self->ecef ();
$self->{debug} and print <<eod;
Debug eci - thetag = $thetag
(x, y) = (@data[0, 1])
eod
my $costh = cos ($thetag);
my $sinth = sin ($thetag);
@data[0, 1] = ($data[0] * $costh - $data[1] * $sinth,
$data[0] * $sinth + $data[1] * $costh);
$self->{debug} and print <<eod;
Debug eci - after rotation,
(x, y) = (@data[0, 1])
eod
if ( @data > 5 ) {
@data[3, 4] = ($data[3] * $costh - $data[4] * $sinth,
$data[3] * $sinth + $data[4] * $costh);
$data[3] -= $data[1] * $self->{angularvelocity};
$data[4] += $data[0] * $self->{angularvelocity};
}
# A bunch of digging says this was added by Subversion commit 697,
# which was first released with Astro::Coord::ECI version 0.014. The
# comment says "Handle equinox in conversion between eci and ecef.
# Correctly, I hope." I'm leaving it in for now, but ...
$self->set (equinox_dynamical => $self->dynamical);
return @{$self->{_ECI_cache}{inertial}{eci} = \@data};
}
# This subroutine converts the object's ECI setting to ECEF, and
# both caches and returns the result.
our $EQUINOX_TOLERANCE = 365 * SECSPERDAY;
sub _convert_eci_to_ecef {
my $self = shift;
my $thetag = thetag ($self->universal);
my $dyn = $self->dynamical;
## my $equi = $self->get ('equinox_dynamical') || do {
## $self->set (equinox_dynamical => $dyn); $dyn};
my $equi = $self->{equinox_dynamical} ||= $dyn;
if (abs ($equi - $dyn) > $EQUINOX_TOLERANCE) {
$self->precess_dynamical ($dyn);
}
my @ecef = $self->eci ();
my $costh = cos (- $thetag);
my $sinth = sin (- $thetag);
@ecef[0, 1] = ($ecef[0] * $costh - $ecef[1] * $sinth,
$ecef[0] * $sinth + $ecef[1] * $costh);
if ( @ecef > 3 ) {
@ecef[ 3, 4 ] = ( $ecef[3] * $costh - $ecef[4] * $sinth,
$ecef[3] * $sinth + $ecef[4] * $costh );
$ecef[3] += $ecef[1] * $self->{angularvelocity};
$ecef[4] -= $ecef[0] * $self->{angularvelocity};
}
$self->{_ECI_cache}{fixed}{ecef} = \@ecef;
defined wantarray
or return;
return @ecef;
}
sub _convert_eci_to_heliocentric_ecliptic_cartesian {
my ( $self ) = @_;
my ( $x, $y, $z ) = $self->eci();
# $x = km in the direction of the Vernal Equinox
# $y = km in the direction 90 degrees east of the Vernal Equinox
# $z = km perpendicular to the equator.
# So we need to rotate z toward y by the obliquity of the ecliptic
my $dynamical = $self->dynamical();
my $epsilon = obliquity( $dynamical );
my $sin_epsilon = sin $epsilon;
my $cos_epsilon = cos $epsilon;
my $Z = $z * $cos_epsilon - $y * $sin_epsilon;
my $Y = $z * $sin_epsilon + $y * $cos_epsilon;
my $sun = $self->get( 'sun' )
or croak q{Attribute 'sun' not set};
my ( $X_s, $Y_s, $Z_s ) = $sun->universal( $self->universal() )->eci();
my @hec = ( $x - $X_s, $Y - $Y_s, $Z - $Z_s );
$self->set( equinox_dynamical => $dynamical );
$self->{_ECI_cache}{inertial}{heliocentric_ecliptic_cartesian} = \@hec;
defined wantarray
or return;
return @hec;
}
sub _convert_heliocentric_ecliptic_cartesian_to_eci {
my ( $self ) = @_;
my ( $x, $y, $z ) = $self->heliocentric_ecliptic_cartesian();
# $x = km in the direction of the Vernal Equinox
# $y = km in the direction 90 degrees east of the Vernal Equinox
# $z = km perpendicular to the ecliptic.
# So we need to rotate z toward y by the obliquity of the ecliptic
my $dynamical = $self->dynamical();
my $epsilon = obliquity( $dynamical );
my $sin_epsilon = sin $epsilon;
my $cos_epsilon = cos $epsilon;
my $Z = $z * $cos_epsilon + $y * $sin_epsilon;
my $Y = - $z * $sin_epsilon + $y * $cos_epsilon;
my $sun = $self->get( 'sun' )
or croak q{Attribute 'sun' not set};
my ( $X_s, $Y_s, $Z_s ) = $sun->universal( $self->universal() )->eci();
my @eci = ( $x + $X_s, $Y + $Y_s, $Z + $Z_s );
$self->set( equinox_dynamical => $dynamical );
$self->{_ECI_cache}{inertial}{eci} = \@eci;
defined wantarray
or return;
return @eci;
}
# my @args = _expand_args_default_station( @_ )
#
# This subroutine handles the placing of the contents of the
# 'station' attribute into the argument list of methods that,
# prior to the introduction of the 'station' attribute, formerly
# took two Astro::Coord::ECI objects and computed the position of
# the second as seen from the first.
sub _expand_args_default_station {
my @args = @_;
if ( ! embodies( $args[1], 'Astro::Coord::ECI' ) ) {
unshift @args, $args[0]->get( 'station' );
embodies( $args[0], 'Astro::Coord::ECI' )
or croak 'Station not set';
defined $args[1]->{universal}
and $args[0]->universal( $args[1]->universal() );
}
return @args;
}
# $value = $self->__initial_inertial
#
# Return the initial setting of the inertial attribute. At this
# level we are assumed not inertial until we acquire a position.
# This is not part of the public interface, but may be used by
# subclasses to set an initial value for this read-only attribute.
# Setting the coordinates explicitly will still set the {inertial}
# attribute appropriately.
sub __initial_inertial{ return };
# $value = _local_mean_delta ($coord)
# Calculate the delta from universal to civil time for the object.
# An exception is raised if the coordinates of the object have not
# been set.
sub _local_mean_delta {
return ($_[0]->geodetic ())[1] * SECSPERDAY / TWOPI;
}
=begin comment
# $string = _rad2dms ($angle)
# Convert radians to degrees, minutes, and seconds of arc.
# Used for debugging.
sub _rad2dms {
my $angle = rad2deg (shift);
my $deg = floor ($angle);
$angle = ($angle - $deg) * 60;
my $min = floor ($angle);
$angle = ($angle - $min) * 60;
return "$deg degrees $min minutes $angle seconds of arc";
}
=end comment
=cut
#######################################################################
#
# Package initialization
#
__PACKAGE__->set (ellipsoid => 'WGS84');
%savatr = map {$_ => 1} (keys %static, qw{dynamical id name universal});
# Note that local_mean_time does not get preserved, because it
# changes if the coordinates change.
1;
__END__
=back
=head2 Attributes
This class has the following attributes:
=over
=item almanac_horizon (radians or string)
This attribute specifies the horizon used for almanac calculations, as
radians above or below the plans of the observer. It can also be set to
the following string:
* C<height> adjusts the horizon in proportion to the observer's height
above sea level. This assumes a spherical Earth and an unobstructed
horizon.
When doing almanac calculations, it is the C<almanac_horizon> of the
observer that is used, not the C<almanac_horizon> of the observed body.
The default is C<0>.
=item angularvelocity (radians per second)
This attribute represents the angular velocity of the Earth's surface in
radians per second. The initial value is 7.292114992e-5, which according
to Jean Meeus is the value for 1996.5. He cites the International Earth
Rotation Service's Annual Report for 1996 (Published at the Observatoire
de Paris, 1997).
Subclasses may place appropriate values here, or provide a period()
method.
=item debug (numeric)
This attribute turns on debugging output. The only supported value of
this attribute is 0. That is to say, the author makes no guarantees of
what will happen if you set it to some other value, nor does he
guarantee that this behavior will not change from release to release.
The default is 0.
=item diameter (numeric, kilometers)
This attribute exists to support classes/instances which represent
astronomical bodies. It represents the diameter of the body, and is
used by the azel() method when computing the upper limb of the body.
It has nothing to do with the semimajor attribute, which always refers
to the Earth, and is used to calculate the latitude and longitude of
the body.
The default is 0.
=item edge_of_earths_shadow (numeric, dimensionless)
This attribute represents the location of the edge of the Earth's shadow
cast by the illuminating body, in terms of the apparent radius of that
body. That is, the edge of the umbra is specified by 1, the middle of
the penumbra by 0, and the edge of the penumbra by -1. It was added for
the benefit of the B<Astro::Coord::ECI::TLE> class, on the same dubious
logic that the L<twilight|/twilight> attribute was added.
The default is 1.
=item ellipsoid (string)
This attribute represents the name of the reference ellipsoid to use.
It must be set to one of the known ellipsoids. If you set this,
flattening and semimajor will be set also. See the documentation to
the known_ellipsoid() method for the initially-valid names, and how
to add more.
The default is 'WGS84'.
=item equinox_dynamical (numeric, dynamical time)
This attribute represents the time of the L</Equinox> to which the
coordinate data are referred. Models implemented by subclasses should
set this to the L</Equinox> to which the model is referred. When setting
positions directly the user should also set the desired
equinox_dynamical if conversion between inertial and Earth-fixed
coordinates is of interest. If this is not set, these conversions will
use the current time setting of the object as the L</Equinox>.
In addition, if you have a position specified in earth-fixed coordinates
and convert it to inertial coordinates, this attribute will be set to
the dynamical time of the object.
Unlike most times in this package, B<this attribute is specified in
dynamical time.> If your equinox is universal time $uni, set this
attribute to $uni + dynamical_delta ($uni). The dynamical_delta
subroutine is found in Astro::Coord::ECI::Utils.
The default is undef.
=item flattening (numeric)
This attribute represents the flattening factor of the reference
ellipsoid. If you set the ellipsoid attribute, this attribute will be
set to the flattening factor for the named ellipsoid. If you set this
attribute, the ellipsoid attribute will become undefined.
The default is appropriate to the default ellipsoid.
=item frequency (numeric, Hertz)
This attribute represents the frequency on which the body transmits, in
Hertz. If specified, the C<azel()> method will return the Doppler shift
if velocities are available.
The default is C<undef>.
=item horizon (numeric, radians)
This attribute represents the distance the effective horizon is above
the geometric horizon. It was added for the
B<Astro::Coord::ECI::TLE::Iridium> class, on the same dubious logic
that the L<twilight|/twilight> attribute was added.
The default is the equivalent of 20 degrees.
=item id (string)
This is an informational attribute, and its setting (or lack thereof)
does not affect the functioning of the class. Certain subclasses will
set this when they are instantiated. See the subclass documentation
for details.
=item inertial (boolean, read-only)
This attribute returns true (in the Perl sense) if the object was
most-recently set to inertial coordinates (i.e. eci, ecliptic, or
equatorial) and false otherwise. If coordinates have not been set,
it is undefined (and therefore false).
=item name (string)
This is an informational attribute, and its setting (or lack thereof)
does not affect the functioning of the class. Certain subclasses will
set this when they are instantiated. See the subclass documentation
for details.
=item refraction (boolean)
Setting this attribute to a true value includes refraction in the
calculation of the azel() method. If set to a false value, atmospheric
refraction is ignored.
The default is true (well, 1 actually).
=item specified (string, read-only)
This attribute records the coordinate system in which the coordinates
of the object were set.
=item semimajor (numeric, kilometers)
This attribute represents the semimajor axis of the reference
ellipsoid. If you set the ellipsoid attribute, this attribute will be
set to the semimajor axis for the named ellipsoid. If you set this
attribute, the ellipsoid attribute will become undefined.
For subclasses representing bodies other than the Earth, this attribute
will be set appropriately.
The default is appropriate to the default ellipsoid.
=item station ( Astro::Coord::ECI object )
This attribute represents the default observing location for relative
coordinate systems. It must be set to an C<Astro::Coord::ECI> object, or
C<undef>.
The intent is that all methods which compute positions or events as seen
from a user-specified location should make use of this. In other words,
if you find one that does not, you have found a bug.
Cascading station objects are B<not> supported. That is, if you have
C<Astro::Coord::ECI> objects C<$a>, C<$b>, C<$c> and so on, if you
$a->set( station => $b ),
then C<< $b->set( station => $c ) >> is not supported, nor is C<<
$c->set( station => $a ) >>. Because of the work involved in preventing
this in general, I have decided to rely on just politely recommending
that It Is Better Not. But specific pernicious cases B<will> throw
exceptions, and I reserve the right to do the work and throw exceptions
in all cases if it becomes a problem.
=item twilight (numeric, radians)
This attribute represents the elevation of the center of the Sun's disk
at the beginning and end of twilight.
Some of the usual values are:
civil twilight: -6 degrees
nautical twilight: -12 degrees
astronomical twilight: -18 degrees
B<Note> that the documentation at this point used to say (or at least
imply) that the strings C<'civil'>, C<'nautical'>, and C<'astronomical'>
were acceptable. This has never been the case to my knowledge, so since
I have received no bug reports, I have considered the bug to be in the
documentation, and removed the relevant text as of version 0.051_01.
The default is -6 degrees (or, actually, the equivalent in radians).
=back
=head1 A NOTE ON VELOCITIES
The velocities returned by the methods in this class are a work in
progress. Velocities set are always returned. But conversions are more
problematic. For one thing, I have not come across decent worked
examples, so I am limited to sanity checking.
I consider a velocity to be sane if the difference between the position
at time C<t> and the position at time C<t + 1> lies between the velocity
at time C<t> and the velocity at time C<t + 1>.
In order to understand what velocity conversions are sane, you need to
understand how coordinate conversions work internally. I will try to
make this plain with some lame ASCII art, which depicts the
relationships among those coordinate systems where velocity conversion
appears to be at least mostly sane in the above sense of the word:
ecliptic()
equatorial()
^
|
eci()
^
|
v
ecef() geocentric() geodetic()
|
v
neu() -> enu()
|
v
azel()
If the object's position and velocity were set in one set of units,
other units are obtained by following the arrows. The arrows below
ecef() are one-way because it is not currently possible to set a
position in these units. Similarly, the arrow from C<eci()> to
C<equatorial()> is one-way because there is currently no way to set an
equatorial velocity. There are no arrows to C<geocentric()>,
C<geodetic()> and C<ecliptic()> because these conversions do not support
velocities.
=head1 TERMINOLOGY AND CONVENTIONS
Partly because this module straddles the divide between geography and
astronomy, the establishment of terminology and conventions was a
thorny and in the end somewhat arbitrary process. Because of this,
documentation of salient terms and conventions seemed to be in order.
=head2 Altitude
This term refers to the distance of a location above mean sea level.
Altitude input to and output from this module is in kilometers.
Maps use "elevation" for this quantity, and measure it in meters. But
we're using L</Elevation> for something different, and I needed
consistent units.
=head2 Azimuth
This term refers to distance around the horizon measured clockwise from
North.
Azimuth output from this module is in radians.
Astronomical sources tend to measure from the South, but I chose the
geodetic standard, which seems to be usual in satellite tracking
software.
=head2 Declination
Declination is the angle a point makes with the plane of the equator
projected onto the sky. North declination is positive, south
declination is negative.
Declination input to and output from this module is in radians.
=head2 Dynamical time
A dynamical time is defined theoretically by the motion of astronomical
bodies. In practice, it is seen to be related to Atomic Time (a.k.a.
TAI) by a constant.
There are actually two dynamical times of interest: TT (Terrestrial
Time, a.k.a. TDT for Terrestrial Dynamical Time), which is defined in
terms of the geocentric ephemerides of solar system bodies, and TDB
(Barycentric Dynamical Time), which is defined in terms of the
barycentre (a.k.a "center of mass") of the solar system. The two differ
by the relativistic effects of the motions of the bodies in the Solar
system, and are generally less than 2 milliseconds different. So unless
you are doing high-precision work they can be considered identical, as
Jean Meeus does in "Astronomical Algorithms".
For practical purposes, TT = TAI + 32.184 seconds. If I ever get the
gumption to do a re-implementation (or alternate implementation) of time
in terms of the DateTime object, this will be the definition of
dynamical time. Until then, though, formula 10.2 on page 78 of Jean
Meeus' "Astronomical Algorithms" second edition, Chapter 10 (Dynamical
Time and Universal Time) is used.
Compare and contrast this to L</Universal time>. This explanation leans
heavily on L<http://star-www.rl.ac.uk/star/docs/sun67.htx/node226.html>,
which contains a more fulsome but eminently readable explanation.
=head2 Earth-Centered, Earth-fixed (ECEF) coordinates
This is a Cartesian coodinate system whose origin is the center of the
reference ellipsoid. The X axis passes through 0 degrees L</Latitude>
and 0 degrees L</Longitude>. The Y axis passes through 90 degrees east
L</Latitude> and 0 degrees L</Longitude>, and the Z axis passes through
90 degrees north L</Latitude> (a.k.a the North Pole).
All three axes are input to and output from this module in kilometers.
Also known as L</XYZ coordinates>, e.g. at
L<http://www.ngs.noaa.gov/cgi-bin/xyz_getxyz.prl>
=head2 Earth-Centered Inertial (ECI) coordinates
This is the Cartesian coordinate system in which NORAD's models predict
the position of orbiting bodies. The X axis passes through 0 hours
L</Right Ascension> and 0 degrees L</Declination>. The Y axis passes
through 6 hours L</Right Ascension> and 0 degrees L</Declination>. The
Z axis passes through +90 degrees L</Declination> (a.k.a. the North
Pole). By implication, these coordinates are referred to a given
L</Equinox>.
All three axes are input to and output from this module in kilometers.
=head2 Ecliptic
The Ecliptic is the plane of the Earth's orbit, projected onto the sky.
Ecliptic coordinates are a spherical coordinate system referred to
the ecliptic and expressed in terms of L</Ecliptic latitude> and
L</Ecliptic longitude>. By implication, Ecliptic coordinates are also
referred to a specific L</Equinox>.
=head3 Ecliptic latitude
Ecliptic longitude is the angular distance of a point above the plane
of the Earth's orbit.
Ecliptic latitude is input to and output from this module in radians.
=head3 Ecliptic longitude
Ecliptic longitude is the angular distance of a point east of the point
where the plane of the Earth's orbit intersects the plane of the
equator. This point is also known as the vernal L</Equinox> and the
first point of Ares.
Ecliptic longitude is input to and output from this module in radians.
=head2 Elevation
This term refers to an angular distance above the horizon.
Elevation output from this module is in radians.
This is the prevailing meaning of the term in satellite tracking.
Astronomers use "altitude" for this quantity, and call the
corresponding coordinate system "altazimuth." But we're using
L</Altitude> for something different.
=head2 Equatorial
Equatorial coordinates are a spherical coordinate system referred to
the plane of the equator projected onto the sky. Equatorial coordinates
are specified in L</Right Ascension> and L</Declination>, and implicitly
referred to a given L</Equinox>
=head2 Equinox
The L</Ecliptic>, L</Equatorial>, and L</Earth-Centered Inertial (ECI)
coordinates> are defined in terms of the location of the intersection of
the celestial equator with the L</Ecliptic>. The actual location of this
point changes in time due to precession of the Earth's axis, so each of
these coordinate systems is implicitly qualified by ("referred to"
appears to be the usual terminology) the relevant time. By a process of
association of ideas, this time is referred to as the equinox of the
data.
=head2 Geocentric
When referring to a coordinate system, this term means that the
coordinate system assumes the Earth is spherical.
=head3 Geocentric latitude
Geocentric latitude is the angle that the ray from the center of the
Earth to the location makes with the plane of the equator. North
latitude is positive, south latitude is negative.
Geocentric latitude is input to and output from this module in radians.
=head2 Geodetic
When referring to a coordinate system, this term means that the
coordinate system assumes the Earth is an ellipsoid of revolution
(or an oblate spheroid if you prefer). A number of standard
L</Reference Ellipsoids> have been adopted over the years.
=head3 Geodetic latitude
Geodetic latitude is the latitude found on maps. North latitude is
positive, south latitude is negative.
Geodetic latitude is input to and output from this module in radians.
Technically speaking, Geodetic latitude is the complement of the angle
the plane of the horizon makes with the plane of the equator. In
this software, the plane of the horizon is determined from a
L</Reference Ellipsoid>.
=head2 Latitude
See either L</Ecliptic latitude>, L</Geocentric latitude> or
L</Geodetic latitude>. When used without qualification, L</Geodetic
latitude> is meant.
=head2 Longitude
When unqualified, this term refers to the angular distance East or West
of the standard meridian. East longitude is positive, West longitude is
negative.
Longitude is input to or output from this module in radians.
For L</Ecliptic longitude>, see that entry.
Jean Meeus reports in "Astronomical Algorithms" that the International
Astronomical Union has waffled on the sign convention. I have taken
the geographic convention.
=head2 Obliquity (of the Ecliptic)
This term refers to the angle the plane of the equator makes with the
plane of the Earth's orbit.
Obliquity is output from this module in radians.
=head2 Reference Ellipsoid
This term refers to a specific ellipsoid adopted as a model of the
shape of the Earth. It is defined in terms of the equatorial radius
in kilometers, and a dimensionless flattening factor which is used
to derive the polar radius, and defined such that the flattening
factor of a sphere is zero.
See the documentation on the reference_ellipsoid() method for a list
of reference ellipsoids known to this class.
=head2 Right Ascension
This term refers to the angular distance of a point east of the
intersection of the plane of the Earth's orbit with the plane of
the equator.
Right Ascension is input to and output from this module in radians.
In astronomical literature it is usual to report right ascension
in hours, minutes, and seconds, with 60 seconds in a minute, 60
minutes in an hour, and 24 hours in a circle.
=head2 Universal time
This term can refer to a number of scales, but the two of interest are
UTC (Coordinated Universal Time) and UT1 (Universal Time 1, I presume).
The latter is in effect mean solar time at Greenwich, though its
technical definition differs in detail from GMT (Greenwich Mean Time).
The former is a clock-based time, whose second is the SI second (defined
in terms of atomic clocks), but which is kept within 0.9 seconds of UT1
by the introduction of leap seconds. These are introduced (typically at
midyear or year end) by prior agreement among the various timekeeping
bodies based on observation; there is no formula for computing when a
leap second will be needed, because of irregularities in the Earth's
rotation.
Jean Meeus' "Astronomical Algorithms", second edition, deals with the
relationship between Universal time and L</Dynamical time> in Chapter 10
(pages 77ff). His definition of "Universal time" seems to refer to UT1,
though he does not use the term.
This software considers Universal time to be equivalent to Perl time.
Since we are interested in durations (time since a given epoch, to be
specific), this is technically wrong in most cases, since leap seconds
are not taken into account. But in the case of the bodies modeled by
the Astro::Coord::ECI::TLE object, the epoch is very recent (within a
week or so), so the error introduced is small. It is larger for
astronomical calculations, where the epoch is typically J2000.0, but the
angular motions involved are smaller, so it all evens out. I hope.
Compare and contrast L</Dynamical time>. This explanation leans heavily
on L<http://star-www.rl.ac.uk/star/docs/sun67.htx/node224.html>.
=head2 XYZ coordinates
See L</Earth-Centered, Earth-fixed (ECEF) coordinates>.
=head1 ACKNOWLEDGMENTS
The author wishes to acknowledge and thank the following individuals and
organizations.
Kazimierz Borkowski, whose "Accurate Algorithms to Transform
Geocentric to Geodetic Coordinates", at
L<http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm>,
was used for transforming geocentric to geodetic coordinates.
Dominik Brodowski (L<http://www.brodo.de/>), whose SGP C-lib
(available at L<http://www.brodo.de/space/sgp/>) provided a
reference implementation that I could easily run, and pick
apart to help get B<Astro::Coord::ECI::TLE> working. Dominik based
his work on Dr. Kelso's Pascal implementation.
Felix R. Hoots and Ronald L. Roehric, the authors of "SPACETRACK
REPORT NO. 3 - Models for Propagation of NORAD Element Sets,"
which provided the basis for the Astro::Coord::ECI::TLE module.
T. S. Kelso, who compiled this report and made it available at
L<http://celestrak.com/>, whose "Computers and Satellites" columns
in "Satellite Times" magazine were invaluable to an understanding
and implementation of satellite tracking software, whose support,
encouragement, patience, and willingness to provide answers on arcane
topics were a Godsend, who kindly granted permission to use his
azimuth/elevation algorithm, and whose Pascal implementation of the
NORAD tracking algorithms indirectly provided a reference
implementation for me to use during development.
John A. Magliacane, whose open-source C<Predict> program, available at
L<http://www.qsl.net/kd2bd/predict.html>, gave me a much-needed leg up
on the velocity calculations in the C<azel()> method.
Jean Meeus, whose book "Astronomical Algorithms" (second edition)
formed the basis for this module and the B<Astro::Coord::ECI::Sun>,
B<Astro::Coord::ECI::Moon>, and B<Astro::Coord::ECI::Star> modules,
and without whom it is impossible to get far in computational
astronomy. Any algorithm not explicitly credited above is probably
due to him.
Dr. Meeus' publisher, Willmann-Bell Inc (L<http://www.willbell.com/>),
which kindly and patiently answered my intellectual-property questions.
Goran Gasparovic of MIT, who asked for (and supplied information for)
the ability to display results in apparent equatorial coordinates,
rather than azimuth and elevation.
Dr. Kevin Cornelius of Ouachita Baptist University, whose handout
"Relationships Among Unit Vectors" for his Mathematical Physics class
(PHYS 4053) provided a much surer and more expeditious way to handling
velocity conversion from spherical to Cartesian (and vice versa) than my
own ancient and rickety matrix math.
=head1 BUGS
Functionality involving velocities is B<untested>, and is quite likely
to be wrong.
Bugs can be reported to the author by mail, or through
L<http://rt.cpan.org/>.
=head1 SEE ALSO
The L<Astro|Astro> package by Chris Phillips. This contains three
function-based modules: L<Astro::Coord|Astro::Coord>, which provides
various astronomical coordinate conversions, plus the calculation of
various ephemeris variables; L<Astro::Time|Astro::Time> contains time
and unit conversions, and L<Astro::Misc|Astro::Misc> contains various
calculations unrelated to position and time.
The B<Astro-Coords> package by Tim Jenness. This contains various modules
to do astronomical calculations, and includes coordinate conversion and
calculation of planetary orbits based on orbital elements. Requires
SLALIB from L<http://www.starlink.rl.ac.uk/Software/software_store.htm>.
The L<Astro::Nova|Astro::Nova> module by Steffen Mueller, which wraps
(and bundles) the libnova celestial mechanics, astrometry and
astrodynamics library found at L<http://libnova.sourceforge.net/>.
The L<http://www.heavens-above.com/> web site, which does B<not> use
this code, but does provide similar functionality.
=head1 AUTHOR
Thomas R. Wyant, III (F<wyant at cpan dot org>)
=head1 COPYRIGHT AND LICENSE
Copyright (C) 2005-2018 by Thomas R. Wyant, III
This program is free software; you can redistribute it and/or modify it
under the same terms as Perl 5.10.0. For more details, see the full text
of the licenses in the directory LICENSES.
This program is distributed in the hope that it will be useful, but
without any warranty; without even the implied warranty of
merchantability or fitness for a particular purpose.
=cut
# ex: set textwidth=72 :