.

vec_temp.py

@@ -1,168 +1,209 @@
 import sympy as sp
 import numpy as np
 import matplotlib.pyplot as plt
-from mpl_toolkits.mplot3d import Axes3D
-from matplotlib.widgets import Slider
+from scipy.integrate import dblquad
 
-# ==============================================================================
-# Part 1: Symbolic Calculations (Expanded)
-# ==============================================================================
+# --- RICH IMPORTS ---
+from rich.progress import Progress, BarColumn, TextColumn, TimeRemainingColumn, TimeElapsedColumn, track
+
+# --- NEW: Import joblib for parallel processing ---
+from joblib import Parallel, delayed, cpu_count
+
+import warnings
+
+# Use this in a Jupyter/IPython environment for nicely formatted math output
+sp.init_printing(use_unicode=True) 
+
+# Reduced for faster testing, feel free to increase it back to 100_000
+n_samples = 20_000 
+
+# --- Part 1: Symbolic Definitions (Unchanged) ---
 
-# 1. Define symbolic variables
 u, v, theta, a, b, gamma = sp.symbols('u v theta a b gamma')
 pi = sp.pi
 
-# 2. Define the Manifold Parameterization
-X = (a + gamma * sp.sin(u/4)) * sp.cos(v/2) + (b * sp.cos(theta) * sp.cos(u/4))
-Y = b * sp.cos(theta) * sp.sin(u/4) + b * sp.sin(theta) * ((1 + sp.cos(v/2))/2) * \
-    ((v/(2*pi))*(sp.cos(u/4) - 1) - u/(2*pi) + (u*v)/(4*pi**2) + 1)
-Z = b * sp.sin(theta) * sp.sin(u/4) + b * sp.sin(theta) * ((1 - sp.cos(v/2))/2) * \
-    ((v/(2*pi))*(1 - sp.cos(u/4)) + sp.cos(u/4) - (u*v)/(4*pi**2))
-
-# Position vector
-r = sp.Matrix([X, Y, Z])
-
-# 3. Define the Vector Field
-# Partial derivatives of the position vector
-r_v = sp.diff(r, v)
-r_theta = sp.diff(r, theta)
-# Sum of component-wise derivatives
-B_vec = sp.Matrix([
-    sp.diff(X, v) + sp.diff(X, theta),
-    sp.diff(Y, v) + sp.diff(Y, theta),
-    sp.diff(Z, v) + sp.diff(Z, theta)
+# Position Vector R
+R_vec = sp.Matrix([
+    (a + gamma * sp.sin(u/4)) * sp.cos(v/2) + (b * sp.cos(theta) * sp.cos(u/4)),
+    b * sp.cos(theta) * sp.sin(u/4) + b * sp.sin(theta) * ((1 + sp.cos(v/2))/2) * \
+        ((v/(2*pi))*(sp.cos(u/4) - 1) - u/(2*pi) + (u*v)/(4*pi**2) + 1),
+    b * sp.sin(theta) * sp.sin(u/4) + b * sp.sin(theta) * ((1 - sp.cos(v/2))/2) * \
+        ((v/(2*pi))*(1 - sp.cos(u/4)) + sp.cos(u/4) - (u*v)/(4*pi**2))
 ])
-B_mag = B_vec.norm() # Magnitude of the vector field ||B||
-
-# 4. Calculate the Surface Area Element
-# The cross product gives a vector normal to the surface
-cross_prod = r_v.cross(r_theta)
-# The magnitude of this vector is the differential area element dS
-dS_mag = cross_prod.norm()
-
-# 5. Define the full integrand for our numerical integral
-# We will integrate the magnitude of the vector field over the surface area
-integrand_expr = B_mag * dS_mag
-
-# ==============================================================================
-# Part 2: Numerical Conversion and Interactive Visualization
-# ==============================================================================
-
-# 1. Lambdify all necessary expressions for speed
-args = [u, v, theta, a, b, gamma]
-X_func = sp.lambdify(args, X, 'numpy')
-Y_func = sp.lambdify(args, Y, 'numpy')
-Z_func = sp.lambdify(args, Z, 'numpy')
-integrand_func = sp.lambdify(args, integrand_expr, 'numpy')
-
-# 2. Define constants and parameter ranges
-A_val = 5.0
-B_val = 2.0
-GAMMA_val = 1.5
-u_plot_range = np.linspace(0, 2, 40) # Range of u for the 2D plot
-
-# Create high-resolution grids for calculation and plotting
-v_res, theta_res = 100, 100
-v_range = np.linspace(0, 2 * np.pi, v_res)
-theta_range = np.linspace(0, 2 * np.pi, theta_res)
-v_grid, theta_grid = np.meshgrid(v_range, theta_range)
-dv = v_range[1] - v_range[0]
-dtheta = theta_range[1] - theta_range[0]
-dA = dv * dtheta # Area of one grid cell for integration
-
-# 3. Set up the plotting figure and axes
-fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8), width_ratios=[1.2, 1])
-ax1 = plt.subplot(1, 2, 1, projection='3d')
-ax2 = plt.subplot(1, 2, 2)
-fig.subplots_adjust(bottom=0.25) # Make space for sliders
-
-# 4. Helper function for numerical integration
-def calculate_integral(u_val, R_val):
-    """Calculates the integral for a specific u and R."""
-    # Calculate manifold coordinates
-    X_surf = X_func(u_val, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    Y_surf = Y_func(u_val, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    Z_surf = Z_func(u_val, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    
-    # Condition: find which points are inside the sphere
-    is_inside = (X_surf**2 + Y_surf**2 + Z_surf**2) < R_val**2
-    
-    # Calculate integrand values across the whole surface
-    integrand_values = integrand_func(u_val, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    
-    # Apply the condition: only consider values inside the sphere
-    integrand_masked = np.where(is_inside, integrand_values, 0)
+
+# Vector Field B
+B_vec = R_vec.diff(v) + R_vec.diff(theta)
+
+print("Successfully defined symbolic Manifold and Vector Field B.")
+
+# --- Part 2: Lambdify (Unchanged) ---
+
+print("Converting to numerical functions...")
+
+params = [u, v, theta, a, b, gamma]
+
+f_R = sp.lambdify(params, R_vec, modules=[{'ImmutableDenseMatrix': np.array}, 'numpy'])
+f_B = sp.lambdify(params, B_vec, modules=[{'ImmutableDenseMatrix': np.array}, 'numpy'])
+
+print("Lambdification complete.")
+
+# --- Part 3: Optimized Integrand (Unchanged) ---
+
+a_val, b_val, gamma_val = 1.0, 1.5, 6
+JACOBIAN_TOL = 1e-8
+EPSILON = 1e-5
+
+def integrand_vectorized(v_vals, theta_vals, u_val):
+    v_vals = np.atleast_1d(v_vals)
+    theta_vals = np.atleast_1d(theta_vals)
+    results = np.zeros_like(v_vals)
     
-    # Sum up the values and multiply by the differential area
-    integral_result = np.sum(integrand_masked) * dA
-    return integral_result
-
-# 5. Define the update functions for the plots
-def update_plots(val):
-    """Updates both plots. Called by sliders."""
-    u_current = u_slider.val
-    R_current = R_slider.val
+    for i, (v_val, theta_val) in enumerate(zip(v_vals, theta_vals)):
+        args = (u_val, v_val, theta_val, a_val, b_val, gamma_val)
+        try:
+            B_center = np.array(f_B(*args)).flatten()
+            R_center = np.array(f_R(*args)).flatten()
+            args_u = (u_val + EPSILON, v_val, theta_val, a_val, b_val, gamma_val)
+            args_v = (u_val, v_val + EPSILON, theta_val, a_val, b_val, gamma_val)
+            args_t = (u_val, v_val, theta_val + EPSILON, a_val, b_val, gamma_val)
+            B_u = np.array(f_B(*args_u)).flatten()
+            B_v = np.array(f_B(*args_v)).flatten()
+            B_t = np.array(f_B(*args_t)).flatten()
+            R_u = np.array(f_R(*args_u)).flatten()
+            R_v = np.array(f_R(*args_v)).flatten()
+            R_t = np.array(f_R(*args_t)).flatten()
+            grad_p_B = np.column_stack([(B_u - B_center) / EPSILON, (B_v - B_center) / EPSILON, (B_t - B_center) / EPSILON])
+            J_val = np.column_stack([(R_u - R_center) / EPSILON, (R_v - R_center) / EPSILON, (R_t - R_center) / EPSILON])
+            det = np.linalg.det(J_val)
+            if abs(det) < JACOBIAN_TOL:
+                results[i] = 0.0
+                continue
+            J_inv = np.linalg.inv(J_val)
+            grad_xyz_B = J_inv.T @ grad_p_B
+            div_B = np.trace(grad_xyz_B)
+            cross_prod = np.cross(J_val[:, 1], J_val[:, 2])
+            norm_val = np.linalg.norm(cross_prod)
+            if norm_val < JACOBIAN_TOL:
+                results[i] = 0.0
+                continue
+            result = div_B * norm_val
+            results[i] = result if np.isfinite(result) else 0.0
+        except Exception:
+            results[i] = 0.0
+    return results[0] if len(results) == 1 else results
+
+def integrand(v_val, theta_val, u_val):
+    return integrand_vectorized(v_val, theta_val, u_val)
+
+# --- Part 4: Integration with Parallel Monte Carlo ---
+
+print("\nStarting integration...")
+
+u_values = np.linspace(0, 2, 100)
+
+# --- NEW: Define a worker function for parallel execution ---
+# This function encapsulates the work for a *single* u_val.
+# It takes all necessary parameters and returns the result tuple.
+def worker_integrate(u_val, n_samples):
+    """
+    Performs Monte Carlo integration for a single u_val.
+    This function is designed to be called in parallel.
+    """
+    # Seed the random number generator differently for each process
+    # to ensure they don't all produce the same random numbers.
+    np.random.seed(int(u_val * 1000) + np.random.randint(0, 1000))
     
-    # --- Update 3D Plot (ax1) ---
-    ax1.clear()
+    v_samples = np.random.uniform(0, 2 * np.pi, n_samples)
+    theta_samples = np.random.uniform(0, 2 * np.pi, n_samples)
     
-    # Calculate manifold coordinates for the current u
-    X_surf = X_func(u_current, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    Y_surf = Y_func(u_current, v_grid, theta_grid, A_val, B_val, GAMMA_val)
-    Z_surf = Z_func(u_current, v_grid, theta_grid, A_val, B_val, GAMMA_val)
+    # Using a list comprehension is slightly cleaner here
+    values = [integrand(v, t, u_val) for v, t in zip(v_samples, theta_samples)]
     
-    # Plot the full manifold semi-transparently
-    ax1.plot_surface(X_surf, Y_surf, Z_surf, color='c', alpha=0.2, rstride=5, cstride=5)
+    # Filter out any non-finite results (inf, NaN)
+    finite_values = [v for v in values if np.isfinite(v)]
     
-    # Highlight the part inside the sphere
-    is_inside = (X_surf**2 + Y_surf**2 + Z_surf**2) < R_current**2
-    X_inside = np.where(is_inside, X_surf, np.nan) # Set points outside to NaN so they are not plotted
-    Y_inside = np.where(is_inside, Y_surf, np.nan)
-    Z_inside = np.where(is_inside, Z_surf, np.nan)
-    ax1.plot_surface(X_inside, Y_inside, Z_inside, alpha=0.8, rstride=5, cstride=5)
-
-    # Draw the sphere wireframe
-    phi, psi = np.mgrid[0:2*np.pi:10j, 0:np.pi:10j]
-    x_sphere = R_current * np.sin(psi) * np.cos(phi)
-    y_sphere = R_current * np.sin(psi) * np.sin(phi)
-    z_sphere = R_current * np.cos(psi)
-    ax1.plot_wireframe(x_sphere, y_sphere, z_sphere, color="r", alpha=0.3)
+    if len(finite_values) == 0:
+        return 0.0, 0.0
     
-    # Aesthetics for 3D plot
-    ax1.set_title(f"Manifold at u={u_current:.2f}, Sphere R={R_current:.2f}")
-    ax1.set_xlabel("X"); ax1.set_ylabel("Y"); ax1.set_zlabel("Z")
-    ax1.set_xlim([-8, 8]); ax1.set_ylim([-8, 8]); ax1.set_zlim([-4, 4])
-
-    # --- Update 2D Plot (ax2) ---
-    # This only needs to re-calculate if the radius changes, but it's simpler to update always.
-    print(f"Recalculating integral curve for R = {R_current:.2f}...")
-    integral_values = [calculate_integral(u_val, R_current) for u_val in u_plot_range]
+    # Monte Carlo estimate: (volume) * (average value)
+    volume = (2 * np.pi) * (2 * np.pi)
+    mean_val = np.mean(finite_values)
+    std_val = np.std(finite_values)
     
-    ax2.clear()
-    ax2.plot(u_plot_range, integral_values, marker='o', linestyle='-')
-    ax2.axvline(u_current, color='r', linestyle='--', label=f'Current u={u_current:.2f}')
-    ax2.set_title(f"Integral vs. u (for Sphere Radius R={R_current:.2f})")
-    ax2.set_xlabel("u")
-    ax2.set_ylabel("Numerical Integral of ||B|| over surface inside sphere")
-    ax2.grid(True)
-    ax2.legend()
-    ax2.set_xlim(0, 2)
+    integral = volume * mean_val
+    error = volume * std_val / np.sqrt(len(finite_values))
     
-    fig.canvas.draw_idle()
-
-# 6. Create the slider widgets
-ax_u_slider = fig.add_axes([0.2, 0.1, 0.65, 0.03])
-u_slider = Slider(ax=ax_u_slider, label='u value', valmin=0.0, valmax=2.0, valinit=0.0)
-
-ax_R_slider = fig.add_axes([0.2, 0.05, 0.65, 0.03])
-R_slider = Slider(ax=ax_R_slider, label='Sphere Radius R', valmin=0.1, valmax=8.0, valinit=4.0)
-
-# 7. Register the update function to the sliders
-u_slider.on_changed(update_plots)
-R_slider.on_changed(update_plots)
-
-# 8. Draw the initial plots
-update_plots(None)
-
-# 9. Show the plot window
+    return integral, error
+
+# --- MODIFIED: The main integration loop is now parallel ---
+# Use n_jobs=-1 to use all available CPU cores.
+n_jobs = -1 
+num_cores = cpu_count()
+print(f"Using parallel Monte Carlo integration on {num_cores} cores ({n_samples} samples per u)...")
+
+# Create a list of "delayed" tasks, one for each u_value.
+# 'delayed' is a wrapper that captures the function and its arguments.
+tasks = [delayed(worker_integrate)(u, n_samples) for u in u_values]
+
+# Execute the tasks in parallel. rich.progress.track provides the progress bar.
+# Parallel(...) returns a list of the results from each worker call.
+parallel_results = Parallel(n_jobs=n_jobs)(
+    track(tasks, description="[bold cyan]Integrating over u values")
+)
+
+# Unpack the list of (result, error) tuples into two separate lists.
+integral_results, integral_errors = zip(*parallel_results)
+
+print("\nIntegration finished.")
+print(f"Average integration error: {np.mean(integral_errors):.2e}")
+print(f"Max integration error: {np.max(integral_errors):.2e}")
+
+
+# --- Part 5: Plotting (Unchanged) ---
+
+plt.figure(figsize=(12, 6))
+
+plt.subplot(1, 2, 1)
+plt.plot(u_values, integral_results, marker='o', linestyle='-', 
+         label='∫ (div B) dA (Monte Carlo)', linewidth=2, markersize=4)
+plt.title('Integral of div(B) over Surface vs. u', fontsize=14)
+plt.xlabel('u', fontsize=12)
+plt.ylabel('Integral Value', fontsize=12)
+plt.grid(True, alpha=0.3)
+plt.legend()
+
+plt.subplot(1, 2, 2)
+plt.semilogy(u_values, integral_errors, marker='s', linestyle='-', 
+             color='red', label='MC Error Estimate', markersize=4)
+plt.title('Integration Error vs. u', fontsize=14)
+plt.xlabel('u', fontsize=12)
+plt.ylabel('Error (log scale)', fontsize=12)
+plt.grid(True, alpha=0.3, which='both')
+plt.legend()
+
+plt.tight_layout()
+plt.show()
+
+# --- Part 6: 3D Visualization (Unchanged) ---
+
+print("\nGenerating 3D visualization...")
+
+u_vis = 2 * np.pi
+v_grid, theta_grid = np.meshgrid(np.linspace(0, 2*np.pi, 30), np.linspace(0, 2*np.pi, 30))
+X, Y, Z = f_R(u_vis, v_grid, theta_grid, a_val, b_val, gamma_val)
+
+v_quiver, theta_quiver = np.meshgrid(np.linspace(0, 2*np.pi, 6), np.linspace(0, 2*np.pi, 6))
+Xq, Yq, Zq = f_R(u_vis, v_quiver, theta_quiver, a_val, b_val, gamma_val)
+Bxq, Byq, Bzq = f_B(u_vis, v_quiver, theta_quiver, a_val, b_val, gamma_val)
+
+fig = plt.figure(figsize=(14, 10))
+ax = fig.add_subplot(111, projection='3d')
+surf = ax.plot_surface(X, Y, Z, alpha=0.7, cmap='viridis', rstride=1, cstride=1, edgecolor='none')
+ax.quiver(Xq, Yq, Zq, Bxq, Byq, Bzq, length=0.3, normalize=True, color='red', arrow_length_ratio=0.3, linewidth=1.5)
+ax.set_title(f'Manifold and Vector Field B at u = {u_vis/np.pi:.2f}π', fontsize=16)
+ax.set_xlabel('X', fontsize=12)
+ax.set_ylabel('Y', fontsize=12)
+ax.set_zlabel('Z', fontsize=12)
+ax.view_init(elev=20., azim=-65)
+fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
 plt.show()