;$Id: pspecm_scl2.pro,v 1.10 2018/04/26 16:16:29 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=11,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
;
siz=1.6
@parameters
!p.charsize=1.8
!p.charsize=2.4
!x.margin=[8.6,0.5]
!x.margin=[3.6,0.5]
!y.margin=[2.2,0.2]
;
cwd,run
half='!s!u 1!n!r!s-!r!d 2!n'
circ_sym,0.7,1
file='specm.sav'
restore,file
;
;  determine normalization
;
k1=k & k1[0]=1. & k1=1./k1 & k1[0]=0.
EEM=total(spec1m)
HHM=total(spec2m)
EEGW=total(grav1m)
EEGW2=total(grav2m)
kM=EEM/total(k1*spec1m)
frac_hel=.5*kM*HHM/EEM
print,'polarization,frac_hel=',EEGW2/EEGW,frac_hel
kGW=EEGW/total(k1*grav1m)
kGW2=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
;
!p.title='!6';+run
;!x.title='!8k!6/!8k!6!d0!n'
!x.title='!6'
;!y.title='!8E!6(!8k!6)  and  !8f!6(!8k!6)'
;!y.title='!8E!6(!8k!6)'
!y.title='!6'
xr=[.9,max(k)*1.1]/kM
default,kin,0
default,hdone,1
default,yr_pspecm_scl2,[2d-5,2d1]
print,"yr_pspecm_scl2=[2d-5,2d6]"
ytickf='logticks_exp'
yr=yr_pspecm_scl2
plot_oo,k/kM,spec1m/spec1m0,xr=xr,yr=yr,ytickf=ytickf
if not kin then begin
  oplot,k/kM,0.5*k*abs(spec2m)/spec1m0,li=1,col=188
  oplot,k/kM,+.5*k*spec2m/spec1m0,ps=8,col=122
  oplot,k/kM,-.5*k*spec2m/spec1m0,ps=8,col=55
  oplot,k/kM,kine1m/spec1m0,li=2
endif
;oplot,k/kM,0.5*k1*kine2m/spec1m0,li=1
;oplot,k/kM,+.5*k1*kine2m/spec1m0,ps=8,col=122
;oplot,k/kM,-.5*k1*kine2m/spec1m0,ps=8,col=55
oplot,k/kM,grav1m/grav1m0,li=3
oplot,k/kM,0.5*abs(grav2m)/grav1m0,li=1,col=188
oplot,k/kM,+.5*grav2m/grav1m0,ps=8,col=122
oplot,k/kM,-.5*grav2m/grav1m0,ps=8,col=55
;oplot,k/kM,k*0+1.,li=3
@postproc
;
fo="(a,4f6.1)"
print,'kM,kGW,kGW2=',k0,kM,kGW,kGW2,fo=fo
print,'$mv ,kGW,fo=fo
print,'$mv idl.ps ~/tex/roper/GW/fig/pspecm_scl2_'+run+'.ps
;
;plot_oo,k[1:*],grav1m[1:*]/spec1m[1:*]*k^2*1e4,yr=[1e-4,1e+1]
;plot_oo,k[1:*],grav1m[1:*]/spec1m[1:*]*k^2*30^2,yr=[1e-4,1e+1]
;plot_oo,k[1:*],grav1m[1:*]/spec1m[1:*]*k^2*10^2,yr=[1e-4,1e+1]
END