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.55;
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($theta0, $phi0, $theta1, $phi1);
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));
227 my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
228 cos($dst)-sin($lat0)*sin($phi1));
230 my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
232 $dir1 -= pi2 if $dir1 > pi2;
234 return ($theta1, $phi1, $dir1);
244 Math::Trig - trigonometric functions
258 # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
259 use Math::Trig ':pi';
261 # Import the conversions between cartesian/spherical/cylindrical.
262 use Math::Trig ':radial';
264 # Import the great circle formulas.
265 use Math::Trig ':great_circle';
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.
274 =head1 TRIGONOMETRIC FUNCTIONS
284 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
287 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
289 The arcus (also known as the inverse) functions of the sine, cosine,
292 B<asin>, B<acos>, B<atan>
294 The principal value of the arc tangent of y/x
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.
301 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
303 The hyperbolic sine, cosine, and tangent
305 B<sinh>, B<cosh>, B<tanh>
307 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
308 and cotanh/coth are aliases)
310 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
312 The area (also known as the inverse) functions of the hyperbolic
313 sine, cosine, and tangent
315 B<asinh>, B<acosh>, B<atanh>
317 The area cofunctions of the hyperbolic sine, cosine, and tangent
318 (acsch/acosech and acoth/acotanh are aliases)
320 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
322 The trigonometric constant B<pi> and some of handy multiples
323 of it are also defined.
325 B<pi, pi2, pi4, pip2, pip4>
327 =head2 ERRORS DUE TO DIVISION BY ZERO
329 The following functions
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
350 cot(0): Division by zero.
351 (Because in the definition of cot(0), the divisor sin(0) is 0)
356 atanh(-1): Logarithm of zero.
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.
366 Note that atan2(0, 0) is not well-defined.
368 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
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.
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.
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:
388 should produce something like this (take or leave few last decimals):
390 1.5707963267949-1.31695789692482i
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>.
395 =head1 PLANE ANGLE CONVERSIONS
397 (Plane, 2-dimensional) angles may be converted with the following functions.
403 $radians = deg2rad($degrees);
407 $radians = grad2rad($gradians);
411 $degrees = rad2deg($radians);
415 $degrees = grad2deg($gradians);
419 $gradians = deg2grad($degrees);
423 $gradians = rad2grad($radians);
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:
431 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
432 $negative_degrees = rad2deg($negative_radians, 1);
434 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
441 $radians_wrapped_by_2pi = rad2rad($radians);
445 $degrees_wrapped_by_360 = deg2deg($degrees);
449 $gradians_wrapped_by_400 = grad2grad($gradians);
453 =head1 RADIAL COORDINATE CONVERSIONS
455 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
456 systems, explained shortly in more detail.
458 You can import radial coordinate conversion functions by using the
461 use Math::Trig ':radial';
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);
470 B<All angles are in radians>.
472 =head2 COORDINATE SYSTEMS
474 B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
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).
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>.
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
500 =head2 3-D ANGLE CONVERSIONS
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
509 =item cartesian_to_cylindrical
511 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
513 =item cartesian_to_spherical
515 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
517 =item cylindrical_to_cartesian
519 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
521 =item cylindrical_to_spherical
523 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
525 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
527 =item spherical_to_cartesian
529 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
531 =item spherical_to_cylindrical
533 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
535 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
539 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
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
546 =head2 great_circle_distance
548 You can compute spherical distances, called B<great circle distances>,
549 by importing the great_circle_distance() function:
551 use Math::Trig 'great_circle_distance';
553 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
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
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
569 $distance = great_circle_distance($lon0, pi/2 - $lat0,
570 $lon1, pi/2 - $lat1, $rho);
572 =head2 great_circle_direction
574 The direction you must follow the great circle (also known as I<bearing>)
575 can be computed by the great_circle_direction() function:
577 use Math::Trig 'great_circle_direction';
579 $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
581 =head2 great_circle_bearing
583 Alias 'great_circle_bearing' for 'great_circle_direction' is also available.
585 use Math::Trig 'great_circle_bearing';
587 $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1);
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
593 =head2 great_circle_destination
595 You can inversely compute the destination if you know the
596 starting point, direction, and distance:
598 use Math::Trig 'great_circle_destination';
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.
607 ($thetad, $phid, $dird) =
608 great_circle_destination($theta, $phi, $diro, $distance);
610 or the midpoint if you know the end points:
612 =head2 great_circle_midpoint
614 use Math::Trig 'great_circle_midpoint';
617 great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
619 The great_circle_midpoint() is just a special case of
621 =head2 great_circle_waypoint
623 use Math::Trig 'great_circle_waypoint';
626 great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
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.
636 The thetas, phis, direction, and distance in the above are all in radians.
638 You can import all the great circle formulas by
640 use Math::Trig ':great_circle';
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.)
652 To calculate the distance between London (51.3N 0.5W) and Tokyo
653 (35.7N 139.8E) in kilometers:
655 use Math::Trig qw(great_circle_distance deg2rad);
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.
663 The direction you would have to go from London to Tokyo (in radians,
664 straight north being zero, straight east being pi/2).
666 use Math::Trig qw(great_circle_direction);
668 my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
670 The midpoint between London and Tokyo being
672 use Math::Trig qw(great_circle_midpoint);
674 my @M = great_circle_midpoint(@L, @T);
676 or about 69 N 89 E, in the frozen wastes of Siberia.
678 B<NOTE>: you B<cannot> get from A to B like this:
680 Dist = great_circle_distance(A, B)
681 Dir = great_circle_direction(A, B)
682 C = great_circle_destination(A, Dist, Dir)
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.
689 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
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%.
695 =head2 Real-valued asin and acos
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:
702 print cos(1e-6)**2+sin(1e-6)**2 - 1,"\n";
703 printf "%.20f", cos(1e-6)**2+sin(1e-6)**2,"\n";
705 which will print something like this
707 -1.11022302462516e-16
708 0.99999999999999988898
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.
720 use Math::Trig qw(asin_real);
722 $real_angle = asin_real($input_sin);
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.
730 use Math::Trig qw(acos_real);
732 $real_angle = acos_real($input_cos);
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.
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... ;-)
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.
752 Do not attempt navigation using these formulas.
758 Jarkko Hietaniemi <F<jhi!at!iki.fi>> and
759 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>.
763 This library is free software; you can redistribute it and/or modify
764 it under the same terms as Perl itself.