;$Id: pkpm8192a_hel.pro,v 1.2 2026/06/15 05:54:58 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=11,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
!p.charsize=1.7
!x.margin=[8.6,.5]
!y.margin=[3.2,.3]
;
siz=2.0
xr=[1.,40.]
yr=[4e-4,1e-2]
!x.title='!7n!6!dM!n(!8t!6)'
!y.title='!8v!6!dA!n(!8t!6)'
lstop=1
lstop=0
;
file='/kpm.txt'
ss1=.5
ss2=1.3
fo='(a,2x,9e11.1)'
;
;-----------------------------------------------------------------------------
dir1='8192_k1em3_k1em1_vA9em3_nu5em5_hel'       ;1/2 decay
dir2='8192_k1em3_k1em1_vA9em3_nu1em5u3_k0_hel2' ;1/2 decay with excursion
;
dir3='8192_k1em3_k1em1_vA9em3_nu1em5' 
dir4='8192_k1em3_k1em1_vA9em3_nu1em5u3_k0'
;
;-----------------------------------------------------------------------------
dir=dir1
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
tA=xiM/vA
plot_oo,xiM,vA,xr=xr,yr=yr,/nodata
oplot,xiM,vA,col=122
circ_sym,1.6,0
oplot,[1,1]*xiM[0],[1,1]*vA[0],ps=8
;
dir=dir2
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
tA=xiM/vA
oplot,xiM,vA,col=122,li=2
circ_sym,1.0,0
oplot,[1,1]*xiM[0],[1,1]*vA[0],col=122,ps=8
;
xx=[5.,30.] & oplot,xx,1.3e-2/xx^.5,li=3
xyouts,16.,.004,siz=siz,'!9A!7n!6!s!dM!n!r!u-1/2!n'
;
;-----------------------------------------------------------------------------
dir=dir3
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
tA=xiM/vA
oplot,xiM,vA,col=55
circ_sym,1.6,0
oplot,[1,1]*xiM[0],[1,1]*vA[0],ps=8
;
dir=dir4
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
tA=xiM/vA
oplot,xiM,vA,col=55,li=2
circ_sym,1.0,0
oplot,[1,1]*xiM[0],[1,1]*vA[0],col=55,ps=8
;
xx=[7.,20.] & oplot,xx,1.5e-2/xx,li=3
xyouts,13.,.0013,siz=siz,'!9A!7n!6!s!dM!n!r!u-1!n'
;
print,"$mv idl.ps ~/Overleaf/Mattia/Isochrones/fig/pkpm8192a_hel.eps"
print,"$mv idl.ps ~/Overleaf/Oksana/MHD-in-axion-inflation/fig/pkpm8192a_hel.eps"
;
END