def estimate_dt_n_steps(w0, hamiltonian, n_periods, n_steps_per_period,
                        dE_threshold=1E-9, func=np.nanmax,
                        **integrate_kwargs):
    if not isinstance(w0, PhaseSpacePosition):
        w0 = np.asarray(w0)
        w0 = PhaseSpacePosition.from_w(w0, units=hamiltonian.units)
    dt = _autodetermine_initial_dt(w0, hamiltonian, dE_threshold=dE_threshold,
                                   **integrate_kwargs)
    n_steps = int(round(10000 / dt))
    orbit = hamiltonian.integrate_orbit(w0, dt=dt, n_steps=n_steps,
                                        **integrate_kwargs)
    circ = orbit.circulation()
    if np.any(circ):
        orbit = orbit.align_circulation_with_z(circulation=circ)
        cyl = orbit.represent_as(coord.CylindricalRepresentation)
        R = cyl.rho.value
        phi = cyl.phi.value
        z = cyl.z.value
        T = np.array([peak_to_peak_period(orbit.t, f).value
                      for f in [R, phi, z]])*orbit.t.unit
    else:
        T = np.array([peak_to_peak_period(orbit.t, f).value
                      for f in orbit.pos])*orbit.t.unit
    T = func(T)
    if np.isnan(T):
        raise RuntimeError("Failed to find period.")
    T = T.decompose(hamiltonian.units).value
    dt = T / float(n_steps_per_period)
    n_steps = int(round(n_periods * T / dt))
    if dt == 0. or dt < 1E-13:
        raise ValueError("Timestep is zero or very small!")
    return dt, n_steps