Conic section parameter of the Kepler problem.
+Conic section semi-latus rectum of the Kepler problem.
Alternative initialiser from system and geometry specifications.
-Given the eccentricity \(e\) and the conic section parameter \(p\), +
Given the eccentricity \(e\) and the conic section semi-latus rectum \(p\), \(\(l = \pm \sqrt{mp}\,,\quad E = (e^2-1) \left|E_\text{min}\right|\,.\)\)
@@ -2331,7 +2331,7 @@geometric specifications positive_angular_mom: whether the angular momentum is positive ecc: eccentricity of the conic section parameter: parameter of the conic section
geometric specifications positive_angular_mom: whether the angular momentum is positive ecc: eccentricity of the conic section parameter: semi-latus rectum of the conic section
hamilflow","text":"Generating dataset for physical systems.
"},{"location":"#development","title":"Development","text":"We use poetry to manage our python environment. Use
poetry install\nto install the requirements. Or run
poetry install --with test\nto install the base environment and test environment for development.
If this is the first you clone the repo and committing code, run
pre-commit install\nfirst.
"},{"location":"changelog/","title":"hamilflow Changelog","text":""},{"location":"changelog/#010-2024-03-21","title":"0.1.0 (2024-03-21)","text":""},{"location":"changelog/#feat","title":"Feat","text":"Setting up repository.
"},{"location":"references/","title":"References","text":""},{"location":"references/maths/trigonometrics/","title":"Trigonometrics","text":"Trigonometric functions.
"},{"location":"references/maths/trigonometrics/#hamilflow.maths.trigonometrics.acos_with_shift","title":"acos_with_shift(x, shift=None)","text":"Arccos with shift.
Source code inhamilflow/maths/trigonometrics.py def acos_with_shift(\n x: \"Collection[float] | npt.ArrayLike\",\n shift: \"Collection[float] | npt.ArrayLike | None\" = None,\n) -> \"npt.ArrayLike\":\n \"\"\"Arccos with shift.\"\"\"\n x = np.asarray(x)\n value = np.arccos(x)\n shift = np.asarray(shift)\n period_shift = (div := np.floor(shift)) * 2 * np.pi\n remainder = shift - div\n value = np.where(remainder <= 0.5, value, 2 * np.pi - value) # noqa: PLR2004\n\n return period_shift + value\nMain module for Brownian motion.
"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion","title":"BrownianMotion","text":"Brownian motion describes motion of small particles with stochastic forces applied to them.
The math of Brownian motion can be modeled with Wiener process.
For consistency, we always use \\(\\mathbf x\\) for displacement, and \\(t\\) for steps. The model we are using is
\\[ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} \\]"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion--references","title":"References","text":"1D Brownian Motion
The dimsion of our Brownian motion is specified by the dimension of the initial condition.
To simulate a 1D Browian motion, we define the system and initial condition:
system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n}\n\ninitial_condition = {\n \"x0\": 0\n}\nThe Brownian motion can be simulated using
bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\nbm(n_steps=100)\n2D Brownian Motion
To simulate a 2D Browian motion,
system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n}\n\ninitial_condition = {\n \"x0\": [0, 0]\n}\n\nbm = BrownianMotion(system=system, initial_condition=initial_condition)\n\nbm(n_steps=100)\nParameters:
Name Type Description Defaultsystem Mapping[str, float] the Brownian motion system definition
requiredinitial_condition Mapping[str, Sequence[float] | ArrayLike] | None the initial condition for the simulation
None Source code in hamilflow/models/brownian_motion.py class BrownianMotion:\n r\"\"\"Brownian motion describes motion of small particles with stochastic forces applied to them.\n\n The math of Brownian motion can be modeled\n with Wiener process.\n\n For consistency, we always use\n $\\mathbf x$ for displacement, and\n $t$ for steps. The model we are using is\n\n $$\n \\begin{align}\n \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) +\n \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t})\n \\end{align}\n $$\n\n References\n ----------\n 1. Brownian motion and random walks. [cited 13 Mar 2024].\n Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html\n 2. Contributors to Wikimedia projects. Brownian motion.\n In: Wikipedia [Internet]. 22 Jan 2024 [cited 13 Mar 2024].\n Available: https://en.wikipedia.org/wiki/Brownian_motion\n\n\n !!! example \"1D Brownian Motion\"\n\n The dimsion of our Brownian motion is specified by\n the dimension of the initial condition.\n\n To simulate a 1D Browian motion, we define the system and initial condition:\n\n ```python\n system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n }\n\n initial_condition = {\n \"x0\": 0\n }\n ```\n\n The Brownian motion can be simulated using\n\n ```python\n bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\n bm(n_steps=100)\n ```\n\n !!! example \"2D Brownian Motion\"\n\n To simulate a 2D Browian motion,\n\n ```python\n system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n }\n\n initial_condition = {\n \"x0\": [0, 0]\n }\n\n bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\n bm(n_steps=100)\n ```\n\n :param system: the Brownian motion system definition\n :param initial_condition: the initial condition for the simulation\n\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: (\n Mapping[str, \"Sequence[float] | npt.ArrayLike\"] | None\n ) = None,\n ) -> None:\n initial_condition = initial_condition or {}\n self.system = BrownianMotionSystem.model_validate(system)\n self.initial_condition = BrownianMotionIC.model_validate(initial_condition)\n\n @property\n def dim(self) -> int:\n \"\"\"Dimension of the Brownian motion.\"\"\"\n return np.asarray(self.initial_condition.x0).size\n\n @property\n def _axis_names(self) -> list[str]:\n return [f\"x_{i}\" for i in range(self.dim)]\n\n def _trajectory(self, n_new_steps: int, seed: int) -> \"npt.NDArray[np.float64]\":\n \"\"\"Give the trajectory of the particle.\n\n We first compute the delta displacement in each step.\n With the displacement at each step, we perform a cumsum\n including the initial coordinate to get the displacement at each step.\n\n :param n_new_steps: number of new steps to simulate, excluding the initial step.\n :param seed: seed for the random generator.\n \"\"\"\n step_history = sp.stats.norm.rvs(\n size=(n_new_steps, self.dim) if self.dim > 1 else n_new_steps,\n scale=self.system.gaussian_scale,\n random_state=np.random.RandomState(seed=seed),\n )\n\n step_history = np.concatenate(\n (np.expand_dims(self.initial_condition.x0, axis=0), step_history),\n )\n\n trajectory = np.cumsum(step_history, axis=0)\n\n return trajectory\n\n def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:\n \"\"\"Generate data from a set of interpretable params for this model.\n\n :param n_steps: total number of steps to be simulated, including the inital step.\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n time_steps = np.arange(0, n_steps) * self.system.delta_t\n\n return self(t=time_steps, seed=seed)\n\n def __call__(self, t: TypeTime, seed: int = 42) -> pd.DataFrame:\n \"\"\"Simulate the coordinates of the particle.\n\n :param t: the time sequence to be used to generate data, 1-D array like\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n n_steps = np.array(t).size\n trajectory = self._trajectory(n_new_steps=n_steps - 1, seed=seed)\n\n df = pd.DataFrame(trajectory, columns=self._axis_names)\n df[\"t\"] = t\n\n return df\ndim: int property","text":"Dimension of the Brownian motion.
"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion.__call__","title":"__call__(t, seed=42)","text":"Simulate the coordinates of the particle.
Parameters:
Name Type Description Defaultt TypeTime the time sequence to be used to generate data, 1-D array like
requiredseed int random generator seed for the stochastic process. Use it to reproduce results.
42 Source code in hamilflow/models/brownian_motion.py def __call__(self, t: TypeTime, seed: int = 42) -> pd.DataFrame:\n \"\"\"Simulate the coordinates of the particle.\n\n :param t: the time sequence to be used to generate data, 1-D array like\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n n_steps = np.array(t).size\n trajectory = self._trajectory(n_new_steps=n_steps - 1, seed=seed)\n\n df = pd.DataFrame(trajectory, columns=self._axis_names)\n df[\"t\"] = t\n\n return df\ngenerate_from(n_steps, seed=42)","text":"Generate data from a set of interpretable params for this model.
Parameters:
Name Type Description Defaultn_steps int total number of steps to be simulated, including the inital step.
requiredseed int random generator seed for the stochastic process. Use it to reproduce results.
42 Source code in hamilflow/models/brownian_motion.py def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:\n \"\"\"Generate data from a set of interpretable params for this model.\n\n :param n_steps: total number of steps to be simulated, including the inital step.\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n time_steps = np.arange(0, n_steps) * self.system.delta_t\n\n return self(t=time_steps, seed=seed)\nBrownianMotionIC","text":" Bases: BaseModel
The initial condition for a Brownian motion.
Attributes:
Name Type Descriptionx0 float | Sequence[float] initial displacement of the particle, the diminsion of this initial condition determines the dimension of the model too.
Source code inhamilflow/models/brownian_motion.py class BrownianMotionIC(BaseModel):\n \"\"\"The initial condition for a Brownian motion.\n\n :cvar x0: initial displacement of the particle,\n the diminsion of this initial condition determines\n the dimension of the model too.\n \"\"\"\n\n x0: float | Sequence[float] = Field(default=1.0)\n\n @field_validator(\"x0\")\n @classmethod\n def _check_x0_types(cls, v: float | Sequence[float]) -> float | Sequence[float]:\n if not isinstance(v, float | int | Sequence):\n # TODO I do not think this raise can be reached\n msg = f\"Value of x0 should be int/float/list of int/float: {v=}\"\n raise TypeError(msg)\n\n return v\nBrownianMotionSystem","text":" Bases: BaseModel
Definition of the Brownian Motion system.
For consistency, we always use \\(\\mathbf x\\) for displacement, and \\(t\\) for steps. The model we are using is
\\[ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} \\]"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotionSystem--references","title":"References","text":"Attributes:
Name Type Descriptionsigma float base standard deviation to be used to compute the variance
delta_t float time granunality of the motion
Source code inhamilflow/models/brownian_motion.py class BrownianMotionSystem(BaseModel):\n r\"\"\"Definition of the Brownian Motion system.\n\n For consistency, we always use\n $\\mathbf x$ for displacement, and\n $t$ for steps. The model we are using is\n\n $$\n \\begin{align}\n \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) +\n \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t})\n \\end{align}\n $$\n\n References\n ----------\n 1. Brownian motion and random walks. [cited 13 Mar 2024].\n Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html\n 2. Contributors to Wikimedia projects. Brownian motion.\n In: Wikipedia [Internet]. 22 Jan 2024 [cited 13 Mar 2024].\n Available: https://en.wikipedia.org/wiki/Brownian_motion\n\n :cvar sigma: base standard deviation\n to be used to compute the variance\n :cvar delta_t: time granunality of the motion\n\n \"\"\"\n\n sigma: float = Field(ge=0.0)\n delta_t: float = Field(ge=0.0, default=1.0)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def gaussian_scale(self) -> float:\n \"\"\"The scale (standard deviation) of the Gaussian term in Brownian motion.\"\"\"\n return self.sigma**2 * self.delta_t\ngaussian_scale: float cached property","text":"The scale (standard deviation) of the Gaussian term in Brownian motion.
"},{"location":"references/models/harmonic_oscillator/","title":"Harmonic Oscillator","text":"Main module for undamped and damped hamornic oscillators.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.ComplexSimpleHarmonicOscillator","title":"ComplexSimpleHarmonicOscillator","text":"Generate time series data for a complex simple harmonic oscillator.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the complex harmonic oscillator.
requiredinitial_condition Mapping[str, tuple[float, float]] the initial condition of the complex harmonic oscillator.
required Source code inhamilflow/models/harmonic_oscillator.py class ComplexSimpleHarmonicOscillator:\n r\"\"\"Generate time series data for a complex simple harmonic oscillator.\n\n :param system: all the params that defines the complex harmonic oscillator.\n :param initial_condition: the initial condition of the complex harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, tuple[float, float]],\n ) -> None:\n self.system = HarmonicOscillatorSystem.model_validate(system)\n self.initial_condition = ComplexSimpleHarmonicOscillatorIC.model_validate(\n initial_condition,\n )\n if self.system.type != \"simple\":\n msg = f\"System is not a Simple Harmonic Oscillator: {self.system}\"\n raise ValueError(\n msg,\n )\n\n @cached_property\n def definition(\n self,\n ) -> dict[str, dict[str, float | tuple[float, float]]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n def _z(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to complex simple harmonic oscillators.\n\n $$\n x(t) = x_+ \\exp(-\\mathbb{i} (\\omega t + \\phi_+)) + x_- \\exp(+\\mathbb{i} (\\omega t + \\phi_-)).\n $$\n \"\"\"\n t = np.asarray(t)\n omega = self.system.omega\n x0, phi = self.initial_condition.x0, self.initial_condition.phi\n phases = -omega * t - phi[0], omega * t + phi[1]\n return x0[0] * np.exp(1j * phases[0]) + x0[1] * np.exp(1j * phases[1])\n\n def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time.\n\n :param t: time(s).\n \"\"\"\n t = np.asarray(t)\n data = self._z(t)\n\n return pd.DataFrame({\"t\": t, \"z\": data})\ndefinition: dict[str, dict[str, float | tuple[float, float]]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.ComplexSimpleHarmonicOscillator.__call__","title":"__call__(t)","text":"Generate time series data for the harmonic oscillator.
Returns a list of floats representing the displacement at each time.
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time(s).
required Source code inhamilflow/models/harmonic_oscillator.py def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time.\n\n :param t: time(s).\n \"\"\"\n t = np.asarray(t)\n data = self._z(t)\n\n return pd.DataFrame({\"t\": t, \"z\": data})\nComplexSimpleHarmonicOscillatorIC","text":" Bases: BaseModel
The initial condition for a complex harmonic oscillator.
Attributes:
Name Type Descriptionx0 tuple[float, float] the initial displacements
phi tuple[float, float] initial phases
Source code inhamilflow/models/harmonic_oscillator.py class ComplexSimpleHarmonicOscillatorIC(BaseModel):\n \"\"\"The initial condition for a complex harmonic oscillator.\n\n :cvar x0: the initial displacements\n :cvar phi: initial phases\n \"\"\"\n\n x0: tuple[float, float] = Field()\n phi: tuple[float, float] = Field(default=(0, 0))\nDampedHarmonicOscillator","text":" Bases: HarmonicOscillatorBase
Generate time series data for a damped harmonic oscillator.
The equation for a general un-driven harmonic oscillator is12
\\[ \\frac{\\mathrm d x^2}{\\mathrm d t^2} + 2\\zeta \\omega \\frac{\\mathrm d x}{\\mathrm dt} + \\omega^2 x = 0, \\]where \\(x\\) is the displacement, \\(\\omega\\) is the angular frequency of an undamped oscillator (\\(\\zeta=0\\)), and \\(\\zeta\\) is the damping ratio.
The solution to the above harmonic oscillator is
\\[ x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right) e^{-\\zeta \\omega t}, \\]where
\\[ \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}. \\]To use this generator,
params = {\"omega\": omega, \"zeta\"=0.2}\n\nho = DampedHarmonicOscillator(params=params)\n\ndf = ho(n_periods=1, n_samples_per_period=10)\ndf will be a pandas dataframe with two columns: t and x.
Contributors to Wikimedia projects. Harmonic oscillator. In: Wikipedia [Internet]. 18 Feb 2024 [cited 20 Feb 2024]. Available: https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator \u21a9
Libretexts. 5.3: General Solution for the Damped Harmonic Oscillator. Libretexts. 13 Apr 2021. Available: https://t.ly/cWTIo. Accessed 20 Feb 2024.\u00a0\u21a9
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the harmonic oscillator.
requiredinitial_condition Mapping[str, float] | None the initial condition of the harmonic oscillator.
None Source code in hamilflow/models/harmonic_oscillator.py class DampedHarmonicOscillator(HarmonicOscillatorBase):\n r\"\"\"Generate time series data for a [damped harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n The equation for a general un-driven harmonic oscillator is[^wiki_ho][^libretext_ho]\n\n $$\n \\frac{\\mathrm d x^2}{\\mathrm d t^2} + 2\\zeta \\omega \\frac{\\mathrm d x}{\\mathrm dt} + \\omega^2 x = 0,\n $$\n\n where $x$ is the displacement, $\\omega$ is the angular frequency of an undamped oscillator ($\\zeta=0$),\n and $\\zeta$ is the damping ratio.\n\n [^wiki_ho]: Contributors to Wikimedia projects. Harmonic oscillator.\n In: Wikipedia [Internet]. 18 Feb 2024 [cited 20 Feb 2024].\n Available: https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator\n\n [^libretext_ho]: Libretexts. 5.3: General Solution for the Damped Harmonic Oscillator. Libretexts. 13 Apr 2021.\n Available: https://t.ly/cWTIo. Accessed 20 Feb 2024.\n\n\n The solution to the above harmonic oscillator is\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n\n To use this generator,\n\n ```python\n params = {\"omega\": omega, \"zeta\"=0.2}\n\n ho = DampedHarmonicOscillator(params=params)\n\n df = ho(n_periods=1, n_samples_per_period=10)\n ```\n\n `df` will be a pandas dataframe with two columns: `t` and `x`.\n\n :param system: all the params that defines the harmonic oscillator.\n :param initial_condition: the initial condition of the harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n super().__init__(system, initial_condition)\n if self.system.type == \"simple\":\n msg = (\n f\"System is not a Damped Harmonic Oscillator: {self.system}\\n\"\n f\"This is a simple harmonic oscillator, use `SimpleHarmonicOscillator`.\"\n )\n raise ValueError(\n msg,\n )\n\n def _x_under_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to under damped harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n \"\"\"\n t = np.asarray(t)\n omega_damp = self.system.omega * np.sqrt(1 - self.system.zeta)\n return (\n self.initial_condition.x0 * np.cos(omega_damp * t)\n + (\n self.system.zeta * self.system.omega * self.initial_condition.x0\n + self.initial_condition.v0\n )\n / omega_damp\n * np.sin(omega_damp * t)\n ) * np.exp(-self.system.zeta * self.system.omega * t)\n\n def _x_critical_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to critical damped harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n \"\"\"\n t = np.asarray(t)\n return self.initial_condition.x0 * np.exp(\n -self.system.zeta * self.system.omega * t,\n )\n\n def _x_over_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to over harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cosh(\\Gamma t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Gamma} \\sinh(\\Gamma t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Gamma = \\omega\\sqrt{ \\zeta^2 - 1 }.\n $$\n \"\"\"\n t = np.asarray(t)\n gamma_damp = self.system.omega * np.sqrt(self.system.zeta - 1)\n\n return (\n self.initial_condition.x0 * np.cosh(gamma_damp * t)\n + (\n self.system.zeta * self.system.omega * self.initial_condition.x0\n + self.initial_condition.v0\n )\n / gamma_damp\n * np.sinh(gamma_damp * t)\n ) * np.exp(-self.system.zeta * self.system.omega * t)\n\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to damped harmonic oscillators.\"\"\"\n t = np.asarray(t)\n if self.system.type == \"under_damped\":\n x = self._x_under_damped(t)\n elif self.system.type == \"over_damped\":\n x = self._x_over_damped(t)\n elif self.system.type == \"critical_damped\":\n x = self._x_critical_damped(t)\n else:\n msg = f\"System type is not damped harmonic oscillator: {self.system.type}\"\n raise ValueError(\n msg,\n )\n\n return x\nHarmonicOscillatorBase","text":" Bases: ABC
Base class to generate time series data for a harmonic oscillator.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the harmonic oscillator.
requiredinitial_condition Mapping[str, float] | None the initial condition of the harmonic oscillator.
None Source code in hamilflow/models/harmonic_oscillator.py class HarmonicOscillatorBase(ABC):\n r\"\"\"Base class to generate time series data for a [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n :param system: all the params that defines the harmonic oscillator.\n :param initial_condition: the initial condition of the harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n initial_condition = initial_condition or {}\n self.system = HarmonicOscillatorSystem.model_validate(system)\n self.initial_condition = HarmonicOscillatorIC.model_validate(initial_condition)\n\n @cached_property\n def definition(self) -> dict[str, dict[str, float]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n @abstractmethod\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to simple harmonic oscillators.\"\"\"\n ...\n\n def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time step.\n\n :param n_periods: Number of periods to generate.\n :param n_samples_per_period: Number of samples per period.\n \"\"\"\n time_delta = self.system.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n data = self._x(time_steps)\n\n return pd.DataFrame({\"t\": time_steps, \"x\": data})\ndefinition: dict[str, dict[str, float]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorBase.__call__","title":"__call__(n_periods, n_samples_per_period)","text":"Generate time series data for the harmonic oscillator.
Returns a list of floats representing the displacement at each time step.
Parameters:
Name Type Description Defaultn_periods int Number of periods to generate.
requiredn_samples_per_period int Number of samples per period.
required Source code inhamilflow/models/harmonic_oscillator.py def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time step.\n\n :param n_periods: Number of periods to generate.\n :param n_samples_per_period: Number of samples per period.\n \"\"\"\n time_delta = self.system.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n data = self._x(time_steps)\n\n return pd.DataFrame({\"t\": time_steps, \"x\": data})\nHarmonicOscillatorIC","text":" Bases: BaseModel
The initial condition for a harmonic oscillator.
Attributes:
Name Type Descriptionx0 float the initial displacement
v0 float the initial velocity
phi float initial phase
Source code inhamilflow/models/harmonic_oscillator.py class HarmonicOscillatorIC(BaseModel):\n \"\"\"The initial condition for a harmonic oscillator.\n\n :cvar x0: the initial displacement\n :cvar v0: the initial velocity\n :cvar phi: initial phase\n \"\"\"\n\n x0: float = Field(default=1.0)\n v0: float = Field(default=0.0)\n phi: float = Field(default=0.0)\nHarmonicOscillatorSystem","text":" Bases: BaseModel
The params for the harmonic oscillator.
Attributes:
Name Type Descriptionomega float angular frequency of the harmonic oscillator
zeta float damping ratio
Source code inhamilflow/models/harmonic_oscillator.py class HarmonicOscillatorSystem(BaseModel):\n \"\"\"The params for the harmonic oscillator.\n\n :cvar omega: angular frequency of the harmonic oscillator\n :cvar zeta: damping ratio\n \"\"\"\n\n omega: float = Field()\n zeta: float = Field(default=0.0)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def period(self) -> float:\n \"\"\"Period of the oscillator.\"\"\"\n return 2 * np.pi / self.omega\n\n @computed_field # type: ignore[misc]\n @cached_property\n def frequency(self) -> float:\n \"\"\"Frequency of the oscillator.\"\"\"\n return 1 / self.period\n\n @computed_field # type: ignore[misc]\n @cached_property\n def type(\n self,\n ) -> Literal[\"simple\", \"under_damped\", \"critical_damped\", \"over_damped\"]:\n \"\"\"Which type of harmonic oscillators.\"\"\"\n if self.zeta == 0:\n return \"simple\"\n elif self.zeta < 1:\n return \"under_damped\"\n elif self.zeta == 1:\n return \"critical_damped\"\n else:\n return \"over_damped\"\n\n @field_validator(\"zeta\")\n @classmethod\n def _check_zeta_non_negative(cls, v: float) -> float:\n if v < 0:\n msg = f\"Value of zeta should be positive: {v=}\"\n raise ValueError(msg)\n\n return v\nfrequency: float cached property","text":"Frequency of the oscillator.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorSystem.period","title":"period: float cached property","text":"Period of the oscillator.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorSystem.type","title":"type: Literal['simple', 'under_damped', 'critical_damped', 'over_damped'] cached property","text":"Which type of harmonic oscillators.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.SimpleHarmonicOscillator","title":"SimpleHarmonicOscillator","text":" Bases: HarmonicOscillatorBase
Generate time series data for a simple harmonic oscillator.
In a one dimensional world, a mass \\(m\\), driven by a force \\(F=-kx\\), is described as
\\[ \\begin{align} F &= - k x \\\\ F &= m a \\end{align} \\]The mass behaves like a simple harmonic oscillator.
In general, the solution to a simple harmonic oscillator is
\\[ x(t) = A \\cos(\\omega t + \\phi), \\]where \\(\\omega\\) is the angular frequency, \\(\\phi\\) is the initial phase, and \\(A\\) is the amplitude.
To use this generator,
params = {\"omega\": omega}\n\nho = SimpleHarmonicOscillator(params=params)\n\ndf = ho(n_periods=1, n_samples_per_period=10)\ndf will be a pandas dataframe with two columns: t and x.
hamilflow/models/harmonic_oscillator.py class SimpleHarmonicOscillator(HarmonicOscillatorBase):\n r\"\"\"Generate time series data for a [simple harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n In a one dimensional world, a mass $m$, driven by a force $F=-kx$, is described as\n\n $$\n \\begin{align}\n F &= - k x \\\\\n F &= m a\n \\end{align}\n $$\n\n The mass behaves like a simple harmonic oscillator.\n\n In general, the solution to a simple harmonic oscillator is\n\n $$\n x(t) = A \\cos(\\omega t + \\phi),\n $$\n\n where $\\omega$ is the angular frequency, $\\phi$ is the initial phase, and $A$ is the amplitude.\n\n\n To use this generator,\n\n ```python\n params = {\"omega\": omega}\n\n ho = SimpleHarmonicOscillator(params=params)\n\n df = ho(n_periods=1, n_samples_per_period=10)\n ```\n\n `df` will be a pandas dataframe with two columns: `t` and `x`.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n super().__init__(system, initial_condition)\n if self.system.type != \"simple\":\n msg = f\"System is not a Simple Harmonic Oscillator: {self.system}\"\n raise ValueError(\n msg,\n )\n\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> np.ndarray:\n r\"\"\"Solution to simple harmonic oscillators.\n\n $$\n x(t) = x_0 \\cos(\\omega t + \\phi).\n $$\n \"\"\"\n return self.initial_condition.x0 * np.cos(\n self.system.omega * np.asarray(t) + self.initial_condition.phi,\n )\nMain module for a harmonic oscillator chain.
"},{"location":"references/models/harmonic_oscillator_chain/#hamilflow.models.harmonic_oscillator_chain.HarmonicOscillatorsChain","title":"HarmonicOscillatorsChain","text":"Generate time series data for a coupled harmonic oscillator chain with periodic boundary condition.
A one-dimensional circle of \\(N\\) interacting harmonic oscillators can be described by the Lagrangian action \\(\\(S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\}\\,,\\)\\) where \\(x_N \\coloneqq x_0\\).
This system can be solved in terms of travelling waves, obtained by discrete Fourier transform.
Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy \\(Y_k = Y^*_{-k \\mod N}\\), see Wikipedia.
Parameters:
Name Type Description Defaultomega float frequence parameter
requiredinitial_conditions Sequence[Mapping[str, float | tuple[float, float]]] a sequence of initial conditions on the Fourier modes. The first element in the sequence is that of the zero mode, taking a position and a velocity. Rest of the elements are that of the independent travelling waves, taking two amplitudes and two initial phases.
requiredodd_dof bool The system will have 2 * len(initial_conditions) + int(odd_dof) - 2 degrees of freedom.
hamilflow/models/harmonic_oscillator_chain.py class HarmonicOscillatorsChain:\n r\"\"\"Generate time series data for a coupled harmonic oscillator chain with periodic boundary condition.\n\n A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action\n $$S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\}\\,,$$\n where $x_N \\coloneqq x_0$.\n\n This system can be solved in terms of _travelling waves_, obtained by discrete Fourier transform.\n\n Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy\n $Y_k = Y^*_{-k \\mod N}$, see [Wikipedia](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#DFT_of_real_and_purely_imaginary_signals).\n\n :param omega: frequence parameter\n :param initial_conditions: a sequence of initial conditions on the Fourier modes.\n The first element in the sequence is that of the zero mode, taking a position and a velocity.\n Rest of the elements are that of the independent travelling waves, taking two amplitudes and two initial phases.\n :param odd_dof: The system will have `2 * len(initial_conditions) + int(odd_dof) - 2` degrees of freedom.\n \"\"\"\n\n def __init__(\n self,\n omega: float,\n initial_conditions: Sequence[Mapping[str, float | tuple[float, float]]],\n odd_dof: bool,\n ) -> None:\n self.n_dof = 2 * len(initial_conditions) + odd_dof - 2\n if not odd_dof:\n prefix = \"For even degrees of freedom, \"\n if self.n_dof == 0:\n raise ValueError(prefix + \"at least 1 travelling wave is needed\")\n amp = cast(tuple[float, float], initial_conditions[-1][\"amp\"])\n if amp[0] != amp[1]:\n msg = \"k == N // 2 must have equal positive and negative amplitudes.\"\n raise ValueError(prefix + msg)\n self.omega = omega\n self.odd_dof = odd_dof\n\n self.free_mode = FreeParticle(cast(Mapping[str, float], initial_conditions[0]))\n\n self.independent_csho_modes = [\n self._sho_factory(\n k,\n cast(tuple[float, float], ic[\"amp\"]),\n cast(tuple[float, float] | None, ic.get(\"phi\")),\n )\n for k, ic in enumerate(initial_conditions[1:], 1)\n ]\n\n def _sho_factory(\n self,\n k: int,\n amp: tuple[float, float],\n phi: tuple[float, float] | None = None,\n ) -> ComplexSimpleHarmonicOscillator:\n return ComplexSimpleHarmonicOscillator(\n {\n \"omega\": 2 * self.omega * np.sin(np.pi * k / self.n_dof),\n },\n {\"x0\": amp} | ({\"phi\": phi} if phi else {}),\n )\n\n @cached_property\n def definition(\n self,\n ) -> dict[\n str,\n float\n | dict[str, dict[str, float | list[float]]]\n | list[dict[str, dict[str, float | tuple[float, float]]]],\n ]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"omega\": self.omega,\n \"n_dof\": self.n_dof,\n \"free_mode\": self.free_mode.definition,\n \"independent_csho_modes\": [\n rwm.definition for rwm in self.independent_csho_modes\n ],\n }\n\n def _z(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"tuple[npt.NDArray[np.complex64], npt.NDArray[np.complex64]]\":\n t = np.asarray(t).reshape(-1)\n all_travelling_waves = [self.free_mode._x(t).reshape(1, -1)] # noqa: SLF001\n\n if self.independent_csho_modes:\n independent_cshos = np.asarray(\n [o._z(t) for o in self.independent_csho_modes], # noqa: SLF001\n )\n all_travelling_waves.extend(\n (\n (independent_cshos, independent_cshos[::-1].conj())\n if self.odd_dof\n else (\n independent_cshos[:-1],\n independent_cshos[[-1]],\n independent_cshos[-2::-1].conj(),\n )\n ),\n )\n\n travelling_waves = np.concatenate(all_travelling_waves)\n original_zs = ifft(travelling_waves, axis=0, norm=\"ortho\")\n return original_zs, travelling_waves\n\n def _x(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"tuple[npt.NDArray[np.float64], npt.NDArray[np.complex64]]\":\n original_xs, travelling_waves = self._z(t)\n\n return np.real(original_xs), travelling_waves\n\n def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator chain.\n\n Returns float(s) representing the displacement at the given time(s).\n\n :param t: time.\n \"\"\"\n t = np.asarray(t)\n original_xs, travelling_waves = self._x(t)\n data = { # type: ignore [var-annotated]\n f\"{name}{i}\": cast(\n \"npt.NDArray[np.float64] | npt.NDArray[np.complex64]\",\n values,\n )\n for name, xs in zip(\n (\"x\", \"y\"),\n (original_xs, travelling_waves),\n strict=False,\n )\n for i, values in enumerate(xs) # type: ignore [arg-type]\n }\n\n return pd.DataFrame(data, index=t)\ndefinition: dict[str, float | dict[str, dict[str, float | list[float]]] | list[dict[str, dict[str, float | tuple[float, float]]]]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator_chain/#hamilflow.models.harmonic_oscillator_chain.HarmonicOscillatorsChain.__call__","title":"__call__(t)","text":"Generate time series data for the harmonic oscillator chain.
Returns float(s) representing the displacement at the given time(s).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time.
required Source code inhamilflow/models/harmonic_oscillator_chain.py def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator chain.\n\n Returns float(s) representing the displacement at the given time(s).\n\n :param t: time.\n \"\"\"\n t = np.asarray(t)\n original_xs, travelling_waves = self._x(t)\n data = { # type: ignore [var-annotated]\n f\"{name}{i}\": cast(\n \"npt.NDArray[np.float64] | npt.NDArray[np.complex64]\",\n values,\n )\n for name, xs in zip(\n (\"x\", \"y\"),\n (original_xs, travelling_waves),\n strict=False,\n )\n for i, values in enumerate(xs) # type: ignore [arg-type]\n }\n\n return pd.DataFrame(data, index=t)\nMain module for a pendulum.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum","title":"Pendulum","text":"Generate time series data for a pendulum.
We describe a generic pendulum system by the Lagrangian action $$ S_L[\\theta] = I \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\}\\,, $$ where \\(\\theta\\) is the angle from the vertical to the pendulum; \\(I\\) is the inertia parameter introduced for dimensional reasons, and \\(\\omega_0\\) the frequency parameter.
Details are collected in the tutorial.
Source code inhamilflow/models/pendulum.py class Pendulum:\n r\"\"\"Generate time series data for a pendulum.\n\n We describe a generic pendulum system by the Lagrangian action\n $$\n S_L\\[\\theta\\] = I \\int_{t_0}^{t_1} \\mathbb{d}t\n \\left\\\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\\\}\\,,\n $$\n where $\\theta$ is the _angle_ from the vertical to the pendulum;\n $I$ is the _inertia parameter_ introduced for dimensional reasons,\n and $\\omega_0$ the _frequency parameter_.\n\n Details are collected in the tutorial.\n \"\"\"\n\n def __init__(\n self,\n system: float | Mapping[str, float],\n initial_condition: float | Mapping[str, float],\n ) -> None:\n if isinstance(system, float | int):\n system = {\"omega0\": system}\n if isinstance(initial_condition, float | int):\n initial_condition = {\"theta0\": initial_condition}\n self.system = PendulumSystem.model_validate(system)\n self.initial_condition = PendulumIC.model_validate(initial_condition)\n\n @cached_property\n def definition(self) -> dict[str, dict[str, Any]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n @property\n def omega0(self) -> float:\n \"\"\"Original angular frequency of the system.\"\"\"\n return self.system.omega0\n\n @property\n def _k(self) -> float:\n return self.initial_condition.k\n\n @property\n def _math_m(self) -> float:\n return self._k**2\n\n @cached_property\n def freq(self) -> float:\n r\"\"\"Frequency.\n\n :return: $\\frac{\\pi}{2K(k^2)}\\omega_0$, where\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n return math.pi * self.omega0 / (2 * ellipk(self._math_m))\n\n @cached_property\n def period(self) -> float:\n r\"\"\"Period.\n\n :return: $\\frac{4K(k^2)}{\\omega_0}$, where\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n return 4 * ellipk(self._math_m) / self.omega0\n\n def _math_u(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"npt.NDArray[np.float64]\":\n return self.omega0 * np.asarray(t)\n\n def u(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Give the convenient generalised coordinate $u$, $\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}$.\n\n :param t: time\n :return: $u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)$, where\n $\\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)\n\n return ph\n\n def theta(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Angle $\\theta$.\n\n :param t: time\n :return: $\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)$, where\n $\\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)\n \"\"\"\n _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)\n\n return 2 * np.arcsin(cn / dn * self._k)\n\n def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate the time sequence from more interpretable params.\n\n :param n_periods: number of periods to include\n :param n_samples_per_period: number of samples in each period\n :return: an array that contains all the timesteps\n \"\"\"\n time_delta = self.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n return self(time_steps)\n\n def __call__(self, t: TypeTime) -> pd.DataFrame:\n \"\"\"Generate the variables of the pendulum in time together with the time steps.\n\n :param t: time steps\n :return: values of the variables\n angle `x`, generalized coordinates `u`, and time `t`.\n \"\"\"\n thetas = self.theta(t)\n us = self.u(t)\n\n return pd.DataFrame({\"t\": t, \"x\": thetas, \"u\": us})\ndefinition: dict[str, dict[str, Any]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.freq","title":"freq: float cached property","text":"Frequency.
Returns:
Type Descriptionfloat \\(\\frac{\\pi}{2K(k^2)}\\omega_0\\), where \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.omega0","title":"omega0: float property","text":"Original angular frequency of the system.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.period","title":"period: float cached property","text":"Period.
Returns:
Type Descriptionfloat \\(\\frac{4K(k^2)}{\\omega_0}\\), where \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.__call__","title":"__call__(t)","text":"Generate the variables of the pendulum in time together with the time steps.
Parameters:
Name Type Description Defaultt TypeTime time steps
requiredReturns:
Type DescriptionDataFrame values of the variables angle x, generalized coordinates u, and time t.
hamilflow/models/pendulum.py def __call__(self, t: TypeTime) -> pd.DataFrame:\n \"\"\"Generate the variables of the pendulum in time together with the time steps.\n\n :param t: time steps\n :return: values of the variables\n angle `x`, generalized coordinates `u`, and time `t`.\n \"\"\"\n thetas = self.theta(t)\n us = self.u(t)\n\n return pd.DataFrame({\"t\": t, \"x\": thetas, \"u\": us})\ngenerate_from(n_periods, n_samples_per_period)","text":"Generate the time sequence from more interpretable params.
Parameters:
Name Type Description Defaultn_periods int number of periods to include
requiredn_samples_per_period int number of samples in each period
requiredReturns:
Type DescriptionDataFrame an array that contains all the timesteps
Source code inhamilflow/models/pendulum.py def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate the time sequence from more interpretable params.\n\n :param n_periods: number of periods to include\n :param n_samples_per_period: number of samples in each period\n :return: an array that contains all the timesteps\n \"\"\"\n time_delta = self.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n return self(time_steps)\ntheta(t)","text":"Angle \\(\\theta\\).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time
requiredReturns:
Type DescriptionNDArray[float64] \\(\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)\\), where \\(\\mathrm{cd}(z, k)\\) is a Jacobian elliptic function
Source code inhamilflow/models/pendulum.py def theta(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Angle $\\theta$.\n\n :param t: time\n :return: $\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)$, where\n $\\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)\n \"\"\"\n _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)\n\n return 2 * np.arcsin(cn / dn * self._k)\nu(t)","text":"Give the convenient generalised coordinate \\(u\\), \\(\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}\\).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time
requiredReturns:
Type DescriptionNDArray[float64] \\(u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)\\), where \\(\\mathrm{am}(x, k)\\) is Jacobi's amplitude function, \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
Source code inhamilflow/models/pendulum.py def u(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Give the convenient generalised coordinate $u$, $\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}$.\n\n :param t: time\n :return: $u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)$, where\n $\\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)\n\n return ph\nPendulumIC","text":" Bases: BaseModel
The initial condition for a pendulum.
Parameters:
Name Type Description Defaulttheta0 \\(-\\frac{\\pi}{2} \\le \\theta_0 \\le \\frac{\\pi}{2}\\), the initial angle
required Source code inhamilflow/models/pendulum.py class PendulumIC(BaseModel):\n r\"\"\"The initial condition for a pendulum.\n\n :param theta0: $-\\frac{\\pi}{2} \\le \\theta_0 \\le \\frac{\\pi}{2}$, the\n initial angle\n \"\"\"\n\n theta0: float = Field(ge=-math.pi / 2, le=math.pi / 2, frozen=True)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def k(self) -> float:\n r\"\"\"A convenient number for elliptic functions.\n\n :return: $\\sin\\frac{\\theta_0}{2}$\n \"\"\"\n return math.sin(self.theta0 / 2)\nk: float cached property","text":"A convenient number for elliptic functions.
Returns:
Type Descriptionfloat \\(\\sin\\frac{\\theta_0}{2}\\)
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.PendulumSystem","title":"PendulumSystem","text":" Bases: BaseModel
The params for the pendulum.
Parameters:
Name Type Description Defaultomega0 \\(\\omega_0 \\coloneqq \\sqrt{\\frac{U}{I}} > 0\\), frequency parameter
required Source code inhamilflow/models/pendulum.py class PendulumSystem(BaseModel):\n r\"\"\"The params for the pendulum.\n\n :param omega0: $\\omega_0 \\coloneqq \\sqrt{\\frac{U}{I}} > 0$, frequency\n parameter\n \"\"\"\n\n omega0: float = Field(gt=0.0, frozen=True)\nExact solution of Kepler dynamics.
"},{"location":"references/models/kepler_problem/dynamics/#hamilflow.models.kepler_problem.dynamics.esolve_u_from_tau_parabolic","title":"esolve_u_from_tau_parabolic(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.
Let \\(\\kappa = 1+9\\tau^2\\), $$ u = -1 + \\left(\\frac{3\\tau}{\\kappa^\\frac{3}{2}} + \\frac{1}{\\kappa}\\right)^\\frac{1}{3} + \\left(\\frac{1}{3\\tau \\kappa^\\frac{3}{2}+\\kappa^2}\\right)^\\frac{1}{3}\\,. $$
Parameters:
Name Type Description Defaultecc float eccentricity, ecc = 0 (unchecked, unused)
requiredtau ArrayLike scaled time
requiredReturns:
Type DescriptionNDArray[float64] convenient radial inverse u
Source code inhamilflow/models/kepler_problem/dynamics.py def esolve_u_from_tau_parabolic(\n ecc: float, # noqa: ARG001\n tau: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.\n\n Let $\\kappa = 1+9\\tau^2$,\n $$ u = -1\n + \\left(\\frac{3\\tau}{\\kappa^\\frac{3}{2}} + \\frac{1}{\\kappa}\\right)^\\frac{1}{3}\n + \\left(\\frac{1}{3\\tau \\kappa^\\frac{3}{2}+\\kappa^2}\\right)^\\frac{1}{3}\\,. $$\n\n :param ecc: eccentricity, ecc = 0 (unchecked, unused)\n :param tau: scaled time\n :return: convenient radial inverse u\n \"\"\"\n tau = np.asarray(tau)\n tau_3 = 3 * tau\n term = 1 + tau_3**2 # 1 + 9 * tau**2\n term1_5 = term**1.5 # (1 + 9 * tau**2)**1.5\n second_term = (tau_3 / term1_5 + 1 / term) ** _1_3\n third_term = 1 / (tau_3 * term1_5 + term**2) ** _1_3\n return -1 + second_term + third_term\ntau_of_1_plus_u_hyperbolic(ecc, u)","text":"Expansion for tau of u in the hyperbolic case at \\(u = -1+0\\).
The exact solution has a logarithmic singularity at \\(u = -1\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\frac{e(1+u)}{2(e^2-1)}\\), $$ \\tau = \\ln \\epsilon + \\frac{e}{2\\epsilon} + 1 - \\frac{e^2+2}{e}\\epsilon + \\left(3+\\frac{2}{e^2}\\right)\\epsilon^2 + O(\\epsilon^3)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_1_plus_u_hyperbolic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the hyperbolic case at $u = -1+0$.\n\n The exact solution has a logarithmic singularity at $u = -1$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\frac{e(1+u)}{2(e^2-1)}$,\n $$ \\tau = \\ln \\epsilon + \\frac{e}{2\\epsilon}\n + 1 - \\frac{e^2+2}{e}\\epsilon + \\left(3+\\frac{2}{e^2}\\right)\\epsilon^2\n + O(\\epsilon^3)\\,. $$\n \"\"\"\n cosqr = ecc**2 - 1\n up1 = ecc * (1 + u) / 2 / cosqr\n diverge = np.log(up1) + ecc / 2 / up1\n regular = 1 - (2 + ecc**2) / ecc * up1 + (3 + 2 / ecc**2) * up1**2\n return (diverge + regular) / cosqr**1.5\ntau_of_e_minus_u_elliptic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = +e-0\\).
The exact solution has a removable singularity at \\(u = +e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}\\), $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon - \\frac{1+9e}{24(1+e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_minus_u_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = +e-0$.\n\n The exact solution has a removable singularity at $u = +e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}$,\n $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon\n - \\frac{1+9e}{24(1+e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n emu = np.sqrt(2 * (ecc - u) / ecc)\n return emu / (1 + ecc) ** 2 - emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3\ntau_of_e_minus_u_hyperbolic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = +e-0\\).
The exact solution has a removable singularity at \\(u = +e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}\\), $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon + \\frac{1+9e}{24(1+e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_minus_u_hyperbolic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = +e-0$.\n\n The exact solution has a removable singularity at $u = +e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}$,\n $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon\n + \\frac{1+9e}{24(1+e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n emu = np.sqrt(2 * (ecc - u) / ecc)\n return emu / (1 + ecc) ** 2 + emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3\ntau_of_e_plus_u_elliptic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = -e+0\\).
The exact solution has a removable singularity at \\(u = -e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e+u)}{e}}\\), $$ \\tau = \\frac{\\pi}{(1-e^2)^\\frac{3}{2}} - \\frac{1}{(1-e)^2}\\epsilon - \\frac{1-9e}{24(1-e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_plus_u_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = -e+0$.\n\n The exact solution has a removable singularity at $u = -e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e+u)}{e}}$,\n $$ \\tau = \\frac{\\pi}{(1-e^2)^\\frac{3}{2}}\n - \\frac{1}{(1-e)^2}\\epsilon\n - \\frac{1-9e}{24(1-e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n epu = np.sqrt(2 * (ecc + u) / ecc)\n const = np.pi / (1 - ecc**2) ** 1.5\n return const - epu / (ecc - 1) ** 2 - epu**3 * (1 - 9 * ecc) / 24 / (1 - ecc) ** 3\ntau_of_u_elliptic(ecc, u)","text":"Calculate the scaled time tau from u in the elliptic case.
Parameters:
Name Type Description Defaultecc float eccentricity, 0 < ecc < 1 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_elliptic(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the scaled time tau from u in the elliptic case.\n\n :param ecc: eccentricity, 0 < ecc < 1 (unchecked)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n return _approximate_at_termina(\n ecc,\n u,\n tau_of_u_exact_elliptic,\n tau_of_e_plus_u_elliptic,\n tau_of_e_minus_u_elliptic,\n )\ntau_of_u_exact_elliptic(ecc, u)","text":"Exact solution for tau of u in the elliptic case.
For \\(-e \\le u \\le e\\), $$ \\tau = -\\frac{\\sqrt{e^2-u^2}}{(1-e^2)(1+u)} + \\frac{\\frac{\\pi}{2} - \\arctan\\frac{e^2+u}{\\sqrt{(1-e^2)(e^2-u^2)}}}{(1-e^2)^{\\frac{3}{2}}}\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_exact_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Exact solution for tau of u in the elliptic case.\n\n For $-e \\le u \\le e$,\n $$ \\tau = -\\frac{\\sqrt{e^2-u^2}}{(1-e^2)(1+u)}\n + \\frac{\\frac{\\pi}{2} - \\arctan\\frac{e^2+u}{\\sqrt{(1-e^2)(e^2-u^2)}}}{(1-e^2)^{\\frac{3}{2}}}\\,. $$\n \"\"\"\n cosqr, eusqrt = 1 - ecc**2, np.sqrt(ecc**2 - u**2)\n trig_numer = np.pi / 2 - np.arctan((ecc**2 + u) / np.sqrt(cosqr) / eusqrt)\n return -eusqrt / cosqr / (1 + u) + trig_numer / cosqr**1.5\ntau_of_u_exact_hyperbolic(ecc, u)","text":"Exact solution for tau of u in the hyperbolic case.
For \\(-1 < u \\le e\\), $$ \\tau = \\frac{\\sqrt{e^2-u^2}}{(e^2-1)(1+u)} - \\frac{\\mathrm{artanh}\\frac{e^2+u}{\\sqrt{(e^2-1)(e^2-u^2)}}}{(e^2-1)^{\\frac{3}{2}}}\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_exact_hyperbolic(\n ecc: float,\n u: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Exact solution for tau of u in the hyperbolic case.\n\n For $-1 < u \\le e$,\n $$ \\tau = \\frac{\\sqrt{e^2-u^2}}{(e^2-1)(1+u)}\n - \\frac{\\mathrm{artanh}\\frac{e^2+u}{\\sqrt{(e^2-1)(e^2-u^2)}}}{(e^2-1)^{\\frac{3}{2}}}\\,. $$\n \"\"\"\n u = np.asarray(u)\n cosqr, eusqrt = ecc**2 - 1, np.sqrt(ecc**2 - u**2)\n trig_numer = np.arctanh(np.sqrt(cosqr) * eusqrt / (ecc**2 + u))\n return eusqrt / cosqr / (1 + u) - trig_numer / cosqr**1.5\ntau_of_u_hyperbolic(ecc, u)","text":"Calculate the scaled time tau from u in the hyperbolic case.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 1 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_hyperbolic(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the scaled time tau from u in the hyperbolic case.\n\n :param ecc: eccentricity, ecc > 1 (unchecked)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n return _approximate_at_termina(\n ecc,\n u,\n tau_of_u_exact_hyperbolic,\n tau_of_1_plus_u_hyperbolic,\n tau_of_e_minus_u_hyperbolic,\n )\ntau_of_u_parabolic(ecc, u)","text":"Calculate the scaled time tau from u in the parabolic case.
For \\(-1 < u \\le 1\\), $$ \\tau = \\frac{\\sqrt{1-u}(2+u)}{3(1+u)^\\frac{3}{2}}\\,. $$
Parameters:
Name Type Description Defaultecc float eccentricity, ecc == 1 (unchecked, unused)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_parabolic(\n ecc: float, # noqa: ARG001\n u: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the scaled time tau from u in the parabolic case.\n\n For $-1 < u \\le 1$,\n $$ \\tau = \\frac{\\sqrt{1-u}(2+u)}{3(1+u)^\\frac{3}{2}}\\,. $$\n\n :param ecc: eccentricity, ecc == 1 (unchecked, unused)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n u = np.asarray(u)\n return np.sqrt(1 - u) * (2 + u) / 3 / (1 + u) ** 1.5\ntau_of_u_prime(ecc, u)","text":"Calculate the first derivative of scaled time tau with respect to u.
\\[ \\tau'(u) = -\\frac{1}{(1+u)^2\\sqrt{e^2-u^2}}\\,. \\]Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] the first derivative scaled time tau with respect to u
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_prime(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the first derivative of scaled time tau with respect to u.\n\n $$ \\tau'(u) = -\\frac{1}{(1+u)^2\\sqrt{e^2-u^2}}\\,. $$\n\n :param ecc: eccentricity, ecc >= 0 (unchecked)\n :param u: convenient radial inverse\n :return: the first derivative scaled time tau with respect to u\n \"\"\"\n u = np.asarray(u)\n return -1 / (1 + u) ** 2 / np.sqrt(ecc**2 - u**2)\ntau_of_u_prime2(ecc, u)","text":"Calculate the second derivative of scaled time tau with respect to u.
\\[ \\tau''(u) = \\frac{2e^2 - u - 3u^2}{(1+u)^3 (e^2-u^2)^\\frac{3}{2}} \\]Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] the second derivative scaled time tau with respect to u
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_prime2(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the second derivative of scaled time tau with respect to u.\n\n $$ \\tau''(u) = \\frac{2e^2 - u - 3u^2}{(1+u)^3 (e^2-u^2)^\\frac{3}{2}} $$\n\n :param ecc: eccentricity, ecc >= 0 (unchecked)\n :param u: convenient radial inverse\n :return: the second derivative scaled time tau with respect to u\n \"\"\"\n u = np.asarray(u)\n u2 = u**2\n return (2 * ecc**2 - u - 3 * u2) / (1 + u) ** 3 / (ecc**2 - u2) ** 1.5\nMain module for Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D","title":"Kepler2D","text":"Kepler problem in two dimensional space.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] the Kepler problem system specification
requiredfirst_integrals Mapping[str, float] the first integrals for the system.
required Source code inhamilflow/models/kepler_problem/model.py class Kepler2D:\n \"\"\"Kepler problem in two dimensional space.\n\n :param system: the Kepler problem system specification\n :param first_integrals: the first integrals for the system.\n \"\"\"\n\n def __init__(\n self,\n system: \"Mapping[str, float]\",\n first_integrals: \"Mapping[str, float]\",\n ) -> None:\n self.system = Kepler2DSystem.model_validate(system)\n\n first_integrals = dict(first_integrals)\n ene = first_integrals[\"ene\"]\n minimal_ene = Kepler2D.minimal_ene(first_integrals[\"angular_mom\"], system)\n if ene < minimal_ene:\n msg = f\"Energy {ene} less than minimally allowed {minimal_ene}\"\n raise ValueError(msg)\n\n self.first_integrals = Kepler2DFI.model_validate(first_integrals)\n\n if 0 <= self.ecc < 1:\n self.tau_of_u = partial(tau_of_u_elliptic, self.ecc)\n elif self.ecc == 1:\n self.tau_of_u = partial(tau_of_u_parabolic, self.ecc)\n elif self.ecc > 1:\n self.tau_of_u = partial(tau_of_u_hyperbolic, self.ecc)\n else:\n raise RuntimeError\n\n @classmethod\n def from_geometry(\n cls,\n system: \"Mapping[str, float]\",\n geometries: \"Mapping[str, bool | float]\",\n ) -> \"Self\":\n r\"\"\"Alternative initialiser from system and geometry specifications.\n\n Given the eccentricity $e$ and the conic section parameter $p$,\n $$l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.$$\n\n :param system: the Kepler problem system specification\n :param geometries: geometric specifications\n `positive_angular_mom`: whether the angular momentum is positive\n `ecc`: eccentricity of the conic section\n `parameter`: parameter of the conic section\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n positive_angular_mom = bool(geometries[\"positive_angular_mom\"])\n ecc, parameter = (float(geometries[k]) for k in [\"ecc\", \"parameter\"])\n abs_angular_mom = math.sqrt(mass * parameter * alpha)\n # abs_minimal_ene = alpha / 2 / parameter: numerically unstable\n abs_minimal_ene = abs(cls.minimal_ene(abs_angular_mom, system))\n ene = (ecc**2 - 1) * abs_minimal_ene\n fi = {\n \"ene\": ene,\n \"angular_mom\": (\n abs_angular_mom if positive_angular_mom else -abs_angular_mom\n ),\n }\n return cls(system, fi)\n\n @staticmethod\n def minimal_ene(\n angular_mom: float,\n system: \"Mapping[str, float]\",\n ) -> float:\n r\"\"\"Minimal possible energy from the system specification and an angular momentum.\n\n $$ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. $$\n\n :param angular_mom: angular momentum\n :param system: system specification\n :return: minimal possible energy\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n return -mass * alpha**2 / (2 * angular_mom**2)\n\n @property\n def mass(self) -> float:\n \"\"\"Mass $m$ from the system specification.\"\"\"\n return self.system.mass\n\n @property\n def alpha(self) -> float:\n r\"\"\"Alpha $\\alpha$ from the system specification.\"\"\"\n return self.system.alpha\n\n @property\n def ene(self) -> float:\n \"\"\"Energy $E$ of the Kepler problem.\"\"\"\n return self.first_integrals.ene\n\n @property\n def angular_mom(self) -> float:\n \"\"\"Angular momentum $l$ of the Kepler problem.\"\"\"\n return self.first_integrals.angular_mom\n\n @property\n def t0(self) -> float:\n r\"\"\"t0 $t_0$ of the Kepler problem.\"\"\"\n return self.first_integrals.t0\n\n @property\n def phi0(self) -> float:\n r\"\"\"phi0 $\\phi_0$ of the Kepler problem.\"\"\"\n return self.first_integrals.phi0\n\n @cached_property\n def period(self) -> float:\n r\"\"\"Period $T$ of the Kepler problem.\n\n For $E < 0$,\n $$ T = \\pi \\alpha \\sqrt{-\\frac{m}{2E^3}}\\,. $$\n \"\"\"\n if self.ene >= 0:\n msg = f\"Only systems with energy < 0 have a period, got {self.ene}\"\n raise TypeError(msg)\n return math.pi * self.alpha * math.sqrt(-self.mass / 2 / self.ene**3)\n\n # FIXME is it called parameter in English?\n @cached_property\n def parameter(self) -> float:\n r\"\"\"Conic section parameter of the Kepler problem.\n\n $$ p = \\frac{l^2}{\\alpha m}\\,. $$\n \"\"\"\n return self.angular_mom**2 / self.mass / self.alpha\n\n @cached_property\n def ecc(self) -> float:\n r\"\"\"Conic section eccentricity of the Kepler problem.\n\n $$ e = \\sqrt{1 + \\frac{2El}{\\alpha^2 m}}\\,. $$\n \"\"\"\n return math.sqrt(\n 1 + 2 * self.ene * self.angular_mom**2 / self.mass / self.alpha**2,\n )\n\n @cached_property\n def period_in_tau(self) -> float:\n r\"\"\"Period in the scaled time tau.\n\n $$ T_\\tau = \\frac{2\\pi}{(1-e^2)^\\frac{3}{2}}\\,. $$\n \"\"\"\n if self.ecc >= 1:\n msg = (\n f\"Only systems with 0 <= eccentricity < 1 have a period, got {self.ecc}\"\n )\n raise TypeError(\n msg,\n )\n return 2 * math.pi / (1 - self.ecc**2) ** 1.5\n\n @property\n def t_to_tau_factor(self) -> float:\n r\"\"\"Scale factor from t to tau.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return abs(self.mass * self.alpha**2 / self.angular_mom**3)\n\n def tau(self, t: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the scaled time tau from t.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return (np.asarray(t) - self.t0) * self.t_to_tau_factor\n\n def u_of_tau(self, tau: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n \"\"\"Give the convenient radial inverse u from tau.\"\"\"\n tau = np.asarray(tau)\n if self.ecc == 0:\n return np.zeros(tau.shape)\n else:\n if self.ecc < 1:\n p = self.period_in_tau\n r = tau % p\n tau = np.where(r <= p / 2, r, p - r)\n else:\n tau = np.abs(tau)\n return u_of_tau(self.ecc, tau) # type: ignore [arg-type]\n\n def r_of_u(self, u: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the radial r from u.\n\n $$ r = \\frac{p}{u+1}\\,. $$\n \"\"\"\n return self.parameter / (np.asarray(u) + 1)\n\n def phi_of_u_tau(\n self,\n u: \"Collection[float] | npt.ArrayLike\",\n tau: \"Collection[float] | npt.ArrayLike\",\n ) -> \"npt.ArrayLike\":\n r\"\"\"Give the angular phi from u and tau.\n\n For $e = 0$,\n $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$\n For $e > 0$,\n $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$\n \"\"\"\n u, tau = np.asarray(u), np.asarray(tau)\n if self.ecc == 0:\n phi = 2 * math.pi * tau / self.period_in_tau\n else:\n if self.ecc < 1:\n shift = tau / self.period_in_tau\n else:\n shift = np.where(tau >= 0, 0, -np.pi)\n phi = acos_with_shift(u / self.ecc, shift) # type: ignore [assignment]\n return phi + self.phi0\n\n def __call__(self, t: \"Collection[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Give a DataFrame of tau, u, r and phi from t.\"\"\"\n tau = self.tau(t)\n u = self.u_of_tau(tau)\n r = self.r_of_u(u)\n phi = self.phi_of_u_tau(u, tau)\n\n return pd.DataFrame({\"t\": t, \"tau\": tau, \"u\": u, \"r\": r, \"phi\": phi})\nalpha: float property","text":"Alpha \\(\\alpha\\) from the system specification.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.angular_mom","title":"angular_mom: float property","text":"Angular momentum \\(l\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.ecc","title":"ecc: float cached property","text":"Conic section eccentricity of the Kepler problem.
\\[ e = \\sqrt{1 + \\frac{2El}{\\alpha^2 m}}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.ene","title":"ene: float property","text":"Energy \\(E\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.mass","title":"mass: float property","text":"Mass \\(m\\) from the system specification.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.parameter","title":"parameter: float cached property","text":"Conic section parameter of the Kepler problem.
\\[ p = \\frac{l^2}{\\alpha m}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.period","title":"period: float cached property","text":"Period \\(T\\) of the Kepler problem.
For \\(E < 0\\), $$ T = \\pi \\alpha \\sqrt{-\\frac{m}{2E^3}}\\,. $$
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.period_in_tau","title":"period_in_tau: float cached property","text":"Period in the scaled time tau.
\\[ T_\\tau = \\frac{2\\pi}{(1-e^2)^\\frac{3}{2}}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.phi0","title":"phi0: float property","text":"phi0 \\(\\phi_0\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.t0","title":"t0: float property","text":"t0 \\(t_0\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.t_to_tau_factor","title":"t_to_tau_factor: float property","text":"Scale factor from t to tau.
\\[ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.__call__","title":"__call__(t)","text":"Give a DataFrame of tau, u, r and phi from t.
Source code inhamilflow/models/kepler_problem/model.py def __call__(self, t: \"Collection[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Give a DataFrame of tau, u, r and phi from t.\"\"\"\n tau = self.tau(t)\n u = self.u_of_tau(tau)\n r = self.r_of_u(u)\n phi = self.phi_of_u_tau(u, tau)\n\n return pd.DataFrame({\"t\": t, \"tau\": tau, \"u\": u, \"r\": r, \"phi\": phi})\nfrom_geometry(system, geometries) classmethod","text":"Alternative initialiser from system and geometry specifications.
Given the eccentricity \\(e\\) and the conic section parameter \\(p\\), \\(\\(l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.\\)\\)
Parameters:
Name Type Description Defaultsystem Mapping[str, float] the Kepler problem system specification
requiredgeometries Mapping[str, bool | float] geometric specifications positive_angular_mom: whether the angular momentum is positive ecc: eccentricity of the conic section parameter: parameter of the conic section
hamilflow/models/kepler_problem/model.py @classmethod\ndef from_geometry(\n cls,\n system: \"Mapping[str, float]\",\n geometries: \"Mapping[str, bool | float]\",\n) -> \"Self\":\n r\"\"\"Alternative initialiser from system and geometry specifications.\n\n Given the eccentricity $e$ and the conic section parameter $p$,\n $$l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.$$\n\n :param system: the Kepler problem system specification\n :param geometries: geometric specifications\n `positive_angular_mom`: whether the angular momentum is positive\n `ecc`: eccentricity of the conic section\n `parameter`: parameter of the conic section\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n positive_angular_mom = bool(geometries[\"positive_angular_mom\"])\n ecc, parameter = (float(geometries[k]) for k in [\"ecc\", \"parameter\"])\n abs_angular_mom = math.sqrt(mass * parameter * alpha)\n # abs_minimal_ene = alpha / 2 / parameter: numerically unstable\n abs_minimal_ene = abs(cls.minimal_ene(abs_angular_mom, system))\n ene = (ecc**2 - 1) * abs_minimal_ene\n fi = {\n \"ene\": ene,\n \"angular_mom\": (\n abs_angular_mom if positive_angular_mom else -abs_angular_mom\n ),\n }\n return cls(system, fi)\nminimal_ene(angular_mom, system) staticmethod","text":"Minimal possible energy from the system specification and an angular momentum.
\\[ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. \\]Parameters:
Name Type Description Defaultangular_mom float angular momentum
requiredsystem Mapping[str, float] system specification
requiredReturns:
Type Descriptionfloat minimal possible energy
Source code inhamilflow/models/kepler_problem/model.py @staticmethod\ndef minimal_ene(\n angular_mom: float,\n system: \"Mapping[str, float]\",\n) -> float:\n r\"\"\"Minimal possible energy from the system specification and an angular momentum.\n\n $$ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. $$\n\n :param angular_mom: angular momentum\n :param system: system specification\n :return: minimal possible energy\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n return -mass * alpha**2 / (2 * angular_mom**2)\nphi_of_u_tau(u, tau)","text":"Give the angular phi from u and tau.
For \\(e = 0\\), $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$ For \\(e > 0\\), $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$
Source code inhamilflow/models/kepler_problem/model.py def phi_of_u_tau(\n self,\n u: \"Collection[float] | npt.ArrayLike\",\n tau: \"Collection[float] | npt.ArrayLike\",\n) -> \"npt.ArrayLike\":\n r\"\"\"Give the angular phi from u and tau.\n\n For $e = 0$,\n $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$\n For $e > 0$,\n $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$\n \"\"\"\n u, tau = np.asarray(u), np.asarray(tau)\n if self.ecc == 0:\n phi = 2 * math.pi * tau / self.period_in_tau\n else:\n if self.ecc < 1:\n shift = tau / self.period_in_tau\n else:\n shift = np.where(tau >= 0, 0, -np.pi)\n phi = acos_with_shift(u / self.ecc, shift) # type: ignore [assignment]\n return phi + self.phi0\nr_of_u(u)","text":"Give the radial r from u.
\\[ r = \\frac{p}{u+1}\\,. \\] Source code inhamilflow/models/kepler_problem/model.py def r_of_u(self, u: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the radial r from u.\n\n $$ r = \\frac{p}{u+1}\\,. $$\n \"\"\"\n return self.parameter / (np.asarray(u) + 1)\ntau(t)","text":"Give the scaled time tau from t.
\\[ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. \\] Source code inhamilflow/models/kepler_problem/model.py def tau(self, t: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the scaled time tau from t.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return (np.asarray(t) - self.t0) * self.t_to_tau_factor\nu_of_tau(tau)","text":"Give the convenient radial inverse u from tau.
Source code inhamilflow/models/kepler_problem/model.py def u_of_tau(self, tau: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n \"\"\"Give the convenient radial inverse u from tau.\"\"\"\n tau = np.asarray(tau)\n if self.ecc == 0:\n return np.zeros(tau.shape)\n else:\n if self.ecc < 1:\n p = self.period_in_tau\n r = tau % p\n tau = np.where(r <= p / 2, r, p - r)\n else:\n tau = np.abs(tau)\n return u_of_tau(self.ecc, tau) # type: ignore [arg-type]\nKepler2DFI","text":" Bases: BaseModel
The first integrals for a Kepler problem.
Attributes:
Name Type Descriptionene float the energy \\(E\\)
angular_mom float the angular momentum \\(l\\)
t0 float the time \\(t_0\\) at which the radial position is closest to 0, defaults to 0
phi0 float the angle \\(\\phi_0\\) at which the radial position is closest to 0, defaults to 0
Source code inhamilflow/models/kepler_problem/model.py class Kepler2DFI(BaseModel):\n r\"\"\"The first integrals for a Kepler problem.\n\n :cvar ene: the energy $E$\n :cvar angular_mom: the angular momentum $l$\n :cvar t0: the time $t_0$ at which the radial position is closest to 0, defaults to 0\n :cvar phi0: the angle $\\phi_0$ at which the radial position is closest to 0, defaults to 0\n \"\"\"\n\n ene: float = Field()\n angular_mom: float = Field()\n t0: float = Field(default=0)\n phi0: float = Field(ge=0, lt=2 * math.pi, default=0)\n\n # TODO process angular momentum = 0\n @field_validator(\"angular_mom\")\n @classmethod\n def _angular_mom_non_zero(cls, v: float) -> float:\n if v == 0:\n msg = \"Only non-zero angular momenta are supported\"\n raise NotImplementedError(msg)\n return v\nKepler2DSystem","text":" Bases: BaseModel
Definition of the Kepler problem.
Potential:
\\[ V(r) = - \\frac{\\alpha}{r}. \\]For reference, if an object is orbiting our Sun, the constant \\(\\alpha = G M_{\\odot} ~ 1.327 \\times 10^{20} \\mathrm{m}^3/\\mathrm{s}^2\\) in SI, which is also called 1 TCB, or 1 solar mass parameter. For computational stability, we recommend using TCB as the unit instead of the large SI values.
Units
When specifying the parameters of the system, be ware of the consistency of the units.
Attributes:
Name Type Descriptionalpha float the proportional constant of the potential energy.
mass float the mass of the orbiting object
Source code inhamilflow/models/kepler_problem/model.py class Kepler2DSystem(BaseModel):\n r\"\"\"Definition of the Kepler problem.\n\n Potential:\n\n $$\n V(r) = - \\frac{\\alpha}{r}.\n $$\n\n For reference, if an object is orbiting our Sun, the constant\n $\\alpha = G M_{\\odot} ~ 1.327 \\times 10^{20} \\mathrm{m}^3/\\mathrm{s}^2$ in SI,\n which is also called 1 TCB, or 1 solar mass parameter. For computational stability, we recommend using\n TCB as the unit instead of the large SI values.\n\n !!! note \"Units\"\n\n When specifying the parameters of the system, be ware of the consistency of the units.\n\n :cvar alpha: the proportional constant of the potential energy.\n :cvar mass: the mass of the orbiting object\n \"\"\"\n\n # TODO add repulsive alpha < 0\n alpha: float = Field(gt=0, default=1.0)\n mass: float = Field(gt=0, default=1.0)\nNumerics for the Kepler problem.
"},{"location":"references/models/kepler_problem/numerics/#hamilflow.models.kepler_problem.numerics.nsolve_u_from_tau_bisect","title":"nsolve_u_from_tau_bisect(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 0, ecc != 1
requiredtau ArrayLike scaled time
requiredReturns:
Type Descriptionlist[OptimizeResult] numeric OptimizeResult from scipy
Source code inhamilflow/models/kepler_problem/numerics.py def nsolve_u_from_tau_bisect(ecc: float, tau: \"npt.ArrayLike\") -> list[OptimizeResult]:\n \"\"\"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.\n\n :param ecc: eccentricity, ecc > 0, ecc != 1\n :param tau: scaled time\n :return: numeric OptimizeResult from scipy\n \"\"\"\n tau_s = np.asarray(tau).reshape(-1)\n if 0 < ecc < 1:\n tau_of_u = tau_of_u_elliptic\n elif ecc > 1:\n tau_of_u = tau_of_u_hyperbolic\n else:\n msg = f\"Expect ecc > 0, ecc != 1, got {ecc}\"\n raise ValueError(msg)\n\n def f(u: float, tau: float) -> np.float64:\n return (\n np.finfo(np.float64).max if u == -1 else np.float64(tau_of_u(ecc, u) - tau)\n )\n\n return [toms748(f, max(-1, -ecc), ecc, (ta,), 2, full_output=True) for ta in tau_s]\nnsolve_u_from_tau_newton(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 0, ecc != 1
requiredtau ArrayLike scaled time
requiredReturns:
Type DescriptionOptimizeResult numeric OptimizeResult from scipy
Raises:
Type DescriptionValueError when ecc is invalid
hamilflow/models/kepler_problem/numerics.py def nsolve_u_from_tau_newton(ecc: float, tau: \"npt.ArrayLike\") -> OptimizeResult:\n \"\"\"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.\n\n :param ecc: eccentricity, ecc > 0, ecc != 1\n :param tau: scaled time\n :raises ValueError: when `ecc` is invalid\n :return: numeric OptimizeResult from scipy\n \"\"\"\n tau = np.asarray(tau)\n if 0 < ecc < 1:\n tau_of_u = tau_of_u_elliptic\n u0 = _u0_elliptic(ecc, tau)\n elif ecc > 1:\n tau_of_u = tau_of_u_hyperbolic\n u0 = _u0_hyperbolic(ecc, tau)\n else:\n msg = f\"Expect ecc > 0, ecc != 1, got {ecc}\"\n raise ValueError(msg)\n\n def f(u: float, tau: \"npt.NDArray[np.float64]\") -> \"npt.NDArray[np.float64]\":\n return tau_of_u(ecc, u) - tau\n\n def fprime(\n u: float,\n tau: \"npt.ArrayLike\", # noqa: ARG001\n ) -> \"npt.NDArray[np.float64]\":\n return tau_of_u_prime(ecc, u)\n\n def fprime2(\n u: float,\n tau: \"npt.ArrayLike\", # noqa: ARG001\n ) -> \"npt.NDArray[np.float64]\":\n return tau_of_u_prime2(ecc, u)\n\n return newton(f, u0, fprime, (tau,), fprime2=fprime2, full_output=True, disp=False)\nu_of_tau(ecc, tau, method='bisect')","text":"Calculate the convenient radial inverse u from tau, using numeric methods.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0
requiredtau ArrayLike scaled time
requiredmethod Literal['bisect', 'newton'] \"newton\" or \"bisect\" numeric methods, defaults to \"bisect\"
'bisect' Returns:
Type DescriptionNDArray[float64] convenient radial inverse u
Raises:
Type DescriptionValueError when method is invalid
ValueError when ecc is invalid
hamilflow/models/kepler_problem/numerics.py def u_of_tau(\n ecc: float,\n tau: \"npt.ArrayLike\",\n method: Literal[\"bisect\", \"newton\"] = \"bisect\",\n) -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the convenient radial inverse u from tau, using numeric methods.\n\n :param ecc: eccentricity, ecc >= 0\n :param tau: scaled time\n :param method: \"newton\" or \"bisect\" numeric methods, defaults to \"bisect\"\n :raises ValueError: when `method` is invalid\n :raises ValueError: when `ecc` is invalid\n :return: convenient radial inverse u\n \"\"\"\n tau = np.asarray(tau)\n if ecc == 0:\n return np.zeros(tau.shape)\n elif ecc == 1:\n return esolve_u_from_tau_parabolic(ecc, tau)\n elif ecc > 0:\n match method:\n case \"bisect\":\n return np.array([s[0] for s in nsolve_u_from_tau_bisect(ecc, tau)])\n case \"newton\":\n return nsolve_u_from_tau_newton(ecc, tau)[0]\n case _:\n msg = f\"Expect 'bisect' or 'newton', got {method}\"\n raise ValueError(msg)\n else:\n msg = f\"Expect ecc >= 0, got {ecc}\"\n raise ValueError(msg)\nWe provide some tutorials to help you get started.
For easier reading, we following the following notations whenever possible.
Notation Description \\(\\mathcal L\\) Lagrangian \\(L\\) Angular Momentum \\(V\\) Potential"},{"location":"tutorials/brownian_motion/","title":"Brownian Motion","text":"In\u00a0[1]: Copied!import plotly.express as px\n\nfrom hamilflow.models.brownian_motion import BrownianMotion\nimport plotly.express as px from hamilflow.models.brownian_motion import BrownianMotion In\u00a0[2]: Copied!
bm_1d = BrownianMotion(\n system={\n \"sigma\": 1,\n \"delta_t\": 1,\n },\n initial_condition={\"x0\": 0},\n)\n bm_1d = BrownianMotion( system={ \"sigma\": 1, \"delta_t\": 1, }, initial_condition={\"x0\": 0}, ) Call the model to generate 1000 steps.
In\u00a0[3]: Copied!df_1d = bm_1d.generate_from(n_steps=1000)\ndf_1d = bm_1d.generate_from(n_steps=1000) In\u00a0[4]: Copied!
px.line(df_1d, x=\"t\", y=\"x_0\")\npx.line(df_1d, x=\"t\", y=\"x_0\") In\u00a0[5]: Copied!
bm_2d = BrownianMotion(\n system={\n \"sigma\": 1,\n \"delta_t\": 1,\n },\n initial_condition={\"x0\": [0, 0]},\n)\n bm_2d = BrownianMotion( system={ \"sigma\": 1, \"delta_t\": 1, }, initial_condition={\"x0\": [0, 0]}, ) We call the model to generate 1000 steps.
In\u00a0[6]: Copied!df_2d = bm_2d.generate_from(n_steps=500)\ndf_2d = bm_2d.generate_from(n_steps=500) In\u00a0[7]: Copied!
(\n px.scatter(df_2d, x=\"x_0\", y=\"x_1\", color=\"t\")\n .update_traces(\n mode=\"lines+markers\",\n marker={\n \"size\": 2.5,\n },\n line={\"width\": 1},\n )\n .update_yaxes(\n scaleanchor=\"x\",\n scaleratio=1,\n )\n)\n ( px.scatter(df_2d, x=\"x_0\", y=\"x_1\", color=\"t\") .update_traces( mode=\"lines+markers\", marker={ \"size\": 2.5, }, line={\"width\": 1}, ) .update_yaxes( scaleanchor=\"x\", scaleratio=1, ) ) In\u00a0[\u00a0]: Copied! \n"},{"location":"tutorials/brownian_motion/#brownian-motion","title":"Brownian Motion\u00b6","text":"
Brownian motion describes motion of small particles with stochastic forces applied to them. The math of Brownian motion can be modeled with Wiener process. In this tutorial, we take a simple form of the model and treat the stochastic forces as Gaussian.
For consistency, we always use $\\mathbf x$ for displacement, and $t$ for steps. The model we are using is
$$ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} $$
Read more:
hamilflow implemented a Brownian motion model called BrownianMotion.
Our BrownianMotion model calculates the dimension of the space based on the dimension of the initial condition $x_0$. To create a 2D Brownian motion model, we need the initial condition to be length 2.
import math\n\nimport numpy as np\nfrom plotly import express as px\n\nfrom hamilflow.models.harmonic_oscillator import ComplexSimpleHarmonicOscillator\nimport math import numpy as np from plotly import express as px from hamilflow.models.harmonic_oscillator import ComplexSimpleHarmonicOscillator In\u00a0[2]: Copied!
t = np.linspace(0, 3, 257)\nsystem_specs = {\"omega\": 2 * math.pi}\n t = np.linspace(0, 3, 257) system_specs = {\"omega\": 2 * math.pi} In\u00a0[3]: Copied! csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (1, 0)})\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (1, 0)}) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, ) In\u00a0[4]: Copied! csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (0, 1)})\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (0, 1)}) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, ) In\u00a0[5]: Copied! csho = ComplexSimpleHarmonicOscillator(\n system_specs,\n initial_condition={\"x0\": (math.cos(math.pi / 12), math.sin(math.pi / 12))},\n)\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator( system_specs, initial_condition={\"x0\": (math.cos(math.pi / 12), math.sin(math.pi / 12))}, ) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, )"},{"location":"tutorials/complex_harmonic_oscillator/#complex-harmonic-oscillator","title":"Complex Harmonic Oscillator\u00b6","text":"In this tutorial, we show a few interesting examples of the complex simple harmonic oscillator.
"},{"location":"tutorials/complex_harmonic_oscillator/#basic-setups","title":"Basic setups\u00b6","text":""},{"location":"tutorials/complex_harmonic_oscillator/#positive-frequency-circular-polarised-mode","title":"Positive-frequency, circular-polarised mode\u00b6","text":"Also known as the left-rotating mode.
"},{"location":"tutorials/complex_harmonic_oscillator/#negative-frequency-circular-polarised-mode","title":"Negative-frequency, circular-polarised mode\u00b6","text":"Also known as the right-rotating mode.
"},{"location":"tutorials/complex_harmonic_oscillator/#positive-frequency-elliptic-polarised-mode","title":"Positive-frequency, elliptic-polarised mode\u00b6","text":""},{"location":"tutorials/complex_harmonic_oscillator/#end-of-notebook","title":"End of Notebook\u00b6","text":""},{"location":"tutorials/harmonic_oscillator/","title":"Harmonic Oscillators","text":"In\u00a0[1]: Copied!import pandas as pd\nimport plotly.express as px\n\nfrom hamilflow.models.harmonic_oscillator import (\n DampedHarmonicOscillator,\n SimpleHarmonicOscillator,\n)\nimport pandas as pd import plotly.express as px from hamilflow.models.harmonic_oscillator import ( DampedHarmonicOscillator, SimpleHarmonicOscillator, ) In\u00a0[2]: Copied!
n_periods = 3\nn_samples_per_period = 200\nn_periods = 3 n_samples_per_period = 200 In\u00a0[3]: Copied!
sho_omega = 0.5\n\nsho = SimpleHarmonicOscillator(system={\"omega\": sho_omega})\n sho_omega = 0.5 sho = SimpleHarmonicOscillator(system={\"omega\": sho_omega}) In\u00a0[4]: Copied! df_sho = sho(n_periods=n_periods, n_samples_per_period=n_samples_per_period)\ndf_sho.head()\ndf_sho = sho(n_periods=n_periods, n_samples_per_period=n_samples_per_period) df_sho.head() Out[4]: t x 0 0.000000 1.000000 1 0.062832 0.999507 2 0.125664 0.998027 3 0.188496 0.995562 4 0.251327 0.992115 In\u00a0[5]: Copied!
px.line(\n df_sho,\n x=\"t\",\n y=\"x\",\n title=rf\"Simple Harmonic Oscillator (omega = {sho_omega})\",\n labels={\n \"x\": r\"Displacement $x(t)$\",\n \"t\": r\"$t$\",\n },\n)\n px.line( df_sho, x=\"t\", y=\"x\", title=rf\"Simple Harmonic Oscillator (omega = {sho_omega})\", labels={ \"x\": r\"Displacement $x(t)$\", \"t\": r\"$t$\", }, ) In\u00a0[6]: Copied! dho_systems = {\n \"Underdamped\": {\"omega\": 0.5, \"zeta\": 0.2},\n \"Critical Damped\": {\"omega\": 0.5, \"zeta\": 1},\n \"Overdamped\": {\n \"omega\": 0.5,\n \"zeta\": 1.2,\n },\n}\n\ndfs_dho = []\n\nfor s_name, s in dho_systems.items():\n dfs_dho.append(\n DampedHarmonicOscillator(system=s)(\n n_periods=n_periods,\n n_samples_per_period=n_samples_per_period,\n ).assign(\n system=rf\"{s_name} (omega = {s.get('omega')}, zeta = {s.get('zeta')})\",\n ),\n )\n\nfig = px.line(\n pd.concat(dfs_dho),\n x=\"t\",\n y=\"x\",\n color=\"system\",\n title=r\"Damped Harmonic Oscillator\",\n labels={\n \"x\": r\"Displacement $x(t)$\",\n \"t\": r\"$t$\",\n },\n)\nfig.update_layout(legend={\"yanchor\": \"top\", \"y\": -0.2, \"xanchor\": \"left\", \"x\": 0})\n dho_systems = { \"Underdamped\": {\"omega\": 0.5, \"zeta\": 0.2}, \"Critical Damped\": {\"omega\": 0.5, \"zeta\": 1}, \"Overdamped\": { \"omega\": 0.5, \"zeta\": 1.2, }, } dfs_dho = [] for s_name, s in dho_systems.items(): dfs_dho.append( DampedHarmonicOscillator(system=s)( n_periods=n_periods, n_samples_per_period=n_samples_per_period, ).assign( system=rf\"{s_name} (omega = {s.get('omega')}, zeta = {s.get('zeta')})\", ), ) fig = px.line( pd.concat(dfs_dho), x=\"t\", y=\"x\", color=\"system\", title=r\"Damped Harmonic Oscillator\", labels={ \"x\": r\"Displacement $x(t)$\", \"t\": r\"$t$\", }, ) fig.update_layout(legend={\"yanchor\": \"top\", \"y\": -0.2, \"xanchor\": \"left\", \"x\": 0}) In\u00a0[\u00a0]: Copied! \n"},{"location":"tutorials/harmonic_oscillator/#harmonic-oscillators","title":"Harmonic Oscillators\u00b6","text":"
In this tutorial, we demo how to generate data of harmonic oscillators.
"},{"location":"tutorials/harmonic_oscillator/#simple-harmonic-oscillator","title":"Simple Harmonic Oscillator\u00b6","text":"For an simple harmonic oscillator, the action of a simple harmonic oscillator is
$$S_L[x] = \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} m \\dot x^2 - \\frac{1}{2} m \\omega^2 x^2 \\right\\}\\,,$$
where the least action principle leads to the following equation of motion,
$$ \\ddot x + \\omega^2 x = 0\\,. $$
A simple harmonic oscillator is a periodic motion.
"},{"location":"tutorials/harmonic_oscillator/#damped-harmonic-oscillator","title":"Damped Harmonic Oscillator\u00b6","text":"A damped harmonic oscillator is a simple harmonic oscillator with damping force that is proportional to its velocity,
$$ \\ddot x + \\omega^2 x = - 2\\xi\\omega \\dot x\\,. $$
In this section, we demonstrate three scenarios of a damped harmonic oscillator.
"},{"location":"tutorials/harmonic_oscillator_chain/","title":"Harmonic oscillator circle","text":"In\u00a0[1]: Copied!import math\n\nimport numpy as np\nfrom plotly import express as px\n\nfrom hamilflow.models.harmonic_oscillator_chain import HarmonicOscillatorsChain\nimport math import numpy as np from plotly import express as px from hamilflow.models.harmonic_oscillator_chain import HarmonicOscillatorsChain In\u00a0[2]: Copied!
t = np.linspace(0, 3, 257)\nomega = 2 * math.pi\npattern_x = r\"x\\d+\"\nt = np.linspace(0, 3, 257) omega = 2 * math.pi pattern_x = r\"x\\d+\" In\u00a0[3]: Copied!
ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] # type: ignore [index]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] # type: ignore [index] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[4]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(2 * math.pi, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(2 * math.pi, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[5]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 0)}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 4\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 0)}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 4 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[6]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[7]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, False)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, False) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[8]: Copied! ics = [{\"x0\": 0, \"v0\": 2}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, False)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 2}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, False) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})"},{"location":"tutorials/harmonic_oscillator_chain/#harmonic-oscillator-circle","title":"Harmonic oscillator circle\u00b6","text":"In this tutorial, we show a few interesting examples of the harmonic oscillator circle.
"},{"location":"tutorials/harmonic_oscillator_chain/#basic-setups","title":"Basic setups\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-right-moving-wave","title":"Fundamental right-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-left-moving-wave","title":"Fundamental left-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#faster-right-moving-wave","title":"Faster right-moving wave\u00b6","text":"Also known as the first harmonic.
"},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-stationary-wave-odd-dof","title":"Fundamental stationary wave, odd dof\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-stationary-wave-even-dof","title":"Fundamental stationary wave, even dof\u00b6","text":"There are stationary nodes at $i = 3, 9$.
"},{"location":"tutorials/harmonic_oscillator_chain/#linearly-moving-chain-and-fundamental-right-moving-wave","title":"Linearly moving chain and fundamental right-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#mathematical-physical-discription","title":"Mathematical-physical discription\u00b6","text":"A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action $$S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\\\}\\\\,,$$ where $x_N \\coloneqq x_0$.
This system can be solved in terms of travelling waves, obtained by discrete Fourier transform.
We can complexify the system $$S_L[x_i] = S_L[x_i, \\phi_j] \\equiv S_L[X^\\ast_i, X_j] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot X^\\ast_i \\delta_{ij} \\dot X_j - \\frac{1}{2}m X^\\ast_i A_{ij} X_j\\right\\\\}\\\\,,$$ where $A_{ij} / \\omega^2$ is equal to $(-2)$ if $i=j$, $1$ if $|i-j|=1$ or $|i-j|=N$, and $0$ otherwise; $X_i \\coloneqq x_i \\mathbb{e}^{-\\phi_i}$, $X^\\ast_i \\coloneqq x_i \\mathbb{e}^{+\\phi_i}$.
$A_{ij}$ can be diagonalised by the inverse discrete Fourier transform $$X_i = (F^{-1})_{ik} Y_k = \\frac{1}{\\sqrt{N}}\\sum_k \\mathbb{e}^{i \\frac{2\\mathbb{\\pi}}{N} k\\mathbb{i}} Y_k\\\\,.$$
Calculating gives $$S_L[X^\\ast_i, X_j] = S_L[Y^\\ast_i, Y_j] = \\sum_{k=0}^{N-1} \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot Y^\\ast_k \\dot Y_k - \\frac{1}{2}m \\omega^2\\cdot4\\sin^2\\frac{2\\mathbb{\\pi}k}{N} Y^\\ast_k Y_k\\right\\\\}\\\\,.$$ We can arrive at an action for the Fourier modes $$S_L[Y, Y^*] = \\sum_{k=0}^{N-1} \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot Y_k^* Y_k - \\frac{1}{2}m \\omega^2\\cdot4\\sin^2\\frac{2\\mathbb{\\pi}k}{N} Y_k^* Y_k\\right\\\\}\\\\,.$$
The origional system can then be solved by $N$ complex oscillators.
Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy $Y_k = Y^*_{-k \\mod N}$, see Wikipedia.
"},{"location":"tutorials/harmonic_oscillator_chain/#end-of-notebook","title":"End of Notebook\u00b6","text":""},{"location":"tutorials/kepler_problem/","title":"Motion in a Central Field","text":"In\u00a0[\u00a0]: Copied! In\u00a0[1]: Copied!# system = {\"alpha\": 1.0, \"mass\": 1.0}\n# ic = {\"r_0\": 1.0, \"phi_0\": 0.0, \"drdt_0\": 0.0, \"dphidt_0\": 0.1}\n\n# t = np.linspace(0, 1, 11)\n # system = {\"alpha\": 1.0, \"mass\": 1.0} # ic = {\"r_0\": 1.0, \"phi_0\": 0.0, \"drdt_0\": 0.0, \"dphidt_0\": 0.1} # t = np.linspace(0, 1, 11) In this section, we briefly discuss the derivations of motion in a central field. Please refer to Mechanics: Vol 1 by Landau and Lifshitz for more details[^1].
The Lagrangian for an object in a central field is
$$ \\mathcal L = \\frac{1}{2} m ({\\dot r}^2 + r^2 {\\dot \\phi}^2) - V(r), $$
where $r$ and $\\phi$ are the polar coordinates, $m$ is the mass of the object, and $V(r)$ is the potential energy. The equations of motion are
$$ \\begin{align} \\frac{\\mathrm d r}{\\mathrm d t} &= \\sqrt{ \\frac{2}{m} (E - V(r)) - \\frac{L^2}{m^2 r^2} } \\\\ \\frac{\\mathrm d \\phi}{\\mathrm d t} &= \\frac{L}{m r^2}, \\end{align} $$
where $E$ is the total energy and $L$ is the angular momentum,
$$ \\begin{align} E =& \\frac{1}{2} m \\left(\\left(\\frac{\\mathrm d r}{ \\mathrm dt} \\right)^2 + r^2 \\left( \\frac{\\mathrm d r}{\\mathrm dt} \\right)^2\\right) + V(r) \\\\ L =& m r^2 \\frac{d\\phi}{dt} \\end{align} $$
Both $E$ and $L$ are conserved. We obtain the coordinates as a function of time by solving the two equations.
For a inverse-square force, the potential energy is
$$ V(r) = - \\frac{\\alpha}{r}, $$
where $\\alpha$ is a constant that specifies the strength of the force. For Newtonian gravity, $\\alpha = G m_0$ with $G$ being the gravitational constant.
First of all, we solve $\\phi(r)$,
$$ \\phi = \\cos^{-1}\\left( \\frac{L/r - m\\alpha/L}{2 m E + m \\alpha^2/L^2} \\right). $$
Define
$$ p = \\frac{L^2}{m\\alpha} $$
and
$$ e = \\sqrt{ 1 + \\frac{2 E L^2}{m \\alpha^2}}, $$
we rewrite the solution $\\phi(r)$ as
$$ r = \\frac{p}{1 + e \\cos{\\phi}}. $$
With the above relationship between $r$ and $\\phi$, and $\\frac{\\mathrm d \\phi}{\\mathrm d} = \\frac{L}{mr^2}$, we find that
$$ \\frac{m\\alpha^2}{L^3} \\mathrm d t = \\frac{1}{(1 + e \\cos{\\phi})^2} \\mathrm d \\phi. $$
The integration on the right hand side depends on the domain of $e$.
$$ \\int\\frac{1}{(1 + e \\cos{\\phi})^2} \\mathrm d \\phi = \\begin{cases} \\frac{1}{(1 - e^2)^{3/2}} \\left( 2 \\tan^{-1} \\sqrt{\\frac{1 - e}{1 + e}} \\tan\\frac{\\phi}{2} - \\frac{e\\sqrt{1 - e^2}\\sin\\phi}{1 + e\\cos\\phi} \\right), & \\text{if } e<1 \\\\ \\frac{1}{2}\\tan{\\frac{\\phi}{2}} + \\frac{1}{6}\\tan^3{\\frac{\\phi}{2}}, & \\text{if } e=1 \\\\ \\frac{1}{(e^2 - 1)^{3/2}}\\left( \\frac{e\\sqrt{e^2-1}\\sin\\phi}{1 + e\\cos\\phi} - \\ln \\left( \\frac{\\sqrt{1 + e} + \\sqrt{e-1}\\tan\\frac{\\phi}{2}}{\\sqrt{1 + e} - \\sqrt{e-1}\\tan\\frac{\\phi}{2} } \\right) \\right), & \\text{if } e> 1. \\end{cases} $$
The value of $t(\\phi)$ is easily obtained from the above formulae.
There exists many numerical methods to solve the Kepler orbits as functions of time, $r(t)$ and $\\phi(t)$. For our use case of the solutions, we choose to integrate the equation of motion directly.
References:
In this tutorial, we generate data for objects moving in a central field.
"},{"location":"tutorials/kepler_problem/#usage","title":"Usage\u00b6","text":""},{"location":"tutorials/kepler_problem/#formalism","title":"Formalism\u00b6","text":""},{"location":"tutorials/pendulum/","title":"Pendulum","text":"In\u00a0[1]: Copied!import math\n\nimport plotly.express as px\n\nfrom hamilflow.models.pendulum import Pendulum\nimport math import plotly.express as px from hamilflow.models.pendulum import Pendulum In\u00a0[2]: Copied!
omega0 = 2 * math.pi\ntheta0 = math.pi / 3\n\nn_periods = 2**2\nn_samples_per_period = 2**8\nomega0 = 2 * math.pi theta0 = math.pi / 3 n_periods = 2**2 n_samples_per_period = 2**8 In\u00a0[3]: Copied!
pen = Pendulum(system=omega0, initial_condition=theta0)\npen = Pendulum(system=omega0, initial_condition=theta0) In\u00a0[4]: Copied!
df_pen = pen.generate_from(\n n_periods=n_periods,\n n_samples_per_period=n_samples_per_period,\n)\ndf_pen.head()\ndf_pen = pen.generate_from( n_periods=n_periods, n_samples_per_period=n_samples_per_period, ) df_pen.head() Out[4]: t x u 0 0.000000 1.047198 1.570796 1 0.004192 1.046897 1.593608 2 0.008384 1.045996 1.616424 3 0.012576 1.044494 1.639247 4 0.016768 1.042393 1.662082 In\u00a0[5]: Copied!
df_pen.describe()\ndf_pen.describe() Out[5]: t x u count 1024.000000 1.024000e+03 1024.000000 mean 2.144268 2.320193e-17 14.124895 std 1.239809 7.453532e-01 7.261249 min 0.000000 -1.047198e+00 1.570796 25% 1.072134 -7.494689e-01 7.848279 50% 2.144268 -3.535251e-16 14.125761 75% 3.216402 7.494689e-01 20.403244 max 4.288536 1.047198e+00 26.680726 In\u00a0[6]: Copied!
px.line(\n df_pen,\n x=\"t\",\n y=\"x\",\n title=rf\"Simple Harmonic Oscillator ($\\omega_0 = {omega0:.4f})$\",\n labels={\"x\": r\"Angle $\\theta(t)$\", \"t\": r\"Time $t$\"},\n)\n px.line( df_pen, x=\"t\", y=\"x\", title=rf\"Simple Harmonic Oscillator ($\\omega_0 = {omega0:.4f})$\", labels={\"x\": r\"Angle $\\theta(t)$\", \"t\": r\"Time $t$\"}, )"},{"location":"tutorials/pendulum/#pendulum","title":"Pendulum\u00b6","text":"In this tutorial, we demonstrate how to generate data of a pendulum, and introduce the mathematics of a pendulum.
"},{"location":"tutorials/pendulum/#constants","title":"Constants\u00b6","text":""},{"location":"tutorials/pendulum/#a-pendulum","title":"A pendulum\u00b6","text":""},{"location":"tutorials/pendulum/#data","title":"Data\u00b6","text":""},{"location":"tutorials/pendulum/#plot","title":"Plot\u00b6","text":""},{"location":"tutorials/pendulum/#todo","title":"TODO\u00b6","text":"We describe a generic pendulum system by the Lagrangian action $$ S_L[\\theta] \\equiv \\int_{t_0}^{t_1} \\mathbb{d}t\\,L(\\theta, \\dot\\theta) \\eqqcolon I \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\}\\,, $$ where $L$ is the Lagrangian; $\\theta$ is the angle from the vertical to the pendulum as the generalised position; $I$ is the inertia parameter, $\\omega_0$ the frequency parameter, and we also call $U \\coloneqq I\\omega_0^2$ the potential parameter.
This setup contains both the single and the physical pendula. For a single pendulum, $$ I = m l^2\\,,\\qquad U = mgl\\,, $$ where $m$ is the mass of the pendulum, $l$ is the length of the rod or cord, and $g$ is the gravitational acceleration.
"},{"location":"tutorials/pendulum/#integral-of-motion","title":"Integral of motion\u00b6","text":"The Lagrangian action does not contain time $t$ explicitly. As a result, the system is invariant under a variation of time, or $\\mathbb{\\delta}S / \\mathbb{\\delta}{t} = 0$. This gives an integral of motion $$ \\dot\\theta\\frac{\\partial L}{\\partial \\dot\\theta} - L \\equiv E \\eqqcolon I \\omega_0^2 \\cos\\theta_0\\,, $$ where $\\theta_0$ is the initial angle.
Substitution gives $$ \\left(\\frac{\\mathbb{d}t}{\\mathbb{d}\\theta}\\right)^2 = \\frac{1}{2\\omega_0^2} \\frac{1}{\\cos\\theta - \\cos\\theta_0}\\,. $$
"},{"location":"tutorials/pendulum/#coordinate-transformation","title":"Coordinate transformation\u00b6","text":"For convenience, introduce the coordinate $u$ and the parameter $k$ $$ \\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{k}\\,,\\qquad k \\coloneqq \\sin\\frac{\\theta_0}{2} \\in [-1, 1]\\,. $$ One arrives at $$ \\left(\\frac{\\mathbb{d}t}{\\mathbb{d}u}\\right)^2 = \\frac{1}{\\omega_0^2} \\frac{1}{1-k^2\\sin^2 u}\\,. $$ The square root of the second factor on the right-hand side makes an elliptic integral.
"},{"location":"tutorials/pendulum/#end-of-notebook","title":"End of Notebook\u00b6","text":""}]} \ No newline at end of file +{"config":{"lang":["en"],"separator":"[\\s\\-]+","pipeline":["stopWordFilter"]},"docs":[{"location":"","title":"Documentation forhamilflow","text":"Generating dataset for physical systems.
"},{"location":"#development","title":"Development","text":"We use poetry to manage our python environment. Use
poetry install\nto install the requirements. Or run
poetry install --with test\nto install the base environment and test environment for development.
If this is the first you clone the repo and committing code, run
pre-commit install\nfirst.
"},{"location":"changelog/","title":"hamilflow Changelog","text":""},{"location":"changelog/#010-2024-03-21","title":"0.1.0 (2024-03-21)","text":""},{"location":"changelog/#feat","title":"Feat","text":"Setting up repository.
"},{"location":"references/","title":"References","text":""},{"location":"references/maths/trigonometrics/","title":"Trigonometrics","text":"Trigonometric functions.
"},{"location":"references/maths/trigonometrics/#hamilflow.maths.trigonometrics.acos_with_shift","title":"acos_with_shift(x, shift=None)","text":"Arccos with shift.
Source code inhamilflow/maths/trigonometrics.py def acos_with_shift(\n x: \"Collection[float] | npt.ArrayLike\",\n shift: \"Collection[float] | npt.ArrayLike | None\" = None,\n) -> \"npt.ArrayLike\":\n \"\"\"Arccos with shift.\"\"\"\n x = np.asarray(x)\n value = np.arccos(x)\n shift = np.asarray(shift)\n period_shift = (div := np.floor(shift)) * 2 * np.pi\n remainder = shift - div\n value = np.where(remainder <= 0.5, value, 2 * np.pi - value) # noqa: PLR2004\n\n return period_shift + value\nMain module for Brownian motion.
"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion","title":"BrownianMotion","text":"Brownian motion describes motion of small particles with stochastic forces applied to them.
The math of Brownian motion can be modeled with Wiener process.
For consistency, we always use \\(\\mathbf x\\) for displacement, and \\(t\\) for steps. The model we are using is
\\[ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} \\]"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion--references","title":"References","text":"1D Brownian Motion
The dimsion of our Brownian motion is specified by the dimension of the initial condition.
To simulate a 1D Browian motion, we define the system and initial condition:
system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n}\n\ninitial_condition = {\n \"x0\": 0\n}\nThe Brownian motion can be simulated using
bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\nbm(n_steps=100)\n2D Brownian Motion
To simulate a 2D Browian motion,
system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n}\n\ninitial_condition = {\n \"x0\": [0, 0]\n}\n\nbm = BrownianMotion(system=system, initial_condition=initial_condition)\n\nbm(n_steps=100)\nParameters:
Name Type Description Defaultsystem Mapping[str, float] the Brownian motion system definition
requiredinitial_condition Mapping[str, Sequence[float] | ArrayLike] | None the initial condition for the simulation
None Source code in hamilflow/models/brownian_motion.py class BrownianMotion:\n r\"\"\"Brownian motion describes motion of small particles with stochastic forces applied to them.\n\n The math of Brownian motion can be modeled\n with Wiener process.\n\n For consistency, we always use\n $\\mathbf x$ for displacement, and\n $t$ for steps. The model we are using is\n\n $$\n \\begin{align}\n \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) +\n \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t})\n \\end{align}\n $$\n\n References\n ----------\n 1. Brownian motion and random walks. [cited 13 Mar 2024].\n Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html\n 2. Contributors to Wikimedia projects. Brownian motion.\n In: Wikipedia [Internet]. 22 Jan 2024 [cited 13 Mar 2024].\n Available: https://en.wikipedia.org/wiki/Brownian_motion\n\n\n !!! example \"1D Brownian Motion\"\n\n The dimsion of our Brownian motion is specified by\n the dimension of the initial condition.\n\n To simulate a 1D Browian motion, we define the system and initial condition:\n\n ```python\n system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n }\n\n initial_condition = {\n \"x0\": 0\n }\n ```\n\n The Brownian motion can be simulated using\n\n ```python\n bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\n bm(n_steps=100)\n ```\n\n !!! example \"2D Brownian Motion\"\n\n To simulate a 2D Browian motion,\n\n ```python\n system = {\n \"sigma\": 1,\n \"delta_t\": 1,\n }\n\n initial_condition = {\n \"x0\": [0, 0]\n }\n\n bm = BrownianMotion(system=system, initial_condition=initial_condition)\n\n bm(n_steps=100)\n ```\n\n :param system: the Brownian motion system definition\n :param initial_condition: the initial condition for the simulation\n\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: (\n Mapping[str, \"Sequence[float] | npt.ArrayLike\"] | None\n ) = None,\n ) -> None:\n initial_condition = initial_condition or {}\n self.system = BrownianMotionSystem.model_validate(system)\n self.initial_condition = BrownianMotionIC.model_validate(initial_condition)\n\n @property\n def dim(self) -> int:\n \"\"\"Dimension of the Brownian motion.\"\"\"\n return np.asarray(self.initial_condition.x0).size\n\n @property\n def _axis_names(self) -> list[str]:\n return [f\"x_{i}\" for i in range(self.dim)]\n\n def _trajectory(self, n_new_steps: int, seed: int) -> \"npt.NDArray[np.float64]\":\n \"\"\"Give the trajectory of the particle.\n\n We first compute the delta displacement in each step.\n With the displacement at each step, we perform a cumsum\n including the initial coordinate to get the displacement at each step.\n\n :param n_new_steps: number of new steps to simulate, excluding the initial step.\n :param seed: seed for the random generator.\n \"\"\"\n step_history = sp.stats.norm.rvs(\n size=(n_new_steps, self.dim) if self.dim > 1 else n_new_steps,\n scale=self.system.gaussian_scale,\n random_state=np.random.RandomState(seed=seed),\n )\n\n step_history = np.concatenate(\n (np.expand_dims(self.initial_condition.x0, axis=0), step_history),\n )\n\n trajectory = np.cumsum(step_history, axis=0)\n\n return trajectory\n\n def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:\n \"\"\"Generate data from a set of interpretable params for this model.\n\n :param n_steps: total number of steps to be simulated, including the inital step.\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n time_steps = np.arange(0, n_steps) * self.system.delta_t\n\n return self(t=time_steps, seed=seed)\n\n def __call__(self, t: TypeTime, seed: int = 42) -> pd.DataFrame:\n \"\"\"Simulate the coordinates of the particle.\n\n :param t: the time sequence to be used to generate data, 1-D array like\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n n_steps = np.array(t).size\n trajectory = self._trajectory(n_new_steps=n_steps - 1, seed=seed)\n\n df = pd.DataFrame(trajectory, columns=self._axis_names)\n df[\"t\"] = t\n\n return df\ndim: int property","text":"Dimension of the Brownian motion.
"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotion.__call__","title":"__call__(t, seed=42)","text":"Simulate the coordinates of the particle.
Parameters:
Name Type Description Defaultt TypeTime the time sequence to be used to generate data, 1-D array like
requiredseed int random generator seed for the stochastic process. Use it to reproduce results.
42 Source code in hamilflow/models/brownian_motion.py def __call__(self, t: TypeTime, seed: int = 42) -> pd.DataFrame:\n \"\"\"Simulate the coordinates of the particle.\n\n :param t: the time sequence to be used to generate data, 1-D array like\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n n_steps = np.array(t).size\n trajectory = self._trajectory(n_new_steps=n_steps - 1, seed=seed)\n\n df = pd.DataFrame(trajectory, columns=self._axis_names)\n df[\"t\"] = t\n\n return df\ngenerate_from(n_steps, seed=42)","text":"Generate data from a set of interpretable params for this model.
Parameters:
Name Type Description Defaultn_steps int total number of steps to be simulated, including the inital step.
requiredseed int random generator seed for the stochastic process. Use it to reproduce results.
42 Source code in hamilflow/models/brownian_motion.py def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:\n \"\"\"Generate data from a set of interpretable params for this model.\n\n :param n_steps: total number of steps to be simulated, including the inital step.\n :param seed: random generator seed for the stochastic process.\n Use it to reproduce results.\n \"\"\"\n time_steps = np.arange(0, n_steps) * self.system.delta_t\n\n return self(t=time_steps, seed=seed)\nBrownianMotionIC","text":" Bases: BaseModel
The initial condition for a Brownian motion.
Attributes:
Name Type Descriptionx0 float | Sequence[float] initial displacement of the particle, the diminsion of this initial condition determines the dimension of the model too.
Source code inhamilflow/models/brownian_motion.py class BrownianMotionIC(BaseModel):\n \"\"\"The initial condition for a Brownian motion.\n\n :cvar x0: initial displacement of the particle,\n the diminsion of this initial condition determines\n the dimension of the model too.\n \"\"\"\n\n x0: float | Sequence[float] = Field(default=1.0)\n\n @field_validator(\"x0\")\n @classmethod\n def _check_x0_types(cls, v: float | Sequence[float]) -> float | Sequence[float]:\n if not isinstance(v, float | int | Sequence):\n # TODO I do not think this raise can be reached\n msg = f\"Value of x0 should be int/float/list of int/float: {v=}\"\n raise TypeError(msg)\n\n return v\nBrownianMotionSystem","text":" Bases: BaseModel
Definition of the Brownian Motion system.
For consistency, we always use \\(\\mathbf x\\) for displacement, and \\(t\\) for steps. The model we are using is
\\[ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} \\]"},{"location":"references/models/brownian_motion/#hamilflow.models.brownian_motion.BrownianMotionSystem--references","title":"References","text":"Attributes:
Name Type Descriptionsigma float base standard deviation to be used to compute the variance
delta_t float time granunality of the motion
Source code inhamilflow/models/brownian_motion.py class BrownianMotionSystem(BaseModel):\n r\"\"\"Definition of the Brownian Motion system.\n\n For consistency, we always use\n $\\mathbf x$ for displacement, and\n $t$ for steps. The model we are using is\n\n $$\n \\begin{align}\n \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) +\n \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t})\n \\end{align}\n $$\n\n References\n ----------\n 1. Brownian motion and random walks. [cited 13 Mar 2024].\n Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html\n 2. Contributors to Wikimedia projects. Brownian motion.\n In: Wikipedia [Internet]. 22 Jan 2024 [cited 13 Mar 2024].\n Available: https://en.wikipedia.org/wiki/Brownian_motion\n\n :cvar sigma: base standard deviation\n to be used to compute the variance\n :cvar delta_t: time granunality of the motion\n\n \"\"\"\n\n sigma: float = Field(ge=0.0)\n delta_t: float = Field(ge=0.0, default=1.0)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def gaussian_scale(self) -> float:\n \"\"\"The scale (standard deviation) of the Gaussian term in Brownian motion.\"\"\"\n return self.sigma**2 * self.delta_t\ngaussian_scale: float cached property","text":"The scale (standard deviation) of the Gaussian term in Brownian motion.
"},{"location":"references/models/harmonic_oscillator/","title":"Harmonic Oscillator","text":"Main module for undamped and damped hamornic oscillators.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.ComplexSimpleHarmonicOscillator","title":"ComplexSimpleHarmonicOscillator","text":"Generate time series data for a complex simple harmonic oscillator.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the complex harmonic oscillator.
requiredinitial_condition Mapping[str, tuple[float, float]] the initial condition of the complex harmonic oscillator.
required Source code inhamilflow/models/harmonic_oscillator.py class ComplexSimpleHarmonicOscillator:\n r\"\"\"Generate time series data for a complex simple harmonic oscillator.\n\n :param system: all the params that defines the complex harmonic oscillator.\n :param initial_condition: the initial condition of the complex harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, tuple[float, float]],\n ) -> None:\n self.system = HarmonicOscillatorSystem.model_validate(system)\n self.initial_condition = ComplexSimpleHarmonicOscillatorIC.model_validate(\n initial_condition,\n )\n if self.system.type != \"simple\":\n msg = f\"System is not a Simple Harmonic Oscillator: {self.system}\"\n raise ValueError(\n msg,\n )\n\n @cached_property\n def definition(\n self,\n ) -> dict[str, dict[str, float | tuple[float, float]]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n def _z(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to complex simple harmonic oscillators.\n\n $$\n x(t) = x_+ \\exp(-\\mathbb{i} (\\omega t + \\phi_+)) + x_- \\exp(+\\mathbb{i} (\\omega t + \\phi_-)).\n $$\n \"\"\"\n t = np.asarray(t)\n omega = self.system.omega\n x0, phi = self.initial_condition.x0, self.initial_condition.phi\n phases = -omega * t - phi[0], omega * t + phi[1]\n return x0[0] * np.exp(1j * phases[0]) + x0[1] * np.exp(1j * phases[1])\n\n def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time.\n\n :param t: time(s).\n \"\"\"\n t = np.asarray(t)\n data = self._z(t)\n\n return pd.DataFrame({\"t\": t, \"z\": data})\ndefinition: dict[str, dict[str, float | tuple[float, float]]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.ComplexSimpleHarmonicOscillator.__call__","title":"__call__(t)","text":"Generate time series data for the harmonic oscillator.
Returns a list of floats representing the displacement at each time.
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time(s).
required Source code inhamilflow/models/harmonic_oscillator.py def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time.\n\n :param t: time(s).\n \"\"\"\n t = np.asarray(t)\n data = self._z(t)\n\n return pd.DataFrame({\"t\": t, \"z\": data})\nComplexSimpleHarmonicOscillatorIC","text":" Bases: BaseModel
The initial condition for a complex harmonic oscillator.
Attributes:
Name Type Descriptionx0 tuple[float, float] the initial displacements
phi tuple[float, float] initial phases
Source code inhamilflow/models/harmonic_oscillator.py class ComplexSimpleHarmonicOscillatorIC(BaseModel):\n \"\"\"The initial condition for a complex harmonic oscillator.\n\n :cvar x0: the initial displacements\n :cvar phi: initial phases\n \"\"\"\n\n x0: tuple[float, float] = Field()\n phi: tuple[float, float] = Field(default=(0, 0))\nDampedHarmonicOscillator","text":" Bases: HarmonicOscillatorBase
Generate time series data for a damped harmonic oscillator.
The equation for a general un-driven harmonic oscillator is12
\\[ \\frac{\\mathrm d x^2}{\\mathrm d t^2} + 2\\zeta \\omega \\frac{\\mathrm d x}{\\mathrm dt} + \\omega^2 x = 0, \\]where \\(x\\) is the displacement, \\(\\omega\\) is the angular frequency of an undamped oscillator (\\(\\zeta=0\\)), and \\(\\zeta\\) is the damping ratio.
The solution to the above harmonic oscillator is
\\[ x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right) e^{-\\zeta \\omega t}, \\]where
\\[ \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}. \\]To use this generator,
params = {\"omega\": omega, \"zeta\"=0.2}\n\nho = DampedHarmonicOscillator(params=params)\n\ndf = ho(n_periods=1, n_samples_per_period=10)\ndf will be a pandas dataframe with two columns: t and x.
Contributors to Wikimedia projects. Harmonic oscillator. In: Wikipedia [Internet]. 18 Feb 2024 [cited 20 Feb 2024]. Available: https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator \u21a9
Libretexts. 5.3: General Solution for the Damped Harmonic Oscillator. Libretexts. 13 Apr 2021. Available: https://t.ly/cWTIo. Accessed 20 Feb 2024.\u00a0\u21a9
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the harmonic oscillator.
requiredinitial_condition Mapping[str, float] | None the initial condition of the harmonic oscillator.
None Source code in hamilflow/models/harmonic_oscillator.py class DampedHarmonicOscillator(HarmonicOscillatorBase):\n r\"\"\"Generate time series data for a [damped harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n The equation for a general un-driven harmonic oscillator is[^wiki_ho][^libretext_ho]\n\n $$\n \\frac{\\mathrm d x^2}{\\mathrm d t^2} + 2\\zeta \\omega \\frac{\\mathrm d x}{\\mathrm dt} + \\omega^2 x = 0,\n $$\n\n where $x$ is the displacement, $\\omega$ is the angular frequency of an undamped oscillator ($\\zeta=0$),\n and $\\zeta$ is the damping ratio.\n\n [^wiki_ho]: Contributors to Wikimedia projects. Harmonic oscillator.\n In: Wikipedia [Internet]. 18 Feb 2024 [cited 20 Feb 2024].\n Available: https://en.wikipedia.org/wiki/Harmonic_oscillator#Damped_harmonic_oscillator\n\n [^libretext_ho]: Libretexts. 5.3: General Solution for the Damped Harmonic Oscillator. Libretexts. 13 Apr 2021.\n Available: https://t.ly/cWTIo. Accessed 20 Feb 2024.\n\n\n The solution to the above harmonic oscillator is\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n\n To use this generator,\n\n ```python\n params = {\"omega\": omega, \"zeta\"=0.2}\n\n ho = DampedHarmonicOscillator(params=params)\n\n df = ho(n_periods=1, n_samples_per_period=10)\n ```\n\n `df` will be a pandas dataframe with two columns: `t` and `x`.\n\n :param system: all the params that defines the harmonic oscillator.\n :param initial_condition: the initial condition of the harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n super().__init__(system, initial_condition)\n if self.system.type == \"simple\":\n msg = (\n f\"System is not a Damped Harmonic Oscillator: {self.system}\\n\"\n f\"This is a simple harmonic oscillator, use `SimpleHarmonicOscillator`.\"\n )\n raise ValueError(\n msg,\n )\n\n def _x_under_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to under damped harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n \"\"\"\n t = np.asarray(t)\n omega_damp = self.system.omega * np.sqrt(1 - self.system.zeta)\n return (\n self.initial_condition.x0 * np.cos(omega_damp * t)\n + (\n self.system.zeta * self.system.omega * self.initial_condition.x0\n + self.initial_condition.v0\n )\n / omega_damp\n * np.sin(omega_damp * t)\n ) * np.exp(-self.system.zeta * self.system.omega * t)\n\n def _x_critical_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to critical damped harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cos(\\Omega t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Omega} \\sin(\\Omega t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Omega = \\omega\\sqrt{ 1 - \\zeta^2}.\n $$\n \"\"\"\n t = np.asarray(t)\n return self.initial_condition.x0 * np.exp(\n -self.system.zeta * self.system.omega * t,\n )\n\n def _x_over_damped(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to over harmonic oscillators.\n\n $$\n x(t) = \\left( x_0 \\cosh(\\Gamma t) + \\frac{\\zeta \\omega x_0 + v_0}{\\Gamma} \\sinh(\\Gamma t) \\right)\n e^{-\\zeta \\omega t},\n $$\n\n where\n\n $$\n \\Gamma = \\omega\\sqrt{ \\zeta^2 - 1 }.\n $$\n \"\"\"\n t = np.asarray(t)\n gamma_damp = self.system.omega * np.sqrt(self.system.zeta - 1)\n\n return (\n self.initial_condition.x0 * np.cosh(gamma_damp * t)\n + (\n self.system.zeta * self.system.omega * self.initial_condition.x0\n + self.initial_condition.v0\n )\n / gamma_damp\n * np.sinh(gamma_damp * t)\n ) * np.exp(-self.system.zeta * self.system.omega * t)\n\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to damped harmonic oscillators.\"\"\"\n t = np.asarray(t)\n if self.system.type == \"under_damped\":\n x = self._x_under_damped(t)\n elif self.system.type == \"over_damped\":\n x = self._x_over_damped(t)\n elif self.system.type == \"critical_damped\":\n x = self._x_critical_damped(t)\n else:\n msg = f\"System type is not damped harmonic oscillator: {self.system.type}\"\n raise ValueError(\n msg,\n )\n\n return x\nHarmonicOscillatorBase","text":" Bases: ABC
Base class to generate time series data for a harmonic oscillator.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] all the params that defines the harmonic oscillator.
requiredinitial_condition Mapping[str, float] | None the initial condition of the harmonic oscillator.
None Source code in hamilflow/models/harmonic_oscillator.py class HarmonicOscillatorBase(ABC):\n r\"\"\"Base class to generate time series data for a [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n :param system: all the params that defines the harmonic oscillator.\n :param initial_condition: the initial condition of the harmonic oscillator.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n initial_condition = initial_condition or {}\n self.system = HarmonicOscillatorSystem.model_validate(system)\n self.initial_condition = HarmonicOscillatorIC.model_validate(initial_condition)\n\n @cached_property\n def definition(self) -> dict[str, dict[str, float]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n @abstractmethod\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> npt.ArrayLike:\n r\"\"\"Solution to simple harmonic oscillators.\"\"\"\n ...\n\n def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time step.\n\n :param n_periods: Number of periods to generate.\n :param n_samples_per_period: Number of samples per period.\n \"\"\"\n time_delta = self.system.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n data = self._x(time_steps)\n\n return pd.DataFrame({\"t\": time_steps, \"x\": data})\ndefinition: dict[str, dict[str, float]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorBase.__call__","title":"__call__(n_periods, n_samples_per_period)","text":"Generate time series data for the harmonic oscillator.
Returns a list of floats representing the displacement at each time step.
Parameters:
Name Type Description Defaultn_periods int Number of periods to generate.
requiredn_samples_per_period int Number of samples per period.
required Source code inhamilflow/models/harmonic_oscillator.py def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator.\n\n Returns a list of floats representing the displacement at each time step.\n\n :param n_periods: Number of periods to generate.\n :param n_samples_per_period: Number of samples per period.\n \"\"\"\n time_delta = self.system.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n data = self._x(time_steps)\n\n return pd.DataFrame({\"t\": time_steps, \"x\": data})\nHarmonicOscillatorIC","text":" Bases: BaseModel
The initial condition for a harmonic oscillator.
Attributes:
Name Type Descriptionx0 float the initial displacement
v0 float the initial velocity
phi float initial phase
Source code inhamilflow/models/harmonic_oscillator.py class HarmonicOscillatorIC(BaseModel):\n \"\"\"The initial condition for a harmonic oscillator.\n\n :cvar x0: the initial displacement\n :cvar v0: the initial velocity\n :cvar phi: initial phase\n \"\"\"\n\n x0: float = Field(default=1.0)\n v0: float = Field(default=0.0)\n phi: float = Field(default=0.0)\nHarmonicOscillatorSystem","text":" Bases: BaseModel
The params for the harmonic oscillator.
Attributes:
Name Type Descriptionomega float angular frequency of the harmonic oscillator
zeta float damping ratio
Source code inhamilflow/models/harmonic_oscillator.py class HarmonicOscillatorSystem(BaseModel):\n \"\"\"The params for the harmonic oscillator.\n\n :cvar omega: angular frequency of the harmonic oscillator\n :cvar zeta: damping ratio\n \"\"\"\n\n omega: float = Field()\n zeta: float = Field(default=0.0)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def period(self) -> float:\n \"\"\"Period of the oscillator.\"\"\"\n return 2 * np.pi / self.omega\n\n @computed_field # type: ignore[misc]\n @cached_property\n def frequency(self) -> float:\n \"\"\"Frequency of the oscillator.\"\"\"\n return 1 / self.period\n\n @computed_field # type: ignore[misc]\n @cached_property\n def type(\n self,\n ) -> Literal[\"simple\", \"under_damped\", \"critical_damped\", \"over_damped\"]:\n \"\"\"Which type of harmonic oscillators.\"\"\"\n if self.zeta == 0:\n return \"simple\"\n elif self.zeta < 1:\n return \"under_damped\"\n elif self.zeta == 1:\n return \"critical_damped\"\n else:\n return \"over_damped\"\n\n @field_validator(\"zeta\")\n @classmethod\n def _check_zeta_non_negative(cls, v: float) -> float:\n if v < 0:\n msg = f\"Value of zeta should be positive: {v=}\"\n raise ValueError(msg)\n\n return v\nfrequency: float cached property","text":"Frequency of the oscillator.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorSystem.period","title":"period: float cached property","text":"Period of the oscillator.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.HarmonicOscillatorSystem.type","title":"type: Literal['simple', 'under_damped', 'critical_damped', 'over_damped'] cached property","text":"Which type of harmonic oscillators.
"},{"location":"references/models/harmonic_oscillator/#hamilflow.models.harmonic_oscillator.SimpleHarmonicOscillator","title":"SimpleHarmonicOscillator","text":" Bases: HarmonicOscillatorBase
Generate time series data for a simple harmonic oscillator.
In a one dimensional world, a mass \\(m\\), driven by a force \\(F=-kx\\), is described as
\\[ \\begin{align} F &= - k x \\\\ F &= m a \\end{align} \\]The mass behaves like a simple harmonic oscillator.
In general, the solution to a simple harmonic oscillator is
\\[ x(t) = A \\cos(\\omega t + \\phi), \\]where \\(\\omega\\) is the angular frequency, \\(\\phi\\) is the initial phase, and \\(A\\) is the amplitude.
To use this generator,
params = {\"omega\": omega}\n\nho = SimpleHarmonicOscillator(params=params)\n\ndf = ho(n_periods=1, n_samples_per_period=10)\ndf will be a pandas dataframe with two columns: t and x.
hamilflow/models/harmonic_oscillator.py class SimpleHarmonicOscillator(HarmonicOscillatorBase):\n r\"\"\"Generate time series data for a [simple harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).\n\n In a one dimensional world, a mass $m$, driven by a force $F=-kx$, is described as\n\n $$\n \\begin{align}\n F &= - k x \\\\\n F &= m a\n \\end{align}\n $$\n\n The mass behaves like a simple harmonic oscillator.\n\n In general, the solution to a simple harmonic oscillator is\n\n $$\n x(t) = A \\cos(\\omega t + \\phi),\n $$\n\n where $\\omega$ is the angular frequency, $\\phi$ is the initial phase, and $A$ is the amplitude.\n\n\n To use this generator,\n\n ```python\n params = {\"omega\": omega}\n\n ho = SimpleHarmonicOscillator(params=params)\n\n df = ho(n_periods=1, n_samples_per_period=10)\n ```\n\n `df` will be a pandas dataframe with two columns: `t` and `x`.\n \"\"\"\n\n def __init__(\n self,\n system: Mapping[str, float],\n initial_condition: Mapping[str, float] | None = None,\n ) -> None:\n super().__init__(system, initial_condition)\n if self.system.type != \"simple\":\n msg = f\"System is not a Simple Harmonic Oscillator: {self.system}\"\n raise ValueError(\n msg,\n )\n\n def _x(self, t: \"Sequence[float] | npt.ArrayLike\") -> np.ndarray:\n r\"\"\"Solution to simple harmonic oscillators.\n\n $$\n x(t) = x_0 \\cos(\\omega t + \\phi).\n $$\n \"\"\"\n return self.initial_condition.x0 * np.cos(\n self.system.omega * np.asarray(t) + self.initial_condition.phi,\n )\nMain module for a harmonic oscillator chain.
"},{"location":"references/models/harmonic_oscillator_chain/#hamilflow.models.harmonic_oscillator_chain.HarmonicOscillatorsChain","title":"HarmonicOscillatorsChain","text":"Generate time series data for a coupled harmonic oscillator chain with periodic boundary condition.
A one-dimensional circle of \\(N\\) interacting harmonic oscillators can be described by the Lagrangian action \\(\\(S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\}\\,,\\)\\) where \\(x_N \\coloneqq x_0\\).
This system can be solved in terms of travelling waves, obtained by discrete Fourier transform.
Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy \\(Y_k = Y^*_{-k \\mod N}\\), see Wikipedia.
Parameters:
Name Type Description Defaultomega float frequence parameter
requiredinitial_conditions Sequence[Mapping[str, float | tuple[float, float]]] a sequence of initial conditions on the Fourier modes. The first element in the sequence is that of the zero mode, taking a position and a velocity. Rest of the elements are that of the independent travelling waves, taking two amplitudes and two initial phases.
requiredodd_dof bool The system will have 2 * len(initial_conditions) + int(odd_dof) - 2 degrees of freedom.
hamilflow/models/harmonic_oscillator_chain.py class HarmonicOscillatorsChain:\n r\"\"\"Generate time series data for a coupled harmonic oscillator chain with periodic boundary condition.\n\n A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action\n $$S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\}\\,,$$\n where $x_N \\coloneqq x_0$.\n\n This system can be solved in terms of _travelling waves_, obtained by discrete Fourier transform.\n\n Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy\n $Y_k = Y^*_{-k \\mod N}$, see [Wikipedia](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#DFT_of_real_and_purely_imaginary_signals).\n\n :param omega: frequence parameter\n :param initial_conditions: a sequence of initial conditions on the Fourier modes.\n The first element in the sequence is that of the zero mode, taking a position and a velocity.\n Rest of the elements are that of the independent travelling waves, taking two amplitudes and two initial phases.\n :param odd_dof: The system will have `2 * len(initial_conditions) + int(odd_dof) - 2` degrees of freedom.\n \"\"\"\n\n def __init__(\n self,\n omega: float,\n initial_conditions: Sequence[Mapping[str, float | tuple[float, float]]],\n odd_dof: bool,\n ) -> None:\n self.n_dof = 2 * len(initial_conditions) + odd_dof - 2\n if not odd_dof:\n prefix = \"For even degrees of freedom, \"\n if self.n_dof == 0:\n raise ValueError(prefix + \"at least 1 travelling wave is needed\")\n amp = cast(tuple[float, float], initial_conditions[-1][\"amp\"])\n if amp[0] != amp[1]:\n msg = \"k == N // 2 must have equal positive and negative amplitudes.\"\n raise ValueError(prefix + msg)\n self.omega = omega\n self.odd_dof = odd_dof\n\n self.free_mode = FreeParticle(cast(Mapping[str, float], initial_conditions[0]))\n\n self.independent_csho_modes = [\n self._sho_factory(\n k,\n cast(tuple[float, float], ic[\"amp\"]),\n cast(tuple[float, float] | None, ic.get(\"phi\")),\n )\n for k, ic in enumerate(initial_conditions[1:], 1)\n ]\n\n def _sho_factory(\n self,\n k: int,\n amp: tuple[float, float],\n phi: tuple[float, float] | None = None,\n ) -> ComplexSimpleHarmonicOscillator:\n return ComplexSimpleHarmonicOscillator(\n {\n \"omega\": 2 * self.omega * np.sin(np.pi * k / self.n_dof),\n },\n {\"x0\": amp} | ({\"phi\": phi} if phi else {}),\n )\n\n @cached_property\n def definition(\n self,\n ) -> dict[\n str,\n float\n | dict[str, dict[str, float | list[float]]]\n | list[dict[str, dict[str, float | tuple[float, float]]]],\n ]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"omega\": self.omega,\n \"n_dof\": self.n_dof,\n \"free_mode\": self.free_mode.definition,\n \"independent_csho_modes\": [\n rwm.definition for rwm in self.independent_csho_modes\n ],\n }\n\n def _z(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"tuple[npt.NDArray[np.complex64], npt.NDArray[np.complex64]]\":\n t = np.asarray(t).reshape(-1)\n all_travelling_waves = [self.free_mode._x(t).reshape(1, -1)] # noqa: SLF001\n\n if self.independent_csho_modes:\n independent_cshos = np.asarray(\n [o._z(t) for o in self.independent_csho_modes], # noqa: SLF001\n )\n all_travelling_waves.extend(\n (\n (independent_cshos, independent_cshos[::-1].conj())\n if self.odd_dof\n else (\n independent_cshos[:-1],\n independent_cshos[[-1]],\n independent_cshos[-2::-1].conj(),\n )\n ),\n )\n\n travelling_waves = np.concatenate(all_travelling_waves)\n original_zs = ifft(travelling_waves, axis=0, norm=\"ortho\")\n return original_zs, travelling_waves\n\n def _x(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"tuple[npt.NDArray[np.float64], npt.NDArray[np.complex64]]\":\n original_xs, travelling_waves = self._z(t)\n\n return np.real(original_xs), travelling_waves\n\n def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator chain.\n\n Returns float(s) representing the displacement at the given time(s).\n\n :param t: time.\n \"\"\"\n t = np.asarray(t)\n original_xs, travelling_waves = self._x(t)\n data = { # type: ignore [var-annotated]\n f\"{name}{i}\": cast(\n \"npt.NDArray[np.float64] | npt.NDArray[np.complex64]\",\n values,\n )\n for name, xs in zip(\n (\"x\", \"y\"),\n (original_xs, travelling_waves),\n strict=False,\n )\n for i, values in enumerate(xs) # type: ignore [arg-type]\n }\n\n return pd.DataFrame(data, index=t)\ndefinition: dict[str, float | dict[str, dict[str, float | list[float]]] | list[dict[str, dict[str, float | tuple[float, float]]]]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/harmonic_oscillator_chain/#hamilflow.models.harmonic_oscillator_chain.HarmonicOscillatorsChain.__call__","title":"__call__(t)","text":"Generate time series data for the harmonic oscillator chain.
Returns float(s) representing the displacement at the given time(s).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time.
required Source code inhamilflow/models/harmonic_oscillator_chain.py def __call__(self, t: \"Sequence[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Generate time series data for the harmonic oscillator chain.\n\n Returns float(s) representing the displacement at the given time(s).\n\n :param t: time.\n \"\"\"\n t = np.asarray(t)\n original_xs, travelling_waves = self._x(t)\n data = { # type: ignore [var-annotated]\n f\"{name}{i}\": cast(\n \"npt.NDArray[np.float64] | npt.NDArray[np.complex64]\",\n values,\n )\n for name, xs in zip(\n (\"x\", \"y\"),\n (original_xs, travelling_waves),\n strict=False,\n )\n for i, values in enumerate(xs) # type: ignore [arg-type]\n }\n\n return pd.DataFrame(data, index=t)\nMain module for a pendulum.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum","title":"Pendulum","text":"Generate time series data for a pendulum.
We describe a generic pendulum system by the Lagrangian action $$ S_L[\\theta] = I \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\}\\,, $$ where \\(\\theta\\) is the angle from the vertical to the pendulum; \\(I\\) is the inertia parameter introduced for dimensional reasons, and \\(\\omega_0\\) the frequency parameter.
Details are collected in the tutorial.
Source code inhamilflow/models/pendulum.py class Pendulum:\n r\"\"\"Generate time series data for a pendulum.\n\n We describe a generic pendulum system by the Lagrangian action\n $$\n S_L\\[\\theta\\] = I \\int_{t_0}^{t_1} \\mathbb{d}t\n \\left\\\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\\\}\\,,\n $$\n where $\\theta$ is the _angle_ from the vertical to the pendulum;\n $I$ is the _inertia parameter_ introduced for dimensional reasons,\n and $\\omega_0$ the _frequency parameter_.\n\n Details are collected in the tutorial.\n \"\"\"\n\n def __init__(\n self,\n system: float | Mapping[str, float],\n initial_condition: float | Mapping[str, float],\n ) -> None:\n if isinstance(system, float | int):\n system = {\"omega0\": system}\n if isinstance(initial_condition, float | int):\n initial_condition = {\"theta0\": initial_condition}\n self.system = PendulumSystem.model_validate(system)\n self.initial_condition = PendulumIC.model_validate(initial_condition)\n\n @cached_property\n def definition(self) -> dict[str, dict[str, Any]]:\n \"\"\"Model params and initial conditions defined as a dictionary.\"\"\"\n return {\n \"system\": self.system.model_dump(),\n \"initial_condition\": self.initial_condition.model_dump(),\n }\n\n @property\n def omega0(self) -> float:\n \"\"\"Original angular frequency of the system.\"\"\"\n return self.system.omega0\n\n @property\n def _k(self) -> float:\n return self.initial_condition.k\n\n @property\n def _math_m(self) -> float:\n return self._k**2\n\n @cached_property\n def freq(self) -> float:\n r\"\"\"Frequency.\n\n :return: $\\frac{\\pi}{2K(k^2)}\\omega_0$, where\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n return math.pi * self.omega0 / (2 * ellipk(self._math_m))\n\n @cached_property\n def period(self) -> float:\n r\"\"\"Period.\n\n :return: $\\frac{4K(k^2)}{\\omega_0}$, where\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n return 4 * ellipk(self._math_m) / self.omega0\n\n def _math_u(\n self,\n t: \"Sequence[float] | npt.ArrayLike\",\n ) -> \"npt.NDArray[np.float64]\":\n return self.omega0 * np.asarray(t)\n\n def u(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Give the convenient generalised coordinate $u$, $\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}$.\n\n :param t: time\n :return: $u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)$, where\n $\\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)\n\n return ph\n\n def theta(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Angle $\\theta$.\n\n :param t: time\n :return: $\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)$, where\n $\\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)\n \"\"\"\n _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)\n\n return 2 * np.arcsin(cn / dn * self._k)\n\n def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate the time sequence from more interpretable params.\n\n :param n_periods: number of periods to include\n :param n_samples_per_period: number of samples in each period\n :return: an array that contains all the timesteps\n \"\"\"\n time_delta = self.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n return self(time_steps)\n\n def __call__(self, t: TypeTime) -> pd.DataFrame:\n \"\"\"Generate the variables of the pendulum in time together with the time steps.\n\n :param t: time steps\n :return: values of the variables\n angle `x`, generalized coordinates `u`, and time `t`.\n \"\"\"\n thetas = self.theta(t)\n us = self.u(t)\n\n return pd.DataFrame({\"t\": t, \"x\": thetas, \"u\": us})\ndefinition: dict[str, dict[str, Any]] cached property","text":"Model params and initial conditions defined as a dictionary.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.freq","title":"freq: float cached property","text":"Frequency.
Returns:
Type Descriptionfloat \\(\\frac{\\pi}{2K(k^2)}\\omega_0\\), where \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.omega0","title":"omega0: float property","text":"Original angular frequency of the system.
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.period","title":"period: float cached property","text":"Period.
Returns:
Type Descriptionfloat \\(\\frac{4K(k^2)}{\\omega_0}\\), where \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.Pendulum.__call__","title":"__call__(t)","text":"Generate the variables of the pendulum in time together with the time steps.
Parameters:
Name Type Description Defaultt TypeTime time steps
requiredReturns:
Type DescriptionDataFrame values of the variables angle x, generalized coordinates u, and time t.
hamilflow/models/pendulum.py def __call__(self, t: TypeTime) -> pd.DataFrame:\n \"\"\"Generate the variables of the pendulum in time together with the time steps.\n\n :param t: time steps\n :return: values of the variables\n angle `x`, generalized coordinates `u`, and time `t`.\n \"\"\"\n thetas = self.theta(t)\n us = self.u(t)\n\n return pd.DataFrame({\"t\": t, \"x\": thetas, \"u\": us})\ngenerate_from(n_periods, n_samples_per_period)","text":"Generate the time sequence from more interpretable params.
Parameters:
Name Type Description Defaultn_periods int number of periods to include
requiredn_samples_per_period int number of samples in each period
requiredReturns:
Type DescriptionDataFrame an array that contains all the timesteps
Source code inhamilflow/models/pendulum.py def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:\n \"\"\"Generate the time sequence from more interpretable params.\n\n :param n_periods: number of periods to include\n :param n_samples_per_period: number of samples in each period\n :return: an array that contains all the timesteps\n \"\"\"\n time_delta = self.period / n_samples_per_period\n time_steps = np.arange(0, n_periods * n_samples_per_period) * time_delta\n\n return self(time_steps)\ntheta(t)","text":"Angle \\(\\theta\\).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time
requiredReturns:
Type DescriptionNDArray[float64] \\(\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)\\), where \\(\\mathrm{cd}(z, k)\\) is a Jacobian elliptic function
Source code inhamilflow/models/pendulum.py def theta(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Angle $\\theta$.\n\n :param t: time\n :return: $\\theta(t) = 2\\arcsin\\!\\big(k\\cdot\\mathrm{cd}(\\omega_0 t, k^2)\\big)$, where\n $\\mathrm{cd}(z, k)$ is a [Jacobian elliptic function](https://dlmf.nist.gov/22.2#E8)\n \"\"\"\n _, cn, dn, _ = ellipj(self._math_u(t), self._math_m)\n\n return 2 * np.arcsin(cn / dn * self._k)\nu(t)","text":"Give the convenient generalised coordinate \\(u\\), \\(\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}\\).
Parameters:
Name Type Description Defaultt Sequence[float] | ArrayLike time
requiredReturns:
Type DescriptionNDArray[float64] \\(u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)\\), where \\(\\mathrm{am}(x, k)\\) is Jacobi's amplitude function, \\(K(m)\\) is Legendre's complete elliptic integral of the first kind
Source code inhamilflow/models/pendulum.py def u(self, t: \"Sequence[float] | npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Give the convenient generalised coordinate $u$, $\\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{\\sin\\frac{\\theta_0}{2}}$.\n\n :param t: time\n :return: $u(t) = \\mathrm{am}\\!\\big(\\omega_0 t + K(k^2), k^2\\big)$, where\n $\\mathrm{am}(x, k)$ is [Jacobi's amplitude function](https://dlmf.nist.gov/22.16#E1),\n $K(m)$ is [Legendre's complete elliptic integral of the first kind](https://dlmf.nist.gov/19.2#E8)\n \"\"\"\n _, _, _, ph = ellipj(self._math_u(t) + ellipk(self._math_m), self._math_m)\n\n return ph\nPendulumIC","text":" Bases: BaseModel
The initial condition for a pendulum.
Parameters:
Name Type Description Defaulttheta0 \\(-\\frac{\\pi}{2} \\le \\theta_0 \\le \\frac{\\pi}{2}\\), the initial angle
required Source code inhamilflow/models/pendulum.py class PendulumIC(BaseModel):\n r\"\"\"The initial condition for a pendulum.\n\n :param theta0: $-\\frac{\\pi}{2} \\le \\theta_0 \\le \\frac{\\pi}{2}$, the\n initial angle\n \"\"\"\n\n theta0: float = Field(ge=-math.pi / 2, le=math.pi / 2, frozen=True)\n\n @computed_field # type: ignore[misc]\n @cached_property\n def k(self) -> float:\n r\"\"\"A convenient number for elliptic functions.\n\n :return: $\\sin\\frac{\\theta_0}{2}$\n \"\"\"\n return math.sin(self.theta0 / 2)\nk: float cached property","text":"A convenient number for elliptic functions.
Returns:
Type Descriptionfloat \\(\\sin\\frac{\\theta_0}{2}\\)
"},{"location":"references/models/pendulum/#hamilflow.models.pendulum.PendulumSystem","title":"PendulumSystem","text":" Bases: BaseModel
The params for the pendulum.
Parameters:
Name Type Description Defaultomega0 \\(\\omega_0 \\coloneqq \\sqrt{\\frac{U}{I}} > 0\\), frequency parameter
required Source code inhamilflow/models/pendulum.py class PendulumSystem(BaseModel):\n r\"\"\"The params for the pendulum.\n\n :param omega0: $\\omega_0 \\coloneqq \\sqrt{\\frac{U}{I}} > 0$, frequency\n parameter\n \"\"\"\n\n omega0: float = Field(gt=0.0, frozen=True)\nExact solution of Kepler dynamics.
"},{"location":"references/models/kepler_problem/dynamics/#hamilflow.models.kepler_problem.dynamics.esolve_u_from_tau_parabolic","title":"esolve_u_from_tau_parabolic(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.
Let \\(\\kappa = 1+9\\tau^2\\), $$ u = -1 + \\left(\\frac{3\\tau}{\\kappa^\\frac{3}{2}} + \\frac{1}{\\kappa}\\right)^\\frac{1}{3} + \\left(\\frac{1}{3\\tau \\kappa^\\frac{3}{2}+\\kappa^2}\\right)^\\frac{1}{3}\\,. $$
Parameters:
Name Type Description Defaultecc float eccentricity, ecc = 0 (unchecked, unused)
requiredtau ArrayLike scaled time
requiredReturns:
Type DescriptionNDArray[float64] convenient radial inverse u
Source code inhamilflow/models/kepler_problem/dynamics.py def esolve_u_from_tau_parabolic(\n ecc: float, # noqa: ARG001\n tau: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.\n\n Let $\\kappa = 1+9\\tau^2$,\n $$ u = -1\n + \\left(\\frac{3\\tau}{\\kappa^\\frac{3}{2}} + \\frac{1}{\\kappa}\\right)^\\frac{1}{3}\n + \\left(\\frac{1}{3\\tau \\kappa^\\frac{3}{2}+\\kappa^2}\\right)^\\frac{1}{3}\\,. $$\n\n :param ecc: eccentricity, ecc = 0 (unchecked, unused)\n :param tau: scaled time\n :return: convenient radial inverse u\n \"\"\"\n tau = np.asarray(tau)\n tau_3 = 3 * tau\n term = 1 + tau_3**2 # 1 + 9 * tau**2\n term1_5 = term**1.5 # (1 + 9 * tau**2)**1.5\n second_term = (tau_3 / term1_5 + 1 / term) ** _1_3\n third_term = 1 / (tau_3 * term1_5 + term**2) ** _1_3\n return -1 + second_term + third_term\ntau_of_1_plus_u_hyperbolic(ecc, u)","text":"Expansion for tau of u in the hyperbolic case at \\(u = -1+0\\).
The exact solution has a logarithmic singularity at \\(u = -1\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\frac{e(1+u)}{2(e^2-1)}\\), $$ \\tau = \\ln \\epsilon + \\frac{e}{2\\epsilon} + 1 - \\frac{e^2+2}{e}\\epsilon + \\left(3+\\frac{2}{e^2}\\right)\\epsilon^2 + O(\\epsilon^3)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_1_plus_u_hyperbolic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the hyperbolic case at $u = -1+0$.\n\n The exact solution has a logarithmic singularity at $u = -1$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\frac{e(1+u)}{2(e^2-1)}$,\n $$ \\tau = \\ln \\epsilon + \\frac{e}{2\\epsilon}\n + 1 - \\frac{e^2+2}{e}\\epsilon + \\left(3+\\frac{2}{e^2}\\right)\\epsilon^2\n + O(\\epsilon^3)\\,. $$\n \"\"\"\n cosqr = ecc**2 - 1\n up1 = ecc * (1 + u) / 2 / cosqr\n diverge = np.log(up1) + ecc / 2 / up1\n regular = 1 - (2 + ecc**2) / ecc * up1 + (3 + 2 / ecc**2) * up1**2\n return (diverge + regular) / cosqr**1.5\ntau_of_e_minus_u_elliptic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = +e-0\\).
The exact solution has a removable singularity at \\(u = +e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}\\), $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon - \\frac{1+9e}{24(1+e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_minus_u_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = +e-0$.\n\n The exact solution has a removable singularity at $u = +e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}$,\n $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon\n - \\frac{1+9e}{24(1+e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n emu = np.sqrt(2 * (ecc - u) / ecc)\n return emu / (1 + ecc) ** 2 - emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3\ntau_of_e_minus_u_hyperbolic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = +e-0\\).
The exact solution has a removable singularity at \\(u = +e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}\\), $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon + \\frac{1+9e}{24(1+e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_minus_u_hyperbolic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = +e-0$.\n\n The exact solution has a removable singularity at $u = +e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e-u)}{e}}$,\n $$ \\tau = \\frac{1}{(1+e)^2}\\epsilon\n + \\frac{1+9e}{24(1+e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n emu = np.sqrt(2 * (ecc - u) / ecc)\n return emu / (1 + ecc) ** 2 + emu**3 * (1 + 9 * ecc) / 24 / (1 + ecc) ** 3\ntau_of_e_plus_u_elliptic(ecc, u)","text":"Expansion for tau of u in the elliptic case at \\(u = -e+0\\).
The exact solution has a removable singularity at \\(u = -e\\), hence this expansion helps with numerics.
Let \\(\\epsilon = \\sqrt{\\frac{2(e+u)}{e}}\\), $$ \\tau = \\frac{\\pi}{(1-e^2)^\\frac{3}{2}} - \\frac{1}{(1-e)^2}\\epsilon - \\frac{1-9e}{24(1-e)^3}\\epsilon^3 + O\\left(\\epsilon^5\\right)\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_e_plus_u_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Expansion for tau of u in the elliptic case at $u = -e+0$.\n\n The exact solution has a removable singularity at $u = -e$, hence this\n expansion helps with numerics.\n\n Let $\\epsilon = \\sqrt{\\frac{2(e+u)}{e}}$,\n $$ \\tau = \\frac{\\pi}{(1-e^2)^\\frac{3}{2}}\n - \\frac{1}{(1-e)^2}\\epsilon\n - \\frac{1-9e}{24(1-e)^3}\\epsilon^3\n + O\\left(\\epsilon^5\\right)\\,. $$\n \"\"\"\n epu = np.sqrt(2 * (ecc + u) / ecc)\n const = np.pi / (1 - ecc**2) ** 1.5\n return const - epu / (ecc - 1) ** 2 - epu**3 * (1 - 9 * ecc) / 24 / (1 - ecc) ** 3\ntau_of_u_elliptic(ecc, u)","text":"Calculate the scaled time tau from u in the elliptic case.
Parameters:
Name Type Description Defaultecc float eccentricity, 0 < ecc < 1 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_elliptic(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the scaled time tau from u in the elliptic case.\n\n :param ecc: eccentricity, 0 < ecc < 1 (unchecked)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n return _approximate_at_termina(\n ecc,\n u,\n tau_of_u_exact_elliptic,\n tau_of_e_plus_u_elliptic,\n tau_of_e_minus_u_elliptic,\n )\ntau_of_u_exact_elliptic(ecc, u)","text":"Exact solution for tau of u in the elliptic case.
For \\(-e \\le u \\le e\\), $$ \\tau = -\\frac{\\sqrt{e^2-u^2}}{(1-e^2)(1+u)} + \\frac{\\frac{\\pi}{2} - \\arctan\\frac{e^2+u}{\\sqrt{(1-e^2)(e^2-u^2)}}}{(1-e^2)^{\\frac{3}{2}}}\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_exact_elliptic(\n ecc: float,\n u: \"npt.NDArray[np.float64]\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Exact solution for tau of u in the elliptic case.\n\n For $-e \\le u \\le e$,\n $$ \\tau = -\\frac{\\sqrt{e^2-u^2}}{(1-e^2)(1+u)}\n + \\frac{\\frac{\\pi}{2} - \\arctan\\frac{e^2+u}{\\sqrt{(1-e^2)(e^2-u^2)}}}{(1-e^2)^{\\frac{3}{2}}}\\,. $$\n \"\"\"\n cosqr, eusqrt = 1 - ecc**2, np.sqrt(ecc**2 - u**2)\n trig_numer = np.pi / 2 - np.arctan((ecc**2 + u) / np.sqrt(cosqr) / eusqrt)\n return -eusqrt / cosqr / (1 + u) + trig_numer / cosqr**1.5\ntau_of_u_exact_hyperbolic(ecc, u)","text":"Exact solution for tau of u in the hyperbolic case.
For \\(-1 < u \\le e\\), $$ \\tau = \\frac{\\sqrt{e^2-u^2}}{(e^2-1)(1+u)} - \\frac{\\mathrm{artanh}\\frac{e^2+u}{\\sqrt{(e^2-1)(e^2-u^2)}}}{(e^2-1)^{\\frac{3}{2}}}\\,. $$
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_exact_hyperbolic(\n ecc: float,\n u: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Exact solution for tau of u in the hyperbolic case.\n\n For $-1 < u \\le e$,\n $$ \\tau = \\frac{\\sqrt{e^2-u^2}}{(e^2-1)(1+u)}\n - \\frac{\\mathrm{artanh}\\frac{e^2+u}{\\sqrt{(e^2-1)(e^2-u^2)}}}{(e^2-1)^{\\frac{3}{2}}}\\,. $$\n \"\"\"\n u = np.asarray(u)\n cosqr, eusqrt = ecc**2 - 1, np.sqrt(ecc**2 - u**2)\n trig_numer = np.arctanh(np.sqrt(cosqr) * eusqrt / (ecc**2 + u))\n return eusqrt / cosqr / (1 + u) - trig_numer / cosqr**1.5\ntau_of_u_hyperbolic(ecc, u)","text":"Calculate the scaled time tau from u in the hyperbolic case.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 1 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_hyperbolic(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the scaled time tau from u in the hyperbolic case.\n\n :param ecc: eccentricity, ecc > 1 (unchecked)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n return _approximate_at_termina(\n ecc,\n u,\n tau_of_u_exact_hyperbolic,\n tau_of_1_plus_u_hyperbolic,\n tau_of_e_minus_u_hyperbolic,\n )\ntau_of_u_parabolic(ecc, u)","text":"Calculate the scaled time tau from u in the parabolic case.
For \\(-1 < u \\le 1\\), $$ \\tau = \\frac{\\sqrt{1-u}(2+u)}{3(1+u)^\\frac{3}{2}}\\,. $$
Parameters:
Name Type Description Defaultecc float eccentricity, ecc == 1 (unchecked, unused)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] scaled time tau
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_parabolic(\n ecc: float, # noqa: ARG001\n u: \"npt.ArrayLike\",\n) -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the scaled time tau from u in the parabolic case.\n\n For $-1 < u \\le 1$,\n $$ \\tau = \\frac{\\sqrt{1-u}(2+u)}{3(1+u)^\\frac{3}{2}}\\,. $$\n\n :param ecc: eccentricity, ecc == 1 (unchecked, unused)\n :param u: convenient radial inverse\n :return: scaled time tau\n \"\"\"\n u = np.asarray(u)\n return np.sqrt(1 - u) * (2 + u) / 3 / (1 + u) ** 1.5\ntau_of_u_prime(ecc, u)","text":"Calculate the first derivative of scaled time tau with respect to u.
\\[ \\tau'(u) = -\\frac{1}{(1+u)^2\\sqrt{e^2-u^2}}\\,. \\]Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] the first derivative scaled time tau with respect to u
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_prime(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the first derivative of scaled time tau with respect to u.\n\n $$ \\tau'(u) = -\\frac{1}{(1+u)^2\\sqrt{e^2-u^2}}\\,. $$\n\n :param ecc: eccentricity, ecc >= 0 (unchecked)\n :param u: convenient radial inverse\n :return: the first derivative scaled time tau with respect to u\n \"\"\"\n u = np.asarray(u)\n return -1 / (1 + u) ** 2 / np.sqrt(ecc**2 - u**2)\ntau_of_u_prime2(ecc, u)","text":"Calculate the second derivative of scaled time tau with respect to u.
\\[ \\tau''(u) = \\frac{2e^2 - u - 3u^2}{(1+u)^3 (e^2-u^2)^\\frac{3}{2}} \\]Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0 (unchecked)
requiredu ArrayLike convenient radial inverse
requiredReturns:
Type DescriptionNDArray[float64] the second derivative scaled time tau with respect to u
Source code inhamilflow/models/kepler_problem/dynamics.py def tau_of_u_prime2(ecc: float, u: \"npt.ArrayLike\") -> \"npt.NDArray[np.float64]\":\n r\"\"\"Calculate the second derivative of scaled time tau with respect to u.\n\n $$ \\tau''(u) = \\frac{2e^2 - u - 3u^2}{(1+u)^3 (e^2-u^2)^\\frac{3}{2}} $$\n\n :param ecc: eccentricity, ecc >= 0 (unchecked)\n :param u: convenient radial inverse\n :return: the second derivative scaled time tau with respect to u\n \"\"\"\n u = np.asarray(u)\n u2 = u**2\n return (2 * ecc**2 - u - 3 * u2) / (1 + u) ** 3 / (ecc**2 - u2) ** 1.5\nMain module for Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D","title":"Kepler2D","text":"Kepler problem in two dimensional space.
Parameters:
Name Type Description Defaultsystem Mapping[str, float] the Kepler problem system specification
requiredfirst_integrals Mapping[str, float] the first integrals for the system.
required Source code inhamilflow/models/kepler_problem/model.py class Kepler2D:\n \"\"\"Kepler problem in two dimensional space.\n\n :param system: the Kepler problem system specification\n :param first_integrals: the first integrals for the system.\n \"\"\"\n\n def __init__(\n self,\n system: \"Mapping[str, float]\",\n first_integrals: \"Mapping[str, float]\",\n ) -> None:\n self.system = Kepler2DSystem.model_validate(system)\n\n first_integrals = dict(first_integrals)\n ene = first_integrals[\"ene\"]\n minimal_ene = Kepler2D.minimal_ene(first_integrals[\"angular_mom\"], system)\n if ene < minimal_ene:\n msg = f\"Energy {ene} less than minimally allowed {minimal_ene}\"\n raise ValueError(msg)\n\n self.first_integrals = Kepler2DFI.model_validate(first_integrals)\n\n if 0 <= self.ecc < 1:\n self.tau_of_u = partial(tau_of_u_elliptic, self.ecc)\n elif self.ecc == 1:\n self.tau_of_u = partial(tau_of_u_parabolic, self.ecc)\n elif self.ecc > 1:\n self.tau_of_u = partial(tau_of_u_hyperbolic, self.ecc)\n else:\n raise RuntimeError\n\n @classmethod\n def from_geometry(\n cls,\n system: \"Mapping[str, float]\",\n geometries: \"Mapping[str, bool | float]\",\n ) -> \"Self\":\n r\"\"\"Alternative initialiser from system and geometry specifications.\n\n Given the eccentricity $e$ and the conic section semi-latus rectum $p$,\n $$l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.$$\n\n :param system: the Kepler problem system specification\n :param geometries: geometric specifications\n `positive_angular_mom`: whether the angular momentum is positive\n `ecc`: eccentricity of the conic section\n `parameter`: semi-latus rectum of the conic section\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n positive_angular_mom = bool(geometries[\"positive_angular_mom\"])\n ecc, parameter = (float(geometries[k]) for k in [\"ecc\", \"parameter\"])\n abs_angular_mom = math.sqrt(mass * parameter * alpha)\n # abs_minimal_ene = alpha / 2 / parameter: numerically unstable\n abs_minimal_ene = abs(cls.minimal_ene(abs_angular_mom, system))\n ene = (ecc**2 - 1) * abs_minimal_ene\n fi = {\n \"ene\": ene,\n \"angular_mom\": (\n abs_angular_mom if positive_angular_mom else -abs_angular_mom\n ),\n }\n return cls(system, fi)\n\n @staticmethod\n def minimal_ene(\n angular_mom: float,\n system: \"Mapping[str, float]\",\n ) -> float:\n r\"\"\"Minimal possible energy from the system specification and an angular momentum.\n\n $$ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. $$\n\n :param angular_mom: angular momentum\n :param system: system specification\n :return: minimal possible energy\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n return -mass * alpha**2 / (2 * angular_mom**2)\n\n @property\n def mass(self) -> float:\n \"\"\"Mass $m$ from the system specification.\"\"\"\n return self.system.mass\n\n @property\n def alpha(self) -> float:\n r\"\"\"Alpha $\\alpha$ from the system specification.\"\"\"\n return self.system.alpha\n\n @property\n def ene(self) -> float:\n \"\"\"Energy $E$ of the Kepler problem.\"\"\"\n return self.first_integrals.ene\n\n @property\n def angular_mom(self) -> float:\n \"\"\"Angular momentum $l$ of the Kepler problem.\"\"\"\n return self.first_integrals.angular_mom\n\n @property\n def t0(self) -> float:\n r\"\"\"t0 $t_0$ of the Kepler problem.\"\"\"\n return self.first_integrals.t0\n\n @property\n def phi0(self) -> float:\n r\"\"\"phi0 $\\phi_0$ of the Kepler problem.\"\"\"\n return self.first_integrals.phi0\n\n @cached_property\n def period(self) -> float:\n r\"\"\"Period $T$ of the Kepler problem.\n\n For $E < 0$,\n $$ T = \\pi \\alpha \\sqrt{-\\frac{m}{2E^3}}\\,. $$\n \"\"\"\n if self.ene >= 0:\n msg = f\"Only systems with energy < 0 have a period, got {self.ene}\"\n raise TypeError(msg)\n return math.pi * self.alpha * math.sqrt(-self.mass / 2 / self.ene**3)\n\n # FIXME is it called parameter in English?\n @cached_property\n def parameter(self) -> float:\n r\"\"\"Conic section semi-latus rectum of the Kepler problem.\n\n $$ p = \\frac{l^2}{\\alpha m}\\,. $$\n \"\"\"\n return self.angular_mom**2 / self.mass / self.alpha\n\n @cached_property\n def ecc(self) -> float:\n r\"\"\"Conic section eccentricity of the Kepler problem.\n\n $$ e = \\sqrt{1 + \\frac{2El}{\\alpha^2 m}}\\,. $$\n \"\"\"\n return math.sqrt(\n 1 + 2 * self.ene * self.angular_mom**2 / self.mass / self.alpha**2,\n )\n\n @cached_property\n def period_in_tau(self) -> float:\n r\"\"\"Period in the scaled time tau.\n\n $$ T_\\tau = \\frac{2\\pi}{(1-e^2)^\\frac{3}{2}}\\,. $$\n \"\"\"\n if self.ecc >= 1:\n msg = (\n f\"Only systems with 0 <= eccentricity < 1 have a period, got {self.ecc}\"\n )\n raise TypeError(\n msg,\n )\n return 2 * math.pi / (1 - self.ecc**2) ** 1.5\n\n @property\n def t_to_tau_factor(self) -> float:\n r\"\"\"Scale factor from t to tau.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return abs(self.mass * self.alpha**2 / self.angular_mom**3)\n\n def tau(self, t: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the scaled time tau from t.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return (np.asarray(t) - self.t0) * self.t_to_tau_factor\n\n def u_of_tau(self, tau: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n \"\"\"Give the convenient radial inverse u from tau.\"\"\"\n tau = np.asarray(tau)\n if self.ecc == 0:\n return np.zeros(tau.shape)\n else:\n if self.ecc < 1:\n p = self.period_in_tau\n r = tau % p\n tau = np.where(r <= p / 2, r, p - r)\n else:\n tau = np.abs(tau)\n return u_of_tau(self.ecc, tau) # type: ignore [arg-type]\n\n def r_of_u(self, u: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the radial r from u.\n\n $$ r = \\frac{p}{u+1}\\,. $$\n \"\"\"\n return self.parameter / (np.asarray(u) + 1)\n\n def phi_of_u_tau(\n self,\n u: \"Collection[float] | npt.ArrayLike\",\n tau: \"Collection[float] | npt.ArrayLike\",\n ) -> \"npt.ArrayLike\":\n r\"\"\"Give the angular phi from u and tau.\n\n For $e = 0$,\n $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$\n For $e > 0$,\n $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$\n \"\"\"\n u, tau = np.asarray(u), np.asarray(tau)\n if self.ecc == 0:\n phi = 2 * math.pi * tau / self.period_in_tau\n else:\n if self.ecc < 1:\n shift = tau / self.period_in_tau\n else:\n shift = np.where(tau >= 0, 0, -np.pi)\n phi = acos_with_shift(u / self.ecc, shift) # type: ignore [assignment]\n return phi + self.phi0\n\n def __call__(self, t: \"Collection[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Give a DataFrame of tau, u, r and phi from t.\"\"\"\n tau = self.tau(t)\n u = self.u_of_tau(tau)\n r = self.r_of_u(u)\n phi = self.phi_of_u_tau(u, tau)\n\n return pd.DataFrame({\"t\": t, \"tau\": tau, \"u\": u, \"r\": r, \"phi\": phi})\nalpha: float property","text":"Alpha \\(\\alpha\\) from the system specification.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.angular_mom","title":"angular_mom: float property","text":"Angular momentum \\(l\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.ecc","title":"ecc: float cached property","text":"Conic section eccentricity of the Kepler problem.
\\[ e = \\sqrt{1 + \\frac{2El}{\\alpha^2 m}}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.ene","title":"ene: float property","text":"Energy \\(E\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.mass","title":"mass: float property","text":"Mass \\(m\\) from the system specification.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.parameter","title":"parameter: float cached property","text":"Conic section semi-latus rectum of the Kepler problem.
\\[ p = \\frac{l^2}{\\alpha m}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.period","title":"period: float cached property","text":"Period \\(T\\) of the Kepler problem.
For \\(E < 0\\), $$ T = \\pi \\alpha \\sqrt{-\\frac{m}{2E^3}}\\,. $$
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.period_in_tau","title":"period_in_tau: float cached property","text":"Period in the scaled time tau.
\\[ T_\\tau = \\frac{2\\pi}{(1-e^2)^\\frac{3}{2}}\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.phi0","title":"phi0: float property","text":"phi0 \\(\\phi_0\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.t0","title":"t0: float property","text":"t0 \\(t_0\\) of the Kepler problem.
"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.t_to_tau_factor","title":"t_to_tau_factor: float property","text":"Scale factor from t to tau.
\\[ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. \\]"},{"location":"references/models/kepler_problem/model/#hamilflow.models.kepler_problem.model.Kepler2D.__call__","title":"__call__(t)","text":"Give a DataFrame of tau, u, r and phi from t.
Source code inhamilflow/models/kepler_problem/model.py def __call__(self, t: \"Collection[float] | npt.ArrayLike\") -> pd.DataFrame:\n \"\"\"Give a DataFrame of tau, u, r and phi from t.\"\"\"\n tau = self.tau(t)\n u = self.u_of_tau(tau)\n r = self.r_of_u(u)\n phi = self.phi_of_u_tau(u, tau)\n\n return pd.DataFrame({\"t\": t, \"tau\": tau, \"u\": u, \"r\": r, \"phi\": phi})\nfrom_geometry(system, geometries) classmethod","text":"Alternative initialiser from system and geometry specifications.
Given the eccentricity \\(e\\) and the conic section semi-latus rectum \\(p\\), \\(\\(l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.\\)\\)
Parameters:
Name Type Description Defaultsystem Mapping[str, float] the Kepler problem system specification
requiredgeometries Mapping[str, bool | float] geometric specifications positive_angular_mom: whether the angular momentum is positive ecc: eccentricity of the conic section parameter: semi-latus rectum of the conic section
hamilflow/models/kepler_problem/model.py @classmethod\ndef from_geometry(\n cls,\n system: \"Mapping[str, float]\",\n geometries: \"Mapping[str, bool | float]\",\n) -> \"Self\":\n r\"\"\"Alternative initialiser from system and geometry specifications.\n\n Given the eccentricity $e$ and the conic section semi-latus rectum $p$,\n $$l = \\pm \\sqrt{mp}\\,,\\quad E = (e^2-1) \\left|E_\\text{min}\\right|\\,.$$\n\n :param system: the Kepler problem system specification\n :param geometries: geometric specifications\n `positive_angular_mom`: whether the angular momentum is positive\n `ecc`: eccentricity of the conic section\n `parameter`: semi-latus rectum of the conic section\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n positive_angular_mom = bool(geometries[\"positive_angular_mom\"])\n ecc, parameter = (float(geometries[k]) for k in [\"ecc\", \"parameter\"])\n abs_angular_mom = math.sqrt(mass * parameter * alpha)\n # abs_minimal_ene = alpha / 2 / parameter: numerically unstable\n abs_minimal_ene = abs(cls.minimal_ene(abs_angular_mom, system))\n ene = (ecc**2 - 1) * abs_minimal_ene\n fi = {\n \"ene\": ene,\n \"angular_mom\": (\n abs_angular_mom if positive_angular_mom else -abs_angular_mom\n ),\n }\n return cls(system, fi)\nminimal_ene(angular_mom, system) staticmethod","text":"Minimal possible energy from the system specification and an angular momentum.
\\[ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. \\]Parameters:
Name Type Description Defaultangular_mom float angular momentum
requiredsystem Mapping[str, float] system specification
requiredReturns:
Type Descriptionfloat minimal possible energy
Source code inhamilflow/models/kepler_problem/model.py @staticmethod\ndef minimal_ene(\n angular_mom: float,\n system: \"Mapping[str, float]\",\n) -> float:\n r\"\"\"Minimal possible energy from the system specification and an angular momentum.\n\n $$ E_\\text{min} = -\\frac{m\\alpha^2}{2l^2}\\,. $$\n\n :param angular_mom: angular momentum\n :param system: system specification\n :return: minimal possible energy\n \"\"\"\n mass, alpha = system[\"mass\"], system[\"alpha\"]\n return -mass * alpha**2 / (2 * angular_mom**2)\nphi_of_u_tau(u, tau)","text":"Give the angular phi from u and tau.
For \\(e = 0\\), $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$ For \\(e > 0\\), $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$
Source code inhamilflow/models/kepler_problem/model.py def phi_of_u_tau(\n self,\n u: \"Collection[float] | npt.ArrayLike\",\n tau: \"Collection[float] | npt.ArrayLike\",\n) -> \"npt.ArrayLike\":\n r\"\"\"Give the angular phi from u and tau.\n\n For $e = 0$,\n $$ \\phi - \\phi_0 = 2\\pi \\frac{\\tau}{T_\\tau}\\,; $$\n For $e > 0$,\n $$ \\cos(\\phi - \\phi_0) = \\frac{u}{e}\\,. $$\n \"\"\"\n u, tau = np.asarray(u), np.asarray(tau)\n if self.ecc == 0:\n phi = 2 * math.pi * tau / self.period_in_tau\n else:\n if self.ecc < 1:\n shift = tau / self.period_in_tau\n else:\n shift = np.where(tau >= 0, 0, -np.pi)\n phi = acos_with_shift(u / self.ecc, shift) # type: ignore [assignment]\n return phi + self.phi0\nr_of_u(u)","text":"Give the radial r from u.
\\[ r = \\frac{p}{u+1}\\,. \\] Source code inhamilflow/models/kepler_problem/model.py def r_of_u(self, u: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the radial r from u.\n\n $$ r = \\frac{p}{u+1}\\,. $$\n \"\"\"\n return self.parameter / (np.asarray(u) + 1)\ntau(t)","text":"Give the scaled time tau from t.
\\[ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. \\] Source code inhamilflow/models/kepler_problem/model.py def tau(self, t: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n r\"\"\"Give the scaled time tau from t.\n\n $$ \\tau = \\frac{\\alpha^2 m}{|l|^3} (t-t_0)\\,. $$\n \"\"\"\n return (np.asarray(t) - self.t0) * self.t_to_tau_factor\nu_of_tau(tau)","text":"Give the convenient radial inverse u from tau.
Source code inhamilflow/models/kepler_problem/model.py def u_of_tau(self, tau: \"Collection[float] | npt.ArrayLike\") -> \"npt.ArrayLike\":\n \"\"\"Give the convenient radial inverse u from tau.\"\"\"\n tau = np.asarray(tau)\n if self.ecc == 0:\n return np.zeros(tau.shape)\n else:\n if self.ecc < 1:\n p = self.period_in_tau\n r = tau % p\n tau = np.where(r <= p / 2, r, p - r)\n else:\n tau = np.abs(tau)\n return u_of_tau(self.ecc, tau) # type: ignore [arg-type]\nKepler2DFI","text":" Bases: BaseModel
The first integrals for a Kepler problem.
Attributes:
Name Type Descriptionene float the energy \\(E\\)
angular_mom float the angular momentum \\(l\\)
t0 float the time \\(t_0\\) at which the radial position is closest to 0, defaults to 0
phi0 float the angle \\(\\phi_0\\) at which the radial position is closest to 0, defaults to 0
Source code inhamilflow/models/kepler_problem/model.py class Kepler2DFI(BaseModel):\n r\"\"\"The first integrals for a Kepler problem.\n\n :cvar ene: the energy $E$\n :cvar angular_mom: the angular momentum $l$\n :cvar t0: the time $t_0$ at which the radial position is closest to 0, defaults to 0\n :cvar phi0: the angle $\\phi_0$ at which the radial position is closest to 0, defaults to 0\n \"\"\"\n\n ene: float = Field()\n angular_mom: float = Field()\n t0: float = Field(default=0)\n phi0: float = Field(ge=0, lt=2 * math.pi, default=0)\n\n # TODO process angular momentum = 0\n @field_validator(\"angular_mom\")\n @classmethod\n def _angular_mom_non_zero(cls, v: float) -> float:\n if v == 0:\n msg = \"Only non-zero angular momenta are supported\"\n raise NotImplementedError(msg)\n return v\nKepler2DSystem","text":" Bases: BaseModel
Definition of the Kepler problem.
Potential:
\\[ V(r) = - \\frac{\\alpha}{r}. \\]For reference, if an object is orbiting our Sun, the constant \\(\\alpha = G M_{\\odot} ~ 1.327 \\times 10^{20} \\mathrm{m}^3/\\mathrm{s}^2\\) in SI, which is also called 1 TCB, or 1 solar mass parameter. For computational stability, we recommend using TCB as the unit instead of the large SI values.
Units
When specifying the parameters of the system, be ware of the consistency of the units.
Attributes:
Name Type Descriptionalpha float the proportional constant of the potential energy.
mass float the mass of the orbiting object
Source code inhamilflow/models/kepler_problem/model.py class Kepler2DSystem(BaseModel):\n r\"\"\"Definition of the Kepler problem.\n\n Potential:\n\n $$\n V(r) = - \\frac{\\alpha}{r}.\n $$\n\n For reference, if an object is orbiting our Sun, the constant\n $\\alpha = G M_{\\odot} ~ 1.327 \\times 10^{20} \\mathrm{m}^3/\\mathrm{s}^2$ in SI,\n which is also called 1 TCB, or 1 solar mass parameter. For computational stability, we recommend using\n TCB as the unit instead of the large SI values.\n\n !!! note \"Units\"\n\n When specifying the parameters of the system, be ware of the consistency of the units.\n\n :cvar alpha: the proportional constant of the potential energy.\n :cvar mass: the mass of the orbiting object\n \"\"\"\n\n # TODO add repulsive alpha < 0\n alpha: float = Field(gt=0, default=1.0)\n mass: float = Field(gt=0, default=1.0)\nNumerics for the Kepler problem.
"},{"location":"references/models/kepler_problem/numerics/#hamilflow.models.kepler_problem.numerics.nsolve_u_from_tau_bisect","title":"nsolve_u_from_tau_bisect(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 0, ecc != 1
requiredtau ArrayLike scaled time
requiredReturns:
Type Descriptionlist[OptimizeResult] numeric OptimizeResult from scipy
Source code inhamilflow/models/kepler_problem/numerics.py def nsolve_u_from_tau_bisect(ecc: float, tau: \"npt.ArrayLike\") -> list[OptimizeResult]:\n \"\"\"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.\n\n :param ecc: eccentricity, ecc > 0, ecc != 1\n :param tau: scaled time\n :return: numeric OptimizeResult from scipy\n \"\"\"\n tau_s = np.asarray(tau).reshape(-1)\n if 0 < ecc < 1:\n tau_of_u = tau_of_u_elliptic\n elif ecc > 1:\n tau_of_u = tau_of_u_hyperbolic\n else:\n msg = f\"Expect ecc > 0, ecc != 1, got {ecc}\"\n raise ValueError(msg)\n\n def f(u: float, tau: float) -> np.float64:\n return (\n np.finfo(np.float64).max if u == -1 else np.float64(tau_of_u(ecc, u) - tau)\n )\n\n return [toms748(f, max(-1, -ecc), ecc, (ta,), 2, full_output=True) for ta in tau_s]\nnsolve_u_from_tau_newton(ecc, tau)","text":"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc > 0, ecc != 1
requiredtau ArrayLike scaled time
requiredReturns:
Type DescriptionOptimizeResult numeric OptimizeResult from scipy
Raises:
Type DescriptionValueError when ecc is invalid
hamilflow/models/kepler_problem/numerics.py def nsolve_u_from_tau_newton(ecc: float, tau: \"npt.ArrayLike\") -> OptimizeResult:\n \"\"\"Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.\n\n :param ecc: eccentricity, ecc > 0, ecc != 1\n :param tau: scaled time\n :raises ValueError: when `ecc` is invalid\n :return: numeric OptimizeResult from scipy\n \"\"\"\n tau = np.asarray(tau)\n if 0 < ecc < 1:\n tau_of_u = tau_of_u_elliptic\n u0 = _u0_elliptic(ecc, tau)\n elif ecc > 1:\n tau_of_u = tau_of_u_hyperbolic\n u0 = _u0_hyperbolic(ecc, tau)\n else:\n msg = f\"Expect ecc > 0, ecc != 1, got {ecc}\"\n raise ValueError(msg)\n\n def f(u: float, tau: \"npt.NDArray[np.float64]\") -> \"npt.NDArray[np.float64]\":\n return tau_of_u(ecc, u) - tau\n\n def fprime(\n u: float,\n tau: \"npt.ArrayLike\", # noqa: ARG001\n ) -> \"npt.NDArray[np.float64]\":\n return tau_of_u_prime(ecc, u)\n\n def fprime2(\n u: float,\n tau: \"npt.ArrayLike\", # noqa: ARG001\n ) -> \"npt.NDArray[np.float64]\":\n return tau_of_u_prime2(ecc, u)\n\n return newton(f, u0, fprime, (tau,), fprime2=fprime2, full_output=True, disp=False)\nu_of_tau(ecc, tau, method='bisect')","text":"Calculate the convenient radial inverse u from tau, using numeric methods.
Parameters:
Name Type Description Defaultecc float eccentricity, ecc >= 0
requiredtau ArrayLike scaled time
requiredmethod Literal['bisect', 'newton'] \"newton\" or \"bisect\" numeric methods, defaults to \"bisect\"
'bisect' Returns:
Type DescriptionNDArray[float64] convenient radial inverse u
Raises:
Type DescriptionValueError when method is invalid
ValueError when ecc is invalid
hamilflow/models/kepler_problem/numerics.py def u_of_tau(\n ecc: float,\n tau: \"npt.ArrayLike\",\n method: Literal[\"bisect\", \"newton\"] = \"bisect\",\n) -> \"npt.NDArray[np.float64]\":\n \"\"\"Calculate the convenient radial inverse u from tau, using numeric methods.\n\n :param ecc: eccentricity, ecc >= 0\n :param tau: scaled time\n :param method: \"newton\" or \"bisect\" numeric methods, defaults to \"bisect\"\n :raises ValueError: when `method` is invalid\n :raises ValueError: when `ecc` is invalid\n :return: convenient radial inverse u\n \"\"\"\n tau = np.asarray(tau)\n if ecc == 0:\n return np.zeros(tau.shape)\n elif ecc == 1:\n return esolve_u_from_tau_parabolic(ecc, tau)\n elif ecc > 0:\n match method:\n case \"bisect\":\n return np.array([s[0] for s in nsolve_u_from_tau_bisect(ecc, tau)])\n case \"newton\":\n return nsolve_u_from_tau_newton(ecc, tau)[0]\n case _:\n msg = f\"Expect 'bisect' or 'newton', got {method}\"\n raise ValueError(msg)\n else:\n msg = f\"Expect ecc >= 0, got {ecc}\"\n raise ValueError(msg)\nWe provide some tutorials to help you get started.
For easier reading, we following the following notations whenever possible.
Notation Description \\(\\mathcal L\\) Lagrangian \\(L\\) Angular Momentum \\(V\\) Potential"},{"location":"tutorials/brownian_motion/","title":"Brownian Motion","text":"In\u00a0[1]: Copied!import plotly.express as px\n\nfrom hamilflow.models.brownian_motion import BrownianMotion\nimport plotly.express as px from hamilflow.models.brownian_motion import BrownianMotion In\u00a0[2]: Copied!
bm_1d = BrownianMotion(\n system={\n \"sigma\": 1,\n \"delta_t\": 1,\n },\n initial_condition={\"x0\": 0},\n)\n bm_1d = BrownianMotion( system={ \"sigma\": 1, \"delta_t\": 1, }, initial_condition={\"x0\": 0}, ) Call the model to generate 1000 steps.
In\u00a0[3]: Copied!df_1d = bm_1d.generate_from(n_steps=1000)\ndf_1d = bm_1d.generate_from(n_steps=1000) In\u00a0[4]: Copied!
px.line(df_1d, x=\"t\", y=\"x_0\")\npx.line(df_1d, x=\"t\", y=\"x_0\") In\u00a0[5]: Copied!
bm_2d = BrownianMotion(\n system={\n \"sigma\": 1,\n \"delta_t\": 1,\n },\n initial_condition={\"x0\": [0, 0]},\n)\n bm_2d = BrownianMotion( system={ \"sigma\": 1, \"delta_t\": 1, }, initial_condition={\"x0\": [0, 0]}, ) We call the model to generate 1000 steps.
In\u00a0[6]: Copied!df_2d = bm_2d.generate_from(n_steps=500)\ndf_2d = bm_2d.generate_from(n_steps=500) In\u00a0[7]: Copied!
(\n px.scatter(df_2d, x=\"x_0\", y=\"x_1\", color=\"t\")\n .update_traces(\n mode=\"lines+markers\",\n marker={\n \"size\": 2.5,\n },\n line={\"width\": 1},\n )\n .update_yaxes(\n scaleanchor=\"x\",\n scaleratio=1,\n )\n)\n ( px.scatter(df_2d, x=\"x_0\", y=\"x_1\", color=\"t\") .update_traces( mode=\"lines+markers\", marker={ \"size\": 2.5, }, line={\"width\": 1}, ) .update_yaxes( scaleanchor=\"x\", scaleratio=1, ) ) In\u00a0[\u00a0]: Copied! \n"},{"location":"tutorials/brownian_motion/#brownian-motion","title":"Brownian Motion\u00b6","text":"
Brownian motion describes motion of small particles with stochastic forces applied to them. The math of Brownian motion can be modeled with Wiener process. In this tutorial, we take a simple form of the model and treat the stochastic forces as Gaussian.
For consistency, we always use $\\mathbf x$ for displacement, and $t$ for steps. The model we are using is
$$ \\begin{align} \\mathbf x(t + \\mathrm dt) &= \\mathbf x(t) + \\mathcal{N}(\\mu=0, \\sigma=\\sigma \\sqrt{\\mathrm d t}) \\end{align} $$
Read more:
hamilflow implemented a Brownian motion model called BrownianMotion.
Our BrownianMotion model calculates the dimension of the space based on the dimension of the initial condition $x_0$. To create a 2D Brownian motion model, we need the initial condition to be length 2.
import math\n\nimport numpy as np\nfrom plotly import express as px\n\nfrom hamilflow.models.harmonic_oscillator import ComplexSimpleHarmonicOscillator\nimport math import numpy as np from plotly import express as px from hamilflow.models.harmonic_oscillator import ComplexSimpleHarmonicOscillator In\u00a0[2]: Copied!
t = np.linspace(0, 3, 257)\nsystem_specs = {\"omega\": 2 * math.pi}\n t = np.linspace(0, 3, 257) system_specs = {\"omega\": 2 * math.pi} In\u00a0[3]: Copied! csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (1, 0)})\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (1, 0)}) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, ) In\u00a0[4]: Copied! csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (0, 1)})\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator(system_specs, initial_condition={\"x0\": (0, 1)}) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, ) In\u00a0[5]: Copied! csho = ComplexSimpleHarmonicOscillator(\n system_specs,\n initial_condition={\"x0\": (math.cos(math.pi / 12), math.sin(math.pi / 12))},\n)\n\ndf = csho(t)\n\narr_z = df[\"z\"].to_numpy(copy=False)\n\npx.line_3d(\n x=arr_z.real,\n y=arr_z.imag,\n z=t,\n labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"},\n)\n csho = ComplexSimpleHarmonicOscillator( system_specs, initial_condition={\"x0\": (math.cos(math.pi / 12), math.sin(math.pi / 12))}, ) df = csho(t) arr_z = df[\"z\"].to_numpy(copy=False) px.line_3d( x=arr_z.real, y=arr_z.imag, z=t, labels={\"x\": \"real\", \"y\": \"imag\", \"z\": \"time\"}, )"},{"location":"tutorials/complex_harmonic_oscillator/#complex-harmonic-oscillator","title":"Complex Harmonic Oscillator\u00b6","text":"In this tutorial, we show a few interesting examples of the complex simple harmonic oscillator.
"},{"location":"tutorials/complex_harmonic_oscillator/#basic-setups","title":"Basic setups\u00b6","text":""},{"location":"tutorials/complex_harmonic_oscillator/#positive-frequency-circular-polarised-mode","title":"Positive-frequency, circular-polarised mode\u00b6","text":"Also known as the left-rotating mode.
"},{"location":"tutorials/complex_harmonic_oscillator/#negative-frequency-circular-polarised-mode","title":"Negative-frequency, circular-polarised mode\u00b6","text":"Also known as the right-rotating mode.
"},{"location":"tutorials/complex_harmonic_oscillator/#positive-frequency-elliptic-polarised-mode","title":"Positive-frequency, elliptic-polarised mode\u00b6","text":""},{"location":"tutorials/complex_harmonic_oscillator/#end-of-notebook","title":"End of Notebook\u00b6","text":""},{"location":"tutorials/harmonic_oscillator/","title":"Harmonic Oscillators","text":"In\u00a0[1]: Copied!import pandas as pd\nimport plotly.express as px\n\nfrom hamilflow.models.harmonic_oscillator import (\n DampedHarmonicOscillator,\n SimpleHarmonicOscillator,\n)\nimport pandas as pd import plotly.express as px from hamilflow.models.harmonic_oscillator import ( DampedHarmonicOscillator, SimpleHarmonicOscillator, ) In\u00a0[2]: Copied!
n_periods = 3\nn_samples_per_period = 200\nn_periods = 3 n_samples_per_period = 200 In\u00a0[3]: Copied!
sho_omega = 0.5\n\nsho = SimpleHarmonicOscillator(system={\"omega\": sho_omega})\n sho_omega = 0.5 sho = SimpleHarmonicOscillator(system={\"omega\": sho_omega}) In\u00a0[4]: Copied! df_sho = sho(n_periods=n_periods, n_samples_per_period=n_samples_per_period)\ndf_sho.head()\ndf_sho = sho(n_periods=n_periods, n_samples_per_period=n_samples_per_period) df_sho.head() Out[4]: t x 0 0.000000 1.000000 1 0.062832 0.999507 2 0.125664 0.998027 3 0.188496 0.995562 4 0.251327 0.992115 In\u00a0[5]: Copied!
px.line(\n df_sho,\n x=\"t\",\n y=\"x\",\n title=rf\"Simple Harmonic Oscillator (omega = {sho_omega})\",\n labels={\n \"x\": r\"Displacement $x(t)$\",\n \"t\": r\"$t$\",\n },\n)\n px.line( df_sho, x=\"t\", y=\"x\", title=rf\"Simple Harmonic Oscillator (omega = {sho_omega})\", labels={ \"x\": r\"Displacement $x(t)$\", \"t\": r\"$t$\", }, ) In\u00a0[6]: Copied! dho_systems = {\n \"Underdamped\": {\"omega\": 0.5, \"zeta\": 0.2},\n \"Critical Damped\": {\"omega\": 0.5, \"zeta\": 1},\n \"Overdamped\": {\n \"omega\": 0.5,\n \"zeta\": 1.2,\n },\n}\n\ndfs_dho = []\n\nfor s_name, s in dho_systems.items():\n dfs_dho.append(\n DampedHarmonicOscillator(system=s)(\n n_periods=n_periods,\n n_samples_per_period=n_samples_per_period,\n ).assign(\n system=rf\"{s_name} (omega = {s.get('omega')}, zeta = {s.get('zeta')})\",\n ),\n )\n\nfig = px.line(\n pd.concat(dfs_dho),\n x=\"t\",\n y=\"x\",\n color=\"system\",\n title=r\"Damped Harmonic Oscillator\",\n labels={\n \"x\": r\"Displacement $x(t)$\",\n \"t\": r\"$t$\",\n },\n)\nfig.update_layout(legend={\"yanchor\": \"top\", \"y\": -0.2, \"xanchor\": \"left\", \"x\": 0})\n dho_systems = { \"Underdamped\": {\"omega\": 0.5, \"zeta\": 0.2}, \"Critical Damped\": {\"omega\": 0.5, \"zeta\": 1}, \"Overdamped\": { \"omega\": 0.5, \"zeta\": 1.2, }, } dfs_dho = [] for s_name, s in dho_systems.items(): dfs_dho.append( DampedHarmonicOscillator(system=s)( n_periods=n_periods, n_samples_per_period=n_samples_per_period, ).assign( system=rf\"{s_name} (omega = {s.get('omega')}, zeta = {s.get('zeta')})\", ), ) fig = px.line( pd.concat(dfs_dho), x=\"t\", y=\"x\", color=\"system\", title=r\"Damped Harmonic Oscillator\", labels={ \"x\": r\"Displacement $x(t)$\", \"t\": r\"$t$\", }, ) fig.update_layout(legend={\"yanchor\": \"top\", \"y\": -0.2, \"xanchor\": \"left\", \"x\": 0}) In\u00a0[\u00a0]: Copied! \n"},{"location":"tutorials/harmonic_oscillator/#harmonic-oscillators","title":"Harmonic Oscillators\u00b6","text":"
In this tutorial, we demo how to generate data of harmonic oscillators.
"},{"location":"tutorials/harmonic_oscillator/#simple-harmonic-oscillator","title":"Simple Harmonic Oscillator\u00b6","text":"For an simple harmonic oscillator, the action of a simple harmonic oscillator is
$$S_L[x] = \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} m \\dot x^2 - \\frac{1}{2} m \\omega^2 x^2 \\right\\}\\,,$$
where the least action principle leads to the following equation of motion,
$$ \\ddot x + \\omega^2 x = 0\\,. $$
A simple harmonic oscillator is a periodic motion.
"},{"location":"tutorials/harmonic_oscillator/#damped-harmonic-oscillator","title":"Damped Harmonic Oscillator\u00b6","text":"A damped harmonic oscillator is a simple harmonic oscillator with damping force that is proportional to its velocity,
$$ \\ddot x + \\omega^2 x = - 2\\xi\\omega \\dot x\\,. $$
In this section, we demonstrate three scenarios of a damped harmonic oscillator.
"},{"location":"tutorials/harmonic_oscillator_chain/","title":"Harmonic oscillator circle","text":"In\u00a0[1]: Copied!import math\n\nimport numpy as np\nfrom plotly import express as px\n\nfrom hamilflow.models.harmonic_oscillator_chain import HarmonicOscillatorsChain\nimport math import numpy as np from plotly import express as px from hamilflow.models.harmonic_oscillator_chain import HarmonicOscillatorsChain In\u00a0[2]: Copied!
t = np.linspace(0, 3, 257)\nomega = 2 * math.pi\npattern_x = r\"x\\d+\"\nt = np.linspace(0, 3, 257) omega = 2 * math.pi pattern_x = r\"x\\d+\" In\u00a0[3]: Copied!
ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] # type: ignore [index]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] # type: ignore [index] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[4]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(2 * math.pi, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(2 * math.pi, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[5]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 0)}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 4\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (0, 0)}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 4 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[6]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, True)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, True) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[7]: Copied! ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, False)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 0}, {\"amp\": (1, 1)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, False) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"}) In\u00a0[8]: Copied! ics = [{\"x0\": 0, \"v0\": 2}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5\nhoc = HarmonicOscillatorsChain(omega, ics, False)\n\ndf_res = hoc(t)\ndf_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)]\npx.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})\n ics = [{\"x0\": 0, \"v0\": 2}, {\"amp\": (1, 0)}] + [{\"amp\": (0, 0)}] * 5 hoc = HarmonicOscillatorsChain(omega, ics, False) df_res = hoc(t) df_x_wide = df_res.loc[:, df_res.columns.str.match(pattern_x)] px.imshow(df_x_wide, origin=\"lower\", labels={\"y\": \"t\"})"},{"location":"tutorials/harmonic_oscillator_chain/#harmonic-oscillator-circle","title":"Harmonic oscillator circle\u00b6","text":"In this tutorial, we show a few interesting examples of the harmonic oscillator circle.
"},{"location":"tutorials/harmonic_oscillator_chain/#basic-setups","title":"Basic setups\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-right-moving-wave","title":"Fundamental right-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-left-moving-wave","title":"Fundamental left-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#faster-right-moving-wave","title":"Faster right-moving wave\u00b6","text":"Also known as the first harmonic.
"},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-stationary-wave-odd-dof","title":"Fundamental stationary wave, odd dof\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#fundamental-stationary-wave-even-dof","title":"Fundamental stationary wave, even dof\u00b6","text":"There are stationary nodes at $i = 3, 9$.
"},{"location":"tutorials/harmonic_oscillator_chain/#linearly-moving-chain-and-fundamental-right-moving-wave","title":"Linearly moving chain and fundamental right-moving wave\u00b6","text":""},{"location":"tutorials/harmonic_oscillator_chain/#mathematical-physical-discription","title":"Mathematical-physical discription\u00b6","text":"A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action $$S_L[x_i] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\sum_{i=0}^{N-1} \\frac{1}{2}m \\dot x_i^2 - \\frac{1}{2}m\\omega^2\\left(x_i - x_{i+1}\\right)^2 \\right\\\\}\\\\,,$$ where $x_N \\coloneqq x_0$.
This system can be solved in terms of travelling waves, obtained by discrete Fourier transform.
We can complexify the system $$S_L[x_i] = S_L[x_i, \\phi_j] \\equiv S_L[X^\\ast_i, X_j] = \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot X^\\ast_i \\delta_{ij} \\dot X_j - \\frac{1}{2}m X^\\ast_i A_{ij} X_j\\right\\\\}\\\\,,$$ where $A_{ij} / \\omega^2$ is equal to $(-2)$ if $i=j$, $1$ if $|i-j|=1$ or $|i-j|=N$, and $0$ otherwise; $X_i \\coloneqq x_i \\mathbb{e}^{-\\phi_i}$, $X^\\ast_i \\coloneqq x_i \\mathbb{e}^{+\\phi_i}$.
$A_{ij}$ can be diagonalised by the inverse discrete Fourier transform $$X_i = (F^{-1})_{ik} Y_k = \\frac{1}{\\sqrt{N}}\\sum_k \\mathbb{e}^{i \\frac{2\\mathbb{\\pi}}{N} k\\mathbb{i}} Y_k\\\\,.$$
Calculating gives $$S_L[X^\\ast_i, X_j] = S_L[Y^\\ast_i, Y_j] = \\sum_{k=0}^{N-1} \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot Y^\\ast_k \\dot Y_k - \\frac{1}{2}m \\omega^2\\cdot4\\sin^2\\frac{2\\mathbb{\\pi}k}{N} Y^\\ast_k Y_k\\right\\\\}\\\\,.$$ We can arrive at an action for the Fourier modes $$S_L[Y, Y^*] = \\sum_{k=0}^{N-1} \\int_{t_0}^{t_1}\\mathbb{d} t \\left\\\\{ \\frac{1}{2}m \\dot Y_k^* Y_k - \\frac{1}{2}m \\omega^2\\cdot4\\sin^2\\frac{2\\mathbb{\\pi}k}{N} Y_k^* Y_k\\right\\\\}\\\\,.$$
The origional system can then be solved by $N$ complex oscillators.
Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy $Y_k = Y^*_{-k \\mod N}$, see Wikipedia.
"},{"location":"tutorials/harmonic_oscillator_chain/#end-of-notebook","title":"End of Notebook\u00b6","text":""},{"location":"tutorials/kepler_problem/","title":"Kepler Problem","text":"In\u00a0[1]: Copied!import numpy as np\nimport pandas as pd\nimport plotly.express as px\nimport plotly.graph_objects as go\n\nfrom hamilflow.models.kepler_problem import Kepler2D\nimport numpy as np import pandas as pd import plotly.express as px import plotly.graph_objects as go from hamilflow.models.kepler_problem import Kepler2D
To make it easy to use the Kepler system, we implemented an interface Kepler2D.from_geometry to specify the configuration using the geometry of the orbit. The geometry of the orbit requires three parameters:
positive_angular_mom: whether the angular momentum is positiveecc: the eccentricity of the orbitparameter: the semi-latus rectum of the orbits, for circular orbits, this is the radiussystem = {\n \"alpha\": 1.0,\n \"mass\": 1.0,\n}\n\nk2d = Kepler2D.from_geometry(\n system=system,\n geometries={\n \"positive_angular_mom\": True,\n \"ecc\": 0,\n \"parameter\": 1.0,\n },\n)\n system = { \"alpha\": 1.0, \"mass\": 1.0, } k2d = Kepler2D.from_geometry( system=system, geometries={ \"positive_angular_mom\": True, \"ecc\": 0, \"parameter\": 1.0, }, ) We define the time steps for a Kepler system. In this example, we create 401 time steps, also with some some negative timestamps for a more symmetric trajectory.
In\u00a0[3]: Copied!# t = np.linspace(-20, 20, 401)\nt = np.arange(-20, 20, 0.125)\n# t = np.linspace(-20, 20, 401) t = np.arange(-20, 20, 0.125)
Orbit time series data is generated by calling the object k2d, as the data generation is done in the __call__ method.
df_p_1 = k2d(t=t)\ndf_p_1 # noqa: B018\ndf_p_1 = k2d(t=t) df_p_1 # noqa: B018 Out[4]: t tau u r phi 0 -20.000 -20.000 0.0 1.0 -20.000 1 -19.875 -19.875 0.0 1.0 -19.875 2 -19.750 -19.750 0.0 1.0 -19.750 3 -19.625 -19.625 0.0 1.0 -19.625 4 -19.500 -19.500 0.0 1.0 -19.500 ... ... ... ... ... ... 315 19.375 19.375 0.0 1.0 19.375 316 19.500 19.500 0.0 1.0 19.500 317 19.625 19.625 0.0 1.0 19.625 318 19.750 19.750 0.0 1.0 19.750 319 19.875 19.875 0.0 1.0 19.875
320 rows \u00d7 5 columns
In\u00a0[5]: Copied!def visualize_orbit(dataframe: pd.DataFrame) -> go.Figure:\n \"\"\"Plot out the trajectory in a polar coordinate.\n\n :param dataframe: dataframe containing the trajectory data\n \"\"\"\n df_chart = dataframe.assign(\n phi_degree=dataframe.phi / (2 * np.pi) * 360,\n )\n range_r = 0, df_chart.r.max() * 1.1\n fig = px.scatter_polar(\n df_chart,\n r=\"r\",\n theta=\"phi_degree\",\n hover_data=[\"t\", \"r\", \"phi\"],\n start_angle=0,\n animation_frame=\"t\",\n color_discrete_sequence=[\"red\"],\n range_r=range_r,\n )\n\n trajectory_traces = list(\n px.line_polar(\n df_chart,\n r=\"r\",\n theta=\"phi_degree\",\n hover_data=[\"t\", \"r\", \"phi\"],\n start_angle=0,\n range_r=range_r,\n ).select_traces(),\n )\n\n fig.update_layout(\n polar={\n \"angularaxis_thetaunit\": \"radians\",\n },\n )\n\n fig.add_traces(trajectory_traces)\n\n fig.layout.updatemenus[0].buttons[0].args[1][\"frame\"][\"duration\"] = 50\n fig.layout.updatemenus[0].buttons[0].args[1][\"transition\"][\"duration\"] = 50\n\n return fig\n def visualize_orbit(dataframe: pd.DataFrame) -> go.Figure: \"\"\"Plot out the trajectory in a polar coordinate. :param dataframe: dataframe containing the trajectory data \"\"\" df_chart = dataframe.assign( phi_degree=dataframe.phi / (2 * np.pi) * 360, ) range_r = 0, df_chart.r.max() * 1.1 fig = px.scatter_polar( df_chart, r=\"r\", theta=\"phi_degree\", hover_data=[\"t\", \"r\", \"phi\"], start_angle=0, animation_frame=\"t\", color_discrete_sequence=[\"red\"], range_r=range_r, ) trajectory_traces = list( px.line_polar( df_chart, r=\"r\", theta=\"phi_degree\", hover_data=[\"t\", \"r\", \"phi\"], start_angle=0, range_r=range_r, ).select_traces(), ) fig.update_layout( polar={ \"angularaxis_thetaunit\": \"radians\", }, ) fig.add_traces(trajectory_traces) fig.layout.updatemenus[0].buttons[0].args[1][\"frame\"][\"duration\"] = 50 fig.layout.updatemenus[0].buttons[0].args[1][\"transition\"][\"duration\"] = 50 return fig In\u00a0[6]: Copied! fig = visualize_orbit(df_p_1)\nfig.show()\nfig = visualize_orbit(df_p_1) fig.show()
We also plot out the ellipse, hyperbolic, parabolic orbits.
In\u00a0[7]: Copied!visualize_orbit(\n Kepler2D.from_geometry(\n system=system,\n geometries={\n \"positive_angular_mom\": True,\n \"ecc\": 0.5,\n \"parameter\": 1.0,\n },\n )(t=t),\n).show()\n visualize_orbit( Kepler2D.from_geometry( system=system, geometries={ \"positive_angular_mom\": True, \"ecc\": 0.5, \"parameter\": 1.0, }, )(t=t), ).show() /home/runner/work/hamilflow/hamilflow/hamilflow/models/kepler_problem/dynamics.py:26: RuntimeWarning:\n\ndivide by zero encountered in divide\n\nIn\u00a0[8]: Copied!
visualize_orbit(\n Kepler2D.from_geometry(\n system=system,\n geometries={\n \"positive_angular_mom\": True,\n \"ecc\": 1,\n \"parameter\": 1.0,\n },\n )(t=t),\n).show()\n visualize_orbit( Kepler2D.from_geometry( system=system, geometries={ \"positive_angular_mom\": True, \"ecc\": 1, \"parameter\": 1.0, }, )(t=t), ).show() In\u00a0[9]: Copied! visualize_orbit(\n Kepler2D.from_geometry(\n system=system,\n geometries={\n \"positive_angular_mom\": True,\n \"ecc\": 2,\n \"parameter\": 1.0,\n },\n )(t=np.arange(-5, 5, 0.125)),\n).show()\n visualize_orbit( Kepler2D.from_geometry( system=system, geometries={ \"positive_angular_mom\": True, \"ecc\": 2, \"parameter\": 1.0, }, )(t=np.arange(-5, 5, 0.125)), ).show() /home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1046: RuntimeWarning:\n\noverflow encountered in scalar add\n\n/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1046: RuntimeWarning:\n\ninvalid value encountered in scalar multiply\n\n/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1054: RuntimeWarning:\n\noverflow encountered in scalar add\n\n/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:995: RuntimeWarning:\n\noverflow encountered in divide\n\n
In this section, we briefly discuss the derivations of motion in a central field. Please refer to Mechanics: Vol 1 by Landau and Lifshitz for more details[^1].
The Lagrangian for an object in a central field is
$$ \\mathcal L = \\frac{1}{2} m ({\\dot r}^2 + r^2 {\\dot \\phi}^2) - V(r), $$
where $r$ and $\\phi$ are the polar coordinates, $m$ is the mass of the object, and $V(r)$ is the potential energy. The equations of motion are
$$ \\begin{align} \\frac{\\mathrm d r}{\\mathrm d t} &= \\sqrt{ \\frac{2}{m} (E - V(r)) - \\frac{L^2}{m^2 r^2} } \\\\ \\frac{\\mathrm d \\phi}{\\mathrm d t} &= \\frac{L}{m r^2}, \\end{align} $$
where $E$ is the total energy and $L$ is the angular momentum,
$$ \\begin{align} E =& \\frac{1}{2} m \\left(\\left(\\frac{\\mathrm d r}{ \\mathrm dt} \\right)^2 + r^2 \\left( \\frac{\\mathrm d r}{\\mathrm dt} \\right)^2\\right) + V(r) \\\\ L =& m r^2 \\frac{d\\phi}{dt} \\end{align} $$
Both $E$ and $L$ are conserved. We obtain the coordinates as a function of time by solving the two equations.
For a inverse-square force, the potential energy is
$$ V(r) = - \\frac{\\alpha}{r}, $$
where $\\alpha$ is a constant that specifies the strength of the force. For Newtonian gravity, $\\alpha = G m_0$ with $G$ being the gravitational constant.
First of all, we solve $\\phi(r)$,
$$ \\phi = \\cos^{-1}\\left( \\frac{L/r - m\\alpha/L}{2 m E + m \\alpha^2/L^2} \\right). $$
Define
$$ p = \\frac{L^2}{m\\alpha} $$
and
$$ e = \\sqrt{ 1 + \\frac{2 E L^2}{m \\alpha^2}}, $$
we rewrite the solution $\\phi(r)$ as
$$ r = \\frac{p}{1 + e \\cos{\\phi}}. $$
With the above relationship between $r$ and $\\phi$, and $\\frac{\\mathrm d \\phi}{\\mathrm d} = \\frac{L}{mr^2}$, we find that
$$ \\frac{m\\alpha^2}{L^3} \\mathrm d t = \\frac{1}{(1 + e \\cos{\\phi})^2} \\mathrm d \\phi. $$
The integration on the right hand side depends on the domain of $e$.
$$ \\int\\frac{1}{(1 + e \\cos{\\phi})^2} \\mathrm d \\phi = \\begin{cases} \\frac{1}{(1 - e^2)^{3/2}} \\left( 2 \\tan^{-1} \\sqrt{\\frac{1 - e}{1 + e}} \\tan\\frac{\\phi}{2} - \\frac{e\\sqrt{1 - e^2}\\sin\\phi}{1 + e\\cos\\phi} \\right), & \\text{if } e<1 \\\\ \\frac{1}{2}\\tan{\\frac{\\phi}{2}} + \\frac{1}{6}\\tan^3{\\frac{\\phi}{2}}, & \\text{if } e=1 \\\\ \\frac{1}{(e^2 - 1)^{3/2}}\\left( \\frac{e\\sqrt{e^2-1}\\sin\\phi}{1 + e\\cos\\phi} - \\ln \\left( \\frac{\\sqrt{1 + e} + \\sqrt{e-1}\\tan\\frac{\\phi}{2}}{\\sqrt{1 + e} - \\sqrt{e-1}\\tan\\frac{\\phi}{2} } \\right) \\right), & \\text{if } e> 1. \\end{cases} $$
The value of $t(\\phi)$ is easily obtained from the above formulae.
There exists many numerical methods to solve the Kepler orbits as functions of time, $r(t)$ and $\\phi(t)$. For our use case of the solutions, we choose to integrate the equation of motion directly.
References:
In this tutorial we will generate data for the Kepler problem.
"},{"location":"tutorials/kepler_problem/#usage","title":"Usage\u00b6","text":""},{"location":"tutorials/kepler_problem/#formalism","title":"Formalism\u00b6","text":""},{"location":"tutorials/pendulum/","title":"Pendulum","text":"In\u00a0[1]: Copied!import math\n\nimport plotly.express as px\n\nfrom hamilflow.models.pendulum import Pendulum\nimport math import plotly.express as px from hamilflow.models.pendulum import Pendulum In\u00a0[2]: Copied!
omega0 = 2 * math.pi\ntheta0 = math.pi / 3\n\nn_periods = 2**2\nn_samples_per_period = 2**8\nomega0 = 2 * math.pi theta0 = math.pi / 3 n_periods = 2**2 n_samples_per_period = 2**8 In\u00a0[3]: Copied!
pen = Pendulum(system=omega0, initial_condition=theta0)\npen = Pendulum(system=omega0, initial_condition=theta0) In\u00a0[4]: Copied!
df_pen = pen.generate_from(\n n_periods=n_periods,\n n_samples_per_period=n_samples_per_period,\n)\ndf_pen.head()\ndf_pen = pen.generate_from( n_periods=n_periods, n_samples_per_period=n_samples_per_period, ) df_pen.head() Out[4]: t x u 0 0.000000 1.047198 1.570796 1 0.004192 1.046897 1.593608 2 0.008384 1.045996 1.616424 3 0.012576 1.044494 1.639247 4 0.016768 1.042393 1.662082 In\u00a0[5]: Copied!
df_pen.describe()\ndf_pen.describe() Out[5]: t x u count 1024.000000 1.024000e+03 1024.000000 mean 2.144268 2.320193e-17 14.124895 std 1.239809 7.453532e-01 7.261249 min 0.000000 -1.047198e+00 1.570796 25% 1.072134 -7.494689e-01 7.848279 50% 2.144268 -3.535251e-16 14.125761 75% 3.216402 7.494689e-01 20.403244 max 4.288536 1.047198e+00 26.680726 In\u00a0[6]: Copied!
px.line(\n df_pen,\n x=\"t\",\n y=\"x\",\n title=rf\"Simple Harmonic Oscillator ($\\omega_0 = {omega0:.4f})$\",\n labels={\"x\": r\"Angle $\\theta(t)$\", \"t\": r\"Time $t$\"},\n)\n px.line( df_pen, x=\"t\", y=\"x\", title=rf\"Simple Harmonic Oscillator ($\\omega_0 = {omega0:.4f})$\", labels={\"x\": r\"Angle $\\theta(t)$\", \"t\": r\"Time $t$\"}, )"},{"location":"tutorials/pendulum/#pendulum","title":"Pendulum\u00b6","text":"In this tutorial, we demonstrate how to generate data of a pendulum, and introduce the mathematics of a pendulum.
"},{"location":"tutorials/pendulum/#constants","title":"Constants\u00b6","text":""},{"location":"tutorials/pendulum/#a-pendulum","title":"A pendulum\u00b6","text":""},{"location":"tutorials/pendulum/#data","title":"Data\u00b6","text":""},{"location":"tutorials/pendulum/#plot","title":"Plot\u00b6","text":""},{"location":"tutorials/pendulum/#todo","title":"TODO\u00b6","text":"We describe a generic pendulum system by the Lagrangian action $$ S_L[\\theta] \\equiv \\int_{t_0}^{t_1} \\mathbb{d}t\\,L(\\theta, \\dot\\theta) \\eqqcolon I \\int_{t_0}^{t_1} \\mathbb{d}t \\left\\{\\frac{1}{2} \\dot\\theta^2 + \\omega_0^2 \\cos\\theta \\right\\}\\,, $$ where $L$ is the Lagrangian; $\\theta$ is the angle from the vertical to the pendulum as the generalised position; $I$ is the inertia parameter, $\\omega_0$ the frequency parameter, and we also call $U \\coloneqq I\\omega_0^2$ the potential parameter.
This setup contains both the single and the physical pendula. For a single pendulum, $$ I = m l^2\\,,\\qquad U = mgl\\,, $$ where $m$ is the mass of the pendulum, $l$ is the length of the rod or cord, and $g$ is the gravitational acceleration.
"},{"location":"tutorials/pendulum/#integral-of-motion","title":"Integral of motion\u00b6","text":"The Lagrangian action does not contain time $t$ explicitly. As a result, the system is invariant under a variation of time, or $\\mathbb{\\delta}S / \\mathbb{\\delta}{t} = 0$. This gives an integral of motion $$ \\dot\\theta\\frac{\\partial L}{\\partial \\dot\\theta} - L \\equiv E \\eqqcolon I \\omega_0^2 \\cos\\theta_0\\,, $$ where $\\theta_0$ is the initial angle.
Substitution gives $$ \\left(\\frac{\\mathbb{d}t}{\\mathbb{d}\\theta}\\right)^2 = \\frac{1}{2\\omega_0^2} \\frac{1}{\\cos\\theta - \\cos\\theta_0}\\,. $$
"},{"location":"tutorials/pendulum/#coordinate-transformation","title":"Coordinate transformation\u00b6","text":"For convenience, introduce the coordinate $u$ and the parameter $k$ $$ \\sin u \\coloneqq \\frac{\\sin\\frac{\\theta}{2}}{k}\\,,\\qquad k \\coloneqq \\sin\\frac{\\theta_0}{2} \\in [-1, 1]\\,. $$ One arrives at $$ \\left(\\frac{\\mathbb{d}t}{\\mathbb{d}u}\\right)^2 = \\frac{1}{\\omega_0^2} \\frac{1}{1-k^2\\sin^2 u}\\,. $$ The square root of the second factor on the right-hand side makes an elliptic integral.
"},{"location":"tutorials/pendulum/#end-of-notebook","title":"End of Notebook\u00b6","text":""}]} \ No newline at end of file diff --git a/sitemap.xml b/sitemap.xml index 208fbf7..a0880b5 100644 --- a/sitemap.xml +++ b/sitemap.xml @@ -2,92 +2,92 @@We also plot out the ellipse, hyperbolic, parabolic orbits.
+visualize_orbit(
+ Kepler2D.from_geometry(
+ system=system,
+ geometries={
+ "positive_angular_mom": True,
+ "ecc": 0.5,
+ "parameter": 1.0,
+ },
+ )(t=t),
+).show()
+/home/runner/work/hamilflow/hamilflow/hamilflow/models/kepler_problem/dynamics.py:26: RuntimeWarning: + +divide by zero encountered in divide -# t = np.linspace(0, 1, 11)
visualize_orbit(
+ Kepler2D.from_geometry(
+ system=system,
+ geometries={
+ "positive_angular_mom": True,
+ "ecc": 1,
+ "parameter": 1.0,
+ },
+ )(t=t),
+).show()
+visualize_orbit(
+ Kepler2D.from_geometry(
+ system=system,
+ geometries={
+ "positive_angular_mom": True,
+ "ecc": 2,
+ "parameter": 1.0,
+ },
+ )(t=np.arange(-5, 5, 0.125)),
+).show()
+/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1046: RuntimeWarning: + +overflow encountered in scalar add + +/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1046: RuntimeWarning: + +invalid value encountered in scalar multiply + +/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:1054: RuntimeWarning: + +overflow encountered in scalar add + +/home/runner/.cache/pypoetry/virtualenvs/hamilflow-yj3wPJwu-py3.12/lib/python3.12/site-packages/scipy/optimize/_zeros_py.py:995: RuntimeWarning: + +overflow encountered in divide + ++