function make_link_list,xx,cellx,cellz,mcellx,mcellz,xini,zini
  si=size(xx) ; 'xx' is the position of matched particles.
  help,xx
  print,'xx=',xx
  nn=si[1] ; # of matched particles found   
  print,'nn=',nn
  head=0*indgen(mcellx+2,mcellz+2,/LONG)
  list=0*indgen(nn,/LONG)
  for ip=0,nn-1 do begin  ;XY
	    icell=floor((xx[ip,0]-xini)/cellx)+1
	    help,icell
	    print,'icell=',icell
	    print,'xx[ip,0]',xx[ip,0]
	    kcell=floor((xx[ip,2]-zini)/cellz)+1
	    list[ip] = head[icell,kcell]
	    head[icell,kcell]=ip
  end   ;XY
  help,head
  print,'head=',head 
; apply PBC to head ...
; first the faces 
  ; now the edges 
  head[0,1:mcellz]        = head[mcellx,1:mcellz]

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

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

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

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

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

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

; done ...
  link_list=create_struct('head',head,'list',list)
  return,link_list
end

;----------------------------------------------------
function get_neighbours,ix,iz

	neighbours = 0*indgen(4,2,/LONG) 

	neighbours(0,0) = ix + 1
	neighbours(0,1) = iz

	neighbours(1,0) = ix - 1
	neighbours(1,1) = iz + 1  

	neighbours(2,0) = ix
	neighbours(2,1) = iz + 1

	neighbours(3,0) = ix + 1
	neighbours(3,1) = iz + 1

	return,neighbours
end
;----------------------------------------
function get_length,rpipj,lengths
	si = size(lengths)
	ilength = 0
	for i1 = 0,si[1]-1 do begin
		lst = lengths[i1]-0.05*lengths[i1]
		lend = lengths[i1]+0.05*lengths[i1]
		if ((rpipj gt lst) and (rpipj lt lend)) then ilength = i1+1
		if (ilength ne 0) then break
	end

	return,ilength 
end

;-----------------------------------------
function get_dv,xxi,xxj,vvi,vvj,lengths,xini,xfin,zini,zfin,etaK
  Lxbox = xfin-xini
  Lzbox = zfin-zini
	Lxby2 = 0.5*Lxbox
	Lzby2 = 0.5*Lzbox
  help,xxi
  print,'xxi=',xxi
  si=size(xxi)
  ni=si[1]
  sj=size(xxj)
  nj=si[1]
  si=size(lengths)
  nell=si[1]
  ell_max=max(lengths) ; maximum separation, to determine the "mcellx" and "mcelly"
  mcellx=floor(Lxbox/ell_max) ; number of cells in x direction
  mcellz=floor(Lzbox/ell_max)
  
  cellx=Lxbox/mcellx ;length of each cell
  cellz=Lzbox/mcellz
  link_listi=make_link_list(xxi,cellx,cellz,mcellx,mcellz,xini,zini)
  link_listj=make_link_list(xxj,cellx,cellz,mcellx,mcellz,xini,zini)
; opening files 
	for i1 = 0,nell-1 do begin
		fname = 'rel_vel/rel_vel_'+strtrim(string(lengths[i1]/etaK),1)+'eta'
		openw,i1+1,fname,/append 
	end

;---relative velocities for family i for the particles in same cell 
  for iz=1,mcellz  do begin
      for ix=1,mcellx  do begin
        ; within the same cell 
        ip = link_listi.head[ix,iz]
        while (ip ne 0) do begin
          jp = link_listi.list[ip]
          while (jp ne 0) do begin
		rx = xxi[jp,0]-xxi[ip,0]
		rz = xxi[jp,2]-xxi[ip,2]
		rpipj = sqrt(rx*rx+rz*rz)

		ilength = get_length(rpipj,lengths)
		if (ilength ne 0) then begin
			dvx = vvi[jp,0]-vvi[ip,0]
			dvz = vvi[jp,2]-vvi[ip,2]

			dv = sqrt(dvx*dvx+dvz*dvz)
			dvl = (rx*dvx+rz*dvz)/rpipj
			printf,ilength,dv,dvl
		endif
            jp = link_listi.list[jp]
          endwhile 
          ip = link_listi.list[ip]
        endwhile 
;
      end
  end
;
;---relative velocities for family i for the particles in neighbouring cells 
  for iz=1,mcellz  do begin
      for ix=1,mcellx  do begin

        ip = link_listi.head[ix,iz]
	neighbours = get_neighbours(ix,iz)
	while (ip ne 0) do begin
		for ineib = 0,12 do begin						
			ix2 = neighbours(ineib,0)      
			iz2 = neighbours(ineib,2)
			jp = link_listi.head[ix2,iz2]
			while (jp ne 0) do begin

				rx = xxi[jp,0]-xxi[ip,0]
				q = floor(rx/Lxby2) 	
				p = 0
				if (q eq 1) then p = -1  		
				if (q eq -2) then p = 1  		
				rx = rx + Lxbox*p
				p = 0
				if (q eq 1) then p = -1  		
				if (q eq -2) then p = 1  		
				rz = xxi[jp,2]-xxi[ip,2]
				q = floor(rz/Lzby2) 	
				p = 0
				if (q eq 1) then p = -1  		
				if (q eq -2) then p = 1  		
				rz = rz + Lzbox*p
				rpipj = sqrt(rx*rx+rz*rz)

				ilength = get_length(rpipj,lengths)	 
				if (ilength ne 0) then begin
					dvx = vvi[jp,0]-vvi[ip,0]
					dvz = vvi[jp,2]-vvi[ip,2]

					dv = sqrt(dvx*dvx+dvz*dvz)
					dvl = (rx*dvx+rz*dvz)/rpipj
					printf,ilength,dv,dvl
				endif

				jp = link_listi.list[jp]       
			endwhile     
		end   
		ip = link_listi.list[ip]
	endwhile

end
	end

	for i1 = 0,nell-1 do begin
		close,i1+1 
	end

end


; main program is here ;
default,iread,0
default,vfile,'var'
vfile='VAR8'
if iread EQ 0 then begin
  pc_read_pvar,object=fp,varfile=vfile
  pc_read_ts,object=ts
  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
SPAWN,'rm  rel_vel/*'
SPAWN,'mkdir  rel_vel'

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] ; left edge of the box 
print,'xini=',xini
xfin = grid.x[dim.l2] ; right edge of the box
print,'xfin=',xfin
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  ;XY

pradius=min(fp.a) ;XY
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*[0.5,1.,2.,5.,10] ; seperations
; loop over snapshots

;for ipread = 1,2 do begin
for ipread = 0,0 do begin
	pvfile = 'PVAR'+strtrim(string(ipread),1) 
	pc_read_pvar,object=fp,varfile=pvfile
	; do the following for each species
	;st1=where(fp.a[*] eq pradius[4]) ;XY
        ;st1=where(fp.a[*] eq pradius[0]) ;XY
        st1=where(fp.a[*] eq fp.a[1]) ;XY
        if st1 GE 0 then begin ;xy
		print,'st1=',st1
		print,'mimmax(fp.xx)',minmax(fp.xx)
		x=fp.xx[st1,0]
		print,'x=',x
		print,'minmax(x)',minmax(x)
		z=fp.xx[st1,2]
		vx=fp.vv[st1,0]
		vz=fp.vv[st1,2]
		si=size(x)
		print,'si=',si
		xx1=findgen(si[1],3)
		help,xx1
		print,'xx1=',xx1
		vv1=findgen(si[1],3)
		print,'si[1]-1',si[1]-1
		for i1=0,si[1]-1 do begin ;XY
			;i1=0  ; XY
			xx1[i1,0]=x[i1]
			xx1[i1,2]=z[i1]
			vv1[i1,0]=vx[i1]
			vv1[i1,2]=vz[i1]
		end
		nst1=si[1]
		print,'nst1=',nst1

		junk = get_dv(xx1,xx1,vv1,vv1,lengths,xini,xfin,zini,zfin,etaK)
            endif else begin
                  print,'do not find the matched particles' 
            endelse	
  
end

END