pro titanmap

;  Plots HST basemap for Cassini coverage studies etc.
;  R. D. Lorenz  1999-2001
;  NB some early versions of IDL have a bug in the map_image 
;  routine
;
; the basemap is derived from the HST map published in
; P. H. Smith, M.T. Lemmon, R. D. Lorenz, L.A. Sromovsky, 
; J. Caldwell and M.D. Allison, šTitanšs Surface, Revealed by 
; HST Imagingš, Icarus  vol.119, pp.336-349, 1996
;
; the basemap is made available for planning purposes only, not 
; for intercomparison with other observations, for which 
; I should be contacted at rlorenz@lpl.arizona.edu


; NB when plotting footprints etc., the map x coordinate is longitude
; EAST.  Thus to plot a point at 90 W, 20 N,  you need a command like
; plots, 360.-90, 20, color=255
; the map has been artificially labelled with longitude WEST, as this
; is the Cassini Surfaces Working Group convention  (and indeed that 
; of IAU etc) 




; tit is title string - pick your own!

tit= 'Titan'

;  lim represents the minimum latitude, longitude and maximum latitude,
; longitude for the map.  Change these values to 'zoom' in on a
; particular region

lim=[-90.,0.,90.,360]

; plot_landing is a flag for displaying the Huygens landing sites
; set to zero to suppress

plot_landing=1

;  you can annotate the map with oplot commands, xyouts etc.
; since the basemap uses black-blue-lightblue colors, best
; annotation or footprint colors are yellow/red/white
; i.e.  color=150 .... color=255


; note that you can generate a gif file for inclusion on websites or
; documents after you get things to your satisfaction on-screen
; with the commands
; a=tvrd()
; write_gif,'filename.gif',a

;(More recent versions of IDL do not support gifs - how lame is 
;that? Software capability is supposed to improve with time!  
; You can do write_bmp, or tiff, etc....


; enjoy!


capcol=255


; commented out - size the window to your own needs
;window,0,xsize=760, ysize=560,$
;title='R.D.Lorenz  RADAR Science Planning'


;  Reads in Map array

z=fltarr(120,61)
shi=60
x=fltarr(120,61)

; assumes 'hstmap.dat' in the current directory - change if that will not be
; the case

openr,1,'hstmap.dat'
readf,1,x
close,1
a=fltarr(120,61)
a=x(*,*,0)
for i=0,119 do begin
for j=0,60 do begin
if a(i,j) gt 800 then a(i,j)=-120
endfor
endfor
mm=min(a)
a=a-mm


z=(a/(max(a)))
z=z*110.


; map set command for established convention with 180 longitude at
; center, and cylindrical projection

map_set, 0.,180., /cyl,   title=tit, limit=lim

z=shift(z,shi,0)
z(*,0:8)=0
for i=0.,15 do begin
z(*,i)=z(*,i)*i/15.
z(*,60-i)=z(*,60-i)*i/15.
endfor
;qq=fltarr(61)
;qq(*)=z(119,*)+z(0,*)
;qq=qq/2.
;z(119,*)=qq(*)
;z(0,*)=qq(*)

z=z(*,0:59)
col1=254
col2=240
col3=220
col4=200
col5=180
lons=findgen(120)*(-3.0)
lats=findgen(60)*(3.0)-90.0
lons=180.-lons
loadct,39

newimage=map_image(z,  startx, starty,xsize, ysize, sc=0.04, /compress)
tv,newimage,startx, starty, xsize=xsize*1.20, ysize=ysize


map_grid,color=capcol
ss=1.4
xyouts,175,0,'180W',color=capcol, size=ss
xyouts,264,0,'90W',color=capcol, size=ss
xyouts,50,30,'30N',color=capcol, size=ss
xyouts,50,-30,'30S',color=capcol, size=ss
xyouts,0,0,'360W',color=capcol, size=ss
xyouts,85,0,'270W',color=capcol, size=ss

for i=-90.,90.,2. do begin
plots,360,i,psym=3
endfor

xyouts, 10,-88, 'Basemap : HST 950nm Map from (Smith et al., 1996)', color=capcol, size=0.9

if (plot_landing eq 1) then begin

lat=9.4
lon=156.
th=1
c1=200
laterr=0.26 * 3.
lonerr=2.2 * 3.
aa=findgen(30) * 2.*!pi/28.



lat=-9.1
lon=168.
xx=lon + lonerr*cos(aa)
;xx=360.-xx
yy=lat + laterr*sin(aa)
oplot,xx,yy, thick=th, color=c1

endif


end