C+
C NAME:
C	smei_eclipse
C PURPOSE:
C	Check whether SMEI is seeing a solar eclipse
C	(by Moon) or is in Earth shadow.
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	logical function smei_eclipse(tt,id,fraction)
C INPUTS:
C	tt(2)		integer		UT time
C	id		integer		10: solar eclipse by Moon
C					 3: solar eclipse by Earth
C					    (i.e. in Earth shadow)
C OUTPUTS:
C	smei_eclipse
C			logical		.TRUE.: eclipse in progress
C	fraction	real		fraction of solar disk covered.
C					For eclipses by Moon the fraction
C					will be larger than one for a total.
C					For Earth shadow the fraction is
C					never larger than one.
C CALLS:
C	Time2smei_eph, Time2jpl_eph, Say
C IMPLICIT:
	implicit double precision (A-H,O-Z)
C INCLUDE:
	include		'sun.h'
	include		'planet.h'
	include		'math.h'
C MODIFICATION HISTORY:
C	FEB-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer		tt(2)
	    integer		id
	    real		fraction

	character		cSay*12		/'smei_eclipse'/
	logical			Time2smei_eph
	logical         	Time2jpl_eph
	double precision	rr(6)
	double precision	rsun(3)
	double precision	r_em(3)

	if (.not. Time2smei_eph(tt,rr,.TRUE.)) call Say(cSay,'E','smei','ephemeris error')

	rx = rr(1)
	ry = rr(2)
	rz = rr(3)

	! id = 3 for Earth, 10 for Moon

	if (.not. Time2jpl_eph(tt,0,3,rr,.TRUE.,.FALSE.)) call Say(cSay,'E','jpl','ephemeris error')

	rsun(1) = rr(1)-rx			! Location of Sun relative to spacecraft
	rsun(2) = rr(2)-ry
	rsun(3) = rr(3)-rz

	asun = dsqrt(rsun(1)*rsun(1)+rsun(2)*rsun(2)+rsun(3)*rsun(3))

	if (id .eq. 10) then

	    if (.not. Time2jpl_eph(tt,10,3,rr,.TRUE.,.FALSE.)) call Say(cSay,'E','jpl','ephemeris error')

	    r_em(1) = rr(1)-rx			! Location of Moon relative to spacecraft
	    r_em(2) = rr(2)-ry
	    r_em(3) = rr(3)-rz

	    a_em = dsqrt(r_em(1)*r_em(1)+r_em(2)*r_em(2)+r_em(3)*r_em(3))

	    elo  = (rsun(1)*r_em(1)+rsun(2)*r_em(2)+rsun(3)*r_em(3))/(asun*a_em)
	    a_em = dasin(PLANET__RMOON/a_em)! Angular radius of Moon

	else
	    r_em(1) = -rx			! Location of Earth relative to spacecraft
	    r_em(2) = -ry
	    r_em(3) = -rz

	    a_em = dsqrt(r_em(1)*r_em(1)+r_em(2)*r_em(2)+r_em(3)*r_em(3))

	    elo = (rsun(1)*r_em(1)+rsun(2)*r_em(2)+rsun(3)*r_em(3))/(asun*a_em)
	    a_em = dasin(PLANET__REARTH_MEAN/a_em)! Angular radius of Earth

	end if

	asun = dasin(SUN__R*1d5/asun)	! Angular radius of Sun
	elo = dacos( max(-1d0,min(elo,1d0)) )	! Solar elongation angle Moon/Earth

	smei_eclipse = elo .le. asun+a_em

	if (elo .le. dabs(asun-a_em)) then	! Total or annular eclipse

	    fraction = a_em/asun
	    fraction = fraction*fraction

	else if (elo .ge. asun+a_em) then 	! No eclipse

	    fraction = 0.0

	else					! Partial eclipse

	    asun2 = asun
	    a_em2 = a_em
	    elo2  = elo

	    sunelo = asun2*elo2
	    em_elo = a_em2*elo2

	    asun2 = asun2*asun2
	    a_em2 = a_em2*a_em2
	    elo2  = elo2 *elo2

	    csun = (asun2+elo2-a_em2)/(2d0*sunelo)
	    c_em = (a_em2+elo2-asun2)/(2d0*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 = dacos( max(-1d0,min(csun,1d0)) )
	    p_em = dacos( max(-1d0,min(c_em,1d0)) )

	    ! Fraction of the solar disk covered

	    fraction = (asun2*(psun-csun*dsin(psun))+a_em2*(p_em-c_em*dsin(p_em)))/(MATH__PI*asun2)

	end if

	if (id .eq. 3) fraction = min(fraction,1.0)

	return
	end