Tuesday, August 19, 2025

Advanced Physics Simulator with Time Crystal Model

  Enjoy this Time-Crystal exploration program.

 
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 (for the new functions)
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 (Unchanged) ---
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 (Modified to accept dilation factor) ---
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

# --- UPDATED: 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

    # --- NEW: Apply Time Crystal "dynamic range modifier" ---
    # We use a hypothetical model where the half-life itself is affected by extreme gravity.
    # The `curvature_parameter` serves as a proxy for the 'valence instability' and the
    # potential for the system to "collimate" into a time crystal state.
    curvature_parameter = calculate_gravitational_curvature_parameter(mass_kg, distance_m)
    time_crystal_modifier = _calculate_time_crystal_half_life_modifier(curvature_parameter)
   
    # We multiply the original half-life by this modifier. A modifier < 1 means the
    # "effective" half-life shortens, speeding up the perceived decay rate, which
    # could model the collapse into a new, more stable state.
    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

    # Gravitational redshift means a distant observer sees the blackbody as cooler.
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    observed_temperature_K = local_temperature_K * dilation_factor

    # Use Wien's law with the observed temperature
    if observed_temperature_K <= 0:
        return np.inf # Effectively no peak wavelength observed
   
    peak_wavelength_m = WIEN_DISPLACEMENT_CONSTANT / observed_temperature_K
    return peak_wavelength_m * 1e9 # Convert to nm

# --- NEW: Quantum and Gravitational Bridge Functions ---
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'.
    This is NOT based on established physics but serves as an illustrative model.
    It proposes that a quantum value (e.g., a spin state) changes with gravitational curvature.
   
    For this example, we'll model it as a simple linear relationship for visualization.
    """
    if curvature_parameter >= 1:
        return 0 # At or inside the event horizon, the effect is undefined or at its maximum
   
    # A hypothetical effect: the base value is 'differentiated' by the curvature.
    # We use (1 - curvature_parameter) to get a value that goes from 1 to 0 as curvature increases.
    differentiated_value = base_quantum_value * (1 - curvature_parameter)
    return differentiated_value

def _calculate_time_crystal_half_life_modifier(curvature_parameter):
    """
    Hypothetical function to model the "time crystal" effect.
    As gravitational curvature increases, the half-life is modified, simulating a
    "collapsing" or "colliding" into a new, stable-like state.
   
    This is not based on established physics.
    The function returns a modifier that, when multiplied by the half-life,
    changes its value. A value < 1 shortens the half-life.
    """
    if curvature_parameter >= 1:
        # At or inside the event horizon, the half-life effectively becomes zero.
        # This models the ultimate collapse into a "time crystal" state.
        return 0.0

    # We use a simple, non-linear function to model the effect.
    # As curvature approaches 1, the modifier approaches zero.
    modifier = np.tanh(1 - curvature_parameter)
    return modifier

# --- Utility Functions (Unchanged) ---
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

# --- Analysis Functions (Modified to include new relativistic options) ---

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/ Time Crystal 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]
   
    # Let's use a base quantum value, for example, a spin state (up/down).
    # We'll use 1.0 as a base value to represent an "undifferentiated" quantum state.
    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()

# --- Main Program Execution (Modified) ---
if __name__ == "__main__":
    print("--- Radiation Analysis Program: Interpreting Data Behavior ---")
    print("\nThis program explores two physics concepts in both a standard, flat-spacetime environment and a relativistic one.")
   
    print("1. **Radioactive Decay:** Compares the decay rate of an isotope as seen by a local observer vs. a distant one in a gravitational field, including a hypothetical 'time crystal' effect.")
    print("2. **Blackbody Radiation:** Compares the peak emission wavelength of a blackbody as seen by a local observer vs. a distant one (gravitational redshift).")
    print("3. **Hypothetical Quantum-Gravity Bridge:** A new, conceptual section to explore how a quantum state might be 'differentiated' by a gravitational field.\n")

    while True:
        main_choice = input("Select: (R)elativistic Comparison, (Q)uantum-Gravity, (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: "))
                # We'll generate a range of distances to show the effect
                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

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

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

No comments: