はじめに
SN 1987Aの画像といえば、明るい equatorial ring (ER) と、その外側に見える northern outer ring (NOR) / southern outer ring (SOR) が印象的です。
ただ、HSTで見えているのは空に投影された2D imageです。三つのringは3D空間ではどう置かれていて、Earthから見るとなぜあの形になるのでしょうか。
そこで、Larsson et al. (2019) の図2、15を参考に、Pythonで簡単に3D可視化してみました。outer ringsについては Tziamtzis et al. (2011) も参考にしています。
今回の目的は精密な3D reconstructionではなく、ER / NOR / SORの立体配置をざっくりイメージするためのポンチ絵を作ることです。
文献値を参考にしていますが、一部には単純化した配置も含まれます。定量解析には使わず、説明用・理解用の模式図くらいに見てもらえればと思います。
triple ringの配置
Larsson et al. のsimple modelでは、三つのringを3D空間では円形として扱います。
| Ring | 傾き |
|---|---|
| ER | 43 deg |
| NOR | 45 deg |
| SOR | 38 deg |
Earthから見た奥行きでは、NORが奥側、SORが手前側です。ERは北側が手前に来る向きにしています。
3D空間の円は、中心を $\mathbf{r}_0$、半径を $R$ とすれば、例えば
\mathbf{r}(\phi)
=
\mathbf{r}_0
+
R\left(
\mathbf{e}_1\cos\phi
+
\mathbf{e}_2\sin\phi
\right)
のように置けます。
つまり今回やっていることは、3D空間に三つのcircleを置いて、それをEarth方向から眺めるというだけです。
ERのprojected semi-major axisは約 $0.82''$ としました。outer ringsの半径は
R_{\rm NOR}=2.2R_{\rm ER},
\qquad
R_{\rm SOR}=2.3R_{\rm ER}
としています。
また、Larsson et al. のsimple modelを参考に、SNからouter-ring centerまでの3D距離を
d_{\rm NOR}=1.3\times10^{18}\ {\rm cm},
\qquad
d_{\rm SOR}=1.0\times10^{18}\ {\rm cm}
として配置しました。
左が3Dで見た配置、右がEarthから見た投影です。
2D画像ではellipseに見えますが、3Dではそれぞれ別の位置にあるringだと考えると、かなりイメージしやすくなります。
なお、左図の実線／破線はringの明るさではなく、奥行きを読みやすくするための表示です。SNを通るsky planeを基準に、Earth側にある部分を実線、奥側にある部分を破線で描いています。NOR（赤）は全体として奥側にあるため、破線になっている部分が多めです。
Larsson et al. の模式図にはERとouter ringsを結ぶ構造も描かれていますが、これは仮説的な構造として置かれたもので、今回の可視化には含めていません。
視点を動かしてみる
せっかく3Dにしたので、視点も少し動かしてみました。
動いているのはringではなく視点です。
斜めから3D配置を眺めたあと、最後にEarthから見た方向へ視点を合わせています。右側はEarthから見た投影に固定しています。
途中で「どっちがEarth側だっけ？」となりやすかったので、$+Z$ をEarth方向として表示し、SNを通る基準面も薄く描いています。
視点を動かしながら見ると、2D投影だけでは少し分かりにくいNOR / SORの奥行きや、ERとの立体的な位置関係が見やすくなります。
まとめ
SN 1987Aのtriple ringをPythonで簡単に3D可視化してみました。
今回見たかったのは、
「この三つのringって、3Dではどう置かれているんだっけ？」
というところです。
論文の模式図だけでも大まかな配置は分かりますが、実際に3Dにして視点を動かしてみると、NOR / SORの奥行きやERとの位置関係がかなりイメージしやすくなりました。
あくまで可視化メモ兼ポンチ絵なので、精密な3D reconstructionや定量解析には使わず、そのくらいの軽い距離感で眺めてもらえればと思います。
完全コード
コード全文はこちら
sn1987a_triple_ring.py
#!/usr/bin/env python3
"""Educational 3D/kinematic-geometric visualization of SN 1987A.
This is not a hydrodynamic simulation and not a complete 3D reconstruction.
The three circular rings follow the simple geometry adopted by Larsson et al.
(2019), with inclinations from Tziamtzis et al. (2011). The optional straight
surfaces between the equatorial and outer rings are explicitly hypothetical.
Internal physical coordinates are right-handed and use cm:
+X = East
+Y = North
+Z = toward the observer (Earth)
Astronomical display convention is applied only at plotting time: East is left
and North is up. The observer projection uses the same +Z line of sight drawn
in the 3D panels.
Primary references
------------------
Larsson et al. 2019, ApJ, 886, 147
https://arxiv.org/abs/1910.09582
https://doi.org/10.3847/1538-4357/ab4ff2
Tziamtzis et al. 2011, A&A, 527, A35
https://arxiv.org/abs/1008.3387
https://doi.org/10.1051/0004-6361/201015576
"""
from __future__ import annotations
import argparse
import json
import math
import shutil
from dataclasses import dataclass
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FFMpegWriter, FuncAnimation, PillowWriter
from matplotlib.lines import Line2D
# =============================================================================
# Parameters and provenance
# =============================================================================
# [A] Literature adopted: Larsson et al. (2019), unless noted otherwise.
DISTANCE_KPC = 49.6
ER_INCLINATION_DEG = 43.0
NOR_INCLINATION_DEG = 45.0
SOR_INCLINATION_DEG = 38.0
ER_RADIUS_ARCSEC = 0.82 # 2018 hotspot semi-major axis
NOR_RADIUS_FACTOR = 2.2 # latest F657N image; simple-model choice
SOR_RADIUS_FACTOR = 2.3 # latest F657N image; simple-model choice
NOR_CENTER_DISTANCE_CM = 1.3e18
SOR_CENTER_DISTANCE_CM = 1.0e18
NOR_WEST_SHIFT_CM = 5.5e16
SOR_WEST_SHIFT_CM = 4.0e16
TARGET_MAJOR_AXIS_TILT_DEG = 7.0
# Near/far constraints from Larsson et al. (2019), with the ER orientation
# checked against Tziamtzis et al. (2011), Fig. 5, and earlier ring literature.
ER_NORTHERN_SIDE_IS_NEAR = True
NOR_CENTER_IS_FAR = True
SOR_CENTER_IS_NEAR = True
# [D] Larsson et al. (2019) hypothetical/simple-model layer.
SHOW_CONNECTING_MATERIAL = False
SHOW_INTERMEDIATE_RINGS = False
EJECTA_EPOCH_DAYS = 11500.0
EJECTA_SPEEDS_KMS = (7000.0, 8000.0, 9000.0)
# [C] Visualization-only settings.
SHOW_EARTH_DEPTH_CUES = True
SKY_PLANE_EXTENT_LY = 2.42
SKY_PLANE_FILL_ALPHA = 0.125
RING_SAMPLES = 720
GENERATOR_COUNT = 42
STATIC_DPI = 180
ANIMATION_DPI = 90
ANIMATION_FPS = 10
ORBIT_FRAMES = 45
ALIGN_FRAMES = 45
HOLD_FRAMES = 30
CAMERA_START_ELEV_DEG = 22.0
CAMERA_START_AZIM_DEG = -55.0
CAMERA_ORBIT_END_AZIM_DEG = -155.0
CAMERA_EARTH_ELEV_DEG = 90.0
CAMERA_EARTH_AZIM_DEG = -90.0
FOV_WIDTH_ARCSEC = 5.0
FOV_HEIGHT_ARCSEC = 5.5
# Restrained colors inspired by the paper's schematic, but not encoding
# brightness, temperature, composition, or measurement certainty.
COLOR_ER = "#777777"
COLOR_NOR = "#C84A44"
COLOR_SOR = "#3D6CA8"
COLOR_SN = "#D39A2C"
COLOR_OBSERVER = "#5B6470"
# Physical constants.
PC_CM = 3.085677581491367e18
LIGHT_YEAR_CM = 9.4607304725808e17
ARCSEC_PER_RADIAN = 206264.80624709636
DAY_S = 86400.0
KM_CM = 1.0e5
# Right-handed physical coordinate basis.
EAST_HAT = np.array([1.0, 0.0, 0.0])
NORTH_HAT = np.array([0.0, 1.0, 0.0])
LOS_HAT = np.array([0.0, 0.0, 1.0]) # SN -> observer
OBSERVER_DIRECTION = LOS_HAT.copy()
SCRIPT_DIR = Path(__file__).resolve().parent
DEFAULT_OUTPUT_DIR = (
SCRIPT_DIR if SCRIPT_DIR.name == "outputs" else SCRIPT_DIR / "outputs"
)
@dataclass(frozen=True)
class RingGeometry:
"""A genuine circle embedded in 3D physical space."""
name: str
color: str
inclination_deg: float
radius_cm: float
center_cm: np.ndarray
e1: np.ndarray
e2: np.ndarray
normal: np.ndarray
phi: np.ndarray
points_cm: np.ndarray
explicit_west_shift_cm: float
@dataclass(frozen=True)
class ModelGeometry:
rings: dict[str, RingGeometry]
system_rotation_deg: float
cm_per_arcsec: float
intermediate_rings: dict[str, dict[float, np.ndarray]]
intersection_residuals: dict[str, dict[float, float]]
# =============================================================================
# Coordinate transforms and geometry
# =============================================================================
def rotation_about_north_deg(angle_deg: float) -> np.ndarray:
"""Right-handed active rotation about physical +Y (North)."""
a = np.deg2rad(angle_deg)
c, s = np.cos(a), np.sin(a)
return np.array([[c, 0.0, s], [0.0, 1.0, 0.0], [-s, 0.0, c]])
def angular_scale_cm_per_arcsec(distance_kpc: float = DISTANCE_KPC) -> float:
distance_cm = distance_kpc * 1.0e3 * PC_CM
return distance_cm / ARCSEC_PER_RADIAN
def base_ring_frame(inclination_deg: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Return e1, e2 and normal before the small north-axis rotation.
e1 is initially East. e2 points toward the northern edge and has positive
+Z, making the northern edge of each ring the near edge. The resulting
normal points toward the southern/near polar direction.
"""
i = np.deg2rad(inclination_deg)
e1 = EAST_HAT.copy()
e2 = np.array([0.0, np.cos(i), np.sin(i)])
normal = np.cross(e1, e2)
return e1, e2, normal / np.linalg.norm(normal)
def projected_basis_metrics(e1: np.ndarray, e2: np.ndarray) -> dict[str, float]:
"""Projected semiaxis factors and major-axis angle for a unit circle."""
basis_on_sky = np.column_stack(
(
[np.dot(e1, EAST_HAT), np.dot(e1, NORTH_HAT)],
[np.dot(e2, EAST_HAT), np.dot(e2, NORTH_HAT)],
)
)
u, singular, _ = np.linalg.svd(basis_on_sky)
major_vector = u[:, 0]
if major_vector[0] < 0.0:
major_vector = -major_vector
angle = np.degrees(np.arctan2(major_vector[1], major_vector[0]))
if angle > 90.0:
angle -= 180.0
if angle <= -90.0:
angle += 180.0
return {
"major_factor": float(singular[0]),
"minor_factor": float(singular[1]),
"axis_ratio": float(singular[1] / singular[0]),
"major_axis_tilt_deg_east_to_north": float(angle),
}
def derive_system_rotation_deg(target_tilt_deg: float) -> float:
"""Derive the +Y rotation that produces the requested ER sky tilt.
Larsson et al. state the rotation axis and the observed approximately 7 deg
major-axis tilt, but not a matrix sign convention. The positive solution
is chosen because it matches their sky figures: the eastern end of the
projected major axis is north of the western end (East is displayed left).
"""
e1_0, e2_0, _ = base_ring_frame(ER_INCLINATION_DEG)
lo, hi = 0.0, 20.0
for _ in range(80):
mid = 0.5 * (lo + hi)
rot = rotation_about_north_deg(mid)
tilt = projected_basis_metrics(rot @ e1_0, rot @ e2_0)[
"major_axis_tilt_deg_east_to_north"
]
if tilt < target_tilt_deg:
lo = mid
else:
hi = mid
return 0.5 * (lo + hi)
def generate_ring(
name: str,
color: str,
inclination_deg: float,
radius_cm: float,
center_distance_cm: float,
center_polar_sign: float,
west_shift_cm: float,
system_rotation: np.ndarray,
) -> RingGeometry:
"""Generate r(phi) = center + R(e1 cos(phi) + e2 sin(phi)).
``center_polar_sign`` is +1 for the southern/near normal and -1 for
northern/far. Locating an outer-ring center along its ring normal is the
simplest convention consistent with the cited side-view geometry and with
the reported 3D center distance. This implementation convention is listed
as derived, not as an independently observed vector component.
"""
e1_0, e2_0, normal_0 = base_ring_frame(inclination_deg)
center_0 = center_polar_sign * center_distance_cm * normal_0
center = system_rotation @ center_0
center = center + np.array([-west_shift_cm, 0.0, 0.0])
e1 = system_rotation @ e1_0
e2 = system_rotation @ e2_0
normal = system_rotation @ normal_0
phi = np.linspace(0.0, 2.0 * np.pi, RING_SAMPLES, endpoint=True)
points = (
center[None, :]
+ radius_cm * np.cos(phi)[:, None] * e1[None, :]
+ radius_cm * np.sin(phi)[:, None] * e2[None, :]
)
return RingGeometry(
name=name,
color=color,
inclination_deg=inclination_deg,
radius_cm=radius_cm,
center_cm=center,
e1=e1,
e2=e2,
normal=normal,
phi=phi,
points_cm=points,
explicit_west_shift_cm=west_shift_cm,
)
def sphere_generator_intersection(
inner_points: np.ndarray,
outer_points: np.ndarray,
sphere_radius_cm: float,
) -> tuple[np.ndarray, np.ndarray, float]:
"""Intersect every straight ER-to-OR generatrix with a centered sphere."""
direction = outer_points - inner_points
qa = np.einsum("ij,ij->i", direction, direction)
qb = 2.0 * np.einsum("ij,ij->i", inner_points, direction)
qc = np.einsum("ij,ij->i", inner_points, inner_points) - sphere_radius_cm**2
disc = np.maximum(qb * qb - 4.0 * qa * qc, 0.0)
root_a = (-qb - np.sqrt(disc)) / (2.0 * qa)
root_b = (-qb + np.sqrt(disc)) / (2.0 * qa)
candidates = np.column_stack((root_a, root_b))
valid = (candidates >= -1.0e-10) & (candidates <= 1.0 + 1.0e-10)
selected = np.full(len(inner_points), np.nan)
for idx in range(len(inner_points)):
roots = candidates[idx, valid[idx]]
if roots.size:
selected[idx] = roots[np.argmin(np.abs(roots))]
if np.isnan(selected).any():
missing = int(np.isnan(selected).sum())
raise RuntimeError(f"Sphere did not intersect {missing} connecting generators")
selected = np.clip(selected, 0.0, 1.0)
points = inner_points + selected[:, None] * direction
residual = np.max(np.abs(np.linalg.norm(points, axis=1) - sphere_radius_cm))
return points, selected, float(residual / sphere_radius_cm)
def build_model() -> ModelGeometry:
cm_per_arcsec = angular_scale_cm_per_arcsec()
er_radius_cm = ER_RADIUS_ARCSEC * cm_per_arcsec
system_rotation_deg = derive_system_rotation_deg(TARGET_MAJOR_AXIS_TILT_DEG)
rotation = rotation_about_north_deg(system_rotation_deg)
rings = {
"ER": generate_ring(
"ER",
COLOR_ER,
ER_INCLINATION_DEG,
er_radius_cm,
0.0,
0.0,
0.0,
rotation,
),
"NOR": generate_ring(
"NOR",
COLOR_NOR,
NOR_INCLINATION_DEG,
NOR_RADIUS_FACTOR * er_radius_cm,
NOR_CENTER_DISTANCE_CM,
-1.0,
NOR_WEST_SHIFT_CM,
rotation,
),
"SOR": generate_ring(
"SOR",
COLOR_SOR,
SOR_INCLINATION_DEG,
SOR_RADIUS_FACTOR * er_radius_cm,
SOR_CENTER_DISTANCE_CM,
+1.0,
SOR_WEST_SHIFT_CM,
rotation,
),
}
intermediate: dict[str, dict[float, np.ndarray]] = {"NOR": {}, "SOR": {}}
residuals: dict[str, dict[float, float]] = {"NOR": {}, "SOR": {}}
for outer_name in ("NOR", "SOR"):
for speed_kms in EJECTA_SPEEDS_KMS:
sphere_radius = speed_kms * KM_CM * EJECTA_EPOCH_DAYS * DAY_S
points, _, residual = sphere_generator_intersection(
rings["ER"].points_cm, rings[outer_name].points_cm, sphere_radius
)
intermediate[outer_name][speed_kms] = points
residuals[outer_name][speed_kms] = residual
return ModelGeometry(
rings=rings,
system_rotation_deg=system_rotation_deg,
cm_per_arcsec=cm_per_arcsec,
intermediate_rings=intermediate,
intersection_residuals=residuals,
)
def project_to_sky_arcsec(points_cm: np.ndarray, cm_per_arcsec: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Orthographic observer projection using the physical basis and +Z LOS."""
east = points_cm @ EAST_HAT / cm_per_arcsec
north = points_cm @ NORTH_HAT / cm_per_arcsec
depth_toward_observer = points_cm @ LOS_HAT / cm_per_arcsec
return east, north, depth_toward_observer
# =============================================================================
# Plot helpers
# =============================================================================
def to_light_years(points_cm: np.ndarray) -> np.ndarray:
return points_cm / LIGHT_YEAR_CM
def ease_in_out(value: float) -> float:
value = float(np.clip(value, 0.0, 1.0))
return value * value * (3.0 - 2.0 * value)
def format_3d_axes(ax: plt.Axes, limit_ly: float = 2.55, show_labels: bool = True) -> None:
"""Minimal 3D styling; X is inverted only for East-left display."""
ax.set_xlim(limit_ly, -limit_ly)
ax.set_ylim(-limit_ly, limit_ly)
ax.set_zlim(-limit_ly, limit_ly)
ax.set_box_aspect((1.0, 1.0, 1.0))
ax.grid(False)
for axis in (ax.xaxis, ax.yaxis, ax.zaxis):
axis.pane.set_alpha(0.0)
axis.line.set_color("#B8BDC5")
ticks = (-2.0, 0.0, 2.0)
ax.set_xticks(ticks)
ax.set_yticks(ticks)
ax.set_zticks(ticks)
ax.tick_params(colors="#626A73", labelsize=8, pad=0)
if show_labels:
ax.set_xlabel("East (+X) [ly]", labelpad=5, color="#46505A")
ax.set_ylabel("North (+Y) [ly]", labelpad=5, color="#46505A")
ax.set_zlabel("toward observer (+Z) [ly]", labelpad=5, color="#46505A")
else:
ax.set_axis_off()
def add_coordinate_triad(ax: plt.Axes, length_ly: float = 0.72) -> list:
artists = []
origin = np.zeros(3)
for vector, label, color in (
(EAST_HAT, "E", "#3F596F"),
(NORTH_HAT, "N", "#3F596F"),
(LOS_HAT, "to observer", "#3F596F"),
):
artists.append(
ax.quiver(
*origin,
*(length_ly * vector),
color=color,
linewidth=1.25,
arrow_length_ratio=0.16,
)
)
endpoint = 1.12 * length_ly * vector
artists.append(ax.text(*endpoint, label, color=color, fontsize=8))
return artists
def plot_ring_with_earth_depth_cue(ax: plt.Axes, ring: RingGeometry) -> list:
"""Draw a ring using physical Earth depth, never animation-camera depth.
Points with physical Z >= 0 (toward Earth) are solid. Points with Z < 0
(away from Earth) are dashed. Color and opacity are unchanged, so the cue
does not encode observed brightness or emissivity.
"""
points = to_light_years(ring.points_cm)
near = points.copy()
far = points.copy()
near[points[:, 2] < 0.0] = np.nan
far[points[:, 2] >= 0.0] = np.nan
artists = []
artists.extend(
ax.plot(
far[:, 0],
far[:, 1],
far[:, 2],
color=ring.color,
lw=2.65,
ls=(0, (3.2, 2.2)),
dash_capstyle="round",
)
)
artists.extend(
ax.plot(
near[:, 0],
near[:, 1],
near[:, 2],
color=ring.color,
lw=2.65,
ls="solid",
solid_capstyle="round",
)
)
return artists
def add_sky_plane_and_center_guides(ax: plt.Axes, model: ModelGeometry) -> list:
"""Add visualization-only Earth-depth guides without changing geometry."""
artists = []
extent = SKY_PLANE_EXTENT_LY
plane_fill_color = "#D5DBE0"
plane_line_color = "#8D98A2"
# A low-alpha reference floor only: it is not a physical disk or surface.
plane_x, plane_y = np.meshgrid([-extent, extent], [-extent, extent])
artists.append(
ax.plot_surface(
plane_x,
plane_y,
np.zeros_like(plane_x),
color=plane_fill_color,
alpha=SKY_PLANE_FILL_ALPHA,
shade=False,
linewidth=0.0,
antialiased=False,
)
)
# Only four interior lines (two in each direction) to reveal perspective.
for value in (-extent / 3.0, extent / 3.0):
artists.extend(
ax.plot(
[-extent, extent],
[value, value],
[0.0, 0.0],
color=plane_line_color,
lw=0.58,
alpha=0.28,
)
)
artists.extend(
ax.plot(
[value, value],
[-extent, extent],
[0.0, 0.0],
color=plane_line_color,
lw=0.58,
alpha=0.28,
)
)
# A restrained perimeter makes the finite plotting reference explicit.
border_x = [-extent, extent, extent, -extent, -extent]
border_y = [-extent, -extent, extent, extent, -extent]
artists.extend(
ax.plot(
border_x,
border_y,
np.zeros(5),
color=plane_line_color,
lw=0.78,
alpha=0.44,
)
)
artists.append(
ax.text(
0.45 * extent,
-0.72 * extent,
0.06,
"reference sky plane (Z = 0)",
color="#66717B",
fontsize=7.2,
ha="center",
bbox=dict(boxstyle="round,pad=0.10", facecolor="white", edgecolor="none", alpha=0.64),
)
)
for name, depth_word, label_offset in (
("NOR", "far", (0.10, -0.16, -0.16)),
("SOR", "near", (-0.34, 0.12, 0.16)),
):
ring = model.rings[name]
center = to_light_years(ring.center_cm)
artists.extend(
ax.plot(
[center[0], center[0]],
[center[1], center[1]],
[0.0, center[2]],
color=ring.color,
lw=1.30,
ls=(0, (2.4, 2.0)),
alpha=0.88,
)
)
artists.append(
ax.scatter(
[center[0]],
[center[1]],
[center[2]],
s=28,
facecolor="white",
edgecolor=ring.color,
linewidth=1.0,
depthshade=False,
)
)
artists.append(
ax.scatter(
[center[0]],
[center[1]],
[0.0],
marker="o",
s=20,
facecolor="white",
edgecolor=ring.color,
linewidth=1.05,
alpha=0.90,
depthshade=False,
)
)
artists.append(
ax.text(
center[0] + label_offset[0],
center[1] + label_offset[1],
center[2] + label_offset[2],
f"{name} · {depth_word}",
color=ring.color,
fontsize=7.4,
bbox=dict(boxstyle="round,pad=0.12", facecolor="white", edgecolor="none", alpha=0.76),
)
)
return artists
def add_fixed_earth_depth_annotations(ax: plt.Axes, style_note_y: float = 0.025) -> list:
"""Fixed panel notes that separate physical Earth depth from camera view."""
direction_note = ax.text2D(
0.98,
0.955,
"Physical Z (camera-independent)\n"
"+Z toward Earth / near\n"
"-Z away from Earth / far",
transform=ax.transAxes,
ha="right",
va="top",
fontsize=6.8,
color="#3F4851",
linespacing=1.25,
bbox=dict(boxstyle="round,pad=0.28", facecolor="white", edgecolor="#D7DCE1", lw=0.55, alpha=0.90),
)
style_note = ax.text2D(
0.02,
style_note_y,
"solid / dashed = Earth-depth cue only\n(not brightness)",
transform=ax.transAxes,
ha="left",
va="bottom",
fontsize=7.1,
color="#58626C",
linespacing=1.2,
bbox=dict(boxstyle="round,pad=0.24", facecolor="white", edgecolor="none", alpha=0.82),
)
return [direction_note, style_note]
def add_sky_direction_arrows(ax: plt.Axes) -> None:
arrow = dict(arrowstyle="-|>", color="#37414A", lw=1.1, mutation_scale=9)
ax.annotate("", xy=(0.84, 0.12), xytext=(0.92, 0.12), xycoords="axes fraction", arrowprops=arrow)
ax.text(0.815, 0.105, "E", transform=ax.transAxes, fontsize=9, color="#37414A")
ax.annotate("", xy=(0.92, 0.22), xytext=(0.92, 0.12), xycoords="axes fraction", arrowprops=arrow)
ax.text(0.905, 0.235, "N", transform=ax.transAxes, fontsize=9, color="#37414A")
def plot_connecting_generators_3d(ax: plt.Axes, er: RingGeometry, outer: RingGeometry) -> None:
indices = np.linspace(0, len(er.points_cm) - 1, GENERATOR_COUNT, endpoint=False).astype(int)
for idx in indices:
segment = np.vstack((er.points_cm[idx], outer.points_cm[idx])) / LIGHT_YEAR_CM
ax.plot(
segment[:, 0],
segment[:, 1],
segment[:, 2],
color=outer.color,
ls=(0, (2.2, 3.2)),
lw=0.75,
alpha=0.22,
)
def plot_scene_3d(
ax: plt.Axes,
model: ModelGeometry,
show_connecting: bool,
show_intermediate: bool,
show_axes: bool = True,
show_legend: bool = True,
show_depth_cues: bool = False,
) -> dict[str, object]:
rings = model.rings
for ring in rings.values():
if show_depth_cues:
plot_ring_with_earth_depth_cue(ax, ring)
else:
p = to_light_years(ring.points_cm)
ax.plot(p[:, 0], p[:, 1], p[:, 2], color=ring.color, lw=2.65, solid_capstyle="round")
depth_artists = []
if show_depth_cues:
depth_artists.extend(add_sky_plane_and_center_guides(ax, model))
depth_artists.extend(
add_fixed_earth_depth_annotations(
ax,
style_note_y=0.12 if show_axes else 0.025,
)
)
if show_connecting:
plot_connecting_generators_3d(ax, rings["ER"], rings["NOR"])
plot_connecting_generators_3d(ax, rings["ER"], rings["SOR"])
if show_intermediate:
alphas = (0.45, 0.65, 0.9)
for outer_name in ("NOR", "SOR"):
for alpha, speed in zip(alphas, EJECTA_SPEEDS_KMS):
p = to_light_years(model.intermediate_rings[outer_name][speed])
ax.plot(
p[:, 0],
p[:, 1],
p[:, 2],
color=rings[outer_name].color,
lw=1.35,
ls=(0, (3.2, 2.6)),
alpha=alpha,
)
sn = ax.scatter([0.0], [0.0], [0.0], marker="*", s=95, color=COLOR_SN, edgecolor="white", linewidth=0.7, zorder=20)
observer_distance = 2.38
observer_line = ax.plot(
[0.0, 0.0],
[0.0, 0.0],
[0.0, observer_distance],
color=COLOR_OBSERVER,
lw=1.05,
ls=(0, (3.0, 3.2)),
alpha=0.66,
)[0]
observer = ax.scatter(
[0.0],
[0.0],
[observer_distance],
s=52,
color=COLOR_OBSERVER,
edgecolor="white",
linewidth=0.8,
alpha=0.88,
)
observer_text = ax.text(0.0, 0.0, observer_distance + 0.13, "Earth / observer", color=COLOR_OBSERVER, fontsize=8, ha="center")
observer_arrow = ax.quiver(
0.0,
0.0,
observer_distance - 0.56,
0.0,
0.0,
0.43,
color=COLOR_OBSERVER,
linewidth=1.05,
arrow_length_ratio=0.28,
)
format_3d_axes(ax, show_labels=show_axes)
if show_axes and not show_depth_cues:
add_coordinate_triad(ax)
if show_legend:
handles = [
Line2D([0], [0], color=COLOR_ER, lw=2.7, label="ER"),
Line2D([0], [0], color=COLOR_NOR, lw=2.7, label="NOR (far-side center)"),
Line2D([0], [0], color=COLOR_SOR, lw=2.7, label="SOR (near-side center)"),
]
if show_connecting:
handles.append(
Line2D(
[0],
[0],
color="#7C8792",
lw=1.1,
ls=(0, (2.2, 3.2)),
label="hypothetical connecting material",
)
)
if show_intermediate:
handles.append(
Line2D(
[0],
[0],
color="#7C8792",
lw=1.3,
ls=(0, (3.2, 2.6)),
label="model sphere intersections (IRs)",
)
)
ax.legend(handles=handles, loc="upper left", frameon=False, fontsize=8)
return {
"sn": sn,
"observer_line": observer_line,
"observer": observer,
"observer_text": observer_text,
"observer_arrow": observer_arrow,
"depth_artists": depth_artists,
}
def plot_sky_projection(
ax: plt.Axes,
model: ModelGeometry,
show_connecting: bool,
show_intermediate: bool,
show_legend: bool = True,
) -> None:
rings = model.rings
for ring in rings.values():
east, north, _ = project_to_sky_arcsec(ring.points_cm, model.cm_per_arcsec)
ax.plot(east, north, color=ring.color, lw=2.65, solid_capstyle="round")
if show_connecting:
indices = np.linspace(0, RING_SAMPLES - 1, GENERATOR_COUNT, endpoint=False).astype(int)
for outer_name in ("NOR", "SOR"):
outer = rings[outer_name]
for idx in indices:
segment = np.vstack((rings["ER"].points_cm[idx], outer.points_cm[idx]))
east, north, _ = project_to_sky_arcsec(segment, model.cm_per_arcsec)
ax.plot(east, north, color=outer.color, lw=0.7, ls=(0, (2.2, 3.2)), alpha=0.18)
if show_intermediate:
alphas = (0.42, 0.64, 0.9)
for outer_name in ("NOR", "SOR"):
for alpha, speed in zip(alphas, EJECTA_SPEEDS_KMS):
east, north, _ = project_to_sky_arcsec(
model.intermediate_rings[outer_name][speed], model.cm_per_arcsec
)
ax.plot(
east,
north,
color=rings[outer_name].color,
lw=1.3,
ls=(0, (3.2, 2.6)),
alpha=alpha,
)
ax.scatter([0.0], [0.0], marker="*", s=86, color=COLOR_SN, edgecolor="white", linewidth=0.7, zorder=20)
ax.text(0.08, 0.08, "SN", fontsize=8, color="#775D20")
for name, offset in (("NOR", (-0.07, 0.18)), ("SOR", (-0.07, -0.28)), ("ER", (0.10, -0.08))):
center = rings[name].center_cm
east, north, _ = project_to_sky_arcsec(center[None, :], model.cm_per_arcsec)
ax.text(east[0] + offset[0], north[0] + offset[1], name, color=rings[name].color, fontsize=9, weight="semibold")
# Astronomical display convention: positive East is deliberately on left.
ax.set_xlim(FOV_WIDTH_ARCSEC / 2.0, -FOV_WIDTH_ARCSEC / 2.0)
ax.set_ylim(-FOV_HEIGHT_ARCSEC / 2.0, FOV_HEIGHT_ARCSEC / 2.0)
ax.set_aspect("equal", adjustable="box")
ax.set_xlabel("RA offset (East +) [arcsec]")
ax.set_ylabel("Dec offset (North +) [arcsec]")
ax.grid(color="#D8DDE2", lw=0.55, alpha=0.55)
ax.spines[["top", "right"]].set_visible(False)
add_sky_direction_arrows(ax)
if show_legend:
handles = [
Line2D([0], [0], color=COLOR_ER, lw=2.7, label="ER"),
Line2D([0], [0], color=COLOR_NOR, lw=2.7, label="NOR"),
Line2D([0], [0], color=COLOR_SOR, lw=2.7, label="SOR"),
]
if show_connecting:
handles.append(
Line2D([0], [0], color="#7C8792", lw=1.0, ls=(0, (2.2, 3.2)), label="hypothetical connectors")
)
if show_intermediate:
handles.append(
Line2D([0], [0], color="#7C8792", lw=1.2, ls=(0, (3.2, 2.6)), label="model IRs")
)
ax.legend(handles=handles, loc="upper right", frameon=False, fontsize=8)
# =============================================================================
# Static figures and animation
# =============================================================================
def create_static_figure(
model: ModelGeometry,
output_path: Path,
show_connecting: bool,
show_intermediate: bool,
) -> None:
fig = plt.figure(figsize=(15.0, 7.25), facecolor="white")
grid = fig.add_gridspec(1, 2, width_ratios=(1.10, 1.0), wspace=0.16)
ax3d = fig.add_subplot(grid[0, 0], projection="3d")
ax2d = fig.add_subplot(grid[0, 1])
plot_scene_3d(
ax3d,
model,
show_connecting,
show_intermediate,
show_axes=True,
show_legend=True,
show_depth_cues=SHOW_EARTH_DEPTH_CUES,
)
ax3d.view_init(elev=22.0, azim=-62.0)
ax3d.set_title("3D geometry", pad=8, fontsize=13, weight="semibold")
plot_sky_projection(ax2d, model, show_connecting, show_intermediate, show_legend=True)
ax2d.set_title("Projection as viewed from Earth", pad=12, fontsize=13, weight="semibold")
fig.suptitle("SN 1987A triple-ring system — simple kinematic-geometric model", fontsize=16, y=0.975, weight="semibold")
caveat = (
"Circular rings: literature-adopted simple geometry. "
+ (
"Dashed connectors are hypothetical, not an observed 3D surface."
if show_connecting
else "Hypothetical ER–OR connecting material is not shown (toggle is off)."
)
)
fig.text(0.5, 0.018, caveat, ha="center", va="bottom", fontsize=9, color="#5E6670")
fig.savefig(output_path, dpi=STATIC_DPI, bbox_inches="tight", facecolor="white")
plt.close(fig)
def create_depth_representative_frame(
model: ModelGeometry,
output_path: Path,
show_connecting: bool,
show_intermediate: bool,
) -> None:
"""Save a representative frame using the same layout as the animation."""
fig = plt.figure(figsize=(12.0, 5.6), facecolor="white")
grid = fig.add_gridspec(1, 2, width_ratios=(1.08, 1.0), wspace=0.12)
ax3d = fig.add_subplot(grid[0, 0], projection="3d")
ax2d = fig.add_subplot(grid[0, 1])
plot_scene_3d(
ax3d,
model,
show_connecting,
show_intermediate,
show_axes=False,
show_legend=False,
show_depth_cues=True,
)
ax3d.view_init(elev=22.0, azim=-68.0)
ax3d.text2D(
0.03,
0.955,
"Camera: oblique view",
transform=ax3d.transAxes,
va="top",
fontsize=8.1,
color="#3F4851",
bbox=dict(boxstyle="round,pad=0.28", facecolor="white", edgecolor="none", alpha=0.82),
)
plot_sky_projection(ax2d, model, show_connecting, show_intermediate, show_legend=False)
ax2d.set_title("Fixed Earth-view projection", fontsize=11.5, weight="semibold")
fig.suptitle("SN 1987A: physical Earth-depth cues in the 3D view", y=0.975, fontsize=14, weight="semibold")
fig.text(
0.5,
0.018,
"The camera is arbitrary; near/far is always defined by physical Z relative to the Z=0 sky plane.",
ha="center",
fontsize=8.8,
color="#5E6670",
)
fig.savefig(output_path, dpi=STATIC_DPI, bbox_inches="tight", facecolor="white")
plt.close(fig)
def create_schematic(
model: ModelGeometry,
output_path: Path,
show_connecting: bool,
show_intermediate: bool,
) -> None:
fig = plt.figure(figsize=(9.0, 7.4), facecolor="white")
ax = fig.add_subplot(111, projection="3d")
plot_scene_3d(ax, model, show_connecting, show_intermediate, show_axes=False, show_legend=True)
ax.view_init(elev=18.0, azim=-60.0)
ax.set_title("SN 1987A triple-ring geometry", fontsize=16, weight="semibold", pad=2)
fig.text(
0.5,
0.035,
"Camera view is schematic; rings are 3D circles. Earth lies along physical +Z.",
ha="center",
fontsize=9,
color="#5E6670",
)
if show_connecting:
fig.text(
0.5,
0.012,
"Dashed ER–OR lines: hypothetical connecting material in the Larsson et al. simple model.",
ha="center",
fontsize=8.5,
color="#6A7178",
)
fig.savefig(output_path, dpi=STATIC_DPI, bbox_inches="tight", facecolor="white")
plt.close(fig)
def camera_for_frame(frame: int) -> tuple[float, float, str]:
if frame < ORBIT_FRAMES:
t = frame / max(ORBIT_FRAMES - 1, 1)
eased = ease_in_out(t)
azim = CAMERA_START_AZIM_DEG + (CAMERA_ORBIT_END_AZIM_DEG - CAMERA_START_AZIM_DEG) * eased
elev = CAMERA_START_ELEV_DEG + 6.0 * np.sin(np.pi * t)
return elev, azim, "Camera orbiting\n(rings and physical +Z remain fixed)"
if frame < ORBIT_FRAMES + ALIGN_FRAMES:
t = (frame - ORBIT_FRAMES) / max(ALIGN_FRAMES - 1, 1)
eased = ease_in_out(t)
elev = CAMERA_START_ELEV_DEG + (CAMERA_EARTH_ELEV_DEG - CAMERA_START_ELEV_DEG) * eased
azim = CAMERA_ORBIT_END_AZIM_DEG + (CAMERA_EARTH_AZIM_DEG - CAMERA_ORBIT_END_AZIM_DEG) * eased
return elev, azim, "Camera moving toward Earth LOS (+Z)"
return CAMERA_EARTH_ELEV_DEG, CAMERA_EARTH_AZIM_DEG, "Camera aligned with Earth LOS\nNorth up, East left"
def create_animation(
model: ModelGeometry,
gif_path: Path,
mp4_path: Path,
show_connecting: bool,
show_intermediate: bool,
make_mp4: bool,
) -> bool:
fig = plt.figure(figsize=(12.0, 5.6), facecolor="white")
grid = fig.add_gridspec(1, 2, width_ratios=(1.08, 1.0), wspace=0.12)
ax3d = fig.add_subplot(grid[0, 0], projection="3d")
ax2d = fig.add_subplot(grid[0, 1])
plot_scene_3d(
ax3d,
model,
show_connecting,
show_intermediate,
show_axes=False,
show_legend=False,
show_depth_cues=True,
)
plot_sky_projection(ax2d, model, show_connecting, show_intermediate, show_legend=False)
ax2d.set_title("Fixed Earth-view projection", fontsize=11.5, weight="semibold")
status = ax3d.text2D(
0.02,
0.96,
"",
transform=ax3d.transAxes,
va="top",
fontsize=8.1,
color="#3F4851",
bbox=dict(boxstyle="round,pad=0.3", facecolor="white", edgecolor="none", alpha=0.82),
)
fig.suptitle("SN 1987A: moving the camera from 3D to the observer's view", y=0.975, fontsize=14, weight="semibold")
fig.text(
0.5,
0.018,
"Only the camera moves; all physical ring coordinates remain fixed.",
ha="center",
fontsize=9,
color="#5E6670",
)
total_frames = ORBIT_FRAMES + ALIGN_FRAMES + HOLD_FRAMES
def update(frame: int):
elev, azim, label = camera_for_frame(frame)
ax3d.view_init(elev=elev, azim=azim)
status.set_text(label)
return (status,)
animation = FuncAnimation(fig, update, frames=total_frames, interval=1000.0 / ANIMATION_FPS, blit=False)
animation.save(gif_path, writer=PillowWriter(fps=ANIMATION_FPS), dpi=ANIMATION_DPI)
mp4_written = False
if make_mp4 and shutil.which("ffmpeg"):
try:
animation.save(
mp4_path,
writer=FFMpegWriter(
fps=ANIMATION_FPS,
codec="libx264",
bitrate=2200,
extra_args=[
"-vf",
"pad=ceil(iw/2)*2:ceil(ih/2)*2",
"-pix_fmt",
"yuv420p",
"-movflags",
"+faststart",
],
),
dpi=ANIMATION_DPI,
)
mp4_written = True
except Exception as exc: # GIF remains a valid fallback.
mp4_path.unlink(missing_ok=True)
print(f"Warning: MP4 encoding failed; GIF was retained ({exc})")
else:
# Avoid leaving a stale MP4 from an earlier run with different settings.
mp4_path.unlink(missing_ok=True)
plt.close(fig)
return mp4_written
# =============================================================================
# Validation
# =============================================================================
def add_check(
checks: list[dict[str, object]],
check_id: str,
description: str,
passed: bool,
measured: object,
expected: object,
tolerance: object | None = None,
) -> None:
entry: dict[str, object] = {
"id": check_id,
"description": description,
"status": "PASS" if bool(passed) else "FAIL",
"measured": measured,
"expected": expected,
}
if tolerance is not None:
entry["tolerance"] = tolerance
checks.append(entry)
def validate_model(model: ModelGeometry) -> dict[str, object]:
rings = model.rings
checks: list[dict[str, object]] = []
metrics = {name: projected_basis_metrics(ring.e1, ring.e2) for name, ring in rings.items()}
for name in ("ER", "NOR", "SOR"):
expected_ratio = abs(math.cos(math.radians(rings[name].inclination_deg)))
measured_ratio = metrics[name]["axis_ratio"]
add_check(
checks,
f"{name.lower()}_axis_ratio",
f"{name} projected axis ratio is approximately cos(inclination)",
abs(measured_ratio - expected_ratio) < 0.01,
measured_ratio,
expected_ratio,
0.01,
)
er_major_arcsec = rings["ER"].radius_cm * metrics["ER"]["major_factor"] / model.cm_per_arcsec
add_check(
checks,
"er_projected_radius",
"ER projected semi-major radius",
abs(er_major_arcsec - ER_RADIUS_ARCSEC) < 1.0e-6,
er_major_arcsec,
ER_RADIUS_ARCSEC,
1.0e-6,
)
nor_ratio = rings["NOR"].radius_cm / rings["ER"].radius_cm
sor_ratio = rings["SOR"].radius_cm / rings["ER"].radius_cm
add_check(checks, "nor_radius_ratio", "NOR/ER intrinsic radius ratio", abs(nor_ratio - NOR_RADIUS_FACTOR) < 1.0e-12, nor_ratio, NOR_RADIUS_FACTOR, 1.0e-12)
add_check(checks, "sor_radius_ratio", "SOR/ER intrinsic radius ratio", abs(sor_ratio - SOR_RADIUS_FACTOR) < 1.0e-12, sor_ratio, SOR_RADIUS_FACTOR, 1.0e-12)
nor_depth = float(np.dot(rings["NOR"].center_cm, LOS_HAT))
sor_depth = float(np.dot(rings["SOR"].center_cm, LOS_HAT))
add_check(checks, "nor_far_side", "NOR center is on the far side (-Z)", nor_depth < 0.0, nor_depth, "< 0 cm")
add_check(checks, "sor_near_side", "SOR center is on the near side (+Z)", sor_depth > 0.0, sor_depth, "> 0 cm")
er = rings["ER"]
north_index = int(np.argmax(er.points_cm[:, 1]))
north_edge_depth = float(er.points_cm[north_index, 2])
add_check(checks, "er_northern_edge_near", "ER northern edge is nearer the observer", north_edge_depth > 0.0, north_edge_depth, "> 0 cm")
for name, expected_shift in (("NOR", NOR_WEST_SHIFT_CM), ("SOR", SOR_WEST_SHIFT_CM)):
explicit_east_component = -rings[name].explicit_west_shift_cm
add_check(
checks,
f"{name.lower()}_west_shift_sign",
f"{name} explicit center correction points West (-X)",
explicit_east_component < 0.0,
explicit_east_component,
f"-{expected_shift:.3e} cm along +X/East",
)
measured_tilt = metrics["ER"]["major_axis_tilt_deg_east_to_north"]
add_check(
checks,
"major_axis_tilt",
"ER projected major-axis tilt (east end north of west end)",
abs(measured_tilt - TARGET_MAJOR_AXIS_TILT_DEG) < 0.01,
measured_tilt,
TARGET_MAJOR_AXIS_TILT_DEG,
0.01,
)
er_radius_ly = rings["ER"].radius_cm / LIGHT_YEAR_CM
add_check(
checks,
"physical_angular_scale",
"0.82 arcsec at 49.6 kpc is consistent with the approximate 0.6 ly description",
abs(er_radius_ly - 0.6) / 0.6 < 0.10,
er_radius_ly,
"approximately 0.6 ly",
"10% (the literature value is rounded)",
)
add_check(
checks,
"los_identity",
"3D observer arrow and 2D projection use the same LOS vector",
bool(np.array_equal(OBSERVER_DIRECTION, LOS_HAT)),
OBSERVER_DIRECTION.tolist(),
LOS_HAT.tolist(),
)
add_check(
checks,
"sky_display_convention",
"Plot limits encode North up and East left",
True,
{"x_limits_arcsec": [FOV_WIDTH_ARCSEC / 2.0, -FOV_WIDTH_ARCSEC / 2.0], "y_limits_arcsec": [-FOV_HEIGHT_ARCSEC / 2.0, FOV_HEIGHT_ARCSEC / 2.0]},
"descending East axis; ascending North axis",
)
maximum_intersection_residual = max(
residual for side in model.intersection_residuals.values() for residual in side.values()
)
add_check(
checks,
"intermediate_ring_sphere_intersections",
"Optional IR curves are numerical intersections, not hand-placed rings",
maximum_intersection_residual < 1.0e-12,
maximum_intersection_residual,
"relative sphere-radius residual < 1e-12",
1.0e-12,
)
center_summary = {}
for name in ("NOR", "SOR"):
east, north, depth = project_to_sky_arcsec(rings[name].center_cm[None, :], model.cm_per_arcsec)
center_summary[name] = {
"east_arcsec": float(east[0]),
"north_arcsec": float(north[0]),
"depth_toward_observer_arcsec_equivalent": float(depth[0]),
"center_distance_cm_after_explicit_west_shift": float(np.linalg.norm(rings[name].center_cm)),
}
output = {
"model_scope": "Educational kinematic-geometric visualization; not hydrodynamics and not a complete 3D reconstruction.",
"coordinate_system": {
"handedness": "right-handed",
"+X": "East",
"+Y": "North",
"+Z": "toward observer",
"observer_projection": "orthographic along +Z onto the X-Y sky plane",
"display": "North up, East left (x-axis inverted only when plotted)",
},
"sources": {
"primary": "Larsson et al. 2019, ApJ 886, 147, doi:10.3847/1538-4357/ab4ff2",
"supporting": "Tziamtzis et al. 2011, A&A 527, A35, doi:10.1051/0004-6361/201015576",
},
"derived": {
"cm_per_arcsec": model.cm_per_arcsec,
"er_radius_cm": rings["ER"].radius_cm,
"er_radius_light_year": er_radius_ly,
"nor_radius_arcsec": rings["NOR"].radius_cm / model.cm_per_arcsec,
"sor_radius_arcsec": rings["SOR"].radius_cm / model.cm_per_arcsec,
"system_rotation_about_north_deg": model.system_rotation_deg,
"outer_ring_projected_centers": center_summary,
"ejecta_sphere_radii_cm": {
str(int(speed)): speed * KM_CM * EJECTA_EPOCH_DAYS * DAY_S for speed in EJECTA_SPEEDS_KMS
},
},
"projected_ring_metrics": metrics,
"checks": checks,
"all_passed": all(check["status"] == "PASS" for check in checks),
}
return output
def write_validation(validation: dict[str, object], output_path: Path) -> None:
output_path.write_text(json.dumps(validation, indent=2, ensure_ascii=False) + "\n", encoding="utf-8")
print("\nValidation")
print("=" * 78)
for check in validation["checks"]:
print(f"{check['status']:>4} {check['id']:<42} {check['description']}")
print("=" * 78)
print(f"all_passed = {validation['all_passed']}")
# =============================================================================
# Command-line entry point
# =============================================================================
def parse_args() -> argparse.Namespace:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--output-dir", type=Path, default=DEFAULT_OUTPUT_DIR, help="Directory for generated artifacts")
parser.add_argument(
"--show-connecting-material",
action=argparse.BooleanOptionalAction,
default=SHOW_CONNECTING_MATERIAL,
help="Show the explicitly hypothetical straight ER-to-OR connecting material",
)
parser.add_argument(
"--show-intermediate-rings",
action=argparse.BooleanOptionalAction,
default=SHOW_INTERMEDIATE_RINGS,
help="Show sphere/generator intersections for 7000, 8000 and 9000 km/s ejecta",
)
parser.add_argument("--skip-animation", action="store_true", help="Generate static figures and validation only")
parser.add_argument("--no-mp4", action="store_true", help="Do not attempt MP4 output")
return parser.parse_args()
def main() -> None:
args = parse_args()
output_dir = args.output_dir.resolve()
output_dir.mkdir(parents=True, exist_ok=True)
model = build_model()
print("SN 1987A triple-ring simple geometry")
print(f"Internal coordinates: +X East, +Y North, +Z toward observer; unit = cm")
print(f"Angular scale: {model.cm_per_arcsec:.6e} cm arcsec^-1 at {DISTANCE_KPC:.1f} kpc")
print(f"Derived ER radius: {model.rings['ER'].radius_cm:.6e} cm = {model.rings['ER'].radius_cm / LIGHT_YEAR_CM:.4f} ly")
print(f"Derived +Y system rotation: {model.system_rotation_deg:.6f} deg for a {TARGET_MAJOR_AXIS_TILT_DEG:.1f} deg sky tilt")
print(f"Hypothetical connecting material shown: {args.show_connecting_material}")
print(f"Intermediate model rings shown: {args.show_intermediate_rings}")
validation = validate_model(model)
write_validation(validation, output_dir / "validation.json")
if not validation["all_passed"]:
raise RuntimeError("One or more geometry validations failed; outputs were not finalized")
create_static_figure(
model,
output_dir / "triple_ring_3d_projection.png",
args.show_connecting_material,
args.show_intermediate_rings,
)
create_depth_representative_frame(
model,
output_dir / "triple_ring_camera_depth_frame.png",
args.show_connecting_material,
args.show_intermediate_rings,
)
create_schematic(
model,
output_dir / "triple_ring_schematic.png",
args.show_connecting_material,
args.show_intermediate_rings,
)
mp4_written = False
if not args.skip_animation:
mp4_written = create_animation(
model,
output_dir / "triple_ring_camera.gif",
output_dir / "triple_ring_camera.mp4",
args.show_connecting_material,
args.show_intermediate_rings,
make_mp4=not args.no_mp4,
)
print("\nGenerated")
print(f"- {output_dir / 'triple_ring_3d_projection.png'}")
print(f"- {output_dir / 'triple_ring_camera_depth_frame.png'}")
print(f"- {output_dir / 'triple_ring_schematic.png'}")
print(f"- {output_dir / 'validation.json'}")
if not args.skip_animation:
print(f"- {output_dir / 'triple_ring_camera.gif'}")
if mp4_written:
print(f"- {output_dir / 'triple_ring_camera.mp4'}")
else:
print("- MP4 skipped: ffmpeg was unavailable or --no-mp4 was used")
if __name__ == "__main__":
main()
基本的にはスクリプトを保存して、
python sn1987a_triple_ring.py
と実行すれば、本文で示した静止画とanimationを生成できます。
ギャラリー
今回作った静止画やanimationをまとめています。
https://drive.google.com/drive/folders/1WPtAmV6bZw7AcRhoM09FSLp5ASTT-z8J?usp=sharing
研究紹介や発表資料、説明用のポンチ絵など、使えそうなものがあればご自由にお使いください。
ただし、本文でも書いた通り精密な3D reconstructionや定量解析を目的としたものではないので、その点だけご注意ください。