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');
37 %EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
39 sub pi2 () { 2 * pi } # use constant generates warning
40 sub pip2 () { pi / 2 } # use constant generates warning
41 use constant DR => pi2/360;
42 use constant RD => 360/pi2;
43 use constant DG => 400/360;
44 use constant GD => 360/400;
45 use constant RG => 400/pi2;
46 use constant GR => pi2/400;
49 # Truncating remainder.
53 # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
54 $_[0] - $_[1] * int($_[0] / $_[1]);
61 sub rad2deg ($) { remt(RD * $_[0], 360) }
63 sub deg2rad ($) { remt(DR * $_[0], pi2) }
65 sub grad2deg ($) { remt(GD * $_[0], 360) }
67 sub deg2grad ($) { remt(DG * $_[0], 400) }
69 sub rad2grad ($) { remt(RG * $_[0], 400) }
71 sub grad2rad ($) { remt(GR * $_[0], pi2) }
73 sub cartesian_to_spherical {
74 my ( $x, $y, $z ) = @_;
76 my $rho = sqrt( $x * $x + $y * $y + $z * $z );
80 $rho ? acos( $z / $rho ) : 0 );
83 sub spherical_to_cartesian {
84 my ( $rho, $theta, $phi ) = @_;
86 return ( $rho * cos( $theta ) * sin( $phi ),
87 $rho * sin( $theta ) * sin( $phi ),
91 sub spherical_to_cylindrical {
92 my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
94 return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
97 sub cartesian_to_cylindrical {
98 my ( $x, $y, $z ) = @_;
100 return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
103 sub cylindrical_to_cartesian {
104 my ( $rho, $theta, $z ) = @_;
106 return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
109 sub cylindrical_to_spherical {
110 return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
113 sub great_circle_distance {
114 my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
116 $rho = 1 unless defined $rho; # Default to the unit sphere.
118 my $lat0 = pip2 - $phi0;
119 my $lat1 = pip2 - $phi1;
122 acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
123 sin( $lat0 ) * sin( $lat1 ) );
130 Math::Trig - trigonometric functions
146 C<Math::Trig> defines many trigonometric functions not defined by the
147 core Perl which defines only the C<sin()> and C<cos()>. The constant
148 B<pi> is also defined as are a few convenience functions for angle
151 =head1 TRIGONOMETRIC FUNCTIONS
161 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
164 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
166 The arcus (also known as the inverse) functions of the sine, cosine,
169 B<asin>, B<acos>, B<atan>
171 The principal value of the arc tangent of y/x
175 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
176 and acotan/acot are aliases)
178 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
180 The hyperbolic sine, cosine, and tangent
182 B<sinh>, B<cosh>, B<tanh>
184 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
185 and cotanh/coth are aliases)
187 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
189 The arcus (also known as the inverse) functions of the hyperbolic
190 sine, cosine, and tangent
192 B<asinh>, B<acosh>, B<atanh>
194 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
195 (acsch/acosech and acoth/acotanh are aliases)
197 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
199 The trigonometric constant B<pi> is also defined.
203 =head2 ERRORS DUE TO DIVISION BY ZERO
205 The following functions
222 cannot be computed for all arguments because that would mean dividing
223 by zero or taking logarithm of zero. These situations cause fatal
224 runtime errors looking like this
226 cot(0): Division by zero.
227 (Because in the definition of cot(0), the divisor sin(0) is 0)
232 atanh(-1): Logarithm of zero.
235 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
236 C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
237 C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
238 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
239 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
240 pi>, where I<k> is any integer.
242 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
244 Please note that some of the trigonometric functions can break out
245 from the B<real axis> into the B<complex plane>. For example
246 C<asin(2)> has no definition for plain real numbers but it has
247 definition for complex numbers.
249 In Perl terms this means that supplying the usual Perl numbers (also
250 known as scalars, please see L<perldata>) as input for the
251 trigonometric functions might produce as output results that no more
252 are simple real numbers: instead they are complex numbers.
254 The C<Math::Trig> handles this by using the C<Math::Complex> package
255 which knows how to handle complex numbers, please see L<Math::Complex>
256 for more information. In practice you need not to worry about getting
257 complex numbers as results because the C<Math::Complex> takes care of
258 details like for example how to display complex numbers. For example:
262 should produce something like this (take or leave few last decimals):
264 1.5707963267949-1.31695789692482i
266 That is, a complex number with the real part of approximately C<1.571>
267 and the imaginary part of approximately C<-1.317>.
269 =head1 PLANE ANGLE CONVERSIONS
271 (Plane, 2-dimensional) angles may be converted with the following functions.
273 $radians = deg2rad($degrees);
274 $radians = grad2rad($gradians);
276 $degrees = rad2deg($radians);
277 $degrees = grad2deg($gradians);
279 $gradians = deg2grad($degrees);
280 $gradians = rad2grad($radians);
282 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
284 =head1 RADIAL COORDINATE CONVERSIONS
286 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
287 systems, explained shortly in more detail.
289 You can import radial coordinate conversion functions by using the
292 use Math::Trig ':radial';
294 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
295 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
296 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
297 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
298 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
299 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
301 B<All angles are in radians>.
303 =head2 COORDINATE SYSTEMS
305 B<Cartesian> coordinates are the usual rectangular I<(x, y,
308 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
309 coordinates which define a point in three-dimensional space. They are
310 based on a sphere surface. The radius of the sphere is B<rho>, also
311 known as the I<radial> coordinate. The angle in the I<xy>-plane
312 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
313 coordinate. The angle from the I<z>-axis is B<phi>, also known as the
314 I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and
315 the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
316 pi/2, rho>. In geographical terms I<phi> is latitude (northward
317 positive, southward negative) and I<theta> is longitude (eastward
318 positive, westward negative).
320 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
321 some texts define the I<phi> to start from the horizontal plane, some
322 texts use I<r> in place of I<rho>.
324 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
325 coordinates which define a point in three-dimensional space. They are
326 based on a cylinder surface. The radius of the cylinder is B<rho>,
327 also known as the I<radial> coordinate. The angle in the I<xy>-plane
328 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
329 coordinate. The third coordinate is the I<z>, pointing up from the
332 =head2 3-D ANGLE CONVERSIONS
334 Conversions to and from spherical and cylindrical coordinates are
335 available. Please notice that the conversions are not necessarily
336 reversible because of the equalities like I<pi> angles being equal to
341 =item cartesian_to_cylindrical
343 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
345 =item cartesian_to_spherical
347 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
349 =item cylindrical_to_cartesian
351 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
353 =item cylindrical_to_spherical
355 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
357 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
359 =item spherical_to_cartesian
361 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
363 =item spherical_to_cylindrical
365 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
367 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
371 =head1 GREAT CIRCLE DISTANCES
373 You can compute spherical distances, called B<great circle distances>,
374 by importing the C<great_circle_distance> function:
376 use Math::Trig 'great_circle_distance'
378 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
380 The I<great circle distance> is the shortest distance between two
381 points on a sphere. The distance is in C<$rho> units. The C<$rho> is
382 optional, it defaults to 1 (the unit sphere), therefore the distance
385 If you think geographically the I<theta> are longitudes: zero at the
386 Greenwhich meridian, eastward positive, westward negative--and the
387 I<phi> are latitudes: zero at the North Pole, northward positive,
388 southward negative. B<NOTE>: this formula thinks in mathematics, not
389 geographically: the I<phi> zero is at the North Pole, not at the
390 Equator on the west coast of Africa (Bay of Guinea). You need to
391 subtract your geographical coordinates from I<pi/2> (also known as 90
394 $distance = great_circle_distance($lon0, pi/2 - $lat0,
395 $lon1, pi/2 - $lat1, $rho);
399 To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N
400 139.8E) in kilometers:
402 use Math::Trig qw(great_circle_distance deg2rad);
404 # Notice the 90 - latitude: phi zero is at the North Pole.
405 @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
406 @T = (deg2rad(139.8),deg2rad(90 - 35.7));
408 $km = great_circle_distance(@L, @T, 6378);
410 The answer may be off by few percentages because of the irregular
411 (slightly aspherical) form of the Earth. The used formula
413 lat0 = 90 degrees - phi0
414 lat1 = 90 degrees - phi1
415 d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
416 sin(lat0) * sin(lat1))
418 is also somewhat unreliable for small distances (for locations
419 separated less than about five degrees) because it uses arc cosine
420 which is rather ill-conditioned for values close to zero.
424 Saying C<use Math::Trig;> exports many mathematical routines in the
425 caller environment and even overrides some (C<sin>, C<cos>). This is
426 construed as a feature by the Authors, actually... ;-)
428 The code is not optimized for speed, especially because we use
429 C<Math::Complex> and thus go quite near complex numbers while doing
430 the computations even when the arguments are not. This, however,
431 cannot be completely avoided if we want things like C<asin(2)> to give
432 an answer instead of giving a fatal runtime error.
436 Jarkko Hietaniemi <F<jhi@iki.fi>> and
437 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.