;$Id: pe_comp.pro,v 1.12 2021/02/05 17:50:30 brandenb Exp $
if !d.name eq 'PS' then begin
  device,xsize=18,ysize=5,yoffset=3
  !p.charthick=1.6 & !p.thick=1.6 & !x.thick=1.6 & !y.thick=1.6
end
;
siz=.9
!p.charsize=.8
!x.margin=[8.3,.5]
!y.margin=[3.2,.5]
loplot=1
;
dir1='F1152a_sig0_t11_M4_ramp01b' & li1=1 & col1=122
dir2='F1152a_sig0_t11_M4_ramp02a' & li2=2 & col2=122
dir3='F1152a_sig0_t11_M4_ramp05a' & li3=3 & col3=55
dir4='F1152a_sig0_t11_M4_ramp1a' & li4=4 & col4=122
dir5='F1152a_sig0_t11_M4_ramp2a' & li5=0
;
default,iread,0
if iread eq 0 then begin
  pc_read_ts,o=ts1,datadir='../'+dir1+'/data'
  pc_read_ts,o=ts2,datadir='../'+dir2+'/data'
  pc_read_ts,o=ts3,datadir='../'+dir3+'/data'
  pc_read_ts,o=ts4,datadir='../'+dir4+'/data'
  pc_read_ts,o=ts5,datadir='../'+dir5+'/data'
  iread=1
endif
;
default,kin,0
default, r_pe,[1.0,1.64]
default,taver1,1.15
default,kf,600.
default,xr_pe,1.
circ_sym,1.3,1
fo="(e8.2,e10.2,i3,f6.2,i5,2x,a)"
!p.multi=[0,2,1]
;
t1a=1.14 & t1b=1.27 & t1c=1.2
t2a=1.14 & t2b=1.32 & t2c=1.3
t3a=1.14 & t3b=1.4 & t3c=1.5
t4a=1.14 & t4b=1.5 & t4c=1.7
t5a=1.14 & t5b=1.7 & t5c=2.2
;
;-----------------------------------------------------------------------------
;  panel 1
;
!x.title='!8t!6'
!y.title='!7X!6!dM!n(!8t!6)'
!y.title='!13E!6!dM!n(!8t!6)'
xr_pe=[1.,3.24]
yr=[0.,.043]
yr=[1e-3,.1]
th=6
;
run=dir1
if kin then EEKM=ts1.EEK else EEKM=ts1.EEM
;plot,ts1.t,EEKM,yr=yr,xr=xr_pe,li=li1
plot_io,ts1.t,EEKM,yr=yr,xr=xr_pe,li=li1
imax=(where(EEKM eq max(EEKM)))[0]
n2=n_elements(EEKM)-1
EEKMmax=mean(EEKM[(imax-1)>0:(imax+1)<n2])
tmax=ts1.t[imax]
;oplot,[1,1]*tmax,[1,1]*EEKMmax,col=122,ps=8
good=where(ts1.t ge t1a and ts1.t le t1b) & tgood=ts1.t(good)
p1=linfit(tgood,alog(EEKM(good)))
if loplot then oplot,tgood,exp(tgood*p1[1]+p1[0]),l=1,th=th
;!p.multi=0
;stop
;
loadct,6
run=dir2
if kin then EEKM=ts2.EEK else EEKM=ts2.EEM
oplot,ts2.t,EEKM,li=li2,col=col2
imax=(where(EEKM eq max(EEKM)))[0]
n2=n_elements(EEKM)-1
EEKMmax=mean(EEKM[(imax-1)>0:(imax+1)<n2])
tmax=ts2.t[imax]
;oplot,[1,1]*tmax,[1,1]*EEKMmax,col=55,ps=8
good=where(ts2.t ge t2a and ts2.t le t2b) & tgood=ts2.t(good)
p2=linfit(tgood,alog(EEKM(good)))
if loplot then oplot,tgood,exp(tgood*p2[1]+p2[0]),col=122,th=th
loadct,5
;
run=dir3
if kin then EEKM=ts3.EEK else EEKM=ts3.EEM
oplot,ts3.t,EEKM,li=li3,col=col3
imax=(where(EEKM eq max(EEKM)))[0]
n2=n_elements(EEKM)-1
EEKMmax=mean(EEKM[(imax-1)>0:(imax+1)<n2])
tmax=ts3.t[imax]
;oplot,[1,1]*tmax,[1,1]*EEKMmax,col=155,ps=8
good=where(ts3.t ge t3a and ts3.t le t3b) & tgood=ts3.t(good)
p3=linfit(tgood,alog(EEKM(good)))
if loplot then oplot,tgood,exp(tgood*p3[1]+p3[0]),li=li3,col=col3,th=th
;
run=dir4
if kin then EEKM=ts4.EEK else EEKM=ts4.EEM
oplot,ts4.t,EEKM,li=li4,col=col4
imax=(where(EEKM eq max(EEKM)))[0]
n2=n_elements(EEKM)-1
EEKMmax=mean(EEKM[(imax-1)>0:(imax+1)<n2])
tmax=ts4.t[imax]
;oplot,[1,1]*tmax,[1,1]*EEKMmax,col=155,ps=8
good=where(ts4.t ge t4a and ts4.t le t4b) & tgood=ts4.t(good)
p4=linfit(tgood,alog(EEKM(good)))
if loplot then oplot,tgood,exp(tgood*p4[1]+p4[0]),li=li4,col=col4,th=th
;
run=dir5
if kin then EEKM=ts5.EEK else EEKM=ts5.EEM
oplot,ts5.t,EEKM,li=li5
imax=(where(EEKM eq max(EEKM)))[0]
n2=n_elements(EEKM)-1
EEKMmax=mean(EEKM[(imax-1)>0:(imax+1)<n2])
tmax=ts5.t[imax]
;oplot,[1,1]*tmax,[1,1]*EEKMmax,col=155,ps=8
good=where(ts5.t ge t5a and ts5.t le t5b) & tgood=ts5.t(good)
p5=linfit(tgood,alog(EEKM(good)))
if loplot then oplot,tgood,exp(tgood*p5[1]+p5[0]),li=li5,col=col5,th=th
;
xyouts,1.05,.065,siz=siz,'!6(a)'
xyouts,2.55,.011,siz=siz,'!7s!6=2'
xyouts,1.10,.002,siz=siz,'!7s!6=0.1' & loadct,6
xyouts,1.39,.004,siz=siz*.8,'0.2',col=122 & loadct,5
xyouts,1.59,.005,siz=siz,'0.5',col=55
xyouts,1.92,.007,siz=siz,'!7s!6=1',col=122
;
;-----------------------------------------------------------------------------
;  panel 2
;
!y.title='!13E!6!dGW!n(!8t!6)/10!u-6!n' & s=1e6
;
yr=[0,2.6e-1]
plot,ts1.t,s*ts1.EEGW,xr=xr_pe,yr=yr,li=li1
good=where(ts1.t ge taver1)
EEGWm=mean(ts1.EEGW(good))
;oplot,ts1.t(good),ts1.t(good)*0.+EEGWm,col=122
print,EEKMmax,EEGWm,kin,tmax,kf,run,fo=fo
good=where(ts1.t ge t1c) & tgood=ts1.t(good)
pp1=linfit(tgood,ts1.EEGW(good)) & mm1=mean(ts1.EEGW(good))
;oplot,tgood,s*(tgood*pp1[1]+pp1[0]),li=li1,th=th
;
loadct,6
oplot,ts2.t,s*ts2.EEGW,li=li2,col=col2
good=where(ts2.t ge taver1)
EEGWm=mean(ts2.EEGW(good))
;oplot,ts2.t(good),ts2.t(good)*0.+EEGWm,col=122
print,EEKMmax,EEGWm,kin,tmax,kf,run,fo=fo
good=where(ts2.t ge t2c) & tgood=ts2.t(good)
pp2=linfit(tgood,ts2.EEGW(good)) & mm2=mean(ts2.EEGW(good))
;oplot,tgood,s*(tgood*pp2[1]+pp2[0]),li=li2,col=col2,th=th
loadct,5
;
oplot,ts3.t,s*ts3.EEGW,li=li3,col=col3
good=where(ts3.t ge taver1)
EEGWm=mean(ts3.EEGW(good))
;oplot,ts3.t(good),ts3.t(good)*0.+EEGWm,col=122
print,EEKMmax,EEGWm,kin,tmax,kf,run,fo=fo
good=where(ts3.t ge t3c) & tgood=ts3.t(good)
pp3=linfit(tgood,ts3.EEGW(good)) & mm3=mean(ts3.EEGW(good))
;oplot,tgood,s*(tgood*pp3[1]+pp3[0]),li=li3,col=col3,th=th
;
oplot,ts4.t,s*ts4.EEGW,li=li4,col=col4
good=where(ts4.t ge taver1)
EEGWm=mean(ts4.EEGW(good))
;oplot,ts4.t(good),ts4.t(good)*0.+EEGWm,col=122
print,EEKMmax,EEGWm,kin,tmax,kf,run,fo=fo
good=where(ts4.t ge t4c) & tgood=ts4.t(good)
pp4=linfit(tgood,ts4.EEGW(good)) & mm4=mean(ts4.EEGW(good))
;oplot,tgood,s*(tgood*pp4[1]+pp4[0]),li=li4,col=col4,th=th
;
oplot,ts5.t,s*ts5.EEGW,li=li5
good=where(ts5.t ge taver1)
EEGWm=mean(ts5.EEGW(good))
;oplot,ts5.t(good),ts5.t(good)*0.+EEGWm,col=122
print,EEKMmax,EEGWm,kin,tmax,kf,run,fo=fo
good=where(ts5.t ge t5c) & tgood=ts5.t(good)
pp5=linfit(tgood,ts5.EEGW(good)) & mm5=mean(ts5.EEGW(good))
;oplot,tgood,s*(tgood*pp5[1]+pp5[0]),li=li5,col=col5,th=th
;
xyouts,1.07,.235,siz=siz,'!6(b)'
;
!p.multi=0
print,'$convert idl.ps /D/Print/Tuusula_GW/pe_comp.png'
print,'$mv idl.ps ~/tex/tina/duration/fig/pe_comp.ps'
print,'$mv idl.ps ~/tex/tina/CircPol/fig/pe_comp.ps'
print,'$mv idl.ps ~/tex/tina/Circular-Polarization-paper-PRL/trunk/fig/pe_comp.ps'
;
print,p1
print,p2
print,p3
print,p4
print,p5
print,mm1,pp1
print,mm2,pp2
print,mm3,pp3
print,mm4,pp4
print,mm5,pp5
END