contour4d.pro (Documentation for /home/davin/idl/socware/)

;+

;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.

;

;

;CREATED BY:	J. McFadden  14-02-07

;FILE:  contour4d.pro

;VERSION 1.

;MODIFICATIONS: 

;		

;-

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



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,/mass)



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

    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.)

	zmin=10.^fix(alog10(min(zdat(ind)))+.1)

	zmax=10.^fix(alog10(max(zdat))+.9)

	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,'CONTOUR2D: 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, $

	  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, $

          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], $

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

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

    

endif



!y.omargin = [0,0]



return

end