# -------------------------------------------------------------------------------
# Cold Fusion Simulation: Deuterium-Based Low Energy Nuclear Reactions in a Lattice
# -------------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import (
    e,          # Elementary charge (C)
    epsilon_0,  # Vacuum permittivity (F/m)
    hbar,       # Reduced Planck constant (J·s)
    k,          # Boltzmann constant (J/K)
    physical_constants,
    N_A         # Avogadro's number (1/mol)
)
from scipy.integrate import quad
from scipy.optimize import root_scalar
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.animation import FuncAnimation
import matplotlib.ticker as ticker

# ----------------------------
# 1. Physical Constants
# ----------------------------

# Atomic mass unit (kg)
u = physical_constants['atomic mass constant'][0]

# Masses (kg)
m_p = physical_constants['proton mass'][0]
m_n = physical_constants['neutron mass'][0]
m_D = m_p + m_n         # Deuteron mass
m_T = m_p + 2 * m_n     # Triton mass
m_He3 = 2 * m_p + m_n   # Helium-3 mass

# Reduced masses (kg)
mu_DD = m_D / 2                         # Reduced mass for D-D reaction
mu_DT = (m_D * m_T) / (m_D + m_T)       # Reduced mass for D-T reaction

# Fine-structure constant
alpha = physical_constants['fine-structure constant'][0]

# Charges
Z1 = 1  # Charge number of deuteron
Z2 = 1  # Charge number of reacting nucleus (1 for D, 1 for T)

# -----------------------------------
# 2. Functions for Fusion Reactions
# -----------------------------------

def S_factor(E, reaction='DD'):
    """
    Astrophysical S-factor for fusion reactions.
    E: Energy in joules.
    reaction: 'DD' for deuterium-deuterium, 'DT' for deuterium-tritium.
    Returns S(E) in joule·m².
    """
    # Convert keV·barn to joule·m² (1 keV = 1.60218e-16 J, 1 barn = 1e-28 m²)
    if reaction == 'DD':
        S0_keVbarn = 50  # S(0) for D-D fusion ~ 50 keV·barn
    elif reaction == 'DT':
        S0_keVbarn = 100  # S(0) for D-T fusion ~ 100 keV·barn
    else:
        raise ValueError("Invalid reaction type. Choose 'DD' or 'DT'.")

    S0 = S0_keVbarn * 1.60218e-16 * 1e-28
    return S0

def eta(E, mu):
    """
    Calculate the Sommerfeld parameter.
    E: Energy in joules.
    mu: Reduced mass of the reacting nuclei.
    Returns dimensionless Sommerfeld parameter.
    """
    return (Z1 * Z2 * e**2) / (4 * np.pi * epsilon_0 * hbar) * np.sqrt(mu / (2 * E))

def sigma_fusion(E, mu, reaction='DD'):
    """
    Fusion cross-section for fusion reactions at energy E (joules).
    E: Energy in joules.
    mu: Reduced mass of the reacting nuclei.
    reaction: 'DD' or 'DT'.
    Returns cross-section in m².
    """
    S = S_factor(E, reaction)
    eta_value = eta(E, mu)
    sigma = (S / E) * np.exp(-2 * np.pi * eta_value)
    return sigma

def maxwell_boltzmann(E, T):
    """
    Maxwell-Boltzmann energy distribution at temperature T (kelvin).
    E: Energy in joules.
    Returns probability density in J⁻¹.
    """
    return 2 * np.sqrt(E / np.pi) * (1 / (k * T)**1.5) * np.exp(-E / (k * T))

def reaction_rate(T, mu, Ue=0, reaction='DD'):
    """
    Calculate the velocity-averaged reaction rate <σv>.
    T: Temperature in kelvin.
    mu: Reduced mass of the reacting nuclei.
    Ue: Electron screening potential in joules.
    reaction: 'DD' or 'DT'.
    Returns reaction rate in m³/s.
    """
    def integrand(E):
        E_eff = E + Ue  # Effective energy due to screening
        sigma = sigma_fusion(E_eff, mu, reaction)
        v = np.sqrt(2 * E / mu)
        return sigma * v * maxwell_boltzmann(E, T)
    
    # Integration limits: from 0 to infinity
    integral, _ = quad(integrand, 0, np.inf, limit=1000)
    return integral

def number_density_D(theta):
    """
    Calculate the number density of deuterons in the lattice.
    theta: Deuterium loading ratio (D/Metal).
    Returns number density in m⁻³.
    """
    # For palladium (Pd)
    rho_M = 12.0e3             # Density of Pd (kg/m³)
    M_M = 106.42e-3            # Molar mass of Pd (kg/mol)
    n_M = rho_M * N_A / M_M    # Number density of Pd atoms (atoms/m³)
    n_D = theta * n_M          # Number density of deuterons (atoms/m³)
    return n_D

def total_reaction_rate(T, theta, Ue=0, reaction='DD'):
    """
    Calculate the total reaction rate per unit volume.
    T: Temperature in kelvin.
    theta: Deuterium loading ratio.
    Ue: Electron screening potential in joules.
    reaction: 'DD' or 'DT'.
    Returns reaction rate per unit volume in reactions/m³/s.
    """
    if reaction == 'DD':
        mu = mu_DD
        n_D1 = number_density_D(theta)  # Deuteron density
        n_D2 = n_D1                     # For D-D, both reactants are deuterons
    elif reaction == 'DT':
        mu = mu_DT
        n_D1 = number_density_D(theta)  # Deuteron density
        # Tritium density (assuming tritium occupies the same sites)
        n_D2 = number_density_D(theta)  # For simplicity, using same loading
    else:
        raise ValueError("Invalid reaction type. Choose 'DD' or 'DT'.")
    
    rate = n_D1 * n_D2 * reaction_rate(T, mu, Ue, reaction)
    if reaction == 'DD':
        rate *= 0.5  # Correcting for identical particles
    return rate

# -----------------------------------------------
# 3. User Interaction and Parameter Definitions
# -----------------------------------------------

# Allow the user to set parameters or use default values

# Temperature range (K)
T_min = float(input("Enter minimum temperature in K (default 0.1 K): ") or 0.1)
T_max = float(input("Enter maximum temperature in K (default 1e5 K): ") or 1e5)
num_T_points = int(input("Enter number of temperature points (default 100): ") or 100)

# Deuterium loading ratio
theta_min = float(input("Enter minimum deuterium loading ratio (0-1, default 0.1): ") or 0.1)
theta_max = float(input("Enter maximum deuterium loading ratio (0-1, default 1.0): ") or 1.0)
num_theta_points = int(input("Enter number of loading ratio points (default 50): ") or 50)

# Screening potential range (eV)
Ue_min = float(input("Enter minimum screening potential in eV (default 0 eV): ") or 0)
Ue_max = float(input("Enter maximum screening potential in eV (default 1000 eV): ") or 1000)
num_Ue_points = int(input("Enter number of screening potential points (default 50): ") or 50)

# Reaction choice
reaction_choice = input("Choose reaction ('DD' or 'DT', default 'DD'): ") or 'DD'

# Convert ranges to arrays
temperatures = np.logspace(np.log10(T_min), np.log10(T_max), num_T_points)
thetas = np.linspace(theta_min, theta_max, num_theta_points)
Ue_eV = np.linspace(Ue_min, Ue_max, num_Ue_points)
Ue_J = Ue_eV * e  # Convert eV to joules

# -------------------------------------------------
# 4. Simulation and Visualization
# -------------------------------------------------

# Helper function for pretty axis formatting
def format_axis(ax):
    ax.grid(True, which='both', ls='--')
    ax.set_facecolor('#f0f0f0')
    ax.tick_params(axis='both', which='major')
    return ax

# 4.1. Reaction Rate vs Temperature
def plot_reaction_rate_vs_temperature():
    rates_T = np.array([total_reaction_rate(T, theta_max, Ue=0, reaction=reaction_choice) for T in temperatures])
    
    plt.figure(figsize=(8, 6))
    plt.loglog(temperatures, rates_T, label=f'Reaction: {reaction_choice}')
    plt.xlabel('Temperature (K)')
    plt.ylabel('Reaction Rate per Unit Volume (reactions/m³/s)')
    plt.title('Fusion Reaction Rate vs Temperature')
    plt.legend()
    format_axis(plt.gca())
    plt.show()

# 4.2. Reaction Rate vs Deuterium Loading Ratio
def plot_reaction_rate_vs_loading_ratio():
    rates_theta = np.array([total_reaction_rate(300, theta, Ue=0, reaction=reaction_choice) for theta in thetas])
    
    plt.figure(figsize=(8, 6))
    plt.plot(thetas, rates_theta)
    plt.xlabel('Deuterium Loading Ratio (D/Metal)')
    plt.ylabel('Reaction Rate per Unit Volume (reactions/m³/s)')
    plt.title('Fusion Reaction Rate vs Deuterium Loading Ratio at 300 K')
    format_axis(plt.gca())
    plt.show()

# 4.3. 3D Plot: Reaction Rate vs Temperature and Loading Ratio
def plot_3D_reaction_rate():
    T_grid, theta_grid = np.meshgrid(temperatures, thetas)
    T_flat = T_grid.flatten()
    theta_flat = theta_grid.flatten()

    rates_flat = np.array([total_reaction_rate(T, theta, Ue=0, reaction=reaction_choice) for T, theta in zip(T_flat, theta_flat)])
    rates_grid = rates_flat.reshape(T_grid.shape)

    fig = plt.figure(figsize=(12, 8))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(np.log10(T_grid), theta_grid, np.log10(rates_grid), cmap=cm.viridis)

    ax.set_xlabel('log10(Temperature) (log10 K)')
    ax.set_ylabel('Deuterium Loading Ratio (D/Metal)')
    ax.set_zlabel('log10(Reaction Rate) (log10 reactions/m³/s)')
    ax.set_title('Reaction Rate vs Temperature and Deuterium Loading Ratio')
    fig.colorbar(surf, shrink=0.5, aspect=10, label='log10(Reactions/m³/s)')
    plt.show()

# 4.4. Reaction Rate vs Electron Screening Potential
def plot_reaction_rate_vs_screening_potential():
    rates_Ue = np.array([total_reaction_rate(300, theta_max, Ue, reaction=reaction_choice) for Ue in Ue_J])
    
    plt.figure(figsize=(8, 6))
    plt.plot(Ue_eV, rates_Ue)
    plt.xlabel('Screening Potential Ue (eV)')
    plt.ylabel('Reaction Rate per Unit Volume (reactions/m³/s)')
    plt.title('Reaction Rate vs Electron Screening Potential at 300 K')
    format_axis(plt.gca())
    plt.show()

# 4.5. Animated Plot: Reaction Rate vs Temperature for Various Ue
def animate_reaction_rate_vs_temperature():
    fig, ax = plt.subplots(figsize=(8, 6))
    line, = ax.loglog([], [], lw=2)
    ax.set_xlabel('Temperature (K)')
    ax.set_ylabel('Reaction Rate per Unit Volume (reactions/m³/s)')
    ax.set_title('Fusion Reaction Rate vs Temperature for Various Ue')
    format_axis(ax)

    def init():
        line.set_data([], [])
        return line,

    def update(Ue_eV_value):
        Ue_J_value = Ue_eV_value * e
        rates_T = np.array([total_reaction_rate(T, theta_max, Ue=Ue_J_value, reaction=reaction_choice) for T in temperatures])
        line.set_data(temperatures, rates_T)
        ax.legend([f'Ue = {Ue_eV_value:.1f} eV'])
        return line,

    ani = FuncAnimation(fig, update, frames=Ue_eV, init_func=init, blit=False, interval=200)
    plt.show()

# 4.6. Fusion Cross-Section vs Energy
def plot_cross_section_vs_energy():
    energies_keV = np.logspace(-1, 3, 1000)  # 0.1 keV to 1000 keV
    energies_J = energies_keV * 1.60218e-16   # Convert keV to joules
    mu = mu_DD if reaction_choice == 'DD' else mu_DT
    sigma_values = [sigma_fusion(E, mu, reaction_choice) for E in energies_J]

    plt.figure(figsize=(8, 6))
    plt.loglog(energies_keV, sigma_values)
    plt.xlabel('Energy (keV)')
    plt.ylabel('Fusion Cross-Section (m²)')
    plt.title(f'Fusion Cross-Section vs Energy for {reaction_choice} Reaction')
    format_axis(plt.gca())
    plt.show()

# -------------------------------------------------
# 5. Execute Plots and Simulations Based on User Input
# -------------------------------------------------

print("\nGenerating plots...")

plot_reaction_rate_vs_temperature()
plot_reaction_rate_vs_loading_ratio()
plot_3D_reaction_rate()
plot_reaction_rate_vs_screening_potential()
plot_cross_section_vs_energy()

animate_choice = input("Do you want to see the animated plot? ('yes' or 'no', default 'no'): ") or 'no'
if animate_choice.lower() == 'yes':
    animate_reaction_rate_vs_temperature()

print("Simulation and plotting completed.") do erro fixing