2 # Trigonometric functions, mostly inherited from Math::Complex.
3 # -- Jarkko Hietaniemi, since April 1997
4 # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
13 use Math::Complex 1.51;
14 use Math::Complex qw(:trig :pi);
16 use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
22 my @angcnv = qw(rad2deg rad2grad
26 my @areal = qw(asin_real acos_real);
28 @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
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);
40 great_circle_direction
44 great_circle_destination
47 my @pi = qw(pi pi2 pi4 pip2 pip4);
49 @EXPORT_OK = (@rdlcnv, @greatcircle, @pi, 'Inf');
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
55 %EXPORT_TAGS = ('radial' => [ @rdlcnv ],
56 'great_circle' => [ @greatcircle ],
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 }
67 # Truncating remainder.
71 # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
72 $_[0] - $_[1] * int($_[0] / $_[1]);
79 sub rad2rad($) { _remt($_[0], pi2) }
81 sub deg2deg($) { _remt($_[0], 360) }
83 sub grad2grad($) { _remt($_[0], 400) }
85 sub rad2deg ($;$) { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) }
87 sub deg2rad ($;$) { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) }
89 sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) }
91 sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) }
93 sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) }
95 sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) }
98 # acos and asin functions which always return a real number
102 return 0 if $_[0] >= 1;
103 return pi if $_[0] <= -1;
108 return &pip2 if $_[0] >= 1;
109 return -&pip2 if $_[0] <= -1;
113 sub cartesian_to_spherical {
114 my ( $x, $y, $z ) = @_;
116 my $rho = sqrt( $x * $x + $y * $y + $z * $z );
120 $rho ? acos_real( $z / $rho ) : 0 );
123 sub spherical_to_cartesian {
124 my ( $rho, $theta, $phi ) = @_;
126 return ( $rho * cos( $theta ) * sin( $phi ),
127 $rho * sin( $theta ) * sin( $phi ),
128 $rho * cos( $phi ) );
131 sub spherical_to_cylindrical {
132 my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
134 return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
137 sub cartesian_to_cylindrical {
138 my ( $x, $y, $z ) = @_;
140 return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
143 sub cylindrical_to_cartesian {
144 my ( $rho, $theta, $z ) = @_;
146 return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
149 sub cylindrical_to_spherical {
150 return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
153 sub great_circle_distance {
154 my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
156 $rho = 1 unless defined $rho; # Default to the unit sphere.
158 my $lat0 = pip2 - $phi0;
159 my $lat1 = pip2 - $phi1;
162 acos_real( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
163 sin( $lat0 ) * sin( $lat1 ) );
166 sub great_circle_direction {
167 my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
169 my $distance = &great_circle_distance;
171 my $lat0 = pip2 - $phi0;
172 my $lat1 = pip2 - $phi1;
175 acos_real((sin($lat1) - sin($lat0) * cos($distance)) /
176 (cos($lat0) * sin($distance)));
178 $direction = pi2 - $direction
179 if sin($theta1 - $theta0) < 0;
181 return rad2rad($direction);
184 *great_circle_bearing = \&great_circle_direction;
186 sub great_circle_waypoint {
187 my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
189 $point = 0.5 unless defined $point;
191 my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
193 return undef if $d == pi;
197 return ($theta0, $phi0) if $sd == 0;
199 my $A = sin((1 - $point) * $d) / $sd;
200 my $B = sin( $point * $d) / $sd;
202 my $lat0 = pip2 - $phi0;
203 my $lat1 = pip2 - $phi1;
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);
209 my $theta = atan2($y, $x);
210 my $phi = acos_real($z);
212 return ($theta, $phi);
215 sub great_circle_midpoint {
216 great_circle_waypoint(@_[0..3], 0.5);
219 sub great_circle_destination {
220 my ( $theta0, $phi0, $dir0, $dst ) = @_;
222 my $lat0 = pip2 - $phi0;
224 my $phi1 = asin_real(sin($lat0)*cos($dst) +
225 cos($lat0)*sin($dst)*cos($dir0));
226 my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
227 cos($dst)-sin($lat0)*sin($phi1));
229 my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
231 $dir1 -= pi2 if $dir1 > pi2;
233 return ($theta1, $phi1, $dir1);
243 Math::Trig - trigonometric functions
257 # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
258 use Math::Trig ':pi';
260 # Import the conversions between cartesian/spherical/cylindrical.
261 use Math::Trig ':radial';
263 # Import the great circle formulas.
264 use Math::Trig ':great_circle';
268 C<Math::Trig> defines many trigonometric functions not defined by the
269 core Perl which defines only the C<sin()> and C<cos()>. The constant
270 B<pi> is also defined as are a few convenience functions for angle
271 conversions, and I<great circle formulas> for spherical movement.
273 =head1 TRIGONOMETRIC FUNCTIONS
283 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
286 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
288 The arcus (also known as the inverse) functions of the sine, cosine,
291 B<asin>, B<acos>, B<atan>
293 The principal value of the arc tangent of y/x
297 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
298 and acotan/acot are aliases). Note that atan2(0, 0) is not well-defined.
300 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
302 The hyperbolic sine, cosine, and tangent
304 B<sinh>, B<cosh>, B<tanh>
306 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
307 and cotanh/coth are aliases)
309 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
311 The arcus (also known as the inverse) functions of the hyperbolic
312 sine, cosine, and tangent
314 B<asinh>, B<acosh>, B<atanh>
316 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
317 (acsch/acosech and acoth/acotanh are aliases)
319 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
321 The trigonometric constant B<pi> and some of handy multiples
322 of it are also defined.
324 B<pi, pi2, pi4, pip2, pip4>
326 =head2 ERRORS DUE TO DIVISION BY ZERO
328 The following functions
345 cannot be computed for all arguments because that would mean dividing
346 by zero or taking logarithm of zero. These situations cause fatal
347 runtime errors looking like this
349 cot(0): Division by zero.
350 (Because in the definition of cot(0), the divisor sin(0) is 0)
355 atanh(-1): Logarithm of zero.
358 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
359 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
360 C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
361 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
362 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
363 pi>, where I<k> is any integer.
365 Note that atan2(0, 0) is not well-defined.
367 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
369 Please note that some of the trigonometric functions can break out
370 from the B<real axis> into the B<complex plane>. For example
371 C<asin(2)> has no definition for plain real numbers but it has
372 definition for complex numbers.
374 In Perl terms this means that supplying the usual Perl numbers (also
375 known as scalars, please see L<perldata>) as input for the
376 trigonometric functions might produce as output results that no more
377 are simple real numbers: instead they are complex numbers.
379 The C<Math::Trig> handles this by using the C<Math::Complex> package
380 which knows how to handle complex numbers, please see L<Math::Complex>
381 for more information. In practice you need not to worry about getting
382 complex numbers as results because the C<Math::Complex> takes care of
383 details like for example how to display complex numbers. For example:
387 should produce something like this (take or leave few last decimals):
389 1.5707963267949-1.31695789692482i
391 That is, a complex number with the real part of approximately C<1.571>
392 and the imaginary part of approximately C<-1.317>.
394 =head1 PLANE ANGLE CONVERSIONS
396 (Plane, 2-dimensional) angles may be converted with the following functions.
402 $radians = deg2rad($degrees);
406 $radians = grad2rad($gradians);
410 $degrees = rad2deg($radians);
414 $degrees = grad2deg($gradians);
418 $gradians = deg2grad($degrees);
422 $gradians = rad2grad($radians);
426 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
427 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
428 If you don't want this, supply a true second argument:
430 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
431 $negative_degrees = rad2deg($negative_radians, 1);
433 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
440 $radians_wrapped_by_2pi = rad2rad($radians);
444 $degrees_wrapped_by_360 = deg2deg($degrees);
448 $gradians_wrapped_by_400 = grad2grad($gradians);
452 =head1 RADIAL COORDINATE CONVERSIONS
454 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
455 systems, explained shortly in more detail.
457 You can import radial coordinate conversion functions by using the
460 use Math::Trig ':radial';
462 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
463 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
464 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
465 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
466 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
467 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
469 B<All angles are in radians>.
471 =head2 COORDINATE SYSTEMS
473 B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
475 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
476 coordinates which define a point in three-dimensional space. They are
477 based on a sphere surface. The radius of the sphere is B<rho>, also
478 known as the I<radial> coordinate. The angle in the I<xy>-plane
479 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
480 coordinate. The angle from the I<z>-axis is B<phi>, also known as the
481 I<polar> coordinate. The North Pole is therefore I<0, 0, rho>, and
482 the Gulf of Guinea (think of the missing big chunk of Africa) I<0,
483 pi/2, rho>. In geographical terms I<phi> is latitude (northward
484 positive, southward negative) and I<theta> is longitude (eastward
485 positive, westward negative).
487 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
488 some texts define the I<phi> to start from the horizontal plane, some
489 texts use I<r> in place of I<rho>.
491 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
492 coordinates which define a point in three-dimensional space. They are
493 based on a cylinder surface. The radius of the cylinder is B<rho>,
494 also known as the I<radial> coordinate. The angle in the I<xy>-plane
495 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
496 coordinate. The third coordinate is the I<z>, pointing up from the
499 =head2 3-D ANGLE CONVERSIONS
501 Conversions to and from spherical and cylindrical coordinates are
502 available. Please notice that the conversions are not necessarily
503 reversible because of the equalities like I<pi> angles being equal to
508 =item cartesian_to_cylindrical
510 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
512 =item cartesian_to_spherical
514 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
516 =item cylindrical_to_cartesian
518 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
520 =item cylindrical_to_spherical
522 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
524 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
526 =item spherical_to_cartesian
528 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
530 =item spherical_to_cylindrical
532 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
534 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
538 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
540 A great circle is section of a circle that contains the circle
541 diameter: the shortest distance between two (non-antipodal) points on
542 the spherical surface goes along the great circle connecting those two
545 =head2 great_circle_distance
547 You can compute spherical distances, called B<great circle distances>,
548 by importing the great_circle_distance() function:
550 use Math::Trig 'great_circle_distance';
552 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
554 The I<great circle distance> is the shortest distance between two
555 points on a sphere. The distance is in C<$rho> units. The C<$rho> is
556 optional, it defaults to 1 (the unit sphere), therefore the distance
559 If you think geographically the I<theta> are longitudes: zero at the
560 Greenwhich meridian, eastward positive, westward negative -- and the
561 I<phi> are latitudes: zero at the North Pole, northward positive,
562 southward negative. B<NOTE>: this formula thinks in mathematics, not
563 geographically: the I<phi> zero is at the North Pole, not at the
564 Equator on the west coast of Africa (Bay of Guinea). You need to
565 subtract your geographical coordinates from I<pi/2> (also known as 90
568 $distance = great_circle_distance($lon0, pi/2 - $lat0,
569 $lon1, pi/2 - $lat1, $rho);
571 =head2 great_circle_direction
573 The direction you must follow the great circle (also known as I<bearing>)
574 can be computed by the great_circle_direction() function:
576 use Math::Trig 'great_circle_direction';
578 $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
580 =head2 great_circle_bearing
582 Alias 'great_circle_bearing' for 'great_circle_direction' is also available.
584 use Math::Trig 'great_circle_bearing';
586 $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1);
588 The result of great_circle_direction is in radians, zero indicating
589 straight north, pi or -pi straight south, pi/2 straight west, and
592 You can inversely compute the destination if you know the
593 starting point, direction, and distance:
595 =head2 great_circle_destination
597 use Math::Trig 'great_circle_destination';
599 # thetad and phid are the destination coordinates,
600 # dird is the final direction at the destination.
602 ($thetad, $phid, $dird) =
603 great_circle_destination($theta, $phi, $direction, $distance);
605 or the midpoint if you know the end points:
607 =head2 great_circle_midpoint
609 use Math::Trig 'great_circle_midpoint';
612 great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
614 The great_circle_midpoint() is just a special case of
616 =head2 great_circle_waypoint
618 use Math::Trig 'great_circle_waypoint';
621 great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
623 Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
624 $phi1). Note that antipodal points (where their distance is I<pi>
625 radians) do not have waypoints between them (they would have an an
626 "equator" between them), and therefore C<undef> is returned for
627 antipodal points. If the points are the same and the distance
628 therefore zero and all waypoints therefore identical, the first point
629 (either point) is returned.
631 The thetas, phis, direction, and distance in the above are all in radians.
633 You can import all the great circle formulas by
635 use Math::Trig ':great_circle';
637 Notice that the resulting directions might be somewhat surprising if
638 you are looking at a flat worldmap: in such map projections the great
639 circles quite often do not look like the shortest routes -- but for
640 example the shortest possible routes from Europe or North America to
641 Asia do often cross the polar regions. (The common Mercator projection
642 does B<not> show great circles as straight lines: straight lines in the
643 Mercator projection are lines of constant bearing.)
647 To calculate the distance between London (51.3N 0.5W) and Tokyo
648 (35.7N 139.8E) in kilometers:
650 use Math::Trig qw(great_circle_distance deg2rad);
652 # Notice the 90 - latitude: phi zero is at the North Pole.
653 sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
654 my @L = NESW( -0.5, 51.3);
655 my @T = NESW(139.8, 35.7);
656 my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
658 The direction you would have to go from London to Tokyo (in radians,
659 straight north being zero, straight east being pi/2).
661 use Math::Trig qw(great_circle_direction);
663 my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
665 The midpoint between London and Tokyo being
667 use Math::Trig qw(great_circle_midpoint);
669 my @M = great_circle_midpoint(@L, @T);
671 or about 68.93N 89.16E, in the frozen wastes of Siberia.
673 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
675 The answers may be off by few percentages because of the irregular
676 (slightly aspherical) form of the Earth. The errors are at worst
677 about 0.55%, but generally below 0.3%.
679 =head2 Real-valued asin and acos
681 For small inputs asin() and acos() may return complex numbers even
682 when real numbers would be enough and correct, this happens because of
683 floating-point inaccuracies. You can see these inaccuracies for
684 example by trying theses:
686 print cos(1e-6)**2+sin(1e-6)**2 - 1,"\n";
687 printf "%.20f", cos(1e-6)**2+sin(1e-6)**2,"\n";
689 which will print something like this
691 -1.11022302462516e-16
692 0.99999999999999988898
694 even though the expected results are of course exactly zero and one.
695 The formulas used to compute asin() and acos() are quite sensitive to
696 this, and therefore they might accidentally slip into the complex
697 plane even when they should not. To counter this there are two
698 interfaces that are guaranteed to return a real-valued output.
704 use Math::Trig qw(asin_real);
706 $real_angle = asin_real($input_sin);
708 Return a real-valued arcus sine if the input is between [-1, 1],
709 B<inclusive> the endpoints. For inputs greater than one, pi/2
710 is returned. For inputs less than minus one, -pi/2 is returned.
714 use Math::Trig qw(acos_real);
716 $real_angle = acos_real($input_cos);
718 Return a real-valued arcus cosine if the input is between [-1, 1],
719 B<inclusive> the endpoints. For inputs greater than one, zero
720 is returned. For inputs less than minus one, pi is returned.
726 Saying C<use Math::Trig;> exports many mathematical routines in the
727 caller environment and even overrides some (C<sin>, C<cos>). This is
728 construed as a feature by the Authors, actually... ;-)
730 The code is not optimized for speed, especially because we use
731 C<Math::Complex> and thus go quite near complex numbers while doing
732 the computations even when the arguments are not. This, however,
733 cannot be completely avoided if we want things like C<asin(2)> to give
734 an answer instead of giving a fatal runtime error.
736 Do not attempt navigation using these formulas.
742 Jarkko Hietaniemi <F<jhi!at!iki.fi>> and
743 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>.
747 This library is free software; you can redistribute it and/or modify
748 it under the same terms as Perl itself.