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