;$Id: pspec_sat_B10.pro,v 1.5 2021/01/16 20:41:50 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
;
@parameters
!p.charsize=1.6
!x.margin=[8.2,0.5]
!y.margin=[3.1,0.4]
;!p.multi=[0,2,2]
siz=1.8
si2=1.6
;
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.
ilast=904
;
!x.title='!8k!6'
!y.title='!8E!6!dM!n(!8k!6)  and  !8E!6!dGW!n(!8k!6)'
xr=minmax(k[1:*])
yr=[3e-18,3e-4]
plot_oo,xr,yr,/nodata
i=isat
;i=ilast
kGW1=total(spec1(1:*,i)/k[1:*])/total(spec1(1:*,i))
print,'kGW1=',kGW1
xyouts,7000.,1e-6,'!8E!6!dM!n(!8k!6)',siz=siz,col=122
oplot,k[1:*],s*spec1(1:*,i),col=122 & print,'t(mag)=',t[i]
;
xx=[310.,2200.] & oplot,xx,1e-13*(xx/100.)^5
xyouts,500.,1e-8,'!9A!8k!6!u5!n',siz=siz
;
for j=0,n_elements(jjj)-1 do begin
  oplot,k[1:*],s*spec1(1:*,jjj[j]),col=122,li=1
  print,'i,t=',j,t[jjj[j]]
endfor
;
i=ilast
oplot,k[1:*],sGW*grav1(1:*,i),col=55
oplot,k[1:*],sGW*grah1(1:*,i)*k[1:*]^2,col=55,li=2
for i=ilast-1,ilast do oplot,k[1:*],sGW*grav1(1:*,i),col=55
for i=ilast-1,ilast do oplot,k[1:*],sGW*grah1(1:*,i)*k[1:*]^2,col=55,li=1
print,'t(GW)=',t[ilast]
xx=[200.,3000.] & oplot,xx,3e-13/(xx/100.)^.5 & xyouts,450.,3e-15,'!9A!8k!6!u-0.5!n',siz=siz
xyouts,7000.,4e-13,'!8E!6!dGW!n(!8k!6)',siz=siz,col=55
xyouts,760.,6e-12,'!8k!6!u2!nSp(!8h!6)/6',siz=si2,col=55
;
mu=param.mu5_const
lam=param2.lambda5
eta=param2.eta
lameta2=lam*eta^2
;oplot,xr,xr*0+16.*mu*eta^2,li=1,thick=6
;oplot,xr,     16.*mu*eta^2*(mu/xr)^2,li=1
oplot,xr,xr*0+1.0*mu/lam,li=3
oplot,[1,1]*mu/2.,yr,li=2
;oplot,[1,1]*mu,yr,li=1
;
cwd,run
print,'$mv idl.ps ~/tex/rei/GW/fig/pspec_sat_'+run+'.ps'
!p.multi=0
END