;
; This program produce data for the scatterplot of: 
; u^2 vs b^2
; rho vs b^2
; rho u^2 vs b^2
;
PRO scatterplot,uu,bb,rho_inn,time,nbins=nbins,nx=nx
;
; Set some defaults
;
default,nbins,100
default,nx,256
;
; Set up some arrays to store the scatter plots in
; (There is too much data to store all individual scatters)
;
v2_b2=fltarr(nbins,nbins)
rho_b2=fltarr(nbins,nbins)
rhov2_b2=fltarr(nbins,nbins)
;
; Define the arrays of which we want scatter plots (e.g. v2 and b2)
;
v2=dot2(uu(3:nx+2,3:nx+2,3:nx+2,*))
b2=dot2(bb(3:nx+2,3:nx+2,3:nx+2,*))
rho=rho_inn(3:nx+2,3:nx+2,3:nx+2)
rhov2=rho*v2
;
; What are the max values?
;
maxv2=max(v2)
maxb2=max(b2)
maxrhov2=max(rhov2)
maxrho=max(rho)
minrho=min(rho)
;
; Find where the different arrays fit into the scatter plots
;
v2index    = floor((v2/maxv2)*(nbins-1))
rhov2index = floor((rhov2/maxrhov2)*(nbins-1))
rhoindex   = floor(((rho-minrho)/(maxrho-minrho))*(nbins-1))
b2index    = floor((b2/maxb2)*(nbins-1))
;
; Set some variables needed for the correlation coefficient
;
meanrhov2=mean(rhov2)
meanb2=mean(b2)
meanrho=mean(rho)
print,'meanrho,meanb2,meanrhov2=',meanrho,meanb2,meanrhov2
upper_rhov2=0
lower_rhov2=0
upper_rho=0
lower_rho=0
lower_b2=0
;
; Loop over all points
;
for x1=0,nx-1 do begin
print,'x1=',x1
    for y1=0,nx-1 do begin
        for z1=0,nx-1 do begin
            v2_b2(v2index(x1,y1,z1),b2index(x1,y1,z1))$
              = v2_b2(v2index(x1,y1,z1),b2index(x1,y1,z1))       + 1
            rho_b2(rhoindex(x1,y1,z1),b2index(x1,y1,z1))$
              = rho_b2(rhoindex(x1,y1,z1),b2index(x1,y1,z1))     + 1
            rhov2_b2(rhov2index(x1,y1,z1),b2index(x1,y1,z1))$
              = rhov2_b2(rhov2index(x1,y1,z1),b2index(x1,y1,z1)) + 1
            ;
            upper_rhov2=upper_rhov2 $
              +(rhov2(x1,y1,z1)-meanrhov2)*(b2(x1,y1,z1)-meanb2)
            lower_rhov2=lower_rhov2+(rhov2(x1,y1,z1)-meanrhov2)^2
            ;
            upper_rho=upper_rho $
              +(rho(x1,y1,z1)-meanrho)*(b2(x1,y1,z1)-meanb2)
            lower_rho=lower_rho+(rho(x1,y1,z1)-meanrho)^2
            ;
            lower_b2=lower_b2+(b2(x1,y1,z1)-meanb2)^2     
        end
    end
end

CC_rhov2_b2=upper_rhov2/sqrt(lower_rhov2*lower_b2)
print,'CC_rhov2_b2=',CC_rhov2_b2
CC_rho_b2=upper_rho/sqrt(lower_rho*lower_b2)
print,'CC_rho_b2=',CC_rho_b2

;
; Save data
;
file='scatterplot_t='+string(time,FORMAT='(F7.2)')+'.sav'
save,file=file,v2_b2,rho_b2,rhov2_b2,time,nbins,maxv2,maxb2,maxrhov2,maxrho,minrho,CC_rhov2_b2,CC_rho_b2
;
END