view release on metacpan or  search on metacpan

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

	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.

In list context, it returns the time, and an indicator which is true if
the body passes above the given elevation at that time, and false if it
passes below the given elevation. In scalar context it just returns the
time.

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>

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

Use this method with caution. I have serious doubts about what will
happen when the invocant is north of the Arctic Circle or south of the
Antarctic Circle, and I am a little unsure of its behaviour near the
Equator. In addition, in an early implementation consecutive calls
returned the same extreme over and over. I believe this to be fixed, but
the fix raises the possibility of skipped extrema. B<Caveat user.>

=item ( $time, $type ) = $body->next_elevation_extreme_tod( $elev, $ipper );

This way to call C<next_elevation_extreme_tod()> is equivalent to the
above, but uses the body's C<station> attribute for the C<$coord>
object.

=cut

sub next_elevation_extreme_tod {
    my ( $self, $body, $angle, $upper ) = _expand_args_default_station( @_ );

    ref $self or croak <<'EOD';
Error - The next_elevation_extreme_tod() method may not be called as a
        class method.
EOD

    $body->represents(__PACKAGE__) or croak <<"EOD";
Error - The first argument to next_elevation_extreme_tod() must be a
        subclass of @{[__PACKAGE__]}.
EOD

    return $body->__next_elevation_extreme_tod( $self, $angle, $upper );
}

sub __next_elevation_extreme_tod {
    return;
}

=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 ( $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 -