;+
; NAME:
;	eclipsed_area
; PURPOSE:
;	Calculates fraction of solar disk eclipsed by Moon or Earth
; CATEGORY:
;	gen/idl/toolbox
; CALLING SEQUENCE:
	FUNCTION eclipsed_area, loc_sun, loc_em, $
		degrees		= degrees	, $
		elo			= elo		, $
		ut			= ut		, $
		plotc		= plotc		, $
		earth		= earth		, $
		radius		= radius	, $
		silent		= silent
; INPUTS:
;	loc_sun 		array[3,n] or array[3,2,n]; type: float
;						if array[3,n]; location of Sun in spherical coordinates
;							(RA,dec,distance, or long, lat, distance)
;						if array[3,2,n] then array[3,0,*] specifies the location
;							of the Sun and array[3,1,*] specifies the location
;							of the Moon or Earth. In this case loc_em is ignored.
;	loc_em			array[3,n]; type: float
;						location of Earth or Moon in same coordinates as for the Sun.
; OPTIONAL INPUT PARAMETERS:
;	/degrees		if set, all angles are in degrees (default: radians)
;	/plotc			if set, then results are plotted
;						(angular separation of Sun and Moon or Earth, and fraction of
;						disk eclipsed).
;	ut=ut			array[*]; type: structure
;						time structure used to plot a time axis (only used
;						if /plotc is set).
; OUTPUTS:
;	area			array[n]; type: float
;						fraction of solar disk eclipsed by Moon or Earth.
; OPTIONAL OUTPUT PARAMETERS:
;	elo 			array[*]; type: float
;						angular distance between Sun and Moon or Earth
; INCLUDE:
	@compile_opt.pro
; CALLS:
;	InitVar, ToRadians, ToDegrees, IsType, sphere_distance, BadValue
;	SyncArgs, IsTime, PlotCurve
; PROCEDURE:
; MODIFICATION HISTORY:
;	OCT-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, earth, /key
InitVar, silent, 0
IF silent LE 0 THEN message, /info, 'eclipse of Sun by '+(['Moon','Earth'])[earth]

IF IsType(loc_sun, /undefined) THEN BEGIN
	message, /info, 'no sun locations specified'
	RETURN, -1
ENDIF

InitVar, plotc, /key

rpm = ToRadians(degrees=degrees)
dpm = ToDegrees(degrees=degrees)

CASE n_params() EQ 2 OF
0: BEGIN
	rsun = reform(loc_sun[*,0,*])
	r_em = reform(loc_sun[*,1,*])
END
1: BEGIN
	rsun = reform(loc_sun)
	r_em = reform(loc_em )
END
ENDCASE

elo = sphere_distance(rsun, r_em, degrees=degrees)*rpm

; Angular size (radius) of Sun and Moon or Earth.

IF IsType(radius, /undefined) THEN radius = !earth.R*earth+!earth.RMoon*(1-earth)

rsun = asin( !sun.RAU/reform(rsun[2,*]) )
r_em = asin( radius/(reform(r_em[2,*])*!sun.au*1e8) )

A = BadValue(0.0)
SyncArgs, elo, A

; Separation larger than radius sun+radius moon or earth: no eclipse

tmp = where(elo GE rsun+r_em)
IF tmp[0] NE -1 THEN A[tmp] = 0.0

; Separation less than absolute difference in radii:
; total eclips (A > 1) or annular eclipse.

tmp = where(elo LE abs(r_em-rsun))
IF tmp[0] NE -1 THEN A[tmp] = (r_em[tmp]/rsun[tmp])^2

tmp = where(abs(rsun-r_em) LT elo and elo LT rsun+r_em)
IF tmp[0] NE -1 THEN BEGIN

	; Partial eclipse: two partially overlapping circles with 
	; radii rsun and r_em. Calculate area common to both circles.

	rsun2 = rsun[tmp]
	r_em2 = r_em[tmp]
	elo2  = elo [tmp]

	sunelo = rsun2*elo2
	em_elo = r_em2*elo2

	rsun2 = rsun2*rsun2
	r_em2 = r_em2*r_em2
	elo2  = elo2 *elo2

	csun = (rsun2+elo2-r_em2)/(2.0*sunelo)
	c_em = (r_em2+elo2-rsun2)/(2.0*em_elo) 

	; psun : angle between sun-moon line and sun  to points of intersection
	; pmoon: angle between sun-moon line and moon to points of intersection

	psun = acos( (csun < 1) > (-1) )
	p_em = acos( (c_em < 1) > (-1) )

	A[tmp] = (rsun2*(psun-csun*sin(psun))+r_em2*(p_em-c_em*sin(p_em)))/(!pi*rsun2)
	IF earth THEN A = A < 1

ENDIF

elo = elo/rpm

IF plotc AND n_elements(ut) GT 0 THEN BEGIN

	CASE IsTime(ut) OF
	0: tt = indgen(elo)
	1: tt = ut
	ENDCASE

	PlotCurve, tt, elo*dpm*60, xmargin=[8,8], ytitle='Angular distance (arcmin)'
	PlotCurve, tt, A, /oplot, /newyaxis, yaxis=0, ytitle='Eclipsed fraction of disk'

ENDIF

RETURN, A  &  END