;+
; NAME:
;  naif_phys
; PURPOSE: (one line)
;   calculate physical data for a target
; DESCRIPTION:
;   calculate physical data for a target: sub observer and solar lat
;   and lon, phase angle 
; CATEGORY:
;  Astronomy
; CALLING SEQUENCE:
;  naif_phys
; INPUTS:
;  target - target string or code
;  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:
; PROCEDURE:
; MODIFICATION HISTORY:
;  Written 2006 May 11 Leslie Young, SwRI
;  2007 June 06 LAY - corrected documentation
;  2008 Mar 25 LAY - corrected sign error in subsollon
;  2011 Jan 13 LAY - Added calls to naif_id to all targ to be string
;-

function naif_phys, targ, obset, obs=obs, lonlatalt=lonlatalt


    naifinit

    targ_code = naif_id(targ)

    nphys = n_elements(obset)

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

    obs_targ = oc_naifstate(obset, targ_code, obs=obs, lonlatalt=lonlatalt,ltime=ltime)
;    CSPICE_SPKEZ, targ, obset, 'J2000', 'LT', 399, obs_targ, ltime
    obsdist = sqrt(total(obs_targ[0:2]^2,/double))

    et = obset - ltime

    sun = 10 ; code for sun
    targ_sun = oc_naifstate(et, sun, obs=targ_code, ltime=ltime)
;    CSPICE_SPKEZ, 10, et, 'J2000', 'LT', targ, targ_sun, ltime
    soldist = sqrt(total(targ_sun[0:2]^2,/double))
    
;
;        Look up the radii of the target body from the PCK file.
;
    cspice_bodvar, targ_code, "RADII", radii
    a = radii[0] 
    b = radii[1]
    c = radii[2]

;           Get the transformation between J2000 and body reference frame
    cspice_bodc2n, targ_code, targ_name, found 
    bodyframe = "IAU_" + targ_name
    cspice_pxform, bodyframe, 'J2000', et, r_body2j2000

;           Negate the observer-target vector and rotate the result
;           into body-fixed coordinates.
;   
    obspos = r_body2j2000 # (-obs_targ[0:2])

;        Find the nearest point on the body to the observer.
    cspice_nearpt, obspos, a, b, c, subobs, alt   

;        Find the latitude and longitude of the sub-observer point.
;    cspice_reclat, subobs,rr, subobslon, subobslat
    cspice_recgeo, subobs, a, (a-c)/a, subobslon, subobslat, alt
    subobslon = (2.*!dpi + subobslon) mod (2.*!dpi)

    sunpos = r_body2j2000 # (targ_sun[0:2])

;        Find the nearest point on the body to the observer.
    cspice_nearpt, sunpos, a, b, c, subsol, alt   

;        Find the latitude and longitude of the sub-sun point.
;    cspice_reclat, subsol,rr, subsollon, subsollat
    cspice_recgeo, subsol, a, (a-c)/a, subsollon, subsollat, alt
    subsollon = (2.*!dpi + subsollon) mod (2.*!dpi) ; 2008 mar 25 - had been (2*!dpi - subsollon)

    phase = vang(-obs_targ[0:2], targ_sun[0:2])

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

end