Add great_circle_direction().

p4raw-id: //depot/perl@9504

diff --git a/lib/Math/Trig.pm b/lib/Math/Trig.pm
index b28f150..2a23590 100644 (file)
--- a/lib/Math/Trig.pm
+++ b/lib/Math/Trig.pm
@@ -32,7 +32,7 @@ my @rdlcnv = qw(cartesian_to_cylindrical
                spherical_to_cartesian
                spherical_to_cylindrical);
 
-@EXPORT_OK = (@rdlcnv, 'great_circle_distance');
+@EXPORT_OK = (@rdlcnv, 'great_circle_distance', 'great_circle_direction');
 
 %EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
 
@@ -130,6 +130,20 @@ sub great_circle_distance {
              sin( $lat0 ) * sin( $lat1 ) );
 }
 
+sub great_circle_direction {
+    my ( $theta0, $phi0, $theta1, $phi1 ) = @_;
+
+    my $lat0 = pip2 - $phi0;
+    my $lat1 = pip2 - $phi1;
+
+    my $direction =
+       atan2(sin($theta0 - $theta1) * cos($lat1),
+             cos($lat0) * sin($lat1) -
+             sin($lat0) * cos($lat1) * cos($theta0 - $theta1));
+
+    return rad2rad($direction);
+}
+
 =pod
 
 =head1 NAME
@@ -383,12 +397,12 @@ Notice that when C<$z> is not 0 C<$rho_c> is not equal to C<$rho_s>.
 
 =back
 
-=head1 GREAT CIRCLE DISTANCES
+=head1 GREAT CIRCLE DISTANCES AND DIRECTIONS
 
 You can compute spherical distances, called B<great circle distances>,
-by importing the C<great_circle_distance> function:
+by importing the great_circle_distance() function:
 
-       use Math::Trig 'great_circle_distance'
+  use Math::Trig 'great_circle_distance';
 
   $distance = great_circle_distance($theta0, $phi0, $theta1, $phi1, [, $rho]);
 
@@ -409,10 +423,26 @@ 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:
+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);
 
@@ -422,8 +452,17 @@ To calculate the distance between London (51.3N 0.5W) and Tokyo (35.7N
 
         $km = great_circle_distance(@L, @T, 6378);
 
-The answer may be off by few percentages because of the irregular
-(slightly aspherical) form of the Earth.  The used formula
+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

diff --git a/t/lib/trig.t b/t/lib/trig.t
index 6949622..def7c15 100755 (executable)
--- a/t/lib/trig.t
+++ b/t/lib/trig.t
@@ -30,7 +30,7 @@ sub near ($$;$) {
     $_[1] ? (abs($_[0]/$_[1] - 1) < $e) : abs($_[0]) < $e;
 }
 
-print "1..23\n";
+print "1..26\n";
 
 $x = 0.9;
 print 'not ' unless (near(tan($x), sin($x) / cos($x)));
@@ -176,4 +176,25 @@ use Math::Trig ':radial';
     print "ok 23\n";
 }
 
+{
+    use Math::Trig 'great_circle_direction';
+
+    print 'not '
+       unless (near(great_circle_direction(0, 0, 0, pi/2), pi));
+    print "ok 24\n";
+
+    print 'not '
+       unless (near(great_circle_direction(0, 0, pi, pi), -pi/2));
+    print "ok 25\n";
+
+    # London to Tokyo.
+    my @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
+    my @T = (deg2rad(139.8),deg2rad(90 - 35.7));
+
+    my $rad = great_circle_direction(@L, @T);
+
+    print 'not ' unless (near($rad, -0.546644569997376));
+    print "ok 26\n";
+}
+
 # eof