Tuesday, August 19, 2025

Physics Simulator for finding EOS factors in Nebulae and Cosmic Strings

My love to ChatGPT and Google Gemini for the collaborative response efforts toward this astronomical exploration project!


import numpy as np
import matplotlib.pyplot as plt
import pubchempy as pcp
import requests
import json
import time
import re

# --- Scientific Constants ---
# Nuclear Physics Constants
CURIE_TO_BQ = 3.7e10           # Conversion factor: 1 Curie = 3.7 x 10^10 Becquerels (disintegrations per second)
LN2 = np.log(2)                # Natural logarithm of 2, used in half-life calculations

# Blackbody Radiation Constants
STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8 # Stefan-Boltzmann constant in W/(m^2·K^4)
PLANCK_CONSTANT = 6.62607015e-34           # Planck's constant in J·s
BOLTZMANN_CONSTANT = 1.380649e-23          # Boltzmann constant in J/K
SPEED_OF_LIGHT = 2.99792458e8              # Speed of light in a vacuum in m/s
WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3 # Wien's displacement law constant in m·K

# General Relativity Constants
GRAVITATIONAL_CONSTANT = 6.67430e-11 # G in m³·kg⁻¹·s⁻²
SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2

# --- Manual Nuclear Data for Gamma Emission ---
MANUAL_GAMMA_DATA = {
    "Cesium-137": {
        "half_life_days": 30.08 * 365.25,
        "gamma_energies_MeV": [0.662],
        "gamma_yields_per_decay": [0.851]
    },
    "Cobalt-60": {
        "half_life_days": 1925.21,
        "gamma_energies_MeV": [1.173, 1.332],
        "gamma_yields_per_decay": [0.999, 0.999]
    },
    "Iodine-131": {
        "half_life_days": 8.02,
        "gamma_energies_MeV": [0.364],
        "gamma_yields_per_decay": [0.817]
    },
    "Potassium-40": {
        "half_life_days": 1.251e9 * 365.25,
        "gamma_energies_MeV": [1.4608],
        "gamma_yields_per_decay": [0.1067]
    },
}

# ---------------------------------------------------------------------------
# PubChem Helper Functions
def _get_compound_cid_pubchem(compound_name):
    try:
        compounds = pcp.get_compounds(compound_name, 'name')
        return compounds[0].cid if compounds else None
    except Exception:
        return None

def _fetch_half_life_from_pubchem(cid):
    base_url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/"
    url = f"{base_url}cid/{cid}/JSON/"
    try:
        response = requests.get(url, timeout=5)
        response.raise_for_status()
        data = response.json()
        search_sections = ['Atomic Mass, Half Life, and Decay', 'Chemical and Physical Properties', 'Depositor-Supplied']
        for record in data.get('Record', {}).get('Section', []):
            if record.get('TOCHeading') in search_sections:
                for prop_section in record.get('Section', []):
                    if 'half-life' in prop_section.get('TOCHeading', '').lower():
                        for item in prop_section.get('Information', []):
                            for value_dict in item.get('Value', {}).get('StringWithMarkup', []):
                                text = value_dict.get('String')
                                if text:
                                    match = re.match(r"(\d+\.?\d*(?:e[+-]?\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)", text.lower())
                                    if match:
                                        return float(match.group(1)), match.group(2)
        return None, None
    except Exception:
        return None, None

def _convert_half_life_to_days(value, unit):
    if unit in ['s', 'seconds']: return value / (24 * 3600)
    if unit in ['min', 'minutes']: return value / (24 * 60)
    if unit in ['hr', 'hours']: return value / 24
    if unit in ['d', 'days']: return value
    if unit in ['y', 'years']: return value * 365.25
    return None

# ---------------------------------------------------------------------------
# Core Radiation Calculation Functions
def calculate_decay_constant(half_life_days):
    if half_life_days is None or half_life_days <= 0:
        return 0
    half_life_seconds = half_life_days * 24 * 3600
    return LN2 / half_life_seconds

def calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds):
    if decay_constant == 0:
        return initial_activity_Bq
    return initial_activity_Bq * np.exp(-decay_constant * time_elapsed_seconds)

def calculate_blackbody_spectral_radiance(wavelength_m, temperature_K):
    if temperature_K <= 0 or wavelength_m <= 0:
        return 0.0
    term1 = (2 * PLANCK_CONSTANT * SPEED_OF_LIGHT_SQ) / (wavelength_m**5)
    term2_exp_arg = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (wavelength_m * BOLTZMANN_CONSTANT * temperature_K)
    if term2_exp_arg > 700:
        return 0.0
    spectral_radiance = term1 / (np.exp(term2_exp_arg) - 1)
    return spectral_radiance

# ---------------------------------------------------------------------------
# General Relativity Functions
def calculate_schwarzschild_radius(mass_kg):
    """Calculates the Schwarzschild radius for a given mass."""
    if mass_kg <= 0:
        return 0.0
    return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def calculate_dilation_factor(mass_kg, distance_m):
    """
    Calculates the time dilation factor for an observer at infinity relative to an object in a gravitational field.
    A factor of 1 means no dilation. Values < 1 mean time is dilated (slowed down).
    """
    if mass_kg <= 0 or distance_m <= 0:
        return 1.0 # No dilation for zero mass or distance
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= schwarzschild_radius:
        return 0.0 # At or inside event horizon, time stops for an outside observer
    return np.sqrt(1 - (schwarzschild_radius / distance_m))

def calculate_relativistic_activity(isotope_name, initial_activity_Bq, time_elapsed_days, mass_kg, distance_m):
    """
    Calculates activity as observed by a distant observer, accounting for gravitational time dilation
    and a hypothetical "time crystal" effect.
    """
    half_life_days = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("half_life_days")
    if half_life_days is None:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None

    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    # Apply gravitational time dilation
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    # An observer far away sees time running slower for the object, so the decay process appears slowed.
    effective_time_local_seconds = time_elapsed_seconds * dilation_factor

    # --- Hypothetical: curvature-based modifier ---
    curvature_parameter = calculate_gravitational_curvature_parameter(mass_kg, distance_m)
    time_crystal_modifier = _calculate_time_crystal_half_life_modifier(curvature_parameter)

    effective_half_life_days = half_life_days * time_crystal_modifier
    decay_constant = calculate_decay_constant(effective_half_life_days)
    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, effective_time_local_seconds)

    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])
    total_gamma_yield_per_decay = sum(gamma_yields)
    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def calculate_relativistic_peak_wavelength(local_temperature_K, mass_kg, distance_m):
    """Calculates the observed peak wavelength of a blackbody, accounting for gravitational redshift."""
    if local_temperature_K <= 0:
        return None
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    observed_temperature_K = local_temperature_K * dilation_factor
    if observed_temperature_K <= 0:
        return np.inf
    peak_wavelength_m = WIEN_DISPLACEMENT_CONSTANT / observed_temperature_K
    return peak_wavelength_m * 1e9 # Convert to nm

# ---------------------------------------------------------------------------
# Quantum and Gravitational Bridge Functions (hypothetical)
def calculate_gravitational_curvature_parameter(mass_kg, distance_m):
    """
    A simplified, unitless parameter representing local gravitational curvature.
    Higher values imply greater curvature.
    """
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= 0 or distance_m <= schwarzschild_radius:
        return np.inf
    return schwarzschild_radius / distance_m

def quantum_differentiation_hypothetical(curvature_parameter, base_quantum_value):
    """
    Hypothetical function to illustrate 'quantum mechanical differentiation'.
    NOT based on established physics.
    """
    if curvature_parameter >= 1:
        return 0
    differentiated_value = base_quantum_value * (1 - curvature_parameter)
    return differentiated_value

def _calculate_time_crystal_half_life_modifier(curvature_parameter):
    """
    Hypothetical function to model a 'time crystal' effect driven by curvature.
    Returns a multiplicative modifier on the half-life.
    """
    if curvature_parameter >= 1:
        return 0.0
    modifier = np.tanh(1 - curvature_parameter)
    return modifier

# ---------------------------------------------------------------------------
# EOS & Stability (Nebulas + Cosmic Strings)
PROTON_MASS = 1.67262192369e-27
RADIATION_CONSTANT = 4 * STEFAN_BOLTZMANN_CONSTANT / SPEED_OF_LIGHT  # a_rad
K_B = BOLTZMANN_CONSTANT

# ---------- Gas / Nebula EOS ----------
def sound_speed_ideal_gas(T_K, mu=2.33, gamma=5/3):
    """
    Isentropic sound speed for an ideal gas:
      c_s = sqrt(gamma * k_B * T / (mu * m_p))
    mu ~ 2.33 for molecular clouds (H2 + He).
    """
    return np.sqrt(gamma * K_B * T_K / (mu * PROTON_MASS))

def pressure_ideal_gas(rho, T_K, mu=2.33):
    return rho * K_B * T_K / (mu * PROTON_MASS)

def pressure_polytrope(rho, K, gamma):
    return K * rho**gamma

def sound_speed_polytrope(rho, K, gamma):
    # c_s^2 = dP/dρ = gamma * K * ρ^(gamma-1)
    return np.sqrt(max(0.0, gamma * K * rho**(gamma - 1)))

def pressure_radiation(T_K):
    # P_rad = (1/3) a T^4
    return (1.0/3.0) * RADIATION_CONSTANT * T_K**4

def mixed_barotropic_pressure(rho, T_K, mu=2.33, gamma=5/3):
    # simple additive gas + radiation pressure
    return pressure_ideal_gas(rho, T_K, mu) + pressure_radiation(T_K)

# ---------- Nebula stability ----------
def jeans_length(rho, T_K, mu=2.33, gamma=5/3):
    """
    λ_J = c_s * sqrt(pi / (G * ρ))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma)
    return cs * np.sqrt(np.pi / (GRAVITATIONAL_CONSTANT * rho))

def jeans_mass(rho, T_K, mu=2.33, gamma=5/3):
    """
    M_J ~ (4π/3) * (λ_J/2)^3 * ρ
    """
    lam = jeans_length(rho, T_K, mu, gamma)
    return (4*np.pi/3.0) * (lam/2.0)**3 * rho

def bonnor_ebert_mass(T_K, P_ext, mu=2.33):
    """
    Critical Bonnor–Ebert mass for isothermal sphere:
      M_BE ≈ 1.18 * (c_s^4) / (G^(3/2) * P_ext^(1/2))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma=1.0)  # isothermal -> gamma=1
    return 1.18 * cs**4 / (GRAVITATIONAL_CONSTANT**1.5 * np.sqrt(P_ext))

# ---------- Cosmic strings (coarse-grained) ----------
def string_wiggliness_params(alpha=0.0, mu_base=1.0):
    """
    Phenomenology:
      U = μ (1+α),   T = μ /(1+α)
    α >= 0 encodes small-scale structure ("wiggles").
    """
    U = mu_base * (1.0 + alpha)
    T = mu_base / (1.0 + alpha)
    return U, T

def string_wave_speeds(alpha=0.0):
    """
    For the model above:
      c_T^2 = T/U = 1/(1+α)^2
    Use c_L ~ 1 (Nambu–Goto) as a simple default.
    """
    c_T = 1.0 / (1.0 + alpha)
    c_L = 1.0  # simple Nambu–Goto proxy
    return c_T, c_L

def string_effective_w(v_rms=0.5):
    """
    Effective EOS parameter for a network (toy model):
      w ≈ -1/3 * (1 - v_rms^2)
    """
    v2 = np.clip(v_rms**2, 0.0, 1.0)
    return -1.0/3.0 * (1.0 - v2)

def string_loop_oscillation_time(loop_length_m, alpha=0.0):
    """
    Fundamental oscillation time τ ~ L / (c_T * c).
    """
    c_T, _ = string_wave_speeds(alpha)
    return loop_length_m / (c_T * SPEED_OF_LIGHT)

# ---------- Visualizers ----------
def plot_nebula_stability_grid(rho_vals, T_vals, mu=2.33, gamma=5/3, mass_scale_Msun=1.0):
    """
    Heatmap: M_J / (mass_scale) to see where collapse is favored (M_J smaller than cloud mass).
    """
    from matplotlib.colors import LogNorm
    R, T = np.meshgrid(rho_vals, T_vals)
    MJ = np.vectorize(jeans_mass)(R, T, mu, gamma)
    plt.figure(figsize=(9,6))
    im = plt.pcolormesh(R, T, MJ/(mass_scale_Msun*1.98847e30), norm=LogNorm())
    plt.colorbar(im, label=r"$M_J/M_\odot$")
    plt.xscale('log'); plt.yscale('log')
    plt.xlabel(r"Density $\rho$ (kg m$^{-3}$)")
    plt.ylabel("Temperature (K)")
    plt.title("Nebula Stability Map (Jeans Mass in Solar Masses)")
    plt.tight_layout(); plt.show()

def plot_string_speeds_vs_wiggliness(alpha_vals):
    cT = [string_wave_speeds(a)[0] for a in alpha_vals]
    plt.figure(figsize=(8,5))
    plt.plot(alpha_vals, cT)
    plt.xlabel("Wiggliness α")
    plt.ylabel("Transverse wave speed c_T (units of c)")
    plt.title("Cosmic String Wave Speed vs Wiggliness")
    plt.grid(True); plt.tight_layout(); plt.show()

# ---------------------------------------------------------------------------
# Analysis Functions
def generate_normal_distribution_series(mean, std_dev, num_points=200, spread_factor=3):
    start = mean - spread_factor * std_dev
    end = mean + spread_factor * std_dev
    if start < 0:
        start = 0
    series_values = np.linspace(start, end, num_points)
    pdf_values = (1 / (std_dev * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((series_values - mean) / std_dev)**2)
    return series_values, pdf_values

def get_gamma_curies_single_point(isotope_name, initial_activity_Bq, time_elapsed_days):
    half_life_days = None
    if isotope_name in MANUAL_GAMMA_DATA:
        half_life_days = MANUAL_GAMMA_DATA[isotope_name]["half_life_days"]
    else:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None
    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])
    total_gamma_yield_per_decay = sum(gamma_yields)
    decay_constant = calculate_decay_constant(half_life_days)
    time_elapsed_seconds = time_elapsed_days * 24 * 3600
    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds)
    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def analyze_gamma_curies_comparison(isotope_name, initial_activity_Bq, mean_time_days, std_dev_time_days, mass_kg, distance_m, num_points=200):
    """
    Compares gamma curies decay in flat spacetime vs. near a massive k-body.
    """
    print(f"\n--- Comparing Gamma Curies: Flat Spacetime vs. Relativistic Environment ---")
    print(f"  Isotope: {isotope_name}, Initial Activity: {initial_activity_Bq:.2e} Bq")
    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")
    if dilation_factor < 0.99:
        print(f"  (A dilation factor less than 1 means time runs slower in the gravitational field.)")
    else:
        print(f"  (Dilation factor is near 1, indicating a weak gravitational effect.)")

    time_series_days, _ = generate_normal_distribution_series(mean_time_days, std_dev_time_days, num_points=num_points)

    # Flat Spacetime Calculation
    flat_curies_series = [get_gamma_curies_single_point(isotope_name, initial_activity_Bq, t_days) for t_days in time_series_days]

    # Relativistic Calculation
    relativistic_curies_series = [calculate_relativistic_activity(isotope_name, initial_activity_Bq, t_days, mass_kg, distance_m) for t_days in time_series_days]

    if all(c is None for c in flat_curies_series):
        print(f"  Calculation failed. Please ensure '{isotope_name}' has half-life data.")
        return

    # Plotting
    plt.figure(figsize=(12, 7))
    plt.plot(time_series_days, flat_curies_series, label='Flat Spacetime (Local Observer)')
    plt.plot(time_series_days, relativistic_curies_series, label=f'Relativistic (Distant Observer) w/ Curvature Modifier', linestyle='--')
    plt.title(f'Gamma Curies Decay: {isotope_name} - Flat vs. Relativistic')
    plt.xlabel('Time (days)')
    plt.ylabel('Gamma Curies (Ci)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def simulate_blackbody_wien_comparison(temperatures_K, mass_kg, distance_m, wavelength_range_nm=(1, 10000)):
    """
    Compares blackbody peak wavelength in flat spacetime vs. near a massive k-body.
    """
    print(f"\n--- Comparing Blackbody Peak Wavelength: Flat vs. Relativistic ---")
    print(f"  Analyzing temperatures: {temperatures_K} K")
    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")

    temp_list_sorted = sorted(temperatures_K)
    flat_peak_wl_list_nm = [(WIEN_DISPLACEMENT_CONSTANT / T) * 1e9 for T in temp_list_sorted if T > 0]
    relativistic_peak_wl_list_nm = [calculate_relativistic_peak_wavelength(T, mass_kg, distance_m) for T in temp_list_sorted if T > 0]

    # Plotting
    plt.figure(figsize=(8, 6))
    plt.plot(temp_list_sorted, flat_peak_wl_list_nm, 'o-', label='Flat Spacetime (Local Observer)')
    plt.plot(temp_list_sorted, relativistic_peak_wl_list_nm, 'r--', label='Relativistic (Distant Observer)')
    plt.title('Wien\'s Law: Peak Wavelength vs. Temperature (Flat vs. Relativistic)')
    plt.xlabel('Temperature (K)')
    plt.ylabel('Peak Wavelength (nm)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def explore_quantum_gravitational_differentiation(mass_kg, distance_m_series):
    """
    New function to explore the hypothetical relationship between gravity and a quantum state.
    """
    print("\n--- Exploring a Hypothetical Gravitational-Quantum Bridge ---")
    print(f"  Analyzing Quantum State near a k-body with Mass={mass_kg:.2e} kg")
    curvature_series = [calculate_gravitational_curvature_parameter(mass_kg, d) for d in distance_m_series]
    base_quantum_value = 1.0
    differentiated_value_series = [quantum_differentiation_hypothetical(c, base_quantum_value) for c in curvature_series]

    plt.figure(figsize=(10, 6))
    plt.plot(distance_m_series, differentiated_value_series, 'b-')
    plt.axvline(x=calculate_schwarzschild_radius(mass_kg), color='r', linestyle='--', label='Schwarzschild Radius')
    plt.title('Hypothetical "Quantum Differentiation" vs. Distance from a k-body')
    plt.xlabel('Distance from k-body (m)')
    plt.ylabel('Hypothetical Differentiated Quantum Value (Unitless)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

# ---------------------------------------------------------------------------
# EOS Stability Explorer UI (Nebulas + Cosmic Strings)
def plot_nebula_stability_driver():
    try:
        m_cloud_Msun = float(input("  Cloud mass to compare (M_sun): "))
        mu = float(input("  Mean molecular weight mu (e.g., 2.33): "))
        gamma = float(input("  Adiabatic index gamma (e.g., 1.4 or 1.6667): "))
    except ValueError:
        print("  Invalid input.")
        return
    # Typical molecular cloud conditions in SI (kg/m^3) and K
    rho_vals = np.logspace(-22, -16, 120)  # kg/m^3
    T_vals = np.logspace(1, 3.5, 120)      # 10 K – ~3000 K
    plot_nebula_stability_grid(rho_vals, T_vals, mu=mu, gamma=gamma, mass_scale_Msun=m_cloud_Msun)
    print("  Darker regions (lower M_J/M_sun) imply easier collapse; where M_cloud > M_J, fragmentation likely.")

def plot_cosmic_string_driver():
    try:
        alpha_min = float(input("  Wiggliness α min (e.g., 0): "))
        alpha_max = float(input("  Wiggliness α max (e.g., 3): "))
    except ValueError:
        print("  Invalid input.")
        return
    alphas = np.linspace(alpha_min, alpha_max, 200)
    plot_string_speeds_vs_wiggliness(alphas)
    # quick demo of loop time
    try:
        L = float(input("  Example loop length L (meters): "))
        a_demo = float(input("  Example α (for τ): "))
    except ValueError:
        print("  Invalid input.")
        return
    tau = string_loop_oscillation_time(L, a_demo)
    print(f"  Loop oscillation time τ ≈ {tau:.3e} s at α={a_demo}. Lower c_T (higher α) → longer τ.")

# ---------------------------------------------------------------------------
# Main Program Execution
if __name__ == "__main__":
    print("--- Radiation + EOS Stability Explorer ---")
    print("\nThis program explores:")
    print("  1) Radioactive decay (flat vs relativistic, with a hypothetical curvature modifier).")
    print("  2) Blackbody redshift (flat vs relativistic).")
    print("  3) Hypothetical quantum–gravity differentiation.")
    print("  4) EOS Stability Explorer: Nebulas (Jeans/Bonnor–Ebert context) and Cosmic Strings (wiggliness & wave speeds).")

    while True:
        main_choice = input("Select: (R)elativistic Comparison, (Q)uantum-Gravity, (S)tability Explorer, (E)xit? ").strip().lower()
        if main_choice == 'e':
            break

        elif main_choice == 'r':
            sub_choice = input("Select Relativistic Scenario: (N)uclear Decay or (B)lackbody Radiation? ").strip().lower()
            try:
                # Get common relativistic parameters
                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a star) in kg: "))
                k_body_distance = float(input("  Enter distance from the center of the k-body in meters: "))
            except ValueError:
                print("  Invalid input. Please enter numerical values for mass and distance.")
                continue

            if sub_choice == 'n':
                print("\n--- Relativistic Nuclear Decay Comparison ---")
                isotope_name = input("  Enter radioactive isotope name (e.g., Cesium-137): ").strip()
                if isotope_name not in MANUAL_GAMMA_DATA and _get_compound_cid_pubchem(isotope_name) is None:
                    print(f"  '{isotope_name}' data not found.")
                    continue
                try:
                    initial_activity = float(input("  Enter initial activity in Bq: "))
                    mean_time = float(input("  Enter mean measurement time (days): "))
                    std_dev_time = float(input("  Enter standard deviation for time (days): "))
                except ValueError:
                    print("  Invalid input.")
                    continue
                analyze_gamma_curies_comparison(isotope_name, initial_activity, mean_time, std_dev_time, k_body_mass, k_body_distance)

            elif sub_choice == 'b':
                print("\n--- Relativistic Blackbody Radiation Comparison ---")
                temp_input_str = input("  Enter temperatures (K) separated by commas (e.g., 300, 1000, 5000): ")
                try:
                    temperatures_list = [float(t.strip()) for t in temp_input_str.split(',')]
                except ValueError:
                    print("  Invalid temperature input.")
                    continue
                simulate_blackbody_wien_comparison(temperatures_list, k_body_mass, k_body_distance)
            else:
                print("  Invalid choice. Please enter 'N' or 'B'.")

        elif main_choice == 'q':
            print("\n--- Quantum-Gravitational Bridge Exploration ---")
            try:
                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a black hole) in kg: "))
                max_distance_factor = float(input("  Enter max distance as a multiple of Schwarzschild Radius (e.g., 5): "))
                schwarzschild_radius_val = calculate_schwarzschild_radius(k_body_mass)
                if schwarzschild_radius_val == 0:
                    print("  Invalid mass. Cannot calculate Schwarzschild Radius.")
                    continue
                distance_series = np.linspace(schwarzschild_radius_val + 1e-9, max_distance_factor * schwarzschild_radius_val, 200)
                explore_quantum_gravitational_differentiation(k_body_mass, distance_series)
            except ValueError:
                print("  Invalid input. Please enter numerical values.")
                continue

        elif main_choice == 's':
            print("\n--- EOS Stability Explorer ---")
            sub = input("Choose: (N)ebula map, (C)osmic-strings? ").strip().lower()
            if sub == 'n':
                plot_nebula_stability_driver()
            elif sub == 'c':
                plot_cosmic_string_driver()
            else:
                print("  Invalid choice.")

        else:
            print("Invalid choice. Please enter 'R', 'Q', 'S', or 'E'.")

    print("\nExiting program. Thank you for exploring these concepts!")

No comments: