Геометрическая Волновая Инженерия

Псевдогиперболоиды 3+ порядков на Питоне

—— coding: utf-8 ——

«»»
ФИНАЛЬНЫЙ СКРИПТ.

Построение открытых псевдогиперболоидов n-го порядка по рекурсивной схеме

B_2 = {d},
B_{k+1} = {R_k + r, R_k - r : r in B_k},
N_n = 2^(n-2).

В ЭТОЙ ВЕРСИИ:
1) строится не просто набор ветвей, а ИТОГОВЫЙ ОБЩИЙ ОБЪЁМ;
2) если отдельные псевдогиперболоидные тороиды пересекаются или вложены,
удаляется только их общая часть как дублирование объёма;
3) сами тороиды как порождающие объекты НЕ удаляются;
4) итоговая фигура есть объединение всех объёмов;
5) если тороиды не пересекаются, они показываются отдельно;
6) правило действует для всех порядков и для всех рядов.

Поддерживаются:

  • вертикальный тип;
  • горизонтальный тип;
  • все порядки от 2 до n;
  • 2D меридиональные сечения общего объёма;
  • 3D поверхности общего объёма вращения;
  • стек одинаковых экземпляров на общей оси вращения.
    «»»

from future import annotations

import math
from pathlib import Path
from dataclasses import dataclass
from typing import Dict, List, Sequence, Tuple

import numpy as np
import matplotlib.pyplot as plt

————————————————————

НАСТРОЙКИ ОТОБРАЖЕНИЯ

————————————————————

plt.rcParams.update({
«figure.dpi»: 140,
«savefig.dpi»: 300,
«font.size»: 11,
«axes.grid»: True,
«grid.alpha»: 0.22,
«font.family»: «DejaVu Sans»,
})

SURFACE_ALPHA = 0.62
MERGE_TOL = 1.0e-9
NPHI = 220

В 2D по умолчанию НЕ заливаем области, чтобы не появлялись

искусственные прямые замыкающие отрезки от fill_betweenx.

FILL_2D_AREAS = False

————————————————————

БАЗОВАЯ ГЕОМЕТРИЯ

————————————————————

def c_focus(a: float, b: float) -> float:
return math.sqrt(a * a + b * b)

def L_open(a: float, b: float, R: float) -> float:
«»»Продольный предел существования открытого профиля вертикального типа.»»»
return a * math.sqrt(1.0 + (R / b) ** 2)

def eta_hyperbola(abs_x: np.ndarray, a: float, b: float) -> np.ndarray:
abs_x = np.asarray(abs_x, dtype=float)
val = (abs_x / a) ** 2 — 1.0
val = np.maximum(val, 0.0)
return b * np.sqrt(val)

def rho_vertical(abs_s: np.ndarray, a: float, b: float, R: float) -> np.ndarray:
«»»
Открытый профиль вертикального типа:
rho(|s|) = R — b*sqrt((|s|/a)^2 — 1), a <= |s| <= L.
«»»
out = R — eta_hyperbola(abs_s, a, b)
return np.maximum(out, 0.0)

def x_horizontal(u: np.ndarray, a: float, b: float, R: float) -> np.ndarray:
«»»
Горизонтальный тип:
x_h(u) = a * sqrt(1 + ((R — |u|)/b)^2), |u| <= R.
«»»
u = np.asarray(u, dtype=float)
return a * np.sqrt(1.0 + ((R — np.abs(u)) / b) ** 2)

————————————————————

БАЗОВЫЕ СЕТКИ ДЛЯ ДВУХ ТИПОВ

————————————————————

def _ensure_zero_in_axis(arr: np.ndarray) -> np.ndarray:
arr = np.asarray(arr, dtype=float)
if np.any(np.isclose(arr, 0.0, atol=1e-14)):
return arr
arr2 = np.concatenate([arr, np.array([0.0])])
arr2.sort()
return arr2

def vertical_base_grid(a: float, b: float, R: float, npts: int = 900) -> Dict[str, np.ndarray]:
L = L_open(a, b, R)
s_left = np.linspace(-L, -a, npts)
s_right = np.linspace(a, L, npts)
d_left = rho_vertical(np.abs(s_left), a, b, R)
d_right = rho_vertical(np.abs(s_right), a, b, R)
return {
«axis_left»: s_left,
«axis_right»: s_right,
«d_left»: d_left,
«d_right»: d_right,
«L»: L,
}

def horizontal_base_grid(a: float, b: float, R: float, npts: int = 900) -> Dict[str, np.ndarray]:
u = np.linspace(-R, R, npts)
u = _ensure_zero_in_axis(u)
d = x_horizontal(u, a, b, R)
return {
«axis»: u,
«d»: d,
«L»: R,
}

————————————————————

ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ

————————————————————

def morphology_class(R: float, offsets: Sequence[float]) -> str:
if not offsets:
return «только базовый 2-й порядок»
R_star = max(offsets)
if R > R_star:
return «с внутренним пересечением»
if math.isclose(R, R_star, rel_tol=1e-12, abs_tol=1e-12):
return «симметричный»
return «кольцевой непересекающийся»

def validate_offsets(base_distance: np.ndarray, offsets: Sequence[float]) -> List[str]:
«»»
Не блокирующая проверка. Строить всё равно можно.
Сообщения носят информационный характер.
«»»
msgs = []
current_max = float(np.max(base_distance))
for i, offset in enumerate(offsets, start=1):
if offset < current_max:
msgs.append(
f»На шаге R{i}: R{i} = {offset:.6f} < {current_max:.6f}. «
f»Это означает пересечение / вложение компонент и необходимость брать объединение объёмов.»
)
elif math.isclose(offset, current_max, rel_tol=1e-12, abs_tol=1e-12):
msgs.append(
f»На шаге R{i}: R{i} = {offset:.6f} равно {current_max:.6f}. «
f»Это предельный случай касания компонент.»
)
current_max = offset + current_max
return msgs

def stack_shifts(axis_half_length: float, h: float, m: int) -> np.ndarray:
«»»
Возвращает сдвиги центров m одинаковых экземпляров вдоль общей оси вращения.

Полная длина одного экземпляра вдоль общей оси = 2 * axis_half_length.
Шаг между центрами соседних экземпляров:
    step = 2 * axis_half_length + h
"""
if int(m) != m or m < 1:
    raise ValueError("Параметр m должен быть целым числом >= 1")
m = int(m)
step = 2.0 * float(axis_half_length) + float(h)
return -np.arange(m, dtype=float) * step

def _ensure_dir(path: Path):
path.mkdir(parents=True, exist_ok=True)

def _save_or_show(fig, path: Path | None, show: bool):
if path is not None:
fig.savefig(path, bbox_inches=»tight»)
if show:
plt.show()
else:
plt.close(fig)

————————————————————

РЕКУРСИВНЫЕ РАДИАЛЬНЫЕ ИНТЕРВАЛЫ

————————————————————

def build_recursive_interval_arrays(base_distance: np.ndarray,
offsets: Sequence[float]) -> List[Tuple[np.ndarray, np.ndarray]]:
«»»
Строит интервалы радиального объёма для одного порядка.

Порядок 2:
    [0, d]

Следующий шаг для каждого интервала [lo, hi]:
    [R - hi, R - lo]  и  [R + lo, R + hi]

Это соответствует объединению объёмов порождающих тороидов,
без удаления самих тороидов как объектов построения.
"""
zero = np.zeros_like(base_distance, dtype=float)
intervals = [(zero.copy(), np.asarray(base_distance, dtype=float).copy())]
for Rk in offsets:
    nxt: List[Tuple[np.ndarray, np.ndarray]] = []
    for lo, hi in intervals:
        a = np.maximum(Rk - hi, 0.0)
        b = np.maximum(Rk - lo, 0.0)
        c = np.maximum(Rk + lo, 0.0)
        d = np.maximum(Rk + hi, 0.0)
        nxt.append((a, b))
        nxt.append((c, d))
    intervals = nxt
return intervals

@dataclass
class Piece:
axis: np.ndarray
lo: np.ndarray
hi: np.ndarray

————————————————————

ПОСТРОЕНИЕ ОБЩЕГО ОБЪЁМА ДЛЯ ДАННОГО ПОРЯДКА

————————————————————

def _interp_piece_to_global_axis(piece: Piece, global_axis: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
x = piece.axis
lo = piece.lo
hi = piece.hi
lo_i = np.interp(global_axis, x, lo)
hi_i = np.interp(global_axis, x, hi)
mask = (global_axis >= x.min() — 1e-12) & (global_axis <= x.max() + 1e-12)
lo_i[~mask] = np.nan
hi_i[~mask] = np.nan
return lo_i, hi_i

def _merge_scalar_intervals(intervals: List[Tuple[float, float]], tol: float = MERGE_TOL) -> List[Tuple[float, float]]:
if not intervals:
return []
data = [(float(min(a, b)), float(max(a, b))) for a, b in intervals if np.isfinite(a) and np.isfinite(b) and b > a + tol]
if not data:
return []
data.sort(key=lambda t: (t[0], t[1]))
merged = [list(data[0])]
for a, b in data[1:]:
if a <= merged[-1][1] + tol:
merged[-1][1] = max(merged[-1][1], b)
else:
merged.append([a, b])
return [(a, b) for a, b in merged]

def _build_merged_components_from_pieces(pieces: List[Piece], n_global: int | None = None) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray]]:
if not pieces:
raise ValueError(«Нет компонент для построения общего объёма»)

# Глобальная ось: объединение всех реальных узлов.
all_axis = np.concatenate([p.axis for p in pieces])
global_axis = np.unique(np.round(all_axis, 12))
global_axis.sort()

if n_global is not None and n_global > global_axis.size:
    dense = np.linspace(global_axis.min(), global_axis.max(), n_global)
    global_axis = np.unique(np.concatenate([global_axis, dense]))
    global_axis.sort()

# ВАЖНО: если между двумя рядами / компонентами есть реальный пустой зазор,
# нельзя позволять линии или поверхности соединять края этого зазора.
# Поэтому добавляем контрольные точки в середины всех пустых промежутков
# между проекциями компонент на общую ось. В этих точках интервалов нет,
# и график разрывается сам, без дорисовывания.
cover = sorted((float(np.min(p.axis)), float(np.max(p.axis))) for p in pieces)
merged_cover: List[List[float]] = []
for a0, b0 in cover:
    if not merged_cover or a0 > merged_cover[-1][1] + MERGE_TOL:
        merged_cover.append([a0, b0])
    else:
        merged_cover[-1][1] = max(merged_cover[-1][1], b0)
gap_points = []
for left, right in zip(merged_cover[:-1], merged_cover[1:]):
    if right[0] > left[1] + MERGE_TOL:
        gap_points.append(0.5 * (left[1] + right[0]))
if gap_points:
    global_axis = np.unique(np.concatenate([global_axis, np.array(gap_points, dtype=float)]))
    global_axis.sort()

interpolated: List[Tuple[np.ndarray, np.ndarray]] = [
    _interp_piece_to_global_axis(piece, global_axis) for piece in pieces
]

merged_per_sample: List[List[Tuple[float, float]]] = []
max_count = 0
for j in range(global_axis.size):
    scalars: List[Tuple[float, float]] = []
    for lo_i, hi_i in interpolated:
        lo = lo_i[j]
        hi = hi_i[j]
        if np.isfinite(lo) and np.isfinite(hi) and hi > lo + MERGE_TOL:
            scalars.append((max(0.0, lo), max(0.0, hi)))
    merged = _merge_scalar_intervals(scalars, tol=MERGE_TOL)
    merged_per_sample.append(merged)
    max_count = max(max_count, len(merged))

lo_components = [np.full(global_axis.shape, np.nan, dtype=float) for _ in range(max_count)]
hi_components = [np.full(global_axis.shape, np.nan, dtype=float) for _ in range(max_count)]

for j, merged in enumerate(merged_per_sample):
    for k, (a, b) in enumerate(merged):
        lo_components[k][j] = a
        hi_components[k][j] = b

return global_axis, lo_components, hi_components

def _contiguous_segments(mask: np.ndarray, axis: np.ndarray | None = None) -> List[Tuple[int, int]]:
mask = np.asarray(mask, dtype=bool)
if mask.size == 0:
return []

# Дополнительная защита: если две соседние активные точки разделены
# слишком большим шагом по оси, это означает реальный пустой промежуток
# между компонентами. Его нельзя соединять прямой линией.
if axis is not None and axis.size == mask.size and axis.size >= 3:
    dif = np.diff(axis)
    pos = dif[np.isfinite(dif) & (dif > 0)]
    if pos.size:
        med = float(np.median(pos))
        jump_thr = max(1e-12, 4.0 * med)
        jump_breaks = np.where(dif > jump_thr)[0]
        for j in jump_breaks:
            if 0 <= j < mask.size - 1:
                mask[j] = False
                mask[j + 1] = False

segments: List[Tuple[int, int]] = []
start = None
for i, ok in enumerate(mask):
    if ok and start is None:
        start = i
    if (not ok) and start is not None:
        if i - start >= 2:
            segments.append((start, i))
        start = None
if start is not None and mask.size - start >= 2:
    segments.append((start, mask.size))
return segments

def _split_on_artificial_jumps(radius: np.ndarray) -> List[Tuple[int, int]]:
«»»
Дополнительно разбивает кривую там, где последовательность точек делает
негеометрический скачок из-за смены индекса объединённого интервала.
Это убирает искусственные прямые соединяющие отрезки.
«»»
radius = np.asarray(radius, dtype=float)
n = radius.size
if n < 2: return [(0, n)] if n > 0 else []
dr = np.abs(np.diff(radius))
finite = np.isfinite(dr)
valid = dr[finite]
if valid.size == 0:
return [(0, n)]
med = float(np.median(valid))
mx = float(np.max(valid))
if med <= 1e-12: jump_thr = max(1e-8, 0.05 * mx) else: jump_thr = max(1e-8, 8.0 * med) break_ids = np.where(dr > jump_thr)[0] + 1
if break_ids.size == 0:
return [(0, n)]
segments = []
st = 0
for br in break_ids:
if br — st >= 2:
segments.append((st, int(br)))
st = int(br)
if n — st >= 2:
segments.append((st, n))
return segments

def _plot_boundary_curve(ax, radius: np.ndarray, axis: np.ndarray, lw: float):
for st, en in _split_on_artificial_jumps(radius):
ax.plot(radius[st:en], axis[st:en], lw=lw)

————————————————————

ПОРОЖДАЮЩИЕ КОМПОНЕНТЫ ДЛЯ ОДНОГО ПОРЯДКА

————————————————————

def vertical_order_union_components(a: float,
b: float,
R: float,
offsets: Sequence[float],
npts: int,
m: int,
h: float) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray], float]:
base = vertical_base_grid(a, b, R, npts=npts)
shifts = stack_shifts(float(base[«L»]), h, m)

interval_left = build_recursive_interval_arrays(base["d_left"], offsets)
interval_right = build_recursive_interval_arrays(base["d_right"], offsets)

pieces: List[Piece] = []
for shift in shifts:
    for lo, hi in interval_left:
        pieces.append(Piece(base["axis_left"] + shift, lo, hi))
    for lo, hi in interval_right:
        pieces.append(Piece(base["axis_right"] + shift, lo, hi))

global_axis, lo_components, hi_components = _build_merged_components_from_pieces(pieces, n_global=max(2500, 3 * npts))
return global_axis, lo_components, hi_components, float(base["L"])

def horizontal_order_union_components(a: float,
b: float,
R: float,
offsets: Sequence[float],
npts: int,
m: int,
h: float) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray], float, np.ndarray]:
base = horizontal_base_grid(a, b, R, npts=npts)
shifts = stack_shifts(float(base[«L»]), h, m)

interval_base = build_recursive_interval_arrays(base["d"], offsets)

pieces: List[Piece] = []
for shift in shifts:
    for lo, hi in interval_base:
        pieces.append(Piece(base["axis"] + shift, lo, hi))

global_axis, lo_components, hi_components = _build_merged_components_from_pieces(pieces, n_global=max(2500, 3 * npts))
return global_axis, lo_components, hi_components, float(base["L"]), shifts

————————————————————

ПОВЕРХНОСТИ ВРАЩЕНИЯ

————————————————————

def revolve_about_x(axis_coord: np.ndarray, radius: np.ndarray, nphi: int = NPHI):
phi = np.linspace(0.0, 2.0 * np.pi, nphi)
S, Phi = np.meshgrid(axis_coord, phi, indexing=»ij»)
RR = np.tile(radius[:, None], (1, phi.size))
X = S
Y = RR * np.cos(Phi)
Z = RR * np.sin(Phi)
return X, Y, Z

def revolve_about_y(axis_coord: np.ndarray, radius: np.ndarray, nphi: int = NPHI):
phi = np.linspace(0.0, 2.0 * np.pi, nphi)
U, Phi = np.meshgrid(axis_coord, phi, indexing=»ij»)
RR = np.tile(radius[:, None], (1, phi.size))
X = RR * np.cos(Phi)
Y = U
Z = RR * np.sin(Phi)
return X, Y, Z

@dataclass
class PlotBounds3D:
xmin: float = math.inf
xmax: float = -math.inf
ymin: float = math.inf
ymax: float = -math.inf
zmin: float = math.inf
zmax: float = -math.inf

def update(self, X: np.ndarray, Y: np.ndarray, Z: np.ndarray):
    self.xmin = min(self.xmin, float(np.nanmin(X)))
    self.xmax = max(self.xmax, float(np.nanmax(X)))
    self.ymin = min(self.ymin, float(np.nanmin(Y)))
    self.ymax = max(self.ymax, float(np.nanmax(Y)))
    self.zmin = min(self.zmin, float(np.nanmin(Z)))
    self.zmax = max(self.zmax, float(np.nanmax(Z)))

def set_axes_equal_real_3d(ax, bounds: PlotBounds3D):
xmid = 0.5 * (bounds.xmin + bounds.xmax)
ymid = 0.5 * (bounds.ymin + bounds.ymax)
zmid = 0.5 * (bounds.zmin + bounds.zmax)
half = 0.5 * max(bounds.xmax — bounds.xmin, bounds.ymax — bounds.ymin, bounds.zmax — bounds.zmin)
if half <= 0:
half = 1.0
ax.set_xlim(xmid — half, xmid + half)
ax.set_ylim(ymid — half, ymid + half)
ax.set_zlim(zmid — half, zmid + half)
try:
ax.set_box_aspect((1.0, 1.0, 1.0))
except Exception:
pass
try:
ax.set_proj_type(«ortho»)
except Exception:
pass

def _plot_surface_piece_x(ax, axis_seg: np.ndarray, radius_seg: np.ndarray, bounds: PlotBounds3D, stride: int = 3):
X, Y, Z = revolve_about_x(axis_seg, radius_seg)
ax.plot_surface(
X[::stride, ::stride], Y[::stride, ::stride], Z[::stride, ::stride],
linewidth=0, edgecolor=»none», alpha=SURFACE_ALPHA, antialiased=False, shade=True
)
bounds.update(X, Y, Z)

def _plot_surface_piece_y(ax, axis_seg: np.ndarray, radius_seg: np.ndarray, bounds: PlotBounds3D,
split_positions: Sequence[float] = (), stride: int = 3):
# Для горизонтального типа разбиваем по положениям сдвигов рядов,
# чтобы в центральных зонах не сглаживались изломы.
split_ids = set()
for pos in split_positions:
idx = int(np.argmin(np.abs(axis_seg — pos)))
if 1 <= idx < axis_seg.size — 1 and np.isclose(axis_seg[idx], pos, atol=1e-10):
split_ids.add(idx)
starts = [0] + sorted(split_ids)
ends = sorted(split_ids) + [axis_seg.size — 1]
for st, en in zip(starts, ends):
seg_axis = axis_seg[st:en + 1]
seg_rad = radius_seg[st:en + 1]
if seg_axis.size < 2:
continue
X, Y, Z = revolve_about_y(seg_axis, seg_rad)
ax.plot_surface(
X[::stride, ::stride], Y[::stride, ::stride], Z[::stride, ::stride],
linewidth=0, edgecolor=»none», alpha=SURFACE_ALPHA, antialiased=False, shade=True
)
bounds.update(X, Y, Z)

————————————————————

2D РИСУНКИ ОБЩЕГО ОБЪЁМА

————————————————————

def _plot_component_2d(ax, axis: np.ndarray, lo: np.ndarray, hi: np.ndarray):
«»»
Рисует только реальные границы объединённого объёма.

ВАЖНО:
- скрипт ничего не дорисовывает;
- fill_betweenx / fill_between не используются;
- если после объединения интервалов возникает скачок нумерации компонент,
  кривая разрывается и не соединяется искусственным прямым отрезком.
"""
mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
for st, en in _contiguous_segments(mask, axis):
    yy = axis[st:en]
    lo_seg = lo[st:en]
    hi_seg = hi[st:en]

    _plot_boundary_curve(ax, hi_seg, yy, lw=1.8)
    _plot_boundary_curve(ax, -hi_seg, yy, lw=1.8)

    if np.nanmax(lo_seg) > MERGE_TOL:
        _plot_boundary_curve(ax, lo_seg, yy, lw=1.5)
        _plot_boundary_curve(ax, -lo_seg, yy, lw=1.5)

def plot_union_section_vertical(order: int,
global_axis: np.ndarray,
lo_components: List[np.ndarray],
hi_components: List[np.ndarray],
m: int,
h: float,
outpath: Path | None = None,
show: bool = True):
fig, ax = plt.subplots(figsize=(9, 10))
for lo, hi in zip(lo_components, hi_components):
_plot_component_2d(ax, global_axis, lo, hi)

rmax = 0.0
for hi in hi_components:
    if np.any(np.isfinite(hi)):
        rmax = max(rmax, float(np.nanmax(hi)))
if rmax <= 0:
    rmax = 1.0

ax.set_xlim(-1.05 * rmax, 1.05 * rmax)
ax.set_title(f"Вертикальный тип, порядок {order}: общий объём, 2D сечение, m={m}, h={h}")
ax.set_xlabel("Радиальная координата")
ax.set_ylabel("Координата вдоль общей оси")
ax.set_aspect("equal", adjustable="box")
fig.tight_layout()
_save_or_show(fig, outpath, show)

def plot_union_section_horizontal(order: int,
global_axis: np.ndarray,
lo_components: List[np.ndarray],
hi_components: List[np.ndarray],
m: int,
h: float,
outpath: Path | None = None,
show: bool = True):
fig, ax = plt.subplots(figsize=(10, 8))
for lo, hi in zip(lo_components, hi_components):
_plot_component_2d(ax, global_axis, lo, hi)

rmax = 0.0
for hi in hi_components:
    if np.any(np.isfinite(hi)):
        rmax = max(rmax, float(np.nanmax(hi)))
if rmax <= 0:
    rmax = 1.0

ax.set_xlim(-1.05 * rmax, 1.05 * rmax)
ax.set_title(f"Горизонтальный тип, порядок {order}: общий объём, 2D сечение, m={m}, h={h}")
ax.set_xlabel("Радиальная координата")
ax.set_ylabel("Координата вдоль общей оси")
ax.set_aspect("equal", adjustable="box")
fig.tight_layout()
_save_or_show(fig, outpath, show)

————————————————————

3D ПОВЕРХНОСТИ ОБЩЕГО ОБЪЁМА

————————————————————

def plot_union_surface_vertical(order: int,
global_axis: np.ndarray,
lo_components: List[np.ndarray],
hi_components: List[np.ndarray],
m: int,
h: float,
outpath: Path | None = None,
show: bool = True,
stride: int = 3):
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=»3d»)
bounds = PlotBounds3D()

for lo, hi in zip(lo_components, hi_components):
    mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
    for st, en in _contiguous_segments(mask, global_axis):
        axis_seg = global_axis[st:en]
        lo_seg = lo[st:en]
        hi_seg = hi[st:en]
        _plot_surface_piece_x(ax, axis_seg, hi_seg, bounds, stride=stride)
        if np.nanmax(lo_seg) > MERGE_TOL:
            _plot_surface_piece_x(ax, axis_seg, lo_seg, bounds, stride=stride)

set_axes_equal_real_3d(ax, bounds)
ax.set_title(f"Вертикальный тип, порядок {order}: 3D поверхность общего объёма, m={m}, h={h}")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
fig.tight_layout()
_save_or_show(fig, outpath, show)

def plot_union_surface_horizontal(order: int,
global_axis: np.ndarray,
lo_components: List[np.ndarray],
hi_components: List[np.ndarray],
shifts: np.ndarray,
m: int,
h: float,
outpath: Path | None = None,
show: bool = True,
stride: int = 3):
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=»3d»)
bounds = PlotBounds3D()

for lo, hi in zip(lo_components, hi_components):
    mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
    for st, en in _contiguous_segments(mask, global_axis):
        axis_seg = global_axis[st:en]
        lo_seg = lo[st:en]
        hi_seg = hi[st:en]
        _plot_surface_piece_y(ax, axis_seg, hi_seg, bounds, split_positions=shifts, stride=stride)
        if np.nanmax(lo_seg) > MERGE_TOL:
            _plot_surface_piece_y(ax, axis_seg, lo_seg, bounds, split_positions=shifts, stride=stride)

set_axes_equal_real_3d(ax, bounds)
ax.set_title(f"Горизонтальный тип, порядок {order}: 3D поверхность общего объёма, m={m}, h={h}")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
fig.tight_layout()
_save_or_show(fig, outpath, show)

————————————————————

ОСНОВНАЯ ФУНКЦИЯ

————————————————————

def run_pseudohyperboloids(a: float = 0.5,
b: float = 1.0,
R: float = 15.0,
offsets: Sequence[float] = (14.0, 30.0),
geometry_type: str = «both», # vertical / horizontal / both
mode: str = «all», # section / surface / all
outdir: str = «pseudohyperboloids_notebook_output»,
npts: int = 900,
h: float = -2.0,
m: int = 3,
show: bool = True) -> Dict[str, object]:
«»»
Главная функция.

Строится ИТОГОВЫЙ ОБЩИЙ ОБЪЁМ для каждого порядка от 2 до n,
включая пересечения и вложения компонент.
"""
offsets = list(offsets)
outdir_path = Path(outdir)
_ensure_dir(outdir_path)

if geometry_type not in {"vertical", "horizontal", "both"}:
    raise ValueError("geometry_type должен быть 'vertical', 'horizontal' или 'both'")
if mode not in {"section", "surface", "all"}:
    raise ValueError("mode должен быть 'section', 'surface' или 'all'")
if a <= 0 or b <= 0 or R <= 0:
    raise ValueError("Параметры a, b, R должны быть положительными")
if int(m) != m or m < 1:
    raise ValueError("Параметр m должен быть целым числом >= 1")

result = {
    "a": a,
    "b": b,
    "R": R,
    "offsets": offsets,
    "n_order": len(offsets) + 2,
    "focus_c": c_focus(a, b),
    "morphology": morphology_class(R, offsets),
    "m": int(m),
    "h": float(h),
    "volume_mode": "union_of_all_toroidal_components",
    "saved_files": [],
    "warnings": [],
}

print(f"a = {a}")
print(f"b = {b}")
print(f"R = {R}")
print(f"offsets = {offsets}")
print(f"Порядок n = {len(offsets) + 2}")
print(f"Фокусное расстояние c = {c_focus(a, b):.6f}")
print(f"Морфологический класс = {result['morphology']}")
print(f"m = {m}")
print(f"h = {h}")
print("Режим построения = общий объём как объединение всех компонент")
print("Правило: удаляется только общая перекрывающаяся часть объёма,")
print("а сами псевдогиперболоидные тороиды как порождающие объекты не удаляются.")
print(f"Папка вывода = {outdir_path.resolve()}")

if geometry_type in {"vertical", "both"}:
    base_v = vertical_base_grid(a, b, R, npts=npts)
    result["warnings"].extend([f"vertical-left: {w}" for w in validate_offsets(base_v["d_left"], offsets)])
    result["warnings"].extend([f"vertical-right: {w}" for w in validate_offsets(base_v["d_right"], offsets)])

    for p in range(2, len(offsets) + 3):
        used_offsets = offsets[:max(0, p - 2)]
        global_axis, lo_components, hi_components, _ = vertical_order_union_components(
            a, b, R, used_offsets, npts=npts, m=m, h=h
        )
        if mode in {"section", "all"}:
            path = outdir_path / f"vertical_order_{p:02d}_section_m{m}_h{h}.png"
            plot_union_section_vertical(
                order=p,
                global_axis=global_axis,
                lo_components=lo_components,
                hi_components=hi_components,
                m=m,
                h=h,
                outpath=path,
                show=show,
            )
            result["saved_files"].append(str(path))
        if mode in {"surface", "all"}:
            path = outdir_path / f"vertical_order_{p:02d}_surface_m{m}_h{h}.png"
            plot_union_surface_vertical(
                order=p,
                global_axis=global_axis,
                lo_components=lo_components,
                hi_components=hi_components,
                m=m,
                h=h,
                outpath=path,
                show=show,
            )
            result["saved_files"].append(str(path))

if geometry_type in {"horizontal", "both"}:
    base_h = horizontal_base_grid(a, b, R, npts=npts)
    result["warnings"].extend([f"horizontal: {w}" for w in validate_offsets(base_h["d"], offsets)])

    for p in range(2, len(offsets) + 3):
        used_offsets = offsets[:max(0, p - 2)]
        global_axis, lo_components, hi_components, _, shifts = horizontal_order_union_components(
            a, b, R, used_offsets, npts=npts, m=m, h=h
        )
        if mode in {"section", "all"}:
            path = outdir_path / f"horizontal_order_{p:02d}_section_m{m}_h{h}.png"
            plot_union_section_horizontal(
                order=p,
                global_axis=global_axis,
                lo_components=lo_components,
                hi_components=hi_components,
                m=m,
                h=h,
                outpath=path,
                show=show,
            )
            result["saved_files"].append(str(path))
        if mode in {"surface", "all"}:
            path = outdir_path / f"horizontal_order_{p:02d}_surface_m{m}_h{h}.png"
            plot_union_surface_horizontal(
                order=p,
                global_axis=global_axis,
                lo_components=lo_components,
                hi_components=hi_components,
                shifts=shifts,
                m=m,
                h=h,
                outpath=path,
                show=show,
            )
            result["saved_files"].append(str(path))

if result["warnings"]:
    print("\nПРЕДУПРЕЖДЕНИЯ:")
    for w in result["warnings"]:
        print(" -", w)

print("\nСОХРАНЁННЫЕ ФАЙЛЫ:")
for f in result["saved_files"]:
    print(" -", f)

return result

————————————————————

ПАРАМЕТРЫ ПО УМОЛЧАНИЮ ДЛЯ ЭТОЙ ТЕМЫ

————————————————————

if name == «main«:
a = 0.5
b = 1.0
R = 15.0
offsets = [14.0, 30.0] # R1, R2, …, R_{n-2}
geometry_type = «both» # «vertical», «horizontal», «both»
mode = «all» # «section», «surface», «all»
outdir = «pseudohyperboloids_notebook_output»
npts = 900
m = 3
h = -2.0
show = True

results = run_pseudohyperboloids(
    a=a,
    b=b,
    R=R,
    offsets=offsets,
    geometry_type=geometry_type,
    mode=mode,
    outdir=outdir,
    npts=npts,
    m=m,
    h=h,
    show=show,
)