This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
t/op/sysio.t for MPE/iX
[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
10use strict;
11
12use Math::Complex qw(:trig);
13
14use vars qw($VERSION $PACKAGE
15 @ISA
d54bf66f 16 @EXPORT @EXPORT_OK %EXPORT_TAGS);
5aabfad6
PP
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
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
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
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
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
PP
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
PP
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
PP
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
PP
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
PP
166
167The arcus (also known as the inverse) functions of the sine, cosine,
168and tangent
169
d54bf66f 170B<asin>, B<acos>, B<atan>
5aabfad6
PP
171
172The principal value of the arc tangent of y/x
173
d54bf66f 174B<atan2>(y, x)
5aabfad6
PP
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
PP
180
181The hyperbolic sine, cosine, and tangent
182
d54bf66f 183B<sinh>, B<cosh>, B<tanh>
5aabfad6
PP
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
PP
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
PP
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
PP
199
200The trigonometric constant B<pi> is also defined.
201
d54bf66f 202$pi2 = 2 * B<pi>;
5aabfad6 203
5cd24f17
PP
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
PP
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
PP
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
PP
242
243=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
5aabfad6
PP
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
PP
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
PP
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,
317pi/2, rho>.
318
319B<Beware>: some texts define I<theta> and I<phi> the other way round,
320some texts define the I<phi> to start from the horizontal plane, some
321texts use I<r> in place of I<rho>.
322
323Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional
324coordinates which define a point in three-dimensional space. They are
325based on a cylinder surface. The radius of the cylinder is B<rho>,
326also known as the I<radial> coordinate. The angle in the I<xy>-plane
327(around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
328coordinate. The third coordinate is the I<z>, pointing up from the
329B<theta>-plane.
330
331=head2 3-D ANGLE CONVERSIONS
332
333Conversions to and from spherical and cylindrical coordinates are
334available. Please notice that the conversions are not necessarily
335reversible because of the equalities like I<pi> angles being equal to
336I<-pi> angles.
337
338=over 4
339
340=item cartesian_to_cylindrical
341
342 ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z);
343
344=item cartesian_to_spherical
345
346 ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z);
347
348=item cylindrical_to_cartesian
349
350 ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z);
351
352=item cylindrical_to_spherical
353
354 ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z);
355
356Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>.
357
358=item spherical_to_cartesian
359
360 ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi);
361
362=item spherical_to_cylindrical
363
364 ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi);
365
366Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
367
368=back
369
370=head1 GREAT CIRCLE DISTANCES
371
372You can compute spherical distances, called B<great circle distances>,
373by importing the C<great_circle_distance> function:
374
375 use Math::Trig 'great_circle_distance'
376
377 $distance = great_circle_distance($theta0, $phi0, $theta1, $phi, [, $rho]);
378
379The I<great circle distance> is the shortest distance between two
380points on a sphere. The distance is in C<$rho> units. The C<$rho> is
381optional, it defaults to 1 (the unit sphere), therefore the distance
382defaults to radians.
383
51301382 384=head1 EXAMPLES
d54bf66f
JH
385
386To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N
387139.8E) in kilometers:
388
389 use Math::Trig qw(great_circle_distance deg2rad);
390
391 # Notice the 90 - latitude: phi zero is at the North Pole.
392 @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
393 @T = (deg2rad(139.8),deg2rad(90 - 35.7));
394
395 $km = great_circle_distance(@L, @T, 6378);
396
397The answer may be off by up to 0.3% because of the irregular (slightly
398aspherical) form of the Earth.
399
5cd24f17 400=head1 BUGS
5aabfad6 401
5cd24f17
PP
402Saying C<use Math::Trig;> exports many mathematical routines in the
403caller environment and even overrides some (C<sin>, C<cos>). This is
404construed as a feature by the Authors, actually... ;-)
5aabfad6 405
5cd24f17
PP
406The code is not optimized for speed, especially because we use
407C<Math::Complex> and thus go quite near complex numbers while doing
408the computations even when the arguments are not. This, however,
409cannot be completely avoided if we want things like C<asin(2)> to give
410an answer instead of giving a fatal runtime error.
5aabfad6 411
5cd24f17 412=head1 AUTHORS
5aabfad6 413
ace5de91
GS
414Jarkko Hietaniemi <F<jhi@iki.fi>> and
415Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>.
5aabfad6
PP
416
417=cut
418
419# eof