;if !d.name eq 'PS' then begin
;  device,filename='idl.ps',xsize=18,ysize=18,yoffset=3,/color
;  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
;end
;
;!p.charsize=3.0
;!x.margin=[4.,1.]
;!y.margin=[2.,.5]
;
;mention directory
FILE_MKDIR, 'box_count'
directory='./box_count/'
; Read particle data and map them on the grid
;
pc_read_param,obj=start
boxsize=start.xyz1
boxsize=boxsize(0)
pc_read_dim,obj=dim
nx=dim.nx
pc_read_pdim,obj=pdim
openr,1,'x.txt'
x=fltarr(nx)
readf,1,x
close,1

ivar0=0
deltaivar=1
ivar1=3


for ivar=ivar0,ivar1 do begin
  restore,'./pvar_file/'+'pvar'+str(ivar)+'.sav'
  print,'t=',t
	;restore,'pvar_file/pvar2.sav'
	;map particles on the grid
	;pc_map_particles_on_grid, pvar.xx, pvar.x, pvar.y, pvar.z,$
	pc_map_particles_on_grid, xx, x, x, x,ineargrid=in
	;
	; set some defaults
	;
	default,width,1
	print,'width=',width
	;
	; Make auxillaries
	;
	box=fltarr(dim.nx+2*width,dim.ny+2*width,dim.nz+2*width)

	; Loop over all particles
	;
	print,'Looping over all particles....'
	for ipar=0,pdim.npar-1 do begin
			ix=in(ipar,0)+width & iy=in(ipar,1)+width & iz=in(ipar,2)+width
			box(ix-width:ix+width,iy-width:iy+width,iz-width:iz+width)=$
			box(ix-width:ix+width,iy-width:iy+width,iz-width:iz+width)+1
	end
	;
	; Fold data from ghost zones
	;
	l1=width & l2=l1+dim.nx-1
	m1=width & m2=m1+dim.ny-1
	n1=width & n2=n1+dim.nz-1
	l1i=l1+width-1
	l2i=l2-width+1
	m1i=m1+width-1
	m2i=m2-width+1
	n1i=n1+width-1
	n2i=n2-width+1
	mx=dim.nx+2*width-1
	my=dim.ny+2*width-1
	mz=dim.nz+2*width-1
	print,'Folding ghost zones....'
	if (width gt 0) then begin
			; x-dir
			box(l1:l1i,0:my,0:mz)=box(l1:l1i,0:my,0:mz)+box(l2+1:mx,0:my,0:mz)
			box(l2i:l2,0:my,0:mz)=box(l2i:l2,0:my,0:mz)+box(0:l1-1,0:my,0:mz)
			box(l2+1:mx,0:my,0:mz)=0
			box(0:l1-1,0:my,0:mz)=0
			; y-dir
			box(0:mx,m1:m1i,0:mz)=box(0:mx,m1:m1i,0:mz)+box(0:mx,m2+1:my,0:mz)
			box(0:mx,m2i:m2,0:mz)=box(0:mx,m2i:m2,0:mz)+box(0:mx,0:m1-1,0:mz)
			box(0:mx,m2+1:my,0:mz)=0
			box(0:mx,0:m1-1,0:mz)=0
			; z-dir
			box(0:mx,0:my,n1:n1i)=box(0:mx,0:my,n1:n1i)+box(0:mx,0:my,n2+1:mz)
			box(0:mx,0:my,n2i:n2)=box(0:mx,0:my,n2i:n2)+box(0:mx,0:my,0:n1-1)
			box(0:mx,0:my,n2+1:mz)=0
			box(0:mx,0:my,0:n1-1)=0
	endif
	;
	; Normalize box values and removed ghose cells
	;
	norm=(1+2*width)^3
	box=box(l1:l2,n1:n2,m1:m2)/norm
	;save,file='box10.sav',box,nx
	save,file=directory+'box'+str(ivar)+'.sav',box,nx
	;
	; Print some relevant data
	;
	print,'minmax(box)=',minmax(box)
	print,'total(box)=',total(box)
	print,'npar=',pdim.npar
	;
;	; Plot data
;	;
;	!x.title='!6'
;	!y.title='!6'
;	contour,box(*,*,nx/2),/fill,nlev=60,/iso
endfor
END