<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_spins(self, q, chiA0, chiB0, omega0, PN_approximant, PN_dt, PN_spin0, PN_phase0, omega0_nrsur):
 """ Evolves spins of the component BHs from an initial orbital frequency = omega0 until t=-100 M from the peak of the waveform. If omega0 < omega0_nrsur, use PN to evolve the spins until orbital frequency = omega0. Then evolves further with the NRSur7dq2 waveform model until t=-100M from the peak. Assumes chiA0 and chiB0 are defined in the inertial frame defined at orbital frequency = omega0 as: The z-axis is along the Newtonian orbital angular momentum when the PN orbital frequency = omega0. The x-axis is along the line of separation from the smaller BH to the larger BH at this frequency. The y-axis completes the triad. Returns spins in the coorbital frame at t=-100M, as well as the coprecessing frame quaternion and orbital phase in the coprecessing frame at this time. """


        if omega0 < omega0_nrsur:
            # If omega0 is below the NRSur7dq2 start frequency, we use PN
            # to evolve the spins until orbital frequency = omega0_nrsur.

            # Note that we update omega0_nrsur here with the PN
            # frequency that was closest to the input omega0_nrsur.
            chiA0_nrsur_copr, chiB0_nrsur_copr, quat0_nrsur_copr, \
                phi0_nrsur, omega0_nrsur \
                = evolve_pn_spins(q, chiA0, chiB0, omega0,
                    omega0_nrsur, approximant=PN_approximant,
                    dt=PN_dt, spinO=PN_spin0,
                    phaseO=PN_phase0)
        else:
            # If omega0>= omega0_nrsur, we evolve spins directly with NRSur7dq2
            # waveform model. We set the coprecessing frame quaternion to
            # identity and orbital phase to 0 at omega=omega0, hence the
            # coprecessing frame is the same as the inertial frame here.

            # Note that we update omega0_nrsur here and set it to omega0
            chiA0_nrsur_copr, chiB0_nrsur_copr, quat0_nrsur_copr, \
                phi0_nrsur, omega0_nrsur \
                = chiA0, chiB0, [1,0,0,0], 0, omega0

        # Load NRSur7dq2 if needed
        if self.nrsur is None:
            self._load_NRSur7dq2()

        # evaluate NRSur7dq2 dynamics
        # We set allow_extrapolation=True always since we test param limits
        # independently
        quat, orbphase, chiA_copr, chiB_copr = self.nrsur.get_dynamics(q,
            chiA0_nrsur_copr, chiB0_nrsur_copr, init_quat=quat0_nrsur_copr,
            init_phase=phi0_nrsur, omega_ref=omega0_nrsur,
            allow_extrapolation=True)

        # get data at time node where remnant fits are done
        fitnode_time = -100
        nodeIdx = np.argmin(np.abs(self.nrsur.tds - fitnode_time))
        quat_fitnode = quat.T[nodeIdx]
        orbphase_fitnode = orbphase[nodeIdx]

        # get coorbital frame spins at the time node
        chiA_coorb_fitnode = utils.rotate_in_plane(chiA_copr[nodeIdx],
                orbphase_fitnode)
        chiB_coorb_fitnode = utils.rotate_in_plane(chiB_copr[nodeIdx],
                orbphase_fitnode)

        return chiA_coorb_fitnode, chiB_coorb_fitnode, quat_fitnode, \
            orbphase_fitnode