Приложение А. Обозначения, параметры, базовые формулы и уровни описания
Настоящее приложение фиксирует единый язык обозначений для Томa 1. Его задача — устранить неоднозначность между геометрической образующей, радиальными интервалами, общим внутренним объёмом, рядной системой и будущими волновыми постановками. Все обозначения согласованы с финальным скриптом построения рядных псевдопараболоидов n-го порядка.
А.1. Базовые параметры псевдопараболоидной геометрии
Таблица А.1. Основные параметры базовой параболической геометрии.
| Обозначение | Смысл | Использование в теории |
| f | Фокусное расстояние параболических ветвей. | Задаёт кривизну исходной параболической образующей; должно быть положительным. |
| R | Расстояние от оси вращения / симметрии до фокальной оси образующей. | Главный масштаб псевдопараболоидной ячейки; должно быть положительным. |
| a = R²/(4f) | Характерная полудлина вертикального типа и предельная величина горизонтального типа. | Вычисляется функцией parabola_a(f, R); задаёт границы базового профиля. |
| dᵥ(s) | Базовая радиальная дистанция вертикального типа. | Используется при построении вертикального псевдопараболоида второго порядка. |
| dₕ(u) | Базовая радиальная дистанция горизонтального типа. | Используется при построении горизонтального псевдопараболоида второго порядка. |
| ξ | Обобщённая осевая координата. | Обозначает s для вертикального типа и u для горизонтального типа. |
| ρ | Радиальная координата меридионального сечения. | Используется в неравенствах внутреннего объёма. |
| φ | Угловая координата вращения. | Используется при восстановлении 3D-поверхностей вращения. |
В базовой редакции второго порядка вертикальный и горизонтальный типы не являются переименованием одной и той же координаты. Они имеют разные осевые диапазоны, разные формулы радиальной дистанции и разные поверхности вращения.
dᵥ(s) = R — 2√(f|s|), |s| ≤ a, a = R²/(4f).
dₕ(u) = (R — |u|)²/(4f), |u| ≤ R.
А.2. Порядок, рекурсивные смещения и формальное ветвевое ядро
Таблица А.2. Параметры рекурсии высших порядков.
| Обозначение | Смысл | Комментарий |
| n | Порядок псевдопараболоидной конструкции. | В финальном скрипте n = len(offsets) + 2. |
| Смещение | Список рекурсивных смещений [R₁, R₂, …, Rₙ₋₂]. | Каждый элемент повышает порядок на один уровень. |
| Rₖ | k-е рекурсивное радиальное смещение. | Действует на интервалы предыдущего порядка, а не заново на исходную образующую. |
| Bₙ | Формальное множество ветвей n-го порядка. | Используется как аналитический скелет, но не как окончательная физическая область. |
| Nₙ = 2ⁿ⁻² | Число формальных ветвей. | Не равно числу фактических Merge-компонент. |
| Iₙ(ξ) | Множество допустимых радиальных интервалов при фиксированной осевой координате. | Основной вычислительный мост от ветвей к объёму. |
| Ωₙ | Общий внутренний объём одиночного псевдопараболоида n-го порядка. | Физически значимая расчётная область до рядной сборки. |
B₂ = {d}, Bₖ₊₁ = {Rₖ₋₁ + r, Rₖ₋₁ — r : r ∈ Bₖ}, k ≥ 2.
Nₙ = 2ⁿ⁻².
Ветвевое ядро показывает происхождение границ, но не определяет полностью внутреннюю область. Для расчётной геометрии используется интервальная рекурсия, после которой выполняется Merge.
А.3. Интервальный оператор и Merge
Таблица А.3. Интервальный уровень построения.
| Объект | Формула / правило | Смысл |
| Второй порядок | I₂(ξ) = {[0, d(ξ)]} | Базовая внутренняя радиальная область. |
| Разностная компонента | [max(Rₖ — hi, 0), max(Rₖ — lo, 0)] | Компонента, возникающая при радиальном отражении относительно смещения Rₖ. |
| Суммовая компонента | [Rₖ + lo, Rₖ + hi] | Компонента, возникающая при радиальном переносе наружу. |
| Обрезание снизу | ρ ≥ 0 | Отрицательный радиус невозможен; поэтому используется max(.,0). |
| Merge | слияние пересекающихся или соприкасающихся интервалов | Удаляет только дублирование общей части объёма, не дорисовывая соединения. |
| NaN-разрывы / gap points | разрыв графика в пустых осевых промежутках | Защита от искусственных прямых линий в 2D и 3D. |
C_R([lo, hi]) = [max(R — hi, 0), max(R — lo, 0)] ∪ [R + lo, R + hi].
Iₖ₊₁(ξ) = Merge(⋃ C_{Rₖ₋₁}([lo, hi])), [lo, hi] ∈ Iₖ(ξ).
А.4. Рядность и общий объём Ωₙ,ₘ
Таблица А.4. Параметры рядной системы.
| Обозначение | Смысл | Комментарий |
| m | Число одинаковых экземпляров в рядной системе. | Целое число m ≥ 1. |
| h | Осевая добавка к шагу между центрами соседних экземпляров. | h > 0 — зазор; h = 0 — касание; h < 0 — перекрытие. |
| L | Половина длины одного экземпляра вдоль общей оси. | Для вертикального типа L = a = R²/(4f); для горизонтального типа L = R. |
| step | Шаг между центрами соседних экземпляров. | step = 2L + h. |
| Δⱼ | Осевой сдвиг j-го экземпляра. | Δⱼ = -j(2L + h), j = 0, 1, …, m-1. |
| Ωₙ,ₘ | Общий внутренний объём рядной системы. | Объединение всех сдвинутых экземпляров Ωₙ с последующим Merge. |
step = 2L + h, Δⱼ = -j(2L + h), j = 0, 1, …, m — 1.
Ωₙ,ₘ = ⋃_{j=0}^{m-1}(Ωₙ + Δⱼ e_ξ).
А.5. Уровни описания и границы физической интерпретации
Таблица А.5. Разделение уровней теории.
| Уровень | Основной объект | Что уже установлено | Что не следует утверждать без расчёта |
| Геометрический | dᵥ, dₕ, Iₙ, Ωₙ, Ωₙ,ₘ | Формулы, интервалы, Merge, рядность, 2D/3D-границы. | Нельзя автоматически объявлять локализацию или усиление волн. |
| Вычислительный | финальный Python-скрипт | Воспроизводимые построения, диагностика смещения, контроль разрывов, равный масштаб 3D. | Скрипт не является FEM/FDTD-решателем. |
| Физический | поля, лучи, моды, спектры, Q, потоки | Постановка будущих задач возможна на области Ωₙ,ₘ. | Результаты должны быть получены отдельным моделированием или экспериментом. |
| Инженерный | макет / резонатор / волновод / акустическая полость | Может быть задана геометрическая форма и контрольные параметры. | Нельзя заранее гарантировать универсальную эффективность для всех волн и частот. |
Таким образом, приложения фиксируют не лозунг, а рабочий аппарат: формулы базовой параболической образующей, рекурсивную интервальную геометрию, Merge, рядность и вычислительную реализацию. Это является достаточной конструктивной базой для передачи псевдопараболоидных областей в последующие физические расчёты.
Приложение Б. Финальный вычислительный скрипт построения рядных псевдопараболоидов n-го порядка
В данном приложении приведён полный текст финального Python-скрипта, на основании которого строятся 2D-меридиональные сечения, 3D-поверхности вращения, порядки от 2 до n, вертикальный и горизонтальный типы, рядные системы и Merge общего объёма.
Код приводится как воспроизводимая вычислительная основа Томa 1. При переносе в рабочую среду необходимо сохранить кодировку UTF-8 и зависимости: numpy, matplotlib, pathlib, dataclasses.
Таблица Б.1. Паспорт финального скрипта.
| Параметр | Значение |
| Имя файла | скрипт_рядных_псевдопараболоидов_N_порядка_ФИНАЛ_2D_оси_R_3D_с_ярким_2D_и_сеткой.py |
| Количество строк | 1603 |
| Язык | Python 3 |
| Назначение | Построение рядных псевдопараболоидов n-го порядка с вертикальным и горизонтальным типами. |
| Ключевые функции | vertical_distance, horizontal_distance, build_recursive_interval_arrays, vertical_order_union_components, horizontal_order_union_components, run_pseudoparaboloids. |
| Графические режимы | profile, section, surface, all. |
Полный листинг скрипта
0001: # -*- coding: utf-8 -*-
0002: «»»
0003: ФИНАЛЬНЫЙ СКРИПТ.
0004:
0005: Рядные псевдопараболоиды n-го порядка по архитектуре скрипта
0006: рядных псевдогиперболоидов, но с другой базовой образующей.
0007:
0008: БАЗОВАЯ ГЕОМЕТРИЯ ПСЕВДОПАРАБОЛОИДА
0009: ————————————
0010: Вертикальный тип строится из двух одинаковых зеркальных параболических ветвей,
0011: соединённых в общей вершине на линии фокусов:
0012:
0013: y(x) = sqrt(4 f |x|) = 2 sqrt(f |x|), |x| <= a,
0014: a = R^2 / (4 f).
0015:
0016: Ось вращения проходит параллельно оси x на расстоянии R от линии фокусов.
0017: Радиальная дистанция от оси вращения до образующей:
0018:
0019: d_v(x) = R — 2 sqrt(f |x|), |x| <= a.
0020:
0021: Горизонтальный тип — та же параболическая образующая, повернутая на 90°:
0022:
0023: d_h(u) = (R — |u|)^2 / (4 f), |u| <= R.
0024:
0025: РЕКУРСИЯ ПОРЯДКА
0026: —————-
0027: Порядок 2:
0028:
0029: [0, d]
0030:
0031: Следующий порядок для каждого интервала [lo, hi]:
0032:
0033: [R_k — hi, R_k — lo] и [R_k + lo, R_k + hi],
0034:
0035: с последующим Merge радиальных интервалов. Это означает:
0036: — сами порождающие тороидальные компоненты не удаляются;
0037: — удаляется только дублирование объёма в местах пересечения / вложения;
0038: — итоговая фигура есть объединение всех объёмов;
0039: — если между компонентами есть пустой зазор, скрипт ничего не дорисовывает.
0040:
0041: ПОДДЕРЖИВАЕТСЯ
0042: —————
0043: — вертикальный тип;
0044: — горизонтальный тип;
0045: — все порядки от 2 до n;
0046: — 2D меридиональные сечения общего объёма;
0047: — 3D поверхности общего объёма вращения;
0048: — стек / ряд m одинаковых экземпляров на общей оси;
0049: — параметр h: h > 0 разносит ряды, h = 0 даёт касание, h < 0 даёт перекрытие;
0050: — защита от искусственных прямых соединений в 2D и 3D.
0051: «»»
0052:
0053: from __future__ import annotations
0054:
0055: import math
0056: from pathlib import Path
0057: from dataclasses import dataclass
0058: from typing import Dict, List, Sequence, Tuple
0059:
0060: import numpy as np
0061: import matplotlib.pyplot as plt
0062:
0063: # ————————————————————
0064: # НАСТРОЙКИ ОТОБРАЖЕНИЯ
0065: # ————————————————————
0066: plt.rcParams.update({
0067: «figure.dpi»: 140,
0068: «savefig.dpi»: 300,
0069: «font.size»: 11,
0070: «axes.grid»: True,
0071: «grid.alpha»: 0.22,
0072: «font.family»: «DejaVu Sans»,
0073: })
0074:
0075: SURFACE_ALPHA = 0.42
0076: MERGE_TOL = 1.0e-9
0077: NPHI = 220
0078:
0079: # ВАЖНО: заливки по умолчанию нет. fill_between / fill_betweenx не используются,
0080: # чтобы не появлялись искусственные замыкающие отрезки.
0081: FILL_2D_AREAS = False
0082:
0083:
0084: # ————————————————————
0085: # БАЗОВАЯ ПАРАБОЛИЧЕСКАЯ ГЕОМЕТРИЯ
0086: # ————————————————————
0087: def parabola_a(f: float, R: float) -> float:
0088: «»»Полудлина вертикального типа / экваториальный радиус горизонтального типа.»»»
0089: if f <= 0 or R <= 0:
0090: raise ValueError(«Параметры f и R должны быть положительными»)
0091: return (R * R) / (4.0 * f)
0092:
0093:
0094: def vertical_distance(abs_s: np.ndarray, f: float, R: float) -> np.ndarray:
0095: «»»
0096: Вертикальный псевдопараболоид:
0097: d_v(|s|) = R — 2*sqrt(f*|s|), 0 <= |s| <= a.
0098: «»»
0099: abs_s = np.asarray(abs_s, dtype=float)
0100: d = R — 2.0 * np.sqrt(np.maximum(0.0, f * abs_s))
0101: return np.maximum(d, 0.0)
0102:
0103:
0104: def horizontal_distance(u: np.ndarray, f: float, R: float) -> np.ndarray:
0105: «»»
0106: Горизонтальный псевдопараболоид:
0107: d_h(u) = (R — |u|)^2 / (4*f), |u| <= R.
0108: «»»
0109: u = np.asarray(u, dtype=float)
0110: d = ((R — np.abs(u)) ** 2) / (4.0 * f)
0111: return np.maximum(d, 0.0)
0112:
0113:
0114: def _ensure_zero_in_axis(arr: np.ndarray) -> np.ndarray:
0115: arr = np.asarray(arr, dtype=float)
0116: if np.any(np.isclose(arr, 0.0, atol=1e-14)):
0117: return arr
0118: arr2 = np.concatenate([arr, np.array([0.0])])
0119: arr2.sort()
0120: return arr2
0121:
0122:
0123: def vertical_base_grid(f: float, R: float, npts: int = 900) -> Dict[str, np.ndarray]:
0124: «»»
0125: Базовая сетка вертикального типа.
0126:
0127: Ось фигуры: s in [-a, a].
0128: Радиальный профиль: d_v(s).
0129: «»»
0130: a = parabola_a(f, R)
0131: s = np.linspace(-a, a, npts)
0132: s = _ensure_zero_in_axis(s)
0133: d = vertical_distance(np.abs(s), f, R)
0134: return {
0135: «axis»: s,
0136: «d»: d,
0137: «L»: float(a),
0138: «a»: float(a),
0139: «max_distance»: float(np.nanmax(d)),
0140: }
0141:
0142:
0143: def horizontal_base_grid(f: float, R: float, npts: int = 900) -> Dict[str, np.ndarray]:
0144: «»»
0145: Базовая сетка горизонтального типа.
0146:
0147: Ось фигуры: u in [-R, R].
0148: Радиальный профиль: d_h(u).
0149: «»»
0150: a = parabola_a(f, R)
0151: u = np.linspace(-R, R, npts)
0152: u = _ensure_zero_in_axis(u)
0153: d = horizontal_distance(u, f, R)
0154: return {
0155: «axis»: u,
0156: «d»: d,
0157: «L»: float(R),
0158: «a»: float(a),
0159: «max_distance»: float(np.nanmax(d)),
0160: }
0161:
0162:
0163: # ————————————————————
0164: # ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ
0165: # ————————————————————
0166: def morphology_class_from_base(base_max_distance: float, offsets: Sequence[float]) -> str:
0167: if not offsets:
0168: return «только базовый 2-й порядок»
0169: R_star = max(offsets)
0170: if base_max_distance > R_star:
0171: return «с внутренним пересечением / вложением»
0172: if math.isclose(base_max_distance, R_star, rel_tol=1e-12, abs_tol=1e-12):
0173: return «предельное касание»
0174: return «кольцевой / разнесённый режим до Merge»
0175:
0176:
0177: def validate_offsets(base_distance: np.ndarray, offsets: Sequence[float]) -> List[str]:
0178: «»»
0179: Не блокирующая проверка. Строить всё равно можно.
0180: Сообщения носят информационный характер.
0181: «»»
0182: msgs: List[str] = []
0183: current_max = float(np.max(base_distance))
0184: for i, offset in enumerate(offsets, start=1):
0185: if offset < current_max:
0186: msgs.append(
0187: f»На шаге R{i}: R{i} = {offset:.6f} < {current_max:.6f}. «
0188: f»Это режим пересечения / вложения компонент; нужен Merge общего объёма.»
0189: )
0190: elif math.isclose(offset, current_max, rel_tol=1e-12, abs_tol=1e-12):
0191: msgs.append(
0192: f»На шаге R{i}: R{i} = {offset:.6f} равно {current_max:.6f}. «
0193: f»Это предельное касание компонент.»
0194: )
0195: current_max = offset + current_max
0196: return msgs
0197:
0198:
0199: def stack_shifts(axis_half_length: float, h: float, m: int) -> np.ndarray:
0200: «»»
0201: Сдвиги центров m одинаковых экземпляров вдоль общей оси вращения.
0202:
0203: Полная длина одного экземпляра вдоль общей оси = 2 * axis_half_length.
0204: Шаг между центрами соседних экземпляров:
0205: step = 2 * axis_half_length + h
0206:
0207: h > 0 — зазор;
0208: h = 0 — касание;
0209: h < 0 — перекрытие / Merge.
0210: «»»
0211: if int(m) != m or m < 1:
0212: raise ValueError(«Параметр m должен быть целым числом >= 1»)
0213: m = int(m)
0214: step = 2.0 * float(axis_half_length) + float(h)
0215: return -np.arange(m, dtype=float) * step
0216:
0217:
0218: def _ensure_dir(path: Path):
0219: path.mkdir(parents=True, exist_ok=True)
0220:
0221:
0222: def _save_or_show(fig, path: Path | None, show: bool):
0223: if path is not None:
0224: fig.savefig(path, bbox_inches=»tight»)
0225: if show:
0226: plt.show()
0227: else:
0228: plt.close(fig)
0229:
0230:
0231: # ————————————————————
0232: # РЕКУРСИВНЫЕ РАДИАЛЬНЫЕ ИНТЕРВАЛЫ
0233: # ————————————————————
0234: def build_recursive_interval_arrays(base_distance: np.ndarray,
0235: offsets: Sequence[float]) -> List[Tuple[np.ndarray, np.ndarray]]:
0236: «»»
0237: Радиальные интервалы для одного порядка.
0238:
0239: Порядок 2:
0240: [0, d]
0241:
0242: Следующий шаг для каждого интервала [lo, hi]:
0243: [R_k — hi, R_k — lo] и [R_k + lo, R_k + hi].
0244:
0245: Затем интервалы объединяются в общем Merge. Это не удаляет порождающие
0246: компоненты как объекты построения, а убирает только дублирование объёма.
0247: «»»
0248: zero = np.zeros_like(base_distance, dtype=float)
0249: intervals = [(zero.copy(), np.asarray(base_distance, dtype=float).copy())]
0250: for Rk in offsets:
0251: nxt: List[Tuple[np.ndarray, np.ndarray]] = []
0252: for lo, hi in intervals:
0253: a = np.maximum(Rk — hi, 0.0)
0254: b = np.maximum(Rk — lo, 0.0)
0255: c = np.maximum(Rk + lo, 0.0)
0256: d = np.maximum(Rk + hi, 0.0)
0257: nxt.append((a, b))
0258: nxt.append((c, d))
0259: intervals = nxt
0260: return intervals
0261:
0262:
0263: @dataclass
0264: class Piece:
0265: axis: np.ndarray
0266: lo: np.ndarray
0267: hi: np.ndarray
0268:
0269:
0270: # ————————————————————
0271: # MERGE ОБЩЕГО ОБЪЁМА
0272: # ————————————————————
0273: def _interp_piece_to_global_axis(piece: Piece, global_axis: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
0274: x = piece.axis
0275: lo = piece.lo
0276: hi = piece.hi
0277: lo_i = np.interp(global_axis, x, lo)
0278: hi_i = np.interp(global_axis, x, hi)
0279: mask = (global_axis >= x.min() — 1e-12) & (global_axis <= x.max() + 1e-12)
0280: lo_i[~mask] = np.nan
0281: hi_i[~mask] = np.nan
0282: return lo_i, hi_i
0283:
0284:
0285: def _merge_scalar_intervals(intervals: List[Tuple[float, float]], tol: float = MERGE_TOL) -> List[Tuple[float, float]]:
0286: if not intervals:
0287: return []
0288: data = [
0289: (float(min(a, b)), float(max(a, b)))
0290: for a, b in intervals
0291: if np.isfinite(a) and np.isfinite(b) and b > a + tol
0292: ]
0293: if not data:
0294: return []
0295: data.sort(key=lambda t: (t[0], t[1]))
0296: merged = [list(data[0])]
0297: for a, b in data[1:]:
0298: if a <= merged[-1][1] + tol:
0299: merged[-1][1] = max(merged[-1][1], b)
0300: else:
0301: merged.append([a, b])
0302: return [(a, b) for a, b in merged]
0303:
0304:
0305: def _build_merged_components_from_pieces(pieces: List[Piece],
0306: n_global: int | None = None) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray]]:
0307: if not pieces:
0308: raise ValueError(«Нет компонент для построения общего объёма»)
0309:
0310: all_axis = np.concatenate([p.axis for p in pieces])
0311: global_axis = np.unique(np.round(all_axis, 12))
0312: global_axis.sort()
0313:
0314: if n_global is not None and n_global > global_axis.size:
0315: dense = np.linspace(global_axis.min(), global_axis.max(), n_global)
0316: global_axis = np.unique(np.concatenate([global_axis, dense]))
0317: global_axis.sort()
0318:
0319: # Если между рядами / компонентами есть реальный пустой зазор, нельзя
0320: # соединять края зазора линиями или поверхностями. Добавляем контрольную
0321: # точку в середину каждого пустого промежутка; там интервалов нет, и график
0322: # разорвётся сам.
0323: cover = sorted((float(np.min(p.axis)), float(np.max(p.axis))) for p in pieces)
0324: merged_cover: List[List[float]] = []
0325: for a0, b0 in cover:
0326: if not merged_cover or a0 > merged_cover[-1][1] + MERGE_TOL:
0327: merged_cover.append([a0, b0])
0328: else:
0329: merged_cover[-1][1] = max(merged_cover[-1][1], b0)
0330:
0331: gap_points: List[float] = []
0332: for left, right in zip(merged_cover[:-1], merged_cover[1:]):
0333: if right[0] > left[1] + MERGE_TOL:
0334: gap_points.append(0.5 * (left[1] + right[0]))
0335: if gap_points:
0336: global_axis = np.unique(np.concatenate([global_axis, np.array(gap_points, dtype=float)]))
0337: global_axis.sort()
0338:
0339: interpolated = [_interp_piece_to_global_axis(piece, global_axis) for piece in pieces]
0340:
0341: merged_per_sample: List[List[Tuple[float, float]]] = []
0342: max_count = 0
0343: for j in range(global_axis.size):
0344: scalars: List[Tuple[float, float]] = []
0345: for lo_i, hi_i in interpolated:
0346: lo = lo_i[j]
0347: hi = hi_i[j]
0348: if np.isfinite(lo) and np.isfinite(hi) and hi > lo + MERGE_TOL:
0349: scalars.append((max(0.0, lo), max(0.0, hi)))
0350: merged = _merge_scalar_intervals(scalars, tol=MERGE_TOL)
0351: merged_per_sample.append(merged)
0352: max_count = max(max_count, len(merged))
0353:
0354: lo_components = [np.full(global_axis.shape, np.nan, dtype=float) for _ in range(max_count)]
0355: hi_components = [np.full(global_axis.shape, np.nan, dtype=float) for _ in range(max_count)]
0356:
0357: for j, merged in enumerate(merged_per_sample):
0358: for k, (a0, b0) in enumerate(merged):
0359: lo_components[k][j] = a0
0360: hi_components[k][j] = b0
0361:
0362: return global_axis, lo_components, hi_components
0363:
0364:
0365: def _contiguous_segments(mask: np.ndarray, axis: np.ndarray | None = None) -> List[Tuple[int, int]]:
0366: mask = np.asarray(mask, dtype=bool).copy()
0367: if mask.size == 0:
0368: return []
0369:
0370: # Защита от соединения краёв реального пустого промежутка.
0371: if axis is not None and axis.size == mask.size and axis.size >= 3:
0372: dif = np.diff(axis)
0373: pos = dif[np.isfinite(dif) & (dif > 0)]
0374: if pos.size:
0375: med = float(np.median(pos))
0376: jump_thr = max(1e-12, 4.0 * med)
0377: jump_breaks = np.where(dif > jump_thr)[0]
0378: for j in jump_breaks:
0379: if 0 <= j < mask.size — 1:
0380: mask[j] = False
0381: mask[j + 1] = False
0382:
0383: segments: List[Tuple[int, int]] = []
0384: start = None
0385: for i, ok in enumerate(mask):
0386: if ok and start is None:
0387: start = i
0388: if (not ok) and start is not None:
0389: if i — start >= 2:
0390: segments.append((start, i))
0391: start = None
0392: if start is not None and mask.size — start >= 2:
0393: segments.append((start, mask.size))
0394: return segments
0395:
0396:
0397: def _split_on_artificial_jumps(radius: np.ndarray) -> List[Tuple[int, int]]:
0398: «»»
0399: Разрыв только при настоящем скачке радиуса, возникающем из-за смены индекса
0400: объединённого интервала. Параболические вершины не считаются ошибкой.
0401: «»»
0402: radius = np.asarray(radius, dtype=float)
0403: n = radius.size
0404: if n < 2:
0405: return [(0, n)] if n > 0 else []
0406: dr = np.abs(np.diff(radius))
0407: valid = dr[np.isfinite(dr)]
0408: if valid.size == 0:
0409: return [(0, n)]
0410:
0411: finite_r = radius[np.isfinite(radius)]
0412: if finite_r.size == 0:
0413: return [(0, n)]
0414: r_range = float(np.nanmax(finite_r) — np.nanmin(finite_r))
0415: med = float(np.median(valid))
0416:
0417: # Жёсткий порог: не режем естественную параболическую вершину из-за большой
0418: # производной, но режем скачки от переиндексации Merge-компонент.
0419: jump_thr = max(1e-8, 12.0 * med, 0.08 * r_range)
0420: break_ids = np.where(dr > jump_thr)[0] + 1
0421: if break_ids.size == 0:
0422: return [(0, n)]
0423:
0424: segments: List[Tuple[int, int]] = []
0425: st = 0
0426: for br in break_ids:
0427: if br — st >= 2:
0428: segments.append((st, int(br)))
0429: st = int(br)
0430: if n — st >= 2:
0431: segments.append((st, n))
0432: return segments
0433:
0434:
0435: def _plot_boundary_curve(ax, radius: np.ndarray, axis: np.ndarray, lw: float):
0436: for st, en in _split_on_artificial_jumps(radius):
0437: ax.plot(radius[st:en], axis[st:en], lw=lw)
0438:
0439:
0440: # ————————————————————
0441: # ПОРОЖДАЮЩИЕ КОМПОНЕНТЫ ДЛЯ ОДНОГО ПОРЯДКА
0442: # ————————————————————
0443: def vertical_order_union_components(f: float,
0444: R: float,
0445: offsets: Sequence[float],
0446: npts: int,
0447: m: int,
0448: h: float) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray], float, np.ndarray]:
0449: base = vertical_base_grid(f, R, npts=npts)
0450: shifts = stack_shifts(float(base[«L»]), h, m)
0451: interval_base = build_recursive_interval_arrays(base[«d»], offsets)
0452:
0453: pieces: List[Piece] = []
0454: for shift in shifts:
0455: for lo, hi in interval_base:
0456: pieces.append(Piece(base[«axis»] + shift, lo, hi))
0457:
0458: global_axis, lo_components, hi_components = _build_merged_components_from_pieces(
0459: pieces,
0460: n_global=max(5000, 6 * npts, 1200 * (len(offsets) + 1)),
0461: )
0462: return global_axis, lo_components, hi_components, float(base[«L»]), shifts
0463:
0464:
0465: def horizontal_order_union_components(f: float,
0466: R: float,
0467: offsets: Sequence[float],
0468: npts: int,
0469: m: int,
0470: h: float) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray], float, np.ndarray]:
0471: base = horizontal_base_grid(f, R, npts=npts)
0472: shifts = stack_shifts(float(base[«L»]), h, m)
0473: interval_base = build_recursive_interval_arrays(base[«d»], offsets)
0474:
0475: pieces: List[Piece] = []
0476: for shift in shifts:
0477: for lo, hi in interval_base:
0478: pieces.append(Piece(base[«axis»] + shift, lo, hi))
0479:
0480: global_axis, lo_components, hi_components = _build_merged_components_from_pieces(
0481: pieces,
0482: n_global=max(5000, 6 * npts, 1200 * (len(offsets) + 1)),
0483: )
0484: return global_axis, lo_components, hi_components, float(base[«L»]), shifts
0485:
0486:
0487: # ————————————————————
0488: # ПОВЕРХНОСТИ ВРАЩЕНИЯ
0489: # ————————————————————
0490: def revolve_about_x(axis_coord: np.ndarray, radius: np.ndarray, nphi: int = NPHI):
0491: phi = np.linspace(0.0, 2.0 * np.pi, nphi)
0492: S, Phi = np.meshgrid(axis_coord, phi, indexing=»ij»)
0493: RR = np.tile(radius[:, None], (1, phi.size))
0494: X = S
0495: Y = RR * np.cos(Phi)
0496: Z = RR * np.sin(Phi)
0497: return X, Y, Z
0498:
0499:
0500: def revolve_about_y(axis_coord: np.ndarray, radius: np.ndarray, nphi: int = NPHI):
0501: phi = np.linspace(0.0, 2.0 * np.pi, nphi)
0502: U, Phi = np.meshgrid(axis_coord, phi, indexing=»ij»)
0503: RR = np.tile(radius[:, None], (1, phi.size))
0504: X = RR * np.cos(Phi)
0505: Y = U
0506: Z = RR * np.sin(Phi)
0507: return X, Y, Z
0508:
0509:
0510: @dataclass
0511: class PlotBounds3D:
0512: xmin: float = math.inf
0513: xmax: float = -math.inf
0514: ymin: float = math.inf
0515: ymax: float = -math.inf
0516: zmin: float = math.inf
0517: zmax: float = -math.inf
0518:
0519: def update(self, X: np.ndarray, Y: np.ndarray, Z: np.ndarray):
0520: self.xmin = min(self.xmin, float(np.nanmin(X)))
0521: self.xmax = max(self.xmax, float(np.nanmax(X)))
0522: self.ymin = min(self.ymin, float(np.nanmin(Y)))
0523: self.ymax = max(self.ymax, float(np.nanmax(Y)))
0524: self.zmin = min(self.zmin, float(np.nanmin(Z)))
0525: self.zmax = max(self.zmax, float(np.nanmax(Z)))
0526:
0527: def is_valid(self) -> bool:
0528: return all(math.isfinite(v) for v in (self.xmin, self.xmax, self.ymin, self.ymax, self.zmin, self.zmax))
0529:
0530:
0531: def set_axes_equal_real_3d(ax, bounds: PlotBounds3D):
0532: if not bounds.is_valid():
0533: return
0534: xmid = 0.5 * (bounds.xmin + bounds.xmax)
0535: ymid = 0.5 * (bounds.ymin + bounds.ymax)
0536: zmid = 0.5 * (bounds.zmin + bounds.zmax)
0537: half = 0.5 * max(bounds.xmax — bounds.xmin, bounds.ymax — bounds.ymin, bounds.zmax — bounds.zmin)
0538: if half <= 0:
0539: half = 1.0
0540: ax.set_xlim(xmid — half, xmid + half)
0541: ax.set_ylim(ymid — half, ymid + half)
0542: ax.set_zlim(zmid — half, zmid + half)
0543: try:
0544: ax.set_box_aspect((1.0, 1.0, 1.0))
0545: except Exception:
0546: pass
0547: try:
0548: ax.set_proj_type(«ortho»)
0549: except Exception:
0550: pass
0551:
0552:
0553: def _plot_surface_piece_x(ax,
0554: axis_seg: np.ndarray,
0555: radius_seg: np.ndarray,
0556: bounds: PlotBounds3D,
0557: stride: int = 2):
0558: X, Y, Z = revolve_about_x(axis_seg, radius_seg)
0559: ax.plot_surface(
0560: X[::stride, ::stride], Y[::stride, ::stride], Z[::stride, ::stride],
0561: linewidth=0, edgecolor=»none», alpha=SURFACE_ALPHA, antialiased=False, shade=True,
0562: )
0563: bounds.update(X, Y, Z)
0564:
0565:
0566: def _plot_surface_piece_y(ax,
0567: axis_seg: np.ndarray,
0568: radius_seg: np.ndarray,
0569: bounds: PlotBounds3D,
0570: split_positions: Sequence[float] = (),
0571: stride: int = 2):
0572: # Разбиваем в положениях центров рядов, чтобы не сглаживать реальные вершины.
0573: split_ids = set()
0574: for pos in split_positions:
0575: idx = int(np.argmin(np.abs(axis_seg — pos)))
0576: if 1 <= idx < axis_seg.size — 1 and np.isclose(axis_seg[idx], pos, atol=1e-10):
0577: split_ids.add(idx)
0578: starts = [0] + sorted(split_ids)
0579: ends = sorted(split_ids) + [axis_seg.size — 1]
0580:
0581: for st, en in zip(starts, ends):
0582: seg_axis = axis_seg[st:en + 1]
0583: seg_rad = radius_seg[st:en + 1]
0584: if seg_axis.size < 2:
0585: continue
0586: X, Y, Z = revolve_about_y(seg_axis, seg_rad)
0587: ax.plot_surface(
0588: X[::stride, ::stride], Y[::stride, ::stride], Z[::stride, ::stride],
0589: linewidth=0, edgecolor=»none», alpha=SURFACE_ALPHA, antialiased=False, shade=True,
0590: )
0591: bounds.update(X, Y, Z)
0592:
0593:
0594: def _surface_subsegments(axis_seg: np.ndarray,
0595: radius_seg: np.ndarray,
0596: split_positions: Sequence[float] = ()) -> List[Tuple[np.ndarray, np.ndarray]]:
0597: «»»
0598: Разбивает граничную кривую поверхности на реальные подучастки, чтобы
0599: не сглаживать и не склеивать внутренние параболические образующие при Merge.
0600:
0601: Учитываются:
0602: — искусственные скачки радиуса из-за переиндексации Merge-компонент;
0603: — реальные точки смены ряда / вершины по split_positions.
0604: «»»
0605: axis_seg = np.asarray(axis_seg, dtype=float)
0606: radius_seg = np.asarray(radius_seg, dtype=float)
0607: if axis_seg.size < 2:
0608: return []
0609:
0610: cut_ids = {0, axis_seg.size}
0611:
0612: for st, en in _split_on_artificial_jumps(radius_seg):
0613: cut_ids.add(st)
0614: cut_ids.add(en)
0615:
0616: for pos in split_positions:
0617: idx0 = int(np.argmin(np.abs(axis_seg — pos)))
0618: if 1 <= idx0 < axis_seg.size — 1:
0619: cut_ids.add(idx0)
0620: cut_ids.add(idx0 + 1)
0621:
0622: cuts = sorted(i for i in cut_ids if 0 <= i <= axis_seg.size)
0623: out: List[Tuple[np.ndarray, np.ndarray]] = []
0624: for a0, b0 in zip(cuts[:-1], cuts[1:]):
0625: if b0 — a0 >= 2:
0626: aa = axis_seg[a0:b0]
0627: rr = radius_seg[a0:b0]
0628: if aa.size >= 2 and np.all(np.isfinite(rr)):
0629: out.append((aa, rr))
0630: return out
0631:
0632:
0633: def _plot_surface_boundary_x(ax,
0634: axis_seg: np.ndarray,
0635: radius_seg: np.ndarray,
0636: bounds: PlotBounds3D,
0637: stride: int = 2):
0638: for aa, rr in _surface_subsegments(axis_seg, radius_seg, split_positions=()):
0639: _plot_surface_piece_x(ax, aa, rr, bounds, stride=stride)
0640:
0641:
0642: def _plot_surface_boundary_y(ax,
0643: axis_seg: np.ndarray,
0644: radius_seg: np.ndarray,
0645: bounds: PlotBounds3D,
0646: split_positions: Sequence[float] = (),
0647: stride: int = 2):
0648: for aa, rr in _surface_subsegments(axis_seg, radius_seg, split_positions=split_positions):
0649: _plot_surface_piece_y(ax, aa, rr, bounds, split_positions=split_positions, stride=stride)
0650:
0651:
0652: # ————————————————————
0653: # БОКОВАЯ 2D-ВИЗУАЛИЗАЦИЯ НА 3D
0654: # ————————————————————
0655: def _plot_side_section_on_vertical_3d(ax,
0656: global_axis: np.ndarray,
0657: lo_components: List[np.ndarray],
0658: hi_components: List[np.ndarray],
0659: bounds: PlotBounds3D):
0660: «»»
0661: Показывает яркое меридиональное 2D-сечение вертикального типа В ПЛОСКОСТИ ОСИ
0662: ВРАЩЕНИЯ 3D: Y = 0. Внешние границы выделяются одним ярким цветом,
0663: внутренние — другим, чтобы сразу видеть полный объединённый объём.
0664: «»»
0665: y_plane = 0.0
0666: outer_color = ‘magenta’
0667: inner_color = ‘lime’
0668:
0669: for lo, hi in zip(lo_components, hi_components):
0670: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
0671: for st, en in _contiguous_segments(mask, global_axis):
0672: xx = global_axis[st:en]
0673: lo_seg = lo[st:en]
0674: hi_seg = hi[st:en]
0675:
0676: for a0, b0 in _split_on_artificial_jumps(hi_seg):
0677: x = xx[a0:b0]
0678: z = hi_seg[a0:b0]
0679: if x.size >= 2:
0680: y = np.full_like(x, y_plane)
0681: ax.plot(x, y, z, color=outer_color, linewidth=3.0, alpha=1.0)
0682: ax.plot(x, y, -z, color=outer_color, linewidth=3.0, alpha=1.0)
0683: if np.nanmax(lo_seg) > MERGE_TOL:
0684: for a0, b0 in _split_on_artificial_jumps(lo_seg):
0685: x = xx[a0:b0]
0686: z = lo_seg[a0:b0]
0687: if x.size >= 2:
0688: y = np.full_like(x, y_plane)
0689: ax.plot(x, y, z, color=inner_color, linewidth=2.6, alpha=1.0)
0690: ax.plot(x, y, -z, color=inner_color, linewidth=2.6, alpha=1.0)
0691:
0692:
0693: def _plot_side_section_on_horizontal_3d(ax,
0694: global_axis: np.ndarray,
0695: lo_components: List[np.ndarray],
0696: hi_components: List[np.ndarray],
0697: bounds: PlotBounds3D):
0698: «»»
0699: Показывает яркое меридиональное 2D-сечение горизонтального типа В ПЛОСКОСТИ ОСИ
0700: ВРАЩЕНИЯ 3D: X = 0. Внешние границы выделяются одним ярким цветом,
0701: внутренние — другим, чтобы сразу видеть полный объединённый объём.
0702: «»»
0703: x_plane = 0.0
0704: outer_color = ‘magenta’
0705: inner_color = ‘lime’
0706:
0707: for lo, hi in zip(lo_components, hi_components):
0708: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
0709: for st, en in _contiguous_segments(mask, global_axis):
0710: yy = global_axis[st:en]
0711: lo_seg = lo[st:en]
0712: hi_seg = hi[st:en]
0713:
0714: for a0, b0 in _split_on_artificial_jumps(hi_seg):
0715: y = yy[a0:b0]
0716: z = hi_seg[a0:b0]
0717: if y.size >= 2:
0718: x = np.full_like(y, x_plane)
0719: ax.plot(x, y, z, color=outer_color, linewidth=3.0, alpha=1.0)
0720: ax.plot(x, y, -z, color=outer_color, linewidth=3.0, alpha=1.0)
0721: if np.nanmax(lo_seg) > MERGE_TOL:
0722: for a0, b0 in _split_on_artificial_jumps(lo_seg):
0723: y = yy[a0:b0]
0724: z = lo_seg[a0:b0]
0725: if y.size >= 2:
0726: x = np.full_like(y, x_plane)
0727: ax.plot(x, y, z, color=inner_color, linewidth=2.6, alpha=1.0)
0728: ax.plot(x, y, -z, color=inner_color, linewidth=2.6, alpha=1.0)
0729:
0730:
0731: # ————————————————————
0732: # 2D СЕЧЕНИЯ ОБЩЕГО ОБЪЁМА
0733: # ————————————————————
0734: def _plot_component_2d(ax, axis: np.ndarray, lo: np.ndarray, hi: np.ndarray):
0735: «»»
0736: Рисует только реальные границы объединённого объёма.
0737:
0738: Никаких служебных линий, никаких искусственных замыканий, никаких
0739: fill_between. Если компонент разорван, линия разрывается.
0740: «»»
0741: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
0742: for st, en in _contiguous_segments(mask, axis):
0743: yy = axis[st:en]
0744: lo_seg = lo[st:en]
0745: hi_seg = hi[st:en]
0746:
0747: _plot_boundary_curve(ax, hi_seg, yy, lw=1.8)
0748: _plot_boundary_curve(ax, -hi_seg, yy, lw=1.8)
0749:
0750: if np.nanmax(lo_seg) > MERGE_TOL:
0751: _plot_boundary_curve(ax, lo_seg, yy, lw=1.5)
0752: _plot_boundary_curve(ax, -lo_seg, yy, lw=1.5)
0753:
0754:
0755: def _axis_limits_from_components(axis: np.ndarray, hi_components: List[np.ndarray]) -> Tuple[float, float, float]:
0756: rmax = 0.0
0757: for hi in hi_components:
0758: if np.any(np.isfinite(hi)):
0759: rmax = max(rmax, float(np.nanmax(hi)))
0760: if rmax <= 0:
0761: rmax = 1.0
0762: amin = float(np.nanmin(axis))
0763: amax = float(np.nanmax(axis))
0764: return rmax, amin, amax
0765:
0766:
0767: def build_recursive_boundary_descriptors(offsets: Sequence[float]):
0768: «»»
0769: Строит точные дескрипторы границ радиальных интервалов без Merge.
0770:
0771: Дескриптор границы:
0772: («const», t) -> r(s) = t
0773: («parab», c, sigma) -> r(s) = c + sigma * d_base(s)
0774:
0775: где d_base(s) — базовая параболическая образующая текущего типа.
0776:
0777: ВАЖНО:
0778: — фокусы вычисляются ТОЛЬКО из этих фактических параболических кривых;
0779: — никакой отдельной «рекурсии фокусов» нет;
0780: — любая новая параболическая кривая получается тем же преобразованием,
0781: что и сама граница, и её фокус переносится синхронно с ней.
0782: «»»
0783: intervals = [((«const», 0.0), («parab», 0.0, +1.0))]
0784: for Rk in offsets:
0785: nxt = []
0786: for lo, hi in intervals:
0787: nxt.append((_desc_sub_const(Rk, hi), _desc_sub_const(Rk, lo)))
0788: nxt.append((_desc_add_const(Rk, lo), _desc_add_const(Rk, hi)))
0789: intervals = nxt
0790: return intervals
0791:
0792:
0793: def _desc_add_const(C: float, desc):
0794: if desc[0] == «const»:
0795: return («const», float(C + desc[1]))
0796: _, c, sigma = desc
0797: return («parab», float(C + c), float(sigma))
0798:
0799:
0800: def _desc_sub_const(C: float, desc):
0801: if desc[0] == «const»:
0802: return («const», float(C — desc[1]))
0803: _, c, sigma = desc
0804: return («parab», float(C — c), float(-sigma))
0805:
0806:
0807: def _split_curve_at_axis_positions(x: np.ndarray,
0808: y: np.ndarray,
0809: split_positions: Sequence[float],
0810: tol: float = 1e-10) -> List[Tuple[np.ndarray, np.ndarray]]:
0811: «»»
0812: Разбивает уже построенную 2D-кривую по осевым положениям, где базовая
0813: формула с |.| меняет ветвь. Это не дорисовывание, а только разбиение
0814: уже существующей линии на параболические участки.
0815: «»»
0816: x = np.asarray(x, dtype=float)
0817: y = np.asarray(y, dtype=float)
0818: if x.size < 2:
0819: return []
0820:
0821: cut_ids = {0, x.size}
0822: for pos in split_positions:
0823: idx = int(np.argmin(np.abs(y — pos)))
0824: if 1 <= idx < x.size — 1 and abs(float(y[idx] — pos)) <= max(tol, 1e-7 * max(1.0, np.nanmax(np.abs(y)))):
0825: cut_ids.add(idx)
0826: cut_ids.add(idx + 1)
0827:
0828: cuts = sorted(cut_ids)
0829: out: List[Tuple[np.ndarray, np.ndarray]] = []
0830: for a0, b0 in zip(cuts[:-1], cuts[1:]):
0831: if b0 — a0 >= 8:
0832: out.append((x[a0:b0], y[a0:b0]))
0833: return out
0834:
0835:
0836: def _quadratic_focus_vertical(x: np.ndarray, y: np.ndarray):
0837: «»»
0838: Фокус параболы, заданной фактически построенными точками y = A*x^2+B*x+C.
0839: «»»
0840: x = np.asarray(x, dtype=float)
0841: y = np.asarray(y, dtype=float)
0842: ok = np.isfinite(x) & np.isfinite(y)
0843: x = x[ok]
0844: y = y[ok]
0845: if x.size < 8 or (np.nanmax(x) — np.nanmin(x)) < 1e-9:
0846: return None
0847:
0848: A, B, C = np.polyfit(x, y, 2)
0849: if abs(A) < 1e-12:
0850: return None
0851:
0852: y_fit = A*x*x + B*x + C
0853: scale = max(1.0, float(np.nanmax(y) — np.nanmin(y)))
0854: rmse = float(np.sqrt(np.mean((y — y_fit)**2))) / scale
0855: if rmse > 2.5e-5:
0856: return None
0857:
0858: xv = -B / (2.0*A)
0859: yv = C — B*B / (4.0*A)
0860: yf = yv + 1.0 / (4.0*A)
0861: return float(xv), float(yf)
0862:
0863:
0864: def _quadratic_focus_horizontal(x: np.ndarray, y: np.ndarray):
0865: «»»
0866: Фокус параболы, заданной фактически построенными точками x = A*y^2+B*y+C.
0867: «»»
0868: x = np.asarray(x, dtype=float)
0869: y = np.asarray(y, dtype=float)
0870: ok = np.isfinite(x) & np.isfinite(y)
0871: x = x[ok]
0872: y = y[ok]
0873: if x.size < 8 or (np.nanmax(y) — np.nanmin(y)) < 1e-9:
0874: return None
0875:
0876: A, B, C = np.polyfit(y, x, 2)
0877: if abs(A) < 1e-12:
0878: return None
0879:
0880: x_fit = A*y*y + B*y + C
0881: scale = max(1.0, float(np.nanmax(x) — np.nanmin(x)))
0882: rmse = float(np.sqrt(np.mean((x — x_fit)**2))) / scale
0883: if rmse > 2.5e-5:
0884: return None
0885:
0886: yv = -B / (2.0*A)
0887: xv = C — B*B / (4.0*A)
0888: xf = xv + 1.0 / (4.0*A)
0889: return float(xf), float(yv)
0890:
0891:
0892: def _adaptive_focus_from_drawn_curve(x: np.ndarray,
0893: y: np.ndarray,
0894: orientation: str,
0895: min_points: int = 18,
0896: depth: int = 0,
0897: max_depth: int = 7) -> List[Tuple[float, float]]:
0898: «»»
0899: Вычисляет фокусы только по фактически построенной 2D-линии.
0900:
0901: Если линия не является одной параболой, она делится на участки.
0902: Это нужно для Merge-случаев, когда видимая граница может переходить
0903: с одной параболической кривой на другую.
0904: «»»
0905: x = np.asarray(x, dtype=float)
0906: y = np.asarray(y, dtype=float)
0907:
0908: ok = np.isfinite(x) & np.isfinite(y)
0909: x = x[ok]
0910: y = y[ok]
0911: if x.size < min_points:
0912: return []
0913:
0914: if orientation == «vertical»:
0915: focus = _quadratic_focus_vertical(x, y)
0916: else:
0917: focus = _quadratic_focus_horizontal(x, y)
0918:
0919: if focus is not None:
0920: return [focus]
0921:
0922: if depth >= max_depth or x.size < 2*min_points:
0923: return []
0924:
0925: mid = x.size // 2
0926: return (
0927: _adaptive_focus_from_drawn_curve(x[:mid], y[:mid], orientation, min_points, depth+1, max_depth)
0928: + _adaptive_focus_from_drawn_curve(x[mid:], y[mid:], orientation, min_points, depth+1, max_depth)
0929: )
0930:
0931:
0932: def _collect_drawn_boundary_curves(axis: np.ndarray,
0933: lo_components: List[np.ndarray],
0934: hi_components: List[np.ndarray]) -> List[Tuple[np.ndarray, np.ndarray]]:
0935: «»»
0936: Возвращает ровно те 2D-линии, которые рисует _plot_component_2d:
0937: внешние границы ±hi и внутренние границы ±lo.
0938: «»»
0939: curves: List[Tuple[np.ndarray, np.ndarray]] = []
0940: for lo, hi in zip(lo_components, hi_components):
0941: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
0942: for st, en in _contiguous_segments(mask, axis):
0943: yy = axis[st:en]
0944: lo_seg = lo[st:en]
0945: hi_seg = hi[st:en]
0946:
0947: for rad in (hi_seg, -hi_seg):
0948: for a0, b0 in _split_on_artificial_jumps(rad):
0949: if b0 — a0 >= 8:
0950: curves.append((np.asarray(rad[a0:b0], dtype=float),
0951: np.asarray(yy[a0:b0], dtype=float)))
0952:
0953: if np.nanmax(lo_seg) > MERGE_TOL:
0954: for rad in (lo_seg, -lo_seg):
0955: for a0, b0 in _split_on_artificial_jumps(rad):
0956: if b0 — a0 >= 8:
0957: curves.append((np.asarray(rad[a0:b0], dtype=float),
0958: np.asarray(yy[a0:b0], dtype=float)))
0959: return curves
0960:
0961:
0962: def _focus_points_from_drawn_2d(axis: np.ndarray,
0963: lo_components: List[np.ndarray],
0964: hi_components: List[np.ndarray],
0965: shifts: np.ndarray,
0966: orientation: str) -> List[Tuple[float, float]]:
0967: «»»
0968: Главная функция фокусов.
0969:
0970: 1. Берёт только фактически построенные на 2D рисунке линии.
0971: 2. Делит их на параболические участки.
0972: 3. Вычисляет фокус каждого участка по его собственному уравнению.
0973: «»»
0974: curves = _collect_drawn_boundary_curves(axis, lo_components, hi_components)
0975: points: List[Tuple[float, float]] = []
0976:
0977: for x, y in curves:
0978: # Разбиваем в осевых центрах рядов, где формула с abs меняет ветвь.
0979: pieces = _split_curve_at_axis_positions(x, y, split_positions=shifts)
0980: if not pieces:
0981: pieces = [(x, y)]
0982:
0983: for px, py in pieces:
0984: points.extend(_adaptive_focus_from_drawn_curve(px, py, orientation=orientation))
0985:
0986: # Удаляем дубли.
0987: unique = sorted({(round(float(x), 10), round(float(y), 10)) for x, y in points},
0988: key=lambda p: (p[1], p[0]))
0989: return [(float(x), float(y)) for x, y in unique]
0990:
0991:
0992: def _plot_focus_points_2d(ax,
0993: points: List[Tuple[float, float]],
0994: rmax: float):
0995: if not points:
0996: return
0997:
0998: xs = [p[0] for p in points]
0999: ys = [p[1] for p in points]
1000:
1001: ax.scatter(xs, ys, s=24, marker=’o’, color=’darkgreen’, zorder=7,
1002: label=’Фокусы фактически построенных 2D-парабол’)
1003:
1004: dx = 0.018 * max(rmax, 1.0)
1005: for x, y in points:
1006: if abs(x) <= 1e-12:
1007: ax.text(x + dx, y, ‘F’, fontsize=8, color=’darkgreen’,
1008: va=’center’, ha=’left’)
1009: else:
1010: ha = ‘left’ if x > 0 else ‘right’
1011: x_text = x + dx if x > 0 else x — dx
1012: ax.text(x_text, y, ‘F’, fontsize=8, color=’darkgreen’,
1013: va=’center’, ha=ha)
1014:
1015:
1016:
1017: # ————————————————————
1018: # АННОТАЦИИ ДЛЯ 2D РИСУНКОВ 2-ГО ПОРЯДКА
1019: # ————————————————————
1020: def _draw_common_rotation_axis_label(ax, axis_min: float, axis_max: float):
1021: «»»Центральная пунктирная ось вращения r=0 для 2D.»»»
1022: ax.axvline(0.0, color=’black’, linestyle=’—‘, linewidth=1.1, zorder=2)
1023: ax.text(0.03, 0.97, ‘Ось вращения’,
1024: transform=ax.transAxes, ha=’left’, va=’top’, fontsize=9)
1025:
1026:
1027: def _annotate_vertical_order2_2d(ax,
1028: f: float,
1029: R: float,
1030: a: float,
1031: focus_points: List[Tuple[float, float]],
1032: rmax: float,
1033: axis_min: float,
1034: axis_max: float):
1035: «»»
1036: Вертикальный псевдопараболоид 2-го порядка:
1037: — центральная пунктирная ось вращения;
1038: — пунктирная фокальная ось образующей;
1039: — размер R между этими двумя пунктирными линиями.
1040: «»»
1041: _draw_common_rotation_axis_label(ax, axis_min, axis_max)
1042:
1043: x_focus_axis = float(R)
1044:
1045: ax.axvline(x_focus_axis, color=’black’, linestyle=’—‘, linewidth=1.05, zorder=2)
1046: ax.text(x_focus_axis + 0.035 * max(rmax, 1.0),
1047: axis_min + 0.72 * (axis_max — axis_min),
1048: ‘Фокальная ось\nобразующей’,
1049: ha=’left’, va=’center’, fontsize=9)
1050:
1051: yR = axis_min + 0.18 * (axis_max — axis_min)
1052: ax.annotate(», xy=(x_focus_axis, yR), xytext=(0.0, yR),
1053: arrowprops=dict(arrowstyle='<->’, lw=1.15, color=’black’))
1054: ax.text(0.5 * x_focus_axis, yR + 0.035 * (axis_max — axis_min),
1055: ‘R’, ha=’center’, va=’bottom’, fontsize=12)
1056: ax.text(0.5 * x_focus_axis, yR — 0.025 * (axis_max — axis_min),
1057: ‘расстояние от оси вращения до фокальной оси’,
1058: ha=’center’, va=’top’, fontsize=8)
1059:
1060: near = [(x, y) for x, y in focus_points if abs(x — x_focus_axis) <= 0.08 * max(R, 1.0)]
1061: if near:
1062: near.sort(key=lambda p: p[1])
1063: y_low = near[0][1]
1064: y_high = near[-1][1]
1065: ax.text(x_focus_axis — 0.04 * max(rmax, 1.0), y_high,
1066: ‘+f’, ha=’right’, va=’center’, fontsize=9)
1067: ax.text(x_focus_axis — 0.04 * max(rmax, 1.0), y_low,
1068: ‘−f’, ha=’right’, va=’center’, fontsize=9)
1069:
1070:
1071: def _annotate_horizontal_order2_2d(ax,
1072: f: float,
1073: R: float,
1074: a: float,
1075: focus_points: List[Tuple[float, float]],
1076: rmax: float,
1077: axis_min: float,
1078: axis_max: float):
1079: «»»
1080: Горизонтальный псевдопараболоид 2-го порядка:
1081: — центральная горизонтальная пунктирная ось симметрии y=0;
1082: — горизонтальная фокальная ось образующей y=R;
1083: — размер R между y=0 и y=R.
1084: «»»
1085: # Центральная вертикальная линия остаётся только как геометрическая ось,
1086: # без подписи «вертикальная ось симметрии».
1087: ax.axvline(0.0, color=’black’, linestyle=’—‘, linewidth=1.1, zorder=2)
1088:
1089: ax.axhline(0.0, color=’black’, linestyle=’—‘, linewidth=1.05, zorder=2)
1090: ax.text(0.98, 0.52, ‘центральная горизонтальная\nось симметрии y = 0’,
1091: transform=ax.transAxes, ha=’right’, va=’bottom’, fontsize=8)
1092:
1093: ax.axhline(R, color=’black’, linestyle=’—‘, linewidth=1.05, zorder=2)
1094: ax.text(0.98, 0.88, ‘Фокальная ось\nобразующей y = R’,
1095: transform=ax.transAxes, ha=’right’, va=’top’, fontsize=9)
1096:
1097: xR = -0.82 * max(rmax, a, 1.0)
1098: ax.annotate(», xy=(xR, R), xytext=(xR, 0.0),
1099: arrowprops=dict(arrowstyle='<->’, lw=1.15, color=’black’))
1100: ax.text(xR — 0.03 * max(rmax, 1.0), 0.5 * R,
1101: ‘R’, ha=’right’, va=’center’, fontsize=12)
1102: ax.text(xR + 0.02 * max(rmax, 1.0), 0.5 * R,
1103: ‘расстояние от оси симметрии\nдо оси фокусов’,
1104: ha=’left’, va=’center’, fontsize=8)
1105:
1106: ax.annotate(», xy=(a, 0.0), xytext=(0.0, 0.0),
1107: arrowprops=dict(arrowstyle='<->’, lw=1.0, color=’black’))
1108: ax.text(0.5 * a, 0.05 * max(R, 1.0),
1109: ‘a = R²/(4f)’, ha=’center’, va=’bottom’, fontsize=9)
1110:
1111:
1112: def _hide_3d_service_axes(ax):
1113: «»»
1114: Возвращает на 3D координатную сетку, оси и tick marks, чтобы можно было
1115: визуально контролировать полный объём и внутренние образующие.
1116: «»»
1117: ax.set_axis_on()
1118: ax.grid(True)
1119: ax.set_xlabel(«X»)
1120: ax.set_ylabel(«Y»)
1121: ax.set_zlabel(«Z»)
1122:
1123:
1124: def plot_union_section_vertical(order: int,
1125: global_axis: np.ndarray,
1126: lo_components: List[np.ndarray],
1127: hi_components: List[np.ndarray],
1128: shifts: np.ndarray,
1129: used_offsets: Sequence[float],
1130: f: float,
1131: R: float,
1132: m: int,
1133: h: float,
1134: outpath: Path | None = None,
1135: show: bool = True):
1136: fig, ax = plt.subplots(figsize=(9, 10))
1137: for lo, hi in zip(lo_components, hi_components):
1138: _plot_component_2d(ax, global_axis, lo, hi)
1139:
1140: rmax, _, _ = _axis_limits_from_components(global_axis, hi_components)
1141: focus_points = _focus_points_from_drawn_2d(
1142: axis=global_axis,
1143: lo_components=lo_components,
1144: hi_components=hi_components,
1145: shifts=shifts,
1146: orientation=»vertical»,
1147: )
1148: _plot_focus_points_2d(ax, focus_points, rmax=rmax)
1149:
1150: if focus_points:
1151: rmax = max(rmax, max(abs(x) for x, _ in focus_points))
1152:
1153: _, axis_min, axis_max = _axis_limits_from_components(global_axis, hi_components)
1154: if order == 2:
1155: _annotate_vertical_order2_2d(
1156: ax=ax, f=f, R=R, a=parabola_a(f, R), focus_points=focus_points,
1157: rmax=rmax, axis_min=axis_min, axis_max=axis_max,
1158: )
1159:
1160: ax.set_xlim(-1.05 * rmax, 1.18 * rmax if order == 2 else 1.05 * rmax)
1161: ax.set_title(f»Псевдопараболоид: вертикальный тип, порядок {order}, 2D, f={f}, R={R}, m={m}, h={h}»)
1162: ax.set_xlabel(«Радиальная координата»)
1163: ax.set_ylabel(«Координата вдоль общей оси»)
1164: ax.set_aspect(«equal», adjustable=»box»)
1165: ax.legend(loc=’best’, fontsize=9)
1166: fig.tight_layout()
1167: _save_or_show(fig, outpath, show)
1168:
1169:
1170: def plot_union_section_horizontal(order: int,
1171: global_axis: np.ndarray,
1172: lo_components: List[np.ndarray],
1173: hi_components: List[np.ndarray],
1174: shifts: np.ndarray,
1175: used_offsets: Sequence[float],
1176: f: float,
1177: R: float,
1178: m: int,
1179: h: float,
1180: outpath: Path | None = None,
1181: show: bool = True):
1182: fig, ax = plt.subplots(figsize=(10, 8))
1183: for lo, hi in zip(lo_components, hi_components):
1184: _plot_component_2d(ax, global_axis, lo, hi)
1185:
1186: rmax, _, _ = _axis_limits_from_components(global_axis, hi_components)
1187: focus_points = _focus_points_from_drawn_2d(
1188: axis=global_axis,
1189: lo_components=lo_components,
1190: hi_components=hi_components,
1191: shifts=shifts,
1192: orientation=»horizontal»,
1193: )
1194: _plot_focus_points_2d(ax, focus_points, rmax=rmax)
1195:
1196: if focus_points:
1197: rmax = max(rmax, max(abs(x) for x, _ in focus_points))
1198:
1199: _, axis_min, axis_max = _axis_limits_from_components(global_axis, hi_components)
1200: if order == 2:
1201: _annotate_horizontal_order2_2d(
1202: ax=ax, f=f, R=R, a=parabola_a(f, R), focus_points=focus_points,
1203: rmax=rmax, axis_min=axis_min, axis_max=axis_max,
1204: )
1205:
1206: ax.set_xlim(-1.10 * max(rmax, parabola_a(f, R)), 1.10 * max(rmax, parabola_a(f, R)))
1207: ax.set_title(f»Псевдопараболоид: горизонтальный тип, порядок {order}, 2D, f={f}, R={R}, m={m}, h={h}»)
1208: ax.set_xlabel(«Радиальная координата»)
1209: ax.set_ylabel(«Координата вдоль общей оси»)
1210: ax.set_aspect(«equal», adjustable=»box»)
1211: ax.legend(loc=’best’, fontsize=9)
1212: fig.tight_layout()
1213: _save_or_show(fig, outpath, show)
1214:
1215:
1216: # ————————————————————
1217: # 3D ПОВЕРХНОСТИ ОБЩЕГО ОБЪЁМА
1218: # ————————————————————
1219: def plot_union_surface_vertical(order: int,
1220: global_axis: np.ndarray,
1221: lo_components: List[np.ndarray],
1222: hi_components: List[np.ndarray],
1223: f: float,
1224: R: float,
1225: m: int,
1226: h: float,
1227: outpath: Path | None = None,
1228: show: bool = True,
1229: stride: int = 2):
1230: fig = plt.figure(figsize=(10, 8))
1231: ax = fig.add_subplot(111, projection=»3d»)
1232: bounds = PlotBounds3D()
1233:
1234: for lo, hi in zip(lo_components, hi_components):
1235: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
1236: for st, en in _contiguous_segments(mask, global_axis):
1237: axis_seg = global_axis[st:en]
1238: lo_seg = lo[st:en]
1239: hi_seg = hi[st:en]
1240: _plot_surface_boundary_x(ax, axis_seg, hi_seg, bounds, stride=stride)
1241: if np.nanmax(lo_seg) > MERGE_TOL:
1242: _plot_surface_boundary_x(ax, axis_seg, lo_seg, bounds, stride=stride)
1243:
1244: _plot_side_section_on_vertical_3d(ax, global_axis, lo_components, hi_components, bounds)
1245: set_axes_equal_real_3d(ax, bounds)
1246: _hide_3d_service_axes(ax)
1247: ax.set_title(f»Псевдопараболоид: вертикальный тип, порядок {order}, 3D, f={f}, R={R}, m={m}, h={h}»)
1248: fig.tight_layout()
1249: _save_or_show(fig, outpath, show)
1250:
1251:
1252: def plot_union_surface_horizontal(order: int,
1253: global_axis: np.ndarray,
1254: lo_components: List[np.ndarray],
1255: hi_components: List[np.ndarray],
1256: shifts: np.ndarray,
1257: f: float,
1258: R: float,
1259: m: int,
1260: h: float,
1261: outpath: Path | None = None,
1262: show: bool = True,
1263: stride: int = 2):
1264: fig = plt.figure(figsize=(10, 8))
1265: ax = fig.add_subplot(111, projection=»3d»)
1266: bounds = PlotBounds3D()
1267:
1268: for lo, hi in zip(lo_components, hi_components):
1269: mask = np.isfinite(lo) & np.isfinite(hi) & (hi > lo + MERGE_TOL)
1270: for st, en in _contiguous_segments(mask, global_axis):
1271: axis_seg = global_axis[st:en]
1272: lo_seg = lo[st:en]
1273: hi_seg = hi[st:en]
1274: _plot_surface_boundary_y(ax, axis_seg, hi_seg, bounds, split_positions=shifts, stride=stride)
1275: if np.nanmax(lo_seg) > MERGE_TOL:
1276: _plot_surface_boundary_y(ax, axis_seg, lo_seg, bounds, split_positions=shifts, stride=stride)
1277:
1278: _plot_side_section_on_horizontal_3d(ax, global_axis, lo_components, hi_components, bounds)
1279: set_axes_equal_real_3d(ax, bounds)
1280: _hide_3d_service_axes(ax)
1281: ax.set_title(f»Псевдопараболоид: горизонтальный тип, порядок {order}, 3D, f={f}, R={R}, m={m}, h={h}»)
1282: fig.tight_layout()
1283: _save_or_show(fig, outpath, show)
1284:
1285:
1286: # ————————————————————
1287: # ДОПОЛНИТЕЛЬНАЯ ДИАГНОСТИКА ОБРАЗУЮЩИХ
1288: # ————————————————————
1289: def plot_base_generatrices(f: float,
1290: R: float,
1291: outpath: Path | None = None,
1292: show: bool = True,
1293: npts: int = 1200):
1294: «»»
1295: Исправленные базовые образующие 2-го порядка.
1296:
1297: Показываются только половинки образующих — по одной стороне от оси вращения,
1298: то есть именно зеркальные два сегмента парабол без дублирования второй
1299: половины полного 2D-сечения.
1300: «»»
1301: fig, axes = plt.subplots(1, 2, figsize=(15, 6.5))
1302:
1303: # —————- Вертикальный тип —————-
1304: ax = axes[0]
1305: a = parabola_a(f, R)
1306: s = np.linspace(-a, a, max(1200, npts))
1307: d = vertical_distance(np.abs(s), f, R)
1308:
1309: # Только правая половина образующей: два зеркальных сегмента парабол.
1310: ax.plot(d, s, color=’navy’, linewidth=2.0, label=’Образующая’)
1311:
1312: # Фокусы образующей.
1313: focus_points_v = [(R, -f), (R, +f)]
1314: _plot_focus_points_2d(ax, focus_points_v, rmax=max(R, float(np.nanmax(d))))
1315:
1316: # Оптимизированные служебные линии.
1317: ax.axvline(0.0, color=’black’, linestyle=’—‘, linewidth=1.05, zorder=1)
1318: ax.text(0.03, 0.97, ‘Ось вращения’, transform=ax.transAxes,
1319: ha=’left’, va=’top’, fontsize=9)
1320:
1321: ax.axvline(R, color=’black’, linestyle=’—‘, linewidth=1.0, zorder=1)
1322: ax.text(R + 0.25, 0.72 * a, ‘Фокальная ось\nобразующей’,
1323: ha=’left’, va=’center’, fontsize=9)
1324:
1325: yR = -a — 0.9
1326: ax.annotate(», xy=(R, yR), xytext=(0.0, yR),
1327: arrowprops=dict(arrowstyle='<->’, lw=1.15, color=’black’))
1328: ax.text(0.5 * R, yR + 0.28, ‘R’, ha=’center’, va=’bottom’, fontsize=12)
1329:
1330: ax.text(R — 0.35, +f, ‘+f’, ha=’right’, va=’center’, fontsize=9)
1331: ax.text(R — 0.35, -f, ‘−f’, ha=’right’, va=’center’, fontsize=9)
1332:
1333: ax.set_xlim(-0.8, 1.28 * R)
1334: ax.set_ylim(-a — 1.3, a + 1.0)
1335: ax.set_title(f’Базовые образующие 2-го порядка: вертикальный тип’)
1336: ax.set_xlabel(‘Радиальная координата’)
1337: ax.set_ylabel(‘Координата вдоль общей оси’)
1338: ax.set_aspect(‘equal’, adjustable=’box’)
1339: ax.legend(loc=’lower right’, fontsize=9)
1340:
1341: # —————- Горизонтальный тип —————-
1342: ax = axes[1]
1343: u = np.linspace(-R, R, max(1200, npts))
1344: d = horizontal_distance(u, f, R)
1345: a_h = float(np.nanmax(d))
1346:
1347: # Только правая половина образующей.
1348: ax.plot(d, u, color=’darkred’, linewidth=2.0, label=’Образующая’)
1349:
1350: # Фокусы образующей.
1351: focus_points_h = [(f, -R), (f, +R)]
1352: _plot_focus_points_2d(ax, focus_points_h, rmax=max(a_h, f))
1353:
1354: # Служебные линии по отработанной схеме 2D.
1355: ax.axvline(0.0, color=’black’, linestyle=’—‘, linewidth=1.05, zorder=1)
1356: ax.axhline(0.0, color=’black’, linestyle=’—‘, linewidth=1.0, zorder=1)
1357: ax.text(0.98, 0.52, ‘центральная горизонтальная\nось симметрии y = 0’,
1358: transform=ax.transAxes, ha=’right’, va=’bottom’, fontsize=8)
1359:
1360: ax.axhline(R, color=’black’, linestyle=’—‘, linewidth=1.0, zorder=1)
1361: ax.text(0.98, 0.88, ‘Фокальная ось\nобразующей y = R’,
1362: transform=ax.transAxes, ha=’right’, va=’top’, fontsize=9)
1363:
1364: xR = -0.72
1365: ax.annotate(», xy=(xR, R), xytext=(xR, 0.0),
1366: arrowprops=dict(arrowstyle='<->’, lw=1.15, color=’black’))
1367: ax.text(xR — 0.18, 0.5 * R, ‘R’, ha=’right’, va=’center’, fontsize=12)
1368:
1369: # Размер a показываем компактно вдоль y=0.
1370: ax.annotate(», xy=(a_h, 0.0), xytext=(0.0, 0.0),
1371: arrowprops=dict(arrowstyle='<->’, lw=1.0, color=’black’))
1372: ax.text(0.5 * a_h, 0.36, ‘a = R²/(4f)’, ha=’center’, va=’bottom’, fontsize=9)
1373:
1374: ax.set_xlim(-1.05, 1.18 * a_h)
1375: ax.set_ylim(-1.18 * R, 1.18 * R)
1376: ax.set_title(f’Базовые образующие 2-го порядка: горизонтальный тип’)
1377: ax.set_xlabel(‘Радиальная координата’)
1378: ax.set_ylabel(‘Координата вдоль общей оси’)
1379: ax.set_aspect(‘equal’, adjustable=’box’)
1380: ax.legend(loc=’lower right’, fontsize=9)
1381:
1382: fig.suptitle(f’Псевдопараболоиды 2-го порядка: базовые образующие, f={f}, R={R}’)
1383: fig.tight_layout()
1384: _save_or_show(fig, outpath, show)
1385:
1386:
1387: # ————————————————————
1388: # ОСНОВНАЯ ФУНКЦИЯ
1389: # ————————————————————
1390: def run_pseudoparaboloids(f: float = 2.0,
1391: R: float = 8.0,
1392: offsets: Sequence[float] = (7.0, 14.0),
1393: geometry_type: str = «both», # vertical / horizontal / both
1394: mode: str = «all», # section / surface / all / profile
1395: outdir: str = «pseudoparaboloids_notebook_output»,
1396: npts: int = 900,
1397: h: float = -2.0,
1398: m: int = 3,
1399: show: bool = True,
1400: make_profile: bool = True) -> Dict[str, object]:
1401: «»»
1402: Главная функция.
1403:
1404: Строит итоговый общий объём для каждого порядка от 2 до n.
1405: n = len(offsets) + 2.
1406:
1407: Параметры:
1408: f фокусное расстояние параболических ветвей;
1409: R вертикальная высота / радиус смещения оси вращения;
1410: offsets R1, R2, …, R_{n-2}; задают рекурсивные уровни;
1411: geometry_type «vertical», «horizontal» или «both»;
1412: mode «section», «surface», «all» или «profile»;
1413: m число рядов;
1414: h зазор / касание / перекрытие между рядами;
1415: show показывать графики на экране;
1416: make_profile сохранить отдельную схему образующих.
1417: «»»
1418: offsets = list(offsets)
1419: outdir_path = Path(outdir)
1420: _ensure_dir(outdir_path)
1421:
1422: if geometry_type not in {«vertical», «horizontal», «both»}:
1423: raise ValueError(«geometry_type должен быть ‘vertical’, ‘horizontal’ или ‘both'»)
1424: if mode not in {«section», «surface», «all», «profile»}:
1425: raise ValueError(«mode должен быть ‘section’, ‘surface’, ‘all’ или ‘profile'»)
1426: if f <= 0 or R <= 0:
1427: raise ValueError(«Параметры f и R должны быть положительными»)
1428: if int(m) != m or m < 1:
1429: raise ValueError(«Параметр m должен быть целым числом >= 1»)
1430:
1431: a = parabola_a(f, R)
1432: base_v = vertical_base_grid(f, R, npts=npts)
1433: base_h = horizontal_base_grid(f, R, npts=npts)
1434:
1435: result: Dict[str, object] = {
1436: «f»: float(f),
1437: «R»: float(R),
1438: «a»: float(a),
1439: «offsets»: offsets,
1440: «n_order»: len(offsets) + 2,
1441: «m»: int(m),
1442: «h»: float(h),
1443: «vertical_base_max_distance»: float(base_v[«max_distance»]),
1444: «horizontal_base_max_distance»: float(base_h[«max_distance»]),
1445: «vertical_morphology»: morphology_class_from_base(float(base_v[«max_distance»]), offsets),
1446: «horizontal_morphology»: morphology_class_from_base(float(base_h[«max_distance»]), offsets),
1447: «volume_mode»: «union_of_all_parabolic_toroidal_components»,
1448: «saved_files»: [],
1449: «warnings»: [],
1450: }
1451:
1452: print(«ПСЕВДОПАРАБОЛОИДЫ n-го ПОРЯДКА»)
1453: print(f»f = {f}»)
1454: print(f»R = {R}»)
1455: print(f»a = R^2/(4f) = {a:.6f}»)
1456: print(f»offsets = {offsets}»)
1457: print(f»Порядок n = {len(offsets) + 2}»)
1458: print(f»m = {m}»)
1459: print(f»h = {h}»)
1460: print(f»Вертикальный базовый максимум d = {base_v[‘max_distance’]:.6f}»)
1461: print(f»Горизонтальный базовый максимум d = {base_h[‘max_distance’]:.6f}»)
1462: print(«Режим построения = общий объём как объединение всех компонент»)
1463: print(«Правило: удаляется только общая перекрывающаяся часть объёма;»)
1464: print(«сами порождающие параболические тороиды как объекты построения не удаляются.»)
1465: print(«Искусственные прямые соединения не дорисовываются.»)
1466: print(f»Папка вывода = {outdir_path.resolve()}»)
1467:
1468: if make_profile or mode == «profile»:
1469: path = outdir_path / f»base_generatrices_f{f}_R{R}.png»
1470: plot_base_generatrices(f=f, R=R, outpath=path, show=show, npts=max(1200, npts))
1471: result[«saved_files»].append(str(path))
1472: if mode == «profile»:
1473: return result
1474:
1475: if geometry_type in {«vertical», «both»}:
1476: result[«warnings»].extend([f»vertical: {w}» for w in validate_offsets(base_v[«d»], offsets)])
1477:
1478: for p in range(2, len(offsets) + 3):
1479: used_offsets = offsets[:max(0, p — 2)]
1480: global_axis, lo_components, hi_components, _, shifts = vertical_order_union_components(
1481: f=f, R=R, offsets=used_offsets, npts=npts, m=m, h=h,
1482: )
1483: if mode in {«section», «all»}:
1484: path = outdir_path / f»vertical_order_{p:02d}_section_m{m}_h{h}.png»
1485: plot_union_section_vertical(
1486: order=p,
1487: global_axis=global_axis,
1488: lo_components=lo_components,
1489: hi_components=hi_components,
1490: shifts=shifts,
1491: used_offsets=used_offsets,
1492: f=f,
1493: R=R,
1494: m=m,
1495: h=h,
1496: outpath=path,
1497: show=show,
1498: )
1499: result[«saved_files»].append(str(path))
1500: if mode in {«surface», «all»}:
1501: path = outdir_path / f»vertical_order_{p:02d}_surface_m{m}_h{h}.png»
1502: plot_union_surface_vertical(
1503: order=p,
1504: global_axis=global_axis,
1505: lo_components=lo_components,
1506: hi_components=hi_components,
1507: f=f,
1508: R=R,
1509: m=m,
1510: h=h,
1511: outpath=path,
1512: show=show,
1513: )
1514: result[«saved_files»].append(str(path))
1515:
1516: if geometry_type in {«horizontal», «both»}:
1517: result[«warnings»].extend([f»horizontal: {w}» for w in validate_offsets(base_h[«d»], offsets)])
1518:
1519: for p in range(2, len(offsets) + 3):
1520: used_offsets = offsets[:max(0, p — 2)]
1521: global_axis, lo_components, hi_components, _, shifts = horizontal_order_union_components(
1522: f=f, R=R, offsets=used_offsets, npts=npts, m=m, h=h,
1523: )
1524: if mode in {«section», «all»}:
1525: path = outdir_path / f»horizontal_order_{p:02d}_section_m{m}_h{h}.png»
1526: plot_union_section_horizontal(
1527: order=p,
1528: global_axis=global_axis,
1529: lo_components=lo_components,
1530: hi_components=hi_components,
1531: shifts=shifts,
1532: used_offsets=used_offsets,
1533: f=f,
1534: R=R,
1535: m=m,
1536: h=h,
1537: outpath=path,
1538: show=show,
1539: )
1540: result[«saved_files»].append(str(path))
1541: if mode in {«surface», «all»}:
1542: path = outdir_path / f»horizontal_order_{p:02d}_surface_m{m}_h{h}.png»
1543: plot_union_surface_horizontal(
1544: order=p,
1545: global_axis=global_axis,
1546: lo_components=lo_components,
1547: hi_components=hi_components,
1548: shifts=shifts,
1549: f=f,
1550: R=R,
1551: m=m,
1552: h=h,
1553: outpath=path,
1554: show=show,
1555: )
1556: result[«saved_files»].append(str(path))
1557:
1558: if result[«warnings»]:
1559: print(«\nПРЕДУПРЕЖДЕНИЯ:»)
1560: for w in result[«warnings»]:
1561: print(» -«, w)
1562:
1563: print(«\nСОХРАНЁННЫЕ ФАЙЛЫ:»)
1564: for filename in result[«saved_files»]:
1565: print(» -«, filename)
1566:
1567: return result
1568:
1569:
1570: # ————————————————————
1571: # ПАРАМЕТРЫ ПО УМОЛЧАНИЮ
1572: # ————————————————————
1573: if __name__ == «__main__»:
1574: # Пример из исходного документа псевдопараболоидов:
1575: # f = 2, R = 8, a = R^2/(4f) = 8.
1576: f = 2.0
1577: R = 8.0
1578:
1579: # offsets задают максимальный порядок: n = len(offsets) + 2.
1580: # Для n=4 нужны два смещения: R1, R2.
1581: offsets = [7.0, 14.0]
1582:
1583: geometry_type = «both» # «vertical», «horizontal», «both»
1584: mode = «all» # «section», «surface», «all», «profile»
1585: outdir = «pseudoparaboloids_notebook_output»
1586: npts = 900
1587: m = 3
1588: h = -2.0
1589: show = True
1590:
1591: results = run_pseudoparaboloids(
1592: f=f,
1593: R=R,
1594: offsets=offsets,
1595: geometry_type=geometry_type,
1596: mode=mode,
1597: outdir=outdir,
1598: npts=npts,
1599: m=m,
1600: h=h,
1601: show=show,
1602: make_profile=True,
1603: )