;$Id: p.pro,v 1.5 2008/07/26 18:20:48 brandenb Exp $
;
;  read data, plot current helicity fluxes, and write averages
;
default,dirname,"~/tex/notes/VCflux/idl/"
default,filename,'1em2.dat'
default,kf,1.5
default,t1,0.
default,t2,9e9
;
;  read run parameters
;
pc_read_param,obj=param,/param2
eta=param.eta
nu=param.nu
B_ext=param.B_ext
Shear=abs(param.Sshear)
B0x=B_ext(0)
B0y=B_ext(1)
B0=sqrt(B0x^2+B0y^2)
;
spawn,'touch parameters.pro'
@parameters
;
;  read time series
;
pc_read_ts,obj=ts
tt=ts.t
urms=ts.urms
brms=ts.brms
nt=n_elements(tt)
;
good=where(tt ge t1 and tt le t2)
urms0=mean(urms(good))
it1=min(good)
it2=nt-1
;
;  calculate the 2 contributions to the flux
;  and present results in normalized form
;
if B0 ne 0. then norm1=1./(B0^2*Shear) else norm1=1.
flux1=2.*ts.exjmz*norm1
flux2=ts.dexbmz*norm1
fluxt=flux1+flux2
;
plot,tt,fluxt,yr=[-.2,1]*max(fluxt)
pc_error_range,tt(good),fluxt(good),mean=fluxtm,error=fluxte,/oplot
oplot,tt,accum(fluxt),col=122
oplot,tt,tt*0
;
oplot,tt,flux2,col=55
pc_error_range,tt(good),flux2(good),mean=flux2m,error=flux2e,/oplot
oplot,tt,accum(flux2),col=55
;
;ratio=ts.dexbmz/ts.exjmz & plot,tt,ratio
;
; open output file
;
fo='(2e9.1,5f7.3,4f8.4,3x,a)'
openw,1,'../'+filename
print,nu,eta,urms0,fluxtm,fluxte,flux2m,flux2e
printf,1,nu,eta,Shear,B0,B0x,B0y,urms0,fluxtm,fluxte,flux2m,flux2e,Run,fo=fo
close,1
;
spawn,'cat ../'+filename+' >> '+dirname+filename
;
print,'gt ',dirname
print,'vi ',filename
;
Rm=urms0/(eta*kf)
Sh=Shear/(urms0*kf)
BBeq=B0/urms0
Cg=BBeq^2
print,'Rm,Sh,Cg=',Rm,Sh,BBeq,Cg
;
END