;$Id: pom2.pro,v 1.16 2026/07/31 10:45:53 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=18,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
  col0=0
endif else begin
  col0=255
end
;
@parameters
!p.charsize=1.7
!x.margin=[7.6,.5]
!y.margin=[3.2,.5]
!p.multi=[0,1,2]
pc_read_ts,obj=ts
pc_read_dim,obj=dim
pc_read_param,obj=param1
pc_read_param,/param2,obj=param
ampluu=param1.ampluu
k1=2.*!pi/(param1.xyz1[0]-param1.xyz0[0])
k0=param1.kpeak*k1
lupw=param.lupw_uu  
nu=param.nu
nx=dim.nx
circ_sym,2.9,1
;
;  panel 1
;
!x.title='!6'
!y.title='!6values'
!y.title='!6terms in Eq.(12)'
default,xr_panel1,[0,3.3]
;default,yr_panel1,[-.01,.1]
default,yr_panel1,[-.01,.2]
!x.range=xr_panel1
;
default,t1,1.3
default,t2,2.7
default,t1d,0.5  ;(for determining early peak values and dyn phases)
default,t2d,1.5  ;(for determining early peak values; changed from 0.5 to 1.5)
default,t1e,0.3
default,t2e,1.0
default,t1l,2.5
default,t2l,3.0
good=where(ts.t ge t1 and ts.t le t2)      ;(red, Rgen1)
early=where(ts.t le t2d)               
early2=where(ts.t ge t1e and ts.t le t2e)  ;(red, Rgen2)
dynphase=where(ts.t ge t1d and ts.t le t2d)
late=where(ts.t ge t1l and ts.t le t2l)
;
genmax=2.*nu*max(ts.qsglnrhom(early))  ;(red)
dismax=   nu*max(ts.q2m(early))        ;(green)
dynmax=      max(ts.quxom(early))      ;(blue)
;
igenmax=where(genmax eq 2.*nu*ts.qsglnrhom(early))
idismax=where(dismax eq    nu*ts.q2m(early))
idynmax=where(dynmax eq       ts.quxom(early))
;
plot,ts.t,deriv(ts.t,.5*ts.orms^2),yr=yr_panel1
oplot,ts.t,ts.t*0
oplot,ts.t,2*nu*ts.qsglnrhom,col=122
oplot,ts.t[igenmax]*[1,1],2*nu*ts.qsglnrhom[igenmax]*[1,1],col=122,ps=8
loadct,6
oplot,ts.t,nu*ts.q2m,col=144
oplot,ts.t[idismax]*[1,1],nu*ts.q2m[idismax]*[1,1],col=144,ps=8
loadct,5
oplot,ts.t,ts.quxom,col=55,th=4
oplot,ts.t[idynmax]*[1,1],ts.quxom[idynmax]*[1,1],col=55,th=4,ps=8
oplot,ts.t,ts.quxom-nu*ts.q2m+2*nu*ts.qsglnrhom,col=155,li=2,th=4
;
siz=1.7
xyouts,siz=siz,0.3,.086,'!6(a)'
xyouts,siz=siz,0.76,0.038,'!8T!6!dgen!n',col=122
xyouts,siz=siz,0.9,0.016,'!8T!6!ddyn!n',col=55
loadct,6
xyouts,siz=siz,1.1,0.05,'!8T!6!ddis!n',col=144
loadct,5
xyouts,siz=siz,0.2,0.02,'!6net',col=155
;
;-----------------------------------------------------------------------------
;  panel 2
;
!x.title='!8t!6'
!y.title='!6ratios'
default,loveride_ratm,0    ;(Rgen1)
default,loveride_ratm2,0   ;(Rgen2)
rat=2*ts.qsglnrhom/ts.q2m
ratd=ts.quxom/(nu*ts.q2m)
ratdx=max(ratd(late))
plot,ts.t,rat,yr=[0,1.2]
oplot,ts.t,ratd,li=2,th=3
tgood=ts.t(good)
tearly=ts.t(early)
tdyn=ts.t(dynphase)
ratm=mean(rat(good))
ratm2=max(rat(early2))
ratdm=mean(ratd(dynphase))  ;(dyn1)
it=findex(ratm2,rat(early2),/rev)
tearly2=ts.t(early2)
if loveride_ratm2 then oplot,tearly2[it]*[1,1],[0,0]+ratm2,col=122,ps=1 else oplot,tearly2[it]*[1,1],[0,0]+ratm2,col=122,ps=8
;
it=findex(ratdx,ratd(late),/rev)
nt=n_elements(ratd(late))
it=it<(nt-1)
print,'AXEL: ',it
tlate=ts.t(late)
oplot,tlate[it]*[1,1],[0,0]+ratdx,col=55,ps=8
;
if loveride_ratm then oplot,tgood,tgood*0+ratm,col=122,li=2 else oplot,tgood,tgood*0+ratm,col=122
oplot,tdyn,tdyn*0+ratdm,col=55
xyouts,siz=siz,3.02,1.07,'!6(b)'
xyouts,siz=siz,0.46,0.86,'!8R!6!dgen2!n',col=122
xyouts,siz=siz,1.8,0.80,'!8R!6!dgen1!n',col=122
xyouts,siz=siz,0.8,0.29,'!8R!6!ddyn1!n',col=55
xyouts,siz=siz,2.4,0.40,'!8R!6!ddyn2!n',col=55
;
;  compute epsK, etc, at t=tcrit=2.7
;
tcrit=2.7
its=findex(tcrit,ts.t)
rhom=ts.rhom(its)
knu=(ts.epsK[its]/(rhom*nu^3))^.25
dx=4.*!pi/nx
kNy=!pi/dx
print,ts.epsK[its]/k0,ts.divu2m[its]/k0^2,ts.orms[its]^2/k0^2,knu,kNy
;
wait,.2
cwd,run
!p.multi=0
if loveride_ratm  then ratm =99.999
if loveride_ratm2 then ratm2=99.999
fo='(e7.1,2f7.3,4e10.1,i6,3e10.1,f6.1,i6,f7.2,f7.2,f7.2,f6.1,i4,2x,a)'
openw,1,'res.tmp'
printf,1,nu,ratm,ratm2,ratdm,genmax,dismax,dynmax,nx,ts.epsK[its],ts.divu2m[its]^.5,ts.orms[its],knu,kNy,ampluu[0],ratdx,ts.umax[its],k0,lupw,run,fo=fo
close,1
spawn,'cat res.tmp'
spawn,'cat res.tmp >> ../idl/res.txt'
cwd,run
print,"$mv idl.ps ~/Overleaf/Eva/BE_collapse/fig/pom2_"+run+".eps"
END