docs: def get_tilting(self, oplane):
          '''
          Main procedure
          '''
          surf_atom1, surf_atom2, surf_atom3, surf_atom4 = oplane
  
          # divide surface atoms into groups by distance between them
          compare = [surf_atom2, surf_atom3, surf_atom4]
          distance_map = []
  
          for i in range(0, 3):
              distance_map.append([ compare[i], self.virtual_atoms.get_distance(surf_atom1, compare[i]) ])
  
          distance_map = sorted(distance_map, key=lambda x: x[1])
  
          distance_map_keys = [i[0] for i in distance_map]
          surf_atom3 = distance_map_keys[2]
          surf_atom2 = distance_map_keys[1]
          surf_atom4 = distance_map_keys[0]
  
          if self.virtual_atoms[surf_atom1].z == self.virtual_atoms[surf_atom2].z and \
              self.virtual_atoms[surf_atom2].z == self.virtual_atoms[surf_atom3].z and \
              self.virtual_atoms[surf_atom3].z == self.virtual_atoms[surf_atom4].z:
              # this is done to prevent false zero tilting
              self.virtual_atoms[surf_atom1].z += 1E-10
              self.virtual_atoms[surf_atom2].z += 1E-10
              self.virtual_atoms[surf_atom3].z -= 1E-10
              self.virtual_atoms[surf_atom4].z -= 1E-10
  
          # new axes will be defined simply as vectors standing on 1 - 3 and 2 - 4 (they are moved to the point of origin)
          self.virtual_atoms.append(Atom('X', (self.virtual_atoms[surf_atom1].x - self.virtual_atoms[surf_atom3].x, self.virtual_atoms[surf_atom1].y - self.virtual_atoms[surf_atom3].y, self.virtual_atoms[surf_atom1].z - self.virtual_atoms[surf_atom3].z)))
          self.virtual_atoms.append(Atom('X', (self.virtual_atoms[surf_atom2].x - self.virtual_atoms[surf_atom4].x, self.virtual_atoms[surf_atom2].y - self.virtual_atoms[surf_atom4].y, self.virtual_atoms[surf_atom2].z - self.virtual_atoms[surf_atom4].z)))
          self.virtual_atoms.append(Atom('X', (0, 0, 0)))
  
          # redefine tilted axes
          surf_atom_first = len(self.virtual_atoms)-3
          surf_atom_second = len(self.virtual_atoms)-2
          center = len(self.virtual_atoms)-1
  
          # inverse arbitrary atom
          self.virtual_atoms.append(Atom('X', (-self.virtual_atoms[surf_atom_first].x, -self.virtual_atoms[surf_atom_first].y, -self.virtual_atoms[surf_atom_first].z)))
          inversed_one = len(self.virtual_atoms)-1
  
          # find and add bisectors, silly swapping
          first_bisector = self.get_bisector_point(surf_atom_first, center, surf_atom_second, self.virtual_atoms)
          sec_bisector = self.get_bisector_point(surf_atom_second, center, inversed_one, self.virtual_atoms)
  
          swap = True
          if first_bisector[0] < 0 and sec_bisector[0] < 0:
              swap = False
          if first_bisector[0] < 0:
              first_bisector[0] *= -1
              first_bisector[1] *= -1
              first_bisector[2] *= -1
          if sec_bisector[0] < 0:
              sec_bisector[0] *= -1
              sec_bisector[1] *= -1
              sec_bisector[2] *= -1
          if swap:
              first_bisector, sec_bisector = sec_bisector, first_bisector
  
          swap = False
          if first_bisector[0] < sec_bisector[0] and first_bisector[1] < 0:
              first_bisector[0] *= -1
              first_bisector[1] *= -1
              first_bisector[2] *= -1
              swap = True
          if first_bisector[0] < sec_bisector[0] and first_bisector[1] > 0:
              swap = True
          if first_bisector[0] > sec_bisector[0] and sec_bisector[1] < 0:
              sec_bisector[0] *= -1
              sec_bisector[1] *= -1
              sec_bisector[2] *= -1
          if swap:
              first_bisector, sec_bisector = sec_bisector, first_bisector
  
          self.virtual_atoms.append(Atom('X', (first_bisector[0], first_bisector[1], first_bisector[2])))
          self.virtual_atoms.append(Atom('X', (sec_bisector[0], sec_bisector[1], sec_bisector[2])))
          first_bisector = len(self.virtual_atoms)-2
          sec_bisector = len(self.virtual_atoms)-1
  
          # use vector cross product to define normal which will play Z axis role
          self.virtual_atoms.append(Atom('X', (
              self.virtual_atoms[first_bisector].y*self.virtual_atoms[sec_bisector].z - self.virtual_atoms[first_bisector].z*self.virtual_atoms[sec_bisector].y,
              self.virtual_atoms[first_bisector].z*self.virtual_atoms[sec_bisector].x - self.virtual_atoms[first_bisector].x*self.virtual_atoms[sec_bisector].z,
              self.virtual_atoms[first_bisector].x*self.virtual_atoms[sec_bisector].y - self.virtual_atoms[first_bisector].y*self.virtual_atoms[sec_bisector].x
          )))
          tilt_z = len(self.virtual_atoms)-1
  
          # Euler angles ZYZ
          alpha = math.degrees(math.atan2(self.virtual_atoms[sec_bisector].z, self.virtual_atoms[first_bisector].z))
          beta = math.degrees(math.atan2(math.sqrt(self.virtual_atoms[tilt_z].x**2 + self.virtual_atoms[tilt_z].y**2), self.virtual_atoms[tilt_z].z))
          gamma = math.degrees(math.atan2(self.virtual_atoms[tilt_z].y, -self.virtual_atoms[tilt_z].x))
  
          # angles adjusting
          adjust_angles = [45, 90, 135, 180, 225, 270, 315, 360]
          tilting = [alpha, beta, gamma]
          for i in range(0, 3):
              tilting[i] = abs(tilting[i])
              if tilting[i] in adjust_angles:
                  tilting[i] = 0.0
                  continue
  
              if tilting[i] > self.MAX_TILTING_DEGREE:
                  for checkpoint in adjust_angles:
                      if checkpoint - self.MAX_TILTING_DEGREE < tilting[i] < checkpoint + self.MAX_TILTING_DEGREE:
                          tilting[i] = abs(tilting[i] - checkpoint)
                          break
          return tilting