;$Id: pcomp_last.pro,v 1.7 2021/06/13 08:53:46 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=22,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
end
;
siz=1.0
si2=1.0
!p.charsize=1.9
!x.margin=[8.2,0.5]
;!x.margin=[8.8,6.5]
!y.margin=[3.2,0.2]
;!y.margin=[3.2,3.2]
!p.multi=[0,2,3]
print,"$mv idl.ps ~/GitHub/Yutong/MGW-NANOGrav/Figures/pcomp_last.eps"
;
bar='!20!s!A$!n!r!6'
bar2='!20!s!u$!n!r!6'
xtit0='!6'
xtit1=bar2+'!7x!6'
xtit2='     '+bar2+'!8x!6'
ytit1='!6Sp(!8h!6)'
ytit2='!7D'+bar+'!8t!6'
xout=1.15
thick2=6
;
dir1='M1024e_exp6k4_k1_kf100_delk0_cont' & yr_k1=[2e-15,2e-10] & lev1=2e-5 & om1cut=0.
dir2='M1024e_exp6k4_k1_kf100_delk3' & yr_k2=[1e-15,1e-11] & lev2=8e-6 & om2cut=3.
dir3='M1024e_exp6k4_k1_kf100_delk10' & yr_k3=[1e-15,1e-12] & lev3=5e-6 & om3cut=10.
;
xr_k=[.9,330.]
grah_file='grah_specs.sav'
spec_file='specm.sav'
om_file='last.sav'
xr_fff=[-1.,1.]*!pi
yr_fff=[-1.,1.]*!pi
;
;-----------------------------------------------------------------------------
;  new units
;
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
print,'fact=',fact
;
;Hstar=2.066d10
;a0=1.2545d15
freq=Hstar/(2.*!pi*a0)
H0=3.24d-18
;
yr=3.15576e7    ;(s)
pc=3.08568e18   ;(cm)
clight=3e10
tfac=a0/Hstar/yr
xfac=clight*a0/Hstar/pc
;
xbox=[0.,1.,1.,0.,0.]
ybox=[0.,0.,1.,1.,0.]
;-----------------------------------------------------------------------------
;
lev=grange(-1.,1.,20)*lev1
restore,'../'+dir1+'/'+grah_file
restore,'../'+dir1+'/'+spec_file
restore,'../'+dir1+'/'+om_file
plot_oo,xr_k,yr_k1,xtit=xtit0,ytit=ytit1,xst=9,/nodata
oplot,k,grah1m,li=1,col=55
axis,/xax,xr=xr_k*freq,xtit='!8f!6 [Hz]'
oplot,om_first,ssp
xyouts,xout,4e-15,bar2+'!7x!6!dcut!n=0',siz=siz
ttlast0=max(tt_last)-!pi
contour,clip(fff_last,minmax(lev)),x,tt_last-ttlast0,/fill,lev=lev,xr=xr,yr=yr_fff,xtit=xtit0,ytit=ytit2,xst=9,yst=9
axis,/xax,xr=xr_fff*xfac,xtit='!8x!6 [pc]'
axis,/yax,yr=yr_fff*tfac,ytit='!7D!8t!6 [yr]'
!x.title='!8h !9X!610!u18!n'
!y.title='!6'
xx0=.96 & dxx=.01 & dyy=.1
yy1=.76
colorbar,pos=[xx0,yy1,xx0+dxx,yy1+dyy], range=minmax(1e18/a0*lev),/top,/ver,$
form='(i3)',charsize=1.2,div=4,ytit='!6',col=255
;
lev=grange(-1.,1.,20)*lev2
restore,'../'+dir2+'/'+grah_file
restore,'../'+dir2+'/'+spec_file
restore,'../'+dir2+'/'+om_file
plot_oo,xr_k,yr_k2,xtit=xtit0,ytit=ytit1,/nodata
oplot,k,grah1m,li=1,col=55
omk=sqrt(om2cut^2+k^2)
oplot,om_first,ssp*2.
oplot,omk,omk/k*grah1m,li=2,col=122,thick=thick2
xyouts,xout,2e-15,bar2+'!7x!6!dcut!n=3',siz=siz
ttlast0=max(tt_last)-!pi
contour,clip(fff_last,minmax(lev)),x,tt_last-ttlast0,/fill,lev=lev,yr=yr_fff,xtit=xtit0,ytit=ytit2,yst=9
;axis,/yax,yr=yr_fff*tfac,ytit='!7D!8t!6!dphys!n [yr]'
axis,/yax,yr=yr_fff*tfac,ytit='!7D!8t!6 [yr]'
s=2*!pi/3. & oplot,s*xbox-!pi,s*ybox-!pi,col=255,thick=thick2
yy1=.43
colorbar,pos=[xx0,yy1,xx0+dxx,yy1+dyy], range=minmax(1e18/a0*lev),/top,/ver,$
form='(i2)',charsize=1.2,div=4,ytit='!6',col=255
;
lev=grange(-1.,1.,20)*lev3
restore,'../'+dir3+'/'+grah_file
restore,'../'+dir3+'/'+spec_file
restore,'../'+dir3+'/'+om_file
plot_oo,xr_k,yr_k3,xtit=xtit1,ytit=ytit1,/nodata
oplot,k,grah1m,li=1,col=55
omk=sqrt(om3cut^2+k^2)
oplot,om_first,ssp
oplot,omk,omk/k*grah1m,li=2,col=122,thick=thick2
xyouts,xout,2e-15,bar2+'!7x!6!dcut!n=10',siz=siz
ttlast0=max(tt_last)-!pi
contour,clip(fff_last,minmax(lev)),x,tt_last-ttlast0,/fill,lev=lev,yr=yr_fff,xtit=xtit2,ytit=ytit2,yst=9
;axis,/yax,yr=yr_fff*tfac,ytit='!7D!8t!6!dphys!n [yr]'
axis,/yax,yr=yr_fff*tfac,ytit='!7D!8t!6 [yr]'
s=2*!pi/10. & oplot,s*xbox-!pi,s*ybox-!pi,col=255,thick=thick2
yy1=.10
colorbar,pos=[xx0,yy1,xx0+dxx,yy1+dyy], range=minmax(1e18/a0*lev),/top,/ver,$
form='(i2)',charsize=1.2,div=4,ytit='!6',col=255
;
; %%BoundingBox: 54 85 564 708
; %%PageBoundingBox: 54 85 595 730
; GitHub/Yutong/MGW-NANOGrav
;
END