pro setup_slab
   name='glass_10x10x10'

   fout='mhd_shock.ic'

   gamma=5./3.
   ll=50L

   nl=2L
   rhosetupl=1.0
   Psetupl=1.0

   nr=1L
   rhosetupr=0.125
   Psetupr=0.1

   readnew,name,h,'HEAD'
   readnew,name,xglass,'POS'
   readnew,name,hglass,'HSML'

   ntot=h.npart(0)*ll*(nl^3+nr^3)

   icount=0LL

   x=fltarr(3,ntot)
   v=fltarr(3,ntot)
   b=fltarr(3,ntot)
   id=lindgen(ntot)+1
   u=fltarr(ntot)
   rho=fltarr(ntot)
   hsml=fltarr(ntot)

   xbase=xglass/h.boxsize/nl
   hbase=hglass/h.boxsize/nl

   FOR ii=0,ll-1 DO BEGIN
      FOR ix=0,nl-1 DO BEGIN
         x0=float(ix)/nl+ii
         FOR iy=0,nl-1 DO BEGIN
            y0=float(iy)/nl
            FOR iz=0,nl-1 DO BEGIN
               z0=float(iz)/nl
               FOR i=0,h.npart(0)-1 DO BEGIN
                  x(0,icount)=xbase(0,i)+x0
                  x(1,icount)=xbase(1,i)+y0
                  x(2,icount)=xbase(2,i)+z0
                  v(0,icount)=0.0
                  v(1,icount)=0.0
                  v(2,icount)=0.0
                  b(0,icount)=0.75
                  b(1,icount)=1.0
                  b(2,icount)=0.0
                  u(icount)=Psetupl/(gamma-1)/rhosetupl
                  hsml(icount)=hbase(i)
                  rho(icount)=rhosetupl
                  icount++
               END
            END
         END
      END
   END

   xbase=xglass/h.boxsize/nr
   hbase=hglass/h.boxsize/nr
   FOR ii=0,ll-1 DO BEGIN
      FOR ix=0,nr-1 DO BEGIN
         x0=float(ix)/nr+ii+ll
         FOR iy=0,nr-1 DO BEGIN
            y0=float(iy)/nr
            FOR iz=0,nr-1 DO BEGIN
               z0=float(iz)/nr
               FOR i=0,h.npart(0)-1 DO BEGIN
                  x(0,icount)=xbase(0,i)+x0
                  x(1,icount)=xbase(1,i)+y0
                  x(2,icount)=xbase(2,i)+z0
                  v(0,icount)=0.0
                  v(1,icount)=0.0
                  v(2,icount)=0.0
                  b(0,icount)=0.75
                  b(1,icount)=-1.0
                  b(2,icount)=0.0
                  u(icount)=Psetupr/(gamma-1)/rhosetupr
                  hsml(icount)=hbase(i)
                  rho(icount)=rhosetupr
                  icount++
               END
            END
         END
      END
   END

   print,icount," ==?== ",ntot

   h.massarr(0)=rhosetupl/(nl^3*h.npart(0))
   h.npart(0)=ntot
   h.parttotal(0)=ntot

   write_head,fout,h
   add_block,fout,x,'POS '
   add_block,fout,v,'VEL '
   add_block,fout,id,'ID  '
   add_block,fout,u,'U   '
   add_block,fout,hsml,'HSML'
   add_block,fout,b,'BFLD'

end