From: Steve Peters Date: Thu, 6 Jul 2006 15:18:58 +0000 (+0000) Subject: Math-Complex now dual-lived with Jarkko Hietaniemi as the maintainer. X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=affad850140455cfdbfe87c7519ca10909a8bf23;p=p5sagit%2Fp5-mst-13.2.git Math-Complex now dual-lived with Jarkko Hietaniemi as the maintainer. Also, sync'ing up with the CPAN version of the module. p4raw-id: //depot/perl@28494 --- diff --git a/Porting/Maintainers.pl b/Porting/Maintainers.pl index 427fdef..d671976 100644 --- a/Porting/Maintainers.pl +++ b/Porting/Maintainers.pl @@ -380,6 +380,13 @@ package Maintainers; 'CPAN' => 1, }, + 'Math::Complex' => + { + 'MAINTAINER' => 'jhi', + 'FILES' => q[lib/Math/Complex.pm lib/Math/Trig.pm], + 'CPAN' => 1, + }, + 'Memoize' => { 'MAINTAINER' => 'mjd', diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm index e4b980b..110e8b6 100644 --- a/lib/Math/Complex.pm +++ b/lib/Math/Complex.pm @@ -9,7 +9,7 @@ package Math::Complex; use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf); -$VERSION = 1.35; +$VERSION = 1.36; BEGIN { unless ($^O eq 'unicosmk') { @@ -66,22 +66,25 @@ my @trig = qw( ), @trig); -@EXPORT_OK = qw(decplx); +my @pi = qw(pi pi2 pi4 pip2 pip4); + +@EXPORT_OK = @pi; %EXPORT_TAGS = ( 'trig' => [@trig], + 'pi' => [@pi], ); use overload - '+' => \&plus, - '-' => \&minus, - '*' => \&multiply, - '/' => \÷, - '**' => \&power, - '==' => \&numeq, - '<=>' => \&spaceship, - 'neg' => \&negate, - '~' => \&conjugate, + '+' => \&_plus, + '-' => \&_minus, + '*' => \&_multiply, + '/' => \&_divide, + '**' => \&_power, + '==' => \&_numeq, + '<=>' => \&_spaceship, + 'neg' => \&_negate, + '~' => \&_conjugate, 'abs' => \&abs, 'sqrt' => \&sqrt, 'exp' => \&exp, @@ -90,7 +93,7 @@ use overload 'cos' => \&cos, 'tan' => \&tan, 'atan2' => \&atan2, - qw("" stringify); + '""' => \&_stringify; # # Package "privates" @@ -107,6 +110,8 @@ my $eps = 1e-14; # Epsilon # c_dirty cartesian form not up-to-date # p_dirty polar form not up-to-date # display display format (package's global when not set) +# bn_cartesian +# bnc_dirty # # Die on bad *make() arguments. @@ -183,7 +188,7 @@ sub make { } $im ||= 0; _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/; - $self->set_cartesian([$re, $im ]); + $self->_set_cartesian([$re, $im ]); $self->display_format('cartesian'); return $self; @@ -217,7 +222,7 @@ sub emake { } $theta ||= 0; _cannot_make("theta", $theta) unless $theta =~ /^$gre$/; - $self->set_polar([$rho, $theta]); + $self->_set_polar([$rho, $theta]); $self->display_format('polar'); return $self; @@ -253,11 +258,18 @@ sub cplxe { sub pi () { 4 * CORE::atan2(1, 1) } # -# pit2 +# pi2 # # The full circle # -sub pit2 () { 2 * pi } +sub pi2 () { 2 * pi } + +# +# pi4 +# +# The full circle twice. +# +sub pi4 () { 4 * pi } # # pip2 @@ -267,19 +279,18 @@ sub pit2 () { 2 * pi } sub pip2 () { pi / 2 } # -# deg1 +# pip4 # -# One degree in radians, used in stringify_polar. +# The eighth circle. # - -sub deg1 () { pi / 180 } +sub pip4 () { pi / 4 } # -# uplog10 +# _uplog10 # # Used in log10(). # -sub uplog10 () { 1 / CORE::log(10) } +sub _uplog10 () { 1 / CORE::log(10) } # # i @@ -297,32 +308,32 @@ sub i () { } # -# ip2 +# _ip2 # # Half of i. # -sub ip2 () { i / 2 } +sub _ip2 () { i / 2 } # # Attribute access/set routines # -sub cartesian {$_[0]->{c_dirty} ? - $_[0]->update_cartesian : $_[0]->{'cartesian'}} -sub polar {$_[0]->{p_dirty} ? - $_[0]->update_polar : $_[0]->{'polar'}} +sub _cartesian {$_[0]->{c_dirty} ? + $_[0]->_update_cartesian : $_[0]->{'cartesian'}} +sub _polar {$_[0]->{p_dirty} ? + $_[0]->_update_polar : $_[0]->{'polar'}} -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] } +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 +# ->_update_cartesian # # Recompute and return the cartesian form, given accurate polar form. # -sub update_cartesian { +sub _update_cartesian { my $self = shift; my ($r, $t) = @{$self->{'polar'}}; $self->{c_dirty} = 0; @@ -331,11 +342,11 @@ sub update_cartesian { # # -# ->update_polar +# ->_update_polar # # Recompute and return the polar form, given accurate cartesian form. # -sub update_polar { +sub _update_polar { my $self = shift; my ($x, $y) = @{$self->{'cartesian'}}; $self->{p_dirty} = 0; @@ -345,34 +356,34 @@ sub update_polar { } # -# (plus) +# (_plus) # # Computes z1+z2. # -sub plus { +sub _plus { my ($z1, $z2, $regular) = @_; - my ($re1, $im1) = @{$z1->cartesian}; + my ($re1, $im1) = @{$z1->_cartesian}; $z2 = cplx($z2) unless ref $z2; - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); unless (defined $regular) { - $z1->set_cartesian([$re1 + $re2, $im1 + $im2]); + $z1->_set_cartesian([$re1 + $re2, $im1 + $im2]); return $z1; } return (ref $z1)->make($re1 + $re2, $im1 + $im2); } # -# (minus) +# (_minus) # # Computes z1-z2. # -sub minus { +sub _minus { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = @{$z1->cartesian}; + my ($re1, $im1) = @{$z1->_cartesian}; $z2 = cplx($z2) unless ref $z2; - my ($re2, $im2) = @{$z2->cartesian}; + my ($re2, $im2) = @{$z2->_cartesian}; unless (defined $inverted) { - $z1->set_cartesian([$re1 - $re2, $im1 - $im2]); + $z1->_set_cartesian([$re1 - $re2, $im1 - $im2]); return $z1; } return $inverted ? @@ -382,28 +393,28 @@ sub minus { } # -# (multiply) +# (_multiply) # # Computes z1*z2. # -sub multiply { +sub _multiply { my ($z1, $z2, $regular) = @_; if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { # if both polar better use polar to avoid rounding errors - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = @{$z2->polar}; + my ($r1, $t1) = @{$z1->_polar}; + my ($r2, $t2) = @{$z2->_polar}; my $t = $t1 + $t2; - if ($t > pi()) { $t -= pit2 } - elsif ($t <= -pi()) { $t += pit2 } + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } unless (defined $regular) { - $z1->set_polar([$r1 * $r2, $t]); + $z1->_set_polar([$r1 * $r2, $t]); return $z1; } return (ref $z1)->emake($r1 * $r2, $t); } else { - my ($x1, $y1) = @{$z1->cartesian}; + my ($x1, $y1) = @{$z1->_cartesian}; if (ref $z2) { - my ($x2, $y2) = @{$z2->cartesian}; + my ($x2, $y2) = @{$z2->_cartesian}; return (ref $z1)->make($x1*$x2-$y1*$y2, $x1*$y2+$y1*$x2); } else { return (ref $z1)->make($x1*$z2, $y1*$z2); @@ -433,41 +444,41 @@ sub _divbyzero { } # -# (divide) +# (_divide) # # Computes z1/z2. # -sub divide { +sub _divide { my ($z1, $z2, $inverted) = @_; if ($z1->{p_dirty} == 0 and ref $z2 and $z2->{p_dirty} == 0) { # if both polar better use polar to avoid rounding errors - my ($r1, $t1) = @{$z1->polar}; - my ($r2, $t2) = @{$z2->polar}; + my ($r1, $t1) = @{$z1->_polar}; + my ($r2, $t2) = @{$z2->_polar}; my $t; if ($inverted) { _divbyzero "$z2/0" if ($r1 == 0); $t = $t2 - $t1; - if ($t > pi()) { $t -= pit2 } - elsif ($t <= -pi()) { $t += pit2 } + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } return (ref $z1)->emake($r2 / $r1, $t); } else { _divbyzero "$z1/0" if ($r2 == 0); $t = $t1 - $t2; - if ($t > pi()) { $t -= pit2 } - elsif ($t <= -pi()) { $t += pit2 } + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } return (ref $z1)->emake($r1 / $r2, $t); } } else { my ($d, $x2, $y2); if ($inverted) { - ($x2, $y2) = @{$z1->cartesian}; + ($x2, $y2) = @{$z1->_cartesian}; $d = $x2*$x2 + $y2*$y2; _divbyzero "$z2/0" if $d == 0; return (ref $z1)->make(($x2*$z2)/$d, -($y2*$z2)/$d); } else { - my ($x1, $y1) = @{$z1->cartesian}; + my ($x1, $y1) = @{$z1->_cartesian}; if (ref $z2) { - ($x2, $y2) = @{$z2->cartesian}; + ($x2, $y2) = @{$z2->_cartesian}; $d = $x2*$x2 + $y2*$y2; _divbyzero "$z1/0" if $d == 0; my $u = ($x1*$x2 + $y1*$y2)/$d; @@ -482,11 +493,11 @@ sub divide { } # -# (power) +# (_power) # # Computes z1**z2 = exp(z2 * log z1)). # -sub power { +sub _power { my ($z1, $z2, $inverted) = @_; if ($inverted) { return 1 if $z1 == 0 || $z2 == 1; @@ -500,65 +511,65 @@ sub power { # If both arguments cartesian, return cartesian, else polar. return $z1->{c_dirty} == 0 && (not ref $z2 or $z2->{c_dirty} == 0) ? - cplx(@{$w->cartesian}) : $w; + cplx(@{$w->_cartesian}) : $w; } # -# (spaceship) +# (_spaceship) # # Computes z1 <=> z2. # Sorts on the real part first, then on the imaginary part. Thus 2-4i < 3+8i. # -sub spaceship { +sub _spaceship { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); my $sgn = $inverted ? -1 : 1; return $sgn * ($re1 <=> $re2) if $re1 != $re2; return $sgn * ($im1 <=> $im2); } # -# (numeq) +# (_numeq) # # Computes z1 == z2. # -# (Required in addition to spaceship() because of NaNs.) -sub numeq { +# (Required in addition to _spaceship() because of NaNs.) +sub _numeq { my ($z1, $z2, $inverted) = @_; - my ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + my ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + my ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); return $re1 == $re2 && $im1 == $im2 ? 1 : 0; } # -# (negate) +# (_negate) # # Computes -z. # -sub negate { +sub _negate { my ($z) = @_; if ($z->{c_dirty}) { - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; $t = ($t <= 0) ? $t + pi : $t - pi; return (ref $z)->emake($r, $t); } - my ($re, $im) = @{$z->cartesian}; + my ($re, $im) = @{$z->_cartesian}; return (ref $z)->make(-$re, -$im); } # -# (conjugate) +# (_conjugate) # -# Compute complex's conjugate. +# Compute complex's _conjugate. # -sub conjugate { +sub _conjugate { my ($z) = @_; if ($z->{c_dirty}) { - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; return (ref $z)->emake($r, -$t); } - my ($re, $im) = @{$z->cartesian}; + my ($re, $im) = @{$z->_cartesian}; return (ref $z)->make($re, -$im); } @@ -577,20 +588,20 @@ sub abs { } } if (defined $rho) { - $z->{'polar'} = [ $rho, ${$z->polar}[1] ]; + $z->{'polar'} = [ $rho, ${$z->_polar}[1] ]; $z->{p_dirty} = 0; $z->{c_dirty} = 1; return $rho; } else { - return ${$z->polar}[0]; + return ${$z->_polar}[0]; } } sub _theta { my $theta = $_[0]; - if ($$theta > pi()) { $$theta -= pit2 } - elsif ($$theta <= -pi()) { $$theta += pit2 } + if ($$theta > pi()) { $$theta -= pi2 } + elsif ($$theta <= -pi()) { $$theta += pi2 } } # @@ -603,11 +614,11 @@ sub arg { return $z unless ref $z; if (defined $theta) { _theta(\$theta); - $z->{'polar'} = [ ${$z->polar}[0], $theta ]; + $z->{'polar'} = [ ${$z->_polar}[0], $theta ]; $z->{p_dirty} = 0; $z->{c_dirty} = 1; } else { - $theta = ${$z->polar}[1]; + $theta = ${$z->_polar}[1]; _theta(\$theta); } return $theta; @@ -630,10 +641,10 @@ sub arg { # sub sqrt { my ($z) = @_; - my ($re, $im) = ref $z ? @{$z->cartesian} : ($z, 0); + my ($re, $im) = ref $z ? @{$z->_cartesian} : ($z, 0); return $re < 0 ? cplx(0, CORE::sqrt(-$re)) : CORE::sqrt($re) if $im == 0; - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; return (ref $z)->emake(CORE::sqrt($r), $t/2); } @@ -650,7 +661,7 @@ sub cbrt { -CORE::exp(CORE::log(-$z)/3) : ($z > 0 ? CORE::exp(CORE::log($z)/3): 0) unless ref $z; - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; return 0 if $r == 0; return (ref $z)->emake(CORE::exp(CORE::log($r)/3), $t/3); } @@ -684,8 +695,8 @@ sub root { 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 $theta_inc = pit2 / $n; + @{$z->_polar} : (CORE::abs($z), $z >= 0 ? 0 : pi); + my $theta_inc = pi2 / $n; my $rho = $r ** (1/$n); my $cartesian = ref $z && $z->{c_dirty} == 0; if (@_ == 2) { @@ -695,12 +706,12 @@ sub root { $i++, $theta += $theta_inc) { my $w = cplxe($rho, $theta); # Yes, $cartesian is loop invariant. - push @root, $cartesian ? cplx(@{$w->cartesian}) : $w; + 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 $cartesian ? cplx(@{$w->_cartesian}) : $w; } } @@ -713,11 +724,11 @@ sub Re { my ($z, $Re) = @_; return $z unless ref $z; if (defined $Re) { - $z->{'cartesian'} = [ $Re, ${$z->cartesian}[1] ]; + $z->{'cartesian'} = [ $Re, ${$z->_cartesian}[1] ]; $z->{c_dirty} = 0; $z->{p_dirty} = 1; } else { - return ${$z->cartesian}[0]; + return ${$z->_cartesian}[0]; } } @@ -730,11 +741,11 @@ sub Im { my ($z, $Im) = @_; return 0 unless ref $z; if (defined $Im) { - $z->{'cartesian'} = [ ${$z->cartesian}[0], $Im ]; + $z->{'cartesian'} = [ ${$z->_cartesian}[0], $Im ]; $z->{c_dirty} = 0; $z->{p_dirty} = 1; } else { - return ${$z->cartesian}[1]; + return ${$z->_cartesian}[1]; } } @@ -763,7 +774,7 @@ sub theta { # sub exp { my ($z) = @_; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; return (ref $z)->emake(CORE::exp($x), $y); } @@ -799,10 +810,10 @@ sub log { _logofzero("log") if $z == 0; return $z > 0 ? CORE::log($z) : cplx(CORE::log(-$z), pi); } - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; _logofzero("log") if $r == 0; - if ($t > pi()) { $t -= pit2 } - elsif ($t <= -pi()) { $t += pit2 } + if ($t > pi()) { $t -= pi2 } + elsif ($t <= -pi()) { $t += pi2 } return (ref $z)->make(CORE::log($r), $t); } @@ -820,7 +831,7 @@ sub ln { Math::Complex::log(@_) } # sub log10 { - return Math::Complex::log($_[0]) * uplog10; + return Math::Complex::log($_[0]) * _uplog10; } # @@ -844,7 +855,7 @@ sub logn { sub cos { my ($z) = @_; return CORE::cos($z) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; my $ey = CORE::exp($y); my $sx = CORE::sin($x); my $cx = CORE::cos($x); @@ -861,7 +872,7 @@ sub cos { sub sin { my ($z) = @_; return CORE::sin($z) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; my $ey = CORE::exp($y); my $sx = CORE::sin($x); my $cx = CORE::cos($x); @@ -942,7 +953,7 @@ sub acos { return CORE::atan2(CORE::sqrt(1-$z*$z), $z) if (! ref $z) && CORE::abs($z) <= 1; $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; return 0 if $x == 1 && $y == 0; my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); @@ -967,7 +978,7 @@ sub asin { return CORE::atan2($z, CORE::sqrt(1-$z*$z)) if (! ref $z) && CORE::abs($z) <= 1; $z = cplx($z, 0) unless ref $z; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; return 0 if $x == 0 && $y == 0; my $t1 = CORE::sqrt(($x+1)*($x+1) + $y*$y); my $t2 = CORE::sqrt(($x-1)*($x-1) + $y*$y); @@ -990,12 +1001,12 @@ sub asin { sub atan { my ($z) = @_; return CORE::atan2($z, 1) unless ref $z; - my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0); + my ($x, $y) = ref $z ? @{$z->_cartesian} : ($z, 0); return 0 if $x == 0 && $y == 0; _divbyzero "atan(i)" if ( $z == i); _logofzero "atan(-i)" if (-$z == i); # -i is a bad file test... my $log = &log((i + $z) / (i - $z)); - return ip2 * $log; + return _ip2 * $log; } # @@ -1061,7 +1072,7 @@ sub cosh { $ex = CORE::exp($z); return $ex ? ($ex + 1/$ex)/2 : $Inf; } - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; $ex = CORE::exp($x); my $ex_1 = $ex ? 1 / $ex : $Inf; return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2, @@ -1081,7 +1092,7 @@ sub sinh { $ex = CORE::exp($z); return $ex ? ($ex - 1/$ex)/2 : "-$Inf"; } - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; my $cy = CORE::cos($y); my $sy = CORE::sin($y); $ex = CORE::exp($x); @@ -1162,7 +1173,7 @@ sub acosh { unless (ref $z) { $z = cplx($z, 0); } - my ($re, $im) = @{$z->cartesian}; + my ($re, $im) = @{$z->_cartesian}; if ($im == 0) { return CORE::log($re + CORE::sqrt($re*$re - 1)) if $re >= 1; @@ -1278,11 +1289,11 @@ sub atan2 { my ($z1, $z2, $inverted) = @_; my ($re1, $im1, $re2, $im2); if ($inverted) { - ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); - ($re2, $im2) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); + ($re1, $im1) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); + ($re2, $im2) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); } else { - ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0); - ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0); + ($re1, $im1) = ref $z1 ? @{$z1->_cartesian} : ($z1, 0); + ($re2, $im2) = ref $z2 ? @{$z2->_cartesian} : ($z2, 0); } if ($im1 || $im2) { # In MATLAB the imaginary parts are ignored. @@ -1345,7 +1356,7 @@ sub display_format { } # -# (stringify) +# (_stringify) # # Show nicely formatted complex number under its cartesian or polar form, # depending on the current display format: @@ -1354,25 +1365,25 @@ sub display_format { # . Otherwise, use the generic current default for all complex numbers, # which is a package global variable. # -sub stringify { +sub _stringify { my ($z) = shift; my $style = $z->display_format; $style = $DISPLAY_FORMAT{style} unless defined $style; - return $z->stringify_polar if $style =~ /^p/i; - return $z->stringify_cartesian; + return $z->_stringify_polar if $style =~ /^p/i; + return $z->_stringify_cartesian; } # -# ->stringify_cartesian +# ->_stringify_cartesian # # Stringify as a cartesian representation 'a+bi'. # -sub stringify_cartesian { +sub _stringify_cartesian { my $z = shift; - my ($x, $y) = @{$z->cartesian}; + my ($x, $y) = @{$z->_cartesian}; my ($re, $im); my %format = $z->display_format; @@ -1428,13 +1439,13 @@ sub stringify_cartesian { # -# ->stringify_polar +# ->_stringify_polar # # Stringify as a polar representation '[r,t]'. # -sub stringify_polar { +sub _stringify_polar { my $z = shift; - my ($r, $t) = @{$z->polar}; + my ($r, $t) = @{$z->_polar}; my $theta; my %format = $z->display_format; @@ -1454,7 +1465,7 @@ sub stringify_polar { # Try to identify pi/n and friends. # - $t -= int(CORE::abs($t) / pit2) * pit2; + $t -= int(CORE::abs($t) / pi2) * pi2; if ($format{polar_pretty_print} && $t) { my ($a, $b); @@ -1571,14 +1582,14 @@ which is also expressed by this formula: In other words, it's the projection of the vector onto the I and I axes. Mathematicians call I the I or I and I -the I of the complex number. The I of C will be -noted C. +the I of the complex number. The I of C is +marked here as C. -The polar notation (also known as the trigonometric -representation) is much more handy for performing multiplications and -divisions of complex numbers, whilst the cartesian notation is better -suited for additions and subtractions. Real numbers are on the I -axis, and therefore I is zero or I. +The polar notation (also known as the trigonometric representation) is +much more handy for performing multiplications and divisions of +complex numbers, whilst the cartesian notation is better suited for +additions and subtractions. Real numbers are on the I axis, and +therefore I or I is zero or I. All the common operations that can be performed on a real number have been defined to work on complex numbers as well, and are merely @@ -1675,6 +1686,8 @@ The definition used for complex arguments of atan2() is -i log((x + iy)/sqrt(x*x+y*y)) +Note that atan2(0, 0) is not well-defined. + The following extra operations are supported on both real and complex numbers: @@ -1895,6 +1908,15 @@ Here are some examples: $j->arg(2); # (the last two aka rho, theta) # can be used also as mutators. +=head2 PI + +The constant C and some handy multiples of it (pi2, pi4, +and pip2 (pi/2) and pip4 (pi/4)) are also available if separately +exported: + + use Math::Complex ':pi'; + $third_of_circle = pi2 / 3; + =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO The division (/) and the following functions @@ -1961,10 +1983,9 @@ Whatever it is, it does not manifest itself anywhere else where Perl runs. =head1 AUTHORS -Daniel S. Lewart > - -Original authors Raphael Manfredi > and -Jarkko Hietaniemi > +Daniel S. Lewart > +Jarkko Hietaniemi > +Raphael Manfredi > =cut diff --git a/lib/Math/Complex.t b/lib/Math/Complex.t index b82350f..d14259d 100755 --- a/lib/Math/Complex.t +++ b/lib/Math/Complex.t @@ -497,14 +497,14 @@ sub test { # check the op= works push @script, <cartesian} : (\$z0, 0)); + my \$za = cplx(ref \$z0 ? \@{\$z0->_cartesian} : (\$z0, 0)); - my (\$z1r, \$z1i) = ref \$z1 ? \@{\$z1->cartesian} : (\$z1, 0); + my (\$z1r, \$z1i) = ref \$z1 ? \@{\$z1->_cartesian} : (\$z1, 0); my \$zb = cplx(\$z1r, \$z1i); \$za $op= \$zb; - my (\$zbr, \$zbi) = \@{\$zb->cartesian}; + my (\$zbr, \$zbi) = \@{\$zb->_cartesian}; check($test, '\$z0 $op= \$z1', \$za, \$z$#args, $args); EOB @@ -684,7 +684,7 @@ __END__ (-100,0):(0,10) (16,-30):(5,-3) -&stringify_cartesian +&_stringify_cartesian (-100,0):"-100" (0,1):"i" (4,-3):"4-3i" @@ -693,7 +693,7 @@ __END__ (-2,4):"-2+4i" (-2,-1):"-2-i" -&stringify_polar +&_stringify_polar [-1, 0]:"[1,pi]" [1, pi/3]:"[1,pi/3]" [6, -2*pi/3]:"[6,-2pi/3]" diff --git a/lib/Math/Trig.pm b/lib/Math/Trig.pm index cd17356..37b3b75 100644 --- a/lib/Math/Trig.pm +++ b/lib/Math/Trig.pm @@ -7,17 +7,17 @@ require Exporter; package Math::Trig; -use 5.006; +use 5.005; use strict; -use Math::Complex 1.35; -use Math::Complex qw(:trig); +use Math::Complex 1.36; +use Math::Complex qw(:trig :pi); -our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS); +use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS); @ISA = qw(Exporter); -$VERSION = 1.03; +$VERSION = 1.04; my @angcnv = qw(rad2deg rad2grad deg2rad deg2grad @@ -42,7 +42,7 @@ my @greatcircle = qw( great_circle_destination ); -my @pi = qw(pi2 pip2 pip4); +my @pi = qw(pi pi2 pi4 pip2 pip4); @EXPORT_OK = (@rdlcnv, @greatcircle, @pi); @@ -54,22 +54,18 @@ my @pi = qw(pi2 pip2 pip4); '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 } -sub DG () { 400/360 } -sub GD () { 360/400 } -sub RG () { 400/pi2 } -sub GR () { pi2/400 } +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 ($$) { +sub _remt ($$) { # Oh yes, POSIX::fmod() would be faster. Possibly. If it is available. $_[0] - $_[1] * int($_[0] / $_[1]); } @@ -78,23 +74,23 @@ sub remt ($$) { # Angle conversions. # -sub rad2rad($) { remt($_[0], pi2) } +sub rad2rad($) { _remt($_[0], pi2) } -sub deg2deg($) { remt($_[0], 360) } +sub deg2deg($) { _remt($_[0], 360) } -sub grad2grad($) { remt($_[0], 400) } +sub grad2grad($) { _remt($_[0], 400) } -sub rad2deg ($;$) { my $d = RD * $_[0]; $_[1] ? $d : deg2deg($d) } +sub rad2deg ($;$) { my $d = _RD * $_[0]; $_[1] ? $d : deg2deg($d) } -sub deg2rad ($;$) { my $d = DR * $_[0]; $_[1] ? $d : rad2rad($d) } +sub deg2rad ($;$) { my $d = _DR * $_[0]; $_[1] ? $d : rad2rad($d) } -sub grad2deg ($;$) { my $d = GD * $_[0]; $_[1] ? $d : deg2deg($d) } +sub grad2deg ($;$) { my $d = _GD * $_[0]; $_[1] ? $d : deg2deg($d) } -sub deg2grad ($;$) { my $d = DG * $_[0]; $_[1] ? $d : grad2grad($d) } +sub deg2grad ($;$) { my $d = _DG * $_[0]; $_[1] ? $d : grad2grad($d) } -sub rad2grad ($;$) { my $d = RG * $_[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 grad2rad ($;$) { my $d = _GR * $_[0]; $_[1] ? $d : rad2rad($d) } sub cartesian_to_spherical { my ( $x, $y, $z ) = @_; @@ -229,24 +225,24 @@ Math::Trig - trigonometric functions =head1 SYNOPSIS - use Math::Trig; + use Math::Trig; - $x = tan(0.9); - $y = acos(3.7); - $z = asin(2.4); + $x = tan(0.9); + $y = acos(3.7); + $z = asin(2.4); - $halfpi = pi/2; + $halfpi = pi/2; - $rad = deg2rad(120); + $rad = deg2rad(120); - # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4). - use Math::Trig ':pi'; + # 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 conversions between cartesian/spherical/cylindrical. + use Math::Trig ':radial'; # Import the great circle formulas. - use Math::Trig ':great_circle'; + use Math::Trig ':great_circle'; =head1 DESCRIPTION @@ -280,7 +276,7 @@ The principal value of the arc tangent of y/x B(y, x) The arcus cofunctions of the sine, cosine, and tangent (acosec/acsc -and acotan/acot are aliases) +and acotan/acot are aliases). Note that atan2(0, 0) is not well-defined. B, B, B, B, B @@ -303,48 +299,51 @@ The arcus cofunctions of the hyperbolic sine, cosine, and tangent B, B, B, B, B -The trigonometric constant B is also defined. +The trigonometric constant B and some of handy multiples +of it are also defined. -$pi2 = 2 * B; +B =head2 ERRORS DUE TO DIVISION BY ZERO The following functions - acoth - acsc - acsch - asec - asech - atanh - cot - coth - csc - csch - sec - sech - tan - tanh + 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 ... + 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... + atanh(-1): Logarithm of zero. + Died at... For the C, C, C, C, C, C, C, C, C, the argument cannot be C<0> (zero). For the C, C, the argument cannot be C<1> (one). For the C, C, the argument cannot be C<-1> (minus one). For the C, C, C, C, the argument cannot be I, where I is any integer. atan2(0, 0) is undefined. +pi>, where I is any integer. + +Note that atan2(0, 0) is not well-defined. =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS @@ -364,11 +363,11 @@ for more information. In practice you need not to worry about getting complex numbers as results because the C takes care of details like for example how to display complex numbers. For example: - print asin(2), "\n"; + print asin(2), "\n"; should produce something like this (take or leave few last decimals): - 1.5707963267949-1.31695789692482i + 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>. @@ -377,25 +376,60 @@ and the imaginary part of approximately C<-1.317>. (Plane, 2-dimensional) angles may be converted with the following functions. - $radians = deg2rad($degrees); - $radians = grad2rad($gradians); +=over + +=item deg2rad + + $radians = deg2rad($degrees); + +=item grad2rad + + $radians = grad2rad($gradians); + +=item rad2deg + + $degrees = rad2deg($radians); - $degrees = rad2deg($radians); - $degrees = grad2deg($gradians); +=item grad2deg - $gradians = deg2grad($degrees); - $gradians = rad2grad($radians); + $degrees = grad2deg($gradians); + +=item deg2grad + + $gradians = deg2grad($degrees); + +=item rad2grad + + $gradians = rad2grad($radians); + +=back The full circle is 2 I 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); + $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(). +=over 4 + +=item rad2rad + + $radians_wrapped_by_2pi = rad2rad($radians); + +=item deg2deg + + $degrees_wrapped_by_360 = deg2deg($degrees); + +=item grad2grad + + $gradians_wrapped_by_400 = grad2grad($gradians); + +=back + =head1 RADIAL COORDINATE CONVERSIONS B are the B and the B @@ -454,29 +488,29 @@ I<-pi> angles. =item cartesian_to_cylindrical - ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); + ($rho, $theta, $z) = cartesian_to_cylindrical($x, $y, $z); =item cartesian_to_spherical - ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); + ($rho, $theta, $phi) = cartesian_to_spherical($x, $y, $z); =item cylindrical_to_cartesian - ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); + ($x, $y, $z) = cylindrical_to_cartesian($rho, $theta, $z); =item cylindrical_to_spherical - ($rho_s, $theta, $phi) = cylindrical_to_spherical($rho_c, $theta, $z); + ($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); + ($x, $y, $z) = spherical_to_cartesian($rho, $theta, $phi); =item spherical_to_cylindrical - ($rho_c, $theta, $z) = spherical_to_cylindrical($rho_s, $theta, $phi); + ($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>. @@ -484,6 +518,13 @@ Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>. =head1 GREAT CIRCLE DISTANCES AND DIRECTIONS +A great circle is section of a circle that contains the circle +diameter: the shortest distance between two (non-antipodal) points on +the spherical surface goes along the great circle connecting those two +points. + +=head2 great_circle_distance + You can compute spherical distances, called B, by importing the great_circle_distance() function: @@ -508,6 +549,8 @@ degrees). $distance = great_circle_distance($lon0, pi/2 - $lat0, $lon1, pi/2 - $lat1, $rho); +=head2 great_circle_direction + The direction you must follow the great circle (also known as I) can be computed by the great_circle_direction() function: @@ -515,13 +558,23 @@ can be computed by the great_circle_direction() function: $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. +=head2 great_circle_bearing + +Alias 'great_circle_bearing' for 'great_circle_direction' is also available. + + use Math::Trig 'great_circle_bearing'; + + $direction = great_circle_bearing($theta0, $phi0, $theta1, $phi1); + +The result of great_circle_direction 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: +=head2 great_circle_destination + use Math::Trig 'great_circle_destination'; # thetad and phid are the destination coordinates, @@ -532,6 +585,8 @@ starting point, direction, and distance: or the midpoint if you know the end points: +=head2 great_circle_midpoint + use Math::Trig 'great_circle_midpoint'; ($thetam, $phim) = @@ -539,6 +594,8 @@ or the midpoint if you know the end points: The great_circle_midpoint() is just a special case of +=head2 great_circle_waypoint + use Math::Trig 'great_circle_waypoint'; ($thetai, $phii) = @@ -569,26 +626,26 @@ Asia do often cross the polar regions. 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); + use Math::Trig qw(great_circle_distance deg2rad); - # Notice the 90 - latitude: phi zero is at the North Pole. - 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. + # Notice the 90 - latitude: phi zero is at the North Pole. + 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 (in radians, straight north being zero, straight east being pi/2). - use Math::Trig qw(great_circle_direction); + use Math::Trig qw(great_circle_direction); - my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi. + my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi. The midpoint between London and Tokyo being - use Math::Trig qw(great_circle_midpoint); + use Math::Trig qw(great_circle_midpoint); - my @M = great_circle_midpoint(@L, @T); + my @M = great_circle_midpoint(@L, @T); or about 68.11N 24.74E, in the Finnish Lapland. @@ -614,8 +671,8 @@ Do not attempt navigation using these formulas. =head1 AUTHORS -Jarkko Hietaniemi > and -Raphael Manfredi >. +Jarkko Hietaniemi > and +Raphael Manfredi >. =cut diff --git a/lib/Math/Trig.t b/lib/Math/Trig.t index 4777240..65bb796 100755 --- a/lib/Math/Trig.t +++ b/lib/Math/Trig.t @@ -15,6 +15,19 @@ BEGIN { } } +BEGIN { + eval { require Test::More }; + if ($@) { + # We are willing to lose testing in e.g. 5.00504. + print "1..0 # No Test::More, skipping\n"; + exit(0); + } else { + import Test::More; + } +} + +plan(tests => 69); + use Math::Trig 1.03; my $pip2 = pi / 2; @@ -31,127 +44,105 @@ if ($^O eq 'unicos') { # See lib/Math/Complex.pm and t/lib/complex.t. sub near ($$;$) { my $e = defined $_[2] ? $_[2] : $eps; - print "# near? $_[0] $_[1] $e\n"; - $_[1] ? (abs($_[0]/$_[1] - 1) < $e) : abs($_[0]) < $e; + my $d = $_[1] ? abs($_[0]/$_[1] - 1) : abs($_[0]); + print "# near? $_[0] $_[1] : $d : $e\n"; + $_[1] ? ($d < $e) : abs($_[0]) < $e; } -print "1..49\n"; - $x = 0.9; -print 'not ' unless (near(tan($x), sin($x) / cos($x))); -print "ok 1\n"; +ok(near(tan($x), sin($x) / cos($x))); -print 'not ' unless (near(sinh(2), 3.62686040784702)); -print "ok 2\n"; +ok(near(sinh(2), 3.62686040784702)); -print 'not ' unless (near(acsch(0.1), 2.99822295029797)); -print "ok 3\n"; +ok(near(acsch(0.1), 2.99822295029797)); $x = asin(2); -print 'not ' unless (ref $x eq 'Math::Complex'); -print "ok 4\n"; +is(ref $x, 'Math::Complex'); # avoid using Math::Complex here $x =~ /^([^-]+)(-[^i]+)i$/; ($y, $z) = ($1, $2); -print 'not ' unless (near($y, 1.5707963267949) and - near($z, -1.31695789692482)); -print "ok 5\n"; +ok(near($y, 1.5707963267949)); +ok(near($z, -1.31695789692482)); -print 'not ' unless (near(deg2rad(90), pi/2)); -print "ok 6\n"; +ok(near(deg2rad(90), pi/2)); -print 'not ' unless (near(rad2deg(pi), 180)); -print "ok 7\n"; +ok(near(rad2deg(pi), 180)); 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"; + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($z, 1)); ($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"; + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 1)); ($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"; + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($z, 0)); ($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"; + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 0)); } { 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"; + ok(near($r, sqrt(3))); + ok(near($t, deg2rad(45))); + ok(near($f, atan2(sqrt(2), 1))); ($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"; - + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 1)); + ($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"; + ok(near($r, sqrt(2))); + ok(near($t, deg2rad(45))); + ok(near($f, deg2rad(90))); ($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"; + ok(near($x, 1)); + ok(near($y, 1)); + ok(near($z, 0)); } { 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"; + ok(near($r, 1)); + ok(near($t, 1)); + ok(near($z, 1)); ($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"; + ok(near($r, 1)); + ok(near($t, 1)); + ok(near($z, 1)); } { use Math::Trig 'great_circle_distance'; - print 'not ' - unless (near(great_circle_distance(0, 0, 0, pi/2), pi/2)); - print "ok 18\n"; + ok(near(great_circle_distance(0, 0, 0, pi/2), pi/2)); - print 'not ' - unless (near(great_circle_distance(0, 0, pi, pi), pi)); - print "ok 19\n"; + ok(near(great_circle_distance(0, 0, pi, pi), pi)); # London to Tokyo. my @L = (deg2rad(-0.5), deg2rad(90 - 51.3)); @@ -159,8 +150,7 @@ use Math::Trig ':radial'; my $km = great_circle_distance(@L, @T, 6378); - print 'not ' unless (near($km, 9605.26637021388)); - print "ok 20\n"; + ok(near($km, 9605.26637021388)); } { @@ -169,62 +159,44 @@ use Math::Trig ':radial'; sub frac { $_[0] - int($_[0]) } my $lotta_radians = deg2rad(1E+20, 1); - print "not " unless near($lotta_radians, 1E+20/$R2D); - print "ok 21\n"; + ok(near($lotta_radians, 1E+20/$R2D)); my $negat_degrees = rad2deg(-1E20, 1); - print "not " unless near($negat_degrees, -1E+20*$R2D); - print "ok 22\n"; + ok(near($negat_degrees, -1E+20*$R2D)); my $posit_degrees = rad2deg(-10000, 1); - print "not " unless near($posit_degrees, -10000*$R2D); - print "ok 23\n"; + ok(near($posit_degrees, -10000*$R2D)); } { use Math::Trig 'great_circle_direction'; - print 'not ' - unless (near(great_circle_direction(0, 0, 0, pi/2), pi)); - print "ok 24\n"; + ok(near(great_circle_direction(0, 0, 0, pi/2), pi)); # 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"; +# ok(near(great_circle_direction(0, 0, pi, pi), -pi()/2)); my @London = (deg2rad( -0.167), deg2rad(90 - 51.3)); my @Tokyo = (deg2rad( 139.5), deg2rad(90 - 35.7)); my @Berlin = (deg2rad ( 13.417), deg2rad(90 - 52.533)); my @Paris = (deg2rad ( 2.333), deg2rad(90 - 48.867)); - print 'not ' - unless (near(rad2deg(great_circle_direction(@London, @Tokyo)), - 31.791945393073)); - print "ok 26\n"; + ok(near(rad2deg(great_circle_direction(@London, @Tokyo)), + 31.791945393073)); - print 'not ' - unless (near(rad2deg(great_circle_direction(@Tokyo, @London)), - 336.069766430326)); - print "ok 27\n"; + ok(near(rad2deg(great_circle_direction(@Tokyo, @London)), + 336.069766430326)); - print 'not ' - unless (near(rad2deg(great_circle_direction(@Berlin, @Paris)), - 246.800348034667)); + ok(near(rad2deg(great_circle_direction(@Berlin, @Paris)), + 246.800348034667)); - print "ok 28\n"; - - print 'not ' - unless (near(rad2deg(great_circle_direction(@Paris, @Berlin)), - 58.2079877553156)); - print "ok 29\n"; + ok(near(rad2deg(great_circle_direction(@Paris, @Berlin)), + 58.2079877553156)); use Math::Trig 'great_circle_bearing'; - print 'not ' - unless (near(rad2deg(great_circle_bearing(@Paris, @Berlin)), - 58.2079877553156)); - print "ok 30\n"; + ok(near(rad2deg(great_circle_bearing(@Paris, @Berlin)), + 58.2079877553156)); use Math::Trig 'great_circle_waypoint'; use Math::Trig 'great_circle_midpoint'; @@ -233,50 +205,39 @@ use Math::Trig ':radial'; ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.0); - print 'not ' unless (near($lon, $London[0])); - print "ok 31\n"; + ok(near($lon, $London[0])); - print 'not ' unless (near($lat, $pip2 - $London[1])); - print "ok 32\n"; + ok(near($lat, $pip2 - $London[1])); ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 1.0); - print 'not ' unless (near($lon, $Tokyo[0])); - print "ok 33\n"; + ok(near($lon, $Tokyo[0])); - print 'not ' unless (near($lat, $pip2 - $Tokyo[1])); - print "ok 34\n"; + ok(near($lat, $pip2 - $Tokyo[1])); ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.5); - print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E - print "ok 35\n"; + ok(near($lon, 1.55609593577679)); # 89.1577 E - print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N - print "ok 36\n"; + ok(near($lat, 1.20296099733328)); # 68.9246 N ($lon, $lat) = great_circle_midpoint(@London, @Tokyo); - print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E - print "ok 37\n"; + ok(near($lon, 1.55609593577679)); # 89.1577 E - print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N - print "ok 38\n"; + ok(near($lat, 1.20296099733328)); # 68.9246 N ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.25); - print 'not ' unless (near($lon, 0.516073562850837)); # 29.5688 E - print "ok 39\n"; + ok(near($lon, 0.516073562850837)); # 29.5688 E + + ok(near($lat, 1.170565013391510)); # 67.0684 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"; + ok(near($lon, 2.17494903805952)); # 124.6154 E - print 'not ' unless (near($lat, 0.952987032741305)); # 54.6021 N - print "ok 42\n"; + ok(near($lat, 0.952987032741305)); # 54.6021 N use Math::Trig 'great_circle_destination'; @@ -285,35 +246,28 @@ use Math::Trig ':radial'; ($lon, $lat) = great_circle_destination(@London, $dir1, $dst1); - print 'not ' unless (near($lon, $Tokyo[0])); - print "ok 43\n"; + ok(near($lon, $Tokyo[0])); - print 'not ' unless (near($lat, $pip2 - $Tokyo[1])); - print "ok 44\n"; + ok(near($lat, $pip2 - $Tokyo[1])); 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"; + ok(near($lon, $London[0])); - print 'not ' unless (near($lat, $pip2 - $London[1])); - print "ok 46\n"; + ok(near($lat, $pip2 - $London[1])); my $dir3 = (great_circle_destination(@London, $dir1, $dst1))[2]; - print 'not ' unless (near($dir3, 2.69379263839118)); # about 154.343 deg - print "ok 47\n"; + ok(near($dir3, 2.69379263839118)); # about 154.343 deg my $dir4 = (great_circle_destination(@Tokyo, $dir2, $dst2))[2]; - print 'not ' unless (near($dir4, 3.6993902625701)); # about 211.959 deg - print "ok 48\n"; + ok(near($dir4, 3.6993902625701)); # about 211.959 deg - print 'not ' unless (near($dst1, $dst2)); - print "ok 49\n"; + ok(near($dst1, $dst2)); } # eof