Astro-satpass
view release on metacpan or search on metacpan
lib/Astro/Coord/ECI.pm view on Meta::CPAN
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 $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
lib/Astro/Coord/ECI.pm view on Meta::CPAN
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;
return $self->local_mean_time() + $self->equation_of_time();
}
=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/ [^0-9] /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/ [^0-9] /smx
and croak 'Invalid Maidenhead locator; ',
"'$chr' is not a digit";
$_ += $chr;
lib/Astro/Coord/ECI.pm view on Meta::CPAN
$body->universal ($end);
$self->universal ($end);
return wantarray ? ($end, $above) : $end;
}
=item ( $delta_psi, $delta_epsilon ) = $self->nutation( $time )
This method calculates the nutation in longitude (delta psi) and
obliquity (delta epsilon) for the given B<dynamical> time. If the time
is unspecified or specified as C<undef>, the current B<dynamical> time
of the object is used.
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 22, pages 143ff. Meeus states that it is good to
0.5 seconds of arc.
=cut
sub nutation {
my ( $self, $time ) = @_;
defined $time
or $time = $self->dynamical();
my $T = jcent2000( $time ); # Meeus (22.1)
my $omega = mod2pi( deg2rad( ( ( $T / 450000 + .0020708 ) * $T -
1934.136261 ) * $T + 125.04452 ) );
my $L = mod2pi( deg2rad( 36000.7698 * $T + 280.4665 ) );
my $Lprime = mod2pi( deg2rad( 481267.8813 * $T + 218.3165 ) );
my $delta_psi = deg2rad( ( -17.20 * sin( $omega )
- 1.32 * sin( 2 * $L )
- 0.23 * sin( 2 * $Lprime )
+ 0.21 * sin( 2 * $omega ) ) / 3600 );
my $delta_epsilon = deg2rad( ( 9.20 * cos( $omega )
+ 0.57 * cos( 2 * $L )
+ 0.10 * cos( 2 * $Lprime )
- 0.09 * cos( 2 * $omega ) ) / 3600 );
return ( $delta_psi, $delta_epsilon );
}
=item $epsilon = $self->obliquity( $time )
This method calculates the obliquity of the ecliptic in radians at the
given B<dynamical> time. If the time is unspecified or specified as
C<undef>, the current B<dynamical> time of the object is used.
The algorithm comes from Jean Meeus' "Astronomical Algorithms", 2nd
Edition, Chapter 22, pages 143ff. The conversion from universal to
dynamical time comes from chapter 10, equation 10.2 on page 78.
=cut
use constant E0BASE => (21.446 / 60 + 26) / 60 + 23;
sub obliquity {
my ( $self, $time ) = @_;
if ( looks_like_number( $self ) ) {
( $self, $time ) = ( __PACKAGE__, $self );
}
defined $time
or $time = $self->dynamical();
my $T = jcent2000 ($time); # Meeus (22.1)
my ( undef, $delta_epsilon ) = $self->nutation( $time );
my $epsilon0 = deg2rad( ( ( 0.001813 * $T - 0.00059 ) * $T - 46.8150 )
* $T / 3600 + E0BASE);
return $epsilon0 + $delta_epsilon;
}
=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
lib/Astro/Coord/ECI.pm view on Meta::CPAN
'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_elevation,
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_elevation,
);
{
# 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 )
&& $value >= - PIOVER2 # Not -PIOVER2 to avoid warning under 5.10.1.
&& $value <= PIOVER2
) {
$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;
}
{
my %dflt = (
twilight => CIVIL_TWILIGHT,
);
# Any attribute specified as angle above or below the horizon.
sub _set_elevation {
my ( $self, $name, $value ) = @_;
defined $value
or $value = ( $dflt{$name} || 0 ); # Default
looks_like_number( $value )
and $value >= - PIOVER2 # Not -PIOVER2 to avoid warning under 5.10.1.
and $value <= PIOVER2
or croak "'$value' is an invalid value for '$name'";
$self->{$name} = $value;
return SET_ACTION_NONE;
}
}
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 ) = @_;
embodies( $value, $self->SUN_CLASS() )
or croak sprintf 'The sun attribute must be a %s',
$self->SUN_CLASS();
ref $value
or $value = $value->new();
$self->{$name} = $value;
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;
lib/Astro/Coord/ECI.pm view on Meta::CPAN
# 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;
}
# 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 );
}
} else {
# $sin_theta = $sin_phi = 0;
# $cos_theta = $cos_phi = 1;
# We used to do the above and then drop through and do the
# velocity calculation if needed. But in this case the
# velocity calulation blows up. So for the moment we are
# just going to punt.
}
}
( run in 0.950 second using v1.01-cache-2.11-cpan-8f98c5d2c55 )