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