;$Id: pkft_old.pro,v 1.1 2020/05/31 03:37:29 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=16,ysize=24,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
;  compute various length scales from the spectra and
;  calculate the instantaneous scaling exponents.
;
default,w,0
default,ihydro,0
@parameters
default,iread,0
if iread eq 0 then begin
  power,'_mag','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
  pc_read_param,o=param,/param2
  iread=1
  print,n
endif
;
!p.charsize=1.6
!x.margin=[8.8,.5]
!y.margin=[3.2,.5]
!p.multi=[0,2,2]
fo2='(f5.2)'
cwd,run
si2=1.4
;
eta=param.eta
kIM=t
kTM=t
kM=t
EM=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))
for it=0,nt-1 do begin
  kIM(it)=1./(total(spec1(*,it)*k1)/total(spec1(*,it)))
  kTM(it)=(total(spec1(*,it)*k^2)/total(spec1(*,it)))^.5
  kM(it)=(total(2*eta*spec1(*,it)*k^2)/eta^3)^.25
  EM(it)=total(spec1(*,it))
endfor
t0=t(0)
dt=t-t0+.4
cwd,run
;
;  define uniform ldt coordinate
;
ldt=alog(dt)
dlndt=.2
dlndt=.1
dlndt=.05
nlndt=nint(max(ldt)/dlndt)
lndt=dlndt*findgen(nlndt)
;
;  first plot
;
print
!x.title='!6'
!y.title='!8k!6!df!n/!8k!6!d1!n'
!p.title='!6'
yr=[4,3000]
yr=[1,1000]
yr=[.3,10]
tt=interpol(t,alog(dt),lndt)
lnEM=interpol(alog(EM),alog(dt),lndt)
lnkIM=interpol(alog(kIM),alog(dt),lndt)
pM=-deriv(lndt,lnEM)
qM=-deriv(lndt,lnkIM)
;
plot,lndt,lnEM,xtit='lndt',ytit='lnEM',ps=-1
plot,lndt,lnkIM,xtit='lndt',ytit='lnkIM',ps=-1
plot,lndt,pM,xtit='lndt',ytit='pM',ps=-1
plot,lndt,qM,xtit='lndt',ytit='qM',ps=-1
;
print
print,'mv idl.ps ~/tex/tina/Classes/fig/pkft_new_'+run+'.ps'
print,'mv idl.ps ~/tex/tina/heldecay/fig/pkft_new_'+run+'.ps'
print
;
print,'write savefile'
save,file='pq.sav',lndt,pM,qM
save,file='kE.sav',lndt,lnkIM,lnEM,tt,t,kIM,EM
;
!p.multi=0
!x.title=''      
!y.title=''
END