fn cloudtopvariation(lat: f64, lon: f64) -> f64 {
    let bandamplitude = 3000.0;
    let bandfreq = 3.0;
    let latrad = lat.to_radians();
    let lonrad = lon.to_radians();
    let bands = bandamplitude * (bandfreq * latrad).cos();
    let vortexlat = -22.0f64;
    let vortexlon = 296.0f64;
    let dlat = (lat - vortexlat) / 8.0;
    let dlon = shortestlondiff(lon, vortexlon) / 15.0;
    let vortexd2 = dlat * dlat + dlon * dlon;
    let vortex = 5000.0 * (-vortexd2 * 0.5).exp();
    let smallvortexlat = -54.0;
    let smallvortexlon = 340.0;
    let dlat2 = (lat - smallvortexlat) / 5.0;
    let dlon2 = shortestlondiff(lon, smallvortexlon) / 8.0;
    let smallvortex = 3000.0 * (-(dlat2 * dlat2 + dlon2 * dlon2) * 0.5).exp();
    let polarcap = if lat.abs() > 70.0 {
        1500.0 * ((lat.abs() - 70.0) / 20.0)
    } else {
        0.0
    };
    let zonalwave = 500.0 * (2.0 * lonrad).sin() * (2.0 * latrad).cos();
    bands + vortex + smallvortex + polarcap + zonalwave
}

fn shortestlondiff(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
}

pub fn neptunecloudtop(lat: f64, lon: f64) -> f64 {
    cloudtopvariation(lat, lon)
}

pub struct Heightmap {
    pub resolution: u32,
    pub data: Vec<f64>,
    pub minelevationm: f64,
    pub maxelevationm: f64,
}

impl Heightmap {
    pub fn new(resolution: u32) -> Self {
        let size = (resolution * resolution) as usize;
        Self {
            resolution,
            data: vec![0.0; size],
            minelevationm: -2000.0,
            maxelevationm: 8000.0,
        }
    }

    pub fn sample(&self, latdeg: f64, londeg: f64) -> f64 {
        let u = ((londeg + 180.0) / 360.0).clamp(0.0, 1.0);
        let v = ((latdeg + 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 radiusat(&self, latdeg: f64, londeg: f64) -> f64 {
        crate::NEPTUNEEQUATORIALRADIUS + self.sample(latdeg, londeg)
    }

    pub fn generateneptunecloudtopography(&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] = neptunecloudtop(lat, lon);
            }
        }
    }

    pub fn generateprocedural(
        &mut self,
        seed: u64,
        octaves: u32,
        persistence: f64,
        lacunarity: f64,
    ) {
        let mut rngstate = 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 octave in 0..octaves {
                    rngstate = rngstate
                        .wrapping_mul(6364136223846793005)
                        .wrapping_add(frequency as u64 + octave as u64);
                    let noise = ((rngstate >> 33) as f64 / u32::MAX as f64) * 2.0 - 1.0;
                    value += noise * amplitude;
                    amplitude *= persistence;
                    frequency *= lacunarity;
                }
                let range = self.maxelevationm - self.minelevationm;
                let idx = (y * self.resolution + x) as usize;
                self.data[idx] = self.minelevationm + (value * 0.5 + 0.5) * range;
            }
        }
    }
}
