;+
; Project     :	Multimission
;
; Name        :	GET_SUNSPICE_PE_ANGLE
;
; Purpose     :	Returns the solar PE angle as seen from a spacecraft
;
; Category    :	SUNSPICE, Orbit
;
; Explanation : This routine returns the solar PE angle as seen from a
;               spacecraft or other solar system body.  The solar PE angle is
;               defined as the angle of the projection of the solar North
;               rotational axis onto the plane of the sky, counterclockwise
;               relative to ecliptic north.  As seen from Earth or one of the
;               two STEREO spacecraft, it varies between approximately +/- 7.25
;               degrees over the course of the year.
;
;               A complement to the PE angle is the B0 angle, defined as the
;               tilt of the solar North axis towards or away from Earth, and
;               which varies between +/- 7.25 degrees over the course of the
;               year.  This parameter is returned as the third value returned
;               by the GET_SUNSPICE_LONLAT routine, e.g.
;
;                       B0 = (GET_SUNSPICE_LONLAT(Date, 'A'))[2,*]
;
; Syntax      :	PE = GET_SUNSPICE_PE_ANGLE( DATE, SPACECRAFT )
;
; Examples    :	Date = '2008-05-06T11:30:00'
;               PE_Ahead = GET_SUNSPICE_PE_ANGLE( Date, 'A' , /DEGREES)
;               PE_Earth = GET_SUNSPICE_PE_ANGLE( Date, 'Earth' , /DEGREES)
;
; Inputs      :	DATE       = The date and time.  This can be input in any
;                            format accepted by ANYTIM2UTC, and can also be an
;                            array of values.
;
;               SPACECRAFT = The name or NAIF numeric code of a spacecraft.
;                            See PARSE_SUNSPICE_NAME for more information about
;                            recognized names.  Can also be the name of a solar
;                            system body, e.g. "Earth", "Mars", "Moon", etc.
;
; Opt. Inputs :	None.
;
; Outputs     :	The result of the function is the PE angle.
;
; Opt. Outputs:	None.
;
; Keywords    : DEGREES = If set, then the angle is returned in units of
;                         degrees, rather than radians.
;
;               CORR = Aberration correction.  Default is 'None'.  Other
;                      possible values are:
;
;                       'LT'    Light travel time
;                       'LT+S'  Light travel time plus stellar aberration
;
;               PRECESS = If set, then ecliptic coordinates are precessed from
;                         the J2000 reference frame to the mean ecliptic of
;                         date.  Only used for HAE/GAE.  Default is PRECESS=0.
;
;               ERRMSG = If defined and passed, then any error messages will be
;                        returned to the user in this parameter rather than
;                        depending on the MESSAGE routine in IDL.  If no errors
;                        are encountered, then a null string is returned.  In
;                        order to use this feature, ERRMSG must be defined
;                        first, e.g.
;
;                               ERRMSG = ''
;                               PE = GET_SUNSPICE_PE_ANGLE( ERRMSG=ERRMSG, ... )
;                               IF ERRMSG NE '' THEN ...
;
;               Will also accept any LOAD_SUNSPICE or ANYTIM2UTC keywords.
;
; Calls       :	GET_SUNSPICE_COORD, CSPICE_RECLAT
;
; Common      :	None.
;
; Restrictions:	This procedure works in conjunction with the Icy/CSPICE
;               package, which is implemented as an IDL Dynamically Loadable
;               Module (DLM).  The Icy source code can be downloaded from
;
;                       ftp://naif.jpl.nasa.gov/pub/naif/toolkit/IDL
;
;               Because this uses dynamic frames, it requires Icy/CSPICE
;               version N0058 or higher.
;
;               Valid ephemeris data must exist for the dates being processed.
;
; Side effects:	Will automatically load the SPICE ephemeris files, if not
;               already loaded.
;
; Prev. Hist. :	Converted from GET_SUNSPICE_P0_ANGLE, which was based on code
;               originally written by Bart de Pontieu, LMSAL
;
; History     :	Version 1, 27-Aug-2019, William Thompson, GSFC
;
; Contact     :	WTHOMPSON
;-
;
function get_sunspice_pe_angle, date, spacecraft, degrees=degrees, corr=k_corr, $
                                precess=precess, errmsg=errmsg, _extra=_extra
;
on_error, 2
if n_params() ne 2 then begin
    message = 'Syntax:  PE = GET_SUNSPICE_PE_ANGLE( DATE, SPACECRAFT )'
    goto, handle_error
endif
;
if n_elements(k_corr) eq 1 then corr = k_corr else corr = 'None'
;
;  Convert DATE into UTC, and determine the dimensionality of the output.
;
utc = anytim2utc(date, /ccsds, _extra=_extra)
ndate = n_elements(utc)
sz = size(utc)
if sz[0] eq 0 then result = 0.0d0 else $
  result = make_array(size=size(utc), /double)
;
;  Parse the name of the spacecraft.
;
sc = parse_sunspice_name(spacecraft, _extra=_extra)
;
;  Make sure that the ephemeris files are loaded.
;
message = ''
load_sunspice, sc, errmsg=message, _extra=_extra
if message ne '' then goto, handle_error
;
;  Get the orientation of the Sun axis in the ecliptic coordinate system.
;
if keyword_set(precess) then frame='ECLIPDATE' else frame='ECLIPJ2000'
cspice_str2et, utc, et
cspice_pxform, 'IAU_SUN', frame, et, rota
eclip_north = [0, 0, 1.0d0]
;
;  Step through the times, and determine the Sun-spacecraft direction.
;
catch, error_status
if error_status ne 0 then begin
    message = !error_state.msg
    catch, /cancel
    goto, handle_error
endif
for i=0,ndate-1 do begin
    sun_north = rota[2,*,i]
    cspice_spkezr, 'Sun', et[i], frame, corr, sc, rad, ltime
    rad = rad[0:2]
    rad = rad / sqrt(total(rad^2, 1))
;
;  Calculate the parts of the solar and ecliptic axes which are perpendicular to
;  the Sun-spacecraft direction.  Renormalize.
;
    snproj = sun_north - total(rad*sun_north) * rad
    snproj = snproj / sqrt(total(snproj^2))
;
    jnproj = eclip_north - total(rad*eclip_north) * rad
    jnproj = jnproj / sqrt(total(jnproj^2))
;
;  Split the solar north projected axis into the parts which are parallel and
;  perpendicular to the ecliptic projected axis.
;
    sxproj = total(snproj*jnproj)
    vecproj = f_crossp(eclip_north, rad)
    vecproj = vecproj / sqrt(total(vecproj^2))
    syproj = total(snproj*vecproj)
;
;  Calculate the PE angle.
;
    result[i] = atan(syproj, sxproj)
endfor
catch, /cancel
;
if keyword_set(degrees) then result = result * (180.d0 / !dpi)
return, result
;
;  Error handling point.
;
handle_error:
if n_elements(errmsg) eq 0 then message, message else $
  errmsg = 'GET_SUNSPICE_PE_ANGLE: ' + message
;
end