Struct DiskStabilityIC 

pub struct DiskStabilityIC {
    pub disk_surface_density: Box<dyn Fn(f64) -> f64 + Send + Sync>,
    pub disk_velocity_dispersion: Box<dyn Fn(f64) -> f64 + Send + Sync>,
    pub bulge: Option<Box<dyn IsolatedEquilibrium>>,
    pub halo_potential: Option<Box<dyn Fn([f64; 3]) -> f64 + Send + Sync>>,
    pub perturbation_mode_m: u32,
    pub perturbation_pattern_speed: f64,
    pub perturbation_amplitude: f64,
}

Disk stability IC: f(E, Lz) for axisymmetric disk plus an optional perturbation.

Fields

disk_surface_density: Box<dyn Fn(f64) -> f64 + Send + Sync>

Disk surface density Σ(R) as function of cylindrical radius.

disk_velocity_dispersion: Box<dyn Fn(f64) -> f64 + Send + Sync>

Radial velocity dispersion σ_R(R).

bulge: Option<Box<dyn IsolatedEquilibrium>>

Optional central bulge component.

halo_potential: Option<Box<dyn Fn([f64; 3]) -> f64 + Send + Sync>>

Optional fixed dark matter halo potential Φ_halo(x).

perturbation_mode_m: u32

Azimuthal mode number m (m=2 = bar, m=3 = triangle).

perturbation_pattern_speed: f64

Pattern speed Ω_p in rad/time.

perturbation_amplitude: f64

Relative amplitude δΣ/Σ.

Implementations

impl DiskStabilityIC

pub fn new( disk_surface_density: Box<dyn Fn(f64) -> f64 + Send + Sync>, disk_velocity_dispersion: Box<dyn Fn(f64) -> f64 + Send + Sync>, perturbation_mode_m: u32, perturbation_pattern_speed: f64, perturbation_amplitude: f64, ) -> Self

pub fn toomre_q(&self, radius: f64) -> f64

Toomre Q(R) = σ_R κ / (3.36 G Σ). Q > 1 means locally stable.

Computes the epicyclic frequency κ(R) = √(R dΩ²/dR + 4Ω²) from the combined disk+bulge+halo circular velocity curve.

pub fn sample_on_grid(&self, domain: &Domain) -> PhaseSpaceSnapshot

Sample onto 6D grid: construct f(E, Lz) for disk using the Shu (1969) distribution function, then superpose an azimuthal perturbation mode.

Shu DF: f(E, L_z) = [Σ(R_c) Ω(R_c)] / [π κ(R_c) σ_R²(R_c)] × exp[-(E - E_c(L_z)) / σ_R²(R_c)]

for L_z > 0 (prograde orbits). f = 0 for retrograde (L_z ≤ 0).

The vertical structure uses a sech²(z/z_0) profile with z_0 = σ_z / √(4πGΣ), assuming σ_z ≈ σ_R (isotropic approximation).