mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
680 lines
29 KiB
Python
680 lines
29 KiB
Python
#!/usr/bin/env python3
|
||
"""
|
||
Pulsar Marble-Jar Model — Multiscale Event-Driven Simulation
|
||
|
||
Physical basis:
|
||
A neutron star = 10^57 neutrons (marbles) in a superfluid (jar).
|
||
Vortices are pinned to crustal nuclei. Mutual friction drains rotation.
|
||
Individual unpinning events are invisible; collectively they spin-down the star.
|
||
When pinning fails catastrophically: glitch = phase-transition flash.
|
||
|
||
Timescales:
|
||
Cruise (years) : magnetic braking + vortex creep
|
||
Glitch (seconds) : avalanche, angular momentum transfer
|
||
Recovery (days) : vortices repin, crust relaxes
|
||
|
||
Units:
|
||
Time : years (cruise) / seconds (glitch)
|
||
Frequency : Hz
|
||
Energy : 10^33 J (one "unit")
|
||
Inertia : 10^38 kg m^2
|
||
|
||
No "up/down" — only torsional direction (rotation axis).
|
||
Genus-3 surface embeds the vortex array topology.
|
||
"""
|
||
|
||
import numpy as np
|
||
import matplotlib.pyplot as plt
|
||
from matplotlib.patches import Circle, FancyBboxPatch
|
||
from mpl_toolkits.mplot3d import Axes3D
|
||
from matplotlib.colors import Normalize
|
||
from matplotlib.cm import ScalarMappable
|
||
import json
|
||
import os
|
||
from datetime import datetime
|
||
from dataclasses import dataclass, field
|
||
from typing import List, Tuple
|
||
from scipy.ndimage import gaussian_filter
|
||
|
||
# ─── Physical Constants (realistic, scaled) ──────────────────────────────────
|
||
YEAR_TO_SEC = 3.154e7
|
||
SEC_TO_YEAR = 1.0 / YEAR_TO_SEC
|
||
|
||
I_CRUST = 0.6 # Crust moment of inertia (fraction of total)
|
||
I_SF = 1.0 # Superfluid moment of inertia
|
||
I_TOTAL = I_CRUST + I_SF
|
||
|
||
B_FIELD = 1e12 # Surface B-field in Gauss (scaled)
|
||
R_STAR = 1.0 # Radius in 10 km units
|
||
|
||
TAU_CHAR = 1e6 # Characteristic spin-down age (years)
|
||
TAU_RECOVERY = 30.0 # Post-glitch recovery time (days)
|
||
TAU_CREEP = 500.0 # Vortex creep coupling timescale (years)
|
||
|
||
T_PIN_BASE = 1.0 # Base pinning threshold
|
||
T_PIN_WIDTH = 0.2 # Distribution width (some pins weaker)
|
||
GLITCH_DF = 0.05 # Fraction of pinned vortices that unpin in a glitch
|
||
|
||
# ─── Genus-3 Surface (unchanged topology) ────────────────────────────────────
|
||
def genus3_surface(u, v, R=3.0, r=1.0, p=0.6):
|
||
rho = R + r * np.cos(v) + p * np.cos(3 * u)
|
||
x = rho * np.cos(u)
|
||
y = rho * np.sin(u)
|
||
z = r * np.sin(v) + 0.3 * np.sin(3 * u)
|
||
return x, y, z
|
||
|
||
|
||
def torsional_depth(u, v):
|
||
"""Local 'depth' = torsional curvature concentration. Higher near lobe centers."""
|
||
d0 = np.minimum(np.abs(u - 0), np.abs(u - 2*np.pi))
|
||
d1 = np.abs(u - 2*np.pi/3)
|
||
d2 = np.abs(u - 4*np.pi/3)
|
||
d_min = np.minimum(np.minimum(d0, d1), d2)
|
||
depth = 3.0 / (1.0 + 5.0 * d_min) * (1.0 + 0.3 * np.cos(v))
|
||
return depth
|
||
|
||
|
||
def genus3_mesh(n_u=80, n_v=40):
|
||
u = np.linspace(0, 2*np.pi, n_u)
|
||
v = np.linspace(0, 2*np.pi, n_v)
|
||
U, V = np.meshgrid(u, v)
|
||
X, Y, Z = genus3_surface(U, V)
|
||
return X, Y, Z, U, V
|
||
|
||
|
||
# ─── Vortex ──────────────────────────────────────────────────────────────────
|
||
@dataclass
|
||
class Vortex:
|
||
u: float
|
||
v: float
|
||
circulation: float = 1.0
|
||
pinned: bool = True
|
||
tension: float = 0.0
|
||
unpin_age: float = 0.0 # When did it last unpin?
|
||
repin_time: float = 0.0 # How long until it repins?
|
||
|
||
def pos(self):
|
||
x, y, z = genus3_surface(self.u, self.v)
|
||
return np.array([x, y, z])
|
||
|
||
|
||
# ─── Flash Event ─────────────────────────────────────────────────────────────
|
||
@dataclass
|
||
class Flash:
|
||
time_yr: float # Time in years
|
||
time_sec: float # Time within glitch in seconds
|
||
energy: float # Flash energy released
|
||
n_unpinned: int # How many vortices unpin
|
||
delta_omega_hz: float # Glitch size (spin-up)
|
||
recovery_time_days: float
|
||
tier: int # Which tier boundary crossed
|
||
|
||
|
||
# ─── Pulsar Marble-Jar Model ─────────────────────────────────────────────────
|
||
class PulsarMarbleJar:
|
||
"""Event-driven multiscale pulsar simulation."""
|
||
|
||
def __init__(self, n_vortices=512, f0_hz=10.0):
|
||
# Two-component rotation
|
||
self.f_crust = f0_hz # Crust spin frequency (Hz)
|
||
self.f_sf = f0_hz * 1.0001 # Superfluid slightly ahead (vortex creep lag)
|
||
self.omega_crust = 2 * np.pi * self.f_crust
|
||
self.omega_sf = 2 * np.pi * self.f_sf
|
||
|
||
# Conserved angular momentum (in units where I_total = 1)
|
||
self.L_total = I_CRUST * self.omega_crust + I_SF * self.omega_sf
|
||
|
||
# Vortex array
|
||
self.vortices: List[Vortex] = []
|
||
self._init_vortices(n_vortices)
|
||
|
||
# State
|
||
self.time_yr = 0.0
|
||
self.in_glitch = False
|
||
self.glitch_start_time_yr = 0.0
|
||
self.glitch_clock_sec = 0.0
|
||
self.recovery_clock_days = 0.0
|
||
|
||
# Events
|
||
self.flashes: List[Flash] = []
|
||
|
||
# History (cruise sampling)
|
||
self.history = {
|
||
't_yr': [],
|
||
'f_crust_hz': [],
|
||
'f_sf_hz': [],
|
||
'E_rot': [],
|
||
'L_total': [],
|
||
'n_pinned': [],
|
||
'tension_max': [],
|
||
'flash_count': [],
|
||
'phase': [] # 'cruise', 'glitch', 'recovery'
|
||
}
|
||
|
||
# Tier boundaries for flash classification
|
||
self.tier_thresholds = [0.5, 1.0, 2.0] # Energy thresholds
|
||
|
||
# Cruise sampling interval
|
||
self.cruise_sample_interval = 100.0 # years
|
||
self.last_sample_time = -self.cruise_sample_interval
|
||
|
||
def _init_vortices(self, n):
|
||
"""Initialize vortices with clustering near lobe centers."""
|
||
for i in range(n):
|
||
if np.random.rand() < 0.35:
|
||
lobe = np.random.choice(3)
|
||
u_c = lobe * 2 * np.pi / 3
|
||
u = u_c + np.random.normal(0, 0.25)
|
||
v = np.pi + np.random.normal(0, 0.4)
|
||
else:
|
||
u = np.random.uniform(0, 2*np.pi)
|
||
v = np.random.uniform(0, 2*np.pi)
|
||
|
||
# Each vortex has slightly different pinning strength
|
||
pin_variation = np.random.normal(0, T_PIN_WIDTH)
|
||
|
||
self.vortices.append(Vortex(
|
||
u=u % (2*np.pi),
|
||
v=v % (2*np.pi),
|
||
circulation=1.0 + 0.5 * np.random.poisson(0.3),
|
||
pinned=True,
|
||
tension=0.0,
|
||
repin_time=np.inf
|
||
))
|
||
|
||
def _vortex_density_grid(self, n_grid=32):
|
||
"""Compute smoothed vortex density on toroidal grid."""
|
||
density = np.zeros((n_grid, n_grid))
|
||
for v in self.vortices:
|
||
iu = int((v.u / (2*np.pi)) * n_grid) % n_grid
|
||
iv = int((v.v / (2*np.pi)) * n_grid) % n_grid
|
||
density[iu, iv] += v.circulation
|
||
density = gaussian_filter(density, sigma=1.2, mode='wrap')
|
||
return density
|
||
|
||
def _update_vortex_tensions(self):
|
||
"""Compute tension on each vortex from Magnus force and local density."""
|
||
density = self._vortex_density_grid()
|
||
n_grid = density.shape[0]
|
||
|
||
max_tension = 0.0
|
||
for v in self.vortices:
|
||
iu = int((v.u / (2*np.pi)) * n_grid) % n_grid
|
||
iv = int((v.v / (2*np.pi)) * n_grid) % n_grid
|
||
local_depth = torsional_depth(v.u, v.v)
|
||
local_rho = density[iu, iv]
|
||
|
||
# Magnus force: F_M ~ ρ_s κ × (v_sf - v_crust)
|
||
# Simplified: tension proportional to differential rotation and local depth
|
||
delta_omega = abs(self.omega_sf - self.omega_crust)
|
||
magnus = delta_omega * v.circulation * (1.0 + 0.8 * local_depth)
|
||
repulsion = 0.2 * local_rho ** 1.5
|
||
v.tension = magnus + repulsion
|
||
max_tension = max(max_tension, v.tension)
|
||
|
||
return max_tension
|
||
|
||
def _enforce_L_conservation(self):
|
||
"""Hard correct to ensure L is exactly conserved."""
|
||
L_now = I_CRUST * self.omega_crust + I_SF * self.omega_sf
|
||
delta_L = self.L_total - L_now
|
||
domega = delta_L / I_TOTAL
|
||
self.omega_crust += domega
|
||
self.omega_sf += domega
|
||
self.f_crust = self.omega_crust / (2 * np.pi)
|
||
self.f_sf = self.omega_sf / (2 * np.pi)
|
||
|
||
# ── Cruise Phase (years) ────────────────────────────────────────────────
|
||
def cruise_step(self, dt_yr):
|
||
"""Integrate spin-down between glitches. Large timestep in years."""
|
||
dt = dt_yr
|
||
|
||
# 1. Magnetic dipole braking: dΩ/dt = -Ω / τ_char * (Ω/Ω_0)^2
|
||
# Characteristic age τ_char ~ 10^6 yr for canonical pulsar
|
||
braking = -self.omega_crust / (TAU_CHAR * YEAR_TO_SEC) * (self.omega_crust / (2*np.pi*10))**2
|
||
braking *= YEAR_TO_SEC # convert to rad/s per year
|
||
self.omega_crust += braking * dt
|
||
|
||
# 2. Vortex creep: superfluid slowly couples to crust
|
||
# Mutual friction torque transfers angular momentum over τ_creep
|
||
delta_omega = self.omega_sf - self.omega_crust
|
||
coupling = delta_omega / TAU_CREEP
|
||
self.omega_sf -= coupling * dt * (I_CRUST / I_SF)
|
||
self.omega_crust += coupling * dt
|
||
|
||
# 3. Vortices slowly drift (creep motion)
|
||
for v in self.vortices:
|
||
if v.pinned:
|
||
# Pinned vortices barely move — thermal creep
|
||
du = 0.001 * np.sin(3*v.u) * dt / TAU_CREEP
|
||
dv = 0.0005 * np.cos(v.v) * dt / TAU_CREEP
|
||
else:
|
||
# Unpinned vortices move with superfluid
|
||
du = 0.1 * self.omega_sf * dt * SEC_TO_YEAR
|
||
dv = 0.02 * dt * SEC_TO_YEAR
|
||
# Check repin
|
||
v.unpin_age += dt * YEAR_TO_SEC
|
||
if v.unpin_age > v.repin_time:
|
||
v.pinned = True
|
||
v.unpin_age = 0
|
||
|
||
v.u = (v.u + du) % (2*np.pi)
|
||
v.v = (v.v + dv) % (2*np.pi)
|
||
|
||
self._enforce_L_conservation()
|
||
self.time_yr += dt
|
||
|
||
# 4. Check for glitch trigger
|
||
max_tension = self._update_vortex_tensions()
|
||
if max_tension > T_PIN_BASE and not self.in_glitch:
|
||
# Enter glitch phase
|
||
self.in_glitch = True
|
||
self.glitch_start_time_yr = self.time_yr
|
||
self.glitch_clock_sec = 0.0
|
||
return True # signal: glitch triggered
|
||
|
||
# 5. Sample history
|
||
if self.time_yr - self.last_sample_time >= self.cruise_sample_interval:
|
||
self._record_history('cruise')
|
||
self.last_sample_time = self.time_yr
|
||
|
||
return False
|
||
|
||
# ── Glitch Phase (seconds) ──────────────────────────────────────────────
|
||
def glitch_step(self, dt_sec):
|
||
"""Resolve glitch dynamics at second-scale resolution."""
|
||
dt = dt_sec
|
||
self.glitch_clock_sec += dt
|
||
|
||
# Count how many vortices exceed threshold
|
||
triggered = [v for v in self.vortices if v.pinned and v.tension > T_PIN_BASE]
|
||
|
||
if len(triggered) > 0 and self.glitch_clock_sec < 10.0:
|
||
# Avalanche: unpin a fraction of triggered vortices
|
||
n_to_unpin = max(1, int(GLITCH_DF * len(triggered)))
|
||
np.random.shuffle(triggered)
|
||
unpinned = triggered[:n_to_unpin]
|
||
|
||
# Each unpinned vortex transfers angular momentum from sf → crust
|
||
# Real glitch: ΔΩ/Ω ~ 10^-8 to 10^-6
|
||
# Each simulated vortex = ~10^14 real vortices
|
||
total_dL = 0.0
|
||
for v in unpinned:
|
||
v.pinned = False
|
||
v.unpin_age = 0.0
|
||
v.repin_time = np.random.exponential(TAU_RECOVERY * 24 * 3600) # seconds
|
||
# Small angular momentum transfer per vortex
|
||
dL = v.circulation * self.omega_sf * 1e-6 * I_SF
|
||
total_dL += dL
|
||
|
||
# Apply transfer
|
||
self.omega_crust += total_dL / I_CRUST
|
||
self.omega_sf -= total_dL / I_SF
|
||
|
||
# Flash energy: rotational energy change + thermal dissipation
|
||
E_flash = 0.5 * total_dL * (self.omega_crust - self.omega_sf)
|
||
|
||
# Classify tier
|
||
tier = sum(1 for th in self.tier_thresholds if E_flash > th)
|
||
|
||
# Record flash
|
||
self.flashes.append(Flash(
|
||
time_yr=self.time_yr + self.glitch_clock_sec * SEC_TO_YEAR,
|
||
time_sec=self.glitch_clock_sec,
|
||
energy=E_flash,
|
||
n_unpinned=len(unpinned),
|
||
delta_omega_hz=total_dL / (2*np.pi * I_CRUST),
|
||
recovery_time_days=TAU_RECOVERY,
|
||
tier=tier
|
||
))
|
||
|
||
# Update tensions after redistribution
|
||
self._update_vortex_tensions()
|
||
|
||
# Vortices move rapidly during glitch
|
||
for v in self.vortices:
|
||
if not v.pinned:
|
||
du = 0.5 * self.omega_sf * dt
|
||
dv = 0.1 * np.sin(v.v) * dt
|
||
v.u = (v.u + du) % (2*np.pi)
|
||
v.v = (v.v + dv) % (2*np.pi)
|
||
|
||
self._enforce_L_conservation()
|
||
|
||
# End glitch after a few seconds
|
||
if self.glitch_clock_sec > 5.0:
|
||
self.in_glitch = False
|
||
self.recovery_clock_days = 0.0
|
||
return True # signal: glitch ended, enter recovery
|
||
|
||
return False
|
||
|
||
# ── Recovery Phase (days) ───────────────────────────────────────────────
|
||
def recovery_step(self, dt_days):
|
||
"""Post-glitch relaxation. Vortices slowly repin."""
|
||
self.recovery_clock_days += dt_days
|
||
|
||
# During recovery, crust and superfluid re-equilibrate
|
||
delta_omega = self.omega_sf - self.omega_crust
|
||
# Recovery follows exponential relaxation with timescale τ_recovery
|
||
relax_rate = dt_days / TAU_RECOVERY
|
||
self.omega_crust += delta_omega * relax_rate * (I_SF / I_TOTAL)
|
||
self.omega_sf -= delta_omega * relax_rate * (I_CRUST / I_TOTAL)
|
||
|
||
# Repin vortices whose time is up
|
||
for v in self.vortices:
|
||
if not v.pinned:
|
||
v.unpin_age += dt_days * 24 * 3600
|
||
if v.unpin_age > v.repin_time:
|
||
v.pinned = True
|
||
v.unpin_age = 0
|
||
|
||
self._enforce_L_conservation()
|
||
self.time_yr += dt_days / 365.25
|
||
|
||
# End recovery
|
||
if self.recovery_clock_days > TAU_RECOVERY:
|
||
return True # signal: back to cruise
|
||
|
||
return False
|
||
|
||
def _record_history(self, phase):
|
||
"""Record state to history buffer."""
|
||
self.history['t_yr'].append(self.time_yr)
|
||
self.history['f_crust_hz'].append(self.f_crust)
|
||
self.history['f_sf_hz'].append(self.f_sf)
|
||
E_rot = 0.5 * I_CRUST * self.omega_crust**2 + 0.5 * I_SF * self.omega_sf**2
|
||
self.history['E_rot'].append(E_rot)
|
||
self.history['L_total'].append(self.L_total)
|
||
self.history['n_pinned'].append(sum(1 for v in self.vortices if v.pinned))
|
||
self.history['tension_max'].append(max(v.tension for v in self.vortices))
|
||
self.history['flash_count'].append(len(self.flashes))
|
||
self.history['phase'].append(phase)
|
||
|
||
# ── Main Run Loop ───────────────────────────────────────────────────────
|
||
def run(self, t_max_yr=2e6):
|
||
"""Run full multiscale simulation."""
|
||
print(f"Starting multiscale simulation: 0 → {t_max_yr:.1e} years")
|
||
print(f"Initial f_crust = {self.f_crust:.2f} Hz")
|
||
print(f"Characteristic age = {TAU_CHAR:.1e} years")
|
||
print(f"Vortices: {len(self.vortices)}")
|
||
|
||
next_print = 0.0
|
||
|
||
while self.time_yr < t_max_yr:
|
||
if not self.in_glitch:
|
||
# Cruise or recovery
|
||
if self.recovery_clock_days > 0:
|
||
# In recovery
|
||
done = self.recovery_step(dt_days=1.0)
|
||
if done:
|
||
self.recovery_clock_days = 0.0
|
||
else:
|
||
# In cruise
|
||
glitch_triggered = self.cruise_step(dt_yr=10.0)
|
||
if glitch_triggered:
|
||
print(f" GLITCH TRIGGERED at t = {self.time_yr:.2e} yr")
|
||
print(f" Max tension = {max(v.tension for v in self.vortices):.3f}")
|
||
else:
|
||
# In glitch — resolve at second scale
|
||
done = self.glitch_step(dt_sec=0.1)
|
||
if done:
|
||
print(f" Glitch resolved: Δf = {self.flashes[-1].delta_omega_hz:.2e} Hz")
|
||
print(f" {self.flashes[-1].n_unpinned} vortices unpinned")
|
||
|
||
# Progress print
|
||
if self.time_yr >= next_print:
|
||
phase = 'GLITCH' if self.in_glitch else ('RECOVERY' if self.recovery_clock_days > 0 else 'CRUISE')
|
||
print(f" t={self.time_yr:.2e} yr | f={self.f_crust:.4f} Hz | "
|
||
f"flashes={len(self.flashes)} | phase={phase}")
|
||
next_print += t_max_yr / 20
|
||
|
||
# Final sample
|
||
self._record_history('end')
|
||
|
||
print(f"\n{'='*60}")
|
||
print("SIMULATION COMPLETE")
|
||
print(f" Total time: {self.time_yr:.2e} years")
|
||
print(f" Final f_crust: {self.f_crust:.4f} Hz")
|
||
print(f" Total flashes: {len(self.flashes)}")
|
||
print(f" Final L deviation: {abs(I_CRUST*self.omega_crust + I_SF*self.omega_sf - self.L_total):.2e}")
|
||
|
||
if self.flashes:
|
||
sizes = [f.delta_omega_hz for f in self.flashes]
|
||
print(f" Glitch sizes: min={min(sizes):.2e}, max={max(sizes):.2e} Hz")
|
||
intervals = np.diff([f.time_yr for f in self.flashes])
|
||
if len(intervals) > 0:
|
||
print(f" Inter-glitch: mean={np.mean(intervals):.2e}, med={np.median(intervals):.2e} yr")
|
||
|
||
|
||
# ─── Visualization ───────────────────────────────────────────────────────────
|
||
def visualize(model, out_dir="/home/allaun/Documents/Research Stack/out"):
|
||
os.makedirs(out_dir, exist_ok=True)
|
||
|
||
fig = plt.figure(figsize=(18, 14))
|
||
|
||
# ── Panel A: Spin-down curve (log-log) ─────────────────────────────────
|
||
ax1 = fig.add_subplot(3, 3, 1)
|
||
t = np.array(model.history['t_yr'])
|
||
f = np.array(model.history['f_crust_hz'])
|
||
ax1.loglog(t + 1, f, 'b-', linewidth=1, label='Crust spin')
|
||
ax1.loglog(t + 1, np.array(model.history['f_sf_hz']), 'r--', alpha=0.5, linewidth=0.8, label='Superfluid')
|
||
|
||
# Mark glitches
|
||
for flash in model.flashes:
|
||
ax1.axvline(flash.time_yr + 1, color='orange', alpha=0.4, linewidth=0.8)
|
||
|
||
# Reference: pure dipole braking f ∝ t^(-1/2)
|
||
t_ref = np.logspace(0, np.log10(t.max()+1), 50)
|
||
f_ref = model.f_crust * (t_ref / (t_ref[0] + 1e3)) ** (-0.5)
|
||
ax1.loglog(t_ref, f_ref, 'g:', alpha=0.5, label='∝ t^(-1/2) dipole')
|
||
|
||
ax1.set_xlabel('Time (years)')
|
||
ax1.set_ylabel('Spin frequency (Hz)')
|
||
ax1.set_title('Spin-down over megayears')
|
||
ax1.legend()
|
||
ax1.set_ylim(bottom=0.1)
|
||
|
||
# ── Panel B: Glitch sizes histogram ────────────────────────────────────
|
||
ax2 = fig.add_subplot(3, 3, 2)
|
||
if model.flashes:
|
||
sizes = [f.delta_omega_hz for f in model.flashes]
|
||
ax2.hist(np.log10(sizes), bins=20, color='steelblue', edgecolor='black', alpha=0.7)
|
||
ax2.set_xlabel('log₁₀(Δf) [Hz]')
|
||
ax2.set_ylabel('Count')
|
||
ax2.set_title(f'Glitch size distribution (n={len(sizes)})')
|
||
else:
|
||
ax2.text(0.5, 0.5, 'No glitches', ha='center', va='center', transform=ax2.transAxes)
|
||
|
||
# ── Panel C: Inter-glitch intervals ────────────────────────────────────
|
||
ax3 = fig.add_subplot(3, 3, 3)
|
||
if len(model.flashes) > 1:
|
||
intervals = np.diff([f.time_yr for f in model.flashes])
|
||
ax3.hist(np.log10(intervals + 1), bins=20, color='forestgreen', edgecolor='black', alpha=0.7)
|
||
ax3.set_xlabel('log₁₀(Δt) [years]')
|
||
ax3.set_ylabel('Count')
|
||
ax3.set_title('Inter-glitch waiting times')
|
||
else:
|
||
ax3.text(0.5, 0.5, 'Need >1 glitch', ha='center', va='center', transform=ax3.transAxes)
|
||
|
||
# ── Panel D: Energy budget over time ───────────────────────────────────
|
||
ax4 = fig.add_subplot(3, 3, 4)
|
||
E_rot = np.array(model.history['E_rot'])
|
||
ax4.semilogy(t + 1, E_rot, 'k-', linewidth=1.5, label='Rotational E')
|
||
# Approximate radiated energy = integral of braking torque
|
||
# E_dot = I Ω dΩ/dt ≈ -I Ω² / τ_char
|
||
E_rad = E_rot[0] * (1 - (f / model.history['f_crust_hz'][0])**2)
|
||
ax4.semilogy(t + 1, E_rot[0] - E_rad, 'm--', alpha=0.5, label='Radiated away')
|
||
ax4.set_xlabel('Time (years)')
|
||
ax4.set_ylabel('Energy')
|
||
ax4.set_title('Energy budget')
|
||
ax4.legend()
|
||
|
||
# ── Panel E: Angular momentum ──────────────────────────────────────────
|
||
ax5 = fig.add_subplot(3, 3, 5)
|
||
L = np.array(model.history['L_total'])
|
||
ax5.plot(t, L, 'k-', linewidth=1.5)
|
||
ax5.set_xlabel('Time (years)')
|
||
ax5.set_ylabel('Angular momentum')
|
||
ax5.set_title(f'L conserved (σ={np.std(L):.2e})')
|
||
|
||
# ── Panel F: Vortex tension over time ──────────────────────────────────
|
||
ax6 = fig.add_subplot(3, 3, 6)
|
||
ax6.semilogy(t + 1, model.history['tension_max'], 'purple', linewidth=1)
|
||
ax6.axhline(T_PIN_BASE, color='red', linestyle='--', label='Pin threshold')
|
||
ax6.set_xlabel('Time (years)')
|
||
ax6.set_ylabel('Max tension')
|
||
ax6.set_title('Peak vortex tension')
|
||
ax6.legend()
|
||
|
||
# ── Panel G: 3D Genus-3 with final vortex state ────────────────────────
|
||
ax7 = fig.add_subplot(3, 3, 7, projection='3d')
|
||
X, Y, Z, U, V = genus3_mesh(n_u=60, n_v=30)
|
||
depth_map = torsional_depth(U, V)
|
||
ax7.plot_surface(X, Y, Z, facecolors=plt.cm.RdYlBu_r(depth_map / depth_map.max()),
|
||
alpha=0.3, rstride=2, cstride=2, linewidth=0.1)
|
||
|
||
vtx_pos = np.array([v.pos() for v in model.vortices])
|
||
tensions = np.array([v.tension for v in model.vortices])
|
||
pinned = np.array([v.pinned for v in model.vortices])
|
||
|
||
if len(vtx_pos) > 0:
|
||
# Color by tension, shape by pinned status
|
||
scatter = ax7.scatter(vtx_pos[:, 0], vtx_pos[:, 1], vtx_pos[:, 2],
|
||
c=tensions, cmap='hot', s=20 + 60*tensions/T_PIN_BASE,
|
||
alpha=0.8, edgecolors='black', linewidth=0.2)
|
||
plt.colorbar(scatter, ax=ax7, shrink=0.5, label='Tension')
|
||
|
||
ax7.set_title('Final vortex array on genus-3')
|
||
|
||
# ── Panel H: Pinned vs unpinned over time ──────────────────────────────
|
||
ax8 = fig.add_subplot(3, 3, 8)
|
||
ax8.plot(t, np.array(model.history['n_pinned']), 'b-', label='Pinned')
|
||
ax8.plot(t, len(model.vortices) - np.array(model.history['n_pinned']), 'r--', label='Unpinned')
|
||
ax8.set_xlabel('Time (years)')
|
||
ax8.set_ylabel('Count')
|
||
ax8.set_title('Vortex pinning state')
|
||
ax8.legend()
|
||
|
||
# ── Panel I: Flash timeline with tiers ─────────────────────────────────
|
||
ax9 = fig.add_subplot(3, 3, 9)
|
||
if model.flashes:
|
||
ft = [f.time_yr for f in model.flashes]
|
||
fe = [f.energy for f in model.flashes]
|
||
fc = [f.n_unpinned for f in model.flashes]
|
||
tiers = [f.tier for f in model.flashes]
|
||
|
||
colors = ['lightblue', 'yellow', 'orange', 'red']
|
||
c_list = [colors[min(t, 3)] for t in tiers]
|
||
|
||
ax9.scatter(ft, fe, c=c_list, s=[20 + 5*c for c in fc],
|
||
alpha=0.7, edgecolors='black', linewidth=0.5)
|
||
ax9.set_xlabel('Time (years)')
|
||
ax9.set_ylabel('Flash energy')
|
||
ax9.set_title('Phase transition flashes')
|
||
ax9.set_yscale('log')
|
||
else:
|
||
ax9.text(0.5, 0.5, 'No flashes', ha='center', va='center', transform=ax9.transAxes)
|
||
|
||
plt.tight_layout()
|
||
out_path = os.path.join(out_dir, f"marble_jar_{datetime.now().strftime('%Y%m%d_%H%M%S')}.png")
|
||
plt.savefig(out_path, dpi=150, bbox_inches='tight')
|
||
plt.close()
|
||
print(f"Saved: {out_path}")
|
||
|
||
# Second figure: Zoomed glitch detail (if any glitches)
|
||
if model.flashes:
|
||
fig2, axes = plt.subplots(2, 3, figsize=(15, 8))
|
||
|
||
for idx, ax in enumerate(axes.flat):
|
||
if idx >= min(6, len(model.flashes)):
|
||
ax.axis('off')
|
||
continue
|
||
|
||
flash = model.flashes[idx]
|
||
# Mock detailed profile: exponential recovery
|
||
t_g = np.linspace(0, TAU_RECOVERY, 200)
|
||
delta_f = flash.delta_omega_hz * np.exp(-t_g / TAU_RECOVERY)
|
||
|
||
ax.plot(t_g, delta_f * 1e6, 'b-', linewidth=1.5)
|
||
ax.axvline(flash.time_sec / 86400, color='red', linestyle='--', alpha=0.5, label='Glitch')
|
||
ax.set_xlabel('Days after glitch')
|
||
ax.set_ylabel('Δf (μHz)')
|
||
ax.set_title(f'Glitch {idx+1}: t={flash.time_yr:.2e} yr, tier={flash.tier}')
|
||
ax.set_yscale('log')
|
||
|
||
plt.tight_layout()
|
||
out_path2 = os.path.join(out_dir, f"marble_jar_glitchdetail_{datetime.now().strftime('%Y%m%d_%H%M%S')}.png")
|
||
plt.savefig(out_path2, dpi=150, bbox_inches='tight')
|
||
plt.close()
|
||
print(f"Saved: {out_path2}")
|
||
return out_path, out_path2
|
||
|
||
return out_path, None
|
||
|
||
|
||
# ─── Export ──────────────────────────────────────────────────────────────────
|
||
def export_data(model, out_dir="/home/allaun/Documents/Research Stack/out"):
|
||
os.makedirs(out_dir, exist_ok=True)
|
||
|
||
data = {
|
||
'metadata': {
|
||
'model': 'PulsarMarbleJarMultiscale',
|
||
'topology': 'genus-3',
|
||
'n_vortices': len(model.vortices),
|
||
't_max_yr': model.time_yr,
|
||
'parameters': {
|
||
'I_crust': I_CRUST,
|
||
'I_sf': I_SF,
|
||
'tau_char': TAU_CHAR,
|
||
'tau_creep': TAU_CREEP,
|
||
'tau_recovery': TAU_RECOVERY,
|
||
'T_pin_base': T_PIN_BASE,
|
||
'glitch_df': GLITCH_DF
|
||
}
|
||
},
|
||
'history': {k: [float(x) if isinstance(x, (int, float, np.floating)) else x
|
||
for x in v]
|
||
for k, v in model.history.items()},
|
||
'flashes': [
|
||
{'time_yr': f.time_yr, 'time_sec': f.time_sec,
|
||
'energy': f.energy, 'n_unpinned': f.n_unpinned,
|
||
'delta_omega_hz': f.delta_omega_hz,
|
||
'tier': f.tier}
|
||
for f in model.flashes
|
||
],
|
||
'final_state': {
|
||
'f_crust_hz': model.f_crust,
|
||
'f_sf_hz': model.f_sf,
|
||
'L_total': model.L_total,
|
||
'n_pinned': sum(1 for v in model.vortices if v.pinned),
|
||
'n_flashes': len(model.flashes)
|
||
}
|
||
}
|
||
|
||
out_path = os.path.join(out_dir, f"marble_jar_data_{datetime.now().strftime('%Y%m%d_%H%M%S')}.json")
|
||
with open(out_path, 'w') as f:
|
||
json.dump(data, f, indent=2)
|
||
print(f"Exported: {out_path}")
|
||
return out_path
|
||
|
||
|
||
# ─── Main ────────────────────────────────────────────────────────────────────
|
||
if __name__ == '__main__':
|
||
print("=" * 70)
|
||
print("PULSAR MARBLE-JAR MODEL (Multiscale)")
|
||
print("Cruise: years | Glitch: seconds | Recovery: days")
|
||
print("=" * 70)
|
||
|
||
model = PulsarMarbleJar(n_vortices=512, f0_hz=10.0)
|
||
model.run(t_max_yr=2e6)
|
||
|
||
print("\nGenerating visualizations...")
|
||
paths = visualize(model)
|
||
|
||
print("\nExporting data...")
|
||
export_data(model)
|
||
|
||
print("\n" + "=" * 70)
|
||
print("ALL DONE")
|
||
print("=" * 70)
|