import numpy as np
import netCDF4 as nc

n = 100
data = np.zeros((n,n,n), dtype=float)
for i in range(n):
    x = 15.*((i+0.5)/float(n)-0.5)
    for j in range(n):
        y = 15.*((j+0.5)/float(n)-0.5)
        for k in range(n):
            z = 15.*((k+0.5)/float(n)-0.5)
            data[i][j][k] = ((np.sin(np.sqrt(y*y+x*x)))**2-0.5)/(0.001*(y*y+x*x)+1.)**2 + 0.5 + \
                            ((np.sin(np.sqrt(z*z+y*y)))**2-0.5)/(0.001*(z*z+y*y)+1.)**2 + 0.5

f = nc.Dataset('sineEnvelope.nc', 'w', format='NETCDF4')
f.createDimension('x', n)
f.createDimension('y', n)
f.createDimension('z', n)
rho = f.createVariable('density', 'f4', ('x', 'y', 'z'))
rho[:,:,:] = data
f.close()
