pro setup_nozzle
  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

  print,'Length of the region to fill is ',L

  fout='nozzle_4000x20x20_dv1'

  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.5-20)
   dL = 0.5

   FOR j=0L,Np-1 DO BEGIN
      IF x[0,j] GT L-dL THEN id[j]=0
   END

;Add nozzle in the left part with a length of dL = 4

   x0=0.5
   dL=8
   dmin=0.15
   amp=(0.5-dmin)/2

   FOR j=0L,Np-1 DO BEGIN
      IF x[0,j] GT x0 and x[0,j] LT x0+dL THEN BEGIN
         dy = abs(x[1,j]-0.5)
         dy_nozzle = 0.5 + amp * (sin((x[0,j] - x0)/dL*2*!Pi + !Pi/2) - 1)
         IF dy GT dy_nozzle THEN id[j] = 0
      END
   END

;Mark a wall at the beginning of the slab as boundaries (x=0-0.5) and
;give it a velocity
   dL = 0.5
   dv = -1
   FOR j=0L,Np-1 DO BEGIN
      IF x[0,j] LT dL THEN BEGIN
         id[j]=0
         v[0,j] = dv
      END
   END

; Move everything to the right by L
   FOR j=0L,Np-1 DO BEGIN
      x[0,j] += L
   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   '

;   plot the geometry to see that things are done correctly
   plot,x[0,*],x[1,*],xr=[0,40],psym=3
;   overplot with larger symbols the fixed particles
   jj=where(id[*] EQ 0)
   oplot,x[0,jj],x[1,jj],psym=3,col=128
;   overplot which particles if the walls are moving
   jj=where(id[*] EQ 0 AND v[0,*] NE 0)
   oplot,x[0,jj],x[1,jj],psym=1,col=128


end