Friday, May 30, 2025

Rakshas' law on nuclear automata and quantum relativity with boundary breaking Swarszchild maxima for cosmic string evaluation - 3 Versions! Updated for General Relativity/Quantum Mechanical exploration & Exceeding Schwarzschild radius

Radical! Nearest neighbour tracing in nuclear automata. Below the program is one that incorporates general relativity and now gravity time-dialation ratio/scalar curvature as a k-body factor due to certain chemical stoichastic reporting systems you may choose to use either, this new update also explores a hypothetical relationship between the gravitational curvature parameter and a quantum mechanical property, this might elucidate how such a "bridge" might be conceptually built, even if not physically validated :-)

Rakshas' law on nuclear automata and quantum relativity with boundary breaking Swarszchild maxima for cosmic string evaluation - 3 Versions! Updated for General Relativity/Quantum Mechanical exploration & Exceeding Schwarzschild radius

 

For the second program please search, I postulate 

And thirdly search Cosmic because it takes many lightyears of a cosmic string to expand past swarszchild radius? 

According to AI, NO!?

 It’s the mass‑per‑unit‑length you’d need for a straight line of matter to have a Schwarzschild radius equal to its own length. In other words, if you could somehow pack 6.7 × 10^26 kilograms into every single meter of a rod, then no matter how long or short it was, its gravitational radius would match or exceed its physical length. 


So no known element — even in its most extreme astrophysical form — comes anywhere close. This number is more in the realm of hypothetical cosmic strings or other exotic matter in theoretical physics.
For a cosmic string with that linear density, the spacetime would be so strongly curved that “velocity” measurements would need to be defined relative to the conical geometry around the string. Locally, you could still measure speeds with clocks and rulers, but the gravitational effects would dominate the dynamics.

Here's the short version while I work on a full report:
• Conical geometry typically arises around cosmic strings, where spacetime has a deficit angle — like a wedge removed from a flat sheet and the edges sewn together. This affects how geodesics (paths of particles and light) behave.


• Anyons, which exist in 2D systems and obey fractional statistics, can induce topological effects. In some models, they can be associated with curvature or torsion in spacetime, though this is speculative and mostly explored in quantum gravity or condensed matter analogs.


• Inverting conical geometry would mean adding an excess angle — like inserting a wedge instead of removing one. This could hypothetically create a spacetime with negative curvature around the defect.


• In such a geometry, velocity could still be locally measurable using proper time and distance, but its interpretation from a distant observer might be altered due to the modified geodesic structure and time dilation effects.


more to come....also take into account inverse time in this context, what can be inferred from known publication and experimentally confirmed through this program? 

Check the latest post :) 

 

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

# --- Manual Nuclear Data for Gamma Emission ---
# This dictionary stores known half-lives and gamma emission properties for isotopes.
# IMPORTANT: These values should be verified from reliable nuclear data sources (e.g., NNDC).
# PubChem provides general information, but specific gamma yields are best from dedicated databases.
MANUAL_GAMMA_DATA = {
    "Cesium-137": {
        "half_life_days": 30.08 * 365.25, # Half-life in days (approx 30.08 years)
        "gamma_energies_MeV": [0.662],
        "gamma_yields_per_decay": [0.851] # Gamma yield per decay (from Ba-137m isomer)
    },
    "Cobalt-60": {
        "half_life_days": 1925.21, # Half-life in days
        "gamma_energies_MeV": [1.173, 1.332],
        "gamma_yields_per_decay": [0.999, 0.999] # Approximately 1 gamma per decay at each energy
    },
    "Iodine-131": {
        "half_life_days": 8.02, # Half-life in days
        "gamma_energies_MeV": [0.364],
        "gamma_yields_per_decay": [0.817] # Gamma yield per decay
    },
    "Potassium-40": {
        "half_life_days": 1.251e9 * 365.25, # Half-life in days (~1.251 billion years)
        "gamma_energies_MeV": [1.4608],
        "gamma_yields_per_decay": [0.1067] # ~10.67% of K-40 decays via electron capture
    },
    # Feel free to add more isotopes here with their respective data
}

# --- PubChem Helper Functions (for half-life lookup if not in manual data) ---
# These functions attempt to fetch half-life data from PubChem,
# but their parsing can be fragile due to variations in PubChem's data structure.
def _get_compound_cid_pubchem(compound_name):
    """
    Searches PubChem for a compound's CID (Compound ID) by name.
    Returns the CID if found, otherwise None.
    """
    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):
    """
    Fetches half-life data from PubChem using a compound's CID.
    Attempts to parse numerical value and unit from the returned text.
    Returns (value, unit) tuple if successful, otherwise (None, None).
    """
    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() # Raise an exception for HTTP errors
        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:
                                    # Use regex to find numerical value and unit
                                    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):
    """
    Converts a half-life value from various units to days.
    Returns the half-life in days, or None if the unit is unrecognized.
    """
    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):
    """
    Calculates the radioactive decay constant (lambda) from the half-life.
    Returns the decay constant in inverse seconds (s^-1).
    If half_life_days is non-positive, returns 0 (implying a stable isotope).
    """
    if half_life_days is None or half_life_days <= 0:
        return 0
    half_life_seconds = half_life_days * 24 * 3600 # Convert days to seconds
    return LN2 / half_life_seconds

def calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds):
    """
    Calculates the activity (number of decays per second) of a radioactive source
    at a given time, using the decay law (N = N0 * e^(-lambda*t)).
    Returns activity in Becquerels (Bq).
    """
    if decay_constant == 0:
        return initial_activity_Bq # For stable isotopes or infinite half-life
    return initial_activity_Bq * np.exp(-decay_constant * time_elapsed_seconds)

def get_gamma_curies_single_point(isotope_name, initial_activity_Bq, time_elapsed_days):
    """
    Calculates the total gamma curies emitted by an isotope at a specific time.
    It first tries to get half-life from the `MANUAL_GAMMA_DATA`, then falls back to PubChem.
    Returns the gamma activity in Curies (Ci), or None if half-life data isn't found.
    """
    half_life_days = None
    if isotope_name in MANUAL_GAMMA_DATA:
        half_life_days = MANUAL_GAMMA_DATA[isotope_name]["half_life_days"]
    else:
        # If not in manual data, try PubChem
        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 # Could not determine half-life

    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 calculate_blackbody_spectral_radiance(wavelength_m, temperature_K):
    """
    Calculates the spectral radiance of a blackbody at a given wavelength and temperature,
    using Planck's Law.
    Returns spectral radiance in Watts per steradian per square meter per meter of wavelength (W/(m³·sr)).
    Returns 0.0 if temperature or wavelength are non-positive or if the exponent is too large.
    """
    if temperature_K <= 0 or wavelength_m <= 0:
        return 0.0

    term1 = (2 * PLANCK_CONSTANT * SPEED_OF_LIGHT**2) / (wavelength_m**5)
    term2_exp_arg = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (wavelength_m * BOLTZMANN_CONSTANT * temperature_K)
    
    # Check to prevent numerical overflow for very large exponent arguments
    if term2_exp_arg > 700:
        return 0.0 # Effectively zero emission at this wavelength/temperature
        
    spectral_radiance = term1 / (np.exp(term2_exp_arg) - 1)
    return spectral_radiance

# --- Utility Function for Generating Bell Curve Data ---
def generate_normal_distribution_series(mean, std_dev, num_points=200, spread_factor=3):
    """
    Generates a series of values that follow a normal (bell curve) distribution.
    This is used to create a distribution for a parameter, like measurement time.
    
    Args:
        mean (float): The mean (average) of the distribution.
        std_dev (float): The standard deviation of the distribution.
        num_points (int): The number of data points to generate.
        spread_factor (int): How many standard deviations away from the mean
                              the series should span (e.g., 3 for 3-sigma range).
    
    Returns:
        tuple: A tuple containing:
            - series_values (numpy.array): The generated data points.
            - pdf_values (numpy.array): The probability density function (PDF) values
                                        corresponding to `series_values`.
    """
    start = mean - spread_factor * std_dev
    end = mean + spread_factor * std_dev
    if start < 0:
        start = 0 # Ensure time or other physical quantities are non-negative
    
    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

# --- Analysis Functions for "Nearest Neighbor Boundaries" and "Displacement" ---

def analyze_gamma_curies_nearest_neighbor_boundaries(isotope_name, initial_activity_Bq, 
                                                     mean_time_days, std_dev_time_days, 
                                                     num_points=200, derivative_threshold_abs=None):
    """
    Calculates gamma curies over a series of measurement times sampled from a bell curve.
    It identifies "nearest neighbor boundaries" by highlighting regions where the
    absolute rate of change (derivative) of gamma curies is above a given threshold.
    This represents where the decay process is changing most rapidly.
    
    Args:
        isotope_name (str): The name of the radioactive isotope (e.g., "Cesium-137").
        initial_activity_Bq (float): The initial activity of the isotope in Becquerels (Bq).
        mean_time_days (float): The mean (average) of the bell curve distribution for measurement times in days.
        std_dev_time_days (float): The standard deviation of the bell curve for measurement times in days.
        num_points (int): The number of time points to simulate.
        derivative_threshold_abs (float, optional): The absolute threshold for the derivative (Ci/day).
                                                    If None, the program will auto-calculate a threshold
                                                    (e.g., 10% of the maximum absolute derivative).
    """
    print(f"\n--- Identifying 'Nearest Neighbor Boundaries' in Gamma Curies ---")
    print(f"  Isotope: {isotope_name}, Initial Activity: {initial_activity_Bq:.2e} Bq")
    print(f"  Time Distribution: Mean={mean_time_days:.2f} days, Std Dev={std_dev_time_days:.2f} days")
    print(f"  'Nearest Neighbor Boundaries' are regions where the decay rate changes significantly.")

    # 1. Generate time series based on a normal distribution
    time_series_days, pdf_values = generate_normal_distribution_series(
        mean_time_days, std_dev_time_days, num_points=num_points)
    
    # 2. Calculate gamma curies for each time point
    curies_series = []
    for t_days in time_series_days:
        curies = get_gamma_curies_single_point(isotope_name, initial_activity_Bq, t_days)
        curies_series.append(curies if curies is not None else np.nan) # Store None as NaN
    curies_series = np.array(curies_series) # Convert to NumPy array for numerical operations

    # Handle cases where half-life was not found (curies_series will be all NaN)
    if np.all(np.isnan(curies_series)):
        print(f"  Calculation failed. Please ensure '{isotope_name}' has half-life data.")
        return

    # 3. Calculate the numerical derivative (rate of change) of gamma curies with respect to time
    # np.gradient provides a robust way to estimate the derivative for discrete data.
    derivative_curies = np.gradient(curies_series, time_series_days)

    # 4. Determine or use a threshold for identifying "boundaries"
    if derivative_threshold_abs is None:
        # If no threshold is provided, auto-calculate it as 10% of the maximum absolute derivative
        # This highlights the top 10% steepest changes.
        max_abs_deriv = np.max(np.abs(derivative_curies[np.isfinite(derivative_curies)]))
        derivative_threshold_abs = max_abs_deriv * 0.1
        print(f"  Auto-calculated derivative threshold: {derivative_threshold_abs:.2e} Ci/day")
    else:
        print(f"  Using provided derivative threshold: {derivative_threshold_abs:.2e} Ci/day")

    # 5. Identify points where the absolute derivative exceeds the threshold
    boundary_indices = np.where(np.abs(derivative_curies) > derivative_threshold_abs)
    boundary_times = time_series_days[boundary_indices]
    boundary_curies = curies_series[boundary_indices]

    # --- Plotting Results ---
    plt.figure(figsize=(12, 7))
    plt.plot(time_series_days, curies_series, label='Gamma Curies Decay')
    plt.scatter(boundary_times, boundary_curies, color='red', s=50, 
                label='Nearest Neighbor Boundaries (High Rate of Change)', zorder=5)
    
    plt.title(f'Gamma Curies with "Nearest Neighbor Boundaries" for {isotope_name}')
    plt.xlabel('Time (days)')
    plt.ylabel('Gamma Curies (Ci)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

    # Optional: Plot the derivative itself for more insight
    plt.figure(figsize=(10, 5))
    plt.plot(time_series_days, derivative_curies, label='Rate of Change (dCuries/dt)')
    plt.axhline(y=derivative_threshold_abs, color='grey', linestyle='--', label='Positive Threshold')
    plt.axhline(y=-derivative_threshold_abs, color='grey', linestyle='--', label='Negative Threshold')
    plt.title('Rate of Change of Gamma Curies over Time')
    plt.xlabel('Time (days)')
    plt.ylabel('dCuries/dt (Ci/day)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

    print(f"\n  Identified {len(boundary_times)} data points in 'Nearest Neighbor Boundary' regions.")
    if len(boundary_times) > 0:
        print("  Example Boundary Points (Time, Curies, Derivative):")
        # Print a few examples from the identified boundary points
        for i in range(min(5, len(boundary_times))):
            idx = boundary_indices[0][i] # Get the original index from the series
            print(f"    Time: {time_series_days[idx]:.2f} days, Curies: {curies_series[idx]:.2e} Ci, Derivative: {derivative_curies[idx]:.2e} Ci/day")

def simulate_blackbody_wien_displacement(temperatures_K, wavelength_range_nm=(1, 10000)):
    """
    Simulates blackbody spectral radiance across a range of wavelengths for given temperatures,
    and demonstrates Wien's Displacement Law by showing how the peak emission wavelength
    ("displaces" or shifts) with temperature.
    
    Args:
        temperatures_K (list): A list of temperatures in Kelvin (K) to simulate.
        wavelength_range_nm (tuple): A tuple (min_wavelength, max_wavelength) in nanometers.
    """
    print(f"\n--- Simulating Blackbody Spectral Radiance & Wien's Displacement Law ---")
    print(f"  Analyzing temperatures: {temperatures_K} K")
    
    wavelengths_nm = np.linspace(wavelength_range_nm[0], wavelength_range_nm[1], 500)
    wavelengths_m = wavelengths_nm * 1e-9 # Convert nanometers to meters

    peak_wavelengths_data = [] # To store (Temperature, Peak Wavelength) for Wien's Law plot

    plt.figure(figsize=(12, 7))
    for T in temperatures_K:
        if T <= 0:
            print(f"  Skipping temperature {T}K as it is non-positive.")
            continue
        radiance_series = [calculate_blackbody_spectral_radiance(wl_m, T) for wl_m in wavelengths_m]
        
        # Find the peak wavelength for the current temperature
        max_radiance_idx = np.argmax(radiance_series)
        peak_wl_m = wavelengths_m[max_radiance_idx]
        peak_wl_nm = wavelengths_nm[max_radiance_idx]
        peak_wavelengths_data.append((T, peak_wl_nm))

        plt.plot(wavelengths_nm, radiance_series, label=f'T = {T} K (Peak: {peak_wl_nm:.2f} nm)')
    
    plt.title('Blackbody Spectral Radiance (Planck\'s Law)')
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Spectral Radiance (W/(m³·sr))')
    plt.yscale('log') # Log scale helps visualize broad range of radiance values
    plt.xscale('log') # Log scale for wavelength can also be useful for broad ranges
    plt.grid(True, which="both", ls="--", alpha=0.7)
    plt.legend()
    plt.tight_layout()
    plt.show()

    print("\n--- Wien's Displacement Law Data Series ---")
    print("  This shows how the peak emission wavelength 'displaces' (shifts) to shorter wavelengths as temperature increases.")
    
    temp_list_sorted = [data[0] for data in sorted(peak_wavelengths_data)]
    peak_wl_list = [data[1] for data in sorted(peak_wavelengths_data)]
    
    for T, peak_wl_nm in sorted(peak_wavelengths_data):
        print(f"  Temperature: {T} K, Calculated Peak Wavelength: {peak_wl_nm:.2f} nm")
        
    plt.figure(figsize=(8, 6))
    plt.plot(temp_list_sorted, peak_wl_list, 'o-', label='Calculated Peak Wavelength')
    
    # Calculate and plot the theoretical curve from Wien's Displacement Law
    theoretical_peak_nm = (WIEN_DISPLACEMENT_CONSTANT / np.array(temp_list_sorted)) * 1e9
    plt.plot(temp_list_sorted, theoretical_peak_nm, 'r--', label='Theoretical (Wien\'s Law)')

    plt.title('Wien\'s Displacement Law: Peak Wavelength vs. Temperature')
    plt.xlabel('Temperature (K)')
    plt.ylabel('Peak Wavelength (nm)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

# --- Main Program Execution ---
if __name__ == "__main__":
    print("--- Radiation Analysis Program: Interpreting Data Behavior ---")
    print("\nThis program helps you understand how different scientific concepts can be analyzed using statistical ideas.")
    print("We'll explore two scenarios:")
    print("  1.  **'Nearest Neighbor Boundaries' in Gamma Curies:** This identifies points where the rate of radioactive decay is changing most rapidly, relative to measurement times that follow a bell curve.")
    print("  2.  **Blackbody Wavelength 'Displacement':** This demonstrates how the peak wavelength of light emitted by a hot object shifts as its temperature changes (Wien's Law).\n")
    print("Please note: The terms 'Nearest Neighbor Boundaries' and 'Edge States' are statistical interpretations here, not direct physical phenomena like in quantum mechanics.\n")

    while True:
        main_choice = input("Select: (N)earest Neighbor Boundaries (Gamma Curies), (B)lackbody Wavelength Displacement, or (Q)uit? ").strip().lower()

        if main_choice == 'q':
            break
        elif main_choice == 'n':
            print("\n--- 'Nearest Neighbor Boundaries' for Gamma Curies ---")
            isotope_name = input("  Enter radioactive isotope name (e.g., Cesium-137, Cobalt-60): ").strip()
            # Check if isotope data is available either manually or via PubChem
            if isotope_name not in MANUAL_GAMMA_DATA and _get_compound_cid_pubchem(isotope_name) is None:
                print(f"  '{isotope_name}' data not found. Please check spelling or add it to MANUAL_GAMMA_DATA.")
                continue
            
            try:
                initial_activity = float(input(f"  Enter initial activity of {isotope_name} in Becquerels (Bq): "))
                mean_time = float(input(f"  Enter mean measurement time for the time distribution (days): "))
                std_dev_time = float(input(f"  Enter standard deviation for the measurement time (days): "))
                num_points = int(input(f"  Enter number of data points for the simulation (e.g., 200): "))
                
                deriv_thresh_input = input(f"  Enter absolute derivative threshold for 'boundaries' (Ci/day, leave blank for auto-calc): ")
                derivative_threshold_abs = float(deriv_thresh_input) if deriv_thresh_input else None

            except ValueError:
                print("  Invalid input. Please enter numerical values.")
                continue
            
            analyze_gamma_curies_nearest_neighbor_boundaries(
                isotope_name, initial_activity, mean_time, std_dev_time, num_points, derivative_threshold_abs)
            
            time.sleep(1) # Small pause for readability

        elif main_choice == 'b':
            print("\n--- Blackbody Wavelength 'Displacement' (Wien's Law) ---")
            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. Please enter numbers separated by commas.")
                continue
            
            simulate_blackbody_wien_displacement(temperatures_list)
            
            time.sleep(1) # Small pause for readability

        else:
            print("Invalid choice. Please enter 'N', 'B', or 'Q'.")

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

 

 General Relativism as a K-Body factor for chemical study.

I postulate there is no Universe in a black hole :-o 

 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.
    """
    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

    decay_constant = calculate_decay_constant(half_life_days)
    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    # Apply 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
    
    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

# --- 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)', 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.")
    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!")

 

Exceed Schwarzschild radius? See Cosmic String Blackhole distribution and chemical radius maybe?


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

# --- Constants ---
CURIE_TO_BQ = 3.7e10
LN2 = np.log(2)

STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8
PLANCK_CONSTANT = 6.62607015e-34
BOLTZMANN_CONSTANT = 1.380649e-23
SPEED_OF_LIGHT = 2.99792458e8
WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3

GRAVITATIONAL_CONSTANT = 6.67430e-11
SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2

# --- Helper Functions ---
def get_compound_cid_pubchem(name)
:
"""Return PubChem CID for compound name."""
try:
compounds = pcp.get_compounds(name, 'name')
return compounds[0].cid if compounds else None
except Exception:
return None

def fetch_half_life_pubchem(cid):
"""Fetch half-life value and unit from PubChem JSON data."""
base_url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/cid/{cid}/JSON/"
try:
r = requests.get(base_url, timeout=5)
r.raise_for_status()
data = r.json()
search_sections = ['Atomic Mass, Half Life, and Decay', 'Chemical and Physical Properties', 'Depositor-Supplied']

for section in data.get('Record', {}).get('Section', []):
if section.get('TOCHeading') in search_sections:
for sub_sec in section.get('Section', []):
if 'half-life' in sub_sec.get('TOCHeading', '').lower():
for info in sub_sec.get('Information', []):
for string_dict in info.get('Value', {}).get('StringWithMarkup', []):
text = string_dict.get('String', '')
if text:
m = re.match(r"(\d+\.?\d*(?:e[+-]?
\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)", text.lower())
if m:
return float(m.group(1)), m.group(2)
return None, None
except Exception:
return None, None

def convert_half_life_to_days(
value, unit):
"""Convert half-life to days."""
unit = unit.lower()
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 Calculations ---
def decay_constant(half_life_days)
:
"""Calculate the decay constant from half-life in days."""
if half_life_days is None or half_life_days <= 0:
return 0
seconds = half_life_days * 24 * 3600
return LN2 / seconds

def activity_at_time(initial_
activity_Bq, decay_const, time_seconds):
"""Calculate the activity at a given time."""
if decay_const == 0:
return initial_activity_Bq
return initial_activity_Bq * np.exp(-decay_const * time_seconds)

def schwarzschild_radius(mass_kg):
"""Calculate the Schwarzschild radius of a mass."""
if mass_kg <= 0:
return 0.0
return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def dilation_factor(mass_kg, distance_m):
"""Calculate relativistic time dilation, handle inside and outside Schwarzschild radius."""
if mass_kg <= 0 or distance_m <= 0:
return 1.0
r_s = schwarzschild_radius(mass_kg)

# Inside the event horizon, apply some factor (could be adjusted further).
if distance_m < r_s:
return np.sqrt(1 - (r_s / distance_m)) # Return modified dilation even inside the event horizon.

# Outside event horizon, apply normal time dilation.
return np.sqrt(1 - (r_s / distance_m))

def fetch_isotope_data(isotope_
name):
"""Fetch isotope data (half-life, gamma energies, yields) from PubChem."""
cid = get_compound_cid_pubchem(
isotope_name)
if not cid:
print(f"Isotope {isotope_name} not found in PubChem.")
return None

val, unit = fetch_half_life_pubchem(cid)
if val is None:
return None

half_life_days = convert_half_life_to_days(val, unit)
gamma_energies = [0.662] # Placeholder for actual data
gamma_yields = [0.851] # Placeholder for actual data

return {"half_life_days": half_life_days, "gamma_energies_MeV": gamma_energies, "gamma_yields_per_decay": gamma_yields}

def gamma_curies(isotope, initial_activity_Bq, time_days):
"""Calculate the total gamma radiation (in Curies) for an isotope."""
isotope_data = fetch_isotope_data(isotope)
if isotope_data is None:
return None

half_life = isotope_data["half_life_days"]
decay_const = decay_constant(half_life)
t_seconds = time_days * 24 * 3600
activity = activity_at_time(initial_
activity_Bq, decay_const, t_seconds)
yields = isotope_data["gamma_yields_
per_decay"]
total_yield = sum(yields)
total_gamma_dps = activity * total_yield
return total_gamma_dps / CURIE_TO_BQ

def relativistic_gamma_curies(
isotope, initial_activity_Bq, time_days, mass_kg, distance_m):
"""Calculate gamma radiation (in Curies) under relativistic conditions."""
isotope_data = fetch_isotope_data(isotope)
if isotope_data is None:
return None

half_life = isotope_data["half_life_days"]
decay_const = decay_constant(half_life)
t_seconds = time_days * 24 * 3600
dil_factor = dilation_factor(mass_kg, distance_m)
effective_time = t_seconds * dil_factor

activity = activity_at_time(initial_
activity_Bq, decay_const, effective_time)
yields = isotope_data["gamma_yields_
per_decay"]
total_yield = sum(yields)
total_gamma_dps = activity * total_yield
return total_gamma_dps / CURIE_TO_BQ

# --- Interactive Input Functions ---
def get_float(prompt, positive=False):
"""Get a float from the user, ensuring it's positive if requested."""
while True:
try:
val = float(input(prompt))
if positive and val <= 0:
print("Please enter a positive value.")
continue
return val
except ValueError:
print("Invalid input, please enter a number.")

def get_isotope_name():
"""Prompt the user to enter an isotope name."""
name = input("Enter radioactive isotope name (e.g., Cesium-137): ").strip()
return name

def get_temperature_list():
"""Prompt the user to enter temperatures for blackbody radiation calculations."""
temp_str = input("Enter temperatures in K separated by commas (e.g., 300, 1000, 5000): ")
try:
temps = [float(t.strip()) for t in temp_str.split(",")]
return temps
except Exception:
print("Invalid temperature input.")
return None

# --- Interactive Menus ---
def menu_relativistic_decay():
"""Handle the user input and plot for relativistic nuclear decay comparison."""
print("\n--- Relativistic Nuclear Decay Comparison ---")
isotope = get_isotope_name()
initial_activity = get_float("Enter initial activity (Bq): ", positive=True)
mean_time = get_float("Enter mean measurement time (days): ", positive=True)
std_dev = get_float("Enter std deviation for time (days): ", positive=True)
mass = get_float("Enter mass of k-body (kg): ", positive=True)
distance = get_float("Enter distance from k-body (m): ", positive=True)

from scipy.stats


## Flat space-time

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

# --- Constants ---
CURIE_TO_BQ = 3.7e10
LN2 = np.log(2)

STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8
PLANCK_CONSTANT = 6.62607015e-34
BOLTZMANN_CONSTANT = 1.380649e-23
SPEED_OF_LIGHT = 2.99792458e8
WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3

GRAVITATIONAL_CONSTANT = 6.67430e-11
SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2

# --- Helper Functions ---
def get_compound_cid_pubchem(name)
:
"""Return PubChem CID for compound name."""
try:
compounds = pcp.get_compounds(name, 'name')
return compounds[0].cid if compounds else None
except Exception:
return None

def fetch_half_life_pubchem(cid):
"""Fetch half-life value and unit from PubChem JSON data."""
base_url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/cid/{cid}/JSON/"
try:
r = requests.get(base_url, timeout=5)
r.raise_for_status()
data = r.json()
search_sections = ['Atomic Mass, Half Life, and Decay', 'Chemical and Physical Properties', 'Depositor-Supplied']

for section in data.get('Record', {}).get('Section', []):
if section.get('TOCHeading') in search_sections:
for sub_sec in section.get('Section', []):
if 'half-life' in sub_sec.get('TOCHeading', '').lower():
for info in sub_sec.get('Information', []):
for string_dict in info.get('Value', {}).get('StringWithMarkup', []):
text = string_dict.get('String', '')
if text:
m = re.match(r"(\d+\.?\d*(?:e[+-]?
\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)", text.lower())
if m:
return float(m.group(1)), m.group(2)
return None, None
except Exception:
return None, None

def convert_half_life_to_days(
value, unit):
"""Convert half-life to days."""
unit = unit.lower()
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 Calculations ---
def decay_constant(half_life_days)
:
"""Calculate the decay constant from half-life in days."""
if half_life_days is None or half_life_days <= 0:
return 0
seconds = half_life_days * 24 * 3600
return LN2 / seconds

def activity_at_time(initial_
activity_Bq, decay_const, time_seconds):
"""Calculate the activity at a given time."""
if decay_const == 0:
return initial_activity_Bq
return initial_activity_Bq * np.exp(-decay_const * time_seconds)

def schwarzschild_radius(mass_kg):
"""Calculate the Schwarzschild radius of a mass."""
if mass_kg <= 0:
return 0.0
return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def dilation_factor(mass_kg, distance_m):
"""Calculate relativistic time dilation."""
if mass_kg <= 0 or distance_m <= 0:
return 1.0
r_s = schwarzschild_radius(mass_kg)
if distance_m <= r_s:
return 0.0
return np.sqrt(1 - (r_s / distance_m))

def fetch_isotope_data(isotope_
name):
"""Fetch isotope data (half-life, gamma energies, yields) from PubChem."""
cid = get_compound_cid_pubchem(
isotope_name)
if not cid:
print(f"Isotope {isotope_name} not found in PubChem.")
return None

val, unit = fetch_half_life_pubchem(cid)
if val is None:
return None

half_life_days = convert_half_life_to_days(val, unit)
gamma_energies = [0.662] # Placeholder for actual data
gamma_yields = [0.851] # Placeholder for actual data

return {"half_life_days": half_life_days, "gamma_energies_MeV": gamma_energies, "gamma_yields_per_decay": gamma_yields}

def gamma_curies(isotope, initial_activity_Bq, time_days):
"""Calculate the total gamma radiation (in Curies) for an isotope."""
isotope_data = fetch_isotope_data(isotope)
if isotope_data is None:
return None

half_life = isotope_data["half_life_days"]
decay_const = decay_constant(half_life)
t_seconds = time_days * 24 * 3600
activity = activity_at_time(initial_
activity_Bq, decay_const, t_seconds)
yields = isotope_data["gamma_yields_
per_decay"]
total_yield = sum(yields)
total_gamma_dps = activity * total_yield
return total_gamma_dps / CURIE_TO_BQ

def relativistic_gamma_curies(
isotope, initial_activity_Bq, time_days, mass_kg, distance_m):
"""Calculate gamma radiation (in Curies) under relativistic conditions, but without smoothing."""
isotope_data = fetch_isotope_data(isotope)
if isotope_data is None:
return None

half_life = isotope_data["half_life_days"]
decay_const = decay_constant(half_life)
t_seconds = time_days * 24 * 3600
activity = activity_at_time(initial_
activity_Bq, decay_const, t_seconds)
yields = isotope_data["gamma_yields_
per_decay"]
total_yield = sum(yields)
total_gamma_dps = activity * total_yield
return total_gamma_dps / CURIE_TO_BQ

# --- Interactive Input Functions ---
def get_float(prompt, positive=False):
"""Get a float from the user, ensuring it's positive if requested."""
while True:
try:
val = float(input(prompt))
if positive and val <= 0:
print("Please enter a positive value.")
continue
return val
except ValueError:
print("Invalid input, please enter a number.")

def get_isotope_name():
"""Prompt the user to enter an isotope name."""
name = input("Enter radioactive isotope name (e.g., Cesium-137): ").strip()
return name

def get_temperature_list():
"""Prompt the user to enter temperatures for blackbody radiation calculations."""
temp_str = input("Enter temperatures in K separated by commas (e.g., 300, 1000, 5000): ")
try:
temps = [float(t.strip()) for t in temp_str.split(",")]
return temps
except Exception:
print("Invalid temperature input.")
return None

# --- Interactive Menus ---
def menu_relativistic_decay():
"""Handle the user input and plot for relativistic nuclear decay comparison."""
print("\n--- Nuclear Decay Comparison ---")
isotope = get_isotope_name()
initial_activity = get_float("Enter initial activity (Bq): ", positive=True)
mean_time = get_float("Enter mean measurement time (days): ", positive=True)
std_dev = get_float("Enter std deviation for time (days): ", positive=True)

from scipy.stats import norm
import matplotlib.pyplot as plt
num_points = 200

# Time series with standard deviation analysis
time_series = np.linspace(mean_time - 3*std_dev, mean_time + 3*std_dev, num_points)
time_series = time_series[time_series >= 0] # no negative time

flat_curies = [gamma_curies(isotope, initial_activity, t) for t in time_series]
relativistic_curies = [relativistic_gamma_curies(
isotope, initial_activity, t, 0, 0) for t in time_series]

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(time_series, flat_curies, label='Flat Spacetime Decay', color='blue')
plt.plot(time_series, relativistic_curies, label='Relativistic Decay', color='red', linestyle='--')
plt.title(f"Gamma Radiation vs Time for {isotope}")
plt.xlabel("Time (days)")
plt.ylabel("Gamma Radiation (Curies)")
plt.legend()
plt.grid(True)
plt.show()

# Run the menu for relativistic decay
menu_relativistic_decay()

 

No comments: