This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Unicode ReadMe update for Unicode 3.1.
[perl5.git] / lib / Math / Trig.pm
CommitLineData
5aabfad6
PP
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
PP
5#
6
7require Exporter;
8package Math::Trig;
9
17f410f9 10use 5.005_64;
5aabfad6
PP
11use strict;
12
13use Math::Complex qw(:trig);
14
17f410f9 15our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
5aabfad6
PP
16
17@ISA = qw(Exporter);
18
19$VERSION = 1.00;
20
ace5de91
GS
21my @angcnv = qw(rad2deg rad2grad
22 deg2rad deg2grad
23 grad2rad grad2deg);
5aabfad6
PP
24
25@EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
26 @angcnv);
27
d54bf66f
JH
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
35@EXPORT_OK = (@rdlcnv, 'great_circle_distance');
36
37%EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
38
9db5a202
JH
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
PP
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
JH
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
JH
76sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) }
77
78sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) }
5aabfad6 79
d54bf66f
JH
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
133=pod
134
5aabfad6
PP
135=head1 NAME
136
137Math::Trig - trigonometric functions
138
139=head1 SYNOPSIS
140
141 use Math::Trig;
3cb6de81 142
5aabfad6
PP
143 $x = tan(0.9);
144 $y = acos(3.7);
145 $z = asin(2.4);
3cb6de81 146
5aabfad6
PP
147 $halfpi = pi/2;
148
ace5de91 149 $rad = deg2rad(120);
5aabfad6
PP
150
151=head1 DESCRIPTION
152
153C<Math::Trig> defines many trigonometric functions not defined by the
4ae80833 154core Perl which defines only the C<sin()> and C<cos()>. The constant
5aabfad6
PP
155B<pi> is also defined as are a few convenience functions for angle
156conversions.
157
158=head1 TRIGONOMETRIC FUNCTIONS
159
160The tangent
161
d54bf66f
JH
162=over 4
163
164=item B<tan>
165
166=back
5aabfad6
PP
167
168The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
169are aliases)
170
d54bf66f 171B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
5aabfad6
PP
172
173The arcus (also known as the inverse) functions of the sine, cosine,
174and tangent
175
d54bf66f 176B<asin>, B<acos>, B<atan>
5aabfad6
PP
177
178The principal value of the arc tangent of y/x
179
d54bf66f 180B<atan2>(y, x)
5aabfad6
PP
181
182The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
183and acotan/acot are aliases)
184
d54bf66f 185B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
5aabfad6
PP
186
187The hyperbolic sine, cosine, and tangent
188
d54bf66f 189B<sinh>, B<cosh>, B<tanh>
5aabfad6
PP
190
191The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
192and cotanh/coth are aliases)
193
d54bf66f 194B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
5aabfad6
PP
195
196The arcus (also known as the inverse) functions of the hyperbolic
197sine, cosine, and tangent
198
d54bf66f 199B<asinh>, B<acosh>, B<atanh>
5aabfad6
PP
200
201The arcus cofunctions of the hyperbolic sine, cosine, and tangent
202(acsch/acosech and acoth/acotanh are aliases)
203
d54bf66f 204B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
5aabfad6
PP
205
206The trigonometric constant B<pi> is also defined.
207
d54bf66f 208$pi2 = 2 * B<pi>;
5aabfad6 209
5cd24f17
PP
210=head2 ERRORS DUE TO DIVISION BY ZERO
211
212The following functions
213
d54bf66f 214 acoth
5cd24f17 215 acsc
5cd24f17 216 acsch
d54bf66f
JH
217 asec
218 asech
219 atanh
220 cot
221 coth
222 csc
223 csch
224 sec
225 sech
226 tan
227 tanh
5cd24f17
PP
228
229cannot be computed for all arguments because that would mean dividing
8c03c583
JH
230by zero or taking logarithm of zero. These situations cause fatal
231runtime errors looking like this
5cd24f17
PP
232
233 cot(0): Division by zero.
234 (Because in the definition of cot(0), the divisor sin(0) is 0)
235 Died at ...
236
8c03c583
JH
237or
238
239 atanh(-1): Logarithm of zero.
240 Died at...
241
242For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
243C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the
244C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the
245C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the
246C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
247pi>, where I<k> is any integer.
5cd24f17
PP
248
249=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
5aabfad6
PP
250
251Please note that some of the trigonometric functions can break out
252from the B<real axis> into the B<complex plane>. For example
253C<asin(2)> has no definition for plain real numbers but it has
254definition for complex numbers.
255
256In Perl terms this means that supplying the usual Perl numbers (also
257known as scalars, please see L<perldata>) as input for the
258trigonometric functions might produce as output results that no more
259are simple real numbers: instead they are complex numbers.
260
261The C<Math::Trig> handles this by using the C<Math::Complex> package
262which knows how to handle complex numbers, please see L<Math::Complex>
263for more information. In practice you need not to worry about getting
264complex numbers as results because the C<Math::Complex> takes care of
265details like for example how to display complex numbers. For example:
266
267 print asin(2), "\n";
3cb6de81 268
5aabfad6
PP
269should produce something like this (take or leave few last decimals):
270
271 1.5707963267949-1.31695789692482i
272
5cd24f17
PP
273That is, a complex number with the real part of approximately C<1.571>
274and the imaginary part of approximately C<-1.317>.
5aabfad6 275
d54bf66f 276=head1 PLANE ANGLE CONVERSIONS
5aabfad6
PP
277
278(Plane, 2-dimensional) angles may be converted with the following functions.
279
ace5de91
GS
280 $radians = deg2rad($degrees);
281 $radians = grad2rad($gradians);
3cb6de81 282
ace5de91
GS
283 $degrees = rad2deg($radians);
284 $degrees = grad2deg($gradians);
3cb6de81 285
ace5de91
GS
286 $gradians = deg2grad($degrees);
287 $gradians = rad2grad($radians);
5aabfad6 288
5cd24f17 289The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
9db5a202
JH
290The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
291If you don't want this, supply a true second argument:
292
293 $zillions_of_radians = deg2rad($zillions_of_degrees, 1);
294 $negative_degrees = rad2deg($negative_radians, 1);
295
296You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
297grad2grad().
5aabfad6 298
d54bf66f
JH
299=head1 RADIAL COORDINATE CONVERSIONS
300
301B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
302systems, explained shortly in more detail.
303
304You can import radial coordinate conversion functions by using the
305C<:radial> tag:
306
307 use Math::Trig ':radial';
308
309 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
310 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
311 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
312 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
313 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
314 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
315
316B<All angles are in radians>.
317
318=head2 COORDINATE SYSTEMS
319
320B<Cartesian> coordinates are the usual rectangular I<(x, y,
321z)>-coordinates.
322
323Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
324coordinates which define a point in three-dimensional space. They are
325based on a sphere surface. The radius of the sphere is B<rho>, also
326known as the I<radial> coordinate. The angle in the I<xy>-plane
327(around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
328coordinate. The angle from the I<z>-axis is B<phi>, also known as the
329I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and
330the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
4b0d1da8
JH
331pi/2, rho>. In geographical terms I<phi> is latitude (northward
332positive, southward negative) and I<theta> is longitude (eastward
333positive, westward negative).
d54bf66f 334
4b0d1da8 335B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
d54bf66f
JH
336some texts define the I<phi> to start from the horizontal plane, some
337texts use I<r> in place of I<rho>.
338
339Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
340coordinates which define a point in three-dimensional space. They are
341based on a cylinder surface. The radius of the cylinder is B<rho>,
342also known as the I<radial> coordinate. The angle in the I<xy>-plane
343(around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
344coordinate. The third coordinate is the I<z>, pointing up from the
345B<theta>-plane.
346
347=head2 3-D ANGLE CONVERSIONS
348
349Conversions to and from spherical and cylindrical coordinates are
350available. Please notice that the conversions are not necessarily
351reversible because of the equalities like I<pi> angles being equal to
352I<-pi> angles.
353
354=over 4
355
356=item cartesian_to_cylindrical
357
358 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
359
360=item cartesian_to_spherical
361
362 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
363
364=item cylindrical_to_cartesian
365
366 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
367
368=item cylindrical_to_spherical
369
370 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
371
372Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
373
374=item spherical_to_cartesian
375
376 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
377
378=item spherical_to_cylindrical
379
380 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
381
382Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
383
384=back
385
386=head1 GREAT CIRCLE DISTANCES
387
388You can compute spherical distances, called B<great circle distances>,
389by importing the C<great_circle_distance> function:
390
391 use Math::Trig 'great_circle_distance'
392
4b0d1da8 393 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
d54bf66f
JH
394
395The I<great circle distance> is the shortest distance between two
396points on a sphere. The distance is in C<$rho> units. The C<$rho> is
397optional, it defaults to 1 (the unit sphere), therefore the distance
398defaults to radians.
399
4b0d1da8
JH
400If you think geographically the I<theta> are longitudes: zero at the
401Greenwhich meridian, eastward positive, westward negative--and the
2d06e7d7 402I<phi> are latitudes: zero at the North Pole, northward positive,
4b0d1da8 403southward negative. B<NOTE>: this formula thinks in mathematics, not
2d06e7d7
JH
404geographically: the I<phi> zero is at the North Pole, not at the
405Equator on the west coast of Africa (Bay of Guinea). You need to
406subtract your geographical coordinates from I<pi/2> (also known as 90
407degrees).
4b0d1da8
JH
408
409 $distance = great_circle_distance($lon0, pi/2 - $lat0,
410 $lon1, pi/2 - $lat1, $rho);
411
51301382 412=head1 EXAMPLES
d54bf66f
JH
413
414To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N
415139.8E) in kilometers:
416
417 use Math::Trig qw(great_circle_distance deg2rad);
418
419 # Notice the 90 - latitude: phi zero is at the North Pole.
420 @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
421 @T = (deg2rad(139.8),deg2rad(90 - 35.7));
422
423 $km = great_circle_distance(@L, @T, 6378);
424
4b0d1da8 425The answer may be off by few percentages because of the irregular
41bd693c
JH
426(slightly aspherical) form of the Earth. The used formula
427
428 lat0 = 90 degrees - phi0
429 lat1 = 90 degrees - phi1
430 d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
431 sin(lat0) * sin(lat1))
432
433is also somewhat unreliable for small distances (for locations
434separated less than about five degrees) because it uses arc cosine
435which is rather ill-conditioned for values close to zero.
d54bf66f 436
5cd24f17 437=head1 BUGS
5aabfad6 438
5cd24f17
PP
439Saying C<use Math::Trig;> exports many mathematical routines in the
440caller environment and even overrides some (C<sin>, C<cos>). This is
441construed as a feature by the Authors, actually... ;-)
5aabfad6 442
5cd24f17
PP
443The code is not optimized for speed, especially because we use
444C<Math::Complex> and thus go quite near complex numbers while doing
445the computations even when the arguments are not. This, however,
446cannot be completely avoided if we want things like C<asin(2)> to give
447an answer instead of giving a fatal runtime error.
5aabfad6 448
5cd24f17 449=head1 AUTHORS
5aabfad6 450
ace5de91 451Jarkko Hietaniemi <F<jhi@iki.fi>> and
6e238990 452Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.
5aabfad6
PP
453
454=cut
455
456# eof