;$Id: pstress_etc.pro,v 1.12 2022/01/29 03:34:28 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=19,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
  col0=0
endif else begin
  col0=255
end
;
;  this is to be run in: 2048_1e2_1e4_4e8_1em6aD_k10_bes
;
!p.charsize=1.7
!x.margin=[10.5,1.8]
!y.margin=[3.2,0.3]
siz=.9
;
@parameters
default,lread,1
if lread then begin
  ;power,'_mag','_Str',k=k,spec1=spec1,spec2=stre1,i=n,tt=t,/noplot,/lks
  power,'_Str','_Str',k=k,spec1=stre1,spec2=stre1,i=n,tt=tstr,/noplot,/lks
  power,'_mag','_mag',k=k,spec1=spec1,spec2=spec1,i=n,tt=tmag,/noplot,/lks
  power,'_GWs','_GWh',k=k,spec1=grav1,spec2=grah1,i=n,tt=t,/noplot,/lks
  pc_read_ts,obj=ts
  lread=0
endif
;
;  select times from t array
;
ntt=n_elements(tt)
it_ts=intarr(ntt)
it_str=intarr(ntt)
it_mag=intarr(ntt)
it_gws=intarr(ntt)
for itt=0,ntt-1 do begin
  it_str[itt]=findex(tt[itt],tstr)
  it_mag[itt]=findex(tt[itt],tmag)
  it_gws[itt]=findex(tt[itt],t)
  it_ts[itt]=findex(tt[itt],ts.t)
endfor
print,it_str
print,it_mag
print,it_gws
;
xlab1=12.
xlab1b=600.
xlab2a=1.13
xlab2b=3.5
xrt=[1.,4.1]
xr=[k[1],max(k)]
!p.multi=[0,2,4]
col=[col0,col0,col0,col0,55,122,155]
dot='!s!9!A .!n!6!r'
circ_sym,.9,1
ps=8
;
;---------------------------------------------------------------------------------
!x.title='!6'
!y.title='!8E!6!dM!n(!8k!6)'
yr_mag=[1d-12,1d-3]
plot_oo,xr,yr_mag,/nodata
for itt=0,ntt-1 do oplot,k[1:*],spec1[1:*,it_mag[itt]],col=col[itt]
xyouts,xlab1,1e-4,siz=siz,'!6(a)'
;
!y.title='!13E!6!dM!n'
plot,ts.t,.5*ts.brms^2,xr=xrt
for itt=0,ntt-1 do oplot,[1,1]*ts.t[it_ts[itt]],[1,1]*.5*ts.brms[it_ts[itt]]^2,col=col[itt],ps=ps
xyouts,xlab2a,.0128,siz=siz,'!6(b)'
xx=[.3,3.] & oplot,xx,.1/xx^.6667
;
;---------------------------------------------------------------------------------
!y.title='!6Sp(!8S!6)'
!y.title='!6Sp(!8T!6!dTT!n)'
yr_str=[1d-11,1d-4]
plot_oo,xr,yr_str,/nodata
for itt=0,ntt-1 do oplot,k[1:*],(tstr[it_str[itt]]/6.)^2*stre1[1:*,it_str[itt]],col=col[itt]
for itt=0,ntt-1 do print,'t=',tstr[it_str[itt]]
;for itt=0,ntt-1 do oplot,k[1:*],col=188,1e-10*k[1:*]^2*t[it_str[itt]]^(4./3.)
for itt=0,ntt-1 do print,t[it_str[itt]]
xyouts,xlab1,2e-5,siz=siz,'!6(c)'
circ_sym,.9,0
;oplot,[11.,30.,100.],[1.6e-10,3e-9,6e-8],ps=8
;xx=[10.,180.] & oplot,xx,6.0e-13*xx^2.5,li=3
oplot,[11.,30.,100.],[4e-10,8e-10,1.6e-9],ps=8
xx=[10.,180.] & oplot,xx,7.5e-11*xx^0.7,li=3
circ_sym,.9,1
;
!y.title='!6<!8T!6!s!u2!n!r!dTT!n>'
tot_stre=total(stre1,1)
plot,tstr,(tstr/6.)^2*tot_stre,xr=xrt
for itt=0,ntt-1 do oplot,[1,1]*tstr[it_str[itt]],[1,1]*(tstr[it_str[itt]]/6.)^2*tot_stre[it_str[itt]],col=col[itt],ps=ps
oplot,ts.t,.4*ts.brms^4,li=1
xyouts,xlab2a,.00032,siz=siz,'!6(d)'
;
;---------------------------------------------------------------------------------
!y.title='!6Sp('+dot+'!8h!6)'
;!y.title='!8E!6!dM!n(!8k!6)'
yr_GWs=[1d-21,5d-10]
plot_oo,xr,yr_GWs,/nodata
for itt=0,ntt-1 do oplot,k[1:*],grav1[1:*,it_GWs[itt]],col=col[itt]
oplot,[1,1]*5000.,yr_GWs,li=3
oplot,[1,1]*50.,yr_GWs,li=3
xyouts,siz=siz,56.,1e-20,'!8k!6!d0!n'
xyouts,siz=siz,3000.,5e-12,'!8k!7!dl!n!6'
xyouts,xlab1b,18e-12,siz=siz,'!6(e)'
xx=[14.,120.] & oplot,xx,1.1e-8/xx^1.5,li=3
xyouts,130.,5e-12,siz=siz,'!9A!8k!6!u-1.5!n'
;
!y.title='!13E!6!dGW!n'
plot,ts.t,ts.EEGW,xr=xrt
for itt=0,ntt-1 do oplot,[1,1]*ts.t[it_ts[itt]],[1,1]*ts.EEGW[it_ts[itt]],col=col[itt],ps=ps
xyouts,xlab2a,3e-11,siz=siz,'!6(f)'
;
;---------------------------------------------------------------------------------
!x.title='!8k!6'
!y.title='!6Sp(!8h!6)'
yr_GWh=[1d-25,4d-12]
plot_oo,xr,yr_GWh,/nodata
for itt=0,ntt-1 do oplot,k[1:*],grah1[1:*,it_GWs[itt]],col=col[itt]
oplot,[1,1]*5000.,yr_GWh,li=3
oplot,[1,1]*50.,yr_GWh,li=3
xyouts,siz=siz,56.,1e-24,'!8k!6!d0!n'
xyouts,siz=siz,3000.,1e-14,'!8k!7!dl!n!6'
arrow,28.,2e-23,40.,1e-21,/data
xyouts,siz=siz,14.,1e-24,'!6break'
xyouts,xlab1b,5e-14,siz=siz,'!6(g)'
xx=[14.,220.] & oplot,xx,4e-10/xx^2.0
xyouts,100.,5e-14,siz=siz,'!9A!8k!6!u-2!n'
;
!x.title='!8t!6'
!y.title='!8h!6!s!u2!n!r!drms!n'
plot,ts.t,ts.hrms^2,xr=xrt
for itt=0,ntt-1 do oplot,[1,1]*ts.t[it_ts[itt]],[1,1]*ts.hrms[it_ts[itt]]^2,col=col[itt],ps=ps
xyouts,xlab2a,1.16e-12,siz=siz,'!6(h)'
;
print,'$mv idl.ps ~/GitHub/Axel/GWspectra_from_CME/fig/pstress_etc.eps'
END