;+
; NAME:
;	CvSky_GSM
; PURPOSE:
;	Euler angles for Ecliptic to GSM conversion
; CATEGORY:
;	smei/gen/idl/ephem
; CALLING SEQUENCE:
	FUNCTION CvSky_GSM, UT, degrees=degrees
; INPUTS:
;	UT	array[n]; type: time structure
; OPTIONAL INPUT PARAMETERS:
;	/degrees    if set all angles are in degrees (default: radians)
; OUTPUTS:
;	Result	array[3,n]; type: float
;		Euler angles
; INCLUDE:
	@compile_opt.pro							; On error, return to caller
; CALLS:
;	ToRadians, ToDegrees, SubArray, SuperArray, TimeGet, CvSky, big_eph
; PROCEDURE:
;	Geocentric Solar Magnetospheric (GSM) coordinates
;	x-axis points from Earth to Sun
;	z-axis is the northern magnetic pole, projected in plane perp to x-axis.
;	The difference between the GSM system and the GSE
;	and GSEQ is a rotation about the X-axis.
; MODIFICATION HISTORY:
;	AUG-2002, Paul Hick (UCSD/CASS)
;	JUN-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bugfix. The conversion of UT to years since 2000 somehow
;		lost the subtraction of 2000.
;-

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

B = TimeGet(UT, /jepoch, /scalar)-2000

; Set up geographic coordinates of Earth's magnetic north pole

A = SuperArray(0*B,2,/lead)
A = SubArray(A, element=0, add = 288.44d0-0.04236d0*B)
A = SubArray(A, element=1, add =  79.53d0-0.03560d0*B)

A /= dpm

; Get ecliptic coordinates of Earth's magnetic north pole

A = CvSky(UT,from_geographic=A,/to_ecliptic,degrees=degrees,/silent)

L = SubArray(A,element=0)		; Ecliptic longitude magnetic pole
D = SubArray(A,element=1)		; Ecliptic latitude magnetic pole
								; Ecliptic longitude Sun
S = big_eph(UT, $
		body	= jpl_body(/sun  ,/string)	, $
		center	= jpl_body(/earth,/string)	, $
		/to_ecliptic	, $
		/to_sphere		, $
		/precess		, $
		degrees=degrees	, $
		/onebody		, $
		/silent)

S = SubArray(S,element=0)		; Geocentric ecliptic longitude of Sun

L -= S							; Ecliptic longitude magnetic pole relative to Sun
								; Angle between projection pole in Y-Z plane and ecliptic North
E = atan(sin(L*rpm)/tan(D*rpm))/rpm

right_angle = 90.0d0/dpm

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

A = SubArray(A, element=0, add=S-right_angle)
A = SubArray(A, element=1, add=-E)
A = SubArray(A, element=2, add=right_angle)

RETURN, A  &  END