pro setup_blast
  fin='glass_10x10x10'
  readnew,fin,h,'HEAD'
  readnew,fin,xin,'POS'
  xin[*,*] /= h.BOXSIZE
  Nin = TOTAL(h.npart,/int)
  
  Ndup = 10L
  
  fout='blast_100x100x100'
  
  Np = Ndup * Ndup * Ndup * Nin
  x=fltarr(3,Np)

  i=0L
  FOR ix=0,Ndup-1 DO BEGIN
     FOR iy=0,Ndup-1 DO BEGIN
        FOR iz=0,Ndup-1 DO BEGIN
           FOR j=0L,Nin-1 DO BEGIN
              x[0,i] = (xin[0,j] + ix) / NDup
              x[1,i] = (xin[1,j] + iy) / NDup
              x[2,i] = (xin[2,j] + iz) / NDup
              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
  Pback = 1e-20
  Pexpl = 1

  v[*,*] = 0
  u[*] = Pback/(gamma-1)/rho
  mpart = (rho * h.BOXSIZE^3)/Np
  m[*] = mpart

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

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

; find central region
  rr=sqrt((x[0,*]-0.5)^2+(x[1,*]-0.5)^2+(x[2,*]-0.5)^2)
  ii=SORT(rr)
  n_ini=20
  u[ii[0:n_ini-1]] = Pexpl/(gamma-1)/rho
  
   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,x[0,*],x[1,*],psym=3
   jj=where(u GT min(u))
   oplot,x[0,jj],x[1,jj],psym=4,symsize=2,thick=3,col=128

   
end