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