function directedAngle,vect1,vect2, pole, degrees=degrees
;directed angular distance between vectors

;compute the vector dot product for both points
if size(vect1,/n_dim) eq 1 then begin
  vectorDotProduct =  dotp(vect1,vect2)/(norm(vect1)*norm(vect2))
  ;make sure the dot product is between +/- 1
  vectorDotProduct = -1.0D > vectorDotProduct < 1.0D
  angle = acos(vectorDotProduct)
  cp = crossp(vect1,vect2)
  ;dot product
  dp = dotp(cp,pole)
  if sign(dp) eq -1 then angle = (!pi *2.) - angle
endif else begin
  vectorDotProduct = total(vect1*vect2,2)/(vector_length(vect1,2) * vector_length(vect2,2))
  ;make sure the dot product is between +/- 1
  vectorDotProduct = -1.0D > vectorDotProduct < 1.0D
  angle = acos(vectorDotProduct)
  ;compute the vector cross product
  cp = [[vect1[*,1]*vect2[*,2]-vect2[*,1]*vect1[*,2]], $
        [vect1[*,2]*vect2[*,0]-vect2[*,2]*vect1[*,0]], $
        [vect1[*,0]*vect2[*,1]-vect2[*,0]*vect1[*,1]]]
  dp = acos(total(cp*pole,2)/(vector_length(cp,2) * vector_length(pole,2)))
  ind = where(sign(dp) eq -1, count)
  if count gt 0 then angle[ind] = (!pi *2.) - angle[ind]
endelse


if keyword_set(degrees) then angle=angle/!dtor
return,angle
end