;$Id: ptxi.pro,v 1.1 2026/06/08 06:38:56 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=2.7
!x.margin=[7.8,.5]
!y.margin=[3.2,2.5]
;
;  procedure to determine fit coefficients listed in txi.txt
;
xr=[1.,100]
xr=[.6,200]
yr=[1e-5,1e-1]
;
file='kpm.txt'
a=rtable(file,head=1,8)
@parameters
;
;  i     t        EMag        xiMag         CM          Lu          IA          Cxi         CEE
;
xr=[1,1e7]
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*] & IA=a[6,*]
!p.multi=[0,2,2]
!x.title='!6'
!y.title='!8I!6!dA!n'
plot_oi,t,IA/(vA*xiM)^2,xr=xr
IArel=mean(IA/(vA*xiM)^2)
oplot,t,t*0+IArel,col=122
;
!x.title='!7n!6!dM!n(!8t!6)'
Cxi=min(xiM/(t^.5*IA^.25))
txi=min(xiM^2/(Cxi^2*IA^.5))
yr=[.1,100]
plot_oo,t,xiM,ps=-1,xr=xr,yr=yr
oplot,t,Cxi*IA^.25*(txi+t)^.5,col=122
;
!x.title='!8t!6'
!y.title='!8v!6!dA!n(!8t!6)'
CvA=max(vA*t^.5/IA^.25)
tvA=min(CvA^2*IA^.5/vA^2)
;
yr=[1e-4,1e-1]
plot_oo,t,vA,ps=-1,xr=xr,yr=yr
oplot,t,CvA*IA^.25/(tvA+t)^.5,col=122
;
yr=[1e-8,1e-1]
plot_oo,t,.5*vA^2,ps=-1,xr=xr,yr=yr
CE=.5*CvA^2
oplot,t,CE*IA^.5/(tvA+t),col=122
;
cwd,run
fo='(2f8.3,2i8,f8.2,2x,a,2x,a)'
openw,1,'tmp.txt'
printf,1,Cxi,CvA,nint(txi),nint(tvA),IArel,crun,run,fo=fo
close,1
spawn,'cat tmp.txt'
spawn,'cat tmp.txt >> ../idl/txi.txt'
!p.multi=0
END