;$Id: pkt_few_comp.pro,v 1.7 2022/11/24 06:07:48 brandenb Exp $
if !d.name eq 'PS' then begin
  ysize=11
  thick=1.6
  device,xsize=18,ysize=ysize,yoffset=3
  !p.charthick=thick & !p.thick=thick & !x.thick=thick & !y.thick=thick
end
;
siz=1.2
xstyle=0
@parameters
default,iread,0
default,ihydro,0
default,thick1,3
default,thick2,3
default,fourthird,1.
yr1_few_comp=[7e-10,.01]
xr1_few_comp=[.016,140.]
iread=0
if iread eq 0 then begin
  power,'_mag','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
  iread=1
endif
nt=n_elements(t)
loadct,5
;
default,siz,1.7
!p.multi=0
half='!s!u 1!n!r!s-!r!d 2!n'
!x.title='!8k!6!7n!6!dM!n(!8t!6)'
!y.title='!6(!8k!6!d0!n!7n!6)!u!7b!6!n !8E!6(!8k,t!6)  [!8B!6!s!u2!n!r!d0!n!8k!6!d0!n]'
!x.margin=[8.0,.3]
!y.margin=[3.2,.3]
!p.charsize=siz
tilde='!9!s!aA!n!r!6'
;
;  The multiplication of E(k) by k is because we plot versus k/k0.
;  We then devide by B^2 k0 to make the spectrum dimensionless.
;
i1=i1_few_comp
default,w,.1
default,k0,60.
default,B0,.1
default,vA0,1.
s=k0/(B0^2*k0)
xr=[.98,1.03*max(k)]/k0
tauA=1./(k0*vA0)
print,'tauA=',tauA
Ene0=total(spec1(*,i1))
Hel0=total(spec2(*,i1))
pc_read_param,/param2,o=param
;eta=param.eta
E=t
D=t
II=t
IIK=t
;
default,xstyle,1
if xstyle then begin
  plot_oo,xr,yr1_few,/nodata,xtickf='logticks_exp',ytickf='logticks_exp'
endif else begin
  plot_oo,xr1_few_comp,yr1_few_comp,/nodata
endelse
print,'blu=.35 & sred=.28'
;
icount=1
default,i1,0
default,iit,indgen(nt)
nit=n_elements(iit)
print,'it,t[it],t[it]*eta0*k0^2,t[it]*etat/xi^2,Brms/etat'
for i=i1,nit-1 do begin
  it=iit(i)
;
;  compute 1/k
;
k1=k & k1(0)=1. & k1=1./k1 & k1(0)=0.
EEM=total(spec1(*,it)) & xi=total(spec1(*,it)*k1)/EEM
Hel_unit=k0*vA0^2
two_xiE=2.*total(spec1(*,it)*k1)/Hel0
Hel1=total(abs(spec2(*,it)))/Hel0
Hel2=abs(total(spec2(*,it)))/Hel0
eps0=total(2*k^2*spec1(*,it))
Brms=sqrt(2.*EEM)
;
default,eta0,2e-5
default,sblu,.42
default,sred,.36
col=50+icount*10
col=170-icount*10
oplot,k*xi,s*spec1(*,it)*(k0*xi)^bet,col=col
;
si2=.9
xx=.5 & dx=.4 & yy=5e-5 & default,dy,2.5
legend,xx,dx,yy/dy^icount,0,col=col,siz=si2,'!8t!6 = '+hershey(1.1*t[it],0)
print,icount,yy/dy^icount,col
;
eta0=eta
default,t0,.1
default,rexp,-.43
etat=eta0*(t[it]>t0)^rexp
eta0eff=eta0*t0^rexp
tunitold=eta0*k0^2
tunit=1./(eta0eff*k0^2)
;etat=eta0
; uncomment for Table 2 in the paper
print,icount,it,t[it],t[it]*eta0*k0^2,t[it]*etat/xi^2,Brms/etat,eps0*xi^2/etat^2
print,icount,col
;print,icount,it,t[it],t[it]/tunit,EEM/Ene0,xi,two_xiE,Hel1,Hel2
fo="('  &',f7.2,' &',f8.4,' &',f6.1,' &',f7.2,' &',f6.2,' &',f6.2,' \\')"
fo_wold="('  &',2f9.2,' &',f8.4,' &',f6.1,' &',f7.2,' &',f6.2,' &',f6.2,' \\')"
; uncomment for Table 4 in the paper
;print,t[it]/tunit,EEM/Ene0,k0*xi,two_xiE,Hel1,Hel2,fo=fo
;print,t[it]/tunit,t[it]/tunitold,EEM/Ene0,k0*xi,two_xiE,Hel1,Hel2,fo=fo_wold
icount=icount+1
;wait,w
endfor
;
@postproc_pkt_few_comp
;
;  diagnostic output
;
print,'B0/[B]=',B0/eta0eff
print,'Brms/[B]=',vA0/eta0eff
;
cwd,run
print,'$mv idl.ps ~/tex/mhd/HallNonhel/fig/pkt_few_comp_'+run+'.eps'
!p.multi=0
END