docs: def true_lam_csd(cell, dr=100, z=None):
      """
      Return CSD from membrane currents as function along the coordinates
      of the electrode along z-axis. 
  
  
      Parameters
      ----------
      cell : `LFPy.cell.Cell` or `LFPy.cell.TemplateCell` object.
          Cell.
      dr : float
          Radius of the cylindrical CSD volume.
      z : numpy.ndarray
          Z-coordinates of electrode.
  
  
      Returns
      ----------
      CSD : numpy.ndarray
          Current-source density (in pA * mum^-3).
      
      """
      if type(z) != type(np.ndarray(shape=0)):
          raise ValueError('type(z) should be a numpy.ndarray')
  
      dz = abs(z[1] - z[0])
      CSD = np.zeros((z.size, cell.tvec.size,))
      r_end = np.sqrt(cell.xend**2 + cell.yend**2)
      r_start = np.sqrt(cell.xstart**2 + cell.ystart**2)
      volume = dz * np.pi * dr**2
  
      for i in range(len(z)):
          aa0 = cell.zstart < z[i] + dz/2
          aa1 = cell.zend < z[i] + dz/2
          bb0 = cell.zstart >= z[i] - dz/2
          bb1 = cell.zend >= z[i] - dz/2
          cc0 = r_start < dr
          cc1 = r_end < dr
          ii = aa0 & bb0 & cc0 # startpoint inside volume
          jj = aa1 & bb1 & cc1 # endpoint inside volume
  
          for j in range(cell.zstart.size):
              # Calc fraction of source being inside control volume from 0-1
              # both start and endpoint in volume
              if ii[j] and jj[j]:
                  CSD[i,] = CSD[i, ] + cell.imem[j, ] / volume
              # Startpoint in volume:
              elif ii[j] and not jj[j]: 
                  z0 = cell.zstart[j]
                  r0 = r_start[j]
                  
                  z1 = cell.zend[j]
                  r1 = r_end[j]
  
                  if r0 == r1:
                      # Segment is parallel with z-axis
                      frac = -(z0 - z[i]-dz/2) / (z1 - z0)
                  else:
                      # Not parallel with z-axis                    
                      L2 = (r1 - r0)**2 + (z1 - z0)**2
      
                      z2 = [z[i]-dz/2, z[i]+dz/2, z[i]-dz/2]
                      r2 = [0, 0, dr]
                      
                      z3 = [z[i]-dz/2, z[i]+dz/2, z[i]+dz/2]
                      r3 = [dr, dr, dr]
      
                      P = []
                      for k in range(3):
                          P.append(_PrPz(r0, z0, r1, z1, r2[k], z2[k], r3[k], z3[k]))
                          if P[k][2]:
                              vL2 = (P[k][0] - r0)**2 + (P[k][1] -z0)**2
                              frac = np.sqrt(vL2 / L2)
                  CSD[i,] = CSD[i, ] + frac * cell.imem[j, ] / volume
              # Endpoint in volume:
              elif jj[j] and not ii[j]: 
                  z0 = cell.zstart[j]
                  r0 = r_start[j]
                  
                  z1 = cell.zend[j]
                  r1 = r_end[j]
                  
                  if r0 == r1:
                      # Segment is parallel with z-axis
                      frac = (z1 - z[i]+dz/2) / (z1 - z0)
                  else:
                      # Not parallel with z-axis
                      L2 = (r1 - r0)**2 + (z1 - z0)**2
      
                      z2 = [z[i]-dz/2, z[i]+dz/2, z[i]-dz/2]
                      r2 = [0, 0, dr]
                      
                      z3 = [z[i]-dz/2, z[i]+dz/2, z[i]+dz/2]
                      r3 = [dr, dr, dr]
      
                      P = []
                      for k in range(3):
                          P.append(_PrPz(r0, z0, r1, z1, r2[k], z2[k], r3[k], z3[k]))
                          if P[k][2]:
                              vL2 = (r1 - P[k][0])**2 + (z1 - P[k][1])**2
                              frac = np.sqrt(vL2 / L2)
                  CSD[i,] = CSD[i, ] + frac * cell.imem[j, ] / volume
              else:
                  pass
  
      return CSD