calculating the difference between 2 angles using st_azimuth and dot product

Question

the function below tries to calculate the difference between 2 angles. I use the mathematical reasoning is as follows :

1) Obtain the 2 angles :

angle1 = ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1))
angle2 = ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2))

2) Obtain vectors (of lenght == 1) from 2 angles :

V1 = (cos(angle1), sin(angle1))
V2 = (cos(angle2), sin(angle2))

3) Compute the dot product :

V1 * V2 = V1[x] * V2[x] + V1[y] * V2[y]
=
cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1)

Since the lenght of the vectors are both 1, it follows that :

the angle between the two vectors = ArcCos(dotProduct)/(length(v1)*legnth(v2)

= ArcCos(dotProduct)

The problem in the function below is that the dot product is yielding values greater than 1, which should not be possible if the math reasoning is correct, and that is my question :

is the reasoning correct ?

CREATE OR REPLACE FUNCTION angleDiff(l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
BEGIN           
  SELECT ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1)) into angle1;
  SELECT ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2)) into angle2;  

  RETURN acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1));
END;
$$  LANGUAGE plpgsql

Answer 1

This merges your logic and whuber's logic, except the return is signed:

CREATE OR REPLACE FUNCTION angleDiff (l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
DECLARE diff FLOAT;
BEGIN           
    SELECT ST_Azimuth (ST_StartPoint(l1), ST_EndPoint(l1)) INTO angle1;
    SELECT ST_Azimuth (ST_StartPoint(l2), ST_EndPoint(l2)) INTO angle2;  
    SELECT degrees (angle2 - angle1) INTO diff;
    CASE
    WHEN diff >   180 THEN RETURN diff - 360;
    WHEN diff <= -180 THEN RETURN diff + 360;
    ELSE                   RETURN diff;
    END CASE;
END;
$$  LANGUAGE plpgsql

The return value is clockwise-positive angle [0 to 180] from line1 to line2. If line2 is "clockwise before" line1, the return value is negative (0 to -180). Ignore the sign if you don't care about order.

Answer 2

Ok, there were 2 problems, but the math reasoning is correct :

1) the dot product was

acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle1));
                                                   ^^^^^^^

it should be :

-- should be : 
acos(cos(angle1) * cos(angle2) + sin(angle1) * sin(angle2));
                                                   ^^^^^^^

and also, the dot product can go slightly outside the range of [-1, 1], so it needs to be checked.

The corrected versions :

CREATE OR REPLACE FUNCTION angleDiff(l1 geometry, l2 geometry)
RETURNS FLOAT AS $$
DECLARE angle1 FLOAT;
DECLARE angle2 FLOAT;
DECLARE dotProduct FLOAT;
BEGIN           
  SELECT ST_Azimuth(ST_StartPoint(l1), ST_EndPoint(l1)) into angle1;
  SELECT ST_Azimuth(ST_StartPoint(l2), ST_EndPoint(l2)) into angle2;  

  select (cos(angle1) * cos(angle2) + sin(angle1) * sin(angle2)) into dotProduct;  

  if dotProduct > 1 then return 0;
  elseif dotProduct < -1 then return 180;
  else return degrees(acos(dotProduct));
  end if;

END;
$$  LANGUAGE plpgsql