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