Radical! Nearest neighbour tracing in nuclear automata.
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!")
