;$Id: pspecm_CME_SerB.pro,v 1.4 2022/11/14 15:43:31 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
;
dirA='512_1e2_1e4_4e8_1em4a'
dirB='512_1e2_1e4_4e8_2em4a'
dirC='512_1e2_1e4_4e8_5em4a'
dirD='512_1e2_1e4_4e8_1em3a'
;
yr_pspecm_scl3=[1d-18,9d-12]
;
siz=.7 & !p.charsize=1.
!x.margin=[7.5,0.3]
!y.margin=[3.2,0.3]
!p.multi=[0,2,1]
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
xr1=[1.0e-7,5e-5]
xr2=[2.2e-4,6e-2]
;
;-----------------------------------------------------------------------------
;  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
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
;
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
print,'X ',freq*k
print,'Y ',OmegaGW
;
xx=[1.3e-7,1.0e-5] & oplot,li=3,xx,6.0e-10*xx^0.5
xyouts,2.0e-6,2.0e-12 ,'!9A!8f!6!u1/2!n',siz=siz
;
;-----------------------------------------------------------------------------
yr=[1d-1,1.01d-0]
;
restore,'../'+dirA+specm
grav_tot=grav1m+grah1m/tm^2+mixed
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
polarGW=grah2m/grah1m
plot_oo,freq*k,polarGW,xr=xr1,yr=yr,ytickf=ytickf,/nodata
;
restore,'../'+dirB+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
polarGW=grah2m/grah1m
oplot,freq*k,polarGW,col=55
;
restore,'../'+dirC+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
polarGW=grah2m/grah1m
oplot,freq*k,polarGW,col=155
;
restore,'../'+dirD+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
polarGW=grah2m/grah1m
oplot,freq*k,polarGW,col=122
;
!p.multi=0
print,'$mv idl.ps ~/GitHub/Tina/NSF-2022-AAG-GWs/pspecm_CME_SerB.eps'
END