pro sza_plot,date,latpole=lat00,lonpole=lon00,kiruna=kiruna,takeoff_gmt=takeoff_gmt
;+ this routine does a sza plot for the location
; clidk on map to find the location
; plot appears in second window
; takeoff_gmt is in hours
;-
if n_elements(lat0) eq 0 then begin
   lat00=67.84
   lon00=20.41
   city='Kiruna'
   local_time=1
   endif
   if keyword_set(kiruna) then begin
   lat00=67.84
   lon00=20.41
   city='Kiruna'
   endif
lat_path=[lat00]
lon_path=[lon00]
flttime=[0]

if n_elements(takeoff_gmt) eq 0 then takeoff_gmt=12      
set_plot,'MAC'
airspeed=600 ; km/hr
device,/retain,/cursor_orig
win0=!d.window+1
win1=!d.window+2
win2=!d.window+3
window,win1,xsize=400,ysize=300,title='Local zenith angle window'
window,win2,xsize=400,ysize=300,title='Flight window'
plot,[0],[0],title='Solar/Lunar Zenith Angle '+nicedate(date),xrange=[0,8],yrange=[0,150],xticks=8,xtitle='Flight Hours',$
ytitle='Zenith Angle'
window,win0,xsize=400,ysize=400, title='Map Window'

erase


rr:   ; draw the map of the region
     map_set,lat00,lon00,/lambert,scale=3.e7
     stdcolor
     map_continents
     map_grid
     for j=400,2400,400 do begin
         range_ring,lat00,lon00,j,360,bearing,latp,lonp
         oplot,lonp,latp,col=1
         endfor
     oplot,[lon00],[lat00],psym=1,col=3
     xyouts,lon00,lat00,city
     
     
xx:  ; find a location
     cursor,lon0,lat0,/data,/down
     if finite(lon0) and finite(lat0) then begin
     print,lon0,lat0
     ns=n_elements(lat_path)
     oldlat=lat_path(ns-1)
     oldlon=lon_path(ns-1)
     oplot,[oldlon,lon0],[oldlat,lat0],col=4
     lat_path=[lat_path,lat0]
     lon_path=[lon_path,lon0]
     gr_circ_rte,oldlat,oldlon,lat0,lon0,2,bearing,dist,del,latp,lonp,rd
     flttime=[flttime,flttime(ns-1)+dist/airspeed]
     ; switch to the other window to show local zenith angle
     wset,win1
     hour=findgen(48)*.5
     yr=fix(strmid(date,0,2))
     mon=fix(strmid(date,2,2))
     day=fix(strmid(date,4,2))

     sz=solar_zenith(lat0+hour*0.,lon0+hour*0.,mon=mon,day=day,year=yr,hour=hour)
     szm=solar_zenith(lat0+hour*0.,lon0+hour*0.,mon=mon,day=day,year=yr,hour=hour,/moon)
     plot,hour,sz,title='Solar/Lunar Zenith Angles '+nicedate(date)+' at (lat,lon) '+ $
        str(lat0) + ',' + str(lon0),xtitle='GMT Hours',$
        ytitle='Zenith Angle',yst=1,xst=1,xrange=[0,24],xticks=6,yrange=[0,150],yticks=5
        oplot,[0.,24.],[94,94]
        oplot,hour,szm,linestyle=1
        sz2=(sz-94.)^2
        n=where(min(sz2) gt sz2-1)
        if n(0) ge 0 then $
           for j=0,n_elements(n)-1 do oplot,hour(n(j))*[1,1],[0,120]
           
     ; switch to last window
     
     wset,win2
    szz=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime) 
     szzm=solar_zenith(lat_path,lon_path,mon=mon,day=day,year=yr,hour=takeoff_gmt+flttime,/moon) 
   plot,flttime,szz,title='Solar/Lunar Zenith Angle '+nicedate(date),xrange=[0,8],yrange=[0,150],xticks=8,xtitle='Flight Hours',$
       ytitle='Zenith Angle'
      oplot,flttime,szzm,linestyle=1
        if flttime(n_elements(flttime)-1) gt 8. then goto,quitit

 
 wset,win0
     map_set,lat00,lon00,/lambert,/noerase,scale=3.e7
       goto,xx
      endif 
 print,'out of range'
quitit:
print, 'Way Points
print,'     Time       Lat       Lon  '

for j=0,n_elements(flttime)-1 do print,flttime(j),lat_path(j),lon_path(j)
return
end