;+
; NAME:
;	jpl_mag
; PURPOSE:
;	Calculate apparent magnitude for body in JPL ephemeris
; CATEGORY:
;	gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION jpl_mag, ut, body	, $
		observer	= observer	, $
		mag_body	= mag_body	, $
		degrees 	= degrees	, $
		phase		= phase 	, $
		int_phase	= int_phase
; 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 apparent magnitude is needed
; OPTIONAL INPUT PARAMETERS:
;	observer=observer
;					scalar or array[nobs]; type: integer
;						default: jpl_body(/earth)
;	mag_body		same structure as 'body'; type: float
;						default: ([-26.73,-0.42,-4.40,-3.86,
;									-1.52,-9.40,-8.88,-7.19,-6.87,-1.0,0.21])[body]
;						absolute magnitude of 'body'
;	/degrees		if set the phase angle is output in degrees
;						(default: radians)
;	int_phase=int_phase
;					scalar; type: string
;						name of function used for the phase integral
;						Must be function with as only input argument
;						the phase angle in radians.
; OUTPUTS:
;	Result			array[nbody,nobs,ntime]; type: float
;						(dummy dimensions are omitted)
;						apparent magnitude
; OPTIONAL OUTPUTS:
;	phase			array[nbody,nobs,ntime]; type: float
;						(dummy dimensions are omitted)
;						phase angle as returned by jpl_phase
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	jpl_phase, ToRadians, InitVar, IsType
; PROCEDURE:
;	Visual absolute magnitudes for all bodies in the JPL ephemeris
;	are hardcoded. Override with the 'mag_body' keyword.
;	The apperent magnitude calculation uses the phase integral
;	for an ideal diffuse reflecting sphere by default.
;	Override this using argument 'int_phase'
;
;	The difference between jpl_eph and JPL's HORIZONS system is usally about ~.2-3, whereas the difference between 
;	an astronomical almanac and jpl_eph is usually on the order of ~.4-5.
; MODIFICATION HISTORY:
;	NOV-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

phase = jpl_phase(ut,body,observer=observer,degrees=degrees,sun=sun,pos=pos)
angle = abs(phase)*ToRadians(degrees=degrees)

;radius = [0.0,2439.0,6052.0,!earth.R,3397.0,71492.0,60268.0,25559.0,24764.0,1195.0,!earth.RMoon]
;radius = radius[body]/(!sun.au*1e8)
;albedo = [0.0,0.106,0.65,0.367,0.150,0.52,0.47,0.51,0.41,0.6,0.12]
;albedo = albedo[body]
;magnitude = -26.73-5.0*alog10(sqrt(albedo)*radius)

InitVar, mag_body, ([-26.73,-0.42,-4.40,-3.86,-1.52,-9.40,-8.88,-7.19,-6.87,-1.0,0.21])[body]

CASE IsType(int_phase,/string) OF
0: phase_integral = ((1-angle/!pi)*cos(angle)+sin(angle)/!pi)/1.5
1: phase_integral = call_function(int_phase,angle)
ENDCASE

magnitude = mag_body
IF n_elements(phase) GT 1 THEN BEGIN
	magnitude = SuperArray(magnitude,n_elements(observer)*n_elements(ut))
	magnitude = reform(magnitude,size(phase,/dim))
ENDIF

magnitude = 2.5*alog10(total(sun*sun,1)*total(pos*pos,1)/phase_integral)+magnitude

RETURN, magnitude  &  END