From 7e5f197a77be82f8e343b9d7b685f1987eefd832 Mon Sep 17 00:00:00 2001
From: Jarkko Hietaniemi <jhi@iki.fi>
Date: Sun, 1 Apr 2001 18:43:09 +0000
Subject: [PATCH] Add great_circle_direction().

p4raw-id: //depot/perl@9504
---
 lib/Math/Trig.pm | 55 +++++++++++++++++++++++++++++++++++++++++++++++--------
 t/lib/trig.t     | 23 ++++++++++++++++++++++-
 2 files changed, 69 insertions(+), 9 deletions(-)

diff --git a/lib/Math/Trig.pm b/lib/Math/Trig.pm
index b28f150..2a23590 100644
--- 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
--- 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
-- 
2.7.4