;+
; NAME:
;	CvSky_Heliographic
; PURPOSE:
;	Euler angles for Ecliptic to Heliographic conversion
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION CvSky_Heliographic, UT, degrees=degrees, longitude=longitude
; INPUTS:
;	UT		array[n]; type: time structure
; OPTIONAL INPUT PARAMETERS:
;	longitude	array[n] or scalar; type: float; default: 0
;				heliographic longitude(s). See PROCEDURE.
;	/degrees	if set all angles are in degrees (default: radians)
; OUTPUTS:
;	abc		array[3,n]; type: float
;			Euler angles for rotation from ecliptic to heliographic
;			coordinate frame
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToDegrees, TimeGet, TimeUnit, SuperArray, SubArray
; SEE ALSO:
;	CvSky_IHG
; PROCEDURE:
;	The ecliptic longitude for the ascending node and the inclination of
;	equator to ecliptic are from R. Green, "Spherical astronomy", p. 431
;
;	If keyword 'longitude' is specified than the Euler angles refer
;	to a rotation from an ecliptic reference frame to a
;	heliographic reference frame with longitude measured from the
;	heliographic meridian through 'longitude', instead of the
;	heliographic prime meridian (i.e. heliographic longitude zero).
;	This is useful, for instance, to set up 'earth-centered'
;	heliographic coordinates with heliographic longitude measured
;	from the heliographic longitude of Earth.
; MODIFICATION HISTORY:
;	AUG-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

dpm = ToDegrees(degrees=degrees)

A = TimeGet(UT, /mjd, /scalar)+1780L		; MJD -1780=JD2398220-JD2400000
B = 0.25d0+(A-!sun.siderealp*long(A/!sun.siderealp))/!sun.siderealp

; A[0,*] = -90 deg + ecliptic longitude of ascending node of solar equator on ecliptic
; A[1,*] = inclination of solar equator to ecliptic
; A[2,*] = +90 deg + angle from ascending node to heliographic prime meridian
;			(prime meridian = heliographic longitude zero)

Y = TimeGet(UT,TimeUnit(/year),/scalar)-1850+(TimeGet(UT,/doy,/scalar)-1.0d0)/365.25d0

A = SuperArray(double(0*A),3,/lead)

A = SubArray(A, element=0, add=-16.333333d0+0.5025d0*Y/36.0d0)
A = SubArray(A, element=1, add=  7.25d0)

L = 360.00d0*(B-long(B))		; Angle asc node to heliographic prime meridian
IF IsType(longitude,/defined) THEN BEGIN
	L += longitude*dpm			; Add input heliographic longitude (in degrees)
	L = AngleRange(L,/degrees)	; Keep inside [0,360]
ENDIF
A = SubArray(A, element=2, add=L)

RETURN, A/dpm  &  END