pro setup_slab
   fin1='grid_10x10x10'
   fin2='grid_12x12x12'
  
   readnew,fin1,h1,'HEAD'
   readnew,fin1,xin1,'POS'
   readnew,fin2,h2,'HEAD'
   readnew,fin2,xin2,'POS'
  
   xin1[*,*] /= h1.BOXSIZE
   Nin1 = TOTAL(h1.npart,/int)

   xin2[*,*] /= h2.BOXSIZE
   Nin2 = TOTAL(h2.npart,/int)

   Ndup = 20L
   fout='slab.ic'
  
   Np = Ndup * (Nin1+Nin2)
   x=fltarr(3,Np)

   i=0L
   FOR ix=0,Ndup-1 DO BEGIN
      FOR j=0L,Nin1-1 DO BEGIN
         x[0,i] = xin1[0,j] + ix 
         x[1,i] = xin1[1,j]
         x[2,i] = xin1[2,j]
         i++
      END
   END

   FOR ix=0,Ndup-1 DO BEGIN
      FOR j=0L,Nin2-1 DO BEGIN
         x[0,i] = xin2[0,j] + Ndup + ix 
         x[1,i] = xin2[1,j]
         x[2,i] = xin2[2,j]
         i++
      END
   END

   v = fltarr(3,Np)
   u=fltarr(Np)
   m=fltarr(Np)
   id=lindgen(Np)

   h1.BOXSIZE=1
   h1.MASSARR[*]=0

   m[*]=1.0/Nin1

   rho1 = m[0] * Nin1
   rho2 = m[0] * Nin2
  
   gamma=5./3.
   P1 = 1.0
   P2 = 0.5

   jj=where(x[0,*] LT Ndup)
   u[jj] = P1/(gamma-1)/rho1

   jj=where(x[0,*] GE Ndup)
   u[jj] = P2/(gamma-1)/rho2

   h1.npart[*]=0
   h1.npart[0]=Np

   h1.PARTTOTAL[*]=0
   h1.PARTTOTAL[0]=Np

   write_head,fout,h1
   add_block,fout,x,'POS '
   add_block,fout,v,'VEL '
   add_block,fout,id,'ID  '
   add_block,fout,m,'MASS'
   add_block,fout,u,'U   '

end