| 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 |