@parameters
;
default,ivar1,18
default,ivar2,62
default,ivar3,1
;
pc_read_param,/param2,obj=param
;
;  read through all snapshots
;
for ivar=ivar1,ivar2,ivar3 do begin
  pc_read_var,/bb,/trimall,obj=var,ivar=ivar
  t=var.t
  x=var.x
  y=var.y
  z=var.z
  iread=1
  nx=n_elements(x)
  ny=n_elements(y)
  nz=n_elements(z)
  ;
  ;  define the various fields: EBmodes
  ;
  nk=nx/2
  EEk=findgen(nk,nz)
  BBk=findgen(nk,nz)
  EBk=findgen(nk,nz)
  EBzk=findgen(nk,nz)
  BBzk=findgen(nk,nz)
  ;
  ;  loop over all heights
  ;
  s=1./float(nz)
  P_int=complexarr(nx,ny)
  R_int=complexarr(nx,ny)
  for iz=0,nz-1 do begin
    ;
    Bx=reform(var.bb(*,*,iz,0))
    By=reform(var.bb(*,*,iz,1))
    Bz=reform(var.bb(*,*,iz,2))
    ;
    P=complex(Bx,By)^2
    R=ebpol(P)
    ;
    P_int=P_int+s*P
    R_int=R_int+s*R
    ;
    E=float(R)
    B=imaginary(R)
    ;
    power2d_hel0,E,E,wav=k,spe=tmp,/quiet,/noplo & EEk(*,iz)=tmp
    power2d_hel0,B,B,wav=k,spe=tmp,/quiet,/noplo & BBk(*,iz)=tmp
    power2d_hel0,E,B,wav=k,spe=tmp,/quiet,/noplo & EBk(*,iz)=tmp
    power2d_hel0,E,Bz,wav=k,spe=tmp,/quiet,/noplo & EBzk(*,iz)=tmp
    power2d_hel0,B,Bz,wav=k,spe=tmp,/quiet,/noplo & BBzk(*,iz)=tmp
    ;
  endfor
  E_int=float(R_int)
  B_int=imaginary(R_int)
  power2d_hel0,E_int,E_int,wav=k,spe=EEk_int,/quiet,/noplo
  power2d_hel0,E_int,B_int,wav=k,spe=EBk_int,/quiet,/noplo
  power2d_hel0,B_int,B_int,wav=k,spe=BBk_int,/quiet,/noplo
  ;
  sum=EEk_int+BBk_int
  dif=EEk_int-BBk_int
  Sm=2*EBk_int/sum
  Tm=dif/sum
 ;
  cvar=str(ivar)
  save,file='sav/EBmodes'+cvar+'_new.sav',x,y,z,t,k,EEk,BBk,EBk,EBzk,BBzk,P_int,R_int,EEk_int,EBk_int,BBk_int
  save,file='sav/EBmodes'+cvar+'_int.sav',t,k,EEk_int,EBk_int,BBk_int
endfor
;
END