lib/Math/Trig.pm

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