;$Id: pstress_etc_LowFreq.pro,v 1.11 2022/10/09 03:09:27 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=10,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
cwd,run
;
;  this is to be run in:
;  LowFreq1024sig1_k01_kf10b
;  LowFreq1024sig0_k01_kf10b2
;
!p.charsize=.9
!x.margin=[8.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,'_SCL','_SCL',k=k,spec1=stre2,spec2=stre2,i=n,tt=tst2,/noplot,/lks
  power,'_VCT','_VCT',k=k,spec1=stre3,spec2=stre3,i=n,tt=tst3,/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_st2=intarr(ntt)
it_st3=intarr(ntt)
it_mag=intarr(ntt)
it_gws=intarr(ntt)
print
;print,'tstr=',tstr
;print
for itt=0,ntt-1 do begin
  it_str[itt]=findex(tt[itt],tstr)
  it_st2[itt]=findex(tt[itt],tst2)
  it_st3[itt]=findex(tt[itt],tst3)
  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=.13
xlab1b=25.
xlab2a=1.13
xlab2b=3.5
xrt=[1.,4.1]
xrt=[1.,80.]
xr=[k[1],max(k)]
!p.multi=[0,2,2]
if ntt eq 7 then begin
  col=[col0,col0,col0,col0,55,122,155]
endif else begin
  col=[col0,col0,col0,55,122,155]
endelse
li=indgen(ntt)
;
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-4]
yr_mag=[2d-9,1d-3] ;(new run)
yr_mag=[2d-12,3d-5] ;(new run, sig0-zero)
yr_mag=[2d-12,3d-4] ;(new run, sig1-zero)
plot_oo,xr,yr_mag,/nodata
for itt=0,ntt-1 do oplot,k[1:*],spec1[1:*,it_mag[itt]],col=col[itt],li=li[itt]
xyouts,xlab1,9e-6,siz=siz,'!6(a)'
;xx=[.4,1.1] & oplot,xx,2e-5*xx^5
xx=[.3,3.1] & oplot,xx,1e-9*xx^4,li=3
xyouts,1.,2.3e-10,'!9A!8k!6!u4!n',siz=siz
;
siz=.9
;t(GW)=      1.54821
;t(GW)=      2.07686
;t(GW)=      5.04138
;t(GW)=      10.0449
;t(GW)=      50.0659
;t(GW)=      140.086
xx=4. & dx=3. & yy=3e-8 & dy=5.
legend,xx,dx,yy,li[0],'!8t!6 = 1.6',siz=siz
legend,xx,dx,yy/dy,li[1],'!8t!6 = 2.1',siz=siz
legend,xx,dx,yy/dy^2,li[2],'!8t!6 = 5.0',siz=siz
legend,xx,dx,yy/dy^3,li[3],'!8t!6 = 10',col=col[3],siz=siz
legend,xx,dx,yy/dy^4,li[4],'!8t!6 = 50',col=col[4],siz=siz
legend,xx,dx,yy/dy^5,li[5],'!8t!6 = 140',col=col[5],siz=siz
;---------------------------------------------------------------------------------
; !y.title='!6Sp(!8T!6!s!u2!n!r!dS!n)'
; plot_oo,xr,yr_str,/nodata
; for itt=0,ntt-1 do oplot,k[1:*],(tst2[it_st2[itt]]/6.)^2*stre2[1:*,it_st2[itt]]*k[1:*]^4,col=col[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,'tst2=',tst2[it_st2[itt]]
; xyouts,xlab1,2e-5,siz=siz,'!6(c)'
; print
;
; !y.title='!6Sp(!8T!6!s!u2!n!r!dV!n)'
; plot_oo,xr,yr_str,/nodata
; for itt=0,ntt-1 do oplot,k[1:*],(tst3[it_st3[itt]]/6.)^2*stre3[1:*,it_st3[itt]]*k[1:*]^2,col=col[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,'tst3=',tst3[it_st3[itt]]
; xyouts,xlab1,2e-5,siz=siz,'!6(c)'
; print
;
!y.title='!6Sp(!8T!6!dTT!n)'
yr_str=[5d-10,2d-6]  ;(new run)
yr_str=[2d-13,5d-8]  ;(new run)
yr_str=[2d-13,5d-7]  ;(new run, sig1-zero)
yr_str=[1d-15,1d-5]  ;(paper)
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],li=li[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,'tstr=',tstr[it_str[itt]]
xyouts,xlab1,3e-7,siz=siz,'!6(b)'
print
;
; !y.title='!6<!8T!6!s!u2!n!r!dTT!n>'
; tot_stre=total(stre1,1)
; ;plot,tmag,tot_stre,xr=xrt
; plot_oo,tmag,tot_stre,xr=xrt
; for itt=0,ntt-1 do oplot,[1,1]*tstr[it_str[itt]],[1,1]*tot_stre[it_str[itt]],col=col[itt],ps=ps
; ;oplot,ts.t,4.*ts.brms^4,li=1
; xyouts,xlab2a,.0035,siz=siz,'!6(d)'
; xx=[10.,70.] & oplot,xx,1e-6/xx^(13./3.)
; xx=[10.,70.] & oplot,xx,3e-7/xx^(10./3.)
;
;---------------------------------------------------------------------------------
!x.title='!8k!6'
!y.title='!6Sp('+dot+'!8h!6)'
;!y.title='!8E!6!dM!n(!8k!6)'
yr_GWs=[1d-11,1d-7]  ;(new run)
yr_GWs=[1d-11,1d-8]  ;(new run)
yr_GWs=[1d-11,1d-7]  ;(new run, sig1-zero)
yr_GWs=[1d-10,3d-8]  ;(paper)
plot_oo,xr,yr_GWs,/nodata
for itt=0,ntt-1 do oplot,k[1:*],grav1[1:*,it_GWs[itt]],col=col[itt],li=li[itt]
for itt=0,ntt-1 do print,'t(GW)=',t[it_GWs[itt]]
oplot,[1,1]*5000.,yr_GWs,li=3
oplot,[1,1]*50.,yr_GWs,li=3
;for itt=0,ntt-1 do oplot,2.*!pi*[1,1]/ts.t[it_ts[itt]],yr_GWs,li=3,col=col[itt]
;for itt=0,ntt-1 do oplot,!pi*[1,1]/(reform(ts.t[it_ts[itt]]))[0],yr_GWs,li=3,col=col[itt]
;-xyouts,siz=siz,56.,1e-20,'!8k!6!d0!n'
;-xyouts,siz=siz,3000.,5e-12,'!8k!7!dl!n!6'
xyouts,xlab1b,1.4e-8,siz=siz,'!6(c)'
if run eq 'LowFreq1024sig1_k01_kf10b_rep2' then begin
  xx=[.12,.6] & oplot,xx,1e-8*xx^1.0,li=3
  xyouts,.2,4e-9,'!9A!8k!6',siz=siz
endif else begin
  xx=[.12,.6] & oplot,xx,3e-8*xx^1.5,li=3
  xyouts,.2,6e-9,'!9A!8k!6!u1.5!n',siz=siz
endelse
;
; !y.title='!13E!6!dGW!n'
; plot,ts.t,ts.EEGW,xr=xrt,yr=[0.,1e-6]
; 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-14,1d-7]
plot_oo,xr,yr_GWh,/nodata
for itt=0,ntt-1 do oplot,k[1:*],grah1[1:*,it_GWs[itt]],col=col[itt],li=li[itt]
oplot,[1,1]*5000.,yr_GWh,li=3
oplot,[1,1]*50.,yr_GWh,li=3
;for itt=0,ntt-1 do oplot,!pi*[1,1]/(reform(ts.t[it_ts[itt]]))[0],yr_GWh,li=3,col=col[itt]
;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,1.7e-8,siz=siz,'!6(d)'
;
; !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)'
;
cwd,run
print,'$mv idl.ps ~/GitHub/Axel/GWspectra_from_CME/fig/pstress_etc_LowFreq_sig1_k05.eps'
print,'$mv idl.ps ~/GitHub/Axel/GWspectra_from_CME/fig/pstress_etc_LowFreq_sig1.eps'
print,'$mv idl.ps ~/GitHub/Axel/GWspectra_from_CME/fig/pstress_etc_LowFreq_sig0.eps'
print,'$mv idl.ps ~/GitHub/Axel/GWspectra_from_CME/fig/pstress_etc_LowFreq_'+run+'.eps'
print,'$mv idl.ps ~/GitHub/Ramkishor/GWspectra_from_mhd_turbulence/fig/pstress_etc_LowFreq_sig1_k05.eps'
print,'$mv idl.ps ~/GitHub/Ramkishor/GWspectra_from_mhd_turbulence/fig/pstress_etc_LowFreq_sig1.eps'
print,'$mv idl.ps ~/tex/prop/computing/pdc/22/fig/pstress_etc_LowFreq_sig1.eps'
print,'$mv idl.ps ~/GitHub/Ramkishor/GWspectra_from_mhd_turbulence/fig/pstress_etc_LowFreq_sig0.eps'
print,'$mv idl.ps ~/GitHub/Ramkishor/GWspectra_from_mhd_turbulence/fig/pstress_etc_LowFreq_'+run+'.eps'
END