;$Id: puniversal.pro,v 1.2 2026/05/21 10:25:25 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
  col0=0
endif else begin
  col0=255
endelse
;
;  This routine can also be used to determine table entries for CM and t_0.
;  In that case, we temporarily replace one run, e.g., the yellow one.
;
i1=21 & i2=29  ;(Orang2)
i1=12 & i2=17  ;(White)
i1=29 & i2=33  ;(Red)
i1=2 & i2=4  ;(Blue)
i1=15 & i2=22  ;(Green)
i1=32 & i2=35  ;(Orange)
!p.charsize=1.7
!x.margin=[7.7,.5]
!y.margin=[3.2,.3]
file='/kpm.txt'
siz=1.5
;
!p.multi=0
!x.title='!8t!6'
!y.title='!8C!d!7n!n!6  and  !8C!d!13E!n!6  '
fo='(f5.1,e10.1,2x,a)'
;
;  For Runs B2 and B3, we temporarily exchange them with Run C4u
;
dir1='1024_k1em2_k1em0_vA1em1_nu1em4b' & istart1=2+6*2  ;black
dir2='1024_k1em3_k1em1_vA1em1_nu1em3b' & istart2=2+6*3  ;blue
dir3='1024_k1em3_k1em1_vA7em3_nu1em4b' & istart3=2+6*3  ;red
dir4='8192_k1em2_k1em0_vA1em1_nu1em5' & istart4=2+6 ;green, running, CB=2.9-2.95
dir5='8192_k1em3_k1em1_vA9em3_nu1em5' & istart5=0+6*0 & i1_dir2=21 & i2_dir2=29   ;Orange
dir6='2048_k1em3_k1em1_vA1em1_nu5em4' & istart6=0+6*0 & i1_dir2=21 & i2_dir2=29   ;Yellow (exp Run B2)
dir6='4096_k1em3_k1em1_vA1em1_nu2em4' & istart6=0+6*0 & i1_dir2=21 & i2_dir2=29   ;Yellow (exp Run B3)
dir6='8192_k1em3_k1em1_vA9em3_nu1em5u2' & istart6=0+6*0 & i1_dir2=21 & i2_dir2=29 ;Yellow Run C4u
;
dir=dir1 & i1=istart1 & t1=3e4 & t2=3e5 & col=col0
dir=dir1 & i1=istart1 & t1=8e2 & t2=2e4 & col=col0
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
xr=[3e2,3e6] & yr=[1e-1,1.3e1] & plot_oo,t,Cxi,xr=xr,yr=yr,li=0 & oplot,t,CEE,li=2
;good=where(t ge t1 and t le t2) & tgood=t(good) & tAgood=tA(good)
;p=linfit(tgood,tAgood)
;xx=grange(t1,t2,300,/log)
;oplot,xx,xx*p[1]+p[0],col=col
;CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
;    i     t        EMag        xiMag        CM          Lu           IA         Cxi         CEE
;
dir=dir2 & i1=istart2 & t1=1e3 & t2=2e5 & col=55 ;should be 2
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
oplot,t,Cxi,li=0,col=col & oplot,t,CEE,li=2,col=col
;
dir=dir3 & i1=istart3 & t1=1e4 & t2=1.5e6 & col=122 ;should be 3
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
oplot,t,Cxi,li=0,col=col & oplot,t,CEE,li=2,col=col
;
loadct,6
dir=dir4 & i1=istart4 & t1=6e2 & t2=7e4 & col=144
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
oplot,t,Cxi,li=0,col=col & oplot,t,CEE,li=2,col=col
loadct,5
;
dir=dir5 & i1=istart5 & t1=1e4 & t2=2e5 & col=155
dir=dir5 & i1=istart5 & t1=7e3 & t2=3e5 & col=155
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
oplot,t,Cxi,li=0,col=col & oplot,t,CEE,li=2,col=col
;
dir=dir6 & i1=istart6 & t1=1e3 & t2=2e5 & col=188  ;(for Run B2)
dir=dir6 & i1=istart6 & t1=1e3 & t2=7e4 & col=188  ;(for Run B3)
dir=dir6 & i1=istart6 & t1=7e3 & t2=3e5 & col=188  ;(for Run C4u)
a=rtable('../'+dir+file,head=1,9)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA & Cxi=a[7,*] & CEE=a[8,*] 
oplot,t,Cxi,li=0,col=col & oplot,t,CEE,li=2,col=col
;
oplot,xr,xr*0+9.5,li=3
oplot,xr,xr*0+4.2,li=3
oplot,xr,xr*0+.30,li=3
oplot,xr,xr*0+.19,li=3
;
siz=1.9
xyouts,siz=siz,500,6.0,'!8C!d!13E!n!6'
xyouts,siz=siz,7e5,.23,'!8C!d!7n!n!6'
siz=1.4
xyouts,siz=siz,1e4,.14,'!6Lu=1100'
xyouts,siz=siz,1e4,7.0,'!6Lu=1100'
xyouts,siz=siz,7e5,.36,'!6Lu=70'
xyouts,siz=siz,7e5,2.8,'!6Lu=70'
;
print,"$mv idl.ps ~/Overleaf/Mattia/Isochrones/fig/puniversal.eps"
END