;$Id: pdel_zt.pro,v 1.14 2009/11/21 13:09:40 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,delmax,.05
default,etamax,.01
default,epsmax,.002
default,iwrite,0
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)
default,kf,5.
;
;  give modified alpmax values in parameters.pro file
;
yrdel=[-.5,1.]*delmax
!p.charsize=2
!p.multi=[0,1,2]
!p.multi=[0,1,3]
plot,tt,eta12,yr=yrdel,ytit='!7-d!6'
if itestfield eq 'B11-B22' then begin
  oplot,tt,+eta21,li=2
  del=.5*(eta12-eta21)
endif else begin
  del=.5*eta12
endelse
oplot,tt,urms/3.,li=1
pc_error_range,tt(good),del(good),mean=delm,error=del_error
pc_error_range,tt(good),+eta12(good),/oplot,mean=eta12m,error=eta12_error
pc_error_range,tt(good),+eta21(good),/oplot,mean=eta21m,error=eta21_error
oplot,tt(good),+accum(eta12(good)),col=55
oplot,tt(good),+accum(eta21(good)),col=122
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
;
;  diagonals the same?
;
if itestfield eq 'B11-B22' then begin
  yreps=[-1.,1.]*epsmax
  eps=.5*(eta11-eta22)
  plot,tt,eps,yr=yreps
  pc_error_range,tt(good),eps(good),/oplot,mean=epsm,error=eps_error
  oplot,tt(good),accum(eps(good)),col=55
  oplot,tt,tt-tt,li=3
endif
;
;  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_t=',delm,' +/-',del_error,fo=fo2
print,'+eta12=',eta12m,' +/-',eta12_error,fo=fo2
print,'+eta21=',eta21m,' +/-',eta21_error,fo=fo2
print,'  epsm=',epsm,' +/-',eps_error,fo=fo2
print,'  urms=',urmsm,' +/-',urmsm_error,fo=fo
print
!p.multi=0
;
@alphaeta_rms
fo='(f6.3,2e8.1,i6,12e10.2,"  ",a)'
fo2='(f6.3,2e8.1,i6,6e10.2,"  ",a)'
pc_read_param,/param2,obj=param
if iuu eq 0 then nu=0. else nu=param.nu
Omega=param.Omega
Shear=param.Sshear
;print,Shear,nu,param.etatest,fix(max(tt)),etam,eta_error,eta12m,eta12_error,eta21m,eta21_error,urmsm,urmsm_error,epsm,eps_error,alprms_all,eta21rms_all,run,fo=fo
;
;  write the same to the disk
;
if n_elements(alprms_all) eq 0 then begin
  print,'.r pc_alpha_flucts'
endif
;
;  write the same to the disk
;
fo='(f6.3,2e8.1,i6,10e10.2,8e8.1,"  ",a)'
openw,1,'summary.dat'
printf,1,Shear,nu,param.etatest,fix(max(tt)),etam,eta_error,$
  eta12m,eta12_error,eta21m,eta21_error,urmsm,urmsm_error, epsm,eps_error,$
  alprms_all,alprms_err_all,etatrms_all,etatrms_err_all,$
  eta21rms_all,eta21rms_err_all,eta12rms_all,eta12rms_err_all,run,fo=fo
close,1
spawn,'cat summary.dat'
if iwrite eq 1 then begin
  ;spawn,'cat summary.dat >> $HOME/tex/mhd/alphaeta/idl/kinshear.dat'
  spawn,'cat summary.dat >> $HOME/tex/mhd/alphaeta/idl/kinshear_Pm1.dat'
endif
;
print
print,Omega,nu,param.etatest,fix(max(tt)),etam,eta_error,delm,del_error,urmsm,urmsm_error,run,fo=fo2
print
;
;  save for other post-processing
;
save,file='eta21etc.sav',tt,eta,eta21,eta12,alp22,shear,urmsm,kf
;
END