package Math::Complex;
-our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);
+use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
-$VERSION = 1.34;
+$VERSION = 1.35;
BEGIN {
unless ($^O eq 'unicosmk') {
my %LOGN;
# Regular expression for floating point numbers.
-my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?)))';
+# These days we could use Scalar::Util::lln(), I guess.
+my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
require Exporter;
sqrt log ln
log10 logn cbrt root
cplx cplxe
+ atan2
),
@trig);
+@EXPORT_OK = qw(decplx);
+
%EXPORT_TAGS = (
'trig' => [@trig],
);
# Die on bad *make() arguments.
sub _cannot_make {
- die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
+ die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
}
-sub _remake {
+sub _make {
my $arg = shift;
- my ($made, $p, $q);
+ my ($p, $q);
- if ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
+ if ($arg =~ /^$gre$/) {
+ ($p, $q) = ($1, 0);
+ } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
($p, $q) = ($1 || 0, $2);
- $made = 'cart';
- } elsif ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
+ } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
($p, $q) = ($1, $2 || 0);
- $made = 'exp';
}
- if ($made) {
+ if (defined $p) {
$p =~ s/^\+//;
+ $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
$q =~ s/^\+//;
+ $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
}
- return ($made, $p, $q);
+ return ($p, $q);
+}
+
+sub _emake {
+ my $arg = shift;
+ my ($p, $q);
+
+ if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
+ ($p, $q) = ($1, $2 || 0);
+ } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
+ ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
+ } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
+ ($p, $q) = ($1, 0);
+ } elsif ($arg =~ /^\s*$gre\s*$/) {
+ ($p, $q) = ($1, 0);
+ }
+
+ if (defined $p) {
+ $p =~ s/^\+//;
+ $q =~ s/^\+//;
+ $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
+ $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
+ }
+
+ return ($p, $q);
}
#
# Create a new complex number (cartesian form)
#
sub make {
- my $self = bless {}, shift;
- my ($re, $im) = @_;
- if (@_ == 1) {
- my ($remade, $p, $q) = _remake($re);
- if ($remade) {
- if ($remade eq 'cart') {
- ($re, $im) = ($p, $q);
- } else {
- return (ref $self)->emake($p, $q);
- }
- }
- }
- my $rre = ref $re;
- if ( $rre ) {
- if ( $rre eq ref $self ) {
- $re = Re($re);
- } else {
- _cannot_make("real part", $rre);
- }
- }
- my $rim = ref $im;
- if ( $rim ) {
- if ( $rim eq ref $self ) {
- $im = Im($im);
- } else {
- _cannot_make("imaginary part", $rim);
- }
- }
+ my $self = bless {}, shift;
+ my ($re, $im);
+ if (@_ == 0) {
+ ($re, $im) = (0, 0);
+ } elsif (@_ == 1) {
+ return (ref $self)->emake($_[0])
+ if ($_[0] =~ /^\s*\[/);
+ ($re, $im) = _make($_[0]);
+ } elsif (@_ == 2) {
+ ($re, $im) = @_;
+ }
+ if (defined $re) {
_cannot_make("real part", $re) unless $re =~ /^$gre$/;
- $im ||= 0;
- _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
- $self->{'cartesian'} = [ $re, $im ];
- $self->{c_dirty} = 0;
- $self->{p_dirty} = 1;
- $self->display_format('cartesian');
- return $self;
+ }
+ $im ||= 0;
+ _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
+ $self->set_cartesian([$re, $im ]);
+ $self->display_format('cartesian');
+
+ return $self;
}
#
# Create a new complex number (exponential form)
#
sub emake {
- my $self = bless {}, shift;
- my ($rho, $theta) = @_;
- if (@_ == 1) {
- my ($remade, $p, $q) = _remake($rho);
- if ($remade) {
- if ($remade eq 'exp') {
- ($rho, $theta) = ($p, $q);
- } else {
- return (ref $self)->make($p, $q);
- }
- }
- }
- my $rrh = ref $rho;
- if ( $rrh ) {
- if ( $rrh eq ref $self ) {
- $rho = rho($rho);
- } else {
- _cannot_make("rho", $rrh);
- }
- }
- my $rth = ref $theta;
- if ( $rth ) {
- if ( $rth eq ref $self ) {
- $theta = theta($theta);
- } else {
- _cannot_make("theta", $rth);
- }
- }
+ my $self = bless {}, shift;
+ my ($rho, $theta);
+ if (@_ == 0) {
+ ($rho, $theta) = (0, 0);
+ } elsif (@_ == 1) {
+ return (ref $self)->make($_[0])
+ if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
+ ($rho, $theta) = _emake($_[0]);
+ } elsif (@_ == 2) {
+ ($rho, $theta) = @_;
+ }
+ if (defined $rho && defined $theta) {
if ($rho < 0) {
$rho = -$rho;
$theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
}
+ }
+ if (defined $rho) {
_cannot_make("rho", $rho) unless $rho =~ /^$gre$/;
- $theta ||= 0;
- _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
- $self->{'polar'} = [$rho, $theta];
- $self->{p_dirty} = 0;
- $self->{c_dirty} = 1;
- $self->display_format('polar');
- return $self;
+ }
+ $theta ||= 0;
+ _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
+ $self->set_polar([$rho, $theta]);
+ $self->display_format('polar');
+
+ return $self;
}
sub new { &make } # For backward compatibility only.
sub polar {$_[0]->{p_dirty} ?
$_[0]->update_polar : $_[0]->{'polar'}}
-sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
-sub set_polar { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
+sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
+ $_[0]->{'cartesian'} = $_[1] }
+sub set_polar { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
+ $_[0]->{'polar'} = $_[1] }
#
# ->update_cartesian
# Die on bad root.
#
sub _rootbad {
- my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
+ my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
my @up = caller(1);
# z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
#
sub root {
- my ($z, $n) = @_;
+ my ($z, $n, $k) = @_;
_rootbad($n) if ($n < 1 or int($n) != $n);
my ($r, $t) = ref $z ?
@{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
- my @root;
- my $k;
my $theta_inc = pit2 / $n;
my $rho = $r ** (1/$n);
- my $theta;
my $cartesian = ref $z && $z->{c_dirty} == 0;
- for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
- my $w = cplxe($rho, $theta);
- # Yes, $cartesian is loop invariant.
- push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
+ if (@_ == 2) {
+ my @root;
+ for (my $i = 0, my $theta = $t / $n;
+ $i < $n;
+ $i++, $theta += $theta_inc) {
+ my $w = cplxe($rho, $theta);
+ # Yes, $cartesian is loop invariant.
+ push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
+ }
+ return @root;
+ } elsif (@_ == 3) {
+ my $w = cplxe($rho, $t / $n + $k * $theta_inc);
+ return $cartesian ? cplx(@{$w->cartesian}) : $w;
}
- return @root;
}
#
#
# (atan2)
#
-# Compute atan(z1/z2).
+# Compute atan(z1/z2), minding the right quadrant.
#
sub atan2 {
my ($z1, $z2, $inverted) = @_;
my ($re1, $im1, $re2, $im2);
if ($inverted) {
($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
- ($re2, $im2) = @{$z1->cartesian};
+ ($re2, $im2) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
} else {
- ($re1, $im1) = @{$z1->cartesian};
+ ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
}
- if ($im2 == 0) {
- return CORE::atan2($re1, $re2) if $im1 == 0;
- return ($im1<=>0) * pip2 if $re2 == 0;
+ if ($im1 || $im2) {
+ # In MATLAB the imaginary parts are ignored.
+ # warn "atan2: Imaginary parts ignored";
+ # http://documents.wolfram.com/mathematica/functions/ArcTan
+ # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
+ my $s = $z1 * $z1 + $z2 * $z2;
+ _divbyzero("atan2") if $s == 0;
+ my $i = &i;
+ my $r = $z2 + $z1 * $i;
+ return -$i * &log($r / &sqrt( $s ));
}
- my $w = atan($z1/$z2);
- my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
- $u += pi if $re2 < 0;
- $u -= pit2 if $u > pi;
- return cplx($u, $v);
+ return CORE::atan2($re1, $re2);
}
#
log(z) = log(r1) + i*t
sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
- atan2(z1, z2) = atan(z1/z2)
+ atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
+
+The definition used for complex arguments of atan2() is
+
+ -i log((x + iy)/sqrt(x*x+y*y))
The following extra operations are supported on both real and complex
numbers:
(root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
+You can return the I<k>th root directly by C<root(z, n, k)>,
+indexing starting from I<zero> and ending at I<n - 1>.
+
The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
order to ensure its restriction to real numbers is conform to what you
would expect, the comparison is run on the real part of the complex
2-3i
-3i
[2,3]
+ [2,-3pi/4]
[2]
in which case the appropriate cartesian and exponential components
will be parsed from the string and used to create new complex numbers.
The imaginary component and the theta, respectively, will default to zero.
-=head1 STRINGIFICATION
+The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
+understand the case of no arguments: this means plain zero or (0, 0).
+
+=head1 DISPLAYING
When printed, a complex number is usually shown under its cartesian
style I<a+bi>, but there are legitimate cases where the polar style
-I<[r,t]> is more appropriate.
+I<[r,t]> is more appropriate. The process of converting the complex
+number into a string that can be displayed is known as I<stringification>.
By calling the class method C<Math::Complex::display_format> and
supplying either C<"polar"> or C<"cartesian"> as an argument, you
(where I<n> is a positive integer and I<k> an integer within [-9, +9]),
this is called I<polar pretty-printing>.
+For the reverse of stringifying, see the C<make> and C<emake>.
+
=head2 CHANGED IN PERL 5.6
The C<display_format> class method and the corresponding
C<i> (the imaginary unit). For the C<atan>, C<acoth>, the argument
cannot be C<-i> (the negative imaginary unit). For the C<tan>,
C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
-is any integer.
+is any integer. atan2(0, 0) is undefined, and if the complex arguments
+are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
Note that because we are operating on approximations of real numbers,
these errors can happen when merely `too close' to the singularities
=head1 BUGS
Saying C<use Math::Complex;> exports many mathematical routines in the
-caller environment and even overrides some (C<sqrt>, C<log>).
+caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
This is construed as a feature by the Authors, actually... ;-)
All routines expect to be given real or complex numbers. Don't attempt to
#!./perl
-# $RCSfile: complex.t,v $
#
# Regression tests for the Math::Complex pacakge
# -- Raphael Manfredi since Sep 1996
# -- Daniel S. Lewart since Sep 1997
BEGIN {
- chdir 't' if -d 't';
- @INC = '../lib';
+ if ($ENV{PERL_CORE}) {
+ chdir 't' if -d 't';
+ @INC = '../lib';
+ }
}
use Math::Complex;
$test = 0;
$| = 1;
my @script = (
- 'my ($res, $s0,$s1,$s2,$s3,$s4,$s5,$s6,$s7,$s8,$s9,$s10, $z0,$z1,$z2);' .
+ 'my ($res, $s0,$s1,$s2,$s3,$s4,$s5,$s6,$s7,$s8,$s9,$s10,$z0,$z1,$z2);' .
"\n\n"
);
my $eps = 1e-13;
my $pi = cplx(pi, 0);
my $pii = cplx(0, pi);
my $pip2 = cplx(pi/2, 0);
+my $pip4 = cplx(pi/4, 0);
my $zero = cplx(0, 0);
+my $inf = 9**9**9;
';
push(@script, $constants);
'coth(0)',
'csc(0)',
'csch(0)',
+ 'atan(cplx(0, 1), cplx(1, 0))',
);
test_loz(
$test++;
push @script, <<EOS;
print "# remake 2+3i\n";
- my \$z = cplx('2+3i');
+ \$z = cplx('2+3i');
print "not " unless \$z == Math::Complex->make(2,3);
print "ok $test\n";
EOS
$test++;
push @script, <<EOS;
- print "# remake 3i\n";
- my \$z = Math::Complex->make('3i');
+ print "# make 3i\n";
+ \$z = Math::Complex->make('3i');
print "not " unless \$z == cplx(0,3);
print "ok $test\n";
EOS
$test++;
push @script, <<EOS;
- print "# remake [2,3]\n";
- my \$z = cplxe('[2,3]');
- print "not " unless \$z == Math::Complex->emake(2,3);
+ print "# emake [2,3]\n";
+ \$z = Math::Complex->emake('[2,3]');
+ print "not " unless \$z == cplxe(2,3);
+ print "ok $test\n";
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print "# make (2,3)\n";
+ \$z = Math::Complex->make('(2,3)');
+ print "not " unless \$z == cplx(2,3);
+ print "ok $test\n";
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print "# emake [2,3pi/8]\n";
+ \$z = Math::Complex->emake('[2,3pi/8]');
+ print "not " unless \$z == cplxe(2,3*\$pi/8);
print "ok $test\n";
EOS
$test++;
push @script, <<EOS;
- print "# remake [2]\n";
- my \$z = Math::Complex->emake('[2]');
+ print "# emake [2]\n";
+ \$z = Math::Complex->emake('[2]');
print "not " unless \$z == cplxe(2);
print "ok $test\n";
EOS
}
+sub test_no_args {
+ push @script, <<'EOS';
+{
+ print "# cplx, cplxe, make, emake without arguments\n";
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ my \$z0 = cplx();
+ print ((\$z0->Re() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print ((\$z0->Im() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ my \$z1 = cplxe();
+ print ((\$z1->rho() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print ((\$z1->theta() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ my \$z2 = Math::Complex->make();
+ print ((\$z2->Re() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print ((\$z2->Im() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ my \$z3 = Math::Complex->emake();
+ print ((\$z3->rho() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+ $test++;
+ push @script, <<EOS;
+ print ((\$z3->theta() == 0) ? "ok $test\n" : "not ok $test\n");
+}
+EOS
+}
+
+sub test_atan2 {
+ push @script, <<'EOS';
+print "# atan2() with some real arguments\n";
+EOS
+ my @real = (-1, 0, 1);
+ for my $x (@real) {
+ for my $y (@real) {
+ next if $x == 0 && $y == 0;
+ $test++;
+ push @script, <<EOS;
+print ((Math::Complex::atan2($y, $x) == CORE::atan2($y, $x)) ? "ok $test\n" : "not ok $test\n");
+EOS
+ }
+ }
+ push @script, <<'EOS';
+ print "# atan2() with some complex arguments\n";
+EOS
+ $test++;
+ push @script, <<EOS;
+ print (abs(atan2(0, cplx(0, 1))) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+ $test++;
+ push @script, <<EOS;
+ print (abs(atan2(cplx(0, 1), 0) - \$pip2) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+ $test++;
+ push @script, <<EOS;
+ print (abs(atan2(cplx(0, 1), cplx(0, 1)) - \$pip4) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+ $test++;
+ push @script, <<EOS;
+ print (abs(atan2(cplx(0, 1), cplx(1, 1)) - cplx(0.553574358897045, 0.402359478108525)) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+}
+
+sub test_decplx {
+}
+
test_remake();
+test_no_args();
+
+test_atan2();
+
+test_decplx();
+
print "1..$test\n";
+#print @script, "\n";
eval join '', @script;
die $@ if $@;
||
($expected =~ /^-?\d/ && $got == $expected)
||
+ (abs(Math::Complex->make($got) - Math::Complex->make($expected)) < $eps)
+ ||
(abs($got - $expected) < $eps)
) {
print "ok $test\n";
|'(root(z, 4))[1] ** 4':'z'
|'(root(z, 5))[3] ** 5':'z'
|'(root(z, 8))[7] ** 8':'z'
+|'(root(z, 8, 0)) ** 8':'z'
+|'(root(z, 8, 7)) ** 8':'z'
|'abs(z)':'r'
|'acot(z)':'acotan(z)'
|'acsc(z)':'acosec(z)'
use 5.006;
use strict;
+use Math::Complex 1.35;
use Math::Complex qw(:trig);
our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
@ISA = qw(Exporter);
-$VERSION = 1.02;
+$VERSION = 1.03;
my @angcnv = qw(rad2deg rad2grad
deg2rad deg2grad
spherical_to_cartesian
spherical_to_cylindrical);
-@EXPORT_OK = (@rdlcnv, 'great_circle_distance', 'great_circle_direction');
+my @greatcircle = qw(
+ great_circle_distance
+ great_circle_direction
+ great_circle_bearing
+ great_circle_waypoint
+ great_circle_midpoint
+ great_circle_destination
+ );
-%EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
+my @pi = qw(pi2 pip2 pip4);
+
+@EXPORT_OK = (@rdlcnv, @greatcircle, @pi);
+
+# See e.g. the following pages:
+# http://www.movable-type.co.uk/scripts/LatLong.html
+# http://williams.best.vwh.net/avform.htm
+
+%EXPORT_TAGS = ('radial' => [ @rdlcnv ],
+ 'great_circle' => [ @greatcircle ],
+ 'pi' => [ @pi ]);
sub pi2 () { 2 * pi }
sub pip2 () { pi / 2 }
+sub pip4 () { pi / 4 }
sub DR () { pi2/360 }
sub RD () { 360/pi2 }
return rad2rad($direction);
}
+*great_circle_bearing = \&great_circle_direction;
+
+sub great_circle_waypoint {
+ my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
+
+ $point = 0.5 unless defined $point;
+
+ my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
+
+ return undef if $d == pi;
+
+ my $sd = sin($d);
+
+ return ($theta0, $phi0) if $sd == 0;
+
+ my $A = sin((1 - $point) * $d) / $sd;
+ my $B = sin( $point * $d) / $sd;
+
+ my $lat0 = pip2 - $phi0;
+ my $lat1 = pip2 - $phi1;
+
+ my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
+ my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
+ my $z = $A * sin($lat0) + $B * sin($lat1);
+
+ my $theta = atan2($y, $x);
+ my $phi = atan2($z, sqrt($x*$x + $y*$y));
+
+ return ($theta, $phi);
+}
+
+sub great_circle_midpoint {
+ great_circle_waypoint(@_[0..3], 0.5);
+}
+
+sub great_circle_destination {
+ my ( $theta0, $phi0, $dir0, $dst ) = @_;
+
+ my $lat0 = pip2 - $phi0;
+
+ my $phi1 = asin(sin($lat0)*cos($dst)+cos($lat0)*sin($dst)*cos($dir0));
+ my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
+ cos($dst)-sin($lat0)*sin($phi1));
+
+ my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
+
+ $dir1 -= pi2 if $dir1 > pi2;
+
+ return ($theta1, $phi1, $dir1);
+}
+
1;
__END__
$rad = deg2rad(120);
+ # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
+ use Math::Trig ':pi';
+
+ # Import the conversions between cartesian/spherical/cylindrical.
+ use Math::Trig ':radial';
+
+ # Import the great circle formulas.
+ use Math::Trig ':great_circle';
+
=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.
+conversions, and I<great circle formulas> for spherical movement.
=head1 TRIGONOMETRIC FUNCTIONS
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.
+pi>, where I<k> is any integer. atan2(0, 0) is undefined.
=head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
=head2 COORDINATE SYSTEMS
-B<Cartesian> coordinates are the usual rectangular I<(x, y,
-z)>-coordinates.
+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
$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:
+The direction you must follow the great circle (also known as I<bearing>)
+can be computed by the great_circle_direction() function:
use Math::Trig 'great_circle_direction';
$direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
+(Alias 'great_circle_bearing' is also available.)
The result is in radians, zero indicating straight north, pi or -pi
straight south, pi/2 straight west, and -pi/2 straight east.
+You can inversely compute the destination if you know the
+starting point, direction, and distance:
+
+ use Math::Trig 'great_circle_destination';
+
+ # thetad and phid are the destination coordinates,
+ # dird is the final direction at the destination.
+
+ ($thetad, $phid, $dird) =
+ great_circle_destination($theta, $phi, $direction, $distance);
+
+or the midpoint if you know the end points:
+
+ use Math::Trig 'great_circle_midpoint';
+
+ ($thetam, $phim) =
+ great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
+
+The great_circle_midpoint() is just a special case of
+
+ use Math::Trig 'great_circle_waypoint';
+
+ ($thetai, $phii) =
+ great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
+
+Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
+$phi1). Note that antipodal points (where their distance is I<pi>
+radians) do not have waypoints between them (they would have an an
+"equator" between them), and therefore C<undef> is returned for
+antipodal points. If the points are the same and the distance
+therefore zero and all waypoints therefore identical, the first point
+(either point) is returned.
+
+The thetas, phis, direction, and distance in the above are all in radians.
+
+You can import all the great circle formulas by
+
+ use Math::Trig ':great_circle';
+
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
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);
+ sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
+ my @L = NESW( -0.5, 51.3);
+ my @T = NESW(139.8, 35.7);
+ my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
-The direction you would have to go from London to Tokyo
+The direction you would have to go from London to Tokyo (in radians,
+straight north being zero, straight east being pi/2).
use Math::Trig qw(great_circle_direction);
- $rad = great_circle_direction(@L, @T);
+ my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
-=head2 CAVEAT FOR GREAT CIRCLE FORMULAS
+The midpoint between London and Tokyo being
-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
+ use Math::Trig qw(great_circle_midpoint);
+
+ my @M = great_circle_midpoint(@L, @T);
+
+or about 68.11N 24.74E, in the Finnish Lapland.
- lat0 = 90 degrees - phi0
- lat1 = 90 degrees - phi1
- d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
- sin(lat0) * sin(lat1))
+=head2 CAVEAT FOR GREAT CIRCLE FORMULAS
-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.
+The answers may be off by few percentages because of the irregular
+(slightly aspherical) form of the Earth. The errors are at worst
+about 0.55%, but generally below 0.3%.
=head1 BUGS
cannot be completely avoided if we want things like C<asin(2)> to give
an answer instead of giving a fatal runtime error.
+Do not attempt navigation using these formulas.
+
=head1 AUTHORS
Jarkko Hietaniemi <F<jhi@iki.fi>> and
#
# Regression tests for the Math::Trig package
#
-# The tests are quite modest as the Math::Complex tests exercise
-# these quite vigorously.
+# The tests here are quite modest as the Math::Complex tests exercise
+# these interfaces quite vigorously.
#
# -- Jarkko Hietaniemi, April 1997
BEGIN {
- chdir 't' if -d 't';
- @INC = '../lib';
+ if ($ENV{PERL_CORE}) {
+ chdir 't' if -d 't';
+ @INC = '../lib';
+ }
}
-use Math::Trig;
+use Math::Trig 1.03;
+
+my $pip2 = pi / 2;
use strict;
$_[1] ? (abs($_[0]/$_[1] - 1) < $e) : abs($_[0]) < $e;
}
-print "1..29\n";
+print "1..49\n";
$x = 0.9;
print 'not ' unless (near(tan($x), sin($x) / cos($x)));
unless (near(great_circle_direction(0, 0, 0, pi/2), pi));
print "ok 24\n";
-# Retired test: Relies on atan(0, 0), which is not portable.
+# Retired test: Relies on atan2(0, 0), which is not portable.
# print 'not '
# unless (near(great_circle_direction(0, 0, pi, pi), -pi()/2));
print "ok 25\n";
print 'not '
unless (near(rad2deg(great_circle_direction(@London, @Tokyo)),
31.791945393073));
-
print "ok 26\n";
+
print 'not '
unless (near(rad2deg(great_circle_direction(@Tokyo, @London)),
336.069766430326));
246.800348034667));
print "ok 28\n";
+
print 'not '
unless (near(rad2deg(great_circle_direction(@Paris, @Berlin)),
58.2079877553156));
print "ok 29\n";
+
+ use Math::Trig 'great_circle_bearing';
+
+ print 'not '
+ unless (near(rad2deg(great_circle_bearing(@Paris, @Berlin)),
+ 58.2079877553156));
+ print "ok 30\n";
+
+ use Math::Trig 'great_circle_waypoint';
+ use Math::Trig 'great_circle_midpoint';
+
+ my ($lon, $lat);
+
+ ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.0);
+
+ print 'not ' unless (near($lon, $London[0]));
+ print "ok 31\n";
+
+ print 'not ' unless (near($lat, $pip2 - $London[1]));
+ print "ok 32\n";
+
+ ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 1.0);
+
+ print 'not ' unless (near($lon, $Tokyo[0]));
+ print "ok 33\n";
+
+ print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
+ print "ok 34\n";
+
+ ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.5);
+
+ print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
+ print "ok 35\n";
+
+ print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
+ print "ok 36\n";
+
+ ($lon, $lat) = great_circle_midpoint(@London, @Tokyo);
+
+ print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
+ print "ok 37\n";
+
+ print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
+ print "ok 38\n";
+
+ ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.25);
+
+ print 'not ' unless (near($lon, 0.516073562850837)); # 29.5688 E
+ print "ok 39\n";
+
+ print 'not ' unless (near($lat, 1.170565013391510)); # 67.0684 N
+ print "ok 40\n";
+ ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.75);
+
+ print 'not ' unless (near($lon, 2.17494903805952)); # 124.6154 E
+ print "ok 41\n";
+
+ print 'not ' unless (near($lat, 0.952987032741305)); # 54.6021 N
+ print "ok 42\n";
+
+ use Math::Trig 'great_circle_destination';
+
+ my $dir1 = great_circle_direction(@London, @Tokyo);
+ my $dst1 = great_circle_distance(@London, @Tokyo);
+
+ ($lon, $lat) = great_circle_destination(@London, $dir1, $dst1);
+
+ print 'not ' unless (near($lon, $Tokyo[0]));
+ print "ok 43\n";
+
+ print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
+ print "ok 44\n";
+
+ my $dir2 = great_circle_direction(@Tokyo, @London);
+ my $dst2 = great_circle_distance(@Tokyo, @London);
+
+ ($lon, $lat) = great_circle_destination(@Tokyo, $dir2, $dst2);
+
+ print 'not ' unless (near($lon, $London[0]));
+ print "ok 45\n";
+
+ print 'not ' unless (near($lat, $pip2 - $London[1]));
+ print "ok 46\n";
+
+ my $dir3 = (great_circle_destination(@London, $dir1, $dst1))[2];
+
+ print 'not ' unless (near($dir3, 2.69379263839118)); # about 154.343 deg
+ print "ok 47\n";
+
+ my $dir4 = (great_circle_destination(@Tokyo, $dir2, $dst2))[2];
+
+ print 'not ' unless (near($dir4, 3.6993902625701)); # about 211.959 deg
+ print "ok 48\n";
+
+ print 'not ' unless (near($dst1, $dst2));
+ print "ok 49\n";
}
# eof
sub tan { sin($_[0]) / cos($_[0]) }
+Note that atan2(0, 0) is not well-defined.
+
=item bind SOCKET,NAME
Binds a network address to a socket, just as the bind system call
projects / p5sagit/p5-mst-13.2.git / commitdiff