wallswitch 0.59.6

use crate::{WallSwitchError, WallSwitchResult, get_random_integer};
use clap::ValueEnum;
use image::RgbImage;
use serde::{Deserialize, Serialize};
use std::{io::Error, path::Path, thread};

/// Represents the supported procedural background overlay effects.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize, ValueEnum)]
#[serde(rename_all = "lowercase")]
#[derive(Default)]
pub enum ProceduralEffect {
    #[value(name = "none")]
    #[default]
    None,
    #[value(name = "fractal")]
    JuliaFractal,
    #[value(name = "star")]
    Starfield,
    #[value(name = "aurora")]
    CosmicAurora,
    #[value(name = "random")]
    Random,
}

impl ProceduralEffect {
    /// Resolves the effect to a concrete rendering variant (resolving Random if selected).
    pub fn resolve(self) -> Self {
        match self {
            Self::Random => match get_random_integer(0, 2) {
                0 => Self::JuliaFractal,
                1 => Self::Starfield,
                _ => Self::CosmicAurora,
            },
            concrete => concrete,
        }
    }
}

// ==============================================================================
// FRACTAL GENERATION ENGINE
// ==============================================================================

/// Julia Fractal Generator and Image Post-Processor.
///
/// # Mathematical Concepts
///
/// ## 1. The Julia Set
/// The Julia set for a quadratic polynomial is defined by the behavior of the feedback loop:
/// `z_{n+1} = z_n^2 + c`
/// Where both `z` and `c` are complex numbers. For a given constant `c`, we iterate the function
/// starting at coordinates mapped to `z_0`. If the magnitude of `z` escapes to infinity (exceeds 2.0),
/// the pixel is colored based on the number of iterations it took to escape. If it remains bounded,
/// it belongs to the Julia set.
///
/// ## 2. Viewport Aspect Ratio Mapping and Random Zoom
/// To prevent the fractal from stretching or compressing on non-square screens (e.g., 16:9 or 21:9),
/// we map pixels to the complex plane using the minimum dimension (`min_dim`) as the scaling reference.
/// We also center the origin `(0,0)` at the exact midpoint of the image.
/// By applying a randomized zoom scale, we introduce varying levels of macro and micro-detail depth.
///
/// ## 3. Trigonometric 360-Degree Rotation
/// To rotate the fractal by a random angle `theta` around the origin without degrading performance,
/// we compute the rotation factors (sine and cosine) exactly once before entering the pixel loop:
/// `rx = cx * cos(theta) - cy * sin(theta)`
/// `ry = cx * sin(theta) + cy * cos(theta)`
///
/// ## 4. Continuous Potential (Smooth Coloring)
/// Integer iteration counts result in color banding artifacts (hard edges between color segments).
/// We solve this by applying a continuous potential calculation that factors in the escape magnitude:
/// `smooth_i = i + 1.0 - ln(ln(|z|)) / ln(2)`
/// This interpolates colors continuously, delivering smooth, professional-grade gradients.
///
/// ## 5. Contrast-Preserving Dynamic Halo Blending
/// To prevent the background image from losing its original contrast and intensity, we avoid flat
/// blending. Instead, blending factors are driven dynamically by the local fractal intensity `t`:
/// * When `t = 0` (outside the fractal), the background image is left completely untouched (100% intensity).
/// * When `t > 0` (inside fractal threads), a local shadow factor darkens the background directly beneath
///   the active pixels. This acts as a drop shadow, allowing vibrant neon fractal colors to stand out
///   clearly, even on pure white or busy backgrounds.
///
/// ## 6. Soft Vignette Shading
/// To seamlessly blend the edges and add depth, we apply a vignette filter.
/// The distance from the pixel to the center is calculated and mapped to a darkening multiplier,
/// creating a soft, professional edge fade.
pub struct FractalGenerator {
    /// Real component of the Julia set constant (c)
    pub c_re: f32,
    /// Imaginary component of the Julia set constant (c)
    pub c_im: f32,
    /// Maximum recursion iterations before determining convergence
    pub max_iterations: u32,
    /// RGB multipliers to tint the fractal intensity
    pub color_palette: [f32; 3],
    /// Randomized scale mapping determining the camera zoom level
    pub zoom: f32,
    /// Precalculated cosine of the random rotation angle
    pub cos_angle: f32,
    /// Precalculated sine of the random rotation angle
    pub sin_angle: f32,
}

impl Default for FractalGenerator {
    fn default() -> Self {
        Self {
            c_re: -0.7,
            c_im: 0.27015,
            max_iterations: 255,
            color_palette: [0.0, 1.0, 1.0], // Default Neon Cyan
            zoom: 3.0,
            cos_angle: 1.0, // 0 degrees rotation
            sin_angle: 0.0,
        }
    }
}

impl FractalGenerator {
    /// Creates a new customized Julia fractal generator.
    pub fn new(
        c_re: f32,
        c_im: f32,
        max_iterations: u32,
        color_palette: [f32; 3],
        zoom: f32,
        angle_degrees: f32,
    ) -> Self {
        let radians = angle_degrees.to_radians();
        Self {
            c_re,
            c_im,
            max_iterations,
            color_palette,
            zoom,
            cos_angle: radians.cos(),
            sin_angle: radians.sin(),
        }
    }

    /// Creates a FractalGenerator with a random aesthetic Julia constant,
    /// a randomized color palette, randomized zoom, and a randomized 360-degree rotation.
    pub fn random() -> Self {
        let presets = [
            (-0.7, 0.27015),     // Classic dendrite
            (-0.4, 0.6),         // Classic cloud swirls
            (-0.8, 0.156),       // Detailed spirals
            (-0.7269, 0.1889),   // Lace structures
            (-0.75, 0.11),       // Feathery branches
            (-0.1, 0.651),       // Cosmic dust style
            (-0.70176, -0.3842), // Dragon-like curves (San Marco fractal boundary)
            (0.355, 0.355),      // Spiral galaxy arms
            (-0.4, -0.59),       // Swirling vortexes
            (-0.54, 0.54),       // Ornamental lace borders
            (-0.74543, 0.11301), // Dense filigree patterns
            (-0.835, -0.2321),   // Lightning rods
            (-0.77269, 0.12428), // Coral reefs
            (-0.51251, 0.5213),  // Fine lace filaments
            (0.4, 0.4),          // Symmetric stellar crowns (fine dust)
            (-0.55, 0.55),       // Intricate leaf outlines
            (-0.624, 0.435),     // Crystalline snowflake patterns
            (-0.162, 1.04),      // Towering minarets
            (-0.12, 0.85),       // Flowing plasma plumes
            (-0.742, 0.1345),    // Intricate branching nodes
            (-0.391, -0.587),    // Swirling storm clouds
            (0.0, 0.8),          // Classic symmetric dendritic structure
            (-0.73, 0.21),       // Feathery dendritic lace
            (-0.81, 0.2),        // Spiral galaxy filaments
            (-0.68, 0.34),       // Delicate coral spirals
            (-0.11, 0.83),       // Plasma tendrils
            (-0.76, 0.08),       // Lightning tree branches
            (-0.72, 0.22),       // Dendritic pine branches
        ];

        let palettes = [
            [1.0, 0.0, 0.8], // Hot Pink / Cyberpunk Magenta
            [0.0, 1.0, 1.0], // Electric Cyan
            [1.0, 0.6, 0.0], // Vivid Orange / Gold
            [0.0, 1.0, 0.2], // Laser Green
            [0.6, 0.0, 1.0], // Deep Neon Violet
            [1.0, 0.1, 0.1], // Vibrant Crimson Red
            [1.0, 1.0, 0.0], // Radioactive Yellow
            [0.0, 0.4, 1.0], // Cobalt Blue
            [0.5, 1.0, 0.0], // Lime Glow
            [1.0, 0.0, 0.4], // Electric Rose
            [0.0, 1.0, 0.6], // Mint Neon
            [1.0, 0.4, 0.4], // Soft Coral Flame
            [0.9, 0.9, 1.0], // Bright Starlight White
            [1.0, 0.8, 0.0], // Amber Glow
            [0.4, 0.0, 0.8], // Electric Indigo
        ];

        let c_idx = get_random_integer(0, (presets.len() - 1) as u64) as usize;
        let p_idx = get_random_integer(0, (palettes.len() - 1) as u64) as usize;

        let zoom = get_random_integer(250, 400) as f32 / 100.0;
        let angle_degrees = get_random_integer(0, 359) as f32;
        let radians = angle_degrees.to_radians();

        let (c_re, c_im) = presets[c_idx];
        let color_palette = palettes[p_idx];

        Self {
            c_re,
            c_im,
            max_iterations: 255,
            color_palette,
            zoom,
            cos_angle: radians.cos(),
            sin_angle: radians.sin(),
        }
    }

    /// Dynamically determines the optimal 45-degree rotation step and zoom factor
    /// to maximize the visual area occupied by the fractal within the target aspect ratio.
    pub fn optimize_fit(&mut self, width: u32, height: u32) {
        let w_f = width as f32;
        let h_f = height as f32;
        let min_dim = w_f.min(h_f);

        // 1. Calculate the theoretical bounding radius of the Julia set
        let c_abs = (self.c_re * self.c_re + self.c_im * self.c_im).sqrt();
        let r_bound = (1.0 + (1.0 + 4.0 * c_abs).sqrt()) / 2.0;
        let search_limit = r_bound * 1.2;

        // 2. Perform a single lightweight pre-scan in unrotated complex space
        let steps = 40; // 40x40 grid is highly efficient and sufficient for fitting
        let inv_steps_minus_1 = 1.0 / (steps - 1) as f32;
        let range = 2.0 * search_limit;
        let scan_iterations = self.max_iterations.min(60); // Lower limit for pre-scan speed

        let mut active_points = Vec::with_capacity(steps * steps / 2);

        for step_y in 0..steps {
            let ry = -search_limit + (step_y as f32 * inv_steps_minus_1) * range;
            for step_x in 0..steps {
                let rx = -search_limit + (step_x as f32 * inv_steps_minus_1) * range;

                let mut z_re = rx;
                let mut z_im = ry;
                let mut i = 0;

                while i < scan_iterations {
                    let re2 = z_re * z_re;
                    let im2 = z_im * z_im;
                    if re2 + im2 > 4.0 {
                        break;
                    }
                    z_im = 2.0 * z_re * z_im + self.c_im;
                    z_re = re2 - im2 + self.c_re;
                    i += 1;
                }

                // Collect points that lie on the active boundary of the fractal
                if i > 3 && i < scan_iterations {
                    active_points.push((rx, ry));
                }
            }
        }

        // 3. Find the rotation step (from 0 to 315 deg) that maximizes visual coverage
        if !active_points.is_empty() {
            let mut best_zoom = f32::MAX;
            let mut best_cos = self.cos_angle;
            let mut best_sin = self.sin_angle;

            // Check 8 structural rotation angles (45-degree increments)
            for angle_step in 0..8 {
                let angle_deg = (angle_step * 45) as f32;
                let rad = angle_deg.to_radians();
                let cos_t = rad.cos();
                let sin_t = rad.sin();

                let mut max_cx_abs = 0.0_f32;
                let mut max_cy_abs = 0.0_f32;

                // Project pre-scanned points back to screen space using inverse rotation
                for &(rx, ry) in &active_points {
                    let cx = rx * cos_t + ry * sin_t;
                    let cy = -rx * sin_t + ry * cos_t;

                    max_cx_abs = max_cx_abs.max(cx.abs());
                    max_cy_abs = max_cy_abs.max(cy.abs());
                }

                // Find zoom bounds for both dimensions
                let zoom_x = 2.0 * max_cx_abs * min_dim / w_f;
                let zoom_y = 2.0 * max_cy_abs * min_dim / h_f;
                let required_zoom = zoom_x.max(zoom_y);

                // A smaller zoom scale means the camera can get closer, maximizing screen coverage
                if required_zoom < best_zoom {
                    best_zoom = required_zoom;
                    best_cos = cos_t;
                    best_sin = sin_t;
                }
            }

            self.zoom = best_zoom * 1.10; // 10% padding margin to prevent edge clipping
            self.cos_angle = best_cos;
            self.sin_angle = best_sin;
        } else {
            // Fallback to safe default bounds if no active structure is scanned
            self.zoom = 2.0 * r_bound * 1.10;
        }
    }

    /// Applies the Julia fractal and vignette blending directly on an in-memory RgbImage in parallel.
    pub fn apply_effect_in_memory(&mut self, rgb_img: &mut RgbImage) {
        let (width, height) = rgb_img.dimensions();

        // 2. Dynamically optimize zoom and rotation so that the structure fits boundaries
        self.optimize_fit(width, height);

        let w_f = width as f32;
        let h_f = height as f32;
        let min_dim = w_f.min(h_f);

        let scale = self.zoom / min_dim;

        let cos_angle = self.cos_angle;
        let sin_angle = self.sin_angle;
        let c_re = self.c_re;
        let c_im = self.c_im;
        let max_iterations = self.max_iterations;
        let color_palette = self.color_palette;

        // Partition flat pixel data into a structured vector of mutable row references (DRY Helper)
        let (mut rows, width_usize) = partition_rows(rgb_img);

        let cores = thread::available_parallelism()
            .map(|n| n.get())
            .unwrap_or(4);
        let chunk_size = (rows.len() / cores).max(1);

        // Precalculate coordinate scaling & rotation steps
        let cx_off = w_f / 2.0;
        let cy_off = h_f / 2.0;

        let dx_re = scale * cos_angle;
        let dx_im = scale * sin_angle;
        let dy_re = -scale * sin_angle;
        let dy_im = scale * cos_angle;

        let start_re = -cx_off * dx_re - cy_off * dy_re;
        let start_im = -cx_off * dx_im - cy_off * dy_im;

        // Precalculate vignette and math constants
        let inv_half_w = 2.0 / w_f;
        let inv_half_h = 2.0 / h_f;

        let inv_ln_2 = std::f32::consts::LOG2_E;

        thread::scope(|scope| {
            for chunk in rows.chunks_mut(chunk_size) {
                scope.spawn(|| {
                    for (y, row_data) in chunk.iter_mut() {
                        let y_f = *y as f32;

                        let rx_row = start_re + y_f * dy_re;
                        let ry_row = start_im + y_f * dy_im;

                        let dy_vignette = y_f * inv_half_h - 1.0;
                        let dy_vignette_sq = dy_vignette * dy_vignette;

                        for x in 0..width_usize {
                            let x_f = x as f32;

                            // Direct coordinate mapping with precalculated rotation & scale steps
                            let rx = rx_row + x_f * dx_re;
                            let ry = ry_row + x_f * dx_im;

                            let mut z_re = rx;
                            let mut z_im = ry;

                            let mut i = 0;
                            while i < max_iterations {
                                let re2 = z_re * z_re;
                                let im2 = z_im * z_im;
                                if re2 + im2 > 4.0 {
                                    break;
                                }
                                z_im = 2.0 * z_re * z_im + c_im;
                                z_re = re2 - im2 + c_re;
                                i += 1;
                            }

                            // Continuous potential smooth coloring formula
                            let t = if i < max_iterations {
                                let mag2 = z_re * z_re + z_im * z_im;
                                if mag2 > 4.0 {
                                    let log_zn = mag2.ln() * 0.5;
                                    let nu = (log_zn * inv_ln_2).ln() * inv_ln_2;
                                    let smooth_i = (i as f32 + 1.0 - nu).max(0.0);
                                    (smooth_i / max_iterations as f32).clamp(0.0, 1.0)
                                } else {
                                    i as f32 / max_iterations as f32
                                }
                            } else {
                                1.0
                            };

                            let idx = x * 3;
                            let original_r = row_data[idx];
                            let original_g = row_data[idx + 1];
                            let original_b = row_data[idx + 2];

                            // Contrast-Preserving Dynamic Halo Blending (drop shadow under threads)
                            let shadow_factor = 1.0 - (t * 0.5);
                            let background_r = original_r as f32 * shadow_factor;
                            let background_g = original_g as f32 * shadow_factor;
                            let background_b = original_b as f32 * shadow_factor;

                            let r_fractal = color_palette[0] * t * 255.0;
                            let g_fractal = color_palette[1] * t * 255.0;
                            let b_fractal = color_palette[2] * t * 255.0;

                            let alpha = t.sqrt() * 0.8;
                            let blended_r = (background_r * (1.0 - alpha)) + (r_fractal * alpha);
                            let blended_g = (background_g * (1.0 - alpha)) + (g_fractal * alpha);
                            let blended_b = (background_b * (1.0 - alpha)) + (b_fractal * alpha);

                            // Soft Vignette Calculation
                            let dx_vignette = x_f * inv_half_w - 1.0;
                            let dist = (dx_vignette * dx_vignette + dy_vignette_sq).sqrt();
                            let vignette = (1.0 - dist * 0.4).clamp(0.1, 1.0);

                            row_data[idx] = (blended_r * vignette).clamp(0.0, 255.0) as u8;
                            row_data[idx + 1] = (blended_g * vignette).clamp(0.0, 255.0) as u8;
                            row_data[idx + 2] = (blended_b * vignette).clamp(0.0, 255.0) as u8;
                        }
                    }
                });
            }
        });
    }

    /// Reads an input image, applies the Julia fractal in parallel, and saves to the output path.
    pub fn apply_effect<P: AsRef<Path>>(
        &mut self,
        input_path: P,
        output_path: P,
    ) -> WallSwitchResult<()> {
        let img = image::open(&input_path)
            .map_err(|e| WallSwitchError::UnableToFind(format!("Failed to open image: {e}")))?;

        let mut rgb_img = img.to_rgb8();
        self.apply_effect_in_memory(&mut rgb_img);

        rgb_img
            .save(&output_path)
            .map_err(|e| WallSwitchError::Io(Error::other(e)))?;

        Ok(())
    }
}

// ==============================================================================
// SHARED UTILITY FUNCTIONS (DRY & Harmonious Color Optimization)
// ==============================================================================

/// Partitions an RGB image buffer into mutable row segments for thread-safe parallel processing.
///
/// Rather than repeating buffer striding code across different generators, this unified function
/// converts a flat 1D byte slice into distinct, row-by-row mutable chunks based on a 3-byte RGB layout.
fn partition_rows(rgb_img: &mut RgbImage) -> (Vec<(usize, &mut [u8])>, usize) {
    let (width, _) = rgb_img.dimensions();
    let width_usize = width as usize;
    let row_stride = width_usize * 3;
    let pixels_buffer = rgb_img.as_mut();

    let rows: Vec<(usize, &mut [u8])> = pixels_buffer
        .chunks_exact_mut(row_stride)
        .enumerate()
        .collect();

    (rows, width_usize)
}

// ==============================================================================
// ALTERNATIVE PROCEDURAL OVERLAYS
// ==============================================================================

/// Helper representation of individual stars in a Starfield.
pub struct Star {
    pub x: f32,
    pub y: f32,
    pub radius: f32,
    pub color: [f32; 3],
    pub intensity: f32,
}

/// Cyberpunk Starfield / Bokeh effect generator.
pub struct StarfieldGenerator {
    pub stars: Vec<Star>,
}

impl StarfieldGenerator {
    /// Generates a randomized list of glowing stars based on target monitor dimension limits.
    pub fn new(count: usize, width: u32, height: u32) -> Self {
        let mut stars = Vec::with_capacity(count);

        let palettes = [
            [1.0, 1.0, 1.0], // White
            [0.6, 0.8, 1.0], // Electric Ice Blue
            [1.0, 0.8, 0.4], // Cosmic Gold
            [1.0, 0.4, 0.8], // Ultraviolet Pink
        ];

        for _ in 0..count {
            let x = get_random_integer(0, width as u64) as f32;
            let y = get_random_integer(0, height as u64) as f32;
            let radius = get_random_integer(5, 45) as f32;
            let intensity = get_random_integer(30, 95) as f32 / 100.0;

            let p_idx = get_random_integer(0, (palettes.len() - 1) as u64) as usize;
            let color = palettes[p_idx];

            stars.push(Star {
                x,
                y,
                radius,
                color,
                intensity,
            });
        }

        Self { stars }
    }

    /// Appends smooth glowing star circles directly onto the image in parallel.
    pub fn apply_effect_in_memory(&self, rgb_img: &mut RgbImage) {
        // Using a soft pastel blue-gray default color for blending to simplify code
        let contrast_color = [0.64, 0.75, 0.85];

        let (mut rows, width_usize) = partition_rows(rgb_img);

        let cores = thread::available_parallelism()
            .map(|n| n.get())
            .unwrap_or(4);
        let chunk_size = (rows.len() / cores).max(1);

        thread::scope(|scope| {
            for chunk in rows.chunks_mut(chunk_size) {
                let stars = &self.stars;
                scope.spawn(move || {
                    for (y, row_data) in chunk.iter_mut() {
                        let y_f = *y as f32;

                        // Gather active stars vertically overlapping with this row to optimize calculations
                        let mut active_stars = Vec::with_capacity(16);
                        for star in stars {
                            let dy = star.y - y_f;
                            let limit = star.radius * 2.0;
                            if dy.abs() < limit {
                                let dy_sq = dy * dy;
                                let star_radius_sq = star.radius * star.radius;
                                active_stars.push((star, dy_sq, star_radius_sq));
                            }
                        }

                        for x in 0..width_usize {
                            let x_f = x as f32;

                            let mut r_contrib = 0.0;
                            let mut g_contrib = 0.0;
                            let mut b_contrib = 0.0;
                            let mut total_alpha = 0.0;

                            for &(star, dy_sq, star_radius_sq) in &active_stars {
                                let dx = star.x - x_f;
                                let dist_sq = dx * dx + dy_sq;

                                if dist_sq < star_radius_sq * 4.0 {
                                    let factor = (-dist_sq / (2.0 * star_radius_sq)).exp();
                                    let alpha = factor * star.intensity;

                                    let r_star =
                                        (star.color[0] * 0.25 + contrast_color[0] * 0.75) * alpha;
                                    let g_star =
                                        (star.color[1] * 0.25 + contrast_color[1] * 0.75) * alpha;
                                    let b_star =
                                        (star.color[2] * 0.25 + contrast_color[2] * 0.75) * alpha;

                                    r_contrib += r_star;
                                    g_contrib += g_star;
                                    b_contrib += b_star;
                                    total_alpha += alpha;
                                }
                            }

                            if total_alpha > 0.001 {
                                let idx = x * 3;
                                let original_r = row_data[idx] as f32;
                                let original_g = row_data[idx + 1] as f32;
                                let original_b = row_data[idx + 2] as f32;

                                let alpha_clamp = total_alpha.min(0.95);

                                let blended_r = (original_r * (1.0 - alpha_clamp))
                                    + (r_contrib * 255.0 / total_alpha * alpha_clamp);
                                let blended_g = (original_g * (1.0 - alpha_clamp))
                                    + (g_contrib * 255.0 / total_alpha * alpha_clamp);
                                let blended_b = (original_b * (1.0 - alpha_clamp))
                                    + (b_contrib * 255.0 / total_alpha * alpha_clamp);

                                row_data[idx] = blended_r.clamp(0.0, 255.0) as u8;
                                row_data[idx + 1] = blended_g.clamp(0.0, 255.0) as u8;
                                row_data[idx + 2] = blended_b.clamp(0.0, 255.0) as u8;
                            }
                        }
                    }
                });
            }
        });
    }
}

/// Cosmic Aurora overlay generator.
pub struct AuroraGenerator {
    pub color_palette: [f32; 3],
    pub density: f32,
}

impl AuroraGenerator {
    pub fn random() -> Self {
        let palettes = [
            [0.2, 1.0, 0.5], // Emerald Green / Aurora Classic
            [0.6, 0.0, 1.0], // Deep Cosmic Violet
            [0.0, 0.8, 1.0], // Neon Ice Teal
            [1.0, 0.0, 0.6], // Solar Flare Pink
        ];
        let idx = get_random_integer(0, (palettes.len() - 1) as u64) as usize;
        let density = get_random_integer(4, 8) as f32;

        Self {
            color_palette: palettes[idx],
            density,
        }
    }

    /// Appends smooth cosmic waves directly onto the image in parallel.
    pub fn apply_effect_in_memory(&self, rgb_img: &mut RgbImage) {
        let (width, height) = rgb_img.dimensions();
        let w_f = width as f32;
        let h_f = height as f32;

        // Using the generator's selected color palette directly to simplify code
        let contrast_color = self.color_palette;

        let (mut rows, width_usize) = partition_rows(rgb_img);

        let cores = thread::available_parallelism()
            .map(|n| n.get())
            .unwrap_or(4);
        let chunk_size = (rows.len() / cores).max(1);

        // Precalculate scaling coefficients and row factors to optimize performance
        let inv_w = 1.0 / w_f;
        let inv_h = 1.0 / h_f;
        let density_u = self.density * 1.5 * inv_w;
        let density_v_coeff = self.density * 2.0;
        let density_w_coeff = self.density * inv_w;
        let density_w4_coeff = self.density * 1.2;

        thread::scope(|scope| {
            for chunk in rows.chunks_mut(chunk_size) {
                scope.spawn(move || {
                    for (y, row_data) in chunk.iter_mut() {
                        let y_f = *y as f32;

                        let v = y_f * inv_h;
                        let w2 = (v * density_v_coeff).cos();
                        let v_density = v * self.density;
                        let v_sq = v * v;

                        for x in 0..width_usize {
                            let x_f = x as f32;

                            let u = x_f * inv_w;

                            let w1 = (x_f * density_u).sin();
                            let w3 = (x_f * density_w_coeff + v_density).sin();
                            let w4 = ((u * u + v_sq).sqrt() * density_w4_coeff).cos();

                            let val = (w1 + w2 + w3 + w4) * 0.25;

                            // High-contrast cos-pow3 wave function to generate distinct glow filaments
                            let intensity = (val * std::f32::consts::PI).cos().abs().powf(3.0);

                            let dx = (u - 0.5) * 2.0;
                            let dy = (v - 0.5) * 2.0;
                            let edge_fade =
                                (1.0 - (dx * dx + dy * dy).sqrt() * 0.5).clamp(0.0, 1.0);

                            // Opacity increased to 0.75 to make waves clearly visible
                            let alpha = intensity * edge_fade * 0.75;

                            if alpha > 0.01 {
                                let idx = x * 3;
                                let original_r = row_data[idx] as f32;
                                let original_g = row_data[idx + 1] as f32;
                                let original_b = row_data[idx + 2] as f32;

                                let r_aurora = contrast_color[0] * 255.0;
                                let g_aurora = contrast_color[1] * 255.0;
                                let b_aurora = contrast_color[2] * 255.0;

                                row_data[idx] =
                                    ((original_r * (1.0 - alpha)) + (r_aurora * alpha)) as u8;
                                row_data[idx + 1] =
                                    ((original_g * (1.0 - alpha)) + (g_aurora * alpha)) as u8;
                                row_data[idx + 2] =
                                    ((original_b * (1.0 - alpha)) + (b_aurora * alpha)) as u8;
                            }
                        }
                    }
                });
            }
        });
    }
}