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.44;
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 @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
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);
38 great_circle_direction
42 great_circle_destination
45 my @pi = qw(pi pi2 pi4 pip2 pip4);
47 @EXPORT_OK = (@rdlcnv, @greatcircle, @pi, 'Inf');
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
53 %EXPORT_TAGS = ('radial' => [ @rdlcnv ],
54 'great_circle' => [ @greatcircle ],
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 }
65 # Truncating remainder.
69 # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
70 $_[0] - $_[1] * int($_[0] / $_[1]);
77 sub rad2rad($) { _remt($_[0], pi2) }
79 sub deg2deg($) { _remt($_[0], 360) }
81 sub grad2grad($) { _remt($_[0], 400) }
83 sub rad2deg ($;$) { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) }
85 sub deg2rad ($;$) { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) }
87 sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) }
89 sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) }
91 sub rad2grad ($;$) { my $d = _RG * $_[0]; $_[1] ? $d : grad2grad($d) }
93 sub grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) }
95 sub cartesian_to_spherical {
96 my ( $x, $y, $z ) = @_;
98 my $rho = sqrt( $x * $x + $y * $y + $z * $z );
102 $rho ? acos( $z / $rho ) : 0 );
105 sub spherical_to_cartesian {
106 my ( $rho, $theta, $phi ) = @_;
108 return ( $rho * cos( $theta ) * sin( $phi ),
109 $rho * sin( $theta ) * sin( $phi ),
110 $rho * cos( $phi ) );
113 sub spherical_to_cylindrical {
114 my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
116 return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
119 sub cartesian_to_cylindrical {
120 my ( $x, $y, $z ) = @_;
122 return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
125 sub cylindrical_to_cartesian {
126 my ( $rho, $theta, $z ) = @_;
128 return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
131 sub cylindrical_to_spherical {
132 return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
135 sub great_circle_distance {
136 my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
138 $rho = 1 unless defined $rho; # Default to the unit sphere.
140 my $lat0 = pip2 - $phi0;
141 my $lat1 = pip2 - $phi1;
144 acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
145 sin( $lat0 ) * sin( $lat1 ) );
148 sub great_circle_direction {
149 my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
151 my $distance = &great_circle_distance;
153 my $lat0 = pip2 - $phi0;
154 my $lat1 = pip2 - $phi1;
157 acos((sin($lat1) - sin($lat0) * cos($distance)) /
158 (cos($lat0) * sin($distance)));
160 $direction = pi2 - $direction
161 if sin($theta1 - $theta0) < 0;
163 return rad2rad($direction);
166 *great_circle_bearing = \&great_circle_direction;
168 sub great_circle_waypoint {
169 my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
171 $point = 0.5 unless defined $point;
173 my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
175 return undef if $d == pi;
179 return ($theta0, $phi0) if $sd == 0;
181 my $A = sin((1 - $point) * $d) / $sd;
182 my $B = sin( $point * $d) / $sd;
184 my $lat0 = pip2 - $phi0;
185 my $lat1 = pip2 - $phi1;
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);
191 my $theta = atan2($y, $x);
194 return ($theta, $phi);
197 sub great_circle_midpoint {
198 great_circle_waypoint(@_[0..3], 0.5);
201 sub great_circle_destination {
202 my ( $theta0, $phi0, $dir0, $dst ) = @_;
204 my $lat0 = pip2 - $phi0;
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));
210 my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
212 $dir1 -= pi2 if $dir1 > pi2;
214 return ($theta1, $phi1, $dir1);
224 Math::Trig - trigonometric functions
238 # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
239 use Math::Trig ':pi';
241 # Import the conversions between cartesian/spherical/cylindrical.
242 use Math::Trig ':radial';
244 # Import the great circle formulas.
245 use Math::Trig ':great_circle';
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.
254 =head1 TRIGONOMETRIC FUNCTIONS
264 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
267 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
269 The arcus (also known as the inverse) functions of the sine, cosine,
272 B<asin>, B<acos>, B<atan>
274 The principal value of the arc tangent of y/x
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.
281 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
283 The hyperbolic sine, cosine, and tangent
285 B<sinh>, B<cosh>, B<tanh>
287 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
288 and cotanh/coth are aliases)
290 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
292 The arcus (also known as the inverse) functions of the hyperbolic
293 sine, cosine, and tangent
295 B<asinh>, B<acosh>, B<atanh>
297 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
298 (acsch/acosech and acoth/acotanh are aliases)
300 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
302 The trigonometric constant B<pi> and some of handy multiples
303 of it are also defined.
305 B<pi, pi2, pi4, pip2, pip4>
307 =head2 ERRORS DUE TO DIVISION BY ZERO
309 The following functions
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
330 cot(0): Division by zero.
331 (Because in the definition of cot(0), the divisor sin(0) is 0)
336 atanh(-1): Logarithm of zero.
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.
346 Note that atan2(0, 0) is not well-defined.
348 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
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.
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.
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:
368 should produce something like this (take or leave few last decimals):
370 1.5707963267949-1.31695789692482i
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>.
375 =head1 PLANE ANGLE CONVERSIONS
377 (Plane, 2-dimensional) angles may be converted with the following functions.
383 $radians = deg2rad($degrees);
387 $radians = grad2rad($gradians);
391 $degrees = rad2deg($radians);
395 $degrees = grad2deg($gradians);
399 $gradians = deg2grad($degrees);
403 $gradians = rad2grad($radians);
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:
411 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
412 $negative_degrees = rad2deg($negative_radians, 1);
414 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
421 $radians_wrapped_by_2pi = rad2rad($radians);
425 $degrees_wrapped_by_360 = deg2deg($degrees);
429 $gradians_wrapped_by_400 = grad2grad($gradians);
433 =head1 RADIAL COORDINATE CONVERSIONS
435 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
436 systems, explained shortly in more detail.
438 You can import radial coordinate conversion functions by using the
441 use Math::Trig ':radial';
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);
450 B<All angles are in radians>.
452 =head2 COORDINATE SYSTEMS
454 B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
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).
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>.
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
480 =head2 3-D ANGLE CONVERSIONS
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
489 =item cartesian_to_cylindrical
491 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
493 =item cartesian_to_spherical
495 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
497 =item cylindrical_to_cartesian
499 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
501 =item cylindrical_to_spherical
503 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
505 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
507 =item spherical_to_cartesian
509 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
511 =item spherical_to_cylindrical
513 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
515 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
519 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
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
526 =head2 great_circle_distance
528 You can compute spherical distances, called B<great circle distances>,
529 by importing the great_circle_distance() function:
531 use Math::Trig 'great_circle_distance';
533 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
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
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
549 $distance = great_circle_distance($lon0, pi/2 - $lat0,
550 $lon1, pi/2 - $lat1, $rho);
552 =head2 great_circle_direction
554 The direction you must follow the great circle (also known as I<bearing>)
555 can be computed by the great_circle_direction() function:
557 use Math::Trig 'great_circle_direction';
559 $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
561 =head2 great_circle_bearing
563 Alias 'great_circle_bearing' for 'great_circle_direction' is also available.
565 use Math::Trig 'great_circle_bearing';
567 $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1);
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
573 You can inversely compute the destination if you know the
574 starting point, direction, and distance:
576 =head2 great_circle_destination
578 use Math::Trig 'great_circle_destination';
580 # thetad and phid are the destination coordinates,
581 # dird is the final direction at the destination.
583 ($thetad, $phid, $dird) =
584 great_circle_destination($theta, $phi, $direction, $distance);
586 or the midpoint if you know the end points:
588 =head2 great_circle_midpoint
590 use Math::Trig 'great_circle_midpoint';
593 great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
595 The great_circle_midpoint() is just a special case of
597 =head2 great_circle_waypoint
599 use Math::Trig 'great_circle_waypoint';
602 great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
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.
612 The thetas, phis, direction, and distance in the above are all in radians.
614 You can import all the great circle formulas by
616 use Math::Trig ':great_circle';
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.
626 To calculate the distance between London (51.3N 0.5W) and Tokyo
627 (35.7N 139.8E) in kilometers:
629 use Math::Trig qw(great_circle_distance deg2rad);
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.
637 The direction you would have to go from London to Tokyo (in radians,
638 straight north being zero, straight east being pi/2).
640 use Math::Trig qw(great_circle_direction);
642 my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
644 The midpoint between London and Tokyo being
646 use Math::Trig qw(great_circle_midpoint);
648 my @M = great_circle_midpoint(@L, @T);
650 or about 68.93N 89.16E, in the frozen wastes of Siberia.
652 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
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%.
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... ;-)
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.
670 Do not attempt navigation using these formulas.
676 Jarkko Hietaniemi <F<jhi!at!iki.fi>> and
677 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>.
681 This library is free software; you can redistribute it and/or modify
682 it under the same terms as Perl itself.