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