;
pro pltst,info,animate=animate,ps=ps
;
; Contour and/or image map.data(nlon,nlat) over polar stereographic
;   projection.
;
file = *info.pfile
map = *info.pmap
vector=map.vector
field = map.field
data = *map.data
lons = *map.lons & nlon = n_elements(lons)
lats = *map.lats & nlat = n_elements(lats)
;
; Add periodic points (add a longitude):
;
   newdata = fltarr(nlon+1,nlat)
   newdata[0:nlon-1,*] = data[*,*]
   newdata[nlon,*] = data[0,*]
   data = newdata
   lons = fltarr(nlon+1)
   lons[0:nlon-1] = *map.lons
   lons[nlon] = 180.
   nlon = nlon+1
;
; Check for all missing value:
;
nmissing = 1
text = 1 & if keyword_set(animate) then text = 0
nodata = check_missing(data,float(map.missing_value),field.name,text=text,$
  nmissing=nmissing)
if nmissing gt 0 then print,'Polar maps: Field ',field.name,': Found ',nmissing,$
  ' missing data points (out of ',n_elements(data),')'
if nodata eq 1 then return
missing_value = float(map.missing_value)
;
; Set map position:
;
wxsize = !d.x_size & wysize = !d.y_size
dia = 0.70*wysize
xcen = 0.42*(wxsize) & ycen = 0.53*wysize
xll_dev = xcen-0.5*dia & yll_dev = ycen-0.5*dia
norm = convert_coord([xll_dev,yll_dev],/device,/to_norm)
xll = norm(0) & yll = norm(1)
xur_dev = xcen+0.5*dia & yur_dev = ycen+0.5*dia
norm = convert_coord([xur_dev,yur_dev],/device,/to_norm)
xur = norm(0) & yur = norm(1)
mappos = [xll,yll,xur,yur]
xmid = .5*(xll+xur)
;
; Set up stereographic projection. 
; If fixpolar='slt', then set local noon at top of projection.
; If fixpolar='lon', then set longitude 0 at top of projection.
; Note that lon at center of projection is up (towards top) for
;   S hem, but down (towards bottom) for N hem.
; If magnetic field, add 70 deg rotation from geographic.
;
projcen = fltarr(2) ; lat,lon of projection center
uthrs = float(map.mtime[1])+float(map.mtime[2])/60. ; ut (hrs)
if map.stperim le 0. then begin ; south hemisphere
  hem = 'S'
  projcen[0] = -90.

  mag = 0 & if field.grid_type eq 'magnetic' then mag = 1
  if map.fixpolar eq 'SLT' then begin  ; local noon constant at top of proj
    projcen[1] = fslt(12.,uthrs,dum,'getlon',mag=mag)
  endif else begin
    projcen[1] = 0.             ; longitude 0 constant at top of proj
  endelse

  map_set,projcen[0],projcen[1],/stereo,position=mappos,/noborder,$
    limit=[-90.,-180.,map.stperim,180.]
endif else begin                ; north hemisphere
  hem = 'N'
  projcen[0] = 90.

  mag = 0 & if field.grid_type eq 'magnetic' then mag = 1
  if map.fixpolar eq 'SLT' then begin  ; local noon constant at top of proj
    projcen[1] = fslt(0.,uthrs,dum,'getlon',mag=mag) 
  endif else begin
    projcen[1] = -180.          ; longitude 0 constant at top of proj
  endelse

  map_set,projcen[0],projcen[1],/stereo,position=mappos,/noborder,$
    limit=[map.stperim,-180.,90.,180.]
endelse
;
fminmax,data,framemin,framemax,missing_value
if map.fmin ge map.fmax then begin
  fminmax,data,fmin,fmax,missing_value
  if map.fmin eq map.fmax then begin
    map.fmin=fmin & map.fmax=fmax
    if not keyword_set(animate) then $
      print,'Polar maps reset map.fmin,max=',map.fmin,map.fmax
  endif else $
    if not keyword_set(animate) then $
      print,'Polar maps did NOT reset min,max: map.fmin,max=',$
        map.fmin,map.fmax
endif else begin
  fmin = map.fmin & fmax = map.fmax
; if not keyword_set(animate) then $
;   print,'Polar maps using preset map.fmin,max=',map.fmin,map.fmax
endelse
;
; Prepare image:
;
if (map.plottype eq 'image_only') or (map.plottype eq 'image+contours') then begin
;
; Byte-scale to an image. Make any missing data black:
;
  missing_indices = where(data eq map.missing_value,nmissing)
  if nmissing gt 0 then begin
    ave = total(data(where(data ne map.missing_value)))/n_elements(data)
    data(missing_indices) = ave
    image = bytscl(data,top=!d.table_size-3,min=fmin,max=fmax)+1
    image(missing_indices) = 0
  endif else $
    image = bytscl(data,top=!d.table_size-3,min=fmin,max=fmax)+1
;
; outside_proj = 0 ; color index for areas outside projection
; outside_proj = 255 ; color index for areas outside projection
  outside_proj = !p.background ; color index for areas outside projection
  image = map_patch(image,lons,lats,xstart=startx,ystart=starty,$
    xsize=xsize,ysize=ysize,lat0=lats[0],lat1=lats[nlat-1],$
    missing=outside_proj)
; image = map_image(image,startx,starty,xsize,ysize,$
;   latmin=lats[0],latmax=lats[nlat-1],$
;   lonmin=lons[0],lonmax=lons[nlon-1],/compress)
;
; For color bar (pro clrbar):
;
  barw = .04 & barh = yur-yll
  xbar = xur+.05 & ybar = yll
  barlabs = [strcompress(string(fmin,format="(e8.1)"),/remove_all),$
    strcompress(string(0.5*(fmin+fmax),format="(e8.1)"),/remove_all),$
    strcompress(string(fmax,format="(e8.1)"),/remove_all)]
  chsize_clrbar = 1.2
  laboff_clrbar = .015
endif ; plottype is image_only or image+contours
;
; Field label:
;
if vector.control ne 'only' then begin
  if (map.log10 eq 'density fields' and isdensity(field.name) gt 0) or $
     (map.log10 eq 'this frame only') or $
     (map.log10 eq 'all fields') then begin
    if strlen(field.long_name) gt 0 then fieldlab = 'LOG10 '+field.long_name $
      else fieldlab = 'LOG10 '+field.name
  endif else begin
    if strlen(field.long_name) gt 0 then fieldlab = field.long_name $
      else fieldlab = field.name
  endelse
  ;if strlen(field.long_name) gt 0 then fieldlab = field.long_name $
  ;  else fieldlab = field.name
  if file.diffs le 0 then begin
    if strlen(field.units) gt 0 then $
      fieldlab = fieldlab + ' (' + strcompress(field.units,/remove_all) + ')'
  endif else fieldlab = fieldlab + ' (' + field.difftype + ' ' + $
    strcompress(field.units,/remove_all) + ')'
endif else begin
  fieldlab= vector.lonvector + '+' + vector.latvector
endelse
;
; Top label (under field lab):
;
uthrs = float(map.mtime[1])+float(map.mtime[2])/60.
if map.plotz eq 0 then toplab = string(format="('DAY = ',i3,'  UT = ',f5.2)",map.mtime[0],uthrs)
if field.nlev gt 0 and map.plotz eq 0 then $
  toplab = toplab + string(format="('  ',a,' = ',f6.2)",$
    field.levshortname,map.zvert)
;
;  Change vertical level label for geopotential altitude plotting
;
if map.plotz ge 1 then toplab = string(format="('DAY=',i3,' UT=',f5.2)",map.mtime[0],uthrs)
if field.nlev gt 0 and map.plotz eq 1 then $
  toplab = toplab + string(format="(' ',a,'=',f6.2,a)",$
    'Geopotential Altitude',map.alt,'km')
if field.nlev gt 0 and map.plotz eq 2 then $
  toplab = toplab + string(format="(' ',a,'=',f6.2,a)",$
    'Geometric Altitude',map.alt,'km')
if field.nlev gt 0 and map.plotz eq 3 then $
  toplab = toplab + string(format="(' ',a,'=',f6.2,a)",$
    'Geometric Altitude (ZG)',map.alt,'km')

if map.plotz eq 0 then toplab = toplab + string(format="('  PERIMLAT = ',f5.1)",map.stperim)
if map.plotz eq 1 then toplab = toplab + string(format="(' PERIMLAT=',f5.1)",map.stperim)
;
; Min,max label:
;
if vector.control ne 'only' then begin
  minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",framemin)+','+$
                         string(format="(g10.4)",framemax)
  if map.plottype ne 'image_only' then $
    minmaxlab = minmaxlab + string(format="(' INTERVAL = ',g10.4)",map.cint)
endif else begin
  lonvecdata=*vector.lonvecdata
  latvecdata=*vector.latvecdata
  winddata=windminmax(lonvecdata,latvecdata,vector.missing_value)
  *vector.vmin=winddata[0]
  *vector.vmax=winddata[1]
  *vector.vint=winddata[2]
  minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",*vector.vmin)+','+$
			      string(format="(g10.4)",*vector.vmax)
  ;minmaxlab = minmaxlab + string(format="(' INTERVAL = ',g10.4)",*vector.vint)
endelse
;
; Bottom labels:
;
filelab = strcompress(field.ncfile,/remove_all)
;
; Contour levels were made by mkclevels.pro.
;
clevels = *map.clevels
nlevels = n_elements(clevels)
clinestyle = intarr(nlevels)
clinestyle[*] = 0 ; solid
for i=0,nlevels-1 do if clevels[i] lt 0. then clinestyle[i] = 2 ; dashed
c_charsize = 1.1
 chsize = 1.2
;
; get wind vector data
;
if vector.control ne 'off' then begin
; zonal data
  zdata = get_field_data(info,vector.ilonvec) 
  *vector.lonvecdata=zdata
  cmin=0. & cmax=0. & cint=0.
  fminmax,zdata,fmin,fmax,map.missing_value
  zlevels = mkclevels(fmin,fmax,cmin,cmax,cint)
  vector.lonint=cint
    
; meridonal data
  mdata = get_field_data(info,vector.ilatvec) 
  *vector.latvecdata=mdata
  cmin=0. & cmax=0. & cint=0.
  fminmax,mdata,fmin,fmax,map.missing_value
  mlevels = mkclevels(fmin,fmax,cmin,cmax,cint)
  vector.latint=cint
    
  *map.data=data
  map.vector=vector
  *info.pmap = map
endif

;
; Display image and/or contour:

;
case map.plottype of
;
; Display image only. Use contour,/nodata to draw axes, ticmarks, labels:
;
  'image_only': begin
    if vector.control ne 'only' then begin
      tv,image,startx,starty,xsize=xsize,ysize=ysize
      clrbar,barw,barh,xbar,ybar,barlabs,botcolor=1,topcolor=!d.table_size-2,$
        charsize=chsize_clrbar,laboffset=laboff_clrbar
      if vector.control eq 'on' then begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
      endif
    endif else begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
;   axes:
;      contour,data,lons,lats,/noerase,position=mappos,/nodata,$
;        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
;        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
;        yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
    endelse
    xyouts,xmid,yur+.08,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.05,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.12,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.16,filelab ,/norm,align=0.5,charsize=1.0
  end
;
; Display image with contour overlay:
;
  'image+contours': begin
    if vector.control ne 'only' then begin
      tv,image,startx,starty,xsize=xsize,ysize=ysize
      contour,data,lons,lats,/noerase,position=mappos,$
        /follow,/overplot,xstyle=5,ystyle=5,c_charsize=c_charsize,$
        levels=clevels,c_linestyle=clinestyle,c_colors=[map.clineclr]
      clrbar,barw,barh,xbar,ybar,barlabs,botcolor=1,topcolor=!d.table_size-2,$
        charsize=chsize_clrbar,laboffset=laboff_clrbar
      if vector.control eq 'on' then begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
      endif
    endif else begin
;   axes:
;      contour,data,lons,lats,/noerase,position=mappos,/nodata,$
;        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
;        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
;        yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
    endelse
    xyouts,xmid,yur+.08,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.05,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.12,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.16,filelab  ,/norm,align=0.5,charsize=1.0
  end
;
; Contour only (B&W):
;
  'monochrome_contours': begin
    if vector.control ne 'only' then begin
      contour,data,lons,lats,/noerase,position=mappos,$
        /follow,/overplot,xstyle=5,ystyle=5,c_charsize=c_charsize,$
        levels=clevels,c_linestyle=clinestyle,c_colors=[map.clineclr]
      if vector.control eq 'on' then begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
      endif
    endif else begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
;   axes:
;      contour,data,lons,lats,/noerase,position=mappos,/nodata,$
;        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
;        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
;        yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
    endelse
    xyouts,xmid,yur+.08,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.05,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.11,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.14,filelab  ,/norm,align=0.5,charsize=1.0
  end
;
; Color fill contour:
;
  'colorfill_contours': begin
    if vector.control ne 'only' then begin
      contour,data,lons,lats,/noerase,position=mappos,$
        /follow,/cell_fill,/overplot,xstyle=5,ystyle=5,levels=clevels
      contour,data,lons,lats,/noerase,position=mappos,$
        /follow,/overplot,xstyle=5,ystyle=5,c_charsize=c_charsize,$
        levels=clevels,c_linestyle=clinestyle,c_colors=[map.clineclr]
      if vector.control eq 'on' then begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
      endif
    endif else begin
        vmagnitude=gridvec(*vector.lonvecdata,*vector.latvecdata,*vector.lonvecs, $
	  *vector.latvecs,map.projection, $
	  spval=map.missing_value,position=mappos, perimlat=map.stperim, $
          scalemag=vector.vmag,scalelen=vector.vlen,xscale=.9,yscale=.08, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
          vector.vmag=vmagnitude
;   axes:
;      contour,data,lons,lats,/noerase,position=mappos,/nodata,$
;        xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1,ticklen=-.015,$
;        xticks=n_elements(xticlabs)-1,xminor=3,xtickname=xticlabs,$
;        yticks=n_elements(yticlabs)-1,yminor=3,ytickname=yticlabs
    endelse
    xyouts,xmid,yur+.08,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.05,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.13,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.16,filelab  ,/norm,align=0.5,charsize=1.0
  end
  else: print,'>>> Polar maps: unknown plot type ',map.plottype
endcase
;
; Add continental outlines and/or map grid if not magnetic grid:
;
if (map.continents eq 1 or map.grid eq 1) and $
   field.grid_type ne 'magnetic' then begin
  if map.continents eq 1 then map_continents,color=map.continents_color
  if map.grid eq 1 then map_grid,/label,color=map.grid_color
endif else if (map.continents eq 1 or map.grid eq 1) and $
  field.grid_type eq 'magnetic' then begin
  if map.continents eq 1 then begin
    print,'WARNING: cannot plot continents on magnetic grid.'
    map.continents = 0
  endif
  if map.grid eq 1 then begin
    print,'WARNING: cannot plot map grid on magnetic grid.'
    map.grid = 0
  endif
endif
;
; Local time axes and labeling for polar plot:
;
labpol_chsize = 1.0
mag = 0 & if field.grid_type eq 'magnetic' then mag = 1
labpol,mappos,hem,labpol_chsize,map.fixpolar,mag=mag
;
; map.fmin,max and map.continents may have been changed, 
; so update state:
;
map.vector=vector
*info.pmap = map
end