;pc_read_var,obj=var,/trimall,/magic,variables=['bb']
;
;  compute chi_B along line of sight
;
x=var.x
y=var.y
z=var.z
nx=n_elements(x)
ny=n_elements(y)
nz=n_elements(z)
ky=shift(findgen(ny)-(ny-1)/2,-(ny-1)/2) & kys=shift(ky,ny/2-1)
kz=shift(findgen(nz)-(nz-1)/2,-(nz-1)/2) & kzs=shift(kz,nz/2-1)
;
;  line of sight is now assumed to be the z direction
;
sigma=1.9
PPs=complexarr(nz,nx,ny)
PP1s=complexarr(nz,nx,ny)
PPdirs=complexarr(nz,nx,ny)
for m=0,ny-1 do begin
for l=0,nx-1 do begin
  Bx=reform(var.bb(l,m,*,0))
  By=reform(var.bb(l,m,*,1))
  BBdir=complex(Bx,By)
  Bperp=sqrt(Bx^2+By^2)
  chiB=atan(By,Bx)
  chiB1=atan(By/Bx)
  BB=Bperp^sigma*complex(cos(chiB),sin(chiB))
  BB1=Bperp^sigma*complex(cos(chiB1),sin(chiB1))
  ;
  PP=fft(BB,1)/nz
  PP1=fft(BB1,1)/nz
  PPdir=fft(BBdir,1)/nz
  ;
  PPs(*,l,m)=shift(PP,nz/2-1)
  PP1s(*,l,m)=shift(PP1,nz/2-1)
  PPdirs(*,l,m)=shift(PPdir,nz/2-1)
endfor
endfor
;
;  average over x and y, which corresponds here to directions 2 and 3
;
PPsm=total(total(PPs,3),2)/(nx*ny)
PP1sm=total(total(PP1s,3),2)/(nx*ny)
PPdirsm=total(total(PPdirs,3),2)/(nx*ny)
;
plot,kzs,abs(PPsm),xr=[-1,1]*8,yr=[0,.005]
oplot,kzs,abs(PP1sm),col=55,li=3
oplot,kzs,abs(PPdirsm)*10,col=122,li=2
;
plot,kzs,abs(PP1sm),xr=[-1,1]*8,yr=[0,.005]
END