;+
; NAME:
;	jpl_phase
; PURPOSE:
;	Calculate phase angle for body in JPL ephemeris
; CATEGORY:
;	gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION jpl_phase, ut, body	, $
		observer	= observer	, $
		degrees 	= degrees	, $
		sun 		= sun		, $
		pos 		= pos
; INPUTS:
;	ut				array[ntime]; type: time structure
;						list of times for which phase angles are needed
;	body			scalar or array[nbody]; type: integer; default: jpl_body(/moon)
;						body for which phase is needed
; OPTIONAL INPUT PARAMETERS:
;	observer=observer
;					scalar or array[nobs]; type: integer; default: jpl_body(/earth)
;						body to used as observer location
;	/degrees		if set the phase angle is returned in degrees
;						(default: radians)
; OUTPUTS:
;	Result			scalar or array[nbody,nobs,ntime]; type: float
;						(dummy dimensions are omitted)
;						phase angles
; OPTIONAL OUTPUTS:
;	sun 			array[3,nbody,nobs,ntime]; type: float
;						(dummy dimensions are omitted)
;						rectangular ecliptic location of Sun in topocentric
;						coordinates of 'body' (AU)
;	pos 			array[3,nbody,nobs,ntime]; type: float
;						(dummy dimensions are omitted)
;						rectangular ecliptic location of observer in topocentric
;						coordinates of 'body' (AU)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, jpl_eph, jpl_body, SuperArray, scalarproduct, ToRadians
;	SubArray
; PROCEDURE:
;	The phase angle is defined as the angle Sun-body-observer
;	The visible fraction of the disk of body as seen by observer
;	is approximately 1-phase/180, so phase=0 is 'full' and
;	phase=180 is 'new'.
;	'waxing' means that the size of the visible disk is increasing,
;	i.e. decreasing phase angle from 180 to 0; 'waning' means
;	the size of the disk is decreasing i.e. phase angle is
;	increasing from 0 to 180.
; MODIFICATION HISTORY:
;	NOV-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, body	 , jpl_body(/moon )
InitVar, observer, jpl_body(/earth)

nbody = n_elements(body)

IF nbody GT 1 THEN BEGIN

	phase = jpl_phase(ut,body[0],observer=observer,degrees=degrees, sun=sun, pos=pos)

	phase = SuperArray(phase,nbody,/lead)
	sun   = SuperArray(sun  ,nbody,after=1)
	pos   = SuperArray(pos  ,nbody,after=1)

	FOR i=0,n_elements(body)-1 DO BEGIN
		tmp_phase = jpl_phase(ut,body[i],observer=observer,degrees=degrees, sun=tmp_sun, pos=tmp_pos)
		phase = SubArray(phase,dim=0,element=i,replace=tmp_phase)
		sun   = SubArray(sun  ,dim=1,element=i,replace=tmp_sun  )
		pos   = SubArray(pos  ,dim=1,element=i,replace=tmp_pos  )
	ENDFOR

	RETURN, phase

ENDIF

; Position of Sun and observer in rectangular topocentric
; coordinates of 'body'.

sun = jpl_eph(ut, jpl_body(/sun), center=body, /to_ecliptic)
pos = jpl_eph(ut, observer		, center=body, /to_ecliptic)

sz6 = size(pos,/dim)
sz3 = [3,2]
IF size(pos,/n_dim) GT 1 THEN sz3 = [sz3,sz6[1:*]]

nb = n_elements(observer)
IF nb NE 1 THEN sun = SuperArray(sun,nb,after=1)
sun = reform(sun,sz6,/overwrite)

pos = reform(pos,sz3,/overwrite)
sun = reform(sun,sz3,/overwrite)
									; Velocity of observer relative to Sun
vob = SubArray(pos,dim=1,element=1)-SubArray(sun,dim=1,element=1)

pos = SubArray(pos,dim=1,element=0) ; Body-to-observer vector
sun = SubArray(sun,dim=1,element=0) ; Body-to-sun vector

; Calculate phase angle from scalar product of sun and pos
; (0 <= phase angle <= 180).

phase = acos(scalarproduct(pos,sun,/angle))/ToRadians(degrees=degrees)

; Negative phase angle means waning disk (increasing phase angle)
; Positive phase angle means waxing disk (decreasing phase angle)

phase *= -(2*(scalarproduct(vob,pos) GE 0)-1)*(2*(scalarproduct(sun,pos) GE 0)-1)

RETURN, phase  &  END