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 % E-mail: david.legland@inrae.fr
0045 % Created: 2005-02-21
0046 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0047 
0048 % test if points are given as matlab spherical coordinates
0049 if size(p1, 2) == 2
0050     [x, y, z] = sph2cart(p1(:,1), p1(:,2), ones(size(p1,1), 1));
0051     p1 = [x y z];
0052     [x, y, z] = sph2cart(p2(:,1), p2(:,2), ones(size(p2,1), 1));
0053     p2 = [x y z];
0054     [x, y, z] = sph2cart(p3(:,1), p3(:,2), ones(size(p3,1), 1));
0055     p3 = [x y z];
0056 end
0057 
0058 % normalize points
0059 p1  = normalizeVector3d(p1);
0060 p2  = normalizeVector3d(p2);
0061 p3  = normalizeVector3d(p3);
0062 
0063 % create the plane tangent to the unit sphere and containing central point
0064 plane = createPlane(p2, p2);
0065 
0066 % project the two other points on the plane
0067 pp1 = planePosition(projPointOnPlane(p1, plane), plane);
0068 pp3 = planePosition(projPointOnPlane(p3, plane), plane);
0069 
0070 % compute angle on the tangent plane
0071 pp2 = zeros(max(size(pp1, 1), size(pp3,1)), 2);
0072 alpha = angle3Points(pp1, pp2, pp3);
0073