; docformat = 'rst'
;+
; :Description:
; Given the angle between two directions (`dphase`) with angle from vertical
; given by `dinc`, `dema`, returns the azimuth difference between them. All inputs
; and the output are in degrees. This function is vectorized, so the 3 input arguments
; can be of any dimension (must be the same for all 3), and that will be the dimension of the
; result.
; This function is complementary to `pp_phase_angle`.
;
;
; :Params:
; dphase : in, required
; The angle between the two directions (phase angle), in degrees.
; dinc : in, required
; The angle between one of the directions and the vertical (incidence angle), in degrees.
; dema : in, required
; The angle between the other direction and the vertical (emission angle), in degrees.
;
; :Examples:
; Taking incidence and emission in the same plane::
;
; print,pp_azdif(30d0,40d0,70d0)
; ;0.0000000
; print,pp_azdif(110d0,40d0,70d0)
; ;180.00001
;
; Taking incidence and emission at the horizon::
;
; print,pp_azdif(110d0,90d0,90d0)
; ;110.00000
;
; :Uses:
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Sep/2009
;-
function pp_azdif,dphase,dinc,dema
compile_opt idl2
;Convert input to radians
ema=dema*!dpi/180d0
inc=dinc*!dpi/180d0
phase=dphase*!dpi/180d0
;Get the necessary sines and cosines
s1=sin(ema)
s2=sin(inc)
c1=cos(ema)
c2=cos(inc)
cp=cos(phase)
;Get the cosine of the azimuth difference
tmp=(cp-(c1*c2))/(s1*s2)
;Treat special cases (incidence or emission is at zenith or nadir)
w=where(s1*s2 eq 0d0,nw)
if (nw gt 0) then tmp[w]=1d0
;Use the complex arc cosine to avoid NaNs due to lack of precision
rdeltaphi=double(acos(complex(tmp,0d0)))
return,rdeltaphi*180d0/!dpi
end