;+
;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
;
;See "conv_units" to change units.
;
; $LastChangedBy: dmitchell $
; $LastChangedDate: 2022-12-26 12:23:43 -0800 (Mon, 26 Dec 2022) $
; $LastChangedRevision: 31379 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_6_0/projects/maven/sta/mvn_sta_functions/contour4d.pro $
;
;CREATED BY:	J. McFadden  14-02-07
;-
pro contour4d,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, $
	tof = tof

if size(/type,tempdat) ne 8 or tempdat.valid eq 0 then begin
    print,'Invalid Data'
    return
endif

; Custom plot margins (temporary fix?)

xomar = !x.omargin
yomar = !y.omargin
!x.omargin = [2,3]
!y.omargin = [0,5]

; Get data

mdat = omni4d(tempdat)
mdat = conv_units(mdat,units)
nenergy = mdat.nenergy
nmass = mdat.nmass

b_ns  = 5.844 					; TOF bins/ns for TOF timing (not sure if this is 5.844 or 5.855 for flight unit) 

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
	if keyword_set(tof) then begin
		ydat=mdat.tof_arr/b_ns 
        	if not keyword_set(ytitle) then ytitle = 'TOF (ns)'
		if not keyword_set(ylin) then ylog=1 else ylog=0
		if ylog then yrange=[3,200] else yrange=[0,200]
	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
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
    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
    str_element,limits,'zticks',index=index
    if index ge 0 then zticks=limits.zticks else zticks=4
    str_element,limits,'charsize',index=index
    if index ge 0 then csize=limits.charsize

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
	if not keyword_set(xmargin) then xmargin=[10,10]
	if not keyword_set(ymargin) then ymargin=[4,2]
    if keyword_set(rotate) then begin
        zdat=transpose(zdat)
        xdat=transpose(xdat)
        ydat=transpose(ydat)
        contour, zdat, ydat, xdat, title=title, xtitle=ytitle, ytitle=xtitle, $
          c_colors=c_colors, levels=levels, yrange=xrange, ystyle=xstyle, $
          ylog=xlog, xrange=yrange, xstyle=ystyle, xlog=ylog, charsize=csize, $
	  fill=fill, xmargin=xmargin, ymargin=ymargin, xticks=yticks, xtickv=ytickv
    endif else begin
        contour, zdat, 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, charsize=csize, $
          fill=fill, xmargin=xmargin, ymargin=ymargin, yticks=yticks, ytickv=ytickv
    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], charsize=csize, $
	ystyle=1, yticks=zticks, ytitle=ztitle, $
	yticklen = 0.02, /norm, ytickformat='(F5.1)',ylog=0

endif

; Restore previous plot margins (temporary fix?)

!x.omargin = xomar
!y.omargin = yomar

return

end