; This routine calculates the normalized mean energy dissipaton rate.
; You should add the value of nu in param.pro
;
;
; Read in time series
;
lepsk='false'
default,file,'time_series.dat'
default,tmax,'infty'
default,tmin,'zero'
pc_read_ts,obj=obj,file=file
@data/index.pro
if (i_epsK ne 0) then lepsK='true'
;
tagnames=tag_names(obj)
lorms='false'
for i=0,n_tags(obj)-1 do begin
    if (tagnames[i] eq 'ORMS') then lorms='true'
end
;
; Some default values
;
default,t1,0.
default,t2,1e9
default,spec_corr,1.
;
; Read input params
;
@param
;
; Read in energy spectrum
;
power,k=k,v1='u',spec1=speu,tt=tt,i=n,/noplot
;
; If time span in time_series.dat is too large one
; might want to reduce what is shown.
;
if (tmax ne 'infty') then begin
    good_abs=where(tt lt tmax)
    tt=tt(good_abs)
end
if (tmin ne 'zero') then begin
    good_abs=where(tt gt tmin)
    tt=tt(good_abs)
end
;
; Where do we want to calculated averages
;
good=where((tt gt t1) and (tt lt t2))
goodts=where((obj.t gt t1) and (obj.t lt t2))
;
; Some auxillary variables
;
k1=double(k(1:*))
spe=speu(1:*,*)
;
; Effective wave number (due to discretisation errors)
;
kNy=max(k1)
keff6=(20.-30*cos(k1*!pi/kNy)+12*cos(2*k1*!pi/kNy)-2.*cos(3*k1*!pi/kNy))/!pi^6*kNy^6.
;
; Set some default values
;
default,hyper6,'true'
default,nu,5e-13
;
; Define some arrays
;
ns_all=(size(tt))[1]
Energy_all=fltarr(ns_all)
time_all=fltarr(ns_all)
um2_all=fltarr(ns_all)
L_all=fltarr(ns_all)
epsilon_all=fltarr(ns_all)
kd=fltarr(ns_all)
;
; Loop over all times
;
for i=0,ns_all-1 do begin
    Energy_all(i)=total(spe(*,i))
    time_all(i)=tt(i)
end
;
; Loop over all times
;
um2_all=2.*Energy_all/3.
urms_all=sqrt(2*Energy_all)
for i=0,ns_all-1 do begin
    L_all(i)=!pi*total(spe(*,i)/k1)/(2*um2_all(i))
    if (hyper6 eq 'true') then begin
        ikd=where(keff6*spe(*,i) eq max(keff6*spe(*,i)))
;        ikd=where(k^6.*spe(*,i) eq max(k^6.*spe(*,i)))
        kd(i)=k1(ikd)
        epsilon_all(i)=2.*nu*int_tabulated(k1,spe(*,i)*keff6)
    endif
    if (hyper6 ne 'true') then begin
        ikd=where(k1^2*spe(*,i) eq max(k1^2*spe(*,i)))
        kd(i)=k1(ikd)
        epsilon_all(i)=2.*nu*int_tabulated(k1,spe(*,i)*k1^2.)
    endif
end
if (lepsK eq 'true') then begin
    epsilon=obj.epsK
    epsilon=interpol(obj.epsK,obj.t,tt)
endif else begin
    epsilon=epsilon_all*spec_corr
end
Ceps_all=epsilon*L_all*um2_all^(-3./2)
;
; If orms is defined; find Re_lambda
;
if (lorms) then begin
;
; Find the Taylor microscale (lambda)
;
    ;lambda=sqrt(5)*obj.urms/obj.orms
;
; Find the Taylor microscale Reynolds number
;
    ;eps_temp=interpol(epsilon,tt,obj.t)/1.7
    ;nueps=eps_temp*(obj.orms)^(-2.)
    ;Re_lambda=obj.urms*lambda/nueps/sqrt(3)
    ;print,'lorms=',lorms
endif else begin
    ;print,'lorms=',lorms
end
;
; Print results
;
print,'-------------------------------------------'
print,'hyper6=',hyper6
print,'nu=',nu
print,'-------------------------------------------'
print,'mean(urms)=',mean(urms_all(good))
print,'mean(um2)=',mean(um2_all(good))
print,'mean(L)=',mean(L_all(good))
print,'mean(epsilon_spectrum)=',mean(epsilon_all(good))
if (lepsK eq 'true') then begin
    print,'mean(epsilon)=',mean(obj.epsK(goodts))
    print,'mean(epsilon_spectrum)/mean(epsilon)=',mean(epsilon_all(good))/mean(obj.epsK(goodts))
endif
print,'mean(Ceps)=',mean(Ceps_all(good))
print,'mean(kd)=',mean(kd(good))
print,'mean(Re_lambda1)=',(mean(kd(good))/1.5)^(2./3)*30.
;if (lorms) then print,'mean(Re_lambda2)=',mean(Re_lambda(good))    
print,'spec_corr=',spec_corr
;
; Plot results
;
!p.multi=[0,2,2]
plot,time_all,um2_all,title='!8u`!u2!n!6'
  oplot,time_all(good),um2_all(good),col=122
plot,time_all,L_all,title='!8L!6'
  oplot,time_all(good),L_all(good),col=122
if (lepsK eq 'true') then begin
    plot,obj.t,obj.epsK,title='!6epsilon'
    oplot,obj.t(goodts),obj.epsK(goodts),col=122
endif else begin
    print,'WARNING! Using epsilon found from the spectrum'
    plot,tt,epsilon,title='!6epsilon'
    oplot,tt(good),epsilon(good),col=122
end
plot,time_all,Ceps_all,title='!8C!6!deps!n'
  oplot,time_all(good),Ceps_all(good),col=122
!p.multi=[0,1,1]
;
END