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 qw(:trig);
15 our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
21 my @angcnv = qw(rad2deg rad2grad
25 @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
28 my @rdlcnv = qw(cartesian_to_cylindrical
29 cartesian_to_spherical
30 cylindrical_to_cartesian
31 cylindrical_to_spherical
32 spherical_to_cartesian
33 spherical_to_cylindrical);
35 @EXPORT_OK = (@rdlcnv, 'great_circle_distance', 'great_circle_direction');
37 %EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
40 sub pip2 () { pi / 2 }
50 # Truncating remainder.
54 # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
55 $_[0] - $_[1] * int($_[0] / $_[1]);
62 sub rad2rad($) { remt($_[0], pi2) }
64 sub deg2deg($) { remt($_[0], 360) }
66 sub grad2grad($) { remt($_[0], 400) }
68 sub rad2deg ($;$) { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) }
70 sub deg2rad ($;$) { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) }
72 sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) }
74 sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) }
76 sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) }
78 sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) }
80 sub cartesian_to_spherical {
81 my ( $x, $y, $z ) = @_;
83 my $rho = sqrt( $x * $x + $y * $y + $z * $z );
87 $rho ? acos( $z / $rho ) : 0 );
90 sub spherical_to_cartesian {
91 my ( $rho, $theta, $phi ) = @_;
93 return ( $rho * cos( $theta ) * sin( $phi ),
94 $rho * sin( $theta ) * sin( $phi ),
98 sub spherical_to_cylindrical {
99 my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
101 return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
104 sub cartesian_to_cylindrical {
105 my ( $x, $y, $z ) = @_;
107 return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
110 sub cylindrical_to_cartesian {
111 my ( $rho, $theta, $z ) = @_;
113 return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
116 sub cylindrical_to_spherical {
117 return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
120 sub great_circle_distance {
121 my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
123 $rho = 1 unless defined $rho; # Default to the unit sphere.
125 my $lat0 = pip2 - $phi0;
126 my $lat1 = pip2 - $phi1;
129 acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
130 sin( $lat0 ) * sin( $lat1 ) );
133 sub great_circle_direction {
134 my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
136 my $distance = &great_circle_distance;
138 my $lat0 = pip2 - $phi0;
139 my $lat1 = pip2 - $phi1;
142 acos((sin($lat1) - sin($lat0) * cos($distance)) /
143 (cos($lat0) * sin($distance)));
145 $direction = pi2 - $direction
146 if sin($theta1 - $theta0) < 0;
148 return rad2rad($direction);
158 Math::Trig - trigonometric functions
174 C<Math::Trig> defines many trigonometric functions not defined by the
175 core Perl which defines only the C<sin()> and C<cos()>. The constant
176 B<pi> is also defined as are a few convenience functions for angle
179 =head1 TRIGONOMETRIC FUNCTIONS
189 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
192 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
194 The arcus (also known as the inverse) functions of the sine, cosine,
197 B<asin>, B<acos>, B<atan>
199 The principal value of the arc tangent of y/x
203 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
204 and acotan/acot are aliases)
206 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
208 The hyperbolic sine, cosine, and tangent
210 B<sinh>, B<cosh>, B<tanh>
212 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
213 and cotanh/coth are aliases)
215 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
217 The arcus (also known as the inverse) functions of the hyperbolic
218 sine, cosine, and tangent
220 B<asinh>, B<acosh>, B<atanh>
222 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
223 (acsch/acosech and acoth/acotanh are aliases)
225 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
227 The trigonometric constant B<pi> is also defined.
231 =head2 ERRORS DUE TO DIVISION BY ZERO
233 The following functions
250 cannot be computed for all arguments because that would mean dividing
251 by zero or taking logarithm of zero. These situations cause fatal
252 runtime errors looking like this
254 cot(0): Division by zero.
255 (Because in the definition of cot(0), the divisor sin(0) is 0)
260 atanh(-1): Logarithm of zero.
263 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
264 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
265 C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
266 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
267 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
268 pi>, where I<k> is any integer.
270 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
272 Please note that some of the trigonometric functions can break out
273 from the B<real axis> into the B<complex plane>. For example
274 C<asin(2)> has no definition for plain real numbers but it has
275 definition for complex numbers.
277 In Perl terms this means that supplying the usual Perl numbers (also
278 known as scalars, please see L<perldata>) as input for the
279 trigonometric functions might produce as output results that no more
280 are simple real numbers: instead they are complex numbers.
282 The C<Math::Trig> handles this by using the C<Math::Complex> package
283 which knows how to handle complex numbers, please see L<Math::Complex>
284 for more information. In practice you need not to worry about getting
285 complex numbers as results because the C<Math::Complex> takes care of
286 details like for example how to display complex numbers. For example:
290 should produce something like this (take or leave few last decimals):
292 1.5707963267949-1.31695789692482i
294 That is, a complex number with the real part of approximately C<1.571>
295 and the imaginary part of approximately C<-1.317>.
297 =head1 PLANE ANGLE CONVERSIONS
299 (Plane, 2-dimensional) angles may be converted with the following functions.
301 $radians = deg2rad($degrees);
302 $radians = grad2rad($gradians);
304 $degrees = rad2deg($radians);
305 $degrees = grad2deg($gradians);
307 $gradians = deg2grad($degrees);
308 $gradians = rad2grad($radians);
310 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
311 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
312 If you don't want this, supply a true second argument:
314 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
315 $negative_degrees = rad2deg($negative_radians, 1);
317 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
320 =head1 RADIAL COORDINATE CONVERSIONS
322 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
323 systems, explained shortly in more detail.
325 You can import radial coordinate conversion functions by using the
328 use Math::Trig ':radial';
330 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
331 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
332 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
333 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
334 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
335 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
337 B<All angles are in radians>.
339 =head2 COORDINATE SYSTEMS
341 B<Cartesian> coordinates are the usual rectangular I<(x, y,
344 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
345 coordinates which define a point in three-dimensional space. They are
346 based on a sphere surface. The radius of the sphere is B<rho>, also
347 known as the I<radial> coordinate. The angle in the I<xy>-plane
348 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
349 coordinate. The angle from the I<z>-axis is B<phi>, also known as the
350 I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and
351 the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
352 pi/2, rho>. In geographical terms I<phi> is latitude (northward
353 positive, southward negative) and I<theta> is longitude (eastward
354 positive, westward negative).
356 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
357 some texts define the I<phi> to start from the horizontal plane, some
358 texts use I<r> in place of I<rho>.
360 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
361 coordinates which define a point in three-dimensional space. They are
362 based on a cylinder surface. The radius of the cylinder is B<rho>,
363 also known as the I<radial> coordinate. The angle in the I<xy>-plane
364 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
365 coordinate. The third coordinate is the I<z>, pointing up from the
368 =head2 3-D ANGLE CONVERSIONS
370 Conversions to and from spherical and cylindrical coordinates are
371 available. Please notice that the conversions are not necessarily
372 reversible because of the equalities like I<pi> angles being equal to
377 =item cartesian_to_cylindrical
379 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
381 =item cartesian_to_spherical
383 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
385 =item cylindrical_to_cartesian
387 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
389 =item cylindrical_to_spherical
391 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
393 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
395 =item spherical_to_cartesian
397 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
399 =item spherical_to_cylindrical
401 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
403 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
407 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
409 You can compute spherical distances, called B<great circle distances>,
410 by importing the great_circle_distance() function:
412 use Math::Trig 'great_circle_distance';
414 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
416 The I<great circle distance> is the shortest distance between two
417 points on a sphere. The distance is in C<$rho> units. The C<$rho> is
418 optional, it defaults to 1 (the unit sphere), therefore the distance
421 If you think geographically the I<theta> are longitudes: zero at the
422 Greenwhich meridian, eastward positive, westward negative--and the
423 I<phi> are latitudes: zero at the North Pole, northward positive,
424 southward negative. B<NOTE>: this formula thinks in mathematics, not
425 geographically: the I<phi> zero is at the North Pole, not at the
426 Equator on the west coast of Africa (Bay of Guinea). You need to
427 subtract your geographical coordinates from I<pi/2> (also known as 90
430 $distance = great_circle_distance($lon0, pi/2 - $lat0,
431 $lon1, pi/2 - $lat1, $rho);
433 The direction you must follow the great circle can be computed by the
434 great_circle_direction() function:
436 use Math::Trig 'great_circle_direction';
438 $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
440 The result is in radians, zero indicating straight north, pi or -pi
441 straight south, pi/2 straight west, and -pi/2 straight east.
443 Notice that the resulting directions might be somewhat surprising if
444 you are looking at a flat worldmap: in such map projections the great
445 circles quite often do not look like the shortest routes-- but for
446 example the shortest possible routes from Europe or North America to
447 Asia do often cross the polar regions.
451 To calculate the distance between London (51.3N 0.5W) and Tokyo
452 (35.7N 139.8E) in kilometers:
454 use Math::Trig qw(great_circle_distance deg2rad);
456 # Notice the 90 - latitude: phi zero is at the North Pole.
457 @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
458 @T = (deg2rad(139.8),deg2rad(90 - 35.7));
460 $km = great_circle_distance(@L, @T, 6378);
462 The direction you would have to go from London to Tokyo
464 use Math::Trig qw(great_circle_direction);
466 $rad = great_circle_direction(@L, @T);
468 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
470 The answers may be off by few percentages because of the irregular
471 (slightly aspherical) form of the Earth. The formula used for
472 grear circle distances
474 lat0 = 90 degrees - phi0
475 lat1 = 90 degrees - phi1
476 d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
477 sin(lat0) * sin(lat1))
479 is also somewhat unreliable for small distances (for locations
480 separated less than about five degrees) because it uses arc cosine
481 which is rather ill-conditioned for values close to zero.
485 Saying C<use Math::Trig;> exports many mathematical routines in the
486 caller environment and even overrides some (C<sin>, C<cos>). This is
487 construed as a feature by the Authors, actually... ;-)
489 The code is not optimized for speed, especially because we use
490 C<Math::Complex> and thus go quite near complex numbers while doing
491 the computations even when the arguments are not. This, however,
492 cannot be completely avoided if we want things like C<asin(2)> to give
493 an answer instead of giving a fatal runtime error.
497 Jarkko Hietaniemi <F<jhi@iki.fi>> and
498 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.