docs: def calculate_matrices(states, Omega=1):
      r"""Calculate the matrices omega_ij, gamma_ij, r_pij.
  
      This function calculates the matrices omega_ij, gamma_ij and r_pij given a
      list of atomic states. The states can be arbitrarily in their fine,
      hyperfine or magnetic detail.
      """
      # We check that all states belong to the same element and the same isotope.
      iso = states[0].isotope
      element = states[0].element
      for state in states[1:]:
          if state.element != element:
              raise ValueError('All states must belong to the same element.')
          if state.isotope != iso:
              raise ValueError('All states must belong to the same isotope.')
  
      # We find the fine states involved in the problem.
      fine_states = find_fine_states(states)
  
      # We find the full magnetic states. The matrices will be first calculated
      # for the complete problem and later reduced to include only the states of
      # interest.
      full_magnetic_states = make_list_of_states(fine_states, 'magnetic',
                                                 verbose=0)
  
      # We calculate the indices corresponding to each sub matrix of fine and
      # hyperfine levels.
  
      # We calculate the frequency differences between states.
      omega_full = calculate_omega_matrix(full_magnetic_states, Omega)
      # We calculate the matrix gamma.
  
      gamma_full = calculate_gamma_matrix(full_magnetic_states, Omega)
      # We calculate the reduced matrix elements
      reduced_matrix_elements = calculate_reduced_matrix_elements(fine_states)
  
      # We calculate the matrices r_-1, r_0, r_1
      r_full = calculate_r_matrices(fine_states, reduced_matrix_elements)
  
      # Reduction to be implemented
      omega = omega_full
      r = r_full
      gamma = gamma_full
      return omega, gamma, r