;-----------------------------------------------------------------------
function lok,step_arr,dir,usemaxstep,lcyl
;
; Do we want to consider this direction?
; We discard a direction if:
;    -stepmax is larger than usemaxstep
;    -we only average data over a cylinder, not a sphere
;     and step_arr(dir,2) ne 0
;
;
; What is the larges step?
;
stepmax=max([step_arr(dir,0),step_arr(dir,1),step_arr(dir,2)])
;
; Is the largest step larger than stepmax?
;
if (stepmax le usemaxstep) then begin
    ;
    ; Are we considering cylinedrical averages only?
    ;
    if (lcyl) then begin
        ;
        ; Is the z-step zero?
        ;
        if (step_arr(dir,2) eq 0) then begin
            ok=-1
        endif else begin
            ok=0
        endelse
    endif else begin
        ok=-1
    endelse
endif else begin
    ok=0
endelse
;
return,ok
;
END
;-----------------------------------------------------------------------
function ymaxval,local_xxx,nr_dir,step_arr,usemaxstep,lcyl,ll,step_length,$
                 derplot
;
ymaxval=-1e37
;
for dir=0,nr_dir-1 do begin
    ok=lok(step_arr,dir,usemaxstep,lcyl)
    ;
    if (ok) then begin
        toplot=local_xxx(*,dir)
        r=ll*step_length(dir)
        if (derplot) then begin
            toplot=deriv(alog(r),alog(local_xxx(*,dir)))
        endif
        if (max(toplot) gt ymaxval) then ymaxval=max(toplot)
    endif
endfor
;
return,ymaxval
;
end
;-----------------------------------------------------------------------
function yminval,local_xxx,nr_dir,step_arr,usemaxstep,lcyl,ll,step_length,$
                 derplot
;
yminval=1e37
;
for dir=0,nr_dir-1 do begin
    ok=lok(step_arr,dir,usemaxstep,lcyl)
    ;
    if (ok) then begin
        toplot=local_xxx(*,dir)
        r=ll*step_length(dir)
        if (derplot) then begin
            toplot=deriv(alog(r),alog(local_xxx(*,dir)))
        endif
        if (min(toplot) lt yminval) then yminval=min(toplot)
    endif
endfor
;
; Check that minimum is larger than (or equal to) zero.
;
if (yminval lt 0) then yminval=0
;
return,yminval
;
end
;-----------------------------------------------------------------------
function aver,ll,step_length,countnumber,toplot,average
;
interpolated=spline(ll*step_length,toplot,ll)
average=(countnumber-1)*average/(countnumber)+interpolated/(countnumber)
;
return,average
;
END
;-----------------------------------------------------------------------

PRO plot_struct_sph,nx,variable=variable,moment=moment,w=w,$
                    max_moment=max_moment,maxstep=maxstep,$
                    usemaxstep=usemaxstep,minstep=minstep, $
                    noprint=noprint,noplot=noplot, $
                    log_plot=log_plot,der_plot=der_plot, $
                    semilog_plot=semilog_plot,cyl=cyl,$
                    xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax, $
                    ll=ll,average=average,debug=debug, $
                    djuidjbi=djuidjbi,read_djuidjbi=read_djuidjbi
;
; nx         : number of grid points
; variable   : Which variable should be plotted?
; moment     : Which moment of the structure function do you want to
;            : show?
; maxstep    : Over how many 'shells' have you averaged?
; max_moment : How many moments have you calculated?
; read_djuidjbi : For Alex's simulations
; 
!p.charsize=1.8
;
; Define some logical variables
;
linplot=-1 
logplot=0 & derplot=0 & semilogplot=0 
lfirst=-1
lcyl=0
if (keyword_set(log_plot))     then logplot=-1
if (keyword_set(semilog_plot)) then semilogplot=-1
if (keyword_set(der_plot))     then derplot=-1
if (logplot or semilogplot)    then linplot=0
if (keyword_set(cyl))          then lcyl=-1
if (logplot and derplot) then begin
    print,'You are not allowed to have both logplot and derplot=T'
    stop
endif
;
; Set some default values
;
default,variable,'SF_flux'
default,moment,1
default,max_moment,6
default,maxstep,2
default,minstep,0
default,w,0.
default,usemaxstep,maxstep
;
; Print values of some variables
;
if (NOT keyword_set(noprint)) then begin
    print,''
    print,'Plotting variables:'
    print,'--------------------------------'
    print,'variable=',variable
    print,'moment=',moment
    print,'max_moment=',max_moment
    print,'maxstep=',maxstep
    print,'minstep=',minstep
    if (keyword_set(cyl)) then print,'Averageing over cylenders'
    print,'--------------------------------'
    print,''
end
;
; Defining arrays
;
lb_nxgrid=fix(alog(nx)/alog(2))
imax=lb_nxgrid*2-1
xxx=fltarr(imax,max_moment,(maxstep+1-minstep)^3)
step_length=fltarr((maxstep+1-minstep)^3)
step_arr=fltarr((maxstep+1-minstep)^3,3)
dist=fltarr(imax)
average=fltarr(imax)
average_lin=fltarr(imax)
;
if (keyword_set(debug)) then begin
    print,'nx=',nx
    print,'imax=',imax
    help,xxx
endif
;
; Find x values of data points (this must be the same
; algorithm as in struct_sph.f90)
;
counter=0
dist(counter)=0
dx=2*!pi/nx
for lb_ll=1,imax-1 do begin
    counter=counter+1
    exp2=lb_ll mod 2
    if (lb_ll eq 1) then exp2=0
    exp1=fix((lb_ll)/2)-exp2
    i=(2^exp1)*(3^exp2)
    dist(counter)=i*dx
end
;
;  Read data file
;
datatopdir='data'
close,1
file=datatopdir+'/'+variable+'.dat'
if (NOT keyword_set(noprint)) then print,'Reading ',file
openr,1,file
readf,1,t
readf,1,nr_dir
readf,1,step_length
readf,1,step_arr
readf,1,xxx
if (keyword_set(read_djuidjbi)) then  readf,1,djuidjbi
close,1
;
; Print time of this snapshot
;
if (NOT keyword_set(noprint)) then print,'time=',t
;
; Some abreviations
;
local_xxx=reform(xxx(1:*,moment-1,0:nr_dir-1))
ll=dist(1:*)
;
; Set max and min values
;
default,xmin,min(ll) 
default,xmax,max(ll)
default,ymax,ymaxval(local_xxx,nr_dir,step_arr, $
                     usemaxstep,lcyl,ll,step_length,derplot)
default,ymin,yminval(local_xxx,nr_dir,step_arr, $
                     usemaxstep,lcyl,ll,step_length,derplot)
!x.range=[xmin,xmax]
!y.range=[ymin,ymax]
;
; Initializing
;
average=0
average_lin=0
countnumber=0
;
; Set x and y labels
;
mm=string(moment,FORMAT='(I1)')
!x.title='!8r!6'
if (linplot) then !y.title='!8SF!6!d'+mm+'!n'
if (logplot) then !y.title='!8SF!6!d'+mm+'!n'
if (derplot) then !y.title='dlog(!8SF!6!d'+mm+'!n)/dlog(!8r!6)'
;
; Loop over all directions
;
for dir=0,nr_dir-1 do begin
    ;
    ; Do we want to consider this direction?
    ; We discard a direction if:
    ;    -stepmax is larger than usemaxstep
    ;    -we only average data over a cylinder, not a sphere
    ;     and step_arr(dir,2) ne 0
    ;
    ok=lok(step_arr,dir,usemaxstep,lcyl)
    ;
    if (ok) then begin
        countnumber=countnumber+1
        toplot=local_xxx(*,dir)
        r=ll*step_length(dir)
        if (derplot) then toplot=deriv(alog(r),alog(local_xxx(*,dir)))
        ;
        ; Do we want to plot the data
        ;
        if (NOT keyword_set(noplot)) then begin
            ;
            ; Is it the first time we want to plot?
            ;
            if (lfirst) then begin
                if (linplot) then plot,r,toplot
                if (logplot) then plot_oo,r,toplot
                if (semilogplot) then plot_oi,r,toplot
                oplot,r,toplot,psym=1,col=122
                lfirst=0
            endif else begin
                wait,w
                oplot,r,toplot
                oplot,r,toplot,psym=1,col=122
            endelse
            ;
            ; Add colors to SF's in the coordinate directions
            ;
            if (fix(total(step_arr(dir,*)) eq 1)) then begin                
                if (fix(step_arr(dir,0)) eq 1) then oplot,r,toplot,li=2,col=122
                if (fix(step_arr(dir,1)) eq 1) then oplot,r,toplot,li=2,col=47
                if (fix(step_arr(dir,2)) eq 1) then oplot,r,toplot,li=2,col=150
            endif                
        endif
        ;
        ; Find average
        ;
        average=aver(ll,step_length(dir),countnumber,toplot,average)
        average_lin=aver(ll,step_length(dir),countnumber,local_xxx(*,dir),average_lin)
    endif
endfor
;
; Plot average
;
oplot,ll,average,col=47,thick=3
if (derplot) then oplot,ll,deriv(alog(r),alog(average_lin)),col=122,thick=3
;
print,'The number of directions is ',countnumber
;
!x.title=''
!y.title=''
!x.range=''
!y.range=''
;
END