;****f* Coordinates/vector_angle
;
; NAME
;    vector_angle - Find the angle between two cartesian vectors
;
; SYNOPSIS
;    angle = vector_angle(X, Y, A)
;
; INPUTS
;    X             Vector 1
;    Y             Vector 2
;    A             Orientation axis (optional)
;
; OUTPUT
;    Angle       Angle between vectors
;
; AUTHOR
;   Met Office, Exeter, UK.
;   Any comments on this software should be given via the ROM SAF
;   Helpdesk at http://www.romsaf.org
;
; COPYRIGHT
;   (c) EUMETSAT. All rights reserved.
;   For further details please refer to the file COPYRIGHT
;   which you should have received as part of this distribution.
;
;****
  
pro SUB_vector_angle, X, Y, A=A, angle
  
;   vector_product = (/ X(2)*Y(3) - X(3)*Y(2),  &
;                       X(3)*Y(1) - X(1)*Y(3),  &
;                       X(1)*Y(2) - X(2)*Y(1) /)

  if keyword_set(A) then begin
     n = A
  endif else begin
     n = CROSSP(X, Y)
  endelse
  
  nn = total(n * n) ; dot_product

  if nn eq 0 then begin
     angle = 0.0
  endif else begin
     n /= sqrt(nn)

     X_0 = CROSSP(n, X)
     X_1 = X - total(n * X) * n
     Y_1 = Y - total(n * Y) * n

     angle = 180.D/!pi*atan(total(X_0*Y_1), total(X_1*Y_1))
  endelse
      
end