docs: def createCellsDensity (self):
          ''' Create population cells based on density'''
          from .. import sim
  
          cells = []
          shape = sim.net.params.shape
          sizeX = sim.net.params.sizeX
          sizeY = sim.net.params.sizeY
          sizeZ = sim.net.params.sizeZ
          
          # calculate volume
          if shape == 'cuboid':
              volume = sizeY/1e3 * sizeX/1e3 * sizeZ/1e3  
          elif shape == 'cylinder':
              volume = sizeY/1e3 * sizeX/1e3/2 * sizeZ/1e3/2 * pi
          elif shape == 'ellipsoid':
              volume = sizeY/1e3/2.0 * sizeX/1e3/2.0 * sizeZ/1e3/2.0 * pi * 4.0 / 3.0
  
          for coord in ['x', 'y', 'z']:
              if coord+'Range' in self.tags:  # if user provided absolute range, convert to normalized
                  self.tags[coord+'normRange'] = [point / getattr(sim.net.params, 'size'+coord.upper()) for point in self.tags[coord+'Range']]
              if coord+'normRange' in self.tags:  # if normalized range, rescale volume
                  minv = self.tags[coord+'normRange'][0] 
                  maxv = self.tags[coord+'normRange'][1] 
                  volume = volume * (maxv-minv)
  
          funcLocs = None  # start with no locations as a function of density function
          if isinstance(self.tags['density'], basestring): # check if density is given as a function 
              if shape == 'cuboid':  # only available for cuboids
                  strFunc = self.tags['density']  # string containing function
                  strVars = [var for var in ['xnorm', 'ynorm', 'znorm'] if var in strFunc]  # get list of variables used 
                  if not len(strVars) == 1:
                      print('Error: density function (%s) for population %s does not include "xnorm", "ynorm" or "znorm"'%(strFunc,self.tags['pop']))
                      return
                  coordFunc = strVars[0] 
                  lambdaStr = 'lambda ' + coordFunc +': ' + strFunc # convert to lambda function 
                  densityFunc = eval(lambdaStr)
                  minRange = self.tags[coordFunc+'Range'][0]
                  maxRange = self.tags[coordFunc+'Range'][1]
  
                  interval = 0.001  # interval of location values to evaluate func in order to find the max cell density
                  maxDensity = max(list(map(densityFunc, (np.arange(minRange, maxRange, interval)))))  # max cell density 
                  maxCells = volume * maxDensity  # max number of cells based on max value of density func 
                  
                  self.rand.Random123(int(maxDensity), sim.net.lastGid, sim.cfg.seeds['loc'])
                  locsAll = minRange + ((maxRange-minRange)) * np.array([self.rand.uniform(0, 1) for i in range(int(maxCells))])  # random location values 
                  locsProb = np.array(list(map(densityFunc, locsAll))) / maxDensity  # calculate normalized density for each location value (used to prune)
                  allrands = np.array([self.rand.uniform(0, 1) for i in range(len(locsProb))])  # create an array of random numbers for checking each location pos 
                  
                  makethiscell = locsProb>allrands  # perform test to see whether or not this cell should be included (pruning based on density func)
                  funcLocs = [locsAll[i] for i in range(len(locsAll)) if i in np.array(makethiscell.nonzero()[0],dtype='int')] # keep only subset of yfuncLocs based on density func
                  self.tags['numCells'] = len(funcLocs)  # final number of cells after pruning of location values based on density func
                  if sim.cfg.verbose: print('Volume=%.2f, maxDensity=%.2f, maxCells=%.0f, numCells=%.0f'%(volume, maxDensity, maxCells, self.tags['numCells']))
              else:
                  print('Error: Density functions are only implemented for cuboid shaped networks')
                  exit(0)
          else:  # NO ynorm-dep
              self.tags['numCells'] = int(self.tags['density'] * volume)  # = density (cells/mm^3) * volume (mm^3)
  
          # calculate locations of cells 
          self.rand.Random123(self.tags['numCells'], sim.net.lastGid, sim.cfg.seeds['loc'])
          self.rand.uniform(0, 1)
          vec = h.Vector(self.tags['numCells']*3)
          vec.setrand(self.rand)
          randLocs = np.array(vec).reshape(self.tags['numCells'], 3)  # create random x,y,z locations
  
          if sim.net.params.shape == 'cylinder':
              # Use the x,z random vales 
              rho = randLocs[:,0] # use x rand value as the radius rho in the interval [0, 1)
              phi = 2 * pi * randLocs[:,2] # use z rand value as the angle phi in the interval [0, 2*pi) 
              x = (1 + sqrt(rho) * cos(phi))/2.0
              z = (1 + sqrt(rho) * sin(phi))/2.0
              randLocs[:,0] = x
              randLocs[:,2] = z
      
          elif sim.net.params.shape == 'ellipsoid':
              # Use the x,y,z random vales 
              rho = np.power(randLocs[:,0], 1.0/3.0) # use x rand value as the radius rho in the interval [0, 1); cuberoot
              phi = 2 * pi * randLocs[:,1] # use y rand value as the angle phi in the interval [0, 2*pi) 
              costheta = (2 * randLocs[:,2]) - 1 # use z rand value as cos(theta) in the interval [-1, 1); ensures uniform dist 
              theta = arccos(costheta)  # obtain theta from cos(theta)
              x = (1 + rho * cos(phi) * sin(theta))/2.0
              y = (1 + rho * sin(phi) * sin(theta))/2.0
              z = (1 + rho * cos(theta))/2.0 
              randLocs[:,0] = x
              randLocs[:,1] = y
              randLocs[:,2] = z
  
          for icoord, coord in enumerate(['x', 'y', 'z']):
              if coord+'normRange' in self.tags:  # if normalized range, rescale random locations
                  minv = self.tags[coord+'normRange'][0] 
                  maxv = self.tags[coord+'normRange'][1] 
                  randLocs[:,icoord] = randLocs[:,icoord] * (maxv-minv) + minv
              if funcLocs and coordFunc == coord+'norm':  # if locations for this coordinate calculated using density function
                  randLocs[:,icoord] = funcLocs
  
          if sim.cfg.verbose and not funcLocs: print('Volume=%.4f, density=%.2f, numCells=%.0f'%(volume, self.tags['density'], self.tags['numCells']))
  
          for i in self._distributeCells(self.tags['numCells'])[sim.rank]:
              gid = sim.net.lastGid+i
              self.cellGids.append(gid)  # add gid list of cells belonging to this population - not needed?
              cellTags = {k: v for (k, v) in self.tags.items() if k in sim.net.params.popTagsCopiedToCells}  # copy all pop tags to cell tags, except those that are pop-specific
              cellTags['pop'] = self.tags['pop']
              cellTags['xnorm'] = randLocs[i,0]  # calculate x location (um)
              cellTags['ynorm'] = randLocs[i,1]  # calculate y location (um)
              cellTags['znorm'] = randLocs[i,2]  # calculate z location (um)
              cellTags['x'] = sizeX * randLocs[i,0]  # calculate x location (um)
              cellTags['y'] = sizeY * randLocs[i,1]  # calculate y location (um)
              cellTags['z'] = sizeZ * randLocs[i,2]  # calculate z location (um)
              cells.append(self.cellModelClass(gid, cellTags)) # instantiate Cell object
              if sim.cfg.verbose: 
                  print(('Cell %d/%d (gid=%d) of pop %s, pos=(%2.f, %2.f, %2.f), on node %d, '%(i, self.tags['numCells']-1, gid, self.tags['pop'],cellTags['x'], cellTags['y'], cellTags['z'], sim.rank)))
          sim.net.lastGid = sim.net.lastGid + self.tags['numCells'] 
          return cells