; docformat = 'rst'
;
; NAME: 
;    subObserver_Coord
; PURPOSE:
;    Return the planetocentric coordinates of the sub-observer point
;+
; :Description:
;    Return the planetocentric coordinates of the sub-observer point
;
; :Categories:
;    Coordinates, Ephemeris
;    
; :Params:
;    V: in, required, type=float
;        The Ecliptic rectangular coordinates to observer (x,y,z)
;    spin: in, required, type=structure
;        Contain spin information (see readSpin documentation)
;    ep: in, required, type=double/string
;        The epoch at which the ephemeris should be computed
;
; :Returns: A two-element array containing the sub-observer coordinates
;
; :Keywords:
;    PA: out, optional, type=float, default=0
;        The target pole angle (deg, positive North to East)
;    LTC: in, optional, type=boolean, default=0
;        Trigger Light-Time Correction
;
; :Uses:
;     framevector_ec2ver, epRead
;
; :Examples:
;    Compute the sub-Earth point of Ceres on J2000.0
;      IDL> epoch ='2000-01-01T00:00:00.0'
;      IDL> xyz = ssoPositions(1,epoch,type='aster', /geocentric)
;      IDL> spin= readSpin('Ceres.spin')
;      IDL> sep = subObserver_coord( -xyz, spin, epoch)
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in November 2013, B. Carry (IMCCE)
;      2014 Nov. - B.Carry (IMCCE) - Header adapted to rst documentation
;      2015 Oct. - B.Carry (OCA) - Added idl2 compile option
;      2018 June - B.Carry (OCA) - Computation of Lon/Lat simplified
;-
function subObserver_Coord, V, spin, ep, PA=PA, LTC=LTC
  COMPILE_OPT hidden, idl2
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---    I -- Initialization And Input Verification              -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--I.1-- Input Vector  -----------------------------------------------------------------
  vLen = sqrt( total( V^2 ) )

 ;--I.2-- Input Spin Structure  ---------------------------------------------------------
  spinFields = tag_names(spin)
  spinRequis = ['LON', 'LAT', 'PD', 'JD0', 'PHI0']
  nbRequis = n_elements(spinRequis)
  for k=0, nbRequis-1 do begin
    pField=where( strCmp(spinRequis[k], spinFields, /FOLD) )
    if pField[0] eq -1 then begin
      message, 'Field '+strTrim(spinRequis[k],2)+' missing in spin structure'
      return, -1
    endif
  endfor

 ;--I.3-- Input Epoch -------------------------------------------------------------
  jd=epRead(ep)

 ;--I.4-- Physical Constants  -----------------------------------------------------
  AUtoKM = 149597870.691d0   ;-km
  cKM_S = 299792.458D        ;-km/s
    cAU_S = cKM_S / AUtoKM   ;-au/s
    cAU_D = cAU_S * 3600.D * 24.D 


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---   II -- Compute the Model Geometry                         -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;

 ;--II.1-- Light-Time Correction  ---------------------------------------------------
  if keyword_set(LTC) then LTC = vLen/cAU_D else LTC = 0D


 ;--II.2-- Convert V into Shape Model Frame  ----------------------------------------
  V_ast = framevector_ec2ver( V, spin, JD-LTC )
  vLen = sqrt( total( V_ast^2 ) )


;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
;---  TAG ---  III -- Coordinates and Pole Angle                         -----------------------;
;-----------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------;
  
 ;--III.1-- Planetocentric Longitude  ------------------------------------------------
  sopLon = (atan( v_ast[1], v_Ast[0] )/!DTOR + 360) mod 360

 ;--III.2-- Planetocentric Latitude  -------------------------------------------------
  sopLat = asin( v_ast[2] / vLen )/!DTOR

 ;--III.3-- Pole Angle  (direct sense)  ----------------------------------------------
  spinLon = spin.lon*!DTOR
  spinLat = spin.lat*!DTOR

  targLon = atan( -V[1], -V[0] )
  targLat = asin( -V[2] / vLen )

  PNum = -cos(spinLat)*sin(targLon-spinLon)
  PDen =  sin(spinLat)*cos(targLat)-cos(spinLat)*sin(targLat)*cos(targLon-spinLon)
  PA = (atan(PNum,PDen) /!DTOR + 360) mod 360

 ;--III.4-- Return Sub-Observer's Coordinates  ----------------------------------------
  return, {lon:sopLon, lat:sopLat}

end