Astro-satpass
view release on metacpan or search on metacpan
lib/Astro/Coord/ECI.pm view on Meta::CPAN
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 => CIVIL_TWILIGHT,
);
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();
ref $static{sun}
and $self->set( sun => $static{sun}->clone() );
@args and $self->set( @args );
exists $self->{almanac_horizon}
or $self->set( almanac_horizon => 0 );
$self->{debug} and do {
require Data::Dumper;
local $Data::Dumper::Terse = 1;
print "Debug - Instantiated ", Data::Dumper::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 the Math Forum,
C<https://www.nctm.org/tmf/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.
Unfortunately, as of early 2021 the National Council of Teachers of
Mathematics restricted the Dr. Math content to their members, but an
annotated and expanded version of the article on haversines is available
at
L<https://www.themathdoctors.org/distances-on-earth-2-the-haversine-formula/>.
If you want the original article, you can feed the URL
C<http://mathforum.org/library/drmath/view/51879.html> to the Wayback
Machine.
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 $angle = $coord->angular_radius();
This method computes the angular radius of the C<$coord> body in
radians, as of the currently-set time. If the C<station> attribute is
set, the range will be computed from it via C<< $self->azel() >>. If
not, the range from the center of the earth will be computed via
C<< $self->ecliptic() >>. If the C<diameter> attribute is not set,
C<undef> is returned.
This method is not supported for terrestrial objects.
=cut
lib/Astro/Coord/ECI.pm view on Meta::CPAN
=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 Clone::clone.
=cut
sub clone {
my ( $self ) = @_;
return Clone::clone( $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
C<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 ();
lib/Astro/Coord/ECI.pm view on Meta::CPAN
$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 $seconds = $self->equation_of_time( $time );
This method returns the equation of time at the given B<dynamical>
time. If the time is C<undef>, the invocant's dynamical time is used.
The algorithm is from W. S. Smart's "Text-Book on Spherical Astronomy",
as reported in Jean Meeus' "Astronomical Algorithms", 2nd Edition,
Chapter 28, page 185.
=cut
sub equation_of_time {
my ( $self, $time ) = @_;
if ( looks_like_number( $self ) ) {
( $self, $time ) = ( __PACKAGE__, $self );
}
defined $time
or $time = $self->dynamical();
my $epsilon = $self->obliquity( $time );
my $y = tan($epsilon / 2);
$y *= $y;
# The following algorithm is from Meeus, chapter 25, page, 163 ff.
my $T = jcent2000( $time ); # Meeus (25.1)
my $L0 = mod2pi( deg2rad( (.0003032 * $T + 36000.76983 ) * $T
+ 280.46646 ) ); # Meeus (25.2)
my $M = mod2pi( deg2rad( ( ( -.0001537 ) * $T + 35999.05029 )
* $T + 357.52911 ) ); # Meeus (25.3)
my $e = ( -0.0000001267 * $T - 0.000042037 ) * $T
+ 0.016708634; # Meeus (25.4)
my $E = $y * sin( 2 * $L0 ) - 2 * $e * sin( $M ) +
4 * $e * $y * sin( $M ) * cos( 2 * $L0 ) -
$y * $y * .5 * sin( 4 * $L0 ) -
1.25 * $e * $e * sin( 2 * $M ); # Meeus (28.3)
return $E * SECSPERDAY / TWOPI; # The formula gives radians.
}
=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.
lib/Astro/Coord/ECI.pm view on Meta::CPAN
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 C<http://star-www.rl.ac.uk/star/docs/sun67.htx/node226.html>,
which contains a more fulsome but eminently readable explanation.
Unfortunately that link has gone the way of the dodo, and I have been
unable to find an adequate replacement.
=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<https://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
lib/Astro/Coord/ECI.pm view on Meta::CPAN
=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 C<http://star-www.rl.ac.uk/star/docs/sun67.htx/node224.html>.
Unfortunately that link has gone the way of the dodo, and I have been
unable to find an adequate replacement.
=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<https://www.uni-saarland.de/lehrstuhl/brodowski/team/prof-dr-dominik-brodowski-llm-upenn.html>), whose SGP C-lib
(available at L<https://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.org/>, 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<https://www.qsl.net/kd2bd/predict.html>, gave me a much-needed leg up
on the velocity calculations in the C<azel()> method.
( run in 0.743 second using v1.01-cache-2.11-cpan-d7f47b0818f )