;$Id: pspecm_scl3n.pro,v 1.11 2019/07/03 19:30:32 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
;
;  mv idl.ps ~/tex/roper/GW/fig/pspecm_scl2.ps
;  mv idl.ps ~/tex/prop/computing/pdc/18/fig/pspecm_scl3n.ps
;
siz=1.3
@parameters
!p.charsize=1.7
;!p.charsize=1.8
!x.margin=[7.6,0.5]
!y.margin=[3.2,0.2]
!p.multi=[0,1,2]
pc_read_param,obj=param
kf=param.kpeak_aa
fac_h=1.263e-18
;
cwd,run
half='!s!u 1!n!r!s-!r!d 2!n'
file='specm.sav'
restore,file
;
;  determine normalization
;
k1=k & k1[0]=1. & k1=1./k1 & k1[0]=0.
EEM=total(spec1m)
kM=EEM/total(k1*spec1m)
kGW=sqrt(total(k^2*grav1m)/total(grav1m))
if cstress_prefactor eq '1' then begin
  EEGW0=EEM^2/(8.*kM^2)
endif else begin
  EEGW0=32.*!pi^2*EEM^2/kM^2
endelse
if lhalf_factor_in_GW then EEGW0=.5*EEGW0
spec1m0=EEM/kM
grav1m0=EEGW0/kM
print,spec1m0,kM
;kM=k
;
;grav1m=s*grav1m
;
default,kin,0
default,hdone,1
default,fac16pi,1.
default,maglabel,'mag'
default,green_lines,1
default,plot_polarization,1
default,yr_pspecm_scl2,[2d-5,2d1]
ytickf='logticks_exp'
yr=yr_pspecm_scl3
;
!p.title='!6';+run
!x.title='!6'
!y.title='!6!8h!6!s!d0!n!r!u2!n!7X!6!dGW!n(!8f!6), !6!8h!6!s!d0!n!r!u2!n!7X!6!d'+maglabel+'!n(!8f!6)'
freq=3.3e-3*2.*!pi/100.
freq=3.3e-3*2.*!pi/4500.
freq=3e-6
;
fac=4.1e-5
facGW=1./(32.*!pi)
xr=[8e-5,8e-1]
xr=[8e-6,8e-2]
circ_sym,0.7,1
plot_oo,freq*k/kM,fac16pi*facGW*fac*k*grav1m,xr=xr,yr=yr,ytickf=ytickf
if plot_polarization then begin
oplot,freq*k/kM,fac16pi*facGW*fac*k*grav2m,col=122,ps=8
oplot,freq*k/kM,-fac16pi*facGW*fac*k*grav2m,col=55,ps=8
endif
;
circ_sym,0.3,0
oplot,freq*k/kM,fac*k*spec1m
if plot_polarization then begin
oplot,freq*k/kM,.5*fac*k^2*spec2m,col=122,ps=8
oplot,freq*k/kM,-.5*fac*k^2*spec2m,col=55,ps=8
endif
;
xyouts,6e-4,1.4e-10,'!6GW',siz=siz
xyouts,2e-3,1.5e-7,'!6'+maglabel+'',siz=siz
;
xyouts,3e-7,2e-8,label
;
;@postproc
;
;xx=[3e-7,3e-3] & oplot,xx,1e-10/(xx/.001)
;xx=[3.2e-3,8e-2] & oplot,xx,3e-10*(xx/.01)^2
;xx=[3e-7,3e-3] & oplot,xx,1e-12/(xx/.001)
;xx=[3.2e-3,8e-2] & oplot,xx,3e-12*(xx/.01)^2
;
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,1.4e-3] & oplot,xx,1e-9*(xx/.001)^3.,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
;
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=2
;
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
;
file=dir+'Maggiore.rtf'
m1=rtable(file,2,head=1)
;
file=dir+'Caprini-C1.rtf'
a1=rtable(file,2,head=1)
oplot,10^a1[0,*],10^a1[1,*]
;
file=dir+'Caprini-C3.rtf'
a3=rtable(file,2,head=1)
oplot,10^a3[0,*],10^a3[1,*]
;
;  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=[1e-23,1e-15]
;
fach=1.263e-18/(freq*k/kM)
fach2=1.263e-18/(freq)
circ_sym,0.7,1
plot_oo,freq*k/kM,fach*sqrt(facGW*fac*k*grav1m),xr=xr,yr=yr,ytickf=ytickf
if plot_polarization then begin
oplot,freq*k/kM,fach*sqrt(facGW*fac*k*grav2m),col=122,ps=8
oplot,freq*k/kM,-fach*sqrt(facGW*fac*k*grav2m),col=55,ps=8
endif
;
;oplot,freq*k/kM,8e-16*sqrt(k*grah1m)/tm^2,co=188
alp=3./(8.*!pi)
oplot,freq*k/kM,alp*8e-16*sqrt(k*grah1m),co=188,li=2,thick=5
;
xyouts,3e-7,2e-8,label
;
;@postproc
;
;xx=[3e-7,3e-3] & oplot,xx,(fac_h/xx)*sqrt(1e-10/(xx/.001))
;xx=[3.2e-3,8e-2] & oplot,xx,(fac_h/xx)*sqrt(3e-10*(xx/.01)^2)
;xx=[3e-7,3e-3] & oplot,xx,(fac_h/xx)*sqrt(1e-12/(xx/.001))
;xx=[3.2e-3,8e-2] & oplot,xx,(fac_h/xx)*sqrt(3e-12*(xx/.01)^2)
;
;  continous limits
;
oplot,10^c1[0,*],10^c1[1,*],li=2
oplot,10^c2[0,*],10^c2[1,*],li=3
;oplot,10^m1[0,*],10^m1[1,*],li=1
oplot,10^a1[0,*],fac_h*sqrt(10^a1[1,*])/10^a1[0,*]
oplot,10^a3[0,*],fac_h*sqrt(10^a3[1,*])/10^a3[0,*]
;
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
;
fo="(a,3f6.1)"
print,'kM,kGW=',k0,kM,kGW,fo=fo
print,'$mv idl.ps ~/tex/roper/GW/fig/pspecm_scl3_'+run+'.ps
;
END