cometsfactory 0.0.3

Comet factory — classify, build and catalogue comets of any type: short-period, long-period, Halley-type, sungrazer, interstellar, main-belt comet, centaur-transition, and extinct.
Documentation
use crate::config::parameters::{AU, PI};

/// Dust tail orientation angle (rad) — the angle between the
/// extended Sun-comet radius vector and the tail direction.
/// For a syndyne with β, aberration angle ≈ β × a_rad × t² / (2 v r).
/// Simplified: θ_tail ≈ atan(v_orbit / v_rad_pressure).
pub fn dust_tail_angle(beta: f64, r_au: f64) -> f64 {
    let radial_accel = beta * 5.93e-3 / (r_au * r_au);
    let v_orbit = 29784.0 / r_au.sqrt();
    if v_orbit < 1.0 {
        return 0.0;
    }
    (radial_accel / v_orbit).atan()
}

/// Ion tail aberration angle (rad) due to solar wind speed and orbital velocity.
/// α = atan(v_orbit / v_sw)
pub fn ion_tail_aberration(r_au: f64, solar_wind_velocity: f64) -> f64 {
    let v_orbit = 29784.0 / r_au.max(0.01).sqrt();
    (v_orbit / solar_wind_velocity.max(1.0)).atan()
}

/// Striae spacing estimate (m) in dust tail from synchrone analysis.
/// Δs ≈ v_ejection × Δt
pub fn striae_spacing(ejection_velocity: f64, time_interval_s: f64) -> f64 {
    ejection_velocity * time_interval_s
}

/// Apparent tail length (degrees) as seen from Earth.
/// θ = atan(L_physical / (Δ × AU))
pub fn apparent_tail_length_deg(physical_length_m: f64, delta_au: f64) -> f64 {
    let d = delta_au * AU;
    if d < 1.0 {
        return 0.0;
    }
    (physical_length_m / d).atan().to_degrees()
}

/// Maximum visible tail surface brightness relative to sky background.
/// Simplified model: ∝ Afρ / (r² Δ²).
pub fn tail_surface_brightness(afrho_cm: f64, r_au: f64, delta_au: f64) -> f64 {
    let denom = r_au * r_au * delta_au * delta_au;
    if denom < 1.0e-10 {
        return 0.0;
    }
    afrho_cm / denom
}

/// Disconnection event probability: ion tail can disconnect during
/// solar wind sector boundary crossings. Returns a 0–1 heuristic.
pub fn disconnection_probability(solar_wind_velocity: f64, gas_production_q: f64) -> f64 {
    let sw_factor: f64 = if solar_wind_velocity > 600_000.0 {
        0.5
    } else {
        0.1
    };
    let q_factor: f64 = if gas_production_q > 1.0e28 { 0.3 } else { 0.1 };
    (sw_factor + q_factor).min(1.0)
}

/// Anti-tail visibility: whether conditions favor an anti-tail
/// (Earth near comet's orbital plane, past perihelion, large dust particles).
/// Returns the projected anti-tail angle (rad).
pub fn anti_tail_angle(comet_inclination_deg: f64, earth_lat_from_plane_deg: f64) -> f64 {
    let i_rad = comet_inclination_deg.to_radians();
    let lat_rad = earth_lat_from_plane_deg.to_radians();
    let projected = i_rad.sin() * lat_rad.cos();
    projected.abs().asin()
}

/// Angular jet feature extent (rad) from nucleus rotation.
/// θ_jet ≈ atan(v_jet / v_expansion)
pub fn jet_feature_extent(jet_velocity: f64, expansion_velocity: f64) -> f64 {
    if expansion_velocity < 1.0 {
        return PI / 2.0;
    }
    (jet_velocity / expansion_velocity).atan()
}