|
Astro::Coord::ECI - Manipulate geocentric coordinates |
Astro::Coord::ECI - Manipulate geocentric coordinates
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);
# 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";
# Instantiate a satellite. The parser returns a list.
my ($tle) = Astro::Coord::ECI::TLE->parse (<<eod);
Pretend this text contains NORAD orbital element data for
a satellite. If it contains more than one set of data, you
get more than one object back.
eod
# Figure out where the satellite is in the observer's sky.
# Note that setting the time runs the model.
my ($right_asc, $decl, $rng) =
$loc->equatorial ($tle->universal ($time));
# Find the satellite's ECI coordinates at the given time.
# It is not necessary to run the model again, unless you
# want coordinates for a different time.
my ($X, $Y, $Z) = $tle->eci ();
print <<eod;
Satellite position
equatorial coordinates:
right ascension: @{[rad2deg ($right_asc)]} degr
declination: @{[rad2deg ($decl)]} degr
range: $rng km
ECI coordinates:
X: $X
Y: $Y
Z: $Z
eod
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 except for the Julian Day routines, which convert Perl time into things like Julian Day or days since Julian 2000.0 (i.e. the aforementioned support routines), and for a few support routines which take or return dynamical times. Time is specified and returned in universal time unless the documentation of the method explicitly states otherwise. In practice, universal time is approximated by normal Perl time.
Known subclasses include Astro::Coord::ECI::Moon to predict the position of the Moon, Astro::Coord::ECI::Star to predict the position of a star, or anything else that can be considered fixed on the celestial sphere, Astro::Coord::ECI::Sun to predict the position of the Sun, and Astro::Coord::ECI::TLE to predict the position of a satellite given the NORAD orbital parameters.
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.
The following methods should be considered public.
set() method once the object has been instantiated.
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 'upper' argument is true, the calculation will be of the upper limb of the object, using the 'diameter' attribute of the $coord2 object.
As a side effect, the time of the $coord object may be set from the $coord object.
If the refraction attribute of the $coord object is
true, the elevation will be corrected for atmospheric refraction using
the correct_for_refraction() method.
The basic algorithm comes from T. S. Kelso's ``Computers and Satellites'' column in ``Satellite Times'', November/December 1995, titled ``Orbital Coordinate Systems, Part II'' and available at http://celestrak.com/columns/v02n02/. If the object represents fixed coordinates, the author's algorithm is used, but the author confesses needing to refer to Dr. Kelso's work to get the signs right.
It's really just a wrapper for Storable::dclone.
This method can also be called as a class method. It is really only
exposed for testing purposes (hence the cumbersome name). The azel()
method calls it for you if the 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.
The algorithm is simple enough to be the author's.
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.
The algorithm comes from Jean Meeus' ``Astronomical Algorithms'', 2nd Edition, Chapter 10, pages 78ff.
This method sets the coordinates represented by the object in terms of Earth-Centered, Earth-fixed (ECEF) coordinates, with x being latitude 0 longitude 0, y being latitude 0 longitude 90 degrees east, and z being latitude 90 degrees north. The velocities are optional, and will default to the rotational velocity at the point being set. The object itself is returned.
This method can also be called as a class method, in which case it instantiates the desired object.
Caveat: Velocities are also returned, but should not at this point be taken seriously unless they were originally set by the same method that is returning them, since I have not at this point got the velocity transforms worked out.
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 http://celestrak.com/columns/v02n01/.
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()
Caveat: Velocities are also returned, but should not at this point be taken seriously unless they were originally set by the same method that is returning them, since I have not at this point got the velocity transforms worked out.
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.
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.
universal() or
dynamical() methods).
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.
universal() or
dynamical() methods).
As a side effect, the time of the $coord object is set from the $coord object.
If the refraction attribute of the $coord object is
true, the coordinates will be corrected for atmospheric refraction using
the correct_for_refraction() method.
The algorithm is lifted pretty much verbatim from the Calculate_RADec subroutine of SGP4 Pascal Library Version 2.65 by T. S. Kelso, available at http://celestrak.com/software/tskelso-sw.asp. Dr. Kelso credits the algorithm to ``Methods of Orbit Determination'' by Pedro Ramon Escobal, pages 401 - 402.
This method can also be called as a class method, in which case it instantiates the desired object.
This method should not be used with map coordinates because map latitude is 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.
The ellipsoid argument is the name of a 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.
This is the method that should be used with map coordinates.
The ellipsoid argument is the name of a 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 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.
If called as a class method, it returns the current default values.
See Attributes for a list of the attributes you can get.
period() method, this method
just returns two pi divided by the period. Otherwise it returns the
contents of the angularvelocity attribute.
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.
The optional $want argument should be specified as true (i.e. 1) if you want the next passage above the observer, or as false (i.e. 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 $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 $body represents the sun), or false if below (i.e. local midnight if the $body represents the sun). If called in scalar context, you get the time only.
The current time of both $coord and $body object are left at the returned time.
The algorithm is by successive approximation. It will croak if the period of the $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.
NOTE that side effects of setting the time (i.e. in subclasses which
define the time_set() method) will take place as a result of calling
this method.
The algorithm comes from Jean Meeus, ``Astronomical Algorithms'', 2nd Edition, Chapter 21, pages 134ff (a.k.a. ``the rigorous method'').
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.
reference_ellipsoid($name)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.
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.
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 Attributes for a list of the attributes you can set.
This method can also be called as a class method, in which case it instantiates the desired object.
This class has the following attributes:
Subclasses may place appropriate values here, or provide a period()
method.
The default is 0.
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.
known_ellipsoid() method for the initially-valid names, and how
to add more.
The default is 'WGS84'.
The default is appropriate to the default ellipsoid.
azel() method. If set to a false value, atmospheric
refraction is ignored.
The default is true (well, 1 actually).
For subclasses representing bodies other than the Earth, this attribute will be set appropriately.
The default is appropriate to the default ellipsoid.
Some of the usual values are:
civil twilight: -6 degrees nautical twilight: -12 degrees astronomical twilight: -18 degrees
The default is -6 degrees (or, actually, the equivalent in radians).
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.
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 Elevation for something different, and I needed consistent units.
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.
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.
This is a Cartesian coodinate system whose origin is the center of the reference ellipsoid. The X axis passes through 0 degrees Latitude and 0 degrees Longitude. The Y axis passes through 90 degrees east Latitude and 0 degrees Longitude, and the Z axis passes through 90 degrees north Latitude (a.k.a the North Pole).
All three axes are input to and output from this module in kilometers.
Also known as XYZ coordinates, e.g. at http://www.ngs.noaa.gov/cgi-bin/xyz_getxyz.prl
This is the Cartesian coordinate system in which NORAD's models predict the position of orbiting bodies. The X axis passes through 0 hours Right Ascension and 0 degrees Declination. The Y axis passes through 6 hours Right Ascension and 0 degrees Declination. The Z axis passes through +90 degrees Declination (a.k.a. the North Pole).
All three axes are input to and output from this module in kilometers.
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 Ecliptic latitude and Ecliptic longitude.
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.
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 equinox and the first point of Ares.
Ecliptic longitude is input to and output from this module in radians.
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 Altitude for something different.
Equatorial coordinates are a spherical coordinate system referred to the plane of the equator projected onto the sky. Equatorial coordinates are specified in Right Ascension and Declination.
When referring to a coordinate system, this term means that the coordinate system assumes the Earth is spherical.
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.
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 Reference Ellipsoids have been adopted over the years.
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 Reference Ellipsoid.
See either Ecliptic latitude, Geocentric latitude or Geodetic latitude. When used without qualification, Geodetic latitude is meant.
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 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.
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.
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.
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.
See Earth-Centered, Earth-fixed (ECEF) coordinates.
The author wishes to acknowledge the following individuals and organizations.
Kazimierz Borkowski, whose ``Accurate Algorithms to Transform Geocentric to Geodetic Coordinates'', at http://www.astro.uni.torun.pl/~kb/Papers/geod/Geod-BG.htm, was used for transforming geocentric to geodetic coordinates.
Dominik Brodowski (http://www.brodo.de/), whose SGP C-lib (available at http://www.brodo.de/space/sgp/) provided a reference implementation that I could easily run, and pick apart to help get 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 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.
Jean Meeus, whose book ``Astronomical Algorithms'' (second edition) formed the basis for this module and the Astro::Coord::ECI::Sun, Astro::Coord::ECI::Moon, and 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 (http://www.willbell.com/), which kindly and patiently answered my intellectual-property questions.
Functionality involving velocities is untested, and is quite likely to be wrong.
Bugs can be reported to the author by mail, or through http://rt.cpan.org/.
The Astro package by Chris Phillips. This contains three function-based modules: Astro::Coord, which provides various astronomical coordinate conversions, plus the calculation of various ephemeris variables; Astro::Time contains time and unit conversions, and Astro::Misc contains various calculations unrelated to position and time.
The 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 http://www.starlink.rl.ac.uk/Software/software_store.htm.
Thomas R. Wyant, III (wyant at cpan dot org)
Copyright 2005, 2006 by Thomas R. Wyant, III (wyant at cpan dot org). All rights reserved.
This module is free software; you can use it, redistribute it and/or modify it under the same terms as Perl itself. Please see http://perldoc.perl.org/index-licence.html for the current licenses.
This software is provided without any warranty of any kind, express or implied. The author will not be liable for any damages of any sort relating in any way to this software.
|
Astro::Coord::ECI - Manipulate geocentric coordinates |