; docformat = 'rst'
;
; NAME:
;    tiltECEQ
; PURPOSE:
;    Compute the angle between the Ecliptic and the Equatorial poles as
;    projected on the plane of the sky (Equatorial pole is the reference)
;
;+
; :Description:
;    Compute the angle between the Ecliptic and the Equatorial poles as
;    projected on the plane of the sky (Equatorial pole is the reference)
;    
; :Categories:
;    Coordinates, Convertion
;    
; :Params:
;    ra: in, required, type=float
;        The Right Ascencion Coordinates (DECIMAL DEGREE)
;    dec: in, required, type=float
;        The Declination Coordinates (DECIMAL DEGREE)
;
; :Returns: The angle between the EQuatorial and ECliptic poles measured
;    on the plane of the sky, counter-clockwise, from the Equatorial pole
;
; :Keywords:
;   degree: in, optional, type=boolean, default=1
;        If set, the angle is reported in degrees
;   radian:  in, optional, type=boolean, default=0
;        If set, the angle is reported in radians
;   sexagesimal:  in, optional, type=boolean, default=0
;        If set, the input RA/DEC are treated as sexagesimal format
;
; :Examples:
;    Prints EC/EQ tilt at vernal point and at several positions::
;        IDL> print, tiltECEQ( 0., 0. )
;        IDL> print, tiltECEQ( 120., 0. )
;        IDL> print, tiltECEQ( 180., +30. )
;
; :Uses:
;
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Written in October 2013, B. Carry (IMCCE)
;      2013-Nov.: B.Carry (IMCCE) - Header in RST
;      2014-Aug.: B.Carry (IMCCE) - Removed dependence on cart2sph
;      2016 Dec.: B.Carry (OCA) - idl2 added
;-
function tiltEcEq, ra, dec, degree=degree, radian=radian, sexagesimal=sexagesimal
 ;--------------------------------------------------------------------------------
 ;--I-- Input verification and Initialization  -----------------------------------
  COMPILE_OPT hidden, idl2

 ;--I.1-- Check Presence of Arguments
  if N_params() LT 2 then begin
    message, /ioError, 'Syntax -  TILT =  tiltEcEq( ra, dec [,/degree, /radian,..])'
    return, -1
 endif
  
 ;--I.2-- Set Output Format
  if not keyword_set(degree) and not keyword_set(radian) then degree=1


 ;--I.3-- Convert Input RA/DEC in Radians 
  if not keyword_set(sexagesimal) then begin
    iRA = !DTOR * ra 
    iDEC= !DTOR * dec
  endif else begin
    iRA = 15. * !DTOR * ten(ra) 
    iDEC=       !DTOR * ten(dec)
  endelse


 ;--------------------------------------------------------------------------------
 ;--II-- Computation of Poles Tilt  ----------------------------------------------

 ;--II.1-- Plane of the Sky (i.e., perpendicular to RA/DEC direction)
  xFP = [ -sin(iDEC)*cos(iRA), -sin(iDEC)*sin(iRA), cos(iDEC) ]
  yFP = [ sin(iRA), -cos(iRA), 0 ]

 ;--II.2-- Earth's obliquity
  epsilon = 23. + 26./60. + 21.4119/3600.
      
 ;--II.3-- Ecliptic Pole in Equatorial
  T3D, /RESET
  T3D, rotate=[epsilon, 0d0, 0d0 ]
  poleEC = Vert_T3D( [[0.,0.,1.],[0.,0.,1.]] )
    poleEC = poleEC[*,0]
    
 ;--II.4-- Projection of Ecliptic Pole on the Plane of the Sky
  poleEC_vec1 = transpose(xFP) # poleEC
  poleEC_vec2 = transpose(yFP) # poleEC

 ;--II.5-- Angle from Equatorial to Ecliptic pole
  tilt = cv_coord( from_rect=[poleEC_vec1, poleEC_vec2], /to_polar, /deg )
    tilt=tilt[0] mod 360.
    if tilt lt 0 then tilt+=360.

 ;--II.6-- Output in DEGREE or RADIAN
  if keyword_set(degree) then return, tilt $
                         else return, tilt*!DTOR

end