; This routine calculates the dissipaton (epsilon) in 
; different ways.
; You should add the value of nu in param.pro
;
!p.charsize=1.4
;
lepsK='false'
lepsM='false'
;
; Read index.pro
;
@data/index.pro
if (i_epsK ne 0) then lepsK='true'
if (i_epsM ne 0) then lepsM='true'
;
; Set some default values
;
default,hyper6,'true'
default,nu,5e-13
default,hydro,'true'
default,fak,1.
default,fakmag,1.
;
; Read input params
;
@param
;
; Read in energy spectrum
;
if (hydro eq 'true') then begin
    power,k=k,v1='u',spec1=spec1,tt=tt,i=n,/noplot
    k1=k(1:*)
    speu=spec1(1:*,*)
endif else begin
    power,k=k,v1='u',v2='b',spec1=spec1,spec2=spec2,tt=tt,i=n,/noplot
    k1=k(1:*)
    speu=spec1(1:*,*)
    speb=spec2(1:*,*)
end
;
; 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.
;
; Read in time series
;
pc_read_ts,obj=obj
;
; 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))
;
; Define some arrays
;
ns_all=(size(tt))[1]
time_all=fltarr(ns_all)
epsilon_kin=fltarr(ns_all)
epsilon_mag=fltarr(ns_all)
print,'-------------------------------------------'
;
; Loop over all times
;
for i=0,ns_all-1 do begin
    time_all(i)=tt(i)
    ;
    ; Pure hydro?
    ;
    if (hydro eq 'true') then begin
        if (hyper6 eq 'true') then begin
            epsilon_kin(i)=2.*nu*int_tabulated(k1,speu(*,i)*keff6)
        endif else begin
            epsilon_kin(i)=2.*nu*int_tabulated(k1,speu(*,i)*k1^2.)
        end
    ;
    ; MHD?
    ;
    endif else begin
        if (hyper6 eq 'true') then begin
            epsilon_kin(i)=2.*nu*int_tabulated(k1,speu(*,i)*keff6)
            epsilon_mag(i)=2.*nu*int_tabulated(k1,speb(*,i)*keff6)
        endif else begin
            epsilon_kin(i)=2.*nu*int_tabulated(k1,speu(*,i)*k1^2.)
            epsilon_mag(i)=2.*nu*int_tabulated(k1,speb(*,i)*k1^2.)
        end
    end
end
;
; Print results
;
print,'mean(epsilon_spectrum_kin)=',mean(epsilon_kin(good))
if (lepsK eq 'true') then begin
    print,'mean(epsilon_kin)=',mean(obj.epsK(goodts))
    print,'mean(epsilon_spectrum_kin)/mean(epsilon_kin)=', $
      mean(epsilon_kin(good))/mean(obj.epsK(goodts))
endif
;
; MHD case?
;
if (hydro ne 'true') then begin
    print,'-------------------------------------------'
    print,'mean(epsilon_spectrum_mag)=',mean(epsilon_mag(good))
    if (lepsM eq 'true') then begin
        print,'mean(epsilon_mag)=',mean(obj.epsM(goodts))
        print,'mean(epsilon_spectrum_mag)/mean(epsilon_mag)=', $
          mean(epsilon_mag(good))/mean(obj.epsM(goodts))
    endif
end
print,'-------------------------------------------'
print,'hyper6=',hyper6
print,'nu=',nu
print,'fak=',fak
print,'fakmag=',fakmag
print,'printing epsK=', lepsK
print,'printing epsM=', lepsM
print,'-------------------------------------------'
;
; Plot results
;
if (hydro eq 'true') then begin
    !p.multi=[0,1,1]
    plot,tt,epsilon_kin,li=2,title='!6epsilon_kin'
    oplot,tt(good),epsilon_kin(good),col=122,li=2	
    oplot,obj.t,obj.epsK*fak
    oplot,obj.t(goodts),obj.epsK(goodts)*fak,col=122
endif else begin
    !p.multi=[0,1,2]
    plot,tt,epsilon_kin,li=2,title='!6epsilon_kin'
    oplot,tt(good),epsilon_kin(good),col=122,li=2	
    oplot,obj.t,obj.epsK*fak
    oplot,obj.t(goodts),obj.epsK(goodts)*fak,col=122
    ;
    plot,tt,epsilon_mag,li=2,title='!6epsilon_mag'
    oplot,tt(good),epsilon_mag(good),col=122,li=2
    oplot,obj.t,obj.epsM*fakmag
    oplot,obj.t(goodts),obj.epsM(goodts)*fakmag,col=122
end
;
!p.multi=[0,1,1]
;
END