# LENRARA Python Simulation Suite with GUI (WORK IN PROGRESS)
# The following is a set of computational function variables to help us simulate Hydrogen Fusion Processes in Condensed Matter Nuclear Science
# License GNU 3.0 Conscious Energy 2024

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.figure import Figure
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation  # Import FuncAnimation
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from scipy.integrate import odeint, solve_ivp
import tkinter as tk
from tkinter import ttk
from tkinter import messagebox  # Import messagebox

# Constants
k_B = 1.380649e-23  # Boltzmann constant, J/K
hbar = 1.0545718e-34  # Reduced Planck constant, J·s
e = 1.60217662e-19  # Elementary charge, C
c = 3.0e8  # Speed of light, m/s
Z_1 = 1  # Atomic number for hydrogen
Z_2 = 1  # Atomic number for hydrogen

# Thermodynamic properties
delta_H = -40e3  # Enthalpy of formation, J/mol H2
delta_S = 130  # Entropy, J/mol·K H2

def create_fcc_lattice(a, n):
    """Create FCC lattice."""
    positions = []
    basis_vectors = [(0, 0, 0), (0.5, 0.5, 0), (0.5, 0, 0.5), (0, 0.5, 0.5)]
    for i in range(n):
        for j in range(n):
            for k_pos in range(n):
                for vec in basis_vectors:
                    pos = (i + vec[0], j + vec[1], k_pos + vec[2])
                    positions.append(pos)
    return np.array(positions) * a

def vant_hoff(delta_H, delta_S, T):
    """Calculate equilibrium pressure using Van't Hoff equation."""
    R = 8.314  # J/(mol·K)
    return np.exp((delta_H - T * delta_S) / (R * T))

def fusion_cross_section(T):
    """Calculate fusion cross-section."""
    term1 = (e**2 / (4 * np.pi * hbar * c))**2
    term2 = (1 / (k_B * T))**2
    term3 = np.exp(-3 * np.pi / (4 * np.sqrt(2)) * (Z_1 * Z_2 * e**2 / (hbar * c))**2 / (k_B * T))
    return 1e-24 * term1 * term2 * term3

def fusion_rate(n1, n2, sigma, v_rel):
    """Calculate fusion rate."""
    return n1 * n2 * sigma * v_rel

def simulate_fusion(T, n1, n2, v_rel):
    """Simulate LENR fusion event."""
    sigma = fusion_cross_section(T)
    rate = fusion_rate(n1, n2, sigma, v_rel)
    return 1 / rate, sigma

def excess_heat_generation(time, P_excess):
    """Calculate excess heat generation."""
    return np.trapz(P_excess, time)

def transmutation_model(y, t, k):
    """Transmutation pathways model."""
    return -k * y

def glow_discharge_power(I, V):
    """Calculate glow discharge power."""
    return I * V

def neutron_production_rate(sigma, rho_e, rho_p, E_e):
    """Calculate neutron production rate."""
    return hbar * sigma * rho_e * rho_p * E_e

def tsc_fusion_rate(S, G, T):
    """Calculate TSC fusion rate."""
    return S * G * T

def neutron_capture(y, t, n_density, capture_rate):
    """Modeling neutron capture and decay."""
    N, D, He, gamma_photons = y
    dNdt = -capture_rate * n_density * N
    dDdt = capture_rate * n_density * N - capture_rate * D
    dHedt = capture_rate * D
    dGammadt = capture_rate * D
    return [dNdt, dDdt, dHedt, dGammadt]

def rayleigh_plesset(t, R, P0, Pv, rho, gamma, mu, sigma):
    """Modeling bubble dynamics under sonication."""
    R_val, R_dot = R
    # Avoid division by zero
    if R_val == 0:
        R_val = 1e-9
    P = P0 - (2 * sigma / R_val) - (4 * mu * R_dot / R_val)
    P_bubble = Pv + (2 * sigma / R_val)
    R_ddot = (P - P_bubble) / (rho * R_val) - (3 / 2) * (R_dot ** 2 / R_val)
    return [R_dot, R_ddot]

def sonication_simulation(nx, nt, dx, dt, source_position, c=1500):
    """Simulate sonication using FDTD."""
    # Initialize pressure fields
    p_current = np.zeros(nx)
    p_previous = np.zeros(nx)
    
    # List to store pressure snapshots
    pressure_history = []
    
    # Time-stepping loop
    for t in range(nt):
        # Simple wave propagation (this is a placeholder; replace with your actual simulation logic)
        p_next = np.zeros(nx)
        for i in range(1, nx-1):
            p_next[i] = (2 * p_current[i] - p_previous[i] +
                         (c * dt / dx)**2 * (p_current[i+1] - 2*p_current[i] + p_current[i-1]))
        
        # Add source (e.g., a sine wave source)
        p_next[source_position] += np.sin(2 * np.pi * 50 * t * dt)
        
        # Update fields
        p_previous, p_current = p_current, p_next
        
        # Store the current pressure field
        pressure_history.append(p_current.copy())
    
    return pressure_history

# GUI Application
class LENRSimulationApp:
    def __init__(self, master):
        self.master = master
        master.title("LENRARA Python Simulation Suite")
        self.create_widgets()

    def create_widgets(self):
        # Create Notebook (Tabs)
        self.tabControl = ttk.Notebook(self.master)

        # Create tabs
        self.tab1 = ttk.Frame(self.tabControl)
        self.tabControl.add(self.tab1, text='Equilibrium Pressure')

        self.tab2 = ttk.Frame(self.tabControl)
        self.tabControl.add(self.tab2, text='LENR Event Simulation')

        self.tab4 = ttk.Frame(self.tabControl)
        self.tabControl.add(self.tab4, text='Transmutation Calculation')

        self.tab7 = ttk.Frame(self.tabControl)
        self.tabControl.add(self.tab7, text='TSC Fusion Rate')

        self.tab8 = ttk.Frame(self.tabControl)
        self.tabControl.add(self.tab8, text='Sonication Simulation')

        self.tabControl.pack(expand=1, fill="both")

        # Create individual tabs
        self.create_equilibrium_pressure_tab()
        self.create_lenr_event_simulation_tab()
        self.create_excess_heat_generation_tab()
        self.create_transmutation_calculation_tab()
        self.create_glow_discharge_tab()
        self.create_neutron_production_rate_tab()
        self.create_tsc_fusion_rate_tab()
        self.create_sonication_simulation_tab()
        self.create_neutron_capture_tab()
        self.create_bubble_dynamics_tab()

    def create_equilibrium_pressure_tab(self):
        # Equilibrium Pressure Tab
        frame = self.tab1

        ttk.Label(frame, text="Temperature (K):").grid(column=0, row=0, padx=5, pady=5)
        self.ep_temperature_entry = ttk.Entry(frame)
        self.ep_temperature_entry.insert(0, "300")
        self.ep_temperature_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_equilibrium_pressure).grid(
            column=0, row=1, padx=5, pady=5, columnspan=2
        )

        self.ep_result_label = ttk.Label(frame, text="Equilibrium Pressure: ")
        self.ep_result_label.grid(column=0, row=2, padx=5, pady=5, columnspan=2)

    def calculate_equilibrium_pressure(self):
        try:
            T = float(self.ep_temperature_entry.get())
            if T <= 0:
                raise ValueError("Temperature must be positive.")
            equilibrium_pressure = vant_hoff(delta_H, delta_S, T)
            self.ep_result_label.config(text=f"Equilibrium Pressure: {equilibrium_pressure:.2e} Pa")
        except ValueError:
            self.ep_result_label.config(text="Invalid input. Please enter a positive numeric value.")

    def create_lenr_event_simulation_tab(self):
        # LENR Event Simulation Tab
        frame = self.tab2

        ttk.Label(frame, text="Temperature (K):").grid(column=0, row=0, padx=5, pady=5)
        self.les_temperature_entry = ttk.Entry(frame)
        self.les_temperature_entry.insert(0, "300")
        self.les_temperature_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Number Density n1 (1/m³):").grid(column=0, row=1, padx=5, pady=5)
        self.les_n1_entry = ttk.Entry(frame)
        self.les_n1_entry.insert(0, "1e28")
        self.les_n1_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Label(frame, text="Number Density n2 (1/m³):").grid(column=0, row=2, padx=5, pady=5)
        self.les_n2_entry = ttk.Entry(frame)
        self.les_n2_entry.insert(0, "1e28")
        self.les_n2_entry.grid(column=1, row=2, padx=5, pady=5)

        ttk.Label(frame, text="Relative Velocity (m/s):").grid(column=0, row=3, padx=5, pady=5)
        self.les_vrel_entry = ttk.Entry(frame)
        self.les_vrel_entry.insert(0, "1e6")
        self.les_vrel_entry.grid(column=1, row=3, padx=5, pady=5)

        ttk.Button(frame, text="Simulate", command=self.simulate_lenr_event).grid(
            column=0, row=4, padx=5, pady=5, columnspan=2
        )

        # Plot area
        self.les_fig = Figure(figsize=(6, 5), dpi=100)
        self.les_canvas = FigureCanvasTkAgg(self.les_fig, master=frame)
        self.les_canvas.get_tk_widget().grid(column=0, row=5, columnspan=2)

        # Control buttons for animation
        control_frame = ttk.Frame(frame)
        control_frame.grid(column=0, row=6, columnspan=2)
        self.pause_button = ttk.Button(control_frame, text="Pause", command=self.pause_animation)
        self.pause_button.grid(column=0, row=0, padx=5, pady=5)
        self.resume_button = ttk.Button(control_frame, text="Resume", command=self.resume_animation)
        self.resume_button.grid(column=1, row=0, padx=5, pady=5)
        self.speed_scale = ttk.Scale(
            control_frame, from_=0.1, to=5.0, value=1.0, orient='horizontal', command=self.update_speed
        )
        self.speed_scale.grid(column=2, row=0, padx=5, pady=5)
        ttk.Label(control_frame, text="Speed").grid(column=2, row=1)

    def simulate_lenr_event(self):
        try:
            a = 1.0  # Lattice parameter
            n_cells = 3  # Number of unit cells along each axis
            fcc_lattice = create_fcc_lattice(a, n_cells)

            T = float(self.les_temperature_entry.get())
            if T <= 0:
                raise ValueError("Temperature must be positive.")
            n1 = float(self.les_n1_entry.get())
            n2 = float(self.les_n2_entry.get())
            v_rel = float(self.les_vrel_entry.get())
            if n1 < 0 or n2 < 0 or v_rel < 0:
                raise ValueError("Number densities and relative velocity must be non-negative.")

            time_to_fusion, sigma = simulate_fusion(T, n1, n2, v_rel)
            decay_positions = fcc_lattice.copy()
            decay_times = np.random.exponential(scale=time_to_fusion, size=decay_positions.shape[0])

            # Animation variables
            self.current_time = 0
            self.max_time = np.max(decay_times)
            self.decay_positions = decay_positions
            self.decay_times = decay_times
            self.n_cells = n_cells  # Store n_cells for use in update_animation

            # Set up the plot
            self.les_fig.clf()
            self.ax = self.les_fig.add_subplot(111, projection='3d')
            self.ax.set_xlim([0, n_cells])
            self.ax.set_ylim([0, n_cells])
            self.ax.set_zlim([0, n_cells])
            self.ax.set_title(f'LENR Event Simulation\nTime to Fusion: {time_to_fusion:.2e} s\nFusion Cross Section: {sigma:.2e} m²')

            # Start the animation
            self.is_paused = False
            self.animation_speed = 1.0  # Default speed
            self.anim = FuncAnimation(
                self.les_fig,
                self.update_animation,
                frames=np.linspace(0, self.max_time, 200),
                interval=50,  # milliseconds between frames
                blit=False,
                repeat=False,
            )
            self.les_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")
        except Exception as e:
            messagebox.showerror("Error", f"An unexpected error occurred: {e}")

    def update_animation(self, frame):
        if self.is_paused:
            return
        self.current_time = frame * self.animation_speed
        if self.current_time > self.max_time:
            self.current_time = self.max_time

        # Filter positions where decay time is greater than current time
        remaining_indices = np.where(self.decay_times > self.current_time)[0]
        self.ax.cla()
        self.ax.set_xlim([0, self.n_cells])
        self.ax.set_ylim([0, self.n_cells])
        self.ax.set_zlim([0, self.n_cells])
        self.ax.set_title(f'LENR Event Simulation\nTime: {self.current_time:.2e} s')
        self.ax.scatter(
            self.decay_positions[remaining_indices, 0],
            self.decay_positions[remaining_indices, 1],
            self.decay_positions[remaining_indices, 2],
            c='r',
            s=20,
        )
        self.les_canvas.draw()

    def pause_animation(self):
        self.is_paused = True

    def resume_animation(self):
        self.is_paused = False

    def update_speed(self, val):
        try:
            self.animation_speed = float(val)
        except ValueError:
            self.animation_speed = 1.0  # Reset to default if invalid

    def create_excess_heat_generation_tab(self):
        # Excess Heat Generation Tab
        frame = self.tab3

        ttk.Label(frame, text="Time Duration (s):").grid(column=0, row=0, padx=5, pady=5)
        self.ehg_time_entry = ttk.Entry(frame)
        self.ehg_time_entry.insert(0, "100")
        self.ehg_time_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_excess_heat_generation).grid(
            column=0, row=1, padx=5, pady=5, columnspan=2
        )

        self.ehg_result_label = ttk.Label(frame, text="Excess Heat Generated: ")
        self.ehg_result_label.grid(column=0, row=2, padx=5, pady=5, columnspan=2)

        # Plot area
        self.ehg_fig = Figure(figsize=(6, 5), dpi=100)
        self.ehg_canvas = FigureCanvasTkAgg(self.ehg_fig, master=frame)
        self.ehg_canvas.get_tk_widget().grid(column=0, row=3, columnspan=2)

    def calculate_excess_heat_generation(self):
        try:
            duration = float(self.ehg_time_entry.get())
            if duration <= 0:
                raise ValueError("Duration must be positive.")
            time = np.linspace(0, duration, 1000)
            P_excess = np.sin(time / 10) + 1  # Arbitrary excess power pattern
            Q_excess = excess_heat_generation(time, P_excess)
            self.ehg_result_label.config(text=f"Excess Heat Generated: {Q_excess:.2e} J")

            # Plotting
            self.ehg_fig.clf()
            ax = self.ehg_fig.add_subplot(111)
            ax.plot(time, P_excess)
            ax.set_xlabel('Time [s]')
            ax.set_ylabel('Excess Power [W]')
            ax.set_title('Excess Power over Time')
            self.ehg_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")

    def create_transmutation_calculation_tab(self):
        # Transmutation Calculation Tab
        frame = self.tab4

        ttk.Label(frame, text="Time Duration (s):").grid(column=0, row=0, padx=5, pady=5)
        self.tc_time_entry = ttk.Entry(frame)
        self.tc_time_entry.insert(0, "100")
        self.tc_time_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Rate Constant k (1/s):").grid(column=0, row=1, padx=5, pady=5)
        self.tc_k_entry = ttk.Entry(frame)
        self.tc_k_entry.insert(0, "0.1")
        self.tc_k_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Label(frame, text="Initial Concentration y0 (1/m³):").grid(column=0, row=2, padx=5, pady=5)
        self.tc_y0_entry = ttk.Entry(frame)
        self.tc_y0_entry.insert(0, "1e23")
        self.tc_y0_entry.grid(column=1, row=2, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_transmutation).grid(
            column=0, row=3, padx=5, pady=5, columnspan=2
        )

        self.tc_result_label = ttk.Label(frame, text="Final Concentration: ")
        self.tc_result_label.grid(column=0, row=4, padx=5, pady=5, columnspan=2)

        # Plot area
        self.tc_fig = Figure(figsize=(6, 5), dpi=100)
        self.tc_canvas = FigureCanvasTkAgg(self.tc_fig, master=frame)
        self.tc_canvas.get_tk_widget().grid(column=0, row=5, columnspan=2)

    def calculate_transmutation(self):
        try:
            duration = float(self.tc_time_entry.get())
            k = float(self.tc_k_entry.get())
            y0 = float(self.tc_y0_entry.get())
            if duration <= 0 or k < 0 or y0 < 0:
                raise ValueError("Duration must be positive and rates/concentrations cannot be negative.")
            t = np.linspace(0, duration, 1000)
            transmutation_result = odeint(transmutation_model, y0, t, args=(k,))
            final_concentration = transmutation_result[-1][0]
            self.tc_result_label.config(text=f"Final Concentration: {final_concentration:.2e} 1/m³")

            # Plotting
            self.tc_fig.clf()
            ax = self.tc_fig.add_subplot(111)
            ax.plot(t, transmutation_result[:, 0], label='Neutron (N)')
            ax.plot(t, transmutation_result[:, 1], label='Deuterium (D)')
            ax.plot(t, transmutation_result[:, 2], label='Helium-3 (He-3)')
            ax.plot(t, transmutation_result[:, 3], label='Gamma photons')
            ax.set_xlabel('Time [s]')
            ax.set_ylabel('Concentration [1/m³]')
            ax.legend()
            ax.set_title('Transmutation over Time')
            self.tc_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")
        except Exception as e:
            messagebox.showerror("Error", f"An unexpected error occurred: {e}")

    def create_glow_discharge_tab(self):
        # Glow Discharge Tab
        frame = self.tab5

        ttk.Label(frame, text="Current I (A):").grid(column=0, row=0, padx=5, pady=5)
        self.gd_current_entry = ttk.Entry(frame)
        self.gd_current_entry.insert(0, "1.0")
        self.gd_current_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Voltage V (V):").grid(column=0, row=1, padx=5, pady=5)
        self.gd_voltage_entry = ttk.Entry(frame)
        self.gd_voltage_entry.insert(0, "100")
        self.gd_voltage_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_glow_discharge).grid(
            column=0, row=2, padx=5, pady=5, columnspan=2
        )

        self.gd_result_label = ttk.Label(frame, text="Glow Discharge Power: ")
        self.gd_result_label.grid(column=0, row=3, padx=5, pady=5, columnspan=2)

    def calculate_glow_discharge(self):
        try:
            I = float(self.gd_current_entry.get())
            V = float(self.gd_voltage_entry.get())
            if I < 0 or V < 0:
                raise ValueError("Current and Voltage must be non-negative.")
            P_glow = glow_discharge_power(I, V)
            self.gd_result_label.config(text=f"Glow Discharge Power: {P_glow:.2e} W")
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")

    def create_neutron_production_rate_tab(self):
        # Neutron Production Rate Tab
        frame = self.tab6

        ttk.Label(frame, text="Cross-section σ (m²):").grid(column=0, row=0, padx=5, pady=5)
        self.npr_sigma_entry = ttk.Entry(frame)
        self.npr_sigma_entry.insert(0, "1e-28")
        self.npr_sigma_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Electron Density ρe (1/m³):").grid(column=0, row=1, padx=5, pady=5)
        self.npr_rhoe_entry = ttk.Entry(frame)
        self.npr_rhoe_entry.insert(0, "1e28")
        self.npr_rhoe_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Label(frame, text="Proton Density ρp (1/m³):").grid(column=0, row=2, padx=5, pady=5)
        self.npr_rhop_entry = ttk.Entry(frame)
        self.npr_rhop_entry.insert(0, "1e28")
        self.npr_rhop_entry.grid(column=1, row=2, padx=5, pady=5)

        ttk.Label(frame, text="Electron Energy Ee (eV):").grid(column=0, row=3, padx=5, pady=5)
        self.npr_Ee_entry = ttk.Entry(frame)
        self.npr_Ee_entry.insert(0, "1e6")
        self.npr_Ee_entry.grid(column=1, row=3, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_neutron_production_rate).grid(
            column=0, row=4, padx=5, pady=5, columnspan=2
        )

        self.npr_result_label = ttk.Label(frame, text="Neutron Production Rate: ")
        self.npr_result_label.grid(column=0, row=5, padx=5, pady=5, columnspan=2)

    def calculate_neutron_production_rate(self):
        try:
            sigma = float(self.npr_sigma_entry.get())
            rho_e = float(self.npr_rhoe_entry.get())
            rho_p = float(self.npr_rhop_entry.get())
            E_e = float(self.npr_Ee_entry.get())
            if sigma < 0 or rho_e < 0 or rho_p < 0 or E_e < 0:
                raise ValueError("All physical parameters must be non-negative.")
            R_n = neutron_production_rate(sigma, rho_e, rho_p, E_e)
            self.npr_result_label.config(text=f"Neutron Production Rate: {R_n:.2e} neutrons/s")
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")

    def create_tsc_fusion_rate_tab(self):
        # TSC Fusion Rate Tab
        frame = self.tab7

        ttk.Label(frame, text="Symmetry Factor S:").grid(column=0, row=0, padx=5, pady=5)
        self.tfr_S_entry = ttk.Entry(frame)
        self.tfr_S_entry.insert(0, "1.0")
        self.tfr_S_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Gibbs Free Energy G (J/mol):").grid(column=0, row=1, padx=5, pady=5)
        self.tfr_G_entry = ttk.Entry(frame)
        self.tfr_G_entry.insert(0, "1.0")
        self.tfr_G_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Label(frame, text="Temperature T (K):").grid(column=0, row=2, padx=5, pady=5)
        self.tfr_T_entry = ttk.Entry(frame)
        self.tfr_T_entry.insert(0, "300")
        self.tfr_T_entry.grid(column=1, row=2, padx=5, pady=5)

        ttk.Button(frame, text="Calculate", command=self.calculate_tsc_fusion_rate).grid(
            column=0, row=3, padx=5, pady=5, columnspan=2
        )

        self.tfr_result_label = ttk.Label(frame, text="TSC Fusion Rate: ")
        self.tfr_result_label.grid(column=0, row=4, padx=5, pady=5, columnspan=2)

    def calculate_tsc_fusion_rate(self):
        try:
            S = float(self.tfr_S_entry.get())
            G = float(self.tfr_G_entry.get())
            T = float(self.tfr_T_entry.get())
            if T <= 0:
                raise ValueError("Temperature must be positive.")
            R_TSC = tsc_fusion_rate(S, G, T)
            self.tfr_result_label.config(text=f"TSC Fusion Rate: {R_TSC:.2e} reactions/s")
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")

    def create_sonication_simulation_tab(self):
        # Sonication Simulation Tab
        frame = self.tab8

        ttk.Label(frame, text="Number of spatial steps nx:").grid(column=0, row=0, padx=5, pady=5)
        self.ss_nx_entry = ttk.Entry(frame)
        self.ss_nx_entry.insert(0, "200")
        self.ss_nx_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Number of time steps nt:").grid(column=0, row=1, padx=5, pady=5)
        self.ss_nt_entry = ttk.Entry(frame)
        self.ss_nt_entry.insert(0, "750")
        self.ss_nt_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Button(frame, text="Simulate", command=self.simulate_sonication).grid(
            column=0, row=2, padx=5, pady=5, columnspan=2
        )

        # Plot area
        self.ss_fig = Figure(figsize=(6, 5), dpi=100)
        self.ss_canvas = FigureCanvasTkAgg(self.ss_fig, master=frame)
        self.ss_canvas.get_tk_widget().grid(column=0, row=3, columnspan=2)

        # Animation Control Buttons
        control_frame = ttk.Frame(frame)
        control_frame.grid(column=0, row=4, columnspan=2, pady=5)

        self.ss_pause_button = ttk.Button(control_frame, text="Pause", command=self.pause_ss_animation)
        self.ss_pause_button.grid(column=0, row=0, padx=5)

        self.ss_resume_button = ttk.Button(control_frame, text="Resume", command=self.resume_ss_animation)
        self.ss_resume_button.grid(column=1, row=0, padx=5)

        self.ss_speed_scale = ttk.Scale(
            control_frame, from_=0.1, to=5.0, value=1.0, orient='horizontal', command=self.update_ss_speed
        )
        self.ss_speed_scale.grid(column=2, row=0, padx=5)
        ttk.Label(control_frame, text="Speed").grid(column=2, row=1)

    def simulate_sonication(self):
        try:
            nx = int(self.ss_nx_entry.get())
            nt = int(self.ss_nt_entry.get())
            if nx <= 0 or nt <= 0:
                raise ValueError("Number of spatial and time steps must be positive integers.")
            dx = 0.01  # Spatial step (m)
            c_sound = 1500  # Speed of sound in m/s (adjust as necessary)
            dt = dx / (2 * c_sound)  # Time step (s), ensuring stability with CFL condition
            source_position = nx // 2

            # Run the simulation and get pressure history
            pressure_history = sonication_simulation(nx, nt, dx, dt, source_position, c=c_sound)
            x = np.arange(nx) * dx

            # Set up the plot for animation
            self.ss_fig.clf()
            self.ss_ax = self.ss_fig.add_subplot(111)
            self.ss_ax.set_xlabel('Position [m]')
            self.ss_ax.set_ylabel('Pressure')
            self.ss_ax.set_title('Pressure Field during Sonication Simulation')
            self.ss_ax.set_xlim(0, nx * dx)
            pressure_min = np.min(pressure_history)
            pressure_max = np.max(pressure_history)
            # Avoid zero range
            if pressure_max - pressure_min == 0:
                pressure_max += 1
            self.ss_ax.set_ylim(pressure_min, pressure_max)
            self.ss_line, = self.ss_ax.plot(x, pressure_history[0], color='b')

            # Animation variables
            self.ss_pause = False
            self.ss_animation_speed = 1.0  # Default speed factor

            # Animation update function
            def update_ss(frame):
                if self.ss_pause:
                    return self.ss_line,
                if frame >= nt:
                    return self.ss_line,
                self.ss_line.set_ydata(pressure_history[frame])
                self.ss_ax.set_title(f'Pressure Field at Time Step {frame+1}/{nt}')
                return self.ss_line,

            # Create the animation
            self.ss_anim = FuncAnimation(
                self.ss_fig,
                update_ss,
                frames=nt,
                interval=50 / self.ss_animation_speed,  # Adjust interval based on speed
                blit=True,
                repeat=False
            )

            self.ss_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")

    def pause_ss_animation(self):
        self.ss_pause = True

    def resume_ss_animation(self):
        self.ss_pause = False

    def update_ss_speed(self, val):
        try:
            speed = float(val)
            if speed <= 0:
                raise ValueError("Speed must be positive.")
            self.ss_animation_speed = speed
            if hasattr(self, 'ss_anim'):
                self.ss_anim.event_source.interval = 50 / self.ss_animation_speed
        except ValueError:
            self.ss_animation_speed = 1.0  # Reset to default if invalid

    def create_neutron_capture_tab(self):
        # Neutron Capture and Decay Tab
        frame = self.tab9

        ttk.Label(frame, text="Time Duration (s):").grid(column=0, row=0, padx=5, pady=5)
        self.nc_time_entry = ttk.Entry(frame)
        self.nc_time_entry.insert(0, "100")
        self.nc_time_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Label(frame, text="Neutron Density (1/m³):").grid(column=0, row=1, padx=5, pady=5)
        self.nc_density_entry = ttk.Entry(frame)
        self.nc_density_entry.insert(0, "1e28")
        self.nc_density_entry.grid(column=1, row=1, padx=5, pady=5)

        ttk.Label(frame, text="Capture Rate (1/s):").grid(column=0, row=2, padx=5, pady=5)
        self.nc_rate_entry = ttk.Entry(frame)
        self.nc_rate_entry.insert(0, "1e-3")
        self.nc_rate_entry.grid(column=1, row=2, padx=5, pady=5)

        ttk.Button(frame, text="Simulate", command=self.simulate_neutron_capture).grid(
            column=0, row=3, padx=5, pady=5, columnspan=2
        )

        # Plot area
        self.nc_fig = Figure(figsize=(6, 5), dpi=100)
        self.nc_canvas = FigureCanvasTkAgg(self.nc_fig, master=frame)
        self.nc_canvas.get_tk_widget().grid(column=0, row=4, columnspan=2)

    def simulate_neutron_capture(self):
        try:
            duration = float(self.nc_time_entry.get())
            n_density = float(self.nc_density_entry.get())
            capture_rate = float(self.nc_rate_entry.get())
            if duration <= 0 or n_density < 0 or capture_rate < 0:
                raise ValueError("Time duration must be positive and densities/rates cannot be negative.")

            y0 = [1e23, 0, 0, 0]  # Initial concentrations: N, D, He, Gamma photons
            t = np.linspace(0, duration, 1000)  # Time array

            sol = odeint(neutron_capture, y0, t, args=(n_density, capture_rate))

            # Plotting
            self.nc_fig.clf()
            ax = self.nc_fig.add_subplot(111)
            ax.plot(t, sol[:, 0], label='Neutron (N)')
            ax.plot(t, sol[:, 1], label='Deuterium (D)')
            ax.plot(t, sol[:, 2], label='Helium-3 (He-3)')
            ax.plot(t, sol[:, 3], label='Gamma photons')
            ax.set_xlabel('Time [s]')
            ax.set_ylabel('Concentration [1/m³]')
            ax.legend()
            ax.set_title('Neutron Capture and Decay')
            self.nc_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")
        except Exception as e:
            messagebox.showerror("Error", f"An unexpected error occurred: {e}")

    def create_bubble_dynamics_tab(self):
        # Bubble Dynamics Tab
        frame = self.tab10

        ttk.Label(frame, text="Time Duration (s):").grid(column=0, row=0, padx=5, pady=5)
        self.bd_time_entry = ttk.Entry(frame)
        self.bd_time_entry.insert(0, "0.01")
        self.bd_time_entry.grid(column=1, row=0, padx=5, pady=5)

        ttk.Button(frame, text="Simulate", command=self.simulate_bubble_dynamics).grid(
            column=0, row=1, padx=5, pady=5, columnspan=2
        )

        # Plot area
        self.bd_fig = Figure(figsize=(6, 5), dpi=100)
        self.bd_canvas = FigureCanvasTkAgg(self.bd_fig, master=frame)
        self.bd_canvas.get_tk_widget().grid(column=0, row=2, columnspan=2)

    def simulate_bubble_dynamics(self):
        try:
            duration = float(self.bd_time_entry.get())
            if duration <= 0:
                raise ValueError("Time duration must be positive.")
            R0 = [1e-6, 0]  # Initial radius and velocity (meters, meters/second)
            t_span = (0, duration)
            t_eval = np.linspace(*t_span, 1000)
            P0 = 101325  # Ambient pressure (Pa)
            Pv = 2338  # Vapor pressure of water at 20°C (Pa)
            rho = 1000  # Density of water (kg/m^3)
            gamma_val = 1.4  # Polytropic index
            mu = 0.001  # Viscosity of water (Pa·s)
            sigma = 0.0728  # Surface tension of water (N/m)

            sol = solve_ivp(
                rayleigh_plesset, t_span, R0, args=(P0, Pv, rho, gamma_val, mu, sigma), t_eval=t_eval, method='RK45'
            )

            # Plotting
            self.bd_fig.clf()
            ax = self.bd_fig.add_subplot(111)
            ax.plot(sol.t, sol.y[0])
            ax.set_xlabel('Time [s]')
            ax.set_ylabel('Bubble Radius [m]')
            ax.set_title('Bubble Dynamics under Sonication')
            ax.grid(True)
            self.bd_canvas.draw()
        except ValueError as ve:
            messagebox.showerror("Invalid input", f"Invalid input: {ve}")
        except Exception as e:
            messagebox.showerror("Error", f"An error occurred during simulation: {e}")

def main():
    root = tk.Tk()
    app = LENRSimulationApp(root)
    root.mainloop()

if __name__ == "__main__":
    main() the required things are defined in the gui tabs so remove the unwanted code parts automatically