lib/perl5/5.6.1/Math/Trig.pm

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