lib/Math/Trig.pm

#
# Trigonometric functions, mostly inherited from Math::Complex.
# -- Jarkko Hietaniemi, since April 1997
# -- Raphael Manfredi, September 1996 (indirectly: because of Math::Complex)
#

require Exporter;
package Math::Trig;

use 5.006;
use strict;

use Math::Complex qw(:trig);

our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);

@ISA = qw(Exporter);

$VERSION = 1.02;

my @angcnv = qw(rad2deg rad2grad
                deg2rad deg2grad
                grad2rad grad2deg);

@EXPORT = (@{$Math::Complex::EXPORT_TAGS{'trig'}},
           @angcnv);

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', 'great_circle_direction');

%EXPORT_TAGS = ('radial' => [ @rdlcnv ]);

sub pi2  () { 2 * pi }
sub pip2 () { pi / 2 }

sub DR  () { pi2/360 }
sub RD  () { 360/pi2 }
sub DG  () { 400/360 }
sub GD  () { 360/400 }
sub RG  () { 400/pi2 }
sub GR  () { pi2/400 }

#
# Truncating remainder.
#

sub remt ($$) {
    # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available.
    $_[0] - $_[1] * int($_[0] / $_[1]);
}

#
# Angle conversions.
#

sub rad2rad($)     { remt($_[0], pi2) }

sub deg2deg($)     { remt($_[0], 360) }

sub grad2grad($)   { remt($_[0], 400) }

sub rad2deg ($;$)  { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) }

sub deg2rad ($;$)  { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) }

sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) }

sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) }

sub rad2grad ($;$) { my $d = RG * $_[0]; $_[1] ? $d : grad2grad($d) }

sub grad2rad ($;$) { my $d = GR * $_[0]; $_[1] ? $d : rad2rad($d) }

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 ) );
}

sub great_circle_direction {
    my ( $theta0, $phi0, $theta1, $phi1 ) = @_;

    my $distance = &great_circle_distance;

    my $lat0 = pip2 - $phi0;
    my $lat1 = pip2 - $phi1;

    my $direction =
        acos((sin($lat1) - sin($lat0) * cos($distance)) /
             (cos($lat0) * sin($distance)));

    $direction = pi2 - $direction
        if sin($theta1 - $theta0) < 0;

    return rad2rad($direction);
}

1;

__END__
=pod

=head1 NAME

Math::Trig - trigonometric functions

=head1 SYNOPSIS

        use Math::Trig;

        $x = tan(0.9);
        $y = acos(3.7);
        $z = asin(2.4);

        $halfpi = pi/2;

        $rad = deg2rad(120);

=head1 DESCRIPTION

C<Math::Trig> defines many trigonometric functions not defined by the
core Perl which defines only the C<sin()> and C<cos()>.  The constant
B<pi> is also defined as are a few convenience functions for angle
conversions.

=head1 TRIGONOMETRIC FUNCTIONS

The tangent

=over 4

=item B<tan>

=back

The cofunctions of the sine, cosine, and tangent (cosec/csc and cotan/cot
are aliases)

B<csc>, B<cosec>, B<sec>, B<sec>, B<cot>, B<cotan>

The arcus (also known as the inverse) functions of the sine, cosine,
and tangent

B<asin>, B<acos>, B<atan>

The principal value of the arc tangent of y/x

B<atan2>(y, x)

The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc
and acotan/acot are aliases)

B<acsc>, B<acosec>, B<asec>, B<acot>, B<acotan>

The hyperbolic sine, cosine, and tangent

B<sinh>, B<cosh>, B<tanh>

The cofunctions of the hyperbolic sine, cosine, and tangent (cosech/csch
and cotanh/coth are aliases)

B<csch>, B<cosech>, B<sech>, B<coth>, B<cotanh>

The arcus (also known as the inverse) functions of the hyperbolic
sine, cosine, and tangent

B<asinh>, B<acosh>, B<atanh>

The arcus cofunctions of the hyperbolic sine, cosine, and tangent
(acsch/acosech and acoth/acotanh are aliases)

B<acsch>, B<acosech>, B<asech>, B<acoth>, B<acotanh>

The trigonometric constant B<pi> is also defined.

$pi2 = 2 * B<pi>;

=head2 ERRORS DUE TO DIVISION BY ZERO

The following functions

        acoth
        acsc
        acsch
        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
runtime errors looking like this

        cot(0): Division by zero.
        (Because in the definition of cot(0), the divisor sin(0) is 0)
        Died at ...

or

        atanh(-1): Logarithm of zero.
        Died at...

For the C<csc>, C<cot>, C<asec>, C<acsc>, C<acot>, C<csch>, C<coth>,
C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
pi>, where I<k> is any integer.

=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS

Please note that some of the trigonometric functions can break out
from the B<real axis> into the B<complex plane>. For example
C<asin(2)> has no definition for plain real numbers but it has
definition for complex numbers.

In Perl terms this means that supplying the usual Perl numbers (also
known as scalars, please see L<perldata>) as input for the
trigonometric functions might produce as output results that no more
are simple real numbers: instead they are complex numbers.

The C<Math::Trig> handles this by using the C<Math::Complex> package
which knows how to handle complex numbers, please see L<Math::Complex>
for more information. In practice you need not to worry about getting
complex numbers as results because the C<Math::Complex> takes care of
details like for example how to display complex numbers. For example:

        print asin(2), "\n";

should produce something like this (take or leave few last decimals):

        1.5707963267949-1.31695789692482i

That is, a complex number with the real part of approximately C<1.571>
and the imaginary part of approximately C<-1.317>.

=head1 PLANE ANGLE CONVERSIONS

(Plane, 2-dimensional) angles may be converted with the following functions.

        $radians  = deg2rad($degrees);
        $radians  = grad2rad($gradians);

        $degrees  = rad2deg($radians);
        $degrees  = grad2deg($gradians);

        $gradians = deg2grad($degrees);
        $gradians = rad2grad($radians);

The full circle is 2 I<pi> radians or I<360> degrees or I<400> gradians.
The result is by default wrapped to be inside the [0, {2pi,360,400}[ circle.
If you don't want this, supply a true second argument:

        $zillions_of_radians  = deg2rad($zillions_of_degrees, 1);
        $negative_degrees     = rad2deg($negative_radians, 1);

You can also do the wrapping explicitly by rad2rad(), deg2deg(), and
grad2grad().

=head1 RADIAL COORDINATE CONVERSIONS

B<Radial coordinate systems> are the B<spherical> and the B<cylindrical>
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<All angles are in radians>.

=head2 COORDINATE SYSTEMS

B<Cartesian> 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<rho>, also
known as the I<radial> coordinate.  The angle in the I<xy>-plane
(around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
coordinate.  The angle from the I<z>-axis is B<phi>, also known as the
I<polar> 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>.  In geographical terms I<phi> is latitude (northward
positive, southward negative) and I<theta> is longitude (eastward
positive, westward negative).

B<BEWARE>: some texts define I<theta> and I<phi> the other way round,
some texts define the I<phi> to start from the horizontal plane, some
texts use I<r> in place of I<rho>.

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<rho>,
also known as the I<radial> coordinate.  The angle in the I<xy>-plane
(around the I<z>-axis) is B<theta>, also known as the I<azimuthal>
coordinate.  The third coordinate is the I<z>, pointing up from the
B<theta>-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<pi> 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 AND DIRECTIONS

You can compute spherical distances, called B<great circle distances>,
by importing the great_circle_distance() function:

  use Math::Trig 'great_circle_distance';

  $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);

The I<great circle distance> 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.

If you think geographically the I<theta> are longitudes: zero at the
Greenwhich meridian, eastward positive, westward negative--and the
I<phi> are latitudes: zero at the North Pole, northward positive,
southward negative.  B<NOTE>: this formula thinks in mathematics, not
geographically: the I<phi> zero is at the North Pole, not at the
Equator on the west coast of Africa (Bay of Guinea).  You need to
subtract your geographical coordinates from I<pi/2> (also known as 90
degrees).

  $distance = great_circle_distance($lon0, pi/2 - $lat0,
                                    $lon1, pi/2 - $lat1, $rho);

The direction you must follow the great circle can be computed by the
great_circle_direction() function:

  use Math::Trig 'great_circle_direction';

  $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);

The result is in radians, zero indicating straight north, pi or -pi
straight south, pi/2 straight west, and -pi/2 straight east.

Notice that the resulting directions might be somewhat surprising if
you are looking at a flat worldmap: in such map projections the great
circles quite often do not look like the shortest routes-- but for
example the shortest possible routes from Europe or North America to
Asia do often cross the polar regions.

=head1 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 direction you would have to go from London to Tokyo

        use Math::Trig qw(great_circle_direction);

        $rad = great_circle_direction(@L, @T);

=head2 CAVEAT FOR GREAT CIRCLE FORMULAS

The answers may be off by few percentages because of the irregular
(slightly aspherical) form of the Earth.  The formula used for
grear circle distances

        lat0 = 90 degrees - phi0
        lat1 = 90 degrees - phi1
        d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
                       sin(lat0) * sin(lat1))

is also somewhat unreliable for small distances (for locations
separated less than about five degrees) because it uses arc cosine
which is rather ill-conditioned for values close to zero.

=head1 BUGS

Saying C<use Math::Trig;> exports many mathematical routines in the
caller environment and even overrides some (C<sin>, C<cos>).  This is
construed as a feature by the Authors, actually... ;-)

The code is not optimized for speed, especially because we use
C<Math::Complex> and thus go quite near complex numbers while doing
the computations even when the arguments are not. This, however,
cannot be completely avoided if we want things like C<asin(2)> to give
an answer instead of giving a fatal runtime error.

=head1 AUTHORS

Jarkko Hietaniemi <F<jhi@iki.fi>> and 
Raphael Manfredi <F<Raphael_Manfredi@pobox.com>>.

=cut

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