;$Id: pfit_all.pro,v 1.6 2026/05/24 20:14:57 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=[6.4,.5]
!y.margin=[3.2,.3]
file='/kpm.txt'
siz=1.5
;
!p.multi=0
!x.title='!8t!6'
!y.title='!8t!6!dA!n(!8t!6)'
fo='(f5.1,e10.1,2x,a)'
fo2='(a8,2e9.1,f7.3)'
openw,1,'fit.tmp'
;
;  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=1e3 & t2=2e4 & col=col0
dir=dir1 & i1=istart1 & t1=5e2 & t2=2e4 & col=col0
dir=dir1 & i1=istart1 & t1=2e2 & t2=2e4 & col=col0
dir=dir1 & i1=istart1 & t1=2e2 & t2=1e4 & col=col0
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
xr=[1e2,3e6] & yr=[3e1,3e5] & plot_oo,t,tA,xr=xr,yr=yr,li=1
;xr=[0,2.3e6] & yr=[0,2.3e5] & plot,t,tA,xr=xr,yr=yr,li=1
;xr=[0,1.3e6] & yr=[0,2.3e5] & plot,t,tA,xr=xr,yr=yr,li=1
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,t(good),t(good)*p[1]+p[0],col=col,ps=-1
;printf,1,'black: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
;
dir=dir2 & i1=istart2 & t1=8e2 & t2=1.4e5 & col=55 ;should be 2
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
oplot,t,tA,col=col,li=1
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,t(good),t(good)*p[1]+p[0],col=col,ps=-1
;printf,1,'blue: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
;
dir=dir3 & i1=istart3 & t1=1e4 & t2=1.0e6 & col=122 ;should be 3
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
oplot,t,tA,col=col,li=1
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,t(good),t(good)*p[1]+p[0],col=col,ps=-1
;printf,1,'red: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
;
loadct,6
dir=dir4 & i1=istart4 & t1=2e2 & t2=8e4 & col=144
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
oplot,t,tA,col=col,li=1
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
oplot,t(good),t(good)*p[1]+p[0],col=col,ps=-1
printf,1,'green: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
loadct,5
;
dir=dir5 & i1=istart5 & t1=7e3 & t2=5e5 & col=155
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
oplot,t,tA,col=col,li=1
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,t(good),t(good)*p[1]+p[0],col=col,ps=-1
;printf,1,'orange: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
print,max(t)
;
dir=dir6 & i1=istart6 & t1=5e2 & t2=2e5 & col=188  ;(for Run B2)
dir=dir6 & i1=istart6 & t1=7e2 & t2=7e4 & col=188  ;(for Run B3)
dir=dir6 & i1=istart6 & t1=7e3 & t2=5e5 & col=188  ;(for Run C4u)
a=rtable('../'+dir+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & tA=xiM/vA
oplot,t,tA,col=col,li=1
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,t(good),t(good)*p[1]+p[0],col=col,ps=-1
;printf,1,'yellow: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
;printf,1,'B2: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
;printf,1,'B3: ',t1,t2,mean(((tA(good)-(t(good)*p[1]+p[0]))/tA(good))^2)^.5,fo=fo2
CM=1./p[1] & print,CM,CM*p[0],dir,fo=fo
;
close,1
spawn,'cat fit.tmp'
spawn,'cat fit.tmp >> ../idl/fit.txt'
print,"$mv idl.ps ~/Overleaf/Mattia/Isochrones/fig/pfit_"+dir+".eps"
END