@parameters
;
ldust=1
ldust=0
if ldust then dir='dust/' else dir='nodust/'
if ldust then prefix='dust_' else prefix=''
;
yr=[-1,1]
xr=[1,100]
!x.title='!8k !6[Mm!u-1!n]'
!y.title='!8c!dS!n!6(!8k!6)'
!y.title='!8c!dA!n!6(!8k!6)'
;
npz=ipz2-ipz1+1
tt=fltarr(npz)
Sm=fltarr(npz,nk)
Tm=fltarr(npz,nk)
EEk=fltarr(npz,nk)
EBk=fltarr(npz,nk)
BBk=fltarr(npz,nk)
plot,xr,yr,/nodata
;
;  perform R_tot integration
;
for ipz=ipz1,ipz2 do begin
  cpz=str(ipz)
  ;restore,dir+'R_int/'+prefix+'R_int_'+cpz+'.sav'
  restore,dir+'/'+prefix+'R_int_'+cpz+'.sav'
  if ipz eq ipz1 then R_tot=R_int else R_tot=R_tot+R_int
  print,'ipz=',ipz1,ipz,ipz2
endfor
;
;  compute pdf based on "old" fmax and nh values
;  that we also used for the individual slices.
;
E_tot=float(R_tot)
B_tot=imaginary(R_tot)
;
E_tot4=sample_2d(E_tot,4,4)
B_tot4=sample_2d(B_tot,4,4)
;
pdf,E_tot,xx,yyE,ndim=nh,fmin=-fmax,fmax=fmax
pdf,B_tot,xx,yyB,ndim=nh,fmin=-fmax,fmax=fmax
;
;  rms values
;
EBrms=sqrt(mean(E_tot^2+B_tot^2))
;
;  compute Tm and Sm
;
power2d_hel0,E_tot,E_tot,wav=k,spe=EEk_tot,/quiet,/noplo
power2d_hel0,E_tot,B_tot,wav=k,spe=EBk_tot,/quiet,/noplo
power2d_hel0,B_tot,B_tot,wav=k,spe=BBk_tot,/quiet,/noplo
;
sum=EEk_tot+BBk_tot
dif=EEk_tot-BBk_tot
Sm_tot=2*EBk_tot/sum
Tm_tot=dif/sum
;
save,file=dir+'EBmodes_EB_tot_histo.sav',xx,yyE,yyB,EBrms
save,file=dir+'EBmodes_EB_tot.sav',E_tot,B_tot
save,file=dir+'EBmodes_EB_tot4.sav',E_tot4,B_tot4
save,file=dir+'EBmodes_ST_tot.sav',k,Sm_tot,Tm_tot
END