;$Id: pkpm8192a_pq.pro,v 1.3 2026/05/21 11:42:28 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=12,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
!p.charsize=1.7
!x.margin=[9.4,.5]
!y.margin=[3.2,.3]
;
xr=[0.,1.] & yr=[2.,0.]
xr=[0.,.66] & yr=[0.,1.4]
!x.title='!8q!6(!8t!6)'
!y.title='!8p!6(!8t!6)'
plot,xr,xr*2,li=3,yr=yr
;oplot,xr,1.*xr,li=1
;oplot,xr,2.5*xr,li=1
;oplot,xr,3.*xr,li=1
;oplot,xr,5.*xr,li=1
oplot,xr,2.*(1.-xr)
lstop=1
lstop=0
;
file='/kpm.txt'
ss1=.5
ss2=1.3
fo='(9e11.1)'
siz=1.7
;
;-----------------------------------------------------------------------------
dir1='1024_k1em2_k1em0_vA1em1_nu1em4b' & istart1=2+6*2  ;black
dir2='1024_k1em3_k1em1_vA7em3_nu1em4b' & istart2=2+6*3  ;red
dir3='1024_k1em3_k1em1_vA1em1_nu1em3b' & istart3=2+6*3  ;blue
dir4='8192_k1em2_k1em0_vA1em1_nu1em5' & istart4=2+6 ;green, running, CB=2.9-2.95
;dir5='8192_k1em3_k1em1_vA9em3_nu1em5' & istart5=2+6*3  ;Orange
dir5='8192_k1em3_k1em1_vA9em3_nu1em5u2' & istart5=0+6*0 & i1_dir2=21 & i2_dir2=29  ;Orange
;
;  i   t      EMag       xiMag        CM          Lu
;  & CM=a[4,*] & Lu=a[5,*]
;
default,CM1,10.
default,CM4,22.
dir=dir1 & i1=istart1
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
q=+deriv(alog(t),alog(xiM))
p=-deriv(alog(t),alog(vA^2/2.))
circ_sym,ss1,0 & oplot,q,p,ps=8
tt=every(t[0,i1:*],6) & qq=every(q[0,i1:*],6) & pp=every(p[0,i1:*],6) & print,tt,fo=fo
circ_sym,ss2,1 & oplot,qq,pp,ps=8
if lstop then stop
;
dir=dir2 & i1=istart2
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
q=+deriv(alog(t),alog(xiM))
p=-deriv(alog(t),alog(vA^2/2.))
circ_sym,ss1,0 & oplot,q,p,ps=8,col=122
tt=every(t[0,i1:*],6) & qq=every(q[0,i1:*],6) & pp=every(p[0,i1:*],6) & print,tt,fo=fo
circ_sym,ss2,1 & oplot,qq,pp,ps=8
if lstop then stop
;
dir=dir3 & i1=istart3
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
q=+deriv(alog(t),alog(xiM))
p=-deriv(alog(t),alog(vA^2/2.))
circ_sym,ss1,0 & oplot,q,p,ps=8,col=55
tt=every(t[0,i1:*],6) & qq=every(q[0,i1:*],6) & pp=every(p[0,i1:*],6) & print,tt,fo=fo
circ_sym,ss2,1 & oplot,qq,pp,ps=8,col=55
xyouts,100.,7e-4,siz=siz,'!8t!6=10!u6!n',col=55
if lstop then stop
;
dir=dir4 & i1=istart4
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
loadct,6
q=+deriv(alog(t),alog(xiM))
p=-deriv(alog(t),alog(vA^2/2.))
circ_sym,ss1,0 & oplot,q,p,ps=8,col=144
tt=every(t[0,i1:*],6) & qq=every(q[0,i1:*],6) & pp=every(p[0,i1:*],6) & print,tt,fo=fo
circ_sym,ss2,1 & oplot,qq,pp,ps=8,col=144
loadct,5
;
dir=dir5 & i1=istart5
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
q=+deriv(alog(t),alog(xiM))
p=-deriv(alog(t),alog(vA^2/2.))
circ_sym,ss1,0 & oplot,q,p,ps=8,col=155
tt=every(t[0,i1:*],6) & qq=every(q[0,i1:*],6) & pp=every(p[0,i1:*],6) & print,tt,fo=fo
circ_sym,ss2,1 & oplot,qq,pp,ps=8,col=155
;
print,"$mv idl.ps ~/Overleaf/Mattia/Isochrones/fig/pkpm8192a_pq.eps"
END