;NAME:		get_mfi_spec
;PURPOSE:	To make tplot variables which will make spectrograms of	
;			3dp data vs. magnetic field direction vs. time
;			for each energy level
;
;INPUT:		DATA: A string of the data type the procedure should get (eg
;			'ph2','pl','el','eh','so','sf')
;
;KEYWORDS:	B3:    Set this keyword to load the 3-second magnetic
;			field data.
;		ZLOG:  This keyword will cause the tplot structures to
;			plot the z-axis of the spectrogram logrithmically.
;		UNITS: Set this keyword to a string containing the output
;			units, such as eflux, flux, df, ncounts, rate,
;			or nrate (default is counts).
;
;OUTPUTS:	Tplot variables named 'B_data_n', where n is the 
;			energy level the variable corresponds
;			to and data is the data type (eg. ph2).
;
;EXAMPLE:	load_3dp_data,'96-03-27/10:',1
;		get_mfi_spec,'so',/zlog,/B3,units='dflux'
;		tplot,['B_so_0','B_so_8']
;
;CREATED BY:	Mike Moyer  96-11-05
;UPDATED BY:	Arjun Raj 97-07-28 (fixed the unit conversions)


pro get_mfi_spec,get_dat,zlog=zlog,units=units,B3=B3,noload = noload

starting_t=systime(1)


;*******ADDED BY ARJUN RAJ******
if not keyword_set(units) then units = 'df'
;*********************************


;setting default values
if not keyword_set(get_dat) then begin
	dat_name='ph2'
	get_dat='get_ph2' 
endif else begin
	dat_name=get_dat
	get_dat='get_'+get_dat
endelse
	
if not keyword_set(zlog) then zlog=0

s=call_function(get_dat,/times)
n_points=n_elements(s)
if n_points le 1 then begin
	print,'Sorry, No Data.'
	return
endif

t=dblarr(2)
 t(0)=s(0)
 t(1)=s(n_points-1)

;*****
;this gets the data for the time interval into the array named data,
;each element of which is the 3d structure for that time interval
;*****
d=call_function(get_dat,t(0))
while d.valid eq 0 do begin
	d=call_function(get_dat,t(0),/advance)
endwhile
data=replicate(d,n_points)
tags=n_tags(d)
n_energies=d.nenergy
n_bins=d.nbins

for i=0,n_points-1 do begin
	q=call_function(get_dat,(s(i)))

	if q.valid eq 0 then begin
		q=d
		q.data=!values.f_nan
		q.time=s(i)
		q.end_time=s(i)+(s(i)-s(i-1))
	endif

	if keyword_set(units) then q=conv_units(q,units)
	for j=0,tags-1 do begin 
		data(i).(j)=q.(j)
	endfor
endfor

;*****

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

if n_bad ne 0 then data(*).geom(bad_bins) = !values.f_nan
	
;to get normalized B field coords for each data point in e1,e2,e3 coords

if keyword_set(B3) then begin
	if not keyword_set(noload) then get_mfi_3s
	;get_data,'B3_gse',data=bold2
endif else begin
	if not keyword_set(noload) then get_mfi
	;get_data,'Bexp',data=bold2
endelse
;*****BAD WAY******
;good_b=where((bold2.x ge t(0)) and (bold2.x le t(1)),n_good_b)
;bold={x:bold2.x(good_b),y:bold2.y(good_b,*)}
;b={x:s,y:fltarr(n_elements(s),3)}
;for i=0,2 do begin
;	b.y(*,0)=interpol(bold.y(*,0),bold.x,b.x)
;	b.y(*,1)=interpol(bold.y(*,1),bold.x,b.x)
;	b.y(*,2)=interpol(bold.y(*,2),bold.x,b.x)
;endfor
;cart_to_sphere,b.y(*,0),b.y(*,1),b.y(*,2),r,bth,bph

;sphere_to_cart,1,bth,bph,bx,by,bz

;bdir=fltarr(n_points,3)
;bdir(*,0)=bx(*)
;bdir(*,1)=by(*)
;bdir(*,2)=bz(*)
;*******END BAD WAY

;********GOOD WAY***********

if keyword_set(B3) then bused= 'B3_gse' else bused = 'Bexp'
print,'B info used is ',bused

store_data,'sometimes',data = {x:s}
interpolate,'sometimes',bused,'newB'
get_data,'newB',data = newB
bdir = newB.y

;******END GOOD WAY*******

;now bx,by,bz are arrays of normalized B coords and bdir is a 
;n_points by 3 array of the B coords
;stop

;let's get the angle between B and the bin (indexed by energy,
;bin_index, and time point).
;How?  Let's get the cartesian coords of each 
;bin(energy,bin_index) then make this a n_points vector, then 
;use angl to take the angle between this (npoints,nenergy,nbin,3)
;vector and the B(npoints) vector


;okay, I figured out how to get a matrix of angles
;we know bdir is the (npoints,3) vec of B direction
;and we make dxyz the (nbins,3) vec for bin direction (for energy=i)
;we define

w=fltarr(n_points,3,n_energies,n_good)
angles=fltarr(n_energies,n_good,n_points)

;now we set every (k,*,i,*) point in w to the xyz coords of that bin
;for the respective energy and bin #

;also get the data/geom into a (n_energies,n_good,n_points) array

dvalues=fltarr(n_energies,n_good,n_points)

for i=0,n_energies-1 do begin
	sphere_to_cart,1.,d.theta(i,good_bins),d.phi(i,good_bins),$
		dx,dy,dz,vec=dxyz
	dxyz=transpose(dxyz)
;stop
	for j=0,n_good-1 do begin
		;for l=0,2 do begin
		;	w(*,l,i,j)=dxyz(j,l)
		;endfor
;stop
		
		;angles(i,j,*)=angl(w(*,*,i,j),bdir)   ;array for angles
;stop
		angles(i,j,*) = angl(bdir,reform(dxyz(j,*)))
		;if units eq 'counts' then $
		;	dvalues(i,j,*)=data(*).data(i,good_bins(j)) / data(*).geom(good_bins(j)) $
		;else $
			dvalues(i,j,*)=data(*).data(i,good_bins(j)) 
				 ;array for data

	endfor
endfor

;so now angles contains the angle between the bin and the B field for 
;everything, and is indexed by (energy,bin#,time point)

;so now our data is loaded in dvalues(n_energies,n_good,n_points) with
;corresponding angles loaded into angles(n_energies,n_good,n_points).


;now all we have to do is sort the data by angle and plot data/geom
;vs. angle vs. time! 


for i=0,n_energies-1 do begin
	istr=strcompress(string(i),/remove_all)
	
	q={x:data(*).time,v:fltarr(n_points,n_good+2),$
		y:fltarr(n_points,n_good+2),energy:data(0).energy(i),$
		ytitle:'',ztitle:'',$
		yrange:[0,180],ystyle:1,zlog:zlog,spec:1}
	
	for k=0,n_points-1 do begin

	;ok, here's where we should find each individual theta point and 
	;index q.v with a (n_points,n_good) subscript
	
		sort_idx=sort(angles(i,*,k)) ;gets angles for each time
			;and energy into monotonically increasing order
		q.v(k,1:n_good)=angles(i,sort_idx,k)
		q.y(k,1:n_good)=dvalues(i,sort_idx,k)
		
		;this is just to get bins at theta=0 and 180
		q.v(k,0)=0.
		q.v(k,n_good+1)=180.
		q.y(k,0)=dvalues(i,sort_idx(0),k)
		q.y(k,n_good+1)=dvalues(i,sort_idx(n_good-1),k)
		
	endfor
	
	
	q.ytitle=dat_name+'!c'+strcompress(string((data(0).energy(i,0)/1000.),$
		format='(g8.2)'),/remove_all)+'keV'
	if zlog ne 0 then q.ztitle=data(0).units_name+'!c(log)' else $
		q.ztitle=data(0).units_name
	
	store_data,'B_'+dat_name+'_'+istr,data=q
	
endfor

print,strcompress(string(systime(1)-starting_t),/remove_all)+' seconds'

end