;$Id: pdelall_zt.pro,v 1.1 2007/09/22 17:56:30 brandenb Exp $
;
;  plot alpa and eta results from time series
;
pc_read_ts,o=o
;
;  set minimum time after which averaging begins
;
t2=1e30
spawn,'touch parameters.pro'
@parameters
@data/index.pro
@data/testfield_info.dat
default,mumax,.05
default,delmax,.05
default,etamax,.01
default,run,''
;
;  introduce abbreviations
;
tt=o.t
urms=o.urms
;
;  del tensor
;
alp11=o.alp11
alp21=o.alp21
;
;  eta tensor as in EMF_i = ... -eta_ij J_j
;  so these new etas are also referred to as eta^*
;  writing this furthermore as EMF = delta x J - etat*J,
;  with delta=(0,0,delta), so we have delta x J = (-Jy, Jx, 0)*delta,
;  eta*_12 = -eta*_21 = delta
;
eta12=-o.eta11
eta22=-o.eta21
;
;  read extra fields
;
if itestfield eq 'B11-B22' then begin
  alp12=o.alp12
  alp22=o.alp22
  eta11=+o.eta12
  eta21=+o.eta22
endif
;
;  range of time where to do the analysis
;
tmin=min(tt)
tmax=max(tt)
default,t1,(tmin+tmax)/2.
good=where(tt gt t1 and tt lt t2)
kf=5.
;
;  give modified alpmax values in parameters.pro file
;
!p.charsize=2
!p.multi=[0,1,3]
;
;  delta
;
yrdel=[-.5,1.]*delmax
del=.5*(eta12-eta21)
plot,tt,del,yr=yrdel,ytit='!7-d!6'
pc_error_range,tt(good),del(good),mean=delm,error=del_error
oplot,tt(good),accum(del(good)),col=55
oplot,tt,tt-tt,li=3
;
;  mu
;
yrmu=[-.5,1.]*mumax
mu=.5*(eta12+eta21)
plot,tt,mu,yr=yrmu,ytit='!7-d!6'
pc_error_range,tt(good),mu(good),mean=mum,error=mu_error
oplot,tt(good),accum(mu(good)),col=55
oplot,tt,tt-tt,li=3
;
;  give modified etamax values in parameters.pro file
;
yreta=[-.2,1.]*etamax
plot,tt,eta22,yr=yreta
if itestfield eq 'B11-B22' then begin
  oplot,tt,eta11,li=2
  eta=.5*(eta11+eta22)
endif else begin
  eta=eta22
endelse
oplot,tt,urms/(3.*kf),li=1
pc_error_range,tt(good),eta(good),/oplot,mean=etam,error=eta_error
oplot,tt(good),accum(eta(good)),col=55
oplot,tt,tt-tt,li=3
;oplot,tt,etat,col=122,li=0
;
pc_error_range,tt(good),urms(good),/oplot,mean=urmsm,error=urmsm_error
;
;  print on screen
;
print
fo='(2(a,f6.4))'
fo2='(2(a,f8.5))'
fo2='(2(a,e9.2))'
print,' eta_t=',etam,' +/-',eta_error,fo=fo2
print,'   del=',delm,' +/-',del_error,fo=fo2
print,'    mu=',mum,' +/-',mu_error,fo=fo2
print,'  urms=',urmsm,' +/-',urmsm_error,fo=fo
print
!p.multi=0
;
fo='(2e8.1,i6,8e10.2,"  ",a)'
fo2='(2e8.1,i6,6e10.2,"  ",a)'
pc_read_param,/param2,obj=param
if iuu eq 0 then nu=0. else nu=param.nu
print,nu,param.etatest,fix(max(tt)),etam,eta_error,eta12m,eta12_error,eta21m,eta21_error,urmsm,urmsm_error,run,fo=fo
print
print,nu,param.etatest,fix(max(tt)),etam,eta_error,delm,del_error,urmsm,urmsm_error,run,fo=fo2
print
END