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