;
; This routine is designed to calculate a histogram of the rotation
; measure through a plasma. I should have been integrated over angles,
; but since we have isotropy and homogenity it is equivalent to
; integrate over position and only consider the x-direction.
;
; nils: coded 25 Oct-04 (adapted from bdotl.pro)
;
;
; To make this work you need to read in the data and then calculate
; the magnetic field:
; .r start
; .r rall
; bb=curl(aa)
; rho=exp(lnrho)
;
PRO bdotl_hist,bb,rho,nx,time,dx=dx,direction1=direction1, $
          direction2=direction2,noplot=noplot, $
          integral=integral,hist=hist,nbins=nbins
;
; bb        : input: Magnetic field
; nx        : input: Number of mesh points in each direction (has to be cubic)
; time      : input: time of this snapshot (for use when writing results)
; dx        : input: mesh size (default assumes Lx=2pi)
; noplot    : input: Do we want to plot the results?
; nbins     : input: How many bins do we want in the histogram?
; direction1: input: Along what direction do we integrate (default is x) 
; direction2: input: Which component of B do we integrate (default is x) 
;                    must be the same as direction1 for line integral
;
; integral  : output: Array containing all the line integrals
; hist      : output: Array containing the histogram
;
; Set some default variables
;
default,dx,6.28/nx
default,nbins,100
default,direction1,'x'
default,direction2,'x'
;
; Set some auxillary variables
;
brms=rms(bb)
kf=1.5
if (direction1 eq 'x') then dir1=[1,0,0]
if (direction1 eq 'y') then dir1=[0,1,0]
if (direction1 eq 'z') then dir1=[0,0,1]
if (direction2 eq 'x') then dir2=[1,0,0]
if (direction2 eq 'y') then dir2=[0,1,0]
if (direction2 eq 'z') then dir2=[0,0,1]
;
; Initialize
;
integral=fltarr(nx,nx)
x1=3
;
; Loop over all pencils
;
for y1=3,nx+2 do begin
    for z1=3,nx+2 do begin
        if (direction1 eq 'x') then sp=[x1,y1,z1]
        if (direction1 eq 'y') then sp=[y1,x1,z1]
        if (direction1 eq 'z') then sp=[y1,z1,x1]
        ;
        ; Begin integration along path
        ;
        int=0
        for i=0,nx-1 do begin
            b1=bb(sp[0]+i*dir1(0),sp[1]+i*dir1(1),sp[2]+i*dir1(2),*)
            r1=rho(sp[0]+i*dir1(0),sp[1]+i*dir1(1),sp[2]+i*dir1(2))
            int=int+(dir2(0)*b1(0)+dir2(1)*b1(1)+dir2(2)*b1(2))*r1*dx
        end
        integral[y1-3,z1-3]=int
    end
end
;
; Find 2 < [RM(r=0)]^2 >
;
mean_2RM02=2*mean(integral^2.)
;
; How many bins do we have
;
hist=intarr(nbins)
Bdl=indgen(nbins)
;
; Find maximum and minimum values (needed when setting up the histogram)
;
maxvalue=max(integral)
minvalue=min(integral)
;
; What is the span in the histogram
;
span=maxvalue-minvalue
;
; Make array with Bdotl values for each bin
;
Bdl=Bdl*span/(nbins-1)+minvalue
print,'minvalue,maxvalue=',minvalue,maxvalue
;
; Loop over all pencils and put them into the histogram
;
for i=0,nx-1 do begin
    for j=0,nx-1 do begin
        index=fix((integral(i,j)-minvalue)/span*(nbins-1))
        hist(index)=hist(index)+1
    end
end
;
; Plot results
;
xvar=Bdl*kf/brms
yvar=hist
if (not keyword_set(noplot)) then begin
    plot,xvar,yvar
end
;
; Calculate structure functions
;
SF2=fltarr(nx/2)
for i=0,nx/2-1 do begin
    for j=0,nx-1 do begin
        for r=0,nx/2-1 do begin
            SF2(r)=SF2(r)+(integral(i+r,j)-integral(i,j))^2
        end
    end
end
for i=0,nx-1 do begin
    for j=0,nx/2-1 do begin
        for r=0,nx/2-1 do begin
            SF2(r)=SF2(r)+(integral(i,j+r)-integral(i,j))^2
        end
    end
end
;
; Save results
;
file='histogram_t='+string(time,FORMAT='(F7.2)') $
  +'_'+direction1+direction2+'.sav'
save,file=file,xvar,yvar,time,SF2,mean_2RM02
;
END