; Diagnostic utility for evaluating spin model performance.
;
; Input: one or two tplot variable names
; Output: A new tplot variable (specified by output keyword)
; Units: Default to radians, convert to degrees with /degrees keyword.
; 
; Input data is assumed to be in DSL coordinates.  If
; two variables are passed, they are assumed to have the
; identical sample count and time tags. No checks are 
; performed...use with caution!
;
; For the one-argument case: output will be the angle between
; the DSL-X axis and each sample vector (after projecting onto
; the spin plane). The angles will be in the interval [0,360) 
; degrees (or equivalent in radians)
;
; For the two-argument case: both arguments are projected onto
; the spin plane, and the output is the angle between the two
; projected vectors.  The result will be in the interval [0, 180]
; degrees (or equivalent in radians).


pro spinplane_angle,v1, v2, output=output, degrees=degrees

; Unit conversion factor

if keyword_set(degrees) then begin
  convert=180.0D/!DPI
endif else begin
  convert=1.0D
endelse

; Get data for first variable
get_data,v1,data=d1
; Project it into the spin plane
proj1=d1.y[*,0:1]

if n_elements(v2) GT 0 then begin
   ; If a second positional argument is present, calculate the
   ; angle between the two sets of vectors after projecting them
   ; onto the spin plane.
   get_data,v2,data=d2
   len1=sqrt(total(proj1*proj1,2))  ; vector lengths for first arg
   proj2=d2.y[*,0:1]
   len2=sqrt(total(proj2*proj2,2))  ; vector lengths for second arg
   dot=total(proj1*proj2,2)         ; 2-d dot product
   costheta=dot/(len1*len2)         ; obtain cos(angle) from dot product
   a=where(costheta LT -1.0D,count)    ; clamp costheta range to [-1, 1]
   if (count GT 0) then costheta[a]=-1.0D
   a=where(costheta GT 1.0D,count)
   if (count GT 0) then costheta[a]=1.0D
   angle=acos(costheta)*convert     ; unit conversion
endif else begin
   ; If only one argument is given, project the vectors onto the spin
   ; plane and return the angles measured from the X axis (in the range
   ; [0, 360) degrees).
   angle=atan(proj1[*,1],proj1[*,0])*convert
endelse

; Shift output range to [0, 360) degrees (or radians equivalent)

a=where(angle LT 0.0D,count)
if (count GT 0 ) then begin
    angle[a] = angle[a] + convert*2.0D*!DPI
endif

; Store result in new tplot variable

store_data,output,data={x:d1.x, y:angle}
end