;$Id: pkf3.pro,v 1.1 2012/02/04 19:27:55 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=12,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
;  mv idl.ps ~/tex/tina/Decay/fig/pkf.ps
;
@parameters
default,iread,0
if iread eq 0 then begin
  power,'_kin','_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot & print,n
  power,'hel_kin','hel_mag',k=k,spec1=spec1h,spec2=spec2h,i=n,tt=t,/noplot & print,n
  pc_read_param,o=param,/param2
  eta=param.eta
  nu=param.nu
  iread=1
endif
;
!p.charsize=1.7
!x.margin=[8.8,.5]
!y.margin=[3.2,.2]
!p.multi=[0,1,3]
;
kI=t
kT=t
kD=t
kM=t
EK=t
EM=t
HK=t
HM=t
nt=n_elements(t)
;
;  compute 1/k
;
k1=k
k1(0)=1.
k1=1./k1
k1(0)=0.
;
;  loop over all times
;
E=t & for it=0,nt-1 do E(it)=total(spec1(*,it))
epsK=t & for it=0,nt-1 do epsK(it)=2.*nu*total(k^2*spec1(*,it))
;
s=k^1.667
yr=[1e-3,1.]
nit=n_elements(iit3)
for i=0,nit-1 do begin
  it=iit3(i)
  plot_oo,k,s*spec2(*,it),li=0,thick=.5,xr=[1,256],yr=yr
  oplot,k,s*spec1(*,it)*k/2.,li=1,thick=.5
  kI(it)=1./(total(spec2(*,it)*k1)/total(spec2(*,it)))
  kT(it)=(total(spec2(*,it)*k^2)/total(spec2(*,it)))^.5
  kD(it)=(total(2*(nu*spec1(*,it)+eta*spec2(*,it))*k^2)/(nu+eta)^3)^.25
  kM(it)=(total(2*eta*spec2(*,it)*k^2)/eta^3)^.25
  oplot,[1,1]*kI(it),yr
  oplot,[1,1]*kT(it),yr,li=1
  oplot,[1,1]*kM(it),yr,li=2
  print,kM(it)
endfor
END