nx=60 & ny=60 & nz=60

; Rotation in the yz-plane, with angle alpha
alpha=0.

rho=fltarr(nx,ny,nz)
uu=fltarr(nx,ny,nz,3)
amplrho=10.
omega=1.0
radius_blob=0.25
xx=-0.5+findgen(nx)/float(nx-1)
yy=-0.5+findgen(ny)/float(ny-1)
zz=-0.5+findgen(nz)/float(nz-1)

for i=0,nx-1 do begin
for j=0,ny-1 do begin
for k=0,nz-1 do begin

rho(i,j,k)=1.0+amplrho*exp(-(xx(i)^2+yy(j)^2+zz(k)^2)/radius_blob^2)
rr=sqrt(xx(i)^2+yy(j)^2+zz(k)^2)
th=atan(xx(i)/zz(k))
ph=atan(-yy(j)/xx(i))
vphi=omega*rr*sin(th)

if (rr lt radius_blob) then begin
  uu(i,j,k,0)=vphi*sin(ph)
  uu(i,j,k,1)=vphi*cos(ph)*cos(alpha)
  uu(i,j,k,2)=vphi*cos(ph)*sin(alpha)
endif

endfor
endfor
endfor

; alpha=0, rotation in xy-plane vx,vy nonzero, vz zero
vel,uu(*,*,nz/2,0),uu(*,*,nz/2,1)

end