;$Id: pspecm_wkf6_600_varyT.pro,v 1.1 2023/04/15 01:26:22 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=5,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
end
;
dir1='M512sig1_k6_ramp1a'
dir2='M512sig1_k6_ramp1a_f001'
dir4='M512sig1_k6_ramp1c_f002'
dir5='M512sig1_k6_ramp1b_f007'
;
siz=1.9 & !p.charsize=1.8
!x.margin=[7.5,0.3]
!y.margin=[3.2,0.3]
!p.multi=0
fac_h=1.263e-18
fac2=2.
half='!s!u 1!n!r!s-!r!d 2!n'
specm='/specm.sav'
label=''
;
;  determine normalization
;
default,maglabel,'mag'
default,green_lines,0
default,plot_polarization,1
ytickf='logticks_exp'
yr=yr_pspecm_scl3
yr=[1e-12,1e-5]
;xr1=[1e-9,3e-7]
;xr2=[1e-4,6e-2]
xr1=[1.4e-9,3e-7]
xr2=[2.2e-4,6e-2]
;
;-----------------------------------------------------------------------------
;  kf=6 runs
;
dirA='M512sig1_kf6_ramp1b_f1'
dirB='M512sig1_kf6_ramp1a_f2'
dirC='M512sig1_kf6_ramp1c_f3'
dirD='M512sig1_kf6_ramp1d_f5'
;
;-----------------------------------------------------------------------------
;  QCD scaling
;
T_in100GeV=1. & gstar_in100=1.
T_in100GeV=1.5d-3 & gstar_in100=.15
T_in100GeV=8.0d-3 & gstar_in100=.15
gS_in100=1.
;gS_in100=.0391
tend=1.
Hstar=2.066d10*T_in100GeV^2*gstar_in100^.5
H0=3.241d-18
a0=1.254d15*T_in100GeV*gS_in100^(1./3.)
;a0=1d12
fact=(Hstar/H0)^2*(tend/a0)^4
;
;Hstar=2.066d10
;a0=1.2545d15
freq=Hstar/(2.*!pi*a0)
H0=3.24d-18
;changed from H0 = 100 km/s/Mpc to 70 km/s/Mpc
;H0=2.268545d-18
rhocrit_fac=3./(8.*!pi)
;fac_h=1.263e-18
;changed H0 from 100 to 70 km/s/Mpc, and computed h_c
;from formula 
fac_h=sqrt(3./2.)*H0/!pi
fac=(Hstar/H0)^2/rhocrit_fac*(1./a0)^4
facGW=1./(16.*!pi)
circ_sym,0.7,1
mixed=0.
;
print,'fact=',fact
print,'fac=',fac
;
restore,'../'+dirA+specm
grav_tot=grav1m+grah1m/tm^2+mixed
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
plot_oo,freq*k,OmegaGW,xr=xr1,yr=yr,ytickf=ytickf,/nodata
;
;  NANOGrav
;
loadct,0
emma='../../../emma/GW/idl/'
up=rtable(emma+'data/GreenUpper2021.csv',2)
dn=rtable(emma+'data/GreenLower2021.csv',2)
xxx=10.^[reform(dn[0,*]),reverse(reform(up[0,*])),dn[0,0]]
yyy=10.^[reform(dn[1,*]),reverse(reform(up[1,*])),dn[1,0]]
polyfill,xxx,yyy<(.9*yr[1]),col=222
loadct,5
;
xx1=10.^reform(dn[0,*])
xx2=10.^reform(up[0,*])
yy1=10.^reform(dn[1,*])
yy2=10.^reform(up[1,*])
;
th=5
fmax=1e-7
oplot,freq*k,OmegaGW
good=where(freq*k le fmax)
pbla=linfit(freq*k(good),OmegaGW(good))
oplot,xx1,xx1*pbla[1]+pbla[0],th=th
;
restore,'../'+dirB+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55
good=where(freq*k le fmax and freq*k ne 0.)
pblu=linfit(alog(freq*k(good)),alog(OmegaGW(good)))
oplot,xx1,exp(alog(xx1)*(pblu[1])[0]+(pblu[0])[0]),th=th,col=55
;
restore,'../'+dirC+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=155
good=where(freq*k le fmax and freq*k ne 0.)
pora=linfit(alog(freq*k(good)),alog(OmegaGW(good)))
oplot,xx1,exp(alog(xx1)*(pora[1])[0]+(pora[0])[0]),th=th,col=155
;
restore,'../'+dirD+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=122
good=where(freq*k le fmax and freq*k ne 0.)
pred=linfit(alog(freq*k(good)),alog(OmegaGW(good)))
oplot,xx1,exp(alog(xx1)*(pred[1])[0]+(pred[0])[0]),th=th,col=122
;
xx=[1.3e-9,1.0e-8] & oplot,li=3,xx,5.0e+6*xx^1.6
xx=[1.8e-9,1.3e-8] & oplot,li=3,xx,2.0e+3*xx^1.6
xyouts,2.0e-9,2.0e-7 ,'!9A!8f!6!u1.6!n',siz=siz
xyouts,6.0e-9,3.0e-11,'!9A!8f!6!u1.6!n',siz=siz
;
print,'$mv idl.ps ~/GitHub/Tina/GWs-BBN/fig/pspecm_wkf6_600.eps'
;print,'$convert idl.ps /D/Print/pspecm_scl5_comp2_single.png'
END