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