pro setup_test2
  fin='glass_10x10x10'
  readnew,fin,h,'HEAD'
  readnew,fin,xin,'POS'
  xin[*,*] /= h.BOXSIZE
  Nin = TOTAL(h.npart,/int)
  
  Ndup = 40L
  Ndup2 = 2L
  L = Ndup / Ndup2

  ; Size of the nozzle
  a=0.175
  
  fout='glass_400x20x20'
  
  Np = Ndup * Ndup2 * Ndup2 * Nin
  x=fltarr(3,Np)

  i=0L
  FOR ix=0,Ndup-1 DO BEGIN
     FOR iy=0,Ndup2-1 DO BEGIN
        FOR iz=0,Ndup2-1 DO BEGIN
           FOR j=0L,Nin-1 DO BEGIN
              x[0,i] = (xin[0,j] + ix) / NDup2
              x[1,i] = (xin[1,j] + iy) / NDup2
              x[2,i] = (xin[2,j] + iz) / NDup2
              i++
           END
        END
     END
  END
  
  v = fltarr(3,Np)
  u=fltarr(Np)
  m=fltarr(Np)
  id=lindgen(Np)

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

  gamma=5./3.
  rho=0.125
  P = 1.0
  
  u[*] = P/(gamma-1)/rho
  mpart = (rho * Ndup2^3 * h.BOXSIZE^3)/Np
  m[*] = mpart

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

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

;Mark a wall at the end of the slab as boundaries (x=19-20)  
   dv = -1.0
   v[0,*] = dv
   dL = 1.0
   dL2 = 6.0

; Split the wall in two parts, where one part is at rest and has a hole
   FOR j=0L,Np-1 DO BEGIN
      IF x[0,j] GT L-dL AND x[0,j] LT L-dL/2 THEN BEGIN
         id[j]=0
         m[j] = mpart * 5
      END
      IF x[0,j] GT L-dL/2 THEN BEGIN
         rr=sqrt((x[1,j]-0.5)^2 + (x[2,j]-0.5)^2)
         IF rr GT a THEN BEGIN
            id[j]=0
            m[j] = mpart * 5
            v[0,j] = 0
         END
      END
   END

;Add nozzle between 10 and 13 as boundary and set it velocity to zero
   x00=8
   x0=12
   x1=15
   
   FOR j=0L,Np-1 DO BEGIN
      IF x[0,j] GT x0 and x[0,j] LT x1 THEN BEGIN
         rr=sqrt((x[1,j]-0.5)^2 + (x[2,j]-0.5)^2)
         d = (x[0,j] - x0) / (x1 - x0) 
         rmax = a + d * (0.5-a)
         IF rr GT rmax THEN BEGIN
            id[j] = 0
            m[j] = mpart * 5
            v[0,j] = 0
         END
      END
      IF x[0,j] LE x0 AND x[0,j] GT x00 THEN BEGIN
         rr=sqrt((x[1,j]-0.5)^2 + (x[2,j]-0.5)^2)
         IF rr GT  a THEN BEGIN
            id[j] = 0
            m[j] = mpart * 5
            v[0,j] = 0
         END ELSE BEGIN
            u[j] = P/(gamma-1)/rho / 17. * 12.
         END
      END
      IF x[0,j] LT x00 THEN BEGIN
         v[0,j] = 0
      END
     
   END
    
   write_head,fout,h
   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   '

   jj=where(abs(x[2,*]-0.5) LT 0.1)
   plot,x[0,jj],x[1,jj],psym=3
   jj=where(abs(x[2,*]-0.5) LT 0.1 AND id[*] EQ 0)
   oplot,x[0,jj],x[1,jj],psym=3,col=128

   
end