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 $lat0 = pip2 - $phi0;
137 my $lat1 = pip2 - $phi1;
140 atan2(sin($theta0 - $theta1) * cos($lat1),
141 cos($lat0) * sin($lat1) -
142 sin($lat0) * cos($lat1) * cos($theta0 - $theta1));
144 return rad2rad($direction);
154 Math::Trig - trigonometric functions
170 C<Math::Trig> defines many trigonometric functions not defined by the
171 core Perl which defines only the C<sin()> and C<cos()>. The constant
172 B<pi> is also defined as are a few convenience functions for angle
175 =head1 TRIGONOMETRIC FUNCTIONS
185 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
188 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
190 The arcus (also known as the inverse) functions of the sine, cosine,
193 B<asin>, B<acos>, B<atan>
195 The principal value of the arc tangent of y/x
199 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
200 and acotan/acot are aliases)
202 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
204 The hyperbolic sine, cosine, and tangent
206 B<sinh>, B<cosh>, B<tanh>
208 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
209 and cotanh/coth are aliases)
211 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
213 The arcus (also known as the inverse) functions of the hyperbolic
214 sine, cosine, and tangent
216 B<asinh>, B<acosh>, B<atanh>
218 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
219 (acsch/acosech and acoth/acotanh are aliases)
221 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
223 The trigonometric constant B<pi> is also defined.
227 =head2 ERRORS DUE TO DIVISION BY ZERO
229 The following functions
246 cannot be computed for all arguments because that would mean dividing
247 by zero or taking logarithm of zero. These situations cause fatal
248 runtime errors looking like this
250 cot(0): Division by zero.
251 (Because in the definition of cot(0), the divisor sin(0) is 0)
256 atanh(-1): Logarithm of zero.
259 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
260 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
261 C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
262 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
263 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
264 pi>, where I<k> is any integer.
266 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
268 Please note that some of the trigonometric functions can break out
269 from the B<real axis> into the B<complex plane>. For example
270 C<asin(2)> has no definition for plain real numbers but it has
271 definition for complex numbers.
273 In Perl terms this means that supplying the usual Perl numbers (also
274 known as scalars, please see L<perldata>) as input for the
275 trigonometric functions might produce as output results that no more
276 are simple real numbers: instead they are complex numbers.
278 The C<Math::Trig> handles this by using the C<Math::Complex> package
279 which knows how to handle complex numbers, please see L<Math::Complex>
280 for more information. In practice you need not to worry about getting
281 complex numbers as results because the C<Math::Complex> takes care of
282 details like for example how to display complex numbers. For example:
286 should produce something like this (take or leave few last decimals):
288 1.5707963267949-1.31695789692482i
290 That is, a complex number with the real part of approximately C<1.571>
291 and the imaginary part of approximately C<-1.317>.
293 =head1 PLANE ANGLE CONVERSIONS
295 (Plane, 2-dimensional) angles may be converted with the following functions.
297 $radians = deg2rad($degrees);
298 $radians = grad2rad($gradians);
300 $degrees = rad2deg($radians);
301 $degrees = grad2deg($gradians);
303 $gradians = deg2grad($degrees);
304 $gradians = rad2grad($radians);
306 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
307 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
308 If you don't want this, supply a true second argument:
310 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
311 $negative_degrees = rad2deg($negative_radians, 1);
313 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
316 =head1 RADIAL COORDINATE CONVERSIONS
318 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
319 systems, explained shortly in more detail.
321 You can import radial coordinate conversion functions by using the
324 use Math::Trig ':radial';
326 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
327 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
328 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
329 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
330 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
331 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
333 B<All angles are in radians>.
335 =head2 COORDINATE SYSTEMS
337 B<Cartesian> coordinates are the usual rectangular I<(x, y,
340 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
341 coordinates which define a point in three-dimensional space. They are
342 based on a sphere surface. The radius of the sphere is B<rho>, also
343 known as the I<radial> coordinate. The angle in the I<xy>-plane
344 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
345 coordinate. The angle from the I<z>-axis is B<phi>, also known as the
346 I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and
347 the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
348 pi/2, rho>. In geographical terms I<phi> is latitude (northward
349 positive, southward negative) and I<theta> is longitude (eastward
350 positive, westward negative).
352 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
353 some texts define the I<phi> to start from the horizontal plane, some
354 texts use I<r> in place of I<rho>.
356 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
357 coordinates which define a point in three-dimensional space. They are
358 based on a cylinder surface. The radius of the cylinder is B<rho>,
359 also known as the I<radial> coordinate. The angle in the I<xy>-plane
360 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
361 coordinate. The third coordinate is the I<z>, pointing up from the
364 =head2 3-D ANGLE CONVERSIONS
366 Conversions to and from spherical and cylindrical coordinates are
367 available. Please notice that the conversions are not necessarily
368 reversible because of the equalities like I<pi> angles being equal to
373 =item cartesian_to_cylindrical
375 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
377 =item cartesian_to_spherical
379 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
381 =item cylindrical_to_cartesian
383 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
385 =item cylindrical_to_spherical
387 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
389 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
391 =item spherical_to_cartesian
393 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
395 =item spherical_to_cylindrical
397 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
399 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
403 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
405 You can compute spherical distances, called B<great circle distances>,
406 by importing the great_circle_distance() function:
408 use Math::Trig 'great_circle_distance';
410 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
412 The I<great circle distance> is the shortest distance between two
413 points on a sphere. The distance is in C<$rho> units. The C<$rho> is
414 optional, it defaults to 1 (the unit sphere), therefore the distance
417 If you think geographically the I<theta> are longitudes: zero at the
418 Greenwhich meridian, eastward positive, westward negative--and the
419 I<phi> are latitudes: zero at the North Pole, northward positive,
420 southward negative. B<NOTE>: this formula thinks in mathematics, not
421 geographically: the I<phi> zero is at the North Pole, not at the
422 Equator on the west coast of Africa (Bay of Guinea). You need to
423 subtract your geographical coordinates from I<pi/2> (also known as 90
426 $distance = great_circle_distance($lon0, pi/2 - $lat0,
427 $lon1, pi/2 - $lat1, $rho);
429 The direction you must follow the great circle can be computed by the
430 great_circle_direction() function:
432 use Math::Trig 'great_circle_direction';
434 $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
436 The result is in radians, zero indicating straight north, pi or -pi
437 straight south, pi/2 straight west, and -pi/2 straight east.
439 Notice that the resulting directions might be somewhat surprising if
440 you are looking at a flat worldmap: in such map projections the great
441 circles quite often do not look like the shortest routes-- but for
442 example the shortest possible routes from Europe or North America to
443 Asia do often cross the polar regions.
447 To calculate the distance between London (51.3N 0.5W) and Tokyo
448 (35.7N 139.8E) in kilometers:
450 use Math::Trig qw(great_circle_distance deg2rad);
452 # Notice the 90 - latitude: phi zero is at the North Pole.
453 @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
454 @T = (deg2rad(139.8),deg2rad(90 - 35.7));
456 $km = great_circle_distance(@L, @T, 6378);
458 The direction you would have to go from London to Tokyo
460 use Math::Trig qw(great_circle_direction);
462 $rad = great_circle_direction(@L, @T);
464 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
466 The answers may be off by few percentages because of the irregular
467 (slightly aspherical) form of the Earth. The formula used for
468 grear circle distances
470 lat0 = 90 degrees - phi0
471 lat1 = 90 degrees - phi1
472 d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
473 sin(lat0) * sin(lat1))
475 is also somewhat unreliable for small distances (for locations
476 separated less than about five degrees) because it uses arc cosine
477 which is rather ill-conditioned for values close to zero.
481 Saying C<use Math::Trig;> exports many mathematical routines in the
482 caller environment and even overrides some (C<sin>, C<cos>). This is
483 construed as a feature by the Authors, actually... ;-)
485 The code is not optimized for speed, especially because we use
486 C<Math::Complex> and thus go quite near complex numbers while doing
487 the computations even when the arguments are not. This, however,
488 cannot be completely avoided if we want things like C<asin(2)> to give
489 an answer instead of giving a fatal runtime error.
493 Jarkko Hietaniemi <F<jhi@iki.fi>> and
494 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.