use crate::physics::orbit::MERCURY_RADIUS;

fn shortest_lon_diff(a: f64, b: f64) -> f64 {
    let mut d = a - b;
    if d > 180.0 {
        d -= 360.0;
    }
    if d < -180.0 {
        d += 360.0;
    }
    d
}

fn basin_depth(lat: f64, lon: f64) -> f64 {
    const BASINS: &[(f64, f64, f64, f64)] = &[
        (30.5, 170.0, 12.0, -2000.0),
        (-33.0, -67.0, 5.5, -1500.0),
        (27.0, -57.0, 2.5, -1200.0),
        (-20.0, -164.0, 5.0, -1500.0),
        (-16.0, -165.0, 4.0, -1000.0),
        (-4.0, -151.0, 3.5, -900.0),
    ];
    let mut depth = 0.0_f64;
    for &(clat, clon, radius_deg, max_d) in BASINS {
        let dlat = (lat - clat) / radius_deg;
        let dlon = shortest_lon_diff(lon, clon) / radius_deg;
        let d2 = dlat * dlat + dlon * dlon;
        if d2 < 4.0 {
            let basin_profile = max_d * (-d2 * 0.8).exp();
            let rim = -max_d * 0.3 * (-((d2.sqrt() - 1.0).powi(2)) * 3.0).exp();
            depth = depth.min(basin_profile + rim);
        }
    }
    depth
}

fn crater_rim_height(lat: f64, lon: f64) -> f64 {
    const CRATERS: &[(f64, f64, f64, f64)] = &[
        (58.4, -30.0, 1.5, 1200.0),
        (-51.0, 180.0, 1.0, 800.0),
        (12.0, 137.0, 2.0, 900.0),
        (44.0, -134.0, 1.5, 700.0),
        (-3.0, -14.0, 0.8, 600.0),
        (7.0, -167.0, 1.0, 1000.0),
    ];
    let mut elev = 0.0_f64;
    for &(clat, clon, radius_deg, rim_h) in CRATERS {
        let dlat = (lat - clat) / radius_deg;
        let dlon = shortest_lon_diff(lon, clon) / radius_deg;
        let d2 = dlat * dlat + dlon * dlon;
        let rim = rim_h * (-((d2.sqrt() - 1.0).powi(2)) * 4.0).exp();
        elev = elev.max(rim);
    }
    elev
}

fn scarp_height(lat: f64, lon: f64) -> f64 {
    const SCARPS: &[(f64, f64, f64, f64, f64, f64)] = &[
        (-53.0, -38.0, 5.0, 0.5, 0.3, 2000.0),
        (-3.0, -100.0, 4.5, 0.5, 1.2, 1000.0),
        (26.0, -65.0, 6.0, 0.6, 0.8, 3000.0),
        (50.0, -32.0, 3.0, 0.4, 0.5, 1500.0),
    ];
    let mut h = 0.0_f64;
    for &(clat, clon, len, width, az, height) in SCARPS {
        let dlat = lat - clat;
        let dlon = shortest_lon_diff(lon, clon);
        let along = dlat * az.cos() + dlon * az.sin();
        let perp = -dlat * az.sin() + dlon * az.cos();
        if along.abs() < len / 2.0 {
            let scarp = height * (-(perp / width).powi(2)).exp();
            h = h.max(scarp);
        }
    }
    h
}

fn smooth_plains_offset(lat: f64, lon: f64) -> f64 {
    const PLAINS: &[(f64, f64, f64, f64)] = &[
        (65.0, 30.0, 20.0, -500.0),
        (70.0, -30.0, 15.0, -400.0),
        (30.5, 170.0, 8.0, -300.0),
    ];
    let mut offset = 0.0_f64;
    for &(clat, clon, radius_deg, depth) in PLAINS {
        let dlat = (lat - clat) / radius_deg;
        let dlon = shortest_lon_diff(lon, clon) / radius_deg;
        let d2 = dlat * dlat + dlon * dlon;
        let contribution = depth * (-d2 * 0.5).exp();
        offset = offset.min(contribution);
    }
    offset
}

pub fn mercury_elevation(lat: f64, lon: f64) -> f64 {
    let basin = basin_depth(lat, lon);
    let crater = crater_rim_height(lat, lon);
    let scarp = scarp_height(lat, lon);
    let plains = smooth_plains_offset(lat, lon);
    let base = 0.0;
    let elev = base + basin + crater + scarp + plains;
    elev.clamp(crate::LOWEST_POINT_M, crate::HIGHEST_POINT_M)
}

pub struct Heightmap {
    pub resolution: u32,
    pub data: Vec<f64>,
    pub min_elevation_m: f64,
    pub max_elevation_m: f64,
}

impl Heightmap {
    pub fn new(resolution: u32) -> Self {
        let size = (resolution * resolution) as usize;
        Self {
            resolution,
            data: vec![0.0; size],
            min_elevation_m: crate::LOWEST_POINT_M,
            max_elevation_m: crate::HIGHEST_POINT_M,
        }
    }

    pub fn sample(&self, lat_deg: f64, lon_deg: f64) -> f64 {
        let u = ((lon_deg + 180.0) / 360.0).clamp(0.0, 1.0);
        let v = ((lat_deg + 90.0) / 180.0).clamp(0.0, 1.0);
        let x = (u * (self.resolution - 1) as f64) as usize;
        let y = (v * (self.resolution - 1) as f64) as usize;
        let idx = y * self.resolution as usize + x;
        self.data.get(idx).copied().unwrap_or(0.0)
    }

    pub fn radius_at(&self, lat_deg: f64, lon_deg: f64) -> f64 {
        MERCURY_RADIUS + self.sample(lat_deg, lon_deg)
    }

    pub fn generate_mercury_topography(&mut self) {
        for y in 0..self.resolution {
            for x in 0..self.resolution {
                let lon = (x as f64 / (self.resolution - 1) as f64) * 360.0 - 180.0;
                let lat = (y as f64 / (self.resolution - 1) as f64) * 180.0 - 90.0;
                let idx = (y * self.resolution + x) as usize;
                self.data[idx] = mercury_elevation(lat, lon);
            }
        }
    }

    pub fn generate_procedural(
        &mut self,
        seed: u64,
        octaves: u32,
        persistence: f64,
        lacunarity: f64,
    ) {
        let mut rng_state = seed;
        for y in 0..self.resolution {
            for x in 0..self.resolution {
                let mut amplitude = 1.0;
                let mut frequency = 1.0;
                let mut value = 0.0;
                for _ in 0..octaves {
                    rng_state = rng_state
                        .wrapping_mul(6364136223846793005)
                        .wrapping_add(frequency as u64);
                    let noise = ((rng_state >> 33) as f64 / u32::MAX as f64) * 2.0 - 1.0;
                    value += noise * amplitude;
                    amplitude *= persistence;
                    frequency *= lacunarity;
                }
                let elevation = value.clamp(-1.0, 1.0);
                let mapped = if elevation < 0.0 {
                    elevation * self.min_elevation_m.abs()
                } else {
                    elevation * self.max_elevation_m
                };
                let idx = (y * self.resolution + x) as usize;
                self.data[idx] = mapped;
            }
        }
    }
}
