code: 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