;$Id: ppxyaver_parametric_combine.pro,v 1.1 2020/08/02 17:24:41 brandenb Exp $
if !d.name eq 'PS' then begin
  ysize=11
  default,thick,4
  device,xsize=18,ysize=ysize,yoffset=3
  !p.charthick=thick & !p.thick=thick & !x.thick=thick & !y.thick=thick
end
;
siz=1.7
!p.multi=0
!x.margin=[6.4,6.5]
!y.margin=[3.2,2.5]
!p.charsize=siz
tilde='!9!s!aA!n!r!6'
;
file='xyaver.sav'
dir2='yHf1em3_t2em5_k180b'
restore,file
eta0=2e-5 & k0=180.
eps0=eta0^3*k0^2
B20=eta0^2
;
nz=n_elements(z)
zeta=(1.-(z-!pi)/2.9)^4
col=grange(180,40,nz)
!x.title='!8B!6!drms!n/[!8B!6]'
!y.title='!7e!6/[!7e!6]'
;xr=minmax(sqrt(b2mz/B20))
;yr=minmax(epsMmz/eps0)
xr=[40,2e4]
yr=[1e1,2e9]
plot_oo,sqrt(b2mz[*,0]/B20),xr=xr,epsMmz[*,0]/eps0,yr=yr,/nodata,xst=9,yst=9
for i=10,nz-1,100 do begin
  ;oplot,sqrt(b2mz[i,*]/B20)/zeta[i]^.33333333,epsMmz[i,*]/eps0,col=col[i]
  oplot,sqrt(b2mz[i,*]/B20),epsMmz[i,*]/eps0,col=col[i]
endfor
;xx=[100.,1000.] & oplot,xx,xx^3/2.9e4
xx=[270.,2000.] & oplot,xx,xx^5/8e10
xyouts,3e3,6e6,'!6bottom',col=166,siz=siz*.8
xyouts,540,3e6,'!6top',col=55,siz=siz*.8
xyouts,380,3e1,'!6slope=5',siz=siz*.8
;xyouts,1e3,1e4,'!6slope=3'
;
xyouts,700.,1.7e8,'time',siz=siz*.8
arrow,/data,3000.,1e9,800.,3e7,thick=4
loadct,0
x0=2.28e2
y0=4.0e2
xx=[x0,18000.]& oplot,xx,xx*0+y0,col=200,li=3
yy=[y0,3e9]& oplot,yy*0+x0,yy,col=200,li=3
xyouts,80,8e2,'!6hel',siz=siz*.7,col=188
xyouts,7e3,2e8,'!6nonhel',siz=siz*.7,col=188
loadct,5
;
;  read second file XX
;
restore,'../'+dir2+'/'+file
for i=700,nz-1,100 do begin
  oplot,sqrt(b2mz[i,*]/B20),epsMmz[i,*]/eps0,col=col[i],l=2
print,'AXEL'
endfor
;
!x.title='    !8B!6!drms!n [G]'
!y.title='!8L!6 [erg s!u-1!n]'
epsunit=2.5d11
epsunit=2.5d9
Bfac=2.5d12 & EMfac=Bfac^2 ;(value from my paper)
Vol=1d12*100
axis,yax=1,yr=yr*epsunit*Vol*1d7
axis,xax=1,xr=xr*sqrt(EMfac)
print,'x0=',x0,x0*sqrt(EMfac)
;
default,texdir,'HallMHD'
print,'$mv idl.ps ~/tex/mhd/'+texdir+'/fig/ppxyaver_parametric.ps'
END