;+
; NAME:
;  pluto_phys
; PURPOSE: (one line)
;   calculate physical data for Pluto
; DESCRIPTION:
;   calculate physical data for a target: sub observer and solar lat
;   and lon, phase angle 
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  pluto_phys
; INPUTS:
;  target - target string or code (not used, assumed to be Pluto,
;           included for consistency)
;  obset - ET (TDB s after J2000) at the observer
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  obs - observer code 
;  lonlatlat - [east long in rad, geodetic lat in rad, alt in km]
;  
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;   a structure with
;     subobslat = sub-observer latitude (radian)
;     subobslon = sub-observer longitude (radian)
;     subsollat = sub-solar latitude (radian)
;     subsollon = sub-solar longitude (radian)
;     phase     = observer-target-sun angle (radian)
;     obsdist   = observer-target distance (km)
;     soldist   = sun-target distance (km)
;  
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
; RESTRICTIONS:
;  DOES NOT RETURN OBSERVER INFORMATION ****** 
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2010 Dec 1 Leslie Young, SwRI
;  2011-08-21 LAY. (1) corrected JD of J2000.  
;     (2) used osc. elems for date of last periapse
;  2013-09-25 return orbital_elements
;     Seems to use the OLD osc. elems.
;  2015-07-26 optionally use pck10 ra, dec
;  2016-04-11 move into astronomy IDL library
;-

function pluto_phys, targ, obset, obs=obs, lonlatalt=lonlatalt, trueanom=trueanom, $
                     xyz=xyz, orbital_elem = orbital_elem, pck10=pck10, _extra=extra

    nphys = n_elements(obset)

    if nphys gt 1 then begin
        phys1 = pluto_phys(targ, obset[0], obs=obs, lonlatalt=lonlatalt, pck10=pck10, _extra=extra)
        phys = replicate(phys1, nphys)
        for i=1,nphys-1 do begin
            phys[i] = pluto_phys(targ, obset[i], obs=obs, lonlatalt=lonlatalt, pck10=pck10, _extra=extra)
        end
        return, phys
    end

    rad_per_deg = !dpi/180.d
    km_per_au = 149597870.691          ; km

    ; ORBITAL ELEMENTS.
; PREVIOUS VALUES
    EC= 2.543748319688926E-01 ; eccentricity in deg
    QR= 2.967363028012835E+01 
    IN= 2.343924587525028E+01 ; Inclination w.r.t xy-plane (degrees) 
    OM= 4.395974009773578E+01 ; Longitude of Ascending Node (degrees)
    W = 1.839936675515462E+02 ; Argument of Perifocus (degrees)
    Tp= 2.447899596660452E+06 ; Julian Date of last periapsis
    N = 3.925806519630605E-03 ; Mean motion (deg/day)
    MA= 2.434747256451740E+01 
    TA= 4.091371916521526E+01
    A = 3.979698051030697E+01 ; semimajor axis.
    AD= 4.992033074048560E+01 
    PR= 9.170090227316496E+04

; NEW VALUE - osc. elements generated near to the time of periapse
 ;2447774.894652778 = A.D. 1989-Sep-05 09:28:18.0000 (CT)
;EC= 2.473683243533040D-01 
;QR= 2.965557896088158D+01 
;IN= 1.714940314326779D+01
;OM= 1.102377353401541D+02 
;W = 1.136062519357873D+02 
;Tp=  2447774.894654024858
;N = 3.984907119697048D-03 
;MA= 3.599999999950302D+02 
;TA= 3.599999999914992D+02
;A = 3.940251243797323D+01 
;AD= 4.914944591506487D+01 
;PR= 9.034087600700943D+04

    orbital_elem = [EC, IN, OM, W, TP, N, A]
    
    ; RA AND DEC OF POLE
    rap = 313.02 * rad_per_deg
    decp = 9.09 * rad_per_deg
    if keyword_set(pck10) then begin
       rap = 132.993 * rad_per_deg
       decp = -6.163 * rad_per_deg
    endif

    ; PRIME MERIDIAN
    w0 = 236.77 * rad_per_deg
    w1 = -56.3623195 * rad_per_deg

    days_since_J2000 = obset / (24.d * 3600.d)

    ; Get the heliocentric position of Pluto
    jd2000 = 2451545.0000d  ; Julian date of J2000
    jd = days_since_J2000 + jd2000
    xyz = ephemlay(orbital_elem, jd, trueanom=trueanom, _extra=extra) ; J2000 XYZ in AU
    r = sqrt(total(xyz^2.d))
    ra = (atan(xyz[1], xyz[0])  + 2.d*!pi) mod (2.d * !Pi)
    dec = asin(xyz[2]/r)
    
    ; get the sub solar latitude and longitude
    w = w0 + w1 * days_since_J2000
    sinlat = -sin(decp)* sin(dec) - cos(decp) * cos(dec) * cos(rap - ra)
    lat = asin(sinlat)
    coslat_sink = -cos(decp) * sin(dec) + sin(decp)*cos(dec)*cos(rap-ra)
    coslat_cosk = cos(dec)*sin(rap - ra)
    k = atan(coslat_sink, coslat_cosk)
    lon = (k - w + 2*!pi) mod (2*!pi)

    subsollat = lat
    subsollon = lon
    soldist = r * km_per_au

    subobslat = 0.
    subobslon = 0.
    phase = 0.
    obsdist = 0.

    return, { subobslat:subobslat, $
              subobslon: subobslon, $
              subsollat: subsollat, $
              subsollon: subsollon, $
              phase: phase, $
              obsdist:obsdist, $
              soldist:soldist  $
            }

end