;+
;       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, f0=f0, g0=g0

naifinit

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


;set predicted time of the occultation
;et = ev.et

londeg = globe.lonlat_targ[0]/deg
latdeg = globe.lonlat_targ[1]/deg

;plot globe of Earth as seen from Pluto at occultation time
!p.background = !d.n_colors-1
!p.color = 0
loadct, 0
map_set, latdeg, londeg, 0., /CONTINENTS, /GRID, /ORTHOGRAPHIC, E_CONTINENTS={FILL:0}, /ISOTROPIC


;Calculate the sub solar point on Earth as seen from Pluto
lonsoldeg = globe.lonlat_sun[0]/deg
latsoldeg = globe.lonlat_sun[1]/deg

if not keyword_set(sunang) then sunang = -12.
deltadeg = 90.+sunAng
delta = deltadeg*deg
theta = findgen(360.+1)*!pi/180
x1 = cos(theta) * sin(delta)
y1 = sin(theta) * sin(delta)
z1 = cos(delta) + (0*theta)

sinlatsol = sin(latsoldeg * deg)
coslatsol = cos(latsoldeg * deg)
rotmat1 = [[-sinlatsol, 0, -coslatsol],$
           [0         , 1,         0], $
           [ coslatsol, 0, -sinlatsol] ]
sinlonsol = sin(lonsoldeg * deg)
coslonsol = cos(lonsoldeg * deg)
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
latdarkdeg = asin(zdark)/deg
londarkdeg = atan(ydark,xdark)/deg
polyfill, londarkdeg, latdarkdeg, color='808080'xl

map_set, latdeg, londeg, 0., /CONTINENTS, /GRID, /ORTHOGRAPHIC, $
   E_CONTINENTS={FILL:0, usa:1, countries:1}, $
  /ISOTROPIC,/noerase



; Adding occultation shadow path

if not keyword_set(f0) then f0 = 0.
if not keyword_set(g0) then g0 = 0.
missDistkm = globe.miss_km
cspice_bodvrd, 'Earth', 'RADII', 3, radii
missDistEarthRadii = missDistkm/radii[0]
radiusEarthRadii = globe.path_km/radii[0]
pa = globe.pa

y = missDistEarthRadii*COS(pa) + g0/radii[0]
x = missDistEarthRadii*SIN(pa) + f0/radii[0]
offy = cos(pa)/radii[0]
offx = sin(pa)/radii[0]
tracky = -2. * sin(pa) * [-1,1]
trackx =  2. * cos(pa) * [-1,1]


;incpt = missDistEarthRadii*(COS(pa) + TAN(pa)*SIN(pa))
;ypt1 = -TAN(pa) + incpt + g0/radii[0]
;ypt2 = TAN(pa) + incpt + g0/radii[0]

xtsave = !x.type
!x.type=0

plots, x + trackx, y + tracky
 

;plots, [[-1, ypt2], [1, ypt1]], thick=2
;plots, [[-1, ypt2+radiusEarthRadii], [1, ypt1+radiusEarthRadii]], thick=2
;plots, [[-1, ypt2-radiusEarthRadii], [1, ypt1-radiusEarthRadii]], thick=2
;plots, cos(theta), sin(theta)
;;plots, [x,y], ps=4


; add error to occ path
if keyword_set(error) then begin
errradiusEarthRadii = globe.path_km/radii[0] + globe.misserr_km/radii[0]

plots, [[-1, ypt2], [1, ypt1]], thick=2
plots, [[-1, ypt2+errradiusEarthRadii], [1, ypt1+errradiusEarthRadii]], $
  thick=2,line=2
plots, [[-1, ypt2-errradiusEarthRadii], [1, ypt1-errradiusEarthRadii]], $
  thick=2,line=2
plots, [x,y], ps=4
end
!x.type = xtsave

end