function make_link_list,xx,cellx,celly,cellz,mcellx,mcelly,mcellz,xini,yini,zini
  si=size(xx)
  help,xx
  nn=si[1]
  head=0*indgen(mcellx+2,mcelly+2,mcellz+2,/LONG) 
  list=0*indgen(nn,/LONG)
  for ip=0,nn-1 do begin
    icell=floor((xx[ip,0]-xini)/cellx)+1
    jcell=floor((xx[ip,1]-yini)/celly)+1
    kcell=floor((xx[ip,2]-zini)/cellz)+1
    list[ip] = head[icell,jcell,kcell]
    head[icell,jcell,kcell]=ip
  end
; apply PBC to head ...
; first the faces 
  head[0,1:mcelly,1:mcellz]=head[mcellx,1:mcelly,1:mcellz]
  head[mcellx+1,1:mcelly,1:mcellz]=head[1,1:mcelly,1:mcellz]
;
  head[1:mcellx,0,1:mcellz]=head[1:mcellx,mcelly,1:mcellz]
  head[1:mcellx,mcelly+1,1:mcellz]=head[1:mcellx,1,1:mcellz]
;
  head[1:mcellx,1:mcelly,0]=head[1:mcellx,1:mcelly,mcellz]
  head[1:mcellx,1:mcelly,mcellz+1]=head[1:mcellx,1:mcelly,1]

; now the edges 
  head[0,0,1:mcellz]      =head[mcellx,mcelly,1:mcellz]
  head[0,mcelly+1,1:mcellz]=head[mcellx,1,1:mcellz]

  head[mcellx+1,0,1:mcellz]      =head[1,mcelly,1:mcellz]
  head[mcellx+1,mcelly+1,1:mcellz]=head[1,1,1:mcellz]

  head[0,1:mcelly,0]      =head[mcellx,1:mcelly,mcellz]
  head[0,1:mcelly,mcellz+1]=head[mcellx,1:mcelly,1]

  head[mcellx+1,1:mcelly,0]      =head[1,1:mcelly,mcellz]
  head[mcellx+1,1:mcelly,mcellz+1]=head[1,1:mcelly,1]

  head[1:mcellx,0,0]      =head[1:mcellx,mcelly,mcellz]
  head[1:mcellx,0,mcellz+1]=head[1:mcellx,mcelly,1]

  head[1:mcellx,mcelly+1,0]      =head[1:mcellx,1,mcellz]
  head[1:mcellx,mcelly+1,mcellz+1]=head[1:mcellx,1,1]

; and the corners
  head[0,0,0]              =head[mcellx,mcelly,mcellz]
  head[0,0,mcellz+1]       =head[mcellx,mcelly,1]
  head[0,mcelly+1,0]       =head[mcellx,1,mcellz]
  head[0,mcelly+1,mcellz+1]=head[mcellx,1,1]

  head[mcellx+1,0,0]              =head[1,mcelly,mcellz]
  head[mcellx+1,0,mcellz+1]       =head[1,mcelly,1]
  head[mcellx+1,mcelly+1,0]       =head[1,1,mcellz]
  head[mcellx+1,mcelly+1,mcellz+1]=head[1,1,1]

; done ...
  link_list=create_struct('head',head,'list',list)
  return,link_list
end
;----------------------------------------
function get_dv,xxi,xxj,vvi,vvj,lengths,xini,xfin,yini,yfin,zini,zfin
  Lxbox = xfin-xini
  Lybox = yfin-yini
  Lzbox = zfin-zini
  si=size(xxi)
  ni=si[1]
  sj=size(xxj)
  nj=si[1]
  si=size(lengths)
  nell=si[1]
  ell_max=max(lengths)
  mcellx=floor(Lxbox/ell_max)
  mcelly=floor(Lybox/ell_max)
  mcellz=floor(Lzbox/ell_max)
  cellx=Lxbox/mcellx
  celly=Lybox/mcelly
  cellz=Lzbox/mcellz
  link_listi=make_link_list(xxi,cellx,celly,cellz,mcellx,mcelly,mcellz,xini,yini,zini)
  link_listj=make_link_list(xxj,cellx,celly,cellz,mcellx,mcelly,mcellz,xini,yini,zini)
;---relative velocities for family i for the particles in same cell 
  for iz=1,mcellz  do begin
    for iy=1,mcelly  do begin
      for ix=1,mcellx  do begin
        ; within the same cell 
        ip = link_listi.head[ix,iy,iz]
        while (ip ne 0) do begin
          jp = link_listi.list[ip]
          while (jp ne 0) do begin
;            print,ip,"->",jp
            jp = link_listi.list[jp]
          endwhile 
          ip = link_listi.list[ip]
        endwhile 
;;
      end
    end
  end
;
;---relative velocities for family i for the particles in neighbouring cells 
;  for iz=1,mcellz  do begin
;    for iy=1,mcelly  do begin
;      for ix=1,mcellx  do begin
;
;        ip = link_listi.head[ix,iy,iz]
;				while (ip ne 0) do begin
;					for ineib = 0,12
;						ix2 = neighbours(ix,iy,iz,ineib,1)      
;						iy2 = neighbours(ix,iy,iz,ineib,2)
;						iz2 = neighbours(ix,iy,iz,ineib,3)
;						jp = link_listi.head[ix2,iy2,iz2]
;						while (jp ne 0) do begin
;							;print,ip,"->",jp
;							jp =link_listi.list[jp]       
;						endwhile     
;					end   
;					ip = link_listi.list[ip]
;				endwhile
;;
;			end
;		end
;	end

; create_structure("name",value) 
end


; main program is here ;
default,iread,0
default,vfile,'var'
vfile='VAR8'
if iread EQ 0 then begin
  pc_read_var,object=ff,varfile=vfile
  pc_read_ts,object=ts
  pc_read_pvar,object=fp
  pc_read_param,/param2,object=prun
  pc_read_param,object=pstart
  pc_read_dim,object=dim
  pc_read_grid,object=grid
  pc_read_kf,kf
  iread=1
end
iavg=100
si=size(ts.t)
itmax=si[1]-1
urms=mean(ts.urms[itmax-iavg:itmax])
epsK=mean(ts.epsK[itmax-iavg:itmax])
etaK=((prun.nu^3)/epsK)^(1./4.)
tauK=sqrt(prun.nu/epsK)
xini = grid.x[dim.l1]
xfin = grid.x[dim.l2] 
yini = grid.y[dim.m1]
yfin = grid.y[dim.m2] 
zini = grid.z[dim.n1]
zfin = grid.z[dim.n2] 
Lbox=grid.x[dim.l2]-grid.x[dim.l1]
tauL=1./(urms*kf)
Re=urms/(kf*prun.nu)
print,'urms,epsK,etaK,tauK,Lbox,tauL,Re='
print,urms,epsK,etaK,tauK,Lbox,tauL,Re
pradius= pstart.ap0
si=size(pradius)
nsp=si[1]
mp=(4./3.)*pstart.rhopmat*!pi*pradius^3
taup=(2./9.)*(pstart.rhopmat*pradius^2)/prun.nu
St=taup/tauK
print,'St='
print,St
print,nsp
lengths=etaK*[10.,20.,30.,40.,50.]
; do the following for each species
st1=where(fp.a[*] eq pradius[0])
x=fp.xx[st1,0]
y=fp.xx[st1,1]
z=fp.xx[st1,2]
vx=fp.vv[st1,0]
vy=fp.vv[st1,1]
vz=fp.vv[st1,2]
si=size(x)
xx1=findgen(si[1],3)
vv1=findgen(si[1],3)
for i1=0,si[1]-1 do begin
  xx1[i1,0]=x[i1]
  xx1[i1,1]=y[i1]
  xx1[i1,2]=z[i1]
  vv1[i1,0]=vx[i1]
  vv1[i1,1]=vy[i1]
  vv1[i1,2]=vz[i1]
end
nst1=si[1]
print,'nst1=',nst1
junk= get_dv(xx1,xx1,vv1,vv1,lengths,xini,xfin,yini,yfin,zini,zfin)
; function : input xxi,vvi,xxj,vvj,length[nell],Lbox
;           output is a structure that contains 
;           nell number of arrays, dv1 to dvnell
;           each of which contains the relative
;           velocities for distance length[nell] 
; function : get_dv 
;             
END