From: Jarkko Hietaniemi Date: Mon, 29 Jun 1998 16:28:53 +0000 (+0300) Subject: add patch to integrate Math::Trig::Radial into Math::Trig X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=d54bf66f316184872ac87c4594ab8c34b6dedc9e;p=p5sagit%2Fp5-mst-13.2.git add patch to integrate Math::Trig::Radial into Math::Trig Message-Id: <199806291328.QAA16916@alpha.hut.fi> Subject: [PATCH] 5.004_68 (or 5.004_04): radial trig p4raw-id: //depot/perl@1267 --- diff --git a/MANIFEST b/MANIFEST index 5c1b5ba..6dd7c59 100644 --- a/MANIFEST +++ b/MANIFEST @@ -481,7 +481,6 @@ lib/Math/BigFloat.pm An arbitrary precision floating-point arithmetic package lib/Math/BigInt.pm An arbitrary precision integer arithmetic package lib/Math/Complex.pm A Complex package lib/Math/Trig.pm A simple interface to complex trigonometry -lib/Math/Trig/Radial.pm Spherical and cylindrical trigonometry lib/Net/Ping.pm Hello, anybody home? lib/Net/hostent.pm By-name interface to Perl's builtin gethost* lib/Net/netent.pm By-name interface to Perl's builtin getnet* diff --git a/lib/Math/Trig.pm b/lib/Math/Trig.pm index a1cbb07..7192d76 100644 --- a/lib/Math/Trig.pm +++ b/lib/Math/Trig.pm @@ -1,6 +1,6 @@ # # Trigonometric functions, mostly inherited from Math::Complex. -# -- Jarkko Hietaniemi, April 1997 +# -- Jarkko Hietaniemi, since April 1997 # -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex) # @@ -13,7 +13,7 @@ use Math::Complex qw(:trig); use vars qw($VERSION $PACKAGE @ISA - @EXPORT); + @EXPORT @EXPORT_OK %EXPORT_TAGS); @ISA = qw(Exporter); @@ -26,13 +26,25 @@ my @angcnv = qw(rad2deg rad2grad @EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}}, @angcnv); -use constant pi2 => 2 * pi; -use constant DR => pi2/360; -use constant RD => 360/pi2; -use constant DG => 400/360; -use constant GD => 360/400; -use constant RG => 400/pi2; -use constant GR => pi2/400; +my @rdlcnv = qw(cartesian_to_cylindrical + cartesian_to_spherical + cylindrical_to_cartesian + cylindrical_to_spherical + spherical_to_cartesian + spherical_to_cylindrical); + +@EXPORT_OK = (@rdlcnv, 'great_circle_distance'); + +%EXPORT_TAGS = ('radial' => [ @rdlcnv ]); + +use constant pi2 => 2 * pi; +use constant pip2 => pi / 2; +use constant DR => pi2/360; +use constant RD => 360/pi2; +use constant DG => 400/360; +use constant GD => 360/400; +use constant RG => 400/pi2; +use constant GR => pi2/400; # # Truncating remainder. @@ -59,6 +71,61 @@ sub rad2grad ($) { remt(RG * $_[0], 400) } sub grad2rad ($) { remt(GR * $_[0], pi2) } +sub cartesian_to_spherical { + my ( $x, $y, $z ) = @_; + + my $rho = sqrt( $x * $x + $y * $y + $z * $z ); + + return ( $rho, + atan2( $y, $x ), + $rho ? acos( $z / $rho ) : 0 ); +} + +sub spherical_to_cartesian { + my ( $rho, $theta, $phi ) = @_; + + return ( $rho * cos( $theta ) * sin( $phi ), + $rho * sin( $theta ) * sin( $phi ), + $rho * cos( $phi ) ); +} + +sub spherical_to_cylindrical { + my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); + + return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); +} + +sub cartesian_to_cylindrical { + my ( $x, $y, $z ) = @_; + + return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); +} + +sub cylindrical_to_cartesian { + my ( $rho, $theta, $z ) = @_; + + return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); +} + +sub cylindrical_to_spherical { + return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); +} + +sub great_circle_distance { + my ( $theta0, $phi0, $theta1, $phi1, $rho ) = @_; + + $rho = 1 unless defined $rho; # Default to the unit sphere. + + my $lat0 = pip2 - $phi0; + my $lat1 = pip2 - $phi1; + + return $rho * + acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + + sin( $lat0 ) * sin( $lat1 ) ); +} + +=pod + =head1 NAME Math::Trig - trigonometric functions @@ -86,68 +153,72 @@ conversions. The tangent - tan +=over 4 + +=item B + +=back The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot are aliases) - csc cosec sec cot cotan +B, B, B, B, B, B The arcus (also known as the inverse) functions of the sine, cosine, and tangent - asin acos atan +B, B, B The principal value of the arc tangent of y/x - atan2(y, x) +B(y, x) The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc and acotan/acot are aliases) - acsc acosec asec acot acotan +B, B, B, B, B The hyperbolic sine, cosine, and tangent - sinh cosh tanh +B, B, B The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch and cotanh/coth are aliases) - csch cosech sech coth cotanh +B, B, B, B, B The arcus (also known as the inverse) functions of the hyperbolic sine, cosine, and tangent - asinh acosh atanh +B, B, B The arcus cofunctions of the hyperbolic sine, cosine, and tangent (acsch/acosech and acoth/acotanh are aliases) - acsch acosech asech acoth acotanh +B, B, B, B, B The trigonometric constant B is also defined. - $pi2 = 2 * pi; +$pi2 = 2 * B; =head2 ERRORS DUE TO DIVISION BY ZERO The following functions - tan - sec - csc - cot - asec + acoth acsc - tanh - sech - csch - coth - atanh - asech acsch - acoth + asec + asech + atanh + cot + coth + csc + csch + sec + sech + tan + tanh cannot be computed for all arguments because that would mean dividing by zero or taking logarithm of zero. These situations cause fatal @@ -196,7 +267,7 @@ should produce something like this (take or leave few last decimals): That is, a complex number with the real part of approximately C<1.571> and the imaginary part of approximately C<-1.317>. -=head1 ANGLE CONVERSIONS +=head1 PLANE ANGLE CONVERSIONS (Plane, 2-dimensional) angles may be converted with the following functions. @@ -211,6 +282,121 @@ and the imaginary part of approximately C<-1.317>. The full circle is 2 I radians or I<360> degrees or I<400> gradians. +=head1 RADIAL COORDINATE CONVERSIONS + +B are the B and the B +systems, explained shortly in more detail. + +You can import radial coordinate conversion functions by using the +C<:radial> tag: + + use Math::Trig ':radial'; + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +B. + +=head2 COORDINATE SYSTEMS + +B coordinates are the usual rectangular I<(x, y, +z)>-coordinates. + +Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a sphere surface. The radius of the sphere is B, also +known as the I coordinate. The angle in the I-plane +(around the I-axis) is B, also known as the I +coordinate. The angle from the I-axis is B, also known as the +I coordinate. The `North Pole' is therefore I<0, 0, rho>, and +the `Bay of Guinea' (think of the missing big chunk of Africa) I<0, +pi/2, rho>. + +B: some texts define I and I the other way round, +some texts define the I to start from the horizontal plane, some +texts use I in place of I. + +Cylindrical coordinates, I<(rho, theta, z)>, are three-dimensional +coordinates which define a point in three-dimensional space. They are +based on a cylinder surface. The radius of the cylinder is B, +also known as the I coordinate. The angle in the I-plane +(around the I-axis) is B, also known as the I +coordinate. The third coordinate is the I, pointing up from the +B-plane. + +=head2 3-D ANGLE CONVERSIONS + +Conversions to and from spherical and cylindrical coordinates are +available. Please notice that the conversions are not necessarily +reversible because of the equalities like I angles being equal to +I<-pi> angles. + +=over 4 + +=item cartesian_to_cylindrical + + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + +=item cartesian_to_spherical + + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + +=item cylindrical_to_cartesian + + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + +=item cylindrical_to_spherical + + ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + +Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>. + +=item spherical_to_cartesian + + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); + +=item spherical_to_cylindrical + + ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + +Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. + +=back + +=head1 GREAT CIRCLE DISTANCES + +You can compute spherical distances, called B, +by importing the C function: + + use Math::Trig 'great_circle_distance' + + $distance = great_circle_distance($theta0, $phi0, $theta1, $phi, [, $rho]); + +The I is the shortest distance between two +points on a sphere. The distance is in C<$rho> units. The C<$rho> is +optional, it defaults to 1 (the unit sphere), therefore the distance +defaults to radians. + +=head EXAMPLES + +To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N +139.8E) in kilometers: + + use Math::Trig qw(great_circle_distance deg2rad); + + # Notice the 90 - latitude: phi zero is at the North Pole. + @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); + @T = (deg2rad(139.8),deg2rad(90 - 35.7)); + + $km = great_circle_distance(@L, @T, 6378); + +The answer may be off by up to 0.3% because of the irregular (slightly +aspherical) form of the Earth. + =head1 BUGS Saying C exports many mathematical routines in the diff --git a/lib/Math/Trig/Radial.pm b/lib/Math/Trig/Radial.pm deleted file mode 100644 index 0001cb7..0000000 --- a/lib/Math/Trig/Radial.pm +++ /dev/null @@ -1,193 +0,0 @@ -package Math::Trig::Radial; - -use strict; -use vars qw(@ISA @EXPORT); -@ISA = qw(Exporter); - -@EXPORT = - qw( - cartesian_to_cylindrical - cartesian_to_spherical - cylindrical_to_cartesian - cylindrical_to_spherical - spherical_to_cartesian - spherical_to_cylindrical - great_circle_distance - ); - -use Math::Trig; - -sub pip2 { pi/2 } - -=pod - -=head1 NAME - -Math::Trig::Radial - spherical and cylindrical trigonometry - -=head1 SYNOPSIS - - use Math::Trig::Radial; - - ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); - ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); - ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); - ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); - ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); - ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); - -=head1 DESCRIPTION - -This module contains a few basic spherical and cylindrical -trigonometric formulas. B, if needed -use C angle unit conversions. - -=head2 COORDINATE SYSTEMS - -B coordinates are the usual rectangular I-coordinates. - -Spherical coordinates are three-dimensional coordinates which define a -point in three-dimensional space. They are based on a sphere surface. -The radius of the sphere is B, also known as the I -coordinate. The angle in the I-plane (around the I-axis) is -B, also known as the I coordinate. The angle from -the I-axis is B, also known as the I coordinate. The -`North Pole' is therefore I<0, 0, rho>, and the `Bay of Guinea' (think -Africa) I<0, pi/2, rho>. - -Cylindrical coordinates are three-dimensional coordinates which define -a point in three-dimensional space. They are based on a cylinder -surface. The radius of the cylinder is B, also known as the -I coordinate. The angle in the I-plane (around the -I-axis) is B, also known as the I coordinate. -The third coordinate is the I. - -=head2 CONVERSIONS - -Conversions to and from spherical and cylindrical coordinates are -available. Please notice that the conversions are not necessarily -reversible because of the equalities like I angles equals I<-pi> -angles. - -=over 4 - -=item cartesian_to_cylindrical - - ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); - -=item cartesian_to_spherical - - ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); - -=item cylindrical_to_cartesian - - ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); - -=item cylindrical_to_spherical - - ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); - -Notice that when C<$z> is not 0 C<$rho_s> is not equal to C<$rho_c>. - -=item spherical_to_cartesian - - ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); - -=item spherical_to_cylindrical - - ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); - -Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. - -=back - -=head2 GREAT CIRCLE DISTANCE - - $distance = great_circle_distance($x0, $y0, $z0, $x1, $y1, $z1 [, $rho]); - -The I is the shortest distance between two -points on a sphere. The distance is in C<$rho> units. The C<$rho> is -optional, it defaults to 1 (the unit sphere), therefore the distance -defaults to radians. The coordinates C<$x0> ... C<$z1> are in -cartesian coordinates. - -=head EXAMPLES - -To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N -139.8E) in kilometers: - - use Math::Trig::Radial; - use Math::Trig; - - my @L = spherical_to_cartesian(1, map { deg2rad $_ } qw(51.3 -0.5)); - my @T = spherical_to_cartesian(1, map { deg2rad $_ } qw(35.7 139.8)); - - $km = great_circle_distance(@L, @T, 6378); - -The answer may be off by up to 0.3% because of the irregular (slightly -aspherical) form of the Earth. - -=head2 AUTHOR - -Jarkko Hietaniemi Fjhi@iki.fiE> - -=cut - -sub cartesian_to_spherical { - my ( $x, $y, $z ) = @_; - - my $rho = sqrt( $x * $x + $y * $y + $z * $z ); - - return ( $rho, - atan2( $y, $x ), - $rho ? acos( $z / $rho ) : 0 ); -} - -sub spherical_to_cartesian { - my ( $rho, $theta, $phi ) = @_; - - return ( $rho * cos( $theta ) * sin( $phi ), - $rho * sin( $theta ) * sin( $phi ), - $rho * cos( $phi ) ); -} - -sub spherical_to_cylindrical { - my ( $x, $y, $z ) = spherical_to_cartesian( @_ ); - - return ( sqrt( $x * $x + $y * $y ), $_[1], $z ); -} - -sub cartesian_to_cylindrical { - my ( $x, $y, $z ) = @_; - - return ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z ); -} - -sub cylindrical_to_cartesian { - my ( $rho, $theta, $z ) = @_; - - return ( $rho * cos( $theta ), $rho * sin( $theta ), $z ); -} - -sub cylindrical_to_spherical { - return ( cartesian_to_spherical( cylindrical_to_cartesian( @_ ) ) ); -} - -sub great_circle_distance { - my ( $x0, $y0, $z0, $x1, $y1, $z1, $rho ) = @_; - - $rho = 1 unless defined $rho; # Default to the unit sphere. - - my ( $r0, $theta0, $phi0 ) = cartesian_to_spherical( $x0, $y0, $z0 ); - my ( $r1, $theta1, $phi1 ) = cartesian_to_spherical( $x1, $y1, $z1 ); - - my $lat0 = pip2 - $phi0; - my $lat1 = pip2 - $phi1; - - return $rho * - acos(cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) + - sin( $lat0 ) * sin( $lat1 ) ); -} - -1; - diff --git a/t/lib/trig.t b/t/lib/trig.t index c2bc2a8..0914711 100755 --- a/t/lib/trig.t +++ b/t/lib/trig.t @@ -25,7 +25,7 @@ sub near ($$;$) { abs($_[0] - $_[1]) < (defined $_[2] ? $_[2] : $eps); } -print "1..7\n"; +print "1..20\n"; $x = 0.9; print 'not ' unless (near(tan($x), sin($x) / cos($x))); @@ -54,4 +54,103 @@ print "ok 6\n"; print 'not ' unless (near(rad2deg(pi), 180)); print "ok 7\n"; +use Math::Trig ':radial'; + +{ + my ($r,$t,$z) = cartesian_to_cylindrical(1,1,1); + + print 'not ' unless (near($r, sqrt(2))) and + (near($t, deg2rad(45))) and + (near($z, 1)); + print "ok 8\n"; + + ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z); + + print 'not ' unless (near($x, 1)) and + (near($y, 1)) and + (near($z, 1)); + print "ok 9\n"; + + ($r,$t,$z) = cartesian_to_cylindrical(1,1,0); + + print 'not ' unless (near($r, sqrt(2))) and + (near($t, deg2rad(45))) and + (near($z, 0)); + print "ok 10\n"; + + ($x,$y,$z) = cylindrical_to_cartesian($r, $t, $z); + + print 'not ' unless (near($x, 1)) and + (near($y, 1)) and + (near($z, 0)); + print "ok 11\n"; +} + +{ + my ($r,$t,$f) = cartesian_to_spherical(1,1,1); + + print 'not ' unless (near($r, sqrt(3))) and + (near($t, deg2rad(45))) and + (near($f, atan2(sqrt(2), 1))); + print "ok 12\n"; + + ($x,$y,$z) = spherical_to_cartesian($r, $t, $f); + + print 'not ' unless (near($x, 1)) and + (near($y, 1)) and + (near($z, 1)); + print "ok 13\n"; + + ($r,$t,$f) = cartesian_to_spherical(1,1,0); + + print 'not ' unless (near($r, sqrt(2))) and + (near($t, deg2rad(45))) and + (near($f, deg2rad(90))); + print "ok 14\n"; + + ($x,$y,$z) = spherical_to_cartesian($r, $t, $f); + + print 'not ' unless (near($x, 1)) and + (near($y, 1)) and + (near($z, 0)); + print "ok 15\n"; +} + +{ + my ($r,$t,$z) = cylindrical_to_spherical(spherical_to_cylindrical(1,1,1)); + + print 'not ' unless (near($r, 1)) and + (near($t, 1)) and + (near($z, 1)); + print "ok 16\n"; + + ($r,$t,$z) = spherical_to_cylindrical(cylindrical_to_spherical(1,1,1)); + + print 'not ' unless (near($r, 1)) and + (near($t, 1)) and + (near($z, 1)); + print "ok 17\n"; +} + +{ + use Math::Trig 'great_circle_distance'; + + print 'not ' + unless (near(great_circle_distance(0, 0, 0, pi/2), pi/2)); + print "ok 18\n"; + + print 'not ' + unless (near(great_circle_distance(0, 0, pi, pi), pi)); + print "ok 19\n"; + + # London to Tokyo. + my @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); + my @T = (deg2rad(139.8),deg2rad(90 - 35.7)); + + my $km = great_circle_distance(@L, @T, 6378); + + print 'not ' unless (near($km, 9605.26637021388)); + print "ok 20\n"; +} + # eof