;+
;NAME:			slice2d
;PURPOSE:		creates a velocity or energy spectrogram
;			with v perp and v para as x and y, and the
;			specified units as z (color axis).
;			DOESN'T fold data points; generates true 2D slice
;			through the B-V plane
;CALL:			ex: slice2d,get_el('20:31'),[keywords]
;KEYWORDS:	ANGLE: the two angle limits of the wedge selected for plot (default [-20,20]).
;           THIRDDIRLIM: the two velocity limits of the wedge. If selected, the ANGLE keyword would be invalid.
;			FINISHED: makes the output publication quality when using ps.
;			ZHEIGHT: the zrange to use.  Def is 100 km/s
;			XRANGE: vector specifying the xrange
;			RANGE: vector specifying the color range
;			UNITS: specifies the units ('eflux','df',etc.) (Def. is 'df')
;			NOZLOG: specifies a linear Z axis
;			THEBDATA: specifies b data to use (def is B3_gse)
;			VAR_LA: vector of tplot variables to show on plot
;			POSITION: positions the plot using a 4-vector
;			ERANGE: specifies the energy range to be used
;			NOFILL: doesn't fill the contour plot with colors
;			NLINES: says how many lines to use if using NOFILL (def 60, max 60)
;			NOOLINES: suppresses the black contour lines
;			NUMOLINES: how many black contour lines (def. 20, max 60)
;			SHOWDATA: plots all the data points over the contour
;			PLOTENERGY: plots using energy instead of velocity
;			VEL: tplot variable containing the velocity data
;			     (default is calculated with v_3d)
;			NOGRID: forces no triangulation
;			NOCROSS: suppresses cross section line plots
;			RESOLUTION: resolution of the mesh (default is 51)
;			NOSMOOTH: suppresses smoothing
;			NOSUN: suppresses the sun direction line
;			RMBINS: removes the sun noise by cutting out certain bins
;			THETA: specifies the theta range for RMBINS (def 20)
;			PHI: specifies the phi range for RMBINS (def 40)
;			NR: removes background noise from ph using noise_remove
;			NOISELEVEL: background level in eflux
;			BOTTOM: level to set as min eflux for background. def. is 0.
;			SR, RS, RM2: removes the sun noise using subtraction
;				REQUIRES write_ph.doc to run
;				Note: automatically sets /nosmooth
;				for smoothing, set /smooth
;			NLOW: used with rm2.  Sets bottom of eflux noise level
;				def. 1e4
;			M: marks the tplot at the current time
;			NOVELLINE: suppresses the velocity line
;			LOGLINES: makes the cross sections plot with log
;			CUT_PARA: the value of vpara to make the 1d cut of vperp
;			CUT_PERP: the value of vperp to make the 1d cut of vpara
;			VEL2: takes a 3-vector velocity and puts it on the plot
;LAST MODIFIED:		1-26-99
;CREATED BY:		Arjun Raj
;EXAMPLES:
;			slice2d,get_el('21:00')
;			displays 2d cut using el
;			slice2d,get_el('21:00'),/nofill,nlines = 30
;			displays 2d cut without filling in the plot with color, using 30 lines
;			slice2d,get_el('23:00'),numolines = 40
;			displays 2d cut with 40 black overplot lines
;			slice2d,get_ph('21:00'),thebdata = 'Bexp'
;			displays 2d cut with Bexp used as b data
;			slice2d,get_ph(/adv), units = 'eflux'
;			advances time one step, uses eflux units
;			slice2d,get_el(/adv),erange = [25,900]
;			restricts energies to cut out photoelectrons
;			slice2d,get_el(/adv),vel = 'Vconv',/subtract
;			transforms into velocity frame, using specified velocity
;			slice2d,get_ph(/adv),/plotenergy,/loglines,/subtract
;			transforms into bulk frame, plots versus energy, and uses a log scale
;REMARKS:		when calling with phb and rm2, use file='write_phb.doc'
;			also, set the noiselevel to 1e5.  This gives the best
;			results
;-


;Yes, I know the program is a mess...
;Sorry.

pro slice2d_themis_longer_esa,sc,typ,current_time,timeinterval,nofix = nofix,nsteps = nsteps,xrange = xrange,$
	range = range, units = units,nozlog = nozlog,zlog = zlog,thebdata = thebdata,$
	b3 = b3,position = position,erange = erange,nofill = nofill,var_label = var_label,$
	nlines = nlines,showdata = showdata,plotenergy = plotenergy,velocity=vel,nosubtract = nosubtract,$
	nosmooth=nosmooth, smooth=smooth,nocross = nocross,cross = cross,nogrid = nogrid,grid = grid,$
	resolution = resolution, pos2 = pos2,nosun = nosun,sundir = sundir,olines = olines,$
	noolines = noolines,numolines = numolines,setpos = setpos,leavezero = leavezero,$
	rmbins = rmbins,nr = nr,noiselevel = noiselevel,bottom = bottom,zheight = zheight,$
	theta = theta,phi = phi,m=m,rm2=rm2,nlow = nlow,angle = angle,phb = phb,filename = filename,$
	novelline = novelline,subtract = subtract,double = double,outfile = outfile,cut_perp = cut_perp,$
	cut_para=cut_para,rs = rs,outcuts = outcuts,sr = sr,loglines = loglines,oldlog = oldlog,$
	logplot = logplot,finished = finished,plotlabel = plotlabel,vel2 = vel2,cut_bulk_vel=cut_bulk_vel,_EXTRA = e,$
	rotation=rotation,filetype=filetype,outputfilename=outputfile,ThirdDirlim=ThirdDirlim,species=species


!p.charsize=1

if not keyword_set(nocross) then begin
    if rotation eq 'BV' then begin
        cross = 0
    endif else begin
        cross =0
    endelse
endif else begin
    cross = 0
endelse
if keyword_set(phb) then filename = 'write_phb.doc'

for i=0,timeinterval/3-1 do begin
  if species eq 'ion' then begin
    thedata2=CALL_FUNCTION('get_th'+sc+'_pei'+typ,current_time+i*3.0)
;    thedata2s=CALL_FUNCTION('thm_sst_psif',current_time+i*3.0,probe=sc)
  endif
  if species eq 'ele' then begin
    thedata2=CALL_FUNCTION('get_th'+sc+'_pee'+typ,current_time+i*3.0)
;    thedata2s=CALL_FUNCTION('thm_sst_psef',current_time+i*3.0,probe=sc)
  endif

  if i eq 0 then begin
    thedata3 = thedata2
;    thedata3s = thedata2s
  endif else begin
    thedata3 = [thedata3,thedata2]
;    thedata3s = [thedata3s,thedata2s]
  endelse

endfor
inumber=i-1

; %%% LOOP START HERE  %%%

for in=0,inumber do begin
thedata=thedata3(in)

bins_2d=fltarr(thedata.nenergy,thedata.nbins)
for i=0,thedata.nbins-1 do begin
;    bins_2d(*,i)=thedata.bins(i)
    bins_2d(*,i)=thedata.bins(*,i)
endfor

if keyword_set(rmbins) then begin
	print, 'Removing bins (bin_remove)'
	thedata = bin_remove(thedata,theta = theta,phi = phi)
endif ;else thedata = thedata2

if keyword_set(sr) then rm2 = 1
if keyword_set(rs) then rm2 = 1
if keyword_set(nofill) then noolines = 1

if keyword_set(rm2) then begin
	print, 'Removing bins (bin_remove2)'
	leavezero = 1
	if not keyword_set(smooth) then nosmooth = 1
	;nr = 1
	load_ph,new,filename = filename
	thedata = bin_remove2(thedata,theta = theta,phi = phi,new= new,nlow = nlow)
endif ;else thedata = thedata2


if keyword_set(nr) then begin
	print,'Removing Noise'
;print, noiselevel
	thedata = noise_remove(thedata,nlevel = noiselevel,bottom = bottom)
	leavezero = 1
endif

;Tai's addition: Always have the leavezero=1 option.
	leavezero = 1

if keyword_set(m) then $
	new_time,'cut2d',thedata.time

numperrow=4

;MODIFICATIONS TO MAKE COMMAND LINE SMALLER

if not keyword_set(zheight) and not keyword_set(angle) then angle = [-1,1]

if not keyword_set(nozlog) then zlog = 1
if not keyword_set(nogrid) then grid = 1

if not keyword_set(nosmooth) then smooth = 1
if not keyword_set(noolines) then begin
	if keyword_set(numolines) then olines = numolines else olines = 20
	endif
if not keyword_set(thebdata) then thebdata = 'B3_gse'
if keyword_set(b3) then thebdata = 'B3_gse'
print, 'B field used is '+thebdata
if not keyword_set(subtract) then nosubtract = 1
if not keyword_set(nosun) then sundir = 1

if keyword_set(zlog) then print,'zl'
if keyword_set(grid) then print,'grid'
if keyword_set(cross) then print,'cross'
if keyword_set(smooth) then print,'smooth'
if keyword_set(olines) then print,'olines'

if not keyword_set(nsteps) then nsteps = 18
if not keyword_set(units) then units = 'df'
if not keyword_set(nlines) then nlines = 60

if not keyword_set(rotation) then rotation='BV'


if filetype eq 'ps' then begin
    SET_PLOT, 'PS'
    DEVICE, FILENAME=outputfile+'.ps',/color,bits_per_pixel=8
endif

;END MODIFICATiONS


perpsym = byte(94)
perpsymbol = string(perpsym)

parasym = byte(47)
parasymbol = string(parasym) + string(parasym)



if keyword_set(finished) and units eq 'df' then thedata.data = thedata.data / 1000.  ;changing units to
									;s^-3 m^-6

if keyword_set(finished) and keyword_set(plotlabel) and !d.name eq 'PS' then begin
;	device,/bold
;	xyouts, 0.0,.95,plotlabel+'!N!7',/normal
endif


if !d.name eq 'PS' then loadct,39

if not keyword_set(resolution) then resolution = 51
if resolution mod 2 eq 0 then resolution = resolution + 1

oldplot = !p.multi

if keyword_set(cross) then begin ;and  !d.name ne 'PS' then begin
	!p.multi = [0,2,1]
	grid = 1
endif

;if not keyword_set(vel) then vel = 'v_3d_ph'

if not keyword_set(position) then begin
	x_size = !d.x_size & y_size = !d.y_size
	xsize = .77
	yoffset = 0.
	d=1.
	if keyword_set(cross) then begin
		yoffset = yoffset + .5
		xsize = xsize/2.+.13/1.5
		y_size = y_size/2.
		x_size = x_size/2.
		d = .5
		if y_size le x_size then $
			pos2 = [.13*d+.05,.03+.13*d,.05+.13*d + xsize * y_size/x_size,.13*d + xsize+.03] else $
			pos2 = [.13*d+.05,.03+.13*d,.05+.13*d + xsize,.13*d + xsize *x_size/y_size+.03]

	endif
	if y_size le x_size then $
		position = [.13*d+.05,.13*d+yoffset,.05+.13*d + xsize * y_size/x_size,.13*d + xsize + yoffset] else $
		position = [.13*d+.05,.13*d+yoffset,.05+.13*d + xsize,.13*d + xsize *x_size/y_size + yoffset]
endif else begin
	if not keyword_set(pos2) then begin
		pos2 = position
		pos2(0) = position(0)
		pos2(2) = position(2)
		pos2(3) = position(1)-.08
		pos2(1) = .1
	endif
endelse


if keyword_set(var_label) and not keyword_set(setpos) then begin; and !d.name ne 'PS' then begin
	if keyword_set(cross) then begin
		pos2(1) = pos2(1) + .04
		pos2(3) = pos2(3) + .04
	endif
	position(1) = position(1) + .02
	position(3) = position(3) + .02
endif


;theonecnt = thedata
thedata = conv_units(thedata,units)


;stop

;thedata.data(*,120)=0


;theonecnt = conv_units(theonecnt,'counts')
;for i = 0,theonecnt.nenergy-1 do theonecnt.data(i,*) = 1
;theonecnt = conv_units(theonecnt,units)
;if theonecnt.units_name eq 'Counts' then theonecnt.data(*,*) = 1.

;**********************************************
;bad_bins=where((thedata.dphi eq 0) or (thedata.dtheta eq 0) or $
;	((thedata.data(0,*) eq 0.) and (thedata.theta(0,*) eq 0.) and $
;	(thedata.phi(0,*) eq 180.)),n_bad)
;good_bins=where(((thedata.dphi ne 0) and (thedata.dtheta ne 0)) and not $
;	((thedata.data(0,*) eq 0.) and (thedata.theta(0,*) eq 0.) and $
;	(thedata.phi(0,*) eq 180.)),n_good)

;if n_bad ne 0 then print,'There are bad bins'


if thedata.valid ne 1 then begin
	print,'Not valid data'
	return
endif

;bad120 = where(good_bins eq 120,count)
;if count eq 1 and thedata.data_name eq 'Pesa High' then begin
;	print, 'Fixing bad 120 bin'
;	if n_bad eq 0 then bad_bins = [120] else bad_bins = [bad_bins,120]
;	good_bins = good_bins(where(good_bins ne 120))
;	n_bad = n_bad + 1
;	n_good = n_good -1
;endif



;*****************************

;In order to find out how many particles there are at all the different locations,
;we must transform the data into cartesian coordinates.


totalx = fltarr(1) & totaly = fltarr(1) & totalz = fltarr(1)
ncounts = fltarr(1)

if not keyword_set(erange) then begin
	erange = [thedata.energy(thedata.nenergy-1,0),thedata.energy(0,0)]
	erange = [min(thedata.energy), max(thedata.energy)]
	eindex = indgen(thedata.nenergy)
endif else begin
	eindex = where(thedata.energy(*,0) ge erange(0) and thedata.energy(*,0) le erange(1))
	erange = [min(thedata.energy(eindex,0)),max(thedata.energy(eindex,0))]
endelse


;stop

mass = thedata.mass / 6.2508206e24

for i = 0, thedata.nenergy-1 do begin
	currbins = where(bins_2d(i,*) ne 0 and thedata.energy(i,*) le erange(1) and thedata.energy(i,*) ge erange(0) and finite(thedata.data(i,*)) eq 1,nbins)
	if nbins ne 0 then begin
;		print, i
		x = fltarr(nbins) & y = fltarr(nbins) & z = fltarr(nbins)
		sphere_to_cart,1,reform(thedata.theta(i,currbins)),reform(thedata.phi(i,currbins)),x,y,z
		if not keyword_set(plotenergy) then begin
			totalx = [totalx, x * reform(sqrt(2*1.6e-19*thedata.energy(i,currbins)/mass))]
			totaly = [totaly, y * reform(sqrt(2*1.6e-19*thedata.energy(i,currbins)/mass))]
			totalz = [totalz, z * reform(sqrt(2*1.6e-19*thedata.energy(i,currbins)/mass))]
		endif else begin
			totalx = [totalx, x * reform(thedata.energy(i,currbins))]
			totaly = [totaly, y * reform(thedata.energy(i,currbins))]
			totalz = [totalz, z * reform(thedata.energy(i,currbins))]
		endelse

		ncounts = [ncounts,reform(thedata.data(i, currbins))]
	endif
endfor

totalx = totalx(1:*)
totaly = totaly(1:*)
totalz = totalz(1:*)
ncounts = ncounts(1:*)

if in eq 0 then begin
  ncounts_t = ncounts
endif else begin
  ncounts_t = ncounts_t+ncounts
endelse

endfor
ncounts=ncounts_t/(inumber+1)

;  %%%  LOOP ENDS HERE  %%%

;*****HERES SOMETHING NEW
;sto
newdata = {dir:fltarr(n_elements(totalx),3), n:fltarr(n_elements(totalx))}

newdata.dir(*,0) = totalx
newdata.dir(*,1) = totaly
newdata.dir(*,2) = totalz
newdata.n = ncounts

;stop


;**********************************************

;get the magnetic field into a variable

get_data,thebdata,data = mgf



;************EXPERIMENTAL INTERPOLATION FIX************
;get_data,thebdata,data = bdata
;index = where(bdata.x le thedata.time + 600 and bdata.x ge thedata.time - 600)
;store_data,thebdata+'cut',data={x:bdata.x(index),y:bdata.y(index,*)}
;********

store_data,'time',data = {x:thedata.time+thedata.integ_t*.5}
;print, thedata.integ_t, ' Thedata.integ_t'
;interpolate,'time',thebdata+'cut','Bfield'
;get_data,'Bfield',data = mgf
;bfield = fltarr(3)
;bfield[0] = mgf.y(0,0)
;bfield[1] = mgf.y(0,1)
;bfield[2] = mgf.y(0,2)

bfield = dat_avg(thebdata, thedata3(0).time, thedata.end_time)

;print, 'BFIELD is ',bfield

;print, 'All data interpolated to ' + time_string(mgf.x)



if keyword_set(nosubtract) then print,'No velocity transform' else begin
	if keyword_set(vel) then print,'Velocity used for subtraction is '+vel else print, 'Velocity used for subtraction is V_3D'
endelse


if keyword_set(vel) then begin
	print,'Using '+vel+' for velocity vector'

;	get_data,vel,data = dummy, index = theindex
;	if theindex eq 0 then begin
;		Print, 'Loading velocity data....'
;		get_3dt,'v_3d','ph',/nr,/rm2
;	endif

;	interpolate,'time',vel,'value'
;	get_data,'value',data = thevalue
;	thevel = 1000.* reform(thevalue.y)

	thevel = 1000. * dat_avg(vel, thedata3(0).time, thedata.end_time)

;	print, thevel

	if keyword_set(plotenergy) then factor = sqrt(total(thevel(*)^2))*mass/2./1.6e-19 else factor = 1.
endif else begin
	print, 'Calculating V with v_3d...'
	thevel = 1000. * v_3d(thedata)
	thevel = 0.01 * j_3d(thedata)/n_3d(thedata)
	for in=0,inumber do begin
	  if in eq 0 then begin
		flux=j_3d(thedata3(0))
		density=n_3d(thedata3(0))
	  endif else begin
	    flux=flux+j_3d(thedata3(in))
	    density=density+n_3d(thedata3(in))
	  endelse
	endfor
	thevel = 0.01 * flux/density
	if keyword_set(plotenergy) then factor = sqrt(total(thevel(*)^2))*mass/2./1.6e-19 else factor = 1.
endelse


if not keyword_set(nosubtract) then begin
	newdata.dir(*,0) = newdata.dir(*,0) - thevel(0)*factor
	newdata.dir(*,1) = newdata.dir(*,1) - thevel(1)*factor
	newdata.dir(*,2) = newdata.dir(*,2) - thevel(2)*factor
endif else begin
	newdata.dir(*,0) = newdata.dir(*,0)
	newdata.dir(*,1) = newdata.dir(*,1)
	newdata.dir(*,2) = newdata.dir(*,2)
endelse





;**************NOW CONVERT TO THE DATA SET REQUIRED*****************


if rotation eq 'BV' then rot=cal_rot(bfield,thevel)
if rotation eq 'BE' then rot=cal_rot(bfield,crossp(bfield,thevel))
if rotation eq 'xy' then rot=cal_rot([1,0,0],[0,1,0])
if rotation eq 'xz' then rot=cal_rot([1,0,0],[0,0,1])
if rotation eq 'yz' then rot=cal_rot([0,1,0],[0,0,1])
if rotation eq 'xvel' then rot=cal_rot([1,0,0],thevel)
if rotation eq 'perp' then begin
    rot=cal_rot(crossp(bfield,crossp(bfield,thevel)),crossp(bfield,thevel))
endif
if rotation eq 'perp_yz' then begin
    rot=cal_rot(CROSSP(CROSSP(bfield,[0,1,0]),bfield),CROSSP(CROSSP(bfield,[0,0,1]),bfield))
endif
if rotation eq 'perp_xy' then begin
    rot=cal_rot(CROSSP(CROSSP(bfield,[1,0,0]),bfield),CROSSP(CROSSP(bfield,[0,1,0]),bfield))
endif
if rotation eq 'perp_xz' then begin
    rot=cal_rot(CROSSP(CROSSP(bfield,[1,0,0]),bfield),CROSSP(CROSSP(bfield,[0,0,1]),bfield))
endif

newdata.dir = newdata.dir#rot
if keyword_set(plotenergy) then factor=1. else factor = 1000.
;vperp = (newdata.dir(*,1)^2 + newdata.dir(*,2)^2)^.5*newdata.dir(*,1)/abs(newdata.dir(*,1))/factor
vperp = newdata.dir(*,1)/factor
vpara = newdata.dir(*,0)/factor
vperp2= newdata.dir(*,2)/factor
zdata = newdata.n



;new 2dslice stuff
;****************

;test only
;index = where(zdata le 0.000002 or zdata ge 0.000004, count)
;if count ne 0 then begin
;    vperp=vperp(index)
;    vpara=vpara(index)
;    zdata=zdata(index)
;endif

if not keyword_set(angle) then begin
	zmag = newdata.dir(*,2)/factor
	index = where(abs(zmag) le zheight,count)
	if count ne 0 then begin
		vperp = vperp(index)
		vpara = vpara(index)
	    vperp2= vperp2(index)
		zdata = zdata(index)
	endif else begin
		message,'NO DATA POINTS AT THAT ZHEIGHT!'
		return
	endelse
	print, 'zheight = ',zheight
endif else begin

if keyword_set(ThirdDirlim) then angle = [-90.,90.]


	if angle(0) eq -1 then angle = [-20.,20.]

	zmag = vperp2

	r = sqrt(vpara^2 + vperp^2+vperp2^2)

	eachangle = asin(zmag/r)
	angle1=min(angle)
	angle2=max(angle)

	index = where(eachangle/!dtor le angle2 and eachangle/!dtor ge angle1,count)
	if count ne 0 then begin
		vperp = vperp(index)
		vpara = vpara(index)
		vperp2=vperp2(index)
		zdata = zdata(index)
	endif else begin
		message,'NO DATA POINTS AT THAT ANGLE!'
		return
	endelse
	print, 'angle = ',angle
endelse

if keyword_set(ThirdDirlim) then begin
    third = vperp2
    index = where(third le max(ThirdDirlim) and third ge min(ThirdDirlim))
    if count ne 0 then begin
        vperp = vperp(index)
        vpara = vpara(index)
        vperp2=vperp2(index)
        zdata = zdata(index)
    endif
endif

;**********************

if keyword_set(sundir) then begin
	sund = [1,0,0]
	sund = sund#rot
	vperpsun = (sund(1)^2 + sund(2)^2)^.5*sund(1)/abs(sund(1))
	vparasun = sund(0)
endif

if not keyword_set(vel) then veldir = v_3d(thedata) else veldir = thevel/1000.
veldir = veldir#rot

;EXPERIMENTAL GET RID OF 0 THING*************
if not keyword_set(leavezero) then begin
	index = where(zdata ne 0)
	vperp = vperp(index)
	vpara = vpara(index)
	vperp2=vperp2(index)
	zdata = zdata(index)
endif else print, 'Zeros left in plot'

;MAKE SURE THERE ARE NO NEGATIVE VALUES!! ***********
index2 = where(zdata lt 0., count)
if count ne 0 then print,'THERE ARE NEGATIVE DATA VALUES'

index = where(zdata ge 0,count)
if count ne 0 then begin
	vperp = vperp(index)
	vperp2= vperp2(index)
	vpara = vpara(index)
	zdata = zdata(index)
endif


;stop

vperp=vperp(sort(vpara))
zdata=zdata(sort(vpara))
vperp2=vperp2(sort(vpara))
vpara=vpara(sort(vpara))

uni2=uniq(vpara)
uni1=[0,uni2(0:n_elements(uni2)-2)+1]

kk=0
for i=0,n_elements(uni2)-1 do begin
    vperpi=vperp(uni1(i):uni2(i))
    vparai=vpara(uni1(i):uni2(i))
    zdatai=zdata(uni1(i):uni2(i))

    vparai=vparai(sort(vperpi))
    zdatai=zdatai(sort(vperpi))
	vperpi=vperpi(sort(vperpi))

	index2=uniq(vperpi)
	if n_elements(index2) eq 1 then begin
	    index1=0
	endif else begin
	    index1=[0,index2(0:n_elements(index2)-2)+1]
	endelse

    for j=0,n_elements(index2)-1 do begin
        vperp(kk)=vperpi(index1(j))
        vpara(kk)=vparai(index1(j))
        if index1(j) eq index2(j) then begin
            zdata(kk)=zdatai(index1(j))
        endif else begin
            zdata_moment=moment(zdatai(index1(j):index2(j)))
            zdata(kk)=zdata_moment(0)
        endelse
        kk=kk+1
    endfor
endfor
vperp=vperp(0:kk-1)
vpara=vpara(0:kk-1)
zdata=zdata(0:kk-1)


;******************NOW TO PLOT THE DATA********************

if not keyword_set(xrange) then begin
	themax = max(abs([vperp,vpara]))
	xrange = [-1*themax,themax]
endif else themax = max(abs(xrange))


if not keyword_set(range) then begin
	if not keyword_set(xrange) then begin
		maximum = max(zdata)
		minimum = min(zdata(where(zdata ne 0)))
	endif else begin
		maximum = max(zdata(where(abs(vperp) le themax and abs(vpara) le themax)))
		minimum = min(zdata(where(zdata ne 0 and abs(vperp) le themax and abs(vpara) le themax)))
	endelse
endif else begin
	maximum = range(1)
	minimum = range(0)
endelse



if keyword_set(zlog) then $
	thelevels = 10.^(indgen(nlines)/float(nlines)*(alog10(maximum) - alog10(minimum)) + alog10(minimum)) $
else $
	thelevels = (indgen(nlines)/float(nlines)*(maximum-minimum)+minimum)
;**********EXTRA STUFF FOR THE CONTOUR LINE OVERPLOTS************
if keyword_set(olines) then begin
	if keyword_set(zlog) then $
		thelevels2 = 10.^(indgen(olines)/float(olines)*(alog10(maximum) - alog10(minimum)) + alog10(minimum)) $
	else $
		thelevels2 = (indgen(olines)/float(olines)*(maximum-minimum)+minimum)

endif
;**********END EXTRA STUFF FOR LINE OVERPLOTS (MORE LATER)*************************************


thecolors = round((indgen(nlines)+1)*(!d.table_size-9)/nlines)+7

if not keyword_set(nofill) then fill = 1 else fill = 0

if not keyword_set(finished) then begin
	if not keyword_set(plotenergy) then begin
	    if rotation eq 'BV' then begin
	        xtitle = 'V Para (km/sec)'
	        ytitle = 'V Perp (km/sec)'
	    endif
	    if rotation eq 'xy' then begin
	        xtitle = 'Vx (km/sec)'
	        ytitle = 'Vy (km/sec)'
	    endif
	    if rotation eq 'xz' then begin
	        xtitle = 'Vx (km/sec)'
	        ytitle = 'Vz (km/sec)'
	    endif
	    if rotation eq 'yz' then begin
	        xtitle = 'Vy (km/sec)'
	        ytitle = 'Vz (km/sec)'
	    endif
	    if rotation eq 'xvel' then begin
	        xtitle = 'Vx (km/sec)'
	        ytitle = 'Vyz (km/sec)'
	    endif
	endif else begin
		xtitle = 'E Para (eV)'
		ytitle = 'E Perp (eV)'
	endelse
endif else begin
	if not keyword_set(plotenergy) then begin
		xtitle = 'V!19!D'+parasymbol+'!N!7 (km/sec)'
		ytitle = 'V!19!D'+perpsymbol+'!N!7 (km/sec)'
	endif else begin
		xtitle = 'E!19!D'+parasymbol+'!N!7 (eV)'
		ytitle = 'E!19!D'+perpsymbol+'!N!7 (eV)'
	endelse
endelse


if keyword_set(grid) then begin
;stop
	x= findgen(resolution)/(resolution-1)*(xrange(1)-xrange(0)) + xrange(0)
	spacing = (xrange(1)-xrange(0))/(resolution-1)
	triangulate,vpara,vperp,tr,b

; test
tr1=tr(0,*)
tr2=tr(1,*)
tr3=tr(2,*)
index = where((vpara(tr1)+vpara(tr2)+vpara(tr3))^2+(vperp(tr1)+vperp(tr2)+vperp(tr3))^2 gt min(vpara^2+vperp^2), count)
if count ne 0 then begin
    tr=tr(*,index)
endif


	thesurf = trigrid(vpara,vperp,zdata,tr,[spacing,spacing], [xrange(0),xrange(0),xrange(1),xrange(1)],xgrid = xg,ygrid = yg )
	if keyword_set(smooth) then thesurf = smooth(thesurf,3)
	if n_elements(xg) mod 2 ne 1 then print,'The line plots are invalid',n_elements(xg)
	;print,n_elements(xg)
;**************************************************
;********EXPERIMENTAL THINGS HERE************
;stop

if keyword_set(logplot) then begin

	vpara2 = vpara
	vperp2 = vperp

	magnitude = .5 * alog(vpara^2 + vperp^2)

	vpara2 = vpara / sqrt(vpara^2 + vperp^2)
	vperp2 = vperp / sqrt(vpara^2 + vperp^2)

	vpara2 = vpara2 * magnitude
	vperp2 = vperp2 * magnitude

;CONVERT BACK AS A CHECK
;magnitude = exp(sqrt(vpara2^2 + vperp2^2) )
;vpara3 = vpara2 / sqrt(vpara2^2 + vperp2^2)
;vperp3 = vperp2 / sqrt(vpara2^2 + vperp2^2)
;vpara4 = vpara3 * magnitude
;vperp4 = vperp3 * magnitude

	xrangeold = xrange

	xrange(0) = -alog(abs(xrange(0)))
	xrange(1) = alog(xrange(1))

;stop


	x= findgen(resolution)/(resolution-1)*(xrange(1)-xrange(0)) + xrange(0)
	spacing = (xrange(1)-xrange(0))/(resolution-1)
	triangulate,vpara2,vperp2,tr,b
	thesurf = trigrid(vpara2,vperp2,zdata,tr,[spacing,spacing], [xrange(0),xrange(0),xrange(1),xrange(1)],xgrid = xg,ygrid = yg )
	if keyword_set(smooth) then thesurf = smooth(thesurf,3)
	if n_elements(xg) mod 2 ne 1 then print,'The line plots are invalid',n_elements(xg)
	;print,n_elements(xg)

	xrange = xrangeold

	indexminus = where(xg lt 0.)
	indexplus = where(xg gt 0.)

	xg(indexminus) = -exp(abs(xg(indexminus)))
	xg(indexplus) = exp(xg(indexplus))
	yg(indexminus) = -exp(abs(yg(indexminus)))
	yg(indexplus) = exp(yg(indexplus))

endif
;stop
;********************************************************
;********************************************************

 ;thedata.data_name+' '+time_string(thedata.time)
	timetitle = thedata.data_name+' '+time_string(thedata3(0).time) + '->' + strmid(time_string(thedata.end_time),11,8)
	;if keyword_set(finished) and keyword_set(plotlabel) then timetitle = '!B
	;if keyword_set(finished) and keyword_set(plotlabel) then timetitle = '!B
	contour,thesurf,xg,yg,$
		/closed,levels=thelevels,c_color = thecolors,fill=fill,$
		title = timetitle, $
		ystyle = 1,$
		ticklen = -0.01,$
		xstyle = 1,$
		xrange = xrange,$
		yrange = xrange,$
		xtitle = xtitle,$
		ytitle = ytitle,position = position
	if keyword_set(olines) then begin
		if !d.name eq 'PS' then somecol = !p.color else somecol = 0
		contour, thesurf,xg,yg,/closed,levels = thelevels2,ystyle = 1+4, $
			xstyle = 1+4,xrange = xrange, yrange = xrange, ticklen = 0,/noerase,position = position,col = somecol
	endif
endif else begin
	contour,zdata,vpara,vperp,/irregular,$
		/closed,levels=thelevels,c_color = thecolors,fill=fill,$
		title = timetitle, $
		ystyle = 1,$
		ticklen = -0.01,$
		xstyle = 1,$
		xrange = xrange,$
		yrange = xrange,$
		xtitle = xtitle,$
		ytitle = ytitle,position = position
	if keyword_set(olines) then begin
		if !d.name eq 'PS' then somecol = !p.color else somecol = 0
		contour, zdata,vpara,vperp,/irregular,/closed, levels = thelevels2, $
			ystyle = 1+4, xstyle = 1+4, ticklen = 0, xrange = xrange, yrange = xrange, position=position,/noerase, col = somecol
	endif
endelse

if not keyword_set(cut_para) then cut_para = 0.
if not keyword_set(cut_perp) then cut_perp = 0.

if keyword_set(cut_bulk_vel) then begin
cut_para= veldir(0)
cut_perp= veldir(1)
endif

;oplot, vpara, vperp, PSYM=1

oplot,[cut_para,cut_para],xrange,linestyle = 2,thick = 2
oplot,xrange,[cut_perp,cut_perp],linestyle = 2,thick = 2
;oplot,[0,0],xrange,linestyle = 1
;oplot,xrange,[0,0],linestyle = 1

if keyword_set(sundir) then oplot,[0,vparasun*max(xrange)],[0,vperpsun*max(xrange)]

if keyword_set(vel2) then begin

;	stop
	vel2 = vel2#rot
	vperpvel2 = (vel2(1)^2 + vel2(2)^2)^.5*vel2(1)/abs(vel2(1))
	vparavel2 = vel2(0)

	bbbb=findgen(36)*(!pi*2/32.)
	usersym,1.5*cos(bbbb),1.5*sin(bbbb),/fill

;	oplot,[vparavel2],[vperpvel2],psym = 8,col= !d.table_size - 10,symsize =1
	oplot,[vparavel2],[vperpvel2],psym = 8,col= 2,symsize =1
endif


if not keyword_set(novelline) then oplot,[0,veldir(0)],[0,veldir(1)],col= !d.table_size-9


if not keyword_set(plotenergy) then begin
	circy=sin(findgen(360)*!dtor)*sqrt(2.*1.6e-19*erange(0)/mass)/1000.
	circx=cos(findgen(360)*!dtor)*sqrt(2.*1.6e-19*erange(0)/mass)/1000.  ;sqrt(2*1.6e-19*energy(i)/mass)
	oplot,circx,circy,thick = 2

	circy=sin(findgen(360)*!dtor)*sqrt(2.*1.6e-19*erange(1)/mass)/1000.
	circx=cos(findgen(360)*!dtor)*sqrt(2.*1.6e-19*erange(1)/mass)/1000.  ;sqrt(2*1.6e-19*energy(i)/mass)
	oplot,circx,circy,thick = 2
endif else begin
	circy=sin(findgen(360)*!dtor)*erange(0)
	circx=cos(findgen(360)*!dtor)*erange(0)  ;sqrt(2*1.6e-19*energy(i)/mass)
	oplot,circx,circy,thick = 2

	circy=sin(findgen(360)*!dtor)*erange(1)
	circx=cos(findgen(360)*!dtor)*erange(1)  ;sqrt(2*1.6e-19*energy(i)/mass)
	oplot,circx,circy,thick = 2
endelse


if not keyword_set(finished) then thetitle = units_string(thedata.units_name) $
else thetitle = units_finalstring(thedata.units_name)

if keyword_set(plotlabel) then xyouts, 0.05,.95,plotlabel+'!N!7',/normal,charsize = 1.5


;if keyword_set(zlog) then thetitle = thetitle + ' (log)'

draw_color_scale,range=[minimum,maximum],log = zlog,yticks=10,title =thetitle


if keyword_set(showdata) then oplot,vpara,vperp,psym=1


if keyword_set(var_label) then begin
	get_data,'time',data = theparticulartime
	;print, 'Vars interp. to '+time_string(theparticulartime.x)

	if keyword_set(vel) then velname = vel else velname = 'V_3D'

	if keyword_set(nosubtract) then outstring = velname + ' used (no trans)' else outstring = vel + 'used (trans)'
	for i = 0, n_elements(var_label)-1 do begin
		print,var_label(i)

;		get_data,var_label(i),data = bdata
;		index = where(bdata.x le thedata.time + 600 and bdata.x ge thedata.time - 600,count)
;		if count ne 0 then begin
;			store_data,var_label(i)+'cut',data={x:bdata.x(index),y:bdata.y(index,*)}
;			interpolate,'time',var_label(i)+'cut','value'
;			get_data,'value',data = thevalue
;			outstring = outstring + '  ' + var_label(i) +'= '+ strtrim(string(format = '(G11.4)',thevalue.y(0)),2)
;		endif

		thevalue = dat_avg(var_label(i), thedata.time, thedata.end_time)

		outstring = outstring + '  ' + var_label(i) +'= '+ strtrim(string(format = '(G11.4)',thevalue),2)

		if i mod numperrow eq 1 then outstring = outstring + '!c'
	endfor
	xyouts,.13,.06,outstring,/normal
endif

if cross eq 1 then begin
	n_elem = n_elements(thesurf(*,0))



	if not keyword_set(cut_perp) then perpval = n_elem/2 else begin
		ind = where(xg ge cut_perp)
		if (xg(ind(0)) - cut_perp) le (cut_perp - xg(ind(0)-1) ) then perpval = ind(0) else perpval = ind(0)-1
	endelse

	if not keyword_set(cut_para) then paraval = n_elem/2 else begin
		ind = where(xg ge cut_para)
		if (xg(ind(0)) - cut_para) le (cut_para - xg(ind(0)-1) ) then paraval = ind(0) else paraval = ind(0)-1
	endelse

;	if keyword_set(zlog) then thetitle = thetitle + ' (log)'

	if keyword_set(plotenergy) then begin
		xtitle = 'Energy (eV)'
		vore = 'E'
	endif else begin
		xtitle = 'Velocity (km/sec)'
		vore = 'V'
	endelse


;HERE COMES SOME NEW COLOR STUFF
	if !d.name eq 'PS' then thecolors = round((indgen(4)+1)*(!d.table_size-9)/4)+7 else begin
		thecolors=indgen(4)
		thecolors = thecolors + 3
	endelse

	if keyword_set(double) then begin
		;first plot vpara on the + side
		plot,xg,[reverse(thesurf(n_elem/2:*,perpval)),thesurf(n_elem/2+1:*,perpval)],$
			xstyle = 1, ystyle =1,$
			xrange = xrange,yrange = [minimum,maximum],ylog = zlog, $
			title = 'Cross Sections',xtitle = xtitle,ytitle = thetitle,$
			position = pos2
		;overplot vpara on the minus side
		oplot,xg,[thesurf(0:n_elem/2,perpval),reverse(thesurf(0:n_elem/2-1,perpval))],color = thecolors(0)
		;now vperp on the + side
		oplot,xg,[reverse(reform(thesurf(paraval,n_elem/2+1:*))),reform(thesurf(paraval,n_elem/2:*))],color = thecolors(1)
		;and now vper on the - side
		oplot,xg,[reform(thesurf(paraval,0:n_elem/2)),reverse(reform(thesurf(paraval,0:n_elem/2-1)))],color = thecolors(2)
	endif else begin
		if not keyword_set(loglines) then begin
			;stop
			;first plot vpara on the + side
			plot,xg(n_elem/2:*),thesurf(n_elem/2:*,perpval),$
				xstyle = 1, ystyle =1,$
				xrange = xrange,yrange = [minimum,maximum],ylog = zlog, $
				title = 'Cross Sections',xtitle = xtitle,ytitle = thetitle,$
				position = pos2
			;overplot vpara on the minus side
			oplot,xg(0:n_elem/2),thesurf(0:n_elem/2,perpval),color = thecolors(0)
			;now vperp on the + side
			oplot,xg(n_elem/2:*),reform(thesurf(paraval,n_elem/2:*)),color = thecolors(1)
			;and now vper on the - side
			oplot,xg(0:n_elem/2),reform(thesurf(paraval,0:n_elem/2)),color = thecolors(2)
		endif else begin
			;*****PUT IN LOGLINES STUFF HERE********
;********EXPERIMENTAL THINGS HERE************
;stop

if not keyword_set(oldlog) then begin

vpara2 = vpara
vperp2 = vperp

magnitude = .5 * alog(vpara^2 + vperp^2)

vpara2 = vpara / sqrt(vpara^2 + vperp^2)
vperp2 = vperp / sqrt(vpara^2 + vperp^2)

vpara2 = vpara2 * magnitude
vperp2 = vperp2 * magnitude

;CONVERT BACK AS A CHECK
;magnitude = exp(sqrt(vpara2^2 + vperp2^2) )
;vpara3 = vpara2 / sqrt(vpara2^2 + vperp2^2)
;vperp3 = vperp2 / sqrt(vpara2^2 + vperp2^2)
;vpara4 = vpara3 * magnitude
;vperp4 = vperp3 * magnitude

xrangeold = xrange

xrange(0) = -alog(abs(xrange(0)))
xrange(1) = alog(xrange(1))

;stop


	x= findgen(resolution)/(resolution-1)*(xrange(1)-xrange(0)) + xrange(0)
	spacing = (xrange(1)-xrange(0))/(resolution-1)
	triangulate,vpara2,vperp2,tr,b
	thesurf = trigrid(vpara2,vperp2,zdata,tr,[spacing,spacing], [xrange(0),xrange(0),xrange(1),xrange(1)],xgrid = xg,ygrid = yg )
	if keyword_set(smooth) then thesurf = smooth(thesurf,3)
	if n_elements(xg) mod 2 ne 1 then print,'The line plots are invalid',n_elements(xg)
	;print,n_elements(xg)

xrange = xrangeold

indexminus = where(xg lt 0.)
indexplus = where(xg gt 0.)

xg(indexminus) = -exp(abs(xg(indexminus)))
xg(indexplus) = exp(xg(indexplus))

endif
;stop
;********************************************************
			;stop
			;plot,xg,[reverse(thesurf(n_elem/2:*,perpval)),thesurf(n_elem/2+1:*,perpval)],$
			plot,xg(n_elem/2:*),thesurf(n_elem/2:*,perpval),$
				xstyle = 1, ystyle =1,$
				xrange = [min(thedata.energy(*,0)),xrange[1]],yrange = [minimum,maximum],ylog = zlog,/xlog, $
				title = 'Cross Sections',xtitle = xtitle,ytitle = thetitle,$
				position = pos2
			;vpara on the minus side
			oplot,xg(n_elem/2:*), reverse(thesurf(0:n_elem/2, perpval)), color = thecolors(0)
			;vperp on the + side
			oplot,xg(n_elem/2:*), reform(thesurf(paraval,n_elem/2:*)),color = thecolors(1)
			;vperp on the - side
			oplot,xg(n_elem/2:*), reverse(reform(thesurf(paraval, 0:n_elem/2))), color = thecolors(2)


		endelse
	endelse



	;put a dotted line
	oplot,[0,0],[minimum,maximum],linestyle = 1
	if not keyword_set(plotenergy) then begin
		oplot,[sqrt(2.*1.6e-19*erange(0)/mass)/1000.,sqrt(2.*1.6e-19*erange(0)/mass)/1000.],[minimum,maximum],linestyle = 5
		oplot,-[sqrt(2.*1.6e-19*erange(0)/mass)/1000.,sqrt(2.*1.6e-19*erange(0)/mass)/1000.],[minimum,maximum],linestyle = 5
		oplot,[sqrt(2.*1.6e-19*erange(1)/mass)/1000.,sqrt(2.*1.6e-19*erange(1)/mass)/1000.],[minimum,maximum],linestyle = 5
		oplot,-[sqrt(2.*1.6e-19*erange(1)/mass)/1000.,sqrt(2.*1.6e-19*erange(1)/mass)/1000.],[minimum,maximum],linestyle = 5
		if keyword_set(onecnt) then begin
			oplot,sqrt(2.*1.6e-19*theonecnt.energy(*,0)/mass)/1000.,theonecnt.data(*,0),color = thecolors(3),linestyle = 3
			oplot,-sqrt(2.*1.6e-19*theonecnt.energy(*,0)/mass)/1000.,theonecnt.data(*,0),color = thecolors(3),linestyle = 3
		endif
	endif else begin
		oplot,[erange(0),erange(0)],[minimum,maximum],linestyle = 5
		oplot,-[erange(0),erange(0)],[minimum,maximum],linestyle = 5
		oplot,[erange(1),erange(1)],[minimum,maximum],linestyle = 5
		oplot,-[erange(1),erange(1)],[minimum,maximum],linestyle = 5
		if keyword_set(onecnt) then begin
			oplot,theonecnt.energy(*,0),theonecnt.data(*,0),color = thecolors(6),linestyle = 3
			oplot,-theonecnt.energy(*,0),theonecnt.data(*,0),color = thecolors(6),linestyle = 3
		endif
	endelse

	;now put the titles on the side of the graph
	positions = -findgen(5)*(pos2(3)-pos2(1))/5 + pos2(3)-.03
	xyouts,pos2(2) + .03,positions(0),vore+' para (+ side)',/norm,charsize = 1.01;.5
	xyouts,pos2(2) + .03,positions(1),vore+' para (- side)',/norm,color = thecolors(0),charsize = 1.01;.5
	xyouts,pos2(2) + .03,positions(2),vore+' perp (+ side)',/norm,color = thecolors(1),charsize = 1.01;.5
	xyouts,pos2(2) + .03,positions(3),vore+' perp (- side)',/norm,color = thecolors(2),charsize = 1.01;.5
	if keyword_set(onecnt) then xyouts,pos2(2) + .03,positions(4),'One count',/norm,color = thecolors(3),charsize = .5
endif

;stop


if keyword_set(outfile) then begin
	openw, thefile, outfile, /get_lun
	printf, thefile, time_string(thedata.time)

	if not keyword_set(cut_perp) then perpval = n_elem/2 else begin
		ind = where(xg ge cut_perp)
		if (xg(ind(0)) - cut_perp) le (cut_perp - xg(ind(0)-1) ) then perpval = ind(0) else perpval = ind(0)-1
	endelse

	if not keyword_set(cut_para) then paraval = n_elem/2 else begin
		ind = where(xg ge cut_para)
		if (xg(ind(0)) - cut_para) le (cut_para - xg(ind(0)-1) ) then paraval = ind(0) else paraval = ind(0)-1
	endelse


	filedata = fltarr(3,n_elements(xg))

	filedata(0,*) = xg
	filedata(1,*) = thesurf(*,perpval)
	filedata(2,*) = thesurf(paraval,*)

	printf, thefile, filedata
	close,/all

endif

;********EXTRA PART*********

;	if not keyword_set(cut_perp) then perpval = n_elem/2 else begin
;		ind = where(xg ge cut_perp)
;		if (xg(ind(0)) - cut_perp) le (cut_perp - xg(ind(0)-1) ) then perpval = ind(0) else perpval = ind(0)-1
;	endelse

;	if not keyword_set(cut_para) then paraval = n_elem/2 else begin
;		ind = where(xg ge cut_para)
;		if (xg(ind(0)) - cut_para) le (cut_para - xg(ind(0)-1) ) then paraval = ind(0) else paraval = ind(0)-1
;	endelse


;	filedata = fltarr(3,n_elements(xg))

;	filedata(0,*) = xg
;	filedata(1,*) = thesurf(*,perpval)
;	filedata(2,*) = thesurf(paraval,*)
;outcuts = filedata

;********END EXTRA PART*******

if filetype eq 'ps' then begin
    DEVICE, /CLOSE
endif

if !d.name ne 'PS' then !p.multi = oldplot


end