;$Id: pkt2_comp.pro,v 1.8 2020/06/20 05:24:04 brandenb Exp $
if !d.name eq 'PS' then begin
  ysize=10
  default,thick,4
  device,xsize=18,ysize=ysize,yoffset=3
  !p.charthick=thick & !p.thick=thick & !x.thick=thick & !y.thick=thick
end
;
siz=1.7
xstyle=0
@parameters
default,iread,0
default,ihydro,0
default,thick1,3
default,thick2,3
default,fourthird,1.
  iread=0
if iread eq 0 then begin
  power,'_mag','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
help,t
  power,'_Lor','hel_Lor',k=k,spec1=spec1h,spec2=spec2h,i=n,tt=t,/noplot
help,t
  power,'_2Lor','hel_2Lor',k=k,spec1=spec1c,spec2=spec2c,i=n,tt=t,/noplot
help,t
  iread=1
endif
nt=n_elements(t)
;
default,siz,1.7
!p.multi=0
half='!s!u 1!n!r!s-!r!d 2!n'
!y.title='!8E!6!dM!n!6(!8k,t!6)'
!x.title='!8k!6/!8k!6!d0!n'
!y.title='!6'
!x.margin=[4.6,.3]
!y.margin=[3.2,.3]
!p.charsize=siz
tilde='!9!s!aA!n!r!6'
;
default,w,.1
default,k0,60.
default,vA0,1.
s=1.
xr=[.98,1.03*max(k)]/k0
tauA=1./(k0*vA0)
print,'tauA=',tauA
pc_read_param,/param2,o=param
;eta=param.eta
E=t
D=t
II=t
IIK=t
;
default,helical,0
default,xstyle,1
if xstyle then begin
  plot_oo,xr,yr1_comp,/nodata,xtickf='logticks_exp',ytickf='logticks_exp'
endif else begin
  plot_oo,xr,yr1_comp,/nodata
endelse
;
;  compute 1/k
;
k1=k & k1(0)=1. & k1=1./k1 & k1(0)=0.
EEM=total(spec1(*,it)) & xi=total(spec1(*,it)*k1)/EEM
eps=eta*total(2*k^2*spec1(*,it))
epsH=eta*total(2*k^2*spec2(*,it))
;
comp1=k^(7./3.)/eps^(2./3.)
comp1b=k^(10./3.)/epsH^(2./3.) ;(not correct)
comp1b=k^(9./3.)/epsH^(2./3.)
if helical then begin
  oplot,k/k0,comp1b*abs(s*spec2(*,it)*k/2.),thick=thick1,col=55,li=1
  oplot,k/k0,comp1b*s*spec2(*,it)*k/2.,thick=thick1,col=55
endif
oplot,k/k0,comp1*s*spec1(*,it),thick=thick2,col=122
oplot,[1,1]/(k0*xi),yr1_comp,li=3
;xyouts,.0070,6e-13,'!8E!6(!8k!6)',col=122
;xyouts,.0150,8e-14 ,'!8H!6(!8k!6)',col=55
;
;  plot II
;
comp2=k^(1./3.)/(eta*eps^(2./3.))
comp2b=k^(3./3.)/(eta*epsH^(2./3.))
oplot,k/k0,comp2*abs(spec1h(*,it)),li=1,thick=thick2
oplot,k/k0,-comp2*   spec1h(*,it), li=0,thick=thick2
;
if helical then begin
  loadct,6
  s3=k/2.
  oplot,k/k0,comp2b*abs(s3*spec2h(*,it)),li=1,thick=thick2,col=122
  oplot,k/k0,comp2b* (-s3)*spec2h(*,it), li=0,thick=thick2,col=122
endif
loadct,5
oplot,xr,xr*0+ym1_comp,li=3
;
siz=1.2
default,xx1,3.0 & dx1=.7*xx1
default,xx2,17. & dx2=.7*xx2
default,comp_label_fac,1.
default,ly1,.84
default,ly2,.59
legend,siz=siz,xx1,dx1,comp_label_fac*10^(ly1),col=122,0,'!8k!6!u7/3!n!8E!6(!8k!6)/!7e!6!u2/3!n'
legend,siz=siz,xx2,dx2,comp_label_fac*10^(ly1),        0,'!9+!8k!6!u1/3!n!8T!dE!n!6(!8k!6)/!7ge!6!u2/3!n'
if helical then begin
  legend,siz=siz,xx1,dx1,comp_label_fac*10^(ly2),col=55 ,0,'!8k!6!u4!n!8H!6(!8k!6)/!7e!6!s!dH!n!r!u2/3!n'      & loadct,6
  legend,siz=siz,xx2,dx2,comp_label_fac*10^(ly2),col=122,0,'!9+!8k!6!u!n!8T!dH!n!6(!8k!6)/!7ge!6!s!dH!n!r!u2/3!n' & loadct,5
  ;oplot,xr,xr*0+ym2_comp,li=3
endif
print,'it,t=',it,t[it]
;
default,xout,.56
default,yout,12.
default,label,'(a)'
xyouts,xout,yout,label
;
cwd,run
print,'$mv idl.ps ~/tex/mhd/HallMHD/fig/pkt_comp_'+run+'.ps'
!p.multi=0
END