""" TLT-002: 2D Carbon Lattice Simulation ====================================== First-principles 2D interference pattern generation. Variables: f (frequency), t (decoherence parameter) No other parameters. No nudging. The physics determines the output. Layer 1: FDTD wave equation solver (single pulsed source) Layer 2: N-wave plane wave interference (all 2D symmetries) Calibration target: graphene (a = 2.46 Angstroms, honeycomb, P6/mmm) """ import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib.colors import Normalize from pathlib import Path import json # ============================================================================ # CONSTANTS (exact values, no approximations) # ============================================================================ C_LIGHT = 299_792_458.0 # m/s (exact, SI definition) H_PLANCK = 6.62607015e-34 # J·s (exact, 2019 SI) CARBON_AMU = 12.011 # amu (standard atomic weight) AMU_KG = 1.66053906660e-27 # kg per amu CARBON_MASS = CARBON_AMU * AMU_KG CARBON_COMPTON_FREQ = CARBON_MASS * C_LIGHT**2 / H_PLANCK # Hz CARBON_COMPTON_WAVELENGTH = C_LIGHT / CARBON_COMPTON_FREQ # m GRAPHENE_LATTICE_CONSTANT = 2.46e-10 # m (2.46 Angstroms) OUTPUT_DIR = Path("/root/rhetroiluso/project_prometheus/time_ledger_theory/" "tlt results/unaudited/2D_lattice_images") OUTPUT_DIR.mkdir(parents=True, exist_ok=True) def print_constants(): """Print all physical constants used.""" print("=" * 70) print("PHYSICAL CONSTANTS") print("=" * 70) print(f" Speed of light c = {C_LIGHT:.0f} m/s") print(f" Planck constant h = {H_PLANCK:.5e} J·s") print(f" Carbon mass = {CARBON_MASS:.5e} kg") print(f" Carbon Compton freq = {CARBON_COMPTON_FREQ:.4e} Hz") print(f" Carbon Compton wvl = {CARBON_COMPTON_WAVELENGTH:.4e} m") print(f" Graphene lattice a = {GRAPHENE_LATTICE_CONSTANT:.2e} m") print(f" Ratio a/lambda_C = {GRAPHENE_LATTICE_CONSTANT/CARBON_COMPTON_WAVELENGTH:.4e}") print("=" * 70) # ============================================================================ # LAYER 2: N-WAVE PLANE WAVE INTERFERENCE # ============================================================================ # This is analytically exact. No numerical approximation. # The interference of N coherent plane waves at equal angular spacing # produces specific 2D patterns determined ONLY by N and wavelength. # # psi(x,y) = sum_{n=0}^{N-1} exp(i * k * (x*cos(theta_n) + y*sin(theta_n))) # I(x,y) = |psi(x,y)|^2 # # where theta_n = 2*pi*n/N, k = 2*pi/lambda # ============================================================================ def generate_nwave_interference(N, wavelength, grid_extent, resolution=1000): """ Generate 2D interference pattern from N plane waves at equal angles. Parameters: N: number of plane waves (2, 3, 4, 6, etc.) wavelength: wavelength of each wave (sets scale) grid_extent: half-width of the computation domain (in same units as wavelength) resolution: number of grid points per axis Returns: X, Y: coordinate meshgrids intensity: |psi|^2 at each point psi: complex wave field """ x = np.linspace(-grid_extent, grid_extent, resolution) y = np.linspace(-grid_extent, grid_extent, resolution) X, Y = np.meshgrid(x, y) k = 2.0 * np.pi / wavelength psi = np.zeros_like(X, dtype=complex) for n in range(N): theta = 2.0 * np.pi * n / N kx = k * np.cos(theta) ky = k * np.sin(theta) psi += np.exp(1j * (kx * X + ky * Y)) intensity = np.abs(psi) ** 2 return X, Y, intensity, psi def plot_interference_pattern(X, Y, intensity, N, title_extra="", filename=None, cmap='inferno'): """Plot and save a 2D interference pattern as PNG.""" fig, ax = plt.subplots(1, 1, figsize=(10, 10)) im = ax.pcolormesh(X, Y, intensity, cmap=cmap, shading='auto') cbar = fig.colorbar(im, ax=ax, shrink=0.8, label='Intensity |ψ|²') ax.set_xlabel('x / λ', fontsize=14) ax.set_ylabel('y / λ', fontsize=14) ax.set_title(f'{N}-Wave Interference Pattern{title_extra}', fontsize=16) ax.set_aspect('equal') ax.tick_params(labelsize=12) plt.tight_layout() if filename: filepath = OUTPUT_DIR / filename fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) def run_layer2_all_symmetries(): """ Generate interference patterns for all natural 2D symmetries. N=2 (stripes), N=3 (hexagonal), N=4 (square), N=6 (hexagonal alt). Uses normalized wavelength = 1.0. """ print("\n" + "=" * 70) print("LAYER 2: N-WAVE PLANE WAVE INTERFERENCE") print("Testing all natural 2D symmetries: N = 2, 3, 4, 6") print("Wavelength = 1.0 (normalized)") print("=" * 70) wavelength = 1.0 grid_extent = 5.0 # show ~10 wavelengths across resolution = 1200 symmetry_names = { 2: "stripes (1D standing wave)", 3: "hexagonal (120° symmetry)", 4: "square (90° symmetry)", 6: "hexagonal alternate (60° symmetry)" } results = {} for N in [2, 3, 4, 6]: print(f"\n Computing N={N}: {symmetry_names[N]}...") X, Y, intensity, psi = generate_nwave_interference( N, wavelength, grid_extent, resolution ) # Measure the pattern characteristics I_max = np.max(intensity) I_min = np.min(intensity) I_range = I_max - I_min # Find peak positions along x-axis (y=0) to measure spacing mid_row = resolution // 2 x_slice = X[mid_row, :] i_slice = intensity[mid_row, :] # Find peaks (local maxima) from scipy.signal import find_peaks peaks, properties = find_peaks(i_slice, height=I_max * 0.8, distance=resolution // 50) if len(peaks) >= 2: peak_positions = x_slice[peaks] spacings = np.diff(peak_positions) avg_spacing = np.mean(spacings) else: avg_spacing = None results[N] = { 'I_max': float(I_max), 'I_min': float(I_min), 'avg_spacing_normalized': float(avg_spacing) if avg_spacing else None, 'num_peaks_found': len(peaks), 'symmetry': symmetry_names[N] } print(f" Intensity range: [{I_min:.2f}, {I_max:.2f}]") print(f" Peaks found on x-axis: {len(peaks)}") if avg_spacing: print(f" Avg peak spacing: {avg_spacing:.4f} λ") # Plot the interference pattern filename = f"carbon_N{N}_{symmetry_names[N].split()[0]}.png" plot_interference_pattern( X, Y, intensity, N, title_extra=f"\n{symmetry_names[N]} — λ = 1.0 (normalized)", filename=filename ) # Also generate a zoomed-in view showing the unit cell fig, ax = plt.subplots(1, 1, figsize=(10, 10)) zoom = 2.0 # show fewer wavelengths for detail mask_x = (X[0, :] >= -zoom) & (X[0, :] <= zoom) mask_y = (Y[:, 0] >= -zoom) & (Y[:, 0] <= zoom) X_z = X[np.ix_(mask_y, mask_x)] Y_z = Y[np.ix_(mask_y, mask_x)] I_z = intensity[np.ix_(mask_y, mask_x)] im = ax.pcolormesh(X_z, Y_z, I_z, cmap='inferno', shading='auto') cbar = fig.colorbar(im, ax=ax, shrink=0.8, label='Intensity |ψ|²') ax.set_xlabel('x / λ', fontsize=14) ax.set_ylabel('y / λ', fontsize=14) ax.set_title(f'{N}-Wave Interference (zoomed)\n{symmetry_names[N]}', fontsize=16) ax.set_aspect('equal') ax.tick_params(labelsize=12) plt.tight_layout() filepath = OUTPUT_DIR / f"carbon_N{N}_{symmetry_names[N].split()[0]}_zoom.png" fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) return results # ============================================================================ # LAYER 1: FDTD 2D WAVE EQUATION SOLVER # ============================================================================ # Solves: d^2(psi)/dt^2 = c^2 * (d^2(psi)/dx^2 + d^2(psi)/dy^2) # using finite-difference time-domain method. # # Source: single point, pulsing at frequency f with decoherence gaps. # Boundary: absorbing (Mur first-order ABC). # Output: accumulated intensity and instantaneous snapshots. # ============================================================================ def fdtd_2d_pulsed(f_norm, t_decohere_norm, n_cycles=30, grid_size=200, courant=0.5): """ 2D FDTD wave equation solver with pulsed source. Uses normalized units: - Spatial unit = 1 (grid spacing dx = 1) - Time unit = dt = courant * dx / (c_norm * sqrt(2)) - c_norm = 1 (normalized speed) - f_norm = normalized frequency - t_decohere_norm = decoherence period in normalized time units Parameters: f_norm: normalized source frequency (cycles per time unit) t_decohere_norm: decoherence gap duration (in time units) n_cycles: number of pulse-rest cycles to simulate grid_size: number of grid points per axis courant: Courant number (stability requires <= 1/sqrt(2) ~ 0.707) Returns: intensity: time-averaged |psi|^2 snapshot: final instantaneous wavefield x, y: coordinate arrays """ nx = ny = grid_size dx = dy = 1.0 c_norm = 1.0 dt = courant * dx / (c_norm * np.sqrt(2.0)) # Wave field arrays psi = np.zeros((nx, ny)) psi_prev = np.zeros((nx, ny)) psi_old_boundary = np.zeros((nx, ny)) # for Mur ABC intensity_acc = np.zeros((nx, ny)) n_steps_acc = 0 # Source position (center) sx, sy = nx // 2, ny // 2 # Courant factor squared r2 = (c_norm * dt / dx) ** 2 # Time parameters period = 1.0 / f_norm if f_norm > 0 else 1.0 cycle_duration = period + t_decohere_norm total_time = n_cycles * cycle_duration n_steps = int(total_time / dt) print(f" FDTD parameters:") print(f" Grid: {nx}x{ny}, dx={dx}") print(f" dt={dt:.4e}, Courant={courant}") print(f" f_norm={f_norm}, period={period:.4f}") print(f" t_decohere={t_decohere_norm:.4f}") print(f" Cycles={n_cycles}, total steps={n_steps}") for step in range(n_steps): t = step * dt # Save old field for Mur ABC psi_old_boundary[:] = psi_prev # Interior update: standard 2D wave equation stencil psi_next = np.zeros_like(psi) psi_next[1:-1, 1:-1] = ( 2.0 * psi[1:-1, 1:-1] - psi_prev[1:-1, 1:-1] + r2 * ( psi[2:, 1:-1] + psi[:-2, 1:-1] + psi[1:-1, 2:] + psi[1:-1, :-2] - 4.0 * psi[1:-1, 1:-1] ) ) # Pulsed source: active during pulse phase, off during decoherence cycle_pos = t % cycle_duration if cycle_pos < period: # Soft source (additive, doesn't overwrite field) psi_next[sx, sy] += np.sin(2.0 * np.pi * f_norm * t) # Mur first-order absorbing boundary conditions # Left boundary (x=0) psi_next[0, :] = psi_old_boundary[1, :] + ( (c_norm * dt - dx) / (c_norm * dt + dx) ) * (psi_next[1, :] - psi[0, :]) # Right boundary (x=nx-1) psi_next[-1, :] = psi_old_boundary[-2, :] + ( (c_norm * dt - dx) / (c_norm * dt + dx) ) * (psi_next[-2, :] - psi[-1, :]) # Bottom boundary (y=0) psi_next[:, 0] = psi_old_boundary[:, 1] + ( (c_norm * dt - dy) / (c_norm * dt + dy) ) * (psi_next[:, 1] - psi[:, 0]) # Top boundary (y=ny-1) psi_next[:, -1] = psi_old_boundary[:, -2] + ( (c_norm * dt - dy) / (c_norm * dt + dy) ) * (psi_next[:, -2] - psi[:, -1]) # Accumulate intensity (skip initial transient) if step > n_steps // 4: intensity_acc += psi_next ** 2 n_steps_acc += 1 # Advance time step psi_prev[:] = psi psi[:] = psi_next # Normalize accumulated intensity if n_steps_acc > 0: intensity_acc /= n_steps_acc x = np.arange(nx) * dx - (nx // 2) * dx y = np.arange(ny) * dy - (ny // 2) * dy return intensity_acc, psi, x, y def run_layer1_fdtd(): """ Run FDTD simulations with pulsed source at various decoherence parameters. """ print("\n" + "=" * 70) print("LAYER 1: FDTD WAVE EQUATION (single pulsed source)") print("Solving: d²ψ/dt² = c²∇²ψ") print("Source: single point, pulsing at frequency f") print("=" * 70) f_norm = 0.05 # normalized frequency (lower = longer wavelength, clearer patterns) # Test several decoherence parameters # t_d = 0 means continuous source (reference) # t_d > 0 introduces gaps between pulses decohere_values = [0.0, 0.5, 1.0, 2.0, 5.0, 10.0] results = {} for t_d in decohere_values: period = 1.0 / f_norm label = f"t_d={t_d:.1f}" if t_d > 0 else "continuous" ratio = t_d / period if period > 0 else 0 print(f"\n Running FDTD: {label} (t_d/T = {ratio:.2f})...") intensity, snapshot, x, y = fdtd_2d_pulsed( f_norm=f_norm, t_decohere_norm=t_d, n_cycles=40, grid_size=300, courant=0.5 ) results[label] = { 't_decohere': t_d, 'ratio_td_T': ratio, 'I_max': float(np.max(intensity)), 'I_mean': float(np.mean(intensity)) } # Plot accumulated intensity fig, axes = plt.subplots(1, 2, figsize=(20, 9)) # Intensity map X, Y = np.meshgrid(x, y) im0 = axes[0].pcolormesh(X, Y, intensity, cmap='inferno', shading='auto') fig.colorbar(im0, ax=axes[0], shrink=0.8, label='Time-averaged |ψ|²') axes[0].set_title(f'Accumulated Intensity\nf={f_norm}, {label}', fontsize=14) axes[0].set_xlabel('x', fontsize=12) axes[0].set_ylabel('y', fontsize=12) axes[0].set_aspect('equal') # Final snapshot vmax = max(abs(np.min(snapshot)), abs(np.max(snapshot))) im1 = axes[1].pcolormesh(X, Y, snapshot, cmap='RdBu_r', shading='auto', vmin=-vmax, vmax=vmax) fig.colorbar(im1, ax=axes[1], shrink=0.8, label='ψ (instantaneous)') axes[1].set_title(f'Instantaneous Wavefield\nf={f_norm}, {label}', fontsize=14) axes[1].set_xlabel('x', fontsize=12) axes[1].set_ylabel('y', fontsize=12) axes[1].set_aspect('equal') plt.tight_layout() safe_label = label.replace("=", "").replace(".", "p") filepath = OUTPUT_DIR / f"FDTD_carbon_{safe_label}.png" fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) # Radial profile analysis fig, ax = plt.subplots(figsize=(12, 6)) # Compute radial average center_x, center_y = len(x) // 2, len(y) // 2 R = np.sqrt((X - x[center_x])**2 + (Y - y[center_y])**2) r_bins = np.linspace(0, np.max(R) * 0.7, 200) r_centers = 0.5 * (r_bins[:-1] + r_bins[1:]) radial_mean = np.zeros(len(r_centers)) for i in range(len(r_centers)): mask = (R >= r_bins[i]) & (R < r_bins[i+1]) if np.any(mask): radial_mean[i] = np.mean(intensity[mask]) ax.plot(r_centers, radial_mean, 'k-', linewidth=1.5) ax.set_xlabel('Radial distance from source', fontsize=14) ax.set_ylabel('Radially-averaged intensity', fontsize=14) ax.set_title(f'Radial Profile — {label}\n' f'(radial symmetry = concentric rings; deviation = structure)', fontsize=14) ax.grid(True, alpha=0.3) plt.tight_layout() filepath = OUTPUT_DIR / f"FDTD_carbon_{safe_label}_radial.png" fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) return results # ============================================================================ # GRAPHENE OVERLAY COMPARISON # ============================================================================ def generate_graphene_overlay(): """ Generate the known graphene lattice for visual comparison. Uses the actual honeycomb geometry with lattice constant a = 1.0 (normalized to match the interference pattern scale). """ print("\n" + "=" * 70) print("GRAPHENE REFERENCE LATTICE") print("Generating known honeycomb geometry for comparison") print("=" * 70) a = 1.0 # normalized lattice constant # Honeycomb lattice vectors a1 = a * np.array([1.0, 0.0]) a2 = a * np.array([0.5, np.sqrt(3.0) / 2.0]) # Basis positions (two atoms per unit cell) basis = [ np.array([0.0, 0.0]), np.array([0.5, np.sqrt(3.0) / 6.0]) * a ] # Actually for honeycomb: basis vectors and atom positions # Lattice vectors for hexagonal: a1 = a * np.array([1.0, 0.0]) a2 = a * np.array([0.5, np.sqrt(3.0) / 2.0]) # Two atoms per unit cell: d1 = np.array([0.0, 0.0]) d2 = (a1 + a2) / 3.0 # Generate lattice points n_range = 8 atoms = [] for i in range(-n_range, n_range + 1): for j in range(-n_range, n_range + 1): R = i * a1 + j * a2 atoms.append(R + d1) atoms.append(R + d2) atoms = np.array(atoms) # Plot honeycomb fig, ax = plt.subplots(figsize=(10, 10)) # Draw bonds (connect atoms within bond distance) bond_length = a / np.sqrt(3.0) for i in range(len(atoms)): for j in range(i + 1, len(atoms)): dist = np.linalg.norm(atoms[i] - atoms[j]) if dist < bond_length * 1.1: ax.plot([atoms[i, 0], atoms[j, 0]], [atoms[i, 1], atoms[j, 1]], 'k-', linewidth=1.0, alpha=0.5) ax.scatter(atoms[:, 0], atoms[:, 1], c='black', s=30, zorder=5) ax.set_xlim(-4, 4) ax.set_ylim(-4, 4) ax.set_aspect('equal') ax.set_xlabel('x / a', fontsize=14) ax.set_ylabel('y / a', fontsize=14) ax.set_title('Graphene Honeycomb Lattice (reference)\n' 'a = 2.46 Å, space group P6/mmm', fontsize=16) ax.grid(True, alpha=0.2) plt.tight_layout() filepath = OUTPUT_DIR / "graphene_reference_lattice.png" fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) # Now overlay with N=3 interference pattern print(" Generating N=3 interference overlay with graphene...") # The N=3 interference pattern has its own characteristic spacing. # We set lambda so the interference maxima align with graphene atoms. # For 3 waves at 120°, the intensity maxima form a hexagonal lattice # with lattice constant a_hex = 2*lambda / 3. # To match graphene's a=1.0 (normalized): lambda = 3*a/2 = 1.5 # BUT: we don't assume this. We show both side by side and overlay. wavelength = 1.0 # keep normalized grid_extent = 4.0 resolution = 800 X, Y, intensity, _ = generate_nwave_interference( 3, wavelength, grid_extent, resolution ) fig, axes = plt.subplots(1, 3, figsize=(30, 9)) # Panel 1: N=3 interference pattern im0 = axes[0].pcolormesh(X, Y, intensity, cmap='inferno', shading='auto') fig.colorbar(im0, ax=axes[0], shrink=0.8) axes[0].set_title('N=3 Interference Pattern\n(3 plane waves at 120°)', fontsize=14) axes[0].set_aspect('equal') axes[0].set_xlabel('x / λ') axes[0].set_ylabel('y / λ') # Panel 2: Graphene lattice for i in range(len(atoms)): for j in range(i + 1, len(atoms)): dist = np.linalg.norm(atoms[i] - atoms[j]) if dist < bond_length * 1.1: axes[1].plot([atoms[i, 0], atoms[j, 0]], [atoms[i, 1], atoms[j, 1]], 'k-', linewidth=1.0, alpha=0.5) axes[1].scatter(atoms[:, 0], atoms[:, 1], c='black', s=30, zorder=5) axes[1].set_xlim(-4, 4) axes[1].set_ylim(-4, 4) axes[1].set_aspect('equal') axes[1].set_title('Graphene Honeycomb\n(a = 2.46 Å reference)', fontsize=14) axes[1].set_xlabel('x / a') axes[1].set_ylabel('y / a') # Panel 3: Overlay — N=3 pattern with graphene atoms # Scale the interference pattern to match graphene lattice constant # N=3 hexagonal maxima spacing: measure from the pattern # For 3 waves, maxima positions form hexagonal lattice im2 = axes[2].pcolormesh(X, Y, intensity, cmap='inferno', shading='auto', alpha=0.8) # Overlay graphene atoms scaled to interference pattern coordinates # The N=3 pattern maxima are at triangular lattice points with # spacing a_tri = 2*lambda / 3 a_tri = 2.0 * wavelength / 3.0 scale = a_tri # scale graphene coords to match atoms_scaled = atoms * scale mask = (np.abs(atoms_scaled[:, 0]) < grid_extent) & \ (np.abs(atoms_scaled[:, 1]) < grid_extent) axes[2].scatter(atoms_scaled[mask, 0], atoms_scaled[mask, 1], c='cyan', s=40, zorder=5, edgecolors='white', linewidths=0.5, alpha=0.9) axes[2].set_xlim(-grid_extent, grid_extent) axes[2].set_ylim(-grid_extent, grid_extent) axes[2].set_aspect('equal') axes[2].set_title('Overlay: N=3 pattern + graphene atoms\n' '(cyan dots = graphene atom positions)', fontsize=14) axes[2].set_xlabel('x / λ') axes[2].set_ylabel('y / λ') fig.colorbar(im2, ax=axes[2], shrink=0.8) plt.tight_layout() filepath = OUTPUT_DIR / "carbon_N3_vs_graphene_overlay.png" fig.savefig(filepath, dpi=150, bbox_inches='tight') print(f" Saved: {filepath}") plt.close(fig) return atoms # ============================================================================ # MAIN # ============================================================================ def main(): print_constants() # Layer 2: All symmetries (fast, analytically exact) layer2_results = run_layer2_all_symmetries() # Graphene reference and overlay generate_graphene_overlay() # Layer 1: FDTD (computationally intensive) layer1_results = run_layer1_fdtd() # Save all numerical results all_results = { 'constants': { 'carbon_compton_freq_Hz': CARBON_COMPTON_FREQ, 'carbon_compton_wavelength_m': CARBON_COMPTON_WAVELENGTH, 'graphene_lattice_constant_m': GRAPHENE_LATTICE_CONSTANT, 'ratio_a_over_lambda_C': GRAPHENE_LATTICE_CONSTANT / CARBON_COMPTON_WAVELENGTH, }, 'layer2_nwave': {str(k): v for k, v in layer2_results.items()}, 'layer1_fdtd': layer1_results, } data_path = OUTPUT_DIR / "carbon_2D_lattice_data.txt" with open(data_path, 'w') as f: f.write("=" * 70 + "\n") f.write("TLT-002: 2D CARBON LATTICE TEST — NUMERICAL RESULTS\n") f.write("=" * 70 + "\n\n") f.write("PHYSICAL CONSTANTS\n") f.write("-" * 40 + "\n") for k, v in all_results['constants'].items(): f.write(f" {k}: {v}\n") f.write("\n\nLAYER 2: N-WAVE INTERFERENCE RESULTS\n") f.write("-" * 40 + "\n") for n_key, data in all_results['layer2_nwave'].items(): f.write(f"\n N = {n_key}:\n") for k, v in data.items(): f.write(f" {k}: {v}\n") f.write("\n\nLAYER 1: FDTD RESULTS\n") f.write("-" * 40 + "\n") for label, data in all_results['layer1_fdtd'].items(): f.write(f"\n {label}:\n") for k, v in data.items(): f.write(f" {k}: {v}\n") print(f"\n Data saved: {data_path}") print("\n" + "=" * 70) print("SIMULATION COMPLETE") print(f"All outputs in: {OUTPUT_DIR}") print("=" * 70) if __name__ == "__main__": main()