Commit | Line | Data |
5aabfad6 |
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 |
5 | # |
6 | |
7 | require Exporter; |
8 | package Math::Trig; |
9 | |
3b825e41 |
10 | use 5.006; |
5aabfad6 |
11 | use strict; |
12 | |
13 | use Math::Complex qw(:trig); |
14 | |
17f410f9 |
15 | our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS); |
5aabfad6 |
16 | |
17 | @ISA = qw(Exporter); |
18 | |
ea0630ea |
19 | $VERSION = 1.02; |
5aabfad6 |
20 | |
ace5de91 |
21 | my @angcnv = qw(rad2deg rad2grad |
d139edd6 |
22 | deg2rad deg2grad |
23 | grad2rad grad2deg); |
5aabfad6 |
24 | |
25 | @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, |
26 | @angcnv); |
27 | |
d54bf66f |
28 | my @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 |
36 | |
37 | %EXPORT_TAGS = ('radial' => [ @rdlcnv ]); |
38 | |
9db5a202 |
39 | sub pi2 () { 2 * pi } |
40 | sub pip2 () { pi / 2 } |
41 | |
42 | sub DR () { pi2/360 } |
43 | sub RD () { 360/pi2 } |
44 | sub DG () { 400/360 } |
45 | sub GD () { 360/400 } |
46 | sub RG () { 400/pi2 } |
47 | sub GR () { pi2/400 } |
5aabfad6 |
48 | |
49 | # |
50 | # Truncating remainder. |
51 | # |
52 | |
53 | sub 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 |
62 | sub rad2rad($) { remt($_[0], pi2) } |
63 | |
64 | sub deg2deg($) { remt($_[0], 360) } |
65 | |
66 | sub grad2grad($) { remt($_[0], 400) } |
5aabfad6 |
67 | |
9db5a202 |
68 | sub rad2deg ($;$) { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) } |
5aabfad6 |
69 | |
9db5a202 |
70 | sub deg2rad ($;$) { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) } |
5aabfad6 |
71 | |
9db5a202 |
72 | sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) } |
5aabfad6 |
73 | |
9db5a202 |
74 | sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) } |
5aabfad6 |
75 | |
9db5a202 |
76 | sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) } |
77 | |
78 | sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) } |
5aabfad6 |
79 | |
d54bf66f |
80 | sub 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 | |
90 | sub 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 | |
98 | sub spherical_to_cylindrical { |
99 | my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); |
100 | |
101 | return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); |
102 | } |
103 | |
104 | sub cartesian_to_cylindrical { |
105 | my ( $x, $y, $z ) = @_; |
106 | |
107 | return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); |
108 | } |
109 | |
110 | sub cylindrical_to_cartesian { |
111 | my ( $rho, $theta, $z ) = @_; |
112 | |
113 | return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); |
114 | } |
115 | |
116 | sub cylindrical_to_spherical { |
117 | return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); |
118 | } |
119 | |
120 | sub 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 |
133 | sub great_circle_direction { |
134 | my ( $theta0, $phi0, $theta1, $phi1 ) = @_; |
135 | |
d139edd6 |
136 | my $distance = &great_circle_distance; |
137 | |
7e5f197a |
138 | my $lat0 = pip2 - $phi0; |
139 | my $lat1 = pip2 - $phi1; |
140 | |
141 | my $direction = |
d139edd6 |
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 |
147 | |
148 | return rad2rad($direction); |
149 | } |
150 | |
ea0630ea |
151 | 1; |
152 | |
153 | __END__ |
d54bf66f |
154 | =pod |
155 | |
5aabfad6 |
156 | =head1 NAME |
157 | |
158 | Math::Trig - trigonometric functions |
159 | |
160 | =head1 SYNOPSIS |
161 | |
162 | use Math::Trig; |
3cb6de81 |
163 | |
5aabfad6 |
164 | $x = tan(0.9); |
165 | $y = acos(3.7); |
166 | $z = asin(2.4); |
3cb6de81 |
167 | |
5aabfad6 |
168 | $halfpi = pi/2; |
169 | |
ace5de91 |
170 | $rad = deg2rad(120); |
5aabfad6 |
171 | |
172 | =head1 DESCRIPTION |
173 | |
174 | C<Math::Trig> defines many trigonometric functions not defined by the |
4ae80833 |
175 | core Perl which defines only the C<sin()> and C<cos()>. The constant |
5aabfad6 |
176 | B<pi> is also defined as are a few convenience functions for angle |
177 | conversions. |
178 | |
179 | =head1 TRIGONOMETRIC FUNCTIONS |
180 | |
181 | The tangent |
182 | |
d54bf66f |
183 | =over 4 |
184 | |
185 | =item B<tan> |
186 | |
187 | =back |
5aabfad6 |
188 | |
189 | The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot |
190 | are aliases) |
191 | |
d54bf66f |
192 | B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan> |
5aabfad6 |
193 | |
194 | The arcus (also known as the inverse) functions of the sine, cosine, |
195 | and tangent |
196 | |
d54bf66f |
197 | B<asin>, B<acos>, B<atan> |
5aabfad6 |
198 | |
199 | The principal value of the arc tangent of y/x |
200 | |
d54bf66f |
201 | B<atan2>(y, x) |
5aabfad6 |
202 | |
203 | The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc |
204 | and acotan/acot are aliases) |
205 | |
d54bf66f |
206 | B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan> |
5aabfad6 |
207 | |
208 | The hyperbolic sine, cosine, and tangent |
209 | |
d54bf66f |
210 | B<sinh>, B<cosh>, B<tanh> |
5aabfad6 |
211 | |
212 | The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch |
213 | and cotanh/coth are aliases) |
214 | |
d54bf66f |
215 | B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh> |
5aabfad6 |
216 | |
217 | The arcus (also known as the inverse) functions of the hyperbolic |
218 | sine, cosine, and tangent |
219 | |
d54bf66f |
220 | B<asinh>, B<acosh>, B<atanh> |
5aabfad6 |
221 | |
222 | The arcus cofunctions of the hyperbolic sine, cosine, and tangent |
223 | (acsch/acosech and acoth/acotanh are aliases) |
224 | |
d54bf66f |
225 | B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh> |
5aabfad6 |
226 | |
227 | The trigonometric constant B<pi> is also defined. |
228 | |
d54bf66f |
229 | $pi2 = 2 * B<pi>; |
5aabfad6 |
230 | |
5cd24f17 |
231 | =head2 ERRORS DUE TO DIVISION BY ZERO |
232 | |
233 | The following functions |
234 | |
d54bf66f |
235 | acoth |
5cd24f17 |
236 | acsc |
5cd24f17 |
237 | acsch |
d54bf66f |
238 | asec |
239 | asech |
240 | atanh |
241 | cot |
242 | coth |
243 | csc |
244 | csch |
245 | sec |
246 | sech |
247 | tan |
248 | tanh |
5cd24f17 |
249 | |
250 | cannot be computed for all arguments because that would mean dividing |
8c03c583 |
251 | by zero or taking logarithm of zero. These situations cause fatal |
252 | runtime errors looking like this |
5cd24f17 |
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 |
258 | or |
259 | |
260 | atanh(-1): Logarithm of zero. |
261 | Died at... |
262 | |
263 | For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>, |
264 | C<asech>, C<acsch>, the argument cannot be C<0> (zero). For the |
265 | C<atanh>, C<acoth>, the argument cannot be C<1> (one). For the |
266 | C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one). For the |
267 | C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k * |
268 | pi>, where I<k> is any integer. |
5cd24f17 |
269 | |
270 | =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS |
5aabfad6 |
271 | |
272 | Please note that some of the trigonometric functions can break out |
273 | from the B<real axis> into the B<complex plane>. For example |
274 | C<asin(2)> has no definition for plain real numbers but it has |
275 | definition for complex numbers. |
276 | |
277 | In Perl terms this means that supplying the usual Perl numbers (also |
278 | known as scalars, please see L<perldata>) as input for the |
279 | trigonometric functions might produce as output results that no more |
280 | are simple real numbers: instead they are complex numbers. |
281 | |
282 | The C<Math::Trig> handles this by using the C<Math::Complex> package |
283 | which knows how to handle complex numbers, please see L<Math::Complex> |
284 | for more information. In practice you need not to worry about getting |
285 | complex numbers as results because the C<Math::Complex> takes care of |
286 | details like for example how to display complex numbers. For example: |
287 | |
288 | print asin(2), "\n"; |
3cb6de81 |
289 | |
5aabfad6 |
290 | should produce something like this (take or leave few last decimals): |
291 | |
292 | 1.5707963267949-1.31695789692482i |
293 | |
5cd24f17 |
294 | That is, a complex number with the real part of approximately C<1.571> |
295 | and the imaginary part of approximately C<-1.317>. |
5aabfad6 |
296 | |
d54bf66f |
297 | =head1 PLANE ANGLE CONVERSIONS |
5aabfad6 |
298 | |
299 | (Plane, 2-dimensional) angles may be converted with the following functions. |
300 | |
ace5de91 |
301 | $radians = deg2rad($degrees); |
302 | $radians = grad2rad($gradians); |
3cb6de81 |
303 | |
ace5de91 |
304 | $degrees = rad2deg($radians); |
305 | $degrees = grad2deg($gradians); |
3cb6de81 |
306 | |
ace5de91 |
307 | $gradians = deg2grad($degrees); |
308 | $gradians = rad2grad($radians); |
5aabfad6 |
309 | |
5cd24f17 |
310 | The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians. |
9db5a202 |
311 | The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle. |
312 | If 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 | |
317 | You can also do the wrapping explicitly by rad2rad(), deg2deg(), and |
318 | grad2grad(). |
5aabfad6 |
319 | |
d54bf66f |
320 | =head1 RADIAL COORDINATE CONVERSIONS |
321 | |
322 | B<Radial coordinate systems> are the B<spherical> and the B<cylindrical> |
323 | systems, explained shortly in more detail. |
324 | |
325 | You can import radial coordinate conversion functions by using the |
326 | C<: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 | |
337 | B<All angles are in radians>. |
338 | |
339 | =head2 COORDINATE SYSTEMS |
340 | |
341 | B<Cartesian> coordinates are the usual rectangular I<(x, y, |
342 | z)>-coordinates. |
343 | |
344 | Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional |
345 | coordinates which define a point in three-dimensional space. They are |
346 | based on a sphere surface. The radius of the sphere is B<rho>, also |
347 | known 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> |
349 | coordinate. The angle from the I<z>-axis is B<phi>, also known as the |
350 | I<polar> coordinate. The `North Pole' is therefore I<0, 0, rho>, and |
351 | the `Bay of Guinea' (think of the missing big chunk of Africa) I<0, |
4b0d1da8 |
352 | pi/2, rho>. In geographical terms I<phi> is latitude (northward |
353 | positive, southward negative) and I<theta> is longitude (eastward |
354 | positive, westward negative). |
d54bf66f |
355 | |
4b0d1da8 |
356 | B<BEWARE>: some texts define I<theta> and I<phi> the other way round, |
d54bf66f |
357 | some texts define the I<phi> to start from the horizontal plane, some |
358 | texts use I<r> in place of I<rho>. |
359 | |
360 | Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional |
361 | coordinates which define a point in three-dimensional space. They are |
362 | based on a cylinder surface. The radius of the cylinder is B<rho>, |
363 | also 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> |
365 | coordinate. The third coordinate is the I<z>, pointing up from the |
366 | B<theta>-plane. |
367 | |
368 | =head2 3-D ANGLE CONVERSIONS |
369 | |
370 | Conversions to and from spherical and cylindrical coordinates are |
371 | available. Please notice that the conversions are not necessarily |
372 | reversible because of the equalities like I<pi> angles being equal to |
373 | I<-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 | |
393 | Notice 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 | |
403 | Notice 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 |
408 | |
409 | You can compute spherical distances, called B<great circle distances>, |
7e5f197a |
410 | by 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 |
415 | |
416 | The I<great circle distance> is the shortest distance between two |
417 | points on a sphere. The distance is in C<$rho> units. The C<$rho> is |
418 | optional, it defaults to 1 (the unit sphere), therefore the distance |
419 | defaults to radians. |
420 | |
4b0d1da8 |
421 | If you think geographically the I<theta> are longitudes: zero at the |
422 | Greenwhich meridian, eastward positive, westward negative--and the |
2d06e7d7 |
423 | I<phi> are latitudes: zero at the North Pole, northward positive, |
4b0d1da8 |
424 | southward negative. B<NOTE>: this formula thinks in mathematics, not |
2d06e7d7 |
425 | geographically: the I<phi> zero is at the North Pole, not at the |
426 | Equator on the west coast of Africa (Bay of Guinea). You need to |
427 | subtract your geographical coordinates from I<pi/2> (also known as 90 |
428 | degrees). |
4b0d1da8 |
429 | |
430 | $distance = great_circle_distance($lon0, pi/2 - $lat0, |
431 | $lon1, pi/2 - $lat1, $rho); |
432 | |
7e5f197a |
433 | The direction you must follow the great circle can be computed by the |
434 | great_circle_direction() function: |
435 | |
436 | use Math::Trig 'great_circle_direction'; |
437 | |
438 | $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1); |
439 | |
440 | The result is in radians, zero indicating straight north, pi or -pi |
441 | straight south, pi/2 straight west, and -pi/2 straight east. |
442 | |
443 | Notice that the resulting directions might be somewhat surprising if |
444 | you are looking at a flat worldmap: in such map projections the great |
445 | circles quite often do not look like the shortest routes-- but for |
446 | example the shortest possible routes from Europe or North America to |
447 | Asia do often cross the polar regions. |
448 | |
51301382 |
449 | =head1 EXAMPLES |
d54bf66f |
450 | |
7e5f197a |
451 | To calculate the distance between London (51.3N 0.5W) and Tokyo |
452 | (35.7N 139.8E) in kilometers: |
d54bf66f |
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 |
462 | The 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 | |
470 | The answers may be off by few percentages because of the irregular |
471 | (slightly aspherical) form of the Earth. The formula used for |
472 | grear circle distances |
41bd693c |
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 | |
479 | is also somewhat unreliable for small distances (for locations |
480 | separated less than about five degrees) because it uses arc cosine |
481 | which is rather ill-conditioned for values close to zero. |
d54bf66f |
482 | |
5cd24f17 |
483 | =head1 BUGS |
5aabfad6 |
484 | |
5cd24f17 |
485 | Saying C<use Math::Trig;> exports many mathematical routines in the |
486 | caller environment and even overrides some (C<sin>, C<cos>). This is |
487 | construed as a feature by the Authors, actually... ;-) |
5aabfad6 |
488 | |
5cd24f17 |
489 | The code is not optimized for speed, especially because we use |
490 | C<Math::Complex> and thus go quite near complex numbers while doing |
491 | the computations even when the arguments are not. This, however, |
492 | cannot be completely avoided if we want things like C<asin(2)> to give |
493 | an answer instead of giving a fatal runtime error. |
5aabfad6 |
494 | |
5cd24f17 |
495 | =head1 AUTHORS |
5aabfad6 |
496 | |
ace5de91 |
497 | Jarkko Hietaniemi <F<jhi@iki.fi>> and |
6e238990 |
498 | Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>. |
5aabfad6 |
499 | |
500 | =cut |
501 | |
502 | # eof |