sphericalAngle

SPHERICALANGLE Compute angle between points on the sphere.

   ALPHA = sphericalAngle(P1, P2, P3)
   Computes angle (P1, P2, P3), i.e. the angle, measured at point P2,
   between the direction (P2, P1) and the direction (P2, P3).
   The result is given in radians, between 0 and 2*PI.

   Points are given either as [x y z] (there will be normalized to lie on
   the unit sphere), or as [phi theta], with phi being the longitude in [0
   2*PI] and theta being the elevation on horizontal [-pi/2 pi/2].


   NOTE: 
   this is an 'oriented' version of the angle computation, that is, the
   result of sphericalAngle(P1, P2, P3) equals
   2*pi-sphericalAngle(P3,P2,P1). To have the more classical relation
   (with results given betwen 0 and PI), it suffices to take the minimum
   of angle and 2*pi-angle.
   
   Examples
     % Use inputs as cartesian coordinates   
     p1 = [0 1 0];
     p2 = [1 0 0];
     p3 = [0 0 1];
     alpha = sphericalAngle(p1, p2, p3)
     alpha =
         1.5708

     % Use inputs as spherical coordinates   
     sph1 = [.1 0];
     sph2 = [0 0];
     sph3 = [0 .1];
     alphas = sphericalAngle(sph1, sph2, sph3)
     alphas =
         1.5708
 

   See also:
   geom3d, angles3d, spheres, sph2cart

0001 function alpha = sphericalAngle(p1, p2, p3)
0002 %SPHERICALANGLE Compute angle between points on the sphere.
0003 %
0004 %   ALPHA = sphericalAngle(P1, P2, P3)
0005 %   Computes angle (P1, P2, P3), i.e. the angle, measured at point P2,
0006 %   between the direction (P2, P1) and the direction (P2, P3).
0007 %   The result is given in radians, between 0 and 2*PI.
0008 %
0009 %   Points are given either as [x y z] (there will be normalized to lie on
0010 %   the unit sphere), or as [phi theta], with phi being the longitude in [0
0011 %   2*PI] and theta being the elevation on horizontal [-pi/2 pi/2].
0012 %
0013 %
0014 %   NOTE:
0015 %   this is an 'oriented' version of the angle computation, that is, the
0016 %   result of sphericalAngle(P1, P2, P3) equals
0017 %   2*pi-sphericalAngle(P3,P2,P1). To have the more classical relation
0018 %   (with results given betwen 0 and PI), it suffices to take the minimum
0019 %   of angle and 2*pi-angle.
0020 %
0021 %   Examples
0022 %     % Use inputs as cartesian coordinates
0023 %     p1 = [0 1 0];
0024 %     p2 = [1 0 0];
0025 %     p3 = [0 0 1];
0026 %     alpha = sphericalAngle(p1, p2, p3)
0027 %     alpha =
0028 %         1.5708
0029 %
0030 %     % Use inputs as spherical coordinates
0031 %     sph1 = [.1 0];
0032 %     sph2 = [0 0];
0033 %     sph3 = [0 .1];
0034 %     alphas = sphericalAngle(sph1, sph2, sph3)
0035 %     alphas =
0036 %         1.5708
0037 %
0038 %
0039 %   See also:
0040 %   geom3d, angles3d, spheres, sph2cart
0041 
0042 %   ---------
0043 %   author : David Legland
0044 %   INRA - TPV URPOI - BIA IMASTE
0045 %   created the 21/02/2005.
0046 %
0047 
0048 %   HISTORY
0049 %   23-05-2006 fix bug for points with angle from center > pi/2
0050 %   05-06-2013 fix bug for points given as spherical coordinates, better
0051 %       support for multiple inputs
0052 
0053 % test if points are given as matlab spherical coordinates
0054 if size(p1, 2) == 2
0055     [x, y, z] = sph2cart(p1(:,1), p1(:,2), ones(size(p1,1), 1));
0056     p1 = [x y z];
0057     [x, y, z] = sph2cart(p2(:,1), p2(:,2), ones(size(p2,1), 1));
0058     p2 = [x y z];
0059     [x, y, z] = sph2cart(p3(:,1), p3(:,2), ones(size(p3,1), 1));
0060     p3 = [x y z];
0061 end
0062 
0063 % normalize points
0064 p1  = normalizeVector3d(p1);
0065 p2  = normalizeVector3d(p2);
0066 p3  = normalizeVector3d(p3);
0067 
0068 % create the plane tangent to the unit sphere and containing central point
0069 plane = createPlane(p2, p2);
0070 
0071 % project the two other points on the plane
0072 pp1 = planePosition(projPointOnPlane(p1, plane), plane);
0073 pp3 = planePosition(projPointOnPlane(p3, plane), plane);
0074 
0075 % compute angle on the tangent plane
0076 pp2 = zeros(max(size(pp1, 1), size(pp3,1)), 2);
0077 alpha = angle3Points(pp1, pp2, pp3);
0078