;-------------------------------------------------------------------------
pro addmapgrid,map,field
;
; 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
  print,'WARNING: cannot plot continents or map grid on magnetic grid.'
  map.continents = 0
  map.grid = 0
endif
end
;-----------------------------------------------------------------------
;
pro pltce,info,animate=animate
;
file = *info.pfile
map = *info.pmap
fields = *info.fields
field = map.field
data = *map.data
lons = *map.lons & nlon = n_elements(lons)
lats = *map.lats & nlat = n_elements(lats)
vector=map.vector
;
;print,' '
;print,'pltce:'
;print,'field.name = ',field.name
;print,'map.projection = ',map.projection
;print,'map.zvert=',map.zvert
;print,'map.mtime=',map.mtime
;print,'map.missing_value=',map.missing_value
;print,'map.stperim=',map.stperim
;print,'nlon=',nlon,' lons=' & print,lons
;print,'nlat=',nlat,' lats=' & print,lats
;print,'data min,max = ',min(data),',',max(data)
;print,format="('data[*,18]=',/,(6e12.4))",data[*,18]
;
; Write binary file for debug:
;file='for_rsi.dat'
;openw,lu,file,/get_lun
;writeu,lu,nlon
;writeu,lu,lons
;writeu,lu,nlat
;writeu,lu,lats
;writeu,lu,data
;close,lu
;free_lun,lu
;print,'Wrote binary file ',file,' lu=',lu
;
; Check for all missing value:
;
nmissing = 1 ; do not init this to 0 because check_missing uses "if keyword_set"
text = 1 & if keyword_set(animate) then text = 0
nodata = check_missing(data,float(map.missing_value),field.name,text=text,$
  nmissing=nmissing)
if nodata eq 1 then return
missing_value = float(map.missing_value)
if nmissing gt 0 then print,'Global Maps: Field ',field.name,': Found ',$
  nmissing,' missing data points (out of ',n_elements(data),')'
;
; Set map position (normalized coords):
;
xll = .15 & yll = .35 
xur = .95 & yur = .85
mappos = [xll,yll,xur,yur]
dlat = lats[1]-lats[0]
;
; Geographic field lats start at +/- 87.5, whereas magnetic fields go to +/- 90.
;
if (lats[0]-dlat/2. eq -90.) and (lats[nlat-1]+dlat/2. eq +90.) then begin
  poledeg = (dlat/2.)/180.*(yur-yll) ; distance between top lat and pole
  conpos = [xll,yll+poledeg,xur,yur-poledeg]
; conpos = [xll,yll+poledeg,xur,yur-poledeg*2.5]
endif else begin
  conpos = [xll,yll,xur,yur]
endelse
xmid = .5*(xll+xur)
;
; Set up map projection:
;
; If censlt >= 0, then it is user-specified local time to fix at center
; of projection.
;
uthrs = float(map.mtime[1])+float(map.mtime[2])/60.
censlt = map.censlt ; local time (hrs) of center of lon axis of projection
;
plat=0 & plon=map.cenlon & rot = 0.
if(map.fixpolar eq 'LON') then censlt=-1 
if censlt ge 0. then begin
  mag = 0 & if field.grid_type eq 'magnetic' then mag = 1
  plon = fslt(censlt,uthrs,plon,'getlon',mag=mag)
endif
map_set,plat,plon,rot,/cylindrical,position=mappos,/noborder
;
; Titles, labels, tic marks, etc:
;
if field.grid_type eq 'geographic' then begin
  xtitle = 'LONGITUDE (DEGREES GEOGRAPHIC)'
  if censlt ge 0. then xtitle = 'SOLAR LOCAL TIME (HOURS)' 
  ytitle = 'LATITUDE (DEGREES GEOGRAPHIC)'
endif else begin
  xtitle = 'LONGITUDE (DEGREES MAGNETIC)'
  if censlt ge 0. then xtitle = 'MAG LOCAL TIME (HOURS)' 
  ytitle = 'LATITUDE (DEGREES MAGNETIC)'
endelse
xticlabs = ['-180','-120','-60','0','60','120','180']
if censlt ge 0. then begin
  sltlabs = fltarr(7)
  sltlabs[3] = censlt
  for i=2,0,-1 do begin
    sltlabs[i] = sltlabs[i+1]-4.
    if(sltlabs[i] ge 24.) then sltlabs[i]=sltlabs[i]-24
    if(sltlabs[i] lt 0.) then sltlabs[i]=sltlabs[i]+24
  endfor
  for i=4,6    do begin
    sltlabs[i] = sltlabs[i-1]+4.
    if(sltlabs[i] ge 24.) then sltlabs[i]=sltlabs[i]-24
    if(sltlabs[i] lt 0.) then sltlabs[i]=sltlabs[i]+24
  endfor
  xticlabs = strarr(7)
  for i=0,6 do begin
	xticlabs[i] = string(format="(f5.1)",sltlabs[i])
	xticlabs[i] = strcompress(xticlabs[i],/remove_all)
  endfor
; print,'pltce: censlt=',censlt,' xticlabs=' & print,xticlabs
endif
if(map.fixpolar eq 'LON') then begin
  lonlabs=fltarr(7)
  lonlabs[3]=map.cenlon
  for i=2,0,-1 do begin 
    lonlabs[i] = lonlabs[i+1]-60.
    if(lonlabs[i] gt 180.) then lonlabs[i]=lonlabs[i]-360
    if(lonlabs[i] lt -180.) then lonlabs[i]=lonlabs[i]+360
  endfor
  for i=4,6    do begin
    lonlabs[i] = lonlabs[i-1]+60.
    if(lonlabs[i] gt 180.) then lonlabs[i]=lonlabs[i]-360
    if(lonlabs[i] lt -180.) then lonlabs[i]=lonlabs[i]+360
  endfor
  xticlabs = strarr(7)
  for i=0,6 do begin
	xticlabs[i] = string(format="(f6.1)",lonlabs[i])
	xticlabs[i] = strcompress(xticlabs[i],/remove_all)
  endfor
endif

yticlabs = ['-90','-60','-30','0','30','60','90']
;
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 map.fmin=fmin & map.fmax=fmax
endif else begin
  fmin = map.fmin & fmax = map.fmax
endelse
;
; Byte-scale to an image. Make any missing data black:
;
if (map.plottype eq 'image_only') or (map.plottype eq 'image+contours') then begin
  missing_indices = where(data eq map.missing_value,nmissing)
  if nmissing gt 0 then begin
    print,'Field ',field.name,' CE map: Found ',nmissing,' missing data points.'
    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
;
; Scale image up and warp image to map projection:
;
 image = map_patch(image,lons,lats,xstart=startx,ystart=starty,$
   xsize=xsize,ysize=ysize) ; dashed black line bug
;
; Set up color bar (pro clrbar):
;
  barw = .90*(xur-xll) & barh = .03
  xbar = xmid-.5*barw & ybar = yll-barh-.1
  cblabels = [strcompress(string(fmin,format="(g10.4)"),/remove_all),$
    strcompress(string(0.5*(fmin+fmax),format="(g10.4)"),/remove_all),$
    strcompress(string(fmax,format="(g10.4)"),/remove_all)]
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 label):
;
uthrs = float(map.mtime[1])+float(map.mtime[2])/60.
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 or geometric altitude plotting
;
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')
;
; 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 label:
;
 botlab1 = 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
;
;Rotate data in longitude if necessary (i.e., fixed slt has been requested):
;
if(map.fixpolar eq 'SLT') then begin
  ix0 = ixfind_nearest(lons,0.)
  ixcen = ixfind_nearest(lons,plon)
  if ixcen lt ix0 then ishift = ix0-ixcen else ishift = -(ixcen-ix0)
  data = shift(data,ishift,0)

; print,format=$
;   "('fix=',a,' ut=',f6.2,' censlt=',f8.2,' plon=',f8.2,' ix0=',i3,' ixcen=',i3,' ishift=',i3)",$
;   map.fixpolar,uthrs,censlt,plon,ix0,ixcen,ishift

endif
if(map.fixpolar eq 'LON') then begin
  ix0 = ixfind_nearest(lons,0.)
  ixcen = ixfind_nearest(lons,plon)
  if ixcen lt ix0 then ishift = ix0-ixcen else ishift = -(ixcen-ix0)
  data = shift(data,ishift,0)
endif
;
; get wind vector data
;
if vector.control ne 'off' then begin
; zonal data
  zdata = get_field_data(info,vector.ilonvec) 
  zdata = shift(zdata,ishift,0)
  *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) 
  mdata = shift(mdata,ishift,0)
  *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.
; Use position=mappos when drawing axes and labels, and use
;   position=conpos when drawing contours.
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
      addmapgrid,map,field
      contour,data,lons,lats,/noerase,/nodata,position=mappos,$
        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
      clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize
      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=xmid,yscale=yll-.3, $
	  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=xmid,yscale=yll-.3, $
	  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+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab   ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.24,botlab1  ,/norm,align=0.5,charsize=chsize
  end
;
; Display image with contour overlay:
;
  'image+contours': begin
    if vector.control ne 'only' then begin
      tv,image,startx,starty,xsize=xsize,ysize=ysize
      addmapgrid,map,field
;   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
;   contours:
      contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
        /follow,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=xmid,yscale=yll-.3, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
	  vector.vmag=vmagnitude
      endif
      clrbar,barw,barh,xbar,ybar,cblabels,charsize=chsize
    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=xmid,yscale=yll-.3, $
	  color=vector.carrows, $
	  xarrowlabel=vector.lonvector, yarrowlabel=vector.latvector)
        if vector.vecanimate eq 'on' then $
	  vector.vmag=vmagnitude
    endelse
      xyouts,xmid,yur+.06,fieldlab ,/norm,align=0.5,charsize=chsize
      xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
      xyouts,xmid,yll-.21,minmaxlab,/norm,align=0.5,charsize=chsize
      xyouts,xmid,yll-.25,botlab1,/norm,align=0.5,charsize=chsize
  end
;
; Contour only (B&W):
;
  'monochrome_contours': begin
    if vector.control ne 'only' then begin
      addmapgrid,map,field
      contour,data,lons,lats,/noerase,/nodata,position=mappos,$
        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
      contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
        /follow,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=xmid,yscale=yll-.3, $
	  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=xmid,yscale=yll-.3, $
	  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+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.11,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.15,botlab1,/norm,align=0.5,charsize=chsize
  end
;
; Color fill contour:
;
  'colorfill_contours': begin
    if vector.control ne 'only' then begin
      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
      contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
        c_charsize=c_charsize,levels=clevels,/follow,/fill
      contour,data,lons,lats,/noerase,position=conpos,xstyle=5,ystyle=5,$
        c_charsize=c_charsize,levels=clevels,/follow,c_linestyle=clinestyle,$
        c_colors=[map.clineclr]
      map_set,/cylindrical,position=mappos,/noborder,/noerase
      addmapgrid,map,field
      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=xmid,yscale=yll-.3, $
	  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=xmid,yscale=yll-.3, $
	  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+.06,fieldlab ,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yur+.03,toplab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.11,minmaxlab,/norm,align=0.5,charsize=chsize
    xyouts,xmid,yll-.15,botlab1,/norm,align=0.5,charsize=chsize
  end
  else: print,'>>> Global CE maps: unknown plot type ',map.plottype
endcase
;
; map.fmin,max and map.continents may have been changed, 
; so update state:
;
map.vector=vector
*info.pmap = map
end