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