;$Id: pkpm_vs_lam.pro,v 1.21 2026/10/27 06:10:30 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.5
!x.margin=[8.8,6.4]
!y.margin=[3.2,2.6]
;
tilde='!9!s!aA!n!r!6'
xr=[.002,.8] & yr=[6e-4,6e-1]
!x.title='!7n!6!dM!n(!8t!6)'
!y.title='!8v!6!dA!n(!8t!6)'
!y.title='!8B!6(!8t!6)'
;
file='/kpm.txt'
;
;-----------------------------------------------------------------------------
dir1='2048_k1_e1em6_k60c_h1em14'  ;black  (dt=4e-4)
dir2='2048_k1_e1em6_k180a_h1em14' ;red
dir3='2048_k1_e1em6_k180b_h1em14' ;blue
dir4='2048_k1_e1em6_k60a_h1em14'  ;orange
;dir5='2048_k1_e1em7_k60a_h1em14'  ;brown  (dt=3e-5) --> to suppress
;dir6='2048_k1_e1em8_k60a_h2em13'  ;green  (dt=3e-5)
dir6='2048_k1_e1em8_k60a_h1em14'  ;green  (dt=3e-5)
;
;  i   t      EMag       xiMag        CM          Lu
;  & CM=a[4,*] & Lu=a[5,*]
;
tcrit=100.
a=rtable('../'+dir1+file,head=1,6)
t=reform(a[1,*]) & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
plot_oo,xiM,vA,ps=1,xr=xr,yr=yr,xst=9,yst=9
circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8
;
a=rtable('../'+dir2+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,xiM,vA,ps=1,col=122
circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8,col=122
;
a=rtable('../'+dir3+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,xiM,vA,ps=1,col=55
circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8,col=55
print,'blue: max(t),tcrit=',max(t),tcrit,t(it)
;
a=rtable('../'+dir4+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,xiM,vA,ps=1,col=155
circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8,col=155
print,'orange: max(t),tcrit=',max(t),tcrit,t(it)
;
loadct,6
;a=rtable('../'+dir5+file,head=1,6)
;t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
;oplot,xiM,vA,ps=1,col=88
;circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8,col=88
;print,'yellow: max(t),tcrit=',max(t),tcrit,t(it)
;
a=rtable('../'+dir6+file,head=1,6)
t=a[1,*] & vA=sqrt(2.*a[2,*]) & xiM=a[3,*]
oplot,xiM,vA,ps=1,col=144
circ_sym,1.3,1 & it=findex(tcrit,t) & oplot,xiM[it]*[1,1],vA[it]*[1,1],ps=8,col=144
loadct,5
print,'green: max(t),tcrit=',max(t),tcrit,t(it)
;
xx=[.06,.4] & oplot,xx,21*xx^2/tcrit,li=3
xx=[.06,.8] & oplot,xx,1*xx^2/tcrit,li=1
xx=[.002,.5] & oplot,xx,6.7e-4/xx,li=0
xyouts,.025,2e-3,'!8B!6=21!7n!6!s!dM!n!r!u2!n/!8t!6'
xyouts,.160,8e-4,'!7n!6!s!dM!n!r!u2!n/!8t!6'
;
!x.title='                !8B!6!drms!n [G]'
!x.title='!8B!6!drms!n [G]'
!y.title='!8L!6 [erg s!u-1!n]'
!x.title='!7n!6!dM!n(!8t!6) [km]'
!y.title='!8B!6(!8t!6) [G]'
;axis,yax=1,yr=yr*EMfac*Vol
epsunit=2.5d9
xfact=100.
yfact=1.6e14*xfact*30.
axis,yax=1,yr=yr*yfact
axis,xax=1,xr=xr*xfact
;
;-----------------------------------------------------------------------------
;
print,"$mv idl.ps ~/Overleaf/Vinay/Helicity-Gauge/figs/pkpm_vs_lam.eps"
print,"$mv idl.ps ~/Overleaf/Vinay/Helicity-Gauge/figs/pkpm_vs_lam_eta.eps"
END