<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def evolve_pn_spins(q, chiA0, chiB0, omega0, omegaTimesM_final, approximant='SpinTaylorT4', dt=0.1, spinO=7, phaseO=7):
 """ Evolves PN spins from a starting orbital frequency and spins to a final frequency. Inputs: q: Mass ratio (q>=1) chiA0: Dimless spin of BhA at initial freq. chiB0: Dimless spin of BhB at initial freq. omega0: Initial orbital frequency in dimless units. omegaTimesM_final: Final orbital frequency in dimless units. approximant: 'SpinTaylorT1/T2/T4'. Default: 'SpinTaylorT4'. dt: Dimless step time for evolution. Default: 0.1 . spinO: Twice PN order of spin effects. Default: 5 . phaseO: Twice PN order in phase. Default: 8 . Outputs (all are time series):
 chiA_end_copr: Spin of BhA at final frequency, in coprecessing frame. chiB_end_copr: Spin of BhB at final frequency, in coprecessing frame. q_copr_end: Coprecessing frame quaternion at final frequency. phi_end: Orbital phase in the coprecessing frame at final frequency. omegaTimesM_end Dimensionless final frequency. Should agree with omegaTimesM_final. The inertial frame is assumed to be aligned to the coorbital frame at orbital frequency = omega0. chiA0 and chiB0 are the inertial/coorbital frame spins at omega0. """

    omega, phi, chiA, chiB, lNhat, e1 = lal_spin_evloution_wrapper(approximant,
            q, omega0, chiA0, chiB0, dt, spinO, phaseO)

    # Compute omega, inertial spins, angular momentum direction and orbital
    # phase when omega = omegaTimesM_final
    end_idx = np.argmin(np.abs(omega - omegaTimesM_final))
    omegaTimesM_end = omega[end_idx]
    chiA_end = chiA[end_idx]
    chiB_end = chiB[end_idx]
    lNhat_end = lNhat[end_idx]
    phi_end = phi[end_idx]

    # Align the z-direction along orbital angular momentum direction
    # at end_idx. This moves us in to the coprecessing frame.
    q_copr_end = _utils.alignVec_quat(lNhat_end)
    chiA_end_copr = _utils.transformTimeDependentVector(
            np.array([q_copr_end]).T,
            np.array([chiA_end]).T, inverse=1).T[0]
    chiB_end_copr = _utils.transformTimeDependentVector(
            np.array([q_copr_end]).T,
            np.array([chiB_end]).T, inverse=1).T[0]

    return chiA_end_copr, chiB_end_copr, q_copr_end, phi_end, omegaTimesM_end