;+
;FUNCTION:	p_3d_new(dat,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins)
;INPUT:	
;	dat:	structure,	2d data structure filled by get_eesa_surv, get_eesa_burst, etc.
;KEYWORDS
;	ENERGY:	fltarr(2),	optional, min,max energy range for integration
;	ERANGE:	fltarr(2),	optional, min,max energy bin numbers for integration
;	EBINS:	bytarr(na),	optional, energy bins array for integration
;					0,1=exclude,include,  
;					na = dat.nenergy
;	ANGLE:	fltarr(2,2),	optional, angle range for integration
;				theta min,max (0,0),(1,0) -90<theta<90 
;				phi   min,max (0,1),(1,1)   0<phi<360 
;	ARANGE:	fltarr(2),	optional, min,max angle bin numbers for integration
;	BINS:	bytarr(nb),	optional, angle bins array for integration
;					0,1=exclude,include,  
;					nb = dat.ntheta
;	BINS:	bytarr(na,nb),	optional, energy/angle bins array for integration
;					0,1=exclude,include
;PURPOSE:
;	Returns the pressure tensor, [Pxx,Pyy,Pzz,Pxy,Pxz,Pyz], eV/cm^3 
;NOTES:	
;	Function normally called by "get_3dt" or "get_2dt" to
;	generate time series data for "tplot.pro".
;
;CREATED BY:
;	J.McFadden	00-2-24	
;LAST MODIFICATION:
;	J.McFadden	05-2-8		Fixed diagonalization
;	J.McFadden	06-2-23		changed the s/c pot calculation to the same as n_2d_new.pro
;	J.McFadden	09-4-29		added a diagonalization to Pxx and Pyy to maximize difference
;-
function p_3d_new,dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra

p3dxx = 0. & p3dyy = 0. & p3dzz = 0. & p3dxy = 0. & p3dxz = 0. & p3dyz = 0.

if dat2.valid eq 0 then begin
  dprint, 'Invalid Data'
  return, [p3dxx,p3dyy,p3dzz,p3dxy,p3dxz,p3dyz]
endif

mass  = dat2.mass
momen = m_3d_new(dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra)
flux=j_3d_new(dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra)
density=n_3d_new(dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra)

p3dxx = (momen(0)-mass*flux(0)*flux(0)/density/1.e10)
p3dyy = (momen(1)-mass*flux(1)*flux(1)/density/1.e10)
p3dzz = (momen(2)-mass*flux(2)*flux(2)/density/1.e10)
p3dxy = (momen(3)-mass*flux(0)*flux(1)/density/1.e10)
p3dxz = (momen(4)-mass*flux(0)*flux(2)/density/1.e10)
p3dyz = (momen(5)-mass*flux(1)*flux(2)/density/1.e10)

; Rotate the tensor about the magnetic field to diagonalize
; This should give a result that diagonalizes the pressure tensor about dat2.magf
; where magf is in s/c coordinates -- ie same as the dat2.theta and dat2.phi coordinates.

; First form the pressure tensor
	p = [[p3dxx,p3dxy,p3dxz],[p3dxy,p3dyy,p3dyz],[p3dxz,p3dyz,p3dzz]]

if finite(total(dat2.magf)) && (total(dat2.magf*dat2.magf) gt 0.) then begin
	bx=dat2.magf(0)
	by=dat2.magf(1)
	bz=dat2.magf(2)
  ; Rotate p about Z-axis by the angle between X and the projection of B on the XY plane
	ph=atan(by,bx)
	rot_ph=([[cos(ph),-sin(ph),0],[sin(ph),cos(ph),0],[0,0,1]])
	p = rot_ph#p#transpose(rot_ph)
  ; Then rotate p about Y-axis by the angle between Bz and B 
	th=!pi/2.-atan(bz,(bx^2+by^2)^.5)
	rot_th=[[cos(th),0,sin(th)],[0,1,0],[-sin(th),0,cos(th)]]
	p = rot_th#p#transpose(rot_th)
  ; Finally diagonalize Pxx and Pyy
	l1 = (p[0,0]+p[1,1] + (p[0,0]^2+p[1,1]^2-2.*p[0,0]*p[1,1]+4.*p[0,1]^2)^.5)/2.
	l2 = (p[0,0]+p[1,1] - (p[0,0]^2+p[1,1]^2-2.*p[0,0]*p[1,1]+4.*p[0,1]^2)^.5)/2.
	if l1 ne l2 then ph=acos(((p[0,0]*l1-p[1,1]*l2)/(l1^2-l2^2))^.5) else ph=0.	; ph is the rotation angle to diagonalize
	rot_ph=([[cos(ph),-sin(ph),0.],[sin(ph),cos(ph),0.],[0.,0.,1.]])
	p = rot_ph#p#transpose(rot_ph)
endif
	return, [p(0,0),p(1,1),p(2,2),p(0,1),p(0,2),p(1,2)]

;	Pressure is in units of eV/cm**3,  [p3dxx,p3dyy,p3dzz,p3dxy,p3dxz,p3dyz]

end

