;$Id: pkt.pro,v 1.8 2020/06/13 19:36:51 brandenb Exp $
if !d.name eq 'PS' then begin
  default,ysize,8
  default,thick,4
  device,xsize=18,ysize=ysize,yoffset=3
  !p.charthick=thick & !p.thick=thick & !x.thick=thick & !y.thick=thick
end
;
@parameters
default,iread,0
default,ihydro,0
default,thick1,3
default,thick2,3
default,fourthird,1.
if iread eq 0 then begin
  power,'_mag','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
  power,'_Lor','hel_Lor',k=k,spec1=spec1h,spec2=spec2h,i=n,tt=t,/noplot
  iread=1
endif
nt=n_elements(t)
;
default,siz,1.7
!p.multi=[0,2,1]
half='!s!u 1!n!r!s-!r!d 2!n'
!x.title='!8k!6/!8k!6!d*!n'
!y.title='!8E!6(!8k,t!6)'
!x.margin=[7.4,1.3]
!y.margin=[3.2,.4]
!p.charsize=siz
tilde='!9!s!aA!n!r!6'
;
;yr1=[2e-14,5e-4]
;yr2=[2e-14,5e-2]
;
default,w,.1
default,k0,60.
default,vA0,1.
s=1./(vA0^2*k0) & default,yr1,[2e-9,2e-2] & default,xr,[.9,576]/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,xstyle,1
if xstyle then begin
  plot_oo,k/k0,s*spec2(*,0),xr=xr,yr=yr1,/nodata,xtickf='logticks_exp',ytickf='logticks_exp'
endif else begin
  plot_oo,k/k0,s*spec2(*,0),xr=xr,yr=yr1,/nodata
endelse
;
iit=indgen(nt)
li1=indgen(nt)
li2=indgen(nt)
it0=iit[0]
nit=n_elements(iit)
for i=0,nit-1 do begin
  it=iit(i)
  oplot,k/k0,s*spec1(*,it),li=li1[i],thick=thick2,col=55
  oplot,k/k0,s*spec2(*,it)*k/2.,li=li2[i],thick=thick1,col=122
  ;Rm=sqrt(2.*mean(spec2(*,it)))/(eta*k0)
  ;print,'i,it,t[it],Rm=',i,it,t[it],Rm
  ;
  ;  compute 1/k
  ;
  k1=k & k1(0)=1. & k1=1./k1 & k1(0)=0.
  EEM=total(spec2(*,it)) & xi=total(spec2(*,it)*k1)/EEM
  wait,w
endfor
;close,1
;if ihydro eq 0 then begin
;  it=n-2 & oplot,k/k0,s*spec1(*,it)/fourthird,li=li1[nit-1],thick=8,col=55
;  it=n-2 & oplot,k/k0,s*spec2(*,it),li=li2[nit-1],thick=8,col=122
;endif else begin
;  oplot,k/k0,s*spec1(*,it),li=0,thick=8,col=55
;endelse
;
siz=1.7
default,xlab,5.
default,ylab,1e-6
;if label ne '' then xyouts,xlab,ylab,label,siz=siz
;
;  plot II
;
default,xstyle,1
if xstyle then begin
  plot_oo,k/k0,s*spec2(*,0),xr=xr,yr=yr2,/nodata,xtickf='logticks_exp',ytickf='logticks_exp'
endif else begin
  plot_oo,k/k0,s*spec2(*,0),xr=xr,yr=yr2,/nodata
endelse
;
iit=indgen(nt)
li1=indgen(nt)
li2=indgen(nt)
it0=iit[0]
nit=n_elements(iit)
for i=0,nit-1 do begin
  it=iit(i)
  oplot,k/k0,abs(s*spec1h(*,it)),li=1,thick=thick2
  oplot,k/k0,    s*spec1h(*,it), li=0,thick=thick2
  wait,w
endfor
;
;it=n-2 & oplot,k/k0,s*spec1(*,it)*k/2.,li=2,thick=8
;
;  overplot slope
;
;print,'@postproc_pkt'
;@postproc_pkt
;xx=[10.,40.] & oplot,xx/k0,1e-10*xx^2,col=122
;xx=[10.,40.] & oplot,xx/k0,1e-11*xx^4,col=122
;
save,file='spec.sav',t,k,spec1,spec2,spec1h,spec2h
spawn,'cvs add -kb spec.sav'
;
if texdir eq 'inflation-hel' then ext='_all' else ext=''
;
cwd,run
default,texdir,'heldecay'
print,'$mv idl.ps ~/tex/mhd/HallMHD/fig/pkt_'+run+ext+'.ps'
!p.multi=0
END