;$Id: pkpm1024a_t.pro,v 1.10 2026/06/07 11:40:08 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=20,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
!p.multi=[0,1,2]
!p.charsize=1.7
!x.margin=[8.4,.5]
!y.margin=[3.2,.4]
oplot_fit=0
;
xr=[1e1,2.3e6]
siz=1.7
si2=1.4
;
file='/kpm.txt'
;
dir1='1024_k1em2_k1em0_vA1em1_nu1em4b'  ;RunA1
dir2='1024_k1em3_k1em1_vA1em1_nu1em3b'  ;RunB1
dir3='1024_k1em3_k1em1_vA7em3_nu1em4b'  ;RunC1
;
tA=0.
s=sqrt(10.)
;
;-----------------------------------------------------------------------------
!x.title='!6'
!y.title='!7n!6!dM!n(!8t!6)'
yr=[7e-1,1.4e2]
a=rtable('../'+dir1+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
plot_oo,t+tA,xiM,xr=xr,yr=yr
xyouts,14.,77.,siz=siz,'!6(a)'
;
a=rtable('../'+dir2+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,t,xiM,col=55
oplot,t,xiM/s,col=55,li=2
xx=5e5 & arrow,/data,col=55,xx,103.,xx,26.
xyouts,6e5,43.,siz=si2,col=55,'!9S!610'
;
a=rtable('../'+dir3+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,t,xiM,col=122
;
xx=[8e2,1e6] & oplot,xx,.028*xx^.5,li=1
xx=[1e1,1e5] & oplot,xx,7.5*xx^.0,li=1,col=122
xx=[1e1,1e4] & oplot,xx,2.3*xx^.0,li=1,col=55
;
circ_sym,1.0,1
xx=7e2 & oplot,[1,1]*xx,[1,1]*.028*xx^.5,ps=8
xx=7e4 & oplot,[1,1]*xx,[1,1]*.028*xx^.5,ps=8,col=122
xx=7e3 & oplot,[1,1]*xx,[1,1]*.028*xx^.5,ps=8,col=55
;
if oplot_fit then begin
  xx=grange(1e1,2e6,100,/log)
  oplot,xx,0.7*(1.+xx/7e2)^.5,li=1
  oplot,xx,7.5*(1.+xx/7e4)^.5,li=1,col=122
  oplot,xx,2.3*(1.+xx/7e3)^.5,li=1,col=55
endif
;-----------------------------------------------------------------------------
fact=sqrt(.1)
!x.title='!8t!6'
!y.title='!8v!6!dA!n(!8t!6)'
yr=[2e-4,2e-2]
a=rtable('../'+dir1+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
plot_oo,t+tA,vA,xr=xr,yr=yr
xyouts,8e5,.012,siz=siz,'!6(b)'
;
a=rtable('../'+dir2+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,t,vA,col=55
oplot,t,vA/s,col=55,li=2
xx=5e5 & arrow,/data,col=50,xx,2.0e-3,xx,5e-4
xyouts,6e5,5e-4,siz=si2,col=55,'!9S!610'
;
a=rtable('../'+dir3+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,t,vA,col=122
;
xx=[3e2,1e6] & oplot,xx,.28/xx^.5,li=1
xx=[1e1,4e2] & oplot,xx,.01700/xx^.0,li=1
xx=[1e1,1e5] & oplot,xx,.00116/xx^.0,li=1,col=122
xx=[1e1,5e3] & oplot,xx,.00540/xx^.0,li=1,col=55
;
circ_sym,1.0,1
xx=2.8e2 & oplot,[1,1]*xx,[1,1]*.28/xx^.5,ps=8
xx=5.8e4 & oplot,[1,1]*xx,[1,1]*.28/xx^.5,ps=8,col=122
xx=2.6e3 & oplot,[1,1]*xx,[1,1]*.28/xx^.5,ps=8,col=55
;
if oplot_fit then begin
  xx=grange(1e1,2e6,100,/log)
  oplot,xx,.01700/(1.+xx/2.8e2)^.5,li=1
  oplot,xx,.00116/(1.+xx/6.3e4)^.5,li=1,col=122
  oplot,xx,.00540/(1.+xx/2.6e3)^.5,li=1,col=55
endif
;
!p.multi=0
print,"$mv idl.ps ~/Overleaf/Mattia/Isochrones/fig/pkpm1024a_t.eps"
END