<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def euler_matrix(phi, theta=None, psi=None):
 """Rotation matrix in 2 and 3 dimensions. Its rows represent the canonical unit vectors as seen from the rotated system while the columns are the rotated unit vectors as seen from the canonical system. Parameters phi : float or `array-like` Either 2D counter-clockwise rotation angle (in radians) or first Euler angle. theta, psi : float or `array-like`, optional Second and third Euler angles in radians. If both are ``None``, a 2D rotation matrix is computed. Otherwise a 3D rotation is computed, where the default ``None`` is equivalent to ``0.0``. The rotation is performed in "ZXZ" rotation order, see the Wikipedia article `Euler angles`_. Returns ------- mat : `numpy.ndarray` Rotation matrix corresponding to the given angles. The returned array has shape ``(ndim, ndim)`` if all angles represent single parameters, with ``ndim == 2`` for ``phi`` only and ``ndim == 3`` for 2 or 3 Euler angles. If any of the angle parameters is an array, the shape of the returned array is ``broadcast(phi, theta, psi).shape + (ndim, ndim)``. References .. _Euler angles: https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix """

    if theta is None and psi is None:
        squeeze_out = (np.shape(phi) == ())
        ndim = 2
        phi = np.array(phi, dtype=float, copy=False, ndmin=1)
        theta = psi = 0.0
    else:
        # `None` broadcasts like a scalar
        squeeze_out = (np.broadcast(phi, theta, psi).shape == ())
        ndim = 3
        phi = np.array(phi, dtype=float, copy=False, ndmin=1)
        if theta is None:
            theta = 0.0
        if psi is None:
            psi = 0.0
        theta = np.array(theta, dtype=float, copy=False, ndmin=1)
        psi = np.array(psi, dtype=float, copy=False, ndmin=1)
        ndim = 3

    cph = np.cos(phi)
    sph = np.sin(phi)
    cth = np.cos(theta)
    sth = np.sin(theta)
    cps = np.cos(psi)
    sps = np.sin(psi)

    if ndim == 2:
        mat = np.array([[cph, -sph],
                        [sph, cph]])
    else:
        mat = np.array([
            [cph * cps - sph * cth * sps,
             -cph * sps - sph * cth * cps,
             sph * sth],
            [sph * cps + cph * cth * sps,
             -sph * sps + cph * cth * cps,
             -cph * sth],
            [sth * sps + 0 * cph,
             sth * cps + 0 * cph,
             cth + 0 * (cph + cps)]])  # Make sure all components broadcast

    if squeeze_out:
        return mat.squeeze()
    else:
        # Move the `(ndim, ndim)` axes to the end
        extra_dims = len(np.broadcast(phi, theta, psi).shape)
        newaxes = list(range(2, 2 + extra_dims)) + [0, 1]
        return np.transpose(mat, newaxes)