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