; This is a routine to find the longitudinal structure functions 
; calculated in different directions.
;
; .r start
; .r rall
;
; Some default variables
;
maxmom=3
maxstep=1
minstep=0
;
; Find total number of directions (should remove repeating directions
; later e.g (1,1,0) and (2,2,0)) 
;
ndir=(maxstep+1)^3.
;
; Define arrays
;
SF=fltarr(nx/2,maxmom+1,ndir)
separation=fltarr(nx/2,ndir)
rvec=fltarr(3,ndir)
;
; Some auxillary variables
;
uuu=uu(l1:l2,m1:m2,n1:n2,*)
uu1=uuu(*,*,*,0)
uu2=uu1
;
; Initialize
;
SF(*,*,*)=0
idir=0
;
; Loop over all directions
;
for xstep=minstep,maxstep do begin
for ystep=minstep,maxstep do begin
for zstep=minstep,maxstep do begin
    if (xstep+ystep+zstep ne 0) then begin
        print,xstep,ystep,zstep
        ;
        ; What is the length of each step
        ;
        rlength=sqrt(xstep^2+ystep^2+zstep^2)
        ;
        ; Find longitudinal component
        ;
        uu1(*,*,*)= $  ;u(x)
          +uuu(*,*,*,0)*xstep/rlength $
          +uuu(*,*,*,1)*ystep/rlength $
          +uuu(*,*,*,2)*zstep/rlength
        ;
        ; Loop over all separations
        ;
        for rr=1,nx/2 do begin
            rvec(*,idir)=[xstep,ystep,zstep]
            separation(rr-1,idir)=rr*rlength*dx
            uu2=shift(uu1,xstep*rr,ystep*rr,zstep*rr) ;u(x+r)
            dval=abs(uu2-uu1) ; |u(x+r)-u(x)|
            ;
            ; Find different moments
            ;
            for mom=1,maxmom do begin
                SF(rr-1,mom,idir)=mean(dval^mom) ;<|u(x+r)-u(x)|^p>
                ;SF(rr-1,mom,idir)=mean(dval(0:nx/2-1,0:ny/2-1,0:nz/2-1)^mom)
            endfor
        endfor
    endif
    idir=idir+1
endfor
endfor
endfor
;
; Save results
;
file='sf_test.sav'
save,file=file,SF,separation,maxmom,rvec,ndir
;
; Plot results
;
print,'--------------------Plotting results-----------------'
;
moment=1 & ymin=0 & ymax=0.1
i=1
w=0.
plot,separation(*,i),SF(*,moment,i),yr=[ymin,ymax],xr=[0,3.14]
oplot,separation(*,i),SF(*,moment,i),li=2,col=150; rvec=[0,0,1]
print,rvec(*,i)
for i=2,ndir-1 do begin
    wait,w
    print,rvec(*,i)
    oplot,separation(*,i),SF(*,moment,i)
    if ((rvec(0,i) eq 1)and(rvec(1,i) eq 0)and(rvec(2,i) eq 0)) then begin
        oplot,separation(*,i),SF(*,moment,i),li=2,col=122; rvec=[1,0,0]
    endif
    if ((rvec(0,i) eq 0)and(rvec(1,i) eq 1)and(rvec(2,i) eq 0)) then begin
        oplot,separation(*,i),SF(*,moment,i),li=2,col=47; rvec=[0,1,0]
    endif
endfor
;
END