docs: def bdplotter(task, **kwargs):
      '''
      bdplotter is based on the fact that phonon DOS/bands and
      electron DOS/bands are the objects of the same kind.
      1) DOS is formatted precomputed / smeared according to a normal distribution
      2) bands are formatted precomputed / interpolated through natural cubic spline function
      '''
      if task == 'bands': # CRYSTAL, "VASP", EXCITING
  
          results = []
  
          if 'precomputed' in kwargs:
              for n in range(len(kwargs['precomputed']['ticks'])):
                  if kwargs['precomputed']['ticks'][n][1] == 'GAMMA': kwargs['precomputed']['ticks'][n][1] = '&#915;'
  
              for stripe in kwargs['precomputed']['stripes']:
                  results.append({'color':'#000000', 'data':[], 'ticks':kwargs['precomputed']['ticks']})
                  for n, val in enumerate(stripe):
                      results[-1]['data'].append([ kwargs['precomputed']['abscissa'][n], val])
  
          else:
              if not 'order' in kwargs: order = sorted( kwargs['values'].keys() ) # TODO
              else: order = kwargs['order']
  
              nullstand = '0 0 0'
              if not '0 0 0' in kwargs['values']: # possible case when there is shifted Gamma point
                  nullstand = order[0]
  
              # reduce k if too much
              if len(order)>20:
                  red_order = []
                  for i in range(0, len(order), int(math.floor(len(order)/10))):
                      red_order.append(order[i])
                  order = red_order
  
              for N in range(len( kwargs['values'][nullstand] )):
                  # interpolate for each curve throughout the BZ
                  results.append({'color':'#000000', 'data':[], 'ticks':[]})
                  d = 0.0
                  x = []
                  y = []
                  bz_vec_ref = [0, 0, 0]
  
                  for bz in order:
                      y.append( kwargs['values'][bz][N] )
                      bz_coords = map(frac2float, bz.split() )
                      bz_vec_cur = dot( bz_coords, linalg.inv( kwargs['xyz_matrix'] ).transpose() )
                      bz_vec_dir = map(sum, zip(bz_vec_cur, bz_vec_ref))
                      bz_vec_ref = bz_vec_cur
                      d += linalg.norm( bz_vec_dir )
                      x.append(d)
                      results[-1]['ticks'].append( [d, bz.replace(' ', '')] )
  
                  # end in nullstand point (normally, Gamma)
                  #y.append(kwargs['values'][nullstand][N])
                  #if d == 0: d+=0.5
                  #else: d += linalg.norm( bz_vec_ref )
                  #x.append(d)
                  #results[-1]['ticks'].append( [d, nullstand.replace(' ', '')] )
  
                  divider = 10 if len(order)<10 else 1.5
                  step = (max(x)-min(x)) / len(kwargs['values']) / divider
  
                  xnew = arange(min(x), max(x)+step/2, step).tolist()
                  ynew = []
                  f = NaturalCubicSpline( array(x), array(y) )
                  for i in xnew:
                      results[-1]['data'].append([ round(  i, 3  ), round(  f(i), 3  ) ]) # round to reduce output
  
          return results
  
      elif task == 'dos': # CRYSTAL, VASP, EXCITING
  
          results = []
  
          if 'precomputed' in kwargs:
              total_dos = [[i, kwargs['precomputed']['total'][n]] for n, i in enumerate(kwargs['precomputed']['x'])]
  
          else:
              tdos = TotalDos( kwargs['eigenvalues'], sigma=kwargs['sigma'] )
              tdos.set_draw_area(omega_min=kwargs['omega_min'], omega_max=kwargs['omega_max'], omega_pitch=kwargs['omega_pitch'])
              total_dos = tdos.calculate()
  
          results.append({'label':'total', 'color': '#000000', 'data': total_dos})
  
          if 'precomputed' in kwargs:
              partial_doses = []
              for k in kwargs['precomputed'].keys():
                  if k in ['x', 'total']: continue
                  partial_doses.append({ 'label': k, 'data': [[i, kwargs['precomputed'][k][n]] for n, i in enumerate(kwargs['precomputed']['x'])] })
  
          elif 'impacts' in kwargs and 'atomtypes' in kwargs:
              # get the order of atoms to evaluate their partial impact
              labels = {}
              types = []
              index, subtractor = 0, 0
              for k, atom in enumerate(kwargs['atomtypes']): # determine the order of atoms for the partial impact of every type
                  if atom not in labels:
                      #if atom == 'X' and not calc.phonons: # artificial GHOST case for phonons, decrease atomic index
                      #    subtractor += 1
                      #    continue
                      labels[atom] = index
                      types.append([k+1-subtractor])
                      index += 1
                  else:
                      types[ labels[atom] ].append(k+1-subtractor)
  
              pdos = PartialDos( kwargs['eigenvalues'], kwargs['impacts'], sigma=kwargs['sigma'] )
              pdos.set_draw_area(omega_min=kwargs['omega_min'], omega_max=kwargs['omega_max'], omega_pitch=kwargs['omega_pitch'])
              partial_doses = pdos.calculate( types, labels )
  
              # add colors to partials
              for i in range(len(partial_doses)):
                  if partial_doses[i]['label'] == 'X': color = '#000000'
                  elif partial_doses[i]['label'] == 'H': color = '#CCCCCC'
                  else:
                      try: color = jmol_to_hex( jmol_colors[ chemical_symbols.index(partial_doses[i]['label']) ] )
                      except ValueError: color = '#FFCC66'
                  partial_doses[i].update({'color': color})
              results.extend(partial_doses)
  
          return results