pro setup_jet
  fin='glass_10x10x10'
  readnew,fin,h,'HEAD'
  readnew,fin,xin,'POS'
  xin[*,*] /= h.BOXSIZE
  Nin = TOTAL(h.npart,/int)
  
  Ndup1 = 10L
  Ndup2 = 30L
  Gap = 4L
;  Gap = 0L
  
  fout='jet_100x300x10_light'
  
  Np_max = Nin * ndup1 * Ndup2
  x = fltarr(3,Np_max)
  v = fltarr(3,Np_max)
  u = fltarr(Np_max)
  m = fltarr(Np_max)
  id = lindgen(Np_max)

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

  gamma=5./3.
  rho = 1.0
  P = 1.0

  u_ini =  P/(gamma-1)/rho
  m_ini =  rho *  (h.BOXSIZE^3 / 1000 / Ndup1^2)
  dv = -1.0
;  dv = 0
  
  m[*] = m_ini
  u[*] = u_ini
  
  ; Lower 1/3
  i=0L
  FOR ix=0,Ndup1-1 DO BEGIN
     FOR iy=0,Ndup1-1 DO BEGIN
        FOR j=0L,Nin-1 DO BEGIN
           x[0,i] = (xin[0,j] + ix) / Ndup1 
           x[1,i] = (xin[1,j] + iy) / Ndup1
           x[2,i] = xin[2,j] / Ndup1
           ; Bottom wall
           IF x[1,i] LT 1. / Ndup1 THEN BEGIN
              m[i] = m_ini * 4
              u[i] = 1.1 * u_ini
              id[i] = 0
           END
           i++
        END
     END
  END

  ndup_thin=5L
 ; Upper 2/3 (hot air)
  FOR ix=0,ndup_thin-1 DO BEGIN
     FOR iy=0,2*ndup_thin-1 DO BEGIN
        FOR j=0L,Nin-1 DO BEGIN
           xx = (xin[0,j] + ix) / ndup_thin
           yy = (xin[1,j] + iy) / ndup_thin
           zz = xin[2,j] 
           IF (yy LE Gap/float(Ndup1)) OR (yy GT Gap/float(Ndup1) AND ABS(xx-0.5) GT 0.15) THEN BEGIN 
              x[0,i] = xx
              x[1,i] = yy + 1
              x[2,i] = zz / Ndup1
              u[i] = u_ini * Ndup1 / (ndup_thin/float(Ndup1)) 
              m[i] = m_ini / ndup_thin
              i++
           END
        END
     END
  END

  ; Upper 2/3 (tube walls)
  FOR iy=Gap,Ndup2-Ndup1-2 DO BEGIN
     FOR j=0L,Nin-1 DO BEGIN
        xx = xin[0,j] / Ndup1 + 0.35
        IF xx LT 0.44 THEN BEGIN
           x[0,i] = xx
           x[1,i] = (xin[1,j] + iy) / Ndup1 + 1
           x[2,i] = xin[2,j] / Ndup1
           id[i] = 0
           m[i] = m_ini * 4
           u[i] = 1.1 * u_ini
           i++
        END
        xx = xin[0,j] / Ndup1 + 0.55
        IF xx GT 0.56 THEN BEGIN
           x[0,i] = xx
           x[1,i] = (xin[1,j] + iy) / Ndup1 + 1
           x[2,i] = xin[2,j] / Ndup1
           id[i] = 0
           m[i] = m_ini * 4
           u[i] = 1.1 * u_ini
           i++
        END
     END
  END

  ; Upper 2/3 (top of tube, including moving part)
  FOR ix=0,2 DO BEGIN
     FOR j=0L,Nin-1 DO BEGIN
        x[0,i] = (xin[0,j] +ix)/ Ndup1 + 0.35
        x[1,i] = xin[1,j] / Ndup1 + 2.9
        x[2,i] = xin[2,j] / Ndup1
        id[i] = 0
        m[i] = m_ini * 8
        u[i] = 2 * u_ini
        IF x[0,i] GT 0.44 AND x[0,i] LE 0.56 THEN v[1,i] = dv
        i++
     END
  END

  ; Upper 2/3 (moving jet)
  FOR ix = 0,1 DO BEGIN
     FOR iz = 0,1 DO BEGIN
        FOR iy = 0,ROUND(2*(Ndup2-Ndup1-Gap-1) / 1.2) - 1 DO BEGIN
           FOR j=0L,Nin-1 DO BEGIN
              zz = (xin[2,j] + iz) / 2 / Ndup1 * 1.2
              IF zz LT 1.0 / Ndup1 THEN BEGIN
                 x[0,i] = (xin[0,j] + ix) / 2 / Ndup1 * 1.2 + 0.44
                 x[1,i] = (xin[1,j] + iy) / 2 / Ndup1 * 1.2 + 1 + Gap / float(Ndup1)
                 x[2,i] = zz
                 v[1,i] = dv
                 u[i] = u_ini / 8 * 1.2^3
                 i++
              END
           END
        END
     END
  END


  IF i GT Np_max THEN stop
  Np=i
  x=x[*,0:Np-1]
  v=v[*,0:Np-1]
  u=u[0:Np-1]
  m=m[0:Np-1]
  id=id[0:Np-1]
  
  h.npart[*]=0
  h.npart[0]=Np

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

;Upper part light fluid (hotter and less dnese
;  jj=WHERE(x[1,*] GT Ndup1)
;  u[jj] = u_ini * 5
;  m[jj] = m_ini / 5
  
;Solid bottom
;  jj=WHERE(x[1,*] LT dL)
;  id[jj] = 0
;  u[jj] *= 1.05

;Solid tube on top (don't forget to make it heavy again)
;  jj=WHERE((abs(x[0,*]-Ndup1/2-dL) LT dL/2 OR abs(x[0,*]-Ndup1/2+dL) LT dL/2) AND x[1,*] GT Ndup1+2*dL)
;  id[jj] = 0
;  u[jj] = u_ini * 1.05
;  m[jj] = m_ini

;Moving Wall 
;  jj=WHERE(abs(x[0,*]-Ndup1/2) LT dL/2 AND x[1,*] GT Ndup2-dL)
;  id[jj] = 0
;  u[jj] = u_ini * 1.05
;  m[jj] = m_ini
; v[1,jj] = dv

;Moving Jet 
;  jj=WHERE(abs(x[0,*]-Ndup1/2) LT dL/2 AND x[1,*] LE Ndup2-dL AND x[1,*] GE Ndup1 + 2*dL)
;  u[jj] = u_ini / 5
;  m[jj] = m_ini * 5
;  v[1,jj] = dv


; Bring positions to [0,1]
   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,*],x[1,*],psym=3
   jj=where(id[*] EQ 0)
   oplot,x[0,jj],x[1,jj],psym=3,col=128

stop
   
end