;+
;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 error due to statistics, dp, eV/cm^3, corrects for spacecraft potential if dat.sc_pot exists, diagonalizes 
;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 dp_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)
dmomen = dm_3d_new(dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra)
dflux=dj_3d_new(dat2,ENERGY=en,ERANGE=er,EBINS=ebins,ANGLE=an,ARANGE=ar,BINS=bins,_extra=_extra)
ddensity=dn_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)

dp3dxx = ((dmomen[0])^2 + (2.*mass/1.e10*flux(0)*dflux(0)/density)^2. + (mass/1.e10*flux(0)*flux(0)*ddensity/density^2)^2)^.5
dp3dyy = ((dmomen[1])^2 + (2.*mass/1.e10*flux(1)*dflux(1)/density)^2. + (mass/1.e10*flux(1)*flux(1)*ddensity/density^2)^2)^.5
dp3dzz = ((dmomen[2])^2 + (2.*mass/1.e10*flux(2)*dflux(2)/density)^2. + (mass/1.e10*flux(2)*flux(2)*ddensity/density^2)^2)^.5
dp3dxy = ((dmomen[3])^2 + (2.*mass/1.e10*flux(0)*dflux(1)/density)^2. + (mass/1.e10*flux(0)*flux(1)*ddensity/density^2)^2)^.5
dp3dxz = ((dmomen[4])^2 + (2.*mass/1.e10*flux(0)*dflux(2)/density)^2. + (mass/1.e10*flux(0)*flux(2)*ddensity/density^2)^2)^.5
dp3dyz = ((dmomen[5])^2 + (2.*mass/1.e10*flux(1)*dflux(2)/density)^2. + (mass/1.e10*flux(1)*flux(2)*ddensity/density^2)^2)^.5

; 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]]
	dp = [[dp3dxx,dp3dxy,dp3dxz],[dp3dxy,dp3dyy,dp3dyz],[dp3dxz,dp3dyz,dp3dzz]]

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 and dp 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)
	dp = rot_ph#dp#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)
	dp = rot_th#dp#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.]])
	dp = rot_ph#dp#transpose(rot_ph)
endif
	return, [dp(0,0),dp(1,1),dp(2,2),dp(0,1),dp(0,2),dp(1,2)]

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

end