This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math::Complex and Math::Trig updates (Re: [perl #37117] Math::Complex atan2 bug)
[perl5.git] / lib / Math / Trig.pm
1 #
2 # Trigonometric functions, mostly inherited from Math::Complex.
3 # -- Jarkko Hietaniemi, since April 1997
4 # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
5 #
6
7 require Exporter;
8 package Math::Trig;
9
10 use 5.006;
11 use strict;
12
13 use Math::Complex 1.35;
14 use Math::Complex qw(:trig);
15
16 our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
17
18 @ISA = qw(Exporter);
19
20 $VERSION = 1.03;
21
22 my @angcnv = qw(rad2deg rad2grad
23                 deg2rad deg2grad
24                 grad2rad grad2deg);
25
26 @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
27            @angcnv);
28
29 my @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 my @greatcircle = qw(
37                      great_circle_distance
38                      great_circle_direction
39                      great_circle_bearing
40                      great_circle_waypoint
41                      great_circle_midpoint
42                      great_circle_destination
43                     );
44
45 my @pi = qw(pi2 pip2 pip4);
46
47 @EXPORT_OK = (@rdlcnv, @greatcircle, @pi);
48
49 # See e.g. the following pages:
50 # http://www.movable-type.co.uk/scripts/LatLong.html
51 # http://williams.best.vwh.net/avform.htm
52
53 %EXPORT_TAGS = ('radial' => [ @rdlcnv ],
54                 'great_circle' => [ @greatcircle ],
55                 'pi'     => [ @pi ]);
56
57 sub pi2  () { 2 * pi }
58 sub pip2 () { pi / 2 }
59 sub pip4 () { pi / 4 }
60
61 sub DR  () { pi2/360 }
62 sub RD  () { 360/pi2 }
63 sub DG  () { 400/360 }
64 sub GD  () { 360/400 }
65 sub RG  () { 400/pi2 }
66 sub GR  () { pi2/400 }
67
68 #
69 # Truncating remainder.
70 #
71
72 sub remt ($$) {
73     # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
74     $_[0] - $_[1] * int($_[0] / $_[1]);
75 }
76
77 #
78 # Angle conversions.
79 #
80
81 sub rad2rad($)     { remt($_[0], pi2) }
82
83 sub deg2deg($)     { remt($_[0], 360) }
84
85 sub grad2grad($)   { remt($_[0], 400) }
86
87 sub rad2deg ($;$)  { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) }
88
89 sub deg2rad ($;$)  { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) }
90
91 sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) }
92
93 sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) }
94
95 sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) }
96
97 sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) }
98
99 sub cartesian_to_spherical {
100     my ( $x, $y, $z ) = @_;
101
102     my $rho = sqrt( $x * $x + $y * $y + $z * $z );
103
104     return ( $rho,
105              atan2( $y, $x ),
106              $rho ? acos( $z / $rho ) : 0 );
107 }
108
109 sub spherical_to_cartesian {
110     my ( $rho, $theta, $phi ) = @_;
111
112     return ( $rho * cos( $theta ) * sin( $phi ),
113              $rho * sin( $theta ) * sin( $phi ),
114              $rho * cos( $phi   ) );
115 }
116
117 sub spherical_to_cylindrical {
118     my ( $x, $y, $z ) = spherical_to_cartesian( @_ );
119
120     return ( sqrt( $x * $x + $y * $y ), $_[1], $z );
121 }
122
123 sub cartesian_to_cylindrical {
124     my ( $x, $y, $z ) = @_;
125
126     return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
127 }
128
129 sub cylindrical_to_cartesian {
130     my ( $rho, $theta, $z ) = @_;
131
132     return ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
133 }
134
135 sub cylindrical_to_spherical {
136     return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) );
137 }
138
139 sub great_circle_distance {
140     my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_;
141
142     $rho = 1 unless defined $rho; # Default to the unit sphere.
143
144     my $lat0 = pip2 - $phi0;
145     my $lat1 = pip2 - $phi1;
146
147     return $rho *
148         acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
149              sin( $lat0 ) * sin( $lat1 ) );
150 }
151
152 sub great_circle_direction {
153     my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
154
155     my $distance = &great_circle_distance;
156
157     my $lat0 = pip2 - $phi0;
158     my $lat1 = pip2 - $phi1;
159
160     my $direction =
161         acos((sin($lat1) - sin($lat0) * cos($distance)) /
162              (cos($lat0) * sin($distance)));
163
164     $direction = pi2 - $direction
165         if sin($theta1 - $theta0) < 0;
166
167     return rad2rad($direction);
168 }
169
170 *great_circle_bearing = \&great_circle_direction;
171
172 sub great_circle_waypoint {
173     my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
174
175     $point = 0.5 unless defined $point;
176
177     my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
178
179     return undef if $d == pi;
180
181     my $sd = sin($d);
182
183     return ($theta0, $phi0) if $sd == 0;
184
185     my $A = sin((1 - $point) * $d) / $sd;
186     my $B = sin(     $point  * $d) / $sd;
187
188     my $lat0 = pip2 - $phi0;
189     my $lat1 = pip2 - $phi1;
190
191     my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
192     my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
193     my $z = $A * sin($lat0)                + $B * sin($lat1);
194
195     my $theta = atan2($y, $x);
196     my $phi   = atan2($z, sqrt($x*$x + $y*$y));
197     
198     return ($theta, $phi);
199 }
200
201 sub great_circle_midpoint {
202     great_circle_waypoint(@_[0..3], 0.5);
203 }
204
205 sub great_circle_destination {
206     my ( $theta0, $phi0, $dir0, $dst ) = @_;
207
208     my $lat0 = pip2 - $phi0;
209
210     my $phi1   = asin(sin($lat0)*cos($dst)+cos($lat0)*sin($dst)*cos($dir0));
211     my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
212                                  cos($dst)-sin($lat0)*sin($phi1));
213
214     my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
215
216     $dir1 -= pi2 if $dir1 > pi2;
217
218     return ($theta1, $phi1, $dir1);
219 }
220
221 1;
222
223 __END__
224 =pod
225
226 =head1 NAME
227
228 Math::Trig - trigonometric functions
229
230 =head1 SYNOPSIS
231
232         use Math::Trig;
233
234         $x = tan(0.9);
235         $y = acos(3.7);
236         $z = asin(2.4);
237
238         $halfpi = pi/2;
239
240         $rad = deg2rad(120);
241
242         # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
243         use Math::Trig ':pi';
244
245         # Import the conversions between cartesian/spherical/cylindrical.
246         use Math::Trig ':radial';
247
248         # Import the great circle formulas.
249         use Math::Trig ':great_circle';
250
251 =head1 DESCRIPTION
252
253 C<Math::Trig> defines many trigonometric functions not defined by the
254 core Perl which defines only the C<sin()> and C<cos()>.  The constant
255 B<pi> is also defined as are a few convenience functions for angle
256 conversions, and I<great circle formulas> for spherical movement.
257
258 =head1 TRIGONOMETRIC FUNCTIONS
259
260 The tangent
261
262 =over 4
263
264 =item B<tan>
265
266 =back
267
268 The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
269 are aliases)
270
271 B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>
272
273 The arcus (also known as the inverse) functions of the sine, cosine,
274 and tangent
275
276 B<asin>, B<acos>, B<atan>
277
278 The principal value of the arc tangent of y/x
279
280 B<atan2>(y, x)
281
282 The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
283 and acotan/acot are aliases)
284
285 B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>
286
287 The hyperbolic sine, cosine, and tangent
288
289 B<sinh>, B<cosh>, B<tanh>
290
291 The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
292 and cotanh/coth are aliases)
293
294 B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>
295
296 The arcus (also known as the inverse) functions of the hyperbolic
297 sine, cosine, and tangent
298
299 B<asinh>, B<acosh>, B<atanh>
300
301 The arcus cofunctions of the hyperbolic sine, cosine, and tangent
302 (acsch/acosech and acoth/acotanh are aliases)
303
304 B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>
305
306 The trigonometric constant B<pi> is also defined.
307
308 $pi2 = 2 * B<pi>;
309
310 =head2 ERRORS DUE TO DIVISION BY ZERO
311
312 The following functions
313
314         acoth
315         acsc
316         acsch
317         asec
318         asech
319         atanh
320         cot
321         coth
322         csc
323         csch
324         sec
325         sech
326         tan
327         tanh
328
329 cannot be computed for all arguments because that would mean dividing
330 by zero or taking logarithm of zero. These situations cause fatal
331 runtime errors looking like this
332
333         cot(0): Division by zero.
334         (Because in the definition of cot(0), the divisor sin(0) is 0)
335         Died at ...
336
337 or
338
339         atanh(-1): Logarithm of zero.
340         Died at...
341
342 For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
343 C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
344 C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
345 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
346 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
347 pi>, where I<k> is any integer.  atan2(0, 0) is undefined.
348
349 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
350
351 Please note that some of the trigonometric functions can break out
352 from the B<real axis> into the B<complex plane>. For example
353 C<asin(2)> has no definition for plain real numbers but it has
354 definition for complex numbers.
355
356 In Perl terms this means that supplying the usual Perl numbers (also
357 known as scalars, please see L<perldata>) as input for the
358 trigonometric functions might produce as output results that no more
359 are simple real numbers: instead they are complex numbers.
360
361 The C<Math::Trig> handles this by using the C<Math::Complex> package
362 which knows how to handle complex numbers, please see L<Math::Complex>
363 for more information. In practice you need not to worry about getting
364 complex numbers as results because the C<Math::Complex> takes care of
365 details like for example how to display complex numbers. For example:
366
367         print asin(2), "\n";
368
369 should produce something like this (take or leave few last decimals):
370
371         1.5707963267949-1.31695789692482i
372
373 That is, a complex number with the real part of approximately C<1.571>
374 and the imaginary part of approximately C<-1.317>.
375
376 =head1 PLANE ANGLE CONVERSIONS
377
378 (Plane, 2-dimensional) angles may be converted with the following functions.
379
380         $radians  = deg2rad($degrees);
381         $radians  = grad2rad($gradians);
382
383         $degrees  = rad2deg($radians);
384         $degrees  = grad2deg($gradians);
385
386         $gradians = deg2grad($degrees);
387         $gradians = rad2grad($radians);
388
389 The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
390 The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
391 If you don't want this, supply a true second argument:
392
393         $zillions_of_radians  = deg2rad($zillions_of_degrees, 1);
394         $negative_degrees     = rad2deg($negative_radians, 1);
395
396 You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
397 grad2grad().
398
399 =head1 RADIAL COORDINATE CONVERSIONS
400
401 B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
402 systems, explained shortly in more detail.
403
404 You can import radial coordinate conversion functions by using the
405 C<:radial> tag:
406
407     use Math::Trig ':radial';
408
409     ($rho, $theta, $z)     = cartesian_to_cylindrical($x, $y, $z);
410     ($rho, $theta, $phi)   = cartesian_to_spherical($x, $y, $z);
411     ($x, $y, $z)           = cylindrical_to_cartesian($rho, $theta, $z);
412     ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
413     ($x, $y, $z)           = spherical_to_cartesian($rho, $theta, $phi);
414     ($rho_c, $theta, $z)   = spherical_to_cylindrical($rho_s, $theta, $phi);
415
416 B<All angles are in radians>.
417
418 =head2 COORDINATE SYSTEMS
419
420 B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
421
422 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
423 coordinates which define a point in three-dimensional space.  They are
424 based on a sphere surface.  The radius of the sphere is B<rho>, also
425 known as the I<radial> coordinate.  The angle in the I<xy>-plane
426 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
427 coordinate.  The angle from the I<z>-axis is B<phi>, also known as the
428 I<polar> coordinate.  The `North Pole' is therefore I<0, 0, rho>, and
429 the `Bay of Guinea' (think of the missing big chunk of Africa) I<0,
430 pi/2, rho>.  In geographical terms I<phi> is latitude (northward
431 positive, southward negative) and I<theta> is longitude (eastward
432 positive, westward negative).
433
434 B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
435 some texts define the I<phi> to start from the horizontal plane, some
436 texts use I<r> in place of I<rho>.
437
438 Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
439 coordinates which define a point in three-dimensional space.  They are
440 based on a cylinder surface.  The radius of the cylinder is B<rho>,
441 also known as the I<radial> coordinate.  The angle in the I<xy>-plane
442 (around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
443 coordinate.  The third coordinate is the I<z>, pointing up from the
444 B<theta>-plane.
445
446 =head2 3-D ANGLE CONVERSIONS
447
448 Conversions to and from spherical and cylindrical coordinates are
449 available.  Please notice that the conversions are not necessarily
450 reversible because of the equalities like I<pi> angles being equal to
451 I<-pi> angles.
452
453 =over 4
454
455 =item cartesian_to_cylindrical
456
457         ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
458
459 =item cartesian_to_spherical
460
461         ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
462
463 =item cylindrical_to_cartesian
464
465         ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
466
467 =item cylindrical_to_spherical
468
469         ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
470
471 Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
472
473 =item spherical_to_cartesian
474
475         ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
476
477 =item spherical_to_cylindrical
478
479         ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
480
481 Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
482
483 =back
484
485 =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
486
487 You can compute spherical distances, called B<great circle distances>,
488 by importing the great_circle_distance() function:
489
490   use Math::Trig 'great_circle_distance';
491
492   $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
493
494 The I<great circle distance> is the shortest distance between two
495 points on a sphere.  The distance is in C<$rho> units.  The C<$rho> is
496 optional, it defaults to 1 (the unit sphere), therefore the distance
497 defaults to radians.
498
499 If you think geographically the I<theta> are longitudes: zero at the
500 Greenwhich meridian, eastward positive, westward negative--and the
501 I<phi> are latitudes: zero at the North Pole, northward positive,
502 southward negative.  B<NOTE>: this formula thinks in mathematics, not
503 geographically: the I<phi> zero is at the North Pole, not at the
504 Equator on the west coast of Africa (Bay of Guinea).  You need to
505 subtract your geographical coordinates from I<pi/2> (also known as 90
506 degrees).
507
508   $distance = great_circle_distance($lon0, pi/2 - $lat0,
509                                     $lon1, pi/2 - $lat1, $rho);
510
511 The direction you must follow the great circle (also known as I<bearing>)
512 can be computed by the great_circle_direction() function:
513
514   use Math::Trig 'great_circle_direction';
515
516   $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
517
518 (Alias 'great_circle_bearing' is also available.)
519 The result is in radians, zero indicating straight north, pi or -pi
520 straight south, pi/2 straight west, and -pi/2 straight east.
521
522 You can inversely compute the destination if you know the
523 starting point, direction, and distance:
524
525   use Math::Trig 'great_circle_destination';
526
527   # thetad and phid are the destination coordinates,
528   # dird is the final direction at the destination.
529
530   ($thetad, $phid, $dird) =
531     great_circle_destination($theta, $phi, $direction, $distance);
532
533 or the midpoint if you know the end points:
534
535   use Math::Trig 'great_circle_midpoint';
536
537   ($thetam, $phim) =
538     great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
539
540 The great_circle_midpoint() is just a special case of
541
542   use Math::Trig 'great_circle_waypoint';
543
544   ($thetai, $phii) =
545     great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
546
547 Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
548 $phi1).  Note that antipodal points (where their distance is I<pi>
549 radians) do not have waypoints between them (they would have an an
550 "equator" between them), and therefore C<undef> is returned for
551 antipodal points.  If the points are the same and the distance
552 therefore zero and all waypoints therefore identical, the first point
553 (either point) is returned.
554
555 The thetas, phis, direction, and distance in the above are all in radians.
556
557 You can import all the great circle formulas by
558
559   use Math::Trig ':great_circle';
560
561 Notice that the resulting directions might be somewhat surprising if
562 you are looking at a flat worldmap: in such map projections the great
563 circles quite often do not look like the shortest routes-- but for
564 example the shortest possible routes from Europe or North America to
565 Asia do often cross the polar regions.
566
567 =head1 EXAMPLES
568
569 To calculate the distance between London (51.3N 0.5W) and Tokyo
570 (35.7N 139.8E) in kilometers:
571
572         use Math::Trig qw(great_circle_distance deg2rad);
573
574         # Notice the 90 - latitude: phi zero is at the North Pole.
575         sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
576         my @L = NESW( -0.5, 51.3);
577         my @T = NESW(139.8, 35.7);
578         my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
579
580 The direction you would have to go from London to Tokyo (in radians,
581 straight north being zero, straight east being pi/2).
582
583         use Math::Trig qw(great_circle_direction);
584
585         my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
586
587 The midpoint between London and Tokyo being
588
589         use Math::Trig qw(great_circle_midpoint);
590
591         my @M = great_circle_midpoint(@L, @T);
592
593 or about 68.11N 24.74E, in the Finnish Lapland.
594
595 =head2 CAVEAT FOR GREAT CIRCLE FORMULAS
596
597 The answers may be off by few percentages because of the irregular
598 (slightly aspherical) form of the Earth.  The errors are at worst
599 about 0.55%, but generally below 0.3%.
600
601 =head1 BUGS
602
603 Saying C<use Math::Trig;> exports many mathematical routines in the
604 caller environment and even overrides some (C<sin>, C<cos>).  This is
605 construed as a feature by the Authors, actually... ;-)
606
607 The code is not optimized for speed, especially because we use
608 C<Math::Complex> and thus go quite near complex numbers while doing
609 the computations even when the arguments are not. This, however,
610 cannot be completely avoided if we want things like C<asin(2)> to give
611 an answer instead of giving a fatal runtime error.
612
613 Do not attempt navigation using these formulas.
614
615 =head1 AUTHORS
616
617 Jarkko Hietaniemi <F<jhi@iki.fi>> and 
618 Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.
619
620 =cut
621
622 # eof