;$Id: ppspec_save_B5.pro,v 1.8 2021/01/17 05:18:44 brandenb Exp $
if !d.name eq 'PS' then begin
  ;device,xsize=18,ysize=8,yoffset=3
  device,xsize=18,ysize=10,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
literate=0
@parameters
!p.charsize=1.6
!x.margin=[8.2,0.5]
!y.margin=[3.1,0.4]
!p.multi=0
si2=1.4
siz=1.6
;
file='spec_range.sav'
dir1='512_1e2_1e4_4e8_1em6a' & i1=21 & fac1=1.
dir2='512_1e2_1e4_4e8_2em6a' & i2=5 & fac2=1.
dir3='512_1e2_1e4_4e8_5em6a' & i3=2 & fac3=1.
dir4='512_1e2_1e4_4e8_1em5a' & i4=7 & col4=55 & fac4=1.
dir7='512_1e2_1e4_4e8_1em4a2'& i7=5 & col7=122 & fac7=1.4
;dir8='512_1e2_1e4_4e8_2em4a' & i8=1 & col8=122 & fac8=1.
;dir9='512_1e2_1e4_4e8_5em4a' & i9=1 & col9=122 & fac9=13000.
dir0='512_1e2_1e4_4e8_1em3a' & i0=1 & col0=122 & fac0=1.0
dir5='512_1e2_1e4_4e8_2em5a' & i5=1 & col5=122 & fac5=3.2
dir6='512_1e2_1e4_4e8_5em5a' & i6=1 & col6=155 & fac6=1.
;
!x.title='!8k!6'
!y.title='!8E!6!dGW!n(!8k,t!6)'
xr=[100.,14000.]
yr=[3e-27,3e-22] ;(for fits only)
yr=[3e-27,3e-24] ;(for the paper)
restore,'../'+dir3+'/'+file & i=i3
plot_oo,k[1:*],k(1:*),yr=yr,/nodata,xr=xr
;
;  choices of different lines
;
dir=dir4 & i=i4 & fac=fac4 & kGW2=476.69 & fact=3.845e-19
dir=dir3 & i=i3 & fac=fac3 & kGW2=226.50 & fact=2.958e-20
dir=dir2 & i=i2 & fac=fac2 & kGW2=104.38 & fact=3.955e-21
dir=dir1 & i=i2 & fac=fac2 & kGW2= 72.46 & fact=4.294e-22
dir=dir6 & i=i6 & fac=fac6 & kGW2=2888.0 & fact=3.082e-16
dir=dir7 & i=i7 & fac=fac7 & kGW2=8244.8 & fact=1.3656e-14
dir=dir0 & i=i0 & fac=fac0 & kGW2=50000.0& fact=3.7114572e-11
dir=dir5 & i=i5 & fac=fac5 & kGW2=982.34 & fact=5.328e-18
;
loadct,6
restore,'../'+dir+'/'+file & i=i
spec=fac*grav1_range(*,i)
oplot,k,spec,col=col5
kGW=total(grav1_range(1:*,i))/total(grav1_range(1:*,i)/k[1:*]) & print,'kGW=',kGW,'  ',dir5
loadct,5
mu50=1e4
;
kGW=kGW2
spec_theo=fact*exp(-(k/mu50)^4)*(k/(kGW^2+k^2))^2
oplot,k,spec_theo,li=2
;
if literate then begin
  nx=30 & ny=30 & xfac=.1 & yfac=.1
  nx=30 & ny=30 & xfac=.01 & yfac=.01
  nx=30 & ny=30 & xfac=.0005 & yfac=.0005
  nx=30 & ny=30 & xfac=1. & yfac=800.
  ;
  res=fltarr(nx,ny)
  kGW_=fltarr(nx)
  fact_=fltarr(ny)
  x=grange(-1.,1.,nx)*xfac
  y=grange(-1.,1.,ny)*yfac
  kGW_=kGW2*exp(x)
  fact_=fact*exp(y)
  ;
  for iy=0,ny-1 do begin
  for ix=0,nx-1 do begin
    kGW=kGW_[ix]
    fact=fact_[iy]
    spec_theo=fact*exp(-(k/mu50)^4)*(k/(kGW^2+k^2))^2
    oplot,k,spec_theo,li=2
    print,'kGW=',kGW
    help,spec,spec_theo
    res(ix,iy)=sqrt(mean((1d24*(spec-spec_theo))^2))
    print,res(ix,iy)
  endfor
  endfor
  contour,res,nlev=100,kGW_,fact_
endif
;
loadct,6
xyouts,6800.,1.7e-25,'!6B5',siz=siz,col=122 & loadct,5
xyouts,5000.,1.9e-26,'!8E!6!s!dGW!n!r!umodel!n',siz=siz
;
cwd,run
print,'$mv idl.ps ~/tex/rei/GW/fig/ppspec_save_B5.ps'
END