;$Id: pslopes.pro,v 1.15 2016/01/18 03:27:13 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=22,yoffset=3
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
;
;  mv idl.ps ~/tex/tina/Classes/fig/pslopes.ps
;
!p.charsize=2.7
!x.margin=[6.8,.5]
!y.margin=[3.2,.5]
!p.multi=[0,1,5]
iplot_aver=1 & dtmin=40. & dtmin4=10. & dtmin5=6.
iplot_aver=1 & dtmin=0. & dtmin4=10. & dtmin5=6.
xr=[0,64]
xr=[0,64]*4
siz=1.4
;
;dir1='Kol1152_k4b_120k0'
dir1='Kol1152_k4c_120k0'
dir1='Kol1152d'
dir1='Kol1152_k4a'
;dir1='Kol1152a_rep2'
dir1='Kol1152_k4c' ;(above the line)
dir1='Kol1152e' ;(returning to the line on another point)
dir1='Kol2304_k2a' ;(on the line, but further left at 1.2, 0.4)
dir1='Kol2304_k4a' ;(XX)
dir1='hydro1152k2_30sig0c_nu5em6' ;(k^2, kp=60, Jan 2016)
dir1b='Kol2304_k4b' ;(on the line, at 1.35, 0.35)
dir1b='hydro1152k4_sig1a' ;(on the line, at 1.35, 0.35, Dec 2015)
dir1b='hydro1152k4_120sig1c' ;(Dec 2015)
dir1b='hydro1152k4_120sig1c_nu1em5' ;(Dec 2015)
dir1b='hydro1152k2_120sig1c_nu1em5' ;(k^2, kp=120, Dec 2015)
dir1b='hydro1152k2_60sig1c_nu1em5' ;(k^2, kp=60, Dec 2015)
dir1b='hydro1152k2_30sig1c_nu5em6' ;(k^2, kp=60, Jan 2016)
dir2='MKol1152a5'
dir2='KH1152k2_sig0d'
dir3='../bforced/short2304pm1_kf60b_noforce'
dir4='H1152k4a'
dir5='H288k2G60_sig1c' ;(Tanmay run; test)
dir5='H1152k4a_120k0'
;
;  prepare all data
;
restore,'../'+dir1+'/pkf0.sav'
ievery=5
dt1=every(dt,ievery)
q1=deriv(alog(dt1),alog(1./every(kik,ievery)))
p1=deriv(alog(dt1),-alog(every(Ek,ievery)))
b1=p1/q1-1.
bb1=2./q1-3.
bbb1=4./(2.-p1)-3.
good1=where(dt1 gt dtmin)
q1m=mean(q1(good1))
p1m=mean(p1(good1))
b1m=mean(b1(good1))
bb1m=mean(bb1(good1))
bbb1m=mean(bbb1(good1))
;
;  prepare all data
;
restore,'../'+dir1b+'/pkf0.sav'
ievery=10
dt1b=every(dt,ievery)
q1b=deriv(alog(dt1b),alog(1./every(kik,ievery)))
p1b=deriv(alog(dt1b),-alog(every(Ek,ievery)))
b1b=p1b/q1b-1.
bb1b=2./q1b-3.
bbb1b=4./(2.-p1b)-3.
good1b=where(dt1b gt dtmin)
q1bm=mean(q1b(good1b))
p1bm=mean(p1b(good1b))
b1bm=mean(b1b(good1b))
bb1bm=mean(bb1b(good1b))
bbb1bm=mean(bbb1b(good1b))
;
restore,'../'+dir2+'/pkf.sav'
ievery=10
dt2=every(dt,ievery)
q2=deriv(alog(dt2),alog(1./every(kik,ievery)))
p2=deriv(alog(dt2),-alog(every(Ek,ievery)))
;q2=deriv(alog(dt),alog(1./kik))
;p2=deriv(alog(dt),-alog(Ek))
b2=p2/q2-1.
bb2=2./q2-3.
bbb2=4./(2.-p2)-3.
dt2=dt
good2=where(dt2 gt dtmin)
q2m=mean(q2(good2))
p2m=mean(p2(good2))
b2m=mean(b2(good2))
bb2m=mean(bb2(good2))
bbb2m=mean(bbb2(good2))
;
restore,'../'+dir3+'/pkf.sav'
ievery=2
dt3=every(dt,ievery)
q3=deriv(alog(dt3),alog(1./every(kik,ievery)))
p3=deriv(alog(dt3),-alog(every(Ek,ievery)))
;q3=deriv(alog(dt),alog(1./kik))
;p3=deriv(alog(dt),-alog(Ek))
b3=p3/q3-1.
bb3=2./q3-3.
bbb3=4./(2.-p3)-3.
dt3=dt
good3=where(dt3 gt dtmin)
q3m=mean(q3(good3))
p3m=mean(p3(good3))
b3m=mean(b3(good3))
bb3m=mean(bb3(good3))
bbb3m=mean(bbb3(good3))
;
restore,'../'+dir4+'/pkf.sav'
ievery=1
dt4=every(dt,ievery)
q4=deriv(alog(dt4),alog(1./every(kik,ievery)))
p4=deriv(alog(dt4),-alog(every(Ek,ievery)))
;q4=deriv(alog(dt),alog(1./kik))
;p4=deriv(alog(dt),-alog(Ek))
b4=p4/q4-1.
bb4=2./q4-3.
bbb4=4./(2.-p4)-3.
dt4=dt
good4=where(dt4 gt dtmin4)
q4m=mean(q4(good4))
p4m=mean(p4(good4))
b4m=mean(b4(good4))
bb4m=mean(bb4(good4))
bbb4m=mean(bbb4(good4))
;
restore,'../'+dir5+'/pkf.sav'
q5=deriv(alog(dt),alog(1./kik))
p5=deriv(alog(dt),-alog(Ek))
b5=p5/q5-1.
bb5=2./q5-3.
bbb5=4./(2.-p5)-3.
dt5=dt
good5=where(dt5 gt dtmin5)
q5m=mean(q5(good5))
p5m=mean(p5(good5))
b5m=mean(b5(good5))
bb5m=mean(bb5(good5))
bbb5m=mean(bbb5(good5))
;
;  panel 1
;
!x.title='!6'
!y.title='!8q!6!dK!n(!8t!6)'
plot,dt1,q1,xr=xr,yr=[.3,.8]
oplot,dt2,q2,col=122
oplot,dt3,q3,col=55
oplot,dt4,q4,col=155,li=2
oplot,dt5,q5,col=155
;
if iplot_aver eq 1 then begin
  oplot,dt1(good1),dt1(good1)*0+q1m,col=188
  oplot,dt2(good2),dt2(good2)*0+q2m,col=188
  oplot,dt3(good3),dt3(good3)*0+q3m,col=188
  oplot,dt4(good4),dt4(good4)*0+q4m,col=188
  oplot,dt5(good5),dt5(good5)*0+q5m,col=188
endif
;
;  panel 2
;
!y.title='!8p!6!dK!n(!8t!6)'
plot,dt1,p1,xr=xr,yr=[.5,1.8]
xyouts,30,1.52,'!6(i)',siz=siz
;
oplot,dt2,p2,col=122
xyouts,40,1.20,'!6(ii)',col=122,siz=siz
;
oplot,dt3,p3,col=55
xyouts,47,.80,'!6(iii)',col=55,siz=siz
;
oplot,dt4,p4,col=155,li=2
oplot,dt5,p5,col=155
xyouts,37,.64,'!6(iv)',col=155,siz=siz
;
if iplot_aver eq 1 then begin
  oplot,dt1(good1),dt1(good1)*0+p1m,col=188
  oplot,dt2(good2),dt2(good2)*0+p2m,col=188
  oplot,dt3(good3),dt3(good3)*0+p3m,col=188
  oplot,dt4(good4),dt4(good4)*0+p4m,col=188
  oplot,dt5(good5),dt5(good5)*0+p5m,col=188
endif
;
;  panel 3
;
!x.title='!8t!6'
!y.title='!7b!6!dK!n(!8t!6)'
plot,dt1,b1,xr=xr,yr=[-.5,4.5]
oplot,dt2,b2,col=122
oplot,dt3,b3,col=55
oplot,dt4,b4,col=155,li=2
oplot,dt5,b5,col=155
;
if iplot_aver eq 1 then begin
  oplot,dt1(good1),dt1(good1)*0+b1m,col=188
  oplot,dt2(good2),dt2(good2)*0+b2m,col=188
  oplot,dt3(good3),dt3(good3)*0+b3m,col=188
  oplot,dt4(good4),dt4(good4)*0+b4m,col=188
  oplot,dt5(good5),dt5(good5)*0+b5m,col=188
endif
;
;  panel 4
;
!x.title='!8t!6'
!y.title='!7b!6!dK!n(!8t!6)'
plot,dt1,bb1,xr=xr,yr=[-.5,4.5]
oplot,dt2,bb2,col=122
oplot,dt3,bb3,col=55
oplot,dt4,bb4,col=155,li=2
oplot,dt5,bb5,col=155
;
if iplot_aver eq 1 then begin
  oplot,dt1(good1),dt1(good1)*0+bb1m,col=188
  oplot,dt2(good2),dt2(good2)*0+b2m,col=188
  oplot,dt3(good3),dt3(good3)*0+bb3m,col=188
  oplot,dt4(good4),dt4(good4)*0+bb4m,col=188
  oplot,dt5(good5),dt5(good5)*0+bb5m,col=188
endif
;
;  panel 5
;
!x.title='!8t!6'
!y.title='!7b!6!dK!n(!8t!6)'
plot,dt1,bbb1,xr=xr,yr=[-.5,11.0]
oplot,dt2,bbb2,col=122
oplot,dt3,bbb3,col=55
oplot,dt4,bbb4,col=155,li=2
oplot,dt5,bbb5,col=155
;
if iplot_aver eq 1 then begin
  oplot,dt1(good1),dt1(good1)*0+bbb1m,col=188
  oplot,dt2(good2),dt2(good2)*0+bbb2m,col=188
  oplot,dt3(good3),dt3(good3)*0+bbb3m,col=188
  oplot,dt4(good4),dt4(good4)*0+bbb4m,col=188
  oplot,dt5(good5),dt5(good5)*0+bbb5m,col=188
endif
;
fo='(a,2f7.2,3f7.1)'
print,'q1m,p1m,b1m=',q1m,p1m,b1m,bb1m,bbb1m,fo=fo
print,'q2m,p2m,b2m=',q2m,p2m,b2m,bb2m,bbb2m,fo=fo
print,'q3m,p3m,b3m=',q3m,p3m,b3m,bb3m,bbb3m,fo=fo
print,'q4m,p4m,b4m=',q4m,p4m,b4m,bb4m,bbb4m,fo=fo
print,'q5m,p5m,b5m=',q5m,p5m,b5m,bb5m,bbb5m,fo=fo
;
END