Geo::Ellipsoid - Longitude and latitude calculations using an ellipsoid model.

  use Geo::Ellipsoid;
  $geo = Geo::Ellipsoid->new(ellipsoid=>'NAD27', units=>'degrees');
  @origin = ( 37.619002, -122.374843 );    # SFO
  @dest = ( 33.942536, -118.408074 );      # LAX
  ( $range, $bearing ) = $geo->to( @origin, @dest );
  ($lat,$lon) = $geo->at( @origin, 2000, 45.0 );
  ( $x, $y ) = $geo->displacement( @origin, $lat, $lon );
  @pos = $geo->location( $lat, $lon, $x, $y );

    f = ( semi-major - semi-minor ) / semi-major

    my $geo = Geo::Ellipsoid->new( 
      ellipsoid => 'WGS84', 
      units => 'radians' ,
      distance_units => 'meter',
      longitude => 0,
      bearing => 0,
    );

    $geo->set_units('degrees');

    $geo->set_distance_unit('kilometer');

    $geo->set_distance_unit(201.168);

  Unit          Units per meter
  --------      ---------------
  foot             0.3048
  kilometer     1000.0
  mile          1609.344
  nm            1852.0

    $geo->set_ellipsoid('NAD27');

    $geo->set_custom_ellipsoid( 'sphere', 6378137, 0 );

    $geo->set_longitude_symmetric(1);

    $geo->set_bearing_symmetric(1);

    $Geo::Ellipsoid->set_defaults( 
      units => 'degrees', 
      ellipsoid => 'GRS80',
      distance_units => 'kilometer',
      longitude => 1,
      bearing => 0
    );

    $Geo::Ellipsoid->set_defaults( 
      uni => 'deg', 
      ell => 'GRS80',
      dis => 'kil',
      lon => 1,
      bea => 0
    );

    ( $lat_scale, $lon_scale ) = $geo->scales($lat0);
    $x = $lon_scale * ($lon - $lon0); 
    $y = $lat_scale * ($lat - $lat0); 

    my $dist = $geo->range( $lat1, $lon1, $lat2, $lon2 );
    my $dist = $geo->range( @origin, @destination );

    my $bearing = $geo->bearing( $lat1, $lon1, $lat2, $lon2 );

    my( $lat2, $lon2 ) = $geo->at( $lat1, $lon1, $range, $bearing );

    my( $dist, $theta ) = $geo->to( $lat1, $lon1, $lat2, $lon2 );
    my $dist = $geo->to( $lat1, $lon1, $lat2, $lon2 );

    my( $x, $y ) = $geo->displacement( $lat1, $lon1, $lat2, $lon2 );

	my @loc = $geo->location( $lat, $lon, $x, $y );

    Ellipsoid        Semi-Major Axis (m.)     1/Flattening
    ---------        -------------------     ---------------
    AIRY                 6377563.396         299.3249646
    AIRY-MODIFIED        6377340.189         299.3249646
    AUSTRALIAN           6378160.0           298.25
    BESSEL-1841          6377397.155         299.1528128
    CLARKE-1880          6378249.145         293.465
    EVEREST-1830         6377276.345         290.8017
    EVEREST-MODIFIED     6377304.063         290.8017
    FISHER-1960          6378166.0           298.3
    FISHER-1968          6378150.0           298.3
    GRS80                6378137.0           298.25722210088
    HOUGH-1956           6378270.0           297.0
    HAYFORD              6378388.0           297.0
    IAU76                6378140.0           298.257
    KRASSOVSKY-1938      6378245.0           298.3
    NAD27                6378206.4           294.9786982138
    NWL-9D               6378145.0           298.25
    SOUTHAMERICAN-1969   6378160.0           298.25
    SOVIET-1985          6378136.0           298.257
    WGS72                6378135.0           298.26
    WGS84                6378137.0           298.257223563

    ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/

#	Geo::Ellipsoid
#
#	This package implements an Ellipsoid class to perform latitude 
#	and longitude calculations on the surface of an ellipsoid.
#
#	This is a Perl conversion of existing Fortran code (see
#	ACKNOWLEDGEMENTS) and the author of this class makes no 
#	claims of originality. Nor can he even vouch for the 
#	results of the calculations, although they do seem to 
#	work for him and have been tested against other methods.

package Geo::Ellipsoid;

use warnings;
use strict;
use 5.006_00;

use Scalar::Util 'looks_like_number';
use Math::Trig;
use Carp;

our $VERSION = '1.12';
our $DEBUG = 0;

# class data and constants
our $degrees_per_radian = 180/pi;
our $eps = 1.0e-23;
our $max_loop_count = 20;
our $twopi = 2 * pi;
our $halfpi = pi/2;
our %defaults = ( 
  ellipsoid => 'WGS84', 
  units => 'radians',
  distance_units => 'meter', 
  longitude => 0,
  latitude => 1,	# allows use of _normalize_output
  bearing => 0,
);
our %distance = (
  'foot'      => 0.3048,
  'kilometer' => 1_000,
  'meter'     => 1.0,
  'mile'      => 1_609.344,
  'nm'        => 1_852,
);

# set of ellipsoids that can be used.
# values are 
#  1) a = semi-major (equatorial) radius of Ellipsoid 
#  2) 1/f = reciprocal of flattening (f), the ratio of the semi-minor 
#     (polar) radius to the semi-major (equatorial) axis, or
#     polar radius = equatorial radius * ( 1 - f )

our %ellipsoids = (
    'AIRY'               => [ 6377563.396, 299.3249646     ],
    'AIRY-MODIFIED'      => [ 6377340.189, 299.3249646     ],
    'AUSTRALIAN'         => [ 6378160.0,   298.25          ],
    'BESSEL-1841'        => [ 6377397.155, 299.1528128     ],
    'CLARKE-1880'        => [ 6378249.145, 293.465         ],
    'EVEREST-1830'       => [ 6377276.345, 300.8017        ],
    'EVEREST-MODIFIED'   => [ 6377304.063, 300.8017        ],
    'FISHER-1960'        => [ 6378166.0,   298.3           ],
    'FISHER-1968'        => [ 6378150.0,   298.3           ],
    'GRS80'              => [ 6378137.0,   298.25722210088 ],
    'HOUGH-1956'         => [ 6378270.0,   297.0           ],
    'HAYFORD'            => [ 6378388.0,   297.0           ],
    'IAU76'              => [ 6378140.0,   298.257         ],
    'KRASSOVSKY-1938'    => [ 6378245.0,   298.3           ],
    'NAD27'              => [ 6378206.4,   294.9786982138  ],
    'NWL-9D'             => [ 6378145.0,   298.25          ],
    'SOUTHAMERICAN-1969' => [ 6378160.0,   298.25          ],
    'SOVIET-1985'        => [ 6378136.0,   298.257         ],
    'WGS72'              => [ 6378135.0,   298.26          ],
    'WGS84'              => [ 6378137.0,   298.257223563   ],
);

sub new
{
  my( $class, %args ) = @_;
  my $self = {%defaults};
  print "new: @_\n" if $DEBUG;
  foreach my $key ( keys %args ) {
    my $val = $args{$key};
    if( $key =~ /^ell/i ) {
      $self->{ellipsoid} = uc $args{$key};
    }elsif( $key =~ /^uni/i ) {
      $self->{units} = $args{$key};
    }elsif( $key =~ /^dis/i ) {
      $self->{distance_units} = $args{$key};
    }elsif( $key =~ /^lon/i ) {
      $self->{longitude} = $args{$key};
    }elsif( $key =~ /^bea/i ) {
      $self->{bearing} = $args{$key};
    }else{
      carp("Unknown argument: $key => $args{$key}");
    }
  }
  set_units($self,$self->{units});
  set_ellipsoid($self,$self->{ellipsoid});
  set_distance_unit($self,$self->{distance_units});
  set_longitude_symmetric($self,$self->{longitude});
  set_bearing_symmetric($self,$self->{bearing});
  print 
    "Ellipsoid(units=>$self->{units},distance_units=>" .
    "$self->{distance_units},ellipsoid=>$self->{ellipsoid}," .
    "longitude=>$self->{longitude},bearing=>$self->{bearing})\n" if $DEBUG;
  bless $self,$class;
  return $self;
}

sub set_units
{
  my $self = shift;
  my $units = shift;
  if( $units =~ /deg/i ) {
    $units = 'degrees';
  }elsif( $units =~ /rad/i ) {
    $units = 'radians';
  }else{
    croak("Invalid units specifier '$units' - please use either " .
      "degrees or radians (the default)") unless $units =~ /rad/i;
  }
  $self->{units} = $units;
}

sub set_distance_unit
{
  my $self = shift;
  my $unit = shift;
  print "distance unit = $unit\n" if $DEBUG;

  my $conversion = 0;

  if( defined $unit ) {
  
    my( $key, $val );
    while( ($key,$val) = each %distance ) {
      my $re = substr($key,0,3);
      print "trying ($key,$re,$val)\n" if $DEBUG;
      if( $unit =~ /^$re/i ) {
        $self->{distance_units} = $unit;
        $conversion = $val;	

	# finish iterating to reset 'each' function call
	while( each %distance ) {}
	last;
      }
    }

    if( $conversion == 0 ) {
      if( looks_like_number($unit) ) {
        $conversion = $unit;
      }else{
        carp("Unknown argument to set_distance_unit: $unit\nAssuming meters");
        $conversion = 1.0;
      }
    }
  }else{
    carp("Missing or undefined argument to set_distance_unit: ".
      "$unit\nAssuming meters");
    $conversion = 1.0;
  }
  $self->{conversion} = $conversion;
}

sub set_ellipsoid
{
  my $self = shift;
  my $ellipsoid = uc shift || $defaults{ellipsoid};
  print "  set ellipsoid to $ellipsoid\n" if $DEBUG;
  unless( exists $ellipsoids{$ellipsoid} ) {
    croak("Ellipsoid $ellipsoid does not exist - please use " .
      "set_custom_ellipsoid to use an ellipsoid not in valid set");
  }
  $self->{ellipsoid} = $ellipsoid;
  my( $major, $recip ) = @{$ellipsoids{$ellipsoid}};
  $self->{equatorial} = $major;
  if( $recip == 0 ) {
    carp("Infinite flattening specified by ellipsoid -- assuming a sphere");
    $self->{polar} = $self->{equatorial};
    $self->{flattening} = 0;
    $self->{eccentricity} = 0;
  }else{
    $self->{flattening} = ( 1.0 / $ellipsoids{$ellipsoid}[1]);
    $self->{polar} = $self->{equatorial} * ( 1.0  - $self->{flattening} );
    $self->{eccentricity} = sqrt( 2.0 * $self->{flattening} - 
      ( $self->{flattening} * $self->{flattening} ) );
  }
}

sub set_custom_ellipsoid
{
  my $self = shift;
  my( $name, $major, $recip ) = @_;
  $name = uc $name;
  $recip = 0 unless defined $recip;
  if( $major ) {
    $ellipsoids{$name} = [ $major, $recip ];
  }else{  
    croak("set_custom_ellipsoid called without semi-major radius parameter");
  }
  set_ellipsoid($self,$name);
}

sub set_longitude_symmetric
{
  my( $self, $sym ) = @_;
  # see if argument passed
  if( $#_ > 0 ) {
    # yes -- use value passed
    $self->{longitude} = $sym;
  }else{
    # no -- set to true
    $self->{longitude} = 1;
  }
}

sub set_bearing_symmetric
{
  my( $self, $sym ) = @_;
  # see if argument passed
  if( $#_ > 0 ) {
    # yes -- use value passed
    $self->{bearing} = $sym;
  }else{
    # no -- set to true
    $self->{bearing} = 1;
  }
}

sub set_defaults
{
  my $self = shift;
  my %args = @_;
  foreach my $key ( keys %args ) {
    if( $key =~ /^ell/i ) {
      $defaults{ellipsoid} = uc $args{$key};
    }elsif( $key =~ /^uni/i ) {
      $defaults{units} = $args{$key};
    }elsif( $key =~ /^dis/i ) {
      $defaults{distance_units} = $args{$key};
    }elsif( $key =~ /^lon/i ) {
      $defaults{longitude} = $args{$key};
    }elsif( $key =~ /^bea/i ) {
      $defaults{bearing} = $args{$key};
    }else{  
      croak("Geo::Ellipsoid::set_defaults called with invalid key: $key");
    }
  }
  print "Defaults set to ($defaults{ellipsoid},$defaults{units}\n"
    if $DEBUG;
}

sub scales
{
  my $self = shift;
  my $units = $self->{units};
  my $lat = $_[0];
  if( defined $lat ) {
    $lat /= $degrees_per_radian if( $units eq 'degrees' );
  }else{
    carp("scales() method requires latitude argument; assuming 0");
    $lat = 0;
  }
  
  my $aa = $self->{equatorial};
  my $bb = $self->{polar};
  my $a2 = $aa*$aa;
  my $b2 = $bb*$bb;
  my $d1 = $aa * cos($lat);
  my $d2 = $bb * sin($lat);
  my $d3 = $d1*$d1 + $d2*$d2;
  my $d4 = sqrt($d3);
  my $n1 = $aa * $bb;
  my $latscl = ( $n1 * $n1 ) / ( $d3 * $d4 * $self->{conversion} );
  my $lonscl = ( $aa * $d1 ) / ( $d4 * $self->{conversion} );
  
  if( $DEBUG ) {
    print "lat=$lat, aa=$aa, bb=$bb\nd1=$d1, d2=$d2, d3=$d3, d4=$d4\n";
    print "latscl=$latscl, lonscl=$lonscl\n";
  }

  if( $self->{units} eq 'degrees' ) {
    $latscl /= $degrees_per_radian;
    $lonscl /= $degrees_per_radian;
  }
  return ( $latscl, $lonscl );
}

    
sub range
{
  my $self = shift;
  my @args = _normalize_input($self->{units},@_);
  my($range,$bearing) = _inverse($self,@args);
  print "inverse(@_[1..4]) returns($range,$bearing)\n" if $DEBUG;
  return $range;
}

sub bearing
{
  my $self = shift;
  my $units = $self->{units};
  my @args = _normalize_input($units,@_);
  my($range,$bearing) = _inverse($self,@args);
  print "inverse(@args) returns($range,$bearing)\n" if $DEBUG;
  my $t = $bearing;
  $self->_normalize_output('bearing',$bearing);
  print "_normalize_output($t) returns($bearing)\n" if $DEBUG;
  return $bearing;
}

sub at
{
  my $self = shift;
  my $units = $self->{units};
  my( $lat, $lon, $az ) = _normalize_input($units,@_[0,1,3]);
  my $r = $_[2];
  print "at($lat,$lon,$r,$az)\n" if $DEBUG;
  my( $lat2, $lon2 ) = _forward($self,$lat,$lon,$r,$az);
  print "_forward returns ($lat2,$lon2)\n" if $DEBUG; 
  $self->_normalize_output('longitude',$lon2);
  $self->_normalize_output('latitude',$lat2);
  return ( $lat2, $lon2 );  
}

sub to
{
  my $self = shift;
  my $units = $self->{units};
  my @args = _normalize_input($units,@_);
  print "to($units,@args)\n" if $DEBUG;
  my($range,$bearing) = _inverse($self,@args);
  print "to: inverse(@args) returns($range,$bearing)\n" if $DEBUG;
  #$bearing *= $degrees_per_radian if $units eq 'degrees';
  $self->_normalize_output('bearing',$bearing);
  if( wantarray() ) {
    return ( $range, $bearing );
  }else{
    return $range;
  }
}

sub displacement
{
  my $self = shift;
  print "displacement(",join(',',@_),"\n" if $DEBUG;
  my @args = _normalize_input($self->{units},@_);
  print "call _inverse(@args)\n" if $DEBUG;
  my( $range, $bearing ) = _inverse($self,@args);
  print "disp: _inverse(@args) returns ($range,$bearing)\n" if $DEBUG;
  my $x = $range * sin($bearing);
  my $y = $range * cos($bearing);
  return ($x,$y);
}

sub location
{
  my $self = shift;
  my $units = $self->{units};
  my($lat,$lon,$x,$y) = @_;
  my $range = sqrt( $x*$x+ $y*$y );
  my $bearing = atan2($x,$y);
  $bearing *= $degrees_per_radian if $units eq 'degrees';
  print "location($lat,$lon,$x,$y,$range,$bearing)\n" if $DEBUG;
  return $self->at($lat,$lon,$range,$bearing);
}

########################################################################
#
#      internal functions
#
#	inverse
#
#	Calculate the displacement from origin to destination.
#	The input to this subroutine is 
#	  ( latitude-1, longitude-1, latitude-2, longitude-2 ) in radians.
#
#	Return the results as the list (range,bearing) with range in the
#	current specified distance unit and bearing in radians.

sub _inverse()
{
  my $self = shift;
  my( $lat1, $lon1, $lat2, $lon2 ) = (@_);
  print "_inverse($lat1,$lon1,$lat2,$lon2)\n" if $DEBUG;
  
  my $a = $self->{equatorial};
  my $f = $self->{flattening};
  
  my $r = 1.0 - $f;
  my $tu1 = $r * sin($lat1) / cos($lat1);
  my $tu2 = $r * sin($lat2) / cos($lat2);
  my $cu1 = 1.0 / ( sqrt(($tu1*$tu1) + 1.0) );
  my $su1 = $cu1 * $tu1;
  my $cu2 = 1.0 / ( sqrt( ($tu2*$tu2) + 1.0 ));
  my $s = $cu1 * $cu2;
  my $baz = $s * $tu2;
  my $faz = $baz * $tu1;
  my $dlon = $lon2 - $lon1;
  
  if( $DEBUG ) {
    printf "lat1=%.8f, lon1=%.8f\n", $lat1, $lon1; 
    printf "lat2=%.8f, lon2=%.8f\n", $lat2, $lon2;
    printf "r=%.8f, tu1=%.8f, tu2=%.8f\n", $r, $tu1, $tu2;
    printf "faz=%.8f, dlon=%.8f\n", $faz, $dlon;
  }
  
  my $x = $dlon;
  my $cnt = 0;
  print "enter loop:\n" if $DEBUG;
  my( $c2a, $c, $cx, $cy, $cz, $d, $del, $e, $sx, $sy, $y );
  do {
    printf "  x=%.8f\n", $x if $DEBUG;
    $sx = sin($x);
    $cx = cos($x);
    $tu1 = $cu2*$sx;
    $tu2 = $baz - ($su1*$cu2*$cx);
  
    printf "    sx=%.8f, cx=%.8f, tu1=%.8f, tu2=%.8f\n", 
      $sx, $cx, $tu1, $tu2 if $DEBUG;
  
    $sy = sqrt( $tu1*$tu1 + $tu2*$tu2 );
    $cy = $s*$cx + $faz;
    $y = atan2($sy,$cy);
    my $sa;
    if( $sy == 0.0 ) {
      $sa = 1.0;
    }else{
      $sa = ($s*$sx) / $sy;
    }
  
    printf "    sy=%.8f, cy=%.8f, y=%.8f, sa=%.8f\n", $sy, $cy, $y, $sa
      if $DEBUG;
  
    $c2a = 1.0 - ($sa*$sa);
    $cz = $faz + $faz;
    if( $c2a > 0.0 ) {
      $cz = ((-$cz)/$c2a) + $cy;
    }
    $e = ( 2.0 * $cz * $cz ) - 1.0;
    $c = ( ((( (-3.0 * $c2a) + 4.0)*$f) + 4.0) * $c2a * $f )/16.0;
    $d = $x;
    $x = ( ($e * $cy * $c + $cz) * $sy * $c + $y) * $sa;
    $x = ( 1.0 - $c ) * $x * $f + $dlon;
    $del = $d - $x;
  
    if( $DEBUG ) {
      printf "    c2a=%.8f, cz=%.8f\n", $c2a, $cz;
      printf "    e=%.8f, d=%.8f\n", $e, $d;
      printf "    (d-x)=%.8g\n", $del;
    }
  
  }while( (abs($del) > $eps) && ( ++$cnt <= $max_loop_count ) );
  
  $faz = atan2($tu1,$tu2);
  $baz = atan2($cu1*$sx,($baz*$cx - $su1*$cu2)) + pi;
  $x = sqrt( ((1.0/($r*$r)) -1.0 ) * $c2a+1.0 ) + 1.0;
  $x = ($x-2.0)/$x;
  $c = 1.0 - $x;
  $c = (($x*$x)/4.0 + 1.0)/$c;
  $d = ((0.375*$x*$x) - 1.0)*$x;
  $x = $e*$cy;
  
  if( $DEBUG ) {
    printf "e=%.8f, cy=%.8f, x=%.8f\n", $e, $cy, $x;
    printf "sy=%.8f, c=%.8f, d=%.8f\n", $sy, $c, $d;
    printf "cz=%.8f, a=%.8f, r=%.8f\n", $cz, $a, $r;
  }
  
  $s = 1.0 - $e - $e;
  $s = (((((((( $sy * $sy * 4.0 ) - 3.0) * $s * $cz * $d/6.0) - $x) * 
    $d /4.0) + $cz) * $sy * $d) + $y ) * $c * $a * $r;
  
  printf "s=%.8f\n", $s if $DEBUG;
  
  # adjust azimuth to (0,360) or (-180,180) as specified
  if( $self->{symmetric} ) {
    $faz += $twopi if $faz < -(pi);
    $faz -= $twopi if $faz >= pi;
  }else{
    $faz += $twopi if $faz < 0;
    $faz -= $twopi if $faz >= $twopi;
  }

  # return result
  my @disp = ( ($s/$self->{conversion}), $faz );
  print "disp = (@disp)\n" if $DEBUG;
  return @disp;
}

#	_forward
#
#	Calculate the location (latitue,longitude) of a point
#	given a starting point and a displacement from that
#	point as (range,bearing)
#
sub _forward
{
  my $self = shift;
  my( $lat1, $lon1, $range, $bearing ) = @_;

  if( $DEBUG ) {
    printf "_forward(lat1=%.8f,lon1=%.8f,range=%.8f,bearing=%.8f)\n",
      $lat1, $lon1, $range, $bearing;
  }

  my $eps = 0.5e-13;

  my $a = $self->{equatorial};
  my $f = $self->{flattening};
  my $r = 1.0 - $f;

  my $tu = $r * sin($lat1) / cos($lat1);
  my $faz = $bearing;
  my $s = $self->{conversion} * $range;
  my $sf = sin($faz);
  my $cf = cos($faz);

  my $baz = 0.0;
  $baz = 2.0 * atan2($tu,$cf) if( $cf != 0.0 );
  
  my $cu = 1.0 / sqrt(1.0 + $tu*$tu);
  my $su = $tu * $cu;
  my $sa = $cu * $sf;
  my $c2a = 1.0 - ($sa*$sa);
  my $x = 1.0 + sqrt( (((1.0/($r*$r)) - 1.0 )*$c2a) +1.0);
  $x = ($x-2.0)/$x;
  my $c = 1.0 - $x;
  $c = ((($x*$x)/4.0) + 1.0)/$c;
  my $d = $x * ((0.375*$x*$x)-1.0);
  $tu = (($s/$r)/$a)/$c;
  my $y = $tu;

  if( $DEBUG ) {
    printf "r=%.8f, tu=%.8f, faz=%.8f\n", $r, $tu, $faz;
    printf "baz=%.8f, sf=%.8f, cf=%.8f\n", $baz, $sf, $cf;
    printf "cu=%.8f, su=%.8f, sa=%.8f\n", $cu, $su, $sa; 
    printf "x=%.8f, c=%.8f, y=%.8f\n", $x, $c, $y;
  }

  my( $cy, $cz, $e, $sy );
  do {
    $sy = sin($y);
    $cy = cos($y);
    $cz = cos($baz+$y);
    $e = (2.0*$cz*$cz)-1.0;
    $c = $y;
    $x = $e * $cy;
    $y = (2.0 * $e) - 1.0;
    $y = ((((((((($sy*$sy*4.0)-3.0)*$y*$cz*$d)/6.0)+$x)*$d)/4.0)-$cz)*$sy*$d) +
      $tu;
    } while( abs($y-$c) > $eps );

  $baz = ($cu*$cy*$cf) - ($su*$sy);
  $c = $r*sqrt(($sa*$sa) + ($baz*$baz));
  $d = $su*$cy + $cu*$sy*$cf;
  my $lat2 = atan2($d,$c);
  $c = $cu*$cy - $su*$sy*$cf;
  $x = atan2($sy*$sf,$c);
  $c = (((((-3.0*$c2a)+4.0)*$f)+4.0)*$c2a*$f)/16.0;
  $d = (((($e*$cy*$c) + $cz)*$sy*$c)+$y)*$sa;
  my $lon2 = $lon1 + $x - (1.0-$c)*$d*$f;
  #$baz = atan2($sa,$baz) + pi;

  # return result
  return ($lat2,$lon2);

}

#	_normalize_input
#
#	Normalize a set of input angle values by converting to
#	radians if given in degrees and by converting to the
#	range [0,2pi), i.e. greater than or equal to zero and
#	less than two pi.
#
sub _normalize_input
{
  my $units = shift;
  my @args = @_;
  return map {
    $_ = deg2rad($_) if $units eq 'degrees';
    while( $_ < 0 ) { $_ += $twopi }
    while( $_ >= $twopi ) { $_ -= $twopi }
    $_
  } @args;
}  

#	_normalize_output
#
#	Normalize a set of output angle values by converting to
#	degrees if needed and by converting to the range [-pi,+pi) or
#	[0,2pi) as needed.
#
sub _normalize_output
{
  my $self = shift;
  my $elem = shift;	# 'bearing' or 'longitude'
  # adjust remaining input values by reference
  for ( @_ ) {
    if( $self->{$elem} ) {
      # normalize to range [-pi,pi)
      while( $_ < -(pi) ) { $_ += $twopi }
      while( $_ >= pi ) { $_ -= $twopi }
    }else{
      # normalize to range [0,2*pi)
      while( $_ < 0 ) { $_ += $twopi }
      while( $_ >= $twopi ) { $_ -= $twopi }
    }
    $_ = rad2deg($_) if $self->{units} eq 'degrees';
  }
}  

1; # End of Geo::Ellipsoid