view release on metacpan or  search on metacpan

lib/Astro/Coord/ECI.pm  view on Meta::CPAN

    # 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);

lib/Astro/Coord/ECI.pm  view on Meta::CPAN

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 sun (Astro::Coord::ECI::Sun class or object)

This attribute represents the Sun. This is typically used to determine
ambient lighting conditions and other Sun-related phenomena. Subclasses
may place more stringent restrictions on this attribute.

=item twilight (numeric, radians)

This attribute represents the elevation of the center of the Sun's disk
at the beginning and end of twilight. This must be a number between
-PI/2 and PI/2 inclusive, or C<undef> which sets the value to the
default.

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 SUBCLASSING

This class is and will remain based on a hash reference.

It is recommended that subclasses store their attributes in a separate
hash hung off the main object in a key named after the subclass. That
is, something like

 sub subclass_method {
     my ( $self ) = @_;
     my $attr = $self->{ ref $self };
     ...
 }

From outside the subclass, these attributes can be made accessible by
overriding C<get()> and C<set()>.

Subclasses B<must> initialize themselves by overriding the C<new()>
method.  The override B<must> call C<< $class->SUPER::new() >>, passing
it any arguments that were passed to the override.

=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