;$Id: ppower_all.pro,v 1.71 2023/01/27 01:44:02 brandenb Exp $
@parameters
print,"$sed.csh data/param2.nml"
default,kin,0
default,gdone,1
default,mdone,0
default,hdone,1
default,hiidone,0
default,hijdone,0
default,giidone,0
default,gijdone,0
default,Stress_done,0
default,get_cstress_prefactor,1
default,divide_by_time,0
get_lhalf_factor_in_GW,lhalf_factor_in_GW
pc_read_param,/param2,obj=param2
pc_read_param,obj=param
if get_cstress_prefactor then cstress_prefactor=param2.cstress_prefactor
;dk=2.*!pi/param.Lxyz[0]
;dk=1.
;
default,lvar,1
default,iread,0
if iread eq 0 then begin
if lvar then begin
if kin then begin
  power,'_kin','hel_kin',k=k,spec1=spec1 ,spec2=spec2 ,i=n,tt=t,/noplot,/lks
  ;power,'_kin','_SCL',k=k,spec1=spec1 ,spec2=spec2 ,i=n,tt=t,/noplot,/lks
endif else begin
  ;power,'_mag','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot,/lks
  power,'_mag','_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot,/lks
  ;power,'_sp','hel_mag',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot,/lks
  power,'_kin','hel_kin',k=k,spec1=kine1,spec2=kine2,i=n,tt=t,/noplot,/lks
  ;power,'_TGW','_GWs',k=k,spec1=spec1,spec2=spec2,i=n,tt=t,/noplot,/lks
endelse
endif
if gdone then $
power,'_GWs','hel_GWs',k=k,spec1=grav1,spec2=grav2,i=n,tt=t,/noplot,/lks
if hdone then $
power,'_GWh','hel_GWh',k=k,spec1=grah1,spec2=grah2,i=n,tt=t,/noplot,/lks
;power,'_GWh','hel_GWh',k=k,spec1=grah1,spec2=grah2x,i=n,tt=t,/noplot,/lks
;power,'_SCL','_SCL',k=k,spec1=grah1,spec2=grah2,i=n,tt=t,/noplot,/lks
if hiidone then $
power,'_GWf','hel_GWf',k=k,spec1=graf1,spec2=graf2,i=n,tt=t,/noplot,/lks
if hijdone then $
power,'_GWg','hel_GWg',k=k,spec1=grag1,spec2=grag2,i=n,tt=t,/noplot,/lks
if giidone then $
power,'_GWd','hel_GWd',k=k,spec1=grad1,spec2=grad2,i=n,tt=t,/noplot,/lks
if gijdone then $
power,'_GWe','hel_GWe',k=k,spec1=grae1,spec2=grae2,i=n,tt=t,/noplot,/lks
if Stress_done then $
power,'_Str','hel_Str',k=k,spec1=stre1,spec2=stre2,i=n,tt=t,/noplot,/lks
;power,'_Str','hel_Str',k=k,spec1=stre1,spec2=stre2,i=n,tt=t,/noplot,/lks
;power,'_Str','hel_Str',k=k,spec1=grav1,spec2=grav2,i=n,tt=t,/noplot,/lks
iread=1
endif
;
;  time series
;
;spec1t=total(spec1,1)
if gdone then grav1t=total(grav1,1)
;
;  $Id: ppower_all.pro,v 1.71 2023/01/27 01:44:02 brandenb Exp $
;
;  This routine is used to view all power spectra
;  that were written during the run (controlled by dspec)
;  The plot range yr is set to a default value.
;  The time interval of plotting is "w" (default value)
;
;   5-oct-02/axel: written
;
spawn,'if (! -f param.pro) touch parameters.pro'
default,w,.2
default,t1,.0
default,t2,1e9
default,iover,0
default,istride,1
default,yr,[3e-12,3e-2]
ampl15=3e-10 & ampl4=1e-13
ampl15=3e-6 & ampl4=1e-8
dk=k[1]-k[0]
;
n=n_elements(t)
nt=n
good=where(t ge t1 and t le t2)
i1=min(good) > 0
i2=max(good) < (n-1)
;xr=dk*[k[1],max(k)]
xr=[k[1],max(k)]
;
for i=i1,i2,istride do begin
  if iover eq 0 or i eq 0 then begin
    plot_oo,xr,yr,tit='t='+str(t(i)),/nodata
    ;plot_oo,k,spec1(*,i),xr=xr,yr=yr,tit='t='+str(t(i))
    ;plot_oi,k,spec1(*,i),xr=xr,yr=[-1,1]*10,tit='t='+str(t(i))
    ;plot_oi,k,spec1(*,i),xr=xr,yr=yr,tit='t='+str(t(i))
  endif
  ;xx=[120.,2000.] & oplot,xx,1e-8/xx^.5
  circ_sym,1.3,0
  if divide_by_time then comp=1./t[i] else comp=1.
  if lvar then begin
  oplot,k,comp*spec1(*,i);,th=7
  if not kin then oplot,k,comp*kine1(*,i),col=55,th=5
  ;oplot,k,Ekin (*,i),col=122
  if kin then begin
;   oplot,k,abs(comp*spec2(*,i))/(2.*k),li=1
;   oplot,k,comp*spec2(*,i)/(2.*k),col=122,ps=8
;   oplot,k,-comp*spec2(*,i)/(2.*k),col=55,ps=8
  endif else begin
    ;oplot,k,comp*kine1(*,i),li=1,thick=5,col=55
    ;oplot,k,abs(comp*spec2(*,i))*k/2.,li=1,th=5
    ;oplot,k,comp*spec2(*,i)*k/2.,col=188,ps=8
    ;oplot,k,-comp*spec2(*,i)*k/2.,col=55,ps=8
    ;xx=[800,2e3] & oplot,xx,1e3/xx^2
  endelse
  endif
  ;
  if gdone then begin
  circ_sym,1.3,1
  oplot,k,comp*grav1(*,i),col=188,thick=6
; oplot,k,abs(comp*grav2(*,i)),li=1
; oplot,k,comp*grav2(*,i),col=122,ps=8
; oplot,k,-comp*grav2(*,i),col=55,ps=8
  loadct,6
  ;oplot,k,abs(grah2(*,i))/grah1(*,i),col=122,li=2
  ;-- oplot,k,grah2(*,i)/(grah1(*,i)),col=122
  loadct,5
  oplot,k,k*0.+1.
    ;xx=[10.,200.] & oplot,xx,1e-11/xx^2
    ;xx=[200,6e3] & oplot,xx,1e0/xx^4
    ;xx=[200,6e3] & oplot,xx,1e0/xx^6
    xx=[1.,100.] & oplot,xx,1e-8*xx^2
  endif
  ;
  if hdone then begin
  circ_sym,.8,1
  oplot,k,comp*grah1(*,i)*k^2,col=155,thick=8,li=2
  ;oplot,k,comp*grah1(*,i)*(10.^2+k^2)^1.5/k,col=122,thick=8,li=2
  oplot,k,comp*grah1(*,i),col=155,thick=1,li=2
; oplot,k,abs(comp*grah2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*grah2(*,i),col=122,ps=8
; oplot,k,-comp*grah2(*,i),col=55,ps=8
  endif
  ;
  if mdone then begin
  circ_sym,.8,1
  oplot,k,abs(comp*gram1(*,i))
  oplot,k,+comp*gram1(*,i),ps=8,col=122
  oplot,k,-comp*gram1(*,i),ps=8,col=55
; oplot,k,abs(comp*grah2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*grah2(*,i),col=122,ps=8
; oplot,k,-comp*grah2(*,i),col=55,ps=8
  endif
  ;
  if hiidone then begin
  circ_sym,.8,1
  oplot,k,k^2*comp*graf1(*,i),col=55,thick=2
; oplot,k,abs(comp*graf2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*graf2(*,i),col=122,ps=8
; oplot,k,-comp*graf2(*,i),col=55,ps=8
  endif
  ;
  if hijdone then begin
  circ_sym,.8,1
  oplot,k,k^2*comp*grag1(*,i),col=55,thick=2,li=2
; oplot,k,abs(comp*grag2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*grag2(*,i),col=122,ps=8
; oplot,k,-comp*grag2(*,i),col=55,ps=8
  endif
  if giidone then begin
  circ_sym,.8,1
  ;-- oplot,k,comp*grad1(*,i),col=122,thick=2
; oplot,k,abs(comp*grad2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*grad2(*,i),col=122,ps=8
; oplot,k,-comp*grad2(*,i),col=55,ps=8
  endif
  ;
  if gijdone then begin
  circ_sym,.8,1
  ;-- oplot,k,comp*grae1(*,i),col=122,thick=2,li=2
; oplot,k,abs(comp*grae2(*,i)),li=1,col=155,thick=4
; oplot,k,comp*grae2(*,i),col=122,ps=8
; oplot,k,-comp*grae2(*,i),col=55,ps=8
  endif
  ;
  if Stress_done then begin
  circ_sym,.8,1
  loadct,6
  oplot,k,comp*stre1(*,i),col=122
  ;oplot,k,abs(comp*stre2(*,i)),li=1,col=122,thick=4
  oplot,k,comp*stre2(*,i),col=122,ps=8
  oplot,k,-comp*stre2(*,i),col=55,ps=8
  loadct,5
  ;xx=[1.4,10] & oplot,xx,1e-4*xx^3
  endif
  ;
  ;kk=[1.,5.] & oplot,kk,.02/kk^8
  ;kk=[10.,100.] & oplot,kk,1e-33*kk^5
  ;kk=[1.,30.] & oplot,kk,1e-10*kk^0
  ;kk=[1.,30.] & oplot,kk,1e-10/kk^2
  ;kk=[30.,300.] & oplot,kk,1e-0/kk^3
  ;kk=[20.,170.] & oplot,kk,1e-1/kk^1.6667
  ;kk=[30.,200.] & oplot,kk,7e-0/kk^3.6667
  ;kk=[1.,50.] & oplot,kk,1e-13*kk^2
  ;kk=[1.,50.] & oplot,kk,1e-13*kk^4
  ;kk=[2.,10.] & oplot,kk,1e-4/kk^2
  ;kk=[10.,100.] & oplot,kk,.0002/kk^2
  ;kk=[8.,80.] & oplot,kk,1e-6/kk^(11./3.)
  print,i,t(i)
  wait,w
end
;
;good=where(t gt t1)
;it1=min(good)
;it2=max(good)
it1=i1
it2=i2
print,'it1,it2=',it1,it2
if lvar then begin
;
;  averaging interval
;
spec1m=total(spec1[*,it1:it2],2)/(it2-it1+1)
spec2m=total(spec2[*,it1:it2],2)/(it2-it1+1)
;
;  replace by just initial time (usually commented out)
;
;spec1m=spec1[*,0]
;spec2m=spec2[*,0]
;
tm=total(t[it1:it2])/(it2-it1+1)
if not kin then begin 
  kine1m=total(kine1[*,it1:it2],2)/(it2-it1+1)
  kine2m=total(kine2[*,it1:it2],2)/(it2-it1+1)
endif
endif
grav1m=total(grav1[*,it1:it2],2)/(it2-it1+1)
grav2m=total(grav2[*,it1:it2],2)/(it2-it1+1)
;
print,'save file?'
stop
spawn,'plot,t,grav1[1,*]'
;plot,t,grav1[1,*]
plot,t,spec1[1,*]
oplot,t[it1:it2],spec1[1,it1:it2],col=122
;
stop
if hdone then begin
  grah1m=total(grah1[*,it1:it2],2)/(it2-it1+1)
  grah2m=total(grah2[*,it1:it2],2)/(it2-it1+1)
  if kin then begin 
    save,file='specm.sav',k,spec1m,spec2m,grav1m,grav2m,grah1m,grah2m, $
      lhalf_factor_in_GW,cstress_prefactor,tm
  endif else begin
    save,file='specm.sav',k,spec1m,spec2m,grav1m,grav2m,grah1m,grah2m, $
      kine1m,kine2m,lhalf_factor_in_GW,cstress_prefactor,tm
  endelse
endif else begin
  if mdone then begin
    save,file='specm.sav',k,spec1m,spec2m,grav1m,grav2m,gram1m,gram2m, $
      kine1m,kine2m,lhalf_factor_in_GW,cstress_prefactor,tm
  endif else begin
    save,file='specm.sav',k,spec1m,spec2m,grav1m,grav2m, $
      kine1m,kine2m,lhalf_factor_in_GW,cstress_prefactor,tm
  endelse
endelse
spawn,'cvs add -kb specm.sav'
print,'sum for t>t1? t1,t2=',t[it1],t[it2]
use_ppower_all_saved=0
print,'use_ppower_all_saved=',use_ppower_all_saved
END