; $Id: pspec_comp_kan.pro,v 1.2 2004/05/24 09:29:30 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=20,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
if !d.name eq 'X' then begin
    window,0,xsize=600,ysize=800
end
;
;  Comparing 1-D spectra: simulations, wind tunnel (both decaying) and Kaneda (forced)
;  mv idl.ps ~/tex/hydro/hyper/fig/pspec_comp_kan.ps
;
power,v1='u',spec1=spec,k=k,i=i,/noplot,tt=tt
ns=(size(spec))[2]
del_itt=10
;
!x.title='!8k!6 (!7m!6!d3!n!8t!6)!u1/6!n'
!y.title='!8E!6(!8k,t!6) !8k!6!u5/3!n!7e!6!u-2/3!n (!8t!6/!7s!6)!u5/4!n'
!p.charsize=1.8
!p.multi=[0,1,2]
;
;  the value of eps512 must be consistent with that in hydro.pro
;  in tex/hydro/hyper/idl
;
eps512=0.000271407/1.7
nu=5e-13
fac=k^1.667/eps512^.667
i1=40 & i2=75 & i3=5
dt=tt-tt(0)
dt_exponent=1.25
tau=24. ;(the same as in urms^2/urms0^2 = [(t-t0)/tau]^1.25
;
; Plot 3D spectrum
;
k_comp=k*(nu*dt(i1))^0.1667
spec_comp=fac*spec(*,i1)*(dt(i1)/tau)^dt_exponent
plot_oo,k_comp,spec_comp,xr=[.018,3.4],yr=[.1,5]
for j=i1,i2,i3 do begin
  print,j,tt(j)-tt(0)
  k_comp=k*(nu*dt(j))^0.1667
  spec_comp=fac*spec(*,j)*(dt(j)/tau)^dt_exponent
  oplot,k_comp,spec_comp
end
;
siz=2.
xx=[.06,.4] & oplot,xx,2.0/xx^.03,li=2
xyouts,.34,1.4,'!9A!8k!6!u-0.03',siz=siz
;
; Plot 1D decay simulations
;
E1=0.*k & E_long=0.*k & nk=n_elements(k)
for i=1,nk-2 do begin
  k0 = k[i]
  E1(i)=integral(k(i:*),spec(i:*,i1)/k(i:*))
  E_long(i)=0.5*integral(k(i:*),(1-(k0/k[i:*])^2)*spec(i:*,i1)/k(i:*))
end
k_comp=k*(nu*dt(i1))^0.1667
spec_comp=fac*E_long*(dt(i1)/tau)^dt_exponent
plot_oo,k_comp,spec_comp,xr=[.0018,3.4],yr=[1e-1,0.7]
for j=i1,i2,i3 do begin
  for i=1,nk-2 do begin
      k0 = k[i]
      E1(i)=integral(k(i:*),spec(i:*,j)/k(i:*))
      E_long(i)=0.5*integral(k(i:*),(1-(k0/k[i:*])^2)*spec(i:*,j)/k(i:*))
  end
  k_comp=k*(nu*dt(j))^0.1667
  spec_comp=fac*E_long*(dt(j)/tau)^dt_exponent
  oplot,k_comp,spec_comp
end
;
; Plot 1D Kaneda 
;
a=rtable('$PENCIL_HOME/nils/hydro/hyper/idl/kaneda.dat',2,head=1)
k_kan=reform(a(0,*))
kkan=k_kan
E_kan=reform(a(1,*))
Ekan=E_kan*kkan^(-5./3)
E1=0.*kkan & E_long=0.*kkan & nk=n_elements(kkan)
for i=1,nk-2 do begin
  k0 = kkan[i]
  E1(i)=integral(kkan(i:*),Ekan(i:*)/kkan(i:*))
  E_long(i)=0.5*integral(kkan(i:*),(1-(k0/kkan[i:*])^2)*Ekan(i:*)/kkan(i:*))
end
oplot,kkan*7.,(E_long*kkan^(5./3))*1.6,li=2
;
; Plot 1D wind tunnel data
;
restore,file='$PENCIL_HOME/nils/hydro/hyper/idl/n900kol.sav'
k_wind=kint & speu_wind=Eint
epswind=4.146735103
Re_lambda=733
nuwind=1.63175E-05
oplot,k_wind*7.,(speu_wind*k_wind^(5./3))*1.,li=1
;
xx=0.03 & dx=0.1
yy=-0.7 & dy=-0.1
legend,xx,dx,10^(yy+dy*0.),0,'512!u3!n (1-dim)'
legend,xx,dx,10^(yy+dy*1.),1,'Wind (1-dim)'
legend,xx,dx,10^(yy+dy*2.),2,'Kaneda (1-dim)'
;
END