@parameters
;
;  sample to compare time series form different directories
;
print,"t1_bfield=2.8"
print,"t2_bfield=4.3"
;
default,iread_spec,0
if iread_spec eq 0 then begin
  power,'_mag','_mag',k=k,spec1=spec1,spec2=spec1,i=n,tt=t,/noplot,/lks
  power,'_GWs','_GWs',k=k,spec1=grav1,spec2=grav1,i=n,tt=t,/noplot,/lks
  power,'_GWh','_GWh',k=k,spec1=grah1,spec2=grah1,i=n,tt=t,/noplot,/lks
  ;power,'_kin','_GWh',k=k,spec1=kine1,spec2=grah1,i=n,tt=t,/noplot,/lks
  power,'_Str','hel_Str',k=k,spec1=stre1,spec2=stre2,i=n,tt=t,/noplot,/lks
  ;power,'hel_mag','hel_GWs',k=k,spec1=spec2,spec2=grav2,i=n,tt=t,/noplot,/lks
  ;power,'hel_kin','hel_GWh',k=k,spec1=kine2,spec2=grah2,i=n,tt=t,/noplot,/lks
  nt=n_elements(t)
  iread_spec=1
endif
;
;  for GW field
;
it1=i1
it2=i2
it2=nt-2
;
;  for magnetic field
;
good_bfield=where(t ge t1_bfield and t le t2_bfield)
jt1=min(good_bfield)
jt2=max(good_bfield)
;
print,'it1,it2=',it1,it2
print,'jt1,jt2=',jt1,jt2
tm=total(t[it1:it2])/(it2-it1+1)
spec1m=total(spec1[*,jt1:jt2],2)/(jt2-jt1+1)
;kine1m=total(kine1[*,jt1:jt2],2)/(jt2-jt1+1)
grav1m=total(grav1[*,it1:it2],2)/(it2-it1+1)
grah1m=total(grah1[*,it1:it2],2)/(it2-it1+1)
;spec2m=total(spec2[*,jt1:jt2],2)/(jt2-jt1+1)
;kine2m=total(kine2[*,jt1:jt2],2)/(jt2-jt1+1)
;grav2m=total(grav2[*,it1:it2],2)/(it2-it1+1)
;grah2m=total(grah2[*,it1:it2],2)/(it2-it1+1)
;
spawn,'cvs add -kb specm.sav'
spawn,'plot,t,grav1[1,*]'
;
!p.multi=[0,1,2]
yr=minmax(grav1[1,*])
plot,t,grav1[1,*]
oplot,[1,1]*t[it1],yr,col=122
oplot,[1,1]*t[it2],yr,col=122
;
ikm=5
;yr=minmax(kine1[ikm,*])
;plot,t,kine1[ikm,*],xr=[t[0],t[2*jt2<it2]]
yr=minmax(spec1[ikm,*])
plot,t,spec1[ikm,*],xr=[t[0],t[2*jt2<it2]]
oplot,[1,1]*t[jt1],yr,col=122
oplot,[1,1]*t[jt2],yr,col=122
;
!p.multi=0
;
save,file='specm.sav',k,spec1m,grav1m,grah1m, $ ;kine1m, $
     ;spec2m,grav2m,grah2m,kine2m, $
     tm
     ;lhalf_factor_in_GW,cstress_prefactor,tm
;
print,'sum for t>t1? t1,t2=',t[it1],t[it2]
use_ppower_all_saved=0
;
spawn,'cvs add -kb specm.sav'
END