function aSunView, rr, vv, alfa, degrees=Degrees
on_error, 2

if keyword_set(Degrees) then dpr = !radeg else dpr = 1.

bOne = n_elements(rr) eq 3

case bOne of
0: begin
	r = rr
	v = vv
	a = alfa/dpr
end
1: begin
	r = [[rr],[rr]]
	v = [[vv],[vv]]
	a = [[alfa],[alfa]]/dpr
end
endcase

r	= cv_coord(from_sphere=r, /to_rect, degrees=Degrees)
s	= [1.,1.,1.]#sqrt(total(r*r,1))
r	= r/s								; Unit vector along r

v	= cv_coord(from_sphere=v, /to_rect, degrees=Degrees)
s	= [1.,1.,1.]#sqrt(total(v*v,1))
v	= v/s								; Unit vector along v

rxv	= shift(r,2,0)*shift(v,1,0)-shift(r,1,0)*shift(v,2,0)	; Cross product r x v
s	= [1.,1.,1.]#sqrt(total(rxv*rxv,1))
rxv = rxv/s								; Unit vector along r x v
										; Cross product (r x v) x r
										; = Unit vector perpendicular to r and v
rxvxr = shift(rxv,2,0)*shift(r,1,0)-shift(rxv,1,0)*shift(r,2,0)

; s1 is a unit vectors in the r,v plane with angles alfa relative to r.
; r divides the r,v plane in two halves. s1 lies in the same half-plane as v
; for alfa in [0,180].

s = r*([1.,1.,1.]#cos(a)) + rxvxr*([1.,1.,1.]#sin(a))

; The cross product s x v is parallel to  r x v.
; If (s x v).(r x v) > 0 then s lies in between r and v
; If (s x v).(r x v) < 0 then v lies in between r and s

sxv = shift(s,2,0)*shift(v,1,0)-shift(s,1,0)*shift(v,2,0)	; Cross product s x v

sr = total(sxv*rxv,1)				; Scalar product (s x v).(r x v)
sr = where(sr lt 0)
if sr(0) ne -1 then s(*,sr) = v(*,sr)

if bOne then s = s(*,0)

s = cv_coord(from_rect=s, /to_sphere, degrees=Degrees)
s(0,*) = s(0,*)+2*!pi*dpr*(s(0,*) lt 0)

return, s  &  end