;$Id: pspec_sat.pro,v 1.14 2021/03/10 09:45:16 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=11,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
;  previously used for 512_1e2_1e4_4e8_1em6aD
;  now adapted for 1024_1e2_1e4_4e8_1em6aD_k20
;
@parameters
!p.charsize=1.0
!x.margin=[8.2,0.5]
!y.margin=[3.1,0.4]
!p.multi=[0,2,2]
siz=1.
si2=.9
;
default,regime,2
print,"$sed.csh data/param2.nml"
pc_read_param,obj=param
pc_read_param,obj=param2,/param2
;
default,iread,0
if iread eq 0 then begin
  power,'_mag','_Str',k=k,spec1=spec1,spec2=stre1,i=n,tt=t,/noplot,/lks
  power,'_GWh','_GWs',k=k,spec1=grah1,spec2=grav1,i=n,tt=t,/noplot,/lks
  iread=1
endif
s=1d0/k[1] ;(scaling factor)
sGW=s/6.
;
;  to determine isat (or, better, imax) as the step when EEM is maximum;
;  compute EEM from spec1:
;
dk=k[1]-k[0]
EEM=dk*total(spec1,1)
EEGW=dk*total(grav1,1)
isat=(findex(max(EEM),EEM))[0]
print,'isat=',isat
;
;  plot
;
!x.title='!8k!6'
!y.title='!8E!6!dM!n(!8k!6)'
xr=minmax(k[1:*])
yr=[3e-27,1e-3]
plot_oo,xr,yr,/nodata
i=isat
oplot,k[1:*],s*spec1(1:*,i),li=0,thick=5
kGW1=total(spec1(1:*,i)/k[1:*])/total(spec1(1:*,i))
print,'kGW1(t=tmax)=',kGW1
xyouts,6000.,4e-13,'!8E!6!dM!n(!8k!6)',siz=siz
;for i=6000,ilast do oplot,k[1:*],s*spec1(1:*,i)
;? i=2000 & oplot,k[1:*],s*spec1(1:*,i),li=2
xx=[310.,2200.] & oplot,xx,1e-5/(xx/100.)^2
xyouts,700.,9e-10,'!9A!8k!6!u-2!n',siz=siz
print,'t(mag)=',t[i]
;
it1=0
nt=n_elements(t)
it2=nt-1
it3=100
;for it=isat,it2,it3 do print,it
for it=isat,it2,it3 do oplot,k[1:*],s*spec1(1:*,it),li=1
for it=it1,isat,it3 do oplot,k[1:*],s*spec1(1:*,it),li=2
i=it2-1
oplot,k[1:*],s*spec1(1:*,i),li=3,thick=5,col=122
;
yr=[3e-27,1e-3]
!y.title='!6Sp(!5T!6) (!8k,t!6)'
plot_oo,xr,yr,/nodata
for it=isat,it2,it3 do oplot,k[1:*],s*stre1(1:*,it),li=1
i=isat
oplot,k[1:*],s*stre1(1:*,i),li=0,thick=5
for it=it1,isat,it3 do oplot,k[1:*],s*stre1(1:*,it),li=2
i=it2-1
oplot,k[1:*],s*stre1(1:*,i),li=3,thick=5,col=122
;
yr=[3d-27,1e-10]
!y.title='!8E!6!dGW!n(!8k,t!6)'
plot_oo,xr,yr,/nodata
for it=isat,it2,it3 do oplot,k[1:*],s*grav1(1:*,it),li=1
for it=it1,isat,it3 do oplot,k[1:*],s*grav1(1:*,it),li=2
i=isat
oplot,k[1:*],s*grav1(1:*,i),li=0,thick=5
i=it2
oplot,k[1:*],s*grav1(1:*,i),li=3,thick=5,col=122
oplot,k[1:*],s*grah1(1:*,i)*k[1:*]^2,li=3,thick=5,col=55
xx=[13,40] & oplot,xx,1e-25*xx^2
;
yr=[3d-27,1e-10]
!y.title='!6Sp(!5h!6) (!8k,t!6)'
plot_oo,xr,yr,/nodata
for it=isat,it2,it3 do oplot,k[1:*],s*grah1(1:*,it),li=1
for it=it1,isat,it3 do oplot,k[1:*],s*grah1(1:*,it),li=2
i=isat
oplot,k[1:*],s*grah1(1:*,i),li=0,thick=5
i=it2
oplot,k[1:*],s*grah1(1:*,i),li=3,thick=5,col=122
xx=[13,40] & oplot,xx,1e-25*xx^2
;
cwd,run
print,'$mv idl.ps ~/tex/rei/GW/fig/pspec_sat_'+run+'.ps'
print,'$mv idl.ps ~/tex/notes/GW_regimeI/fig/pspec_sat_'+run+'.ps'
!p.multi=0
END