;$Id: pspecm_hel_delk_kf600.pro,v 1.2 2021/03/25 10:22:40 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=6,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
end
;
dir1='M1152e_exp6k4'
dir2='M1152e_exp6k4_delk10'
dir3='M1152e_exp6k4_delk50'
dir4='M1152e_exp6k4_delk200'
;
xs=.34
ys1=.06
ys2=.6
;
xr=[2e-7,3e-4]
yr_pspecm_scl3=[1d-18,3d-13]
yr_pspechm_scl3=[1d-23,2.2d-18]
;
xr=[2e-4,2e-1]
yr_pspecm_scl3=[2d-18,3d-13]
yr_pspechm_scl3=[1d-26,2.2d-21]
;
thick4=5
siz=1.0
si2=1.0
!p.charsize=1.2
!x.margin=[7.8,0.5]
!y.margin=[3.2,0.2]
!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'
wspecm='/specm.txt'
label=''
;
;  determine normalization
;
default,maglabel,'mag'
default,green_lines,0
default,plot_polarization,1
ytickf='logticks_exp'
yr=yr_pspecm_scl3
;
;  new units
;
T_in100GeV=1.5d-3 & gstar_in100=.15
T_in100GeV=1. & gstar_in100=1.
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
print,'fact=',fact
;
;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.
;
;  1st panel
;
!p.title='!6';+run
!x.title='!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[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
plot_oo,xr,yr,ytickf=ytickf,/nodata
oplot,freq*k,OmegaGW
;
restore,'../'+dir2+specm
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55
;
loadct,6
restore,'../'+dir3+specm
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=122
loadct,5
;
restore,'../'+dir4+specm
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=155
;
;restore,'../'+dir5+specm
;grav_tot=grav1m+grah1m/tm[0]^2+mixed
;OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k,OmegaGW,col=122
;
;xyouts,xs*5.5e-8,ys1*6.0e-14,'!8k!6!df!n=2',siz=siz,col=155
;xyouts,xs*2.0e-7,ys1*1.0e-13,'!66',siz=siz
;xyouts,xs*7.0e-7,ys1*7.3e-13,'!620',siz=siz,col=55
;xyouts,xs*2.1e-6,ys1*4.8e-12,'!660',siz=siz,col=122
;
;xx=xs*[3e-9,1.0e-7] & oplot,li=3,xx,ys1*3e-2*xx
;xyouts,xs*1.4e-8,ys1*2.3e-11,'!9A!8f!6',siz=siz
xout=9.0e-2
xyouts,xout,8.0e-14,'!6(a)',siz=si2
;
;  2nd panel
;
!y.title='!8h!6!dc!n(!8f!6)'
default,kin,0
default,hdone,1
ytickf='logticks_exp'
yr=yr_pspechm_scl3
circ_sym,0.7,1
;
restore,'../'+dir1+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
plot_oo,freq*k,sqrt(k*grah1m)/a0,xr=xr,yr=yr,ytickf=ytickf,/nodata
oplot,freq*k,sqrt(k*grah1m)/a0
;
;  2nd line
;
restore,'../'+dir2+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,sqrt(k*grah1m)/a0,col=55
;
;  3rd line
;
loadct,6
restore,'../'+dir3+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,sqrt(k*grah1m)/a0,col=122
loadct,5
;
;  4th line
;
restore,'../'+dir4+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm[0]^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,sqrt(k*grah1m)/a0,col=155
;
;  5th line
;
;restore,'../'+dir5+specm
;fach=1.263e-18/(freq*k[1:*])
;grav_tot=grav1m+grah1m/tm[0]^2+mixed
;OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k,sqrt(k*grah1m)/a0,col=122
;
;xyouts,xs*3.0e-8,ys2*2.0e-17,'!8k!6!df!n=2',siz=siz,col=155
;xyouts,xs*1.5e-7,ys2*2.3e-17,'!66',siz=siz
;xyouts,xs*7.0e-7,ys2*1.8e-18,'!620',siz=siz,col=55
;xyouts,xs*1.4e-6,ys2*1.6e-17,'!660',siz=siz,col=122
;
;xx=xs*[2.6e-9,1.3e-7] & oplot,li=3,xx,ys2*6e-18/xx^.5
;xyouts,xs*1.9e-7,ys2*1.0e-14,'!9A!8f!6!u-1/2!n',siz=siz
xyouts,xout,6.0e-22,'!6(b)',siz=si2
;
print,"$mv idl.ps ~/GitHub/Yutong/MGW-NANOGrav/Figures/pspecm_hel_delk_kf600.eps"
END