;$Id: pol_analytical.pro,v 1.4 2021/01/21 01:33:31 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
;
dir1='F1152d2_sig1_t11_M2c_double'
dir2='F1152sig07c'
dir3='F1152sig05a'
dir4='F1152sig03a'
dir5='F1152sig01a'

yr=[-1,1]
restore,'../'+dir1+'/specm.sav'
plot_oi,k[1:*],grav2m[1:*]/grav1m[1:*],yr=yr
restore,'../'+dir2+'/specm.sav'
oplot,k[1:*],grav2m[1:*]/grav1m[1:*]
restore,'../'+dir3+'/specm.sav'
oplot,k[1:*],grav2m[1:*]/grav1m[1:*]
restore,'../'+dir4+'/specm.sav'
oplot,k[1:*],grav2m[1:*]/grav1m[1:*]
restore,'../'+dir5+'/specm.sav'
oplot,k[1:*],grav2m[1:*]/grav1m[1:*]

;
; Kraichnan (1973) turbulence model based on Kolmogorov turbulence
; adding broken power law with k^2 spectrum for low wave numbers
; Kahniashvili, Gogoberidze and Ratra (2005) polarization model
;
dir_ext='r_Kolbr_nS_-3.67_nA_-4.67/'

yr=[-1,1]

dir='~/tex/roper/helical/data/'
dir=dir+dir_ext

h=(indgen(10) + 1)/10.
h=[0.1, 0.3, 0.5, 0.7, 1.]
h=string(h,format='(f4.2)')

k_ast=2*!DPI*100
factor = 1.
for i = 0, 4 do begin
    file=dir + 'h=' + h[i] + '.txt'
    c1=rtable(file,2,head=1)
    k2=c1[0,*]*k_ast
    P=c1[1,*]
    if (i eq 0) then begin
        oplot,k2,P*factor,col=122
        ;oplot,k2,P,thick=4,color=60
    endif else begin
        oplot,k2,P*factor,col=122
    endelse
endfor
;oplot,k2,P,thick=4,col=122

!x.title='!8k!6'
!y.title='!13P!6!dGW!n(!8k!6)'


END