;$Id: pkft_new.pro,v 1.12 2019/12/27 10:15:00 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
  if ihydro eq 1 then begin
    power,'_kin','hel_kin',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
    print,n
  endif else begin
    power,'_kin','_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot
    print,n
  endelse
  pc_read_param,o=param,/param2
  iread=1
endif
;
!p.charsize=3.3
!x.margin=[8.8,.5]
!y.margin=[3.2,.5]
!p.multi=[0,3,3]
fo2='(f5.2)'
cwd,run
si2=1.4
;
if ihydro eq 0 then begin
  eta=param.eta
  kIM=t
  kTM=t
  kM=t
  EM=t
endif
nu=param.nu
kIK=t
kTK=t
kK=t
EK=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))
for it=0,nt-1 do begin
  if ihydro eq 0 then begin
    kIM(it)=1./(total(spec2(*,it)*k1)/total(spec2(*,it)))
    kTM(it)=(total(spec2(*,it)*k^2)/total(spec2(*,it)))^.5
    kM(it)=(total(2*eta*spec2(*,it)*k^2)/eta^3)^.25
    EM(it)=total(spec2(*,it))
  endif
  kIK(it)=1./(total(spec1(*,it)*k1)/total(spec1(*,it)))
  kTK(it)=(total(spec1(*,it)*k^2)/total(spec1(*,it)))^.5
  kK(it)=(total(2*nu*spec1(*,it)*k^2)/nu^3)^.25
  EK(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]
default,xr1_pkft_new,[1,200]
default,yr1_pkft_new,[.3,10]
plot_oo,dt,kIK,xr=xr1_pkft_new,yr=yr1_pkft_new
lnkIK=interpol(alog(kIK),alog(dt),lndt)
tt=interpol(t,alog(dt),lndt)
qK=-deriv(lndt,lnkIK)
;oplot,exp(lndt),exp(lnkIK),col=55,ps=5
if ihydro eq 0 then begin
  oplot,dt,kIM,col=122,ps=-1
  oplot,dt,kTM,col=122,li=1
  oplot,dt,kM,col=122,li=2
  lnkIM=interpol(alog(kIM),alog(dt),lndt)
  qM=-deriv(lndt,lnkIM)
  ;oplot,exp(lndt),exp(lnkIM),col=122,ps=5
endif
oplot,dt,kIK,col=55
oplot,dt,kTK,col=55,li=1
;oplot,dt,kK,col=55,li=2
print,'first plot: kM (red), kK (blue)'
;xx=[3.,30.] & oplot,xx,1e2/xx^.6667
;xx=[3.,30.] & oplot,xx,3e0/xx^.5
;xx=[3.,30.] & oplot,xx,3e0/xx^.2
; & oplot,xx,3e0/xx^.2
xx=xr1_pkft_new
default,xfac1_pkft_new,1e2
default,xexp1_pkft_new,.3
oplot,xx,xfac1_pkft_new/xx^xexp1_pkft_new,ps=-1
xyouts,x1_pkft_new,y1_pkft_new,'!6slope= '+str(xexp1_pkft_new,fo=fo2),siz=si2
;
;  second plot
;
!x.title='!8t!6/!8t!6!d0!n'
!y.title='!8E!6!di!n'
!p.title='!6'+run
default,xr2_pkft_new,[1,200]
default,yr2_pkft_new,[3e-4,.1]
plot_oo,dt,EK,xr=xr2_pkft_new,yr=yr2_pkft_new,/nodata
lnEK=interpol(alog(EK),alog(dt),lndt)
pK=-deriv(lndt,lnEK)
if ihydro eq 0 then begin
  oplot,dt,EM,col=122,ps=-1
  lnEM=interpol(alog(EM),alog(dt),lndt)
  pM=-deriv(lndt,lnEM)
endif
oplot,dt,EK,col=55,ps=-1
print,'second plot: EM (red), EK (blue)'
;xx=[3.,30.] & oplot,xx,3e-2/xx^.6667
;xx=[3.,30.] & oplot,xx,3e-2/xx^.2
xx=xr2_pkft_new
default,xfac2_pkft_new,1e2
default,xexp2_pkft_new,.3
oplot,xx,xfac2_pkft_new/xx^xexp2_pkft_new,ps=-1
xyouts,x2_pkft_new,y2_pkft_new,'!6slope= '+str(xexp2_pkft_new,fo=fo2),siz=si2
;
;  third plot
;
!x.title='!8E!6!dM!n'
if ihydro eq 1 then !x.title='!8E!6!dK!n'
!y.title='!8k!6!df!n/!8k!6!d1!n'
!p.title='!6'
default,xr3_pkft_new,[5e-3,.14]
default,yr3_pkft_new,[1,140.]
plot_oo,EK,kIK,xr=xr3_pkft_new,yr=yr3_pkft_new,/nodata
b1K=deriv(lnkIK,lnEK)
if ihydro eq 0 then begin
  oplot,EM,kIM,col=122,ps=-1
  oplot,EM,kIK,col=55
  b1M=deriv(lnkIM,lnEM)
endif else begin
  oplot,EK,kIK,col=55,ps=-1
endelse
print,'third plot: kM vs EM'
xx=[.01,.1]
xx=xr3_pkft_new
default,xfac3_pkft_new,1e2
default,xexp3_pkft_new,1e2
oplot,xx,xfac3_pkft_new*xx^xexp3_pkft_new
xyouts,x3_pkft_new,y3_pkft_new,'!6slope= '+str(xexp3_pkft_new,fo=fo2),siz=si2
;
xx=[0,2]
circ_sym,1.3,1
;
;  straight p vs q
;
!y.title='!8q!6'
!x.title='!8p!6'
;plot,xx,2.*(1.-xx),xr=[0,1],yr=[0,2]
plot,2.*(1.-xx),xx,xr=[0,2],yr=[-.8,1]
bet=0. & oplot,xx,xx/(1.+bet) ;(beta=0 line)
bet=1. & oplot,xx,xx/(1.+bet) ;(beta=1 line)
bet=2. & oplot,xx,xx/(1.+bet) ;(beta=2 line)
bet=3. & oplot,xx,xx/(1.+bet) ;(beta=3 line)
for ilndt=1,nlndt-1 do begin
  oplot,pK[0:ilndt],qK[0:ilndt],ps=8,col=122
  oplot,pK[ilndt-1:ilndt],qK[ilndt-1:ilndt],ps=8,col=55
  if ihydro eq 0 then oplot,pM[0:ilndt],qM[0:ilndt],ps=8,col=55
  wait,w
endfor
;
;  pb vs q
;
!y.title='!8q!6'
!x.title='!8p!6!u(!7b!6)!n'
pbK=2.*(b1K+1.)/(b1K+3.)
;plot,xx,2.*(1.-xx),xr=[0,1.5],yr=[0,2]
plot,2.*(1.-xx),xx,yr=[-.5,1.5],xr=[0,2.4]
bet=0. & oplot,xx,xx/(1.+bet) ;(beta=0 line)
bet=1. & oplot,xx,xx/(1.+bet) ;(beta=1 line)
bet=3. & oplot,xx,xx/(1.+bet) ;(beta=3 line)
for ilndt=1,nlndt-1 do begin
  oplot,pbK[0:ilndt],qK[0:ilndt],ps=8,col=122
  if ihydro eq 0 then begin
    pbM=2.*(b1M+1.)/(b1M+3.)
    oplot,pbM[0:ilndt],qM[0:ilndt],ps=8,col=55
  endif
  wait,w
endfor
;
;  p vs qb
;
!y.title='!8q!6!u(!7b!6)!n'
!x.title='!8p!6'
qbK=2./(b1K+3.)
;plot,xx,2.*(1.-xx),xr=[0,1],yr=[0,2]
plot,2.*(1.-xx),xx,xr=[0,2.4],yr=[-.8,1]
bet=0. & oplot,xx,xx/(1.+bet) ;(beta=0 line)
bet=1. & oplot,xx,xx/(1.+bet) ;(beta=1 line)
bet=3. & oplot,xx,xx/(1.+bet) ;(beta=3 line)
for ilndt=1,nlndt-1 do begin
  oplot,pK[0:ilndt],qbK[0:ilndt],ps=8,col=122
  if ihydro eq 0 then begin
    qbM=2./(b1M+3.)
    oplot,pM[0:ilndt],qbM[0:ilndt],ps=8,col=55
  endif
  wait,w
endfor
;
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'
if ihydro eq 0 then begin
  save,file='pq.sav',lndt,pK,qK,pbK,qbK,b1K,pM,qM,pbM,qbM,b1M,ihydro
  save,file='kE.sav',lndt,lnkIK,lnEK,lnkIM,lnEM,ihydro,tt,t,kIK,kIM,EK,EM
endif else begin
  save,file='pq.sav',lndt,pK,qK,pbK,qbK,b1K,ihydro
  save,file='kE.sav',lndt,lnkIK,lnEK,ihydro,tt,t,kIK,kIM,EK,EM
endelse
;
;  plot Rey and Rm
;
Re=sqrt(2.*EK)/(nu*kIK)
plot,Re,ytit='Re'
;
if ihydro eq 0 then begin
  Re2=sqrt(2.*EK)/(nu*kIM)
  plot,Re2,ytit='Re2'
endif
;
!p.multi=0
!x.title=''      
!y.title=''
plot_oo,t,1./kIM,xr=[.01,500],yr=[.01,.1]
;
END