; ****************************************************************
; ANGLES--computes wind axes rotation needed to get mean w and mean v 
; equal zero.  Rotates phi degrees around the z-axis and then theta degrees 
; around the y-axis.
; Linsey Marr, June 1995
; modified by Ken Davis, August 1995.
; further modified by Ken Davis June, 1996 - original code is ROTATED.PRO.
; ANGLES only computes the rotation angles and does not perform the rotation.
; ****************************************************************
pro angles, u, v, w, theta, phi, badval

; Computes rotation angles from long-term winds (day long?) so that 
;  vmean=0 and wmean=0.
; These rotation angles can then be applied to short-term (hour long?)
;   time series.

; Variables--
; g             array of indices of good data points
; count         number of good data points
; umean         mean wind, using good points only
; vmean
; wmean
; theta         angle of rotation about y-axis
; phi           angle of rotation about z-axis
; r             rotation matrix (note: IDL is column major)
; a             temporary array of product of wind vector and rotation matrix
; badval	marker for bad or missing data points

;get double precision pi
pi=2.*asin(double(1.))

; Get indices of good data points
g = where(u ne badval, gn)
h = where(v ne badval, hn)
k = where(w ne badval, kn)

; Find mean wind
if (gn gt 0) then umean = total(u(g),/double) / gn
if (hn gt 0) then vmean = total(v(h),/double) / hn
if (kn gt 0) then wmean = total(w(k),/double) / kn

if ((gn le 0) or (hn le 0) or (kn le 0)) then begin
  print,'Unable to compute rotation angles - no data for at least one wind dir'
  printf,8,'Unable to compute rotation angles - no data for at least one wind dir'
  theta=badval
  phi=badval
  return
endif else begin

; Find angles to rotate through

;theta = asin(wmean / sqrt(umean^2 + vmean^2 + wmean^2))
;phi = asin(vmean / sqrt(umean^2 + vmean^2))

;wwg, april 19, 2000

if ((umean^2 + vmean^2 + wmean^2) ne 0. and (umean^2 + vmean^2) ne 0.0) then begin
theta = asin(wmean / sqrt(umean^2 + vmean^2 + wmean^2))
phi = asin(vmean / sqrt(umean^2 + vmean^2))
endif else begin
theta=badval
phi=badval
return
endelse   


; Put horizontal rotation into the correction quadrant.
if ((umean lt 0.) and (vmean lt 0.)) then phi=-pi-phi
if ((umean lt 0.) and (vmean gt 0.)) then phi=pi-phi

endelse

return
end