C+
C NAME:
C	ECLIPTIC_HELIOGRAPHIC
C PURPOSE:
C	Converts ecliptic coordinates into heliographic coordinates or v.v.
C CATEGORY:
C	Math: coordinate transformation
C CALLING SEQUENCE:
	subroutine ECLIPTIC_HELIOGRAPHIC(ID,iYr,Doy,PHI,RLAT)
C INPUTS:
C	ID	integer	ID=0: ecliptic     ---> heliographic
C			ID=1: heliographic ---> ecliptic
C	iYr	integer	year of current date
C	Doy	real	day of year, including a fraction for the time of day
C INPUTS/OUTPUTS:
C	PHI	real	ecliptic longitude for equinox of current date/
C			heliographic longitude
C	RLAT	real	ecliptic latitude for equinox of current date/
C			heliographic latitude
C INCLUDE:
	include		'sun.h'
C CALLS:
C	Julian, rotate
C SIDE EFFECTS:
C	ALFA	-90 deg. plus ecliptic longitude of ascending node of solar 
C		equator on ecliptic
C	BETA	inclination of solar equator with respect to ecliptic
C	GAMMA	90 deg. plus angle from ascending node to heliographic prime 
C		meridian 
C PROCEDURE:
C 	See R. Green, Spherical Astronomy, Cambridge UP, 1985, Section 17.7,
C	p. 430.
C MODIFICATION HISTORY:
C	Paul Hick (MPAE,UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	ID
	    integer	iYr
	    real	Doy
	    real	PHI
	    real	RLAT

	double precision	JD
	double precision	JEpoch
	double precision	GAM
	double precision	PERIOD

	PERIOD = SUN__SIDEREALP
	if (ID .ne. 0 .and. ID .ne. 1) stop 'ECLIPTIC_HELIOGRAPHIC: invalid ID'

	ALFA	= -16.333333+0.5025*(iYr-1850+(Doy-1)/365.25)/36.0
	BETA	= SUN__I 				! Euler angles
	call Julian(0,iYr,Doy,JD,JEpoch)		! Get JD number
	JD	= JD-2398220.0d0
	GAM	= 0.25d0+(JD-PERIOD*int(JD/PERIOD))/PERIOD
	GAMMA	= 360.0*(GAM-int(GAM))

	if (ID .eq. 0) then
	    call rotate(ALFA,BETA,GAMMA,PHI,RLAT)	! To heliographic system
	else
	    call rotate(-GAMMA,-BETA,-ALFA,PHI,RLAT)	! To ecliptic system
	end if

	return
	end