subroutine cw_propagator(t0,x0,h,n,tf,xf,report)

    implicit none

    real(wp),intent(in)               :: t0      !! initialize time [sec]
    real(wp),dimension(6),intent(in)  :: x0      !! initial state in RST coordinates [km,km/s]
    real(wp),intent(in)               :: h       !! abs(time step) [sec]
    real(wp),intent(in)               :: n       !! mean motion of target orbit (`sqrt(mu/a**3)`) [1/sec]
    real(wp),intent(in)               :: tf      !! final time [sec]
    real(wp),dimension(6),intent(out) :: xf      !! final state in RST coordinates [km,km/s]
    procedure(report_func),optional   :: report  !! to report each point

    real(wp) :: t,dt,t2
    real(wp),dimension(6) :: x
    logical :: last,export

    export = present(report)

    if (export) call report(t0,x0)  !first point

    if (h==zero) then
        xf = x0
    else

        t = t0
        x = x0
        dt = sign(h,tf-t0)  !time step  (correct sign)
        do
            t2 = t + dt
            last = ((dt>=zero .and. t2>=tf) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tf))         !
            if (last) dt = tf-t                     !
            xf = cw_equations(x,dt,n)  ! propagate
            if (last) exit
            if (export) call report(t2,xf)   !intermediate point
            x = xf
            t = t2
        end do

    end if

    if (export) call report(tf,xf)   !last point

    end subroutine cw_propagator