import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# Настраиваем стиль

plt.style.use(‘seaborn-v0_8-whitegrid’)

plt.rcParams.update({‘font.size’: 11, ‘font.family’: ‘sans-serif’})

# ==============================================================================

# ПАРАМЕТРЫ ФИГУРЫ И ВРАЩЕНИЯ

# ==============================================================================

f = 2.0  # Фокусное расстояние парабол

R = 8.0  # Радиус смещения оси вращения ВВЕРХ

a = (R**2) / (4 * f)  # Вычисляем ширину a, при которой ветвь пересечет ось вращения

                     # (чтобы фигура замкнулась в «полюсах»)

# ==============================================================================

# ШАГ 1: БАЗОВЫЙ ПРОФИЛЬ (НАША ГАЛОЧКА)

# ==============================================================================

# Правая ветвь: y = sqrt(4fx)

x_right = np.linspace(0, a, 200)

y_right = np.sqrt(4 * f * x_right)

# Левая ветвь: y = sqrt(-4fx)

x_left = np.linspace(-a, 0, 200)

y_left = np.sqrt(-4 * f * x_left)

# Объединяем в один массив для вращения

x_full = np.concatenate((x_left, x_right))

y_full = np.concatenate((y_left, y_right))

# ==============================================================================

# ШАГ 2: РАДИУС ВРАЩЕНИЯ

# ==============================================================================

# Ось вращения поднята на R.

# Локальный радиус (расстояние от оси вращения до профиля галочки):

# r(x) = R — y_full

r_x = R — y_full

# ==============================================================================

# ВИЗУАЛИЗАЦИЯ: 2D СХЕМА И 3D ПОВЕРХНОСТЬ

# ==============================================================================

fig = plt.figure(figsize=(16, 7))

# —————————————————————————-

# ГРАФИК 1: 2D Схема Вращения

# —————————————————————————-

ax1 = fig.add_subplot(121)

# Линия фокусов (исходная ось)

ax1.axhline(0, color=’gray’, linestyle=’-.’, linewidth=1.5, label=’Линия фокусов (y=0)’)

# Новая ось вращения (сдвинутая вверх на R)

ax1.axhline(R, color=’red’, linestyle=’—‘, linewidth=2, label=f’Ось вращения (y={R})’)

# Исходная образующая («галочка»)

ax1.plot(x_full, y_full, ‘b-‘, linewidth=3, label=’Образующая (галочка)’)

# Зеркальное отражение профиля сверху (чтобы показать границы 3D фигуры в сечении)

ax1.plot(x_full, 2*R — y_full, ‘b:’, linewidth=2, alpha=0.5, label=’Зеркальный профиль’)

# Отмечаем фокусы и общую точку

ax1.plot(f, 0, ‘g*’, markersize=12, label=’Фокусы (F)’)

ax1.plot(-f, 0, ‘g*’, markersize=12)

ax1.plot(0, 0, ‘ko’, markersize=8, label=’Центральная точка’)

# Оформление 2D

ax1.set_title(‘Схема вращения (2D сечение)’, fontsize=14, fontweight=’bold’)

ax1.set_xlabel(‘Ось X’, fontweight=’bold’)

ax1.set_ylabel(‘Ось Y’, fontweight=’bold’)

ax1.axis(‘equal’)

ax1.legend(loc=’lower center’, fontsize=10, bbox_to_anchor=(0.5, -0.25), ncol=2)

ax1.grid(True, alpha=0.3)

# —————————————————————————-

# ГРАФИК 2: 3D Псевдоэллипсоид 2-го порядка

# —————————————————————————-

ax2 = fig.add_subplot(122, projection=’3d’)

# Создаем 3D сетку углов вращения

phi = np.linspace(0, 2 * np.pi, 150)

X_mesh, Phi_mesh = np.meshgrid(x_full, phi)

# Радиус для каждой точки X

R_mesh = R — np.concatenate((np.sqrt(-4 * f * x_left), np.sqrt(4 * f * x_right)))

# Переход в декартовы координаты 3D

# X остается X, а Y и Z формируют круги радиуса R_mesh

Y_mesh = R_mesh * np.cos(Phi_mesh)

Z_mesh = R_mesh * np.sin(Phi_mesh)

# Построение 3D поверхности

surf = ax2.plot_surface(X_mesh, Y_mesh, Z_mesh,

                        cmap=’coolwarm’, alpha=0.9,

                        edgecolor=’none’, antialiased=True)

# Ось вращения в 3D

ax2.plot([x_full[0], x_full[-1]], [0, 0], [0, 0], color=’red’, linestyle=’—‘, linewidth=3, label=’Ось вращения’)

# Настройка 3D осей

ax2.set_title(‘3D Псевдоэллипсоид 2-го порядка’, fontsize=14, fontweight=’bold’)

ax2.set_xlabel(‘X (вдоль оси)’)

ax2.set_ylabel(‘Y’)

ax2.set_zlabel(‘Z’)

# Выравнивание масштабов в 3D (чтобы фигура не искажалась)

max_range = np.array([X_mesh.max()-X_mesh.min(), Y_mesh.max()-Y_mesh.min(), Z_mesh.max()-Z_mesh.min()]).max() / 2.0

mid_x = (X_mesh.max()+X_mesh.min()) * 0.5

mid_y = (Y_mesh.max()+Y_mesh.min()) * 0.5

mid_z = (Z_mesh.max()+Z_mesh.min()) * 0.5

ax2.set_xlim(mid_x — max_range, mid_x + max_range)

ax2.set_ylim(mid_y — max_range, mid_y + max_range)

ax2.set_zlim(mid_z — max_range, mid_z + max_range)

# Убираем лишние панели для красоты

ax2.xaxis.pane.fill = False

ax2.yaxis.pane.fill = False

ax2.zaxis.pane.fill = False

ax2.legend(loc=’upper right’)

plt.tight_layout()

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# Настраиваем стиль

plt.style.use(‘seaborn-v0_8-darkgrid’)

plt.rcParams.update({‘font.size’: 11, ‘font.family’: ‘sans-serif’})

# ==============================================================================

# ПАРАМЕТРЫ ФИГУРЫ И ВРАЩЕНИЯ

# ==============================================================================

f = 2.0  # Фокусное расстояние парабол

R = 8.0  # Радиус смещения оси вращения

a = (R**2) / (4 * f)  # Ширина резонатора до смыкания на полюсах

# ==============================================================================

# ПРОФИЛЬ ПОВЕРХНОСТИ

# ==============================================================================

# Правая и левая ветви

x_right = np.linspace(0, a, 100)

x_left = np.linspace(-a, 0, 100)

x_full = np.concatenate((x_left, x_right))

# Формируем галочку

y_full = np.concatenate((np.sqrt(-4 * f * x_left), np.sqrt(4 * f * x_right)))

# Радиус поверхности (расстояние от оси вращения до стенки)

r_surf = R — y_full

# ==============================================================================

# 3D МОДЕЛИРОВАНИЕ

# ==============================================================================

fig = plt.figure(figsize=(14, 10))

ax = fig.add_subplot(111, projection=’3d’)

# 1. Сетка для поверхности

phi = np.linspace(0, 2 * np.pi, 100)

X_mesh, Phi_mesh = np.meshgrid(x_full, phi)

# 2. Декартовы координаты поверхности

Y_mesh = (R — np.interp(X_mesh, x_full, y_full)) * np.cos(Phi_mesh)

Z_mesh = (R — np.interp(X_mesh, x_full, y_full)) * np.sin(Phi_mesh)

# Построение поверхности (полупрозрачная, чтобы видеть сквозь нее)

surf = ax.plot_surface(X_mesh, Y_mesh, Z_mesh,

                       color=’lightblue’, alpha=0.4,

                       edgecolor=’royalblue’, linewidth=0.2, antialiased=True)

# 3. ПОСТРОЕНИЕ ВНЕШНИХ ФОКУСНЫХ КОЛЕЦ

# Левое кольцо (x = -f), радиус = R

Y_ring = R * np.cos(phi)

Z_ring = R * np.sin(phi)

X_ring_left = np.full_like(phi, -f)

# Правое кольцо (x = +f), радиус = R

X_ring_right = np.full_like(phi, f)

# Рисуем кольца (толстые, яркие)

ax.plot(X_ring_left, Y_ring, Z_ring, color=’red’, linewidth=4, label=’Левое фокусное кольцо (внешнее)’)

ax.plot(X_ring_right, Y_ring, Z_ring, color=’darkorange’, linewidth=4, label=’Правое фокусное кольцо (внешнее)’)

# Для наглядности соединим центры колец с самими кольцами тонкими лучами (как спицы)

for angle in np.linspace(0, 2*np.pi, 8):

    ax.plot([-f, -f], [0, R*np.cos(angle)], [0, R*np.sin(angle)], color=’red’, linestyle=’:’, alpha=0.5)

    ax.plot([f, f], [0, R*np.cos(angle)], [0, R*np.sin(angle)], color=’darkorange’, linestyle=’:’, alpha=0.5)

# 4. Ось вращения

ax.plot([-a*1.2, a*1.2], [0, 0], [0, 0], color=’black’, linestyle=’-.’, linewidth=2, label=’Центральная ось вращения’)

# Настройка вида

ax.set_title(‘3D Псевдоэллипсоид с Внешними Фокусными Кольцами’, fontsize=16, fontweight=’bold’, pad=20)

ax.set_xlabel(‘Ось X (продольная)’, fontweight=’bold’)

ax.set_ylabel(‘Ось Y’, fontweight=’bold’)

ax.set_zlabel(‘Ось Z’, fontweight=’bold’)

# Выравнивание масштабов 3D

max_range = np.array([X_mesh.max()-X_mesh.min(), Y_mesh.max()-Y_mesh.min(), Z_mesh.max()-Z_mesh.min()]).max() / 2.0

mid_x = 0

mid_y = 0

mid_z = 0

ax.set_xlim(mid_x — max_range, mid_x + max_range)

ax.set_ylim(mid_y — max_range, mid_y + max_range)

ax.set_zlim(mid_z — max_range, mid_z + max_range)

# Убираем серый фон панелей для «чистого» космического вида

ax.xaxis.pane.fill = False

ax.yaxis.pane.fill = False

ax.zaxis.pane.fill = False

ax.legend(loc=’upper right’, fontsize=11, framealpha=0.9)

# Угол обзора (чуть сбоку и сверху, чтобы хорошо видеть кольца)

ax.view_init(elev=20, azim=45)

plt.tight_layout()

plt.show()

Параметры фигуры:

f = 2.0

R = 8.0

a = 8.0

Длина = 16.0

Максимальный радиус = 8.0

Максимальный диаметр = 16.0

import numpy as np

import matplotlib.pyplot as plt

# =========================

# Параметры псевдоэллипсоида ГВИ

# =========================

f = 2.0

R = 8.0

# Вычисляем границу по X

a = R**2 / (4 * f)

# =========================

# 1. Образующая в плоскости XY

# y(x) = sqrt(4 f |x|) = 2*sqrt(f*|x|)

# =========================

x = np.linspace(-a, a, 1000)

y = 2 * np.sqrt(f * np.abs(x))

# Проверка размеров

length = 2 * a

max_radius = R

max_diameter = 2 * R

print(«Параметры фигуры:»)

print(f»f = {f}»)

print(f»R = {R}»)

print(f»a = {a}»)

print(f»Длина = {length}»)

print(f»Максимальный радиус = {max_radius}»)

print(f»Максимальный диаметр = {max_diameter}»)

# =========================

# 2. Построение 2D-сечения

# =========================

plt.figure(figsize=(10, 6))

plt.plot(x, y, label=r’$y(x)=2\sqrt{f|x|}$’, linewidth=2)

plt.axhline(R, linestyle=’—‘, linewidth=1.5, label=’Ось вращения: y = R’)

plt.axvline(-a, linestyle=’:’, linewidth=1)

plt.axvline(a, linestyle=’:’, linewidth=1)

# Фокусы двумерного сечения

plt.scatter([f, -f], [0, 0], s=60, label=’Фокусы’)

plt.text(f, 0.3, ‘F1′, ha=’center’)

plt.text(-f, 0.3, ‘F2′, ha=’center’)

plt.title(‘Образующая псевдоэллипсоида ГВИ в плоскости XY’)

plt.xlabel(‘X’)

plt.ylabel(‘Y’)

plt.axis(‘equal’)

plt.grid(True)

plt.legend()

plt.tight_layout()

# =========================

# 3. Построение 3D-поверхности вращения

# Вращение вокруг прямой y = R, параллельной оси X

# =========================

theta = np.linspace(0, 2 * np.pi, 180)

X, Theta = np.meshgrid(x, theta)

Y_curve = 2 * np.sqrt(f * np.abs(X))

rho = R — Y_curve   # расстояние от образующей до оси вращения

# Параметризация поверхности

Y = R + rho * np.cos(Theta)

Z = rho * np.sin(Theta)

fig = plt.figure(figsize=(12, 9))

ax = fig.add_subplot(111, projection=’3d’)

ax.plot_surface(X, Y, Z, rstride=2, cstride=2, edgecolor=’none’, alpha=0.9)

# Ось вращения

ax.plot(x, np.full_like(x, R), np.zeros_like(x), linestyle=’—‘, linewidth=2)

# Фокальные окружности при x = ±f, радиус R

phi = np.linspace(0, 2 * np.pi, 300)

# Окружность при x = +f

x_f1 = np.full_like(phi, f)

y_f1 = R + R * np.cos(phi)

z_f1 = R * np.sin(phi)

ax.plot(x_f1, y_f1, z_f1, linewidth=2)

# Окружность при x = -f

x_f2 = np.full_like(phi, -f)

y_f2 = R + R * np.cos(phi)

z_f2 = R * np.sin(phi)

ax.plot(x_f2, y_f2, z_f2, linewidth=2)

ax.set_title(‘3D-поверхность псевдоэллипсоида ГВИ’)

ax.set_xlabel(‘X’)

ax.set_ylabel(‘Y’)

ax.set_zlabel(‘Z’)

# Для более наглядной геометрии

ax.set_box_aspect([length, 2 * R, 2 * R])

plt.tight_layout()

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

# Настройка стиля

plt.style.use(‘seaborn-v0_8-darkgrid’)

plt.rcParams.update({‘font.size’: 11})

# ============================================================================

# ВХОДНЫЕ ПАРАМЕТРЫ

# ============================================================================

f = 2.0  # Фокусное расстояние парабол

R = 8.0  # Высота куполов от экватора

a = (R**2) / (4 * f)  # Экваториальный радиус (a = 8.0)

# ============================================================================

# МАТЕМАТИКА ФОРМЫ

# ============================================================================

z_vals = np.linspace(-R, R, 400)

# Радиус сечения r(z) для горизонтального псевдопараболоида

r_vals = (R — np.abs(z_vals))**2 / (4 * f)

# ============================================================================

# ВИЗУАЛИЗАЦИЯ

# ============================================================================

fig = plt.figure(figsize=(18, 6))

# —————————————————————————-

# 1. ОБРАЗУЮЩАЯ (Одна ветвь)

# —————————————————————————-

ax1 = fig.add_subplot(131)

ax1.plot(r_vals, z_vals, ‘b-‘, linewidth=3, label=’Образующая: r(z)’)

ax1.axvline(0, color=’red’, linestyle=’—‘, linewidth=2, label=’Ось вращения (Z)’)

# Геометрические фокусы парабол образующей

ax1.plot(f, R, ‘g*’, markersize=14, label=’Верхний фокус’)

ax1.plot(f, -R, ‘g*’, markersize=14, label=’Нижний фокус’)

ax1.set_title(‘1. Образующая (Параболы)’, fontsize=14, fontweight=’bold’)

ax1.set_xlabel(‘Радиус r (Ось X)’)

ax1.set_ylabel(‘Высота (Ось Z)’)

ax1.set_xlim(-1, a*1.2)

ax1.set_ylim(-R*1.2, R*1.2)

ax1.legend(loc=’lower right’)

ax1.set_aspect(‘equal’)

# —————————————————————————-

# 2. ПОЛНОЕ 2D СЕЧЕНИЕ

# —————————————————————————-

ax2 = fig.add_subplot(132)

ax2.plot(r_vals, z_vals, ‘b-‘, linewidth=2.5)

ax2.plot(-r_vals, z_vals, ‘b-‘, linewidth=2.5)

ax2.axvline(0, color=’red’, linestyle=’—‘, linewidth=1.5, alpha=0.5)

# Экваториальный внутренний аттрактор

ax2.plot(a, 0, ‘yX’, markersize=12, label=’Внутренний острый экватор’)

ax2.plot(-a, 0, ‘yX’, markersize=12)

# ВНЕШНИЕ ФОКУСЫ (4 точки в 2D сечении)

ax2.plot(f, R, ‘g*’, markersize=12, label=’Внешние фокусы’)

ax2.plot(-f, R, ‘g*’, markersize=12)

ax2.plot(f, -R, ‘g*’, markersize=12)

ax2.plot(-f, -R, ‘g*’, markersize=12)

ax2.set_title(‘2. Полное 2D Сечение’, fontsize=14, fontweight=’bold’)

ax2.set_xlabel(‘Ось X’)

ax2.set_ylabel(‘Ось Z’)

ax2.set_xlim(-a*1.2, a*1.2)

ax2.set_ylim(-R*1.2, R*1.2)

ax2.legend(loc=’lower center’, fontsize=9, bbox_to_anchor=(0.5, -0.2), ncol=2)

ax2.set_aspect(‘equal’)

# —————————————————————————-

# 3. 3D МОДЕЛЬ (С ФОКАЛЬНЫМИ КОЛЬЦАМИ)

# —————————————————————————-

ax3 = fig.add_subplot(133, projection=’3d’)

theta = np.linspace(0, 2*np.pi, 100)

Z_mesh, Theta_mesh = np.meshgrid(z_vals, theta)

R_mesh = (R — np.abs(Z_mesh))**2 / (4 * f)

X_mesh = R_mesh * np.cos(Theta_mesh)

Y_mesh = R_mesh * np.sin(Theta_mesh)

# Поверхность диска

surf = ax3.plot_surface(X_mesh, Y_mesh, Z_mesh, cmap=’plasma’, alpha=0.7, edgecolor=’none’)

# ВНУТРЕННИЙ АТТРАКТОР (Экваториальная кромка)

X_eq = a * np.cos(theta); Y_eq = a * np.sin(theta); Z_eq = np.zeros_like(theta)

ax3.plot(X_eq, Y_eq, Z_eq, color=’yellow’, linewidth=3, label=’Внутренний экватор (Аттрактор)’)

# ВНЕШНИЕ ФОКАЛЬНЫЕ КОЛЬЦА («НИМБЫ»)

# Верхнее кольцо (z = R, радиус = f)

X_halo_top = f * np.cos(theta); Y_halo_top = f * np.sin(theta); Z_halo_top = np.full_like(theta, R)

ax3.plot(X_halo_top, Y_halo_top, Z_halo_top, color=’lime’, linewidth=4, label=’Верхнее фокальное кольцо (Внешнее)’)

# Нижнее кольцо (z = -R, радиус = f)

Z_halo_bot = np.full_like(theta, -R)

ax3.plot(X_halo_top, Y_halo_top, Z_halo_bot, color=’lime’, linewidth=4, label=’Нижнее фокальное кольцо (Внешнее)’)

# Центральная ось

ax3.plot([0,0], [0,0], [-R*1.2, R*1.2], ‘r—‘, linewidth=2, alpha=0.5)

ax3.set_title(‘3. Горизонтальный Псевдопараболоид’, fontsize=14, fontweight=’bold’)

ax3.set_xlabel(‘X’)

ax3.set_ylabel(‘Y’)

ax3.set_zlabel(‘Z’)

max_range = a * 1.1

ax3.set_xlim(-max_range, max_range)

ax3.set_ylim(-max_range, max_range)

ax3.set_zlim(-max_range, max_range)

ax3.xaxis.pane.fill = False; ax3.yaxis.pane.fill = False; ax3.zaxis.pane.fill = False

ax3.legend(loc=’lower center’, fontsize=9, bbox_to_anchor=(0.5, -0.15))

plt.tight_layout()

plt.show()

МОНТЕ КАРТО ПЕРЕСЧЁТ

• 200 лучей на каждую точку сетки,
• 50 отражений на луч,
• 2 режима возбуждения:
  • центральный источник,
  • равномерное объёмное возбуждение,
• метрика удержания:
  вероятность p_run_ge5, то есть шанс получить не менее 5 подряд отражений в асимптотической клиновой зоне аттрактора.

#!/usr/bin/env python3

# -*- coding: utf-8 -*-

from __future__ import annotations

import math, os

from typing import Iterable

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from matplotlib.ticker import PercentFormatter

OUT_DIR = «mc_corrected_pseudoparaboloids_v3»

os.makedirs(OUT_DIR, exist_ok=True)

try:

    plt.style.use(«seaborn-v0_8-whitegrid»)

except Exception:

    plt.style.use(«default»)

plt.rcParams.update({«font.size»:11, «figure.figsize»:(10,7), «axes.grid»:True, «grid.alpha»:0.3})

R_BASE=10.0

U0=0.20

N_RAYS=200

N_REFL=50

BURN_IN=8

CENTER_BALL_FRAC=0.03

K_GRID=[0.05,0.10,0.125,0.1667,0.25,0.50]

TARGET_ETAS=[0.95,0.90,0.80,0.70,0.60,0.50]

def a_from(R,f): return R*R/(4.0*f)

def F_vertical(p,R,f):

    x,y,z=p

    return math.hypot(y,z) — (R — math.sqrt(max(0.0,4.0*f*abs(x))))

def F_horizontal(p,R,f):

    x,y,z=p

    return math.hypot(x,y) — ((R-abs(z))**2)/(4.0*f)

def normal_vertical(p,R,f,eps=1e-9):

    x,y,z=p

    rho=max(math.hypot(y,z),eps); ax=max(abs(x),eps)

    grad=np.array([(math.sqrt(f/ax))*(1 if x>=0 else -1), y/rho, z/rho], float)

    return grad/np.linalg.norm(grad)

def normal_horizontal(p,R,f,eps=1e-9):

    x,y,z=p

    rho=max(math.hypot(x,y),eps)

    dz=0.0 if abs(z)<eps else ((R-abs(z))/(2.0*f))*(1 if z>=0 else -1)

    grad=np.array([x/rho,y/rho,dz], float)

    return grad/np.linalg.norm(grad)

def rand_dir(rng):

    v=rng.normal(size=3); return v/np.linalg.norm(v)

def sample_center_ball(R,rng,frac=CENTER_BALL_FRAC):

    r=frac*R*(rng.uniform()**(1/3)); u=rng.uniform(-1,1); phi=rng.uniform(0,2*np.pi)

    s=math.sqrt(max(0.0,1-u*u))

    return np.array([r*s*math.cos(phi), r*s*math.sin(phi), r*u], float)

def sample_uniform_vertical(R,f,rng):

    a=a_from(R,f); x=rng.uniform(-a,a)

    rv=max(R-math.sqrt(max(0.0,4.0*f*abs(x))),0.0)

    rho=rv*math.sqrt(rng.uniform()); phi=rng.uniform(0,2*np.pi)

    return np.array([x, rho*math.cos(phi), rho*math.sin(phi)], float)

def sample_uniform_horizontal(R,f,rng):

    z=rng.uniform(-R,R); rv=((R-abs(z))**2)/(4.0*f)

    rho=rv*math.sqrt(rng.uniform()); phi=rng.uniform(0,2*np.pi)

    return np.array([rho*math.cos(phi), rho*math.sin(phi), z], float)

def hit_boundary(F,p,d,step,tmax,max_steps=30000):

    t=step

    for _ in range(max_steps):

        if F(p+t*d)>=0:

            lo,hi=t-step,t

            for __ in range(28):

                mid=0.5*(lo+hi)

                if F(p+mid*d)<=0: lo=mid

                else: hi=mid

            return 0.5*(lo+hi)

        t+=step

        if t>tmax: return None

    return None

def trace_metric(kind,R,f,regime,seed):

    rng=np.random.default_rng(seed)

    a=a_from(R,f)

    step=max(R,a)/200.0

    tmax=2.5*max(R,a)

    max_runs=[]; refl_fracs=[]

    for _ in range(N_RAYS):

        if regime==’center’:

            p=sample_center_ball(R,rng)

        else:

            p=sample_uniform_vertical(R,f,rng) if kind==’vertical’ else sample_uniform_horizontal(R,f,rng)

        d=rand_dir(rng)

        if kind==’vertical’:

            F=lambda q: F_vertical(q,R,f); normal=lambda q: normal_vertical(q,R,f)

        else:

            F=lambda q: F_horizontal(q,R,f); normal=lambda q: normal_horizontal(q,R,f)

        flags=[]

        for j in range(N_REFL):

            t=hit_boundary(F,p,d,step,tmax)

            if t is None: break

            pt=p+t*d

            if j>=BURN_IN:

                if kind==’vertical’:

                    s=a-abs(pt[0]); active=(s/a)<=U0

                else:

                    rho=math.hypot(pt[0],pt[1]); delta=a-rho; active=(delta/a)<=U0

                flags.append(active)

            n=normal(pt); d=d-2*np.dot(d,n)*n; d=d/np.linalg.norm(d); p=pt-1e-7*n

        if flags:

            refl_fracs.append(sum(flags)/len(flags))

            cur=mr=0

            for fl in flags:

                if fl: cur+=1; mr=max(mr,cur)

                else: cur=0

            max_runs.append(mr)

    max_runs=np.array(max_runs,float); refl_fracs=np.array(refl_fracs,float)

    return {‘refl_frac’:float(np.mean(refl_fracs)), ‘mean_max_run’:float(np.mean(max_runs)), ‘p_run_ge5’:float(np.mean(max_runs>=5)), ‘p_run_ge8’:float(np.mean(max_runs>=8))}

def fit_power(Ks, etas):

    best=None

    for C in np.logspace(-2,4,400):

        for p in np.linspace(0.5,4.5,401):

            pred=1.0/(1.0 + C*(Ks**p))

            err=np.mean((pred-etas)**2)

            if best is None or err<best[0]:

                best=(err,C,p)

    return best

def K_from_eta(eta,C,p): return ((1.0/eta — 1.0)/C)**(1.0/p)

def save_table_png(df, title, filename):

    fig_h=0.8+0.55*len(df)

    fig,ax=plt.subplots(figsize=(12,fig_h)); ax.axis(‘off’)

    table=ax.table(cellText=df.values,colLabels=df.columns,cellLoc=’center’,loc=’center’)

    table.auto_set_font_size(False); table.set_fontsize(10); table.scale(1.0,1.45)

    ax.set_title(title,pad=16)

    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR,filename), dpi=180, bbox_inches=’tight’); plt.close()

def main():

    rows=[]

    seedbase=0

    for regime in [‘center’,’uniform’]:

        for K in K_GRID:

            for kind in [‘vertical’,’horizontal’]:

                seedbase+=1

                f=K*R_BASE

                m=trace_metric(kind,R_BASE,f,regime,1000+seedbase)

                rows.append({‘regime’:regime,’kind’:kind,’K’:K,’f’:f,’a’:a_from(R_BASE,f),**m})

                print(regime,kind,K,m)

    df=pd.DataFrame(rows)

    df.to_csv(os.path.join(OUT_DIR,’mc_raw_long.csv’), index=False, encoding=’utf-8-sig’)

    pivot=df.pivot_table(index=[‘regime’,’K’], columns=’kind’, values=’p_run_ge5′).reset_index()

    pivot.to_csv(os.path.join(OUT_DIR,’mc_raw_pivot.csv’), index=False, encoding=’utf-8-sig’)

    fits={}

    for regime in [‘center’,’uniform’]:

        d=pivot[pivot.regime==regime]

        Ks=d[‘K’].to_numpy()

        for kind in [‘vertical’,’horizontal’]:

            fits[(regime,kind)] = fit_power(Ks, d[kind].to_numpy())

    with open(os.path.join(OUT_DIR,’fit_summary.txt’),’w’,encoding=’utf-8′) as f:

        for regime in [‘center’,’uniform’]:

            for kind in [‘vertical’,’horizontal’]:

                _,C,p=fits[(regime,kind)]

                f.write(f'{regime} {kind}: eta = 1 / (1 + {C:.4f} * K^{p:.4f})\\n’)

    tables=[]

    for regime in [‘center’,’uniform’]:

        _,Cv,pv=fits[(regime,’vertical’)]

        _,Ch,ph=fits[(regime,’horizontal’)]

        rows=[]

        for eta in TARGET_ETAS:

            Kv=K_from_eta(eta,Cv,pv); Kh=K_from_eta(eta,Ch,ph)

            rows.append({

                ‘eta_target’:eta,

                ‘eta_percent’:f'{eta*100:.0f}%’,

                ‘K_vertical’:Kv,

                ‘f_vertical’:R_BASE*Kv,

                ‘a_vertical’:a_from(R_BASE,R_BASE*Kv),

                ‘alpha_vertical_deg’:math.degrees(math.atan(2*Kv)),

                ‘K_horizontal’:Kh,

                ‘f_horizontal’:R_BASE*Kh,

                ‘a_horizontal’:a_from(R_BASE,R_BASE*Kh),

                ‘alpha_horizontal_deg’:math.degrees(math.atan(2*Kh)),

            })

        t=pd.DataFrame(rows).round(4)

        t.to_csv(os.path.join(OUT_DIR,f’mc_refined_{regime}_table.csv’), index=False, encoding=’utf-8-sig’)

        save_table_png(t, f’Monte Carlo refined table — {regime}’, f’mc_refined_{regime}_table.png’)

        tables.append((regime,t))

    Ks = np.array(K_GRID)

    for regime in [‘center’,’uniform’]:

        d=pivot[pivot.regime==regime]

        Kplot=np.logspace(np.log10(Ks.min()*0.8), np.log10(Ks.max()*1.2), 400)

        plt.figure(figsize=(10,6))

        for kind, marker in [(‘vertical’,’o’),(‘horizontal’,’s’)]:

            plt.semilogx(d[‘K’].to_numpy(), d[kind].to_numpy(), marker, label=f'{kind} raw’)

            _,C,p=fits[(regime,kind)]

            plt.semilogx(Kplot, 1/(1+C*Kplot**p), label=f'{kind} fit’)

        plt.gca().yaxis.set_major_formatter(PercentFormatter(1.0))

        plt.ylim(0,1.0); plt.xlabel(‘K = f/R’); plt.ylabel(‘η_MC (p_run_ge5)’)

        plt.title(f’Monte Carlo fit — {regime}’)

        plt.legend(); plt.grid(True,alpha=0.3); plt.tight_layout()

        plt.savefig(os.path.join(OUT_DIR,f’mc_fit_{regime}.png’), dpi=180, bbox_inches=’tight’); plt.close()

if __name__ == ‘__main__’:

    main()