Some vibe coding I've been doing around Lambertian reflectance and how to position lights in front of a surface such that it illuminates it perfectly with least amount of contrast (diff between darket and lightest). IT was covered in three videos so far: https://youtu.be/O_toTbA5PNc https://youtu.be/-O-LlncUULM https://youtu.be/PtKmvTEz1Og

irradiance-optimisation.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from irradiance_studies_3d_corners import calc_irradiance_points
import time
from enum import Enum

class VisualizationMode(Enum):
    OPTIMIZATION = "optimization"  # Original optimization mode
    CONTRAST_MAP = "contrast_map"  # 2D contrast map at fixed Z
    TRAJECTORY_3D = "trajectory_3d"  # 3D trajectory of best points

# Visualization mode control
MODE = VisualizationMode.TRAJECTORY_3D

FIXED_Z = 4.0  # Fixed Z height for contrast map
GRID_POINTS = 100  # Number of points in each dimension for contrast map

def calc_irradiance_batch(point_sources_batch, a_receiver=0.5, sample_points=100):
    """
    Vectorized irradiance computation for multiple 4-point Lambertian configurations.

    Assumes all sources are above the receiver (i.e., no source at Z=0).

    Parameters:
        point_sources_batch : ndarray of shape (N, 4, 4)
            Each entry is a configuration of 4 sources with (x, y, z, intensity)
        a_receiver : float
            Half-size of square receiver (centered at origin in XY plane at Z=0)
        sample_points : int
            Number of sampling points along X and Y axis

    Returns:
        Z_batch : ndarray of shape (N, sample_points, sample_points)
            Irradiance maps for each configuration
    """
    import numpy as np

    N = point_sources_batch.shape[0]
    grid = np.linspace(-a_receiver, a_receiver, sample_points)
    X, Y = np.meshgrid(grid, grid)  # Shape: (H, W)
    H, W = X.shape

    # Expand grid to match source batch
    X = X[None, :, :]  # shape (1, H, W)
    Y = Y[None, :, :]  # shape (1, H, W)

    Z_total = np.zeros((N, H, W))

    for s in range(4):
        xs = point_sources_batch[:, s, 0][:, None, None]  # (N,1,1)
        ys = point_sources_batch[:, s, 1][:, None, None]
        zs = point_sources_batch[:, s, 2][:, None, None]
        intensities = point_sources_batch[:, s, 3][:, None, None]

        dx = X - xs
        dy = Y - ys
        dz = zs

        r = np.sqrt(dx**2 + dy**2 + dz**2)  # Always > 0
        contrib = (intensities * dz / r) / (r**2)
        Z_total += contrib

    Z_total /= np.pi
    return Z_total

# Objective function: contrast ratio = (max - min) / max
def contrast_ratio(params, a_receiver=0.5, sample_points=100):
    x, y, z = params

    # Enforce symmetry: 4 point sources mirrored across both axes
    sources = [
        ( x,  y, z),
        (-x,  y, z),
        ( x, -y, z),
        (-x, -y, z)
    ]

    _, _, Z = calc_irradiance_points(sources, a_receiver, sample_points)
    max_irr = np.max(Z)
    min_irr = np.min(Z)
    return (max_irr - min_irr) / max_irr if max_irr != 0 else 1

if MODE == VisualizationMode.OPTIMIZATION:
    # Original optimization mode
    # Initial guess for (x, y, z)
    initial_guess = [2.0, 2.0, 4.0]

    # Define bounds: x, y ∈ [-3, 3], z ∈ [3, 5]
    bounds = [(-4, 4), (-4, 4), (3, 5)]

    # Perform optimization
    result = minimize(contrast_ratio, initial_guess, bounds=bounds, method='L-BFGS-B')

    # Get optimized symmetric point sources
    x_opt, y_opt, z_opt = result.x
    optimized_sources = [
        ( x_opt,  y_opt, z_opt),
        (-x_opt,  y_opt, z_opt),
        ( x_opt, -y_opt, z_opt),
        (-x_opt, -y_opt, z_opt)
    ]

    # Print results
    print("\nOptimized source coordinates:")
    for i, (x, y, z) in enumerate(optimized_sources):
        print(f"Source {i+1}: (x={x:.3f}, y={y:.3f}, z={z:.3f})")
    print(f"\nContrast ratio: {result.fun * 100:.6f}%")

    # Calculate angle between XY plane and line from origin to optimized point
    # Using first quadrant point (x_opt, y_opt, z_opt)
    horizontal_dist = np.sqrt(x_opt**2 + y_opt**2)
    angle_rad = np.arctan2(z_opt, horizontal_dist)
    angle_deg = np.degrees(angle_rad)
    print(f"\nAngle between XY plane and optimized point: {angle_deg:.2f}°")

    # Plot final irradiance distribution
    X, Y, Z = calc_irradiance_points(optimized_sources, a_receiver=0.5, sample_points=200)
    plt.figure(figsize=(6,5))
    plt.pcolormesh(X, Y, Z, shading='auto', cmap='viridis')
    plt.colorbar(label='Irradiance')
    plt.title(f"Irradiance Map (Contrast: {(result.fun * 100):.6f}%)")
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.gca().set_aspect('equal')
    plt.show()

def get_contrast_map(z):
    # Generate grid of candidate positions (only for one quadrant due to symmetry)
    x_range = np.linspace(0.01, 3.0, GRID_POINTS)
    y_range = np.linspace(0.01, 3.0, GRID_POINTS)
    Xg, Yg = np.meshgrid(x_range, y_range)

    # Flatten to create batch of light configurations
    X_flat = Xg.ravel()
    Y_flat = Yg.ravel()

    FIXED_Z = z
    sources_batch = []

    for x, y in zip(X_flat, Y_flat):
        sources = [
            ( x,  y, FIXED_Z, 1.0),
            (-x,  y, FIXED_Z, 1.0),
            ( x, -y, FIXED_Z, 1.0),
            (-x, -y, FIXED_Z, 1.0)
        ]
        sources_batch.append(sources)

    sources_batch = np.array(sources_batch)  # shape: (N, 4, 4)

    # Vectorized batch irradiance computation
    Z_maps = calc_irradiance_batch(sources_batch, a_receiver=0.5, sample_points=50)  # Try 50 for speed test

    # Compute contrast per irradiance map
    max_irr = np.max(Z_maps, axis=(1, 2))
    min_irr = np.min(Z_maps, axis=(1, 2))
    contrast = (max_irr - min_irr) / max_irr

    # Reshape back to grid
    contrast_grid = contrast.reshape(GRID_POINTS, GRID_POINTS)

    # Create full grid by mirroring across both axes
    full_contrast = np.zeros((GRID_POINTS*2, GRID_POINTS*2))
    # First quadrant (original)
    full_contrast[GRID_POINTS:, GRID_POINTS:] = contrast_grid
    # Second quadrant (mirror across Y)
    full_contrast[GRID_POINTS:, :GRID_POINTS] = np.fliplr(contrast_grid)
    # Third quadrant (mirror across both)
    full_contrast[:GRID_POINTS, :GRID_POINTS] = np.flipud(np.fliplr(contrast_grid))
    # Fourth quadrant (mirror across X)
    full_contrast[:GRID_POINTS, GRID_POINTS:] = np.flipud(contrast_grid)

    # Create full coordinate grid
    full_x = np.linspace(-3, 3, GRID_POINTS*2)
    full_y = np.linspace(-3, 3, GRID_POINTS*2)
    X_full, Y_full = np.meshgrid(full_x, full_y)

    return X_full, Y_full, full_contrast, contrast_grid, Xg, Yg

if MODE == VisualizationMode.TRAJECTORY_3D:
    # Calculate best points across Z heights
    z_heights = np.arange(1, 5, 0.2)
    best_points = []
    contrast_values = []
    
    print("\nCalculating best points across Z heights...")
    start_time = time.time()
    
    for z in z_heights:
        _, _, _, contrast_grid, Xg, Yg = get_contrast_map(z)
        best_idx = np.unravel_index(np.argmin(contrast_grid), contrast_grid.shape)
        best_x, best_y = Xg[best_idx], Yg[best_idx]
        best_points.append((best_x, best_y, z))
        contrast_values.append(contrast_grid[best_idx])
        print(f"Z={z:.1f}: Best point at ({best_x:.3f}, {best_y:.3f}) with contrast {contrast_grid[best_idx]:.6f}")
    
    duration = time.time() - start_time
    print(f"\nCalculation time: {duration:.2f} seconds")
    print(f"Number of Z heights: {len(z_heights)}")
    print(f"Average time per height: {duration/len(z_heights):.2f} seconds")
    
    # Plot 3D trajectory
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Extract coordinates
    x_coords = [p[0] for p in best_points]
    y_coords = [p[1] for p in best_points]
    z_coords = [p[2] for p in best_points]
    
    # Normalize contrast values for point sizes (scale between 20 and 200)
    contrast_array = np.array(contrast_values)
    min_contrast = np.min(contrast_array)
    max_contrast = np.max(contrast_array)
    point_sizes = 20 + 180 * (contrast_array - min_contrast) / (max_contrast - min_contrast)
    
    # Plot points and connecting lines
    ax.plot(x_coords, y_coords, z_coords, 'b-', label='Best point trajectory')
    scatter = ax.scatter(x_coords, y_coords, z_coords, c=contrast_array, 
                        cmap='viridis', marker='o', s=point_sizes)
    
    # Plot mirrored points for symmetry
    ax.plot(-np.array(x_coords), y_coords, z_coords, 'b-')
    ax.plot(x_coords, -np.array(y_coords), z_coords, 'b-')
    ax.plot(-np.array(x_coords), -np.array(y_coords), z_coords, 'b-')
    
    # Plot mirrored points with same sizes
    ax.scatter(-np.array(x_coords), y_coords, z_coords, c=contrast_array, 
              cmap='viridis', marker='o', s=point_sizes)
    ax.scatter(x_coords, -np.array(y_coords), z_coords, c=contrast_array, 
              cmap='viridis', marker='o', s=point_sizes)
    ax.scatter(-np.array(x_coords), -np.array(y_coords), z_coords, c=contrast_array, 
              cmap='viridis', marker='o', s=point_sizes)
    
    # Add colorbar
    cbar = plt.colorbar(scatter)
    cbar.set_label('Contrast Ratio')
    
    # Add labels and title
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    ax.set_title('Best Point Trajectory Across Z Heights\n(Point size ∝ Contrast)')
    
    # Set equal aspect ratio
    ax.set_box_aspect([1,1,1])
    
    plt.show()

if MODE == VisualizationMode.CONTRAST_MAP:
    # Start timing
    start_time = time.time()

    X_full, Y_full, full_contrast, contrast_grid, Xg, Yg = get_contrast_map(FIXED_Z)

    # End timing and print stats
    duration = time.time() - start_time
    print(f"\nCalculation time: {duration:.2f} seconds")
    print(f"Grid size: {GRID_POINTS}x{GRID_POINTS} = {GRID_POINTS * GRID_POINTS} points")
    print(f"Average time per point: {duration / (GRID_POINTS * GRID_POINTS) * 1000:.2f} ms")

    # Plot contrast map using mirrored results
    plt.figure(figsize=(8, 6))
    plt.pcolormesh(X_full, Y_full, full_contrast, shading='auto', cmap='viridis')
    plt.colorbar(label='Contrast Ratio')
    plt.title(f'Contrast Map at Z={FIXED_Z}')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.gca().set_aspect('equal')

    # Mark the best (lowest contrast) point and its mirrors
    best_idx = np.unravel_index(np.argmin(contrast_grid), contrast_grid.shape)
    best_x, best_y = Xg[best_idx], Yg[best_idx]
    plt.plot(best_x, best_y, 'r*', markersize=10, label=f'Best: ({best_x:.3f}, {best_y:.3f})')
    plt.plot(-best_x, best_y, 'r*', markersize=10)
    plt.plot(best_x, -best_y, 'r*', markersize=10)
    plt.plot(-best_x, -best_y, 'r*', markersize=10)
    plt.legend()
    plt.show()

irradiance-studies-3d-corners.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, CheckButtons

def calc_irradiance_points(point_sources, a_receiver=0.5, sample_points=100):
    """
    Calculate irradiance on a square XY-plane receiver at Z=0 from Lambertian point sources.

    Parameters:
        point_sources : list of tuples
            Each tuple is (x, y, z [, intensity]) specifying a source's position and optional intensity.
            If intensity is omitted, it defaults to 1.0.
        a_receiver : float
            Half-width of square receiver centered at origin (default = 0.5, i.e., 1x1 square).
        sample_points : int
            Number of points per axis in the sampling grid (higher = finer detail).

    Returns:
        X, Y : 2D arrays containing coordinates of sampling grid on receiver plane
        Z : 2D array containing irradiance values at each sample point
    """
    # Create uniform grid of (x, y) sample points across the receiver
    x = np.linspace(-a_receiver, a_receiver, sample_points)
    y = np.linspace(-a_receiver, a_receiver, sample_points)
    X, Y = np.meshgrid(x, y)  # X and Y shape: (sample_points, sample_points)

    # Initialize irradiance array (same shape as grid)
    Z = np.zeros_like(X)

    # Loop over all point sources
    for source in point_sources:
        # Unpack source position and optional intensity
        if len(source) == 3:
            x_s, y_s, z_s = source
            intensity = 1.0  # Default intensity
        else:
            x_s, y_s, z_s, intensity = source

        # Compute distances from source to each point on the receiver plane
        dx = X - x_s     # horizontal distance in X direction
        dy = Y - y_s     # horizontal distance in Y direction
        dz = z_s         # vertical height from receiver to source (constant)

        # Euclidean distance from each receiver point to the source
        r = np.sqrt(dx**2 + dy**2 + dz**2)

        # Avoid division by zero (i.e., directly under a source)
        mask = r > 0

        # Lambertian irradiance contribution: (cosθ / r²) = (dz / r) / r²
        Z[mask] += intensity * (dz / r[mask]) / (r[mask]**2)

    # Normalize by Lambertian factor (1 / π)
    Z /= np.pi

    return X, Y, Z

def calculate_irradiance(x_receiver: np.ndarray, x_c: float, z_c: float, line_length: float) -> np.ndarray:
    """
    Calculate irradiance at each point of the receiver surface from two point sources.
    
    Parameters:
        x_receiver : array
            x-coordinates of receiver surface points
        x_c, z_c : float
            Coordinates of point source center
        line_length : float
            Distance between the two point sources (will place points at ±line_length/2 from center)
    
    Returns:
        array of irradiance values at each receiver point
    """
    # Calculate positions of the two point sources
    y1 = -line_length/2
    y2 = line_length/2
    
    # Ensure x_receiver is 1D
    x_r = np.asarray(x_receiver).flatten()
    
    # Calculate distances for first point
    dx1 = x_r - x_c
    dy1 = y1
    dz1 = z_c
    r1 = np.sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1)
    
    # Calculate distances for second point
    dx2 = x_r - x_c
    dy2 = y2
    dz2 = z_c
    r2 = np.sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2)
    
    # Calculate irradiance from both points
    irradiance1 = np.zeros_like(r1)
    irradiance2 = np.zeros_like(r2)
    
    # Avoid division by zero
    np.divide(dz1, r1, out=irradiance1, where=r1 > 0)
    np.divide(irradiance1, r1*r1, out=irradiance1, where=r1 > 0)
    
    np.divide(dz2, r2, out=irradiance2, where=r2 > 0)
    np.divide(irradiance2, r2*r2, out=irradiance2, where=r2 > 0)
    
    # Sum contributions from both points and normalize
    return (irradiance1 + irradiance2) / (2 * np.pi)

def plot_top_down(x_c, z_c, x_c_mirrored, z_c_mirrored, a_receiver, ax,
                  show_both, auto_scale=True, irradiance_max=1.0, irradiance_min=0.0, line_length=0.0):
    """
    Plot top-down view of the receiver plane and irradiance from two point sources
    at the ends of a former line source.

    Parameters:
        x_c, z_c : float
            Coordinates of primary point source center
        x_c_mirrored, z_c_mirrored : float
            Coordinates of mirrored point source center
        a_receiver : float
            Half-width of receiver surface (receiver is square)
        ax : matplotlib.axes.Axes
            The axes to plot on
        show_both : bool
            Whether to show both sets of sources
        auto_scale : bool
            Whether to automatically scale irradiance axis
        irradiance_max, irradiance_min : float
            Manual irradiance scale limits
        line_length : float
            Distance between two point sources (centered around Y=0)
    """
    ax.clear()
    
    # Calculate point source positions
    y_offsets = np.array([-line_length/2, line_length/2])
    point_sources = []
    
    # Add primary point sources
    for y_l in y_offsets:
        point_sources.append((x_c, y_l, z_c))
    
    # Add mirrored point sources if show_both is True
    if show_both:
        for y_l in y_offsets:
            point_sources.append((x_c_mirrored, y_l, z_c_mirrored))
    
    # Print point source coordinates
    print("\nPoint source coordinates:")
    for i, (x, y, z) in enumerate(point_sources):
        print(f"Source {i+1}: (x={x:.3f}, y={y:.3f}, z={z:.3f})")
    
    # Calculate irradiance distribution
    X, Y, Z = calc_irradiance_points(point_sources, a_receiver)

    im = ax.pcolormesh(X, Y, Z, shading='auto', cmap='viridis')

    if auto_scale:
        min_irr = max(irradiance_min, np.min(Z))
        max_irr = np.max(Z)
        padding = (max_irr - min_irr) * 0.1
        im.set_clim(min_irr - padding, max_irr + padding)
    else:
        im.set_clim(irradiance_min, irradiance_max)

    # Update colorbar
    if not hasattr(ax, 'colorbar'):
        ax.colorbar = plt.colorbar(im, ax=ax, label='Irradiance')
    else:
        ax.colorbar.mappable = im
        ax.colorbar.update_normal(im)
        ax.colorbar.ax.figure.canvas.draw_idle()

    # Plot point sources
    ax.scatter([x_c, x_c], y_offsets, c='red', s=100, label='Point sources')
    if show_both:
        ax.scatter([x_c_mirrored, x_c_mirrored], y_offsets, c='red', s=100)

    ax.plot([-a_receiver, a_receiver, a_receiver, -a_receiver, -a_receiver],
            [-a_receiver, -a_receiver, a_receiver, a_receiver, -a_receiver],
            'b-', linewidth=2, label='Receiver plane')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    
    # Calculate and show the difference in the title
    if auto_scale:
        min_irr = max(irradiance_min, np.min(Z))
        max_irr = np.max(Z)
        diff = (max_irr - min_irr) / max_irr if max_irr > 0 else 0
    else:
        diff = (irradiance_max - irradiance_min) / irradiance_max if irradiance_max > 0 else 0
    ax.set_title(f'Top View (contrast ratio: {(diff * 100):.6f}%)')
    
    ax.grid(True)
    ax.set_aspect('equal')
    ax.set_xlim(-a_receiver*1.2, a_receiver*1.2)
    ax.set_ylim(-a_receiver*1.2, a_receiver*1.2)

    plt.subplots_adjust(right=0.85)


def plot_surfaces(angle_deg, a_receiver, distance, ax_side, ax_top, show_legend=True, show_both=True, xy_scale=1.0, irradiance_max=1.0, irradiance_min=0.0, auto_scale=True, line_length=0.0):
    """
    Plot side view and top view of the system.
    
    Parameters:
        angle_deg : float
            Angle between surfaces in degrees (positive means emitter tilted up)
        a_receiver : float
            Half-width of receiver surface
        distance : float
            Distance between centers of emitter and receiver
        ax_side : matplotlib.axes.Axes
            The axes for side view
        ax_top : matplotlib.axes.Axes
            The axes for top view
        show_legend : bool
            Whether to show the legend
        show_both : bool
            Whether to show both point sources
        xy_scale : float
            Scale factor for X and Y axes in side view
        irradiance_max : float
            Maximum irradiance value for manual scale
        irradiance_min : float
            Minimum irradiance value for manual scale
        auto_scale : bool
            Whether to automatically scale irradiance axis to show full range
        line_length : float
            Distance between the two point sources
    """
    # Clear the axes
    ax_side.clear()
    
    # Remove old twin axis if it exists
    if hasattr(ax_side, 'twin_axis'):
        try:
            ax_side.twin_axis.remove()
        except KeyError:
            pass  # Ignore if the axis was already removed
    
    # Create twin axis for irradiance
    ax2 = ax_side.twinx()
    ax_side.twin_axis = ax2  # Store reference to twin axis
    
    # Convert angles to radians
    angle_rad = np.radians(angle_deg)
    ref_angle_rad = np.radians(90+angle_deg)
    
    # Receiver surface
    x_receiver = np.linspace(-a_receiver, a_receiver, 100)
    z_receiver = np.zeros_like(x_receiver)  # z is 0 for receiver

    # Center of receiver is at (0, 0)
    # Calculate point source positions (original and mirrored)
    x_c = distance * np.cos(ref_angle_rad)
    z_c = distance * np.sin(ref_angle_rad)  # This is now z instead of y
    x_c_mirrored = -x_c  # Mirror over X-axis
    z_c_mirrored = z_c   # Z position stays the same

    # Calculate irradiance from both point sources
    irradiance1 = calculate_irradiance(x_receiver, x_c, z_c, line_length)
    irradiance2 = calculate_irradiance(x_receiver, x_c_mirrored, z_c_mirrored, line_length)  # Mirrored point source
    irradiance = irradiance1 + irradiance2 if show_both else irradiance1
    
    # Plot side view
    ax_side.plot(x_receiver, z_receiver, 'b-', linewidth=3, label='Receiver surface')
    
    # Plot original point sources (as points in side view)
    ax_side.scatter([x_c, x_c], [z_c, z_c], c='red', s=100, label=f'Point sources ({angle_deg:.1f}°)')
    
    # Plot mirrored point sources if show_both is True
    if show_both:
        ax_side.scatter([x_c_mirrored, x_c_mirrored], [z_c_mirrored, z_c_mirrored], c='red', s=100)

    # Plot the reference lines between centers
    ax_side.plot([0, x_c], [0, z_c], 'k--', label=f'Reference line ({90+angle_deg:.1f}°)')
    if show_both:
        ax_side.plot([0, x_c_mirrored], [0, z_c_mirrored], 'k--')

    # Plot center points
    if show_both:
        ax_side.scatter([0, x_c, x_c_mirrored], [0, z_c, z_c_mirrored], 
                       c=['blue', 'red', 'red'], s=100,
                       label='Center points (receiver, emitters)')
    else:
        ax_side.scatter([0, x_c], [0, z_c], 
                       c=['blue', 'red'], s=100,
                       label='Center points (receiver, emitter)')

    # Plot irradiance
    ax2.plot(x_receiver, irradiance, 'g-', label='Irradiance', alpha=0.7)
    ax2.set_ylabel('Irradiance (arbitrary units)', color='g')
    ax2.tick_params(axis='y', labelcolor='g')

    # Get axes dimensions in pixels
    bbox = ax_side.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    width_pixels = bbox.width * fig.dpi
    height_pixels = bbox.height * fig.dpi

    xy_ratio = width_pixels / height_pixels
    
    # Set fixed ranges for X and Z axes
    ax_side.set_xlim(-5 / xy_scale * xy_ratio/2, 5 / xy_scale * xy_ratio/2)
    ax_side.set_ylim(0, 5 / xy_scale)
    
    # Set irradiance axis limits based on auto_scale
    if auto_scale:
        # Use actual min and max values with padding
        min_irradiance = max(irradiance_min, np.min(irradiance))
        max_irradiance = np.max(irradiance)
        padding = (max_irradiance - min_irradiance) * 0.1  # 10% padding
        ax2.set_ylim(min_irradiance - padding, max_irradiance + padding)
    else:
        ax2.set_ylim(irradiance_min, irradiance_max)
    
    # Move irradiance axis to the right
    ax2.spines['right'].set_position(('outward', 0))

    ax_side.set_xlabel('X')
    ax_side.set_ylabel('Z')
    ax_side.grid(True)
    
    # Handle legend visibility
    if show_legend:
        # Combine legends
        lines1, labels1 = ax_side.get_legend_handles_labels()
        lines2, labels2 = ax2.get_legend_handles_labels()
        ax_side.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
    else:
        # Remove legends from both axes
        if ax_side.get_legend() is not None:
            ax_side.get_legend().remove()
        if ax2.get_legend() is not None:
            ax2.get_legend().remove()
    
    ax_side.set_title(f'Side View - Point sources at {angle_deg:.1f}°')
    
    # Plot top view
    plot_top_down(x_c, z_c, x_c_mirrored, z_c_mirrored, a_receiver, ax_top, show_both, auto_scale, irradiance_max, irradiance_min, line_length)

def create_controls():
    """Create or recreate all controls"""
    global angle_slider, a_receiver_slider, distance_slider, xy_scale_slider, irradiance_max_slider, irradiance_min_slider, line_length_slider
    global legend_checkbox, both_emitters_checkbox, auto_scale_checkbox, reset_button
    
    # Create axes
    ax_side = plt.axes([0.1, 0.5, 0.35, 0.4])
    ax_top = plt.axes([0.55, 0.5, 0.35, 0.4])
    
    # Create sliders
    ax_angle = plt.axes([0.25, 0.35, 0.65, 0.03])
    ax_receiver = plt.axes([0.25, 0.30, 0.65, 0.03])
    ax_distance = plt.axes([0.25, 0.25, 0.65, 0.03])
    ax_xy_scale = plt.axes([0.25, 0.20, 0.65, 0.03])
    ax_irradiance_scale = plt.axes([0.25, 0.15, 0.65, 0.03])
    ax_irradiance_min = plt.axes([0.25, 0.10, 0.65, 0.03])
    ax_line_length = plt.axes([0.25, 0.05, 0.65, 0.03])
    
    angle_slider = Slider(ax_angle, 'Emitter Angle (deg)', -90, 90, valinit=angle_init)
    a_receiver_slider = Slider(ax_receiver, 'Receiver half-width', 0.1, 2.0, valinit=a_receiver_init)
    distance_slider = Slider(ax_distance, 'Distance between centers', 0.5, 8.0, valinit=distance_init)
    xy_scale_slider = Slider(ax_xy_scale, 'X/Y Scale', 0.1, 10.0, valinit=xy_scale_init)
    irradiance_max_slider = Slider(ax_irradiance_scale, 'Max Irradiance', 0.01, 1.0, valinit=irradiance_max_init)
    irradiance_min_slider = Slider(ax_irradiance_min, 'Min Irradiance', 0.0, 0.15, valinit=irradiance_min_init)
    line_length_slider = Slider(ax_line_length, 'Line Length', 0.0, 6.0, valinit=line_length_init)
    
    # Create checkboxes
    ax_legend = plt.axes([0.8, 0.025, 0.1, 0.04])
    legend_checkbox = CheckButtons(ax_legend, ['Show Legend'], [show_legend_init])
    
    ax_both = plt.axes([0.65, 0.025, 0.1, 0.04])
    both_emitters_checkbox = CheckButtons(ax_both, ['Both Emitters'], [show_both_init])
    
    ax_auto_scale = plt.axes([0.5, 0.025, 0.1, 0.04])
    auto_scale_checkbox = CheckButtons(ax_auto_scale, ['Auto Scale'], [auto_scale_init])
    
    # Create reset button
    reset_ax = plt.axes([0.35, 0.025, 0.1, 0.04])
    reset_button = Button(reset_ax, 'Reset', color='lightgoldenrodyellow', hovercolor='0.975')
    
    return ax_side, ax_top

def update(val):
    """Update the plot when sliders change"""
    # Get current values from all controls
    angle = angle_slider.val
    a_rec = a_receiver_slider.val
    dist = distance_slider.val
    xy_scale_val = xy_scale_slider.val
    irr_max_val = irradiance_max_slider.val
    irr_min_val = irradiance_min_slider.val
    line_length = line_length_slider.val
    show_legend = legend_checkbox.get_status()[0]
    show_both = both_emitters_checkbox.get_status()[0]
    auto_scale = auto_scale_checkbox.get_status()[0]
    
    # Disable/enable irradiance scale slider based on auto_scale
    irradiance_max_slider.active = not auto_scale
    irradiance_min_slider.active = not auto_scale
    
    # Clear both plot areas
    ax_side.clear()
    ax_top.clear()
    
    # Plot the surfaces
    plot_surfaces(angle, a_rec, dist, ax_side, ax_top, 
                 show_legend=show_legend, 
                 show_both=show_both,
                 xy_scale=xy_scale_val,
                 irradiance_max=irr_max_val,
                 irradiance_min=irr_min_val,
                 auto_scale=auto_scale,
                 line_length=line_length)
    
    # Force redraw
    fig.canvas.draw_idle()

def reset(event):
    """Reset all controls to initial values"""
    global angle_slider, a_receiver_slider, distance_slider, xy_scale_slider, irradiance_max_slider, irradiance_min_slider, line_length_slider
    global legend_checkbox, both_emitters_checkbox, auto_scale_checkbox
    
    angle_slider.reset()
    a_receiver_slider.reset()
    distance_slider.reset()
    xy_scale_slider.reset()
    irradiance_max_slider.reset()
    irradiance_min_slider.reset()
    line_length_slider.reset()
    legend_checkbox.set_active(0) if not show_legend_init else None
    both_emitters_checkbox.set_active(0) if not show_both_init else None
    auto_scale_checkbox.set_active(0) if not auto_scale_init else None
    update(None)  # Force an update after reset

if __name__ == '__main__':
    # Create the figure
    fig = plt.figure(figsize=(15, 8))

    # Initial values
    angle_init = -30
    a_receiver_init = 1.0
    distance_init = 4.0
    xy_scale_init = 3.5
    irradiance_max_init = 0.039
    irradiance_min_init = 0.031
    line_length_init = 4
    show_legend_init = False
    show_both_init = True
    auto_scale_init = True

    # Global variables for sliders and controls
    angle_slider = None
    a_receiver_slider = None
    distance_slider = None
    xy_scale_slider = None
    irradiance_max_slider = None
    irradiance_min_slider = None
    line_length_slider = None
    legend_checkbox = None
    both_emitters_checkbox = None
    auto_scale_checkbox = None
    reset_button = None

    # Create initial controls
    ax_side, ax_top = create_controls()
    
    # Register the update function with each slider and checkbox
    angle_slider.on_changed(update)
    a_receiver_slider.on_changed(update)
    distance_slider.on_changed(update)
    xy_scale_slider.on_changed(update)
    irradiance_max_slider.on_changed(update)
    irradiance_min_slider.on_changed(update)
    line_length_slider.on_changed(update)
    legend_checkbox.on_clicked(update)
    both_emitters_checkbox.on_clicked(update)
    auto_scale_checkbox.on_clicked(update)
    reset_button.on_clicked(reset)
    
    # Initial plot
    plot_surfaces(angle_init, a_receiver_init, distance_init, ax_side, ax_top, 
                 show_legend=show_legend_init, 
                 show_both=show_both_init,
                 xy_scale=xy_scale_init,
                 irradiance_max=irradiance_max_init,
                 irradiance_min=irradiance_min_init,
                 auto_scale=auto_scale_init,
                 line_length=line_length_init)
    
    # Set initial state of irradiance scale slider
    irradiance_max_slider.active = not auto_scale_init
    irradiance_min_slider.active = not auto_scale_init
    
    plt.show()

irradiance-studies-3d-lines.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, CheckButtons

def calculate_irradiance(x_receiver: np.ndarray, x_c: float, z_c: float, line_length: float) -> np.ndarray:
    """
    Calculate irradiance at each point of the receiver surface from a line source.

    Parameters:
        x_receiver : array
            x-coordinates of receiver surface points
        x_c, z_c : float
            Coordinates of line source center
        line_length : float
            Length of the line source (total length, will extend ±line_length/2 from center)

    Returns:
        array of irradiance values at each receiver point
    """
    n_points = 25
    y_points = np.linspace(-line_length/2, line_length/2, n_points)

    x_r = np.asarray(x_receiver).flatten()[:, np.newaxis]  # Shape: (n_receiver, 1)
    y_l = y_points[np.newaxis, :]                          # Shape: (1, n_points)

    dx = x_r - x_c                                          # (n_receiver, n_points)
    dy = -y_l                                               # (1, n_points), broadcasted
    dz = z_c                                                # Scalar

    r = np.sqrt(dx*dx + dy*dy + dz*dz)                      # Total distance from source point to receiver

    irradiance = np.zeros_like(r)
    np.divide(dz, r, out=irradiance, where=r > 0)
    np.divide(irradiance, r*r, out=irradiance, where=r > 0)

    return np.sum(irradiance, axis=1) / (np.pi * n_points)


def plot_top_down(x_c, z_c, x_c_mirrored, z_c_mirrored, a_receiver, ax, show_both, auto_scale=True, irradiance_max=1.0, irradiance_min=0.0, line_length=0.0):
    """
    Plot top-down view of the receiver plane and line sources.
    Line sources are kept at Y=0, only their X positions vary.
    
    Parameters:
        x_c, z_c : float
            Coordinates of line source center
        x_c_mirrored, z_c_mirrored : float
            Coordinates of mirrored line source center
        a_receiver : float
            Half-width of receiver surface
        ax : matplotlib.axes.Axes
            The axes to plot on
        show_both : bool
            Whether to show both line sources
        auto_scale : bool
            Whether to automatically scale irradiance axis to show full range
        irradiance_max : float
            Maximum irradiance value for manual scale
        irradiance_min : float
            Minimum irradiance value for manual scale
        line_length : float
            Length of the line source
    """
    # Clear the axes
    ax.clear()
    
    sample_points = 100

    # Create a grid of points on the receiver plane
    x = np.linspace(-a_receiver, a_receiver, sample_points)
    y = np.linspace(-a_receiver, a_receiver, sample_points)
    X, Y = np.meshgrid(x, y)
    
    # Create points along the line source
    y_points = np.linspace(-line_length/2, line_length/2, sample_points)
    
    # Reshape arrays for broadcasting
    X_flat = X.flatten()[:, np.newaxis]  # Shape: (n_points^2, 1)
    Y_flat = Y.flatten()[:, np.newaxis]  # Shape: (n_points^2, 1)
    y_l = y_points[np.newaxis, :]        # Shape: (1, n_points)
    
    # Calculate distance vectors for all combinations
    dx1 = X_flat - x_c                    # Shape: (n_points^2, n_points)
    dy1 = Y_flat - y_l                    # Shape: (n_points^2, n_points)
    dz1 = z_c                             # Scalar
    
    r1 = np.sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1)  # Shape: (n_points^2, n_points)
    
    # Calculate irradiance contribution from first line source
    mask1 = r1 > 0
    Z1 = np.zeros_like(r1)
    Z1[mask1] = (dz1 / r1[mask1]) / (r1[mask1] * r1[mask1])
    
    # Sum contributions from all points along the line
    Z = np.sum(Z1, axis=1)
    
    if show_both:
        # Calculate for second line source
        dx2 = X_flat - x_c_mirrored
        dy2 = Y_flat - y_l
        dz2 = z_c_mirrored
        
        r2 = np.sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2)
        
        mask2 = r2 > 0
        Z2 = np.zeros_like(r2)
        Z2[mask2] = (dz2 / r2[mask2]) / (r2[mask2] * r2[mask2])
        
        Z += np.sum(Z2, axis=1)
    
    # Normalize by 1/π factor from Lambertian source and number of points
    Z /= (np.pi * sample_points)
    
    # Reshape back to 2D grid
    Z = Z.reshape(X.shape)
    
    # Plot the receiver plane with irradiance
    im = ax.pcolormesh(X, Y, Z, shading='auto', cmap='viridis')
    
    # Set colorbar limits based on auto_scale
    if auto_scale:
        # Use actual min and max values with padding
        min_irradiance = max(irradiance_min, np.min(Z))
        max_irradiance = np.max(Z)
        padding = (max_irradiance - min_irradiance) * 0.1  # 10% padding
        im.set_clim(min_irradiance - padding, max_irradiance + padding)
    else:
        im.set_clim(irradiance_min, irradiance_max)
    
    # Add or update colorbar
    if not hasattr(ax, 'colorbar'):
        ax.colorbar = plt.colorbar(im, ax=ax, label='Irradiance')
    else:
        ax.colorbar.mappable = im
        ax.colorbar.update_normal(im)
        ax.colorbar.ax.figure.canvas.draw_idle()
    
    # Plot line sources
    ax.plot([x_c, x_c], [-line_length/2, line_length/2], 'r-', linewidth=2, label='Line source')
    if show_both:
        ax.plot([x_c_mirrored, x_c_mirrored], [-line_length/2, line_length/2], 'r-', linewidth=2)
    
    # Plot receiver outline
    ax.plot([-a_receiver, a_receiver, a_receiver, -a_receiver, -a_receiver],
            [-a_receiver, -a_receiver, a_receiver, a_receiver, -a_receiver],
            'b-', linewidth=2, label='Receiver plane')
    
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    
    # Calculate and show the difference in the title
    if auto_scale:
        min_irr = max(irradiance_min, np.min(Z))
        max_irr = np.max(Z)
        diff = max_irr - min_irr
    else:
        diff = irradiance_max - irradiance_min
    ax.set_title(f'Top View (diff: {diff:.6f})')

    ax.grid(True)
    
    # Set fixed aspect ratio and limits
    ax.set_aspect('equal')
    ax.set_xlim(-a_receiver*1.2, a_receiver*1.2)
    ax.set_ylim(-a_receiver*1.2, a_receiver*1.2)
    
    # Adjust layout to prevent shrinking
    plt.subplots_adjust(right=0.85)  # Make room for colorbar

def plot_surfaces(angle_deg, a_receiver, distance, ax_side, ax_top, show_legend=True, show_both=True, xy_scale=1.0, irradiance_max=1.0, irradiance_min=0.0, auto_scale=True, line_length=0.0):
    """
    Plot side view and top view of the system.
    
    Parameters:
        angle_deg : float
            Angle between surfaces in degrees (positive means emitter tilted up)
        a_receiver : float
            Half-width of receiver surface
        distance : float
            Distance between centers of emitter and receiver
        ax_side : matplotlib.axes.Axes
            The axes for side view
        ax_top : matplotlib.axes.Axes
            The axes for top view
        show_legend : bool
            Whether to show the legend
        show_both : bool
            Whether to show both line sources
        xy_scale : float
            Scale factor for X and Y axes in side view
        irradiance_max : float
            Maximum irradiance value for manual scale
        irradiance_min : float
            Minimum irradiance value for manual scale
        auto_scale : bool
            Whether to automatically scale irradiance axis to show full range
        line_length : float
            Length of the line source
    """
    # Clear the axes
    ax_side.clear()
    
    # Remove old twin axis if it exists
    if hasattr(ax_side, 'twin_axis'):
        try:
            ax_side.twin_axis.remove()
        except KeyError:
            pass  # Ignore if the axis was already removed
    
    # Create twin axis for irradiance
    ax2 = ax_side.twinx()
    ax_side.twin_axis = ax2  # Store reference to twin axis
    
    # Convert angles to radians
    angle_rad = np.radians(angle_deg)
    ref_angle_rad = np.radians(90+angle_deg)
    
    # Receiver surface
    x_receiver = np.linspace(-a_receiver, a_receiver, 100)
    z_receiver = np.zeros_like(x_receiver)  # z is 0 for receiver

    # Center of receiver is at (0, 0)
    # Calculate line source positions (original and mirrored)
    x_c = distance * np.cos(ref_angle_rad)
    z_c = distance * np.sin(ref_angle_rad)  # This is now z instead of y
    x_c_mirrored = -x_c  # Mirror over X-axis
    z_c_mirrored = z_c   # Z position stays the same

    # Calculate irradiance from both line sources
    irradiance1 = calculate_irradiance(x_receiver, x_c, z_c, line_length)
    irradiance2 = calculate_irradiance(x_receiver, x_c_mirrored, z_c_mirrored, line_length)  # Mirrored line source
    irradiance = irradiance1 + irradiance2 if show_both else irradiance1
    
    # Plot side view
    ax_side.plot(x_receiver, z_receiver, 'b-', linewidth=3, label='Receiver surface')
    
    # Plot original line source (as a point in side view since it extends in Y direction)
    ax_side.scatter([x_c], [z_c], c='red', s=100, label=f'Line source ({angle_deg:.1f}°)')
    
    # Plot mirrored line source if show_both is True
    if show_both:
        ax_side.scatter([x_c_mirrored], [z_c_mirrored], c='red', s=100)

    # Plot the reference lines between centers
    ax_side.plot([0, x_c], [0, z_c], 'k--', label=f'Reference line ({90+angle_deg:.1f}°)')
    if show_both:
        ax_side.plot([0, x_c_mirrored], [0, z_c_mirrored], 'k--')

    # Plot center points
    if show_both:
        ax_side.scatter([0, x_c, x_c_mirrored], [0, z_c, z_c_mirrored], 
                       c=['blue', 'red', 'red'], s=100,
                       label='Center points (receiver, emitters)')
    else:
        ax_side.scatter([0, x_c], [0, z_c], 
                       c=['blue', 'red'], s=100,
                       label='Center points (receiver, emitter)')

    # Plot irradiance
    ax2.plot(x_receiver, irradiance, 'g-', label='Irradiance', alpha=0.7)
    ax2.set_ylabel('Irradiance (arbitrary units)', color='g')
    ax2.tick_params(axis='y', labelcolor='g')

    # Get axes dimensions in pixels
    bbox = ax_side.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    width_pixels = bbox.width * fig.dpi
    height_pixels = bbox.height * fig.dpi

    xy_ratio = width_pixels / height_pixels
    
    # Set fixed ranges for X and Z axes
    ax_side.set_xlim(-5 / xy_scale * xy_ratio/2, 5 / xy_scale * xy_ratio/2)
    ax_side.set_ylim(0, 5 / xy_scale)
    
    # Set irradiance axis limits based on auto_scale
    if auto_scale:
        # Use actual min and max values with padding
        min_irradiance = max(irradiance_min, np.min(irradiance))
        max_irradiance = np.max(irradiance)
        ax2.set_ylim(min_irradiance, max_irradiance)
    else:
        ax2.set_ylim(irradiance_min, irradiance_max)
    
    # Move irradiance axis to the right
    ax2.spines['right'].set_position(('outward', 0))

    ax_side.set_xlabel('X')
    ax_side.set_ylabel('Z')
    ax_side.grid(True)
    
    # Handle legend visibility
    if show_legend:
        # Combine legends
        lines1, labels1 = ax_side.get_legend_handles_labels()
        lines2, labels2 = ax2.get_legend_handles_labels()
        ax_side.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
    else:
        # Remove legends from both axes
        if ax_side.get_legend() is not None:
            ax_side.get_legend().remove()
        if ax2.get_legend() is not None:
            ax2.get_legend().remove()
    
    ax_side.set_title(f'Side View - Line source at {angle_deg:.1f}°')
    
    # Plot top view
    plot_top_down(x_c, z_c, x_c_mirrored, z_c_mirrored, a_receiver, ax_top, show_both, auto_scale, irradiance_max, irradiance_min, line_length)

# Create the figure
fig = plt.figure(figsize=(15, 8))

# Initial values
angle_init = -27.5
a_receiver_init = 1.0
distance_init = 4.0
xy_scale_init = 3.5
irradiance_max_init = 0.039
irradiance_min_init = 0.031
line_length_init = 0.5
show_legend_init = False
show_both_init = True
auto_scale_init = False

# Global variables for sliders and controls
angle_slider = None
a_receiver_slider = None
distance_slider = None
xy_scale_slider = None
irradiance_max_slider = None
irradiance_min_slider = None
line_length_slider = None
legend_checkbox = None
both_emitters_checkbox = None
auto_scale_checkbox = None
reset_button = None

def create_controls():
    """Create or recreate all controls"""
    global angle_slider, a_receiver_slider, distance_slider, xy_scale_slider, irradiance_max_slider, irradiance_min_slider, line_length_slider
    global legend_checkbox, both_emitters_checkbox, auto_scale_checkbox, reset_button
    
    # Create axes
    ax_side = plt.axes([0.1, 0.5, 0.35, 0.4])
    ax_top = plt.axes([0.55, 0.5, 0.35, 0.4])
    
    # Create sliders
    ax_angle = plt.axes([0.25, 0.35, 0.65, 0.03])
    ax_receiver = plt.axes([0.25, 0.30, 0.65, 0.03])
    ax_distance = plt.axes([0.25, 0.25, 0.65, 0.03])
    ax_xy_scale = plt.axes([0.25, 0.20, 0.65, 0.03])
    ax_irradiance_scale = plt.axes([0.25, 0.15, 0.65, 0.03])
    ax_irradiance_min = plt.axes([0.25, 0.10, 0.65, 0.03])
    ax_line_length = plt.axes([0.25, 0.05, 0.65, 0.03])
    
    angle_slider = Slider(ax_angle, 'Emitter Angle (deg)', -90, 90, valinit=angle_init)
    a_receiver_slider = Slider(ax_receiver, 'Receiver half-width', 0.1, 2.0, valinit=a_receiver_init)
    distance_slider = Slider(ax_distance, 'Distance between centers', 0.5, 8.0, valinit=distance_init)
    xy_scale_slider = Slider(ax_xy_scale, 'X/Y Scale', 0.1, 10.0, valinit=xy_scale_init)
    irradiance_max_slider = Slider(ax_irradiance_scale, 'Max Irradiance', 0.01, 1.0, valinit=irradiance_max_init)
    irradiance_min_slider = Slider(ax_irradiance_min, 'Min Irradiance', 0.0, 0.15, valinit=irradiance_min_init)
    line_length_slider = Slider(ax_line_length, 'Line Length', 0.0, 2.0, valinit=line_length_init)
    
    # Create checkboxes
    ax_legend = plt.axes([0.8, 0.025, 0.1, 0.04])
    legend_checkbox = CheckButtons(ax_legend, ['Show Legend'], [show_legend_init])
    
    ax_both = plt.axes([0.65, 0.025, 0.1, 0.04])
    both_emitters_checkbox = CheckButtons(ax_both, ['Both Emitters'], [show_both_init])
    
    ax_auto_scale = plt.axes([0.5, 0.025, 0.1, 0.04])
    auto_scale_checkbox = CheckButtons(ax_auto_scale, ['Auto Scale'], [auto_scale_init])
    
    # Create reset button
    reset_ax = plt.axes([0.35, 0.025, 0.1, 0.04])
    reset_button = Button(reset_ax, 'Reset', color='lightgoldenrodyellow', hovercolor='0.975')
    
    return ax_side, ax_top

def update(val):
    """Update the plot when sliders change"""
    # Get current values from all controls
    angle = angle_slider.val
    a_rec = a_receiver_slider.val
    dist = distance_slider.val
    xy_scale_val = xy_scale_slider.val
    irr_max_val = irradiance_max_slider.val
    irr_min_val = irradiance_min_slider.val
    line_length = line_length_slider.val
    show_legend = legend_checkbox.get_status()[0]
    show_both = both_emitters_checkbox.get_status()[0]
    auto_scale = auto_scale_checkbox.get_status()[0]
    
    # Disable/enable irradiance scale slider based on auto_scale
    irradiance_max_slider.active = not auto_scale
    irradiance_min_slider.active = not auto_scale
    
    # Clear both plot areas
    ax_side.clear()
    ax_top.clear()
    
    # Plot the surfaces
    plot_surfaces(angle, a_rec, dist, ax_side, ax_top, 
                 show_legend=show_legend, 
                 show_both=show_both,
                 xy_scale=xy_scale_val,
                 irradiance_max=irr_max_val,
                 irradiance_min=irr_min_val,
                 auto_scale=auto_scale,
                 line_length=line_length)
    
    # Force redraw
    fig.canvas.draw_idle()

def reset(event):
    """Reset all controls to initial values"""
    global angle_slider, a_receiver_slider, distance_slider, xy_scale_slider, irradiance_max_slider, irradiance_min_slider, line_length_slider
    global legend_checkbox, both_emitters_checkbox, auto_scale_checkbox
    
    angle_slider.reset()
    a_receiver_slider.reset()
    distance_slider.reset()
    xy_scale_slider.reset()
    irradiance_max_slider.reset()
    irradiance_min_slider.reset()
    line_length_slider.reset()
    legend_checkbox.set_active(0) if not show_legend_init else None
    both_emitters_checkbox.set_active(0) if not show_both_init else None
    auto_scale_checkbox.set_active(0) if not auto_scale_init else None
    update(None)  # Force an update after reset

# Create initial controls
ax_side, ax_top = create_controls()

# Register the update function with each slider and checkbox
angle_slider.on_changed(update)
a_receiver_slider.on_changed(update)
distance_slider.on_changed(update)
xy_scale_slider.on_changed(update)
irradiance_max_slider.on_changed(update)
irradiance_min_slider.on_changed(update)
line_length_slider.on_changed(update)
legend_checkbox.on_clicked(update)
both_emitters_checkbox.on_clicked(update)
auto_scale_checkbox.on_clicked(update)
reset_button.on_clicked(reset)

# Initial plot
plot_surfaces(angle_init, a_receiver_init, distance_init, ax_side, ax_top, 
             show_legend=show_legend_init, 
             show_both=show_both_init,
             xy_scale=xy_scale_init,
             irradiance_max=irradiance_max_init,
             irradiance_min=irradiance_min_init,
             auto_scale=auto_scale_init,
             line_length=line_length_init)

# Set initial state of irradiance scale slider
irradiance_max_slider.active = not auto_scale_init
irradiance_min_slider.active = not auto_scale_init

plt.show()

irradiance-studies.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button, CheckButtons

def calculate_irradiance(x_receiver: np.ndarray, x_c: float, y_c: float, angle_rad: float, a_emitter: float) -> np.ndarray:
    """
    Calculate irradiance at each point of the receiver surface using numerical integration.
    This version properly handles tilted emitters by computing the full 2D integral.
    
    Parameters:
        x_receiver : array
            x-coordinates of receiver surface points
        x_c, y_c : float
            Coordinates of emitter center
        angle_rad : float
            Angle of emitter in radians (positive means tilted up)
        a_emitter : float
            Half-width of emitter surface
    
    Returns:
        array of irradiance values at each receiver point
    """
    # Vector along emitter (both x and y components needed for tilted emitter)
    v_x = a_emitter * np.cos(angle_rad)
    v_y = a_emitter * np.sin(angle_rad)
    
    # Emitter endpoints
    x_emitter_start = x_c - v_x
    y_emitter_start = y_c - v_y
    x_emitter_end = x_c + v_x
    y_emitter_end = y_c + v_y
    
    # Calculate irradiance at each receiver point
    irradiance = np.zeros_like(x_receiver)
    
    # Number of points to sample along emitter for numerical integration
    n_samples = 100
    emitter_x = np.linspace(x_emitter_start, x_emitter_end, n_samples)
    emitter_y = np.linspace(y_emitter_start, y_emitter_end, n_samples)
    
    # For each receiver point
    for i, x_r in enumerate(x_receiver):
        # For each emitter point
        for j in range(n_samples):
            # Distance vector from emitter point to receiver point
            dx = x_r - emitter_x[j]
            dy = -emitter_y[j]  # y_r is always 0
            r = np.sqrt(dx*dx + dy*dy)
            
            if r > 0:  # Avoid division by zero
                # Calculate angles for Lambertian cosine terms
                # θi: angle between emitter normal and direction to receiver
                # θo: angle between receiver normal and direction to emitter
                
                # Emitter normal is perpendicular to emitter direction
                emitter_normal_x = -np.sin(angle_rad)
                emitter_normal_y = np.cos(angle_rad)
                
                # Direction to receiver point (normalized)
                dir_x = dx / r
                dir_y = dy / r
                
                # Calculate cosines of angles
                cos_theta_i = abs(emitter_normal_x * dir_x + emitter_normal_y * dir_y)
                cos_theta_o = abs(dir_y)  # Receiver normal is (0,1)
                
                # Add contribution from this emitter point
                # Using inverse square law and both cosine terms
                irradiance[i] += (cos_theta_i * cos_theta_o) / (r * r)
    
    # Normalize by number of samples and emitter length
    # The 1/π factor comes from the Lambertian source
    irradiance *= (a_emitter / n_samples) / np.pi
    
    return irradiance

def plot_surfaces(angle_deg, a_emitter, a_receiver, distance, ax, show_legend=True, show_both=True, xy_scale=1.0, irradiance_scale=1.0):
    """
    Plot emitter and receiver surfaces with specified angle and sizes.
    Includes a second emitter mirrored over the Y-axis.
    
    Parameters:
        angle_deg : float
            Angle between surfaces in degrees (positive means emitter tilted up)
        a_emitter : float
            Half-width of emitter surface
        a_receiver : float
            Half-width of receiver surface
        distance : float
            Distance between centers of emitter and receiver
        ax : matplotlib.axes.Axes
            The axes to plot on
        show_legend : bool
            Whether to show the legend
        show_both : bool
            Whether to show both emitters or just one
        xy_scale : float
            Scale factor for X and Y axes (maintains aspect ratio)
        irradiance_scale : float
            Scale factor for irradiance axis
    """
    # Clear the axes
    ax.clear()
    
    # Remove old twin axis if it exists
    if hasattr(ax, 'twin_axis'):
        ax.twin_axis.remove()
    
    # Create twin axis for irradiance
    ax2 = ax.twinx()
    ax.twin_axis = ax2  # Store reference to twin axis
    
    # Convert angles to radians
    angle_rad = np.radians(angle_deg)
    ref_angle_rad = np.radians(90+angle_deg)
    
    # Receiver surface
    x_receiver = np.linspace(-a_receiver, a_receiver, 100)
    y_receiver = np.zeros_like(x_receiver)

    # Vector along emitter
    v_x = a_emitter * np.cos(angle_rad)
    v_y = a_emitter * np.sin(angle_rad)

    # Center of receiver is at (0, 0)
    # Calculate emitter center positions (original and mirrored)
    x_c = distance * np.cos(ref_angle_rad)
    y_c = distance * np.sin(ref_angle_rad)
    x_c_mirrored = -x_c  # Mirror over Y-axis
    y_c_mirrored = y_c   # Y position stays the same

    # Emitter endpoints (original)
    x_emitter_start = x_c - v_x
    y_emitter_start = y_c - v_y
    x_emitter_end = x_c + v_x
    y_emitter_end = y_c + v_y

    # Emitter endpoints (mirrored)
    x_emitter_start_mirrored = x_c_mirrored + v_x  # Note the + instead of - due to mirroring
    y_emitter_start_mirrored = y_c_mirrored - v_y
    x_emitter_end_mirrored = x_c_mirrored - v_x    # Note the - instead of + due to mirroring
    y_emitter_end_mirrored = y_c_mirrored + v_y

    # Calculate irradiance from both emitters
    irradiance1 = calculate_irradiance(x_receiver, x_c, y_c, angle_rad, a_emitter)
    irradiance2 = calculate_irradiance(x_receiver, x_c_mirrored, y_c_mirrored, -angle_rad, a_emitter)  # Negative angle for mirroring
    irradiance = irradiance1 + irradiance2 if show_both else irradiance1
    
    # Plot surfaces
    ax.plot(x_receiver, y_receiver, 'b-', linewidth=3, label='Receiver surface')
    
    # Plot original emitter
    ax.plot([x_emitter_start, x_emitter_end], [y_emitter_start, y_emitter_end],
            'r-', linewidth=3, label=f'Emitter surface ({angle_deg:.1f}°)')
    
    # Plot mirrored emitter if show_both is True
    if show_both:
        ax.plot([x_emitter_start_mirrored, x_emitter_end_mirrored], 
                [y_emitter_start_mirrored, y_emitter_end_mirrored],
                'r-', linewidth=3)

    # Plot the reference lines between centers
    ax.plot([0, x_c], [0, y_c], 'k--', label=f'Reference line ({90+angle_deg:.1f}°)')
    if show_both:
        ax.plot([0, x_c_mirrored], [0, y_c_mirrored], 'k--')

    # Plot center points
    if show_both:
        ax.scatter([0, x_c, x_c_mirrored], [0, y_c, y_c_mirrored], 
                   c=['blue', 'red', 'red'], s=100,
                   label='Center points (receiver, emitters)')
    else:
        ax.scatter([0, x_c], [0, y_c], 
                   c=['blue', 'red'], s=100,
                   label='Center points (receiver, emitter)')

    # Plot irradiance
    ax2.plot(x_receiver, irradiance, 'g-', label='Irradiance', alpha=0.7)
    ax2.set_ylabel('Irradiance (arbitrary units)', color='g')
    ax2.tick_params(axis='y', labelcolor='g')

    # Calculate appropriate limits to show both surfaces
    x_min = min(-a_receiver, x_emitter_start, x_emitter_start_mirrored if show_both else float('inf')) - 0.5
    x_max = max(a_receiver, x_emitter_end, x_emitter_end_mirrored if show_both else float('-inf')) + 0.5
    y_max = max(0, y_emitter_start, y_emitter_end, 
                y_emitter_start_mirrored if show_both else 0, 
                y_emitter_end_mirrored if show_both else 0) + 0.5

    # Get axes dimensions in pixels
    bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    width_pixels = bbox.width * fig.dpi
    height_pixels = bbox.height * fig.dpi

    xy_ratio = width_pixels / height_pixels

    # Apply manual scale factor while maintaining aspect ratio
    # scale_factor = min(base_x_scale, base_y_scale) / xy_scale
    
    # Set fixed ranges for X and Y axes
    ax.set_xlim(-5 / xy_scale * xy_ratio/2, 5 / xy_scale * xy_ratio/2)
    ax.set_ylim(0, 5 / xy_scale)
    
    # Set fixed range for irradiance axis
    ax2.set_ylim(0, 1 / irradiance_scale)
    
    # Move irradiance axis to the right
    ax2.spines['right'].set_position(('outward', 0))

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.grid(True)
    
    # Handle legend visibility
    if show_legend:
        # Combine legends
        lines1, labels1 = ax.get_legend_handles_labels()
        lines2, labels2 = ax2.get_legend_handles_labels()
        ax.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
    else:
        # Remove legends from both axes
        if ax.get_legend() is not None:
            ax.get_legend().remove()
        if ax2.get_legend() is not None:
            ax2.get_legend().remove()
    
    ax.set_title(f'Emitter at {angle_deg:.1f}°')

def update(val):
    """Update the plot when sliders change"""
    angle = angle_slider.val
    a_emit = a_emitter_slider.val
    a_rec = a_receiver_slider.val
    dist = distance_slider.val
    xy_scale_val = xy_scale_slider.val
    irr_scale_val = irradiance_scale_slider.val
    show_legend = legend_checkbox.get_status()[0]  # Get first checkbox state
    show_both = both_emitters_checkbox.get_status()[0]  # Get second checkbox state
    plot_surfaces(angle, a_emit, a_rec, dist, ax, 
                 show_legend=show_legend, 
                 show_both=show_both,
                 xy_scale=xy_scale_val,
                 irradiance_scale=irr_scale_val)
    fig.canvas.draw_idle()

# Create the figure and axes
fig, ax = plt.subplots(figsize=(10, 8))
plt.subplots_adjust(bottom=0.5)  # Make room for additional sliders

# Initial values
angle_init = -45
a_emitter_init = 1.0
a_receiver_init = 1.0
distance_init = 2.0
xy_scale_init = 1.0
irradiance_scale_init = 1.0

# Create sliders
ax_angle = plt.axes([0.25, 0.35, 0.65, 0.03])
ax_emitter = plt.axes([0.25, 0.30, 0.65, 0.03])
ax_receiver = plt.axes([0.25, 0.25, 0.65, 0.03])
ax_distance = plt.axes([0.25, 0.20, 0.65, 0.03])
ax_xy_scale = plt.axes([0.25, 0.15, 0.65, 0.03])
ax_irradiance_scale = plt.axes([0.25, 0.10, 0.65, 0.03])

angle_slider = Slider(ax_angle, 'Emitter Angle (deg)', -90, 90, valinit=angle_init)
a_emitter_slider = Slider(ax_emitter, 'Emitter half-width', 0.1, 2.0, valinit=a_emitter_init)
a_receiver_slider = Slider(ax_receiver, 'Receiver half-width', 0.1, 2.0, valinit=a_receiver_init)
distance_slider = Slider(ax_distance, 'Distance between centers', 0.5, 8.0, valinit=distance_init)
xy_scale_slider = Slider(ax_xy_scale, 'X/Y Scale', 0.1, 10.0, valinit=xy_scale_init)
irradiance_scale_slider = Slider(ax_irradiance_scale, 'Irradiance Scale', 0.1, 60.0, valinit=irradiance_scale_init)

# Create checkboxes for legend and both emitters
ax_legend = plt.axes([0.8, 0.025, 0.1, 0.04])
legend_checkbox = CheckButtons(ax_legend, ['Show Legend'], [True])

ax_both = plt.axes([0.65, 0.025, 0.1, 0.04])
both_emitters_checkbox = CheckButtons(ax_both, ['Both Emitters'], [True])

# Register the update function with each slider
angle_slider.on_changed(update)
a_emitter_slider.on_changed(update)
a_receiver_slider.on_changed(update)
distance_slider.on_changed(update)
xy_scale_slider.on_changed(update)
irradiance_scale_slider.on_changed(update)
legend_checkbox.on_clicked(update)
both_emitters_checkbox.on_clicked(update)

# Add a reset button
reset_ax = plt.axes([0.5, 0.025, 0.1, 0.04])
reset_button = Button(reset_ax, 'Reset', color='lightgoldenrodyellow', hovercolor='0.975')

def reset(event):
    angle_slider.reset()
    a_emitter_slider.reset()
    a_receiver_slider.reset()
    distance_slider.reset()
    xy_scale_slider.reset()
    irradiance_scale_slider.reset()

reset_button.on_clicked(reset)

# Initial plot
plot_surfaces(angle_init, a_emitter_init, a_receiver_init, distance_init, ax)

plt.show()