t/eci.t - metacpan.org

package main;

use strict;
use warnings;

use lib qw{ inc };

use Astro::Coord::ECI;
use Astro::Coord::ECI::Utils qw{ :time deg2rad PERL2000 rad2deg };
use My::Module::Test qw{ :tolerance velocity_sanity };
use POSIX qw{strftime floor};
use Test::More 0.88;

use constant EQUATORIALRADIUS => 6378.14;	# Meeus page 82.
use constant TIMFMT => '%d-%b-%Y %H:%M:%S';

Astro::Coord::ECI->set (debug => 0);

# universal time
# Tests: universal()

# We just make sure we get the same thing back.

{
    my $want = timegm ( 0, 0, 0, 1, 0, 100 );
    my $got = Astro::Coord::ECI->universal( $want )->universal();

    cmp_ok $got, '==', $want,
	'Univeral time round trip: Jan 1 2000';

    $want = timegm ( 0, 0, 0, 1, 0, 105 );
    $got = Astro::Coord::ECI->universal( $want )->universal();

    cmp_ok $got, '==', $want,
	'Univeral time round trip: Jan 1 2005';
}


# universal time -> dynamical time
# Tests: dynamical()

{
    my $univ = timegm( 0, 0, 0, 1, 0, 100 );
    my $dyn = floor(
	Astro::Coord::ECI->universal( $univ )->dynamical + .5 );

    cmp_ok $dyn, '==', $univ + 65,
	'Universal to dynamical time: Jan 1 2000';

    $univ = timegm( 0, 0, 0, 1, 0, 105 );
    $dyn = floor(
	Astro::Coord::ECI->universal( $univ )->dynamical + .5 );

    cmp_ok $dyn, '==', $univ + 72,
	'Universal to dynamical time: Jan 1 2005';
}


# dynamical time -> universal time
# tests: dynamical()

{
    my $dyn = timegm( 0, 0, 0, 1, 0, 100 );
    my $univ = floor(
	Astro::Coord::ECI->dynamical( $dyn
	)->universal() + .5 );

    cmp_ok $univ, '==', $dyn - 65,
	'Dynamical to universal time: Jan 1 2000 dynamical';

    $dyn = timegm( 0, 0, 0, 1, 0, 105 );
    $univ = floor(
	Astro::Coord::ECI->dynamical( $dyn
	)->universal() + .5 );

    cmp_ok $univ, '==', $dyn - 72,
	'Dynamical to universal time: Jan 1 2005 dynamical';
}


# ecef
# Tests: ecef()

# All we do here is be sure we get back what we put in.

{
    my ( $X, $Y, $Z ) = Astro::Coord::ECI->ecef( 3000, 4000, 5000 )->ecef();

    cmp_ok $X, '==', 3000, 'ECEF round-trip: X';

    cmp_ok $Y, '==', 4000, 'ECEF round-trip: Y';

    cmp_ok $Z, '==', 5000, 'ECEF round-trip: Z';
}


# geodetic -> geocentric
# Tests: geodetic()

# Meeus, page 82, example 11a

# Both TopoZone and Google say the observatory is latitude 34 deg 13'33"
# N (=   degrees), longitude 118 deg 03'25" W (= -118.056944444444
# degrees).  The test uses Meeus' latitude of 33 deg 21'22" N (since
# that's what Meeus himself uses) but the TopoZone/Google value for
# longitude, since longitude does not affect the calculation, but my
# Procrustean validation expects it.

# We also test the antpodal (sort of) point, since finding a bug in my
# implementation of Borkowski's algorithm when running on a point in the
# southern hemisphere. No, this particular test doesn't use that
# algorithm, but once bitten, twice shy.

{
    my ( $phiprime, undef, $rho ) =
	Astro::Coord::ECI->new( ellipsoid => 'IAU76' )->
	geodetic( .58217396455, -2.060487233536, 1.706 )->geocentric();
    my $rhosinphiprime = $rho / EQUATORIALRADIUS * sin ($phiprime);
    my $rhocosphiprime = $rho / EQUATORIALRADIUS * cos ($phiprime);

    cmp_ok sprintf( '%.6f', $rhosinphiprime ), '==', .546861,
	'geodetic to geocentric: rho * sin( phiprime )';

    cmp_ok sprintf( '%.6f', $rhocosphiprime ), '==', .836339,
	'geodetic to geocentric: rho * cos( phiprime )';

    ( $phiprime, undef, $rho ) =
	Astro::Coord::ECI->new( ellipsoid => 'IAU76' )->
	geodetic( -.58217396455, 2.060487233536, 1.706 )->geocentric();
    $rhosinphiprime = $rho / EQUATORIALRADIUS * sin ($phiprime);
    $rhocosphiprime = $rho / EQUATORIALRADIUS * cos ($phiprime);

    cmp_ok sprintf( '%.6f', $rhosinphiprime ), '==', -.546861,
	'geodetic to geocentric: rho * sin( phiprime )';

    cmp_ok sprintf( '%.6f', $rhocosphiprime ), '==', .836339,
	'geodetic to geocentric: rho * cos( phiprime )';
}


# geocentric -> geodetic
# Tests: geodetic()

# Borkowski
# For this, we just invert Meeus' example.

# We also test the antpodal point, since finding a bug in my
# implementation of Borkowski's algorithm when running on a point in the
# southern hemisphere.

{
    my ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'IAU76' )->
	geocentric(
	    0.579094339305825, -2.060487233536, 6373.41803380646,
	)->geodetic();

    tolerance_frac $lat, .58217396455, 1e-6,
	'Geocentric to geodetic: latitude';

    tolerance_frac $long, -2.060487233536, 1e-6,
	'Geocentric to geodetic: longitude';

    tolerance_frac $elev, 1.706, 1e-3,
	'Geocentric to geodetic: height above sea level';

#	[IAU76 => -0.579094339305825, 1.08110542005979, 6373.41803380646,
#		[-.58217396455, 1e-6], [1.08110542005979, 1e-6], [1.706, 1e-3]],
#	) {

    ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'IAU76' )->
	geocentric(
	    -0.579094339305825, 1.08110542005979, 6373.41803380646,
	)->geodetic();

    tolerance_frac $lat, -.58217396455, 1e-6,
	'Geocentric to geodetic: latitude';

    tolerance_frac $long, 1.08110542005979, 1e-6,
	'Geocentric to geodetic: longitude';

    tolerance_frac $elev, 1.706, 1e-3,
	'Geocentric to geodetic: height above sea level';
}

# geodetic -> Earth-Centered, Earth-Fixed
# Tests: geocentric() (and geodetic())

# Continuing the above example, but ecef coordinates. Book answer from
# http://www.ngs.noaa.gov/cgi-bin/xyz_getxyz.prl is
#    x                y             z       Elipsoid
# -2508975.4549 -4707403.8939  3487953.2711 GRS80

note <<'EOD';

In the following twelve tests the tolerance is degraded because the book
solution is calculated using a different, and apparently simpler model
attributed to Escobal, "Methods of Orbit Determination", 1965, Wiley &
Sons, Inc., pp. 27-29.

EOD

{
    my ( $x, $y, $z ) =
	Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	geodetic( .58217396455, -2.060487233536, 1.706 )->ecef();

    tolerance_frac $x, -2508.9754549, 1e-5,
	'Geodetic to ECEF: X';

    tolerance_frac $y, -4707.4038939, 1e-5,
	'Geodetic to ECEF: Y';

    tolerance_frac $z, 3487.9532711, 1e-5,
	'Geodetic to ECEF: Z';

    ( $x, $y, $z ) =
	Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	geodetic( -.58217396455, 1.08110542005979, 1.706 )->ecef();

    tolerance_frac $x, 2508.9754549, 1e-5,
	'Geodetic to ECEF: X';

    tolerance_frac $y, 4707.4038939, 1e-5,
	'Geodetic to ECEF: Y';

    tolerance_frac $z, -3487.9532711, 1e-5,
	'Geodetic to ECEF: Z';
}


# Earth-Centered, Earth-Fixed -> geodetic
# Tests: geocentric() (and geodetic())

# Continuing the above example, but ecef coordinates. We use the book
# solution of the opposite test as our input, and vice versa.

{
    my ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	ecef( -2508.9754549, -4707.4038939, 3487.9532711 )->geodetic();

    tolerance_frac $lat, .58217396455, 1e-5,
	'ECEF to Geodetic: latitude';

    tolerance_frac $long, -2.060487233536, 1e-5,
	'ECEF to Geodetic: longitude';

    tolerance_frac $elev, 1.706, 1e-5,
	'ECEF to Geodetic: height above sea level';

    ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	ecef( 2508.9754549, 4707.4038939, -3487.9532711 )->geodetic();

    tolerance_frac $lat, -.58217396455, 1e-5,
	'ECEF to Geodetic: latitude';

    tolerance_frac $long, 1.08110542005979, 1e-5,
	'ECEF to Geodetic: longitude';

    tolerance_frac $elev, 1.706, 1e-5,
	'ECEF to Geodetic: height above sea level';
}


# geodetic -> eci
# Tests: eci() (and geodetic() and geocentric())

# Standard is from http://celestrak.com/columns/v02n03/ (Kelso)

{
    my $time = timegm( 0, 0, 9, 1, 9, 95 );

    my ( $x, $y, $z ) =
	Astro::Coord::ECI->new( ellipsoid => 'WGS72' )->
	geodetic( deg2rad( 40 ), deg2rad( -75 ), 0 )->eci( $time );

    tolerance_frac $x, 1703.295, 1e-6, 'Geodetic to ECI: X';

    tolerance_frac $y, 4586.650, 1e-6, 'Geodetic to ECI: Y';

    tolerance_frac $z, 4077.984, 1e-6, 'Geodetic to ECI: Z';

    $time = timegm( 0, 0, 9, 1, 9, 95 );

    ( $x, $y, $z ) =
	Astro::Coord::ECI->new( ellipsoid => 'WGS72' )->
	geodetic( deg2rad( -40 ), deg2rad( 105 ), 0 )->eci( $time );

    tolerance_frac $x, -1703.295, 1e-6, 'Geodetic to ECI: X';

    tolerance_frac $y, -4586.650, 1e-6, 'Geodetic to ECI: Y';

    tolerance_frac $z, -4077.984, 1e-6, 'Geodetic to ECI: Z';
}


# eci -> geodetic
# Tests: eci() (and geodetic() and geocentric())

# This is the reverse of the previous test.

{
    my $time = timegm( 0, 0, 9, 1, 9, 95 );

    my ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'WGS72' )->
	eci( 1703.295, 4586.650, 4077.984, $time )->geodetic();

    tolerance_frac $lat, deg2rad( 40 ), 1e-6,
	'ECI to geodetic: latitude';

    tolerance_frac $long, deg2rad( -75 ), 1e-6,
	'ECI to geodetic: longitude';

    $elev += EQUATORIALRADIUS;
    tolerance_frac $elev, EQUATORIALRADIUS, 1e-6,
	'ECI to geodetic: distance from center';

    $time = timegm (0, 0, 9, 1, 9, 95);

    ( $lat, $long, $elev ) =
	Astro::Coord::ECI->new( ellipsoid => 'WGS72' )->
	eci( -1703.295, -4586.650, -4077.984, $time )->geodetic();

    tolerance_frac $lat, deg2rad( -40 ), 1e-6,
	'ECI to geodetic: latitude';

    tolerance_frac $long, deg2rad( 105 ), 1e-6,
	'ECI to geodetic: longitude';

    $elev += EQUATORIALRADIUS;
    tolerance_frac $elev, EQUATORIALRADIUS, 1e-6,
	'ECI to geodetic: distance from center';
}


# azel
# Tests: azel() (and geodetic(), geocentric(), and eci())

# Book solution from
# http://www.satcom.co.uk/article.asp?article=1

note <<'EOD';

In the following three tests the tolerance is degraded because the
book solution is calculated by http://www.satcom.co.uk/article.asp?article=1
which apparently assumes an exactly synchronous orbit. Their exact
altitude assumption is undocumented, as is their algorithm. So the
tests are really more of a sanity check.

EOD

{
    my $time = timegm (0, 0, 5, 27, 7, 105);
    my $sta = Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	geodetic( deg2rad( 38 ), deg2rad( -80 ), 1 );
    my $sat = Astro::Coord::ECI->new( ellipsoid => 'GRS80' )->
	universal( $time )->
	geodetic( deg2rad( 0 ), deg2rad( -75 ), 35800 );

    my ( $azm, $elev, $rng ) = $sta->azel( $sat );

    tolerance_frac $azm, deg2rad( 171.906 ), 1e-3,
	'Azimuth for observer';

    tolerance_frac $elev, deg2rad( 45.682 ), 1e-3,
	'Elevation for observer';

    tolerance_frac $rng, 37355.457, 1e-3,
	'Range for observer';

    # Same as above tests, but with the station in the 'station'
    # attribute.

    $sat->set( station => $sta );

    ( $azm, $elev, $rng ) = $sat->azel();

    tolerance_frac $azm, deg2rad( 171.906 ), 1e-3,
	'Azimuth for observer';

    tolerance_frac $elev, deg2rad( 45.682 ), 1e-3,
	'Elevation for observer';

    tolerance_frac $rng, 37355.457, 1e-3,
	'Range for observer';

    # enu
    # Tests: enu()

    # From the azel() test above, converted as described at
    # http://geostarslib.sourceforge.net/main.html#conv

    my ( $East, $North, $Up ) = $sat->enu();

    tolerance $East, 3675, 10,
	'East for observer';

    tolerance $North, -25840, 10,
	'North for observer';

    tolerance $Up, 26746, 10,
	'Up for observer';
}

# atmospheric refraction.
# Tests: correct_for_refraction()

# Based on Meeus' Example 16.a.

{
    my $got = Astro::Coord::ECI->
	correct_for_refraction( deg2rad( .5541 ) );

    tolerance_frac $got, deg2rad( 57.864 / 60 ), 1e-4,
	'Correction for atmospheric refraction';
}


# Angle between two points as seen from a third.
# Tests: angle()

{
    my $A = Astro::Coord::ECI->ecef( 0, 0, 0 );
    my $B = Astro::Coord::ECI->ecef( 1, 0, 0 );
    my $C = Astro::Coord::ECI->ecef( 0, 1, 0 );

    my $got = $A->angle ($B, $C);

    tolerance $got, deg2rad( 90 ), 1e-6,
	'Angle between two points as seen from a third';
}


# Precession of equinoxes.
# Tests: precess()

# Based on Meeus' example 21.b.

use constant LIGHTYEAR2KILOMETER => 9.4607e12;

{
    note q{Precession - no 'station' attribute};

    my $alpha0 = 41.054063;
    my $delta0 = 49.227750;
    my $rho = 36.64;
    my $t0 = PERL2000;
    my $alphae = deg2rad( 41.547214 );
    my $deltae = deg2rad( 49.348483 );
    my $time = timegm( 0, 0, 0, 13, 10, 128 ) + .19 * 86400;

    my $eci = Astro::Coord::ECI->dynamical( $t0 )->equatorial(
	deg2rad( $alpha0 ), deg2rad( $delta0 ),
	$rho *  LIGHTYEAR2KILOMETER )->set( equinox_dynamical => $t0 );
    my $utim = Astro::Coord::ECI->dynamical( $time )->universal();
    my ( $alpha, $delta ) = $eci->precess( $utim )->equatorial();
    my $tolerance = 1e-6;

    tolerance $alpha, $alphae, $tolerance,
	'Precession of equinoxes: right ascension';

    tolerance $delta, $deltae, $tolerance,
	'Precession of equinoxes: declination';
}

{
    note q{Precession - with 'station' attribute};

    my $alpha0 = 41.054063;
    my $delta0 = 49.227750;
    my $rho = 36.64;
    my $t0 = PERL2000;
    my $alphae = deg2rad( 41.547214 );
    my $deltae = deg2rad( 49.348483 );
    my $time = timegm( 0, 0, 0, 13, 10, 128 ) + .19 * 86400;

    my $sta = Astro::Coord::ECI->dynamical( $t0 )->eci( 0, 0, 0 )->set(
	equinox_dynamical => $t0 );
    my $eci = Astro::Coord::ECI->dynamical( $t0 )->equatorial(
	deg2rad( $alpha0 ), deg2rad( $delta0 ),
	$rho *  LIGHTYEAR2KILOMETER )->set(
	equinox_dynamical => $t0,
	station	=> $sta,
    );
    my $utim = Astro::Coord::ECI->dynamical( $time )->universal();
    my ( $alpha, $delta ) = $eci->precess( $utim )->equatorial();
    my $tolerance = 1e-6;

    tolerance $alpha, $alphae, $tolerance,
	'Precession of equinoxes: right ascension';

    tolerance $delta, $deltae, $tolerance,
	'Precession of equinoxes: declination';

    cmp_ok $sta->get( 'equinox_dynamical' ), '==', $time,
	'Station object was precessed';
}


# Right ascension/declination to ecliptic lat/lon
# Tests: ecliptic() (and obliquity())

# Based on Meeus' example 13.a, page 95.

# Meeus' example involves the star Pollux. We use an arbitrary (and much
# too small) rho, because it doesn't come into the conversion anyway.
# The time matters because it figures in to the obliquity of the
# ecliptic. Unfortunately Meeus didn't give us the time in his example,
# only the obliquity. The time used in the example was chosen because it
# gave the desired obliquity value of 23.4392911 degrees.

{
    my $time = timegm( 36, 27, 2, 30, 6, 109 );

    my ( $lat, $long ) = Astro::Coord::ECI->equatorial(
	deg2rad( 116.328942 ), deg2rad( 28.026183 ), 1e12, $time )->ecliptic();

    tolerance_frac $lat, deg2rad( 6.684170 ), 1e-6,
	'Equatorial to Ecliptic: latitude';


    tolerance_frac $long, deg2rad( 113.215630 ), 1e-6,
	'Equatorial to Ecliptic: longitude';
}


# Ecliptic lat/lon to right ascension/declination
# Tests: ecliptic() (and obliquity())

# Based on inverting the above test.

{
    my $time = timegm( 36, 27, 2, 30, 6, 109 );

    my ( $ra, $dec ) = Astro::Coord::ECI->ecliptic(
	deg2rad( 6.684170 ), deg2rad( 113.215630 ), 1e12, $time )->equatorial();

    tolerance_frac $ra, deg2rad( 116.328942 ), 1e-6,
	'Ecliptic to Equatorial: Right ascension';

    tolerance_frac $dec, deg2rad( 28.026183 ), 1e-6,
	'Ecliptic to Equatorial: Declination';
}

use constant ASTRONOMICAL_UNIT => 149_597_870; # Meeus, Appendix 1, pg 407

# Ecliptic lat/long to ECI
# Tests: equatorial() (and ecliptic())

# This test is based on Meeus' example 26.a.

{
    my $time = timegm( 0, 0, 0, 13, 9, 92 );
    my $lat = .62 / 3600;
    my $lon = 199.907347;
    my $rho = .99760775 * ASTRONOMICAL_UNIT;
    my $expx = -0.9379952 * ASTRONOMICAL_UNIT;
    my $expy = -0.3116544 * ASTRONOMICAL_UNIT;
    my $expz = -0.1351215 * ASTRONOMICAL_UNIT;

    my ( $x, $y, $z ) = Astro::Coord::ECI->dynamical( $time )->
	ecliptic( deg2rad( $lat ), deg2rad( $lon ), $rho
	)->eci();
    my $tolerance = 1e-5;

    tolerance_frac $x, $expx, $tolerance, 'Ecliptic to ECI: X';

    tolerance_frac $y, $expy, $tolerance, 'Ecliptic to ECI: Y';

    tolerance_frac $z, $expz, $tolerance, 'Ecliptic to ECI: Z';
}

# universal time to local mean time
# Tests: local_mean_time()

# This test is based on http://www.statoids.com/tconcept.html

{
    my $time = timegm( 0, 0, 0, 1, 0, 101 );
    my $lat = 29/60 + 40;
    my $lon = -( 8/60 + 86 );
    my $offset = -( ( 5 * 60 + 44 ) * 60 + 32 );

    my $got = floor( Astro::Coord::ECI->geodetic(
	    deg2rad( $lat ), deg2rad( $lon ), 0 )->universal(
	    $time )->local_mean_time() );
    my $want = $time + $offset;

    cmp_ok $got, '==', $want, 'Universal time to Local mean time';
}


# local mean time to universal time
# Tests: local_mean_time()

# This test is the inverse of the previous one.

{
    my $time = timegm( 28, 15, 18, 31, 11, 100 );
    my $lat = 29/60 + 40;
    my $lon = -( 8/60 + 86 );
    my $offset = -( ( 5 * 60 + 44 ) * 60 + 32 );

    my $got = floor( Astro::Coord::ECI->geodetic(
	    deg2rad( $lat ), deg2rad( $lon ), 0 )->local_mean_time(
	    $time )->universal() + .5 );
    my $want = $time - $offset;

    cmp_ok $got, '==', $want, 'Local mean time to universal time';
}


# Equatorial coordinates relative to observer.

# I don't have a book solution for this, but if I turn off atmospheric
# refraction, you should get the same result as if you simply subtract
# the ECI coordinates of the observer from the ECI coordinates of the
# body.


{	# Begin local symbol block;

    note 'Equatorial relative to location, specified explicitly';

    my $time = time ();
    my $station = Astro::Coord::ECI->geodetic(
	deg2rad (  38.898741 ),
	deg2rad ( -77.037684 ),
	0.01668
    );
    $station->set( refraction => 0 );
    my $body = Astro::Coord::ECI->eci(
	2328.97048951, -5995.22076416, 1719.97067261,
	2.91207230, -0.98341546, -7.09081703, $time);
    my @staeci = $station->eci( $time );
    my @bdyeci = $body->eci();
    my @want = Astro::Coord::ECI->eci(
	    ( map {$bdyeci[$_] - $staeci[$_]} 0 .. 5 ), $time )->
	    equatorial();
    my @got = $station->equatorial ($body);

    tolerance $got[0], $want[0], 0.000001,
	'Right ascension relative to explicit location';

    tolerance $got[1], $want[1], 0.000001,
	'Declination relative to explicit location';

    tolerance $got[2], $want[2], 0.000001,
	'Right ascension relative to explicit location';

}	# End local symbol block.

{	# Begin local symbol block;

    note 'Equatorial relative to location specified in station attribute';

    my $time = time ();
    my $station = Astro::Coord::ECI->geodetic(
	deg2rad (  38.898741 ),
	deg2rad ( -77.037684 ),
	0.01668
    );
    $station->set( refraction => 0 );
    my $body = Astro::Coord::ECI->eci(
	2328.97048951, -5995.22076416, 1719.97067261,
	2.91207230, -0.98341546, -7.09081703, $time);
    my @staeci = $station->eci( $time );
    my @bdyeci = $body->eci();
    my @want = Astro::Coord::ECI->eci(
	    ( map {$bdyeci[$_] - $staeci[$_]} 0 .. 5 ), $time )->
	    equatorial();

    $body->set( station => $station );
    my @got = $body->equatorial_apparent();

    tolerance $got[0], $want[0], 0.000001,
	'Right ascension relative to station attribute';

    tolerance $got[1], $want[1], 0.000001,
	'Declination relative to station attribute';

    tolerance $got[2], $want[2], 0.000001,
	'Right ascension relative to station attribute';

}	# End local symbol block.

# represents().

is( Astro::Coord::ECI->represents(), 'Astro::Coord::ECI',
    'Astro::Coord::ECI->represents() returns itself' );

ok( Astro::Coord::ECI->represents( 'Astro::Coord::ECI' ),
    q{Astro::Coord::ECI->represents('Astro::Coord::ECI') is true} );

ok( ! Astro::Coord::ECI->represents( 'Astro::Coord::ECI::TLE' ),
    q{Astro::Coord::ECI->represents('Astro::Coord::ECI::TLE') is false} );

{
    # Maidenhead Locator System. Reference implementation is at
    # http://home.arcor.de/waldemar.kebsch/The_Makrothen_Contest/fmaidenhead.html
    my $got;
    my ( $lat, $lon ) = ( 38.898748, -77.037684 );

    my $sta = Astro::Coord::ECI->new()->geodetic(
	deg2rad( $lat ),
	deg2rad( $lon ),
	0,
    );


    ( $got ) = $sta->maidenhead( 3 );
    is $got, 'FM18lv', "Maidenhead precision 3 for $lat, $lon is 'FM18lv'";

    ( $got ) = $sta->maidenhead( 2 );
    is $got, 'FM18', "Maidenhead precision 2 for $lat, $lon is 'FM18'";

    ( $got ) = $sta->maidenhead( 1 );
    is $got, 'FM', "Maidenhead precision 1 for $lat, $lon is 'FM'";

}

# Velocity sanity tests
{
    my $time = timegm( 0, 0, 12, 1, 3, 112 );

    my $body = Astro::Coord::ECI->new(
	name => 'eci coordinates',
    )->eci(
	1000, 1000, 1000, 0, 0, 0, $time );

    my $sta = Astro::Coord::ECI->new(
	name => 'White House',
    )->geodetic(
	deg2rad( 38.8987 ),
	deg2rad( -77.0377 ),
	17 / 1000,
    );

    $body->set( station => $sta );

    velocity_sanity ecef => $body;

    velocity_sanity neu => $body->universal( $time );

    velocity_sanity enu => $body->universal( $time );

    velocity_sanity azel => $body->universal( $time ), $sta;

    velocity_sanity azel => $body->universal( $time );

    velocity_sanity eci => $sta->universal( $time );

    {
	$body->set( frequency => 1_000_000 );
	my @azel = $body->azel();
	if ( @azel > 6 ) {
	    # This is just to be sure it is calculated consistently; I
	    # have no idea whether I am getting the right answer.
	    tolerance $azel[6], .25, 0.01, 'Doppler shift';
	} else {
	    fail 'No Doppler calculated';
	}
	$body->set( frequency => undef );
    }

    # I would love to have tested the internal routines before this, but
    # the critical thing is velocity, and I have no test data. So
    # assuming that the sanity tests are correct, we go back the other
    # way and see how we do.

    my @want = $body->universal( $time )->neu();
    my @got = Astro::Coord::ECI::_convert_spherical_to_cartesian(
	$body->azel() );
    my @coord = qw{ X Y Z X_dot Y_dot Z_dot };
    foreach my $inx ( 0 .. 5 ) {
	tolerance $got[$inx], $want[$inx], 1e-12, $coord[$inx];
    }
}

done_testing;

# need to test:
#    dip
#    get (class, object)
#    reference_ellipsoid
#    set (class, object with and without resetting the object)

1;

# ex: set textwidth=72 :
	Global
`s`	Focus search bar
`?`	Bring up this help dialog
	GitHub
`g` `p`	Go to pull requests
`g` `i`	go to github issues (only if github is preferred repository)
	POD
`g` `a`	Go to author
`g` `c`	Go to changes
`g` `i`	Go to issues
`g` `d`	Go to dist
`g` `r`	Go to repository/SCM
`g` `s`	Go to source
`g` `b`	Go to file browse
	Search terms
module: (e.g. module:Plugin)
distribution: (e.g. distribution:Dancer auth)
author: (e.g. author:SONGMU Redis)
version: (e.g. version:1.00)