;$Id: pkft.pro,v 1.11 2014/07/19 08:21:04 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/pkft_short2048pm1b3_noforce.ps
;  mv idl.ps ~/tex/tina/Decay/fig/pkft_short2048pm1b4_noforce.ps
;  mv idl.ps ~/tex/tina/Decay/fig/pkft_short2048pm100b4_noforce.ps
;  mv idl.ps ~/tex/tina/Decay/fig/pkft_short2304pm1_kf60b_noforce.ps
;  mv idl.ps ~/tex/tina/Decay/fig/pkft_short2304pm100_kf60b_noforce.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
  iread=1
endif
eta=param.eta
nu=param.nu ;*t^1.5
;
!p.charsize=1.7
!x.margin=[8.8,.5]
!y.margin=[3.2,.2]
!p.multi=[0,2,2]
cwd,run
;
kI=t
kIK=t
kT=t
kD=t
kK=t
kM=t
RR=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(it)*total(k^2*spec1(*,it))
epsK=t & for it=0,nt-1 do epsK(it)=2.*nu*total(k^2*spec1(*,it))
for it=0,nt-1 do begin
  kI(it)=1./(total(spec2(*,it)*k1)/total(spec2(*,it)))
  kIK(it)=1./(total(spec1(*,it)*k1)/total(spec1(*,it)))
  kT(it)=(total(spec2(*,it)*k^2)/total(spec2(*,it)))^.5
  ;kD(it)=(total(2*(nu(it)*spec1(*,it)+eta*spec2(*,it))*k^2)/(nu(it)+eta)^3)^.25
  ;kK(it)=(total(2*nu(it)*spec1(*,it)*k^2)/nu(it)^3)^.25
  kD(it)=(total(2*(nu*spec1(*,it)+eta*spec2(*,it))*k^2)/(nu+eta)^3)^.25
  kK(it)=(total(2*nu*spec1(*,it)*k^2)/nu^3)^.25
  kM(it)=(total(2*eta*spec2(*,it)*k^2)/eta^3)^.25
  ;RR(it)=(total(2.*(spec1(*,it)+spec2(*,it))))^.5/((nu(it)+eta)*kT(it))
  RR(it)=(total(2.*(spec1(*,it)+spec2(*,it))))^.5/((nu+eta)*kT(it))
  EK(it)=total(spec1(*,it))
  EM(it)=total(spec2(*,it))
  HK(it)=total(spec1h(*,it))
  HM(it)=total(spec2h(*,it))
endfor
t0=t(0)
dt=t-t0+1.
!x.range=[1,250]
!x.range=[1,200]
cwd,run
;
!x.title='!8t!6/!8t!6!d0!n'
!y.title='!8k!6!df!n/!8k!6!d1!n'
plot_oo,dt,kI,yr=[1,200],tit='!6'+run
oplot,dt,kT,li=1
;oplot,dt,kD,li=2
oplot,dt,.01*kM,li=2
oplot,dt,.1*kK,li=3,col=55
xx=[8.,200.] & oplot,xx,40./xx^.5,col=122
;
;  slopes
;
plot,dt,deriv(alog(dt),alog(kI)),yst=0,ytit='slope'
oplot,dt,deriv(alog(dt),alog(kM)),li=1
;
!y.title='!8E!6!dM!n  and  !8E!6!dK!n'
plot_oo,dt,EM,yr=[2e-5,5e-2]
oplot,dt,EK,li=2
xx=[.8,200.] & oplot,xx,.04/xx^1.0,col=122
save,file='pkf0.sav',dt,kI,EM,EK
;
;  slopes
;
plot,dt,deriv(alog(dt),alog(EM)),yst=0,ytit='slope',yr=[-1.4,-.4]
oplot,dt,deriv(alog(dt),alog(EK)),li=2
;
;   plot,dt,(2*EM)^.5/(eta*kI)
;   plot,dt,(2*EK)^.5/(nu*kI)
;
!p.multi=0
stop
;
!p.multi=[0,1,2]
!y.title='!8R!6'
plot_oo,dt,RR,yr=yr1
oplot,dt,dt*0+Rm,col=122
;
!y.title='!8m!6'
slope1=-deriv(alog(dt),alog(kI))
slope2=-deriv(alog(dt),alog(EM))
slope3=-deriv(alog(dt),alog(EM+EK))
slopeh=-deriv(alog(dt),alog(HM))
plot_oi,dt,slope1,yr=yr2
oplot,dt,slope2,li=1
oplot,dt,slope3,col=122
oplot,dt,slopeh,col=55
oplot,dt,dt*0+mm
oplot,dt,dt*0+nn,col=122
oplot,dt,dt*0+ss,col=188,li=1
oplot,dt,dt*0+nn-mm,col=55
;
openw,1,'results.tmp'
fo='(f5.1,3f5.2,2x,a)'
print,Rm,mm,nn,nn-mm,run,fo=fo
printf,1,Rm,mm,nn,nn-mm,run,fo=fo
save,file='pkf.sav',dt,kI,kIK,kT,kD,kK,kM,EM,EK,HM,HK
close,1
spawn,'cat results.tmp'
spawn,'cat results.tmp >> ../idl/results.dat'
!p.multi=0
print,'plot,dt,kI'
END