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