;$Id: pspecm_scl5_comp.pro,v 1.26 2020/08/06 17:06:00 roper Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=16,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
dir0='E288e_t11_km1b'    ;(acB)
dir1='E288e_t11_M4e'     ;(ac1)
dir1='E1152f_t11_M3a_single'     ;(ac1, 2.67e-16, L=0.12, eta=5e-5, F=3, dt=8e-5)
dir1='E1152g2_t11_M4a_single'     ;(ac1, 2.9e-18, L=0.31, eta=5e-4, F=10, dt=7e-5)
dir1='E1152g2_t11_M5b_single'     ;(ac1, 8.7e-18, L=0.31, eta=1e-3, F=20, dt=3e-5)
dir1='E1152g2_t11_M5c_single'     ;(ac1,        , L=0.31, eta=1e-3, F=50, dt=7e-5)
dir1='E1152e_t11_M4a_single'     ;(ac1, 5.0e-13, L=0.031, eta=5e-5, F=3, dt=5e-6)
dir1='E1152e_t11_M4f_double_rep' ;(ac1, but with slightly more viscosity)
dir1='E1152e_t11_M4d_double'     ;(ac1, consistent with table)
dir2='F1152d2_sig1_t11_M2c' ;(hel1, consistent with table, but redo in double)
dir2='F1152d2_sig1_t11_M2c_double' ;(hel1)
;dir2='F288e_sig1_t11_M4' ;(hel1, old)
;dir2='F1152d2_sig1_t11_M4' ;(hel1)
;dir2='F1152e_sig1_t11_M4' ;(hel1)
;dir2='F1152b_sig1_t11_M4' ;(hel1)
;dir3='F288a_sig1_t11b'   ;(hel2)
;dir3='F1152a_sig1_t11b'   ;(hel2, ?)
dir3='F1152a_sig1_t11d_double'        ;(hel2)
dir4='F1152b_sig1'        ;(hel3 single)
dir5='F1152a_sig0_t11b'   ;(noh2)
dir5='F1152b_sig0_t11_M4'   ;(noh1, preferred)
;
yr_pspecm_scl3=[1d-16,5d-8]
yr_pspechm_scl3=[2d-24,6.2d-16]
;
yr_pspecm_scl3=[1d-20,9d-10]
yr_pspechm_scl3=[2d-27,6.2d-20]
;
yr_pspecm_scl3=[1d-28,2d-9]
yr_pspechm_scl3=[2d-30,6.2d-20]
yr_pspechm_scl3=[1d-31,6.2d-20]
yr_pspechm_scl3=[5d-29,3.2d-20]
;
yr_pspecm_scl3=[5d-25,2d-9]
yr_pspechm_scl3=[2d-28,3.2d-20]
;
thick4=5
siz=1.5
!p.charsize=1.7
!x.margin=[7.3,0.5]
!y.margin=[3.2,0.2]
!p.multi=[0,1,2]
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
;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)
xr=[2.4e-4,4.0e-2]
xr=[2.4e-5,8.0e-2]
xr=[2.4e-4,16.0e-2]
circ_sym,0.7,1
mixed=0.
;
;  1st panel
;
!p.title='!6';+run
!x.title='!6'
!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
ikmax1=500 ;(to suppress the flat part)
plot_oo,freq*k[1:ikmax1],OmegaGW[1:ikmax1],xr=xr,yr=yr,ytickf=ytickf
k1=k
k1[1:*]=1./k[1:*]
print,'AXEL_dir1: ',total(OmegaGW*k1)
;
restore,'../'+dir0+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k,OmegaGW,li=2,col=155
;
restore,'../'+dir2+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=122
print,'freq,k[1],OmegaGW[1]=',freq,k[1],OmegaGW[1]
k1=k
k1[1:*]=1./k[1:*]
print,'freq,k[1],100*total(k1*OmegaGW[1])=',freq,k[1],100*total(k1*OmegaGW[1])
;help,k
;
restore,'../'+dir3+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55
print,dir3,facGW*fac*minmax(k[1:*])
help,k
;
restore,'../'+dir4+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
ikmax=220 ;(to suppress the flat part)
oplot,freq*k[1:ikmax],OmegaGW[1:ikmax],col=55,li=1
;
restore,'../'+dir5+specm
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
oplot,freq*k,OmegaGW,col=55,li=2
;
;  first plot
;
xx=4e-4 & dx=4e-4
legend,xx,dx,1e-19,1,siz=siz,'!6hel3',col=55
legend,xx,dx,1e-21,0,siz=siz,'!6hel2',col=55
legend,xx,dx,1e-23,2,siz=siz,'!6noh1',col=55
;
xyouts,5.0e-2,1.0e-17,'!6ac1',siz=siz
xyouts,4.6e-2,1.7e-22,'!6hel1',siz=siz,col=122
;xyouts,2.8e-3,4.3e-16,'!6hel2',siz=siz,col=55
;
if green_lines then begin
  loadct,6
  xx=6*[1.6e-3,0.7e-2] & oplot,xx,6*4.0e-7/(xx/.001)^0.667,col=122
  xx=6*[2.0e-3,1.0e-2] & oplot,xx,2.0e-7/(xx/.001)^2.667,col=122
  xyouts,.05,1e-7,'!9A!8f!6!u-2/3!n',col=122,siz=siz
  xyouts,.04,5e-13,'!9A!8f!6!u-8/3!n',col=122,siz=siz
  ;
  xx=[5e-4,1.4e-3] & oplot,xx,2e-7*(xx/.001)^5.,col=122
  xx=[5e-4,2.4e-3] & oplot,xx,1e-9*(xx/.001)^2.2,col=122
  xyouts,4e-4,5e-8,'!9A!8f!6!u5!n',col=122,siz=siz
  xyouts,2.7e-4,4e-11,'!9A!8f!6!u2.2!n',col=122,siz=siz
  loadct,5
endif
;
;  Caprini+16, Cornish, Maggiore00 points
;
loadct,6
dir='~/tex/roper/PencilGW/data/'
file=dir+'Cornish.rtf'
c1=rtable(file,2,head=1)
;oplot,10^c1[0,*],((10^c1[1,*]/fac_h)*10^c1[0,*])^2,li=4,col=122
;
file=dir+'Cornish20Yrs.dat'
c2=rtable(file,2,head=1)
;oplot,10^c2[0,*],((10^c2[1,*]/fac_h)*10^c2[0,*])^2,li=3,col=122
;
file=dir+'Maggiore.rtf'
m1=rtable(file,2,head=1)
;
file=dir+'Caprini-C3.rtf'
a3=rtable(file,2,head=1)
;oplot,10^a3[0,*],10^a3[1,*],li=2,col=122
;
file=dir+'Caprini-C1.rtf'
a1=rtable(file,2,head=1)
;oplot,10^a1[0,*],10^a1[1,*],li=1,col=122
;
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,3.0e-2,1.0e-13,'!6(i)',siz=siz,col=122
;xyouts,1.0e-2,1.5e-11,'!6(ii)',siz=siz,col=122
;xyouts,9e-3,3.4e-10,'!6(iii)',siz=siz,col=122
;xyouts,4e-3,3.8e-10,'!6(iv)',siz=siz,col=122
;
loadct,5
;
;  overplot
;
xx=[3e-4,8e-4] & oplot,li=3,xx,1e-11*xx
xx=[7e-3,3e-2] & oplot,li=3,xx,5e-23/xx^(5.)
xyouts,6.0e-4,3.0e-16,'!9A!8f!6',siz=siz*.8
xyouts,1.5e-2,2.8e-13,'!9A!8f!6!u-5!n',siz=siz*.8
;
;  2nd panel
;
!y.title='!8h!6!dc!n(!8f!6)'
!x.title='!8f!6 [Hz]'
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^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;plot_oo,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),xr=xr,yr=yr,ytickf=ytickf
plot_oo,freq*k[1:ikmax1],(sqrt(k*grah1m)/a0)[1:ikmax1],xr=xr,yr=yr,ytickf=ytickf
;
;  0th line
;
restore,'../'+dir0+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),col=155
;oplot,freq*k,8e-16*sqrt(k*grah1m),co=155,li=1
;
;  2nd line
;
restore,'../'+dir2+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),col=122
oplot,freq*k,sqrt(k*grah1m)/a0,col=122
;
;  3rd line
;
restore,'../'+dir3+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),col=55
oplot,freq*k,sqrt(k*grah1m)/a0,col=55
;oplot,freq*k,8e-16*sqrt(k*grah1m),co=55,li=1
;
;  4th line
;
restore,'../'+dir4+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),col=55,li=1
oplot,freq*k,sqrt(k*grah1m)/a0,col=55,li=1
;oplot,freq*k,8e-16*sqrt(k*grah1m),co=55,li=1
;
;  5th line
;
restore,'../'+dir5+specm
fach=1.263e-18/(freq*k[1:*])
grav_tot=grav1m+grah1m/tm^2+mixed
OmegaGW=facGW*fac*k*grav_tot
;oplot,freq*k[1:*],fach*sqrt(OmegaGW[1:*]),col=55,li=2
oplot,freq*k,sqrt(k*grah1m)/a0,col=55,li=2
;oplot,freq*k,8e-16*sqrt(k*grah1m),co=55,li=1
;
;  second plot
;
xx=4e-4 & dx=4e-4 & yy=2e-25
legend,xx,dx,yy,1,siz=siz,'!6hel3',col=55
legend,xx,dx,yy/1e1,0,siz=siz,'!6hel2',col=55
legend,xx,dx,yy/1e2,2,siz=siz,'!6noh1',col=55
;
xyouts,3.6e-3,3.0e-21,'!6ac1',siz=siz
xyouts,7.5e-2,3.0e-26,'!6hel1',siz=siz,col=122
;xyouts,2.8e-3,4.3e-24,'!6hel2',siz=siz,col=55
;
;  continous limits
;
loadct,6
;oplot,10^c1[0,*],10^c1[1,*],li=4,col=122
;oplot,10^c2[0,*],10^c2[1,*],li=3,col=122
;oplot,10^a3[0,*],fac_h*sqrt(10^a3[1,*])/10^a3[0,*],li=2,col=122
;oplot,10^a1[0,*],fac_h*sqrt(10^a1[1,*])/10^a1[0,*],li=1,col=122
oplot,10^n1[0,*],10^n1[1,*],li=3,col=122,thick=thick4
;
;xyouts,3.6e-2,2.0e-23,'!6(i)',siz=siz,col=122
;xyouts,1.2e-2,6.3e-22,'!6(ii)',siz=siz,col=122
;xyouts,9.0e-3,1.9e-21,'!6(iii)',siz=siz,col=122
;xyouts,5.0e-3,5.3e-21,'!6(iv)',siz=siz,col=122
;
loadct,5
;
;  overplot lines
;
xx=[3.0e-4,1e-3] & oplot,li=3,xx,4e-24/xx^.5
xx=[7.0e-3,3e-2] & oplot,li=3,xx,1e-29/xx^(7./2.)
xyouts,5.0e-4,2.0e-23,'!9A!8f!6!u-1/2!n',siz=siz*.8
xyouts,1.5e-2,4.2e-23,'!9A!8f!6!u-7/2!n',siz=siz*.8
;
if green_lines then begin
  loadct,6
  xx=6*[2.0e-3,1.8e-2] & oplot,xx,5.0e-19/(xx/.001)^(7./3.),col=122
  xyouts,.02,2e-23,'!9A!8f!6!u-7/3!n',col=122,siz=siz
  ;
  xx=[5e-4,2.4e-3] & oplot,xx,7e-21*(xx/.001)^(.1),col=122
  xyouts,2.7e-4,1e-21,'!9A!8f!6!u0.1!n',col=122,siz=siz
  loadct,5
endif
;
print,'$mv idl.ps ~/tex/roper/GW/fig/pspecm_scl5_comp.ps'
END