#! /usr/bin/env python

import math

def __main__():
  nx = 50
  ny = 50
  nz = 50
  print_metadata(nx, ny, nz)
  for i in range(nx):
    for j in range(ny):
      for k in range(nz):
        rho = density_func(i, j, k)
        print rho, '   ',

## a function to generate some sample densities for mayavi to visualize
def density_func(x, y, z):
  off1 = (5, 20, 25)
  off2 = (30, 10, 5)
  f1 = math.exp(-((x-off1[0])*(x-off1[0])/2.0 + (y-off1[1])*(y-off1[0])
                  + (z-off1[2])*(z-off1[2]))/50.0)
  f2 = math.exp(-((x-off2[0])*(x-off2[0]) + (y-off2[1])*(y-off2[0]) + (z-off2[2])*(z-off2[2]))/40.0)
  return f1 + f2

def print_metadata(nx, ny, nz):
  print '# vtk DataFile Version 2.0 --- %s' % 'mayavi sample'
  print 'Volume example'
  print 'ASCII'
  print 'DATASET STRUCTURED_POINTS'
  print 'DIMENSIONS %d %d %d' % (nx, ny, nz)
  print 'ASPECT_RATIO %g %g %g' % (1.0/nx, 1.0/ny, 1.0/nz)
  print 'ORIGIN 0 0 0'
  print 'POINT_DATA %d' % (int) (nx*ny*nz)
  print 'SCALARS volume_scalars double 1'
  print 'LOOKUP_TABLE default'

__main__()