|
1 | 1 | import datetime |
| 2 | +from dataclasses import dataclass, field |
2 | 3 |
|
3 | 4 | import numpy as np |
4 | 5 | import matplotlib.pyplot as plt |
5 | 6 | import cartopy.crs as ccrs |
6 | 7 |
|
7 | 8 |
|
8 | 9 | from cmocean import cm |
9 | | - |
10 | 10 | from cartopy.feature.nightshade import Nightshade |
11 | 11 | from joblib import Parallel, delayed |
| 12 | +from matplotlib.colors import Normalize |
| 13 | +from mpl_toolkits.mplot3d.art3d import Line3DCollection |
12 | 14 | from scipy.integrate import solve_ivp |
13 | 15 | from scipy.special import j0 |
14 | 16 | from tqdm import tqdm |
15 | 17 |
|
| 18 | + |
16 | 19 | from matplotloom import Loom |
17 | 20 |
|
18 | 21 | def test_sine_wave(test_output_dir): |
@@ -252,3 +255,71 @@ def plot_frame(day_of_year, loom, frame_number): |
252 | 255 |
|
253 | 256 | assert (test_output_dir / "test_night_time_shading.mp4").is_file() |
254 | 257 | assert (test_output_dir / "test_night_time_shading.mp4").stat().st_size > 0 |
| 258 | + |
| 259 | +def test_lorenz(test_output_dir): |
| 260 | + @dataclass |
| 261 | + class Lorenz: |
| 262 | + dt: float = 0.01 |
| 263 | + sigma: float = 10.0 |
| 264 | + rho: float = 28.0 |
| 265 | + beta: float = 8.0 / 3.0 |
| 266 | + x: float = 1.0 |
| 267 | + y: float = 1.0 |
| 268 | + z: float = 1.0 |
| 269 | + |
| 270 | + def step(self): |
| 271 | + dx = self.sigma * (self.y - self.x) |
| 272 | + dy = self.x * (self.rho - self.z) - self.y |
| 273 | + dz = self.x * self.y - self.beta * self.z |
| 274 | + self.x += dx * self.dt |
| 275 | + self.y += dy * self.dt |
| 276 | + self.z += dz * self.dt |
| 277 | + |
| 278 | + @property |
| 279 | + def position(self) -> tuple[float, float, float]: |
| 280 | + return self.x, self.y, self.z |
| 281 | + |
| 282 | + @dataclass |
| 283 | + class LorenzPlotter: |
| 284 | + steps_per_frame: int = 20 |
| 285 | + attractor = Lorenz() |
| 286 | + points: list[tuple[float, float, float]] = field(default_factory=list) |
| 287 | + |
| 288 | + def initialize(self, steps: int): |
| 289 | + self.points = [self.attractor.position] |
| 290 | + for _ in range(steps): |
| 291 | + self.attractor.step() |
| 292 | + self.points.append(self.attractor.position) |
| 293 | + |
| 294 | + @property |
| 295 | + def frames(self) -> list[int]: |
| 296 | + return list(range(1, len(self.points) // self.steps_per_frame)) |
| 297 | + |
| 298 | + def get_frame(self, i: int, loom: Loom): |
| 299 | + fig, ax = plt.subplots(figsize=(12, 8), subplot_kw={'projection': '3d'}) |
| 300 | + points = np.array(self.points[: i * self.steps_per_frame]) |
| 301 | + xs, ys, zs = points.T |
| 302 | + segments = np.array([points[:-1], points[1:]]).transpose(1, 0, 2) |
| 303 | + norm = Normalize(vmin=0, vmax=len(xs)) |
| 304 | + colors = plt.get_cmap('inferno')(norm(np.arange(len(xs) - 1))) |
| 305 | + lc = Line3DCollection(segments, colors=colors, linewidth=0.5) |
| 306 | + ax.add_collection3d(lc) |
| 307 | + ax.set_xlim(-30, 30) |
| 308 | + ax.set_ylim(-30, 30) |
| 309 | + ax.set_zlim(0, 50) |
| 310 | + ax.view_init( |
| 311 | + azim=(np.pi * 1.7 + 0.8 * np.sin(2.0 * np.pi * i * self.steps_per_frame / len(self.frames) / 10)) |
| 312 | + * 180.0 |
| 313 | + / np.pi |
| 314 | + ) |
| 315 | + ax.set_axis_off() |
| 316 | + ax.grid(visible=False) |
| 317 | + loom.save_frame(fig, i - 1) |
| 318 | + |
| 319 | + with Loom(test_output_dir / 'test_lorenz.mp4', fps=60, parallel=True) as loom: |
| 320 | + attractor = LorenzPlotter() |
| 321 | + attractor.initialize(1000) |
| 322 | + Parallel(n_jobs=-1)(delayed(attractor.get_frame)(i, loom) for i in attractor.frames[:5]) |
| 323 | + |
| 324 | + assert (test_output_dir / "test_lorenz.mp4").is_file() |
| 325 | + assert (test_output_dir / "test_lorenz.mp4").stat().st_size > 0 |
0 commit comments