;+
;NAME:			get_v_spec
;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).
;			Does NOT tranform into bulk flow frame (unlike get_v_spec_t)
;CALL:			ex: get_v_spec,get_el('20:31'),[keywords]
;KEYWORDS:		XRANGE: vector specifying the xrange
;			NOFIX: tells the program not to rebin the data
;			       with fixangdata.pro
;			NSTEPS: tells the program how many angle steps to use
;			        when it calls fixangdata.pro
;			RANGE: vector specifying the color range
;			VAR_LA: vector of tplot variables to show on plot
;			UNITS: specifies the units ('eflux','df',etc.)
;			ZLOG: specifies a logarithmic Z axis
;			B3: uses B3 data
;			POSITION: positions the plot using a 4-vector
;			ERANGE: specifies the energy range to be used
;			NOFILL: doesn't fill the contour plot
;			NLINES: says how many lines to use if using NOFILL
;			SHOWDATA: plots all the data points over the contour
;			PLOTENERGY: plots using energy instead of velocity
;LAST MODIFIED:		8-15-97
;CREATED BY:		Arjun Raj
;-





pro get_v_spec, thedata,nofix = nofix,nsteps = nsteps,xrange = xrange,range = range, units = units,zlog = zlog,b3 = b3,position = position,erange = erange,nofill = nofill,var_label = var_label,nlines = nlines,showdata = showdata,plotenergy = plotenergy,_EXTRA = e


if not keyword_set(position) then begin
	if !d.y_size le !d.x_size then $
		position = [.13,.13,.13 + .77 * !d.y_size/!d.x_size,.13 + .77] else $
		position = [.13,.13,.13 + .77,.13 + .77 *!d.x_size/!d.y_size]
endif

if keyword_set(var_label) then begin
	position(1) = position(1) + .06
	position(3) = position(3) + .06
endif

if not keyword_set(nsteps) then nsteps = 18
if not keyword_set(units) then units = 'eflux'
if not keyword_set(nlines) then nlines = 60
thedata = conv_units(thedata,units)
;**********************************************
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 print,'Not valid data'

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

;get the magnetic field into a variable
if keyword_set(B3) then begin
	get_data,'B3_gse',data = mgf,index = theindex
	if theindex eq 0 then begin
		get_mfi_3s
		get_data,'B3_gse',data = mgf
	endif
	Bname = 'B3_gse'
endif else begin
	get_data,'Bexp',data = mgf,index = theindex
	if theindex eq 0 then begin
		get_mfi
		get_data,'Bexp',data = mgf
	endif
	Bname = 'Bexp'
endelse

store_data,'time',data = {x:thedata.time}
interpolate,'time',Bname,'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)

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

data = {dir:fltarr(n_good,3),energy:reverse(thedata.energy(*,0)),n:fltarr(thedata.nenergy,n_good)}

x = fltarr(n_good) & y = fltarr(n_good) & z = fltarr(n_good)

sphere_to_cart,1,thedata.theta(0,good_bins),thedata.phi(0,good_bins),x,y,z
data.dir(*,0) = x & data.dir(*,1) = y & data.dir(*,2) = z

if units eq 'counts' then $
	for i = 0,thedata.nenergy - 1 do $
		data.n(i,*) = thedata.data(thedata.nenergy-1-i,good_bins)/thedata.geom $
 else $
	for i = 0,thedata.nenergy - 1 do $
		data.n(i,*) = thedata.data(thedata.nenergy-1-i,good_bins)

;now the variable data contains both a tag with the directions associated with all the bins
;as well as all the counts.

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

;toplot = {ang:fltarr(n_good),energy:data.energy(eindex),n:fltarr(n_elements(eindex),n_good)}
toplot = {ang:fltarr(n_good),energy:data.energy,n:fltarr(n_elements(data.energy),n_good)}

angles = angl(data.dir,bfield)
angindex = sort(angles)
toplot.ang = angles(angindex)


;toplot.n(*,*)=data.n(eindex,*)
toplot.n(*,*)=data.n(*,*)
toplot.n(*,*)=toplot.n(*,angindex)

if not keyword_set(nofix) then newdata = fixangdata(toplot,nsteps,/noextra) else newdata = toplot

uneindex = where(newdata.energy lt erange(0) or data.energy gt erange(1),counts)
if counts ne 0 then newdata.n(uneindex,*)=0

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

if strpos(thedata.data_name, 'Eesa') ne -1 then mass = 9.1e-31 $
	else mass = 1.67e-27


vperp = fltarr(1)
vpara = fltarr(1)
zdata = fltarr(1)

for i = 0, n_elements(newdata.energy)-1 do begin
	if not keyword_set(plotenergy) then begin
		add_perp = sqrt(2*1.6e-19 * newdata.energy(i) / mass)*sin(!dtor * newdata.ang)
		add_para = sqrt(2*1.6e-19 * newdata.energy(i) / mass)*cos(!dtor * newdata.ang)
	endif else begin
		add_perp = newdata.energy(i)*sin(!dtor*newdata.ang)
		add_para = newdata.energy(i)*cos(!dtor*newdata.ang)
	endelse
	add_z = newdata.n(i,*)
	vperp = [vperp,add_perp]
	vpara = [vpara,add_para]	
	zdata = [zdata,add_z(*)]
	vperp = [vperp,-1*add_perp]
	vpara = [vpara,add_para]
	zdata = [zdata,add_z(*)]

endfor

vperp = vperp(1:*)
vpara = vpara(1:*)
zdata = zdata(1:*)

if not keyword_set(plotenergy) then begin
	vperp = vperp/1000.
	vpara = vpara/1000.
endif

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


add_z = fltarr(n_elements(newdata.ang)*2)
add_z(*) = 0
if not keyword_set(energy) then begin
	add_para = sqrt(2*1.6e-19 * (newdata.energy(0)-5) / mass)*cos(!dtor * newdata.ang)/1000.
	add_perp = sqrt(2*1.6e-19 * (newdata.energy(0)-5) / mass)*sin(!dtor * newdata.ang)/1000.
endif else begin
	add_para = (newdata.energy(0)-2)*cos(!dtor*newdata.ang)
	add_perp = (newdata.energy(0)-2)*cos(!dtor*newdata.ang)
endelse
add_para = [add_para,add_para]
add_perp = [add_perp,-add_perp]


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)

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(plotenergy) then begin
	xtitle = 'V Para (km/sec)'
	ytitle = 'V Perp (km/sec)'
endif else begin
	xtitle = 'E Para (eV)'
	ytitle = 'E Perp (eV)'
endelse

contour,[zdata,add_z],[vpara,add_para],[vperp,add_perp],/irregular,$
		/closed,levels=thelevels,c_color = thecolors,fill=fill,$
		title = thedata.data_name+' '+time_string(thedata.time) ,$
		ystyle = 1,$
		ticklen = -0.01,$
		xstyle = 1,$
		xrange = xrange,$
		yrange = xrange,$
		xtitle = xtitle,$
		ytitle = ytitle,position = position

thetitle = thedata.units_name
if keyword_set(zlog) then thetitle = thetitle + ' (log)'
draw_color_scale,range=[minimum,maximum],log = zlog,yticks=10,title =thetitle

;wi,5
;print,'interpolating...'
;;z = tri_surf(zdata,vpara,vperp,bounds=[-themax,-themax,themax,themax])
;z = tri_surf(zdata,vpara,vperp,nx=100,ny=100);

;x = findgen(100)*2*themax/99-themax
;y = findgen(100)*2*themax/99-themax

;help,/str,x
;help,/str,y
;help,/str,z

;stop
;contour,z,x,y,$
;		/closed,levels=thelevels,c_color = thecolors,fill=fill,$
;		title = thedata.data_name+' '+time_string(thedata.time) ,$
;		ystyle = 1,$
;		ticklen = -0.01,$
;		xstyle = 1,$
;		xrange = xrange,$
;		yrange = xrange,$
;		xtitle = 'V Para (km/sec)',$
;		ytitle = 'V Perp (km/sec)',position = position

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


if keyword_set(var_label) then begin
	outstring = ''
	for i = 0, n_elements(var_label)-1 do begin
		print,var_label(i)
		interpolate,'time',var_label(i),'value'
		get_data,var_label(i),data = thename
		get_data,'value',data = thevalue
		outstring = outstring + '  ' + var_label(i) +'= '+ strtrim(string(format = '(G11.4)',thevalue.y(0)),2)
		if i mod 2 eq 1 then outstring = outstring + '!c' 
	endfor
	xyouts,.13,.10,outstring,/normal
endif


end