Re: Copious warnings from Sys::Syslog
[p5sagit/p5-mst-13.2.git] / lib / Math / Trig.pm
1 #
2 # Trigonometric functions, mostly inherited from Math::Complex.
3 # -- Jarkko Hietaniemi, since April 1997
4 # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
5 #
6
7 require Exporter;
8 package Math::Trig;
9
10 use 5.006;
11 use strict;
12
13 use Math::Complex qw(:trig);
14
15 our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
16
17 @ISA = qw(Exporter);
18
19 $VERSION = 1.01;
20
21 my @angcnv = qw(rad2deg rad2grad
22              deg2rad deg2grad
23              grad2rad grad2deg);
24
25 @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
26            @angcnv);
27
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);
34
35 @EXPORT_OK = (@rdlcnv, 'great_circle_distance', 'great_circle_direction');
36
37 %EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
38
39 sub pi2  () { 2 * pi }
40 sub pip2 () { pi / 2 }
41
42 sub DR  () { pi2/360 }
43 sub RD  () { 360/pi2 }
44 sub DG  () { 400/360 }
45 sub GD  () { 360/400 }
46 sub RG  () { 400/pi2 }
47 sub GR  () { pi2/400 }
48
49 #
50 # Truncating remainder.
51 #
52
53 sub remt ($$) {
54     # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
55     $_[0] - $_[1] * int($_[0] / $_[1]);
56 }
57
58 #
59 # Angle conversions.
60 #
61
62 sub rad2rad($)     { remt($_[0], pi2) }
63
64 sub deg2deg($)     { remt($_[0], 360) }
65
66 sub grad2grad($)   { remt($_[0], 400) }
67
68 sub rad2deg ($;$)  { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) }
69
70 sub deg2rad ($;$)  { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) }
71
72 sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) }
73
74 sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) }
75
76 sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) }
77
78 sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) }
79
80 sub cartesian_to_spherical {
81     my ( $x, $y, $z ) = @_;
82
83     my $rho = sqrt( $x * $x + $y * $y + $z * $z );
84
85     return ( $rho,
86              atan2( $y, $x ),
87              $rho ? acos( $z / $rho ) : 0 );
88 }
89
90 sub spherical_to_cartesian {
91     my ( $rho, $theta, $phi ) = @_;
92
93     return ( $rho * cos( $theta ) * sin( $phi ),
94              $rho * sin( $theta ) * sin( $phi ),
95              $rho * cos( $phi   ) );
96 }
97
98 sub spherical_to_cylindrical {
99     my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
100
101     return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
102 }
103
104 sub cartesian_to_cylindrical {
105     my ( $x, $y, $z ) = @_;
106
107     return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
108 }
109
110 sub cylindrical_to_cartesian {
111     my ( $rho, $theta, $z ) = @_;
112
113     return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
114 }
115
116 sub cylindrical_to_spherical {
117     return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
118 }
119
120 sub great_circle_distance {
121     my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
122
123     $rho = 1 unless defined $rho; # Default to the unit sphere.
124
125     my $lat0 = pip2 - $phi0;
126     my $lat1 = pip2 - $phi1;
127
128     return $rho *
129         acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
130              sin( $lat0 ) * sin( $lat1 ) );
131 }
132
133 sub great_circle_direction {
134     my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
135
136     my $lat0 = pip2 - $phi0;
137     my $lat1 = pip2 - $phi1;
138
139     my $direction =
140         atan2(sin($theta0 - $theta1) * cos($lat1),
141               cos($lat0) * sin($lat1) -
142               sin($lat0) * cos($lat1) * cos($theta0 - $theta1));
143
144     return rad2rad($direction);
145 }
146
147 =pod
148
149 =head1 NAME
150
151 Math::Trig - trigonometric functions
152
153 =head1 SYNOPSIS
154
155         use Math::Trig;
156
157         $x = tan(0.9);
158         $y = acos(3.7);
159         $z = asin(2.4);
160
161         $halfpi = pi/2;
162
163         $rad = deg2rad(120);
164
165 =head1 DESCRIPTION
166
167 C<Math::Trig> defines many trigonometric functions not defined by the
168 core Perl which defines only the C<sin()> and C<cos()>.  The constant
169 B<pi> is also defined as are a few convenience functions for angle
170 conversions.
171
172 =head1 TRIGONOMETRIC FUNCTIONS
173
174 The tangent
175
176 =over 4
177
178 =item B<tan>
179
180 =back
181
182 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
183 are aliases)
184
185 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
186
187 The arcus (also known as the inverse) functions of the sine, cosine,
188 and tangent
189
190 B<asin>, B<acos>, B<atan>
191
192 The principal value of the arc tangent of y/x
193
194 B<atan2>(y, x)
195
196 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
197 and acotan/acot are aliases)
198
199 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
200
201 The hyperbolic sine, cosine, and tangent
202
203 B<sinh>, B<cosh>, B<tanh>
204
205 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
206 and cotanh/coth are aliases)
207
208 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
209
210 The arcus (also known as the inverse) functions of the hyperbolic
211 sine, cosine, and tangent
212
213 B<asinh>, B<acosh>, B<atanh>
214
215 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
216 (acsch/acosech and acoth/acotanh are aliases)
217
218 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
219
220 The trigonometric constant B<pi> is also defined.
221
222 $pi2 = 2 * B<pi>;
223
224 =head2 ERRORS DUE TO DIVISION BY ZERO
225
226 The following functions
227
228         acoth
229         acsc
230         acsch
231         asec
232         asech
233         atanh
234         cot
235         coth
236         csc
237         csch
238         sec
239         sech
240         tan
241         tanh
242
243 cannot be computed for all arguments because that would mean dividing
244 by zero or taking logarithm of zero. These situations cause fatal
245 runtime errors looking like this
246
247         cot(0): Division by zero.
248         (Because in the definition of cot(0), the divisor sin(0) is 0)
249         Died at ...
250
251 or
252
253         atanh(-1): Logarithm of zero.
254         Died at...
255
256 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
257 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
258 C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
259 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
260 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
261 pi>, where I<k> is any integer.
262
263 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
264
265 Please note that some of the trigonometric functions can break out
266 from the B<real axis> into the B<complex plane>. For example
267 C<asin(2)> has no definition for plain real numbers but it has
268 definition for complex numbers.
269
270 In Perl terms this means that supplying the usual Perl numbers (also
271 known as scalars, please see L<perldata>) as input for the
272 trigonometric functions might produce as output results that no more
273 are simple real numbers: instead they are complex numbers.
274
275 The C<Math::Trig> handles this by using the C<Math::Complex> package
276 which knows how to handle complex numbers, please see L<Math::Complex>
277 for more information. In practice you need not to worry about getting
278 complex numbers as results because the C<Math::Complex> takes care of
279 details like for example how to display complex numbers. For example:
280
281         print asin(2), "\n";
282
283 should produce something like this (take or leave few last decimals):
284
285         1.5707963267949-1.31695789692482i
286
287 That is, a complex number with the real part of approximately C<1.571>
288 and the imaginary part of approximately C<-1.317>.
289
290 =head1 PLANE ANGLE CONVERSIONS
291
292 (Plane, 2-dimensional) angles may be converted with the following functions.
293
294         $radians  = deg2rad($degrees);
295         $radians  = grad2rad($gradians);
296
297         $degrees  = rad2deg($radians);
298         $degrees  = grad2deg($gradians);
299
300         $gradians = deg2grad($degrees);
301         $gradians = rad2grad($radians);
302
303 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
304 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
305 If you don't want this, supply a true second argument:
306
307         $zillions_of_radians  = deg2rad($zillions_of_degrees, 1);
308         $negative_degrees     = rad2deg($negative_radians, 1);
309
310 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
311 grad2grad().
312
313 =head1 RADIAL COORDINATE CONVERSIONS
314
315 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
316 systems, explained shortly in more detail.
317
318 You can import radial coordinate conversion functions by using the
319 C<:radial> tag:
320
321     use Math::Trig ':radial';
322
323     ($rho, $theta, $z)     = cartesian_to_cylindrical($x, $y, $z);
324     ($rho, $theta, $phi)   = cartesian_to_spherical($x, $y, $z);
325     ($x, $y, $z)           = cylindrical_to_cartesian($rho, $theta, $z);
326     ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
327     ($x, $y, $z)           = spherical_to_cartesian($rho, $theta, $phi);
328     ($rho_c, $theta, $z)   = spherical_to_cylindrical($rho_s, $theta, $phi);
329
330 B<All angles are in radians>.
331
332 =head2 COORDINATE SYSTEMS
333
334 B<Cartesian> coordinates are the usual rectangular I<(x, y,
335 z)>-coordinates.
336
337 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
338 coordinates which define a point in three-dimensional space.  They are
339 based on a sphere surface.  The radius of the sphere is B<rho>, also
340 known as the I<radial> coordinate.  The angle in the I<xy>-plane
341 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
342 coordinate.  The angle from the I<z>-axis is B<phi>, also known as the
343 I<polar> coordinate.  The `North Pole' is therefore I<0, 0, rho>, and
344 the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
345 pi/2, rho>.  In geographical terms I<phi> is latitude (northward
346 positive, southward negative) and I<theta> is longitude (eastward
347 positive, westward negative).
348
349 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
350 some texts define the I<phi> to start from the horizontal plane, some
351 texts use I<r> in place of I<rho>.
352
353 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
354 coordinates which define a point in three-dimensional space.  They are
355 based on a cylinder surface.  The radius of the cylinder is B<rho>,
356 also known as the I<radial> coordinate.  The angle in the I<xy>-plane
357 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
358 coordinate.  The third coordinate is the I<z>, pointing up from the
359 B<theta>-plane.
360
361 =head2 3-D ANGLE CONVERSIONS
362
363 Conversions to and from spherical and cylindrical coordinates are
364 available.  Please notice that the conversions are not necessarily
365 reversible because of the equalities like I<pi> angles being equal to
366 I<-pi> angles.
367
368 =over 4
369
370 =item cartesian_to_cylindrical
371
372         ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
373
374 =item cartesian_to_spherical
375
376         ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
377
378 =item cylindrical_to_cartesian
379
380         ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
381
382 =item cylindrical_to_spherical
383
384         ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
385
386 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
387
388 =item spherical_to_cartesian
389
390         ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
391
392 =item spherical_to_cylindrical
393
394         ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
395
396 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
397
398 =back
399
400 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
401
402 You can compute spherical distances, called B<great circle distances>,
403 by importing the great_circle_distance() function:
404
405   use Math::Trig 'great_circle_distance';
406
407   $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
408
409 The I<great circle distance> is the shortest distance between two
410 points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
411 optional, it defaults to 1 (the unit sphere), therefore the distance
412 defaults to radians.
413
414 If you think geographically the I<theta> are longitudes: zero at the
415 Greenwhich meridian, eastward positive, westward negative--and the
416 I<phi> are latitudes: zero at the North Pole, northward positive,
417 southward negative.  B<NOTE>: this formula thinks in mathematics, not
418 geographically: the I<phi> zero is at the North Pole, not at the
419 Equator on the west coast of Africa (Bay of Guinea).  You need to
420 subtract your geographical coordinates from I<pi/2> (also known as 90
421 degrees).
422
423   $distance = great_circle_distance($lon0, pi/2 - $lat0,
424                                     $lon1, pi/2 - $lat1, $rho);
425
426 The direction you must follow the great circle can be computed by the
427 great_circle_direction() function:
428
429   use Math::Trig 'great_circle_direction';
430
431   $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
432
433 The result is in radians, zero indicating straight north, pi or -pi
434 straight south, pi/2 straight west, and -pi/2 straight east.
435
436 Notice that the resulting directions might be somewhat surprising if
437 you are looking at a flat worldmap: in such map projections the great
438 circles quite often do not look like the shortest routes-- but for
439 example the shortest possible routes from Europe or North America to
440 Asia do often cross the polar regions.
441
442 =head1 EXAMPLES
443
444 To calculate the distance between London (51.3N 0.5W) and Tokyo
445 (35.7N 139.8E) in kilometers:
446
447         use Math::Trig qw(great_circle_distance deg2rad);
448
449         # Notice the 90 - latitude: phi zero is at the North Pole.
450         @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
451         @T = (deg2rad(139.8),deg2rad(90 - 35.7));
452
453         $km = great_circle_distance(@L, @T, 6378);
454
455 The direction you would have to go from London to Tokyo
456
457         use Math::Trig qw(great_circle_direction);
458
459         $rad = great_circle_direction(@L, @T);
460
461 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
462
463 The answers may be off by few percentages because of the irregular
464 (slightly aspherical) form of the Earth.  The formula used for
465 grear circle distances
466
467         lat0 = 90 degrees - phi0
468         lat1 = 90 degrees - phi1
469         d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
470                        sin(lat0) * sin(lat1))
471
472 is also somewhat unreliable for small distances (for locations
473 separated less than about five degrees) because it uses arc cosine
474 which is rather ill-conditioned for values close to zero.
475
476 =head1 BUGS
477
478 Saying C<use Math::Trig;> exports many mathematical routines in the
479 caller environment and even overrides some (C<sin>, C<cos>).  This is
480 construed as a feature by the Authors, actually... ;-)
481
482 The code is not optimized for speed, especially because we use
483 C<Math::Complex> and thus go quite near complex numbers while doing
484 the computations even when the arguments are not. This, however,
485 cannot be completely avoided if we want things like C<asin(2)> to give
486 an answer instead of giving a fatal runtime error.
487
488 =head1 AUTHORS
489
490 Jarkko Hietaniemi <F<jhi@iki.fi>> and 
491 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.
492
493 =cut
494
495 # eof