code: def bode_phase_numerical_data(
      system,
      initial_exp=-5,
      final_exp=5,
      freq_unit="rad/sec",
      phase_unit="rad",
      phase_unwrap=True,
      **kwargs,
  ):
      """
      Returns the numerical data of the Bode phase plot of the system.
      It is internally used by ``bode_phase_plot`` to get the data
      for plotting Bode phase plot. Users can use this data to further
      analyse the dynamics of the system or plot using a different
      backend/plotting-module.
  
      Parameters
      ==========
  
      system : SISOLinearTimeInvariant
          The system for which the Bode phase plot data is to be computed.
      initial_exp : Number, optional
          The initial exponent of 10 of the semilog plot. Defaults to -5.
      final_exp : Number, optional
          The final exponent of 10 of the semilog plot. Defaults to 5.
      freq_unit : string, optional
          User can choose between ``'rad/sec'`` (radians/second) and '``'Hz'`` (Hertz) as frequency units.
      phase_unit : string, optional
          User can choose between ``'rad'`` (radians) and ``'deg'`` (degree) as phase units.
      phase_unwrap : bool, optional
          Set to ``True`` by default.
  
      Returns
      =======
  
      tuple : (x, y)
          x = x-axis values of the Bode phase plot.
          y = y-axis values of the Bode phase plot.
  
      Raises
      ======
  
      NotImplementedError
          When a SISO LTI system is not passed.
  
          When time delay terms are present in the system.
  
      ValueError
          When more than one free symbol is present in the system.
          The only variable in the transfer function should be
          the variable of the Laplace transform.
  
          When incorrect frequency or phase units are given as input.
  
      Examples
      ========
  
      >>> from sympy.abc import s
      >>> from sympy.physics.control.lti import TransferFunction
      >>> from sympy.physics.control.control_plots import bode_phase_numerical_data
      >>> tf1 = TransferFunction(s**2 + 1, s**4 + 4*s**3 + 6*s**2 + 5*s + 2, s)
      >>> bode_phase_numerical_data(tf1)   # doctest: +SKIP
      ([1e-05, 1.4472354033813751e-05, 2.035581932165858e-05,..., 47577.3248186011, 67884.09326036123, 100000.0],
      [-2.5000000000291665e-05, -3.6180885085e-05, -5.08895483066e-05,...,-3.1415085799262523, -3.14155265358979])
  
      See Also
      ========
  
      bode_magnitude_plot, bode_phase_numerical_data
  
      """
      _check_system(system)
      expr = system.to_expr()
      freq_units = ("rad/sec", "Hz")
      phase_units = ("rad", "deg")
      if freq_unit not in freq_units:
          raise ValueError(
              'Only "rad/sec" and "Hz" are accepted frequency units.'
          )
      if phase_unit not in phase_units:
          raise ValueError(
              'Only "rad" and "deg" are accepted phase units.'
          )
  
      _w = Dummy("w", real=True)
      if freq_unit == "Hz":
          repl = I * _w * 2 * pi
      else:
          repl = I * _w
      w_expr = expr.subs({system.var: repl})
  
      if phase_unit == "deg":
          phase = arg(w_expr) * 180 / pi
      else:
          phase = arg(w_expr)
  
      x, y = LineOver1DRangeSeries(
          phase,
          (_w, 10**initial_exp, 10**final_exp),
          xscale="log",
          **kwargs,
      ).get_points()
  
      half = None
      if phase_unwrap:
          if phase_unit == "rad":
              half = pi
          elif phase_unit == "deg":
              half = 180
      if half:
          unit = 2 * half
          for i in range(1, len(y)):
              diff = y[i] - y[i - 1]
              if diff > half:  # Jump from -half to half
                  y[i] = y[i] - unit
              elif diff < -half:  # Jump from half to -half
                  y[i] = y[i] + unit
  
      return x, y
  
imports: from sympy.core.numbers import I, pi
  from sympy.functions.elementary.exponential import exp, log
  from sympy.polys.partfrac import apart
  from sympy.core.symbol import Dummy
  from sympy.external import import_module
  from sympy.functions import arg, Abs
  from sympy.integrals.laplace import _fast_inverse_laplace
  from sympy.physics.control.lti import SISOLinearTimeInvariant
  from sympy.plotting.series import LineOver1DRangeSeries
  from sympy.polys.polytools import Poly
  from sympy.printing.latex import latex
label_desc: Use image node_id 12 for calling a global function with example usage: bode_phase_numerical_data(system, initial_exp, final_exp, freq_unit, phase_unit, phase_unwrap) and returns: x,  y