;+
;       NAME:
;               oc_plotglobe.Pro
;
;       PURPOSE:
;		Make a globe of the Earth for an occultation
;
;	DESCRIPTION:
;		Constructs a globe as seen from the occulting body at the time of the 
;		occultation.  Also, includes shadow where the sun is below a specified number
;		of degrees below the horizon.
;
;       CALLING SEQUENCE:
;              oc_plotglobe, ev, globe, sunAng=sunang
;
; INPUTS:
;		body   - string with the occultation body name: "Pluto"
;		utc    - utc string giving the time of the occultation, string
;       sunAng - the angle of the sun relative to the horizon for the shadow, deg
;       missDist    - the closest approach distance of the occultation, arcsec
;       paDeg    - the position angle at the closest approach distance, deg
;		radius  - the radius of the occulting body, km
;
;       OUTPUTS:
;
;       REVISON HISTORY:
;       05-April-05 CBO SwRI
;-

pro	oc_plotglobe, globe, sunang=sunang, error=error

naifinit

deg = !dpi / 180.d
theta = findgen(360.*10+1)*deg/10
cos = cos(theta)
sin = sin(theta)


;-------------------------------------------------------------
;plot globe of Earth as seen from Pluto at occultation time

londeg = globe.lonlat_targ[0]/deg
latdeg = globe.lonlat_targ[1]/deg
!p.background = !d.n_colors-1
!p.color = 0
loadct, 0, /silent
map_set, latdeg, londeg, 0., /CONTINENTS, /GRID, /ORTHOGRAPHIC, E_CONTINENTS={FILL:0}, /ISOTROPIC


;-------------------------------------------------------------
; Plot darkness shaded

;The sub solar point on Earth as seen from Pluto
lonsol = globe.lonlat_sun[0]
latsol = globe.lonlat_sun[1]
sinlatsol = sin(latsol)
coslatsol = cos(latsol)
sinlonsol = sin(lonsol)
coslonsol = cos(lonsol)

if not keyword_set(sunang) then sunang = -12. * deg
delta = !pi/2.+sunAng

; x1, y1, z1 describe a cone an angle delta from z axis
x1 = cos(theta) * sin(delta) 
y1 = sin(theta) * sin(delta)
z1 = cos(delta) + (0*theta)

; Rotate x1, y1, z1 to be a cone an angle delta from
; sub-solar point
rotmat1 = [[-sinlatsol, 0, -coslatsol],$
           [0         , 1,         0], $
           [ coslatsol, 0, -sinlatsol] ]

rotmat2 = [ [coslonsol, -sinlonsol, 0], $
            [sinlonsol,  coslonsol, 0], $
            [   0,            0,    1] ]

xyz1 = [ [x1], [y1],[z1] ]           
xyz2 = rotmat1 ## xyz1
x2=xyz2[*,0]
y2=xyz2[*,1]
z2=xyz2[*,2]

xyz3 = rotmat2 ## xyz2
x3=xyz3[*,0]
y3=xyz3[*,1]
z3=xyz3[*,2]

xdark=x3
ydark=y3
zdark=z3

; convert xyz to lat and lon (in degrees for plotting)
latdarkdeg = asin(zdark)/deg
londarkdeg = atan(ydark,xdark)/deg
polyfill, londarkdeg, latdarkdeg, color='808080'xl

;-------------------------------------------------------------
; Replot map over sun shading

map_set, latdeg, londeg, 0., /CONTINENTS, /GRID, /ORTHOGRAPHIC, $
   E_CONTINENTS={FILL:0, usa:1, countries:1}, $
  /ISOTROPIC,/noerase
; prepare to plot w/ x, y in units of Earth radii
xtsave = !x.type
!x.type=0
; plot a circle around the Earth limb
plots, cos(theta), sin(theta)
; return to plotting x, y in units longitude and latitude
!x.type = xtsave

;-------------------------------------------------------------
; Adding occultation shadow path

; Get the Earth radius
cspice_bodvrd, 'Earth', 'RADII', 3, radii

; offset from Earth center, also converts km to Earth radii
offx = sin(globe.pa)/radii[0]
offy = cos(globe.pa)/radii[0]  ; 

; points parallel to track path
trackx =  2. * cos(globe.pa) * [-1,0,1]
tracky = -2. * sin(globe.pa) * [-1,0,1]

; the center track
ctrx = globe.miss_km*offx + trackx
ctry = globe.miss_km*offy + tracky

; prepare to plot w/ x, y in units of Earth radii
xtsave = !x.type
!x.type=0


plots, ctrx, ctry, thick=2, clip=[-1,-1,1,1], noclip=0
for i=-1,1,2 do $
	plots, $
		ctrx + i*globe.path_km*offx, $
		ctry + i*globe.path_km*offy, $
		thick=2, clip=[-1,-1,1,1], noclip=0

; add error to occ path
if keyword_set(error) then begin

  for i=-1,1,2 do $
	  plots, $
	  	ctrx + i*(globe.path_km + globe.misserr_km)*offx, $
	  	ctry + i*(globe.path_km + globe.misserr_km)*offy, $
	  	thick=2,line=2, clip=[-1,-1,1,1], noclip=0
end
!x.type = xtsave

end