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