
n = 30;

min = -2;
max = 2;
del = (max - min)/(n-1);

grid = linspace(min,max,n);

printf("# vtk DataFile Version 2.0\n");
printf("Flow near sphere\n");
printf("ASCII\n");

printf("DATASET STRUCTURED_POINTS\n");
printf("DIMENSIONS  %d %d %d\n", n, n, n);
printf("ORIGIN    %.3g %.3g %.3g\n",min,min,min);
printf("SPACING   %.3g %.3g %.3g\n",del,del,del);

printf("POINT_DATA   %d\n",n^3);
printf("VECTORS flow float\n");
      
for z = grid
  for y = grid
    for x = grid
      r = sqrt(x^2 + y^2 + z^2);
      if r>1
	v = [1 + 0.5/r^3 - 1.5*x^2/r^5, -1.5*x*y/r^5, -1.5*x*z/r^5];
      else
	v = [0,0,0];
      endif
      printf(" %g", v);
      printf("\n");
    endfor
  endfor
endfor

printf("SCALARS potential float\n");
printf("LOOKUP_TABLE default\n");
     
for z = grid
  for y = grid
    for x = grid
      r = sqrt(x^2 + y^2 + z^2);
      if r>1
	p = x*(1 + 0.5/r^3);
      else
	p = 0;
      endif
      printf(" %g", p);
      printf("\n");
    endfor
  endfor
endfor
