<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def frontogenesis(thta, u, v, dx, dy, dim_order='yx'):
 r"""Calculate the 2D kinematic frontogenesis of a temperature field. The implementation is a form of the Petterssen Frontogenesis and uses the formula outlined in [Bluestein1993]_ pg.248-253. .. math:: F=\frac{1}{2}\left|\nabla \theta\right|[D cos(2\beta)-\delta] * :math:`F` is 2D kinematic frontogenesis * :math:`\theta` is potential temperature * :math:`D` is the total deformation * :math:`\beta` is the angle between the axis of dilitation and the isentropes * :math:`\delta` is the divergence Parameters thta : (M, N) ndarray Potential temperature u : (M, N) ndarray x component of the wind v : (M, N) ndarray y component of the wind dx : float or ndarray The grid spacing(s) in the x-direction. If an array, there should be one item less than the size of `u` along the applicable axis. dy : float or ndarray The grid spacing(s) in the y-direction. If an array, there should be one item less than the size of `u` along the applicable axis. Returns ------- (M, N) ndarray 2D Frontogenesis in [temperature units]/m/s Notes ----- If inputs have more than two dimensions, they are assumed to have either leading dimensions of (x, y) or trailing dimensions of (y, x), depending on the value of ``dim_order``. Conversion factor to go from [temperature units]/m/s to [temperature units/100km/3h] :math:`1.08e4*1.e5` """

    # Get gradients of potential temperature in both x and y
    ddy_thta = first_derivative(thta, delta=dy, axis=-2)
    ddx_thta = first_derivative(thta, delta=dx, axis=-1)

    # Compute the magnitude of the potential temperature gradient
    mag_thta = np.sqrt(ddx_thta**2 + ddy_thta**2)

    # Get the shearing, stretching, and total deformation of the wind field
    shrd = shearing_deformation(u, v, dx, dy, dim_order=dim_order)
    strd = stretching_deformation(u, v, dx, dy, dim_order=dim_order)
    tdef = total_deformation(u, v, dx, dy, dim_order=dim_order)

    # Get the divergence of the wind field
    div = divergence(u, v, dx, dy, dim_order=dim_order)

    # Compute the angle (beta) between the wind field and the gradient of potential temperature
    psi = 0.5 * np.arctan2(shrd, strd)
    beta = np.arcsin((-ddx_thta * np.cos(psi) - ddy_thta * np.sin(psi)) / mag_thta)

    return 0.5 * mag_thta * (tdef * np.cos(2 * beta) - div)