[Summary view] [Print] [Text view]

   1  #
   2  # Trigonometric functions, mostly inherited from Math::Complex.
   3  # -- Jarkko Hietaniemi, since April 1997
   4  # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
   5  #
   6  
   7  require Exporter;
   8  package Math::Trig;
   9  
  10  use 5.005;
  11  use strict;
  12  
  13  use Math::Complex 1.36;
  14  use Math::Complex qw(:trig :pi);
  15  
  16  use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
  17  
  18  @ISA = qw(Exporter);
  19  
  20  $VERSION = 1.04;
  21  
  22  my @angcnv = qw(rad2deg rad2grad
  23          deg2rad deg2grad
  24          grad2rad grad2deg);
  25  
  26  @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
  27         @angcnv);
  28  
  29  my @rdlcnv = qw(cartesian_to_cylindrical
  30          cartesian_to_spherical
  31          cylindrical_to_cartesian
  32          cylindrical_to_spherical
  33          spherical_to_cartesian
  34          spherical_to_cylindrical);
  35  
  36  my @greatcircle = qw(
  37               great_circle_distance
  38               great_circle_direction
  39               great_circle_bearing
  40               great_circle_waypoint
  41               great_circle_midpoint
  42               great_circle_destination
  43              );
  44  
  45  my @pi = qw(pi pi2 pi4 pip2 pip4);
  46  
  47  @EXPORT_OK = (@rdlcnv, @greatcircle, @pi);
  48  
  49  # See e.g. the following pages:
  50  # http://www.movable-type.co.uk/scripts/LatLong.html
  51  # http://williams.best.vwh.net/avform.htm
  52  
  53  %EXPORT_TAGS = ('radial' => [ @rdlcnv ],
  54              'great_circle' => [ @greatcircle ],
  55              'pi'     => [ @pi ]);
  56  
  57  sub _DR  () { pi2/360 }
  58  sub _RD  () { 360/pi2 }
  59  sub _DG  () { 400/360 }
  60  sub _GD  () { 360/400 }
  61  sub _RG  () { 400/pi2 }
  62  sub _GR  () { pi2/400 }
  63  
  64  #
  65  # Truncating remainder.
  66  #
  67  
  68  sub _remt ($$) {
  69      # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
  70      $_[0] - $_[1] * int($_[0] / $_[1]);
  71  }
  72  
  73  #
  74  # Angle conversions.
  75  #
  76  
  77  sub rad2rad($)     { _remt($_[0], pi2) }
  78  
  79  sub deg2deg($)     { _remt($_[0], 360) }
  80  
  81  sub grad2grad($)   { _remt($_[0], 400) }
  82  
  83  sub rad2deg ($;$)  { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) }
  84  
  85  sub deg2rad ($;$)  { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) }
  86  
  87  sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) }
  88  
  89  sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) }
  90  
  91  sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) }
  92  
  93  sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) }
  94  
  95  sub cartesian_to_spherical {
  96      my ( $x, $y, $z ) = @_;
  97  
  98      my $rho = sqrt( $x * $x + $y * $y + $z * $z );
  99  
 100      return ( $rho,
 101               atan2( $y, $x ),
 102               $rho ? acos( $z / $rho ) : 0 );
 103  }
 104  
 105  sub spherical_to_cartesian {
 106      my ( $rho, $theta, $phi ) = @_;
 107  
 108      return ( $rho * cos( $theta ) * sin( $phi ),
 109               $rho * sin( $theta ) * sin( $phi ),
 110               $rho * cos( $phi   ) );
 111  }
 112  
 113  sub spherical_to_cylindrical {
 114      my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
 115  
 116      return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
 117  }
 118  
 119  sub cartesian_to_cylindrical {
 120      my ( $x, $y, $z ) = @_;
 121  
 122      return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
 123  }
 124  
 125  sub cylindrical_to_cartesian {
 126      my ( $rho, $theta, $z ) = @_;
 127  
 128      return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
 129  }
 130  
 131  sub cylindrical_to_spherical {
 132      return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
 133  }
 134  
 135  sub great_circle_distance {
 136      my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
 137  
 138      $rho = 1 unless defined $rho; # Default to the unit sphere.
 139  
 140      my $lat0 = pip2 - $phi0;
 141      my $lat1 = pip2 - $phi1;
 142  
 143      return $rho *
 144          acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
 145               sin( $lat0 ) * sin( $lat1 ) );
 146  }
 147  
 148  sub great_circle_direction {
 149      my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
 150  
 151      my $distance = &great_circle_distance;
 152  
 153      my $lat0 = pip2 - $phi0;
 154      my $lat1 = pip2 - $phi1;
 155  
 156      my $direction =
 157      acos((sin($lat1) - sin($lat0) * cos($distance)) /
 158           (cos($lat0) * sin($distance)));
 159  
 160      $direction = pi2 - $direction
 161      if sin($theta1 - $theta0) < 0;
 162  
 163      return rad2rad($direction);
 164  }
 165  
 166  *great_circle_bearing = \&great_circle_direction;
 167  
 168  sub great_circle_waypoint {
 169      my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
 170  
 171      $point = 0.5 unless defined $point;
 172  
 173      my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
 174  
 175      return undef if $d == pi;
 176  
 177      my $sd = sin($d);
 178  
 179      return ($theta0, $phi0) if $sd == 0;
 180  
 181      my $A = sin((1 - $point) * $d) / $sd;
 182      my $B = sin(     $point  * $d) / $sd;
 183  
 184      my $lat0 = pip2 - $phi0;
 185      my $lat1 = pip2 - $phi1;
 186  
 187      my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
 188      my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
 189      my $z = $A * sin($lat0)                + $B * sin($lat1);
 190  
 191      my $theta = atan2($y, $x);
 192      my $phi   = acos($z);
 193      
 194      return ($theta, $phi);
 195  }
 196  
 197  sub great_circle_midpoint {
 198      great_circle_waypoint(@_[0..3], 0.5);
 199  }
 200  
 201  sub great_circle_destination {
 202      my ( $theta0, $phi0, $dir0, $dst ) = @_;
 203  
 204      my $lat0 = pip2 - $phi0;
 205  
 206      my $phi1   = asin(sin($lat0)*cos($dst)+cos($lat0)*sin($dst)*cos($dir0));
 207      my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
 208                   cos($dst)-sin($lat0)*sin($phi1));
 209  
 210      my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
 211  
 212      $dir1 -= pi2 if $dir1 > pi2;
 213  
 214      return ($theta1, $phi1, $dir1);
 215  }
 216  
 217  1;
 218  
 219  __END__
 220  =pod
 221  
 222  =head1 NAME
 223  
 224  Math::Trig - trigonometric functions
 225  
 226  =head1 SYNOPSIS
 227  
 228      use Math::Trig;
 229  
 230      $x = tan(0.9);
 231      $y = acos(3.7);
 232      $z = asin(2.4);
 233  
 234      $halfpi = pi/2;
 235  
 236      $rad = deg2rad(120);
 237  
 238      # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
 239      use Math::Trig ':pi';
 240  
 241      # Import the conversions between cartesian/spherical/cylindrical.
 242      use Math::Trig ':radial';
 243  
 244          # Import the great circle formulas.
 245      use Math::Trig ':great_circle';
 246  
 247  =head1 DESCRIPTION
 248  
 249  C<Math::Trig> defines many trigonometric functions not defined by the
 250  core Perl which defines only the C<sin()> and C<cos()>.  The constant
 251  B<pi> is also defined as are a few convenience functions for angle
 252  conversions, and I<great circle formulas> for spherical movement.
 253  
 254  =head1 TRIGONOMETRIC FUNCTIONS
 255  
 256  The tangent
 257  
 258  =over 4
 259  
 260  =item B<tan>
 261  
 262  =back
 263  
 264  The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
 265  are aliases)
 266  
 267  B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
 268  
 269  The arcus (also known as the inverse) functions of the sine, cosine,
 270  and tangent
 271  
 272  B<asin>, B<acos>, B<atan>
 273  
 274  The principal value of the arc tangent of y/x
 275  
 276  B<atan2>(y, x)
 277  
 278  The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
 279  and acotan/acot are aliases).  Note that atan2(0, 0) is not well-defined.
 280  
 281  B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
 282  
 283  The hyperbolic sine, cosine, and tangent
 284  
 285  B<sinh>, B<cosh>, B<tanh>
 286  
 287  The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
 288  and cotanh/coth are aliases)
 289  
 290  B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
 291  
 292  The arcus (also known as the inverse) functions of the hyperbolic
 293  sine, cosine, and tangent
 294  
 295  B<asinh>, B<acosh>, B<atanh>
 296  
 297  The arcus cofunctions of the hyperbolic sine, cosine, and tangent
 298  (acsch/acosech and acoth/acotanh are aliases)
 299  
 300  B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
 301  
 302  The trigonometric constant B<pi> and some of handy multiples
 303  of it are also defined.
 304  
 305  B<pi, pi2, pi4, pip2, pip4>
 306  
 307  =head2 ERRORS DUE TO DIVISION BY ZERO
 308  
 309  The following functions
 310  
 311      acoth
 312      acsc
 313      acsch
 314      asec
 315      asech
 316      atanh
 317      cot
 318      coth
 319      csc
 320      csch
 321      sec
 322      sech
 323      tan
 324      tanh
 325  
 326  cannot be computed for all arguments because that would mean dividing
 327  by zero or taking logarithm of zero. These situations cause fatal
 328  runtime errors looking like this
 329  
 330      cot(0): Division by zero.
 331      (Because in the definition of cot(0), the divisor sin(0) is 0)
 332      Died at ...
 333  
 334  or
 335  
 336      atanh(-1): Logarithm of zero.
 337      Died at...
 338  
 339  For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
 340  C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
 341  C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
 342  C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
 343  C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
 344  pi>, where I<k> is any integer.
 345  
 346  Note that atan2(0, 0) is not well-defined.
 347  
 348  =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
 349  
 350  Please note that some of the trigonometric functions can break out
 351  from the B<real axis> into the B<complex plane>. For example
 352  C<asin(2)> has no definition for plain real numbers but it has
 353  definition for complex numbers.
 354  
 355  In Perl terms this means that supplying the usual Perl numbers (also
 356  known as scalars, please see L<perldata>) as input for the
 357  trigonometric functions might produce as output results that no more
 358  are simple real numbers: instead they are complex numbers.
 359  
 360  The C<Math::Trig> handles this by using the C<Math::Complex> package
 361  which knows how to handle complex numbers, please see L<Math::Complex>
 362  for more information. In practice you need not to worry about getting
 363  complex numbers as results because the C<Math::Complex> takes care of
 364  details like for example how to display complex numbers. For example:
 365  
 366      print asin(2), "\n";
 367  
 368  should produce something like this (take or leave few last decimals):
 369  
 370      1.5707963267949-1.31695789692482i
 371  
 372  That is, a complex number with the real part of approximately C<1.571>
 373  and the imaginary part of approximately C<-1.317>.
 374  
 375  =head1 PLANE ANGLE CONVERSIONS
 376  
 377  (Plane, 2-dimensional) angles may be converted with the following functions.
 378  
 379  =over
 380  
 381  =item deg2rad
 382  
 383      $radians  = deg2rad($degrees);
 384  
 385  =item grad2rad
 386  
 387      $radians  = grad2rad($gradians);
 388  
 389  =item rad2deg
 390  
 391      $degrees  = rad2deg($radians);
 392  
 393  =item grad2deg
 394  
 395      $degrees  = grad2deg($gradians);
 396  
 397  =item deg2grad
 398  
 399      $gradians = deg2grad($degrees);
 400  
 401  =item rad2grad
 402  
 403      $gradians = rad2grad($radians);
 404  
 405  =back
 406  
 407  The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
 408  The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
 409  If you don't want this, supply a true second argument:
 410  
 411      $zillions_of_radians  = deg2rad($zillions_of_degrees, 1);
 412      $negative_degrees     = rad2deg($negative_radians, 1);
 413  
 414  You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
 415  grad2grad().
 416  
 417  =over 4
 418  
 419  =item rad2rad
 420  
 421      $radians_wrapped_by_2pi = rad2rad($radians);
 422  
 423  =item deg2deg
 424  
 425      $degrees_wrapped_by_360 = deg2deg($degrees);
 426  
 427  =item grad2grad
 428  
 429      $gradians_wrapped_by_400 = grad2grad($gradians);
 430  
 431  =back
 432  
 433  =head1 RADIAL COORDINATE CONVERSIONS
 434  
 435  B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
 436  systems, explained shortly in more detail.
 437  
 438  You can import radial coordinate conversion functions by using the
 439  C<:radial> tag:
 440  
 441      use Math::Trig ':radial';
 442  
 443      ($rho, $theta, $z)     = cartesian_to_cylindrical($x, $y, $z);
 444      ($rho, $theta, $phi)   = cartesian_to_spherical($x, $y, $z);
 445      ($x, $y, $z)           = cylindrical_to_cartesian($rho, $theta, $z);
 446      ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
 447      ($x, $y, $z)           = spherical_to_cartesian($rho, $theta, $phi);
 448      ($rho_c, $theta, $z)   = spherical_to_cylindrical($rho_s, $theta, $phi);
 449  
 450  B<All angles are in radians>.
 451  
 452  =head2 COORDINATE SYSTEMS
 453  
 454  B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
 455  
 456  Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
 457  coordinates which define a point in three-dimensional space.  They are
 458  based on a sphere surface.  The radius of the sphere is B<rho>, also
 459  known as the I<radial> coordinate.  The angle in the I<xy>-plane
 460  (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
 461  coordinate.  The angle from the I<z>-axis is B<phi>, also known as the
 462  I<polar> coordinate.  The North Pole is therefore I<0, 0, rho>, and
 463  the Gulf of Guinea (think of the missing big chunk of Africa) I<0,
 464  pi/2, rho>.  In geographical terms I<phi> is latitude (northward
 465  positive, southward negative) and I<theta> is longitude (eastward
 466  positive, westward negative).
 467  
 468  B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
 469  some texts define the I<phi> to start from the horizontal plane, some
 470  texts use I<r> in place of I<rho>.
 471  
 472  Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
 473  coordinates which define a point in three-dimensional space.  They are
 474  based on a cylinder surface.  The radius of the cylinder is B<rho>,
 475  also known as the I<radial> coordinate.  The angle in the I<xy>-plane
 476  (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
 477  coordinate.  The third coordinate is the I<z>, pointing up from the
 478  B<theta>-plane.
 479  
 480  =head2 3-D ANGLE CONVERSIONS
 481  
 482  Conversions to and from spherical and cylindrical coordinates are
 483  available.  Please notice that the conversions are not necessarily
 484  reversible because of the equalities like I<pi> angles being equal to
 485  I<-pi> angles.
 486  
 487  =over 4
 488  
 489  =item cartesian_to_cylindrical
 490  
 491      ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
 492  
 493  =item cartesian_to_spherical
 494  
 495      ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
 496  
 497  =item cylindrical_to_cartesian
 498  
 499      ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
 500  
 501  =item cylindrical_to_spherical
 502  
 503      ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
 504  
 505  Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
 506  
 507  =item spherical_to_cartesian
 508  
 509      ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
 510  
 511  =item spherical_to_cylindrical
 512  
 513      ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
 514  
 515  Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
 516  
 517  =back
 518  
 519  =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
 520  
 521  A great circle is section of a circle that contains the circle
 522  diameter: the shortest distance between two (non-antipodal) points on
 523  the spherical surface goes along the great circle connecting those two
 524  points.
 525  
 526  =head2 great_circle_distance
 527  
 528  You can compute spherical distances, called B<great circle distances>,
 529  by importing the great_circle_distance() function:
 530  
 531    use Math::Trig 'great_circle_distance';
 532  
 533    $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
 534  
 535  The I<great circle distance> is the shortest distance between two
 536  points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
 537  optional, it defaults to 1 (the unit sphere), therefore the distance
 538  defaults to radians.
 539  
 540  If you think geographically the I<theta> are longitudes: zero at the
 541  Greenwhich meridian, eastward positive, westward negative--and the
 542  I<phi> are latitudes: zero at the North Pole, northward positive,
 543  southward negative.  B<NOTE>: this formula thinks in mathematics, not
 544  geographically: the I<phi> zero is at the North Pole, not at the
 545  Equator on the west coast of Africa (Bay of Guinea).  You need to
 546  subtract your geographical coordinates from I<pi/2> (also known as 90
 547  degrees).
 548  
 549    $distance = great_circle_distance($lon0, pi/2 - $lat0,
 550                                      $lon1, pi/2 - $lat1, $rho);
 551  
 552  =head2 great_circle_direction
 553  
 554  The direction you must follow the great circle (also known as I<bearing>)
 555  can be computed by the great_circle_direction() function:
 556  
 557    use Math::Trig 'great_circle_direction';
 558  
 559    $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
 560  
 561  =head2 great_circle_bearing
 562  
 563  Alias 'great_circle_bearing' for 'great_circle_direction' is also available.
 564  
 565    use Math::Trig 'great_circle_bearing';
 566  
 567    $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1);
 568  
 569  The result of great_circle_direction is in radians, zero indicating
 570  straight north, pi or -pi straight south, pi/2 straight west, and
 571  -pi/2 straight east.
 572  
 573  You can inversely compute the destination if you know the
 574  starting point, direction, and distance:
 575  
 576  =head2 great_circle_destination
 577  
 578    use Math::Trig 'great_circle_destination';
 579  
 580    # thetad and phid are the destination coordinates,
 581    # dird is the final direction at the destination.
 582  
 583    ($thetad, $phid, $dird) =
 584      great_circle_destination($theta, $phi, $direction, $distance);
 585  
 586  or the midpoint if you know the end points:
 587  
 588  =head2 great_circle_midpoint
 589  
 590    use Math::Trig 'great_circle_midpoint';
 591  
 592    ($thetam, $phim) =
 593      great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
 594  
 595  The great_circle_midpoint() is just a special case of
 596  
 597  =head2 great_circle_waypoint
 598  
 599    use Math::Trig 'great_circle_waypoint';
 600  
 601    ($thetai, $phii) =
 602      great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
 603  
 604  Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
 605  $phi1).  Note that antipodal points (where their distance is I<pi>
 606  radians) do not have waypoints between them (they would have an an
 607  "equator" between them), and therefore C<undef> is returned for
 608  antipodal points.  If the points are the same and the distance
 609  therefore zero and all waypoints therefore identical, the first point
 610  (either point) is returned.
 611  
 612  The thetas, phis, direction, and distance in the above are all in radians.
 613  
 614  You can import all the great circle formulas by
 615  
 616    use Math::Trig ':great_circle';
 617  
 618  Notice that the resulting directions might be somewhat surprising if
 619  you are looking at a flat worldmap: in such map projections the great
 620  circles quite often do not look like the shortest routes-- but for
 621  example the shortest possible routes from Europe or North America to
 622  Asia do often cross the polar regions.
 623  
 624  =head1 EXAMPLES
 625  
 626  To calculate the distance between London (51.3N 0.5W) and Tokyo
 627  (35.7N 139.8E) in kilometers:
 628  
 629      use Math::Trig qw(great_circle_distance deg2rad);
 630  
 631      # Notice the 90 - latitude: phi zero is at the North Pole.
 632      sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
 633      my @L = NESW( -0.5, 51.3);
 634      my @T = NESW(139.8, 35.7);
 635      my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
 636  
 637  The direction you would have to go from London to Tokyo (in radians,
 638  straight north being zero, straight east being pi/2).
 639  
 640      use Math::Trig qw(great_circle_direction);
 641  
 642      my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
 643  
 644  The midpoint between London and Tokyo being
 645  
 646      use Math::Trig qw(great_circle_midpoint);
 647  
 648      my @M = great_circle_midpoint(@L, @T);
 649  
 650  or about 89.16N 68.93E, practically at the North Pole.
 651  
 652  =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
 653  
 654  The answers may be off by few percentages because of the irregular
 655  (slightly aspherical) form of the Earth.  The errors are at worst
 656  about 0.55%, but generally below 0.3%.
 657  
 658  =head1 BUGS
 659  
 660  Saying C<use Math::Trig;> exports many mathematical routines in the
 661  caller environment and even overrides some (C<sin>, C<cos>).  This is
 662  construed as a feature by the Authors, actually... ;-)
 663  
 664  The code is not optimized for speed, especially because we use
 665  C<Math::Complex> and thus go quite near complex numbers while doing
 666  the computations even when the arguments are not. This, however,
 667  cannot be completely avoided if we want things like C<asin(2)> to give
 668  an answer instead of giving a fatal runtime error.
 669  
 670  Do not attempt navigation using these formulas.
 671  
 672  =head1 AUTHORS
 673  
 674  Jarkko Hietaniemi <F<jhi!at!iki.fi>> and 
 675  Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>.
 676  
 677  =cut
 678  
 679  # eof