<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def from_spherical_coords(theta_phi, phi=None):
 """Return the quaternion corresponding to these spherical coordinates Assumes the spherical coordinates correspond to the quaternion R via R = exp(phi*z/2) * exp(theta*y/2) The angles naturally must be in radians for this to make any sense. Note that this quaternion rotates `z` onto the point with the given spherical coordinates, but also rotates `x` and `y` onto the usual basis vectors (theta and phi, respectively) at that point. Parameters theta_phi: float or array of floats This argument may either contain an array with last dimension of size 2, where those two elements describe the (theta, phi) values in radians for each point; or it may contain just the theta values in radians, in which case the next argument must also be given. phi: None, float, or array of floats If this array is given, it must be able to broadcast against the first argument. Returns ------- R: quaternion array If the second argument is not given to this function, the shape will be the same as the input shape except for the last dimension, which will be removed. If the second argument is given, this output array will have the shape resulting from broadcasting the two input arrays against each other. """

    # Figure out the input angles from either type of input
    if phi is None:
        theta_phi = np.asarray(theta_phi, dtype=np.double)
        theta = theta_phi[..., 0]
        phi  = theta_phi[..., 1]
    else:
        theta = np.asarray(theta_phi, dtype=np.double)
        phi = np.asarray(phi, dtype=np.double)

    # Set up the output array
    R = np.empty(np.broadcast(theta, phi).shape + (4,), dtype=np.double)

    # Compute the actual values of the quaternion components
    R[..., 0] =  np.cos(phi/2)*np.cos(theta/2)  # scalar quaternion components
    R[..., 1] = -np.sin(phi/2)*np.sin(theta/2)  # x quaternion components
    R[..., 2] =  np.cos(phi/2)*np.sin(theta/2)  # y quaternion components
    R[..., 3] =  np.sin(phi/2)*np.cos(theta/2)  # z quaternion components

    return as_quat_array(R)