pro close_packed_lattice_new

!p.multi=[0,3,1]

lz=1L

n=12L
r=0.6

rho_target=0.125

fout="packed_12x12x12"


dx=2*r
dy=sqrt(3.)*r

x0=[0.,r]

nx=fix((n/r)/4)*2L
ny=fix((n/r)/4)*2L

ylatice=0

n2d=(nx+4L)*(ny+4L)
x2d=fltarr(2,n2d)

i=0L
for iy=-1,ny+2 do begin
   for ix=-1,nx+2 do begin
      x2d(0,i)=x0[ylatice]+ix*dx
      x2d(1,i)=iy*dy
      i=i+1
   end
   ylatice=1-ylatice
end

dz_z=r*2.
dx_z=3.*r/3.
dy_z=sqrt(3)*r/3.
nz=fix((n*lz/dz_z)/3*n*n/nx/ny)

ntot=n2d*nz*3L
x=fltarr(3,ntot)
col=intarr(ntot)

for i=0L,nz-1 do begin
   for j=0L,2 do begin
      iz=i*3+j
      x(0,iz*n2d:(iz+1)*n2d-1)=x2d(0,*)+j*dx_z + dx_z / 2
      x(1,iz*n2d:(iz+1)*n2d-1)=x2d(1,*)+j*dy_z + dy_z / 4
      x(2,iz*n2d:(iz+1)*n2d-1)=iz*dz_z + dz_z / 2
      col(iz*n2d:(iz+1)*n2d-1)=j
   end
end

xbox=[0,1,1,0,0]
ybox=[0,0,1,1,0]

plot,xbox,ybox,xr=[-0.1,1.1],yr=[-0.1,1.1],xst=1,yst=1,xtit='x',ytit='y'
oplot,x(0,*)/n,x(1,*)/n,psym=4,col=256-64

plot,xbox,ybox,xr=[-0.1,1.1],yr=[-0.1,1.1],xst=1,yst=1,xtit='z',ytit='y'
oplot,x(2,*)/n,x(1,*)/n,psym=4,col=256-64

plot,xbox,ybox,xr=[-0.1,1.1],yr=[-0.1,1.1],xst=1,yst=1,xtit='x',ytit='z'
oplot,x(0,*)/n,x(2,*)/n,psym=4,col=256-64

print,n,n,n*lz
print,nx*dx,ny*dy,nz*3*dz_z

eps=0.01

x(0,*)=x(0,*)*n/(nx*dx) + r/2
x(1,*)=x(1,*)*n/(ny*dy) + r/2
x(2,*)=x(2,*)*n*lz/(nz*dz_z*3) + r/2

jj=where((x(0,*) GE 0) AND (x(0,*) LT n) AND $
         (x(1,*) GE 0) AND (x(1,*) LT n) AND $
         (x(2,*) GE 0) AND (x(2,*) LT n*LZ))

x=x(*,jj)/n
nfinal=N_ELEMENTS(jj)
pos=fltarr(3,nfinal)
pos(0,*)=x(2,*)
pos(1,*)=x(1,*)
pos(2,*)=x(0,*)
vel=fltarr(3,nfinal)
vel(*,*)=0
vel(0,*)=0.1
mass=fltarr(nfinal)
id=lindgen(nfinal)+1
u=fltarr(nfinal)
u(*)=1.0

print,N_ELEMENTS(jj),n^3*LZ

m=1.0/n^3*rho_target
mass(*)=m

print,"mass =",m

rho=N_ELEMENTS(jj)*m/lz

print,"rho =",rho,rho_target

npart=lonarr(6)
massarr=dblarr(6)
time=0.0D
redshift=0.0D
flag_sfr=0L
flag_feedback=0L
partTotal=lonarr(6)
flag_cooling=0L
num_files=1L
BoxSize=1.0D
Omega0=0.0D
OmegaLambda=0.0D
HubbleParam=0.7D
flag_stellarage=0L ;          /*!< flags whether the file contains formation times of star particles */
flag_metals=0L ;              /*!< flags whether the file contains metallicity values for gas and star particles */ 
npartTotalHighWord=lonarr(6)    ;
Labels=bytarr(2,15)
bytesleft=256-6*4 - 6*8 - 8 - 8 - 2*4 - 6*4 - 2*4 - 4*8 - 2*4 - 6*4 - 2*15
la=bytarr(bytesleft)

   h = { head , npart:npart,$
	 massarr:massarr,$
	 time:time,$
	 redshift:redshift,$
	 flag_sfr:flag_sfr,$
	 flag_feedback:flag_feedback,$
	 partTotal:partTotal,$
	 flag_cooling:flag_cooling,$
	 num_files:num_files,$
	 BoxSize:BoxSize,$
	 Omega0:Omega0,$
	 OmegaLambda:OmegaLambda,$
         HubbleParam:HubbleParam,$
         flag_stellarage:flag_stellarage,$
         flag_metals:flag_metals,$
         npartTotalHighWord:npartTotalHighWord,$
         Labels:Labels,$
         la:la}

   h.npart(0)=nfinal
   h.partTotal(0)=nfinal

   write_head,fout,h
   add_block,fout,pos,'POS '
   add_block,fout,vel,'VEL '
   add_block,fout,id,'ID  '
   add_block,fout,mass,'MASS'
   add_block,fout,u,'U   '

stop
end