;+
;.r /Users/cmfowler/IDL/STATIC_routines/Generic/contour4d_edit.pro
;Original routine from Jim McFadden. CMF edited some keywords, and now uses it in plot_c6, etc.
;
;PROCEDURE: contour4d,data
;PURPOSE:
; Produces contour plots of energy-mass distributions from 4D data structures.
;INPUTS:
;   data   - structure containing 4d data  (obtained from get_mvn_?() routine)
;   e.g. "get_mvn_c6, get_mvn_ce, etc."
;KEYWORDS:
; LIMITS - A structure containing limits and display options.
;             see: "options", "xlim" and "ylim", to change limits
; UNITS  - convert to given data units before plotting
; TITLE - Title to be plotted,
;   - set title=' ' for NO Title!
;   - set title='1' for just the time in the title
; XTITLE - xtitle to be plotted,
;   - set xtitle=' ' for NO xtitle!, default determined by VEL keyword
; YTITLE - ytitle to be plotted,
;   - set ytitle=' ' for NO ytitle!, default=data.units_name
; ZTITLE - ztitle to be plotted,
;   - set ztitle=' ' for NO ytitle!, default='Log!D10!N('+y_units+')'
; RETRACE - set to number of retrace steps removed,
;   - default set to 0
; VEL - If set, x-axis is velocity km/s  -- Default is Energy (eV)
;
; NCONT - Number of contours to be plotted, default = 8
; LEVELS - Explicit contour levels, default levels spaced down
;          from max by 10^.5
; FILL - If set, contours are filled with solid color or gray scale
; BW - If set, contours are white, no affect on fill plots
; PURE - If set, 6 pure colors are cycled through for levels
; ROTATE - Exchanges x and y axes for non-polar plots
; LABEL - Labels the contour levels
;       XMARGIN - Change xmargin from default
; YMARGIN - Change ymargin from default
; POINTS - adds data points to plot
;top_c, bottom_c: top and bottom indices to use in the colorbar. Default if not set is bottom=0, top=254.
;See "conv_units" to change units.
;
;
;CREATED BY:  J. McFadden  14-02-07
;FILE:  contour4d.pro
;VERSION 1.
;MODIFICATIONS:
;
;-

pro contour4d_edit,tempdat,   $

  LIMITS = limits, $

  UNITS = units,   $

  TITLE = title, $

  YTITLE = ytitle, $

  XTITLE = xtitle, $

  ZTITLE = ztitle, $

  RETRACE = retrace,   $

  VEL = vel, $

  NCONT = ncont, $

  LEVELS = levels, $

  FILL = fill, $

  BW = bw, $

  PURE = pure, $

  ROTATE = rotate, $

  LABEL = label, $

  XMARGIN = xmargin, $

  YMARGIN = ymargin, $

  POINTS = points, $
  mass = mass, $
  zrange = zrange, $
  ylin = ylin, $
  xlin = xlin, $
  twt = twt, charsize=charsize, top_c=top_c, bottom_c = bottom_c, position=position, charthick=charthick, thick=thick, xrange=xrange

  ;$%$%$ CMF : specify colorbar indice range:
  if size(bottom_c, /type) eq 0 then bottom_c = 0  ;default is full colorbar
  if size(top_c, /type) eq 0 then top_c =254

  if size(/type,tempdat) ne 8 or tempdat.valid eq 0 then begin

    print,'Invalid Data'

    return

  endif



  !y.omargin =[2,3]               ; temporary fix
  !x.omargin =[0,5]               ; temporary fix

  mdat = omni4d(tempdat)



  mdat = conv_units(mdat,units)

  nenergy = mdat.nenergy

  nmass = mdat.nmass



  if not keyword_set(title) then begin

    title = mdat.project_name+'  '+mdat.data_name+' ' + $

      mdat.units_name

    ;    title = title + '!C'+time_string((mdat.time+mdat.end_time)/2.)

    title = title + '!C'+time_string(mdat.time)+'-'+strmid(time_string(mdat.end_time),11,8)

  endif else if title eq '1' then begin

    title = '' + time_string((mdat.time+mdat.end_time)/2.)

  endif


  if keyword_set(mass) then begin
    ydat=mdat.mass_arr
    if not keyword_set(ytitle) then ytitle = 'Mass/q (amu/q)'
    if not keyword_set(ylin) then ylog=1 else ylog=0

    yrange=[.5,max(ydat)+1]

  endif else begin
    ydat=replicate(1.,mdat.nenergy)#findgen(mdat.nmass)
    if not keyword_set(ytitle) then ytitle = 'Mass Bin'
    ylog=0

    yrange=[0,max(ydat)+1]

  endelse
  ystyle=1

  yticks = 0

  ytickv = 0


  zdat = mdat.data

  if keyword_set(twt) then zdat = zdat/mdat.twt_arr



  if keyword_set(vel) then begin

    ; str_element,limits,'velocity',value=vel

    xdat = velocity(mdat.energy,mdat.mass*mdat.mass_arr)

    if not keyword_set(xtitle) then xtitle = 'Velocity (km/s)'

  endif else begin

    xdat = mdat.energy

    if not keyword_set(xtitle) then xtitle = 'Energy/q  (eV/q)'

  endelse



  ;if keyword_set(retrace) then begin

  ;       if retrace gt 0 then zdat(0:retrace-1,*)=0. $

  ;                else zdat(nenergy+retrace:nenergy-1,*)=0.

  ;endif



  ; Set plot limits



  xmin=.9*min(xdat(*,0))-1. > .1

  xmax=1.1*max(xdat(*,0))
  if keyword_set(vel) then xmax = xmax > 100.

  if keyword_set(vel) then xmin = 1
  if size(xrange,/type) eq 0 then xrange=[xmin,xmax]

  xstyle=1

  if not keyword_set(xlin) then xlog=1 else xlog=0

  if not keyword_set(zrange) then begin

    ind=where(zdat ne 0.,zcount)

    if zcount gt 1 then zmin=10.^fix(alog10(min(zdat(ind)))+.1) else zmin=1

    if zcount gt 1 then zmax=10.^fix(alog10(max(zdat))+.9) > 10.*zmin else zmax=10

    zrange=[zmin,zmax]
  endif else begin
    zmax=max(zrange)
    Zmin=min(zrange)
  endelse

  zlog=1


  if keyword_set(limits) then begin



    str_element,limits,'xrange',index=index

    if index ge 0 then xrange=limits.xrange

    str_element,limits,'xstyle',index=index

    if index ge 0 then xstyle=limits.xstyle

    str_element,limits,'xlog',index=index

    if index ge 0 then xlog=limits.xlog



    str_element,limits,'yrange',index=index

    if index ge 0 then begin

      yrange=limits.yrange

      yticks=0

      ytickv=0

    endif

    str_element,limits,'ystyle',index=index

    if index ge 0 then ystyle=limits.ystyle

    str_element,limits,'ylog',index=index

    if index ge 0 then ylog=limits.ylog



    str_element,limits,'zrange',index=index

    if index ge 0 then begin

      zrange=limits.zrange

      if n_elements(limits.zrange) GT 1 then begin

        zmin = limits.zrange(0)

        zmax = limits.zrange(1)

      endif

    endif

    str_element,limits,'zlog',index=index

    if index ge 0 then zlog=limits.zlog



  endif


  if not keyword_set(units) then str_element,limits,'units',value=units



  if not keyword_set(levels) then begin

    if keyword_set(ncont) then begin

      if (ncont GT 29 or ncont LT 1) then begin

        print,'contour4d: Keyword NCONT must be in the range 1 to 29.'

        print,'Setting NCONT = 25'

        ncont=25

      endif

    endif

    if not keyword_set(ncont) then ncont=25
    if zlog eq 1 then begin

      levels=zrange(1)*10.^(-(findgen(ncont)/(ncont-1) * (alog10(zmax) - $

        alog10(zmin))))

      levels=reverse(levels)

    endif else begin

      levels=(zrange(1)-zrange(0))*findgen(ncont)/ncont + zrange(0)


    endelse
  endif else ncont = dimen1(levels)



  if keyword_set(bw) and not keyword_set(fill) then c_colors=[1000] else $

    begin

    if keyword_set(pure) then c_colors=bytescale(pure=ncont) else $

      c_colors=bytescale(findgen(ncont))

  endelse

  ;Remove colors that lie below lower limit; set ones above to max:
  zdat2 = zdat
  if size(zrange,/type) ne 0 then begin
    iRM = where(zdat2 eq 0, niRM)
    if niRM eq 0 then zdat2[iRM] = !values.f_nan
    ;   iRM = where(zdat2 lt zrange[0], niRM)
    ;   if niRM gt 0 then zdat2[iRM] = !values.f_nan ;zrange[0]

    ;   iRM = where(zdat2 gt zrange[1], niRM)
    ;   if niRM gt 0 then zdat2[iRM] = !values.f_nan ;zrange[1]
  endif

  ;$%$%$ CMF - overwrite colors here to take into account qualcolors:
  c_colors = mvn_sta_bytscl2(findgen(ncont), bottom=bottom_c, top=top_c)


  if not keyword_set(xmargin) then xmargin=[10,10]

  if not keyword_set(ymargin) then ymargin=[4,2]

  if keyword_set(rotate) then begin

    zdat2=transpose(zdat2)   ;NOTE zdat2 here after dealing with min/max above

    xdat=transpose(xdat)

    ydat=transpose(ydat)

    contour, zdat2, ydat, xdat, title=title, xtitle=ytitle, ytitle=xtitle, $

      c_colors=c_colors, levels=levels, yrange=xrange, ystyle=xstyle, $

      ylog=xlog, xrange=xrange, xstyle=ystyle, xlog=ylog, $

      fill=fill, xmargin=xmargin, ymargin=ymargin, xticks=yticks, xtickv=ytickv, charsize=charsize, position=position, $
      charthick=charthick, thick=thick

  endif else begin
    contour, zdat2, xdat, ydat, title=title, xtitle=xtitle, ytitle=ytitle, $

      c_colors=c_colors, levels=levels, yrange=yrange, ystyle=ystyle, $

      ylog=ylog, xrange=xrange, xstyle=xstyle, xlog=xlog, $

      fill=fill, xmargin=xmargin, ymargin=ymargin, yticks=yticks, ytickv=ytickv, charsize=charsize, position=position, $
      charthick=charthick, thick=thick

  endelse

  if keyword_set(points) then begin

    if keyword_set(xrange) and keyword_set(yrange) then begin

      ind=where(xdat gt xrange(0) and xdat lt xrange(1) and ydat gt yrange(0) and ydat lt yrange(1))

      xyouts,xdat(ind),ydat(ind),'.'
    endif else if keyword_set(xrange) then begin

      ind=where(xdat gt xrange(0) and xdat lt xrange(1))

      xyouts,xdat(ind),ydat(ind),'.'

    endif else if keyword_set(yrange) then begin

      ind=where(ydat gt yrange(0) and ydat lt yrange(1))

      xyouts,xdat(ind),ydat(ind),'.'

    endif else xyouts,xdat,ydat,'.'

  endif



  time_stamp



  if keyword_set(label) then begin



    bar_levels = [[c_colors],[c_colors]]



    bar_width = 0.03*(!x.window(1)-!x.window(0))

    bar_start = 0.01

    bar_pos1_norm = [!x.window(1)+bar_start, !y.window(0)]

    bar_height = !y.window(1) - !y.window(0)

    bar_pos1_dev = convert_coord(bar_pos1_norm(0), bar_pos1_norm(1), $

      /norm, /to_device)

    bar_pos2_dev = convert_coord(bar_pos1_norm(0) + bar_width, $

      bar_pos1_norm(1) + bar_height, /norm, $

      /to_device)

    x_size = -(bar_pos1_dev(0) - bar_pos2_dev(0))

    y_size = -(bar_pos1_dev(1) - bar_pos2_dev(1))

    bar_pos2_norm = convert_coord(bar_pos1_dev(0)+x_size, $

      bar_pos1_dev(1)+y_size, /device, $

      /to_norm)

    bar_levels = transpose(congrid(bar_levels,y_size,x_size))



    if !d.name EQ 'PS' then bar_levels = $

      transpose(congrid([[c_colors],[c_colors]],200,50))



    tv, bar_levels, bar_pos1_norm(0), bar_pos1_norm(1), $

      xsize = bar_pos2_norm(0) - bar_pos1_norm(0), $

      ysize = bar_pos2_norm(1) - bar_pos1_norm(1), /norm



    bar_min = alog10(levels(0))

    bar_max = alog10(levels(ncont-1))



    y_units = strlowcase(units)

    if keyword_set(twt) then begin

      case y_units of

        'counts':y_units = 'counts/tofbin'

        'eflux':y_units = 'eV/(cm!E2!N-s-sr-eV-tofbin)'

        'rate':y_units = '#/sec-tofbin'

        'crate':y_units = 'Deadtime Corrected #/sec-tofbin'

        'flux':y_units = '#/(cm!E2!N-s-sr-eV-tofbin)'

        'df':y_units = '#/(cm!E3!N-(km/sec)!E3!N/tofbin'

        else:y_units = units

      endcase

      if not keyword_set(ztitle) then ztitle='Log!D10!N('+y_units+')'

      if label eq 3 then ztitle=' '
    endif else begin
      case y_units of

        'counts':y_units = 'counts'

        'eflux':y_units = 'eV/(cm!E2!N-s-sr-eV)'

        'rate':y_units = '#/sec'

        'crate':y_units = 'Deadtime Corrected #/sec'

        'flux':y_units = '#/(cm!E2!N-s-sr-eV)'

        'df':y_units = '#/(cm!E3!N-(km/sec)!E3!N)'

        else:y_units = units

      endcase

      if not keyword_set(ztitle) then ztitle='Log!D10!N('+y_units+')'

      if label eq 3 then ztitle=' '
    endelse


    axis, bar_pos2_norm(0),0,0, yaxis=1, $
      yrange=[bar_min,bar_max], $

      ystyle=1, yticks=4, ytitle=ztitle, $

      yticklen = 0.02, /norm, ytickformat='(F5.1)',ylog=0, charsize=charsize, charthick=charthick



  endif



  !y.omargin = [0,0]


  return

end