;$Id: pspecm_wkf6.pro,v 1.2 2021/09/21 02:37:28 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=12,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
dir1='M512sig1_k6_ramp1a'
dir2='M512sig1_k6_ramp1a_f001'
dir4='M512sig1_k6_ramp1c_f002'
dir5='M512sig1_k6_ramp1b_f007'
;
yr_pspecm_scl3=[1d-16,2d-5]
;
siz=.8 & !p.charsize=1.
siz=1.3 & !p.charsize=1.7
!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
;
Hstar=2.066d10
a0=1.2545d15
freq=Hstar/(2.*!pi*a0)
H0=3.24d-18
rhocrit_fac=3./(8.*!pi)
fac=(Hstar/H0)^2/rhocrit_fac*(1./a0)^4
facGW=1./(16.*!pi)
xr=[1e-9,6e-2]
circ_sym,0.7,1
mixed=0.
;
!x.title='!6frequency !8f!6 [Hz]'
!y.title='!6!8h!6!s!d0!n!r!u2!n!7X!6!dGW!n(!8f!6)'
;
restore,'../'+dir1+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
plot_oo,freq*k,OmegaGW,xr=xr,yr=yr,ytickf=ytickf,/nodata
oplot,freq*k,OmegaGW
;
restore,'../'+dir2+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55
;
restore,'../'+dir4+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=155
;
restore,'../'+dir5+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=122
;
xx=[2.8e-4,1.3e-3] & oplot,li=3,xx,5.0e-10*xx
xx=[2e-4,1e-3] & oplot,li=3,xx,2.0e+0*xx^2.5
xyouts,2.0e-4,1.9e-8,'!9A!8f!6!u2.5!n',siz=siz
xyouts,5.0e-4,1.0e-13,'!9A!8f!6',siz=siz
;
;  Cornish points
;
loadct,6
dir='~/tex/roper/PencilGW/data/'
file=dir+'Robson_Cornish.rtf'
n1=rtable(file,2,head=1)
oplot,10^n1[0,*],((10^n1[1,*]/fac_h)*10^n1[0,*])^2,li=3,col=122,thick=thick4
xyouts,0.9e-2,4.2e-11,'!6  LISA!climit',siz=siz*.8,col=122
loadct,5
;
;-----------------------------------------------------------------------------
;  kf=6 runs
;
dirA='M512sig1_kf6_ramp1b_f1'
dirB='M512sig1_kf6_ramp1a_f2'
dirC='M512sig1_kf6_ramp1c_f3'
dirD='M512sig1_kf6_ramp1d_f5'
;
restore,'../'+dirA+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,li=1
;
restore,'../'+dirB+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,li=1,col=55
;
restore,'../'+dirC+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,li=1,col=155
;
restore,'../'+dirD+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,li=1,col=122
;
print,'fac (first)=',fac
;-----------------------------------------------------------------------------
;  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
;
;  QCD scaling
;
T_in100GeV=1. & gstar_in100=1.
T_in100GeV=1.5d-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
oplot,freq*k,OmegaGW
;
restore,'../'+dirB+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55
;
restore,'../'+dirC+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=155
;
restore,'../'+dirD+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=122
;
;xx=[3e-9,2e-3] & oplot,li=3,xx,2.0e+5*xx^2
xx=[1.7e-9,1.2e-8] & oplot,li=3,xx,6.0e+7*xx^1.7
xyouts,2.0e-9,5.0e-7,'!9A!8f!6!u1.7!n',siz=siz
;
print,'$mv idl.ps ~/tex/tina/BBN/fig/pspecm_wkf6.eps'
;print,'$convert idl.ps /D/Print/pspecm_scl5_comp2_single.png'
END